More on quaternions - affines

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

@@ -177,7 +177,7 @@ public class ComboMatch {
 		
 		
 		boolean use_saved_collection = true; // false;
-		boolean save_collection =     true;
+		boolean save_collection =     false; // true;
 		boolean process_correlation = true; // use false to save new version of data
 		int     num_tries_fit =       10;
 		boolean update_match =        true; // use false to save new version of data
@@ -210,7 +210,7 @@ public class ComboMatch {
 			pattern_match = false;
 		}
 		
-		GenericJTabbedDialog gd = new GenericJTabbedDialog("Set image pair",1200,800);
+		GenericJTabbedDialog gd = new GenericJTabbedDialog("Set image pair",1200,900);
 //		gd.addChoice      ("Files list/data path (w/o extension):", FILES_LISTS_PATHS, FILES_LISTS_PATHS[default_list_choice]);
 		gd.addChoice      ("Files list/data path (w/o extension):", FILES_LISTS_PATHS, omtch_img_set);
 		gd.addCheckbox    ("Use saved maps collection ", use_saved_collection, " (if available). If false - use files list.");
@@ -433,7 +433,7 @@ public class ComboMatch {
 					};
 			
 			
-			GenericJTabbedDialog gdo = new GenericJTabbedDialog("Set image pair",1200,900);
+			GenericJTabbedDialog gdo = new GenericJTabbedDialog("Set image pair",1200,1000);
 			gdo.addNumericField("Correlation size",          corr_size,       0,4,"",
 					"Correlation size, power of 2");
 			gdo.addNumericField("Extracted size",            extr_size,  0,4,"",
@@ -858,6 +858,12 @@ public class ComboMatch {
 		}
         if (process_correlation || render_match || pattern_match || test_multi_lma) {
 //        	int [] gpu_pair;
+			boolean removeTilt =   false;
+			boolean removeRot =    false;
+			boolean removeScale =  false;
+			boolean removeOffset = false;
+			boolean max_is_scale = false;
+
         	if (gpu_spair == null) {
         		ArrayList<Point> pairs_list = new ArrayList<Point>(); 
         		for (OrthoMap map : maps_collection.ortho_maps) {
@@ -905,7 +911,7 @@ public class ComboMatch {
         		boolean flt_show_alt =       clt_parameters.imp.flt_show_alt;      //  true;
         		boolean flt_update_config = false;
         		String  flt_extra_line =    "--- select a single image ---";
-				GenericJTabbedDialog gdf = new GenericJTabbedDialog("Select pairs filter/display",800,500);
+				GenericJTabbedDialog gdf = new GenericJTabbedDialog("Select pairs filter/display",800,600);
 				gdf.addCheckbox    ("Filter pairs",                flt_list, "Filter available pairs.");
 				gdf.addCheckbox    ("Keep undefined pairs only",   flt_undef_only, "Keep only undefined pairs (affines== null).");
 				gdf.addNumericField("Minimal scene overlap (0..1)",flt_min_overlap,  3,7,"", "Minimal overlap of the scenes to keep (0-no overlap, 1.0 - smaller scene is inside the parger one.");
@@ -999,15 +1005,27 @@ public class ComboMatch {
 						false,             // boolean          use_tab,
 						flt_extra_line);   // String           extra_line)
 				
-				GenericJTabbedDialog gdc = new GenericJTabbedDialog("Select image pair",1200,100);
+				GenericJTabbedDialog gdc = new GenericJTabbedDialog("Select image pair",1200,300);
 				int num_choice_lines = 50;
 				gdc.addChoice("Image pair:",
 						choices_all,
 						choices_all[choices_all.length - 1], // -1],
 						"Select processed image pair or request a single image selection", num_choice_lines);
+				gdc.addCheckbox    ("Remove tilt",             removeTilt,   "Remove tilts from the pairwise affine transform.");
+				gdc.addCheckbox    ("Remove Rotation",         removeRot,    "Remove rotation from the pairwise affine transform.");
+				gdc.addCheckbox    ("Remove Scale",            removeScale,  "Remove scale from the pairwise affine transform.");
+				gdc.addCheckbox    ("Remove Translation",      removeOffset, "Remove translation from the pairwise affine transform.");
+				gdc.addCheckbox    ("Max eigenvalue is scale", max_is_scale, "Maximal eigenvalue is scale (false - minimal is).");
+				//        			boolean max_is_scale = false;
+
 				gdc.showDialog();
 				if (gdc.wasCanceled()) return false;
 				int pair=      gdc.getNextChoiceIndex();
+				removeTilt =   gdc.getNextBoolean();
+				removeRot =    gdc.getNextBoolean();
+				removeScale =  gdc.getNextBoolean();
+				removeOffset = gdc.getNextBoolean();
+				max_is_scale = gdc.getNextBoolean();
 				if (pair >= (choices_all.length -1)) {
 					int default_choice = 0;
 					int num_scene_lines = 50;
@@ -1118,7 +1136,36 @@ public class ComboMatch {
         			
         			
         		} else {
-        			affine1 = pairwiseOrthoMatch.getAffine();
+        			affine1 = pairwiseOrthoMatch.getAffine().clone();
+        			if (debugLevel > -4) {
+            			System.out.println("removeTilt="+  removeTilt);
+            			System.out.println("removeRot="+   removeRot);
+            			System.out.println("removeScale="+ removeScale);
+            			System.out.println("removeOffset="+removeOffset);
+            			System.out.println("max_is_scale="+max_is_scale);
+            			if (removeTilt || removeRot || removeScale || removeOffset) {
+                			SingularValueDecomposition svd_affine_pair =  SingularValueDecomposition.singularValueDecomposeScaleTiltGamma(affine1, true); // y_down_ccw); // boolean y_down_ccw)
+                			System.out.println("Full (original) affine transform:");
+            	        	System.out.println(QuatUtils.affinesToString(affine1, "full_affine"));
+            	        	System.out.println("full_affine_svd" + svd_affine_pair.toString(true)); // use_degrees));
+            			}
+        			}
+        			affine1 = SingularValueDecomposition.removeTiltRotScale(
+        					affine1,       // double [][] A,
+        					removeTilt,    // boolean removeTilt,
+        					removeRot,     // boolean removeRot,
+        					removeScale,   // boolean removeScale,
+        					removeOffset,  // boolean removeOffset)
+        					max_is_scale); //boolean max_is_scale)
+        			if (debugLevel > -4) {
+            			if (removeTilt || removeRot || removeScale || removeOffset) {
+                			SingularValueDecomposition svd_affine_pair =  SingularValueDecomposition.singularValueDecomposeScaleTiltGamma(affine1, true); // y_down_ccw); // boolean y_down_ccw)
+                			System.out.println("Final differential affine transform:");
+            	        	System.out.println(QuatUtils.affinesToString(affine1, "affine1    "));
+            	        	System.out.println("affine1_svd    " + svd_affine_pair.toString(true)); // use_degrees));
+            			}
+        			}
+        			
         			double [][][] affines = {affine0,affine1};
         			int [] zooms = {initial_zoom, min_zoom_lev, 1000,1000}; // make automatic
 //        			double scale = 2.0; // scale vectors when warping;
@@ -1239,6 +1286,14 @@ public class ComboMatch {
         				}
         				if (render_match) {
         					String title=String.format("multi_%03d-%03d_%s-%s_zoom%d_%d",gpu_pair[0],gpu_pair[1],gpu_spair[0],gpu_spair[1],min_zoom_lev,zoom_lev);
+        					if (removeTilt || removeRot || removeScale || removeOffset) {
+        						title += "_no-";
+        						if (removeTilt)   title += "T";
+        						if (removeRot)    title += "R";
+        						if (removeScale)  title += "S";
+        						if (removeOffset) title += "O";
+        						if (max_is_scale) title += "X"; 
+        					}
         					// Avoid renderMulti() - it duplicates code renderMultiDouble()
         					int eq_mode = 2; // calculate
         					ImagePlus imp_img_pair = 	maps_collection.renderMulti (
@@ -2077,7 +2132,7 @@ adjusted affines[1] for a pair: 1694564291_293695/1694564778_589341
 					{2139.833, 1474.5},
 					{ 461.5,   1493.5}}};
 		boolean set_coords = true; // false;
-		GenericJTabbedDialog gd = new GenericJTabbedDialog("Set image pair",1090,900);
+		GenericJTabbedDialog gd = new GenericJTabbedDialog("Select processing moder",1090,900);
 		gd.addStringField ("First image path",  image_paths[0], 120, "First image full path");		
 		gd.addStringField ("Second image path", image_paths[1], 120, "Second image full path");		
 		gd.addCheckbox    ("Set image markers to following coordinates", set_coords, "Delete markers and set them according to data.");

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

@@ -3732,7 +3732,21 @@ public class OrthoMap implements Comparable <OrthoMap>, Serializable{
 			final double []  weight,
 			final int        width,
 			final double []  xy0) {
-		int debug = -1;
+		return getPlane(
+				   data, // final double []   data,
+				mask,  // final boolean [] mask,
+				weight, //final double []  weight,
+				width, // final int        width,
+				xy0, // final double []  xy0,
+				null); // String dbg_title)
+	}
+	public static double [] getPlane(
+			final double []   data,
+			final boolean [] mask,
+			final double []  weight,
+			final int        width,
+			final double []  xy0,
+			String dbg_title) {
 		final double [][][] mdatai = new double [data.length][][];
 		AtomicInteger anum_good = new AtomicInteger(0);
 		final Thread[] threads = ImageDtt.newThreadArray();
@@ -3770,7 +3784,7 @@ public class OrthoMap implements Comparable <OrthoMap>, Serializable{
 				null,       // damping,    // double [] damping, null OK
 				-1);        // debug level
 		if (approx2d != null) {
-			if (debug > 0) {
+			if (dbg_title != null) {
 				double [] plane=new double[data.length];
 				double [] diff=new double[data.length];
 				double [] dweight = (weight != null) ? weight: new double[data.length];
@@ -3794,7 +3808,7 @@ public class OrthoMap implements Comparable <OrthoMap>, Serializable{
 						width,
 						data.length/width,
 						true,
-						"plane_approximation",
+						dbg_title,
 						new String[] {"data","approx","diff","masked","weight","mask"});
 			}
 			

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

@@ -5204,7 +5204,7 @@ public class OrthoMapsCollection  implements Serializable{
 		
 		boolean flt_update_config =  false;
 		
-		GenericJTabbedDialog gdf = new GenericJTabbedDialog("Select pairs filter/display",800,1100);
+		GenericJTabbedDialog gdf = new GenericJTabbedDialog("Select pairs filter/display",800,1200);
 		gdf.addMessage("Filter pairs parameters");
 		gdf.addNumericField("Minimal scene overlap (0..1)",flt_min_overlap,  3,7,"", "Minimal overlap of the scenes to keep (0-no overlap, 1.0 - smaller scene is inside the parger one.");
 		gdf.addNumericField("Maximal scene overlap (0..1)",flt_max_overlap,  3,7,"", "Maximal overlap of the scenes to keep (0-no overlap, 1.0 - smaller scene is inside the parger one.");
@@ -6342,7 +6342,7 @@ public class OrthoMapsCollection  implements Serializable{
 		boolean flt_update_config =  false;
 		boolean select_pairs =       false;
 		
-		GenericJTabbedDialog gdf = new GenericJTabbedDialog("Select pairs filter/display",800,1100);
+		GenericJTabbedDialog gdf = new GenericJTabbedDialog("Select pairs filter/display",800,1200);
 		gdf.addMessage("Filter pairs parameters");
 		gdf.addNumericField("Minimal scene overlap (0..1)",flt_min_overlap,  3,7,"", "Minimal overlap of the scenes to keep (0-no overlap, 1.0 - smaller scene is inside the parger one.");
 		gdf.addNumericField("Maximal scene overlap (0..1)",flt_max_overlap,  3,7,"", "Maximal overlap of the scenes to keep (0-no overlap, 1.0 - smaller scene is inside the parger one.");
@@ -6581,7 +6581,7 @@ public class OrthoMapsCollection  implements Serializable{
 		boolean flt_update_config =  false;
 		
 		
-		GenericJTabbedDialog gdf = new GenericJTabbedDialog("Select pairs filter/display",800,1100);
+		GenericJTabbedDialog gdf = new GenericJTabbedDialog("Select pairs filter/display",800,1200);
 		gdf.addMessage("Filter pairs parameters");
 		gdf.addNumericField("Minimal scene overlap (0..1)",flt_min_overlap,  3,7,"", "Minimal overlap of the scenes to keep (0-no overlap, 1.0 - smaller scene is inside the parger one.");
 		gdf.addNumericField("Maximal scene overlap (0..1)",flt_max_overlap,  3,7,"", "Maximal overlap of the scenes to keep (0-no overlap, 1.0 - smaller scene is inside the parger one.");

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

@@ -45,7 +45,8 @@ public final class QuatUtils { // static class
 		String fmt_deg = "[%12.9f, %12.9f,%12.9f, %12.9f], tilt=%12.7f\u00B0, dir=%12.7f\u00B0, scale=%12.10f";
 		double s = degrees ? (180/Math.PI):1;
 		String fmt=degrees ? fmt_deg : fmt_rad;
-		return String.format(fmt, quat[0], quat[1],quat[2],quat[3], s*2*Math.acos(quat[0]), s*Math.atan2(quat[2],quat[1]), scale);
+		String rslt=String.format(fmt, quat[0], quat[1],quat[2],quat[3], s*2*Math.acos(quat[0]), s*Math.atan2(quat[2],quat[1]), scale);
+		return rslt;
 	}
 	
 	/**
@@ -300,7 +301,10 @@ public final class QuatUtils { // static class
 		if (t == 0) {
 			return UNIT_QUAT;
 		}
-		double [] axis = {ty/t,-tx/t,0}; // pi/2 CW
+		tx/=t;
+		ty/=t;
+		
+		double [] axis = {ty,-tx,0}; // pi/2 CW
 		
 		/* sin (A/2) = +-sqrt((1 - cos(A))/2)
 		 * cos (A/2) = +-sqrt((1 + cos(A))/2)
@@ -313,6 +317,64 @@ public final class QuatUtils { // static class
 		return quat;
 	}
 	
+	public static String tiltToString (
+			double [] txy,
+			boolean	y_down_ccw,
+			boolean degrees) {
+		String fmt_rad = " tilt= %10.7f, dir=%10.7f";
+		String fmt_deg = " tilt= %10.5f\u00B0, dir=%10.5f\u00B0";
+		double k = degrees ? (180/Math.PI):1;
+		String fmt=degrees ? fmt_deg : fmt_rad;
+
+		String s = String.format(" = [%12.9f, %12.9f, %15.9f]", txy[0], txy[1], txy[2]);
+		double tx = txy[0], ty = y_down_ccw?-txy[1]:txy[1]; // invert y
+		double t2 = tx*tx + ty*ty;
+		double t = Math.sqrt(t2);
+		if (t == 0) {
+			return s;
+		}
+		tx/=t;
+		ty/=t;
+		
+		double [] axis = {ty,-tx,0}; // pi/2 CW
+		double dir = Math.atan2(axis[1],axis[0]);
+		double tilt = Math.atan(t);
+		
+		s+= String.format(fmt, k*tilt, k*dir);
+		return s;
+	}
+	
+	
+	public static double [] sceneRelLocalGround(
+			double []   txy,
+			double [][] affine, // used for rot and scale Now can be null
+			boolean	    y_down_ccw) { // boolean y_down_ccw) for tilts and affine
+		double [] qtilts = tiltToQuaternion(txy,y_down_ccw);
+		boolean inv_aff = true; // invert affine
+		double scale = 1.0;
+		double rot = 0;
+		if (affine != null) {
+			if (inv_aff) {
+				affine = matInverse2x2(affine);
+			}
+			SingularValueDecomposition svd= SingularValueDecomposition.singularValueDecomposeScaleTiltBeta(
+					affine,
+					y_down_ccw);
+			scale = svd.getMaxScale();
+			rot = svd.getRotAngle(); // TODO: check sign !
+		}		
+		
+		// Now beta,rot correspond to Y - up
+		double chalfrot = Math.cos(rot/2); 
+		double shalfrot = Math.sin(rot/2); 
+		double [] qrot = {chalfrot,0,0,shalfrot}; // rotation around vertical axis
+		double [] qinv_tilts = invert(qtilts); //scene to ground
+		double [] qrt=  multiply(qrot,qinv_tilts);
+//		double [] quat= scale(multiply(qrot,qinv_tilts),scale);
+		double [] quat= scale(multiply(qinv_tilts,qrot),scale);
+		return    quat;
+	}	
+	
 	/**
 	 * Remove rotation around Z-axis, keep only tilt around axis in XY plane
 	 * @param quat_in input rotation that may include rotation around Z (vertical) axis
@@ -563,7 +625,7 @@ q11 - q22 = c3
 	 */
 	public static double [][] affineToQuatScaled (
 			double [][] affine,
-			boolean     stretch,
+			boolean     stretch,// false
 			boolean     y_down_ccw) {
 		affine = matInverse2x2(affine);
 		

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

@@ -29,10 +29,16 @@ public class SingularValueDecomposition {
 		
 	}
 	// A = getR1() * getW() getR2() = getR1() * getW() * getIR1() * getRot()
+	public static double [][] getUnity(){
+		return new double [][] {{1,0},{0,1}};
+	}
 	public double [][] getR1(){	return rotMatrix(beta);}
 	public double [][] getR2(){	return rotMatrix(gamma);}
-	public double [][] getIR1(){return rotMatrix(-beta);}
+	public double [][] getIR1(){return rotMatrix(-beta);} // never used
 	public double [][] getRot(){return rotMatrix(rot);}
+	public double [][] getW(){
+		return new double [][] {{w1,0},{0,w1}};
+	}
 	public double getTiltAngle() {
 		if (ratio <= 1.0) return Math.acos(ratio);
 		else return Math.acos(1/ratio);
@@ -46,19 +52,16 @@ public class SingularValueDecomposition {
 	public String toString(boolean degrees) {
 //System.out.println("svd_affine_pair=    ["+svd_affine_pair.scale+  ","+svd_affine_pair.getTiltAngle()+ ","+svd_affine_pair.gamma+ ","+svd_affine_pair.rot+
 // "] tilt="+(svd_affine_pair.getTiltAngle()*180/Math.PI)+ "\u00B0, dir="+(svd_affine_pair.gamma*180/Math.PI)+"\u00B0");
-		String fmt_rad = "scale=%10.8f,tilt= %10.7f, gamma=%10.7f, beta=%10.7f, rot=%10.7f, w1=%10.8f, w2=%10.8f, ratio=%10.8f";
-		String fmt_deg = "scale=%10.8f,tilt= %10.5f\u00B0, gamma=%10.5f\u00B0, beta=%10.5f\u00B0, rot=%10.5f\u00B0, w1=%10.8f, w2=%10.8f, ratio=%10.8f";
+		String fmt_rad = " scale=%10.8f,tilt= %10.7f, gamma=%10.7f, beta=%10.7f, rot=%10.7f, w1=%10.8f, w2=%10.8f, ratio=%10.8f";
+		String fmt_deg = " scale=%10.8f,tilt= %10.5f\u00B0, gamma=%10.5f\u00B0, beta=%10.5f\u00B0, rot=%10.5f\u00B0, w1=%10.8f, w2=%10.8f, ratio=%10.8f";
 		double s = degrees ? (180/Math.PI):1;
 		String fmt=degrees ? fmt_deg : fmt_rad;
 		return String.format(fmt, scale, s*getTiltAngle(), s*gamma, s*beta, s*rot, w1, w2, ratio);
 	}
 	
-	public double [][] getW(){
-		return new double [][] {{w1,0},{0,w1}};
-	}
 	
 	/**
-	 * Decomposing linear transform into rotation
+	 * Decomposing linear transform into rotations and scales. OK to use full 2x3 affine
 	 * R1={{cos(beta),sin(beta)},{-sin(beta), cos(beta)}} ,
 	 * non-uniform scale transform W={{w1,0}{0,w2}}, and rotation
 	 * R2={{cos(gamma),sin(gamma)},{-sin(gamma), cos(gamma)}}
@@ -91,6 +94,64 @@ public class SingularValueDecomposition {
 		return svd;
 	}
 
+	
+	public static double [][] removeTiltRotScale(
+			double [][] A,
+			boolean removeTilt,
+			boolean removeRot,
+			boolean removeScale,
+			boolean removeOffset,
+			boolean max_is_scale) {
+		SingularValueDecomposition svd = singularValueDecompose(A);
+		double [][] AR = removeTiltRotScale(
+				svd, // SingularValueDecomposition svd,
+				removeTilt,   // boolean removeTilt,
+				removeRot,    // boolean removeRot,
+				removeScale, // boolean removeScale)
+				max_is_scale); // boolean max_is_scale)
+		if (A[0].length < 3) {
+			return AR;
+		}
+		A=A.clone();
+		if (removeOffset) {
+			for (int i = 0; i < 2; i++) {
+				A[i][2] = 0;	
+			}
+		}
+		for (int i = 0; i < 2; i++) {
+			for (int j = 0; j < 2; j++) {
+				A[i][j] = AR[i][j];	
+			}
+		}
+		return A;
+	}
+	
+	public static double [][] removeTiltRotScale(
+			SingularValueDecomposition svd,
+			boolean removeTilt,
+			boolean removeRot,
+			boolean removeScale,
+			boolean max_is_scale) { 
+		double [][] R1 =  svd.getR1();
+		double [][] IR1 = svd.getIR1();
+//		double [][] W =   svd.getW();
+		double scale = max_is_scale?Math.max(svd.w1,svd.w2):Math.min(svd.w1,svd.w2);
+		double w1 =svd.w1, w2 = svd.w2;
+		if (removeScale) {
+			w1 /= scale;
+			w2 /= scale;
+			scale = 1.0;
+		}
+		if (removeTilt) {
+			w1 = scale;
+			w2 = scale;
+		}
+		double [][]W = {{w1,0},{0, w2}};
+		double [][] R =   removeRot?   getUnity() : svd.getRot();
+		double [][] A = QuatUtils.matMult2x2(QuatUtils.matMult2x2(QuatUtils.matMult2x2(R1, W), IR1), R);
+		return A;
+	}	
+	
 	public static double [] getMinMaxEigenValues(double [][] A) {
 		double a00=A[0][0],a01=A[0][1],a10=A[1][0],a11=A[1][1];
 		double b00=(a00+a11)/2; // , b11 = b00;