for functions to Differentials_, added multithreaded version

src/main/java/epfl/Differentials_.java

@@ -3,6 +3,8 @@ package epfl;
 | Version: April 6, 2014
 \=====================================================================*/
 
+import java.util.concurrent.atomic.AtomicInteger;
+
 /*=====================================================================
 | Philippe Thevenaz
 | EPFL/STI/IMT-LS/LIB/BM.4.137
@@ -77,6 +79,12 @@ public final static int GRADIENT_X = 6;
 public final static int GRADIENT_Y = 7;
 public final static int GRADIENT_XY = 8;
 public final static int NORMALIZE_LMAX = 9;
+public final static int LOCAL_MIN = 10;  // local minimum within radius
+public final static int LOCAL_MAX = 11;  // local maximum within radius
+public final static int LOCAL_DIFF = 12; // local maximum - local minimum within radius
+
+
+private final static int THREADS_MAX = 100;
 
 /*....................................................................
 	private variables
@@ -103,7 +111,10 @@ private static String[] operations = {
 	"Gradient X",
 	"Gradient Y",
 	"Directional derivative",
-	"Normalize by local maximum"
+	"Normalize by local maximum",
+	"Local min",
+	"Local max",
+	"Local diff"
 };
 
 /*....................................................................
@@ -147,17 +158,19 @@ public void run (
 	GenericDialog dialog = new GenericDialog("Differentials");
 	dialog.addChoice("Operation:", operations, operations[operation]);
 	dialog.addMessage("--- optional parameters ---");
-	dialog.addNumericField("Multiply d/dx by",       dir_x, 3,5,"");
-	dialog.addNumericField("Multiply d/dy by",       dir_y, 3,5,"");
-	dialog.addNumericField("Normalization radius",   max_norm_rad, 3,5,"");
+	dialog.addNumericField("Multiply_d/dx_by",       dir_x, 3,5,"");
+	dialog.addNumericField("Multiply_d/dy_by",       dir_y, 3,5,"");
+	dialog.addNumericField("Normalization_radius",   max_norm_rad, 3,5,"");
 
 	dialog.showDialog();
 	if (dialog.wasCanceled()) {
 		return;
 	}
 	final String operationString = dialog.getNextChoice();
+    System.out.println("operationString='"+operationString+"'");
+	
 	for (int choice = 0; (choice < operations.length); choice++) {
-		if (operationString == operations[choice]) {
+		if (operationString.equals(operations[choice])) {
 			operation = choice;
 			break;
 		}
@@ -165,7 +178,7 @@ public void run (
 	dir_x =        (float) dialog.getNextNumber();
 	dir_y =        (float) dialog.getNextNumber();
 	max_norm_rad = dialog.getNextNumber();
-
+    System.out.println("operation="+operation+", dir_x="+dir_x+", dir_y="+dir_y+", max_norm_rad="+max_norm_rad);
 	imp.startTiming();
 	if (1 < imp.getStackSize()) {
 		if (!(imp.getProcessor().getPixels() instanceof float[])) {
@@ -357,6 +370,156 @@ private void doIt (
 		throw new IllegalArgumentException("Float image required");
 	}
 	switch (operation) {
+	case LOCAL_DIFF:
+		{
+			final int ir = (int) Math.ceil(max_norm_rad);
+			final double r2 = max_norm_rad * max_norm_rad;
+			final float[] floatPixels =  (float[])ip.getPixels();
+			final double [] pixels = new double[floatPixels.length];
+			for (int i = 0; i < pixels.length; i++) {
+				pixels[i] = floatPixels[i];
+			}
+			final Thread[] threads = newThreadArray(THREADS_MAX);
+			final AtomicInteger ai = new AtomicInteger(0);
+			for (int y00 = 0; (y00 < height); y00++) {
+				final int y0 = y00;
+				int ymin = Math.max(y0-ir, 0);
+				int ymax = Math.min(y0+ir, height-1);
+				ai.set(0);
+				for (int ithread = 0; ithread < threads.length; ithread++) {
+					threads[ithread] = new Thread() {
+						public void run() {
+							for (int x0 = ai.getAndIncrement(); x0 < width; x0 = ai.getAndIncrement()) {
+								int k = y0*width + x0;
+								int xmin = Math.max(x0-ir, 0);
+								int xmax = Math.min(x0+ir, width-1);
+								double d = pixels[k];
+								double dmax = d;
+								double dmin = d;
+								for (int y = ymin; y <= ymax; y++) {
+									for (int x = xmin; x <= xmax; x++) {
+										int dx = x-x0;
+										int dy = y-y0;
+										if ((dx*dx+dy*dy) > r2) {
+											continue;
+										}
+										if (pixels[y*width+x] > dmax) {
+											dmax=pixels[y*width+x];
+										} else if (pixels[y*width+x] < dmin) {
+											dmin=pixels[y*width+x];
+										}
+
+									}
+								}
+								floatPixels[k] = (float) (dmax - dmin);
+							}
+						}
+					};
+				}		      
+				startAndJoin(threads);
+				stepProgressBar();
+			}
+		}
+		break;
+	
+	
+	
+		case LOCAL_MIN:
+		{
+			final int ir = (int) Math.ceil(max_norm_rad);
+			final double r2 = max_norm_rad * max_norm_rad;
+			final float[] floatPixels =  (float[])ip.getPixels();
+			final double [] pixels = new double[floatPixels.length];
+			for (int i = 0; i < pixels.length; i++) {
+				pixels[i] = floatPixels[i];
+			}
+			final Thread[] threads = newThreadArray(THREADS_MAX);
+			final AtomicInteger ai = new AtomicInteger(0);
+			for (int y00 = 0; (y00 < height); y00++) {
+				final int y0 = y00;
+				int ymin = Math.max(y0-ir, 0);
+				int ymax = Math.min(y0+ir, height-1);
+				ai.set(0);
+				for (int ithread = 0; ithread < threads.length; ithread++) {
+					threads[ithread] = new Thread() {
+						public void run() {
+							for (int x0 = ai.getAndIncrement(); x0 < width; x0 = ai.getAndIncrement()) {
+								int k = y0*width + x0;
+								int xmin = Math.max(x0-ir, 0);
+								int xmax = Math.min(x0+ir, width-1);
+								double d = pixels[k];
+								double dmin = d;
+								for (int y = ymin; y <= ymax; y++) {
+									for (int x = xmin; x <= xmax; x++) {
+										int dx = x-x0;
+										int dy = y-y0;
+										if ((dx*dx+dy*dy) > r2) {
+											continue;
+										}
+										if (pixels[y*width+x] < dmin) {
+											dmin=pixels[y*width+x];
+										}
+									}
+								}
+								floatPixels[k] = (float) dmin;
+							}
+						}
+					};
+				}		      
+				startAndJoin(threads);
+				stepProgressBar();
+			}
+		}
+		break;
+		
+		case LOCAL_MAX:
+		{
+			final int ir = (int) Math.ceil(max_norm_rad);
+			final double r2 = max_norm_rad * max_norm_rad;
+			final float[] floatPixels =  (float[])ip.getPixels();
+			final double [] pixels = new double[floatPixels.length];
+			for (int i = 0; i < pixels.length; i++) {
+				pixels[i] = floatPixels[i];
+			}
+			final Thread[] threads = newThreadArray(THREADS_MAX);
+			final AtomicInteger ai = new AtomicInteger(0);
+			for (int y00 = 0; (y00 < height); y00++) {
+				final int y0 = y00;
+				int ymin = Math.max(y0-ir, 0);
+				int ymax = Math.min(y0+ir, height-1);
+				ai.set(0);
+				for (int ithread = 0; ithread < threads.length; ithread++) {
+					threads[ithread] = new Thread() {
+						public void run() {
+							for (int x0 = ai.getAndIncrement(); x0 < width; x0 = ai.getAndIncrement()) {
+								int k = y0*width + x0;
+								int xmin = Math.max(x0-ir, 0);
+								int xmax = Math.min(x0+ir, width-1);
+								double d = pixels[k];
+								double dmax = d;
+								for (int y = ymin; y <= ymax; y++) {
+									for (int x = xmin; x <= xmax; x++) {
+										int dx = x-x0;
+										int dy = y-y0;
+										if ((dx*dx+dy*dy) > r2) {
+											continue;
+										}
+										if (pixels[y*width+x] > dmax) {
+											dmax=pixels[y*width+x];
+										}
+									}
+								}
+								floatPixels[k] = (float) dmax;
+							}
+						}
+					};
+				}		      
+				startAndJoin(threads);
+				stepProgressBar();
+			}
+		}
+		break;
+		
 		case NORMALIZE_LMAX:
 			{
 				int ir = (int) Math.ceil(max_norm_rad);
@@ -771,6 +934,9 @@ private void setupProgressBar (
 	completed = 0;
 	lastTime = System.currentTimeMillis();
 	switch (operation) {
+		case LOCAL_MIN:
+		case LOCAL_MAX:
+		case LOCAL_DIFF:
 		case NORMALIZE_LMAX:
 			processDuration = stackSize * (width + 1 * height);
 			break;
@@ -849,4 +1015,36 @@ private void symmetricFirMirrorOnBounds (
 	}
 } /* end symmetricFirMirrorOnBounds */
 
+/* Create a Thread[] array as large as the number of processors available.
+ * From Stephan Preibisch's Multithreading.java class. See:
+ * http://repo.or.cz/w/trakem2.git?a=blob;f=mpi/fruitfly/general/MultiThreading.java;hb=HEAD
+ */
+public static Thread[] newThreadArray(int maxCPUs) { // USED in lwir
+	int n_cpus = Runtime.getRuntime().availableProcessors();
+	if (n_cpus>maxCPUs)n_cpus=maxCPUs;
+	return new Thread[n_cpus];
+}
+/* Start all given threads and wait on each of them until all are done.
+ * From Stephan Preibisch's Multithreading.java class. See:
+ * http://repo.or.cz/w/trakem2.git?a=blob;f=mpi/fruitfly/general/MultiThreading.java;hb=HEAD
+ */
+public static void startAndJoin(Thread[] threads) // USED in lwir
+{
+	for (int ithread = 0; ithread < threads.length; ++ithread)
+	{
+		threads[ithread].setPriority(Thread.NORM_PRIORITY);
+		threads[ithread].start();
+	}
+
+	try
+	{
+		for (int ithread = 0; ithread < threads.length; ++ithread)
+			threads[ithread].join();
+	} catch (InterruptedException ie)
+	{
+		throw new RuntimeException(ie);
+	}
+}
+
+
 } /* end class Differentials_ */