Automatic matching of multiple scene sets

src/main/java/com/elphel/imagej/calibration/CalibrationFileManagement.java

@@ -8,6 +8,7 @@ import javax.swing.JFileChooser;
 import javax.swing.filechooser.FileFilter;
 
 import ij.IJ;
+import ij.Prefs;
 import ij.io.OpenDialog;
 
 
@@ -187,16 +188,23 @@ public class CalibrationFileManagement {
 	  }
 	  
 	  public static void saveStringToFile (String path,String data, boolean append){
+		  int dir_sep = path.lastIndexOf(Prefs.getFileSeparator());
+		  if (dir_sep >=0) { // create directory if it does not exist
+			  File dir = new File(path.substring(0, dir_sep));
+			  if (!dir.exists()) {
+				  dir.mkdirs();
+			  }
+		  }
 		  BufferedWriter writer = null;
 		  try  {
 			  writer = new BufferedWriter( new FileWriter( path, append));
 			  writer.write( data);
 
 		  }  catch ( IOException e)  {
-				String msg = e.getMessage();
-				if (msg==null || msg.equals(""))  msg = ""+e;
-				IJ.showMessage("Error",msg);
-				throw new IllegalArgumentException (msg);
+			  String msg = e.getMessage();
+			  if (msg==null || msg.equals(""))  msg = ""+e;
+			  IJ.showMessage("Error",msg);
+			  throw new IllegalArgumentException (msg);
 		  }
 		  finally  {
 			  try {

src/main/java/com/elphel/imagej/orthomosaic/OrthoMap.java

@@ -191,6 +191,13 @@ public class OrthoMap implements Comparable <OrthoMap>, Serializable{
 			PairwiseOrthoMatch match) {
 		pairwise_matches.put(name, match);
 	}
+
+	public void unsetMatch(
+			String name) {
+		pairwise_matches.remove(name);
+	}
+
+	
 	public PairwiseOrthoMatch getMatch(String name) {
 		return pairwise_matches.get(name);
 	}

src/main/java/com/elphel/imagej/orthomosaic/OrthoPairLMA.java

@@ -27,6 +27,7 @@
 package com.elphel.imagej.orthomosaic;
 
 import java.awt.Rectangle;
+import java.util.Arrays;
 import java.util.concurrent.atomic.AtomicInteger;
 
 import com.elphel.imagej.common.ShowDoubleFloatArrays;
@@ -105,6 +106,8 @@ public class OrthoPairLMA {
 			int           min_good_tiles,
 			double        max_std,      // maximal standard deviation to limit center area  
 			double        min_std_rad,  // minimal radius of the central area (if less - fail)
+			double        tile_rad,
+			int           min_tiles_rad,
 			final int     debug_level) {
 		tile_centers = centers;
 		this.width = width;
@@ -127,16 +130,17 @@ public class OrthoPairLMA {
 		origin = new double[2];
 		if (origin_center) {
 			origin = new double [] {
-					woi.x + 0.5 * woi.width  * GPUTileProcessor.DTT_SIZE,
-					woi.y + 0.5 * woi.height * GPUTileProcessor.DTT_SIZE};
+					(woi.x + 0.5 * woi.width)  * GPUTileProcessor.DTT_SIZE,
+					(woi.y + 0.5 * woi.height) * GPUTileProcessor.DTT_SIZE};
 		}
 		N = woi.width * woi.height;
 		parameters_vector = new double [] {1,0,0,1,0,0};
 		if (min_std_rad > 0) {
-			int         min_tiles = 4;
+//			int         min_tiles = 4;
 			getCenterRadius(
-					max_std,       // final double      max_std,      // maximal standard deviation to limit center area  
-					min_tiles,     // final int         min_tiles,
+					max_std,       // final double      max_std,      // maximal standard deviation to limit center area
+					tile_rad,      // min_std_rad,   // final double      min_radius,
+					min_tiles_rad, // min_tiles,     // final int         min_tiles,
 					vector_XYS,    // final double [][] vector_XYS,
 					weights_extra, // final double []   weights_extra, // null or additional weights (such as elevation-based)
 					centers);      // final double [][] centers)
@@ -251,6 +255,11 @@ public class OrthoPairLMA {
 	 * In that case try to adjust only the central area first, then increase that area
 	 * in next iterations. 
 	 * @param max_std       maximal standard deviation (average for X and Y) inside center area
+	 * @param min_radius    minimal central zone radius that has to have >= min_tiles.
+	 *                      Should be increased (with min_tiles) for high resolution images.
+	 *                      Also there can be a bush right in the center - need to handle it too.
+	 *                      Maybe for the final it should be a fraction of the minimal overlap
+	 *                      dimension? 
 	 * @param min_tiles     minimal tiles in the center zone
 	 * @param vector_XYS    2D correlation-measured X,Y, and strength
 	 * @param weights_extra null or optional additional weights of the samples
@@ -258,9 +267,10 @@ public class OrthoPairLMA {
 	 * @return maximal radius from the center (in pixels) where standard deviation of the inside samples is
 	 *         below max_std. Returns Double.POSITIVE_INFINITY if std is not reached
 	 */
-	private double getCenterRadius(
+	public double getCenterRadius(
 	        final double      max_std,      // maximal standard deviation to limit center area  
 //	        final double      min_std_rad,  // minimal radius of the central area (if less - fail)
+	        final double      min_radius,
 	        final int         min_tiles,
 			final double [][] vector_XYS,
 			final double []   weights_extra, // null or additional weights (such as elevation-based)
@@ -276,7 +286,53 @@ public class OrthoPairLMA {
 		final double [][] sy_arr =  new double [threads.length][rad_length];
 		final double [][] sy2_arr = new double [threads.length][rad_length];
 		final int    [][] sn_arr =  new int [threads.length][rad_length];
-		
+		boolean dbg = false;
+		if (dbg) {
+			String [] titles = {"fx","fy","vw","w","cent-x","cent-y","dx","dy","r_t","irt"};
+			final double [][] dbg_img = new double [titles.length][woi.width*woi.height];
+			for (int i = 0; i < dbg_img.length; i++) {
+				Arrays.fill(dbg_img[i], Double.NaN);
+			}
+			for (int ithread = 0; ithread < threads.length; ithread++) {
+				threads[ithread] = new Thread() {
+					public void run() {
+						for (int iTile = ai.getAndIncrement(); iTile < N; iTile = ai.getAndIncrement()) {
+							int tileX = iTile % woi.width + woi.x;
+							int tileY = iTile / woi.width + woi.y;
+							int aTile = tileY * width + tileX;
+							if ((vector_XYS[aTile] != null) && (centers[aTile] != null)) {
+								double w = vector_XYS[aTile][2];
+								if (weights_extra != null) w *= weights_extra[aTile]; 
+								if (Double.isNaN(w)) w = 0;
+								double dx = centers[aTile][0] - origin[0];
+								double dy = centers[aTile][1] - origin[1];
+								double r_t = Math.sqrt(dx*dx+dy*dy)/GPUTileProcessor.DTT_SIZE; // radius in tiles
+								int irt = (int) Math.round(r_t);
+								dbg_img[0][iTile] = vector_XYS[aTile][0];
+								dbg_img[1][iTile] = vector_XYS[aTile][1];
+								dbg_img[2][iTile] = vector_XYS[aTile][2];
+								dbg_img[3][iTile] = w;
+								dbg_img[4][iTile] = centers[aTile][0];
+								dbg_img[5][iTile] = centers[aTile][1];
+								dbg_img[6][iTile] = dx;
+								dbg_img[7][iTile] = dy;
+								dbg_img[8][iTile] = r_t;
+								dbg_img[9][iTile] = irt;
+							}
+						}
+					}
+				};
+			}		      
+			ImageDtt.startAndJoin(threads);
+			ai.set(0);
+			ShowDoubleFloatArrays.showArrays(
+					dbg_img,
+					woi.width,
+					woi.height,
+					true,
+					"getCenterRadius",
+					titles);
+		}
 		for (int ithread = 0; ithread < threads.length; ithread++) {
 			threads[ithread] = new Thread() {
 				public void run() {
@@ -335,6 +391,7 @@ public class OrthoPairLMA {
 		ImageDtt.startAndJoin(threads);
 		double sc0=0, scx=0,scx2=0,scy=0,scy2=0,std_prev=0, std=0;
 		int    scn=0;
+		int    scnc = 0;
 		center_radius = Double.POSITIVE_INFINITY; 
 		for (int iRad=0; iRad < rad_length; iRad++) {
 			sc0+= s0 [iRad];
@@ -343,7 +400,10 @@ public class OrthoPairLMA {
 			scy+= sy [iRad];
 			scy2+=sy2[iRad];
 			scn+= sn [iRad];
-			if ((scn > min_tiles) && (sc0 > 0)) {// for one tile gets scrt() of a small negative error
+			if (iRad <= ((int) Math.round(min_radius))) {
+				scnc = scn;
+			}
+			if ((scn > min_tiles) && (sc0 > 0)) {// for one tile gets sqcrt() of a small negative error
 				std = Math.sqrt((scx2*sc0 - scx*scx + scy2*sc0 - scy*scy)/2)/sc0;
 			}
 			if ((scn >= min_tiles) && (std >= max_std)) {
@@ -356,6 +416,9 @@ public class OrthoPairLMA {
 			}
 			std_prev = std;
 		}
+		if ((min_tiles > 0) && (min_radius > 0) && (scnc < min_tiles)) {
+			center_radius = 0;
+		}
 		return center_radius;
 	}
 	

src/main/java/com/elphel/imagej/orthomosaic/PairwiseOrthoMatch.java

@@ -11,6 +11,7 @@ public class PairwiseOrthoMatch implements Serializable {
 	public transient double [][] jtj =    new double [6][6];
 	public int zoom_lev;
 	public double rms = Double.NaN;
+	public transient int [] nxy = null; // not saved, just to communicate for logging 
 	public PairwiseOrthoMatch() {
 		
 	}
@@ -19,9 +20,10 @@ public class PairwiseOrthoMatch implements Serializable {
 			double [][] jtj,
 			double rms,
 			int zoom_lev) {
-		this.affine = affine;
-		this.jtj = jtj;
-		this.zoom_lev=zoom_lev;
+		this.affine =  affine;
+		this.jtj =     jtj;
+		this.zoom_lev= zoom_lev;
+		this.rms =     rms;
 	}
 	public PairwiseOrthoMatch clone() {
 		double [][] affine = {this.affine[0].clone(),this.affine[1].clone()};
@@ -34,6 +36,9 @@ public class PairwiseOrthoMatch implements Serializable {
 				jtj,
 				this.rms,
 				this.zoom_lev);
+		if (nxy != null) {
+			pom.nxy = nxy.clone();
+		}
 		return pom;
 	}
 	public PairwiseOrthoMatch getInverse(double [] rd) {