working on expansion + bug fixes

src/main/java/AlignmentCorrection.java

-import java.util.ArrayList;
-
-import Jama.Matrix;
-import ij.IJ;
-import ij.ImagePlus;
-import ij.ImageStack;
-import ij.WindowManager;
-
 /**
  **
  ** AlignmentCorrection - try to apply minor adjustments to the misaligned camera
@@ -29,6 +21,13 @@ import ij.WindowManager;
  ** -----------------------------------------------------------------------------**
  **
  */
+import java.util.ArrayList;
+
+import Jama.Matrix;
+import ij.IJ;
+import ij.ImagePlus;
+import ij.ImageStack;
+import ij.WindowManager;
 
 public class AlignmentCorrection {
 	static int NUM_SLICES =      10; // disp, strength, dx0, dy0, dx1, dy1, dx2, dy2, dx3, dy3)

src/main/java/EyesisCorrectionParameters.java

@@ -2070,6 +2070,9 @@ public class EyesisCorrectionParameters {
   		public double     ex_strength      = 0.18;  // minimal 4-corr strength to trust tile
   		public double     ex_nstrength     = 0.4;   // minimal 4-corr strength divided by channel diff for new (border) tiles
   		
+  		public boolean    ex_over_bgnd     = false; // Allow expansion over previously identified background (infinity)
+  		public double     ex_min_over      = 1.0;   // When expanding over background, disregard lower disparity 
+  		
   		public double     pt_super_trust   = 1.6;   // If strength exceeds ex_strength * super_trust, do not apply ex_nstrength and plate_ds
   		public boolean    pt_keep_raw_fg   = true;  // Do not replace raw tiles by the plates, if raw is closer (like poles)
   		public double     pt_scale_pre     = 1.5;   // Scale plates strength before comparing to raw strength
@@ -2417,7 +2420,7 @@ public class EyesisCorrectionParameters {
 		public boolean    tsLoopMulti          = true;   // Repeat multi-choice assignment while succeeding 
 		public boolean    tsReset              = false;  // Reset tiles to surfaces assignment
 		public boolean    tsShow               = false;  // Show results of tiles to surfaces assignment
-		public int        tsNumClust           = 50;     // Number of clusters to keep
+		public int        tsNumClust           = 500;     // Number of clusters to keep
 
 		public int        tsConsensMode        = 7;      // Which assignments to match +1 - combo, +2 grown single, +4 plane seeds 
 		public int        tsConsensAgree       = 1;      // Minimal number of assignments to agree
@@ -2624,6 +2627,9 @@ public class EyesisCorrectionParameters {
 			properties.setProperty(prefix+"ex_strength",      this.ex_strength +"");
 			properties.setProperty(prefix+"ex_nstrength",     this.ex_nstrength +"");
 			
+			properties.setProperty(prefix+"ex_over_bgnd",     this.ex_over_bgnd+"");
+			properties.setProperty(prefix+"ex_min_over",      this.ex_min_over +"");
+
 			properties.setProperty(prefix+"pt_super_trust",   this.pt_super_trust +"");
 			properties.setProperty(prefix+"pt_keep_raw_fg",   this.pt_keep_raw_fg+"");
 			properties.setProperty(prefix+"pt_scale_pre",     this.pt_scale_pre +"");
@@ -3141,6 +3147,9 @@ public class EyesisCorrectionParameters {
   			if (properties.getProperty(prefix+"ex_strength")!=null)       this.ex_strength=Double.parseDouble(properties.getProperty(prefix+"ex_strength"));
   			if (properties.getProperty(prefix+"ex_nstrength")!=null)      this.ex_nstrength=Double.parseDouble(properties.getProperty(prefix+"ex_nstrength"));
 
+  			if (properties.getProperty(prefix+"ex_over_bgnd")!=null)      this.ex_over_bgnd=Boolean.parseBoolean(properties.getProperty(prefix+"ex_over_bgnd"));
+  			if (properties.getProperty(prefix+"ex_min_over")!=null)       this.ex_min_over=Double.parseDouble(properties.getProperty(prefix+"ex_min_over"));
+
   			if (properties.getProperty(prefix+"pt_super_trust")!=null)    this.pt_super_trust=Double.parseDouble(properties.getProperty(prefix+"pt_super_trust"));
   			if (properties.getProperty(prefix+"pt_keep_raw_fg")!=null)    this.pt_keep_raw_fg=Boolean.parseBoolean(properties.getProperty(prefix+"pt_keep_raw_fg"));
   			if (properties.getProperty(prefix+"pt_scale_pre")!=null)      this.pt_scale_pre=Double.parseDouble(properties.getProperty(prefix+"pt_scale_pre"));
@@ -3680,6 +3689,9 @@ public class EyesisCorrectionParameters {
   			gd.addNumericField("Minimal 4-corr strength to trust tile",                                        this.ex_strength,  3);
   			gd.addNumericField("Minimal 4-corr strength divided by channel diff for new (border) tiles",       this.ex_nstrength,  3);
   			
+  			gd.addCheckbox    ("Allow expansion over previously identified background (infinity)",             this.ex_over_bgnd);
+  			gd.addNumericField("When expanding over background, disregard lower disparity ",                   this.ex_min_over,  3);
+
   			gd.addMessage     ("********* Plates filetering when building initial z-map *********");
   			gd.addNumericField("If strength exceeds ex_strength * super_trust, do not apply ex_nstrength and plate_ds", this.pt_super_trust,  3);
   			gd.addCheckbox    ("Do not replace raw tiles by the plates, if raw is closer (like poles)",        this.pt_keep_raw_fg);
@@ -3803,10 +3815,10 @@ public class EyesisCorrectionParameters {
   			gd.addNumericField("Consider merging initial planes if disparity difference below",                this.stSmallDiff,  6);
   			gd.addNumericField("Consider merging initial planes if jumps between ratio above",                 this.stHighMix,  6);
 
-  			gd.addNumericField("Outlayer tiles weaker than this may be replaced from neighbors",               this.outlayerStrength,  6);
-  			gd.addNumericField("Replace weak outlayer tiles that do not have neighbors within this disparity difference", this.outlayerDiff,  6);
-  			gd.addNumericField("Replace weak outlayer tiles that have higher disparity than weighted average", this.outlayerDiffPos,  6);
-  			gd.addNumericField("Replace weak outlayer tiles that have lower disparity than weighted average",  this.outlayerDiffNeg,  6);
+  			gd.addNumericField("Outlier tiles weaker than this may be replaced from neighbors",               this.outlayerStrength,  6);
+  			gd.addNumericField("Replace weak outlier tiles that do not have neighbors within this disparity difference", this.outlayerDiff,  6);
+  			gd.addNumericField("Replace weak outlier tiles that have higher disparity than weighted average", this.outlayerDiffPos,  6);
+  			gd.addNumericField("Replace weak outlier tiles that have lower disparity than weighted average",  this.outlayerDiffNeg,  6);
 
   			gd.addCheckbox    ("Combine with all previous after refine pass",                                  this.combine_refine);
   			gd.addNumericField("Disregard weaker tiles when combining scans",                                  this.combine_min_strength,  6);
@@ -4222,6 +4234,9 @@ public class EyesisCorrectionParameters {
   			this.ex_strength=           gd.getNextNumber();
   			this.ex_nstrength=          gd.getNextNumber();
 
+  			this.ex_over_bgnd=          gd.getNextBoolean();
+  			this.ex_min_over=           gd.getNextNumber();
+  			
   			this.pt_super_trust=        gd.getNextNumber();
   			this.pt_keep_raw_fg=        gd.getNextBoolean();
   			this.pt_scale_pre=          gd.getNextNumber();

src/main/java/LinkPlanes.java

@@ -3879,7 +3879,7 @@ public class LinkPlanes {
 					weak_pairs,       // final boolean [][]              weak_fg_pairs, 
 					1.0,              // double                    relax,
 					plWeakFgRelax,    // final double                    relax_weak_fg,
-					debugLevel,                  // final int        debugLevel)
+					1+ debugLevel,                  // final int        debugLevel)
 					dbg_tileX,
 					dbg_tileY);
 			// Calculate costs of merging planes of the same supertile and remove those that are exceed threshold

src/main/java/MeasuredLayers.java

@@ -337,7 +337,7 @@ public class MeasuredLayers {
 	 * @return [smplSide * smplSide]array of weights, 1.0 in the center
 	 */
 	
-	public double [] getSampleWindow(int smplSide, boolean all1){
+	static double [] getSampleWindow(int smplSide, boolean all1){
 		double [] weights = new double [smplSide * smplSide];
 		for (int sy = 0; sy < smplSide; sy++){
 			for (int sx = 0; sx < smplSide; sx++){
@@ -898,6 +898,9 @@ public class MeasuredLayers {
 			boolean    null_if_none,
 			int        debugLevel)
 	{
+		boolean use_new = true; // false;
+		
+		if (use_new) {
 		return getDisparityStrength (
 				num_layer,      // int num_layer,
 				stX,            // int stX,
@@ -914,6 +917,25 @@ public class MeasuredLayers {
 				0.2,            // double     max_rel_tilt, //  = 0.2; // (pix / disparity) per tile
 				null_if_none,   // boolean    null_if_none,
 				debugLevel);    // int        debugLevel)
+		} else {
+			return getDisparityStrength_old (
+					num_layer,      // int num_layer,
+					stX,            // int stX,
+					stY,            // int stY,
+					sel_in,         // boolean [] sel_in,
+//					null,           // double []  tiltXY, // null - free with limit on both absolute (2.0?) and relative (0.2) values 
+					strength_floor, // double     strength_floor,
+					strength_pow,   // double     strength_pow,
+					smplSide,       // int        smplSide, //        = 2;      // Sample size (side of a square)
+					smplNum,        // int        smplNum, //         = 3;      // Number after removing worst (should be >1)
+					smplRms,        // double     smplRms, //         = 0.1;    // Maximal RMS of the remaining tiles in a sample
+					smplWnd,        // boolean    smplWnd, //
+//					2.0,            // double     max_abs_tilt, //  = 2.0; // pix per tile
+//					0.2,            // double     max_rel_tilt, //  = 0.2; // (pix / disparity) per tile
+					null_if_none,   // boolean    null_if_none,
+					debugLevel);    // int        debugLevel)
+			
+		}
 	}
 	
 	public double[][] getDisparityStrength (
@@ -921,7 +943,9 @@ public class MeasuredLayers {
 			int stX,
 			int stY,
 			boolean [] sel_in,
-			double []  tiltXY, // null - free with limit on both absolute (2.0?) and relative (0.2) values 
+			double [][] atiltXY, // null - free with limit on both absolute (2.0?) and relative (0.2) values
+			                     // if it has just one element - apply to all tiles, otherwise it is supposed
+			                     // to be per-tile of a supertile array
 			double     strength_floor,
 			double     strength_pow,
 			int        smplSide, //        = 2;      // Sample size (side of a square)
@@ -1018,7 +1042,7 @@ public class MeasuredLayers {
 					if (num_in_sample >= smplNum){ // try, remove worst
 						sample_loop:
 						{
-						boolean en_tilt = (tiltXY == null);
+						boolean en_tilt = (atiltXY == null);
 						if (en_tilt) { // make sure there are enough samples and not all of them are on the same line
 							if (!notColinear(smpl_sel,smplSide)){
 								en_tilt = false;
@@ -1128,7 +1152,10 @@ public class MeasuredLayers {
 										}
 									}
 									if (iworst < 0){
-										System.out.println("**** this is a BUG in getDisparityStrength() can not find the worst sample ****");
+										if (debugLevel > 0) {
+											System.out.println("**** this may be BUG in getDisparityStrength() can not find the worst sample  - all tiles fit perfectly ****");
+										}
+										// this can happen if some samples are the same and all the pixels fit exactly - use all of them
 										break;
 									}
 									// remove worst sample
@@ -1158,6 +1185,15 @@ public class MeasuredLayers {
 							
 						} else { // fixed tilt
 							// tilt around center
+							double [] tiltXY= {0.0,0.0};
+							if (atiltXY != null){ // if null ( free but failed co-linear test - keep both tilts == 0.0
+								if (atiltXY.length == 1){
+									tiltXY = atiltXY[0]; // common for all tiles (such as constant disparity)
+								} else if (atiltXY[indx] != null){
+									tiltXY = atiltXY[indx];
+								}
+							}
+							
 							if ((tiltXY != null) && (tiltXY[0] != 0.0) && (tiltXY[1] != 0.0)){
 								for (int sy = 0; sy < smplSide; sy++){
 									for (int sx = 0; sx < smplSide; sx++){

src/main/java/SuperTiles.java

@@ -1436,13 +1436,17 @@ public class SuperTiles{
 		// TODO: Remove when promoting PlaneData
 		final TilePlanes tpl = new TilePlanes(tileSize,superTileSize, geometryCorrection);
 		final double [][][][] plane_disp_strength = new double [nStiles][][][];
+		// DEBUG feature:
+		final double [][] zero_tilts = {{0.0,0.0}}; // set to null for float
 		for (int ithread = 0; ithread < threads.length; ithread++) {
 			threads[ithread] = new Thread() {
 				public void run() {
+					double [][][] plane_tilts = null; // used only forworld_plane_norm != null
+					// to get tile disparities needed to calculate tilts
 					for (int nsTile = ai.getAndIncrement(); nsTile < nStiles; nsTile = ai.getAndIncrement()) {
 //						int dl =  (nsTile == debug_stile) ? 3 : 0;
                         int dl = ((debugLevel > 1) && (nsTile == debug_stile)) ? 3: debugLevel;
-						if (dl > 1){
+						if (dl > 2){
 							System.out.println("getPlaneDispStrengths(): nsTile="+nsTile);
 						}
 						int stileY = nsTile / stilesX;  
@@ -1457,7 +1461,48 @@ public class SuperTiles{
 								measuredLayers,     // MeasuredLayers measuredLayers,
 								plPreferDisparity);   // boolean preferDisparity)
 
-						
+						if (smplMode && (world_plane_norm != null)) {
+							plane_tilts = new double [measuredLayers.getNumLayers()][][];
+							for (int ml = 0; ml < plane_tilts.length; ml++) if ((stMeasSel & ( 1 << ml)) != 0){
+								double [][] tile_disp_strengths = measuredLayers.getDisparityStrength(
+										ml,             // int num_layer,
+										stileX,         // int stX,
+										stileY,         // int stY,
+										null,           // boolean [] sel_in,
+										strength_floor, // double strength_floor,
+										strength_pow,   // double strength_pow,
+										true);          // boolean null_if_none);
+								if (tile_disp_strengths == null){
+									plane_tilts[ml] = 	zero_tilts;								
+								} else {
+									plane_tilts[ml] = pd0.getDisparityTilts(
+											world_plane_norm,     // double []     world_normal_xyz,
+											tile_disp_strengths,  // double [][]   tile_disp_strengths,
+											debugLevel);         // int           debugLevel);
+									if(dl > 2){
+										if (plane_tilts[ml] != null){
+											for (int k = 0; k < 2; k++) {
+												System.out.println((k > 0) ? "tilt Y, pix/tile" : "tilt X, pix/tile");
+												for (int i = 0; i < 2 * superTileSize; i++){
+													System.out.print(i+",");
+													for (int j = 0; j < 2 * superTileSize; j++){
+														int it = j + 2 * superTileSize * i;
+														if (plane_tilts[ml][it] !=null){
+															System.out.print(String.format("%f7.4", plane_tilts[ml][it][k]));
+														}
+														if (j < (2 * superTileSize - 1)){
+															System.out.print(", ");
+														}
+													}
+													System.out.println();
+												}
+												System.out.println();
+											}
+										}
+									}
+								}
+							}
+						}						
 
 						plane_disp_strength[nsTile] = new double[measuredLayers.getNumLayers()][][];
 						
@@ -1469,13 +1514,18 @@ public class SuperTiles{
 										stileX,         // int stX,
 										stileY,         // int stY,
 										null,           // boolean [] sel_in,
+//										((world_plane_norm == null)? zero_tilts: null),      // double []  tiltXY, // null - free with limit on both absolute (2.0?) and relative (0.2) values 
+										((world_plane_norm == null)? zero_tilts: plane_tilts[ml]),      // double []  tiltXY, // null - free with limit on both absolute (2.0?) and relative (0.2) values 
 										strength_floor, // double strength_floor,
 										strength_pow,   // double strength_pow,
 										smplSide,       // int        smplSide, // = 2;   // Sample size (side of a square)
 										smplNum,        //int        smplNum,   // = 3;   // Number after removing worst (should be >1)
 										smplRms,        //double     smplRms,   // = 0.1; // Maximal RMS of the remaining tiles in a sample
-										true,           // boolean null_if_none);
+//										true,           // boolean null_if_none);
 										smplWnd,        // boolean         smplWnd, //
+										2.0,            // double     max_abs_tilt, //  = 2.0; // pix per tile
+										0.2,            // double     max_rel_tilt, //  = 0.2; // (pix / disparity) per tile
+										true,           // boolean null_if_none);
 										dl);
 							} else {
 								plane_disp_strength[nsTile][ml] =  measuredLayers.getDisparityStrength(
@@ -1505,7 +1555,7 @@ public class SuperTiles{
 								if (dl > 0) {
 									System.out.println("Plane tilted disparity for stileX = "+stileX+" stileY="+stileY+", average disparity "+(sd/sw));
 								}
-								if (dl>1) {
+								if (dl>2) {
 									String [] dbg_titles = showSupertileSeparationTitles( plane_disp_strength[nsTile], null);
 									double [][] dbg_img = showSupertileSeparation(false, plane_disp_strength[nsTile], null); // plane_sels);
 									showDoubleFloatArrays sdfa_instance = new showDoubleFloatArrays();
@@ -1519,7 +1569,7 @@ public class SuperTiles{
 										null, // boolean [][]  tile_sel, // null - do not use, {} use all (will be modified)
 										plane_disp_strength[nsTile], // double [][][] disp_str, // calculate just once if null 
 										dl); // 1); // int           debugLevel);
-								if (dl>1) {
+								if (dl>2) {
 									String [] dbg_titles = showSupertileSeparationTitles( plane_disp_strength[nsTile], null);
 									double [][] dbg_img = showSupertileSeparation(false, plane_disp_strength[nsTile], null); // plane_sels);
 									showDoubleFloatArrays sdfa_instance = new showDoubleFloatArrays();
@@ -2623,7 +2673,7 @@ public class SuperTiles{
 				highMix,    //final double     highMix,    // stHighMix         = 0.4;   // Consider merging initial planes if jumps between ratio above
 				world_hor, // final double []  world_hor, // horizontal plane normal (default [0.0, 1.0, 0.0])
 				show_histograms, // final boolean    show_histograms,
-				debugLevel+0, // final int        debugLevel,
+				debugLevel+1, // final int        debugLevel,
 				dbg_X, // final int        dbg_X,
 				dbg_Y); // final int        dbg_Y)
 

src/main/java/TilePlanes.java

@@ -972,6 +972,82 @@ public class TilePlanes {
 			return true;
 		}
 		
+		/**
+		 * This method is used to filter sample plates by averaging subset of the square
+		 * sample and removing outliers. Currently only constant disparity and horizontal
+		 * surfaces are used, this method is used for horizontal ones to find tilts
+		 * d_disp/d_tx, d_disp/d_ty measured in pixels per tile
+		 * @param world_normal_xyz world space normal vector, currently only "up" - (0,1,0) is used
+		 * @param tile_disp_strengths [0] - unfiltered disparities for the supertile,
+		 *                            [1] - unfiltered correlation strengths for the
+		 *                                  supertile (just 0/non-0)
+		 * @param debugLevel
+		 * @return per tile arrays of either nulls or tilt-x, tilt-y pairs
+		 */
+		
+		public double [][] getDisparityTilts(
+				double []     world_normal_xyz,
+				double [][]   tile_disp_strengths,
+				int           debugLevel)
+		{
+			final int stSize2 = 2* stSize;
+			final int stCenter = (stSize2 + 1) * superTileSize / 2; 
+			final double [][] tile_tilts = new double [stSize2*stSize2][];
+			final double [] apy_vector = {0.0, 0.0, 1.0};
+			final Matrix py_vector = new Matrix(apy_vector,3); 
+			invalidateCalculated();			
+			double [] px_py = getCenterPxPy();
+			Matrix normal_row = new Matrix(world_normal_xyz,1); // 1x3
+			Matrix normal_col = new Matrix(world_normal_xyz,3); // 3x1
+			
+			// find world coordinates of the center of tile intersection with the plane
+			for (int lTile = 0; lTile < tile_disp_strengths[1].length; lTile++) if (tile_disp_strengths[1][lTile] > 0.0){
+				int tY = lTile / stSize2;
+				int tX = lTile % stSize2;
+				double px = px_py[0] + tX - stCenter;
+				double py = px_py[1] + tY - stCenter;
+				double disp = tile_disp_strengths[0][lTile];
+				// get world coordinates for each tile as determined individually
+				Matrix t_xyz = new Matrix(geometryCorrection.getWorldCoordinates(
+						px,
+						px,
+						disp,
+						this.correctDistortions),3);
+//				double n_by_w =  normal_row.times(t_xyz).get(0, 0);
+				// take pixel vector parallel to py and convert it to the world space
+				Matrix jacobian =  new Matrix(geometryCorrection.getWorldJacobian(
+						px,
+						py,
+						disp,
+						this.correctDistortions,
+						(debugLevel > 2)
+						));
+				Matrix inv_jacobian = jacobian.inverse();
+				Matrix wpy_vector = jacobian.times(py_vector); // 3 rows, 1 column py vector in the real world
+				// get 2 in-plane vectors in the real world
+				Matrix  wv1 = cross3d(normal_col, wpy_vector); 
+				Matrix  wv2 = cross3d(normal_col, wv1);
+				
+				// convert then to pixel space
+				Matrix pv1 = inv_jacobian.times(wv1); 
+				Matrix pv2 = inv_jacobian.times(wv2);
+				
+				// Get plane normal in the pixel space
+				Matrix pn = cross3d(pv1, pv2);
+				
+				// convert to the two tilts (d/dpx, d/dpy). We will need in pix/tile, not pix/pix
+				double [] txty = {
+						-pn.get(1, 0)/pn.get(0, 0)*stSize,
+						-pn.get(2, 0)/pn.get(0, 0)*stSize};
+				tile_tilts[lTile] = txty;
+			}
+			
+//			if (debugLevel > 1) {
+//				System.out.println("st_xyz = {"+st_xyz.get(0, 0)+","+st_xyz.get(1, 0)+","+st_xyz.get(2, 0)+"}"+" ="+n_by_w);
+//			}
+			return tile_tilts;
+			
+		}
 		/**
 		 * Tilt disparity values around the supertile center (this.zxy) so constant disparity in the output
 		 * corresponds to the real world plane parallel to the provided one. Used to discriminate tiles by
@@ -995,7 +1071,7 @@ public class TilePlanes {
 		{
 			final int stSize2 = 2* stSize;
 			invalidateCalculated();			
-			double [] px_py = getCenterPxPy();
+			double [] px_py = getCenterPxPy(); // tile center
 			// find world coordinates of the center of tile intersection with the plane
 			Matrix st_xyz = new Matrix(geometryCorrection.getWorldCoordinates(
 						px_py[0],
@@ -2251,7 +2327,8 @@ public class TilePlanes {
 			double val = values[0];
 			if ((dispNorm > 0.0) && (zxy[0] > dispNorm)){
 				double k = dispNorm / zxy[0];
-				val *= k*k;
+//				val *= k*k;
+				val *= k ; // reducing correction for the large disparities (making it sqrt() )
 			}
 			return val;
 		}
@@ -3671,6 +3748,13 @@ public class TilePlanes {
 				int debugLevel){
 			return getWorldXYZ(this.correctDistortions, debugLevel);
 		}
+		/**
+		 * Get vector from the camera perpendicular to the plane converted to the real world 
+		 * @param correct_distortions
+		 * @param debugLevel
+		 * @return
+		 */
+		
 		public double [] getWorldXYZ(
 				boolean correct_distortions,
 				int debugLevel)