added new surface generation, working on resolving 'impossible triangle' conflicts

src/main/java/EyesisCorrectionParameters.java

@@ -2171,7 +2171,8 @@ public class EyesisCorrectionParameters {
   		public double     plMaxEigen           =   0.3;  // Maximal eigenvalue of a plane 
   		public boolean    plDbgMerge           =   true; // Combine 'other' plane with current
   		public double     plWorstWorsening     =   3.0;  // Worst case worsening after merge
-  		public double     plOKMergeEigen       =   0.01; // If result of the merged planes is below, OK to bypass worst worsening
+  		public double     plOKMergeEigen       =   0.03; // If result of the merged planes is below, OK to bypass worst worsening
+  		public double     plMaxWorldSin2       =   0.1;  // Maximal sine squared of the world angle between planes to merge. Set to >= 1.0 to disable
   		public double     plWeakWorsening      =   1.0;  // Relax merge requirements for weaker planes
   		public boolean    plMutualOnly         =   true; // keep only mutual links, remove weakest if conflict
   		public boolean    plFillSquares        =   true; // Add diagonals to full squares
@@ -2217,12 +2218,13 @@ public class EyesisCorrectionParameters {
   		public double     msScaleProj          =   1.5;  // Scale projection of the plane ellipsoid 
   		public double     msFractUni           =   0.3;  // Spread this fraction of the ellipsoid weight among extended (double) supertile
   		
-		public boolean    tsNoEdge             = true;   // Do not assigned tiles to thesurface edges (not having all 8 neighbors)
+		public boolean    tsNoEdge             = true;   // Do not assign tiles to the surface edges (not having all 8 neighbors)
+		public boolean    tsUseCenter          = true;   // Only assign outside of 8x8 center if no suitable alternative
   		public double     tsMaxDiff            = 0.3;    // Maximal disparity difference when assigning tiles
 		public double     tsMinDiffOther       = 0.35;   // Minimal disparity difference to be considered as a competitor surface
 		public double     tsMinStrength        = 0.05;   // Minimal tile correlation strength to be assigned
 		public double     tsMaxStrength        = 10.0;   // Maximal tile correlation strength to be assigned
-		public double     tsMinSurface         = 0.0001; // Minimal surface strength at the tile location
+		public double     tsMinSurface         = 0.001;  // Minimal surface strength at the tile location
 		public int        tsMoveDirs           = 3;      // Allowed tile disparity correction: 1 increase, 2 - decrease, 3 - both directions
 		public double     tsSurfStrPow         = 0.0;    // Raise surface strengths ratio to this power when comparing candidates
 		public double     tsAddStrength        = 0.01;   // Add to strengths when calculating pull of assigned tiles
@@ -2517,6 +2519,7 @@ public class EyesisCorrectionParameters {
 			properties.setProperty(prefix+"plDbgMerge",       this.plDbgMerge+"");
 			properties.setProperty(prefix+"plWorstWorsening", this.plWorstWorsening +"");
 			properties.setProperty(prefix+"plOKMergeEigen",   this.plOKMergeEigen +"");
+			properties.setProperty(prefix+"plMaxWorldSin2",   this.plMaxWorldSin2 +"");
 			properties.setProperty(prefix+"plWeakWorsening",  this.plWeakWorsening +"");
 			properties.setProperty(prefix+"plMutualOnly",     this.plMutualOnly+"");
 			properties.setProperty(prefix+"plFillSquares",    this.plFillSquares+"");
@@ -2557,6 +2560,7 @@ public class EyesisCorrectionParameters {
 			properties.setProperty(prefix+"msFractUni",       this.msFractUni +"");
 
 			properties.setProperty(prefix+"tsNoEdge",         this.tsNoEdge+"");
+			properties.setProperty(prefix+"tsUseCenter",      this.tsUseCenter+"");
 			properties.setProperty(prefix+"tsMaxDiff",        this.tsMaxDiff +"");
 			properties.setProperty(prefix+"tsMinDiffOther",   this.tsMinDiffOther +"");
 			properties.setProperty(prefix+"tsMinStrength",    this.tsMinStrength +"");
@@ -2842,6 +2846,7 @@ public class EyesisCorrectionParameters {
   			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+"plOKMergeEigen")!=null)    this.plOKMergeEigen=Double.parseDouble(properties.getProperty(prefix+"plOKMergeEigen"));
+  			if (properties.getProperty(prefix+"plMaxWorldSin2")!=null)    this.plMaxWorldSin2=Double.parseDouble(properties.getProperty(prefix+"plMaxWorldSin2"));
   			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"));
 
@@ -2885,6 +2890,7 @@ public class EyesisCorrectionParameters {
 
 
   			if (properties.getProperty(prefix+"tsNoEdge")!=null)          this.tsNoEdge=Boolean.parseBoolean(properties.getProperty(prefix+"tsNoEdge"));
+  			if (properties.getProperty(prefix+"tsUseCenter")!=null)       this.tsUseCenter=Boolean.parseBoolean(properties.getProperty(prefix+"tsUseCenter"));
   			if (properties.getProperty(prefix+"tsMaxDiff")!=null)         this.tsMaxDiff=Double.parseDouble(properties.getProperty(prefix+"tsMaxDiff"));
   			if (properties.getProperty(prefix+"tsMinDiffOther")!=null)    this.tsMinDiffOther=Double.parseDouble(properties.getProperty(prefix+"tsMinDiffOther"));
   			if (properties.getProperty(prefix+"tsMinStrength")!=null)     this.tsMinStrength=Double.parseDouble(properties.getProperty(prefix+"tsMinStrength"));
@@ -3196,6 +3202,7 @@ public class EyesisCorrectionParameters {
   			gd.addCheckbox    ("Combine 'other' plane with the current (unused)",                              this.plDbgMerge);
   			gd.addNumericField("Worst case worsening after merge",                                             this.plWorstWorsening,  6);
   			gd.addNumericField("If result of the merged planes is below, OK to bypass worst worsening",        this.plOKMergeEigen,  6);
+  			gd.addNumericField("Maximal sine squared of the world angle between planes to merge. Set to >= 1.0 to disable", this.plMaxWorldSin2,  6);
   			gd.addNumericField("Relax merge requirements for weaker planes",                                   this.plWeakWorsening,  6);
   			gd.addCheckbox    ("Keep only mutual links, remove weakest if conflict",                           this.plMutualOnly);
 
@@ -3239,8 +3246,9 @@ public class EyesisCorrectionParameters {
   			gd.addNumericField("Spread this fraction of the ellipsoid weight among extended (double) supertile",this.msFractUni,  6);
 
   			gd.addMessage     ("--- Tiles assignment ---");
-  			
-  			gd.addCheckbox    ("Do not assigned tiles to the surface edges (not having all 8 neighbors)",         this.tsNoEdge);
+  			 
+  			gd.addCheckbox    ("Do not assign tiles to the surface edges (not having all 8 neighbors)",           this.tsNoEdge);
+  			gd.addCheckbox    ("Only assign outside of 8x8 center if no suitable alternative",                    this.tsUseCenter);
   			gd.addNumericField("Maximal disparity difference when assigning tiles",                               this.tsMaxDiff,  6);
   			gd.addNumericField("Minimal disparity difference to be considered as a competitor surface",           this.tsMinDiffOther,  6);
   			gd.addNumericField("Minimal tile correlation strength to be assigned",                                this.tsMinStrength,  6);
@@ -3534,6 +3542,7 @@ public class EyesisCorrectionParameters {
   			this.plDbgMerge=            gd.getNextBoolean();
   			this.plWorstWorsening=      gd.getNextNumber();
   			this.plOKMergeEigen=        gd.getNextNumber();
+  			this.plMaxWorldSin2=        gd.getNextNumber();
   			this.plWeakWorsening=       gd.getNextNumber();
   			this.plMutualOnly=          gd.getNextBoolean();
 
@@ -3576,6 +3585,7 @@ public class EyesisCorrectionParameters {
   			this.msFractUni=            gd.getNextNumber();
 
   			this.tsNoEdge=              gd.getNextBoolean();
+  			this.tsUseCenter=           gd.getNextBoolean();
   			this.tsMaxDiff=             gd.getNextNumber();
   			this.tsMinDiffOther=        gd.getNextNumber();
   			this.tsMinStrength=         gd.getNextNumber();
@@ -3626,7 +3636,8 @@ public class EyesisCorrectionParameters {
   		
   		public boolean showTsDialog() {
   			GenericDialog gd = new GenericDialog("Set CLT tiles to surfaces assignment parameters");
-  			gd.addCheckbox    ("Do not assigned tiles to the surface edges (not having all 8 neighbors)",         this.tsNoEdge);
+  			gd.addCheckbox    ("Do not assign tiles to the surface edges (not having all 8 neighbors)",           this.tsNoEdge);
+  			gd.addCheckbox    ("Only assign outside of 8x8 center if no suitable alternative",                    this.tsUseCenter);
   			gd.addNumericField("Maximal disparity difference when assigning tiles",                               this.tsMaxDiff,  6);
   			gd.addNumericField("Minimal disparity difference to be considered as a competitor surface",           this.tsMinDiffOther,  6);
   			gd.addNumericField("Minimal tile correlation strength to be assigned",                                this.tsMinStrength,  6);
@@ -3659,6 +3670,7 @@ public class EyesisCorrectionParameters {
   			gd.showDialog();
   			if (gd.wasCanceled()) return false;
   			this.tsNoEdge=              gd.getNextBoolean();
+  			this.tsUseCenter=           gd.getNextBoolean();
   			this.tsMaxDiff=             gd.getNextNumber();
   			this.tsMinDiffOther=        gd.getNextNumber();
   			this.tsMinStrength=         gd.getNextNumber();

src/main/java/QuadCLT.java

@@ -5233,10 +5233,11 @@ public class QuadCLT {
         			  tp.clt_3d_passes.get(scanIndex), // CLTPass3d   scan,
         			  "X3D-"+scanIndex);
 			  }
-			  boolean showTri = ((scanIndex < next_pass + 1) && clt_parameters.show_triangles) ||(scanIndex == 49) || (scanIndex == 54);
+			  boolean showTri = ((scanIndex < next_pass + 1) && clt_parameters.show_triangles) ||((scanIndex - next_pass) == 73);
 			  generateClusterX3d(
 					  x3dOutput,
 					  texturePath,
+					  "shape_id-"+(scanIndex - next_pass), // id
 					  scan.bounds,
 					  scan.selected,
 					  scan_disparity, // scan.disparity_map[ImageDtt.DISPARITY_INDEX_CM],
@@ -5286,7 +5287,8 @@ public class QuadCLT {
 		  tp.showScan(
 				  tp.clt_3d_passes.get(next_pass-1),   // CLTPass3d   scan,
 				  "after_pass2-"+(next_pass-1)); //String title)
-		  tp.thirdPassSetup( // prepare tile tasks for the second pass based on the previous one(s)
+//		  tp.thirdPassSetup( // prepare tile tasks for the second pass based on the previous one(s)
+		  tp.thirdPassSetupSurf( // prepare tile tasks for the second pass based on the previous one(s)
 				  clt_parameters,
 				  clt_parameters.bgnd_range, // double            disparity_far, 
 				  clt_parameters.grow_disp_max, // other_range, //double            disparity_near,   //
@@ -5323,11 +5325,6 @@ public class QuadCLT {
 					  threadsMax,  // maximal number of threads to launch                         
 					  updateStatus,
 					  debugLevel);
-/**			  
-			  if ((scanIndex == 49) || (scanIndex == 54) ) tp.showScan(
-					  scan, // tp.clt_3d_passes.get(scanIndex), // CLTPass3d   scan,
-					  "MEASURED-"+scanIndex);
-*/					  
 		  }
 
 		  // TEMPORARY EXIT
@@ -5353,13 +5350,14 @@ public class QuadCLT {
 
 			  CLTPass3d scan = tp.clt_3d_passes.get(scanIndex);
 
+/*			  
 			  if ((scanIndex == 73) ) {
 				  tp.showScan(
 						  tp.clt_3d_passes.get(scanIndex), // CLTPass3d   scan,
 						  "SELECTED-"+scanIndex);
 			  }
 
-
+*/
 
 			  // TODO: use new updated disparity, for now just what was forced for the picture
 			  double [] scan_disparity = new double [tilesX * tilesY];
@@ -5383,17 +5381,21 @@ public class QuadCLT {
 					  scan_disparity[i] = sdw;
 				  }
 			  }
-
+/*
 			  if ((scanIndex == 73)) {
 				  tp.showScan(
 						  tp.clt_3d_passes.get(scanIndex), // CLTPass3d   scan,
 						  "X3D-"+scanIndex);
 			  }
-			  boolean showTri = ((scanIndex < next_pass + 1) && clt_parameters.show_triangles) ||(scanIndex == 73);
+*/			  
+//			  boolean showTri = ((scanIndex < next_pass + 1) && clt_parameters.show_triangles) ||(scanIndex < 3);
+			  boolean showTri = ((scanIndex < next_pass + 1) && clt_parameters.show_triangles) ||((scanIndex - next_pass) == 73);
+			  
 //			  boolean showTri = ((scanIndex < next_pass + 1) && clt_parameters.show_triangles) ||(scanIndex == 49) || (scanIndex == 54);
 			  generateClusterX3d(
 					  x3dOutput,
 					  texturePath,
+					  "shape_id-"+(scanIndex - next_pass), // id
 					  scan.bounds,
 					  scan.selected,
 					  scan_disparity, // scan.disparity_map[ImageDtt.DISPARITY_INDEX_CM],
@@ -5429,6 +5431,7 @@ public class QuadCLT {
 	  public void generateClusterX3d(
 			  X3dOutput  x3dOutput,
 			  String     texturePath,
+			  String     id,
 			  Rectangle  bounds,
 			  boolean [] selected,
 			  double []  disparity,
@@ -5486,6 +5489,7 @@ public class QuadCLT {
 
 		  x3dOutput.addCluster(
 				  texturePath,
+				  id,
 				  texCoord,
 				  worldXYZ,
 				  triangles);

src/main/java/SuperTiles.java

@@ -21,6 +21,7 @@
  ** -----------------------------------------------------------------------------**
  **
  */
+import java.awt.Point;
 import java.util.ArrayList;
 import java.util.Collections;
 import java.util.Comparator;
@@ -192,6 +193,19 @@ public class SuperTiles{
 	{
 		this.tileSurface = tileSurface;
 	}
+	public void setTileSurface(
+			GeometryCorrection   geometryCorrection)
+	{
+		this.tileSurface = new TileSurface(
+				tileProcessor.getTileSize(),      // int tileSize,
+				tileProcessor.getSuperTileSize(), // int superTileSize,
+				tileProcessor.getTilesX(),        // int tilesX,
+				tileProcessor.getTilesY(),        // int tilesY,
+				geometryCorrection,               // GeometryCorrection geometryCorrection,
+				tileProcessor.threadsMax);        // int threadsMax); 
+	}
+
+	
 	public TileSurface getTileSurface()
 	{
 		return tileSurface;
@@ -1510,10 +1524,11 @@ public class SuperTiles{
 	}
 	
 	public double [][] showSupertileSeparation(
+			boolean        useWorld,
 			double [][][]  disp_strength,
 			boolean [][][] selections)
 	{
-		return 	showSupertileSeparation(disp_strength, selections, null);
+		return 	showSupertileSeparation(useWorld, disp_strength, selections, null);
 	}
 
 	public String [] showSupertileSeparationTitles(
@@ -1548,6 +1563,7 @@ public class SuperTiles{
 	}
 
 	public double [][] showSupertileSeparation(
+			boolean useWorld,
 			double [][][]  disp_strength,
 			boolean [][][] selections,
 			TilePlanes.PlaneData [] planes
@@ -1588,13 +1604,17 @@ public class SuperTiles{
 		}
 		for (int npd = 0; npd < num_pd; npd++) if (planes[npd +LOWEST_PLANE(planes.length)] != null){
 			double [][] ellipsoids = planes[npd +LOWEST_PLANE(planes.length)].getDoublePlaneDisparityStrength(
+					useWorld,
 					null, // double [] window,
 					true, // boolean   use_sel,
 					true, // boolean   divide_by_area,
 					1.5, // double   scale_projection,
 					1); // int       debugLevel)
 			
-			data [num_pm + 3 * num_ml + 0 * num_pd + npd] = planes[npd +LOWEST_PLANE(planes.length)].getDoublePlaneDisparity(false);
+			data [num_pm + 3 * num_ml + 0 * num_pd + npd] = planes[npd +LOWEST_PLANE(planes.length)].getDoublePlaneDisparity(
+					useWorld,
+					false);
+			
 			boolean [][] msel = planes[npd +LOWEST_PLANE(planes.length)].getMeasSelection();
 			int ntiles = data [num_pm + 3 * num_ml + 0 * num_pd + npd].length; 
 			data [num_pm + 3 * num_ml + 1 * num_pd + npd] = new double [ntiles];
@@ -1616,6 +1636,7 @@ public class SuperTiles{
 	
 	// calculate "tilted" disparity, so planes parallel to the same world plane would have the same disparity
 	// also produces non-tilted, if world_plane_norm == null
+	// Protecting from behind the horizon - set strength of all tiles in the negative disparity area to 0
 	public double [][][][] getPlaneDispStrengths(
 			final double []  world_plane_norm, // real world normal vector to a suggested plane family (0,1,0) for horizontal planes
 			final int        stMeasSel, //            = 1;      // Select measurements for supertiles : +1 - combo, +2 - quad +4 - hor +8 - vert
@@ -1712,17 +1733,27 @@ public class SuperTiles{
 								if (dl > 0) {
 									System.out.println("Plane tilted disparity for stileX = "+stileX+" stileY="+stileY+", average disparity "+(sd/sw));
 								}
+								if (dl>1) {
+									String [] dbg_titles = showSupertileSeparationTitles( plane_disp_strength[nsTile], null);
+									double [][] dbg_img = showSupertileSeparation(false, plane_disp_strength[nsTile], null); // plane_sels);
+									showDoubleFloatArrays sdfa_instance = new showDoubleFloatArrays();
+									sdfa_instance.showArrays(dbg_img, 2 * superTileSize, 2* superTileSize, true, "plane_separation_original_disp"+nsTile,dbg_titles);
+									dbg_img = showSupertileSeparation(true, plane_disp_strength[nsTile], null); // plane_sels);
+									sdfa_instance.showArrays(dbg_img, 2 * superTileSize, 2* superTileSize, true, "plane_separation_original_world"+nsTile,dbg_titles);
+								}
 								plane_disp_strength[nsTile] = pd0.getDisparityToPlane(
 										world_plane_norm, // double []     world_normal_xyz,
 										sd / sw, // average disparity // double        disp_center,
 										null, // boolean [][]  tile_sel, // null - do not use, {} use all (will be modified)
 										plane_disp_strength[nsTile], // double [][][] disp_str, // calculate just once if null 
-										1); // int           debugLevel);
+										dl); // 1); // int           debugLevel);
 								if (dl>1) {
 									String [] dbg_titles = showSupertileSeparationTitles( plane_disp_strength[nsTile], null);
-									double [][] dbg_img = showSupertileSeparation(plane_disp_strength[nsTile], null); // plane_sels);
+									double [][] dbg_img = showSupertileSeparation(false, plane_disp_strength[nsTile], null); // plane_sels);
 									showDoubleFloatArrays sdfa_instance = new showDoubleFloatArrays();
-									sdfa_instance.showArrays(dbg_img, 2 * superTileSize, 2* superTileSize, true, "plane_separation"+nsTile,dbg_titles);
+									sdfa_instance.showArrays(dbg_img, 2 * superTileSize, 2* superTileSize, true, "plane_separation_disp"+nsTile,dbg_titles);
+									dbg_img = showSupertileSeparation(true, plane_disp_strength[nsTile], null); // plane_sels);
+									sdfa_instance.showArrays(dbg_img, 2 * superTileSize, 2* superTileSize, true, "plane_separation_world"+nsTile,dbg_titles);
 								}
 							} else {
 								plane_disp_strength[nsTile] = null;
@@ -1988,9 +2019,11 @@ public class SuperTiles{
 
 							if (dl > 3) {
 								String [] dbg_titles = showSupertileSeparationTitles( disp_strength, plane_sels);
-								double [][] dbg_img = showSupertileSeparation(disp_strength, plane_sels);
+								double [][] dbg_img = showSupertileSeparation(false,disp_strength, plane_sels);
 								showDoubleFloatArrays sdfa_instance = new showDoubleFloatArrays();
-								sdfa_instance.showArrays(dbg_img, 2 * superTileSize, 2* superTileSize, true, "disp_clusterize"+nsTile+"-"+debugLevel,dbg_titles);
+								sdfa_instance.showArrays(dbg_img, 2 * superTileSize, 2* superTileSize, true, "disp_clusterize_disp"+nsTile+"-"+debugLevel,dbg_titles);
+								dbg_img = showSupertileSeparation(true, disp_strength, plane_sels);
+								sdfa_instance.showArrays(dbg_img, 2 * superTileSize, 2* superTileSize, true, "disp_clusterize_world"+nsTile+"-"+debugLevel,dbg_titles);
 							}
 							plane_selections[nsTile] = plane_sels;
 						}
@@ -2302,6 +2335,7 @@ public class SuperTiles{
 			final double     plTargetEigen, //        =   0.1;  // Remove outliers until main axis eigenvalue (possibly scaled by plDispNorm) gets below
 			final double     plFractOutliers, //      =   0.3;  // Maximal fraction of outliers to remove
 			final int        plMaxOutliers, //        =    20;  // Maximal number of outliers to remove
+//			final double     plVertWors,    //        =    1.5  // if rotating plane vertical does not increase 'eigenvalue' more, use vertical  
 			final boolean    plPreferDisparity, // Always start with disparity-most axis (false - lowest eigenvalue)
 			final GeometryCorrection geometryCorrection,
 			final boolean    correct_distortions,
@@ -2445,9 +2479,11 @@ public class SuperTiles{
 							}
 							if (dl > 2) {
 								String [] dbg_titles = showSupertileSeparationTitles( disp_strength[nsTile], plane_selections[nsTile], result_planes[nsTile]);
-								double [][] dbg_img =  showSupertileSeparation(disp_strength[nsTile], plane_selections[nsTile], result_planes[nsTile]);
+								double [][] dbg_img =  showSupertileSeparation(false, disp_strength[nsTile], plane_selections[nsTile], result_planes[nsTile]);
 								showDoubleFloatArrays sdfa_instance = new showDoubleFloatArrays();
-								sdfa_instance.showArrays(dbg_img, 2 * superTileSize, 2* superTileSize, true, "create_planes-"+nsTile+"-"+debugLevel,dbg_titles);
+								sdfa_instance.showArrays(dbg_img, 2 * superTileSize, 2* superTileSize, true, "create_planes_disp-"+nsTile+"-"+debugLevel,dbg_titles);
+								dbg_img =  showSupertileSeparation(true, disp_strength[nsTile], plane_selections[nsTile], result_planes[nsTile]);
+								sdfa_instance.showArrays(dbg_img, 2 * superTileSize, 2* superTileSize, true, "create_planes_world-"+nsTile+"-"+debugLevel,dbg_titles);
 							}
 						}
 					}
@@ -2528,7 +2564,7 @@ public class SuperTiles{
 
 		// get per-tile disparity strength again (may consider using non-filtered data here)
 		
-		double [][][][] disp_strength = getPlaneDispStrengths(
+		double [][][][] disp_strength = getPlaneDispStrengths( // here use actual disparity, not tilted
 				null,                // final double []  world_plane_norm, // real world normal vector to a suggested plane family (0,1,0) for horizontal planes
 				stMeasSel,           // final int        stMeasSel, //            = 1;      // Select measurements for supertiles : +1 - combo, +2 - quad +4 - hor +8 - vert
 
@@ -2568,25 +2604,6 @@ public class SuperTiles{
 				dbg_X,               // final int        dbg_X,
 				dbg_Y);              // final int        dbg_Y)
 		this.planes = new_planes; // save as "measured" (as opposed to "smoothed" by neighbors) planes
-/*		
-		TileSurface tileSurface = new TileSurface(
-				tileProcessor.getTileSize(),      // int tileSize,
-				tileProcessor.getSuperTileSize(), // int superTileSize,
-				tileProcessor.getTilesX(),        // int tilesX,
-				tileProcessor.getTilesY(),        // int tilesY,
-				geometryCorrection,               // GeometryCorrection geometryCorrection,
-				tileProcessor.threadsMax);        // int threadsMax); 
-
-		
-		TileSurface.TileData [][] tileData = tileSurface.createTileShells (
-				msUseSel,            // final boolean                   use_sel,
-				msDivideByArea,      // final boolean                   divide_by_area,
-				msScaleProj,         // final double                    scale_projection,
-				this.planes,         // final TilePlanes.PlaneData [][] planes,
-				debugLevel,          // final int        debugLevel,
-				dbg_X,               // final int        dbg_X,
-				dbg_Y);              // final int        dbg_Y)
-*/		
 		
 	}
 	
@@ -2751,9 +2768,11 @@ public class SuperTiles{
 						// plane_disp_strength
 						if (dl > 2) {
 							String [] dbg_titles = showSupertileSeparationTitles(plane_disp_strength[nsTile], null);
-							double [][] dbg_img = showSupertileSeparation(plane_disp_strength[nsTile], null);
+							double [][] dbg_img = showSupertileSeparation(false,plane_disp_strength[nsTile], null);
 							showDoubleFloatArrays sdfa_instance = new showDoubleFloatArrays();
-							sdfa_instance.showArrays(dbg_img, 2 * superTileSize, 2* superTileSize, true, "HOR_SEP_"+nsTile,dbg_titles);
+							sdfa_instance.showArrays(dbg_img, 2 * superTileSize, 2* superTileSize, true, "HOR_SEP_DISP"+nsTile,dbg_titles);
+							dbg_img = showSupertileSeparation(true, plane_disp_strength[nsTile], null);
+							sdfa_instance.showArrays(dbg_img, 2 * superTileSize, 2* superTileSize, true, "HOR_SEP_WORLD"+nsTile,dbg_titles);
 						}
 						
 						ArrayList<TilePlanes.PlaneData> st_planes = new ArrayList<TilePlanes.PlaneData>();
@@ -2945,9 +2964,11 @@ public class SuperTiles{
 						
 						if (dl > 2) {
 							String [] dbg_titles = showSupertileSeparationTitles( disp_strength, plane_sels);
-							double [][] dbg_img = showSupertileSeparation(disp_strength, plane_sels);
+							double [][] dbg_img = showSupertileSeparation(false,disp_strength, plane_sels);
 							showDoubleFloatArrays sdfa_instance = new showDoubleFloatArrays();
-							sdfa_instance.showArrays(dbg_img, 2 * superTileSize, 2* superTileSize, true, "initial_separation"+nsTile,dbg_titles);
+							sdfa_instance.showArrays(dbg_img, 2 * superTileSize, 2* superTileSize, true, "initial_separation_disp"+nsTile,dbg_titles);
+							dbg_img = showSupertileSeparation(true, disp_strength, plane_sels);
+							sdfa_instance.showArrays(dbg_img, 2 * superTileSize, 2* superTileSize, true, "initial_separation_world"+nsTile,dbg_titles);
 						}
 
 						if (dl > 2) {
@@ -2966,9 +2987,11 @@ public class SuperTiles{
 										disp_strength, // double [][][] disp_str, // calculate just once if null 
 										1); // int           debugLevel);
 								String [] dbg_titles = showSupertileSeparationTitles( hor_disp_strength, plane_sels);
-								double [][] dbg_img = showSupertileSeparation(hor_disp_strength, plane_sels);
+								double [][] dbg_img = showSupertileSeparation(false, hor_disp_strength, plane_sels);
 								showDoubleFloatArrays sdfa_instance = new showDoubleFloatArrays();
-								sdfa_instance.showArrays(dbg_img, 2 * superTileSize, 2* superTileSize, true, "hor_separation"+nsTile,dbg_titles);
+								sdfa_instance.showArrays(dbg_img, 2 * superTileSize, 2* superTileSize, true, "hor_separation_disp"+nsTile,dbg_titles);
+								dbg_img = showSupertileSeparation(true, hor_disp_strength, plane_sels);
+								sdfa_instance.showArrays(dbg_img, 2 * superTileSize, 2* superTileSize, true, "hor_separation_world"+nsTile,dbg_titles);
 							}
 						}
 						
@@ -3656,12 +3679,305 @@ public class SuperTiles{
 		return num_added;
 
 	}
+	class Conflict{
+		int start_layer;
+		int end_layer;
+		boolean [] start_dirs = new boolean[8]; // 4 last - opposite direction
+		Conflict(int start_layer, int end_layer, int start_dir)
+		{
+			if (end_layer > start_layer) {
+				this.start_layer =               start_layer;
+				this.end_layer =                 end_layer;
+				this.start_dirs[start_dir] = true;
+			} else {
+				this.start_layer =               end_layer;
+				this.end_layer =                 start_layer;
+				this.start_dirs[start_dir + 4] = true;
+
+			}
+		}
+		boolean combine (Conflict other_conflict)
+		{
+			if ((other_conflict.start_layer == this.start_layer) && (other_conflict.end_layer == this.end_layer)) {
+				for (int i = 0; i < start_dirs.length; i++) start_dirs[i] |= other_conflict.start_dirs[i];
+				return true;
+			}
+			return false;
+		}
+		int [] toArray()
+		{
+			int dirs_bits = 0;
+			for (int i = 0; i < start_dirs.length; i++){
+				if (start_dirs[i]) dirs_bits |= (1 << i);
+			}
+			int [] rslt = {start_layer, end_layer, dirs_bits};
+			return rslt;
+		}
+	}
+	
+	
+	public int [][] detectTriangularTileConflicts(
+			int nsTile0,
+			TileSurface.TileNeibs tnSurface)
+	{
+		ArrayList<Conflict> conflicts_list= new ArrayList<Conflict>();
+		if ( planes[nsTile0] != null) {
+			for (int np0 = 0; np0 < planes[nsTile0].length; np0++){
+				if (planes[nsTile0][np0] != null) {
+					for (int dir = 0; dir < 8; dir +=2){
+						int np1 = planes[nsTile0][np0].getNeibBest(dir);
+						if (np1 >=0) {
+							int nsTile1 = tnSurface.getNeibIndex(nsTile0, dir);
+							if (nsTile1 >= 0){
+								int dir1 = (dir + 3) % 8;
+								int np2 = planes[nsTile1][np1].getNeibBest(dir1);
+								if (np2 >= 0){
+									int nsTile2 = tnSurface.getNeibIndex(nsTile1, dir1);
+									if (nsTile2 >= 0){
+										int dir2 = (dir1 + 3) % 8;
+										int np3 = planes[nsTile2][np2].getNeibBest(dir2);
+										if ((np3 >= 0) && (np3 != np0)){
+											int cls = conflicts_list.size();
+											Conflict conflict = new Conflict(np0, np3, dir/2);
+											if ((cls == 0) || !conflicts_list.get(cls - 1).combine(conflict)) {
+												conflicts_list.add(conflict);
+											}
+										}
+									}
+								}
+							}
+						}
+					}
+				}
+			}
+		}
+		if(conflicts_list.isEmpty()){
+			return null;
+		}
+		int [][] conflicts = new int [conflicts_list.size()][];
+		int indx=0;
+		for (Conflict conflict:conflicts_list){
+			conflicts[indx++] = conflict.toArray();
+		}
+		return conflicts;
+	}
+	/**
+	 * Find "triangular" conflicts after running selectNeighborPlanesMutual.
+	 * Such conflicts happen is when starting N (or other ortho direction) from some node,
+	 * then turning 135 degrees right, and then 135 right again it will get to the different
+	 * layer than started (and all 3 connections exist. 
+	 * @param debugLevel
+	 * @return 3-d array, first index - tile number, 2-nd - tile conflict number. Each conflict
+	 * is stored as {start layer, end layer, bitmask of start directions} (+1 - N, +2 - E, +4 - S, +8 - W)
+	 */
+
+	
+	public int [][][] detectTriangularConflicts(
+			final int debugLevel)
+	{
+		final int tilesX =        tileProcessor.getTilesX();
+		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 nStiles =       stilesX * stilesY;
+		final int [][][] conflicts = new int [planes.length][][];
+		final TileSurface.TileNeibs tnSurface = tileSurface.new TileNeibs(stilesX, stilesY);
+		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 nsTile0 = ai.getAndIncrement(); nsTile0 < nStiles; nsTile0 = ai.getAndIncrement()) {
+						if ( planes[nsTile0] != null) {
+							conflicts[nsTile0] =  detectTriangularTileConflicts(
+									nsTile0,
+									tnSurface);
+						}
+					}
+				}
+			};
+		}
+		ImageDtt.startAndJoin(threads);
+		if (debugLevel > -1){
+			int [] num_conflicts = new int [8];
+			for (int nsTile = 0; nsTile < conflicts.length; nsTile++) if (conflicts[nsTile] != null){
+				for (int nc = 0; nc < conflicts[nsTile].length; nc++){
+					int np = 0;
+					for (int i = 0; i < 8; i++){
+						if ((conflicts[nsTile][nc][2] & (1 << i)) != 0) {
+							np ++;
+						}
+					}
+					num_conflicts[np]++;
+					if (debugLevel > 0) {
+						System.out.println("detectTriangularConflicts() nsTile ="+nsTile+
+								" nl1 = "+ conflicts[nsTile][nc][0]+
+								" nl2 = "+ conflicts[nsTile][nc][1]+
+								" bits = "+ conflicts[nsTile][nc][2]);
+					}
+				}
+			}
+			System.out.print("detectTriangularConflicts() detected conflicts:");
+			for (int i = 0; i < 8; i++){
+				if (num_conflicts[i] >0) System.out.print(" type "+i+": "+num_conflicts[i]);
+			}
+			System.out.println();
+		}
+		return conflicts;
+	}
+	
+	/**
+	 * Find out which of the layers to use to resolve the conflict. Multiple conflict may need to
+	 * be pre-filtered or split to prevent using simultaneously different neighbors for the same
+	 * direction
+	 * @param nsTile supertile index
+	 * @param nl1 first conflicting layer (lower of the too
+	 * @param nl2 second conflicting layer
+	 * @param dir_mask direction bit mask: bit 0: nl1->N->SE->W->nl2, bit 1: nl1->E->SW->N->nl2,...,
+	 * bit 4 : nl1->E->NW->S->nl2 (backwards of bit1), bit 5: nl1->S->NE->W->nl2 (backwards of bit1)
+	 *  	 * @param rquality maximal degradation by merging (does not depend on the total weight)
+	 * @param okMergeEigen if result eigenvalue of the merged planes is below, OK to bypass worst worsening
+	 * @param maxWorldSin2 maximal square of the sine of the angle between the planes to allow merge (>= 1.0 - disable)
+	 * @param maxEigen maximal eigenvalue of each of the merged planes
+	 * @param tnSurface TileNeibs instance to navigate tile index and control array borders
+	 * @param preferDisparity - the first eigenvalue/vector is the most disparity-like
+	 *                          (false - smallest eigenvalue)
+	 * @param debugLevel
+	 * @return winning layer or -1 if not possible
+	 * TODO: add more parameters - some thresholds 
+	 */
+	public int resolveTriangularConflict(
+			int nsTile,
+			int nl1,
+			int nl2,
+			int dir_mask,
+			double rquality,
+			double weakWorsening,
+			double okMergeEigen,
+			double maxWorldSin2,
+			double dispNorm,
+			double maxEigen, // maximal eigenvalue of planes to consider
+			TileSurface.TileNeibs tnSurface,
+			boolean preferDisparity,
+			int debugLevel)
+	{
+		ArrayList<Point> surf_list = new ArrayList<Point>();
+		for (int dir = 0; dir < 8; dir +=2){
+			if ((dir_mask & (1 << (dir / 2))) != 0){
+				int nsTile1 = tnSurface.getNeibIndex(nsTile, dir);
+				int nsTile2 = tnSurface.getNeibIndex(nsTile, (dir + 2) % 8);
+				surf_list.add(new Point(nsTile1,nl1));
+				surf_list.add(new Point(nsTile2,nl2));
+			}
+			if ((dir_mask & (1 << ((dir / 2) + 4))) != 0){
+				int nsTile1 = tnSurface.getNeibIndex(nsTile, dir);
+				int nsTile2 = tnSurface.getNeibIndex(nsTile, (dir + 2) % 8);
+				surf_list.add(new Point(nsTile1,nl2));
+				surf_list.add(new Point(nsTile2,nl1));
+			}
+		}
+		TilePlanes.PlaneData merged_plane = null;
+		TilePlanes.PlaneData this_plane1 = planes[nsTile][nl1];
+		TilePlanes.PlaneData this_plane2 = planes[nsTile][nl2];
+		for (Point p : surf_list) {
+			TilePlanes.PlaneData other_plane = this_plane1.getPlaneToThis(  // layer here does not matter
+					planes[p.x][p.y],
+					debugLevel - 1); // debugLevel);
+			if (merged_plane == null){
+				merged_plane = other_plane; 
+			} else {
+				merged_plane = merged_plane.mergePlaneToThis(
+						other_plane, // PlaneData otherPd,
+						1.0,         // double    scale_other,
+						false,       // boolean   ignore_weights,
+						true, // boolean   sum_weights,
+						preferDisparity, 
+						debugLevel - 1); // int       debugLevel)
+			}
+		}
+		TilePlanes.PlaneData [] these = {this_plane1, this_plane2};
+		TilePlanes.PlaneData [] merged = new TilePlanes.PlaneData [2];
+
+		// now find which is better and if any of them fits at all. Print other criteria also - maybe use them later
+		double [] merge_ev = new double[2];
+		boolean [] fit = {true,true}; // initially both fit
+		double [] this_rq = new double[2];
+		for (int np = 0; np < 2; np++){
+			merged[np] = these[np].mergePlaneToThis(
+					merged_plane, // PlaneData otherPd,
+					1.0,          // double    scale_other,
+					false,        // boolean   ignore_weights,
+					true,         // boolean   sum_weights,
+					preferDisparity, 
+					debugLevel - 1); // int       debugLevel)
+
+			// verify result plane smallest eigenvalue is OK  
+			if ((maxEigen != 0.0) &&
+					(merged[np].getValue() > corrMaxEigen(
+							maxEigen,
+							dispNorm,
+							merged[np]))) {
+				fit[np] = false;
+				continue;
+			}
+			double w1 = merged_plane.getWeight();
+			double w2 = merged[np].getWeight();
+			this_rq[np] = mergeRQuality(
+					merged_plane.getValue(), // double L1,
+					merged[np].getValue(), // double L2,
+					merged[np].getValue(), // double L,
+					w1, // double w1,
+					w2); // double w2)
+			double this_rq_norm = this_rq[np];
+			if ((w1 + w2) < weakWorsening) this_rq_norm *= (w1 + w2) / weakWorsening; // forgive more for weak planes
+
+			// verify that eigenvalue worsening by the merge is OK
+			if (    (this_rq_norm > rquality) &&
+					(merged[np].getValue() > okMergeEigen)){
+				fit[np] = false;
+				continue;
+			}
+			// verify
+			if ((maxWorldSin2 < 1.0) && (merged_plane.getWorldSin2(merged[np]) <=maxWorldSin2)) {
+				fit[np] = false;
+				continue;
+			}
+			this_rq[np] /= (w1 + w2); // for comparison reduce this value for stronger planes 
+			if (debugLevel >0) {
+				System.out.println("nsTile="+nsTile+":"+np+", nsTile="+nsTile+":"+np+", this_rq="+this_rq[np]+
+						" w1="+w1+" w2="+w2+
+						" L1="+merged_plane.getValue()+" L2="+these[np].getValue()+" L="+merged[np].getValue());
+
+				System.out.println("nsTile="+nsTile+":"+np+", world sin2 ="+
+						merged_plane.getWorldSin2(these[np]));
+
+				System.out.println("nsTile="+nsTile+":"+np+
+						", world dist this="+ Math.sqrt(merged_plane.getWorldPlaneDist2(these[np]))+
+						", world dist other="+Math.sqrt(these[np].getWorldPlaneDist2(merged_plane))+
+						", world dist sum="+Math.sqrt(merged_plane.getWorldPlaneDist2(these[np])+
+								these[np].getWorldPlaneDist2(merged_plane)));
+			}
+		}
+		if (fit[0] || fit[1]){
+			if (!fit[0]) return nl2;
+			if (!fit[1]) return nl1;
+			
+			if (this_rq[0] < this_rq[1]) return nl1;
+			else return nl2;
+		}
+
+		return -1;
+	}
+	
+	
 	
 	/**
 	 * Find mutual links between multi-layer planes for supertiles. requires that for each plane there are calculated smalles eigenvalues
 	 * for merging with each plane for each of 8 neighbors 
 	 * @param rquality maximal degradation by merging (does not depend on the total weight)
 	 * @param okMergeEigen if result eigenvalue of the merged planes is below, OK to bypass worst worsening
+	 * @param maxWorldSin2 maximal square of the sine of the angle between the planes to allow merge (>= 1.0 - disable)
 	 * @param maxEigen maximal eigenvalue of each of the merged planes
 	 * @param minWeight minimal weight of each of the planes
 	 * @param debugLevel debug level
@@ -3672,6 +3988,7 @@ public class SuperTiles{
 			final double rquality,
 			final double weakWorsening,
 			final double okMergeEigen,
+			final double maxWorldSin2,
 			final double dispNorm,
 			final double maxEigen, // maximal eigenvalue of planes to consider
 			final double minWeight, // minimal pain weight to consider
@@ -3781,12 +4098,14 @@ public class SuperTiles{
 																w2); // double w2)
 														double this_rq_norm = this_rq;
 														if ((w1 + w2) < weakWorsening) this_rq_norm *= (w1 + w2) / weakWorsening; // forgive more for weak planes 
-														if ((this_rq_norm <= rquality) ||(merge_ev[np] <= okMergeEigen)) {   
-															this_rq /= (w1 + w2); // for comparision reduce this value for stronger planes 
-															if (Double.isNaN(best_rqual) || (this_rq < best_rqual)){ // OK if Double.isNaN(this_rq[np])
-																best_rqual = this_rq;
-																best_pair[0]= np0;
-																best_pair[1]= np;
+														if ((this_rq_norm <= rquality) ||(merge_ev[np] <= okMergeEigen)) {
+															if ((maxWorldSin2 >= 1.0) || (planes[nsTile0][np0].getWorldSin2(planes[nsTile][np]) <=maxWorldSin2)) {
+																this_rq /= (w1 + w2); // for comparison reduce this value for stronger planes 
+																if (Double.isNaN(best_rqual) || (this_rq < best_rqual)){ // OK if Double.isNaN(this_rq[np])
+																	best_rqual = this_rq;
+																	best_pair[0]= np0;
+																	best_pair[1]= np;
+																}
 															}
 														}
 														if (debugLevel > 0){
@@ -3801,6 +4120,13 @@ public class SuperTiles{
 															System.out.println("nsTile0="+nsTile0+":"+np0+", nsTile="+nsTile+":"+np+", this_rq="+this_rq+
 																	" w1="+w1+" w2="+w2+
 																	" L1="+planes[nsTile0][np0].getValue()+" L2="+planes[nsTile][np].getValue()+" L="+merge_ev[np]);
+															System.out.println("nsTile0="+nsTile0+":"+np0+", nsTile="+nsTile+":"+np+", world sin2 ="+
+																	planes[nsTile0][np0].getWorldSin2(planes[nsTile][np]));
+															System.out.println("nsTile0="+nsTile0+":"+np0+", nsTile="+nsTile+":"+np+
+																	", world dist this="+ Math.sqrt(planes[nsTile0][np0].getWorldPlaneDist2(planes[nsTile][np]))+
+																	", world dist other="+Math.sqrt(planes[nsTile][np].getWorldPlaneDist2(planes[nsTile0][np0]))+
+																	", world dist sum="+Math.sqrt(planes[nsTile0][np0].getWorldPlaneDist2(planes[nsTile][np])+
+																			planes[nsTile][np].getWorldPlaneDist2(planes[nsTile0][np0])));
 														}
 													}															
 												}														
@@ -3822,7 +4148,8 @@ public class SuperTiles{
 											planes[nsTile0][np].setNeibBest(dir,-1);
 										}
 									}
-									for (int np = 0; np < other_matched.length; np++) if (planes[nsTile][np] != null){
+//									for (int np = 0; np < other_matched.length; np++) if (planes[nsTile][np] != null){
+									for (int np = 0; np < other_matched.length; np++) if ((planes[nsTile][np] != null) && (planes[nsTile][np].getWeight() > 0.0)){ // disregard 0-weight planes
 										if (!other_matched[np]){
 											planes[nsTile][np].setNeibBest(dir + 4,-1);
 										}
@@ -3833,9 +4160,10 @@ public class SuperTiles{
 											System.out.println("nsTile0="+nsTile0+":"+np+" : ["+
 											bn[0]+","+bn[1]+","+bn[2]+","+bn[3]+","+bn[4]+","+bn[5]+","+bn[6]+","+bn[7]+"]");
 										}
-										for (int np = 0; np < other_matched.length; np++) if (planes[nsTile][np] != null){
+//										for (int np = 0; np < other_matched.length; np++) if (planes[nsTile][np] != null){
+										for (int np = 0; np < other_matched.length; np++) if ((planes[nsTile][np] != null) && (planes[nsTile][np].getWeight() > 0.0)){ // disregard 0-weight planes
 											int [] bn = planes[nsTile][np].getNeibBest();
-											System.out.println("nsTile="+nsTile+":"+np+" : ["+
+											System.out.println("nsTile="+nsTile+":"+np+" best neighbors : ["+
 											bn[0]+","+bn[1]+","+bn[2]+","+bn[3]+","+bn[4]+","+bn[5]+","+bn[6]+","+bn[7]+"]");
 										}
 									}
@@ -3893,7 +4221,7 @@ public class SuperTiles{
 				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];
+						double disp = planes[i][j].getSinglePlaneDisparity(false)[centerIndex];
 						/*
 								if (disp > dispNorm) {
 									eigVal *= dispNorm / disp;
@@ -3924,7 +4252,9 @@ 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++) if (planes[nsTile][np] != null){
+//					for (int np = 0; np < planes[nsTile].length; np++) if (planes[nsTile][np] != null){
+					for (int np = 0; np < planes[nsTile].length; np++) if ((planes[nsTile][np] != null) && (planes[nsTile][np].getWeight() > 0.0)){ // disregard 0-weight planes
+					
 						if (planes[nsTile][np].getWeight() >= minWeight){
 							//&& (planes[nsTile][np].getValue() < maxEigen)){
 							double eigVal = planes[nsTile][np].getValue();
@@ -3938,7 +4268,7 @@ public class SuperTiles{
 									maxEigen,
 									dispNorm,
 									disp)) {
-								plane = planes[nsTile][np].getPlaneDisparity(use_NaN);
+								plane = planes[nsTile][np].getSinglePlaneDisparity(use_NaN);
 								// add connections to neighbors if available
 								int [] neib_best = planes[nsTile][np].getNeibBest();
 								if (neib_best != null){
@@ -4048,7 +4378,7 @@ public class SuperTiles{
 						for (int np = 0; np < planes[nsTile].length; np++){
 							if ((shells[nsTile][np] == (nShell + 1)) && ((nlayer == 0) || (num_rendered < nlayer)) ) {
 								num_rendered ++;
-								plane = planes[nsTile][np].getPlaneDisparity(use_NaN);
+								plane = planes[nsTile][np].getSinglePlaneDisparity(use_NaN);
 								int [] neib_best = planes[nsTile][np].getNeibBest();
 								if (fuse && (neib_best != null)){
 									if (debug) {
@@ -4057,7 +4387,7 @@ public class SuperTiles{
 										dbg_this = new double [dbg_titles.length][];
 										dbg_other = new double [dbg_titles.length][];
 										dbg_img3[0] = planes[nsTile][np].getTriplePlaneDisparity();
-										dbg_img[0] = planes[nsTile][np].getPlaneDisparity(false);
+										dbg_img[0] = planes[nsTile][np].getSinglePlaneDisparity(false);
 									}
 
 
@@ -4092,7 +4422,7 @@ public class SuperTiles{
 										sdfa_instance.showArrays(dbg_this, superTileSize, superTileSize, true, "weights_this_"+nsTile,dbg_titles);
 										sdfa_instance.showArrays(dbg_other, superTileSize, superTileSize, true, "weights_other_"+nsTile,dbg_titles);
 									}
-									double [] plane_other = planes[nsTile][np].getPlaneDisparity(false);
+									double [] plane_other = planes[nsTile][np].getSinglePlaneDisparity(false);
 									for (int i = 0; i < weight_this.length; i++){
 										plane_combo[i] += weight_this[i] * plane_other[i];
 										double w = weight_this[i] + weight_other[i];
@@ -4214,7 +4544,7 @@ public class SuperTiles{
 						for (int np = 0; np < planes[nsTile].length; np++){
 							if ((shells[nsTile][np] == (nShell + 1)) && ((nlayer == 0) || (num_rendered < nlayer)) ) {
 								num_rendered ++;
-								plane = planes[nsTile][np].getPlaneDisparity(use_NaN);
+								plane = planes[nsTile][np].getSinglePlaneDisparity(use_NaN);
 								int [] neib_best = planes[nsTile][np].getNeibBest();
 								if (fuse && (neib_best != null)){
 									if (debug) {
@@ -4223,7 +4553,7 @@ public class SuperTiles{
 										dbg_this = new double [dbg_titles.length][];
 										dbg_other = new double [dbg_titles.length][];
 										dbg_img3[0] = planes[nsTile][np].getTriplePlaneDisparity();
-										dbg_img[0] = planes[nsTile][np].getPlaneDisparity(false);
+										dbg_img[0] = planes[nsTile][np].getSinglePlaneDisparity(false);
 									}
 
 
@@ -4258,7 +4588,7 @@ public class SuperTiles{
 										sdfa_instance.showArrays(dbg_this, superTileSize, superTileSize, true, "weights_this_"+nsTile,dbg_titles);
 										sdfa_instance.showArrays(dbg_other, superTileSize, superTileSize, true, "weights_other_"+nsTile,dbg_titles);
 									}
-									double [] plane_other = planes[nsTile][np].getPlaneDisparity(false);
+									double [] plane_other = planes[nsTile][np].getSinglePlaneDisparity(false);
 									for (int i = 0; i < weight_this.length; i++){
 										plane_combo[i] += weight_this[i] * plane_other[i];
 										double w = weight_this[i] + weight_other[i];
@@ -4437,11 +4767,11 @@ public class SuperTiles{
 									System.out.println("Bug?  new_planes["+nsTile0+"]["+np0+"] == null");
 									continue;
 								}
-								if (dl > 0) dbg_img[ 0] = this_new_plane.getPlaneDisparity(false);
-								if (dl > 0) dbg_img[ 1] = measured_planes[nsTile0][np0].getPlaneDisparity(false);
+								if (dl > 0) dbg_img[ 0] = this_new_plane.getSinglePlaneDisparity(false);
+								if (dl > 0) dbg_img[ 1] = measured_planes[nsTile0][np0].getSinglePlaneDisparity(false);
 
 								int [] neibs = this_new_plane.getNeibBest();
-								this_new_plane.setWeight(0.0);
+								this_new_plane.setWeight(0.0); //
 								double num_merged = 0.0; // double to add fractional pull weight of the center
 								for (int dir = 0; dir < neibs.length; dir++){
 									if (neibs[dir] >= 0) {
@@ -4451,7 +4781,7 @@ public class SuperTiles{
 										TilePlanes.PlaneData other_plane = this_new_plane.getPlaneToThis(
 												mod_planes[nsTile][neibs[dir]], // neighbor, previous value
 												dl - 1); // debugLevel);
-										if (dl > 0) dbg_img[ 2 + dir] = other_plane.getPlaneDisparity(false);
+										if (dl > 0) dbg_img[ 2 + dir] = other_plane.getSinglePlaneDisparity(false);
 										if ((other_plane != null) && ((other_plane.getValue() <= maxValue) || (maxValue == 0))) {
 											if (this_new_plane.getWeight() > 0.0){
 												this_new_plane = this_new_plane.mergePlaneToThis(
@@ -4469,7 +4799,7 @@ public class SuperTiles{
 												this_new_plane = other_plane; 
 											}
 											new_planes[nsTile0][np0] = this_new_plane;
-											if (dl > 0) dbg_img[10 + dir] = this_new_plane.getPlaneDisparity(false);
+											if (dl > 0) dbg_img[10 + dir] = this_new_plane.getSinglePlaneDisparity(false);
 										} else if (stopBadNeib){
 											this_new_plane = null;
 											break;
@@ -4489,7 +4819,9 @@ public class SuperTiles{
 											(measured_planes != null) &&
 											(measured_planes[nsTile0] != null) &&
 											(measured_planes[nsTile0][np0] != null) &&
-											((measured_planes[nsTile0][np0].getValue() < maxValue) || (maxValue == 0))
+//											((measured_planes[nsTile0][np0].getValue() < maxValue) || (maxValue == 0))
+											((measured_planes[nsTile0][np0].getValue() < maxValue) || (maxValue == 0) || 
+													(this_new_plane == null) || (this_new_plane.getWeight() == 0.0)) // keep measured if impossible to merge
 											){ // merge with "measured"
 
 										if (this_new_plane.getWeight() > 0.0){
@@ -4508,7 +4840,7 @@ public class SuperTiles{
 											num_merged = 1.0;	
 										}
 										new_planes[nsTile0][np0] = this_new_plane;
-										if (dl > 0) dbg_img[18] =  this_new_plane.getPlaneDisparity(false);
+										if (dl > 0) dbg_img[18] =  this_new_plane.getSinglePlaneDisparity(false);
 									}
 									if ((num_merged > 0.0) && (this_new_plane != null)){
 										this_new_plane.scaleWeight(1.0/num_merged);
@@ -4522,9 +4854,9 @@ public class SuperTiles{
 									// calculate largest disparity difference between old and new plane
 									if (rslt_diffs != null){ // filter out outliers here?
 										// get plane for both old and new, calc rms of diff
-										double [] oldPlane = mod_planes[nsTile0][np0].getPlaneDisparity(
+										double [] oldPlane = mod_planes[nsTile0][np0].getSinglePlaneDisparity(
 												false); // use_NaN)
-										double [] newPlane = new_planes[nsTile0][np0].getPlaneDisparity(
+										double [] newPlane = new_planes[nsTile0][np0].getSinglePlaneDisparity(
 												false); // use_NaN)
 										double s = 0.0;
 										for (int i = 0; i < oldPlane.length; i++){
@@ -4546,7 +4878,7 @@ public class SuperTiles{
 										}
 										if (debugLevel > -1) {
 											//										if ((s > 5.0) || (dl > 0)){
-											if (s > 5.0){
+											if ((s > 5.0) || (dl > 0)){
 												System.out.println("planesSmoothStep() nsTile0="+nsTile0+
 														" num_merged="+num_merged + " rms = "+s+
 														" new_weight = "+new_planes[nsTile0][np0].getWeight()+
@@ -4563,6 +4895,9 @@ public class SuperTiles{
 										showDoubleFloatArrays sdfa_instance = new showDoubleFloatArrays(); // just for debugging?
 										sdfa_instance.showArrays(dbg_img, superTileSize, superTileSize, true, "smooth_step_x"+stx0+"_y"+sty0, titles);
 									}
+								} else { // if (this_new_plane != null)
+									this_new_plane = mod_planes[nsTile0][np0].clone();
+									new_planes[nsTile0][np0] = this_new_plane;
 								}
 							}
 						}								
@@ -5160,7 +5495,6 @@ public class SuperTiles{
 														TilePlanes.PlaneData other_plane = plane_triads[nSplit][nis].getPlaneToThis(
 																neib_planes[nsTile][np],
 																dl-1); // debugLevel);
-														// if (dl > 0) dbg_img[ 2 + dir] = other_plane.getPlaneDisparity(false);
 														if (other_plane != null){
 															if (plane_triads[nSplit][nis].getWeight() > 0.0){
 																plane_triads[nSplit][nis] = plane_triads[nSplit][nis].mergePlaneToThis(
@@ -5173,7 +5507,6 @@ public class SuperTiles{
 															} else {
 																plane_triads[nSplit][nis] = other_plane; 
 															}
-															//if (dl > 0) dbg_img[10 + dir] = this_new_plane.getPlaneDisparity(false);
 														} else {
 															plane_triads[nSplit][nis] = null;
 															break;

src/main/java/SurfaceData.java

+/**
+ **
+ ** SurfaceData - represents rectangualar area with per-tile disparity
+ **
+ ** Copyright (C) 2017 Elphel, Inc.
+ **
+ ** -----------------------------------------------------------------------------**
+ **  
+ **  SurfaceData.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 SurfaceData {
+	private int []     window;
+	private int        tilesX;
+	private int        tilesY;
+	private double []  disparity;
+	private boolean [] selected; // should include border
+	private boolean [] border;
+	private int []     sizes;
+	private int []     clusters;
+
+	public  SurfaceData(int [] tilesWH, int [] window)
+	{
+		this.tilesX = tilesWH[0];
+		this.tilesY = tilesWH[1];
+		this.window = window;
+		int len = window[2]*window[3];
+		this.disparity = new double [len];
+		this.selected =  new boolean [len];
+		this.border =    new boolean [len];
+	}
+	public SurfaceData(int tilesX, int tilesY, int x0, int y0, int width,int height)
+	{
+		this.tilesX = tilesX;
+		this.tilesY = tilesY;
+		int [] window = {x0,y0,width,height}; 
+		this.window = window;
+		int len = window[2]*window[3];
+		this.disparity = new double [len];
+		this.selected =  new boolean [len];
+		this.border =    new boolean [len];
+	}
+	
+	public void setCluster(int [] clusters)
+	{
+		this.clusters = clusters;
+	}
+	public int []  getCluster()
+	{
+		return this.clusters;
+	}
+	
+	public int getWidth()
+	{
+		return window[2];
+	}
+	public int getHeight()
+	{
+		return window[3];
+	}
+	public int getLeft()
+	{
+		return window[0];
+	}
+	public int getTop()
+	{
+		return window[1];
+	}
+	public int [] getWindow()
+	{
+		return window;
+	}
+	
+	
+	public void setDisparity(double [] disparity)
+	{
+		this.disparity = disparity;
+	}
+	public void setDisparity(int indx, double disparity)
+	{
+		this.disparity[indx] = disparity;
+	}
+	public double [] getDisparity()
+	{
+		return this.disparity;
+	}
+	public double getDisparity(int indx)
+	{
+		return this.disparity[indx];
+	}
+
+	
+	public void setSelected(boolean [] selected)
+	{
+		this.selected = selected;
+		this.sizes = null;
+	}
+	public void setSelected(int indx, boolean selected)
+	{
+		this.selected[indx] = selected;
+		this.sizes = null;
+	}
+	public boolean [] getSelected()
+	{
+		return this.selected;
+	}
+	public boolean getSelected(int indx)
+	{
+		return this.selected[indx];
+	}
+
+	
+	public void setBorder(boolean [] border)
+	{
+		this.border = border;
+		this.sizes = null;
+	}
+	public void setBorder(int indx, boolean border)
+	{
+		this.border[indx] = border;
+		this.sizes = null;
+	}
+	public boolean [] getBorder()
+	{
+		return this.border;
+	}
+	public boolean getBorder(int indx)
+	{
+		return this.border[indx];
+	}
+	
+	int getImageIndex(int indx)
+	{
+		int x = (indx % window[2]) + window[0];
+		int y = (indx / window[2]) + window[1];
+		if ((y <0) || (x < 0) || (y >= tilesY) || (x >= tilesX)){
+			return -1;
+		}
+		return  y * tilesX + x;
+	}
+	
+	public int getSize(
+			boolean include_border)
+	{
+		if (this.sizes == null){
+			int [] sizes = {0,0};
+			this.sizes = sizes;
+			for (int i = 0; i < selected.length; i++){
+				if (selected[i]) {
+					sizes[1]++;
+					if (!border[i]){
+						sizes[0]++;
+					}
+				}
+			}
+		}
+		return this.sizes[include_border ? 1 : 0];
+
+	}
+	public double [] getDisparityNaN(
+			boolean include_border)
+	{
+		double [] data = new double [tilesX*tilesY];
+		for (int i = 0; i < data.length; i++){
+			data[i] = Double.NaN;
+		}
+		for (int i = 0; i < disparity.length; i++){
+			if (selected[i] && (include_border || !border[i])){
+				data[getImageIndex(i)] = disparity[i];
+			}
+		}
+		return data;
+	}
+
+}

src/main/java/TilePlanes.java

@@ -53,7 +53,7 @@ public class TilePlanes {
 		double []   values =     null; // [3] -eigenvalues
 		int         num_points = 0;
 		double      weight =     0.0;
-		double []   center_xyz = null; // center of this supertile this plane center in world coordinates
+		double []   center_xyz = null; // center of this this "plane" (ellipsoid) center in world coordinates
 		double []   world_xyz =  null; // world coordinates of the nearest point of the plane, in meters
 		double []   world_v1 =   null; // world in-plane vector, corresponding to vectors[1] 
 		double []   world_v2 =   null; // world in-plane vector, corresponding to vectors[1]
@@ -566,6 +566,8 @@ public class TilePlanes {
 		 * Tilt disparity values around the supertile center (this.zxy) so constant disparity in the output
 		 * corresponds to the real world plane parallel to the provided one. Used to discriminate tiles by
 		 * the effective disparity value (disparity in the center of the supertile of the parallel plane)
+		 * Adding protection from behind the horizon - areas where disparity is negative zero the
+		 * result strengths
 		 * @param world_normal_xyz real world 3d vector of the plane normal (0.0, 1.0, 0.0 - horizontal)
 		 * @param disp_center dispariy in at the center of the supertile (to rotate around)
 		 * @param tile_sel multi-layer tile selection (or null to use all available tiles)
@@ -620,15 +622,21 @@ public class TilePlanes {
 							// then calculate disparity from z of that point
 							// inner product of transposed
 							double n_by_p = normal_row.times(w_xyz).get(0, 0);
-							double z = st_xyz.get(2, 0)*n_by_p / n_by_w;
 							if (disp_str[ml][1][indx] > 0){ // do not bother with zero-strength
+								double z;
+								if ((n_by_p * n_by_w) > 0.0) { 
+									z = st_xyz.get(2, 0)*n_by_p / n_by_w;
+									// convert z to disparity
+									eff_disp_str[ml][0][indx] = geometryCorrection.getDisparityFromZ (-z);
+								} else {
+									z = 0.0;
+									eff_disp_str[ml][1][indx] = 0.0; // behind the horizon
+								}
 								if (debugLevel > 1) {
 									System.out.println("dy = "+dy+", dx=" + dx+ " {"+w_xyz.get(0, 0)+","+w_xyz.get(1, 0)+","+w_xyz.get(2, 0)+"}"+" z="+z+" n_by_p = "+n_by_p
 											+" disp = "+disp_str[ml][0][indx]+" px = "+(px_py[0] + x)+" py = "+(px_py[1] + y));
 								}
 							}
-							// convert z to disparity
-							eff_disp_str[ml][0][indx] = geometryCorrection.getDisparityFromZ (-z);
 						}
 					}
 				}
@@ -1268,7 +1276,7 @@ public class TilePlanes {
 		 * @param useNaN replace unselected tiles with Double.NaN
 		 * @return array of disparity values for the plane (not including overlapped areas)
 		 */
-		public double[] getPlaneDisparity(
+		public double[] getSinglePlaneDisparity( // getPlaneDisparity(
 				boolean useNaN)
 		{
 			double [] disparities = new double[superTileSize*superTileSize];
@@ -1294,18 +1302,29 @@ public class TilePlanes {
 			return disparities;
 		}
 
+		public double[] getDoublePlaneDisparity(
+				boolean useNaN)
+		{
+			return getDoublePlaneDisparity(
+					true,
+					useNaN);
+		}
+		
 		/**
 		 * Get disparity values for the tiles of this overlapping supertile as [2*superTileSize * 2*superTileSize] array
+		 * @param useWorld calculate disparity in the real world (false - just px, py, disparity plane)
 		 * @param useNaN replace unselected tiles with Double.NaN
 		 * @return array of disparity values for the plane (not including overlapped areas)
 		 */
 		public double[] getDoublePlaneDisparity(
+				boolean useWorld,
 				boolean useNaN)
 		{
 			double [] disparities = new double[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 [] pxyc = getCenterPxPy(); // center of this supertile, not plane center 
 			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
@@ -1313,7 +1332,16 @@ public class TilePlanes {
 				for (int sx = -superTileSize; sx < superTileSize; sx++){
 					double x = tileSize * (sx + 0.5) + 0.5 - zxy[1];
 					if (plane_sel[indx] || !useNaN ||  (plane_sel==null)){
-						disparities[indx] = zxy[0] - (normal[1] * x + normal[2] * y)/normal[0];
+						if (useWorld) {
+							disparities[indx] = geometryCorrection.getPlaneDisparity( // disparity (at this center) for crossing other supertile plane
+									getWorldXYZ(this.correctDistortions), // will calculate if not yet done so. Should it use otherPd, not pd? and then clone later?
+									x + pxyc[0] + zxy[1],
+									y + pxyc[1] + zxy[2],
+									this.correctDistortions);
+						} else {
+							disparities[indx] = zxy[0] - (normal[1] * x + normal[2] * y)/normal[0];
+						}
+						
 					} else {
 						disparities[indx] = Double.NaN;
 					}
@@ -1327,7 +1355,18 @@ public class TilePlanes {
 				int     dir,
 				boolean useNaN)
 		{
-			double [] plane_all = getDoublePlaneDisparity(useNaN);
+			return getDoublePlaneDisparity(
+					true,
+					dir,
+					useNaN);
+		}
+
+		public double[] getDoublePlaneDisparity(
+				boolean useWorld,
+				int     dir,
+				boolean useNaN)
+		{
+			double [] plane_all = getDoublePlaneDisparity(useWorld,useNaN);
 			double [] plane = new double [superTileSize * superTileSize];
 			int [] start_index = {
 					superTileSize/2,                            // N    4
@@ -1345,18 +1384,33 @@ public class TilePlanes {
 		}
 
 		public double[] getTriplePlaneDisparity()
+		{
+			return getTriplePlaneDisparity(true);
+		}
+
+		public double[] getTriplePlaneDisparity(
+				boolean   useWorld)
 		{
 			double [] disparities = new double[9*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 [] pxyc = getCenterPxPy(); // center of this supertile, not plane center 
 			for (int sy = -3 * superTileSize / 2; sy < 3* superTileSize / 2; 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 = -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];
+					if (useWorld) {
+						disparities[indx] = geometryCorrection.getPlaneDisparity( // disparity (at this center) for crossing other supertile plane
+								getWorldXYZ(this.correctDistortions), // will calculate if not yet done so. Should it use otherPd, not pd? and then clone later?
+								x + pxyc[0] + zxy[1],
+								y + pxyc[1] + zxy[2],
+								this.correctDistortions);
+					} else {
+						disparities[indx] = zxy[0] - (normal[1] * x + normal[2] * y)/normal[0];
+					}
 					indx++;
 				}
 			}
@@ -1366,7 +1420,14 @@ public class TilePlanes {
 		public double[] getTriplePlaneDisparity(
 				int     dir)
 		{
-			double [] plane_all = getTriplePlaneDisparity();
+			return getTriplePlaneDisparity(true, dir);
+		}
+
+		public double[] getTriplePlaneDisparity(
+				boolean   useWorld,
+				int     dir)
+		{
+			double [] plane_all = getTriplePlaneDisparity(useWorld);
 			double [] plane = new double [superTileSize * superTileSize];
 			int [] start_index = {
 					superTileSize,                               // N    8
@@ -1383,19 +1444,36 @@ public class TilePlanes {
 			return plane;
 		}
 
+		public double[][] getDoublePlaneDisparityStrength(
+				double [] window,
+				boolean   use_sel,
+				boolean   divide_by_area,
+				double    scale_projection,
+				int       debugLevel)
+		{		return getDoublePlaneDisparityStrength(
+				true,
+				window,
+				use_sel,
+				divide_by_area,
+				scale_projection,
+				debugLevel);
+		}		
 		/**
 		 * 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
 		 * Sharp weights - when selecting the best match - use exponent of (delta_disp) ^2 ?
-		 * Or divide weight by ellipse arae?
+		 * Or divide weight by ellipse area?
+		 * @param useWorld calculate disparity in the real world (false - just px, py, disparity plane)
 		 * @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 divide_by_area divide weights by ellipsoid area
 		 * @param scale_projection use plane ellipsoid projection for weight: 0 - do not use, > 0 linearly scale ellipsoid 
 		 * @return a pair of arrays {disparity, strength}, each [2*superTileSize * 2*superTileSize]
 		 */
+		// obsolete, convert to another version ?
 		public double[][] getDoublePlaneDisparityStrength(
+				boolean   useWorld,
 				double [] window,
 				boolean   use_sel,
 				boolean   divide_by_area,
@@ -1435,15 +1513,32 @@ public class TilePlanes {
 					}
 				}
 			}
-
+			double [] pxyc = getCenterPxPy(); // center of this supertile, not plane center
+			if (debugLevel > 0) {
+				double world_disp = geometryCorrection.getPlaneDisparity( // disparity (at this center) for crossing other supertile plane
+						getWorldXYZ(this.correctDistortions), // will calculate if not yet done so. Should it use otherPd, not pd? and then clone later?
+						zxy[1] + pxyc[0],
+						zxy[2] + pxyc[1],
+						this.correctDistortions);
+				System.out.println("getDoublePlaneDisparityStrength(): zxy = ["+zxy[0]+", "+zxy[1]+", "+zxy[2]+"], disp = "+world_disp);
+			}
 			int indx = 0;
 			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];
+				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];
+					if (useWorld) {
+						disp_strength[0][indx] = geometryCorrection.getPlaneDisparity( // disparity (at this center) for crossing other supertile plane
+								getWorldXYZ(this.correctDistortions), // will calculate if not yet done so. Should it use otherPd, not pd? and then clone later?
+								x + pxyc[0] + zxy[1],
+								y + pxyc[1] + zxy[2],
+								this.correctDistortions);
+					} else {
+						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;
@@ -1464,12 +1559,33 @@ public class TilePlanes {
 			return disp_strength;
 		}
 
+		public double[][] getDoublePlaneDisparityStrength(
+				double [] window,
+				int       dir,
+				boolean   use_sel,
+				boolean   divide_by_area,
+				double    scale_projection,
+				double    fraction_uni,
+				int       debugLevel)
+		{		
+			return getDoublePlaneDisparityStrength(
+					true,
+					window,
+					dir,
+					use_sel,
+					divide_by_area,
+					scale_projection,
+					fraction_uni,
+					debugLevel);
+		}
+		
 		/**
 		 * 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
 		 * Sharp weights - when selecting the best match - use exponent of (delta_disp) ^2 ?
 		 * Or divide weight by ellipse area?
+		 * @param useWorld calculate disparity in the real world (false - just px, py, disparity plane)
 		 * @param window null or window function as [2*superTileSize * 2*superTileSize] array
 		 * @param dir - source tile shift from the target: -1 center, 0 - N, 1 - NE
 		 * @param use_sel use plane selection (this.sel_mask) to select only some part of the plane
@@ -1481,6 +1597,7 @@ public class TilePlanes {
 		 * 
 		 */
 		public double[][] getDoublePlaneDisparityStrength(
+				boolean   useWorld,
 				double [] window,
 				int       dir,
 				boolean   use_sel,
@@ -1522,9 +1639,7 @@ public class TilePlanes {
 					}
 				}
 			}
-//			int ss1 = superTileSize / 2;
 			int ss2 = superTileSize;
-//			int ss3 = 3 *ss1;
 			int ss4 = 2 * superTileSize;
 			
 			
@@ -1535,27 +1650,30 @@ public class TilePlanes {
 					{ss2, ss4,   0, ss2, -ss2,  ss2 },  // NE	
 					{  0, ss4,   0, ss2,    0,  ss2 },  // E	
 					{  0, ss2,   0, ss2,  ss2,  ss2 },  // SE	
-//					{  0, ss2,   0, ss2,  ss2,    0 },  // S	
 					{  0, ss2,   0, ss4,  ss2,    0 },  // S	
 					{  0, ss2, ss2, ss4,  ss2, -ss2 },  // SW	
 					{  0, ss4, ss2, ss4,    0, -ss2 },  // W	
-//					{ss2, ss4,   0, ss4, -ss2, -ss2 }}; // NW
 					{ss2, ss4, ss2, ss4, -ss2, -ss2 }}; // NW
 			int dir1 = dir + 1;
-
-//			for (int sy = -superTileSize; sy < superTileSize; sy++){
+			double [] pxyc = getCenterPxPy(); // center of this supertile, not plane center 
 			for (int iy = offsets[dir1][0]; iy < offsets[dir1][1]; iy++){
 				// 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 * (iy - ss2 + 0.5) + 0.5  - zxy[2];
 				int oy = iy + offsets[dir1][4]; //vert index in the result tile
-//				for (int sx = -superTileSize; sx < superTileSize; sx++){
 				for (int ix = offsets[dir1][2]; ix < offsets[dir1][3]; ix++){
 					double x = tileSize * (ix - ss2 + 0.5) + 0.5 - zxy[1];
-//					int indx = ss2 * oy + ix + offsets[dir1][5];
 					int indx = ss4 * oy + ix + offsets[dir1][5];
 					int indx_i = iy * ss4 + ix; // input index
-					disp_strength[0][indx] = zxy[0] - (normal[1] * x + normal[2] * y)/normal[0];
+					if (useWorld) {
+						disp_strength[0][indx] = geometryCorrection.getPlaneDisparity( // disparity (at this center) for crossing other supertile plane
+								getWorldXYZ(this.correctDistortions), // will calculate if not yet done so. Should it use otherPd, not pd? and then clone later?
+								x + pxyc[0] + zxy[1],
+								y + pxyc[1] + zxy[2],
+								this.correctDistortions);
+					} else {
+						disp_strength[0][indx] = zxy[0] - (normal[1] * x + normal[2] * y)/normal[0];
+					}
 					double w = weight;
 					if ((w > 0.0) && (scale_projection > 0.0)){
 						double [] xy = {x,y};
@@ -1576,13 +1694,30 @@ public class TilePlanes {
 			return disp_strength;
 		}
 
-
+		public double[][] getSinglePlaneDisparityStrength(
+				double [] window,
+				int       dir,
+				boolean   use_sel,
+				boolean   divide_by_area,
+				double    scale_projection,
+				int       debugLevel)
+		{
+			return getSinglePlaneDisparityStrength(
+					true,
+					window,
+					dir,
+					use_sel,
+					divide_by_area,
+					scale_projection,
+					debugLevel);
+		}
 		/**
 		 * Get disparity values for the tiles of this overlapping supertile as [superTileSize * superTileSize] array
 		 * and weights combined from provided window function, optional selection and using ellipsoid projection on the
 		 * px, py plane (constant disparity
 		 * Sharp weights - when selecting the best match - use exponent of (delta_disp) ^2 ?
 		 * Or divide weight by ellipse arae?
+		 * @param useWorld calculate disparity in the real world (false - just px, py, disparity plane)
 		 * @param window null or window function as [2*superTileSize * 2*superTileSize] array
 		 * @param dir - source tile shift from the targer: -1 center, 0 - N, 1 - NE
 		 * @param use_sel use plane selection (this.sel_mask) to select only some part of the plane
@@ -1593,6 +1728,7 @@ public class TilePlanes {
 		 * 
 		 */
 		public double[][] getSinglePlaneDisparityStrength(
+				boolean   useWorld,
 				double [] window,
 				int       dir,
 				boolean   use_sel,
@@ -1651,19 +1787,26 @@ public class TilePlanes {
 					{ss1, ss3, ss3, ss4, -ss1, -ss3 },  // W	
 					{ss3, ss4, ss3, ss4, -ss3, -ss3 }}; // NW
 			int dir1 = dir + 1;
-
-//			for (int sy = -superTileSize; sy < superTileSize; sy++){
+			double [] pxyc = getCenterPxPy(); // center of this supertile, not plane center 
 			for (int iy = offsets[dir1][0]; iy < offsets[dir1][1]; iy++){
 				// 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 * (iy - ss2 + 0.5) + 0.5  - zxy[2];
 				int oy = iy + offsets[dir1][4]; //vert index in the result tile
-//				for (int sx = -superTileSize; sx < superTileSize; sx++){
 				for (int ix = offsets[dir1][2]; ix < offsets[dir1][3]; ix++){
 					double x = tileSize * (ix - ss2 + 0.5) + 0.5 - zxy[1];
 					int indx = ss2 * oy + ix + offsets[dir1][5];
 					int indx_i = iy * ss4 + ix; // ss2;
-					disp_strength[0][indx] = zxy[0] - (normal[1] * x + normal[2] * y)/normal[0];
+					if (useWorld) {
+						disp_strength[0][indx] = geometryCorrection.getPlaneDisparity( // disparity (at this center) for crossing other supertile plane
+								getWorldXYZ(this.correctDistortions), // will calculate if not yet done so. Should it use otherPd, not pd? and then clone later?
+								x + pxyc[0] + zxy[1],
+								y + pxyc[1] + zxy[2],
+								this.correctDistortions);
+					} else {
+						disp_strength[0][indx] = zxy[0] - (normal[1] * x + normal[2] * y)/normal[0];
+					}
+					
 					double w = weight;
 					if (window != null) w *= window[indx_i];
 					if (use_sel && (sel_mask != null) && !(sel_mask[indx_i])) w = 0.0;
@@ -1703,7 +1846,54 @@ public class TilePlanes {
 			return new Matrix (ar);
 		}
 
+		/** 
+		 * Get sin squared of the angle between planes in the real world
+		 * @param otherPd other plane data
+		 * @param correct_distortions true if the lens distortions should be corrected 
+		 * @return sine squared of the angle between normals to the planes
+		 */
+		public double getWorldSin2(
+				PlaneData otherPd,
+				boolean correct_distortions)
+		{
+			Matrix this_wv =  new Matrix(this.getWorldXYZ(correct_distortions, 0),3);
+			Matrix other_wv = new Matrix(otherPd.getWorldXYZ(correct_distortions, 0),3);
+			Matrix cp = cross3d(this_wv, other_wv);
+			double cp2 = cp.transpose().times(cp).get(0, 0);
+			double this_wv2 = this_wv.transpose().times(this_wv).get(0, 0);
+			double other_wv2 = other_wv.transpose().times(other_wv).get(0, 0);
+			return cp2/(this_wv2 * other_wv2);
+		}
+		public double getWorldSin2(
+				PlaneData otherPd)
+		{
+			return getWorldSin2(otherPd, this.correctDistortions);
+		}
 
+		/**
+		 * Get distance squared from the other plane and the center of the current "plane" (ellipsoid) 
+		 * @param otherPd other plane data
+		 * @param correct_distortions true if the lens distortions should be corrected 
+		 * @return distance squared from other plane to the (ellipsoid) center of this one (in meters)
+		 */
+		public double getWorldPlaneDist2(
+				PlaneData otherPd,
+				boolean correct_distortions)
+		{
+			Matrix this_center =  new Matrix(this.getCenterXYZ(correct_distortions, 0),3);
+			Matrix other_wv =     new Matrix(otherPd.getWorldXYZ(correct_distortions, 0),3);
+			double w_dot_w = other_wv.transpose().times(other_wv).get(0, 0);
+			double w_dot_p = other_wv.transpose().times(this_center).get(0, 0);
+			return (w_dot_w - w_dot_p)*(w_dot_w - w_dot_p)/w_dot_w;
+		}
+
+		public double getWorldPlaneDist2(
+				PlaneData otherPd)
+		{
+			return getWorldPlaneDist2(otherPd, this.correctDistortions);
+		}
+		
+		
 		/**
 		 * 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)
@@ -2047,6 +2237,26 @@ public class TilePlanes {
 				return world_v1;
 			}
 		}
+		public double [] getCenterXYZ(
+				boolean correct_distortions,
+				int debugLevel)
+		{
+			double delta = 0.0001;
+			if (center_xyz != null) return center_xyz;
+			setCorrectDistortions(correct_distortions);
+			// get pixel coordinates of the plane origin point
+			double [] px_py = getCenterPxPy();
+			
+			double px = px_py[0] + zxy[1];
+			double py = px_py[1] + zxy[2];
+			double disp =  zxy[0];
+			center_xyz = geometryCorrection.getWorldCoordinates(
+					px,
+					py,
+					disp,
+					this.correctDistortions);
+			return center_xyz;
+		}
 
 		public double [] getWorldXYZ(
 				boolean correct_distortions,
@@ -2056,14 +2266,7 @@ public class TilePlanes {
 			if (world_xyz != null) return world_xyz;
 			setCorrectDistortions(correct_distortions);
 			// get pixel coordinates of the plane origin point
-//			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];
 			double [] px_py = getCenterPxPy();
-//			if ((px_py[0] == 1760) && (px_py[1] == 1056)){ // 27, 15
-//				System.out.println("getWorldXYZ, px_py = {"+px_py[0]+","+px_py[1]+"}");
-//				debugLevel = 2;
-//			}
-			
 			double px = px_py[0] + zxy[1];
 			double py = px_py[1] + zxy[2];
 			double disp =  zxy[0];
@@ -2075,15 +2278,9 @@ public class TilePlanes {
 			Matrix xyz = new Matrix(center_xyz, 3); // column matrix
 			Matrix dpxpy = new Matrix(vectors[0],3); // 3 rows, 1 column
 			if (debugLevel > 0){
-//				if (sTileXY[0] == 27 ){
-//					System.out.println("STOP");
-//					
-//				}
 				System.out.println("getWorldXYZ("+sTileXY[0]+","+sTileXY[1]+"), correctDistortions="+correctDistortions+", xyz= {"+
 						xyz.get(0, 0)+","+xyz.get(1, 0)+","+xyz.get(2, 0)+"}, weight = "+getWeight());
-//				xyz.print(10, 6); // w,d
 				
-//				double [] dpxpy = geometryCorrection.getImageCoordinates(xyz.getColumnPackedCopy(),this.correctDistortions);
 				System.out.println("getWorldXYZ("+sTileXY[0]+","+sTileXY[1]+"): disp, px, py="+disp+","+px+","+py); // + ", reversed:"+dpxpy[0]+","+dpxpy[1]+","+dpxpy[2]);
 				
 				Matrix xyz1 = new Matrix(geometryCorrection.getWorldCoordinates(
@@ -2112,8 +2309,6 @@ public class TilePlanes {
 			world_v2 = v2.getColumnPackedCopy();
 			
 			Matrix norm_xyz =  cross3d(v1,v2);
-//			norm_xyz = norm_xyz.times(1.0/norm_xyz.normF()); // unity normal vector;
-//			norm_xyz = jacobian.times(dpxpy); // plane normal vector in world xyz
 			if (debugLevel > 0){
 				System.out.println("getWorldXYZ("+sTileXY[0]+","+sTileXY[1]+"): norm_xyz={"+
 						norm_xyz.get(0, 0)+", "+norm_xyz.get(1, 0)+", "+norm_xyz.get(2, 0)+"}, (dpxpy={"+
@@ -2145,7 +2340,6 @@ public class TilePlanes {
 				System.out.println("+getWorldXYZ("+sTileXY[0]+","+sTileXY[1]+"): unit plane normal={"+
 						norm_xyz.get(0, 0)+", "+norm_xyz.get(1, 0)+", "+norm_xyz.get(2, 0)+"})");
 				double dotprod = xyz.transpose().times(norm_xyz).get(0,0);
-//				Matrix wn = norm_xyz.times(-dotprod);
 				Matrix wn = norm_xyz.times(dotprod);
 				System.out.println(":getWorldXYZ("+sTileXY[0]+","+sTileXY[1]+"): xyz.transpose().times(norm_xyz).get(0,0) ="+dotprod);
 				System.out.println("?getWorldXYZ("+sTileXY[0]+","+sTileXY[1]+"):  plane normal={"+
@@ -2161,7 +2355,6 @@ public class TilePlanes {
 			
 			// convert plane normal vector to world coordinates
 			//world_xyz
-//			world_xyz = norm_xyz.times(-(xyz.transpose().times(norm_xyz).get(0,0))).getColumnPackedCopy();
 			world_xyz = norm_xyz.times((xyz.transpose().times(norm_xyz).get(0,0))).getColumnPackedCopy();
 			return world_xyz;
 		}

src/main/java/TileProcessor.java

@@ -2949,6 +2949,7 @@ public class TileProcessor {
 			final boolean     updateStatus,
 			final int         debugLevel)
 	{
+		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);
 		SuperTiles st = scan_prev.getSuperTiles();
@@ -2957,6 +2958,10 @@ public class TileProcessor {
 			return false;
 		}
 
+		tileSurface.testSimpleConnected(
+				clt_parameters.tileX,
+				clt_parameters.tileY);
+		
 		double [][][]  dispStrength = st.getDisparityStrengths(
 				clt_parameters.stMeasSel); // int        stMeasSel) //            = 1;      // Select measurements for supertiles : +1 - combo, +2 - quad +4 - hor +8 - vert)
 
@@ -2972,7 +2977,8 @@ public class TileProcessor {
         // assign tiles that do not depend on other assigned tiles - single pass
 		if (clt_parameters.tsEnSingle) {
 			int [] stats=  tileSurface.assignTilesToSurfaces(
-					clt_parameters.tsNoEdge ,         // final boolean       noEdge,
+					clt_parameters.tsNoEdge ,       // final boolean       noEdge,
+					clt_parameters.tsUseCenter,     // final boolean       useCenter,
 					clt_parameters.tsMaxDiff,       //final double        maxDiff,
 					clt_parameters.tsMinDiffOther,  //final double        minDiffOther, // should be >= maxDiff
 					clt_parameters.tsMinStrength,   //final double        minStrength,
@@ -3000,7 +3006,8 @@ public class TileProcessor {
 		if (clt_parameters.tsEnMulti) {
 			for (int nTry = 0; nTry < 100; nTry++) {
 				int [] stats=  tileSurface.assignTilesToSurfaces(
-						clt_parameters.tsNoEdge ,         // final boolean       noEdge,
+						clt_parameters.tsNoEdge ,       // final boolean       noEdge,
+						clt_parameters.tsUseCenter,     // final boolean       useCenter,
 						clt_parameters.tsMaxDiff,       //final double        maxDiff,
 						clt_parameters.tsMinDiffOther,  //final double        minDiffOther, // should be >= maxDiff
 						clt_parameters.tsMinStrength,   //final double        minStrength,
@@ -3109,9 +3116,9 @@ public class TileProcessor {
 
 		tileSurface.showClusterStats(
 				cluster_stats1, // int [][] cluster_stats,
-				clt_parameters.tsNumClust); // int max_clusters){
+				-1); //clt_parameters.tsNumClust); // int max_clusters){
 		
-		if (clt_parameters.tsShow){
+		if (clt_parameters.tsShow && (debugLevel > -1)){
 			tileSurface.showClusters(
 					"clusters_individual", // String title,
 					cluster_stats1, // int [][] cluster_stats,
@@ -3133,8 +3140,15 @@ public class TileProcessor {
 
 		tileSurface.showClusterStats(
 				cluster_stats1a, // int [][] cluster_stats,
-				clt_parameters.tsNumClust); // int max_clusters){
+				-1); //clt_parameters.tsNumClust); // int max_clusters){
 
+		if (clt_parameters.tsShow && (debugLevel > -1)){
+			tileSurface.showClusters(
+					"clusters_individual_fixed", // String title,
+					cluster_stats1a, // int [][] cluster_stats,
+					clt_parameters.tsNumClust, // int max_clusters
+					clusters1a); // int [][] clusters); // final double [][][] dispStrength)
+		}
 		
 		
 		
@@ -3159,9 +3173,9 @@ public class TileProcessor {
 
 		tileSurface.showClusterStats(
 				cluster_stats2, // int [][] cluster_stats,
-				clt_parameters.tsNumClust); // int max_clusters){
+				-1); //clt_parameters.tsNumClust); // int max_clusters){
 		
-		if (clt_parameters.tsShow){
+		if (clt_parameters.tsShow && (debugLevel > -1)){
 			tileSurface.showClusters(
 					"clusters_merged", // String title,
 					cluster_stats2, // int [][] cluster_stats,
@@ -3184,9 +3198,9 @@ public class TileProcessor {
 
 		tileSurface.showClusterStats(
 				cluster_stats2a, // int [][] cluster_stats,
-				clt_parameters.tsNumClust); // int max_clusters){
+				-1); //clt_parameters.tsNumClust); // int max_clusters){
 
-		if (clt_parameters.tsShow){
+		if (clt_parameters.tsShow && (debugLevel > -1)){
 			tileSurface.showClusters(
 					"clusters_merged_split", // String title,
 					cluster_stats2a, // int [][] cluster_stats,
@@ -3194,13 +3208,28 @@ public class TileProcessor {
 					clusters2a); // int [][] clusters); // final double [][][] dispStrength)
 		}
 		
-		tileSurface.mergeAndGrow( // TODO: add result
+		tileSurface.mergeAndGrow( // T ODO: add result
 				assigned_sel, // final boolean [][] sel_in,
 				0, // final int           debugLevel,
 				clt_parameters.tileX,
 				clt_parameters.tileY);
 		
-		
+		if (clt_parameters.tsShow){
+			int numToShow = clt_parameters.tsNumClust;
+			int num_surf = tileSurface.getSurfaceDataLength();
+			if (numToShow > num_surf) numToShow = num_surf;
+			String [] titles = new String[2*numToShow];
+			double [][] img_data = new double [2*numToShow][];
+			for (int i = 0; i < numToShow; i++){
+				titles[2 * i    ] = "disp_"+i;
+				titles[2 * i + 1] = "pure_disp_"+i;
+				img_data[2 * i    ] = tileSurface.getSurfaceData(i).getDisparityNaN(true);
+				img_data[2 * i + 1] = tileSurface.getSurfaceData(i).getDisparityNaN(false);
+				
+			}
+			showDoubleFloatArrays sdfa_instance = new showDoubleFloatArrays(); // just for debugging?
+			sdfa_instance.showArrays(img_data, tilesX, tilesY, true, "final_surfaces",titles);
+		}		
 		return true;
 	}	
 	
@@ -3253,8 +3282,8 @@ public class TileProcessor {
 			dbg_hist[1] = scan_prev.getSuperTiles().showDisparityHistogram();
 		}
 
-//		SuperTiles st = 
-				scan_prev.setSuperTiles(
+		//		SuperTiles st = 
+		scan_prev.setSuperTiles(
 				clt_parameters.stStepNear,       // double     step_disparity,
 				clt_parameters.stStepFar,        // double     step_near,
 				clt_parameters.stStepThreshold,  // double     step_threshold,
@@ -3280,7 +3309,7 @@ public class TileProcessor {
 			showDoubleFloatArrays sdfa_instance = new showDoubleFloatArrays(); // just for debugging?
 			sdfa_instance.showArrays(dbg_hist, hist_width0, hist_height0, true, "disparity_supertiles_histograms",dbg_st_titles);
 		}
-		
+
 		if (clt_parameters.stSmplMode){ // just temporary making it last
 			scan_prev.setSuperTiles(
 					clt_parameters.stStepNear,       // double     step_disparity,
@@ -3297,12 +3326,23 @@ public class TileProcessor {
 					clt_parameters.stSmplRms,   // Maximal RMS of the remaining tiles in a sample
 					clt_parameters.stMeasSel); // bitmask of the selected measurements for supertiles : +1 - combo, +2 - quad +4 - hor +8 - vert
 		}
-		
-		
-		
+
+
+
 		SuperTiles st = scan_prev.getSuperTiles();
+		/*
+		TileSurface tileSurface = new TileSurface(
+				st.tileProcessor.getTileSize(),      // int tileSize,
+				st.tileProcessor.getSuperTileSize(), // int superTileSize,
+				st.tileProcessor.getTilesX(),        // int tilesX,
+				st.tileProcessor.getTilesY(),        // int tilesY,
+				geometryCorrection,               // GeometryCorrection geometryCorrection,
+				st.tileProcessor.threadsMax);        // int threadsMax);
+				*/ 
+		st.setTileSurface(geometryCorrection); // tileSurface);
 		
-// moved here
+
+		// moved here
 		if (clt_parameters.dbg_migrate) {
 			// separate each supertile data into clusters, trying both horizontal and perpendicular to view planes
 			double []  world_hor = {0.0, 1.0, 0.0};
@@ -3332,8 +3372,8 @@ public class TileProcessor {
 					clt_parameters.msUseSel,        // final boolean    msUseSel,        // final boolean                   use_sel,
 					clt_parameters.msDivideByArea,  // final boolean    msDivideByArea,  // final boolean                   divide_by_area,
 					clt_parameters.msScaleProj,     //final double     msScaleProj,     // final double                    scale_projection,
-					
-					
+
+
 					clt_parameters.stSmallDiff,     //       = 0.4;   // Consider merging initial planes if disparity difference below
 					clt_parameters.stHighMix,       // stHighMix         = 0.4;   // Consider merging initial planes if jumps between ratio above
 					world_hor,                      // final double []  world_hor, // horizontal plane normal (default [0.0, 1.0, 0.0])
@@ -3365,31 +3405,8 @@ public class TileProcessor {
 					0, // -1, // debugLevel,                  // final int        debugLevel)
 					clt_parameters.tileX,
 					clt_parameters.tileY);
-/*			
-			st.processPlanes3(
-					null, // final boolean [] selected, // or null
-					0.3, // final double     min_disp,
-					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
-					clt_parameters.plFractOutliers, //      =   0.3;  // Maximal fraction of outliers to remove
-					clt_parameters.plMaxOutliers, //        =    20;  // Maximal number of outliers to remove\
-					clt_parameters.plPreferDisparity,
-					geometryCorrection,
-					clt_parameters.correct_distortions,
-//					false, // clt_parameters.dbg_migrate && clt_parameters.stSmplMode , // final boolean                    smplMode, //        = true;   // Use sample mode (false - regular tile mode)
-					clt_parameters.stSmplMode , // final boolean                    smplMode, //        = true;   // Use sample mode (false - regular tile mode)
-					clt_parameters.stSmplSide , // final int                        smplSide, //        = 2;      // Sample size (side of a square)
-					clt_parameters.stSmplNum , // final int                        smplNum, //         = 3;      // Number after removing worst
-					clt_parameters.stSmplRms , // final double                     smplRms, //         = 0.1;    // Maximal RMS of the remaining tiles in a sample
-					
-					0, // -1, // debugLevel,                  // final int        debugLevel)
-					clt_parameters.tileX,
-					clt_parameters.tileY);
-*/					
 		}
-		
+
 		showDoubleFloatArrays sdfa_instance = null;
 		if (debugLevel > -1) sdfa_instance = new showDoubleFloatArrays(); // just for debugging?
 
@@ -3401,15 +3418,19 @@ public class TileProcessor {
 
 		st.selectNeighborPlanesMutual(
 				clt_parameters.plWorstWorsening, // final double worst_worsening,
-				clt_parameters.plOKMergeEigen,   // final double okMergeEigen,
 				clt_parameters.plWeakWorsening,  // final double worst_worsening,
+				clt_parameters.plOKMergeEigen,   // final double okMergeEigen,
+				clt_parameters.plMaxWorldSin2,   // final double maxWorldSin2,
+
 				clt_parameters.plDispNorm,
 				clt_parameters.plMaxEigen,
 				clt_parameters.plMinStrength,
 				0, // final int debugLevel)
 				clt_parameters.tileX,
 				clt_parameters.tileY);
-		
+		st.detectTriangularConflicts(
+				1); // final int debugLevel)
+
 		if (clt_parameters.plSplitApply) {
 			while (true) {
 				int num_added = 0;
@@ -3424,29 +3445,29 @@ public class TileProcessor {
 				if (num_added == 0) break;
 			}
 		}
-		
+
 		double [][][]  dispStrength = st.getDisparityStrengths(
 				clt_parameters.stMeasSel); // int        stMeasSel) //            = 1;      // Select measurements for supertiles : +1 - combo, +2 - quad +4 - hor +8 - vert)
 		boolean [][] tileSel =  st.getMeasurementSelections(
 				clt_parameters.stMeasSel); // int        stMeasSel) //            = 1;      // Select measurements for supertiles : +1 - combo, +2 - quad +4 - hor +8 - vert)
 
-		
-		
-		
+
+
+
 		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
+						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); 
 
-				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,
@@ -3490,14 +3511,18 @@ public class TileProcessor {
 
 			st.selectNeighborPlanesMutual(
 					clt_parameters.plWorstWorsening, // final double worst_worsening,
-					clt_parameters.plOKMergeEigen,   // final double okMergeEigen,
 					clt_parameters.plWeakWorsening,  // final double worst_worsening,
+					clt_parameters.plOKMergeEigen,   // final double okMergeEigen,
+					clt_parameters.plMaxWorldSin2,   // final double maxWorldSin2,
 					clt_parameters.plDispNorm,
 					clt_parameters.plMaxEigen,
 					clt_parameters.plMinStrength,
 					0, // final int debugLevel)
 					clt_parameters.tileX,
 					clt_parameters.tileY);
+			st.detectTriangularConflicts(
+					1); // final int debugLevel)
+			
 		}
 
 		while (true) {
@@ -3516,8 +3541,8 @@ public class TileProcessor {
 			}
 			if (num_added == 0) break;
 		}
-		
-		
+
+
 		TilePlanes.PlaneData [][] planes_mod = null;
 
 		// smooth planes (by averaging with neighbors and the "measured" one with variable "pull")  
@@ -3543,7 +3568,7 @@ public class TileProcessor {
 				clt_parameters.plMaxEigen, //  maxEigen,
 				clt_parameters.plMinStrength, //  minWeight,
 				st.getPlanesMod());
-		
+
 		if (clt_parameters.show_planes){
 			double [] split_lines = st.showSplitLines(
 					split_planes,
@@ -3554,20 +3579,20 @@ public class TileProcessor {
 			int [] wh = st.getShowPlanesWidthHeight();
 			double [][] plane_data_nonan = st.getShowPlanes(
 					(planes_mod != null) ? st.getPlanesMod():st.getPlanes(),
-					clt_parameters.plMinStrength, //  minWeight,
-					clt_parameters.plMaxEigen, //  maxEigen,
-					clt_parameters.plDispNorm,
-					false, //boolean use_NaN)
-					0.0,
-					10.0);
+							clt_parameters.plMinStrength, //  minWeight,
+							clt_parameters.plMaxEigen, //  maxEigen,
+							clt_parameters.plDispNorm,
+							false, //boolean use_NaN)
+							0.0,
+							10.0);
 			double [][] plane_data_nan = st.getShowPlanes(
 					(planes_mod != null) ? st.getPlanesMod():st.getPlanes(),
-					clt_parameters.plMinStrength, //  minWeight,
-					clt_parameters.plMaxEigen, //  maxEigen,
-					clt_parameters.plDispNorm,
-					true, //boolean use_NaN)
-					0.0,
-					10.0);
+							clt_parameters.plMinStrength, //  minWeight,
+							clt_parameters.plMaxEigen, //  maxEigen,
+							clt_parameters.plDispNorm,
+							true, //boolean use_NaN)
+							0.0,
+							10.0);
 			double [][] plane_data = new double [plane_data_nonan.length + plane_data_nan.length + 2][];
 			int indx = 0;
 			for (int i = 0; i < plane_data_nonan.length; i++){
@@ -3586,10 +3611,10 @@ public class TileProcessor {
 				if (Double.isNaN(plane_data[indx][i])) plane_data[indx][i] = 0.0; 
 				if (plane_data[indx-1][i] > 0) plane_data[indx][i] = Double.NaN; 
 			}
-			
+
 			sdfa_instance.showArrays(plane_data, wh[0], wh[1], true, "plane_data");
-			
-			
+
+/*
 			TileSurface tileSurface = new TileSurface(
 					st.tileProcessor.getTileSize(),      // int tileSize,
 					st.tileProcessor.getSuperTileSize(), // int superTileSize,
@@ -3598,8 +3623,8 @@ public class TileProcessor {
 					geometryCorrection,               // GeometryCorrection geometryCorrection,
 					st.tileProcessor.threadsMax);        // int threadsMax); 
 			st.setTileSurface(tileSurface);
-			
-					tileSurface.createTileShells (
+*/
+			st.tileSurface.createTileShells (
 					clt_parameters.msUseSel,            // final boolean                   use_sel,
 					clt_parameters.msDivideByArea,      // final boolean                   divide_by_area,
 					clt_parameters.msScaleProj,         // final double                    scale_projection,
@@ -3609,21 +3634,21 @@ public class TileProcessor {
 					clt_parameters.tileX,
 					clt_parameters.tileY);
 
-			tileSurface.InitTilesAssignment(
+			st.tileSurface.InitTilesAssignment(
 					true,
 					dispStrength, // final double [][][]                            dispStrength,
 					tileSel,      // final boolean [][]                             tileSel,
-	                debugLevel);    // final int                                      debugLevel,
-			
-			
+					debugLevel);    // final int                                      debugLevel,
+
+
 			if (debugLevel > -10){
 				return; // just cut off the rest
 			}
-			
-//*******************************************************************************
-			
-			
-	        // detect connected "shells" by running wave algorithm over multi-plane supertiles, create integer array (same structure as planes and selected_planes
+
+			//*******************************************************************************
+
+
+			// detect connected "shells" by running wave algorithm over multi-plane supertiles, create integer array (same structure as planes and selected_planes
 			// Each element [per supertile][per disparity plane] is either 0 (not used) or unique shell index (started from 1)
 			int [][] shells = st.createSupertileShells(
 					selected_planes, // boolean[][]               selection, // may be null
@@ -3649,7 +3674,7 @@ public class TileProcessor {
 					10.0);               // double arrow_white)
 			// save surfaces with SuperTiles instance. They can be used to snap to for the per-tile disparity maps.
 			st.setSurfaces(surfaces);
-			
+
 			plane_data = st.getShowShells(
 					0, // int  nlayer, // over multi-layer - do not render more than nlayer on top of each other
 					st.getPlanesMod(),   // TilePlanes.PlaneData [][] planes,
@@ -3660,7 +3685,7 @@ public class TileProcessor {
 					false,               // boolean use_NaN,
 					10.0,                 // double arrow_dark,
 					10.0);               // double arrow_white)
-			
+
 			sdfa_instance.showArrays(plane_data, wh[0], wh[1], true, "shells_all");
 
 			plane_data = st.getShowShells(
@@ -3673,9 +3698,9 @@ public class TileProcessor {
 					false,               // boolean use_NaN,
 					10.0,                 // double arrow_dark,
 					10.0);               // double arrow_white)
-			
+
 			sdfa_instance.showArrays(plane_data, wh[0], wh[1], true, "shells2_all");
-			
+
 			plane_data = st.getShowShells(
 					0, // int  nlayer, // over multi-layer - do not render more than nlayer on top of each other
 					st.getPlanesMod(),   // TilePlanes.PlaneData [][] planes,
@@ -3700,7 +3725,7 @@ public class TileProcessor {
 					10.0);               // double arrow_white)
 			sdfa_instance.showArrays(plane_data, wh[0], wh[1], true, "shells2_all_arrows");
 
-			
+
 			plane_data = st.getShowShells(
 					0, // int  nlayer, // over multi-layer - do not render more than nlayer on top of each other
 					st.getPlanesMod(),   // TilePlanes.PlaneData [][] planes,
@@ -3724,7 +3749,7 @@ public class TileProcessor {
 					10.0,                 // double arrow_dark,
 					10.0);               // double arrow_white)
 			sdfa_instance.showArrays(plane_data, wh[0], wh[1], true, "shells2");
-// also test snap here			
+			// also test snap here			
 			int [] snap_surf =  st.snapSort(
 					st.cltPass3d.getDisparity(),                  // final double []  disparity,
 					st.cltPass3d.getStrength(),                   // final double []  strength,
@@ -3737,7 +3762,7 @@ public class TileProcessor {
 					clt_parameters.plSnapDispWeight, //  final double     snapDispWeight, //    =  .5,   // Maximal disparity diff. by weight product to snap to plane
 					clt_parameters.plSnapZeroMode,   // final int        snapZeroMode,   // Zero strength snap mode: 0: no special treatment, 1 - strongest, 2 - farthest
 					1); // final int        debugLevel)
-			
+
 			double [] snap_disparity =  st.snapDisparity(
 					st.cltPass3d.getDisparity(),                  // final double []  disparity,
 					st.cltPass3d.getStrength(),                   // final double []  strength,
@@ -3774,7 +3799,7 @@ public class TileProcessor {
 			snap_img[4] = snap_disparity_masked;
 			snap_img[6] = strength_masked;
 			sdfa_instance.showArrays(snap_img, tilesX, tilesY, true, "snap",snap_titles);
-			
+
 
 			boolean [] these_tiles = scan_prev.getSelected();
 
@@ -3821,7 +3846,7 @@ public class TileProcessor {
 					tilesY*clt_parameters.transform_size,
 					true,
 					"surf-neighbors");		
-			
+
 		}
 	}
 
@@ -4060,57 +4085,34 @@ public class TileProcessor {
 			final int         debugLevel)
 	{
 		CLTPass3d scan_prev = clt_3d_passes.get(clt_3d_passes.size() -1); // get last one
-///		double [] dbg_orig_disparity =  null;
-///		double [] dbg_with_super_disp = null;
-///		double [] dbg_outlayers =       null;
 		showDoubleFloatArrays sdfa_instance = null;
 		if (debugLevel > -1) sdfa_instance = new showDoubleFloatArrays(); // just for debugging?
 		boolean [] these_tiles = scan_prev.getSelected();
-//		double [] this_disparity     = scan_prev.getDisparity(); // returns a copy of the FPGA-generated disparity combined with the target one
-//		double [] this_strength =      scan_prev.getStrength(); // cloned, can be modified/ read back
 		double [] this_disparity     = scan_prev.getDisparity().clone(); // returns a copy of the FPGA-generated disparity combined with the target one
 		double [] this_strength =      scan_prev.getStrength().clone(); // cloned, can be modified/ read back
-		
+
 		scan_prev.fixNaNDisparity(
 				null, // border, // boolean [] select,   // which tiles to correct (null - all)
 				this_disparity, // scan_prev.getDisparity(), // double [] disparity,
 				this_strength); // scan_prev.getStrength()); // double [] strength)
-		
-///		dbg_orig_disparity = this_disparity.clone();
-		
+
+
 		SuperTiles st = scan_prev.getSuperTiles();
 
-//		if (clt_parameters.plSnapDispAny >= 0.0) {
-			double [] snap_disparity =  st.snapDisparity(
-					this_disparity,                  // final double []  disparity,
-					this_strength,                   // final double []  strength,
-					null, // this_selection,         // final boolean [] selection,       // can be null
-					clt_parameters.plDispNorm,       //  final double     dispNorm, // plDispNorm
-					clt_parameters.plSnapDispAny,    //final double     snapDispAny,    //    =  .2,   // Maximal (scaled by plDispNorm) disparity difference to snap to plane at any strength
-					clt_parameters.plSnapStrengthAny,//final double     snapStrengthAny,    //    =  .2,   //  Maximal strength to fit any distance (if does not fit otherwise - treat as zero strength
-					clt_parameters.plSnapNegAny,     //final double     snapNegAny,     //    =  .2,   // Maximal negative disparity difference from the best match
-					clt_parameters.plSnapDispMax,    //  final double     snapDispMax,    //    =  .5,   // Maximal (scaled by plDispNorm) disparity difference to snap to plane at low strength
-					clt_parameters.plSnapDispWeight, //  final double     snapDispWeight, //    =  .5,   // Maximal disparity diff. by weight product to snap to plane
-					clt_parameters.plSnapZeroMode,   // final int        snapZeroMode,   // Zero strength snap mode: 0: no special treatment, 1 - strongest, 2 - farthest
-					1); // final int        debugLevel)
-			System.arraycopy(snap_disparity,  0, this_disparity, 0, snap_disparity.length);
-//		}
-/**		
-// Replace weak outlayers
-		if (clt_parameters.replaceWeakOutlayers) {
-//			boolean[] outlayers = 
-					scan_prev.replaceWeakOutlayers(
-					null, // final boolean [] selection,
-					clt_parameters.outlayerStrength , //final double weakStrength,    // strength to be considered weak, subject to this replacement
-					clt_parameters.outlayerDiff, // final double maxDiff)
-					clt_parameters.outlayerDiffPos, // final double maxDiff)
-					clt_parameters.outlayerDiffNeg, // final double maxDiff)
-					0.5 * disparity_far,
-					2.0 * disparity_near,
-					debugLevel);
+		double [] snap_disparity =  st.snapDisparity(
+				this_disparity,                  // final double []  disparity,
+				this_strength,                   // final double []  strength,
+				null, // this_selection,         // final boolean [] selection,       // can be null
+				clt_parameters.plDispNorm,       //  final double     dispNorm, // plDispNorm
+				clt_parameters.plSnapDispAny,    //final double     snapDispAny,    //    =  .2,   // Maximal (scaled by plDispNorm) disparity difference to snap to plane at any strength
+				clt_parameters.plSnapStrengthAny,//final double     snapStrengthAny,    //    =  .2,   //  Maximal strength to fit any distance (if does not fit otherwise - treat as zero strength
+				clt_parameters.plSnapNegAny,     //final double     snapNegAny,     //    =  .2,   // Maximal negative disparity difference from the best match
+				clt_parameters.plSnapDispMax,    //  final double     snapDispMax,    //    =  .5,   // Maximal (scaled by plDispNorm) disparity difference to snap to plane at low strength
+				clt_parameters.plSnapDispWeight, //  final double     snapDispWeight, //    =  .5,   // Maximal disparity diff. by weight product to snap to plane
+				clt_parameters.plSnapZeroMode,   // final int        snapZeroMode,   // Zero strength snap mode: 0: no special treatment, 1 - strongest, 2 - farthest
+				1); // final int        debugLevel)
+		System.arraycopy(snap_disparity,  0, this_disparity, 0, snap_disparity.length);
 
-		}
-*/		
 		int [] snap_surf =  st.snapSort(
 				this_disparity,                  // final double []  disparity,
 				this_strength,                   // final double []  strength,
@@ -4123,8 +4125,6 @@ public class TileProcessor {
 				clt_parameters.plSnapDispWeight, //  final double     snapDispWeight, //    =  .5,   // Maximal disparity diff. by weight product to snap to plane
 				clt_parameters.plSnapZeroMode,   // final int        snapZeroMode,   // Zero strength snap mode: 0: no special treatment, 1 - strongest, 2 - farthest
 				1); // final int        debugLevel)
-		
-		//
 
 		DisparityProcessor dp = new DisparityProcessor(this, clt_parameters.transform_size * geometryCorrection.getScaleDzDx());
 		boolean [] grown = these_tiles.clone();
@@ -4134,7 +4134,7 @@ public class TileProcessor {
 				null);     // boolean [] prohibit)
 		boolean [] border = grown.clone();
 		for (int i = 0; i < border.length; i++) border[i] &= !these_tiles[i];
-		
+
 
 		int num_unmatched = st.getNumNotMatched(
 				snap_surf,
@@ -4148,12 +4148,11 @@ public class TileProcessor {
 		int [][] neib_surf =   new int [num_surf][];
 		double [][] dbg_neibs = new double [num_surf][];
 		for (int ns = 0 ; ns < num_surf; ns ++){
-			
 			neib_surf[ns] = st.getNeighbors( // creates neighbors mask from bitmask
 					these_tiles, // final boolean [] selected, // or null
 					snap_surf,   // final int []     snap_surf, // use this size - it matches image, not supertiles
 					ns);        // final int        surf_index)
-			
+
 			dbg_neibs[ns]= dp.dbgShowNeighbors(
 					these_tiles, // grown, // these_tiles, 
 					neib_surf[ns], // _orig, // int [] neighbors,
@@ -4161,22 +4160,20 @@ public class TileProcessor {
 					-1.0, // Double.NaN, // double bgnd,
 					10.0); // double fgnd)
 		}
-		
+
 		sdfa_instance.showArrays(dbg_neibs,
-				                 tilesX*clt_parameters.transform_size,
-				                 tilesY*clt_parameters.transform_size,
-				                 true,
-				                 "surf-neighbors");		
+				tilesX*clt_parameters.transform_size,
+				tilesY*clt_parameters.transform_size,
+				true,
+				"surf-neighbors");		
 		int numScans0 = 0;
 		for (int ns = 0 ; ns < num_surf; ns ++){
-//		for (int ns = 1 ; (ns < num_surf) && (ns < 2); ns ++){
+			//		for (int ns = 1 ; (ns < num_surf) && (ns < 2); ns ++){
 			boolean [] surface_selection = new boolean [these_tiles.length];
 			for (int i = 0; i < surface_selection.length; i++) {
 				surface_selection[i] = these_tiles[i] && (neib_surf[ns][i] >= 0);
 			}
-			
-			
-			
+
 			grown = surface_selection.clone();
 			growTiles(
 					2,          // grow tile selection by 1 over non-background tiles 1: 4 directions, 2 - 8 directions, 3 - 8 by 1, 4 by 1 more
@@ -4184,7 +4181,7 @@ public class TileProcessor {
 					null);     // boolean [] prohibit)
 			border = grown.clone();
 			for (int i = 0; i < border.length; i++) border[i] &= !these_tiles[i];
-			
+
 			int [] enum_clusters = enumerateClusters(
 					true, // boolean diag_en,
 					surface_selection); // grown); // these_tiles); // boolean [] tiles_src)
@@ -4196,13 +4193,76 @@ public class TileProcessor {
 					new_disparity0, // this_disparity, //  [] disparity_in, masked ok too
 					enum_clusters, // int []    clusters_in,
 					(ns == 1)?2:0); // debugLevel);
-			
+
 			if (debugLevel > -1){
 				System.out.println("thirdPassSetup(): surface=" + ns + " created "+ numScans0+ " FPGA passes.");
 			}
-			
+
 		}		
+	}
+	
+	
+	
+	public void thirdPassSetupSurf( // prepare tile tasks for the second pass based on the previous one(s)
+			EyesisCorrectionParameters.CLTParameters           clt_parameters,
+			double            disparity_far,    //
+			double            disparity_near,   // 
+			GeometryCorrection geometryCorrection,
+			final int         threadsMax,  // maximal number of threads to launch                         
+			final boolean     updateStatus,
+			final int         debugLevel)
+	{
+		CLTPass3d scan_prev = clt_3d_passes.get(clt_3d_passes.size() -1); // get last one
+		SuperTiles st = scan_prev.getSuperTiles();
+		TileSurface tileSurface = st.getTileSurface();
+		if (tileSurface == null){
+			System.out.println("tileSurface does not exist, aborting");
+			return;
+		}
 		
+		int num_surf_proc = clt_parameters.tsNumClust;
+		int num_surf = tileSurface.getSurfaceDataLength();
+		if (num_surf_proc > num_surf) num_surf_proc = num_surf;
+		int op = ImageDtt.setImgMask(0, 0xf);
+		op =     ImageDtt.setPairMask(op,0xf);
+		op =     ImageDtt.setForcedDisparity(op,true);
+		int numClusters = 0; // needed?
+		for (int ns = 0 ; ns < num_surf_proc; ns ++){
+			if (ns == 73){
+				System.out.println("thirdPassSetupSurf() ns = "+ns);
+			}
+			// return to original dimensions
+			double [][] disparityTask = new double [tilesY][tilesX];
+			int [][]    tile_op =       new int [tilesY][tilesX];
+			boolean [] borderTiles =    new boolean[tilesY*tilesX]; // to zero alpha in the images
+			for (int ty = 0; ty < tilesY; ty++) for (int tx = 0; tx <tilesX; tx++){
+				int indx =  tilesX * ty + tx;
+				disparityTask[ty][tx] = 0.0;
+				tile_op[ty][tx] =       0;
+				borderTiles[indx] =     false;
+			}
+			SurfaceData sd = tileSurface.getSurfaceData(ns);
+			boolean [] sel = sd.getSelected();
+			for (int i = 0; i < sel.length; i++) if (sel[i]){
+				int indx = sd.getImageIndex(i);
+				int tx = indx % tilesX; 
+				int ty = indx / tilesX;
+				tile_op[ty][tx] = op;
+				disparityTask[ty][tx] = sd.getDisparity(i);
+				borderTiles[indx] = sd.getBorder(i);
+			}
+			// Create FPGA task for this cluster
+
+			CLTPass3d scan_next =new CLTPass3d(this);
+			scan_next.disparity =     disparityTask;
+			scan_next.tile_op =       tile_op;
+			scan_next.border_tiles =  borderTiles;
+			clt_3d_passes.add(scan_next);
+			numClusters++;
+		}
+		if (debugLevel > -1){
+			System.out.println("thirdPassSetupSurf():  created "+ numClusters+ " FPGA passes.");
+		}
 	}
 	
 	

src/main/java/TileSurface.java

@@ -44,6 +44,7 @@ public class TileSurface {
 		private int threadsMax = 100;
 		private int [][] tileLayers = null;
 		private TileData [][] tileData = null;
+		private SurfaceData [] surfaceData = null; // result surfaces
 		
 		static int STAT_UNASSIGNED = 0; // index of number of unassigned tiles
 		static int STAT_ASSIGNED =   1; // index of number of assigned tiles
@@ -521,7 +522,7 @@ public class TileSurface {
 				final boolean                   divide_by_area,
 				final double                    scale_projection,
 				final double                    fraction_uni,
-				final TilePlanes.PlaneData [][] planes,
+				final TilePlanes.PlaneData [][] planes, // why plane here has weighth == 0.0 ?
 				final int                       debugLevel,
 				final int                       dbg_X,
 				final int                       dbg_Y)
@@ -543,7 +544,7 @@ public class TileSurface {
 								}
 								double [][][] disp_strength = new double [planes[nsTile].length][][];
 								for (int np = 0; np < disp_strength.length; np++){
-									if (planes[nsTile][np] != null){
+									if ((planes[nsTile][np] != null) && (planes[nsTile][np].getWeight() > 0.0)){ // disregard 0-weight planes
 										disp_strength[np] = planes[nsTile][np].getDoublePlaneDisparityStrength(
 												getWindow(),      // double [] window,
 												-1,              // int dir (-1 - center, 0- N, 1 - NE, .. 7 - NW
@@ -553,6 +554,7 @@ public class TileSurface {
 												fraction_uni, //				double    fraction_uni,
 												debugLevel-1);   // int       debugLevel)
 										// multiply disparities by strengths to calculate weighted averages
+										double [] disp0 = disp_strength[np][0].clone(); // to use where strength == 0 
 										for (int i = 0; i < disp_strength[np][1].length; i++){
 											disp_strength[np][0][i] *= disp_strength[np][1][i]; 
 										}
@@ -584,6 +586,8 @@ public class TileSurface {
 										for (int i = 0; i < disp_strength[np][1].length; i++){ 
 											if (disp_strength[np][1][i] > 0.0){
 												disp_strength[np][0][i] /= disp_strength[np][1][i]; 
+											} else {
+												disp_strength[np][0][i] = disp0[i]; // may have discontinuity
 											}
 										}
 										if (use_sel){ // zero out selection after averaging, apply to this tile
@@ -677,7 +681,8 @@ public class TileSurface {
 								}
 								corners[nsTile] = new int [planes[nsTile].length][][];
 								for (int np = 0; np < planes[nsTile].length; np++){
-									if (planes[nsTile][np] != null){
+//									if (planes[nsTile][np] != null){
+									if ((planes[nsTile][np] != null) && (planes[nsTile][np].getWeight() > 0.0)){ // disregard 0-weight planes
 										int [] neibs = planes[nsTile][np].getNeibBest();
 										if (neibs == null) {
 											System.out.println("getSupertilesTopology(), nsTile = "+nsTile+" neibs= null");
@@ -786,7 +791,8 @@ public class TileSurface {
 								}
 								meshes[nsTile] = new int [planes[nsTile].length][][][];
 								for (int np = 0; np < planes[nsTile].length; np++){
-									if (planes[nsTile][np] != null){
+//									if (planes[nsTile][np] != null){
+									if ((planes[nsTile][np] != null) && (planes[nsTile][np].getWeight() > 0.0)){ // disregard 0-weight planes
 										int [][] pre_mesh = new int [len_st2][2];
 										for (int i = 0; i < len_st2; i ++){
 											pre_mesh[i][0] = nsTile;
@@ -1037,7 +1043,9 @@ public class TileSurface {
 								}
 								int stileY = nsTile / stilesX;  
 								int stileX = nsTile % stilesX;
-								for (int np = 0; np < planes[nsTile].length; np++) if (planes[nsTile][np] != null){
+//								for (int np = 0; np < planes[nsTile].length; np++) if (planes[nsTile][np] != null){
+								for (int np = 0; np < planes[nsTile].length; np++) if ((planes[nsTile][np] != null) && (planes[nsTile][np].getWeight() > 0.0)){
+									
 									int [][][] src_mesh = lappingMeshes[nsTile][np];
 									double [][] disp_strength = fusedSupertilePlanes[nsTile][np];
 									TileData [] dual_mesh = new TileData [len2]; // full overlapping dual-sized mesh
@@ -1629,7 +1637,7 @@ public class TileSurface {
 				final int                       dbg_Y)
 		{
 
-			double [][][][] fused_planes = fuseSupertilePlanes (
+			double [][][][] fused_planes = fuseSupertilePlanes ( // some planes have zero weight that still participate - should they be eliminated?
 					use_sel,           // final boolean                   use_sel,
 					divide_by_area,    // final boolean                   divide_by_area,
 					scale_projection,  // final double                    scale_projection,
@@ -1755,6 +1763,26 @@ public class TileSurface {
 			return surfaceTilesX * (nTile /  imageTilesX) + (nTile %  imageTilesX);
 		}
 
+		/**
+		 * Convert from image tile index to the surface supertile index (surface tiles are all
+		 * full superTileSize),
+		 * TODO: update/remove if surface grid will be trimmed to fit image
+		 * Currently there are 324 tiles horizontally in the image and 328 in the surfaces
+		 * @param nTile image tile index in scan order
+		 * @return surface tile index in scan order
+		 */
+		public int getSurfaceSuperTileIndex(
+				int nTile)
+		{
+			// calculate index in tileData (has different dimensions - TODO: trim?
+//			int surfaceTilesX = stilesX * superTileSize;
+//			return surfaceTilesX * (nTile /  imageTilesX) + (nTile %  imageTilesX);
+			return stilesX * (nTile /  imageTilesX / superTileSize) + (nTile %  imageTilesX)/superTileSize;
+		}
+
+
+		
+		
 		/**
 		 * Convert from  surface tile index (surface tiles are all full superTileSize) to
 		 * the image tile index
@@ -1776,6 +1804,314 @@ public class TileSurface {
 			return imageTilesX * ty + tx;
 		}
 		
+		/**
+		 * Grow around surface tile, without going back. May still produce multi-layer result
+		 * that will need to be filtered. Used to fill large flat gaps
+		 * @param nsTile0
+		 * @param nl0
+		 * @return
+		 */
+		public boolean [][] growSimpleConnected(
+				int nsTile0,
+				int nl0)
+		{
+			int tilesX = stilesX * superTileSize;
+			if ((tileData[nsTile0] == null) || (tileData[nsTile0][nl0] == null)){
+				return null; // error
+			}
+			TileNeibs tnSurface = new TileNeibs(stilesX * superTileSize, stilesY * superTileSize);
+			boolean [][] simple_surf = new boolean [tileData.length][];
+			int x0 = nsTile0 % tilesX;
+			int y0 = nsTile0 / tilesX;
+			ArrayList<Point> wave_list = new ArrayList<Point>();
+			simple_surf[nsTile0] = new boolean [tileData[nsTile0].length];
+			simple_surf[nsTile0][nl0] = true;
+			wave_list.add(new Point(nsTile0, nl0));
+			while (!wave_list.isEmpty()){
+				Point p = wave_list.remove(0);
+				int nsTile = p.x;
+				int nl = p.y;
+				int x = nsTile % tilesX;
+				int y = nsTile / tilesX;
+				int dirs = 0xff;
+				if      (x > x0) dirs &= 0x1f;
+				else if (x < x0) dirs &= 0xf1;
+				if      (y > y0) dirs &= 0x7c;
+				else if (y < y0) dirs &= 0xc7;
+				int [] neibs = tileData[nsTile][nl].getNeighbors();
+				for (int dir = 0; dir < neibs.length; dir ++) {
+					int nl1 = neibs[dir];
+					if ((nl1 >= 0) && ((dirs & (1 << dir)) != 0)){
+						int nsTile1 = tnSurface.getNeibIndex(nsTile, dir);
+						if ((nsTile1 >= 0) && ((simple_surf[nsTile1] == null) || !simple_surf[nsTile1][nl1])){
+							if (simple_surf[nsTile1] == null) {
+								simple_surf[nsTile1] = new boolean [tileData[nsTile1].length];
+							}
+							simple_surf[nsTile1][nl1] = true;
+							wave_list.add(new Point(nsTile1, nl1));
+						}
+					}
+				}
+			}
+			return simple_surf;
+		}
+
+		public int [][] growSimpleIConnected(
+				int nsTile0,
+				int nl0)
+		{
+			int tilesX = stilesX * superTileSize;
+			if ((tileData[nsTile0] == null) || (tileData[nsTile0][nl0] == null)){
+				return null; // error
+			}
+			TileNeibs tnSurface = new TileNeibs(stilesX * superTileSize, stilesY * superTileSize);
+			int [][] simple_surf = new int [tileData.length][];
+			int x0 = nsTile0 % tilesX;
+			int y0 = nsTile0 / tilesX;
+			ArrayList<Point> wave_list = new ArrayList<Point>();
+			simple_surf[nsTile0] = new int [tileData[nsTile0].length];
+			simple_surf[nsTile0][nl0] = 1;
+			wave_list.add(new Point(nsTile0, nl0));
+			while (!wave_list.isEmpty()){
+				Point p = wave_list.remove(0);
+				int nsTile = p.x;
+				int nl = p.y;
+				int x = nsTile % tilesX;
+				int y = nsTile / tilesX;
+				int dirs = 0xff;
+				if      (x > x0) dirs &= 0x1f;
+				else if (x < x0) dirs &= 0xf1;
+				if      (y > y0) dirs &= 0x7c;
+				else if (y < y0) dirs &= 0xc7;
+				int [] neibs = tileData[nsTile][nl].getNeighbors();
+				for (int dir = 0; dir < neibs.length; dir ++) {
+					int nl1 = neibs[dir];
+					if ((nl1 >= 0) && ((dirs & (1 << dir)) != 0)){
+						int nsTile1 = tnSurface.getNeibIndex(nsTile, dir);
+						if ((nsTile1 >= 0) && ((simple_surf[nsTile1] == null) || (simple_surf[nsTile1][nl1] == 0))){
+							if (simple_surf[nsTile1] == null) {
+								simple_surf[nsTile1] = new int [tileData[nsTile1].length];
+							}
+							simple_surf[nsTile1][nl1] = simple_surf[nsTile][nl]+1;
+							wave_list.add(new Point(nsTile1, nl1));
+						}
+					}
+				}
+			}
+			return simple_surf;
+		}
+		
+		
+		
+		public double [][] getShowSimpleConnected(
+				int          nsTile0,
+				boolean [][] simple_surf)
+		{
+			int numSurfaces = 1;
+			for (int nsTile = 0; nsTile < simple_surf.length; nsTile++){
+				if (simple_surf[nsTile] != null){
+					int num_filled = 0;
+					for (int i = 0; i < simple_surf[nsTile].length; i++) if (simple_surf[nsTile][i]) num_filled++;
+					if (num_filled > numSurfaces) numSurfaces = num_filled; 
+				}
+			}
+			double [][] img_data = new double [numSurfaces][tileData.length];
+			for (int nsTile = 0; nsTile < tileData.length; nsTile++){
+				if (simple_surf[nsTile] != null){
+					int ns = 0;
+					for (int nl = 0; nl < simple_surf[nsTile].length; nl++) if (simple_surf[nsTile][nl]){
+						img_data[ns][nsTile] = tileData[nsTile][nl].getDisparity();
+						ns++;
+					}
+					for (; ns < numSurfaces; ns++) img_data[ns][nsTile] = Double.NaN;
+				} else {
+					for (int ns = 0; ns < numSurfaces; ns++) img_data[ns][nsTile] = Double.NaN;
+				}
+			}
+			// mark start point:
+			for (int ns = 0; ns < numSurfaces; ns++) img_data[ns][nsTile0] = Double.NaN;
+			return img_data;
+		}
+		public double [][] getShowSimpleConnected(
+				int          nsTile0,
+				int [][] simple_surf)
+		{
+			int numSurfaces = 1;
+			for (int nsTile = 0; nsTile < simple_surf.length; nsTile++){
+				if (simple_surf[nsTile] != null){
+					int num_filled = 0;
+					for (int i = 0; i < simple_surf[nsTile].length; i++) if (simple_surf[nsTile][i] !=0) num_filled++;
+					if (num_filled > numSurfaces) numSurfaces = num_filled; 
+				}
+			}
+			double [][] img_data = new double [numSurfaces][tileData.length];
+			for (int nsTile = 0; nsTile < tileData.length; nsTile++){
+				if (simple_surf[nsTile] != null){
+					int ns = 0;
+					for (int nl = 0; nl < simple_surf[nsTile].length; nl++) if (simple_surf[nsTile][nl] != 0){
+						img_data[ns][nsTile] = tileData[nsTile][nl].getDisparity();
+						ns++;
+					}
+					for (; ns < numSurfaces; ns++) img_data[ns][nsTile] = Double.NaN;
+				} else {
+					for (int ns = 0; ns < numSurfaces; ns++) img_data[ns][nsTile] = Double.NaN;
+				}
+			}
+			// mark start point:
+			for (int ns = 0; ns < numSurfaces; ns++) img_data[ns][nsTile0] = Double.NaN;
+			return img_data;
+		}
+
+		public double [][] getShowSimpleConnectedDistanceLayer(
+				boolean show_distance, // false - layer number
+				int [][] simple_surf)
+		{
+			int numSurfaces = 1;
+			for (int nsTile = 0; nsTile < simple_surf.length; nsTile++){
+				if (simple_surf[nsTile] != null){
+					int num_filled = 0;
+					for (int i = 0; i < simple_surf[nsTile].length; i++) if (simple_surf[nsTile][i] !=0) num_filled++;
+					if (num_filled > numSurfaces) numSurfaces = num_filled; 
+				}
+			}
+			double [][] img_data = new double [numSurfaces][tileData.length];
+			for (int nsTile = 0; nsTile < simple_surf.length; nsTile++){
+				if (simple_surf[nsTile] != null){
+					int ns = 0;
+					for (int nl = 0; nl < simple_surf[nsTile].length; nl++) if (simple_surf[nsTile][nl] != 0){
+						
+						img_data[ns][nsTile] = show_distance? simple_surf[nsTile][nl]: nl;
+						ns++;
+					}
+					for (; ns < numSurfaces; ns++) img_data[ns][nsTile] = Double.NaN;
+				} else {
+					for (int ns = 0; ns < numSurfaces; ns++) img_data[ns][nsTile] = Double.NaN;
+				}
+			}
+			return img_data;
+		}
+		
+
+		
+		public void testSimpleConnected(
+				int tileX,
+				int tileY)
+		{
+			testSimpleConnected(stilesX * superTileSize * tileY+ tileX);
+			printSurfaceConnections(
+					tileX - 7, // int left,
+					tileY - 7, // int top,
+					15,        // int width,
+					15,        // int height,
+					1,         // int extraH,
+					1);        // int extraV)
+			
+		}
+		
+		
+		public void testSimpleConnected(
+				int nsTile)
+		{
+			showDoubleFloatArrays sdfa_instance = new showDoubleFloatArrays();
+			double [][][] data =  new double [tileData[nsTile].length][][];
+			double [][][] dist =  new double [tileData[nsTile].length][][];
+			double [][][] layer = new double [tileData[nsTile].length][][];
+			int [][][] simple_surf = new int [tileData[nsTile].length][][];
+			int num_layers = 0;
+			for (int nl = 0; nl < tileData[nsTile].length; nl++) if (tileData[nsTile][nl] != null) {
+				simple_surf[nl] = growSimpleIConnected(nsTile, nl);
+				data[nl] =  getShowSimpleConnected(nsTile, simple_surf[nl]);
+				dist[nl] =  getShowSimpleConnectedDistanceLayer(true,simple_surf[nl]);
+				layer[nl] = getShowSimpleConnectedDistanceLayer(false,simple_surf[nl]);
+				num_layers += data[nl].length; 
+			}
+			double [][] img_data = new double [num_layers * 3][];
+			String [] titles = new String [num_layers * 3];
+			int indx = 0;
+			for (int nl = 0; nl < data.length; nl++) if (data[nl] != null) {
+				for (int i = 0; i < data[nl].length; i++){
+					titles[indx + 0 * num_layers] = ""+nl+((data[nl].length>1)?(":"+i):"");
+					titles[indx + 1 * num_layers] = "dist"+nl+((data[nl].length>1)?(":"+i):"");
+					titles[indx + 2 * num_layers] = "layer"+nl+((data[nl].length>1)?(":"+i):"");
+					img_data[indx + 0 * num_layers]= data[nl][i];
+					img_data[indx + 1 * num_layers]= dist[nl][i];
+					img_data[indx + 2 * num_layers]= layer[nl][i];
+					indx++;
+				}
+			}
+			sdfa_instance.showArrays(img_data,   stilesX * superTileSize,  stilesY * superTileSize, true, "simple_"+nsTile, titles);
+		}
+		
+		public void printSurfaceConnections(
+				int left,
+				int top,
+				int width,
+				int height,
+				int extraH,
+				int extraV)
+		{
+			int tilesX = stilesX * superTileSize;
+			int numSurfaces = 0;
+	 			for (int tileY = top; tileY < (top + height); tileY++) {
+				for (int tileX = left; tileX < (left + width); tileX++) {
+					int nsTile = tileY * tilesX + tileX;
+					if (tileData[nsTile] != null){
+						if (tileData[nsTile].length > numSurfaces) numSurfaces = tileData[nsTile].length; 
+					}
+				}
+			}
+ 			String hor_gap = "";
+ 			for (int i = 0; i < extraH; i++) hor_gap += " ";
+ 			String vert_gap = "";
+ 			for (int i = 0; i < extraV; i++) vert_gap += "\n";
+ 			for (int ns = 0; ns < numSurfaces; ns++){
+ 				System.out.println("\n ===== Surface number:"+ns+" =====");
+ 	 			for (int tileY = top; tileY < (top + height); tileY++) {
+ 	 				for (int l = 0; l < 3; l++){
+ 	 					for (int tileX = left; tileX < (left + width); tileX++) {
+ 	 						int nsTile = tileY * tilesX + tileX;
+ 	 						if ((tileData[nsTile]!= null) && (tileData[nsTile].length > ns) && (tileData[nsTile][ns] != null)){
+ 	 							int [] neibs = tileData[nsTile][ns].getNeighbors();
+ 	 							int master_st = tileData[nsTile][ns].getDbgNsTile();
+ 	 							switch (l){
+ 	 							case 0: System.out.print(" "+
+ 	 									((neibs[7] >= 0)?neibs[7]:" ")+
+ 	 									((neibs[0] >= 0)?neibs[0]:" ")+
+ 	 									((neibs[1] >= 0)?neibs[1]:" ")+" ");
+ 	 							break;
+ 	 							case 1: System.out.print(""+
+ 	 									((neibs[6] >= 0)?neibs[6]:" ")+
+ 	 									((master_st % stilesX) % 10)+
+ 	 									("*")+
+ 	 									((master_st / stilesX) % 10)+
+ 	 									((neibs[2] >= 0)?neibs[2]:" "));
+ 	 							break;
+ 	 							case 2: System.out.print(" "+
+ 	 									((neibs[5] >= 0)?neibs[5]:" ")+
+ 	 									((neibs[4] >= 0)?neibs[4]:" ")+
+ 	 									((neibs[3] >= 0)?neibs[3]:" ")+" ");
+ 	 							break;
+ 	 							}
+ 	 						} else {
+ 	 							switch (l){
+ 	 							case 0:
+ 	 							case 2: System.out.print("     "); break;
+ 	 							case 1: System.out.print("  .  "); break;
+ 	 							}
+ 	 						}
+ 	 						if (tileX < (left + width -1)){
+ 	 							System.out.print(hor_gap);
+ 	 						} else {
+ 	 							System.out.println();
+ 	 						}
+ 	 					} 	 					
+ 	 				}
+ 	 				if (tileY < (top + height -1)){
+ 	 					System.out.print(vert_gap);
+ 	 				}
+ 	 			}
+ 			}
+		}
 		
 		/**
 		 * Assign tiles to a certain disparity surface if there is only one surface candidate
@@ -1795,7 +2131,7 @@ public class TileSurface {
 		 * @param dbg_Y debug tile Y coordinate
 		 * @return 
 		 */
-		public int [] assignTilesToSingleCandidate(
+		public int [] assignTilesToSingleCandidate( // not used
 				final boolean       noEdge,
 				final double        maxDiff,
 				final double        minDiffOther,
@@ -2171,6 +2507,7 @@ public class TileSurface {
 		/**
 		 * Assign tiles to a certain disparity surface if there is only one surface candidate
 		 * @param noEdge do not assign tiles to the surface edges (can not add border later)
+		 * @param useCenter only assign outside of 8x8 center if no suitable alternative
 		 * @param maxDiff maximal (normalized) disparity difference
 		 * @param minDiffOther minimal disparity difference to closest 2-nd place candidate
 		 * @param minStrength minimal processed (floor subtracted) correlation strength of the candidate 
@@ -2197,6 +2534,7 @@ public class TileSurface {
 		
 		public int [] assignTilesToSurfaces(
 				final boolean       noEdge,
+				final boolean       useCenter,
 				final double        maxDiff,
 				final double        minDiffOther, // should be >= maxDiff
 				final double        minStrength,
@@ -2272,6 +2610,7 @@ public class TileSurface {
 									} else {
 										// calculate index in tileData (has different dimensions - TODO: trim?
 										int nSurfTile = getSurfaceTileIndex(nTile);
+										int nSurfSuperTile = getSurfaceSuperTileIndex(nTile);
 										if ((tileData[nSurfTile] == null) || (tileData[nSurfTile].length == 0)){
 											stats_all[numThread][NO_SURF] ++;
 											tileLayers[fml][nTile] = IMPOSSIBLE;
@@ -2280,7 +2619,11 @@ public class TileSurface {
 											int num_fit = 0;
 											int num_fit_other = 0;
 											int fit = -1;
+											int num_fit_center = 0;
+											int num_fit_other_center = 0;
+											int fit_center = -1;
 											boolean [] bad_surface = new boolean [tileData[nSurfTile].length]; // has less than 8 neighbors if controlled
+											boolean [] center_surface = new boolean [tileData[nSurfTile].length]; // is inside the center part
 											for (int ns = 0; ns < tileData[nSurfTile].length; ns++){
 												double surf_disp_diff = getNormDispFromSurface (
 														dispStrength[fml][0][nTile], // double disp_tile,
@@ -2303,9 +2646,29 @@ public class TileSurface {
 														if (Math.abs(surf_disp_diff) <= minDiffOther){
 															num_fit_other ++;
 														}
+														// separately calculate for center of the tiles getDbgNsTile
+														TileData td = tileData[nSurfTile][ns];
+														if (tileData[nSurfTile][ns].getDbgNsTile() == nSurfSuperTile){
+															center_surface[ns] = true;
+															if (Math.abs(surf_disp_diff) <= maxDiff){
+																fit_center = ns;   // no rating for fit "quality" here
+																num_fit_center ++;
+															}
+															if (Math.abs(surf_disp_diff) <= minDiffOther){
+																num_fit_other_center ++;
+															}
+														}
 													}
 												}
 											}
+											if (useCenter && (num_fit_center > 0)) { // use only center if possible
+												fit =           fit_center;
+												num_fit =       num_fit_center;
+												num_fit_other = num_fit_other_center;
+												for (int i = 0; i < center_surface.length; i++){
+													if (!center_surface[i]) bad_surface[i] = true;
+												}
+											}
 											if (num_fit < 1){
 												stats_all[numThread][TOO_FAR] ++;
 											} else if ((num_fit == 1) && (num_fit_other <= 1)){ // assign
@@ -2946,6 +3309,160 @@ public class TileSurface {
 			return grown_cluster;
 		}
 		
+		public void setSurfaceData(
+				final int [][] cluster_groups,
+				final int [][] clusters_pure, // each tile/layer belongs to a single pure cluster
+				final int      debugLevel,
+				final int      dbg_X,
+				final int      dbg_Y)
+		{
+			if (debugLevel >-1){
+				System.out.println("setSurfaceData(): groups: "+cluster_groups.length);
+			}
+			final SurfaceData [] sdata = new SurfaceData[cluster_groups.length];
+			final Thread[] threads = ImageDtt.newThreadArray(threadsMax);
+			final AtomicInteger ai = new AtomicInteger(0);
+			final int [] tilesWH = {imageTilesX, imageTilesY};
+			for (int ithread = 0; ithread < threads.length; ithread++) {
+				threads[ithread] = new Thread() {
+					public void run() {
+						for (int iGClust = ai.getAndIncrement(); iGClust < sdata.length; iGClust = ai.getAndIncrement()) {
+							int dl = ((debugLevel > -1) && (iGClust == dbg_X)) ? 1:0; 
+						
+							int [] window = getClusterBBox(
+									cluster_groups[iGClust], // final int []   nClust, // 1-based
+									1,                       // final int      border,
+									clusters_pure);          // final int [][] clusters)
+							
+							int [][] tile_indices =  getClusterBBoxIndices(
+									window,                  // final int [] window,
+									0);                      // border); // final int    border) // maybe 0, actual value just saves time
+							int [] img_tile_indices = tile_indices[0]; 
+							int [] surf_tile_indices = tile_indices[1];
+							boolean [][] combined = new boolean [img_tile_indices.length][];
+							if (dl > 0){
+								String dbg_s = "";
+								for (int indxClust = 0; indxClust < cluster_groups[iGClust].length; indxClust++){
+									dbg_s += " " +cluster_groups[iGClust][indxClust];
+								}
+								System.out.println("setSurfaceData(): iGClust = "+iGClust+" window = [ "+window[0]+", "+window[1]+", "+window[2]+", "+window[3]+" ],  clusters: ["+dbg_s+"]");
+							}
+							for (int indxClust = 0; indxClust< cluster_groups[iGClust].length; indxClust++){
+								if (dl > 0){
+									System.out.println("setSurfaceData():  cluster_groups["+iGClust+"]["+indxClust+"] = "+cluster_groups[iGClust][indxClust]);
+								}
+								boolean [][] extracted_cluster = extractCluster(
+										cluster_groups[iGClust][indxClust],   // final int      numToSplit, // 1-based
+										clusters_pure,       // final int [][] clusters,
+										window,              // final int []   window,
+										img_tile_indices,    // final int []   bbox_indices,				
+										debugLevel,          // final int      debugLevel,
+										dbg_X,               // final int      dbg_X,
+										dbg_Y);              // final int      dbg_Y)
+								for (int neTile = 0; neTile < combined.length; neTile++){
+									if (extracted_cluster[neTile] != null){
+										if (combined[neTile] == null){
+											combined[neTile] = extracted_cluster[neTile].clone();
+										} else {
+											for (int nl = 0; nl < extracted_cluster[neTile].length; nl++){
+												combined[neTile][nl] |= extracted_cluster[neTile][nl];
+											}
+										}
+									}
+								}
+							}
+							boolean [] pure_sel = new boolean [combined.length];
+							for (int neTile = 0; neTile < combined.length; neTile++){
+								if (combined[neTile] != null){
+									for (int nl = 0; nl < combined[neTile].length; nl++){
+										if (combined[neTile][nl]) {
+											pure_sel[neTile] = true;
+											break;
+										}
+									}									
+								}
+							}
+							boolean [][] grown_combined = growClusterOnce( // used from inside threads
+									combined,   // final boolean [][] cluster_local_pure,
+									window,            // final int []   window,
+									img_tile_indices,  // final int []   bbox_indices,				
+									surf_tile_indices, // final int []       surf_indices,				
+									0, // dl, // debugLevel,       // final int      debugLevel,
+									dbg_X,            // final int      dbg_X,
+									dbg_Y);           // final int      dbg_Y)
+							boolean [] selected = new boolean [combined.length];
+							boolean [] border  =  new boolean [combined.length];
+							for (int neTile = 0; neTile < grown_combined.length; neTile++){
+								if (grown_combined[neTile] != null){
+									for (int nl = 0; nl < grown_combined[neTile].length; nl++){
+										if (grown_combined[neTile][nl]) {
+											selected[neTile] = true;
+											if (!pure_sel[neTile]){
+												border[neTile] = true;
+											}
+											break;
+										}
+									}									
+								}
+							}
+							sdata[iGClust] = new 	SurfaceData(
+									tilesWH,
+									window);
+							sdata[iGClust].setSelected(selected);
+							sdata[iGClust].setBorder(border);
+							sdata[iGClust].setCluster(cluster_groups[iGClust]);
+							for (int neTile = 0; neTile < grown_combined.length; neTile++){
+								if (selected[neTile]){
+									for (int nl = 0; nl < grown_combined[neTile].length; nl++){
+										if (grown_combined[neTile][nl]){
+											sdata[iGClust].setDisparity(
+													neTile,
+													tileData[surf_tile_indices[neTile]][nl].getDisparity());
+										}
+									}
+								}
+							}
+							if (dl > 0){
+								System.out.println("setSurfaceData(): sdata ["+iGClust+"] prepared");
+							}
+						}
+					}
+				};
+			}		      
+			ImageDtt.startAndJoin(threads);
+			ArrayList<SurfaceData> sdata_list = new ArrayList<SurfaceData>();
+			for (int i = 0; i < sdata.length; i++){
+				sdata_list.add(sdata[i]);
+			}
+			Collections.sort(sdata_list, new Comparator<SurfaceData>() {
+				@Override
+				public int compare(SurfaceData lhs, SurfaceData rhs) {
+					// -1 - less than, 1 - greater than, 0 - equal, all inverted for descending
+					return (lhs.getSize(true) > rhs.getSize(true)) ? -1 : (lhs.getSize(true) < rhs.getSize(true)) ? 1 : 0;
+				}
+			});
+			this.surfaceData = 	sdata_list.toArray(new SurfaceData[0] );
+			if (debugLevel >-1){
+				if (debugLevel > 0){
+					for (int i = 0; i < this.surfaceData.length; i++) {
+						System.out.print(i+": "+this.surfaceData[i].getSize(true)+" ("+this.surfaceData[i].getSize(false)+") [");
+						for (int j = 0; j < this.surfaceData[i].getCluster().length; j++){
+							System.out.print(this.surfaceData[i].getCluster()[j] +", ");
+						}
+						System.out.println("]");
+					}
+				}
+				System.out.println("setSurfaceData(): DONE");
+			}
+ 		}
+		
+		public SurfaceData getSurfaceData(int indx){
+			return this.surfaceData[indx];
+		}
+
+		public int getSurfaceDataLength(){
+			return this.surfaceData.length;
+		}
 		
 		public int [][] mergeNoConflict(
 				final int [][] matchedGrown,
@@ -2955,7 +3472,7 @@ public class TileSurface {
 				final int      dbg_X,
 				final int      dbg_Y)
 		{
-			if (debugLevel >-1){
+			if (debugLevel > 0){
 				System.out.println("mergeNoConflict(): groups: "+matchedGrown.length);
 			}
 			final int num_grown = matchedGrown.length;
@@ -2976,7 +3493,7 @@ public class TileSurface {
 							} else { //>= 2 clusters
 								int nGClust = iGClust + 1;
 								int dl = ((debugLevel > -1) && (num_subs > 3)) ? 3: 0;
-								if (dl > 0){
+								if (dl > 3){
 									System.out.println("mergeNoConflict(): nGClust: "+nGClust+" num_subs= "+num_subs);
 								}
 								
@@ -2998,7 +3515,6 @@ public class TileSurface {
 								for (int nSub = 0; nSub < num_subs; nSub++){
 									subs_pure[nSub] = extractCluster(
 											matchedGrown[iGClust][nSub],   // final int      numToSplit, // 1-based
-											0, // border,     // final int      border,
 											clusters_pure,    // final int [][] clusters,
 											window,           // final int []   window,
 											img_tile_indices, // final int []   bbox_indices,				
@@ -3098,7 +3614,7 @@ public class TileSurface {
 									pre_merged_subs[iGClust][ng] = new int [groups_list.get(ng).size()];
 									int nc = 0;
 									for (Integer cl: groups_list.get(ng)){
-										pre_merged_subs[iGClust][ng][nc++] = cl;
+										pre_merged_subs[iGClust][ng][nc++] = matchedGrown[iGClust][cl]; // cl;
 									}
 								}
 							}
@@ -3123,7 +3639,7 @@ public class TileSurface {
 				}
 			}
 			
-			if (debugLevel > -1) {
+			if (debugLevel > 0) {
 				for (int i = 0; i < num_new_groups; i++){
 					System.out.print("mergeNoConflict(): "+ (i+1)+" [");
 					for (int j = 0; j < merged_subs[i].length; j++) {
@@ -3256,7 +3772,7 @@ public class TileSurface {
 				};
 			}		      
 			ImageDtt.startAndJoin(threads);
-			if (debugLevel > -1) {
+			if (debugLevel > 0) {
 				for (int i = 0; i < num_grown; i++){
 					System.out.print("matchPureGrown(): "+ (i+1)+" [");
 					for (int j = 0; j < grown_subs[i].length; j++) {
@@ -3271,17 +3787,17 @@ public class TileSurface {
 		
 		
 		public int [] getClusterBBox(
-				final int      numToSplit, // 1-based
+				final int      nClust, // 1-based
 				final int      border,
 				final int [][] clusters)
 		{
-			if (numToSplit == 74) { //177){
-				System.out.println("getClusterBBox() numToSplit="+numToSplit);
-			}
+//			if (numToSplit == 74) { //177){
+//				System.out.println("getClusterBBox() numToSplit="+numToSplit);
+//			}
 			int x_max = -1, y_max = -1, x_min = imageTilesX, y_min = imageTilesY;
 			for (int nTile = 0; nTile < clusters.length; nTile++) if (clusters[nTile] != null){
 				for (int nl = 0; nl < clusters[nTile].length; nl++){
-					if (clusters[nTile][nl] == numToSplit){
+					if (clusters[nTile][nl] == nClust){
 						int tX = nTile % imageTilesX;
 						
 						int tY = nTile / imageTilesX;
@@ -3299,6 +3815,41 @@ public class TileSurface {
 			int [] bbox = {x0, y0, width, height};
 			return bbox;
 		}		
+		public int [] getClusterBBox(
+				final int []   nClust, // 1-based
+				final int      border,
+				final int [][] clusters)
+		{
+			int x_max = -1, y_max = -1, x_min = imageTilesX, y_min = imageTilesY;
+			for (int nTile = 0; nTile < clusters.length; nTile++) if (clusters[nTile] != null){
+				for (int nl = 0; nl < clusters[nTile].length; nl++){
+					boolean sel = false;
+					for (int i = 0; i < nClust.length; i++){
+						if (clusters[nTile][nl] == nClust[i]){
+							sel = true;
+							break;
+						}
+					}
+					if (sel){
+						int tX = nTile % imageTilesX;
+						int tY = nTile / imageTilesX;
+						if (tX < x_min) x_min = tX;
+						if (tX > x_max) x_max = tX;
+						if (tY < y_min) y_min = tY;
+						if (tY > y_max) y_max = tY;
+					}
+				}
+			}
+			int x0 = x_min - border; // may be negative
+			int y0 = y_min - border; // may be negative
+			int width = x_max +1 + border - x0; 
+			int height = y_max +1 + border - y0;
+			int [] bbox = {x0, y0, width, height};
+			return bbox;
+		}		
+
+		
+		
 		
 		public int [][] getClusterBBoxIndices(
 				final int [] window,
@@ -3334,7 +3885,6 @@ public class TileSurface {
 		
 		public boolean [][] extractCluster(
 				final int      numToSplit, // 1-based
-				final int      border,
 				final int [][] clusters,
 				final int []   window,
 				final int []   bbox_indices,				
@@ -3407,7 +3957,6 @@ public class TileSurface {
 
 							boolean [][] extracted_cluster = extractCluster(
 									numToSplit + 1,   // final int      numToSplit, // 1-based
-									0, // border,       // final int      border,
 									clusters,     // final int [][] clusters,
 									window,       // final int []   window,
 									split_indices[numToSplit][0], // final int []   bbox_indices,				
@@ -3629,6 +4178,63 @@ public class TileSurface {
 			return split_cluster;
 		}
 		
+		public boolean [][] growSelectionOnce(
+//				int grow,
+				final boolean [][] sel_in,
+				final int          debugLevel,
+				final int          dbg_X,
+				final int          dbg_Y)
+		{
+			final boolean [][] selection =  new boolean [sel_in.length][];
+			final Thread[] threads = ImageDtt.newThreadArray(threadsMax);
+			final AtomicInteger ai = new AtomicInteger(0);
+			final TileNeibs tnImage =   new TileNeibs(imageTilesX, imageTilesY);
+			//			final TileNeibs tnSurface = new TileNeibs(stilesX * superTileSize, stilesY * superTileSize);
+			ai.set(0);
+			// make sure even empty tiles have surface selection arrays defined
+			for (int ithread = 0; ithread < threads.length; ithread++) {
+				threads[ithread] = new Thread() {
+					public void run() {
+						for (int nTile = ai.getAndIncrement(); nTile < sel_in.length; nTile = ai.getAndIncrement()) {
+							boolean not_empty = false;
+							boolean [] stack = null;
+							if (sel_in[nTile] != null) {
+								stack = sel_in[nTile].clone(); 
+							}
+							int nSurfTile = getSurfaceTileIndex(nTile);
+							if (tileData[nSurfTile] != null) {
+								if (stack == null) {
+									stack = new boolean[tileData[nSurfTile].length];
+								}
+								for (int nl = 0; nl < stack.length; nl++){
+									if (stack[nl]) {
+										not_empty = true;
+									} else {
+										int [] neibs = tileData[nSurfTile][nl].getNeighbors();
+										for (int dir = 0; dir < neibs.length; dir++){
+											if (neibs[dir] >= 0){
+												int nTile1 = tnImage.getNeibIndex(nTile, dir);
+												if ((nTile1 >= 0) && (sel_in[nTile1] != null) && sel_in[nTile1][neibs[dir]]){
+													stack[nl] = true;
+													not_empty = true;
+												}
+											}
+										}
+									}
+								}								
+							}
+							if (not_empty){
+								selection[nTile] = stack;
+							}
+						}
+					}
+				};
+			}		      
+			ImageDtt.startAndJoin(threads);
+			return selection;
+		}
+
+		
 		
 		public boolean [][] growSelection(
 				int grow,
@@ -3747,8 +4353,8 @@ public class TileSurface {
 					dbg_X,
 					dbg_Y);
 
-			boolean [][] grown_sel = growSelection(
-					2,              // int grow,
+			boolean [][] grown_sel = growSelectionOnce(
+//					2,              // int grow,
 					assigned_sel,   // final boolean [][] sel_in,
 					debugLevel,
 					dbg_X,
@@ -3772,7 +4378,14 @@ public class TileSurface {
 					debugLevel,
 					dbg_X,
 					dbg_Y);
- 					
+
+			setSurfaceData(
+					merged_no_conflict, // final int [][] cluster_groups,
+					clusters_pure, // final int [][] clusters_pure, // each tile/layer belongs to a single pure cluster
+					0, // final int      debugLevel,
+					dbg_X,
+					dbg_Y);
+			
 			System.out.println("mergeAndGrow() done");
 			
 

src/main/java/X3dOutput.java

@@ -110,7 +110,8 @@ public class X3dOutput {
 	}
 	
 	public void addCluster(
-			String url,
+			String      url,
+			String      id,
 			double [][] texCoord,
 			double [][] coordinate,
 			int [][] triangles)
@@ -156,6 +157,7 @@ public class X3dOutput {
 		
         Element el_shape =        x3dDoc.createElement("Shape");
 		el_Scene.appendChild(el_shape);
+		el_shape.setAttribute("id",id);
 		
         Element el_appearance =   x3dDoc.createElement("Appearance");
 		el_shape.appendChild(el_appearance);