Fixed MB correction for very small blurs, added thin wire over sky to

sky detection

src/main/java/com/elphel/imagej/gpu/GpuQuad.java

@@ -3822,6 +3822,8 @@ public class GpuQuad{ // quad camera description
 			final boolean []          valid_tiles,            
 			final int                 threadsMax)  // maximal number of threads to launch
 	{
+		final double min_sub = 1e-12;
+		final double min_len = 0.1; // pix
 		int num_pairs = Correlation2d.getNumPairs(num_cams);
 		//change to fixed 511?
 		final int task_code = ((1 << num_pairs)-1) << GPUTileProcessor.TASK_CORR_BITS; //  correlation only
@@ -3875,20 +3877,32 @@ public class GpuQuad{ // quad camera description
 						dx /= mb_len; // unit vector
 						dy /= mb_len;
 						mb_len *= mb_tau; // now in pixels
-						double mb_offs = 1.0; // try 1 pixel. Maybe adjust for non-ortho, e.g. sqrt(2) for diagonal?
-						double min_offs = mb_len_scale * mb_len;
-						if (mb_offs < min_offs) {
-							mb_offs = min_offs;
-						}
-						dx *= mb_offs;
-						dy *= mb_offs;
-						double [] centerXY_sub = {centerXY[0]+dx,centerXY[1]+dy};
-						tp_task_sub.setCenterXY(centerXY_sub);
-						double exp_offs = Math.exp(-mb_offs/mb_len);
-						double gain =    1.0/(1.0 - exp_offs);
-						double gain_sub = -gain * exp_offs;
-						tp_task.setScale(gain);
-						tp_task_sub.setScale(gain_sub);
+						/*
+						double [] centerXY_sub = centerXY;
+						if (mb_len < min_len) {
+							tp_task.setScale(1.0);
+							tp_task_sub.task = 0; // disable
+							tp_task.setScale(-min_sub);
+						} else {
+						*/
+							double mb_offs = 1.0; // try 1 pixel. Maybe adjust for non-ortho, e.g. sqrt(2) for diagonal?
+							double min_offs = mb_len_scale * mb_len;
+							if (mb_offs < min_offs) {
+								mb_offs = min_offs;
+							}
+							dx *= mb_offs;
+							dy *= mb_offs;
+							double [] centerXY_sub = new double[] {centerXY[0]+dx,centerXY[1]+dy};
+							tp_task_sub.setCenterXY(centerXY_sub);
+							double exp_offs = Math.exp(-mb_offs/mb_len);
+							double gain =    1.0/(1.0 - exp_offs);
+							double gain_sub = -gain * exp_offs;
+							if (gain_sub > -min_sub) {
+								gain_sub = -min_sub;
+							}
+							tp_task.setScale(gain);
+							tp_task_sub.setScale(gain_sub);
+//						}
 						boolean bad_margins = false;
 						if (calcPortsCoordinatesAndDerivatives) { // for non-GPU?
 							double [][] disp_dist = new double[quad_main][]; // used to correct 3D correlations (not yet used here)

src/main/java/com/elphel/imagej/tileprocessor/IntersceneMatchParameters.java

@@ -152,6 +152,12 @@ public class IntersceneMatchParameters {
 	public  double  sky_bottom_override =           -300; // maximal average sky value to override lowest_sky_row test
 	public  int     sky_override_shrink =             10; // shrink detected sky before finding hottest tile there
 	
+	// Increase fom if there are enough LMA-defined similar tiles around
+	public  double  disp_boost_min  =                  0.5;
+	public  double  disp_boost_diff  =                 0.35;
+	public  int     disp_boost_neibs  =                2;
+	public  double  disp_boost_amount  =               2.0;
+	
 	
 	// Some "AGC" to adjust how much to discard
 	public  int     margin =                       1;      // do not use tiles if their centers are closer to the image edge
@@ -545,15 +551,18 @@ public class IntersceneMatchParameters {
 		gd.addNumericField("Hottest sky tile to override lowest row",this.sky_bottom_override, 5,7,"",
 				"If the detected sky is all cold enough, bypass lowest row test, allow to raise camera.");
 		gd.addNumericField("Shrink before finding hottest sky",      this.sky_override_shrink, 0,3,"",
-				"Shrink detected sky before looking for the hottest skyt tile (blurred skyline in wet atmosphere).");
-		
+				"Shrink detected sky before looking for the hottest sky tile (blurred skyline in wet atmosphere).");
+		
+		gd.addMessage  ("Boost FOM if there are enough LMA-defined good neighbors around (Thin wires over the sky)");
+		gd.addNumericField("Minimal disparity",                      this.disp_boost_min, 5,7,"pix",
+				"Minimal LMA-defined disparity to boost FOM.");
+		gd.addNumericField("Disparity difference",                   this.disp_boost_diff, 5,7,"pix",
+				"Maximal disparity difference to neighbor to count.");
+		gd.addNumericField("Number of good neighbors",               this.disp_boost_neibs, 0,3,"",
+				"Number of neighbors (of 8) to have small disparity difference to boost FOM.");
+		gd.addNumericField("Boost amount",                           this.disp_boost_amount, 5,7,"x",
+				"Multiply FOM by this value if number of neighbors is exactly minimal. Sacale proportional to the total number of neighbors.");
 		gd.addTab         ("Inter-Match", "Parameters for full-resolution scene matching");
-//		gd.addTab("Interscene Equalization","Equalization of the interscene correlation confidence to improve camera X,Y,Z matching");
-
-
-		
-		
-		
 		
 		
 		gd.addMessage  ("Interscene match parameters");
@@ -931,6 +940,11 @@ public class IntersceneMatchParameters {
 		this.sky_bottom_override =      gd.getNextNumber();    
 		this.sky_override_shrink =(int) gd.getNextNumber();
 		
+		this.disp_boost_min =           gd.getNextNumber();    
+		this.disp_boost_diff =          gd.getNextNumber();    
+		this.disp_boost_neibs =   (int) gd.getNextNumber();
+		this.disp_boost_amount =        gd.getNextNumber();    
+		
 		this.margin =             (int) gd.getNextNumber();
 		this.sensor_mask_inter=   (int) gd.getNextNumber();
 		this.use_partial =              gd.getNextBoolean();
@@ -1206,6 +1220,11 @@ public class IntersceneMatchParameters {
 		properties.setProperty(prefix+"sky_bottom_override",  this.sky_bottom_override+""); // double
 		properties.setProperty(prefix+"sky_override_shrink",  this.sky_override_shrink+""); // int
 		
+		properties.setProperty(prefix+"disp_boost_min",       this.disp_boost_min+"");      // double
+		properties.setProperty(prefix+"disp_boost_diff",      this.disp_boost_diff+"");     // int
+		properties.setProperty(prefix+"disp_boost_neibs",     this.disp_boost_neibs+"");    // double
+		properties.setProperty(prefix+"disp_boost_amount",    this.disp_boost_amount+"");   // double
+		
 		properties.setProperty(prefix+"margin",               this.margin+"");              // int
 		properties.setProperty(prefix+"sensor_mask_inter",    this.sensor_mask_inter+"");   // int
 		properties.setProperty(prefix+"use_partial",          this.use_partial+"");         // boolean
@@ -1433,6 +1452,11 @@ public class IntersceneMatchParameters {
 		if (properties.getProperty(prefix+"lowest_sky_row")!=null)       this.lowest_sky_row=Integer.parseInt(properties.getProperty(prefix+"lowest_sky_row"));
 		if (properties.getProperty(prefix+"sky_bottom_override")!=null)  this.sky_bottom_override=Double.parseDouble(properties.getProperty(prefix+"sky_bottom_override"));
 		if (properties.getProperty(prefix+"sky_override_shrink")!=null)  this.sky_override_shrink=Integer.parseInt(properties.getProperty(prefix+"sky_override_shrink"));
+
+		if (properties.getProperty(prefix+"disp_boost_min")!=null)       this.disp_boost_min=Double.parseDouble(properties.getProperty(prefix+"disp_boost_min"));
+		if (properties.getProperty(prefix+"disp_boost_diff")!=null)      this.disp_boost_diff=Double.parseDouble(properties.getProperty(prefix+"disp_boost_diff"));
+		if (properties.getProperty(prefix+"disp_boost_neibs")!=null)     this.disp_boost_neibs=Integer.parseInt(properties.getProperty(prefix+"disp_boost_neibs"));
+		if (properties.getProperty(prefix+"disp_boost_amount")!=null)    this.disp_boost_amount=Double.parseDouble(properties.getProperty(prefix+"disp_boost_amount"));
 		
 		if (properties.getProperty(prefix+"margin")!=null)               this.margin=Integer.parseInt(properties.getProperty(prefix+"margin"));
 		if (properties.getProperty(prefix+"sensor_mask_inter")!=null)    this.sensor_mask_inter=Integer.parseInt(properties.getProperty(prefix+"sensor_mask_inter"));
@@ -1678,6 +1702,11 @@ public class IntersceneMatchParameters {
 		imp.lowest_sky_row =        this.lowest_sky_row;
 		imp.sky_bottom_override =   this.sky_bottom_override;
 		imp.sky_override_shrink =   this.sky_override_shrink;
+
+		imp.disp_boost_min =        this.disp_boost_min;
+		imp.disp_boost_diff =       this.disp_boost_diff;
+		imp.disp_boost_neibs =      this.disp_boost_neibs;
+		imp.disp_boost_amount =     this.disp_boost_amount;
 		
 		imp.margin                = this.margin;
 		imp.sensor_mask_inter     = this.sensor_mask_inter;

src/main/java/com/elphel/imagej/tileprocessor/OpticalFlow.java

@@ -2802,7 +2802,7 @@ public class OpticalFlow {
 		for (int ithread = 0; ithread < threads.length; ithread++) {
 			threads[ithread] = new Thread() {
 				public void run() {
-					for (int nTile = ai.getAndIncrement(); nTile < tiles; nTile = ai.getAndIncrement()) if (!Double.isNaN(disparity_ref[nTile])) {
+					for (int nTile = ai.getAndIncrement(); nTile < tiles; nTile = ai.getAndIncrement()) if (!Double.isNaN(disparity_ref[nTile])) {// null !!!
 						double disparity = disparity_ref[nTile];
 						int tileY = nTile / tilesX;  
 						int tileX = nTile % tilesX;
@@ -4139,6 +4139,11 @@ public class OpticalFlow {
  		int    lowest_sky_row =      clt_parameters.imp.lowest_sky_row;
  		double sky_bottom_override = clt_parameters.imp.sky_bottom_override;
  		int    sky_override_shrink = clt_parameters.imp.sky_override_shrink;
+		double disp_boost_min =      clt_parameters.imp.disp_boost_min; //  0.5;
+		double disp_boost_diff  =    clt_parameters.imp.disp_boost_diff; // 0.35;
+		int    disp_boost_neibs  =   clt_parameters.imp.disp_boost_neibs; // 2;
+		double disp_boost_amount  =  clt_parameters.imp.disp_boost_amount; // 2.0;
+
 		
 		boolean [] ref_blue_sky = null; // turn off "lma" in the ML output  
 		
@@ -4278,7 +4283,11 @@ public class OpticalFlow {
 							lowest_sky_row,        //  =   50;// appears that low - invalid, remove completely
 							sky_bottom_override,   // double sky_temp_override,     // really cold average seed - ignore lowest_sky_row filter
 							sky_override_shrink,   // int    shrink_for_temp,       // shrink before finding hottest sky
-							sky_highest_min,       //  =   100; // lowest absolute value should not be higher (requires photometric) 
+							sky_highest_min,       //  =   100; // lowest absolute value should not be higher (requires photometric)
+							disp_boost_min,        // double disp_boost_min,    //  = 0.5;
+							disp_boost_diff,       //double disp_boost_diff,   //  = 0.35;
+							disp_boost_neibs,   //int    disp_boost_neibs,  //  = 2;
+							disp_boost_amount,  //double disp_boost_amount, //  = 2.0;
 							dsi[TwoQuadCLT.DSI_STRENGTH_AUX], // double [] strength,
 							dsi[TwoQuadCLT.DSI_SPREAD_AUX], // double [] spread,
 							dsi[TwoQuadCLT.DSI_DISPARITY_AUX_LMA], // double [] spread,
@@ -4332,7 +4341,11 @@ public class OpticalFlow {
 							lowest_sky_row,        //  =     50;// appears that low - invalid, remove completely
 							sky_bottom_override,   // double sky_temp_override,     // really cold average seed - ignore lowest_sky_row filter
 							sky_override_shrink,   // int    shrink_for_temp,       // shrink before finding hottest sky
-							sky_highest_min,       //  =   100; // lowest absolute value should not be higher (requires photometric) 
+							sky_highest_min,       //  =   100; // lowest absolute value should not be higher (requires photometric)
+							disp_boost_min,        // double disp_boost_min,    //  = 0.5;
+							disp_boost_diff,       //double disp_boost_diff,   //  = 0.35;
+							disp_boost_neibs,      //int    disp_boost_neibs,  //  = 2;
+							disp_boost_amount,     //double disp_boost_amount, //  = 2.0;
 							dsi[TwoQuadCLT.DSI_STRENGTH_AUX],      // double [] strength,
 							dsi[TwoQuadCLT.DSI_SPREAD_AUX],        // double [] spread,
 							dsi[TwoQuadCLT.DSI_DISPARITY_AUX_LMA], //double [] disp_lma,
@@ -4456,6 +4469,7 @@ public class OpticalFlow {
 							last_diff);
 				}
 				// Refine with LMA
+				
 				scenes_xyzatr[scene_index] = adjustPairsLMAInterscene(
 						clt_parameters,      // CLTParameters  clt_parameters,
 						quadCLTs[ref_index], // QuadCLT reference_QuadCLT,
@@ -4595,6 +4609,7 @@ public class OpticalFlow {
 						clt_parameters,     // CLTParameters  clt_parameters,
 						reliable_ref,       // boolean []     reliable_ref, // null or bitmask of reliable reference tiles
 						quadCLTs,           // QuadCLT []     quadCLTs,
+						!batch_mode, // boolean        test_motion_blur,
 						debugLevel) ;       // int            debugLevel)
 				// should update earliest_scene
 				if ((ref_index - earliest_scene + 1) < min_num_scenes) {
@@ -14066,9 +14081,10 @@ public double[][] correlateIntersceneDebug( // only uses GPU and quad
 			CLTParameters  clt_parameters,
 			boolean []     reliable_ref, // null or bitmask of reliable reference tiles			
     		QuadCLT []     quadCLTs,
+    		boolean        test_motion_blur,
 			int            debugLevel)
 	{  
-		boolean test_motion_blur = true;
+//		boolean test_motion_blur = true;//false
 		
 		boolean mb_en =       clt_parameters.imp.mb_en;
 		double  mb_tau =      clt_parameters.imp.mb_tau;      // 0.008; // time constant, sec
@@ -14831,6 +14847,11 @@ public double[][] correlateIntersceneDebug( // only uses GPU and quad
 			double         max_rms,
 			int            debug_level)
 	{
+	    boolean use_lma_dsi =          clt_parameters.imp.use_lma_dsi;
+		if (ref_disparity == null) { 
+			ref_disparity = reference_QuadClt.getDLS()[use_lma_dsi?1:0];
+		}
+
 		// TODO: set reference as new version of adjustPairsLMAInterscene() assumes it set
 		int        margin = clt_parameters.imp.margin;
 		double [][] pXpYD_ref = transformToScenePxPyD( // full size - [tilesX*tilesY], some nulls

src/main/java/com/elphel/imagej/tileprocessor/QuadCLTCPU.java

@@ -457,12 +457,19 @@ public class QuadCLTCPU {
 			int    lowest_sky_row,        //  =     50;// appears that low - invalid, remove completely
 			double sky_temp_override,     // really cold average seed - ignore lowest_sky_row filter
 			int shrink_for_temp,          // shrink before finding hottest sky
-			double sky_highest_min,       //  =   -50; // 100; // lowest absolute value should not be higher (requires photometric) 
+			double sky_highest_min,       //  =   -50; // 100; // lowest absolute value should not be higher (requires photometric)
+			
+			double disp_boost_min,    //  = 0.5;
+			double disp_boost_diff,   //  = 0.35;
+			int    disp_boost_neibs,  //  = 2;
+			double disp_boost_amount, //  = 2.0;
+			
 			int    width,
 			double [] strength,
 			double [] spread,
 			double [] disparity,
 			double [] avg_val,
+			
 			int debugLevel) { // >0 to show
 		if ((strength == null) || (spread==null)) {
 			return null;
@@ -556,11 +563,36 @@ public class QuadCLTCPU {
 		String [] dbg_in_titles =    {"fom", "strength", "spread", "disparity", "avg_val", "tscale"};
 		String [] dbg_titles = {"sky", "seed", "max", "shrank","full_shrank","neck_shrank","reexpand"};
 		
-		if (debugLevel>0) {
-			double [] fom = new double[strength.length];
-			for (int i = 0; i < fom.length; i++) {
-				fom[i] =  Math.max(strength[i], min_strength) * spread[i] * temp_scales[i];
+		double [][] dbg_img = (debugLevel>0) ? new double [dbg_titles.length][strength.length]:null;
+		TileNeibs tn = new TileNeibs(width,strength.length/width);
+
+		boolean [] sky_tiles =      new  boolean [strength.length];
+		boolean [] prohibit_tiles = new  boolean [strength.length];
+		double [] fom = new double[strength.length];
+		for (int i = 0; i < sky_tiles.length; i++) {
+			fom[i] = Math.max(strength[i], min_strength) * spread[i];
+			if (temp_scales != null) {
+				fom[i] *=temp_scales[i];
 			}
+		}		
+		if (disp_boost_amount > 0) {
+			for (int nTile = 0; nTile < sky_tiles.length; nTile++) if (disparity[nTile] >= disp_boost_min){
+				int num_neibs = 0;
+				for (int dir = 0; dir < 8; dir++) {
+					int nTile1 = tn.getNeibIndex(nTile, dir);
+					if ((nTile1 >= 0) && (Math.abs(disparity[nTile] - disparity[nTile1]) <=disp_boost_diff)) {
+						num_neibs++;
+					}
+				}
+				if (num_neibs >= disp_boost_neibs ) {
+					fom[nTile] *= disp_boost_amount * num_neibs / disp_boost_neibs;
+				}
+			}			
+		}
+		
+		
+		
+		if (debugLevel>0) {
 			(new ShowDoubleFloatArrays()).showArrays(
 					new double[][] {fom, strength, spread, disparity, avg_val, temp_scales},
 					width,
@@ -570,20 +602,14 @@ public class QuadCLTCPU {
 					dbg_in_titles); //	dsrbg_titles);
 
 		}
-		double [][] dbg_img = (debugLevel>0) ? new double [dbg_titles.length][strength.length]:null;
 		
-		boolean [] sky_tiles =      new  boolean [strength.length];
-		boolean [] prohibit_tiles = new  boolean [strength.length];
 		for (int i = 0; i < sky_tiles.length; i++) {
-			double d = Math.max(strength[i], min_strength) * spread[i];
-			if (temp_scales != null) {
-				d *=temp_scales[i];
-			}
-			prohibit_tiles[i] = (d >= sky_lim);
-			int row = i/width;
-//			sky_tiles[i] =    (row < seed_rows) &&  (d < sky_seed) && !(disparity[i] > disparity_seed);
-			sky_tiles[i] = (d < sky_seed) && !(disparity[i] > disparity_seed);
+			prohibit_tiles[i] = (fom[i] >= sky_lim);
+			sky_tiles[i] = (fom[i] < sky_seed) && !(disparity[i] > disparity_seed);
 		}
+		
+		
+		
 		//seed_rows
 		if (dbg_img != null) {
 			for (int i = 0; i < sky_tiles.length; i++) {
@@ -591,7 +617,6 @@ public class QuadCLTCPU {
 				dbg_img[2][i] = prohibit_tiles[i]? 0 : 1;
 			}			
 		}
-		TileNeibs tn = new TileNeibs(width,sky_tiles.length/width);
 		tn.shrinkSelection(
 				sky_shrink,    // int        shrink,           // grow tile selection by 1 over non-background tiles 1: 4 directions, 2 - 8 directions, 3 - 8 by 1, 4 by 1 more
 				sky_tiles,     // boolean [] tiles,
@@ -723,7 +748,11 @@ public class QuadCLTCPU {
 			int    lowest_sky_row,        //  =     50;// appears that low - invalid, remove completely
 			double sky_temp_override,     // really cold average seed - ignore lowest_sky_row filter
 			int    shrink_for_temp,       // shrink before finding hottest sky
-			double sky_highest_max,       //  =   100; // lowest absolute value should not be higher (requires photometric) 
+			double sky_highest_max,       //  =   100; // lowest absolute value should not be higher (requires photometric)
+			double disp_boost_min,    //  = 0.5;
+			double disp_boost_diff,   //  = 0.35;
+			int    disp_boost_neibs,  //  = 2;
+			double disp_boost_amount, //  = 2.0;
 			double [] strength,
 			double [] spread,
 			double [] disp_lma,
@@ -746,6 +775,10 @@ public class QuadCLTCPU {
 				sky_temp_override,     // double sky_temp_override,     // really cold average seed - ignore lowest_sky_row filter
 				shrink_for_temp,       // int shrink_for_temp,          // shrink before finding hottest sky
 				sky_highest_max,       //  =   100; // lowest absolute value should not be higher (requires photometric) 
+				disp_boost_min,    //  = 0.5;
+				disp_boost_diff,   //  = 0.35;
+				disp_boost_neibs,  //  = 2;
+				disp_boost_amount, //  = 2.0;
 				width,
 				strength,
 				spread,