removing foreground bridges

src/main/java/Conflicts.java

@@ -284,7 +284,7 @@ public class Conflicts {
 		final int stilesY =       (tilesY + superTileSize -1)/superTileSize;
 		final int nStiles =       stilesX * stilesY;
 		final int [][][] conflicts = new int [st.getPlanes().length][][];
-		final TileSurface.TileNeibs tnSurface = st.getTileSurface().new TileNeibs(stilesX, stilesY);
+		final TileNeibs tnSurface = new TileNeibs(stilesX, stilesY);
 		final Thread[] threads = ImageDtt.newThreadArray(st.getTileProcessor().threadsMax);
 		final AtomicInteger ai = new AtomicInteger(0);
 		for (int ithread = 0; ithread < threads.length; ithread++) {
@@ -343,7 +343,7 @@ public class Conflicts {
 			HashMap<Integer,Integer> replacement_tiles, // null is OK
 			int [][][]    replacement_neibs,               // null OK if  replacement_tiles == null
 			double [][][] replacement_val_weights,
-			TileSurface.TileNeibs tnSurface)
+			TileNeibs tnSurface)
 	{
 		TilePlanes.PlaneData [][] planes = st.getPlanes();
 		// generate conflicts if not provided
@@ -405,7 +405,7 @@ public class Conflicts {
 			HashMap<Integer,Integer> replacement_tiles, // null is OK
 			int [][][]    replacement_neibs,               // null OK if  replacement_tiles == null
 			ConnectionCosts connectionCosts,
-			TileSurface.TileNeibs tnSurface)
+			TileNeibs tnSurface)
 	{
 		TilePlanes.PlaneData [][] planes = st.getPlanes();
 		// generate conflicts if not provided
@@ -465,7 +465,7 @@ public class Conflicts {
 			int nsTile0,
 			HashMap<Integer,Integer> replacement_tiles, // null is OK - will use only planes data
 			int [][][] replacement_neibs,               // null OK if  replacement_tiles == null
-			TileSurface.TileNeibs tnSurface)
+			TileNeibs tnSurface)
 	{
 		TilePlanes.PlaneData [][] planes = st.getPlanes();
 		ArrayList<Conflict> conflicts_list= new ArrayList<Conflict>();

src/main/java/ConnectionCosts.java

@@ -28,7 +28,7 @@ import java.util.HashSet;
 public class ConnectionCosts {
 	TilePlanes.PlaneData [][] planes =  null;
 	boolean preferDisparity = false;
-	TileSurface.TileNeibs tnSurface;
+	TileNeibs tnSurface;
 	double         orthoWeight;
 	double         diagonalWeight;
 	double         starPwr;    // Divide cost by number of connections to this power
@@ -57,7 +57,7 @@ public class ConnectionCosts {
 			double         starValPwr, //  Raise value of each tile before averaging
 			int            steps,
 			TilePlanes.PlaneData [][] planes,
-			TileSurface.TileNeibs tnSurface,
+			TileNeibs tnSurface,
 			boolean preferDisparity)
 	{
 		this.planes =          planes;
@@ -481,7 +481,7 @@ public class ConnectionCosts {
 			double diagonalWeight,
 			double starPwr,    // Divide cost by number of connections to this power
 			double   starWeightPwr,    // Use this power of tile weight when calculating connection cost
-			TileSurface.TileNeibs tnSurface,
+			TileNeibs tnSurface,
 			boolean preferDisparity,
 			int    debugLevel)
 	{
@@ -506,7 +506,7 @@ public class ConnectionCosts {
 			double diagonalWeight,
 			double starPwr,    // Divide cost by number of connections to this power
 			double   starWeightPwr,    // Use this power of tile weight when calculating connection cost
-			TileSurface.TileNeibs tnSurface,
+			TileNeibs tnSurface,
 			boolean preferDisparity,
 			int    debugLevel)
 	{
@@ -566,7 +566,7 @@ public class ConnectionCosts {
 			double   diagonalWeight,
 			double   starPwr,    // Divide cost by number of connections to this power
 			double   starWeightPwr,    // Use this power of tile weight when calculating connection cost
-			TileSurface.TileNeibs tnSurface,
+			TileNeibs tnSurface,
 			boolean preferDisparity,
 			int     debugLevel)
 	{
@@ -594,7 +594,7 @@ public class ConnectionCosts {
 			double   diagonalWeight,
 			double   starPwr,    // Divide cost by number of connections to this power
 			double   starWeightPwr,    // Use this power of tile weight when calculating connection cost
-			TileSurface.TileNeibs tnSurface,
+			TileNeibs tnSurface,
 			boolean preferDisparity,
 			int     debugLevel)
 	{

src/main/java/LinkPlanes.java

@@ -813,6 +813,17 @@ public class LinkPlanes {
 		ImageDtt.startAndJoin(threads);
 	}
 
+	/**
+	 * Merge the supertile planes with agreeing neighbors, non-exclusively (no considering other planes
+	 * of the same supertile. Start with the best fit, then goes to lower quality, until the individual
+	 * merge quality falls below scaled quality of the best, pre-set minimum or the merged plane becomes
+	 * too thick
+	 * Separately calculates merged weighted plane and with equal weights of the neighbors
+	 * @param planes array of plane instances for the same supertile
+	 * @param debugLevel
+	 * @param dbg_X
+	 * @param dbg_Y
+	 */
 	public void setNonExclusive(
 			final TilePlanes.PlaneData [][] planes,
 			final int     debugLevel,
@@ -834,7 +845,7 @@ public class LinkPlanes {
 		//				final int debug_stile = 9 * stilesX + 26;
 		final int debug_stile = dbg_Y * stilesX + dbg_X;
 
-		final Thread[] threads = ImageDtt.newThreadArray(st.tileProcessor.threadsMax);
+		final Thread[] threads = ImageDtt.newThreadArray((debugLevel > 1)? 1 : st.tileProcessor.threadsMax);
 		final AtomicInteger ai = new AtomicInteger(0);
 
 		for (int ithread = 0; ithread < threads.length; ithread++) {
@@ -2426,7 +2437,7 @@ public class LinkPlanes {
 		final int stilesX =       (tilesX + superTileSize -1)/superTileSize;  
 		final int stilesY =       (tilesY + superTileSize -1)/superTileSize;
 		final int nStiles =       stilesX * stilesY; 
-		final TileSurface.TileNeibs tnSurface = st.tileSurface.new TileNeibs(stilesX, stilesY);
+		final TileNeibs tnSurface = new TileNeibs(stilesX, stilesY);
 		final Thread[] threads = ImageDtt.newThreadArray(st.tileProcessor.threadsMax);
 		final AtomicInteger ai = new AtomicInteger(0);
 		for (int ithread = 0; ithread < threads.length; ithread++) {
@@ -2441,7 +2452,7 @@ public class LinkPlanes {
 							starValPwr,     //double          starValPwr, //  Raise value of each tile before averaging
 							steps,          // int            steps,
 							planes,         // TilePlanes.PlaneData [][] planes,
-							tnSurface,      // TileSurface.TileNeibs tnSurface,
+							tnSurface,      // TileNeibs tnSurface,
 							preferDisparity); // boolean preferDisparity)
 					int [] mod_supertiles = new int[1];
 					for (int nsTile = ai.getAndIncrement(); nsTile < nStiles; nsTile = ai.getAndIncrement()) {
@@ -2484,7 +2495,7 @@ public class LinkPlanes {
 		//				final int tileSize =      tileProcessor.getTileSize();
 		final int stilesX =       (tilesX + superTileSize -1)/superTileSize;  
 		final int stilesY =       (tilesY + superTileSize -1)/superTileSize;
-		final TileSurface.TileNeibs tnSurface = st.tileSurface.new TileNeibs(stilesX, stilesY);
+		final TileNeibs tnSurface = new TileNeibs(stilesX, stilesY);
 		final Thread[] threads = ImageDtt.newThreadArray(st.tileProcessor.threadsMax);
 		final AtomicInteger ai = new AtomicInteger(0);
 		for (int ithread = 0; ithread < threads.length; ithread++) {
@@ -2499,7 +2510,7 @@ public class LinkPlanes {
 							starValPwr,      //double     starValPwr, //  Raise value of each tile before averaging
 							steps,          // int            steps,
 							planes,         // TilePlanes.PlaneData [][] planes,
-							tnSurface,      // TileSurface.TileNeibs tnSurface,
+							tnSurface,      // TileNeibs tnSurface,
 							preferDisparity); // boolean preferDisparity)
 					int [] supertiles = new int[1];
 					for (int isTile = ai.getAndIncrement(); isTile < mod_supertiles.length; isTile = ai.getAndIncrement()) {

src/main/java/SuperTiles.java

@@ -1025,147 +1025,6 @@ public class SuperTiles{
 		return bgTileDispStrength;
 	}
 
-	public void processPlanes(
-
-			final boolean [] selected, // or null
-			final double     min_disp,
-			final boolean    invert_disp, // use 1/disparity
-			final double     plDispNorm,
-			final int        debugLevel)
-	{
-
-		final int tilesX =        tileProcessor.getTilesX();
-		final int tilesY =        tileProcessor.getTilesY();
-		final int superTileSize = tileProcessor.getSuperTileSize();
-		final int tileSize =      tileProcessor.getTileSize();
-		final double  [] disparity = cltPass3d.getDisparity();
-		final double  [] strength =  cltPass3d.getStrength();
-		final int stilesX = (tilesX + superTileSize -1)/superTileSize;  
-		final int stilesY = (tilesY + superTileSize -1)/superTileSize;
-		final int nStiles = stilesX * stilesY; 
-		final Thread[] threads = ImageDtt.newThreadArray(tileProcessor.threadsMax);
-		final AtomicInteger ai = new AtomicInteger(0);
-		final int st_start = -superTileSize/2;
-		final int superTileSize2 = 2 * superTileSize;
-		final double [][] lapWeight = getLapWeights();
-		final int len2 = superTileSize2*superTileSize2;
-		final double []  double_zero =  new double  [len2];
-		final boolean [] boolean_zero = new boolean [len2];
-
-		final int debug_stile = 18 * stilesX + 25; 
-
-
-		for (int ithread = 0; ithread < threads.length; ithread++) {
-			threads[ithread] = new Thread() {
-				public void run() {
-					TilePlanes tpl = new TilePlanes(tileSize,superTileSize);
-					double [] stDisparity = new double [superTileSize2*superTileSize2];
-					double [] stStrength =  new double [superTileSize2*superTileSize2];
-					boolean [] stSel = new boolean [superTileSize2*superTileSize2];
-					for (int nsTile = ai.getAndIncrement(); nsTile < nStiles; nsTile = ai.getAndIncrement()) {
-						int stileY = nsTile / stilesX;  
-						int stileX = nsTile % stilesX;
-						double sw = 0.0; // sum weights
-						double [] hist = new double [numBins];
-						int tY0 = stileY * superTileSize + st_start;
-						int tX0 = stileX * superTileSize + st_start;
-						System.arraycopy(double_zero,  0, stDisparity, 0, len2);
-						System.arraycopy(double_zero,  0, stStrength,  0, len2);
-						System.arraycopy(boolean_zero, 0, stSel,       0, len2);
-						for (int tY = 0; tY < superTileSize2; tY++){
-							int tileY = tY0 +tY;
-							if ((tileY >= 0) && (tileY < tilesY)) {
-								for (int tX = 0; tX < superTileSize2; tX++){
-									int tileX = tX0 +tX;
-									if ((tileX >= 0) && (tileX < tilesX)) {
-										int indx = tileY*tilesX + tileX;
-										double d = disparity[indx];
-										if (!Double.isNaN(d) && (d >= min_disp) &&((selected == null) || selected[indx])){
-											if (invert_disp){
-												d = 1.0/d;
-											}
-											double w = strength[indx] - strength_floor;
-											if (w > 0.0){
-												if (strength_pow != 1.0) w = Math.pow(w, strength_pow);
-												w *= lapWeight[tY][tX];
-												int indx_out =  tY * superTileSize2 + tX;
-												stDisparity[indx_out] = d;
-												stStrength[indx_out] =  w;
-												stSel[indx_out] =       true;
-												sw +=w;
-											}
-										}
-									}
-								}
-							}
-						}
-						if (sw >0){
-							//									int dl = ((nsTile >= debug_stile-1) && (nsTile <= debug_stile+1) ) ? 1 : 0;
-							//									int dl = (stileY == 17) ? 1 : 0;
-							int dl = (stileY >= 0) ? 1 : 0;
-							double [][][] rslt = tpl.getCovar(
-									stDisparity,
-									stStrength,
-									stSel,
-									plDispNorm,
-									0); // dl); // debugLevel);
-							if (dl > 0) {
-								int       numPoints =  (int) rslt[2][0][2];
-								double    kz =   rslt[2][0][1];
-								double    swc =  rslt[2][0][0];
-								double [] szxy = rslt[2][1];
-								double [][] eig_val =  rslt[0];
-								double [][] eig_vect = rslt[1];
-								if (swc > 1.0) {
-									System.out.println("Processing planes, nsTile="+nsTile+", stileX="+stileX+", stileY="+stileY+", numPoints="+numPoints+
-											", kz = "+kz+", sw = "+sw+ ", swc = "+swc+ ", center=["+szxy[0]+","+szxy[1]+","+szxy[2]+"]"+
-											", eig_val = {"+eig_val[0][0]+","+eig_val[1][1]+","+eig_val[2][2]+"}"+
-											", eig_vect[0] = {"+eig_vect[0][0]+","+eig_vect[1][0]+","+eig_vect[2][0]+"}");
-								}
-							}
-
-							/*
-		double [][][] rslt = {
-				eig.getD().getArray(),
-				eig.getV().getArray(),
-				{
-					{sw,kz},
-					{swz, swx, swy}}};
-		return rslt;
-
-							 */
-
-
-							if (dl > 1) {
-								System.out.println("Processing planes with average disparity");
-								double ss = 0.0, sd = 0.0;
-								for (int i = 0; i < stDisparity.length; i++){
-									if (!Double.isNaN(stDisparity[i]) && stSel[i]){
-										ss += stStrength[i];
-										sd += stStrength[i] * stDisparity[i];
-									}
-								}
-								if (ss > 0) {
-									sd /= ss; 
-								}
-								for (int i = 0; i < stDisparity.length; i++){
-									stDisparity[i] = sd;
-								}
-								tpl.getCovar(
-										stDisparity,
-										stStrength,
-										stSel,
-										plDispNorm,
-										dl); // debugLevel);
-							}
-						}
-					}
-				}
-			};
-		}		      
-		ImageDtt.startAndJoin(threads);
-	}
-
 
 	public double corrMaxEigen(
 			double maxEigen,
@@ -1766,12 +1625,37 @@ public class SuperTiles{
 		return plane_selections;
 	}
 
-	// use both horizontal and const disparity tiles to create tile clusters
-	// Add max_diff (maximal disparity difference while extracting initial tile selection) and max_tries (2..3) parameters
-	
-	// Add separate method to create + remove outliers from all planes (2 different ones)?
-	// TODO later re-assign pixels according to existing plane parameters
-	// Sort plane data by center (plane or supertile) disparity
+	/**
+	 * Use both horizontal and const disparity tiles to create tile clusters
+	 *Add max_diff (maximal disparity difference while extracting initial tile selection) and max_tries (2..3) parameters
+	 * Add separate method to create + remove outliers from all planes (2 different ones)?
+	 * TODO later re-assign pixels according to existing plane parameters
+	 * Sort plane data by center (plane or supertile) disparity
+	 * 	
+	 * @param growSelection
+	 * @param stMeasSel
+	 * @param plDispNorm
+	 * @param plMinPoints
+	 * @param plPreferDisparity
+	 * @param geometryCorrection
+	 * @param correct_distortions
+	 * @param smplMode
+	 * @param smplSide
+	 * @param smplNum
+	 * @param smplRms
+	 * @param bin_blur_hor
+	 * @param bin_blur_vert
+	 * @param max_diff_hor
+	 * @param max_diff_vert
+	 * @param max_tries
+	 * @param smallDiff
+	 * @param highMix
+	 * @param world_hor
+	 * @param debugLevel
+	 * @param dbg_X
+	 * @param dbg_Y
+	 * @return
+	 */
 
 	public boolean [][][][]  initialDiscriminateTiles(
 			final int        growSelection,                     // grow initial selection before processing 
@@ -2058,15 +1942,36 @@ public class SuperTiles{
 		return planes_selections;
 	}
 	
-	public TilePlanes.PlaneData [][] createPlanesFromSelectionsOld(
-			final boolean [][][][] plane_selections, //  = new boolean [nStiles][][][]; // num_tiles
-			final double  [][][][] disp_strength,			
+	/**
+	 * Re-assign tiles to the planes according to fitted planes. Scans each known plane parallel with specified
+	 * size steps and specified number of steps, maximizing overall fit quality.
+	 * Total weight of the fitted tiles? RMS? 
+	 * @param planes per-supertile, per plane array of plane data instances. Should have
+	 *        nonexclusiveStar and nonexclusiveStarEq calculated
+	 * @param plDispNorm - normalization of the measured disparity precision - closer objects will have disparity
+	 *  difference proportionally reduced
+	 * @param plMinPoints minimal number of points in a plane
+	 * @param plPreferDisparity kept for historical reasons - when true, select disparity-most vector even if it has higher eigenvalue
+	 * @param geometryCorrection GeometryCorrection instance to use
+	 * @param correct_distortions correct geometrical distortions when converting to/from world coordinates 
+	 * @param smplMode sample mode
+	 * @param smplSide sample side for averaging/filtering tile disparities
+	 * @param smplNum number of best samples used fro averaging
+	 * @param smplRms maximal sample disparity rms to consider sample valid
+	 * @param max_disp_diff maximal disparity difference from the planes 
+	 * @param disp_step disparity scan step when searching for optimal fit
+	 * @param steps_per_plane number of disparity steps to try around initial plane
+	 * @param mode which of the calculated planes to use - 0 - weighted, 1 - equalized, 2 heaviest of the two, weighted combination
+	 * @param debugLevel debug level
+	 * @param dbg_X supertile X for elevated debug level
+	 * @param dbg_Y supertile X for elevated debug level
+	 * @return per-tile, per plane, per measurement layer (type of correlation - combo, 4, hor, vert), per tile selections of
+	 * filtered tiles
+	 */
+	public boolean [][][][]  refineDiscriminateTiles(
+			final TilePlanes.PlaneData [][] planes, 
 			final double     plDispNorm,
 			final int        plMinPoints, //          =     5;  // Minimal number of points for plane detection
-			final double     plTargetEigen, //        =   0.1;  // Remove outliers until main axis eigenvalue (possibly scaled by plDispNorm) gets below
-			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,
@@ -2076,182 +1981,20 @@ public class SuperTiles{
 			final int        smplNum, //         = 3;      // Number after removing worst
 			final double     smplRms, //         = 0.1;    // Maximal RMS of the remaining tiles in a sample
 			
+			final double     max_disp_diff,    // maximal disparity difference from the plane to consider tile 
+			final double     disp_step,        // parallel move known planes around original know value for the best overall fit 
+			final int        steps_per_plane,  // number of steps for each plane to search for the best fit
+			final int        mode, // 0 - weighted, 1 - equalized, 2 - best, 3 - combined
+
 			final int        debugLevel,
 			final int        dbg_X,
 			final int        dbg_Y)
 	{
-		final int tilesX =        tileProcessor.getTilesX();
-		final int tilesY =        tileProcessor.getTilesY();
-		final int superTileSize = tileProcessor.getSuperTileSize();
-		final int tileSize =      tileProcessor.getTileSize();
-
-		final int stilesX = (tilesX + superTileSize -1)/superTileSize;  
-		final int stilesY = (tilesY + superTileSize -1)/superTileSize;
-		final int nStiles = stilesX * stilesY; 
-		final Thread[] threads = ImageDtt.newThreadArray(tileProcessor.threadsMax);
-		final AtomicInteger ai = new AtomicInteger(0);
-		final TilePlanes.PlaneData [][] result_planes = new TilePlanes.PlaneData[nStiles][];
-		//		this.planes = new TilePlanes.PlaneData[nStiles][];
-		final int debug_stile = (debugLevel > -1)? (dbg_Y * stilesX + dbg_X):-1;
-		// TODO: Remove when promoting PlaneData
-		final TilePlanes tpl = new TilePlanes(tileSize,superTileSize, geometryCorrection);
-
-		for (int ithread = 0; ithread < threads.length; ithread++) {
-			threads[ithread] = new Thread() {
-				public void run() {
-					for (int nsTile = ai.getAndIncrement(); nsTile < nStiles; nsTile = ai.getAndIncrement()) {
-						if (nsTile == debug_stile){
-							System.out.println("createPlanesFromSelections(): nsTile="+nsTile);
-						}
-						if (plane_selections[nsTile] != null) {
-							
-							
-							
-							int stileY = nsTile / stilesX;  
-							int stileX = nsTile % stilesX;
-							int [] sTiles = {stileX, stileY};
-							int dl =  (nsTile == debug_stile) ? 3 : 0;
+		// create a list of usable planes according to the mode
 		
 		
 		
-							result_planes[nsTile] = null;
-							// first make a plane from all tiles
-							ArrayList<TilePlanes.PlaneData> st_planes = new ArrayList<TilePlanes.PlaneData>();
-							TilePlanes.PlaneData pd0 = tpl.new  PlaneData (
-									sTiles, // int [] sTileXY, 
-									tileSize, // int tileSize,
-									geometryCorrection, // GeometryCorrection   geometryCorrection,
-									correct_distortions,
-									measuredLayers,     // MeasuredLayers measuredLayers,
-									plPreferDisparity);   // boolean preferDisparity)
-
-							
-							// iterate through all plane selections
-							for (int ps = 0; ps < plane_selections[nsTile].length; ps++) {
-								TilePlanes.PlaneData pd = pd0.clone(); 
-								boolean OK = (pd.getPlaneFromMeas(
-										plane_selections[nsTile][ps], // tile_sel,       // boolean [][] tile_sel, // null - do not use, {} use all (will be modified)
-										disp_strength[nsTile],
-										Double.NaN,    // double       disp_far, // minimal disparity to select (or NaN)
-										Double.NaN,    // double       disp_near, // maximal disparity to select (or NaN)
-										plDispNorm,    // 0.0,            // plDispNorm,  // double       dispNorm,   //  Normalize disparities to the average if above
-										0.0,            // double       min_weight,
-										plMinPoints,    // int          min_tiles,
-										strength_floor, // 
-										strength_pow,   // double       strength_pow,
-										// update !
-										smplMode,
-										smplSide,
-										smplNum,
-										smplRms,
-										dl) != null);            // int          debugLevel)
-								if (OK) {
-									if (dl > 0) {
-										if (pd.getWeight() > 1.0) {
-											System.out.println("Processing subplane["+nsTile+"]["+ps+"]"+
-													", stileX="+stileX+
-													", stileY="+stileY+
-													", numPoints="+ pd.getNumPoints()+
-													", swc = "+pd.getWeight()+
-													", center=["+pd.getZxy()[0]+","+pd.getZxy()[1]+","+pd.getZxy()[2]+"]"+
-													", eig_val = {"+pd.getValues()[0]+","+pd.getValues()[1]+","+pd.getValues()[2]+"}"+
-													", eig_vect[0] = {"+pd.getVector()[0]+","+pd.getVector()[1]+","+pd.getVector()[2]+"}");
-										}
-									}
-									// now try to remove outliers
-									int max_outliers = (int) Math.round(pd.getNumPoints() * plFractOutliers);
-									if (max_outliers > plMaxOutliers) max_outliers = plMaxOutliers;
-									double targetV = plTargetEigen;
-									double z0 = pd.getZxy()[0];
-									if ((plDispNorm > 0.0) && (z0 > plDispNorm)) {
-										double dd = (plDispNorm + z0)/ plDispNorm; // > 1
-										targetV *= dd * dd; // > original
-									}
-									if (pd.getValues()[0] > targetV) {
-										OK = pd.removeOutliers( // getPlaneFromMeas should already have run
-												disp_strength[nsTile], 
-												targetV,      // double     targetEigen, // target eigenvalue for primary axis (is disparity-dependent, so is non-constant)
-												max_outliers, // int        maxRemoved,  // maximal number of tiles to remove (not a constant)
-												dl); // int        debugLevel)
-										if (!OK) {
-											continue;
-										}
-										if (dl > 0) {
-											if (pd.getWeight() > 1.0) {
-												System.out.println("Removed outliers["+nsTile+"]["+ps+"]"+
-														", stileX="+stileX+
-														", stileY="+stileY+
-														", numPoints="+ pd.getNumPoints()+
-														", swc = "+pd.getWeight()+
-														", center=["+pd.getZxy()[0]+","+pd.getZxy()[1]+","+pd.getZxy()[2]+"]"+
-														", eig_val = {"+pd.getValues()[0]+","+pd.getValues()[1]+","+pd.getValues()[2]+"}"+
-														", eig_vect[0] = {"+pd.getVector()[0]+","+pd.getVector()[1]+","+pd.getVector()[2]+"}");
-											}
-										}
-									}
-									double [] norm_xyz = pd.getWorldXYZ(
-											correct_distortions);
-									st_planes.add(pd);
-									if (dl > 0) {
-										System.out.println("World normal["+nsTile+"]["+ps+"] = {"+
-												norm_xyz[0]+", "+norm_xyz[1]+", "+norm_xyz[2]+"}");
-
-									}
-									// calculate the world planes too									
-//									if (debugLevel > -1){
-										pd.getWorldPlaneFromMeas(
-												plane_selections[nsTile][ps], // tile_sel,       // boolean [][] tile_sel, // null - do not use, {} use all (will be modified)
-												disp_strength[nsTile],
-												Double.NaN,    // double       disp_far, // minimal disparity to select (or NaN)
-												Double.NaN,    // double       disp_near, // maximal disparity to select (or NaN)
-												plDispNorm,    // 0.0,            // plDispNorm,  // double       dispNorm,   //  Normalize disparities to the average if above
-												0.0,            // double       min_weight,
-												plMinPoints,    // int          min_tiles,
-												strength_floor, // 
-												strength_pow,   // double       strength_pow,
-												// update !
-												smplMode,
-												smplSide,
-												smplNum,
-												smplRms,
-												dl);
-//									}
-								}
-							}
-							if (st_planes.size() > 0){
-								// sort planes by increasing disparity (tile center or plane center ? ) Using plane center
-								Collections.sort(st_planes, new Comparator<TilePlanes.PlaneData>() {
-									@Override
-									public int compare(TilePlanes.PlaneData lhs, TilePlanes.PlaneData rhs) {
-										// -1 - less than, 1 - greater than, 0 - equal
-										return (rhs.getZxy()[0] > lhs.getZxy()[0]) ? -1 : (rhs.getZxy()[0] < lhs.getZxy()[0] ) ? 1 : 0;
-									}
-								});
-								if (LOWEST_PLANE(2) > 0) st_planes.add(0, st_planes.get(0)); // insert dummy at pos 0;
-								result_planes[nsTile] = st_planes.toArray(new TilePlanes.PlaneData[0] );
-								if (LOWEST_PLANE(2) > 0) result_planes[nsTile][0] = null; // remove dummy
-								if (dl >0){
-									System.out.println("createPlanesFromSelections(): nsTile="+nsTile);
-								}
-							}
-							if (dl > 2) {
-								String [] dbg_titles = showSupertileSeparationTitles( 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_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);
-							}
-							
-							
-							
-						}
-					}
-				}
-			};
-		}		      
-		ImageDtt.startAndJoin(threads);
-		return result_planes;		
+		return null;
 	}
 	
 	public TilePlanes.PlaneData [][] createPlanesFromSelections(
@@ -2284,7 +2027,9 @@ public class SuperTiles{
 		final int stilesX = (tilesX + superTileSize -1)/superTileSize;  
 		final int stilesY = (tilesY + superTileSize -1)/superTileSize;
 		final int nStiles = stilesX * stilesY; 
-		final Thread[] threads = ImageDtt.newThreadArray(tileProcessor.threadsMax);
+//		final Thread[] threads = ImageDtt.newThreadArray(tileProcessor.threadsMax);
+		final Thread[] threads = ImageDtt.newThreadArray((debugLevel > 1)? 1 : tileProcessor.threadsMax);
+		
 		final AtomicInteger ai = new AtomicInteger(0);
 		final TilePlanes.PlaneData [][] result_planes = new TilePlanes.PlaneData[nStiles][];
 		//		this.planes = new TilePlanes.PlaneData[nStiles][];
@@ -2336,6 +2081,52 @@ public class SuperTiles{
 											dl);                       // int          debugLevel);
 
 							if ((st_planes != null) && (!st_planes.isEmpty())){
+								if (dl > 0){
+									System.out.println("======= createPlanesFromSelections(): nsTile="+nsTile+" detecting bridges ==========");
+								}
+								boolean [][][] split_sels = pd0.filterBridges(
+										st_planes, // ArrayList<PlaneData> tilePlanes,
+										3, // int max_grow, // make configurable?
+										3, // int max_grow_far,
+										dl); // int debugLevel)
+								if (split_sels !=null){
+									if (dl > -1){
+										System.out.println("======= createPlanesFromSelections(): nsTile="+nsTile+" removing bridges ==========");
+									}
+									if (dl > 2) {
+										String [] dbg_titles = showSupertileSeparationTitles( 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, "pre_bridge_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, "pre_bridge_world-"+nsTile+"-"+debugLevel,dbg_titles);
+									}
+									
+									
+									
+									st_planes = pd0.createTilePlanesFromSelections(
+											"" + nsTile+" bridges ",  // String        suffix,
+											split_sels,               // boolean [][][] plane_selections, //  = new boolean [nStiles][][][]; // num_tiles
+											disp_strength[nsTile],    // double  [][][] disp_strength,			
+											plDispNorm,               // double       dispNorm,   //  Normalize disparities to the average if above
+											plMinPoints,              // int          min_tiles,
+											plTargetEigen,            // double       plTargetEigen, //        =   0.1;  // Remove outliers until main axis eigenvalue (possibly scaled by plDispNorm) gets below
+											0.0,                      // double       plFractOutliers, //      =   0.3;  // Maximal fraction of outliers to remove
+											0,                        // plMaxOutliers,            // int          plMaxOutliers, //        =    20;  // Maximal number of outliers to remove
+											strength_floor,           // double       strength_floor,
+											strength_pow,             // double       strength_pow,
+											correct_distortions,      // boolean      correct_distortions,
+											smplMode,                 // boolean      smplMode, //        = true;   // Use sample mode (false - regular tile mode)
+											smplSide,                 // int          smplSide, //        = 2;      // Sample size (side of a square)
+											smplNum,                  //  int          smplNum, //         = 3;      // Number after removing worst
+											smplRms,                  //  double       smplRms, //         = 0.1;    // Maximal RMS of the remaining tiles in a sample
+											dl);                       // int          debugLevel);
+								}
+								
+							}
+							
+							if ((st_planes != null) && (!st_planes.isEmpty())){
+								
 								if (LOWEST_PLANE(2) > 0) st_planes.add(0, st_planes.get(0)); // insert dummy at pos 0;
 								result_planes[nsTile] = st_planes.toArray(new TilePlanes.PlaneData[0] );
 								if (LOWEST_PLANE(2) > 0) result_planes[nsTile][0] = null; // remove dummy
@@ -2465,7 +2256,7 @@ public class SuperTiles{
 				smplNum,             // final int        smplNum, //         = 3;      // Number after removing worst
 				smplRms,             // final double     smplRms, //         = 0.1;    // Maximal RMS of the remaining tiles in a sample
 
-				debugLevel+1,          // final int        debugLevel,
+				debugLevel + 2, // 1,          // final int        debugLevel,
 				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
@@ -2488,132 +2279,6 @@ public class SuperTiles{
 		return wh;
 	}
 
-	// TODO: Obsolete (never used)?
-	TilePlanes.PlaneData [][] getNeibPlanes(
-			final int     dir,     // 0: get from up (N), 1:from NE, ... 7 - from NW
-			final boolean preferDisparity,
-			final boolean dbgMerge, // Combine 'other' plane with current
-			final int     dbg_X,
-			final int     dbg_Y
-			)
-	{
-		final int tilesX =        tileProcessor.getTilesX();
-		final int tilesY =        tileProcessor.getTilesY();
-		final int superTileSize = tileProcessor.getSuperTileSize();
-		final int tileSize =      tileProcessor.getTileSize();
-		final int stilesX = (tilesX + superTileSize -1)/superTileSize;  
-		final int stilesY = (tilesY + superTileSize -1)/superTileSize;
-		final double [] nan_plane = new double [superTileSize*superTileSize];
-		for (int i = 0; i < nan_plane.length; i++) nan_plane[i] = Double.NaN;
-		final int [][] dirsYX = {{-1, 0},{-1,1},{0,1},{1,1},{1,0},{1,-1},{0,-1},{-1,-1}};
-		TilePlanes.PlaneData [][] neib_planes = new TilePlanes.PlaneData[stilesY * stilesX][];
-		TilePlanes tpl = null; // new TilePlanes(tileSize,superTileSize, geometryCorrection);
-		//				final int debug_stile = 20 * stilesX + 27;
-		//				final int debug_stile = 17 * stilesX + 27;
-		final int debug_stile = dbg_Y * stilesX + dbg_X;
-		//				final int debug_stile = -1;
-
-		for (int sty0 = 0; sty0 < stilesY; sty0++){
-			int sty = sty0 + dirsYX[dir][0];
-			for (int stx0 = 0; stx0 < stilesX; stx0++){
-				int nsTile0 = sty0 * stilesX + stx0; // result tile
-				int stx = stx0 + dirsYX[dir][1];
-				int nsTile = sty * stilesX + stx; // from where to get
-				if ((sty >= 0) && (sty < stilesY) && (stx >= 0) && (stx < stilesX) && (planes[nsTile] != null)) {
-					neib_planes[nsTile0] = new TilePlanes.PlaneData[planes[nsTile].length];
-					for (int np = 0; np < planes[nsTile].length; np++){
-						GeometryCorrection geometryCorrection = planes[nsTile][np].getGeometryCorrection();
-						boolean correct_distortions = planes[nsTile][np].getCorrectDistortions();
-						int [] sTileXY0 = {stx0,sty0}; 
-						if (tpl == null) {
-							tpl = new TilePlanes(tileSize,superTileSize, geometryCorrection);
-						}
-						TilePlanes.PlaneData pd = tpl.new PlaneData(
-								sTileXY0,
-								tileSize,
-								superTileSize,
-								geometryCorrection,
-								correct_distortions);
-						/*
-								TilePlanes.PlaneData pd = planes[nsTile0][np].clone();
-								pd.setS
-						 */
-						if (nsTile0 == debug_stile) {
-							System.out.println("getNeibPlanes(): nsTile0="+nsTile0+", nsTile="+nsTile+", np = "+np+" dir = "+ dir ); 
-						}
-						TilePlanes.PlaneData other_pd = pd.getPlaneToThis(
-								planes[nsTile][np], // PlaneData otherPd,
-								(nsTile0 == debug_stile)? 1:0); // int       debugLevel)
-						int num_this = (planes[nsTile0] == null)? 0:(planes[nsTile0].length); 
-						int best_index = 0;
-						if (dbgMerge && (num_this > 0)) {
-							double best_value = Double.NaN;
-							if (num_this > 1){
-								// find the best candidate for merge
-								for (int indx = 1; indx < planes[nsTile0].length; indx++){
-									if (nsTile0 == debug_stile) {
-										System.out.println("this with this, same weight:");
-										TilePlanes.PlaneData this_this_pd = planes[nsTile0][indx].mergePlaneToThis(
-												planes[nsTile0][indx], // PlaneData otherPd,
-												1.0,      // double    scale_other,
-												1.0, // double     starWeightPwr,    // Use this power of tile weight when calculating connection cost
-												true,     // boolean   ignore_weights,
-												true, // boolean   sum_weights,
-												preferDisparity,
-												(nsTile0 == debug_stile)? 1:0); // int       debugLevel)
-										System.out.println("other with other, same weight:");
-										TilePlanes.PlaneData other_other_pd = other_pd.mergePlaneToThis(
-												other_pd, // PlaneData otherPd,
-												1.0,      // double    scale_other,
-												1.0, // double     starWeightPwr,    // Use this power of tile weight when calculating connection cost
-												true,     // boolean   ignore_weights,
-												true, // boolean   sum_weights,
-												preferDisparity,
-												(nsTile0 == debug_stile)? 1:0); // int       debugLevel)
-										System.out.println("other with this, same weight:");
-									}
-									TilePlanes.PlaneData merged_pd = planes[nsTile0][indx].mergePlaneToThis(
-											other_pd, // PlaneData otherPd,
-											1.0,      // double    scale_other,
-											1.0, // double     starWeightPwr,    // Use this power of tile weight when calculating connection cost
-											true,     // boolean   ignore_weights,
-											true, // boolean   sum_weights,
-											preferDisparity,
-											(nsTile0 == debug_stile)? 1:0); // int       debugLevel)
-									if (!(merged_pd.getValue() > best_value)) { // Double.isNaN(best_value) will work too
-										best_value = merged_pd.getValue();
-										best_index = indx;
-									}
-								}
-							}
-							// now merge with weights with the best plane of this supertile
-							if (nsTile0 == debug_stile) {
-								System.out.println("this with this, proportional weight, np(other):"+np+" best fit(this):"+best_index);
-							}									
-							TilePlanes.PlaneData merged_pd = planes[nsTile0][best_index].mergePlaneToThis(
-									other_pd, // PlaneData otherPd,
-									1.0,      // double    scale_other,
-									1.0, // double     starWeightPwr,    // Use this power of tile weight when calculating connection cost
-									false,    // boolean   ignore_weights,
-									true, // boolean   sum_weights,
-									preferDisparity,
-									(nsTile0 == debug_stile)? 1:0); // int       debugLevel)
-							if (merged_pd != null) {
-								merged_pd.scaleWeight(0.5);
-							}
-							neib_planes[nsTile0][np] = merged_pd;
-						} else  {
-							neib_planes[nsTile0][np] = other_pd;
-						}
-					}
-				} else {
-					neib_planes[nsTile0] = null;
-				}
-			}
-		}
-		return neib_planes;
-	}
-	
 	public int fillSquares()
 	{
 		final int tilesX =        tileProcessor.getTilesX();
@@ -2807,7 +2472,7 @@ public class SuperTiles{
 		final int stilesX =       (tilesX + superTileSize -1)/superTileSize;
 		final int stilesY =       (tilesY + superTileSize -1)/superTileSize;
 		final int dbgTile = dbg_Y * stilesX + dbg_X;
-		final TileSurface.TileNeibs tnSurface = tileSurface.new TileNeibs(stilesX, stilesY);
+		final TileNeibs tnSurface = new TileNeibs(stilesX, stilesY);
 		int [] rslt = {0,0};
 		for (int nsTile = 0; nsTile < conflicts.length; nsTile++) if (conflicts[nsTile] != null) {
 			// conflicts may disappear after being fixed, recheck for null
@@ -2827,7 +2492,7 @@ public class SuperTiles{
 								start_dir4,      // int start_dir4,
 								right,           // boolean right,
 								reverse,         // boolean reverse,
-								tnSurface,       // TileSurface.TileNeibs tnSurface,
+								tnSurface,       // TileNeibs tnSurface,
 								conflicts,       // int [][][] conflicts,
 								conflict_stats,  // Conflicts conflict_stats, // to be updated after applying resolution
 								starSteps, // How far to look around when calculationg connection cost
@@ -2857,7 +2522,7 @@ public class SuperTiles{
 			int start_dir4,
 			boolean right,
 			boolean reverse,
-			TileSurface.TileNeibs tnSurface,
+			TileNeibs tnSurface,
 			int [][][] conflicts,
 			Conflicts conflict_stats, // to be updated after applying resolution
 			int        starSteps, // How far to look around when calculationg connection cost
@@ -2979,7 +2644,7 @@ public class SuperTiles{
 					nsTiles[isTile],   // int nsTile0,
 					replacement_tiles, //HashMap<Integer,Integer> replacement_tiles, //
 					neibs_prev, // int [][][] replacement_neibs,
-					tnSurface); // TileSurface.TileNeibs tnSurface)
+					tnSurface); // TileNeibs tnSurface)
 		}
 		if (debugLevel > 1) {
 			System.out.println("Calculating original conflicts");
@@ -3040,7 +2705,7 @@ public class SuperTiles{
 					nsTiles[isTile],   // int nsTile0,
 					replacement_tiles, //HashMap<Integer,Integer> replacement_tiles, //
 					neibs, // neibs, // int [][][] replacement_neibs,
-					tnSurface); // TileSurface.TileNeibs tnSurface)
+					tnSurface); // TileNeibs tnSurface)
 		}
 		Conflicts new_conflicts_stats = new Conflicts(
 				new_conflicts,
@@ -3118,7 +2783,7 @@ public class SuperTiles{
 			int nl1,
 			int nl2,
 			int dir_mask,
-			TileSurface.TileNeibs tnSurface,
+			TileNeibs tnSurface,
 			int [][][] conflicts,
 			Conflicts conflict_stats, // to be updated after applying resolution
 			int        starSteps, // How far to look around when calculationg connection cost
@@ -3220,7 +2885,7 @@ public class SuperTiles{
 					nsTiles[isTile],   // int nsTile0,
 					replacement_tiles, //HashMap<Integer,Integer> replacement_tiles, //
 					neibs_prev, // int [][][] replacement_neibs,
-					tnSurface); // TileSurface.TileNeibs tnSurface)
+					tnSurface); // TileNeibs tnSurface)
 		}
 		if (debugLevel > 1) {
 			System.out.println("Involved supertiles:");
@@ -3292,7 +2957,7 @@ public class SuperTiles{
 						nsTiles[isTile],   // int nsTile0,
 						replacement_tiles, //HashMap<Integer,Integer> replacement_tiles, //
 						neibs_vars[variant], // neibs, // int [][][] replacement_neibs,
-						tnSurface); // TileSurface.TileNeibs tnSurface)
+						tnSurface); // TileNeibs tnSurface)
 			}
 			variant_conflicts_stats[variant] = new Conflicts(
 					variant_conflicts[variant],
@@ -3414,7 +3079,7 @@ public class SuperTiles{
 		final int stilesX =       (tilesX + superTileSize -1)/superTileSize;
 		final int stilesY =       (tilesY + superTileSize -1)/superTileSize;
 		final int dbgTile = dbg_Y * stilesX + dbg_X;
-		final TileSurface.TileNeibs tnSurface = tileSurface.new TileNeibs(stilesX, stilesY);
+		final TileNeibs tnSurface = new TileNeibs(stilesX, stilesY);
 		int [] rslt = {0,0};
 		for (int nsTile = 0; nsTile < conflicts.length; nsTile++) if (conflicts[nsTile] != null) {
 			// conflicts may disappear after being fixed, recheck for null
@@ -3473,7 +3138,7 @@ public class SuperTiles{
 		final int stilesX =       (tilesX + superTileSize -1)/superTileSize;
 		final int stilesY =       (tilesY + superTileSize -1)/superTileSize;
 		final int dbgTile = dbg_Y * stilesX + dbg_X;
-		final TileSurface.TileNeibs tnSurface = tileSurface.new TileNeibs(stilesX, stilesY);
+		final TileNeibs tnSurface = new TileNeibs(stilesX, stilesY);
 		int [] rslt = {0,0};
 		for (int nsTile = 0; nsTile < conflicts.length; nsTile++) if (conflicts[nsTile] != null) {
 			// conflicts may disappear after being fixed, recheck for null
@@ -3511,7 +3176,7 @@ public class SuperTiles{
 			int nsTile,
 			int nl1,
 			int nl2,
-			TileSurface.TileNeibs tnSurface,
+			TileNeibs tnSurface,
 			int [][][] conflicts,
 			Conflicts conflict_stats, // to be updated after applying resolution
 			int        starSteps, // How far to look around when calculationg connection cost
@@ -3619,7 +3284,7 @@ public class SuperTiles{
 					nsTiles[isTile],   // int nsTile0,
 					replacement_tiles, //HashMap<Integer,Integer> replacement_tiles, //
 					neibs_prev, // int [][][] replacement_neibs,
-					tnSurface); // TileSurface.TileNeibs tnSurface)
+					tnSurface); // TileNeibs tnSurface)
 		}
 		if (debugLevel > 1) {
 			System.out.println("Involved supertiles:");
@@ -3649,7 +3314,7 @@ public class SuperTiles{
 						nsTiles[isTile],   // int nsTile0,
 						replacement_tiles, //HashMap<Integer,Integer> replacement_tiles, //
 						neibs_vars[variant], // neibs, // int [][][] replacement_neibs,
-						tnSurface); // TileSurface.TileNeibs tnSurface)
+						tnSurface); // TileNeibs tnSurface)
 			}
 			variant_conflicts_stats[variant] = new Conflicts(
 					variant_conflicts[variant],
@@ -3809,7 +3474,7 @@ public class SuperTiles{
 			int nl1,
 			int nl2,
 			LinkPlanes lp,
-			TileSurface.TileNeibs tnSurface,
+			TileNeibs tnSurface,
 			int [][][] conflicts,
 			Conflicts conflict_stats, // to be updated after applying resolution
 			int        starSteps, // How far to look around when calculationg connection cost
@@ -3922,7 +3587,7 @@ public class SuperTiles{
 					nsTiles[isTile],   // int nsTile0,
 					replacement_tiles, // HashMap<Integer,Integer> replacement_tiles, //
 					neibs_prev, // int [][][] replacement_neibs,
-					tnSurface); // TileSurface.TileNeibs tnSurface)
+					tnSurface); // TileNeibs tnSurface)
 // calculate cost of conflicts			
 			
 			conflicts_old_cost += iconflicts.getConflictsCost(
@@ -3932,7 +3597,7 @@ public class SuperTiles{
 					replacement_tiles,     // HashMap<Integer,Integer> replacement_tiles, // null is OK
 					neibs_prev,            // int [][][]    replacement_neibs,               // null OK if  replacement_tiles == null
 					connectionCosts,       // ConnectionCosts connectionCosts,
-					tnSurface); // TileSurface.TileNeibs tnSurface)
+					tnSurface); // TileNeibs tnSurface)
 
 			
 		}
@@ -3971,7 +3636,7 @@ public class SuperTiles{
 						nsTiles[isTile],   // int nsTile0,
 						replacement_tiles, //HashMap<Integer,Integer> replacement_tiles, //
 						neibs_vars[variant], // neibs, // int [][][] replacement_neibs,
-						tnSurface); // TileSurface.TileNeibs tnSurface)
+						tnSurface); // TileNeibs tnSurface)
 				// connectionCosts now contains last calculated	val/weight pairs for broader array of tile data
 				conflicts_var_cost[variant] += iconflicts.getConflictsCost(
 						nsTiles[isTile],                    // int           nsTile0,
@@ -3980,7 +3645,7 @@ public class SuperTiles{
 						replacement_tiles,                  // HashMap<Integer,Integer> replacement_tiles, // null is OK
 						neibs_vars[variant],                // int [][][]    replacement_neibs,               // null OK if  replacement_tiles == null
 						connectionCosts,                    // ConnectionCosts connectionCosts,
-						tnSurface);                         // TileSurface.TileNeibs tnSurface)
+						tnSurface);                         // TileNeibs tnSurface)
 			}
 			variant_conflicts_stats[variant] = new Conflicts(
 					variant_conflicts[variant],
@@ -4159,7 +3824,7 @@ public class SuperTiles{
 		final int stilesX =       (tilesX + superTileSize -1)/superTileSize;  
 		final int stilesY =       (tilesY + superTileSize -1)/superTileSize;
 		final int nStiles =       stilesX * stilesY; 
-		final TileSurface.TileNeibs tnSurface = tileSurface.new TileNeibs(stilesX, stilesY);
+		final TileNeibs tnSurface = 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++) {
@@ -4174,7 +3839,7 @@ public class SuperTiles{
 							starValPwr,     //double          starValPwr, //  Raise value of each tile before averaging
 							steps,          // int            steps,
 							planes,         // TilePlanes.PlaneData [][] planes,
-							tnSurface,      // TileSurface.TileNeibs tnSurface,
+							tnSurface,      // TileNeibs tnSurface,
 							preferDisparity); // boolean preferDisparity)
 					int [] mod_supertiles = new int[1];
 					for (int nsTile = ai.getAndIncrement(); nsTile < nStiles; nsTile = ai.getAndIncrement()) {
@@ -4214,7 +3879,7 @@ public class SuperTiles{
 		//				final int tileSize =      tileProcessor.getTileSize();
 		final int stilesX =       (tilesX + superTileSize -1)/superTileSize;  
 		final int stilesY =       (tilesY + superTileSize -1)/superTileSize;
-		final TileSurface.TileNeibs tnSurface = tileSurface.new TileNeibs(stilesX, stilesY);
+		final TileNeibs tnSurface = 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++) {
@@ -4229,7 +3894,7 @@ public class SuperTiles{
 							starValPwr,      //double     starValPwr, //  Raise value of each tile before averaging
 							steps,          // int            steps,
 							planes,         // TilePlanes.PlaneData [][] planes,
-							tnSurface,      // TileSurface.TileNeibs tnSurface,
+							tnSurface,      // TileNeibs tnSurface,
 							preferDisparity); // boolean preferDisparity)
 					int [] supertiles = new int[1];
 					for (int isTile = ai.getAndIncrement(); isTile < mod_supertiles.length; isTile = ai.getAndIncrement()) {
@@ -4273,7 +3938,7 @@ public class SuperTiles{
 			int nl1,
 			int nl2,
 			int dir_mask,
-			TileSurface.TileNeibs tnSurface,
+			TileNeibs tnSurface,
 			int        debugLevel)
 	{
 		int [][] layers0 = new int[5][2]; // 4 - center, 0 - N, 1 - E, 2 - S, 3 - w. 4 has 2, each other 0,1 or 2 layers
@@ -4684,7 +4349,7 @@ public class SuperTiles{
 		final int stilesX =       (tilesX + superTileSize -1)/superTileSize;
 		final int stilesY =       (tilesY + superTileSize -1)/superTileSize;
 		final int dbgTile = dbg_Y * stilesX + dbg_X;
-		final TileSurface.TileNeibs tnSurface = tileSurface.new TileNeibs(stilesX, stilesY);
+		final TileNeibs tnSurface = new TileNeibs(stilesX, stilesY);
 		final Conflicts iconflicts = new Conflicts(this); 
 
 		int [] rslt = {0,0};
@@ -4870,12 +4535,12 @@ public class SuperTiles{
 												nsTiles[isTile],   // int nsTile0,
 												replacement_tiles, //HashMap<Integer,Integer> replacement_tiles, //
 												neibs_prev, // int [][][] replacement_neibs,
-												tnSurface); // TileSurface.TileNeibs tnSurface)
+												tnSurface); // TileNeibs tnSurface)
 										conflicts_new[isTile] = iconflicts.detectTriangularTileConflicts(
 												nsTiles[isTile],   // int nsTile0,
 												replacement_tiles, //HashMap<Integer,Integer> replacement_tiles, //
 												neibs, // int [][][] replacement_neibs,
-												tnSurface); // TileSurface.TileNeibs tnSurface)
+												tnSurface); // TileNeibs tnSurface)
 									}
 									if (debugLevel > 1) {
 										System.out.println("Involved supertiles:");
@@ -4974,7 +4639,7 @@ public class SuperTiles{
 			int dir4,
 			double maxEigen, // maximal eigenvalue of planes to consider
 			double     starWeightPwr,    // Use this power of tile weight when calculating connection cost			
-			TileSurface.TileNeibs tnSurface,
+			TileNeibs tnSurface,
 			boolean preferDisparity,
 			int debugLevel)
 	{
@@ -5067,7 +4732,7 @@ public class SuperTiles{
 			double diagonalWeight,
 			double starPwr,    // Divide cost by number of connections to this power
 			double     starWeightPwr,    // Use this power of tile weight when calculating connection cost			
-			TileSurface.TileNeibs tnSurface,
+			TileNeibs tnSurface,
 			boolean preferDisparity,
 			int    debugLevel)
 	{
@@ -5106,7 +4771,7 @@ public class SuperTiles{
 	
 	public int [] getInvolvedSupertiles(
 			int [] mod_supertiles,
-			TileSurface.TileNeibs tnSurface)
+			TileNeibs tnSurface)
 	{
 		HashSet<Integer> stiles_set = new HashSet<Integer>();
 		for (int i = 0; i < mod_supertiles.length; i++){
@@ -5157,7 +4822,7 @@ public class SuperTiles{
 			double         diagonalWeight,
 			double         starPwr,    // Divide cost by number of connections to this power
 			double     starWeightPwr,    // Use this power of tile weight when calculating connection cost
-			TileSurface.TileNeibs tnSurface,
+			TileNeibs tnSurface,
 			boolean preferDisparity,
 			int    debugLevel)
 	{
@@ -6156,7 +5821,7 @@ public class SuperTiles{
 										smplSide,
 										smplNum,
 										smplRms,
-										dl);            // int          debugLevel)
+										dl - 2);            // int          debugLevel)
 								if (disp_strength == null) {
 									System.out.println("=== BUG in applyMergePlanes(): failed to getPlaneFromMeas() for merged planes");  
 									break;
@@ -6177,7 +5842,7 @@ public class SuperTiles{
 										smplSide,
 										smplNum,
 										smplRms,
-										dl);            // int          debugLevel)
+										dl - 2);            // int          debugLevel)
 								
 								
 								

src/main/java/TileNeibs.java

+import java.util.ArrayList;
+import java.util.Arrays;
+
+/**
+ ** TileNeibs - handles walking inside rectangular area
+ **
+ ** Copyright (C) 2017 Elphel, Inc.
+ **
+ ** -----------------------------------------------------------------------------**
+ **  
+ **  TileNeibs.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 TileNeibs{
+	int sizeX;
+	int sizeY;
+	public int dirs = 8;
+	public TileNeibs(int size){
+		this.sizeX = size;
+		this.sizeY = size;
+	}
+	public TileNeibs(int sizeX, int sizeY){
+		this.sizeX = sizeX;
+		this.sizeY = sizeY;
+	}
+	public int numNeibs() // TODO: make configurable to 
+	{
+		return dirs;
+	}
+	public int opposite(int dir){
+		return (dir + dirs / 2) % dirs;
+	}
+
+	/**
+	 * Get x,y pair from index
+	 * @param indx element index
+	 * @return array of {x,y}
+	 */
+	
+	int [] getXY(int indx)
+	{
+		int [] xy = {indx % sizeX ,indx / sizeX};
+		return xy;
+	}
+	
+	/**
+	 * Get element index from x and y
+	 * @param x horizontal position
+	 * @param y vertical position
+	 * @return element linescan index
+	 */
+	
+	int getIndex(int x, int y){
+		return y * sizeX + x;
+	}
+	
+	/**
+	 * Get 2d element index after step N, NE, ... NW. Returns -1 if leaving array   
+	 * @param indx start index
+	 * @param dir step direction (CW from up)
+	 * @return new index or -1 if leaving 
+	 */
+	int getNeibIndex(int indx, int dir)
+	{
+		int y = indx / sizeX;
+		int x = indx % sizeX;
+		if (dir < 0) return indx;
+		switch (dir % dirs){
+		case 0: return (y == 0) ?                                    -1 : (indx - sizeX); 
+		case 1: return ((y == 0)           || ( x == (sizeX - 1))) ? -1 : (indx - sizeX + 1); 
+		case 2: return (                      ( x == (sizeX - 1))) ? -1 : (indx        + 1); 
+		case 3: return ((y == (sizeY - 1)) || ( x == (sizeX - 1))) ? -1 : (indx + sizeX + 1); 
+		case 4: return ((y == (sizeY - 1))                       ) ? -1 : (indx + sizeX); 
+		case 5: return ((y == (sizeY - 1)) || ( x == 0))           ? -1 : (indx + sizeX - 1); 
+		case 6: return (                      ( x == 0))           ? -1 : (indx        - 1); 
+		case 7: return ((y == 0)           || ( x == 0))           ? -1 : (indx - sizeX - 1); 
+		default: return indx;
+		}
+	}
+	
+	/**
+	 * Return tile segment for 50% overlap. -1 - center, 0 N, 1 - NE,... 7 - NW
+	 * @param indx element index
+	 * @return which of the 9 areas this element belongs
+	 */
+	int getSegment(int indx)
+	{
+		int s1x = sizeX / 4;
+		int s3x = 3 * sizeX / 4;
+		int s1y = sizeY / 4;
+		int s3y = 3 * sizeY / 4;
+		int x = indx % sizeX;
+		int y = indx / sizeX;
+		boolean up = y < s1y;
+		boolean down = y >= s3y;
+		boolean left = x < s1x;
+		boolean right = x >= s3x;
+		if (up){
+			if (left) return 7;
+			if (right) return 1;
+			return 0;
+		}
+		if (down){
+			if (left) return 5;
+			if (right) return 3;
+			return 4;
+		}
+		if (left) return 6;
+		if (right) return 2;
+		return -1;
+	}
+
+	/**
+	 * Find if the step leaves the center half of all area
+	 * @param indx start point
+	 * @param dir direction
+	 * @return direction to the new tile (assuming 50% overlap) or -1 if did not cross the border
+	 */
+	int leaveOvderlapedCenter(int indx, int dir)
+	{
+		int segm = getSegment(indx);
+		int indx1 = getNeibIndex(indx, dir);
+		if (indx1 < 0 ) {
+			return -1; // should not happen
+		}
+		int segm1 = getSegment(indx1);
+		if (segm1 == segm) return -1;
+		if (segm == -1) return segm1;
+		int [][] dxy = {
+				{ 0, 0},
+				{-1, 0},
+				{-1, 1},
+				{ 0, 1},
+				{ 1, 1},
+				{ 1, 0},
+				{ 1,-1},
+				{ 0,-1},
+				{-1,-1}};
+		int dx = dxy[segm1 + 1][1] - dxy[segm + 1][1]; 
+		int dy = dxy[segm1 + 1][0] - dxy[segm + 1][0];
+		for (int dp1 = 0; dp1 <=8; dp1++) {
+			int sdx = (dx > 0) ? 1: ( (dx < 0) ? -1 : 0);
+			int sdy = (dy > 0) ? 1: ( (dy < 0) ? -1 : 0);
+			if ((dxy[dp1][0] == sdy) && (dxy[dp1][1] == sdx)){
+				return dp1 -1;
+			}
+		}
+		return -1; // should not happen
+	}
+	
+	public void growSelection(
+			int        grow,           // grow tile selection by 1 over non-background tiles 1: 4 directions, 2 - 8 directions, 3 - 8 by 1, 4 by 1 more
+			boolean [] tiles,
+			boolean [] prohibit)
+	{
+		boolean [] src_tiles = tiles.clone(); // just in case
+		// grow
+		boolean hor = true;
+		for (; grow > 0; grow--){
+			boolean single = (grow ==1) && hor;
+			src_tiles = tiles.clone();
+			int num_new = 0;
+			if (hor){
+				for (int tileY = 0; tileY < sizeY; tileY++){
+					for (int tileX = 0; tileX < (sizeX - 1); tileX++){
+						int tindx = tileY * sizeX + tileX;
+						if ((prohibit == null) || (!prohibit[tindx] && !prohibit[tindx + 1])) {
+							if (!tiles[tindx + 1] && src_tiles[tindx]){
+								num_new++;
+							}
+							tiles[tindx + 1] |= src_tiles[tindx];  
+						}
+
+					}
+					for (int tileX = 1; tileX < sizeX; tileX++){
+						int tindx = tileY * sizeX + tileX;
+						if ((prohibit == null) || (!prohibit[tindx] && !prohibit[tindx - 1])) {
+							if (!tiles[tindx - 1] && src_tiles[tindx]){
+								num_new++;
+							}
+							tiles[tindx - 1] |= src_tiles[tindx];
+						}
+					}
+				}
+			}
+			if (!hor || single){ // do vertically, but from previous state
+				for (int tileX = 0; tileX < sizeX; tileX++){
+					for (int tileY = 0; tileY < (sizeY - 1); tileY++){
+						int tindx = tileY * sizeX + tileX;
+						if ((prohibit == null) || (!prohibit[tindx] && !prohibit[tindx + sizeX])) {
+							if (!tiles[tindx + sizeX] && src_tiles[tindx]){
+								num_new++;
+							}
+							tiles[tindx + sizeX] |= src_tiles[tindx];
+						}
+
+					}
+					for (int tileY = 1; tileY < sizeY; tileY++){
+						int tindx = tileY * sizeX + tileX;
+						if ((prohibit == null) || (!prohibit[tindx] && !prohibit[tindx - sizeX])) {
+							if (!tiles[tindx - sizeX] && src_tiles[tindx]){
+								num_new++;
+							}
+							tiles[tindx - sizeX] |= src_tiles[tindx];
+						}
+					}
+				}
+			}
+			hor = !hor;
+			if (num_new == 0){
+				break;
+			}
+		}
+	}
+	
+	public boolean [] boundShape(
+			boolean [] selection,
+			boolean octo)
+	{
+		boolean [] bound_shape = new boolean [selection.length];
+		int min_x=-1, max_x=-1, min_y=-1, max_y=-1; 
+		int min_s=-1, max_s=-1, min_d=-1, max_d=-1; 
+		boolean is_set = false;
+		for (int i = 0; i < selection.length; i++) if (selection[i]){
+			int [] xy = getXY(i);
+			int [] sd = {xy[0]+xy[1],xy[0]-xy[1]};
+			if (!is_set) {
+				min_x = xy[0];
+				max_x = xy[0];
+				min_y = xy[1];
+				max_y = xy[1];
+				if (octo) {
+				min_s = sd[0];
+				max_s = sd[0];
+				min_d = sd[1];
+				max_d = sd[1];
+				}
+				is_set = true;
+			} else {
+				if      (xy[0] < min_x) min_x = xy[0];
+				else if (xy[0] > max_x) max_x = xy[0]; 
+				if      (xy[1] < min_y) min_y = xy[1];
+				else if (xy[1] > max_y) max_y = xy[1];
+				if (octo) {
+					if      (sd[0] < min_s) min_s = sd[0];
+					else if (sd[0] > max_s) max_s = sd[0]; 
+					if      (sd[1] < min_d) min_d = sd[1];
+					else if (sd[1] > max_d) max_d = sd[1];
+				}
+				
+			}
+		}
+		for (int y = min_y; y <= max_y; y++){
+			for (int x = min_x; x <= max_x; x++){
+				if (!octo ||
+				(((x + y) >= min_s) && ((x + y) <= max_s) && ((x - y) >= min_d) && ((x - y) <= max_d))) {
+					bound_shape[getIndex(x,y)] = true; 
+					
+				}
+			}
+		}
+		return bound_shape;
+	}
+	
+	
+	/**
+	 * Enumerate clusters on rectangular area
+	 * @param tiles   selected tiles, size should be sizeX * sizeY
+	 * @param ordered if true, order tiles from largest to smallest5
+	 * @return integer array, where 0 is unused, 1+ cluster it belongs to
+	 */
+	
+	
+	public int [] enumerateClusters(
+			boolean [] tiles,
+			boolean ordered)
+	{
+		int [] waves = new int [tiles.length];
+		for (int i = 0; i < tiles.length; i++) waves[i] = tiles[i] ? 0: -1;
+		ArrayList<Integer> front = new ArrayList<Integer>();
+		int [] enum_clust = new int[tiles.length];
+		int numClust = 0;
+		for (int start_indx = 0; start_indx < tiles.length; start_indx++) if (waves[start_indx] == 0){ // found first pixel of a new cluster
+			numClust ++;
+			Integer ipx = start_indx;
+			Integer ipx1;
+			front.clear();
+			int area = 1;
+			waves[ipx] = area;
+			enum_clust[ipx] = numClust;
+			front.add(ipx);
+			while (!front.isEmpty()) {
+				ipx = front.remove(0);// get oldest element
+				for (int d = 0; d < dirs; d++){
+					ipx1 = 	getNeibIndex(ipx, d);
+					if (ipx1 >= 0){
+						if (waves[ipx1] == 0) {
+							area++;
+							waves[ipx1] = area;
+							enum_clust[ipx1] = numClust;
+							front.add(ipx1);
+						}
+					}
+				}
+			}
+		}
+		if (!ordered) {
+			return enum_clust;
+		}
+		
+		// count cluster
+		int []clustSizes = new int [numClust];
+		for (int i = 0; i < clustSizes.length; i++) clustSizes[i] = 0;
+		for (int i = 0; i < enum_clust.length; i++) if (enum_clust[i] > 0) clustSizes[enum_clust[i]-1]++;
+
+		class Pair implements Comparable<Pair> {
+			public final int index;
+			public final int value;
+
+			public Pair(int index, int value) {
+				this.index = index;
+				this.value = value;
+			}
+
+			@Override
+			public int compareTo(Pair other) {
+				return Integer.valueOf(this.value).compareTo(other.value);
+			}
+		}
+
+		Pair[] pairs = new Pair[numClust];
+		for (int i = 0; i < clustSizes.length; i++) pairs [i] = new Pair (i, clustSizes[i]);
+		Arrays.sort(pairs);
+		int [] revIndex = new int [numClust];
+		for (int i = 0; i < revIndex.length; i++) revIndex [pairs[i].index] = (numClust - i); // array was in accending order
+		int [] enum_clust_ordered = new int[tiles.length];
+		for (int i=0; i < enum_clust_ordered.length; i++){
+			enum_clust_ordered[i] = (enum_clust[i] > 0) ? revIndex[enum_clust[i] - 1] : 0; 
+		}
+		return enum_clust_ordered;  
+	}
+	public int getMax(
+			int [] data)
+	{
+		int mx = data[0];
+		for (int i = 1; i < data.length; i++){
+			if (data[i] > mx){
+				mx= data[i];
+			}
+		}
+		return mx;
+	}
+}

src/main/java/TilePlanes.java

@@ -22,8 +22,10 @@
  **
  */
 import java.util.ArrayList;
+import java.util.Arrays;
 import java.util.Collections;
 import java.util.Comparator;
+import java.util.HashSet;
 
 import Jama.EigenvalueDecomposition;
 import Jama.Matrix;
@@ -2210,13 +2212,13 @@ public class TilePlanes {
 			return plane;
 		}
 
-		public double[][] getDoublePlaneDisparityStrength(
+		public double[][] getDoublePlaneDisparityStrength_old(
 				double [] window,
 				boolean   use_sel,
 				boolean   divide_by_area,
 				double    scale_projection,
 				int       debugLevel)
-		{		return getDoublePlaneDisparityStrength(
+		{		return getDoublePlaneDisparityStrength_old(
 				true,
 				window,
 				use_sel,
@@ -2259,7 +2261,7 @@ public class TilePlanes {
 		 * @return a pair of arrays {disparity, strength}, each [2*superTileSize * 2*superTileSize]
 		 */
 		// obsolete, convert to another version ?
-		public double[][] getDoublePlaneDisparityStrength(
+		public double[][] getDoublePlaneDisparityStrength_old(
 				boolean   useWorld,
 				double [] window,
 				boolean   use_sel,
@@ -3331,14 +3333,11 @@ public class TilePlanes {
 			return world_xyz;
 		}
 		
-		public ArrayList<TilePlanes.PlaneData> createTilePlanesFromSelections(
+		public ArrayList<PlaneData> createTilePlanesFromSelections(
 				String        suffix,
 				boolean [][][] plane_selections, //  = new boolean [nStiles][][][]; // num_tiles
 				double  [][][] disp_strength,			
-				// double       disp_far, // minimal disparity to select (or NaN)
-				// double       disp_near, // maximal disparity to select (or NaN)
 				double       dispNorm,   //  Normalize disparities to the average if above
-//				double       min_weight,
 				int          min_tiles,
 				double       plTargetEigen, //        =   0.1;  // Remove outliers until main axis eigenvalue (possibly scaled by plDispNorm) gets below
 				double       plFractOutliers, //      =   0.3;  // Maximal fraction of outliers to remove
@@ -3354,11 +3353,11 @@ public class TilePlanes {
 		{
 			
 			// first make a plane from all tiles
-			ArrayList<TilePlanes.PlaneData> st_planes = new ArrayList<TilePlanes.PlaneData>();
+			ArrayList<PlaneData> st_planes = new ArrayList<PlaneData>();
 			
 			// iterate through all plane selections
 			for (int ps = 0; ps < plane_selections.length; ps++) {
-				TilePlanes.PlaneData pd = this.clone(); 
+				PlaneData pd = this.clone(); 
 				boolean OK = (pd.getPlaneFromMeas(
 						plane_selections[ps], // tile_sel,       // boolean [][] tile_sel, // null - do not use, {} use all (will be modified)
 						disp_strength,
@@ -3426,7 +3425,6 @@ public class TilePlanes {
 
 					}
 					// calculate the world planes too									
-//					if (debugLevel > -1){
 						pd.getWorldPlaneFromMeas(
 								plane_selections[ps], // tile_sel,       // boolean [][] tile_sel, // null - do not use, {} use all (will be modified)
 								disp_strength,
@@ -3447,9 +3445,9 @@ public class TilePlanes {
 			}
 			if (st_planes.size() > 0){
 				// sort planes by increasing disparity (tile center or plane center ? ) Using plane center
-				Collections.sort(st_planes, new Comparator<TilePlanes.PlaneData>() {
+				Collections.sort(st_planes, new Comparator<PlaneData>() {
 					@Override
-					public int compare(TilePlanes.PlaneData lhs, TilePlanes.PlaneData rhs) {
+					public int compare(PlaneData lhs, PlaneData rhs) {
 						// -1 - less than, 1 - greater than, 0 - equal
 						return (rhs.getZxy()[0] > lhs.getZxy()[0]) ? -1 : (rhs.getZxy()[0] < lhs.getZxy()[0] ) ? 1 : 0;
 					}
@@ -3459,296 +3457,351 @@ public class TilePlanes {
 			return null;		
 		}
 
-		
-		
-		
-	}
-	
-	
-	
-	
-	//TODO: Remove below methods and promote PlaneData (no TilePlanes) after tested
-	
-	/**
-	 * Calculate covariance matrix for a subset of tile data (disparities)
-	 * Subtract weight floor from weight
-	 * @param data    data array - square (2*stSize) * (2*stSize)  
-	 * @param weight  per sample weight (should have floor already subtracted)
-	 * @param select  sample selection
-	 * @return covariance (diagonal) matrix: [0]: disparity, [1]: d<disparity>/dx, [2]: d<disparity>/dy,
-	 */
-	public double [][][] getCovar(
-			double []  data,
-			double []  weight,
-			boolean [] select,
-			double     plDispNorm, //  Normalize disparities to the average if above
-			int        debugLevel){
-		double mindet = 1E-15;
-		int stSize2 = 2 * stSize;
-//		Matrix covar = new Matrix(3,3);
-		double [][] acovar = new double [3][3];
-		int numPoints = 0;
-		double sw =0.0, swz = 0.0, swx = 0.0, swy = 0.0;
-
-		for (int indx = 0; indx < data.length; indx++){
-			if (select[indx] && (weight[indx] > 0)){
-				numPoints++;
-				double w = weight[indx];
-				double d = data[indx];
-				// referencing samples to centers of pixels
-				double x = ((indx % stSize2) - stSize + 0.5) * tileSize + 0.5; // in pixels, not in tiles
-				double y = ((indx / stSize2) - stSize + 0.5) * tileSize + 0.5;
-				sw  += w;
-				swz += w * d;
-				swx += w * x;
-				swy += w * y;
+		public boolean [][][]  refineDiscriminateTiles(
+				final PlaneData  [] planes, 
+				final int        stMeasSel, //            = 1;      // Select measurements for supertiles : +1 - combo, +2 - quad +4 - hor +8 - vert
+				final double     plDispNorm,
+				final int        plMinPoints, //          =     5;  // Minimal number of points for plane detection
+				final GeometryCorrection geometryCorrection,
+				final boolean    correct_distortions,
+
+				final boolean    smplMode, //        = true;   // Use sample mode (false - regular tile mode)
+				final int        smplSide, //        = 2;      // Sample size (side of a square)
+				final int        smplNum, //         = 3;      // Number after removing worst
+				final double     smplRms, //         = 0.1;    // Maximal RMS of the remaining tiles in a sample
+				
+				final double     max_disp_diff,    // maximal disparity difference from the plane to consider tile 
+				final double     disp_range,       // parallel move known planes around original know value for the best overall fit
+				final int        amplitude_steps,  // number of steps (each direction) for each plane to search for the best fit (0 - single, 1 - 1 each side)
+				final int        num_variants,     // total number of variants to try (protect from too many planes) 
+				final int        mode, // 0 - weighted, 1 - equalized, 2 - best, 3 - combined
+				final int        debugLevel)
+		{
+			if (planes == null) return null;
+			// create a list of usable planes according to the mode
+			ArrayList<PlaneData> tilePlanes = new ArrayList<PlaneData>();
+			for (int np = 0; np < planes.length; np++) if (planes[np] != null){
+				PlaneData weighted_pd = planes[np].getNonexclusiveStar();
+				PlaneData equal_pd =    planes[np].getNonexclusiveStarEq();
+				switch (mode){
+				case 0:
+					if (weighted_pd == null)throw new IllegalArgumentException ("refineDiscriminateTiles(): getNonexclusiveStar() retgurned null");
+					tilePlanes.add(weighted_pd);
+					break;
+				case 1:
+					if (equal_pd == null)throw new IllegalArgumentException ("refineDiscriminateTiles(): getNonexclusiveStarEq() retgurned null");
+					tilePlanes.add(equal_pd);
+					break;
+				case 2:
+					if (weighted_pd == null)throw new IllegalArgumentException ("refineDiscriminateTiles(): getNonexclusiveStar() retgurned null");
+					if (equal_pd == null)throw new IllegalArgumentException ("refineDiscriminateTiles(): getNonexclusiveStarEq() retgurned null");
+					if (weighted_pd.getWeight() > equal_pd.getWeight()) tilePlanes.add(weighted_pd);
+					else                                                tilePlanes.add(equal_pd);
+					break;
+				case 3:
+					if (weighted_pd == null)throw new IllegalArgumentException ("refineDiscriminateTiles(): getNonexclusiveStar() retgurned null");
+					if (equal_pd == null)throw new IllegalArgumentException ("refineDiscriminateTiles(): getNonexclusiveStarEq() retgurned null");
+					PlaneData combo_pd =  weighted_pd.mergePlaneToThis(
+							equal_pd, // PlaneData otherPd,
+							1.0,      // double    scale_other,
+							1.0,      // double    starWeightPwr,    // Use this power of tile weight when calculating connection cost
+							false,    // boolean   ignore_weights,
+							true,     // boolean   sum_weights,
+							preferDisparity, // boolean   preferDisparity, // Always start with disparity-most axis (false - lowest eigenvalue)
+							debugLevel - 2); // int       debugLevel)
+					break;
+				default:
+		    		throw new IllegalArgumentException ("refineDiscriminateTiles(): invalid mode="+mode);
 				}
 			}
-		if (sw == 0.0) {
+			if (tilePlanes.isEmpty()){
 				return null;
 			}
-		swz /= sw;
-		swx /= sw;
-		swy /= sw;
-		
-		
-		// TODO: scale disparity to make same scale for 3 axes?
-		
-		double kz = ((plDispNorm > 0.0) && (swz > plDispNorm)) ? (plDispNorm / swz) : 1.0; 
-		if (debugLevel > 0){
-			System.out.println("getCovar(): sw = "+sw +", swz = "+swz +", swx = "+swx +", swy = "+swy+" kz = "+kz);
+           // get measured disparity/strength data, filtered, not tilted
+			double [][][] disp_strength = new double[measuredLayers.getNumLayers()][][];
+			for (int ml = 0; ml < disp_strength.length; ml++) if ((stMeasSel & ( 1 << ml)) != 0){
+				if (smplMode) {
+					disp_strength[ml] = measuredLayers.getDisparityStrength( // expensive to calculate (improve removing outlayers
+							ml, // int num_layer,
+							getSTileXY()[0],        // int stX,
+							getSTileXY()[1],        // int stY,
+							null,                   // boolean [] sel_in, - use all 
+							strength_floor,
+							measured_strength_pow,  //
+							smplSide,               // = 2;      // Sample size (side of a square)
+							smplNum,                // = 3;      // Number after removing worst
+							smplRms,                // = 0.1;    // Maximal RMS of the remaining tiles in a sample
+							true);                  // boolean null_if_none)
+				} else {
+					disp_strength[ml] = measuredLayers.getDisparityStrength(
+							ml,                     // int num_layer,
+							getSTileXY()[0],        // int stX,
+							getSTileXY()[1],        // int stY,
+							null,                   // boolean [] sel_in, - use all
+							strength_floor,         //  double strength_floor,
+							measured_strength_pow,  // double strength_pow,
+							true);                  // boolean null_if_none);
 				}
-		for (int indx = 0; indx < data.length; indx++){
-			if (select[indx] && (weight[indx] > 0)){
-				double w = weight[indx] / sw;
-				double d = kz * (data[indx] - swz);
-				double wd = w*d;
-				double x = ((indx % stSize2) - stSize + 0.5) * tileSize + 0.5 - swx;
-				double y = ((indx / stSize2) - stSize + 0.5) * tileSize + 0.5 - swy;
-				acovar [0][0] += wd * d;
-				acovar [0][1] += wd * x;
-				acovar [0][2] += wd * y;
-				acovar [1][1] += w * x * x;
-				acovar [1][2] += w * x * y;
-				acovar [2][2] += w * y * y;
 			}
+			double [] window =	getWindow (2 * superTileSize);
+ 
+			// get all original planes 
+			int num_planes = tilePlanes.size();
+			double [][][] pds = new double [num_planes][][];
+			for (int np = 0; np < pds.length; np++){
+				PlaneData pd = tilePlanes.get(np);
+				pds[np] = pd.getDoublePlaneDisparityStrength(
+						false, // boolean   useWorld,
+						null, // double [] window,
+						-1, // int       dir,
+						false, // boolean   use_sel,
+						false, // boolean   divide_by_area,
+						1.0, // double    scale_projection,
+						0.0, // double    fraction_uni,
+						0); // int       debugLevel)
+			}			
+			
+			// calculate number of the variants for each plane
+			int  extra_vars = (((int) Math.pow(num_variants, 1.0/num_planes)) - 1) / 2; // 0 - single, 1 - 3 (1 each direction), ...
+			if (extra_vars > amplitude_steps) extra_vars = amplitude_steps; //
+			int steps = 2 * extra_vars + 1;
+			double [] disps = new double [steps];
+			disps[0] = 0.0; // center
+			for (int i = 0; i < extra_vars; i++){
+				disps[2 * i + 1] =  (disp_range * (i + 1)) / extra_vars; 
+				disps[2 * i + 2] = -disps[2 * i + 1]; 
+			}
+			int [] state = new int [num_planes]; // variants counter
+//			for (int variant = 0; i < num_variants; variant++ )
+			while (true){
+				// evaluate current variant
+				double [] weights = new double [num_planes];
+				double [] err2 =    new double [num_planes];
+				boolean 
+				
+				
+
+				//  calculate next variant
+				boolean all_done = true;
+				for (int i = 0; i < num_planes; i ++){
+					if (state[i] < (steps - 1)){
+						state[i] ++;
+						for (int j = 0; j < i; j++){
+							state[j] = 0;
+						}
+						all_done = false;
+						break;
 					}
-		acovar [1][0] = acovar [0][1]; 
-		acovar [2][0] = acovar [0][2]; 
-		acovar [2][1] = acovar [1][2];
-		Matrix covar = new Matrix(acovar);
-//		if (Math.abs(covar.det()) < mindet){
-//			debugLevel = 5;
-//		}
-
-		EigenvalueDecomposition eig = covar.eig();
-		if (Double.isNaN(eig.getV().get(0, 0))){
-			System.out.println("getCovar(): Double.isNaN(eig.getV().get(0, 0))");
-			debugLevel = 20;
 				}
-		if (debugLevel > 3){
-			System.out.println("getCovar(): sw = "+sw +", swz = "+swz +", swx = "+swx +", swy = "+swy +", covar.det() = "+covar.det());
-			System.out.println("getCovar(): covarianvce matrix, number of used points:"+numPoints);
-			covar.print(10, 6); // w,d
-			System.out.println("getCovar(): eigenvalues");
-			eig.getD().print(10, 6); // w,d
-			System.out.println("getCovar(): eigenvectors");
-			eig.getV().print(10, 6); // w,d
+				if (all_done){
+					break;
 				}
-		if ((eig.getD().get(0, 0) == 0.0) || (Math.abs(covar.det()) < mindet)) {
-			return null; // testing with zero eigenvalue
-			// Problem with zero eigenvalue is with derivatives and coordinate conversion
+
+			}
+
+			
+			
+			
+			return null;
 		}
-		double [][][] rslt = {
-				eig.getD().getArray(),
-				eig.getV().getArray(),
+		
+		
+		// detect and split/eliminate planes that are use disconnected tiles
+		// and are not separated by higher disparity planes
+		// return planes/selections (no remove outliers!)
+		boolean[][][] filterBridges(
+				ArrayList<PlaneData> tilePlanes,
+				int max_grow_these,
+				int max_grow_far,
+				int debugLevel)
 		{
-					{sw,kz,numPoints},
-					{swz, swx, swy}}};
-		return rslt;
+			double [][] pds = new double [tilePlanes.size()][];
+			double [] max_disp = null;
+			boolean [][] selections = new boolean[pds.length][];
+			final int tsize = 4 * stSize * stSize;
+			for (int np = 0; np < pds.length; np++){
+				PlaneData pd = tilePlanes.get(np);
+				if (pd != null){
+					pds[np] = pd.getDoublePlaneDisparityStrength(
+							false, // boolean   useWorld,
+							null, // double [] window,
+							-1, // int       dir,
+							false, // boolean   use_sel,
+							false, // boolean   divide_by_area,
+							1.0, // double    scale_projection,
+							0.0, // double    fraction_uni,
+							0)[0]; // int       debugLevel)
+					if (max_disp == null){
+						max_disp = pds[np].clone(); 
+					} else {
+						for (int i = 0; i < max_disp.length; i++){
+							max_disp[i] = Math.max(max_disp[i], pds[np][i]);
 						}
-	// TODO: obsolete - remove
-	public PlaneData getPlane(
-			int [] sTileXY,
-			double []  data,
-			double []  weight,
-			boolean [] select, // null OK, will enable all tiles
-			boolean    correctDistortions,
-			boolean    preferDisparity, // Always start with disparity-most axis (false - lowest eigenvalue)
-			int        debugLevel){
-		if (select == null) {
-			select = new boolean [4*stSize];
-			for (int i = 0; i < select.length; i++) select[i] = true;
 					}
-		if (debugLevel > 0){
-			System.out.println("getPlane()");
-		}
-		double [][][] rslt = getCovar(
-				data,
-				weight,
-				select,
-				0.0,
-				debugLevel); // debugLevel1); //0); // debugLevel);
-		if (rslt == null) return null;
-		int       numPoints =  (int) rslt[2][0][2];
-		double    swc =  rslt[2][0][0];
-		double [] szxy = rslt[2][1];
-		double [][] eig_val =  rslt[0];
-		double [][] eig_vect = rslt[1];
-		// find vector most orthogonal to view // (anyway it all works with that assumption), make it first
-		// TODO normalize to local linear scales
-		int oindx = 0;
-		if (preferDisparity) {
-			for (int i = 1; i <3; i++){
-				if (Math.abs(eig_vect[0][i]) > Math.abs(eig_vect[0][oindx])){
-					oindx = i;
+					boolean [][] ms = pd.getMeasSelection(); // should not be null
+					// combine selections from all measurement layers
+					for (int ml = 0; ml < ms.length; ml++) if (ms[ml] != null){
+						if (selections[np] == null) {
+							selections[np] = ms[ml].clone(); 
+						} else {
+							for (int i = 0; i < ms[ml].length; i++) {
+								selections[np][i] |= ms[ml][i];
 							}
 						}
-		} else {
-			for (int i = 1; i < 3 ; i++){
-				if (eig_val[i][i] < eig_val[oindx][oindx]){
-					oindx = i;
 					}
 				}
 			}
-		if (eig_val[oindx][oindx] == 0.0){
-			System.out.println("getPlane(): zero eigenvalue!!");
+			class TileSelections{
+				int        np;   // number of original plane
+				boolean [] mask; // selection mask to apply to measurement layers
+				TileSelections(int nl, boolean [] mask){
+					this.np = nl;
+					this.mask = mask;
 				}
-		
-		/*
-		// Find two other axis - "mostly X" (horizontal) and "mostly Y" (vertical) 
-		int vindx = (oindx == 0)? 1 : 0;
-		int hindx = (oindx == 0)? 2 : ((oindx == 1) ? 2 : 1);
-		if (Math.abs(eig_vect[2][vindx]) < Math.abs(Math.abs(eig_vect[2][hindx]))){
-			int tmp = vindx;
-			vindx = hindx;
-			hindx = tmp;
 			}
-		*/
-		// select 2 other axes for increasing eigenvalues (so v is short axis, h  is the long one)
-		int vindx = (oindx == 0)? 1 : 0;
-		int hindx = (oindx == 0)? 2 : ((oindx == 1) ? 2 : 1);
-		if (eig_val[vindx][vindx] > eig_val[hindx][hindx]){
-			int tmp = vindx;
-			vindx = hindx;
-			hindx = tmp;
+			ArrayList<TileSelections> split_planes = new ArrayList<TileSelections>();
+			TileNeibs tileNeibs = new TileNeibs(2 * stSize, 2 * stSize);
+			HashSet<Integer> old_planes = new HashSet<Integer>();
+			for (int np = 0; np < pds.length; np++){
+				PlaneData pd = tilePlanes.get(np);
+				if (pd != null){
+					old_planes.add(np);
+					// see if all the plane tiles have other planes with higher disparity
+					boolean all_far = true;
+					for (int i = 0; i < max_disp.length; i++)if (selections[np][i]){
+						if (pds[np][i] == max_disp[i]){
+							all_far = false;
+							break;
 						}
-		
-		PlaneData pd = new PlaneData(
-				sTileXY,
-				this.tileSize,
-				this.stSize,
-				this.geometryCorrection,
-				correctDistortions);
-		pd.setZxy(szxy);
-		
-//		pd.setValues(eig_val[oindx][oindx],eig_val[hindx][hindx],eig_val[vindx][vindx]); // eigenvalues [0] - thickness, 2 other to detect skinny (poles)
-		pd.setValues(eig_val[oindx][oindx],eig_val[vindx][vindx],eig_val[hindx][hindx]); // eigenvalues [0] - thickness, 2 other to detect skinny (poles)
-
-		double [][] plane = {
-				{eig_vect[0][oindx],eig_vect[1][oindx],eig_vect[2][oindx]},  // plane normal to camera
-				{eig_vect[0][vindx],eig_vect[1][vindx],eig_vect[2][vindx]},  // "horizontal" axis // to detect skinny planes and poles
-				{eig_vect[0][hindx],eig_vect[1][hindx],eig_vect[2][hindx]}}; //  "vertical"   axis  // to detect skinny planes and poles
-		/*
-		// Make normal be towards camera (positive disparity), next vector - positive in X direction (right), last one - in positive Y (down)
-		for (int v = 0; v <3; v++) {
-			if (plane[v][v] < 0.0) for (int i = 0; i < 3; i ++) plane[v][i] = -plane[v][i];
 					}
-		*/
-		// Make normal be towards camera (positive disparity), next vector - positive in X direction (right)
-		for (int v = 0; v < 2; v++) {
-			if (plane[v][v] < 0.0) for (int i = 0; i < 3; i ++) plane[v][i] = -plane[v][i];
+					if (!all_far){
+						boolean [] grown_sel = selections[np].clone();
+						tileNeibs.growSelection( //
+								max_grow_these,  // int        grow,           // grow tile selection by 1 over non-background tiles 1: 4 directions, 2 - 8 directions, 3 - 8 by 1, 4 by 1 more
+								grown_sel, // boolean [] tiles,
+								null);      // boolean [] prohibit,
+						int [] clusters = tileNeibs.enumerateClusters(
+								grown_sel, // boolean [] tiles,
+								false); // boolean ordered)
+
+						int num_clusters = tileNeibs.getMax(
+								clusters); // int [] data)
+						if (num_clusters > 1){
+							boolean [] dbg_grown_sel =grown_sel.clone();
+							// determine area that should have connections between these clusters, not blocked by known
+							// tiles with lower disparity
+							boolean [] bound_octo= tileNeibs.boundShape(
+									selections[np], // boolean [] selection,
+									true); // boolean octo)
+							// now combine all "offending" selections - known tiles belonging to farther planes
+							boolean [] sel_offend = new boolean[tsize];
+							for (int np_other = 0; np_other < pds.length; np_other++) if ((selections[np_other] != null) && (np_other != np)){
+								for (int i = 0; i < tsize; i++) if (selections[np_other][i]){
+									sel_offend [i] |= (pds[np][i] > pds[np_other][i]);
+								}
+							}
+							// now grow offending selections to fill gaps
+							boolean [] dbg_sel_offend = sel_offend.clone(); 
+							tileNeibs.growSelection( //
+									max_grow_far,  // int        grow,           // grow tile selection by 1 over non-background tiles 1: 4 directions, 2 - 8 directions, 3 - 8 by 1, 4 by 1 more
+									sel_offend, // boolean [] tiles,
+									null);      // boolean [] prohibit,
+							// create prohibited mask - all not in grown selection that are outside bounding octagon or belong to grown offenders
+							
+							boolean [] prohibit = new boolean [tsize];
+							for (int i = 0; i < prohibit.length; i++) {
+								prohibit[i] = !grown_sel[i] && (!bound_octo[i] || sel_offend[i] ); // (pds[np][i] > max_disp[i]); 
+							}
+							tileNeibs.growSelection( //
+									stSize,     // int        grow,           // grow tile selection by 1 over non-background tiles 1: 4 directions, 2 - 8 directions, 3 - 8 by 1, 4 by 1 more
+									grown_sel,  // boolean [] tiles,
+									prohibit);      // boolean [] prohibit,
+							// re-clusterize what is remaining
+							int [] dbg_clusters = clusters.clone();
+							clusters = tileNeibs.enumerateClusters(
+									grown_sel, // boolean [] tiles,
+									false); // boolean ordered)
+							num_clusters = tileNeibs.getMax(
+									clusters); 
+							if (num_clusters > 1){
+								boolean [][] cluster_sels = new boolean [num_clusters][selections[np].length];
+								for (int i = 0; i < selections[np].length; i++) if (selections[np][i]){
+									cluster_sels[clusters[i]-1][i] = true;
+								}
+								for (int nc = 0; nc < num_clusters; nc++){
+									split_planes.add(new TileSelections(np, cluster_sels[nc]));
+								}
+								old_planes.remove(np);
+							}
+							if (debugLevel > 2) {
+								String [] dbg_titles = {"sel","grown", "clust0","bound","offend","off_grown", "prohibit", "clusters","clus_masked"};
+								double [][] dbg_img =  new double [dbg_titles.length][tsize];
+								double lon = 3.0;
+								for (int i = 0; i < tsize; i++) {
+									dbg_img[0][i] = selections[np][i] ? lon: 0.0; 
+									dbg_img[1][i] = dbg_grown_sel[i] ? lon: 0.0; 
+									dbg_img[2][i] = dbg_clusters[i]; 
+									dbg_img[3][i] = bound_octo[i] ? lon: 0.0; 
+									dbg_img[4][i] = dbg_sel_offend[i] ? lon: 0.0; 
+									dbg_img[5][i] = sel_offend[i] ? lon: 0.0; 
+									dbg_img[6][i] = prohibit[i] ? lon: 0.0; 
+									dbg_img[7][i] = clusters[i]; 
+									dbg_img[8][i] = selections[np][i] ? clusters[i] : 0.0; 
 								}
 								
-		// make  direction last vector so px (x) py (.) disp < 0 (left-hand coordinate system) 
-		if (new Matrix(plane).det() > 0){
-			for (int i = 0; i < 3; i ++) plane[2][i] = -plane[2][i];
+								showDoubleFloatArrays sdfa_instance = new showDoubleFloatArrays();
+								sdfa_instance.showArrays(dbg_img, 2 * superTileSize, 2* superTileSize, true, "bridges-"+np+"-"+debugLevel,dbg_titles);
+								sdfa_instance.showArrays(pds,     2 * superTileSize, 2* superTileSize, true, "pds-bridges-"+np+"-"+debugLevel,dbg_titles);
 							}
-		
-		pd.setVectors   (plane);
-		pd.setNumPoints (numPoints);
-		pd.setWeight    (swc);
-		pd.setPlaneSelection(select);
-		return pd;
 						}
-	
-	public PlaneData removePlaneOutliers(
-			PlaneData  pd,     // already found or null 
-			int [] sTileXY,    // may be null if pd is not null 
-			double []  data,
-			double []  weight,
-			boolean [] select, // will be modified
-			boolean    correctDistortions,
-			double     targetEigen, // target eigenvalue for primary axis (is disparity-dependent, so is non-constant)
-			int        maxRemoved,  // maximal number of tiles to remove (not a constant)
-			int        minLeft,     // minimal number of tiles to keep
-			boolean    preferDisparity, // Always start with disparity-most axis (false - lowest eigenvalue)
-			int        debugLevel){
-		int stSize2 = 2 * stSize;
-		if (pd == null) {
-			pd = getPlane(
-					sTileXY, 
-					data,
-					weight,
-					select, // null OK
-					correctDistortions,
-					preferDisparity,
-					debugLevel);
-		} else if (select != null){
-			pd.setPlaneSelection(select);
 					}
-		if (pd == null){
-			return null; // zero eigenvalues
 				}
-		if (maxRemoved > (pd.getNumPoints() - minLeft)) maxRemoved = pd.getNumPoints() - minLeft;
-		int numRemoved = 0;
-		for (; (pd.getValue() > targetEigen) && (numRemoved < maxRemoved); numRemoved++){
-			if (debugLevel > 2){
-				System.out.println("removePlaneOutliers("+sTileXY[0]+":"+sTileXY[1]+"): numRemoved = "+numRemoved+" eigenValue = "+pd.getValue()+" target = "+targetEigen);
 			}
-			// make a plane and find the worst (largest disparity difference) tile
-			// z = -(x*Vx + y*Vy)/Vz
-			double worst_d2 = 0.0;
-			int worst_index =  -1;
-			double [] v = pd.getVector();
-			double [] zxy0 = pd.getZxy();
-			for (int indx = 0; indx < data.length; indx++){
-				if (select[indx] && (weight[indx] > 0)){
-//					double w = weight[indx];
-					double x = ((indx % stSize2) - stSize) - zxy0[1];
-					double y = ((indx / stSize2) - stSize) - zxy0[2];
-					double d = data[indx];
-					d -= zxy0[0];
-					d += (x * v[1]+y*v[2])/v[0];
-					double d2 = d*d;
-					if (d2 > worst_d2){
-						worst_d2 = d2;
-						worst_index = indx;
+			if (!split_planes.isEmpty()){
+				boolean [][][] split_selections = new boolean [old_planes.size() + split_planes.size()][][];
+				int ns = 0;
+				for (Integer np: old_planes){
+					split_selections[ns++] = tilePlanes.get(np).getMeasSelection();
 				}
+				for (TileSelections ts: split_planes){
+					PlaneData pd = tilePlanes.get(ts.np);
+					boolean [][] ms = tilePlanes.get(ts.np).getMeasSelection().clone();
+					for (int ml = 0; ml < ms.length; ml++) if (ms[ml] != null){
+						ms[ml] = ms[ml].clone();
+						for (int i = 0; i < ms[ml].length; i++){
+							ms[ml][i] &= ts.mask[i];
 						}
 					}
-			if (worst_index < 0) {
-				System.out.println("This is a BUG in removePlaneOutliers()");
-				break;
+					split_selections[ns++] = ms;
 				}
-			select[worst_index] = false;
-			if (debugLevel > 2){
-				System.out.println("removePlaneOutliers() worst_index = " + worst_index);
+				return split_selections;
+				
+			}
+			return null; // no changes
 		}
 
-			pd = getPlane(    // re-calculate with point removed
-					pd.getSTileXY(),
-					data,
-					weight,
-					select,
-					correctDistortions,
-					preferDisparity,
-					debugLevel);
-			if (pd == null) {
-				return null;
+		public double [] getWindow (
+				int size)
+		{
+			double [] wnd1d = new double [size];
+			for (int i = 0; i < size/2; i++){
+				wnd1d[i] = 0.5 * (1.0 - Math.cos(2*Math.PI*(i+0.5)/size));
+				wnd1d[size - i -1] = wnd1d[i]; 
+			}
+			double [] wnd = new double [size * size];
+			int indx = 0;
+			for (int i = 0; i < size; i++){
+				for (int j = 0; j < size; j++){
+					wnd[indx++] = wnd1d[i]*wnd1d[j];
 				}
 			}
-		return pd;
+			return wnd;
 		}
 		
+		
+		
+		
+		
+
+	}
 }

src/main/java/TileProcessor.java

@@ -3386,7 +3386,7 @@ public class TileProcessor {
 			LinkPlanes lp = new LinkPlanes (clt_parameters, st);
 			// try to merge multiple times
 			
-			int max_num_merge_try = 5;
+			int max_num_merge_try = 8;
 			TilePlanes.PlaneData [][][] dbg_orig_planes = new TilePlanes.PlaneData [max_num_merge_try][][];
 			
 			for (int num_merge_try = 0; num_merge_try < max_num_merge_try; num_merge_try++){
@@ -3467,20 +3467,6 @@ public class TileProcessor {
 						clt_parameters.tileX,
 						clt_parameters.tileY);
 //				System.out.println("merge_cost_data.length = " + merge_cost_data.length);
-				/*
-//				double  [][][][][][]  merge_cost_data_eq = 
-						lp.costSameTileConnections(
-						true, // final boolean ignore_weights,
-						1.8,  // final double threshold_worst,
-						1000.0, //final double threshold_world_worst,
-						st.planes, // ffinal TilePlanes.PlaneData [][] planes,
-						merge_candidates,       // final int [][][] merge_candidates,
-						plane_nooverlaps, // final boolean [][][]   valid_candidates, // will be updated
-						2,                      // -1, // debugLevel,                  // final int        debugLevel)
-						clt_parameters.tileX,
-						clt_parameters.tileY);
-//				System.out.println("merge_cost_data_eq.length = " + merge_cost_data_eq.length);
-				*/
 				int [][][] merge_groups = lp.extractMergeSameTileGroups(
 						st.planes,              // final TilePlanes.PlaneData [][] planes,
 						merge_candidates,       // final int [][][] merge_candidates,

src/main/java/TileSurface.java

@@ -235,144 +235,6 @@ public class TileSurface {
 			
 		}
 
-		public class TileNeibs{
-			int sizeX;
-			int sizeY;
-			public int dirs = 8;
-			public TileNeibs(int size){
-				this.sizeX = size;
-				this.sizeY = size;
-			}
-			public TileNeibs(int sizeX, int sizeY){
-				this.sizeX = sizeX;
-				this.sizeY = sizeY;
-			}
-			public int numNeibs() // TODO: make configurable to 
-			{
-				return dirs;
-			}
-			public int opposite(int dir){
-				return (dir + dirs / 2) % dirs;
-			}
-
-			/**
-			 * Get x,y pair from index
-			 * @param indx element index
-			 * @return array of {x,y}
-			 */
-			
-			int [] getXY(int indx)
-			{
-				int [] xy = {indx % sizeX ,indx / sizeX};
-				return xy;
-			}
-			
-			/**
-			 * Get element index from x and y
-			 * @param x horizontal position
-			 * @param y vertical position
-			 * @return element linescan index
-			 */
-			
-			int getIndex(int x, int y){
-				return y * sizeX + x;
-			}
-			
-			/**
-			 * Get 2d element index after step N, NE, ... NW. Returns -1 if leaving array   
-			 * @param indx start index
-			 * @param dir step direction (CW from up)
-			 * @return new index or -1 if leaving 
-			 */
-			int getNeibIndex(int indx, int dir)
-			{
-				int y = indx / sizeX;
-				int x = indx % sizeX;
-				if (dir < 0) return indx;
-				switch (dir % dirs){
-				case 0: return (y == 0) ?                                    -1 : (indx - sizeX); 
-				case 1: return ((y == 0)           || ( x == (sizeX - 1))) ? -1 : (indx - sizeX + 1); 
-				case 2: return (                      ( x == (sizeX - 1))) ? -1 : (indx        + 1); 
-				case 3: return ((y == (sizeY - 1)) || ( x == (sizeX - 1))) ? -1 : (indx + sizeX + 1); 
-				case 4: return ((y == (sizeY - 1))                       ) ? -1 : (indx + sizeX); 
-				case 5: return ((y == (sizeY - 1)) || ( x == 0))           ? -1 : (indx + sizeX - 1); 
-				case 6: return (                      ( x == 0))           ? -1 : (indx        - 1); 
-				case 7: return ((y == 0)           || ( x == 0))           ? -1 : (indx - sizeX - 1); 
-				default: return indx;
-				}
-			}
-			
-			/**
-			 * Return tile segment for 50% overlap. -1 - center, 0 N, 1 - NE,... 7 - NW
-			 * @param indx element index
-			 * @return which of the 9 areas this element belongs
-			 */
-			int getSegment(int indx)
-			{
-				int s1x = sizeX / 4;
-				int s3x = 3 * sizeX / 4;
-				int s1y = sizeY / 4;
-				int s3y = 3 * sizeY / 4;
-				int x = indx % sizeX;
-				int y = indx / sizeX;
-				boolean up = y < s1y;
-				boolean down = y >= s3y;
-				boolean left = x < s1x;
-				boolean right = x >= s3x;
-				if (up){
-					if (left) return 7;
-					if (right) return 1;
-					return 0;
-				}
-				if (down){
-					if (left) return 5;
-					if (right) return 3;
-					return 4;
-				}
-				if (left) return 6;
-				if (right) return 2;
-				return -1;
-			}
-
-			/**
-			 * Find if the step leaves the center half of all area
-			 * @param indx start point
-			 * @param dir direction
-			 * @return direction to the new tile (assuming 50% overlap) or -1 if did not cross the border
-			 */
-			int leaveOvderlapedCenter(int indx, int dir)
-			{
-				int segm = getSegment(indx);
-				int indx1 = getNeibIndex(indx, dir);
-				if (indx1 < 0 ) {
-					return -1; // should not happen
-				}
-				int segm1 = getSegment(indx1);
-				if (segm1 == segm) return -1;
-				if (segm == -1) return segm1;
-				int [][] dxy = {
-						{ 0, 0},
-						{-1, 0},
-						{-1, 1},
-						{ 0, 1},
-						{ 1, 1},
-						{ 1, 0},
-						{ 1,-1},
-						{ 0,-1},
-						{-1,-1}};
-				int dx = dxy[segm1 + 1][1] - dxy[segm + 1][1]; 
-				int dy = dxy[segm1 + 1][0] - dxy[segm + 1][0];
-				for (int dp1 = 0; dp1 <=8; dp1++) {
-					int sdx = (dx > 0) ? 1: ( (dx < 0) ? -1 : 0);
-					int sdy = (dy > 0) ? 1: ( (dy < 0) ? -1 : 0);
-					if ((dxy[dp1][0] == sdy) && (dxy[dp1][1] == sdx)){
-						return dp1 -1;
-					}
-				}
-				return -1; // should not happen
-			}
-		}
-
 		public int getNStileDir(
 				int nsTile,
 				int dir)