assigning tiles to surfaces

src/main/java/EyesisCorrectionParameters.java

@@ -2215,6 +2215,18 @@ public class EyesisCorrectionParameters {
   		public double     msScaleProj          =   1.5;  // Scale projection of the plane ellipsoid 
   		public double     msFractUni           =   0.3;  // Spread this fraction of the ellipsoid weight among extended (double) supertile
   		
+  		public double     tsMaxDiff            = 0.3;    // Maximal disparity difference when assigning tiles
+		public double     tsMinDiffOther       = 0.35;   // Minimal disparity difference to be considered as a competitor surface
+		public double     tsMinStrength        = 0.05;   // Minimal tile correlation strength to be assigned
+		public double     tsMaxStrength        = 10.0;   // Maximal tile correlation strength to be assigned
+		public int        tsMoveDirs           = 3;      // Allowed tile disparity correction: 1 increase, 2 - decrease, 3 - both directions
+		public boolean    tsEnMulti            = false;  // Allow assignment when several surfaces fit
+		public double     tsSurfStrPow         = 0.0;    // Raise surface strengths ratio to this power when comparing candidates
+		public double     tsSigma              = 2.0;    // Radius of influence (in tiles) of the previously assigned tiles
+		public double     tsNSigma             = 2.0;    // Maximal relative to radius distance to calculate influence
+		public double     tsMinAdvantage       = 3.0;    // Minimal ratio of the best surface candidate to the next one to make selection
+		public boolean    tsReset              = false;  // Reset tiles to surfaces assignment
+  		
   		public boolean    replaceWeakOutlayers =   true; // false; 
   		
   		public boolean    dbg_migrate =            true; 
@@ -2519,6 +2531,17 @@ public class EyesisCorrectionParameters {
 			properties.setProperty(prefix+"msScaleProj",      this.msScaleProj +"");
 			properties.setProperty(prefix+"msFractUni",       this.msFractUni +"");
 
+			properties.setProperty(prefix+"tsMaxDiff",        this.tsMaxDiff +"");
+			properties.setProperty(prefix+"tsMinDiffOther",   this.tsMinDiffOther +"");
+			properties.setProperty(prefix+"tsMinStrength",    this.tsMinStrength +"");
+			properties.setProperty(prefix+"tsMaxStrength",    this.tsMaxStrength +"");
+  			properties.setProperty(prefix+"tsMoveDirs",       this.tsMoveDirs+"");
+			properties.setProperty(prefix+"tsEnMulti",        this.tsEnMulti+"");
+			properties.setProperty(prefix+"tsSurfStrPow",    this.tsSurfStrPow +"");
+			properties.setProperty(prefix+"tsSigma",         this.tsSigma +"");
+			properties.setProperty(prefix+"tsNSigma",         this.tsNSigma +"");
+			properties.setProperty(prefix+"tsMinAdvantage",  this.tsMinAdvantage +"");
+			
 			properties.setProperty(prefix+"dbg_migrate",            this.dbg_migrate+"");
   			
 			properties.setProperty(prefix+"show_ortho_combine",     this.show_ortho_combine+"");
@@ -2817,6 +2840,18 @@ public class EyesisCorrectionParameters {
   			if (properties.getProperty(prefix+"msFractUni")!=null)        this.msFractUni=Double.parseDouble(properties.getProperty(prefix+"msFractUni"));
 
 
+  			if (properties.getProperty(prefix+"tsMaxDiff")!=null)         this.tsMaxDiff=Double.parseDouble(properties.getProperty(prefix+"tsMaxDiff"));
+  			if (properties.getProperty(prefix+"tsMinDiffOther")!=null)    this.tsMinDiffOther=Double.parseDouble(properties.getProperty(prefix+"tsMinDiffOther"));
+  			if (properties.getProperty(prefix+"tsMinStrength")!=null)     this.tsMinStrength=Double.parseDouble(properties.getProperty(prefix+"tsMinStrength"));
+  			if (properties.getProperty(prefix+"tsMaxStrength")!=null)     this.tsMaxStrength=Double.parseDouble(properties.getProperty(prefix+"tsMaxStrength"));
+  			if (properties.getProperty(prefix+"tsMoveDirs")!=null)        this.tsMoveDirs=Integer.parseInt(properties.getProperty(prefix+"tsMoveDirs"));
+  			if (properties.getProperty(prefix+"tsEnMulti")!=null)         this.tsEnMulti=Boolean.parseBoolean(properties.getProperty(prefix+"tsEnMulti"));
+  			if (properties.getProperty(prefix+"tsSurfStrPow")!=null)     this.tsSurfStrPow=Double.parseDouble(properties.getProperty(prefix+"tsSurfStrPow"));
+  			if (properties.getProperty(prefix+"tsSigma")!=null)          this.tsSigma=Double.parseDouble(properties.getProperty(prefix+"tsSigma"));
+  			if (properties.getProperty(prefix+"tsNSigma")!=null)          this.tsNSigma=Double.parseDouble(properties.getProperty(prefix+"tsNSigma"));
+  			if (properties.getProperty(prefix+"tsMinAdvantage")!=null)   this.tsMinAdvantage=Double.parseDouble(properties.getProperty(prefix+"tsMinAdvantage"));
+  			
+  			
   			if (properties.getProperty(prefix+"dbg_migrate")!=null)       this.dbg_migrate=Boolean.parseBoolean(properties.getProperty(prefix+"dbg_migrate"));
 
   			if (properties.getProperty(prefix+"show_ortho_combine")!=null)     this.show_ortho_combine=Boolean.parseBoolean(properties.getProperty(prefix+"show_ortho_combine"));
@@ -3139,6 +3174,17 @@ public class EyesisCorrectionParameters {
   			gd.addNumericField("Scale projection of the plane ellipsoid",                                      this.msScaleProj,  6);
   			gd.addNumericField("Spread this fraction of the ellipsoid weight among extended (double) supertile",this.msFractUni,  6);
   			
+  			gd.addNumericField("Maximal disparity difference when assigning tiles",                               this.tsMaxDiff,  6);
+  			gd.addNumericField("Minimal disparity difference to be considered as a competitor surface",           this.tsMinDiffOther,  6);
+  			gd.addNumericField("Minimal tile correlation strength to be assigned",                                this.tsMinStrength,  6);
+  			gd.addNumericField("Maximal tile correlation strength to be assigned",                                this.tsMaxStrength,  6);
+  			gd.addNumericField("Allowed tile disparity correction: 1 increase, 2 - decrease, 3 - both directions",this.tsMoveDirs,  0);
+  			gd.addCheckbox    ("Allow assignment when several surfaces fit",                                      this.tsEnMulti);
+  			gd.addNumericField("Raise surface strengths ratio to this power when comparing candidates",           this.tsSurfStrPow,  6);
+  			gd.addNumericField("Radius of influence (in tiles) of the previously assigned tiles",                 this.tsSigma,  6);
+  			gd.addNumericField("Maximal relative to radius distance to calculate influence",                      this.tsNSigma,  6);
+  			gd.addNumericField("Minimal ratio of the best surface candidate to the next one to make selection",   this.tsMinAdvantage,  6);
+  			
   			gd.addCheckbox    ("Test new mode after migration",                                                this.dbg_migrate);
 
   			gd.addMessage     ("--- Other debug images ---");
@@ -3443,6 +3489,18 @@ public class EyesisCorrectionParameters {
   			this.msScaleProj=           gd.getNextNumber();
   			this.msFractUni=            gd.getNextNumber();
 
+  			this.tsMaxDiff=             gd.getNextNumber();
+  			this.tsMinDiffOther=        gd.getNextNumber();
+  			this.tsMinStrength=         gd.getNextNumber();
+  			this.tsMaxStrength=         gd.getNextNumber();
+  			this.tsMoveDirs=      (int) gd.getNextNumber();
+  			this.tsEnMulti=             gd.getNextBoolean();
+  			this.tsSurfStrPow=          gd.getNextNumber();
+  			this.tsSigma=               gd.getNextNumber();
+  			this.tsNSigma=              gd.getNextNumber();
+  			this.tsMinAdvantage=        gd.getNextNumber();
+  			this.tsReset              = false;  // Reset tiles to surfaces assignment
+  			
   			this.dbg_migrate=           gd.getNextBoolean();
 
   			this.show_ortho_combine=    gd.getNextBoolean();
@@ -3458,8 +3516,40 @@ public class EyesisCorrectionParameters {
   			this.show_flaps_dirs=       gd.getNextBoolean();
   			this.show_first_clusters=   gd.getNextBoolean();
   			this.show_planes=           gd.getNextBoolean();
+  			
   			return true;
   		}
+  		
+  		public boolean showTsDialog() {
+  			GenericDialog gd = new GenericDialog("Set CLT parameters");
+  			gd.addNumericField("Maximal disparity difference when assigning tiles",                               this.tsMaxDiff,  6);
+  			gd.addNumericField("Minimal disparity difference to be considered as a competitor surface",           this.tsMinDiffOther,  6);
+  			gd.addNumericField("Minimal tile correlation strength to be assigned",                                this.tsMinStrength,  6);
+  			gd.addNumericField("Maximal tile correlation strength to be assigned",                                this.tsMaxStrength,  6);
+  			gd.addNumericField("Allowed tile disparity correction: 1 increase, 2 - decrease, 3 - both directions",this.tsMoveDirs,  0);
+  			gd.addCheckbox    ("Allow assignment when several surfaces fit",                                      this.tsEnMulti);
+  			gd.addNumericField("Raise surface strengths ratio to this power when comparing candidates",           this.tsSurfStrPow,  6);
+  			gd.addNumericField("Radius of influence (in tiles) of the previously assigned tiles",                 this.tsSigma,  6);
+  			gd.addNumericField("Maximal relative to radius distance to calculate influence",                      this.tsNSigma,  6);
+  			gd.addNumericField("Minimal ratio of the best surface candidate to the next one to make selection",   this.tsMinAdvantage,  6);
+  			gd.addCheckbox    ("Reset tiles to surfaces assignment",                                              false);
+  			WindowTools.addScrollBars(gd);
+  			gd.showDialog();
+  			if (gd.wasCanceled()) return false;
+  			this.tsMaxDiff=             gd.getNextNumber();
+  			this.tsMinDiffOther=        gd.getNextNumber();
+  			this.tsMinStrength=         gd.getNextNumber();
+  			this.tsMaxStrength=         gd.getNextNumber();
+  			this.tsMoveDirs=      (int) gd.getNextNumber();
+  			this.tsEnMulti=             gd.getNextBoolean();
+  			this.tsSurfStrPow=          gd.getNextNumber();
+  			this.tsSigma=               gd.getNextNumber();
+  			this.tsNSigma=              gd.getNextNumber();
+  			this.tsMinAdvantage=        gd.getNextNumber();
+  			this.tsReset=               gd.getNextBoolean();
+  			return true;
+  		}
+  		
     }
 
     public static class DCTParameters {

src/main/java/Eyesis_Correction.java

@@ -507,6 +507,7 @@ private Panel panel1,
 			addButton("CLT reset 3D",              panelClt1, color_stop);
 			addButton("CLT 3D",                    panelClt1, color_conf_process);
 			addButton("CLT planes",                panelClt1, color_conf_process);
+			addButton("CLT ASSIGN",                panelClt1, color_process);
 			addButton("CLT OUT 3D",                panelClt1, color_process);
 						
 			add(panelClt1);
@@ -4782,6 +4783,26 @@ private Panel panel1,
         		UPDATE_STATUS, //final boolean    updateStatus,
         		DEBUG_LEVEL); //final int        debugLevel);
         return;
+    } else if (label.equals("CLT ASSIGN")) {
+    	DEBUG_LEVEL=MASTER_DEBUG_LEVEL;
+    	EYESIS_CORRECTIONS.setDebug(DEBUG_LEVEL);
+        if (QUAD_CLT == null){
+       		System.out.println("QUAD_CLT is null, nothing to show");
+       		return;
+        }
+
+    	if (!CLT_PARAMETERS.showTsDialog()) return;
+        
+        boolean OK = QUAD_CLT.assignCLTPlanes(
+        		CLT_PARAMETERS,  // EyesisCorrectionParameters.DCTParameters           dct_parameters,
+        		THREADS_MAX, //final int          threadsMax,  // maximal number of threads to launch                         
+        		UPDATE_STATUS, //final boolean    updateStatus,
+        		DEBUG_LEVEL); //final int        debugLevel);
+        if (!OK){
+       		System.out.println("Could not assign tiles to surfaces, probably \"CLT planes\" command did not run");
+       		return;
+        }
+        return;
         
     } else if (label.equals("CLT OUT 3D")) {
     	DEBUG_LEVEL=MASTER_DEBUG_LEVEL;
@@ -4791,7 +4812,6 @@ private Panel panel1,
        		return;
         }
         
-        
     	String configPath=null;
     	if (EYESIS_CORRECTIONS.correctionsParameters.saveSettings) {
     		configPath=EYESIS_CORRECTIONS.correctionsParameters.selectResultsDirectory(

src/main/java/QuadCLT.java

@@ -4460,7 +4460,6 @@ public class QuadCLT {
 		  }
 		  return rslt;
 	  }
-//	  public ImagePlus [] cltDisparityScan(
 	  public void showCLTPlanes(
 			  EyesisCorrectionParameters.CLTParameters           clt_parameters,
 			  final int          threadsMax,  // maximal number of threads to launch                         
@@ -4481,9 +4480,33 @@ public class QuadCLT {
 		  			threadsMax,
 		  			updateStatus,
 		  			debugLevel);
-//		  	CLTPass3d last_scan = tp.clt_3d_passes.get(tp.clt_3d_passes.size() -1); // get last one
 		  	
 	  }
+	  
+	  public boolean assignCLTPlanes(
+			  EyesisCorrectionParameters.CLTParameters           clt_parameters,
+			  final int          threadsMax,  // maximal number of threads to launch                         
+			  final boolean    updateStatus,
+			  final int        debugLevel)
+	  {
+		  	if (tp == null){
+		  		System.out.println("showCLTPlanes(): tp is null");
+		  		return false;
+		  	}
+		  	if (tp.clt_3d_passes == null){
+		  		System.out.println("showCLTPlanes(): tp.clt_3d_passes is null");
+		  		return false;
+		  	}
+		  	return tp.assignTilesToSurfaces(
+		  			clt_parameters,
+		  			geometryCorrection,		  			
+		  			threadsMax,
+		  			updateStatus,
+		  			debugLevel);
+		  	
+	  }
+	  
+	  
 	  public void out3d(
 			  EyesisCorrectionParameters.CLTParameters           clt_parameters,
 			  final int          threadsMax,  // maximal number of threads to launch                         

src/main/java/SuperTiles.java

@@ -73,6 +73,9 @@ public class SuperTiles{
 	int [][]    shell_map = null;   // per supertile, per disparity plane - shell index + 1 (0 - none)
 	double [][] surfaces;  // per shell, per tile (linescan order) disparity value or NaN in missing supertiles
 	
+	TileSurface tileSurface = null;	
+	
+
 	/**
 	 * currently lowest plane for each includes all tiles, so do not use it. May change in the future
 	 * @param np total number of planes in a supertile
@@ -185,6 +188,15 @@ public class SuperTiles{
 		}
 	}
 	
+	public void setTileSurface( TileSurface tileSurface)
+	{
+		this.tileSurface = tileSurface;
+	}
+	public TileSurface getTileSurface()
+	{
+		return tileSurface;
+	}
+	
 	public void initFuseCoeff(
 			double scale_diag,
 			boolean debug)

src/main/java/TileProcessor.java

@@ -2942,6 +2942,53 @@ public class TileProcessor {
 	}
 
 
+	public boolean assignTilesToSurfaces(
+			EyesisCorrectionParameters.CLTParameters           clt_parameters,
+			GeometryCorrection geometryCorrection,
+			final int         threadsMax,  // maximal number of threads to launch                         
+			final boolean     updateStatus,
+			final int         debugLevel)
+	{
+		CLTPass3d scan_prev = clt_3d_passes.get(clt_3d_passes.size() -1); // get last one
+//		boolean show_st =    clt_parameters.stShow || (debugLevel > 1);
+		SuperTiles st = scan_prev.getSuperTiles();
+		TileSurface tileSurface = st.getTileSurface();
+		if (tileSurface == null){
+			return false;
+		}
+
+		double [][][]  dispStrength = st.getDisparityStrengths(
+				clt_parameters.stMeasSel); // int        stMeasSel) //            = 1;      // Select measurements for supertiles : +1 - combo, +2 - quad +4 - hor +8 - vert)
+
+		boolean [][] tileSel =  st.getMeasurementSelections(
+				clt_parameters.stMeasSel); // int        stMeasSel) //            = 1;      // Select measurements for supertiles : +1 - combo, +2 - quad +4 - hor +8 - vert)
+
+		tileSurface.InitTilesAssignment(
+				clt_parameters.tsReset,
+				dispStrength, // final double [][][]                            dispStrength,
+				tileSel,      // final boolean [][]                             tileSel,
+                debugLevel);    // final int                                    debugLevel,
+
+		 int [] stats=  tileSurface.assignTilesToSurfaces(
+				clt_parameters.tsMaxDiff,       //final double        maxDiff,
+				clt_parameters.tsMinDiffOther,  //final double        minDiffOther, // should be >= maxDiff
+				clt_parameters.tsMinStrength,   //final double        minStrength,
+				clt_parameters.tsMaxStrength,   //final double        maxStrength,
+				clt_parameters.tsMoveDirs,      //final int           moveDirs, // 1 increase disparity, 2 - decrease disparity, 3 - both directions
+				clt_parameters.tsEnMulti,       //final boolean       enMulti,
+				clt_parameters.tsSurfStrPow,    //final double        surfStrPow, // surface strength power
+				clt_parameters.tsSigma,         //final double        sigma,
+				clt_parameters.tsNSigma,        //final double        nSigma,
+				clt_parameters.tsMinAdvantage,  //final double        minAdvantage,
+				clt_parameters.plDispNorm,      // final double        dispNorm, // disparity normalize (proportionally scale down disparity difference if above
+				dispStrength,                   // final double [][][] dispStrength,
+				0, // -1,                       // debugLevel,                  // final int        debugLevel)
+				clt_parameters.tileX,
+				clt_parameters.tileY);
+		 tileSurface.printStats(stats);
+		return true;
+	}	
+	
 //======================
 	public void showPlanes(
 			EyesisCorrectionParameters.CLTParameters           clt_parameters,
@@ -3174,9 +3221,11 @@ public class TileProcessor {
 				st.tileProcessor.getTilesY(),        // int tilesY,
 				geometryCorrection,               // GeometryCorrection geometryCorrection,
 				st.tileProcessor.threadsMax);        // int threadsMax); 
+		st.setTileSurface(tileSurface);
 		
 		
-		TileSurface.TileData [][] tileData = tileSurface.createTileShells (
+//		TileSurface.TileData [][] tileData = 
+				tileSurface.createTileShells (
 				clt_parameters.msUseSel,            // final boolean                   use_sel,
 				clt_parameters.msDivideByArea,      // final boolean                   divide_by_area,
 				clt_parameters.msScaleProj,         // final double                    scale_projection,
@@ -3186,16 +3235,14 @@ public class TileProcessor {
 				clt_parameters.tileX,
 				clt_parameters.tileY);
 
-		int [][] tiles_layers = tileSurface.sortTilesToSurfaces(
+		tileSurface.InitTilesAssignment(
+				true,
 				dispStrength, // final double [][][]                            dispStrength,
 				tileSel,      // final boolean [][]                             tileSel,
-				tileData,     // final TileData [][]                            tileData_src,
-				// parameters
-				clt_parameters, // final EyesisCorrectionParameters.CLTParameters clt_parameters,
-                debugLevel,    // final int                                      debugLevel,
-				clt_parameters.tileX,
-				clt_parameters.tileY);
-		
+                debugLevel);    // final int                                      debugLevel,
+		if (debugLevel > -10){
+			return; // just cut off the rest
+		}
 		
 
 		
@@ -3572,6 +3619,10 @@ public class TileProcessor {
 		}
 	}
 
+	
+	
+	
+	
 	public void secondPassSetup( // prepare tile tasks for the second pass based on the previous one(s)
 			//			  final double [][][]       image_data, // first index - number of image in a quad
 			EyesisCorrectionParameters.CLTParameters           clt_parameters,

src/main/java/TileSurface.java

+import java.awt.Point;
 import java.util.ArrayList;
 import java.util.Collections;
 import java.util.Comparator;
@@ -40,6 +41,8 @@ public class TileSurface {
 		private int [] ss_dirs8;
 		private double [] window;
 		private int threadsMax = 100;
+		private int [][] tileLayers = null;
+		private TileData [][] tileData = null;
 		
 		static int STAT_UNASSIGNED = 0; // index of number of unassigned tiles
 		static int STAT_ASSIGNED =   1; // index of number of assigned tiles
@@ -1649,7 +1652,7 @@ public class TileSurface {
 					dbg_Y);            // final int                       dbg_Y);
 			
 			showSurfaceDS (tileData, "tileData");			
-			
+			this.tileData = tileData;
 			return tileData;
 		}
 		
@@ -1712,7 +1715,8 @@ public class TileSurface {
 		}
 		
 		/**
-		 * Convert from image tile index to surface tile index (surface tiles are all full superTileSize
+		 * Convert from image tile index to the surface tile index (surface tiles are all
+		 * full superTileSize),
 		 * TODO: update/remove if surface grid will be trimmed to fit image
 		 * Currently there are 324 tiles horizontally in the image and 328 in the surfaces
 		 * @param nTile image tile index in scan order
@@ -1722,9 +1726,31 @@ public class TileSurface {
 				int nTile)
 		{
 			// calculate index in tileData (has different dimensions - TODO: trim?
-			return stilesX * superTileSize * (nTile /  imageTilesX) + (nTile %  imageTilesX);
+			int surfaceTilesX = stilesX * superTileSize;
+			return surfaceTilesX * (nTile /  imageTilesX) + (nTile %  imageTilesX);
+		}
 
+		/**
+		 * Convert from  surface tile index (surface tiles are all full superTileSize) to
+		 * the image tile index
+		 * TODO: update/remove if surface grid will be trimmed to fit image
+		 * Currently there are 324 tiles horizontally in the image and 328 in the surfaces
+		 * @param nSurfTile surface tile index in scan order
+		 * @return image tile index in scan order or -1 if outside of the image tiles
+		 */
+		public int getImageTileIndex(
+				int nSurfTile)
+		{
+			// calculate index in tileData (has different dimensions - TODO: trim?
+			int surfaceTilesX = stilesX * superTileSize;
+			int tx = nSurfTile %  surfaceTilesX;
+			int ty = nSurfTile /  surfaceTilesX;
+			if ((tx > imageTilesX) || (ty > imageTilesY)){
+				return -1; // out of image
 			}
+			return imageTilesX * ty + tx;
+		}
+		
 		
 		/**
 		 * Assign tiles to a certain disparity surface if there is only one surface candidate
@@ -1743,15 +1769,15 @@ public class TileSurface {
 		 * @param dbg_Y debug tile Y coordinate
 		 * @return 
 		 */
-		public int [] assignSingleCandidate(
+		public int [] assignTilesToSingleCandidate(
 				final double        maxDiff,
 				final double        minDiffOther,
 				final double        minStrength,
 				final double        maxStrength,
 				final int           moveDirs, // 1 increase disparity, 2 - decrease disparity, 3 - both directions
 				final double        dispNorm, // disparity normalize (proportionally scale down disparity difference if above
-				final int [][]      tileLayers,
-				final TileData [][] tileData,
+//				final int [][]      tileLayers,
+//				final TileData [][] tileData,
 				final double [][][] dispStrength,
                 final int           debugLevel,
 				final int           dbg_X,
@@ -1831,6 +1857,18 @@ public class TileSurface {
 			return stats_new;
 		}
 		
+		public void printStats(int []stats)
+		{
+			System.out.println("printStats(x):"+
+					" NEW_ASSIGNED = " + stats[NEW_ASSIGNED] +
+					" NO_SURF=" + stats[NO_SURF] +
+					" TOO_WEAK=" + stats[TOO_WEAK] +
+					" TOO_STRONG=" + stats[TOO_STRONG] +
+					" TOO_FAR=" + stats[TOO_FAR] +
+					" NOT_UNIQUE=" + stats[NOT_UNIQUE]);
+
+		}
+		
 		
 		/**
 		 * Assign tiles to a certain disparity surface if there is only one surface candidate
@@ -1855,20 +1893,20 @@ public class TileSurface {
 		 * @return 
 		 */
 		
-		public int [] assignSingleCandidate(
+		public int [] assignTilesToSurfaces(
 				final double        maxDiff,
 				final double        minDiffOther, // should be >= maxDiff
 				final double        minStrength,
 				final double        maxStrength,
 				final int           moveDirs, // 1 increase disparity, 2 - decrease disparity, 3 - both directions
 				final boolean       enMulti,
-				final double        surfStrPow, // surface string power
+				final double        surfStrPow, // surface strength power
 				final double        sigma,
 				final double        nSigma,
 				final double        minAdvantage,
 				final double        dispNorm, // disparity normalize (proportionally scale down disparity difference if above
-				final int [][]      tileLayers,
-				final TileData [][] tileData,
+//				final int [][]      tileLayers, // use this.
+//				final TileData [][] tileData,// use this.
 				final double [][][] dispStrength,
                 final int           debugLevel,
 				final int           dbg_X,
@@ -1889,6 +1927,25 @@ public class TileSurface {
 			final AtomicInteger ai = new AtomicInteger(0);
 			final boolean en_lower =  (moveDirs & 1) != 0;
 			final boolean en_higher = (moveDirs & 2) != 0;
+			final double radius = sigma * nSigma;
+			final double rsigma2 = 1.0 / ( 2.0 * sigma * sigma);
+			final int iradius = (int) Math.round(radius + 0.001);  
+			final int field_size = 2 * iradius + 1;
+			final int center_index = iradius * (field_size + 1);
+			final double cost_start = 1.0;
+			final double cost_ortho = 1.0;
+			final double cost_diag  = Math.sqrt(2.0);
+			final int surfTilesX = stilesX * superTileSize;
+			final int [] ldirs8 = {
+					-field_size,
+					-field_size + 1,
+					1,
+					field_size + 1,
+					field_size,
+					field_size - 1,
+					-1,
+					-field_size - 1}; 
+			
 			for (int ml = 0; ml < tileLayers.length; ml++) if (tileLayers[ml] != null){
 				final int fml = ml;
 				ai_numThread.set(0);
@@ -1938,6 +1995,9 @@ public class TileSurface {
 												stats_all[numThread][NOT_UNIQUE] ++;
 											} else { // multi, enabled
 												int [] candidates =       new int     [num_fit_other];
+												// reversed mapping - from layer to candidate (-1 = not a candidate)
+//												int [] rcandidates =      new int     [tileData[nSurfTile].length];
+//												for (int i = 0; i < rcandidates.length; i++) rcandidates[i] = -1;
 												boolean [] close_enough = new boolean [num_fit_other];
 												num_fit_other = 0;
 												for (int ns = 0; ns < tileData[nSurfTile].length; ns++){
@@ -1948,22 +2008,168 @@ public class TileSurface {
 													if (((surf_disp_diff >= 0) && en_higher) || ((surf_disp_diff <= 0) && en_lower)){
 														if (Math.abs(surf_disp_diff) <= minDiffOther){
 															close_enough[num_fit_other] = (Math.abs(surf_disp_diff) <= maxDiff);
+//															rcandidates[ns] = num_fit_other;
 															candidates[num_fit_other++] = ns;
 														}
 													}
 												}
-												double [][] advantages = new double [num_fit_other][num_fit_other];
+												double [][][] distances = new double [num_fit_other][field_size  * field_size ][];
+												// for each local index get surface tile index
+												int [] surfIndices =  new int [field_size  * field_size];
+												int [] imageIndices = new int [field_size  * field_size];
+												int stx0 = (nTile % surfTilesX) - iradius;
+												int sty0 = (nTile / surfTilesX) - iradius;
+												for (int iy = 0; iy < field_size; iy++){
+													for (int ix = 0; ix < field_size; ix++){
+														int indx = iy * field_size + ix;
+														surfIndices[indx] = (sty0 + iy)*  surfTilesX + (stx0 + ix);
+														imageIndices[indx] = getImageTileIndex(surfIndices[indx]);
+													}
+												}
+												// calculate distances from the center point for each surface with wave
+												// algorithm, first regardless of who is closer.
+												// later, when comparing pairs only use the same side
+												for (int isurf = 0; isurf < num_fit_other; isurf++){
+													ArrayList<Point> lwave = new ArrayList<Point>();
+													Point p0 = new Point(center_index, candidates[isurf]);
+													distances[isurf][p0.x] = new double [tileData[nSurfTile].length];
+													distances[isurf][p0.x][p0.y] = cost_start;
+													lwave.add(p0);
+													// run wave build radius (plus 1.0) along each surface connections,
+													// until next radius >= radius
+													while (!lwave.isEmpty()){
+														p0 = lwave.remove(0);
+														int [] neibs = tileData[surfIndices[p0.x]][p0.y].getNeighbors();
+														// try ortho directions first
+														double new_dist =  distances[isurf][p0.x][p0.y] + cost_ortho;
+														if (new_dist <= (radius + cost_start)) {
+															for (int dir = 0; dir < 8; dir +=2) if ( neibs[dir] >= 0){
+																Point pn = new Point (p0.x + ldirs8[dir], neibs[dir]);
+																if (distances[isurf][pn.x] == null){
+																	distances[isurf][pn.x] = new double [tileData[surfIndices[pn.x]].length];
+																}
+																if ((distances[isurf][pn.x][pn.y] == 0) || (distances[isurf][pn.x][pn.y] > new_dist)){
+																	distances[isurf][pn.x][pn.y] = new_dist;
+																	lwave.add(pn);
+																}
+															}
+														}
+														// try diagonal directions second
+														new_dist =  distances[isurf][p0.x][p0.y] + cost_diag;
+														if (new_dist <= (radius + cost_start)) {
+															for (int dir = 1; dir < 8; dir +=2) if ( neibs[dir] >= 0){
+																Point pn = new Point (p0.x + ldirs8[dir], neibs[dir]);
+																if (distances[isurf][pn.x] == null){
+																	distances[isurf][pn.x] = new double [tileData[surfIndices[pn.x]].length];
+																}
+																if ((distances[isurf][pn.x][pn.y] == 0) || (distances[isurf][pn.x][pn.y] > new_dist)){
+																	distances[isurf][pn.x][pn.y] = new_dist;
+																	lwave.add(pn);
+																}
+															}
+														}
+													}
+												}											
+												
+												// pulls belong to pairs, not individual surfaces (difference when they cross)
+												double [][] surface_pulls = new double [num_fit_other][num_fit_other];
 												// now calculate advantage of each one surface (close_enough) to each other (as a ratio)
 												// and then see if it is above minAdvantage
-												for (int ns1 = 0; ns1 < num_fit_other; ns1++){
-													for (int ns2 = ns1 + 1; ns2 < num_fit_other; ns2++){
-														if (close_enough[ns1] || close_enough[ns2]){
-															
+												for (int other_ml = 0; other_ml < tileLayers.length; other_ml++) if (tileLayers[other_ml] != null){
+													for (int lindx = 0; lindx < imageIndices.length; lindx ++){
+														if (imageIndices[lindx] >= 0) {
+															int nsurf = tileLayers_src[other_ml][imageIndices[lindx]] - 1; // assigned surface number (>=0)
+															if ( nsurf >= 0){
+																double strength = dispStrength[other_ml][1][imageIndices[lindx]];
+																// see if this tile belongs to any of the considered surfaces
+																int num_found = 0;
+																boolean [] on_surface = new boolean [num_fit_other];
+																for (int i = 0; i < num_fit_other; i++) {
+																	if ((distances[i][lindx] != null) && (distances[i][lindx][nsurf] > 0.0)){
+																		num_found++;
+																		on_surface[i] = true;
+																	}
+																}
+																if (num_found > 0) { // tile lies on at least one of the considered surfaces
+																	for (int is1 = 0; is1 < num_fit_other; is1++) if (on_surface[is1]) {
+																		// is2 can be any other candidate, just check it is on the same side
+																		// of is1 as in the center (new tile)
+																		for (int is2 = 0; is2 < num_fit_other; is2++) if (is2 != is1) {
+																			boolean good_pair = true;
+																			if (distances[is2][lindx] != null) { // otherwise OK
+																				for (int i = 0; i < distances[is2][lindx].length; i++){
+																					if (distances[is2][lindx][i] >= 0){
+																						if (	((is2 > is1) && (i < nsurf)) ||
+																								((is2 > is1) && (i < nsurf))) {
+																							good_pair = false; // surfaces cross between
+																							break;
+																						}
+																					}
+																				}
+																			}
+																			if (good_pair){
+																				double r = distances[is1][lindx][nsurf] - cost_start;
+																				// pull to is1 when in pair with is2
+																				surface_pulls[is1][is2] = Math.exp(- r * r * rsigma2) * strength ; 
+																			}
+																		}
+																	}
 																	
+																}
+															}
+														}
+													}
+												}
 												
+												double [][] advantages = new double [num_fit_other][num_fit_other];
+												for (int is1 = 0; is1 < num_fit_other; is1++){
+													for (int is2 = is1 + 1; is2 < num_fit_other; is2++){
+														advantages[is1][is2] = surface_pulls[is1][is2]/surface_pulls[is2][is1];
+														advantages[is2][is1] = 1.0/advantages[is1][is2];
+														if (surfStrPow != 0.0){ // consider surface strength also 
+															double str1 = tileData[nSurfTile][candidates[is1]].getStrength();
+															double str2 = tileData[nSurfTile][candidates[is1]].getStrength();
+															if ((str1 > 0.0) && (str2 > 0.0)){
+																advantages[is1][is2] *= Math.pow(str1/str2, surfStrPow);
+																advantages[is2][is1] = 1.0/advantages[is1][is2];
+															} else if (str1 > 0.0) {
+																advantages[is1][is2] = 2.0 * minAdvantage; // sure will win
+																advantages[is2][is1] = (minAdvantage > 0.0) ? (1.0/advantages[is1][is2]) : 0.0; 
+																//minAdvantage
+															} else if (str2 > 0.0) {
+																advantages[is2][is1] = 2.0 * minAdvantage; // sure will win
+																advantages[is1][is2] = (minAdvantage > 0.0) ? (1.0/advantages[is2][is1]) : 0.0; 
+															} else { // both zero - do nothing about surface strengths
+																
+															}
+														}
+													}
+												}
+												// Now see if we have a winner that is (2 could not satisfy, look for the first:
+												// a) close enough, and
+												// b) sufficient advantage over all other candidates
+												fit = -1;
+												for (int is1 = 0; is1 < num_fit_other; is1++){
+													if (close_enough[is1]){ //
+														boolean is_a_winner = true; 
+														for (int is2 = is1 + 1; is2 < num_fit_other; is2++){
+															if (advantages[is1][is2] < minAdvantage){
+																is_a_winner = false;
+																break;
+															}
+														}														
+														if (is_a_winner){
+															fit = is1;
+															break;
 														}
 													}
 												}
+												if (fit >= 0) {
+													tileLayers[fml][nTile] = fit + 1;
+													stats_all[numThread][NEW_ASSIGNED] ++;
+												} else {
+													stats_all[numThread][NOT_UNIQUE] ++;
+												}
 											}
 										}
 									}
@@ -1983,6 +2189,35 @@ public class TileSurface {
 		}
 
 		
+		public int [] InitTilesAssignment(
+				final boolean            force,
+				final double [][][]      dispStrength,
+				final boolean [][]       tileSel,
+                final int                debugLevel)
+		{
+			if (force || (this.tileLayers == null)) {
+				int [][] tileLayers = new int [tileSel.length][];
+				for (int ml = 0; ml < tileSel.length; ml++){
+					if (tileSel[ml] != null){
+						tileLayers[ml] = new int [tileSel[ml].length];
+						for (int i = 0; i < tileSel[ml].length; i++){
+							tileLayers[ml][i] = tileSel[ml][i] ? 0: -1; // 0 - unassigned, -1 - prohibited	
+						}
+					}
+				}
+				this.tileLayers = tileLayers;
+			}
+			int []stats = getTilesAssignStats(tileLayers);
+			if (debugLevel >= -1) {
+				System.out.println("sortTilesToSurfaces(): using "+stats[STAT_NUM_ML] +" measurement layers"+
+						", number of assigned tiles: "+stats[STAT_ASSIGNED]+
+						", number of unassigned tiles: "+stats[STAT_UNASSIGNED]+
+						", number of prohibited tiles: "+stats[STAT_PROHIBITED]+
+						", number of impossible tiles: "+stats[STAT_IMPOSSIBLE]);
+			}
+			return stats;
+		}
+		
 		
 		
 		public int [][] sortTilesToSurfaces(