debugging plane split

src/main/java/CLTPass3d.java

@@ -533,7 +533,8 @@ public class CLTPass3d{
 				double     max_disparity,
 				double     strength_floor,
 				double     strength_pow,
-				double     stBlurSigma)
+				double     stBlurSigma,
+				int        measSel)
 		{
 			this.superTiles = new SuperTiles(
 					this,
@@ -544,9 +545,18 @@ public class CLTPass3d{
 					max_disparity,
 					strength_floor,
 					strength_pow,
-					stBlurSigma);
+					stBlurSigma,
+					measSel);
 			return this.superTiles;
 		}
+		public double [] showDisparityHistogram(int measSel)
+		{
+			if (this.superTiles == null){
+				return null;
+			}
+			return this.superTiles.showDisparityHistogram(measSel);
+		}
+
 		public double [] showDisparityHistogram()
 		{
 			if (this.superTiles == null){

src/main/java/EyesisCorrectionParameters.java

@@ -2123,6 +2123,9 @@ public class EyesisCorrectionParameters {
   		public double     stUseDisp         = 0.15;  // Use background disparity from supertiles if tile strength is less 
   		public double     stStrengthScale   = 50.0;  // Multiply st strength if used instead of regular strength
   		
+  		public int        stMeasSel         = 1;     // Select measurements for supertiles : +1 - combo, +2 - quad +4 - hor +8 - vert
+  		
+  		
   		public double     outlayerStrength  = 0.3;   // Outlayer tiles weaker than this may be replaced from neighbors
   		public double     outlayerDiff      = 0.4;   // Replace weak outlayer tiles that do not have neighbors within this disparity difference
   		public double     outlayerDiffPos   = 1.0;   // Replace weak outlayer tiles that have higher disparity than weighted average
@@ -2157,16 +2160,28 @@ public class EyesisCorrectionParameters {
   		public boolean    plFillSquares        =   true; // Add diagonals to full squares
   		public boolean    plCutCorners         =   true; // Add ortho to 45-degree corners
 
-
-  		public double     plPull               =  .1;   // Relative weight of original (measured) plane when combing with neighbors
+  		public double     plPull               =  .3;   // Relative weight of original (measured) plane compared to average neighbor pull
+  		                                                // when combing with neighbors
   		public int        plIterations         =  10;   // Maximal number of smoothing iterations for each step
   		public boolean    plStopBad            =  true; // Do not update supertile if any of connected neighbors is not good (false: just skip that neighbor)
   		public int        plPrecision          =  6;    // Maximal step difference (1/power of 10)
   		
   		public double     plSplitPull          =  .5;   // Relative weight of center plane when splitting into pairs
+  		public double     plNormPow            =  .5;   // 0.0: 8 neighbors pull 8 times as 1, 1.0 - same as 1
   		public int        plSplitMinNeib       =  2;    // Minimal number of neighbors to split plane in pairs
   		public double     plSplitMinWeight     =  2.0;  // Minimal weight of split plains to show
-  		public double     plSplitMinQuality    =  1.1;  // Minimal split quality to show
+  		public double     plSplitMinQuality    =  1.1;  // Maximal normalized disparity difference from the plane to consider
+  		
+  		public boolean    plSplitApply         = true;  // Apply plane split to pairs
+  		public boolean    plNonExclusive       = true;  // Allow tiles to belong to both planes of the pair
+  		public boolean    plUseOtherPlanes     = false; // Allow other tiles from the same supertile
+  		public boolean    plAllowParallel      = true;  // Allow parallel shift of the specified planes before adding
+  		public double     plMaxDiff            = 0.3;   // Maximal normalized tile disparity difference from the plane to consider
+  		public double     plOtherDiff          = 1.4;   // Maximal difference of the added tile ratio to the average  disparity difference
+  		
+  		
+  		public boolean    plSplitXY            = true;  // Separate tiles for split planes by X, Y 
+  		public double     plSplitXYTolerance   = 0.2;   // Disparity tolerance when separating by X, Y 
   		
   		
   		public boolean    plFuse               =  true; // Fuse planes together (off for debug only)
@@ -2399,6 +2414,8 @@ public class EyesisCorrectionParameters {
 			properties.setProperty(prefix+"stMinBgFract",     this.stMinBgFract +"");
 			properties.setProperty(prefix+"stUseDisp",        this.stUseDisp +"");
 			properties.setProperty(prefix+"stStrengthScale",  this.stStrengthScale +"");
+  			properties.setProperty(prefix+"stMeasSel",        this.stMeasSel+"");
+			
 			properties.setProperty(prefix+"outlayerStrength", this.outlayerStrength +"");
 			properties.setProperty(prefix+"outlayerDiff",     this.outlayerDiff +"");
 			properties.setProperty(prefix+"outlayerDiffPos",  this.outlayerDiffPos +"");
@@ -2433,6 +2450,7 @@ public class EyesisCorrectionParameters {
 			properties.setProperty(prefix+"plCutCorners",     this.plCutCorners+"");
 
 			properties.setProperty(prefix+"plPull",           this.plPull +"");
+			properties.setProperty(prefix+"plNormPow",        this.plNormPow +"");
   			properties.setProperty(prefix+"plIterations",     this.plIterations+"");
 			properties.setProperty(prefix+"plStopBad",        this.plStopBad+"");
   			properties.setProperty(prefix+"plPrecision",      this.plPrecision+"");
@@ -2441,6 +2459,14 @@ public class EyesisCorrectionParameters {
   			properties.setProperty(prefix+"plSplitMinNeib",   this.plSplitMinNeib+"");
 			properties.setProperty(prefix+"plSplitMinWeight", this.plSplitMinWeight +"");
 			properties.setProperty(prefix+"plSplitMinQuality",this.plSplitMinQuality +"");
+			properties.setProperty(prefix+"plSplitApply",     this.plSplitApply+"");
+			properties.setProperty(prefix+"plNonExclusive",   this.plNonExclusive+"");
+			properties.setProperty(prefix+"plUseOtherPlanes", this.plUseOtherPlanes+"");
+			properties.setProperty(prefix+"plAllowParallel",  this.plAllowParallel+"");
+			properties.setProperty(prefix+"plMaxDiff",        this.plMaxDiff +"");
+			properties.setProperty(prefix+"plOtherDiff",      this.plOtherDiff +"");
+			properties.setProperty(prefix+"plSplitXY",        this.plSplitXY+"");
+			properties.setProperty(prefix+"plSplitXYTolerance",this.plSplitXYTolerance +"");
 
   			properties.setProperty(prefix+"plFuse",           this.plFuse+"");
 			properties.setProperty(prefix+"plKeepOrphans",    this.plKeepOrphans+"");
@@ -2661,6 +2687,8 @@ public class EyesisCorrectionParameters {
   			if (properties.getProperty(prefix+"stMinBgFract")!=null)      this.stMinBgFract=Double.parseDouble(properties.getProperty(prefix+"stMinBgFract"));
   			if (properties.getProperty(prefix+"stUseDisp")!=null)         this.stUseDisp=Double.parseDouble(properties.getProperty(prefix+"stUseDisp"));
   			if (properties.getProperty(prefix+"stStrengthScale")!=null)   this.stStrengthScale=Double.parseDouble(properties.getProperty(prefix+"stStrengthScale"));
+  			if (properties.getProperty(prefix+"stMeasSel")!=null)         this.stMeasSel=Integer.parseInt(properties.getProperty(prefix+"stMeasSel"));
+
   			if (properties.getProperty(prefix+"outlayerStrength")!=null)  this.outlayerStrength=Double.parseDouble(properties.getProperty(prefix+"outlayerStrength"));
   			if (properties.getProperty(prefix+"outlayerDiff")!=null)      this.outlayerDiff=Double.parseDouble(properties.getProperty(prefix+"outlayerDiff"));
   			if (properties.getProperty(prefix+"outlayerDiffPos")!=null)   this.outlayerDiffPos=Double.parseDouble(properties.getProperty(prefix+"outlayerDiffPos"));
@@ -2695,6 +2723,7 @@ public class EyesisCorrectionParameters {
   			if (properties.getProperty(prefix+"plCutCorners")!=null)      this.plCutCorners=Boolean.parseBoolean(properties.getProperty(prefix+"plCutCorners"));
 
   			if (properties.getProperty(prefix+"plPull")!=null)            this.plPull=Double.parseDouble(properties.getProperty(prefix+"plPull"));
+  			if (properties.getProperty(prefix+"plNormPow")!=null)         this.plNormPow=Double.parseDouble(properties.getProperty(prefix+"plNormPow"));
   			if (properties.getProperty(prefix+"plIterations")!=null)      this.plIterations=Integer.parseInt(properties.getProperty(prefix+"plIterations"));
   			if (properties.getProperty(prefix+"plStopBad")!=null)         this.plStopBad=Boolean.parseBoolean(properties.getProperty(prefix+"plStopBad"));
   			if (properties.getProperty(prefix+"plPrecision")!=null)       this.plPrecision=Integer.parseInt(properties.getProperty(prefix+"plPrecision"));
@@ -2703,6 +2732,14 @@ public class EyesisCorrectionParameters {
   			if (properties.getProperty(prefix+"plSplitMinNeib")!=null)    this.plSplitMinNeib=Integer.parseInt(properties.getProperty(prefix+"plSplitMinNeib"));
   			if (properties.getProperty(prefix+"plSplitMinWeight")!=null)  this.plSplitMinWeight=Double.parseDouble(properties.getProperty(prefix+"plSplitMinWeight"));
   			if (properties.getProperty(prefix+"plSplitMinQuality")!=null) this.plSplitMinQuality=Double.parseDouble(properties.getProperty(prefix+"plSplitMinQuality"));
+  			if (properties.getProperty(prefix+"plSplitApply")!=null)      this.plSplitApply=Boolean.parseBoolean(properties.getProperty(prefix+"plSplitApply"));
+  			if (properties.getProperty(prefix+"plNonExclusive")!=null)    this.plNonExclusive=Boolean.parseBoolean(properties.getProperty(prefix+"plNonExclusive"));
+  			if (properties.getProperty(prefix+"plUseOtherPlanes")!=null)  this.plUseOtherPlanes=Boolean.parseBoolean(properties.getProperty(prefix+"plUseOtherPlanes"));
+  			if (properties.getProperty(prefix+"plAllowParallel")!=null)   this.plAllowParallel=Boolean.parseBoolean(properties.getProperty(prefix+"plAllowParallel"));
+  			if (properties.getProperty(prefix+"plMaxDiff")!=null)         this.plMaxDiff=Double.parseDouble(properties.getProperty(prefix+"plMaxDiff"));
+  			if (properties.getProperty(prefix+"plOtherDiff")!=null)       this.plOtherDiff=Double.parseDouble(properties.getProperty(prefix+"plOtherDiff"));
+  			if (properties.getProperty(prefix+"plSplitXY")!=null)         this.plSplitXY=Boolean.parseBoolean(properties.getProperty(prefix+"plSplitXY"));
+  			if (properties.getProperty(prefix+"plSplitXYTolerance")!=null)this.plSplitXYTolerance=Double.parseDouble(properties.getProperty(prefix+"plSplitXYTolerance"));
 
   			if (properties.getProperty(prefix+"plFuse")!=null)            this.plFuse=Boolean.parseBoolean(properties.getProperty(prefix+"plFuse"));
   			if (properties.getProperty(prefix+"plKeepOrphans")!=null)     this.plKeepOrphans=Boolean.parseBoolean(properties.getProperty(prefix+"plKeepOrphans"));
@@ -2944,6 +2981,8 @@ public class EyesisCorrectionParameters {
   			gd.addNumericField("Minimal fraction of the disparity histogram to use as background",             this.stMinBgFract,  6);
   			gd.addNumericField("Use background disparity from supertiles if tile strength is less",            this.stUseDisp,  6);
   			gd.addNumericField("Multiply st strength if used instead of regular strength ",                    this.stStrengthScale,  6);
+  			gd.addNumericField("Select measurements for supertiles : +1 - combo, +2 - quad +4 - hor +8 - vert",this.stMeasSel,  0);
+
   			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);
@@ -2978,15 +3017,24 @@ public class EyesisCorrectionParameters {
   			gd.addCheckbox    ("Add diagonals to full squares",                                                this.plFillSquares);
   			gd.addCheckbox    ("Add ortho to 45-degree corners",                                               this.plCutCorners);
 
-  			gd.addNumericField("Relative weight of original (measured) plane when combing with neighbors",     this.plPull,  6);
+  			gd.addNumericField("Relative (to average neighbor) weight of the measured plane when combing with neighbors",     this.plPull,  6);
+  			gd.addNumericField("0.0: 8 neighbors pull 8 times as 1, 1.0 - same as 1",                          this.plNormPow,  6);
   			gd.addNumericField("Maximal number of smoothing iterations for each step",                         this.plIterations,  0);
   			gd.addCheckbox    ("Do not update supertile if any of connected is not good (false: just skip that neighbor)", this.plStopBad);
   			gd.addNumericField("Maximal step difference (1/power of 10)",                                      this.plPrecision,  0);
 
   			gd.addNumericField("Relative weight of center plane when splitting into pairs",                    this.plSplitPull,  6);
   			gd.addNumericField("Minimal number of neighbors to split plane in pairs",                          this.plSplitMinNeib,  0);
-  			gd.addNumericField(" Minimal weight of split plains to show",                                      this.plSplitMinWeight,  6);
+  			gd.addNumericField("Minimal weight of split plains to show",                                       this.plSplitMinWeight,  6);
   			gd.addNumericField("Minimal split quality to show",                                                this.plSplitMinQuality,  6);
+  			gd.addCheckbox    ("Apply plane split to pairs",                                                   this.plSplitApply);
+  			gd.addCheckbox    ("Allow tiles to belong to both planes of the pair",                             this.plNonExclusive);
+  			gd.addCheckbox    ("Allow other tiles from the same supertile",                                    this.plUseOtherPlanes);
+  			gd.addCheckbox    ("Allow parallel shift of the specified planes before adding",                   this.plAllowParallel);
+  			gd.addNumericField("Maximal normalized tile disparity difference from the plane to consider",      this.plMaxDiff,  6);
+  			gd.addNumericField("Maximal difference of the added tile ratio to the average  disparity difference",this.plOtherDiff,  6);
+  			gd.addCheckbox    ("Separate tiles for split planes by X, Y",                                      this.plSplitXY);
+  			gd.addNumericField("Disparity tolerance when separating by X, Y",                                  this.plSplitXYTolerance,  6);
   			
   			gd.addCheckbox    ("Fuse planes together (off for debug only)",                                    this.plFuse);
   			gd.addCheckbox    ("Keep unconnected supertiles",                                                  this.plKeepOrphans);
@@ -3216,6 +3264,7 @@ public class EyesisCorrectionParameters {
   			this.stMinBgFract=          gd.getNextNumber();
   			this.stUseDisp=             gd.getNextNumber();
   			this.stStrengthScale=       gd.getNextNumber();
+  			this.stMeasSel=       (int) gd.getNextNumber();
   			this.outlayerStrength=      gd.getNextNumber();
   			this.outlayerDiff=          gd.getNextNumber();
   			this.outlayerDiffPos=       gd.getNextNumber();
@@ -3251,6 +3300,7 @@ public class EyesisCorrectionParameters {
   			this.plCutCorners=          gd.getNextBoolean();
 
   			this.plPull=                gd.getNextNumber();
+  			this.plNormPow=             gd.getNextNumber();
   			this.plIterations=    (int) gd.getNextNumber();
   			this.plStopBad=             gd.getNextBoolean();
   			this.plPrecision=     (int) gd.getNextNumber();
@@ -3259,6 +3309,14 @@ public class EyesisCorrectionParameters {
   			this.plSplitMinNeib=  (int) gd.getNextNumber();
   			this.plSplitMinWeight=      gd.getNextNumber();
   			this.plSplitMinQuality=     gd.getNextNumber();
+  			this.plSplitApply=          gd.getNextBoolean();
+  			this.plNonExclusive=        gd.getNextBoolean();
+  			this.plUseOtherPlanes=      gd.getNextBoolean();
+  			this.plAllowParallel=       gd.getNextBoolean();
+  			this.plMaxDiff=             gd.getNextNumber();
+  			this.plOtherDiff=           gd.getNextNumber();
+  			this.plSplitXY=             gd.getNextBoolean();
+  			this.plSplitXYTolerance=    gd.getNextNumber();
 
   			this.plFuse=                gd.getNextBoolean();
   			this.plKeepOrphans=         gd.getNextBoolean();

src/main/java/MeasuredLayers.java

@@ -224,7 +224,7 @@ public class MeasuredLayers {
 	 * @return number of layers (some may be uninitialized
 	 */
 	
-	public int getnumLayers()
+	public int getNumLayers()
 	{
 		if (layers == null){
 			return 0;

src/main/java/SuperTiles.java

@@ -25,8 +25,9 @@ import java.util.ArrayList;
 import java.util.Collections;
 import java.util.Comparator;
 import java.util.concurrent.atomic.AtomicInteger;
-import java.awt.Point;
-import java.util.ArrayList;
+
+//import java.awt.Point;
+//import java.util.ArrayList;
 
 //import TilePlanes.PlaneData;
 
@@ -52,7 +53,7 @@ public class SuperTiles{
 	double [][][] maxMinMax = null;
 	double []   bgDisparity = null;
 	double []   bgStrength = null;
-	
+	int         measSel = 1; // bitmask of the selected measurements for supertiles : +1 - combo, +2 - quad +4 - hor +8 - vert
 	MeasuredLayers measuredLayers = null;
 	
 	CLTPass3d cltPass3d;
@@ -63,6 +64,15 @@ public class SuperTiles{
 	int [][]    shell_map = null;   // per supertile, per disparity plane - shell index + 1 (0 - none)
 	double [][] surfaces;  // per shell, per tile (linescan order) disparity value or NaN in missing supertiles
 
+	/**
+	 * currently lowest plane for each includes all tiles, so do not use it. May change in the future
+	 * @param np total number of planes in a supertile
+	 * @return lowest plane to iterate through
+	 */
+	
+	static int LOWEST_PLANE(int np){ // 
+		return (np >1) ? 1 : 0;
+	}
 
 	public SuperTiles(
 			CLTPass3d cltPass3d,
@@ -73,7 +83,8 @@ public class SuperTiles{
 			double                  max_disparity,
 			double                  strength_floor,
 			double                  strength_pow,
-			double                  stBlurSigma)
+			double                  stBlurSigma,
+			int                     measSel)
 	{
 		this.cltPass3d =           cltPass3d;
 		this.tileProcessor =       cltPass3d.getTileProcessor();
@@ -85,6 +96,7 @@ public class SuperTiles{
 		this.strength_floor = strength_floor;
 		this.strength_pow =   strength_pow;
 		this.stBlurSigma =    stBlurSigma;
+		this.measSel =        measSel;
 		//				this.step_threshold_near = this.step_threshold_far * this.step_far / step_near;
 		this.step_threshold_near = this.step_threshold_far * step_near / this.step_far ;
 		this.bin_far =             this.step_threshold_far / this.step_far; 
@@ -96,18 +108,6 @@ public class SuperTiles{
 		for (bin = 0; bin < numBins; bin ++){
 			bin_centers[bin] = binToDisparity(bin);  
 		}
-
-		getDisparityHistograms(null); // calculate and blur supertiles (for all, not just selected?)
-		if (tileProcessor.globalDebugLevel > -1){
-			System.out.println("SuperTiles(): min_disparity = "+min_disparity+", max_disparity="+max_disparity);
-			System.out.println("SuperTiles(): step_far = "+step_far+", step_near="+step_near);
-			System.out.println("SuperTiles(): numBins = "+numBins+", bin_far="+bin_far+", bin_near = "+bin_near);
-			System.out.println("SuperTiles(): step_threshold_far = "+step_threshold_far+", step_threshold_near="+step_threshold_near);
-			for (int i = 0; i<numBins; i++){
-				System.out.println(i+": "+bin_centers[i]+", "+disparityToBin(bin_centers[i]));
-			}
-			//
-		}
 		initFuseCoeff(0.5, false); // true);
 		
 // Set up MeasuredLayers		
@@ -138,9 +138,21 @@ public class SuperTiles{
 				cltPass3d.getVertStrength(),      // double [] strength,
 				cltPass3d.getSelected());         // boolean [] selection) // may be null
 
+		getDisparityHistograms(measSel); // calculate and blur supertiles (for all, not just selected?)
+		if (tileProcessor.globalDebugLevel > -1){
+			System.out.println("SuperTiles(): min_disparity = "+min_disparity+", max_disparity="+max_disparity);
+			System.out.println("SuperTiles(): step_far = "+step_far+", step_near="+step_near);
+			System.out.println("SuperTiles(): numBins = "+numBins+", bin_far="+bin_far+", bin_near = "+bin_near);
+			System.out.println("SuperTiles(): step_threshold_far = "+step_threshold_far+", step_threshold_near="+step_threshold_near);
+			for (int i = 0; i<numBins; i++){
+				System.out.println(i+": "+bin_centers[i]+", "+disparityToBin(bin_centers[i]));
+			}
+			//
+		}
+		
 		String [] titles = {"d0","s0","d1","s1","d2","s2","d3","s3"};
 		double [][] dbg_img = new double [titles.length][];
-		for (int i = 0; i < measuredLayers.getnumLayers(); i++){
+		for (int i = 0; i < measuredLayers.getNumLayers(); i++){
 			dbg_img[2 * i] =     measuredLayers.getDisparity(i);
 			dbg_img[2 * i + 1] = measuredLayers.getStrength(i);
 		}
@@ -315,6 +327,7 @@ public class SuperTiles{
 		return lapWeight;
 	}
 	// updates disparityHistograms
+/*	
 	public double [][] getDisparityHistograms()
 	{
 		return getDisparityHistograms(cltPass3d.selected, tileProcessor.globalDebugLevel);
@@ -324,13 +337,13 @@ public class SuperTiles{
 	{
 		return getDisparityHistograms(selected, tileProcessor.globalDebugLevel);
 	}
-
-	public double [][] getDisparityHistograms(
+*/
+	public double [][] getDisparityHistogramsOld( // OLD!
 			final boolean [] selected, // or null
 			final int        debugLevel)
 	{
 		if (this.disparityHistograms != null) return this.disparityHistograms;
-		final double step_disparity = step_near; // TODO: implement
+///		final double step_disparity = step_near; // TODO: implement
 		final int tilesX =        tileProcessor.getTilesX();
 		final int tilesY =        tileProcessor.getTilesY();
 		final int superTileSize = tileProcessor.superTileSize;
@@ -403,6 +416,80 @@ public class SuperTiles{
 		return this.disparityHistograms; // dispHist;
 	}
 
+	public double [][] getDisparityHistograms(
+			final int        measSel) // bitmask of the selected measurements for supertiles : +1 - combo, +2 - quad +4 - hor +8 - vert
+//			final int        debugLevel)
+	{
+		if (this.disparityHistograms != null) return this.disparityHistograms;
+//		final double step_disparity = step_near; // TODO: implement
+		final int tilesX =        tileProcessor.getTilesX();
+		final int tilesY =        tileProcessor.getTilesY();
+		final int superTileSize = tileProcessor.superTileSize;
+		final int stilesX = (tilesX + superTileSize -1)/superTileSize;  
+		final int stilesY = (tilesY + superTileSize -1)/superTileSize;
+		final int nStiles = stilesX * stilesY; 
+		final double [][] dispHist =     new double [nStiles][]; // now it will be sparse 
+		final double []   strengthHist = new double [nStiles];
+		final Thread[] threads = ImageDtt.newThreadArray(tileProcessor.threadsMax);
+		final AtomicInteger ai = new AtomicInteger(0);
+		for (int ithread = 0; ithread < threads.length; ithread++) {
+			threads[ithread] = new Thread() {
+				public void run() {
+					for (int nsTile = ai.getAndIncrement(); nsTile < nStiles; nsTile = ai.getAndIncrement()) {
+						int stileY = nsTile / stilesX;  
+						int stileX = nsTile % stilesX;  
+						double sw = 0.0; // sum weights
+						double [] hist = new double [numBins];
+						for (int nl = 0; nl < measuredLayers.getNumLayers(); nl ++) {
+							if ((measSel & (1 << nl)) != 0) {
+								double [][] disp_strength =  measuredLayers.getDisparityStrength(
+										nl,             // 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 (disp_strength != null) {
+									for (int indx = 0; indx < disp_strength[1].length; indx++) {
+										double w = disp_strength[1][indx]; 
+										if ( w > 0.0){
+											double d = disp_strength[0][indx]; 
+											int bin = disparityToBin(d);
+											if ((bin >= 0) && (bin < numBins)){ // maybe collect below min and above max somewhere?
+												hist[bin] += w; // +1]
+												sw +=w;
+											}
+										}
+									}
+								}
+							}
+						}
+						strengthHist[nsTile] = sw / superTileSize / superTileSize; // average strength per tile in the super-tile
+						if (sw > 0){
+							for (int i = 0; i<numBins; i++){
+								hist[i] /= sw; 
+							}
+							dispHist[nsTile] = hist;
+						} else {
+							dispHist[nsTile] = null;
+						}
+					}
+				}
+			};
+		}		      
+		ImageDtt.startAndJoin(threads);
+		this.disparityHistograms = dispHist;
+		this.stStrength =          strengthHist;
+		if (this.stBlurSigma > 0.0) {
+			blurDisparityHistogram(0); // debugLevel);
+		}
+		return this.disparityHistograms; // dispHist;
+	}
+
+	
+	
+	
 	
 
 	public void blurDisparityHistogram( // in-place
@@ -433,7 +520,7 @@ public class SuperTiles{
 
 	public double [][][] getMaxMinMax(){
 		// first find all integer maximums, and if the top is flat - use the middle. If not flat - use 2-nd degree polynomial
-		if (disparityHistograms == null) getDisparityHistograms();
+		if (disparityHistograms == null) getDisparityHistograms(this.measSel);
 		final int globalDebugLevel = tileProcessor.globalDebugLevel;
 		maxMinMax = new double [disparityHistograms.length][][];
 		final Thread[] threads = ImageDtt.newThreadArray(tileProcessor.threadsMax);
@@ -560,11 +647,16 @@ public class SuperTiles{
 		int sTilesX = (tileProcessor.getTilesX() + superTileSize -1)/superTileSize;  
 		return sTilesX * (numBins + 1) + 1;
 	}
-
 	public double [] showDisparityHistogram()
+	{
+		return showDisparityHistogram(this.measSel);
+	}
+	
+	public double [] showDisparityHistogram(
+			final int        measSel) // bitmask of the selected measurements for supertiles : +1 - combo, +2 - quad +4 - hor +8 - vert
 	{
 		if (disparityHistograms == null){
-			getDisparityHistograms(); // calculate and blur with the current settings, specified at instantiation
+			getDisparityHistograms(measSel); // calculate and blur with the current settings, specified at instantiation
 		}
 		return showDisparityHistogram(disparityHistograms);
 	}
@@ -1528,7 +1620,8 @@ public class SuperTiles{
 	public void processPlanes3(
 			final boolean [] selected, // or null
 			final double     min_disp,
-			final boolean    invert_disp, // use 1/disparity
+//			final boolean    invert_disp, // use 1/disparity
+			final int        stMeasSel, //            = 1;      // Select measurements for supertiles : +1 - combo, +2 - quad +4 - hor +8 - vert
 			final double     plDispNorm,
 			final int        plMinPoints, //          =     5;  // Minimal number of points for plane detection
 			final double     plTargetEigen, //        =   0.1;  // Remove outliers until main axis eigenvalue (possibly scaled by plDispNorm) gets below
@@ -1554,7 +1647,17 @@ public class SuperTiles{
 		final AtomicInteger ai = new AtomicInteger(0);
 		this.planes = new TilePlanes.PlaneData[nStiles][];
 		final int debug_stile = (debugLevel > -1)? (dbg_Y * stilesX + dbg_X):-1;
-		final boolean [][] dflt_select = {{}, null, null, null, null}; // use layer 0 (combo) only
+		
+//		final boolean [][] dflt_select = {{}, null, null, null, null}; // use layer 0 (combo) only
+		final boolean [][] dflt_select = new boolean [measuredLayers.getNumLayers()][];
+		for (int i = 0; i < dflt_select.length; i++){
+			if ((stMeasSel & (1 << i)) !=0){
+				dflt_select[i] = new boolean[0];
+			} else {
+				dflt_select[i] = null;
+			}
+		}
+		
 		// TODO: Remove when promoting PlaneData
 		final TilePlanes tpl = new TilePlanes(tileSize,superTileSize, geometryCorrection);
 		
@@ -1569,7 +1672,7 @@ public class SuperTiles{
 		if (debugLevel > -1) {
 			String [] titles = {"d0","s0","d1","s1","d2","s2","d3","s3","s","d"};
 			double [][] dbg_img = new double [titles.length][];
-			for (int i = 0; i < measuredLayers.getnumLayers(); i++){
+			for (int i = 0; i < measuredLayers.getNumLayers(); i++){
 				dbg_img[2 * i] =     measuredLayers.getDisparity(i);
 				dbg_img[2 * i + 1] = measuredLayers.getStrength(i);
 			}
@@ -1578,6 +1681,9 @@ public class SuperTiles{
 			showDoubleFloatArrays sdfa_instance = new showDoubleFloatArrays();
 			sdfa_instance.showArrays(dbg_img,  tileProcessor.getTilesX(), tileProcessor.getTilesY(), true, "measuredLayers",titles);
 		}
+//		if (maxMinMax == null) 
+			getMaxMinMax();
+
 		for (int ithread = 0; ithread < threads.length; ithread++) {
 			threads[ithread] = new Thread() {
 				public void run() {
@@ -1660,9 +1766,14 @@ public class SuperTiles{
 
 							// now try for each of the disparity-separated clusters (only for multi-peak histograms)
 
+
 							double [][] mm = maxMinMax[nsTile];
+							if (mm == null){
+								double [][][] dbg_min_max = maxMinMax;
+								System.out.println("maxMinMax["+nsTile+"] == null");
+							}
 
-							if (mm.length > 1) { // multiple maximums - separate into multiple selections
+							if ((mm!= null) &&  (mm.length > 1)) { // multiple maximums - separate into multiple selections // null pointer
 								for (int m = 0; m < (mm.length +1)/2; m++){
 									double [] far_near = new double [2];
 									if (m == 0) {
@@ -1954,6 +2065,7 @@ public class SuperTiles{
 							for (int np0 = 0; np0 < planes[nsTile0].length; np0++){ // nu
 								//										planes[nsTile0][np0].initNeibBest(); // 
 								TilePlanes.PlaneData this_plane = planes[nsTile0][np0];
+								if (this_plane != null) {
 									this_plane.initMergedValue();
 									for (int dir = 0; dir < 4; dir++){ // just half directions - relations are symmetrical
 										int stx = stx0 + dirsYX[dir][1];
@@ -1969,6 +2081,7 @@ public class SuperTiles{
 
 													this_plane.initMergedValue(dir,other_planes.length); // filled with NaN
 													for (int np = 0; np < other_planes.length; np ++){
+														if (other_planes[np] != null) {
 															TilePlanes.PlaneData other_plane = this_plane.getPlaneToThis(
 																	other_planes[np],
 																	dl-1); // debugLevel);
@@ -1982,7 +2095,7 @@ public class SuperTiles{
 																		dl-1); // int       debugLevel)
 
 																if (merged_pd !=null) { // now always, but may add later
-///															merged_pd.scaleWeight(0.5);
+																	///															merged_pd.scaleWeight(0.5);
 																	this_plane.setNeibMatch(dir, np, merged_pd.getValue()); // smallest eigenValue
 																}
 															}
@@ -1995,6 +2108,8 @@ public class SuperTiles{
 								}
 							}
 						}
+					}
+				}
 			};
 		}		      
 		ImageDtt.startAndJoin(threads);
@@ -2014,6 +2129,7 @@ public class SuperTiles{
 						if ( planes[nsTile0] != null) {
 							for (int np0 = 0; np0 < planes[nsTile0].length; np0++){ // nu
 								TilePlanes.PlaneData this_plane = planes[nsTile0][np0];
+								if (this_plane != null) {
 									for (int dir = 4; dir < 8; dir++){ // other half - copy from opposite
 										int stx = stx0 + dirsYX[dir][1];
 										int sty = sty0 + dirsYX[dir][0];
@@ -2038,6 +2154,7 @@ public class SuperTiles{
 							}
 						}
 					}
+				}
 			};
 		}		      
 		ImageDtt.startAndJoin(threads);
@@ -2081,7 +2198,7 @@ public class SuperTiles{
 										double [] neib_worse = this_plane.getMergedValue(dir);
 										if (neib_worse != null) {
 											int best_index = -1;  
-											int np_min = (neib_worse.length > 1) ? 1:0; // Modify if overall plane will be removed
+											int np_min = LOWEST_PLANE(neib_worse.length);
 											for (int np = np_min; np < neib_worse.length; np++){
 												if (	((worst_worsening == 0.0) || (neib_worse[np] < worst_worsening)) &&
 														!Double.isNaN(neib_worse[np]) &&
@@ -2372,11 +2489,13 @@ public class SuperTiles{
 						if ( planes[nsTile0] != null) {
 							for (int np0 = 0; np0 < planes[nsTile0].length; np0++){ // nu
 								TilePlanes.PlaneData this_plane = planes[nsTile0][np0];
+								if (this_plane != null) {
 									this_plane.initNeibBest();
 								}
 							}
 						}
 					}
+				}
 			};
 		}		      
 		ImageDtt.startAndJoin(threads);
@@ -2402,13 +2521,14 @@ public class SuperTiles{
 
 								int num_other_planes = 0;
 								for (int np0 = np0_min; np0 < planes[nsTile0].length; np0++){
-									if (planes[nsTile0][np0].getMergedValue(dir) != null){
+									if ((planes[nsTile0][np0] != null) && (planes[nsTile0][np0].getMergedValue(dir) != null)){
 										int l = planes[nsTile0][np0].getMergedValue(dir).length;
 										if (l > num_other_planes)num_other_planes = l;
 									}
 								}
 								if (num_other_planes > 0){ // will eliminate bad margins
-									int np_min = (num_other_planes > 1) ? 1:0; // Modify if overall plane will be removed
+									int np_min = LOWEST_PLANE(num_other_planes);
+									
 									boolean [] this_matched = new boolean [planes[nsTile0].length];
 									boolean [] other_matched = new boolean [num_other_planes];
 									int num_pairs = this_matched.length - np0_min;
@@ -2417,7 +2537,7 @@ public class SuperTiles{
 									for (int pair = 0; pair < num_pairs; pair ++){
 										int [] best_pair = {-1,-1};
 										double best_rqual = Double.NaN;
-										for (int np0 = np0_min; np0 < this_matched.length; np0++){
+										for (int np0 = np0_min; np0 < this_matched.length; np0++) if (planes[nsTile0][np0] != null){
 											double [] merge_ev = planes[nsTile0][np0].getMergedValue(dir);
 											if (!this_matched[np0] &&
 													(merge_ev != null) &&
@@ -2429,6 +2549,7 @@ public class SuperTiles{
 													(planes[nsTile0][np0].getWeight() > minWeight)) {
 												for (int np = np_min; np < merge_ev.length; np++){
 													if (!other_matched[np] &&
+															(planes[nsTile][np] != null) &&
 															!Double.isNaN(merge_ev[np]) &&
 															((maxEigen == 0.0) ||
 																	(planes[nsTile][np].getValue() < corrMaxEigen(
@@ -2476,12 +2597,12 @@ public class SuperTiles{
 										planes[nsTile][best_pair[1]].setNeibBest(dir + 4,best_pair[0]);
 									}
 									// disable remaining neighbors
-									for (int np = 0; np < this_matched.length; np++){
+									for (int np = 0; np < this_matched.length; np++) if (planes[nsTile0][np] != null){
 										if (!this_matched[np]){
 											planes[nsTile0][np].setNeibBest(dir,-1);
 										}
 									}
-									for (int np = 0; np < other_matched.length; np++){
+									for (int np = 0; np < other_matched.length; np++) if (planes[nsTile][np] != null){
 										if (!other_matched[np]){
 											planes[nsTile][np].setNeibBest(dir + 4,-1);
 										}
@@ -2537,8 +2658,8 @@ public class SuperTiles{
 			}
 			if ((planes[i] != null) && (planes[i].length > num_planes)){
 				int num_sat = 0;
-				for (int j = 0; j < planes[i].length; j++){
-					if (planes[i][j].getWeight() >= minWeight){
+				for (int j = 0; j < planes[i].length; j++) if (planes[i][j] != null){
+//					if (planes[i][j].getWeight() >= minWeight){ /// Later does not compare too minWeight
 						double eigVal = planes[i][j].getValue(); // ************** generated does not have value ????????????
 						double disp = planes[i][j].getPlaneDisparity(false)[centerIndex];
 						/*
@@ -2551,7 +2672,7 @@ public class SuperTiles{
 								maxEigen,
 								dispNorm,
 								disp)) num_sat ++;
-					}
+//					}
 				}
 				if (num_sat > num_planes) num_planes = num_sat;
 			}
@@ -2571,7 +2692,7 @@ public class SuperTiles{
 					System.out.println("getShowPlanes():"+((planes[nsTile] != null)? planes[nsTile].length : "null"));
 				}
 				if (planes[nsTile] != null) {
-					for (int np = 0; np < planes[nsTile].length; np++){
+					for (int np = 0; np < planes[nsTile].length; np++) if (planes[nsTile][np] != null){
 						if (planes[nsTile][np].getWeight() >= minWeight){
 							//&& (planes[nsTile][np].getValue() < maxEigen)){
 							double eigVal = planes[nsTile][np].getValue();
@@ -2822,6 +2943,8 @@ public class SuperTiles{
 			final double      maxValue, // do not combine with too bad planes with primary eigenvalue above this value ( 0 any OK)
 			final int         num_passes,
 			final boolean     stopBadNeib, // do not update supertile if any of connected neighbors is not good (false: just skip that neighbor)
+			final double      normPow,
+			
 			final double      maxDiff, // maximal change in any of the disparity values
 			final boolean     preferDisparity, // Always start with disparity-most axis (false - lowest eigenvalue)
 			final int         debugLevel,
@@ -2836,6 +2959,7 @@ public class SuperTiles{
 					meas_pull,       //  relative pull of the original (measured) plane with respect to the average of the neighbors
 					maxValue,        // final double maxValue, // do not combine with too bad planes
 					stopBadNeib,     // do not update supertile if any of connected neighbors is not good (false: just skip that neighbor)
+					normPow,         // 0.0 - 8 neighbors pull 8 times more than 1. 1.0 - same as one
 					this.planes,     // final TilePlanes.PlaneData[][] measured_planes,
 					this.planes_mod, // final TilePlanes.PlaneData[][] mod_planes,
 					true,            // final boolean calc_diff,
@@ -2861,6 +2985,7 @@ public class SuperTiles{
 			final double meas_pull,//  relative pull of the original (measured) plane with respect to the average of the neighbors
 			final double maxValue, // do not combine with too bad planes
 			final boolean stopBadNeib, // do not update supertile if any of connected neighbors is not good (false: just skip that neighbor)
+			final double normPow,
 			final TilePlanes.PlaneData[][] measured_planes,
 			final TilePlanes.PlaneData[][] mod_planes,
 			final boolean calc_diff,
@@ -2874,7 +2999,7 @@ public class SuperTiles{
 		final int tilesY =        tileProcessor.getTilesY();
 		final int superTileSize = tileProcessor.getSuperTileSize();
 		final int stilesX = (tilesX + superTileSize -1)/superTileSize;  
-		final int stilesY = (tilesY + superTileSize -1)/superTileSize;
+//		final int stilesY = (tilesY + superTileSize -1)/superTileSize;
 		final int debug_stile = dbg_Y * stilesX + dbg_X;
 		final TilePlanes.PlaneData[][] new_planes = copyPlanes(mod_planes);
 		final Thread[] threads = ImageDtt.newThreadArray(tileProcessor.threadsMax);
@@ -2946,6 +3071,14 @@ public class SuperTiles{
 									}
 								}
 								if (this_new_plane != null) {
+									// average weight over participating directions, so the relative pull
+									// does not depend on number of neighbors
+									if ((num_merged > 0.0) && (this_new_plane != null) && (normPow > 1.0)){
+										double scale = Math.pow(num_merged, normPow);
+										this_new_plane.scaleWeight(1.0/scale);
+										num_merged /=scale;
+									}
+									
 									if (    (meas_pull > 0.0) &&
 											(measured_planes != null) &&
 											(measured_planes[nsTile0] != null) &&
@@ -2962,7 +3095,7 @@ public class SuperTiles{
 													preferDisparity,
 													dl - 1);                           // int       debugLevel)
 											if (this_new_plane != null){
-												num_merged += meas_pull;	
+												num_merged += meas_pull;	// num_merged was 1.0 and weight is averaged over all neighbors 
 											}
 										} else {
 											this_new_plane = measured_planes[nsTile0][np0].clone();
@@ -3155,6 +3288,341 @@ public class SuperTiles{
 		return split_lines;
 	}
 	
+	/**
+	 * Re-generate "planes" (ellipsoids) from the measured data according to provided
+	 * plane pairs for some supertiles.
+	 * Current version destroys connections between neighbor planes, they need to be
+	 * re-calculated
+	 * TODO: consider amplifying difference (angle) between planes before/during/after
+	 * running this method. Like 3d "sharpening" to compensate for being washed out.
+	 * @param planes per-supertile (in linescan order), per-plane data to be replaced
+	 * @param brokenPd per-supertile , per-plane 3 PD instances (two first make the pair).
+	 *        generated by breakPlanesToPairs method. This data will not be modified.
+	 * @param max_diff maximal normalized disparity difference from the plane to consider
+	 * @param other_diff maximal difference of the added tile ratio to the average
+	 *  disparity difference of the exclusively selected tiles
+	 * @param non_exclusive allow tiles to belong to both planes of the pair
+	 * @param use_other_planes allow other tiles from the same supertile
+	 * @param measSel (with use_other_planes) select measurements for supertiles:
+	 *  +1 - combo, +2 - quad +4 - hor +8 - vert
+	 * @param allow_parallel allow parallel shift of the specified planes before adding
+	 *  more tiles (from the same pair and/or other tiles of the supertile (see
+	 *  non_exclusive, use_other_planes
+	 * @param splitXY generate XY separation masks for split planes, so tiles will not
+	 *        be allowed to snap to the other plane in pair
+	 * @param splitXYTolerance XY separation disparity tolerance (there will be some
+	 *  overlapping if > 0.0
+	 * @param disp_far disparity low limit 
+	 * @param disp_near disparity high limit (Double.NaN - no limit)
+	 * @param dispNorm normalize disparity difference (keep relative) for disparities above
+	 * @param min_weight minimal total weight of the plane
+	 * @param min_tiles minimal number of tiles in a plane
+	 * @param targetEigen target primary eigenvalue when removing outliers
+	 * @param fractOutliers maximal fraction of outliers from all tiles to remove
+	 * @param maxOutliers maximal number of outliers to remove
+	 * @param strength_floor subtract from correlation strength (and limit by 0) to use
+	 *  correlation strength as weight
+	 * @param strength_pow raise correlation strength (after subtracting strength_floor) to
+	 *  this power before using as weight
+	 * @param debugLevel debug level
+	 * @param dbg_X supertile horizontal index to show debug information
+	 * @param dbg_Y supertile vertical index to show debug information
+	 * @return number of planes replaced
+	 */
+	// TODO: Mark replaced planes and reduce eigenvalue requirements for them (in case they are small)
+	// TODO:  planesSmooth() add power function to modify pull/neighbors relation (8 pull stronger than 1, but not 8 times stronger
+	public int replaceBrokenPlanes(
+			final TilePlanes.PlaneData[][]   planes,
+			final TilePlanes.PlaneData[][][] brokenPd,
+			final double                     max_diff, // maximal disparity difference (0 - any), will be normalized by dispNorm
+			final double                     other_diff, 
+			final boolean                    non_exclusive,
+			final boolean                    use_other_planes, // TODO:
+			final int                        measSel, // Select measurements for supertiles : +1 - combo, +2 - quad +4 - hor +8 - vert
+			final boolean                    allow_parallel,
+			
+			final boolean                    splitXY,
+			final double                     splitXYTolerance,
+			
+			// parameters to generate ellipsoids			
+			final double                     disp_far, // minimal disparity to select (or NaN)
+			final double                     disp_near, // maximal disparity to select (or NaN)
+			final double                     dispNorm,   //  Normalize disparities to the average if above
+			final double                     min_weight,
+			final int                        min_tiles,
+			// parameters to reduce outliers			
+			final double                     targetEigen,   //     =   0.1;  // Remove outliers until main axis eigenvalue (possibly scaled by plDispNorm) gets below
+			final double                     fractOutliers, //     =   0.3;  // Maximal fraction of outliers to remove
+			final int                        maxOutliers,   //     =   20;  // Maximal number of outliers to remove
+			final double                     strength_floor,
+			final double                     strength_pow,
+			final int                        debugLevel,
+			final int                        dbg_X,
+			final int                        dbg_Y)
+	{
+		final int tilesX =        tileProcessor.getTilesX();
+		final int superTileSize = tileProcessor.getSuperTileSize();
+		final int stilesX = (tilesX + superTileSize -1)/superTileSize;  
+		final int debug_stile = dbg_Y * stilesX + dbg_X;
+		final Thread[] threads = ImageDtt.newThreadArray(tileProcessor.threadsMax);
+		final AtomicInteger ai = new AtomicInteger(0);
+		final AtomicInteger ai_numThread = new AtomicInteger(0);
+		final int [] replaced = new int[threads.length];
+		final int [][] dbg_dirsYX = {{-1, 0},{-1,1},{0,1},{1,1},{1,0},{1,-1},{0,-1},{-1,-1}};
+		
+		for (int ithread = 0; ithread < threads.length; ithread++) {
+			threads[ithread] = new Thread() {
+				public void run() {
+					double [][] dbg_img = null;
+					String [] dbg_titles = null;
+					int numThread = ai_numThread.getAndIncrement(); // unique number of thread to write to rslt_diffs[numThread]
+					for (int nsTile = ai.getAndIncrement(); nsTile < planes.length; nsTile = ai.getAndIncrement()) {
+						int sty = nsTile / stilesX;  
+						int stx = nsTile % stilesX;
+						int dl = ((debugLevel > -1) && (nsTile == debug_stile)) ? 4:0;
+						int num_bplanes = 2; // brokenPd[nsTile].length;
+						int num_meas_layers = measuredLayers.getNumLayers();
+						int [] dbg_indices = {
+								0,
+								num_bplanes,
+								num_bplanes + 3 * num_meas_layers,
+								num_bplanes + 3 * num_meas_layers + 2 * (num_bplanes * num_meas_layers),
+								num_bplanes + 3 * num_meas_layers + 2 * (num_bplanes * num_meas_layers) + num_bplanes,
+						};
+
+						if ( brokenPd[nsTile] != null) {
+							if (dl > 0){
+								System.out.println("replaceBrokenPlanes nsTile="+nsTile);
+								dbg_titles = new String[dbg_indices[4]];
+								for (int i = 0; i < num_bplanes; i++){
+									dbg_titles[i + dbg_indices[0]]="src_plane_"+i;
+								}
+								for (int j = 0; j < num_meas_layers; j++){
+									dbg_titles[3* j + dbg_indices[1] + 0] =     "disp_"+j;
+									dbg_titles[3* j + dbg_indices[1] + 1] =     "mask_disp_"+j;
+									dbg_titles[3* j + dbg_indices[1] + 2] =     "str_"+j;
+								}
+								for (int i = 0; i < num_bplanes; i++){
+									for (int j = 0; j < num_meas_layers; j++){
+										dbg_titles[2 * (num_meas_layers * i +j) + dbg_indices[2] + 0] =  "t"+i+"_m"+j;
+										dbg_titles[2 * (num_meas_layers * i +j) + dbg_indices[2] + 1] =  "ft"+i+"_m"+j;
+									}
+								}
+								for (int i = 0; i < num_bplanes; i++){
+									dbg_titles[i + dbg_indices[3]]="plane_"+i;
+								}
+							}
+							
+							TilePlanes.PlaneData[][] bpd = new TilePlanes.PlaneData[brokenPd[nsTile].length][];
+							for (int np = 0; np < bpd.length; np++){
+								if (brokenPd[nsTile][np] != null){
+									bpd[np] = new TilePlanes.PlaneData[2];
+									bpd[np][0] = brokenPd[nsTile][np][0].clone();
+									bpd[np][1] = brokenPd[nsTile][np][1].clone(); // no need to clone [2]
+								}
+							}
+							if (splitXY) {
+								for (int np = 0; np < bpd.length; np++){
+									if (brokenPd[nsTile][np] != null){
+										// if it will fail (return false) no selection masks, still possible to separate
+										planes[nsTile][np].calcSelMasks(
+												bpd[np][0],
+												bpd[np][1],
+												splitXYTolerance, // double    tolerance,
+												min_tiles, 
+												dl); // int       debugLevel
+									}
+								}
+							}
+							
+							// now work on the pairs in bpd, null them out in case of failure
+							for (int np = 0; np < bpd.length; np++){
+								if (dl > 1) {
+									int ss2 = 2 * superTileSize;
+									dbg_img = new double [dbg_titles.length][];
+									for (int ni = 0; ni < num_bplanes ; ni++) {
+										dbg_img[ni + dbg_indices[0]] =  bpd[np][ni].getDoublePlaneDisparity(false);
+										boolean [] sm = bpd[np][ni].getSelMask();
+										// checkerboard pattern over disabled part of the plane
+										if (sm != null) {
+											for (int i = 0; i < sm.length; i++){
+												if (!sm[i] && ((((i % ss2)+ (i / ss2)) & 1) != 0)){
+													dbg_img[ni + dbg_indices[0]][i] = Double.NaN;
+												}
+											}
+										}
+									}
+									// add connections to neighbors
+									for (int dir = 0; dir <8; dir++){
+										int indx = (ss2 / 2) * (ss2 + 1);
+										int dindx = dbg_dirsYX[dir][0] * ss2 + dbg_dirsYX[dir][1];
+										for (int l = 0; l < (ss2 / 2 - 2); l++){
+											if ( l > 2){ // keep center not modified
+												for (int ni = 0; ni < num_bplanes ; ni++) {
+													if (bpd[np][ni].getNeibBest(dir) >= 0){
+														dbg_img[ni + dbg_indices[0]][indx] = Double.NaN;
+													}
+												}
+											}
+											indx += dindx; 
+										}
+									}
+									// show disparity (all and masked for this plane) and strength
+									for (int ml = 0; ml < num_meas_layers; ml++) {
+										double [][] dbg_disp_str = planes[nsTile][np].getMeasuredLayers().getDisparityStrength(
+												ml,                     // int num_layer,
+												stx,        // int stX,
+												sty,        // int stY,
+												null,       // boolean [] sel_in,
+												strength_floor,         //  double strength_floor,
+												strength_pow,  // double strength_pow,
+												true);                  // boolean null_if_none);
+										if (dbg_disp_str != null){
+											boolean [] dbg_msel = planes[nsTile][np].getMeasSelection(ml);
+											dbg_img[3 * ml + dbg_indices[1] + 0] = dbg_disp_str[0];
+											dbg_img[3 * ml + dbg_indices[1] + 2] = dbg_disp_str[1];
+											dbg_img[3 * ml + dbg_indices[1] + 1] = dbg_disp_str[0].clone();
+											if (dbg_msel != null) {
+												for (int i = 0; i < dbg_msel.length; i++){
+													if (!dbg_msel[i]) dbg_img[3 * ml + dbg_indices[1] + 1][i] = Double.NaN;
+												}
+											}
+										}
+									}
+									
+								} // (dl > 1)
+								
+								if (brokenPd[nsTile][np] != null){
+									// Is it a single set to replace current plane? If yes, use all tiles, not just original selection
+									int np_min = LOWEST_PLANE(planes[nsTile].length);
+									boolean single_plane = (np_min == (planes[nsTile].length - 1));
+									
+									Boolean OK = planes[nsTile][np].splitPlaneTiles ( // now always OK TODO: add layer mask or make use_other_planes int
+											bpd[np],          // PlaneData [] pd_set,
+											single_plane,     // boolean      single_plane,
+											max_diff,         // double       max_diff, // maximal disparity difference (0 - any), will be normalized by dispNorm
+											other_diff,
+											non_exclusive,    // boolean      non_exclusive,
+											use_other_planes, // boolean      use_other_planes,
+											measSel,          // int          measSel, // Select measurements for supertiles : +1 - combo, +2 - quad +4 - hor +8 - vert
+											allow_parallel,   //boolean      allow_parallel,
+											dl); // int          debugLevel)
+									if (!OK) {
+										brokenPd[nsTile][np] = null;
+										continue;
+									}
+									if (dl > 1) {// show initial selection (made of disparity clone - dbg_img[3 * ml + dbg_indices[1] + 0] = dbg_disp_str[0];
+										for (int ni = 0; ni < num_bplanes ; ni++) {
+											boolean [][] msel = new boolean [num_meas_layers][];
+											for (int ml = 0; ml < num_meas_layers; ml++) {
+												msel[ml] = bpd[np][ni].getMeasSelection(ml);
+												if (msel[ml] != null) {
+													dbg_img[2 * (num_meas_layers * ni +ml) + dbg_indices[2] + 0] =  dbg_img[3 * ml + dbg_indices[1] + 0].clone();
+													for (int i = 0; i < msel[ml].length; i++){
+														if (!msel[ml][i]) dbg_img[2 * (num_meas_layers * ni +ml) + dbg_indices[2] + 0][i] = Double.NaN;
+													}
+												}
+											}
+										}
+									}
+									
+									for (int npip = 0; npip < bpd[np].length; npip++) {
+										OK = bpd[np][npip].getPlaneFromMeas(
+												null,           // boolean [][] tile_sel, // null - do not use, {} use all (will be modified)
+												disp_far,       // double       disp_far, // minimal disparity to select (or NaN)
+												disp_near,      // double       disp_near, // maximal disparity to select (or NaN)
+												dispNorm,       // double       dispNorm,   //  Normalize disparities to the average if above
+												min_weight,     // double       min_weight,
+												min_tiles,      // int          min_tiles,
+												strength_floor, // double       strength_floor,
+												strength_pow,   // double       strength_pow,
+												dl);            // int          debugLevel)
+										if (!OK) break;
+										// remove outliers //removeOutliers
+										// now try to remove outliers
+										int max_outliers = (int) Math.round(bpd[np][npip].getNumPoints() * fractOutliers);
+										if (max_outliers > maxOutliers) max_outliers = maxOutliers;
+										double targetV = corrMaxEigen(
+												targetEigen,
+												dispNorm,
+												bpd[np][npip]);
+										if (bpd[np][npip].getValue() > targetV) {
+											OK = bpd[np][npip].removeOutliers( // getPlaneFromMeas should already have run
+													targetV,      // double     targetEigen, // target eigenvalue for primary axis (is disparity-dependent, so is non-constant)
+													max_outliers, // int        maxRemoved,  // maximal number of tiles to remove (not a constant)
+													dl); // int        debugLevel)
+											if (!OK) break;
+										}										
+									}
+									if (!OK) {
+										brokenPd[nsTile][np] = null;
+										continue;
+									}
+									
+									if (dl > 1) {// filtered selection (made of disparity clone - dbg_img[3 * ml + dbg_indices[1] + 0] = dbg_disp_str[0];
+										for (int ni = 0; ni < num_bplanes ; ni++) {
+											boolean [][] msel = new boolean [num_meas_layers][];
+											for (int ml = 0; ml < num_meas_layers; ml++) {
+												msel[ml] = bpd[np][ni].getMeasSelection(ml);
+												if (msel[ml] != null) {
+													dbg_img[2 * (num_meas_layers * ni +ml) + dbg_indices[2] +1] =  dbg_img[3 * ml + dbg_indices[1] + 0].clone();
+													for (int i = 0; i < msel[ml].length; i++){
+														if (!msel[ml][i]) dbg_img[2 * (num_meas_layers * ni +ml) + dbg_indices[2] + 1][i] = Double.NaN;
+													}
+												}
+											}
+										}
+										// show new planes
+										for (int ni = 0; ni < num_bplanes ; ni++) {
+											dbg_img[ni + dbg_indices[3]] =  bpd[np][ni].getDoublePlaneDisparity(false);
+										}
+										// actually show the image
+										
+										if (debugLevel > -1){
+											showDoubleFloatArrays sdfa_instance = new showDoubleFloatArrays();
+											sdfa_instance.showArrays(dbg_img, 2 * superTileSize, 2 * superTileSize, true, "replaceBrokenPlanes"+stx+"_y"+sty+"_p"+np, dbg_titles);
+										}
+									}
+								}
+							}
+							// now see if anything is left
+							int npairs = 0;
+							for (int np = 0; np < bpd.length; np++){
+								if (bpd[np] != null) npairs++;
+							}
+							
+							// Fixing that first plane in multi=-plane is not processed
+							int np_start = LOWEST_PLANE(2); // planes[nsTile].length);
+
+							if (npairs > 0){
+								TilePlanes.PlaneData[] old_planes = planes[nsTile];
+								planes[nsTile] = new TilePlanes.PlaneData[old_planes.length + npairs + np_start];
+								int npr = 0;
+								if (np_start > 0){
+//									planes[nsTile][npr++] = old_planes[0]; 
+									
+									planes[nsTile][npr++] = null; // see if there will be any conflicts, fis or replace with old_planes[0];  
+								}
+								for (int np = 0; np < bpd.length; np++){
+									if (bpd[np] != null) {
+										planes[nsTile][npr++] = bpd[np][0]; 
+										planes[nsTile][npr++] = bpd[np][1]; 
+									} else {
+										planes[nsTile][npr++] = old_planes[np]; 
+									}
+								}
+							}
+							replaced[numThread] += npairs;
+						}								
+					}
+				}
+			};
+		}		      
+		ImageDtt.startAndJoin(threads);
+		int replaced_planes = 0;
+		for (int i = 0; i < replaced.length; i++) replaced_planes += replaced[i];
+		return replaced_planes;
+	}
 	
 	/**
 	 * Create candidate planes to break a single plane in 2 by splitting consecutive connected
@@ -3419,7 +3887,7 @@ public class SuperTiles{
 													// to make each plane
 													// TODO save gain from splitting to planes 
 
-													if (debugLevel > -1){
+													if (debugLevel > 3){
 														showDoubleFloatArrays sdfa_instance = new showDoubleFloatArrays();
 														sdfa_instance.showArrays(dbg_img, 3*superTileSize, 3*superTileSize, true, "breakPlanesToPairs"+stx0+"_y"+sty0+"_p"+np0, titles);
 													}
@@ -3439,8 +3907,6 @@ public class SuperTiles{
 	}
 
 
-
-
 	public boolean [][] selectPlanes(
 			final double dispNorm,
 			final double maxEigen, // maximal eigenvalue of planes to consider
@@ -3460,7 +3926,7 @@ public class SuperTiles{
 				if (planes[nsTile] != null) {
 					selection[nsTile] = new boolean [planes[nsTile].length];
 					for (int np = 0; np < planes[nsTile].length; np++){
-						if (planes[nsTile][np].getWeight() >= minWeight){
+						if ((planes[nsTile][np] != null) && (planes[nsTile][np].getWeight() >= minWeight)){
 							//&& (planes[nsTile][np].getValue() < maxEigen)){
 							double eigVal = planes[nsTile][np].getValue();
 							double disp = planes[nsTile][np].getZxy()[0];

src/main/java/TilePlanes.java

@@ -64,15 +64,13 @@ public class TilePlanes {
 //		Lav = Math.sqrt((L1 * L1 * w1 + L2 * L2 * w2)/(w1 + w2))
 // worsening_12 = (L - Lav) * (w1 + w2) * (w1 + w2) / (Lav * x1 * w2)		
 		
-		
-		
-		
 		int         tileSize;
 		int         superTileSize;
 		int []      sTileXY =    null; // X and Y indices of this superTile in the image
 		
 		MeasuredLayers measuredLayers     =  null;
 		boolean [][] measuredSelection =     null; // [number of layers in measuredLayers][2*superTileSize * 2*superTileSize]
+		boolean   [] sel_mask          =     null; // selection mask - may be used for splitting plane along a line - each half can have mask
 		double       measured_strength_pow = 1.0;
 		double       strength_floor =        0.0;
 		double       min_weight =            0.0;  // minimal weight of the ellipsoid
@@ -87,6 +85,7 @@ public class TilePlanes {
 					this.tileSize,
 					this.superTileSize,
 					this.geometryCorrection);
+			pd.correctDistortions = this.correctDistortions;
 			pd.num_points = this.num_points; 
 			pd.weight =     this.weight; 
 			if (this.plane_sel != null)  pd.plane_sel =  this.plane_sel.clone(); 
@@ -104,14 +103,11 @@ public class TilePlanes {
 				pd.vectors[2] = this.vectors[2].clone();
 			}
 			if (this.measuredLayers != null) pd.measuredLayers = this.measuredLayers;
-			if (this.measuredSelection != null) {
-				pd.measuredSelection = this.measuredSelection.clone();
-				for (int i = 0; i < this.measuredSelection.length; i++){
-					if (this.measuredSelection[i] != null) {
-						pd.measuredSelection[i] = this.measuredSelection[i].clone();
-					}
-				}
-			}
+			
+			pd.setMeasSelection(this.measuredSelection);
+			
+			if (this.sel_mask != null) pd.sel_mask = this.sel_mask.clone();
+			
 			pd.measured_strength_pow = this.measured_strength_pow;
 			pd.strength_floor =        this.strength_floor;
 
@@ -125,6 +121,42 @@ public class TilePlanes {
 			return pd;
 		}
 		
+		public void setSelMask (boolean []sel_mask)
+		{
+			this.sel_mask = sel_mask;
+		}
+		public boolean [] getSelMask ()
+		{
+			return this.sel_mask;
+		}
+		
+		public void setMeasSelection(boolean [][] meas_sel)
+		{
+			if (meas_sel == null) 
+				this.measuredSelection = null;
+			else {
+				this.measuredSelection = meas_sel.clone();
+				for (int i = 0; i < meas_sel.length; i++){
+					if (meas_sel[i] != null) {
+						this.measuredSelection[i] = meas_sel[i].clone();
+					}
+				}
+			}
+		}
+		
+		public boolean [] getMeasSelection(int nl){
+			if (this.measuredSelection == null) {
+				return null;
+			}
+			return 	this.measuredSelection[nl];
+		}
+		
+		
+		public MeasuredLayers getMeasuredLayers()
+		{
+			return this.measuredLayers;
+		}
+		
 		public void copyNeib(
 				PlaneData src,
 				PlaneData dst)
@@ -180,6 +212,290 @@ public class TilePlanes {
 			this.measuredLayers = measuredLayers;
 			this.preferDisparity = preferDisparity;
 		}
+
+		/**
+		 * Create separation masks for two crossing planes. Uses neighbors connections to distinguish between concave and convex 
+		 * @param pd1 first plane data instance to create sel_mak
+		 * @param pd2 second plane data instance to create sel_mak
+		 * @param tolerance disparity tolerance for separation (will be normalized by dispNorm for large disparities)
+		 * @param min_tiles minimal number of tiles in each half
+		 * @param debugLevel debug level
+		 * @return true if OK, false if the planes are not crossing (masks not created)
+		 */
+		public boolean calcSelMasks(
+				PlaneData pd1,
+				PlaneData pd2,
+				double    tolerance,
+				int       min_tiles,
+				int       debugLevel
+				)
+		{
+			int st2 = 2 * superTileSize;
+			int [] dirs = {-st2, -st2 + 1, 1, st2 + 1, st2, st2 - 1, -1, -st2 - 1};
+			double [][] planes = {
+					pd1.getDoublePlaneDisparity(false),
+					pd2.getDoublePlaneDisparity(false)};
+			PlaneData [] pair = {pd1,pd2};
+			int len2 = planes[0].length;
+			boolean [][] sel_masks = new boolean [2][len2];
+			double concave = 0; // positive - concave, negative - convex
+			int cent_indx = (st2/2 -1) * (st2+1);
+			for (int np = 0; np < 2; np++){
+				int [] neibs = pair[np].getNeibBest();
+				for (int dir = 0; dir < dirs.length; dir++) {
+					if (neibs[dir] >= 0) { // neighbor connected
+						int indx = cent_indx + dirs[dir];
+						concave += (planes[np][indx] - planes[np][cent_indx]) - (planes[1 - np][indx] - planes[1 - np][cent_indx]);
+					}
+				}
+			}
+			// swap plane selections for convex plane intersection (edge is the closest)
+			int first = (concave < 0) ? 1 : 0;
+			int [] nums = {0, 0};
+			
+			for (int i = 0; i <len2; i++){
+				double d_av = 0.5 * (planes[0][i] + planes[1][i]);
+				double diff = planes[0][i] - planes[1][i];
+				if ((dispNorm > 0.0) && (d_av > dispNorm)){
+					diff *= dispNorm/d_av;
+				}
+				if (diff > -tolerance) {
+					sel_masks[first][i] = true;
+					nums[first] ++;
+				}
+				if (-diff > -tolerance) {
+					sel_masks[1-first][i] = true;
+					nums[1-first] ++;
+				}
+			}
+			if ((nums[0] < min_tiles) || (nums[1] < min_tiles)){
+				return false;
+			}
+			// apply selections
+			pd1.setSelMask(sel_masks[0]);
+			pd2.setSelMask(sel_masks[1]);
+			return true;
+		}
+		
+		/**
+		 * Spit tiles belonging to this between multiple PlaneData instances
+		 * @param pd_set array of plane data instances
+		 * @param single_plane it is a single plane to split - use all tiles, not just previously
+		 *  selected
+		 * @param max_diff maximal normalized disparity difference from the plane to consider
+		 * @param other_diff maximal difference of the added tile ratio to the average
+		 *  disparity difference of the exclusively selected tiles
+		 * 
+		 * @param non_exclusive allow the same tile data to belong to multiple PD instances
+		 * @param use_other_planes allow the same tile not included in this PD to be used
+		 * @param measSel (with use_other_planes) select measurements for supertiles :
+		 *  +1 - combo, +2 - quad +4 - hor +8 - vert
+		 * @param allow_parallel allow parallel shift of each plane before adding more data
+		 * @param debugLevel debug level
+		 * @return true if OK
+		 */
+		
+		public boolean splitPlaneTiles (
+				PlaneData [] pd_set,
+				boolean      single_plane,
+				double       max_diff, // maximal disparity difference (0 - any), will be normalized by dispNorm
+				double       other_diff, 
+				boolean      non_exclusive,
+				boolean      use_other_planes,
+				int          measSel, // Select measurements for supertiles : +1 - combo, +2 - quad +4 - hor +8 - vert
+				boolean      allow_parallel,
+				int          debugLevel)
+		{
+			if (debugLevel>0){
+				System.out.println("splitPlaneTiles");
+			}
+			
+			double [][] planes = new double [pd_set.length][];
+			boolean [][][] tsel = new boolean [pd_set.length][measuredLayers.getNumLayers()][];
+			boolean [][] sel_masks = new boolean [pd_set.length][];
+			for (int np = 0; np < pd_set.length; np ++) if (pd_set[np] != null){
+				planes[np] = pd_set[np].getDoublePlaneDisparity(false);
+				sel_masks[np] = pd_set[np].getSelMask();
+				for (int nl = 0; nl < measuredLayers.getNumLayers(); nl++){
+					if (measuredSelection[nl] != null){
+						tsel[np][nl] = new boolean [measuredSelection[nl].length];
+					}
+				}
+//				pd_set[np].setMeasSelection(measuredSelection); // copy same selection to all planes
+			}
+			double max_diff2 = max_diff * max_diff;
+			double other_diff2 =other_diff * other_diff; 
+			double [] sd2 = new double[pd_set.length];
+			double [] sd2_av = new double[pd_set.length];
+			double [] sw =  new double[pd_set.length];
+			double [][][] disp_strength = new double[measuredLayers.getNumLayers()][][];
+			// split exclusively, calculate rms for each, then add others if RMS is not increased
+			for (int nl = 0; nl < measuredLayers.getNumLayers(); nl ++){
+				if ((measuredSelection[nl] != null) &&  ((measSel & (1 << nl)) !=0)) {
+					disp_strength[nl] = measuredLayers.getDisparityStrength(
+							nl,                     // int num_layer,
+							getSTileXY()[0],        // int stX,
+							getSTileXY()[1],        // int stY,
+							(single_plane ? null : measuredSelection[nl]),  // boolean [] sel_in,
+							strength_floor,         //  double strength_floor,
+							measured_strength_pow,  // double strength_pow,
+							true);                  // boolean null_if_none);
+					if (disp_strength[nl] != null) {
+						for (int indx = 0; indx < disp_strength[nl][1].length; indx++){
+							double w = disp_strength[nl][1][indx];
+							if (w > 0.0){
+								double d = disp_strength[nl][0][indx];
+								double d2_best = Double.NaN;
+								int np_best = 0;
+								for (int np = 0; np < planes.length; np++) if (planes[np] != null){
+									if ((sel_masks[np] == null) || sel_masks[np][indx]) {
+										double d_av = 0.5 * (d + planes[np][indx]);
+										double d2 = (d - planes[np][indx]);
+										if ((dispNorm > 0.0) && (d_av > dispNorm)){
+											d2 *= dispNorm/d_av;
+										}
+										d2 *= d2;
+										if (!(d2 >= d2_best)){ // so d2_best NaN is OK
+											d2_best = d2;
+											np_best = np;
+										}
+									}
+								}
+								if ((max_diff2 == 0.0) || (d2_best < max_diff2)) {
+									tsel[np_best][nl][indx] = true;
+									sd2[np_best] += w * d2_best;
+									sw[np_best] += w;
+								}
+							}
+						}
+					}
+				}
+			}
+			for (int np = 0; np < planes.length; np++) if (planes[np] != null){
+				if (sw[np] > 0.0){
+					sd2_av[np] = sd2[np] / sw[np];
+				}
+			}
+			double [] sd = new double[pd_set.length]; // will be all 0.0;
+
+			if (allow_parallel && (non_exclusive || use_other_planes)){
+				for (int nl = 0; nl < measuredLayers.getNumLayers(); nl ++){
+					if (disp_strength[nl] != null) {
+						//disp_strength[nl] = measuredLayers.getDisparityStrength(
+						for (int indx = 0; indx < disp_strength[nl][1].length; indx++){
+							double w = disp_strength[nl][1][indx];
+							if (w > 0.0){
+								double d = disp_strength[nl][0][indx];
+								for (int np = 0; np < planes.length; np++) if (planes[np] != null){
+									if (!tsel[np][nl][indx]) { // not already in that plane
+										sd[np] += w * d;
+									}
+								}
+							}
+						}
+					}
+				}
+				for (int np = 0; np < planes.length; np++) if (planes[np] != null){
+					if (sw[np] > 0.0){
+						sd[np] /= sw[np];
+						sd2_av[np] -= sd[np] * sd[np];
+						sd2[np] = sd2_av[np] * sw[np]; // to add more 
+					}
+				}
+				
+			}
+			
+			
+			if (non_exclusive) {
+				for (int nl = 0; nl < measuredLayers.getNumLayers(); nl ++){
+					if ((measuredSelection[nl] != null) &&  ((measSel & (1 << nl)) !=0)) {
+						if (disp_strength[nl] != null) {
+							//disp_strength[nl] = measuredLayers.getDisparityStrength(
+							for (int indx = 0; indx < disp_strength[nl][1].length; indx++){
+								double w = disp_strength[nl][1][indx];
+								if (w > 0.0){
+									double d = disp_strength[nl][0][indx];
+									for (int np = 0; np < planes.length; np++) if (planes[np] != null){
+										if ((sel_masks[np] == null) || sel_masks[np][indx]) {
+											if (!tsel[np][nl][indx]) { // not already in that plane
+												double d_av = 0.5 * (d - planes[np][indx]);
+												double d2 = (d - planes[np][indx]);
+												d2 -= sd[np]; // subtract parallel shift
+												if ((dispNorm > 0.0) && (d_av > dispNorm)){
+													d2 *= dispNorm/d_av;
+												}
+												d2 *= d2;
+												if (d2 <= sd2_av[np] * other_diff2) { // not more than exclusive tile variance 
+													tsel[np][nl][indx] = true;
+													sd2[np] += w * d2;
+													sw[np] += w;
+												}
+											}
+										}
+									}
+								}
+							}					
+						}
+					}
+				}
+			}
+			
+			if (use_other_planes && !single_plane) { // no need if single_plane - it already got all planes it could
+				for (int nl = 0; nl < measuredLayers.getNumLayers(); nl ++){
+					if ((measuredSelection[nl] != null) &&  ((measSel & (1 << nl)) !=0)) {
+						// recalculate for all measure tiles, not just selected in the original PD
+						disp_strength[nl] = measuredLayers.getDisparityStrength(
+								nl,                     // int num_layer,
+								getSTileXY()[0],        // int stX,
+								getSTileXY()[1],        // int stY,
+								null,                   // boolean [] sel_in, null here - all measured data
+								strength_floor,         //  double strength_floor,
+								measured_strength_pow,  // double strength_pow,
+								true);                  // boolean null_if_none);
+						//disp_strength[nl] = measuredLayers.getDisparityStrength(
+						for (int indx = 0; indx < disp_strength[nl][1].length; indx++){
+							double w = disp_strength[nl][1][indx];
+							if (w > 0.0){
+								double d = disp_strength[nl][0][indx];
+								for (int np = 0; np < planes.length; np++) if (planes[np] != null){
+									if ((sel_masks[np] == null) || sel_masks[np][indx]) {
+										if (!tsel[np][nl][indx]) { // not already in that plane
+											double d_av = 0.5 * (d - planes[np][indx]);
+											double d2 = (d - planes[np][indx]);
+											d2 -= sd[np]; // subtract parallel shift
+											if ((dispNorm >0.0) && (d_av > dispNorm)){
+												d2 *= dispNorm/d_av;
+											}
+											d2 *= d2;
+											if (d2 <= sd2_av[np]  * other_diff2) { // not more than exclusive tile variance 
+												tsel[np][nl][indx] = true;
+												sd2[np] += w * d2;
+												sw[np] += w;
+											}
+										}
+									}
+								}
+							}
+						}					
+					}
+				}
+			}
+			
+			// re-calculate variance (just for debug, not needed
+			for (int np = 0; np < planes.length; np++)  if (planes[np] != null){
+				if (sw[np] > 0.0){
+					sd2_av[np] = sd2[np] / sw[np];
+				}
+			}
+			// apply selections to each PD
+			for (int np = 0; np < planes.length; np++)  if (planes[np] != null){
+				pd_set[np].setMeasSelection(tsel[np]);
+			}
+			// need to re-calculate new planes, remove outliers
+			return true;
+		}
+		
+		
 		/**
 		 * Remove outliers from the set of tiles contributing to a single plane ellipsoid
 		 * Should run after getPlaneFromMeas as some parameter4s will be copied from that run
@@ -327,7 +643,7 @@ public class TilePlanes {
 			this.min_weight =             min_weight;
 			this.min_tiles =              min_tiles;
 			this.dispNorm =               dispNorm;
-			if (debugLevel > 1){
+			if (debugLevel > 3){
 				System.out.println("getPlaneFromMeas()");
 			}
 

src/main/java/TileProcessor.java

@@ -2753,7 +2753,8 @@ public class TileProcessor {
 					clt_parameters.stMaxDisparity,   // double     max_disparity,
 					clt_parameters.stFloor,          // double     strength_floor,
 					clt_parameters.stPow,            // double     strength_pow,
-					0.0); // NO BLUR double     stBlurSigma)
+					0.0,// NO BLUR double     stBlurSigma)
+					clt_parameters.stMeasSel); // bitmask of the selected measurements for supertiles : +1 - combo, +2 - quad +4 - hor +8 - vert
 			dbg_hist[0] = scan_prev.showDisparityHistogram();
 
 			scan_prev.setSuperTiles(
@@ -2764,7 +2765,8 @@ public class TileProcessor {
 					clt_parameters.stMaxDisparity,   // double     max_disparity,
 					clt_parameters.stFloor,          // double     strength_floor,
 					clt_parameters.stPow,            // double     strength_pow,
-					clt_parameters.stSigma); // with blur double     stBlurSigma)
+					clt_parameters.stSigma, // with blur double     stBlurSigma)
+					clt_parameters.stMeasSel); // bitmask of the selected measurements for supertiles : +1 - combo, +2 - quad +4 - hor +8 - vert
 			dbg_hist[1] = scan_prev.showDisparityHistogram();
 
 			dbg_hist[2] = scan_prev.showMaxMinMax();
@@ -2942,6 +2944,50 @@ public class TileProcessor {
 	{
 		trimCLTPasses(); // make possible to run this method multiple time - remove extra passes added by it last time		
 		CLTPass3d scan_prev = clt_3d_passes.get(clt_3d_passes.size() -1); // get last one
+
+		boolean show_st =    clt_parameters.stShow || (debugLevel > 1);
+		// recalculate supertiles (may be removed later)
+		//		if (use_supertiles || show_st) {	
+		String [] dbg_st_titles = {"raw", "blurred"+clt_parameters.stSigma,"max-min-max"}; 
+		double [][] dbg_hist = new double[dbg_st_titles.length][];
+		if (show_st) { // otherwise only blured version is needed 
+			scan_prev.setSuperTiles(
+					clt_parameters.stStepNear,       // double     step_disparity,
+					clt_parameters.stStepFar,        // double     step_near,
+					clt_parameters.stStepThreshold,  // double     step_threshold,
+					clt_parameters.stMinDisparity,   // double     min_disparity,
+					clt_parameters.stMaxDisparity,   // double     max_disparity,
+					clt_parameters.stFloor,          // double     strength_floor,
+					clt_parameters.stPow,            // double     strength_pow,
+					0.0, // NO BLUR double     stBlurSigma)
+					clt_parameters.stMeasSel); // bitmask of the selected measurements for supertiles : +1 - combo, +2 - quad +4 - hor +8 - vert
+			dbg_hist[0] = scan_prev.showDisparityHistogram();
+		}
+
+//		SuperTiles st = 
+				scan_prev.setSuperTiles(
+				clt_parameters.stStepNear,       // double     step_disparity,
+				clt_parameters.stStepFar,        // double     step_near,
+				clt_parameters.stStepThreshold,  // double     step_threshold,
+				clt_parameters.stMinDisparity,   // double     min_disparity,
+				clt_parameters.stMaxDisparity,   // double     max_disparity,
+				clt_parameters.stFloor,          // double     strength_floor,
+				clt_parameters.stPow,            // double     strength_pow,
+				clt_parameters.stSigma,          // with blur double     stBlurSigma)
+				clt_parameters.stMeasSel); // bitmask of the selected measurements for supertiles : +1 - combo, +2 - quad +4 - hor +8 - vert
+		if (show_st) { // otherwise only blured version is needed 
+			dbg_hist[1] = scan_prev.showDisparityHistogram();
+			dbg_hist[2] = scan_prev.showMaxMinMax();
+		}
+
+
+		if (show_st){
+			int hist_width0 =  scan_prev.showDisparityHistogramWidth();
+			int hist_height0 = dbg_hist[0].length/hist_width0;
+			showDoubleFloatArrays sdfa_instance = new showDoubleFloatArrays(); // just for debugging?
+			sdfa_instance.showArrays(dbg_hist, hist_width0, hist_height0, true, "disparity_supertiles_histograms",dbg_st_titles);
+		}
+		
 		SuperTiles st = scan_prev.getSuperTiles();
 		
 // moved here
@@ -2949,7 +2995,7 @@ public class TileProcessor {
 		st.processPlanes3(
 				null, // final boolean [] selected, // or null
 				0.3, // final double     min_disp,
-				false, // final boolean    invert_disp, // use 1/disparity
+				clt_parameters.stMeasSel, //            =     1   //Select measurements for supertiles : +1 - combo, +2 - quad +4 - hor +8 - vert
 				clt_parameters.plDispNorm, //           =   2.0;  // Normalize disparities to the average if above
 				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
@@ -3013,6 +3059,84 @@ public class TileProcessor {
 			if (num_added == 0) break;
 		}
 		
+		TilePlanes.PlaneData[][][]       split_planes =   // use original (measured planes. See if smoothed are needed here)
+				st.breakPlanesToPairs(
+				st.getPlanes(), // Mod(),             // final TilePlanes.PlaneData[][] center_planes, // measured_planes,
+				st.getPlanes(), // Mod(),             // final TilePlanes.PlaneData[][] neib_planes,   //mod_planes,
+				clt_parameters.plSplitPull ,          // final double                   center_pull,
+				clt_parameters.plSplitMinNeib ,       // min_neibs, // 2
+				clt_parameters.plSplitMinWeight,      // final double splitMinWeight,  //     =  2.0;  // Minimal weight of split plains to show
+				clt_parameters.plSplitMinQuality,     // final double splitMinQuality, //    =  1.1;  // Minimal split quality to show
+
+				clt_parameters.plPreferDisparity,
+				1,                                   // final int debugLevel)
+				clt_parameters.tileX,
+				clt_parameters.tileY); 
+		
+		if (clt_parameters.plSplitApply) {
+			int numSplitPlanes = st.replaceBrokenPlanes(
+					st.getPlanes(),                 // final TilePlanes.PlaneData[][]   planes,
+					split_planes,                    // final TilePlanes.PlaneData[][][] brokenPd,
+					clt_parameters.plMaxDiff,        //  final double                     max_diff, // maximal disparity difference (0 - any), will be normalized by dispNorm
+					clt_parameters.plOtherDiff,      //  final double                     other_diff, // other_diff maximal difference of the added tile ratio to the average disparity difference
+					clt_parameters.plNonExclusive,   // final boolean                    non_exclusive,
+					clt_parameters.plUseOtherPlanes, // final boolean                    use_other_planes, // TODO:
+					clt_parameters.stMeasSel,        // final int                        measSel, // Select measurements for supertiles : +1 - combo, +2 - quad +4 - hor +8 - vert
+					clt_parameters.plAllowParallel,  // final boolean                    allow_parallel,
+					clt_parameters.plSplitXY,        // final boolean                    splitXY,
+					clt_parameters.plSplitXYTolerance, // final double                   splitXYTolerance,
+					// parameters to generate ellipsoids			
+					0.0, // 3,                             // final double                     disp_far, // minimal disparity to select (or NaN)
+					Double.NaN,                      // final double                     disp_near, // maximal disparity to select (or NaN)
+					clt_parameters.plDispNorm,       // final double                     dispNorm,   //  Normalize disparities to the average if above
+					0.0,                             // final double                     min_weight,
+					clt_parameters.plMinPoints,      // final int                        min_tiles,
+					// parameters to reduce outliers			
+					clt_parameters.plTargetEigen,    // final double                     targetEigen,   //     =   0.1;  // Remove outliers until main axis eigenvalue (possibly scaled by plDispNorm) gets below
+					clt_parameters.plFractOutliers,  // final double                     fractOutliers, //     =   0.3;  // Maximal fraction of outliers to remove
+					clt_parameters.plMaxOutliers,    // final int                        maxOutliers,   //     =   20;  // Maximal number of outliers to remove
+					clt_parameters.stFloor,          // final double                     strength_floor,
+					clt_parameters.stPow,            // final double                     strength_pow,
+					1,                               // final int debugLevel)
+					clt_parameters.tileX,
+					clt_parameters.tileY);
+			if (debugLevel > -1){
+				System.out.println("replaceBrokenPlanes(): Replaced " + numSplitPlanes + " planes by pairs.");
+			}
+			// now re-create connections and best neighbors
+			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);
+			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")  
@@ -3022,7 +3146,8 @@ 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.plMaxEigen,                    // final double      maxValue, // do not combine with too bad planes
 					clt_parameters.plIterations,                  // final int         num_passes,
-					clt_parameters.plStopBad,                     //Do not update supertile if any of connected neighbors is not good (false: just skip that neighbor)
+					clt_parameters.plStopBad,                     // Do not update supertile if any of connected neighbors is not good (false: just skip that neighbor)
+					clt_parameters.plNormPow,                     // 0.0: 8 neighbors pull 8 times as 1, 1.0 - same as 1
 					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)
@@ -3059,6 +3184,7 @@ public class TileProcessor {
 		// save surfaces with SuperTiles instance. They can be used to snap to for the per-tile disparity maps.
 		st.setSurfaces(surfaces);
 
+/*
   		TilePlanes.PlaneData[][][]       split_planes =  
 				st.breakPlanesToPairs(
 				st.getPlanes(), // Mod(),             // final TilePlanes.PlaneData[][] center_planes, // measured_planes,
@@ -3072,7 +3198,7 @@ public class TileProcessor {
 				1,                                   // final int debugLevel)
 				clt_parameters.tileX,
 				clt_parameters.tileY); 
-		
+*/		
 		
 		if (clt_parameters.show_planes){
 			double [] split_lines = st.showSplitLines(
@@ -3398,7 +3524,8 @@ public class TileProcessor {
 					clt_parameters.stMaxDisparity,   // double     max_disparity,
 					clt_parameters.stFloor,          // double     strength_floor,
 					clt_parameters.stPow,            // double     strength_pow,
-					0.0); // NO BLUR double     stBlurSigma)
+					0.0, // NO BLUR double     stBlurSigma)
+					clt_parameters.stMeasSel); // bitmask of the selected measurements for supertiles : +1 - combo, +2 - quad +4 - hor +8 - vert
 			dbg_hist[0] = scan_prev.showDisparityHistogram();
 
 			SuperTiles st = scan_prev.setSuperTiles(
@@ -3409,7 +3536,8 @@ public class TileProcessor {
 					clt_parameters.stMaxDisparity,   // double     max_disparity,
 					clt_parameters.stFloor,          // double     strength_floor,
 					clt_parameters.stPow,            // double     strength_pow,
-					clt_parameters.stSigma); // with blur double     stBlurSigma)
+					clt_parameters.stSigma, // with blur double     stBlurSigma)
+					clt_parameters.stMeasSel); // bitmask of the selected measurements for supertiles : +1 - combo, +2 - quad +4 - hor +8 - vert
 			dbg_hist[1] = scan_prev.showDisparityHistogram();
 
 			dbg_hist[2] = scan_prev.showMaxMinMax();
@@ -3612,84 +3740,6 @@ public class TileProcessor {
 			
 		}		
 		
-/**		
-		int [] neighbors = dp.getNeighbors( // creates neighbors mask from bitmask
-				grown, // these_tiles, // grown, // these_tiles, // boolean [] selected,
-				tilesX);
-		//	int [] neighbors_orig = neighbors.clone();
-		double [] dbg_neib = dp.dbgShowNeighbors(
-				grown, // these_tiles, // grown, // these_tiles,
-				neighbors, // _orig, // int [] neighbors,
-				clt_parameters.transform_size, // int    tile_size,
-				-1.0, // double bgnd,
-				1.0); // double fgnd)
-
-		double [] new_disparity = this_disparity.clone();
-		double [][]dbgDeriv = new double [2][]; // [these_tiles.length];
-		//	sdfa_instance.showArrays(dbg_neib,tilesX*clt_parameters.transform_size, tilesY*clt_parameters.transform_size,"neighbors");
-		dp.smoothDisparity(
-				clt_parameters.tiDispPull,   // final double     dispPull, // clt_parameters.tiDispPull or 0.0
-				3,                           // final int        mask,     // 1 - work on internal elements, 2 - on border elements, 3 - both (internal first);
-				clt_parameters.tiIterations, //  final int        num_passes,
-				Math.pow(10.0,  -clt_parameters.tiPrecision), // final double     maxDiff, // maximal change in any of the disparity values
-				neighbors,                   // final int     [] neighbors, // +1 - up (N), +2 - up-right - NE, ... +0x80 - NW
-				new_disparity,               // final double  [] disparity,          // current disparity value
-				this_disparity,              // final double  [] measured_disparity, // measured disparity
-				this_strength,               // final double  [] strength,
-				null, // this_hor_disparity,          // final double     hor_disparity, // not yet used
-				null, // hor_strength_conv,           // final double     hor_strength, // not yet used
-				these_tiles,                 // grown, // these_tiles,                 // final boolean [] selected,
-				border, // final boolean [] border,
-				clt_parameters,
-				//			dbgDeriv,                    // final double [][] dbgDeriv, //double [2][len] or null;
-				threadsMax,                  // maximal number of threads to launch                         
-				debugLevel);
-		int [][][] clustersNO= dp.extractNonOlerlap(
-				true, // diag_en,
-				neighbors, // +1 - up (N), +2 - up-right - NE, ... +0x80 - NW
-				grown, // these_tiles, // final boolean []  selected, // only inner?
-				border,          // border should be diagonal!
-				threadsMax,      // maximal number of threads to launch                         
-				debugLevel);
-
-		
-		int numScans= 0;
-
-		if (clt_parameters.shUseFlaps) {
-			numScans =  createTileOverlapTasks(
-					clt_parameters.max_clusters, // 50,                                           // int       maxClusters,
-					clt_parameters.shMinArea,                     // int       minClusterArea,
-					new_disparity,                                // [] disparity_in, masked ok too
-					this_strength,  // double [] strength_in,
-					Math.pow(10.0,  -clt_parameters.tiPrecision), // double      maxChange,   // adjust border disparity until change is below this.
-					clt_parameters.shMinStrength,                 // double    minStrength,
-					clustersNO,                                   // int []    clusters_in,
-					disparity_far,
-					disparity_near,
-					clt_parameters.show_shells,
-					debugLevel);
-
-		} else {
-			int [] enum_clusters = enumerateClusters(
-					true, // boolean diag_en,
-					grown); // these_tiles); // boolean [] tiles_src)
-
-			numScans =  createTileTasks(
-					50, // int       maxClusters,
-					0,  // int       minClusterArea,
-					new_disparity, // this_disparity, //  [] disparity_in, masked ok too
-					this_strength,  // double [] strength_in,
-					0.0,            // double    minStrength,
-					enum_clusters, // int []    clusters_in,
-					disparity_far,
-					disparity_near,
-					debugLevel);
-		}
-
-		if (debugLevel > -1){
-			System.out.println("thirdPassSetup(): created "+ numScans+ " FPGA passes.");
-		}
-*/		
 	}
 	
 	
@@ -4281,7 +4331,8 @@ public class TileProcessor {
 					clt_parameters.stMaxDisparity,   // double     max_disparity,
 					clt_parameters.stFloor,          // double     strength_floor,
 					clt_parameters.stPow,            // double     strength_pow,
-					0.0); // NO BLUR double     stBlurSigma)
+					0.0,// NO BLUR double     stBlurSigma)
+					clt_parameters.stMeasSel); // bitmask of the selected measurements for supertiles : +1 - combo, +2 - quad +4 - hor +8 - vert
 			dbg_hist[0] = scan_prev.showDisparityHistogram();
 
 			scan_prev.setSuperTiles(
@@ -4292,7 +4343,8 @@ public class TileProcessor {
 					clt_parameters.stMaxDisparity,   // double     max_disparity,
 					clt_parameters.stFloor,          // double     strength_floor,
 					clt_parameters.stPow,            // double     strength_pow,
-					clt_parameters.stSigma); // with blur double     stBlurSigma)
+					clt_parameters.stSigma, // with blur double     stBlurSigma)
+					clt_parameters.stMeasSel); // bitmask of the selected measurements for supertiles : +1 - combo, +2 - quad +4 - hor +8 - vert
 			dbg_hist[1] = scan_prev.showDisparityHistogram();
 
 			dbg_hist[2] = scan_prev.showMaxMinMax();