handling zero eigenvalues

src/main/java/EyesisCorrectionParameters.java

@@ -2142,6 +2142,7 @@ public class EyesisCorrectionParameters {
   		public double     grow_disp_step =   6.0; // Increase disparity (from maximal tried) if nothing found in that tile // TODO: handle enclosed dips?  
   		public double     grow_min_diff =    0.5; // Grow more only if at least one channel has higher variance from others for the tile  
   		
+  		public boolean    plPreferDisparity    =   false;// Always start with disparity-most axis (false - lowest eigenvalue)
   		public double     plDispNorm           =   3.0;  // Normalize disparities to the average if above (now only for eigenvalue comparison)
   		public int        plMinPoints          =     5;  // Minimal number of points for plane detection
   		public double     plTargetEigen        =   0.1;  // Remove outliers until main axis eigenvalue (possibly scaled by plDispNorm) gets below
@@ -2150,7 +2151,8 @@ public class EyesisCorrectionParameters {
   		public double     plMinStrength        =   0.1;  // Minimal total strength of a plane 
   		public double     plMaxEigen           =   0.3;  // Maximal eigenvalue of a plane 
   		public boolean    plDbgMerge           =   true; // Combine 'other' plane with current
-  		public double     plWorstWorsening     =   1.0;  // Worst case worsening after merge 
+  		public double     plWorstWorsening     =   3.0;  // Worst case worsening after merge
+  		public double     plWeakWorsening      =   1.0;  // Relax merge requirements for weaker planes
   		public boolean    plMutualOnly         =   true; // keep only mutual links, remove weakest if conflict
   		public boolean    plFillSquares        =   true; // Add diagonals to full squares
   		public boolean    plCutCorners         =   true; // Add ortho to 45-degree corners
@@ -2406,6 +2408,7 @@ public class EyesisCorrectionParameters {
 			properties.setProperty(prefix+"grow_disp_step",   this.grow_disp_step +"");
 			properties.setProperty(prefix+"grow_min_diff",    this.grow_min_diff +"");
 
+			properties.setProperty(prefix+"plPreferDisparity",this.plPreferDisparity+"");
 			properties.setProperty(prefix+"plDispNorm",       this.plDispNorm +"");
   			properties.setProperty(prefix+"plMinPoints",      this.plMinPoints+"");
 			properties.setProperty(prefix+"plTargetEigen",    this.plTargetEigen +"");
@@ -2415,6 +2418,7 @@ public class EyesisCorrectionParameters {
 			properties.setProperty(prefix+"plMaxEigen",       this.plMaxEigen +"");
 			properties.setProperty(prefix+"plDbgMerge",       this.plDbgMerge+"");
 			properties.setProperty(prefix+"plWorstWorsening", this.plWorstWorsening +"");
+			properties.setProperty(prefix+"plWeakWorsening",  this.plWeakWorsening +"");
 			properties.setProperty(prefix+"plMutualOnly",     this.plMutualOnly+"");
 
 			properties.setProperty(prefix+"plFillSquares",    this.plFillSquares+"");
@@ -2657,6 +2661,7 @@ public class EyesisCorrectionParameters {
   			if (properties.getProperty(prefix+"grow_disp_step")!=null)    this.grow_disp_step=Double.parseDouble(properties.getProperty(prefix+"grow_disp_step"));
   			if (properties.getProperty(prefix+"grow_min_diff")!=null)     this.grow_min_diff=Double.parseDouble(properties.getProperty(prefix+"grow_min_diff"));
 
+  			if (properties.getProperty(prefix+"plPreferDisparity")!=null) this.plPreferDisparity=Boolean.parseBoolean(properties.getProperty(prefix+"plPreferDisparity"));
   			if (properties.getProperty(prefix+"plDispNorm")!=null)        this.plDispNorm=Double.parseDouble(properties.getProperty(prefix+"plDispNorm"));
   			if (properties.getProperty(prefix+"plMinPoints")!=null)       this.plMinPoints=Integer.parseInt(properties.getProperty(prefix+"plMinPoints"));
   			if (properties.getProperty(prefix+"plTargetEigen")!=null)     this.plTargetEigen=Double.parseDouble(properties.getProperty(prefix+"plTargetEigen"));
@@ -2666,6 +2671,7 @@ public class EyesisCorrectionParameters {
   			if (properties.getProperty(prefix+"plMaxEigen")!=null)        this.plMaxEigen=Double.parseDouble(properties.getProperty(prefix+"plMaxEigen"));
   			if (properties.getProperty(prefix+"plDbgMerge")!=null)        this.plDbgMerge=Boolean.parseBoolean(properties.getProperty(prefix+"plDbgMerge"));
   			if (properties.getProperty(prefix+"plWorstWorsening")!=null)  this.plWorstWorsening=Double.parseDouble(properties.getProperty(prefix+"plWorstWorsening"));
+  			if (properties.getProperty(prefix+"plWeakWorsening")!=null)  this.plWeakWorsening=Double.parseDouble(properties.getProperty(prefix+"plWeakWorsening"));
   			if (properties.getProperty(prefix+"plMutualOnly")!=null)      this.plMutualOnly=Boolean.parseBoolean(properties.getProperty(prefix+"plMutualOnly"));
 
   			if (properties.getProperty(prefix+"plFillSquares")!=null)     this.plFillSquares=Boolean.parseBoolean(properties.getProperty(prefix+"plFillSquares"));
@@ -2930,6 +2936,7 @@ public class EyesisCorrectionParameters {
   			gd.addNumericField("Increase disparity (from maximal tried) if nothing found in that tile",        this.grow_disp_step,  6);
   			gd.addNumericField("Grow more only if at least one channel has higher variance from others for the tile", this.grow_min_diff,  6);
   			gd.addMessage     ("--- Planes detection ---");
+  			gd.addCheckbox    ("Always start with disparity-most axis (false - lowest eigenvalue)",            this.plPreferDisparity);
   			gd.addNumericField("Normalize disparities to the average if above",                                this.plDispNorm,  6);
   			gd.addNumericField("Minimal number of points for plane detection",                                  this.plMinPoints,  0);
   			gd.addNumericField("Remove outliers until main axis eigenvalue (possibly scaled by plDispNorm) gets below", this.plTargetEigen,  6);
@@ -2937,8 +2944,9 @@ public class EyesisCorrectionParameters {
   			gd.addNumericField("Maximal number of outliers to remove",                                         this.plMaxOutliers,  0);
   			gd.addNumericField("Minimal total strength of a plane",                                            this.plMinStrength,  6);
   			gd.addNumericField("Maximal eigenvalue of a plane",                                                this.plMaxEigen,  6);
-  			gd.addCheckbox    ("Combine 'other' plane with the current",                                       this.plDbgMerge);
+  			gd.addCheckbox    ("Combine 'other' plane with the current (unused)",                              this.plDbgMerge);
   			gd.addNumericField("Worst case worsening after merge",                                             this.plWorstWorsening,  6);
+  			gd.addNumericField("Relax merge requirements for weaker planes",                                   this.plWeakWorsening,  6);
   			gd.addCheckbox    ("Keep only mutual links, remove weakest if conflict",                           this.plMutualOnly);
 
   			gd.addCheckbox    ("Add diagonals to full squares",                                                this.plFillSquares);
@@ -3191,6 +3199,7 @@ public class EyesisCorrectionParameters {
   			this.grow_disp_step=        gd.getNextNumber();
   			this.grow_min_diff=         gd.getNextNumber();
 
+  			this.plPreferDisparity=     gd.getNextBoolean();
   			this.plDispNorm=            gd.getNextNumber();
   			this.plMinPoints=     (int) gd.getNextNumber();
   			this.plTargetEigen=         gd.getNextNumber();
@@ -3200,6 +3209,7 @@ public class EyesisCorrectionParameters {
   			this.plMaxEigen=            gd.getNextNumber();
   			this.plDbgMerge=            gd.getNextBoolean();
   			this.plWorstWorsening=      gd.getNextNumber();
+  			this.plWeakWorsening=       gd.getNextNumber();
   			this.plMutualOnly=          gd.getNextBoolean();
 
   			this.plFillSquares=         gd.getNextBoolean();

src/main/java/TilePlanes.java

@@ -452,6 +452,7 @@ public class TilePlanes {
 				PlaneData otherPd,
 				double    scale_other,
 				boolean   ignore_weights,
+				boolean   preferDisparity, // Always start with disparity-most axis (false - lowest eigenvalue)
 				int       debugLevel)
 		{
 			if (debugLevel > 0) {
@@ -474,7 +475,11 @@ public class TilePlanes {
 			double sum_weight =      scale_other *  otherPd.weight + this.weight;
 			double other_fraction = ignore_weights? (scale_other/(scale_other + 1.0)): ((scale_other *  otherPd.weight) / sum_weight); 
 			Matrix common_center =  this_center.times(1.0 - other_fraction).plus(other_center.times(other_fraction));
-			Matrix other_offset =   other_center.minus(this_center); // other center from this center 
+			Matrix other_offset =   other_center.minus(this_center); // other center from this center
+			if ((this.values[0] == 0.0) || (otherPd.values[0] == 0.0)) {
+				System.out.println("Zero eigenvalue");
+				debugLevel = 10;
+			}
 			if (debugLevel > 0) {
 				System.out.println("other_eig_vals");
 				other_eig_vals.print(8, 6);
@@ -495,9 +500,6 @@ public class TilePlanes {
 				System.out.println("other_fraction="+other_fraction);
 			}
 			
-			
-			
-			
 			double [][] acovar = { // covariance matrix of center masses (not yet scaled by weight)
 					{other_offset.get(0,0)*other_offset.get(0,0), other_offset.get(0,0)*other_offset.get(1,0), other_offset.get(0,0)*other_offset.get(2,0)},
 					{other_offset.get(1,0)*other_offset.get(0,0), other_offset.get(1,0)*other_offset.get(1,0), other_offset.get(1,0)*other_offset.get(2,0)},
@@ -525,9 +527,12 @@ public class TilePlanes {
 				System.out.println("covar with other_covar and this_covar");
 				covar.print(8, 6);
 			}
-			
+			if (Double.isNaN(covar.get(0, 0))){
+				System.out.println("covar is NaN !");
+				covar.print(8, 6);
+			}
 			// extract new eigenvalues, eigenvectors
-			EigenvalueDecomposition eig = covar.eig();
+			EigenvalueDecomposition eig = covar.eig(); // verify NaN - it gets stuck
 //			eig.getD().getArray(),
 //			eig.getV().getArray(),
 			if (debugLevel > 0) {
@@ -542,12 +547,30 @@ public class TilePlanes {
 			double [][] eig_vect = eig.getV().getArray();
 			double [][] eig_val = eig.getD().getArray();
 			// make towards camera, left coordinate system
+/*			
 			int oindx = 0;
 			for (int i = 1; i <3; i++){
 				if (Math.abs(eig_vect[0][i]) > Math.abs(eig_vect[0][oindx])){
 					oindx = i;
 				}
 			}
+*/			
+			int oindx = 0;
+			if (preferDisparity) {
+				for (int i = 1; i <3; i++){
+					if (Math.abs(eig_vect[0][i]) > Math.abs(eig_vect[0][oindx])){
+						oindx = i;
+					}
+				}
+			} else {
+				for (int i = 1; i <3; i++){
+					if (eig_val[i][i] < eig_val[oindx][oindx]){
+						oindx = i;
+					}
+				}
+
+			}
+			
 			// select 2 other axes for increasing eigenvalues (so v is short axis, h  is the long one)
 			int vindx = (oindx == 0)? 1 : 0;
 			int hindx = (oindx == 0)? 2 : ((oindx == 1) ? 2 : 1);
@@ -886,6 +909,7 @@ public class TilePlanes {
 			boolean [] select,
 			double     plDispNorm, //  Normalize disparities to the average if above
 			int        debugLevel){
+		double mindet = 1E-15;
 		int stSize2 = 2 * stSize;
 //		Matrix covar = new Matrix(3,3);
 		double [][] acovar = new double [3][3];
@@ -939,8 +963,21 @@ public class TilePlanes {
 		acovar [2][0] = acovar [0][2]; 
 		acovar [2][1] = acovar [1][2];
 		Matrix covar = new Matrix(acovar);
+		if (Math.abs(covar.det()) < mindet){
+			debugLevel = 5;
+		}
+
 		EigenvalueDecomposition eig = covar.eig();
+		if (Double.isNaN(eig.getV().get(0, 0))){
+			System.out.println("getCovar(): Double.isNaN(eig.getV().get(0, 0))");
+			debugLevel = 20;
+		}
+		
+		if (eig.getD().get(0, 0) == 0.0){
+			debugLevel = 10;
+		}
 		if (debugLevel > 0){
+			System.out.println("getCovar(): sw = "+sw +", swz = "+swz +", swx = "+swx +", swy = "+swy +", covar.det() = "+covar.det());
 			System.out.println("getCovar(): covarianvce matrix, number of used points:"+numPoints);
 			covar.print(10, 6); // w,d
 			System.out.println("getCovar(): eigenvalues");
@@ -948,6 +985,10 @@ public class TilePlanes {
 			System.out.println("getCovar(): eigenvectors");
 			eig.getV().print(10, 6); // w,d
 		}
+		if ((eig.getD().get(0, 0) == 0.0) || (Math.abs(covar.det()) < mindet)) {
+			return null; // testing with zero eigenvalue
+			// Problem with zero eigenvalue is with derivatives and coordinate conversion
+		}
 		double [][][] rslt = {
 				eig.getD().getArray(),
 				eig.getV().getArray(),
@@ -962,13 +1003,12 @@ public class TilePlanes {
 			double []  data,
 			double []  weight,
 			boolean [] select, // null OK, will enable all tiles
+			boolean    preferDisparity, // Always start with disparity-most axis (false - lowest eigenvalue)
 			int        debugLevel){
 		if (select == null) {
 			select = new boolean [4*stSize];
 			for (int i = 0; i < select.length; i++) select[i] = true;
 		}
-//		int debugLevel1 = ((sTileXY[0] == 27) && (sTileXY[1] == 20))? 1: 0; 
-//		int debugLevel1 = ((sTileXY[0] == 27) && (sTileXY[1] == 17))? 1: 0; 
 		int debugLevel1 = ((sTileXY[0] == 27) && (sTileXY[1] == 16))? 1: 0; // check why v[0][0] <0  
 		
 		
@@ -980,7 +1020,6 @@ public class TilePlanes {
 				debugLevel1); //0); // debugLevel);
 		if (rslt == null) return null;
 		int       numPoints =  (int) rslt[2][0][2];
-//		double    kz =   rslt[2][0][1]; // == 1.0
 		double    swc =  rslt[2][0][0];
 		double [] szxy = rslt[2][1];
 		double [][] eig_val =  rslt[0];
@@ -988,11 +1027,23 @@ public class TilePlanes {
 		// find vector most orthogonal to view // (anyway it all works with that assumption), make it first
 		// TODO normalize to local linear scales
 		int oindx = 0;
-		for (int i = 1; i <3; i++){
-			if (Math.abs(eig_vect[0][i]) > Math.abs(eig_vect[0][oindx])){
-				oindx = i;
+		if (preferDisparity) {
+			for (int i = 1; i <3; i++){
+				if (Math.abs(eig_vect[0][i]) > Math.abs(eig_vect[0][oindx])){
+					oindx = i;
+				}
 			}
+		} else {
+			for (int i = 1; i < 3 ; i++){
+				if (eig_val[i][i] < eig_val[oindx][oindx]){
+					oindx = i;
+				}
+			}
+		}
+		if (eig_val[oindx][oindx] == 0.0){
+			System.out.println("getPlane(): zero eigenvalue!!");
 		}
+		
 		/*
 		// Find two other axis - "mostly X" (horizontal) and "mostly Y" (vertical) 
 		int vindx = (oindx == 0)? 1 : 0;
@@ -1058,6 +1109,7 @@ public class TilePlanes {
 			double     targetEigen, // target eigenvalue for primary axis (is disparity-dependent, so is non-constant)
 			int        maxRemoved,  // maximal number of tiles to remove (not a constant)
 			int        minLeft,     // minimal number of tiles to keep
+			boolean    preferDisparity, // Always start with disparity-most axis (false - lowest eigenvalue)
 			int        debugLevel){
 		int stSize2 = 2 * stSize;
 		if (pd == null) {
@@ -1066,10 +1118,14 @@ public class TilePlanes {
 					data,
 					weight,
 					select, // null OK
+					preferDisparity,
 					debugLevel);
 		} else if (select != null){
 			pd.setPlaneSelection(select);
 		}
+		if (pd == null){
+			return null; // zero eigenvalues
+		}
 		if (maxRemoved > (pd.getNumPoints() - minLeft)) maxRemoved = pd.getNumPoints() - minLeft;
 		int numRemoved = 0;
 		for (; (pd.getValue() > targetEigen) && (numRemoved < maxRemoved); numRemoved++){
@@ -1107,7 +1163,11 @@ public class TilePlanes {
 					data,
 					weight,
 					select,
+					preferDisparity,
 					debugLevel);
+			if (pd == null) {
+				return null;
+			}
 		}
 		return pd;
 	}

src/main/java/TileProcessor.java

@@ -2952,7 +2952,8 @@ public class TileProcessor {
 				clt_parameters.plMinPoints, //           =     5;  // Minimal number of points for plane detection
 				clt_parameters.plTargetEigen, //         =   0.1;  // Remove outliers until main axis eigenvalue (possibly scaled by plDispNorm) gets below
 				clt_parameters.plFractOutliers, //      =   0.3;  // Maximal fraction of outliers to remove
-				clt_parameters.plMaxOutliers, //        =    20;  // Maximal number of outliers to remove
+				clt_parameters.plMaxOutliers, //        =    20;  // Maximal number of outliers to remove\
+				clt_parameters.plPreferDisparity,
 				geometryCorrection,
 				clt_parameters.correct_distortions,
 				0, // -1, // debugLevel,                  // final int        debugLevel)
@@ -2964,26 +2965,36 @@ public class TileProcessor {
 		if (debugLevel > -1) sdfa_instance = new showDoubleFloatArrays(); // just for debugging?
 
 		st.matchPlanes(
+				clt_parameters.plPreferDisparity,
 				clt_parameters.tileX,
 				clt_parameters.tileY); 
 
 		st.selectNeighborPlanesMutual(
 				clt_parameters.plWorstWorsening, // final double worst_worsening,
+				clt_parameters.plWeakWorsening, // final double worst_worsening,
 				clt_parameters.plDispNorm,
 				clt_parameters.plMaxEigen,
 				clt_parameters.plMinStrength,
 				0, // final int debugLevel)
 				clt_parameters.tileX,
 				clt_parameters.tileY);
-		
-		if (clt_parameters.plFillSquares){
-			st.fillSquares();
-		}
-		if (clt_parameters.plCutCorners){
-			st.cutCorners();
+		while (true) {
+			int num_added = 0;
+			if (clt_parameters.plFillSquares){
+				num_added += st.fillSquares();
+			}
+			if (debugLevel > -1) {
+				System.out.println("after fillSquares() added "+num_added);
+			}
+			if (clt_parameters.plCutCorners){
+				num_added += st.cutCorners();
+			}
+			if (debugLevel > -1) {
+				System.out.println("after plCutCorners() added (cumulative) "+num_added);
+			}
+			if (num_added == 0) break;
 		}
 		
-		
 		TilePlanes.PlaneData [][] planes_mod = null;
 
 		// smooth planes (by averaging with neighbors and the "measured" one with variable "pull")  
@@ -2993,6 +3004,7 @@ public class TileProcessor {
 					clt_parameters.plPull,                        // final double      meas_pull,//  relative pull of the original (measured) plane with respect to the average of the neighbors
 					clt_parameters.plIterations,                  // final int         num_passes,
 					Math.pow(10.0,  -clt_parameters.plPrecision), // final double      maxDiff, // maximal change in any of the disparity values
+					clt_parameters.plPreferDisparity,
 					0, // final int debugLevel)
 					clt_parameters.tileX,
 					clt_parameters.tileY);