tested derivative for the OrthoMultiLMA

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

@@ -166,7 +166,7 @@ public class ComboMatch {
 		int num_tries_fit =           10;
 		boolean update_match =        true; // use false to save new version of data
 		boolean render_match =        false; // true;
-
+		boolean test_multi_lma =      false;
 		boolean pattern_match =       true; // false;
 		
 		boolean bounds_to_indices =   true;
@@ -202,7 +202,8 @@ public class ComboMatch {
 
 		gd.addCheckbox    ("Update match if calculated", update_match, "Will update correlation match for a pair if found.");
 		gd.addCheckbox    ("Render match", render_match, "Render a pair of matched images.");
-		
+		gd.addCheckbox    ("Test multi LMA", test_multi_lma, "Temporary debug.");
+		//
 		gd.addCheckbox    ("Pattern match", pattern_match, "Search for patterns for both images in a pair, first is primary.");
 		
 		gd.addCheckbox    ("Bounds to selected images", bounds_to_indices, "Set combo image bounds to selected images only. False - all images.");
@@ -262,7 +263,7 @@ public class ComboMatch {
 		num_tries_fit =      (int) gd.getNextNumber();
 		update_match=              gd.getNextBoolean();
 		render_match=              gd.getNextBoolean();
-		
+		test_multi_lma=            gd.getNextBoolean();
 		pattern_match=             gd.getNextBoolean();
 		
 		bounds_to_indices=         gd.getNextBoolean();
@@ -774,7 +775,7 @@ public class ComboMatch {
 					orthoMapsCollection_path); // String orthoMapsCollection_path);
 			if (!ok) return false;
 		}
-        if (process_correlation || render_match || pattern_match) {
+        if (process_correlation || render_match || pattern_match || test_multi_lma) {
 //        	int [] gpu_pair;
         	if (gpu_spair == null) {
         		ArrayList<Point> pairs_list = new ArrayList<Point>(); 
@@ -859,24 +860,8 @@ public class ComboMatch {
         	double agl_ratio = max_agl/50.0;
         	double metric_error_adj = metric_error * agl_ratio * agl_ratio; // metric_error settings is good for 50m. Increase for higher Maybe squared?
         	int initial_zoom = max_zoom_lev - 4;  // another algorithm?
-/*
-        	if (GPU_QUAD_AFFINE == null) {
-            	System.out.println("Setting up GPU");
-        		try {
-        			GPU_QUAD_AFFINE = new GpuQuad(//
-        					GPU_TILE_PROCESSOR,   // GPUTileProcessor gpuTileProcessor,
-        					gpu_max_width,        // final int        max_width,
-        					gpu_max_height,       // final int        max_height,
-        					1,                    // final int        num_colors, // normally 1?
-        					clt_parameters.gpu_debug_level);
-        		} catch (Exception e) {
-        			System.out.println("Failed to initialize GpuQuad class");
-        			// TODO Auto-generated catch block
-        			e.printStackTrace();
-        			return false;
-        		} // final int debugLevel);
-        	}
-*/        	
+        	
+        	// Here - always start with unity affine0 and affine1 from possibly matched pair 
     		double [][] affine0 = {{1,0,0},{0,1,0}}; // will always stay the same
     		double [][] affine1 = null;
 
@@ -893,6 +878,13 @@ public class ComboMatch {
 				}
         		
         	} else {
+        		if (test_multi_lma) {
+        			OrthoMultiLMA.testMultiLMA(
+        					clt_parameters,  // CLTParameters       clt_parameters,
+        					maps_collection, // OrthoMapsCollection maps_collection,
+        					gpu_pair);       // int []              indices)
+        			return true;
+        		}
         		if (process_correlation && !use_marked_image) { // match may or may not exist
         			// if match exists - ask if use it. If not - open dialog and start spiral
         			pairwiseOrthoMatch = initialPairAdjust(
@@ -1095,6 +1087,7 @@ public class ComboMatch {
 		boolean log_append =           clt_parameters.imp.pwise_log_append;
 		String log_path =              clt_parameters.imp.pwise_log_path;
 		
+		boolean     pmtch_use_affine =    clt_parameters.imp.pmtch_use_affine;
 		double      max_std =             clt_parameters.imp.pmtch_max_std;    // 1.5;      // maximal standard deviation to limit center area  
 		double      min_std_rad =         clt_parameters.imp.pmtch_min_std_rad;// 2.0;  // minimal radius of the central area (if less - fail)
 		double      rad_fraction =        clt_parameters.imp.pmtch_cent_rad; // center circle radius fraction of 0.5* min(width, height) in tiles
@@ -1126,6 +1119,7 @@ public class ComboMatch {
 		gd.addNumericField("Spiral search debug level",spiral_debug, 0,3,"","Debug level during Spiral search.");
 		
 		gd.addMessage("Parameters, common to all matching, not only spiral");
+		gd.addCheckbox    ("Use scenes' affine",               pmtch_use_affine, "Use known scenes' affine matrices, false - start from scratch (unity) ones.");
 		gd.addNumericField("Central area standard deviation",  max_std,  3,7,"",	"Central area limit by the standard deviation.");
 		gd.addNumericField("Central area minimal radius",      min_std_rad,  3,7,"tile",	"Minimal radius of the central area after all LMA passes.");
 		gd.addNumericField("Central area radius as fraction",  rad_fraction,   3,7,"",	"Central area radius as fraction of half minimal WOI dimension.");
@@ -1154,6 +1148,7 @@ public class ComboMatch {
 		log_path =                 gd.getNextString();
 		spiral_debug =       (int) gd.getNextNumber();
 		search_step=               gd.getNextNumber();
+		pmtch_use_affine=          gd.getNextBoolean();
 		max_std=                   gd.getNextNumber();
 		min_std_rad=               gd.getNextNumber();
 		rad_fraction=              gd.getNextNumber();
@@ -1180,6 +1175,8 @@ public class ComboMatch {
 				clt_parameters,  // CLTParameters    clt_parameters,
 				frac_remove,     // double           frac_remove, //  =        0.25
 				metric_error,    // double           metric_error,
+				pmtch_use_affine, // boolean          pmtch_use_affine,
+
 				max_std,         // double           max_std,     // maximal standard deviation to limit center area  
 				min_std_rad,     // double           min_std_rad, // minimal radius of the central area (if less - fail)
 				rad_fraction,    // double           rad_fraction,
@@ -1198,6 +1195,7 @@ public class ComboMatch {
 				min_overlap,     // int              min_overlap, // 3000
 				ignore_rms,      // boolean          ignore_rms,
 				null,//double []        max_rms_iter, //  = {1.0, 0.6};//
+				1.0, // double           overlap,
 				spiral_debug);     // int              debugLevel){
 		if (log_append && (log_path != null)) { // assuming directory exists
 			StringBuffer sb = new StringBuffer();

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

@@ -828,7 +828,8 @@ public class OrthoMap implements Comparable <OrthoMap>, Serializable{
 	 * image ones. Both are in meters relative to the "vert_meters" point
 	 * @param affine [2][3] array were the last column is made of the offsets (in meters)
 	 * @return [2][3] array that converts from the image metric coordinates relative to the
-	 * "vert_meters" point to the rectified metric coordinates relative to the same point.  
+	 * "vert_meters" point to the rectified metric coordinates relative to the same point.
+	 * 
 	 */
 	public static double [][] invertAffine(double [][] affine){
 		Matrix A = new Matrix (

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

@@ -5,13 +5,16 @@ import java.io.ObjectInputStream;
 import java.io.ObjectOutputStream;
 import java.io.Serializable;
 
+import Jama.Matrix;
+
 public class PairwiseOrthoMatch implements Serializable {
 	private static final long serialVersionUID = 1L;
 	public double [][]           affine = new double[2][3];
 	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 transient int [] nxy = null; // not saved, just to communicate for logging
+	public transient double overlap = 0.0;
 	public PairwiseOrthoMatch() {
 		
 	}
@@ -19,11 +22,13 @@ public class PairwiseOrthoMatch implements Serializable {
 			double [][] affine,
 			double [][] jtj,
 			double rms,
-			int zoom_lev) {
+			int zoom_lev,
+			double overlap) {
 		this.affine =  affine;
 		this.jtj =     jtj;
 		this.zoom_lev= zoom_lev;
 		this.rms =     rms;
+		this.overlap = overlap;
 	}
 	public PairwiseOrthoMatch clone() {
 		double [][] affine = {this.affine[0].clone(),this.affine[1].clone()};
@@ -35,19 +40,63 @@ public class PairwiseOrthoMatch implements Serializable {
 				affine,
 				jtj,
 				this.rms,
-				this.zoom_lev);
+				this.zoom_lev,
+				this.overlap);
 		if (nxy != null) {
 			pom.nxy = nxy.clone();
 		}
 		return pom;
 	}
+	/**
+	 * If this match is calculated from scene0 to scene1, inverse one is from scene1 to scene0  
+	 * @param rd displacement from scene0 reference (vertical) point to scene1 reference point
+	 * (GNSS-derived), referenced as -V in the notes below. (it is calculated as 0 from 1)
+	 * @return inverted 3x2 matrix
+	 * 
+		Xs0, Xs1 - metric coordinates (from top-left) in each of the source images
+		V - image center offset of the second scene from the first in rectified space (GPS-derived)
+		S0, S1 - vertical point offset from the top-left in image coordinates
+		A0,A1 - affine transform matrices (sensor coordinates from rectified coordinates)
+		B0, B1 centers offsets in image coordinates
+		Xr - rectified coordinates. Xr1, Xr2 - just different points
+		
+		Xsi = Ai*(Xr - V) + Bi + Si
+		
+		1) A00 = E, A1, B1 - second image match when the first is unity:
+		
+		Xs0 = Xr1 + S0
+		Xs1 = A1 * (Xr1 - V) + B1 + S1;
+		
+		2) calculate inverted A1,B1 of the first scene relative to the second one. In this case Xs0
+		of the first scene should still correspond to Xs1
+		Xs0 = A0 * (Xr2 + V) + B0 + S0;
+		Xs1 = Xr2 + S1
+		
+		----
+		Xs1 =  A1 * (Xs0 - S0 - V) + B1 + S1;
+		Xs1 - B1 - S1 =  A1 * (Xs0 - S0 - V);
+		A1inv*(Xs1 - B1 - S1) =  Xs0 - S0 - V;
+		Xs0 = A1inv*(Xs1 - B1 - S1) + S0 + V
+		---
+		Xr2 = Xs1 - S1
+		Xs0 = A0 * (Xs1 - S1 + V) + B0 + S0;
+		
+		A1inv*(Xs1 - B1 - S1) + S0 + V = A0 * (Xs1 - S1 + V) + B0 + S0
+		=> A0=A1inv
+		A1inv*(Xs1 - B1 - S1) + S0 + V = A1inv * (Xs1 - S1 + V) + B0 + S0
+		-A1inv*( B1)  + V = A1inv * ( V) + B0
+		A1inv * ( V) + B0 + A1inv*(B1) - V = 0
+		B0 = -((A1inv-E)*V + A1inv*B1)
+	 */
+	
 	public PairwiseOrthoMatch getInverse(double [] rd) {
 		double [][] affine = OrthoMap.invertAffine(getAffine());
 		PairwiseOrthoMatch inverted_match = new PairwiseOrthoMatch(
 				affine, // double [][] affine,
-				null, // double [][] jtj,
+				jtj, // double [][] jtj,
 				rms, // double rms,
-				zoom_lev); // int zoom_lev)
+				zoom_lev, // int zoom_lev)
+				overlap); // 
 		double [] corr = {
 				rd[0] * (affine[0][0]-1.0)+ rd[1]*affine[0][1],
 				rd[0] * affine[1][0]+       rd[1]*(affine[1][1]-1.0)};
@@ -55,14 +104,73 @@ public class PairwiseOrthoMatch implements Serializable {
 		affine[1][2] += corr[1];
 		return inverted_match;	
 	}
-
-	public PairwiseOrthoMatch getInverse() {
-		PairwiseOrthoMatch inverted_match = new PairwiseOrthoMatch(
-				OrthoMap.invertAffine(getAffine()), // double [][] affine,
-				null, // double [][] jtj,
-				rms, // double rms,
-				zoom_lev); // int zoom_lev)
-		return inverted_match;
+	
+	/**
+	 * Create differential PairwiseOrthoMatch instance from 2 affines of the scenes and an offset vector
+	 * @param affine0
+	 * @param affine1
+	 * @param rd
+	 * Continued from getInverse notes:
+	 * =========================== V = -rd
+		Xsi = Ai*(Xr - V) + Bi + Si
+		
+		Xs0 = A0*(Xr - 0) + B0 + S0
+		Xs1 = A1*(Xr - V) + B1 + S1
+		
+		--- difference Ad, Bd 
+		Xs0 = Xr1 + S0
+		Xs1 = Ad * (Xr1 - V) + Bd + S1;
+		
+		
+		Xs0 = A0*(Xr - 0) + B0 + S0
+		Xs0 - (B0 + S0) = A0*Xr
+		Xr = A0inv*(Xs0 - (B0 + S0))
+		Xs1 = A1*((A0inv*(Xs0 - (B0 + S0))) - V) + B1 + S1
+		
+		Xr1 = Xs0-S0
+		
+		Xs1 = Ad * ((Xs0-S0) - V) + Bd + S1;
+		
+		A1*((A0inv*(Xs0 - (B0 + S0))) - V) + B1 + S1 = Ad * ((Xs0-S0) - V) + Bd + S1
+		
+		A1*((A0inv*(Xs0 - (B0 + S0))) - V) + B1 = Ad * ((Xs0-S0) - V) + Bd
+		
+		A1*((A0inv*(Xs0 - (B0 + S0))) - V) + B1 = Ad * ((Xs0-S0) - V) + Bd
+		A1 * A0inv * Xs0 - A1 * A0inv * (B0 + S0) - A1*V + B1 - Ad *Xs0 + Ad * S0 +Ad*V -Bd
+		=> Ad = A1 * A0inv:
+		 - A1 * A0inv * (B0 + S0) - A1*V + B1  + Ad * S0 +Ad*V -Bd =
+		 - Ad * (B0 + S0) - A1*V + B1  + Ad * S0 + Ad*V - Bd =
+		 Ad * (S0 + V - B0 - S0) - A1*V + B1- Bd = 
+		  Ad * (V - B0) - A1 * V + B1 - Bd
+		Bd =  Ad * (V - B0) - A1 * V + B1
+		
+		--- check if A0 = E, B0 = 0
+		Ad = A1inv
+		Bd = A1 * V - A1 *V + B1 = B1
+		--- check if A1 = E, B1 = 0, V1= -V
+		Bd =  A0inv * (V1 - B0) - V1 =
+		(A0inv-E)*V1 -A0inv*B0 = 
+		-(A0inv-E)*V -A0inv*B0 = 
+	 */
+	public PairwiseOrthoMatch (
+			double [][] affine0,
+			double [][] affine1,
+			double [] rd) {
+		Matrix A0 = new Matrix (
+				new double [][] {{affine0[0][0],affine0[0][1]},{affine0[1][0],affine0[1][1]}});
+		Matrix B0 = new Matrix(new double [][] {{affine0[0][2]},{affine0[1][2]}});
+		Matrix A1 = new Matrix (
+				new double [][] {{affine1[0][0],affine1[0][1]},{affine1[1][0],affine1[1][1]}});
+		Matrix B1 = new Matrix(new double [][] {{affine1[0][2]},{affine1[1][2]}});
+		Matrix V =  new Matrix(new double [][] {{-rd[0]},{-rd[1]}});
+		Matrix A =  A1.times(A0.inverse());
+		Matrix B = A.times(V.minus(B0)).minus(A1.times(V)).plus(B1);
+		affine = new double[][] {
+				{A.get(0,0),A.get(0,1), B.get(0,0)},
+				{A.get(1,0),A.get(1,1), B.get(1,0)}};
+//				jtj = null;
+//				rms = Double.NaN; // double rms,
+//				zoom_lev = 0; // int zoom_lev)
 	}
 	
 	public double [][] getAffine(){
@@ -76,6 +184,7 @@ public class PairwiseOrthoMatch implements Serializable {
 				oos.writeObject(jtj[i][j]);
 			}
 		}
+		oos.writeObject(overlap);
 	}
 	private void readObject(ObjectInputStream ois) throws ClassNotFoundException, IOException {
 		ois.defaultReadObject();
@@ -88,6 +197,7 @@ public class PairwiseOrthoMatch implements Serializable {
 				}
 			}
 		}
+		overlap = (Double) ois.readObject();
 	}
 	//private void readObjectNoData() throws ObjectStreamException; // used to modify default values
 }

src/main/java/com/elphel/imagej/tileprocessor/IntersceneMatchParameters.java

@@ -108,6 +108,7 @@ public class IntersceneMatchParameters {
     public double   rln_neib_rstr =     0.35;    // minimal neighbors phase correlation maximums relative to max str
 	
     // Pairwise scenes matching
+    public boolean  pmtch_use_affine =  false;   // when matching pairs, start with known scene affine matrices, false - use unity
     public double   pmtch_frac_remove = 0.1;
 	public double   pmtch_metric_err =  0.05; // 0.02;//  2 cm
 	public double   pmtch_max_std =     1.5;      // maximal standard deviation to limit center area  
@@ -737,6 +738,7 @@ public class IntersceneMatchParameters {
 				"Minimal neighbors phase correlation maximums relative to maximal strength.");
 
 		gd.addMessage  ("Pairwise scenes matching");
+		gd.addCheckbox    ("Use scenes' affine",               this.pmtch_use_affine, "Use known scenes' affine matrices, false - start from scratch (unity) ones.");
 		gd.addNumericField("Remove fraction of worst matches", this.pmtch_frac_remove,  3,7,"",	"When fitting scenes remove this fraction of worst match tiles.");
 		gd.addNumericField("Maximal metric error",             this.pmtch_metric_err,  3,7,"m",	"Maximal tolerable fitting error caused by elevation variations.");
 		gd.addNumericField("Central area standard deviation",  this.pmtch_max_std,  3,7,"",	"Central area limit by the standard deviation.");
@@ -1622,6 +1624,7 @@ public class IntersceneMatchParameters {
 		this.rln_sngl_rstr =                  gd.getNextNumber();
 		this.rln_neib_rstr =                  gd.getNextNumber();
 		
+		this.pmtch_use_affine=                gd.getNextBoolean();
 		this.pmtch_frac_remove  =             gd.getNextNumber();
 		this.pmtch_metric_err   =             gd.getNextNumber();
 		this.pmtch_max_std      =             gd.getNextNumber();
@@ -2120,6 +2123,7 @@ public class IntersceneMatchParameters {
 		properties.setProperty(prefix+"rln_sngl_rstr",                 this.rln_sngl_rstr +"");                // double
 		properties.setProperty(prefix+"rln_neib_rstr",                 this.rln_neib_rstr +"");                // double
 
+		properties.setProperty(prefix+"pmtch_use_affine",              this.pmtch_use_affine +"");             // boolean
 		properties.setProperty(prefix+"pmtch_frac_remove",             this.pmtch_frac_remove +"");            // double
 		properties.setProperty(prefix+"pmtch_metric_err",              this.pmtch_metric_err +"");             // double
 		properties.setProperty(prefix+"pmtch_max_std",                 this.pmtch_max_std +"");                // double
@@ -2582,6 +2586,7 @@ public class IntersceneMatchParameters {
 		if (properties.getProperty(prefix+"rln_sngl_rstr")!=null)      this.rln_sngl_rstr=Double.parseDouble(properties.getProperty(prefix+"rln_sngl_rstr"));
 		if (properties.getProperty(prefix+"rln_neib_rstr")!=null)      this.rln_neib_rstr=Double.parseDouble(properties.getProperty(prefix+"rln_neib_rstr"));
 
+		if (properties.getProperty(prefix+"pmtch_use_affine")!=null)   this.pmtch_use_affine= Boolean.parseBoolean(properties.getProperty(prefix+"pmtch_use_affine"));		
    		if (properties.getProperty(prefix+"pmtch_frac_remove")!=null)  this.pmtch_frac_remove=Double.parseDouble(properties.getProperty(prefix+  "pmtch_frac_remove"));
 		if (properties.getProperty(prefix+"pmtch_metric_err")!=null)   this.pmtch_metric_err= Double.parseDouble(properties.getProperty(prefix+  "pmtch_metric_err"));
 		if (properties.getProperty(prefix+"pmtch_max_std")!=null)      this.pmtch_max_std=    Double.parseDouble(properties.getProperty(prefix+  "pmtch_max_std"));
@@ -3071,6 +3076,7 @@ public class IntersceneMatchParameters {
 		imp.rln_sngl_rstr                 = this.rln_sngl_rstr;
 		imp.rln_neib_rstr	              = this.rln_neib_rstr;
 
+		imp.pmtch_use_affine =              this.pmtch_use_affine;
         imp.pmtch_frac_remove             = this.pmtch_frac_remove;
 		imp.pmtch_metric_err              = this.pmtch_metric_err;
 		imp.pmtch_max_std                 = this.pmtch_max_std;