redesigning shell extraction

src/main/java/EyesisCorrectionParameters.java

@@ -2155,6 +2155,13 @@ public class EyesisCorrectionParameters {
   		
   		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 double     plBlurBinVert        =   1.2;  // Blur disparity histograms for constant disparity clusters by this sigma (in bins)
+  		public double     plBlurBinHor         =   0.8;  // Blur disparity histograms for horizontal clusters by this sigma (in bins)
+  		public double     plMaxDiffVert        =   0.4;  // Maximal normalized disparity difference when initially assigning to vertical plane
+  		public double     plMaxDiffHor         =   0.2;  // Maximal normalized disparity difference when initially assigning to horizontal plane
+  		public int        plInitPasses         =     3;  // Number of initial passes to assign tiles to vert (const disparity) and hor planes
+  		
   		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
   		public double     plFractOutliers      =   0.3;  // Maximal fraction of outliers to remove
@@ -2453,7 +2460,14 @@ public class EyesisCorrectionParameters {
 
 			properties.setProperty(prefix+"plPreferDisparity",this.plPreferDisparity+"");
 			properties.setProperty(prefix+"plDispNorm",       this.plDispNorm +"");
-  			properties.setProperty(prefix+"plMinPoints",      this.plMinPoints+"");
+
+			properties.setProperty(prefix+"plBlurBinVert",    this.plBlurBinVert +"");
+			properties.setProperty(prefix+"plBlurBinHor",     this.plBlurBinHor +"");
+			properties.setProperty(prefix+"plMaxDiffVert",    this.plMaxDiffVert +"");
+			properties.setProperty(prefix+"plMaxDiffHor",     this.plMaxDiffHor +"");
+  			properties.setProperty(prefix+"plInitPasses",     this.plInitPasses+"");
+			
+			properties.setProperty(prefix+"plMinPoints",      this.plMinPoints+"");
 			properties.setProperty(prefix+"plTargetEigen",    this.plTargetEigen +"");
 			properties.setProperty(prefix+"plFractOutliers",  this.plFractOutliers +"");
   			properties.setProperty(prefix+"plMaxOutliers",    this.plMaxOutliers+"");
@@ -2463,10 +2477,8 @@ public class EyesisCorrectionParameters {
 			properties.setProperty(prefix+"plWorstWorsening", this.plWorstWorsening +"");
 			properties.setProperty(prefix+"plWeakWorsening",  this.plWeakWorsening +"");
 			properties.setProperty(prefix+"plMutualOnly",     this.plMutualOnly+"");
-
 			properties.setProperty(prefix+"plFillSquares",    this.plFillSquares+"");
 			properties.setProperty(prefix+"plCutCorners",     this.plCutCorners+"");
-
 			properties.setProperty(prefix+"plPull",           this.plPull +"");
 			properties.setProperty(prefix+"plNormPow",        this.plNormPow +"");
   			properties.setProperty(prefix+"plIterations",     this.plIterations+"");
@@ -2738,6 +2750,13 @@ public class EyesisCorrectionParameters {
 
   			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+"plBlurBinVert")!=null)     this.plBlurBinVert=Double.parseDouble(properties.getProperty(prefix+"plBlurBinVert"));
+  			if (properties.getProperty(prefix+"plBlurBinHor")!=null)      this.plBlurBinHor=Double.parseDouble(properties.getProperty(prefix+"plBlurBinHor"));
+  			if (properties.getProperty(prefix+"plMaxDiffVert")!=null)     this.plMaxDiffVert=Double.parseDouble(properties.getProperty(prefix+"plMaxDiffVert"));
+  			if (properties.getProperty(prefix+"plMaxDiffHor")!=null)      this.plMaxDiffHor=Double.parseDouble(properties.getProperty(prefix+"plMaxDiffHor"));
+  			if (properties.getProperty(prefix+"plInitPasses")!=null)      this.plInitPasses=Integer.parseInt(properties.getProperty(prefix+"plInitPasses"));
+  			
   			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"));
   			if (properties.getProperty(prefix+"plFractOutliers")!=null)   this.plFractOutliers=Double.parseDouble(properties.getProperty(prefix+"plFractOutliers"));
@@ -2746,7 +2765,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+"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"));
@@ -3047,7 +3066,14 @@ public class EyesisCorrectionParameters {
   			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("Blur disparity histograms for constant disparity clusters by this sigma (in bins)",   this.plBlurBinVert,  6);
+  			gd.addNumericField("Blur disparity histograms for horizontal clusters by this sigma (in bins)",           this.plBlurBinHor,  6);
+  			gd.addNumericField("Maximal normalized disparity difference when initially assigning to vertical plane",  this.plMaxDiffVert,  6);
+  			gd.addNumericField("Maximal normalized disparity difference when initially assigning to horizontal plane",this.plMaxDiffHor,  6);
+  			gd.addNumericField("Number of initial passes to assign tiles to vert (const disparity) and hor planes",   this.plInitPasses,  0);
+  			
+  			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);
   			gd.addNumericField("Maximal fraction of outliers to remove",                                       this.plFractOutliers,  6);
   			gd.addNumericField("Maximal number of outliers to remove",                                         this.plMaxOutliers,  0);
@@ -3340,6 +3366,13 @@ public class EyesisCorrectionParameters {
 
   			this.plPreferDisparity=     gd.getNextBoolean();
   			this.plDispNorm=            gd.getNextNumber();
+
+  			this.plBlurBinVert=         gd.getNextNumber();
+  			this.plBlurBinHor=          gd.getNextNumber();
+  			this.plMaxDiffVert=         gd.getNextNumber();
+  			this.plMaxDiffHor=          gd.getNextNumber();
+  			this.plInitPasses=    (int) gd.getNextNumber();
+  			
   			this.plMinPoints=     (int) gd.getNextNumber();
   			this.plTargetEigen=         gd.getNextNumber();
   			this.plFractOutliers=       gd.getNextNumber();

src/main/java/TilePlanes.java

@@ -165,6 +165,10 @@ public class TilePlanes {
 			}
 			return 	this.measuredSelection[nl];
 		}
+
+		public boolean [][] getMeasSelection(){
+			return 	this.measuredSelection;
+		}
 		
 		
 		public MeasuredLayers getMeasuredLayers()
@@ -1006,10 +1010,8 @@ public class TilePlanes {
 				debugLevel = 20;
 			}
 			
-//			if (eig.getD().get(0, 0) == 0.0){
-//				debugLevel = 10;
-//			}
-			if (debugLevel > 0){
+			if (debugLevel > 3){
+//				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:"+num_tiles);
 				covar.print(10, 6); // w,d
@@ -1022,15 +1024,6 @@ public class TilePlanes {
 				return null; // testing with zero eigenvalue
 				// Problem with zero eigenvalue is with derivatives and coordinate conversion
 			}
-/*			
-			double [][][] rslt = {
-					eig.getD().getArray(),
-					eig.getV().getArray(),
-					{
-						{sw,kz,numPoints},
-						{swz, swx, swy}}};
-			return rslt;
-*/
 
 			double [][] eig_val =  eig.getD().getArray(); // rslt[0];
 			double [][] eig_vect = eig.getV().getArray(); // rslt[1];
@@ -1363,13 +1356,13 @@ public class TilePlanes {
 				double y = tileSize * (sy + 0.5) + 0.5  - zxy[2];
 				for (int sx = -3 * superTileSize/2; sx < 3 * superTileSize / 2; sx++){
 					double x = tileSize * (sx + 0.5) + 0.5 - zxy[1];
-						disparities[indx] = zxy[0] - (normal[1] * x + normal[2] * y)/normal[0];
+					disparities[indx] = zxy[0] - (normal[1] * x + normal[2] * y)/normal[0];
 					indx++;
 				}
 			}
 			return disparities;
 		}
-		
+
 		public double[] getTriplePlaneDisparity(
 				int     dir)
 		{
@@ -1389,13 +1382,79 @@ public class TilePlanes {
 			}
 			return plane;
 		}
-		
-		
-//		double px = tileSize*(superTileSize * sTileXY[0] + superTileSize/2) + zxy[1];  // [3] - plane point {disparity, x, y), x=0, y=0 is a 4,4 point of an 8x8 supertile
-//		double py = tileSize*(superTileSize * sTileXY[1] + superTileSize/2) + zxy[2];
-		
-		
-		
+
+		/**
+		 * Get disparity values for the tiles of this overlapping supertile as [2*superTileSize * 2*superTileSize] array
+		 * and weights combined from provided window function, optional selection and using ellipsoid projection on the
+		 * px, py plane (constant disparity
+		 * @param window null or window function as [2*superTileSize * 2*superTileSize] array
+		 * @param use_sel use plane selection (this.sel_mask) to select only some part of the plane
+		 * @param scale_projection use plane ellipsoid projection for weight: 0 - do not use, > 0 linearly scale ellipsoid 
+		 * @return a pair of ar5rays {disparity, strength}, each [2*superTileSize * 2*superTileSize]
+		 */
+		public double[][] getDoublePlaneDisparityStrength(
+				double [] window,
+				boolean   use_sel,
+				double    scale_projection,
+				int       debugLevel)
+		{
+			double [][] disp_strength = new double[2][4*superTileSize*superTileSize];
+			int indx = 0;
+			double [] normal = getVector();
+			double [] zxy =    getZxy(); // {disparity, x center in pixels, y center in pixels (relative to a supertile center)
+			double    weight = getWeight();
+			double k_gauss = 0;
+			Matrix val2d = null, vect2d = null;
+			if (scale_projection > 0.0){
+				double [] vals3d =      getValues();
+				double [][] vectors3d = getVectors();
+				double [][] acovar = new double [2][2];
+				for (int i = 0; i < 2; i++){
+					for (int j = i; j < 2; j++){
+						acovar[i][j] = 0.0;
+						for (int k = 0; k < 3; k++){
+							acovar[i][j] += vals3d[k] * vectors3d[k][i+1] * vectors3d[k][j+1]; // 0 - z, disparity == 0
+						}
+						if (i != j) {
+							acovar[j][i] =acovar[i][j];
+						}
+					}
+				}
+				Matrix covar = new Matrix(acovar); // 2d, x y only
+				EigenvalueDecomposition eig = covar.eig();
+				val2d = eig.getD();
+				vect2d = eig.getV().transpose();
+				k_gauss = 0.5/(scale_projection*scale_projection);
+			}
+
+			for (int sy = -superTileSize; sy < superTileSize; sy++){
+				// adding half-tile and half-pixel to match the center of the pixel. Supertile center is between
+				// pixel 31 and pixel 32 (counting from 0) in both directions
+				double y = tileSize * (sy + 0.5) + 0.5  - zxy[2];
+				for (int sx = -superTileSize; sx < superTileSize; sx++){
+					double x = tileSize * (sx + 0.5) + 0.5 - zxy[1];
+					disp_strength[0][indx] = zxy[0] - (normal[1] * x + normal[2] * y)/normal[0];
+					double w = weight;
+					if (window != null) w *= window[indx];
+					if (use_sel && (sel_mask != null) && !(sel_mask[indx])) w = 0.0;
+					if ((w > 0.0) && (scale_projection > 0.0)){
+						double [] xy = {x,y};
+						Matrix vxy = vect2d.times(new Matrix(xy,2)); // verify if it is correct
+						double r2 = 0; 
+						for (int i = 0; i <2; i++){
+							double d = vxy.get(i,0);
+							r2 += d * d / val2d.get(i, i);
+						}
+						w *= Math.exp(-k_gauss*r2); // verify it is correct size - maybe it should be -0.5*r2 ?
+					}					
+					disp_strength[1][indx] = w;
+					indx++;
+				}
+			}
+			return disp_strength;
+		}
+
+
 		/**
 		 * Cross product of 2 3-d vectors as column matrices
 		 * @param v1
@@ -1412,8 +1471,8 @@ public class TilePlanes {
 					{av1[0][0]*av2[1][0] - av1[1][0]*av2[0][0]}};
 			return new Matrix (ar);
 		}
-		
-		
+
+
 		/**
 		 * Combine 2 Plane instances using centers, eigenvalues eihenvectors and total weights of this and other PlaneData objects
 		 * other plane should already be transformed to the same supertile coordinate system (with getPlaneToThis() method)
@@ -1476,7 +1535,7 @@ public class TilePlanes {
 				other_offset.print(8, 6);
 				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)},
@@ -1493,7 +1552,7 @@ public class TilePlanes {
 				System.out.println("covar");
 				covar.print(8, 6);
 			}			
-			
+
 			covar.plusEquals(other_covar.times(other_fraction));
 			if (debugLevel > 0) {
 				System.out.println("covar with other_covar");
@@ -1510,8 +1569,8 @@ public class TilePlanes {
 			}
 			// extract new eigenvalues, eigenvectors
 			EigenvalueDecomposition eig = covar.eig(); // verify NaN - it gets stuck
-//			eig.getD().getArray(),
-//			eig.getV().getArray(),
+			//			eig.getD().getArray(),
+			//			eig.getV().getArray(),
 			if (debugLevel > 0) {
 				System.out.println("eig.getV()");
 				eig.getV().print(8, 6);
@@ -1519,19 +1578,19 @@ public class TilePlanes {
 				eig.getD().print(8, 6);
 			}
 
-			
-			
+
+
 			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++){
@@ -1547,7 +1606,7 @@ public class TilePlanes {
 				}
 
 			}
-			
+
 			// 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);
@@ -1565,18 +1624,18 @@ public class TilePlanes {
 			for (int v = 0; v < 2; v++) {
 				if (plane[v][v] < 0.0) for (int i = 0; i < 3; i ++) plane[v][i] = -plane[v][i];
 			}
-			
+
 			// make  direction last vector so px (x) py (.) disp < 0 (left-hand coordinate system) 
 			if (new Matrix(plane).det() > 0){
 				for (int i = 0; i < 3; i ++) plane[2][i] = -plane[2][i];
 			}
-			
-			
+
+
 			PlaneData pd = this.clone(); // will copy selections too
 			pd.invalidateCalculated();   // real world vectors
 			pd.setValues(eig_val[oindx][oindx],eig_val[vindx][vindx],eig_val[hindx][hindx]); // eigenvalues [0] - thickness, 2 other to detect skinny (poles)
 			pd.setVectors(plane);
-			
+
 			pd.setZxy(common_center.getColumnPackedCopy()); // set new center
 			// what weight to use? cloned is original weight for this supertile
 			// or use weighted average like below?
@@ -1588,7 +1647,7 @@ public class TilePlanes {
 
 			return pd;
 		}
-		
+
 		
 		/**
 		 * Convert plane data from other supertile to this one (disparity, px, py) for the center of this supertile
@@ -1956,11 +2015,7 @@ public class TilePlanes {
 			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){
+		if (debugLevel > 3){
 			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

src/main/java/TileSurface.java

+/**
+ **
+ ** TileSurface - hadle tile surfaces
+ **
+ ** Copyright (C) 2017 Elphel, Inc.
+ **
+ ** -----------------------------------------------------------------------------**
+ **  
+ **  TileSurface.java is free software: you can redistribute it and/or modify
+ **  it under the terms of the GNU General Public License as published by
+ **  the Free Software Foundation, either version 3 of the License, or
+ **  (at your option) any later version.
+ **
+ **  This program is distributed in the hope that it will be useful,
+ **  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ **  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ **  GNU General Public License for more details.
+ **
+ **  You should have received a copy of the GNU General Public License
+ **  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ ** -----------------------------------------------------------------------------**
+ **
+ */
+
+public class TileSurface {
+//	public 
+		private int tileSize;
+		private int stSize;
+		private int tilesX;
+		private int tilesY;
+		private int [] st_dirs8;
+		
+//		private int nsTilesstSize =   0; // 8;
+		GeometryCorrection   geometryCorrection = null;
+		public TileSurface(
+				int tileSize,
+				int stSize,
+				int tilesX,
+				int tilesY,
+				GeometryCorrection geometryCorrection){
+			this.tileSize = tileSize;
+			this.stSize = stSize;
+			this.geometryCorrection =geometryCorrection;
+			this.tilesX =  tilesX;
+			this.tilesY =  tilesY;
+			int [] dirs =  {-tilesX, -tilesX + 1, 1, tilesX + 1, tilesX, tilesX - 1, -1, -tilesX - 1};
+			this.st_dirs8 = dirs;
+		}
+
+		public class TileData{
+			double disparity;
+			double strength;
+			double enable;
+			int [] neighbors = null;
+		}
+		/**
+		 * Get tile surface number from supertile number, direction (-1 same) and the supertile plane index
+		 * @param nsTile number of the supertile
+		 * @param dir direction -1 - same supertile, 0 - N (up), 1 - NE, .. 7 - NW
+		 * @param np number of the supertile plane
+		 * @param planes array of the per-supertile, per plane plane data (each level can be null)
+		 * @return unique tile plane index
+		 */
+		public int getTileSurfaceNumber (
+				int nsTile,
+				int dir,              // direction, or -1 (same)
+				int np,
+				TilePlanes.PlaneData [][] planes
+				){
+			if ((planes[nsTile] == null) || (planes[nsTile][np] == null)){
+				return -1; // empty supertile or supertile plane
+			}
+			if (dir < 0) {
+				return np;
+			}
+			int tsn = planes[nsTile].length;
+			for (int d = 0; d < dir; d ++){
+				int nsTile1 = nsTile + st_dirs8[d];
+				if (planes[nsTile1] != null){
+					tsn += planes[nsTile1].length;
+				}
+			}
+			return tsn + np; 
+		}
+		/**
+		 * Get supertile direction and the plane number that contributeted to a specific tile surface 
+		 * @param nsTile supertile index
+		 * @param tp tile surface index (generated by getTileSurfaceNumber)
+		 * @param planes array of the per-supertile, per plane plane data (each level can be null)
+		 * @return a pair of {dir, plane index}. dir is -1 for the plane in the same supertile, 0..7 for neighbors
+		 */
+		
+		public int []  getSuperTileDirPlane (
+				int nsTile,
+				int tp,
+				TilePlanes.PlaneData [][] planes
+				)
+		{
+			int num_planes = (planes[nsTile] == null)? 0: planes[nsTile].length;
+			int [] rslt = {-1, tp};
+			if (tp < num_planes) return rslt;
+			tp -= num_planes;
+			for (int d = 0; d < st_dirs8.length; d ++){
+				int nsTile1 = nsTile + st_dirs8[d];
+				num_planes = (planes[nsTile1] == null)? 0: planes[nsTile1].length;
+				if (tp < num_planes){
+					rslt[0] = d;
+					rslt[1] = tp;
+					return rslt;
+				}
+				tp -= num_planes; 
+			}
+			return null; // error - invalid input
+		}
+
+		
+}