Improved filtering

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

@@ -1148,6 +1148,9 @@ public class CLTPass3d{
 			}
 		}
 
+        public void setSecondMax(double [] second_max) { // setting "fake" for the non-measurement scan
+        	this.second_max = second_max;
+        }
         public void setSecondMax() {
         	this.second_max = getSecondMaxDiff(false);
         }

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

@@ -3982,16 +3982,8 @@ public class OpticalFlow {
 		            debug_level+1);
 		}		
 	}
-	/*
-	public class VideoPathParams{
-		public String path;
-		public int    width;
-		public int    height;
-		public int    gap;
-		public double fps;
-		public double real_fps;
-	}
-	*/
+
+	
 	/**
 	 * Build series of poses from just a single (reference) scene.
 	 * @param quadCLT_main
@@ -4142,6 +4134,12 @@ public class OpticalFlow {
 		
 		double  min_ref_str =        clt_parameters.imp.min_ref_str;
 		
+		double sky_seed =       7.0;  // start with product of strength by diff_second below this
+		double sky_lim  =      15.0; // then expand to product of strength by diff_second below this
+		int    sky_shrink =       4;
+		int    sky_expand_extra = 0; // 1?
+		
+		boolean [] ref_blue_sky = null; // turn off "lma" in the ML output  
 		
 		if (reuse_video) { // disable all other options
 			generate_mapped = false;
@@ -4200,6 +4198,7 @@ public class OpticalFlow {
 		}
 		
 		// 1. Reference scene DSI
+		while (quadCLTs[ref_index].getBlueSky() == null) {
 			if (build_ref_dsi) {
 				TwoQuadCLT.copyJP4src( // actually there is no sense to process multiple image sets. Combine with other
 						// processing?
@@ -4232,6 +4231,7 @@ public class OpticalFlow {
 				dsi[TwoQuadCLT.DSI_DISPARITY_AUX] =     aux_last_scan[0];
 				dsi[TwoQuadCLT.DSI_STRENGTH_AUX] =      aux_last_scan[1];
 				dsi[TwoQuadCLT.DSI_DISPARITY_AUX_LMA] = aux_last_scan[2];
+				dsi[TwoQuadCLT.DSI_SPREAD_AUX] =        aux_last_scan[3];
 
 				if (quadCLT_main.correctionsParameters.clt_batch_dsi_cm_strength) {
 					CLTPass3d scan = new CLTPass3d(quadCLTs[ref_index].tp);
@@ -4248,7 +4248,7 @@ public class OpticalFlow {
 							-1.0,           // final double        mismatch_override, // keep tile with large mismatch if there is LMA with really strong correlation
 							threadsMax,     // final int         threadsMax,  // maximal number of threads to launch
 							updateStatus,   // final boolean     updateStatus,
-						debugLevel);    // final int         debugLevel);
+							debugLevel-2);  // final int         debugLevel);
 					dsi[TwoQuadCLT.DSI_STRENGTH_AUX] =      scan.getStrength();
 					if (debugLevel > 1) {
 						quadCLTs[ref_index].tp.showScan(
@@ -4258,6 +4258,16 @@ public class OpticalFlow {
 				} else {
 					dsi[TwoQuadCLT.DSI_STRENGTH_AUX] =      aux_last_scan[1];
 				}
+				int tilesX = quadCLTs[ref_index].getTileProcessor().getTilesX();
+				quadCLTs[ref_index].setBlueSky  (
+						sky_seed,           // double sky_seed, //  =       7.0;  // start with product of strength by diff_second below this
+						sky_lim,            // double sky_lim, //   =      15.0; // then expand to product of strength by diff_second below this
+						sky_shrink,         // int    sky_shrink, //  =       4;
+						sky_expand_extra,   // int    sky_expand_extra, //  = 100; // 1?
+						dsi[TwoQuadCLT.DSI_STRENGTH_AUX], // double [] strength,
+						dsi[TwoQuadCLT.DSI_SPREAD_AUX], // double [] spread,
+						debugLevel);        // int debugLevel)
+
 				quadCLTs[ref_index].saveDSIAll (
 						"-DSI_MAIN", // String suffix, // "-DSI_MAIN"
 						dsi);
@@ -4281,8 +4291,24 @@ public class OpticalFlow {
 				quadCLTs[ref_index].set_orient(0); // reset orientations
 				quadCLTs[ref_index].set_accum(0);  // reset accumulations ("build_interscene") number
 				earliest_scene = 0;  // reset failures, try to use again all scenes
-			
-		} // if (build_ref_dsi) {
+			} else {// if (build_ref_dsi) {
+				// read DSI_MAIN
+				double [][] dsi = quadCLTs[ref_index].readDsiMain();
+				if (dsi[TwoQuadCLT.DSI_SPREAD_AUX] == null) {
+					System.out.println("DSI_MAIN file has old format and does not have spread data, will recalculate.");
+				} else {
+					quadCLTs[ref_index].setBlueSky  (
+							sky_seed,           // double sky_seed, //  =       7.0;  // start with product of strength by diff_second below this
+							sky_lim,            // double sky_lim, //   =      15.0; // then expand to product of strength by diff_second below this
+							sky_shrink,         // int    sky_shrink, //  =       4;
+							sky_expand_extra,   // int    sky_expand_extra, //  = 100; // 1?
+							dsi[TwoQuadCLT.DSI_STRENGTH_AUX], // double [] strength,
+							dsi[TwoQuadCLT.DSI_SPREAD_AUX], // double [] spread,
+							debugLevel);        // int debugLevel)
+				}
+			}
+		} // while (blue_sky == null) 
+		ref_blue_sky = quadCLTs[ref_index].getBlueSky();
 		
 		quadCLTs[ref_index] = (QuadCLT) quadCLT_main.spawnQuadCLT( // restores dsi from "DSI-MAIN"
 				set_channels[ref_index].set_name,
@@ -4290,6 +4316,7 @@ public class OpticalFlow {
 				colorProcParameters, //
 				threadsMax,
 				debugLevel);
+		quadCLTs[ref_index].setBlueSky(ref_blue_sky);
 		
 		quadCLTs[ref_index].setDSRBG(
 				clt_parameters, // CLTParameters  clt_parameters,
@@ -4566,6 +4593,7 @@ public class OpticalFlow {
 						colorProcParameters, //
 						threadsMax,
 						debugLevel);
+				quadCLTs[ref_index].setBlueSky(ref_blue_sky);
 				quadCLTs[ref_index].setDSRBG(
 						clt_parameters, // CLTParameters  clt_parameters,
 						threadsMax,     // int            threadsMax,  // maximal number of threads to launch
@@ -5347,6 +5375,8 @@ public class OpticalFlow {
 							combo_dsn_final,     // double [][]    combo_dsn_final, // dls,
 							quadCLTs[ref_index], // QuadCLT        scene,
 							debugLevel); // int            debugLevel);// > 0
+			
+			
 			intersceneMlExport(
 					clt_parameters,      // CLTParameters        clt_parameters,
 					ers_reference,       // ErsCorrection        ers_reference,
@@ -8132,15 +8162,22 @@ public class OpticalFlow {
 		
 
 		double  disp_ampl = Math.max(Math.abs(disparity_low),Math.abs(disparity_high));
-		//		String suffix =ref_scene.correctionsParameters.mlDirectory; // now "ML32
 		int      indx_ref = scenes.length - 1; // Always added to the end even if out-of order
 		QuadCLT  ref_scene = scenes[indx_ref]; // ordered by increasing timestamps
-//		int num_scenes =  scenes.length; // FIXME: skip unmatched
 		final int num_scenes =  indx_ref - ref_scene.getEarliestScene(scenes) + 1;
 		final int tilesX = ref_scene.getTileProcessor().getTilesX();
 		final int tilesY = ref_scene.getTileProcessor().getTilesY();
 		final int tiles = tilesX * tilesY;
 
+		boolean [] blue_sky = ref_scene.getBlueSky();
+		double  [] payload_blue_sky = new double[tiles];
+		Arrays.fill(payload_blue_sky,Double.NaN);
+		if (blue_sky != null) {
+			for (int i = 0; i < tiles; i++) {
+				payload_blue_sky[i] = blue_sky[i] ? 1.0 : 0.0;
+			}
+		}
+		
 		double  fat_zero_single = clt_parameters.getGpuFatZero(ref_scene.isMonochrome()); // for single scene
 		ImageDtt image_dtt;
 		image_dtt = new ImageDtt(
@@ -8214,7 +8251,9 @@ public class OpticalFlow {
 					combo_dsn_final[COMBO_DSN_INDX_LMA_BG],     // masked copy from BG disparity 
 					combo_dsn_final[COMBO_DSN_INDX_CHANGE_BG],  // increment, BG
 					combo_dsn_final[COMBO_DSN_INDX_DISP_FG],    //cumulative disparity (from CM or POLY), FG == COMBO_DSN_INDX_DISP
-					combo_dsn_final[COMBO_DSN_INDX_DISP_BG_ALL] // cumulative BG disparity (Use FG where no BG is available)
+					combo_dsn_final[COMBO_DSN_INDX_DISP_BG_ALL],// cumulative BG disparity (Use FG where no BG is available)
+					payload_blue_sky                            // detected featureless sky - 1.0, reliable - 0.0, no data - NaN
+					
 			};
 			for (int i = 0; i < payload.length; i++) {
 				add_tile_meta(
@@ -8275,6 +8314,7 @@ public class OpticalFlow {
 			imp_ml.setProperty("metaBGLastDiff",      ""+9);
 			imp_ml.setProperty("metaFGDisparity",     ""+10); //=="metaGTDisparity"
 			imp_ml.setProperty("metaBGDisparityAll",  ""+11);
+			imp_ml.setProperty("metaBlueSky",         ""+12);
 			(new JP46_Reader_camera(false)).encodeProperiesToInfo(imp_ml);			
 			FileSaver fs=new FileSaver(imp_ml);
 			fs.saveAsTiff(file_path);
@@ -8340,8 +8380,9 @@ public class OpticalFlow {
 					combo_dsn_final[COMBO_DSN_INDX_STRENGTH_BG],// background strength
 					combo_dsn_final[COMBO_DSN_INDX_LMA_BG],     // masked copy from BG disparity 
 					combo_dsn_final[COMBO_DSN_INDX_CHANGE_BG],  // increment, BG
-					combo_dsn_final[COMBO_DSN_INDX_DISP_FG],    //cumulative disparity (from CM or POLY), FG == COMBO_DSN_INDX_DISP
-					combo_dsn_final[COMBO_DSN_INDX_DISP_BG_ALL] // cumulative BG disparity (Use FG where no BG is available)
+					combo_dsn_final[COMBO_DSN_INDX_DISP_FG],    // cumulative disparity (from CM or POLY), FG == COMBO_DSN_INDX_DISP
+					combo_dsn_final[COMBO_DSN_INDX_DISP_BG_ALL],// cumulative BG disparity (Use FG where no BG is available)
+					payload_blue_sky                            // detected featureless sky - 1.0, reliable - 0.0, no data - NaN
 			};
 			for (int i = 0; i < payload.length; i++) {
 				add_tile_meta(
@@ -8418,6 +8459,7 @@ public class OpticalFlow {
 			imp_ml.setProperty("metaBGLastDiff",      ""+9);
 			imp_ml.setProperty("metaFGDisparity",     ""+10); //=="metaGTDisparity"
 			imp_ml.setProperty("metaBGDisparityAll",  ""+11);
+			imp_ml.setProperty("metaBlueSky",         ""+12);
 			(new JP46_Reader_camera(false)).encodeProperiesToInfo(imp_ml);			
 			FileSaver fs=new FileSaver(imp_ml);
 			fs.saveAsTiff(file_path);
@@ -11590,7 +11632,6 @@ public double[][] correlateIntersceneDebug( // only uses GPU and quad
 			final int         weak_min_neibs,
 			final double      strong_strength,
 			final double      weak_strength)
-			
 	{
 		
 		final int tilesX = scene.getTileProcessor().getTilesX();
@@ -11598,23 +11639,34 @@ public double[][] correlateIntersceneDebug( // only uses GPU and quad
 		final int transform_size = scene.getTileProcessor().getTileSize();
 		final int tiles =tilesX * tilesY;
 
-		String [] dbg_titles = {"str", "lma", "clean-lma", "disp","-lma","by-lma","-nonlma", "few_weak", "old-disp","old-sngl","weak","filled"}; 
+		String [] dbg_titles = {"str", "lma", "clean-lma", "disp","-sky","-lma","by-lma","-nonlma", "few_weak", "old-disp","old-sngl","weak","filled"}; 
 		double [][] dbg_img = new double [dbg_titles.length][];
 		double [] clean_lma = dls[1].clone();
-		for (int i = 0; i <clean_lma.length; i++) {
+		double [] clean_disparity = dls[0].clone();
+		for (int i = 00; i <clean_lma.length; i++) {
 			if (dls[2][i] < min_strength_lma) {
 				clean_lma[i] = Double.NaN;
 			}
 		}
+		boolean [] blue_sky = scene.getBlueSky();
+		if (blue_sky != null) {
+			for (int i = 0; i < clean_lma.length; i++) if (blue_sky[i]){
+				clean_lma[i] = Double.NaN;
+				clean_disparity[i] = 0.0;
+			}
+		}
+		
 		//Remove crazy LMA high-disparity tiles
 		dbg_img[0] = dls[2].clone();
 		dbg_img[1] = dls[1].clone();
 		dbg_img[2] = clean_lma.clone();
 		dbg_img[3] = dls[0].clone();
+		dbg_img[4] = clean_disparity.clone();
 		double [] disp_outliers = QuadCLT.removeDisparityLMAOutliers( // nothing removed (trying to remove bad LMA)
 				false,                       // final boolean     non_ma,
 //				dls, //final double [][] dls,
-				new double[][] {dls[1], dls[1], clean_lma}, //final double [][] dls,
+//				new double[][] {dls[0],  clean_lma, dls[2]}, //final double [][] dls,
+				new double[][] {clean_disparity,  clean_lma, dls[2]}, //final double [][] dls,
 				outliers_lma_max_strength,   // final double      max_strength,  // do not touch stronger
 				outliers_lma_nth_fromextrem, // final int         nth_fromextrem, // 0 - compare to max/min. 1 - second max/min, ... 
 				outliers_tolerance_absolute, // final double      tolerance_absolute,
@@ -11623,9 +11675,9 @@ public double[][] correlateIntersceneDebug( // only uses GPU and quad
 				outliers_max_iter,           // final int         max_iter,
 				threadsMax,                  // final int         threadsMax,
 				debug_level);                // final int         debug_level)
-		dbg_img[4] = disp_outliers.clone();
+		dbg_img[5] = disp_outliers.clone();
 		disp_outliers = QuadCLT.removeDisparityOutliersByLMA( // removed sky, keeps sky edge near strong objects
-				new double[][] {disp_outliers, dls[1], clean_lma}, //final double [][] dls,
+				new double[][] {disp_outliers, clean_lma, dls[2]}, //final double [][] dls,
 				outliers_from_lma_max_strength,  // final double      max_strength,  // do not touch stronger
 				diff_from_lma_pos,           // final double      diff_from_lma_pos,   // Difference from farthest FG objects (OK to have large, e.g. 100)
 				diff_from_lma_neg,           // final double      diff_from_lma_neg,   // Difference from nearest BG objects (small, as FG are usually more visible)
@@ -11634,11 +11686,11 @@ public double[][] correlateIntersceneDebug( // only uses GPU and quad
 				tilesX,                      // final int         width,               //tilesX
 				threadsMax,                  // final int         threadsMax,
 				debug_level);                // final int         debug_level)
-		dbg_img[5] = disp_outliers.clone();
+		dbg_img[6] = disp_outliers.clone();
 		// mostly filter infinity, clouds, sky
 		disp_outliers = QuadCLT.removeDisparityLMAOutliers( // filter non-lma tiles // removed too few !!!
 				true,                       // final boolean     non_ma,
-				new double[][] {disp_outliers, dls[1], clean_lma}, //final double [][] dls,
+				new double[][] {disp_outliers, clean_lma, dls[2]}, //final double [][] dls,
 				outliers_max_strength,   // final double      max_strength,  // do not touch stronger
 				outliers_nth_fromextrem, // final int         nth_fromextrem, // 0 - compare to max/min. 1 - second max/min, ... 
 				outliers_tolerance_absolute, // final double      tolerance_absolute,
@@ -11647,10 +11699,10 @@ public double[][] correlateIntersceneDebug( // only uses GPU and quad
 				outliers_max_iter,           // final int         max_iter,
 				threadsMax,                  // final int         threadsMax,
 				debug_level);                // final int         debug_level)
-		dbg_img[6] = disp_outliers.clone();
+		dbg_img[7] = disp_outliers.clone();
 		
 		disp_outliers = QuadCLT.removeFewWeak( // filter non-lma tiles // removed too few !!!
-				new double[][] {disp_outliers, dls[1], clean_lma}, //final double [][] dls,
+				new double[][] {disp_outliers, clean_lma, dls[2]}, //final double [][] dls,
 				strong_strength,             // final double      strong,
 				weak_strength,               // final double      weak,
 				weak_min_neibs,              // final int         min_neibs, 
@@ -11661,10 +11713,10 @@ public double[][] correlateIntersceneDebug( // only uses GPU and quad
 				threadsMax,                  // final int         threadsMax,
 				debug_level);                // final int         debug_level)
 		
-		dbg_img[7] = disp_outliers.clone();
+		dbg_img[8] = disp_outliers.clone();
 		// Pre- 2022 filters, some may be obsolete 
 		disp_outliers = QuadCLT.removeDisparityOutliers(
-				new double[][] {disp_outliers, clean_lma}, //final double [][] dls,
+				new double[][] {disp_outliers, dls[2]}, //final double [][] dls,
 				outliers_max_strength,       // final double      max_strength,  // do not touch stronger
 				outliers_nth_fromextrem,     // final int         nth_fromextrem, // 0 - compare to max/min. 1 - second max/min, ... 
 				outliers_tolerance_absolute, // final double      tolerance_absolute,
@@ -11674,10 +11726,10 @@ public double[][] correlateIntersceneDebug( // only uses GPU and quad
 				false,                       // final boolean     fit_completely, // do not add tolerance when replacing
 				threadsMax,                  // final int         threadsMax,
 				debug_level);                // final int         debug_level)
-		dbg_img[8] = disp_outliers.clone();
+		dbg_img[9] = disp_outliers.clone();
 		// remove extreme single-tile outliers (some may be strong - 0.404)
 		disp_outliers = QuadCLT.removeDisparityOutliers(
-				new double[][] {disp_outliers, clean_lma}, //final double [][] dls,
+				new double[][] {disp_outliers,  dls[2]}, //final double [][] dls,
 				outliers_max_strength2,       // final double      max_strength,  // do not touch stronger
 				outliers_nth_fromextrem2,     // final int         nth_fromextrem, // 0 - compare to max/min. 1 - second max/min, ... 
 				outliers_tolerance_absolute2, // final double      tolerance_absolute,
@@ -11687,16 +11739,23 @@ public double[][] correlateIntersceneDebug( // only uses GPU and quad
 				false,                       //  final boolean     fit_completely, // do not add tolerance when replacing
 				threadsMax,                  // final int         threadsMax,
 				debug_level);                // final int         debug_level)
-		dbg_img[9] = disp_outliers.clone();
+		dbg_img[10] = disp_outliers.clone();
 		double [] disp = QuadCLT.blurWeak(
-				new double[][] {disp_outliers, clean_lma}, //final double [][] dls,
+				new double[][] {disp_outliers, dls[2]}, //final double [][] dls,
 				min_strength_blur,    // double      min_strength_blur,
 				min_strength_replace, // double      min_strength_replace,
 				num_blur,             // int         n,
 				tilesX,               // int         width,
 				sigma);               // double      sigma);
-		dbg_img[10] = disp.clone();
+		dbg_img[11] = disp.clone();
 		double [][] ds = {disp,     dls[2]};
+		if (blue_sky != null) { // Temporary, dix - pass blue_sky to fillDisparityStrength()
+			for (int i = 0; i < clean_lma.length; i++) if (blue_sky[i]){
+				ds[0][i] = 0.0;
+				ds[1][i] = 0.0001;
+			}
+		}
+		// ignores results of the last step that produces zero-strength?
 		final double [][] ds_filled = QuadCLT.fillDisparityStrength(
 				ds,                // final double [][] ds0,
 				min_disparity,     // final double      min_disparity,
@@ -11706,8 +11765,15 @@ public double[][] correlateIntersceneDebug( // only uses GPU and quad
 				max_change,        // final double      max_change,
 				tilesX,            // final int         width,
 				threadsMax,        // final int         threadsMax,
-				debug_level); // final int         debug_level)
-		dbg_img[11] = ds_filled[0].clone();
+				debug_level+1); // final int         debug_level)
+		if (blue_sky != null) {
+			for (int i = 0; i < clean_lma.length; i++) if (blue_sky[i]){
+				ds_filled[0][i] = 0.0;
+			}
+		}
+		
+		
+		dbg_img[12] = ds_filled[0].clone();
 		
 		if ((debug_level > 0)) {// && (clt_parameters.ofp.enable_debug_images)) {
 			(new ShowDoubleFloatArrays()).showArrays(
@@ -11751,17 +11817,6 @@ public double[][] correlateIntersceneDebug( // only uses GPU and quad
     			valid_tiles, // final boolean []          valid_tiles,            
     			threadsMax); // final int                 threadsMax)  // maximal number of threads to launch
 		
-		//FIXME:  not clear where tp_tasks was supposed to go? They are not needed, but  valid_tiles - are
-
-		/*
-		scene.getGPU().setInterTasks(
-				pXpYD, // final double [][]         pXpYD, // per-tile array of pX,pY,disparity triplets (or nulls)
-    			scene.getGeometryCorrection(), // final GeometryCorrection  geometryCorrection,
-    			disparity_corr, // final double              disparity_corr,
-    			margin, // final int                 margin,      // do not use tiles if their centers are closer to the edges
-    			valid_tiles, // final boolean []          valid_tiles,            
-    			threadsMax); // final int                 threadsMax)  // maximal number of threads to launch
-    			*/
 		ai.set(0);
 		for (int ithread = 0; ithread < threads.length; ithread++) {
 			threads[ithread] = new Thread() {

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

@@ -482,8 +482,8 @@ public class QuadCLT extends QuadCLTCPU {
 		}
 		return disparity_out;
 	}
-	
-	public static double [][] fillDisparityStrength(
+	@Deprecated
+	public static double [][] fillDisparityStrength0(
 			final double [][] ds0,
 			final double      min_disparity,
 			final double      max_sym_disparity, // lower disparity - num_bottom = 8; 
@@ -494,6 +494,8 @@ public class QuadCLT extends QuadCLTCPU {
 			final int         threadsMax,
 			final int         debug_level)
 	{
+		final double max_disp_diff = 1.00; // no pull from pixels that are this high than average of lowest
+		final int num_bot_avg = 3; // number of lowest neighbors to average
         final double   diagonal_weight = 0.5 * Math.sqrt(2.0); // relative to ortho
 		double wdiag = 0.25 *diagonal_weight / (diagonal_weight + 1.0);
 		double wortho = 0.25 / (diagonal_weight + 1.0);
@@ -509,7 +511,7 @@ public class QuadCLT extends QuadCLTCPU {
 		final Thread[] threads = ImageDtt.newThreadArray(threadsMax);
 		final AtomicInteger ai = new AtomicInteger(0);
 		final AtomicInteger anum_gaps = new AtomicInteger(0);
-		final int dbg_tile = 2191;
+		final int dbg_tile = 89;
 		for (int ithread = 0; ithread < threads.length; ithread++) {
 			threads[ithread] = new Thread() {
 				public void run() {
@@ -566,17 +568,28 @@ public class QuadCLT extends QuadCLTCPU {
 									}
 								}
 							}
+							swd /= sw; // here = not zero;
+
 							System.arraycopy(neibs, 0, neibs_sorted, 0, 8);
 							Arrays.sort(neibs_sorted); // NaNs in the end
 							if (Double.isNaN(neibs_sorted[min_defined - 1])) {
 								continue; // too few defined neighbors
 							}
-							swd /= sw; // here = not zero;
+							double swd_low = 0, sw_low = 0;
+							for (int i = 0; i < num_bot_avg; i++) if (!Double.isNaN(neibs_sorted[i])) {
+								 swd_low += neibs_sorted[i];
+								 sw_low += 1.0;
+							}
+							double max_d = swd_low / sw_low + max_disp_diff;
+							
 							if (!fill_all[0]) {
 								anum_gaps.getAndIncrement();
 							}
 							double max_disp = neibs_sorted[num_bottom - 1];
-//							if (!(ds[0][nTile] > max_sym_disparity)){
+							if (!(max_disp < max_d)) {
+								max_disp = max_d; // this will reduce number of averaged if tries to pull to high 
+							}
+							
 							if (!(swd > max_sym_disparity)){ // here compare to average, not this!
 								max_disp = Double.NaN; // so it will average all
 							}
@@ -598,6 +611,10 @@ public class QuadCLT extends QuadCLTCPU {
 								}
 								amax_diff.accumulate(d2);
 							}
+							if ((debug_level >0) && (nTile == dbg_tile)) {
+								System.out.println("fillDisparityStrength() nTile="+nTile+" swd="+swd);
+							}
+							
 							disp_new[nTile] = swd;		
 						}
 					}
@@ -607,7 +624,7 @@ public class QuadCLT extends QuadCLTCPU {
 			ai.set(0);
 			System.arraycopy(disp_new, 0, ds[0], 0, tiles);
 			if ((debug_level > 0) && fill_all[0]) {
-				System.out.println("fillDisparityStrength() npass="+npass+", change="+Math.sqrt(amax_diff.get())+" ("+max_change+")");
+				System.out.println("fillDisparityStrength() num_gaps="+num_gaps+", npass="+npass+", change="+Math.sqrt(amax_diff.get())+" ("+max_change+")");
 			}
 			if (fill_all[0] && (amax_diff.get() < max_change2)) {
 				break;
@@ -623,109 +640,199 @@ public class QuadCLT extends QuadCLTCPU {
 		return ds;
 	}
 
-	/*
-    public double [][]  getDSRBG (){
-    	return dsrbg;
+	// trying new version
+	public static double [][] fillDisparityStrength(
+			final double [][] ds0,
+			final double      min_disparity,
+			final double      max_sym_disparity, // lower disparity - num_bottom = 8; 
+			final int         num_bottom, // average this number of lowest disparity neighbors (of 8)
+			final int         num_passes,
+			final double      max_change,
+			final int         width,
+			final int         threadsMax,
+			final int         debug_level)
+	{
+		final double max_disp_diff =     1.00; // no pull from pixels that are this high than average of lowest
+		final double max_disp_diff_rel = 0.2; // no pull from pixels that are this high than average of lowest
+		final double anchor_up_abs =     0.2; // if anchor is more than that larger, pull as if it is that latger
+		final double anchor_up_rel =     0.2; // add to anchor_up_abs multiplied by anchor disparity
+		final int    min_float_neibs =   2; // minimal number of float neibs to override fixed;
+		
+        final double   diagonal_weight = 0.5 * Math.sqrt(2.0); // relative to ortho
+		double wdiag = 0.25 *diagonal_weight / (diagonal_weight + 1.0);
+		double wortho = 0.25 / (diagonal_weight + 1.0);
+		final double [] neibw = {wortho, wdiag, wortho, wdiag, wortho, wdiag, wortho, wdiag}; 
+
+		final double max_change2 = max_change * max_change;
+		final double [][] ds = new double [][] {ds0[0].clone(),ds0[1].clone()};
+		final int tiles = ds[0].length;
+		final TileNeibs tn =  new TileNeibs(width, tiles/width);
+		
+		final boolean [] floating =      new boolean[tiles]; // which tiles will change
+		final double [] anchors =        new double [tiles]; // average of the known fixed neighbors
+		final double [] anchor_weights = new double [tiles]; // weights of anchors
+		Arrays.fill(anchors,Double.NaN);
+		
+		final int [] tile_indices = new int [tiles];
+		final Thread[] threads = ImageDtt.newThreadArray(threadsMax);
+		final AtomicInteger ai = new AtomicInteger(0);
+		final AtomicInteger anum_gaps = new AtomicInteger(0);
+		final int dbg_tile = -3379;
+		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 (ds[0][nTile] < min_disparity) {
+							ds[0][nTile] = Double.NaN; // min_disparity;
+						}
+						if (Double.isNaN(ds[0][nTile]) || (ds[1][nTile] <= 0)) { // blue_sky is made small >0 weight
+							ds[0][nTile] = Double.NaN;
+							ds[1][nTile] = 0.0;
+							int indx = anum_gaps.getAndIncrement();
+							tile_indices[indx] = nTile;
+							floating[nTile] = true;
+						}
+					}
+				}
+			};
+		}		      
+		ImageDtt.startAndJoin(threads);
+		ai.set(0);
+		final int num_gaps = anum_gaps.get(); 
+		if (num_gaps == 0) {
+			return ds; // no gaps already
 		}
 		
-    public void setDSRBG(
-			CLTParameters  clt_parameters,
-			int            threadsMax,  // maximal number of threads to launch
-			boolean        updateStatus,
-			int            debugLevel)
-    {
-    	setDSRBG(
-    			this.dsi[is_aux?TwoQuadCLT.DSI_DISPARITY_AUX:TwoQuadCLT.DSI_DISPARITY_MAIN],
-    			this.dsi[is_aux?TwoQuadCLT.DSI_STRENGTH_AUX:TwoQuadCLT.DSI_STRENGTH_MAIN],
-    			this.dsi[is_aux?TwoQuadCLT.DSI_DISPARITY_AUX_LMA:TwoQuadCLT.DSI_DISPARITY_MAIN_LMA],
-    			clt_parameters,
-    			threadsMax,
-    	        updateStatus,
-    			debugLevel);
+		for (int ithread = 0; ithread < threads.length; ithread++) {
+			threads[ithread] = new Thread() {
+				public void run() {
+					double [] neibs = new double[8];
+					for (int indx = ai.getAndIncrement(); indx < num_gaps; indx = ai.getAndIncrement()) {
+						int nTile = tile_indices[indx];
+						if ((debug_level >0) && (nTile == dbg_tile)) {
+							System.out.println("fillDisparityStrength() nTile="+nTile);
+						}
+						Arrays.fill(neibs, Double.NaN);
+						double min_def_disp = Double.NaN;
+						for (int dir = 0; dir < 8; dir++) {
+							int nt_neib = tn.getNeibIndex(nTile, dir);
+							if ((nt_neib >= 0) && !floating[nt_neib]) {
+								neibs[dir] = ds[0][nt_neib];
+								if (!(neibs[dir] >= min_def_disp)) { // true if was NaN or new is smaller
+									min_def_disp = neibs[dir];
+								}
 							}
-    public void setDSRBG(
-    		double [] disparity,
-    		double [] strength,
-    		double [] disparity_lma,
-			CLTParameters  clt_parameters,
-			int            threadsMax,  // maximal number of threads to launch
-			boolean        updateStatus,
-			int            debugLevel)
-    {
-    	double[][] rbg = getTileRBG(
-    			clt_parameters,
-    			disparity,
-    			strength,
-    			disparity_lma,
-    			threadsMax,  // maximal number of threads to launch
-    			updateStatus,
-    			debugLevel);
-    	double [][] dsrbg = {
-    			disparity,
-    			strength,
-    			disparity_lma,
-    			rbg[0],rbg[1],rbg[2]};
-    	this.dsrbg = dsrbg;
 						}
+						if (!Double.isNaN(min_def_disp)) { // no fixed neighbors
+							double max_to_vag = min_def_disp * (1.0 + max_disp_diff_rel) + max_disp_diff;
+							double swd = 0.0, sw = 0.0;
+							for (int dir = 0; dir < 8; dir++) {
+								if (neibs[dir] < max_to_vag) {// NaN OK, will be false
+									sw +=  neibw[dir];
+									swd += neibw[dir] * neibs[dir]; 
+								}
+							}
+							anchors[nTile] = swd / sw; // here = not zero as there will be at least one tile for min_def_disp
+							anchor_weights[nTile] = sw;
+						}
+					}
+				}
+			};
+		}		      
+		ImageDtt.startAndJoin(threads);
 		
-	public double[][] getTileRBG(
-			CLTParameters  clt_parameters,
-			double []      disparity,
-			double []      strength,
-    		double []      disparity_lma,
-			int            threadsMax,  // maximal number of threads to launch
-			boolean        updateStatus,
-			int            debugLevel)
-	{
-		CLTPass3d scan = new CLTPass3d(tp);
-		scan.setCalcDisparityStrength(
-				disparity,
-				strength);
-		if (disparity_lma == null) {
-			scan.resetLMA();
-		} else {
-			scan.setLMA(disparity_lma);
+		ai.set(0);
+		final boolean [] fill_all = {false};
+		final double [] disp_new = ds[0].clone();
+		DoubleAccumulator amax_diff =  new DoubleAccumulator (Double::max, Double.NEGATIVE_INFINITY);
+		for (int npass = 0; npass < num_passes; npass+= fill_all[0]? 1:0 ) { // do not limit initial passes
+			anum_gaps.set(0);
+			amax_diff.reset();
+			for (int ithread = 0; ithread < threads.length; ithread++) {
+				threads[ithread] = new Thread() {
+					public void run() {
+						double [] neibs = new double[8];
+						for (int indx = ai.getAndIncrement(); indx < num_gaps; indx = ai.getAndIncrement()) {
+							int nTile = tile_indices[indx];
+							if ((debug_level >0) && (nTile == dbg_tile)) {
+								System.out.println("fillDisparityStrength() nTile="+nTile);
 							}
-		boolean [] selection = new boolean [disparity.length];
-		for (int i = 0; i < disparity.length; i++) {
-			selection[i] = (!Double.isNaN(disparity[i]) && ((strength == null) || (strength[i] > 0)));
+							if (!fill_all[0] && !Double.isNaN(ds[0][nTile])) {
+								continue; // fill only new
 							}
-		scan.setTileOpDisparity(selection, disparity);
-		// will work only with GPU
-		// reset bayer source, geometry correction/vector
-		//this.new_image_data =     true;
-		QuadCLT savedQuadClt =  gpuQuad.getQuadCLT();
-		if (savedQuadClt != this) {
-			gpuQuad.updateQuadCLT(this);
+							Arrays.fill(neibs, Double.NaN);
+							// average only floating, use anchors for fixed
+							double swd = 0.0, sw = 0.0;
+							double new_val; //  = Double.NaN;
+							int num_floats = 0;
+							for (int dir = 0; dir < 8; dir++) {
+								int nt_neib = tn.getNeibIndex(nTile, dir);
+								if ((nt_neib >= 0) && floating[nt_neib]) {
+									neibs[dir] = ds[0][nt_neib];
+									if (! Double.isNaN(neibs[dir])) {
+										sw +=  neibw[dir];
+										swd += neibw[dir] * neibs[dir];
+										num_floats++;
+									}
+								}
+							}
+							if ((sw > 0) || (anchor_weights[nTile] > 0)) {
+								// (num_floats >= min_float_neibs) to disallow floats BG pull down through just diagonal 
+								if ((sw > 0) && (anchor_weights[nTile] > 0) && (num_floats >= min_float_neibs)) {
+									double avg_flt = swd/sw;
+									double max_fix = avg_flt + anchor_up_abs+ avg_flt*anchor_up_rel;
+									double avg_fix = anchors[nTile];
+									if (avg_fix > max_fix) {
+										avg_fix = max_fix;
+									}
+									new_val = (avg_fix * anchor_weights[nTile] + swd)/(sw + anchor_weights[nTile]);
+								} else if (anchor_weights[nTile] > 0) { // (sw == 0) || (num_floats < min_float_neibs)
+									new_val = anchors[nTile];
+								} else { // OK even if number of float neighbors is <  min_float_neibs - there are no fixed competitors 
+									new_val = swd/sw;
+								}
+								if (fill_all[0]) {
+									double d2 = max_change2 * 2;
+									if (!Double.isNaN(disp_new[nTile])){
+										double d = new_val -  disp_new[nTile];
+										d2 = d * d;
+										if ((debug_level > 0) &&  (d2 > 15)) {
+											System.out.println(
+													"fillDisparityStrength(): nTile="+nTile+
+													" tileX="+(nTile%width)+ " tileY="+(nTile/width)+
+													", d="+d);
+										}
+									}
+									amax_diff.accumulate(d2);
 								} else {
-			savedQuadClt = null;
+									anum_gaps.getAndIncrement();
 								}
-		setPassAvgRBGA( // get image from a single pass, return relative path for x3d // USED in lwir
-				clt_parameters, // CLTParameters           CLTParameters           clt_parameters,,
-				scan,
-				threadsMax,  // maximal number of threads to launch
-				updateStatus,
-				debugLevel);
-		double [][] rgba = scan.getTilesRBGA();
-		if (debugLevel > -1) { // -2) {
-			String title = image_name+"-RBGA";
-			String [] titles = {"R","B","G","A"};
-			(new ShowDoubleFloatArrays()).showArrays(
-					rgba,
-					tp.getTilesX(),
-					tp.getTilesY(),
-					true,
-					title,
-					titles);
+								disp_new[nTile] = new_val;
 							}
-        // Maybe resotore y caller?
-		if (savedQuadClt !=  null) {
-			gpuQuad.updateQuadCLT(savedQuadClt);
 						}
-		return rgba;
 					}
-	
-*/	
+				};
+			}		      
+			ImageDtt.startAndJoin(threads);
+			ai.set(0);
+			System.arraycopy(disp_new, 0, ds[0], 0, tiles);
+			if ((debug_level > 0) && fill_all[0]) {
+				System.out.println("fillDisparityStrength() num_gaps="+num_gaps+", npass="+npass+", change="+Math.sqrt(amax_diff.get())+" ("+max_change+")");
+			}
+			if (fill_all[0] && (amax_diff.get() < max_change2)) {
+				break;
+			}
+			if (anum_gaps.get() == 0) { // no new tiles filled
+				fill_all[0] = true; 
+			}
+			if ((debug_level>0) && (npass == (num_passes-1))){
+				System.out.println("fillDisparityStrength() LAST PASS ! npass="+npass+", change="+Math.sqrt(amax_diff.get())+" ("+max_change+")");
+				System.out.println("fillDisparityStrength() LAST PASS ! npass="+npass+", change="+Math.sqrt(amax_diff.get())+" ("+max_change+")");
+				System.out.println("fillDisparityStrength() LAST PASS ! npass="+npass+", change="+Math.sqrt(amax_diff.get())+" ("+max_change+")");
+			}
+		} // for (int npass = 0; npass < num_passes; npass+= fill_all[0]? 1:0 )
+		return ds;
+	}
 	
 		
 	@Deprecated

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

@@ -95,7 +95,8 @@ public class QuadCLTCPU {
 	public static int             INDEX_DSI_MAIN =  2;
 	public static final String [] FGBG_TITLES_ADJ = {"disparity","strength"};
 //	public static final String [] FGBG_TITLES = {"disparity","strength", "rms","rms-split","fg-disp","fg-str","bg-disp","bg-str"};
-	public static final String [] FGBG_TITLES_AUX = {"disparity","strength", "rms","rms-split","fg-disp","fg-str","bg-disp","bg-str","aux-disp","aux-str"};
+	public static final String [] FGBG_TITLES_AUX = 
+		{"disparity","strength", "rms","rms-split","fg-disp","fg-str","bg-disp","bg-str","aux-disp","aux-str"};
 //	public static final enum      FGBG           {DISPARITY,  STRENGTH,   RMS,  RMS_SPLIT,  FG_DISP,  FG_STR,  BG_DISP,  BG_STR};
 	public static final int       FGBG_DISPARITY = 0;
 	public static final int       FGBG_STRENGTH =  1;
@@ -168,7 +169,7 @@ public class QuadCLTCPU {
     public int                                             num_orient = 0; 
     //number of times scenes are accumulated: 0 - none, 1 - after first orientation, 2 - after second orientation
     public int                                             num_accum =  0;
-    
+    public boolean[] blue_sky =                            null;
     public void inc_orient()        {num_orient++;}
     public void inc_accum()         {num_accum++;}
     public void set_orient(int num) {num_orient = num;}
@@ -408,6 +409,126 @@ public class QuadCLTCPU {
 		return x3d_path;
 	}
 
+	/**
+	 * Discriminate "blue sky" areas with no details at infinity. Such areas
+	 * have both low strength and low pixel value variations between channels
+	 * when calculated for 0 disparity. Assuming that there are no large enough
+	 * featureless areas in the scene itself (not always true, of course).
+	 * So first find tiles with a product of strength and disparity is less than
+	 * sky_seed (<sky_lim), shrink it by sky_shrink (to eliminate small non-infinity
+	 * false positive areas), then expand limited by  sky_lim (expecting a strong
+	 * enough skyline). The (optionally) expand by sky_expand_extra more.
+	 * Expansion by 1 is horizontal/vertical only, by 2 includes diagonals,
+	 * and so on. 
+	 *  
+	 * @param sky_seed minimal value of strength*spread to seed sky areas 
+	 * @param sky_lim maximal value of strength*spread over which sky area will
+	 *                        be expanded
+	 * @param sky_shrink shrink initial sky area to eliminate small non-sky areas. 
+	 * @param sky_expand_extra additionally expand
+	 * @param width number of tiles in a row 
+	 * @param strength 1d array of tile strengths in scanline order
+	 * @param spread 1d array of tile spreads (second maximal difference from
+	 *                  of per-sensor pixel values in a tile to their average
+	 *                  in scanline order. 
+	 * @param debugLevel debug level
+	 * @return boolean 1d array of the pixels belonging to the blue sky.
+	 */
+	public static boolean [] getBlueSky  (
+			double sky_seed, //  =       7.0;  // start with product of strength by diff_second below this
+			double sky_lim, //   =      15.0; // then expand to product of strength by diff_second below this
+			int    sky_shrink, //  =       4;
+			int    sky_expand_extra, //  = 100; // 1?
+			int    width,
+			double [] strength,
+			double [] spread,
+			int debugLevel) {
+		if ((strength == null) || (spread==null)) {
+			return null;
+		}
+		String [] dbg_titles = {"sky", "seed", "max", "shrank"};
+		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 = strength[i] * spread[i];
+			sky_tiles[i] =      (d < sky_seed);
+			prohibit_tiles[i] = (d >= sky_lim);
+		}
+		if (dbg_img != null) {
+			for (int i = 0; i < sky_tiles.length; i++) {
+				dbg_img[1][i] = sky_tiles[i]? 1 : 0;
+				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,
+				null);         // boolean [] prohibit)
+		if (dbg_img != null) {
+			for (int i = 0; i < sky_tiles.length; i++) {
+				dbg_img[3][i] = sky_tiles[i]? 1 : 0;
+			}			
+		}
+		tn.growSelection(
+				2*width ,        // 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,
+				prohibit_tiles); // boolean [] prohibit)
+		if (sky_expand_extra > 0) {
+			tn.growSelection(
+					sky_expand_extra , // 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,
+					null);             // boolean [] prohibit)
+		}
+
+		if (dbg_img != null) {
+			for (int i = 0; i < sky_tiles.length; i++) {
+				dbg_img[0][i] = sky_tiles[i]? 1 : 0;
+			}
+			
+			(new ShowDoubleFloatArrays()).showArrays(
+					dbg_img,
+					width,
+					sky_tiles.length/width,
+					true,
+					"sky_selection",
+					dbg_titles); //	dsrbg_titles);
+		}
+		
+		return sky_tiles;
+	}
+	
+	public boolean [] getBlueSky () {
+		return this.blue_sky;
+	}
+
+	public void setBlueSky (boolean [] blue_sky) {
+		this.blue_sky = blue_sky;
+	}
+	
+	public void setBlueSky (
+			double sky_seed, //  =       7.0;  // start with product of strength by diff_second below this
+			double sky_lim, //   =      15.0; // then expand to product of strength by diff_second below this
+			int    sky_shrink, //  =       4;
+			int    sky_expand_extra, //  = 100; // 1?
+			double [] strength,
+			double [] spread,
+			int debugLevel) {
+		int width = tp.getTilesX();
+		this.blue_sky = getBlueSky  (
+				sky_seed, //  =       7.0;  // start with product of strength by diff_second below this
+				sky_lim, //   =      15.0; // then expand to product of strength by diff_second below this
+				sky_shrink, //  =       4;
+				sky_expand_extra, //  = 100; // 1?
+				width,
+				strength,
+				spread,
+				debugLevel);
+	}
+	
+	
 	public void setDSI(
 			double [][] dsi) {
 		this.dsi = dsi;
@@ -415,7 +536,7 @@ public class QuadCLTCPU {
 
 	public void setDSIFromCombo(
 			double [][] combo_dsi) {
-		this.dsi = new double [TwoQuadCLT.DSI_SLICES.length][];
+		this.dsi = new double [TwoQuadCLT.DSI_SLICES.length][]; // will not have DSI_SPREAD_AUX
 		this.dsi[is_aux?TwoQuadCLT.DSI_DISPARITY_AUX:TwoQuadCLT.DSI_DISPARITY_MAIN] =
 				combo_dsi[OpticalFlow.COMBO_DSN_INDX_DISP_FG];
 		this.dsi[is_aux?TwoQuadCLT.DSI_STRENGTH_AUX:TwoQuadCLT.DSI_STRENGTH_MAIN] =
@@ -473,7 +594,14 @@ public class QuadCLTCPU {
 		}
 		boolean [] reliable = new boolean [strength.length];
 		for (int i = 0; i < reliable.length; i++) {
-			reliable[i] = (strength[i] >= min_strength) && (!needs_lma || !Double.isNaN(disparity_lma[i]));
+			reliable[i] = (strength[i] >= min_strength) &&
+					(!needs_lma || !Double.isNaN(disparity_lma[i]));
+		}
+		boolean [] blue_sky = getBlueSky();
+		if (blue_sky != null) {
+			for (int i = 0; i < reliable.length; i++) {
+				reliable[i] &= !blue_sky[i];
+			}
 		}
 		return reliable;
 	}
@@ -1360,7 +1488,8 @@ public class QuadCLTCPU {
 	{
 		String x3d_path = getX3dDirectory();
 		String title = image_name+TwoQuadCLT.DSI_COMBO_SUFFIX;
-		ImagePlus imp = (new ShowDoubleFloatArrays()).makeArrays(dsi,tp.getTilesX(), tp.getTilesY(),  title, TwoQuadCLT.DSI_SLICES);
+		ImagePlus imp = (new ShowDoubleFloatArrays()).makeArrays(
+				dsi,tp.getTilesX(), tp.getTilesY(),  title, TwoQuadCLT.DSI_SLICES);
 		eyesisCorrections.saveAndShow(
 				imp,      // ImagePlus             imp,
 				x3d_path, // String                path,
@@ -1373,7 +1502,8 @@ public class QuadCLTCPU {
 	public void showDSI(double [][] dsi)
 	{
 		  String title = image_name + TwoQuadCLT.DSI_COMBO_SUFFIX;
-		  (new ShowDoubleFloatArrays()).showArrays(dsi, tp.getTilesX(), tp.getTilesY(), true, title, TwoQuadCLT.DSI_SLICES);
+		  (new ShowDoubleFloatArrays()).showArrays(
+				  dsi, tp.getTilesX(), tp.getTilesY(), true, title, TwoQuadCLT.DSI_SLICES);
 	}
 
 	public void saveDSIMain(){saveDSIMain(this.dsi);}
@@ -1382,10 +1512,15 @@ public class QuadCLTCPU {
 	{
 		String x3d_path = getX3dDirectory();
 		String title = image_name+"-DSI_MAIN";
-		String []   titles =   {TwoQuadCLT.DSI_SLICES[TwoQuadCLT.DSI_DISPARITY_MAIN], TwoQuadCLT.DSI_SLICES[TwoQuadCLT.DSI_STRENGTH_MAIN]};
-		double [][] dsi_main = {dsi[TwoQuadCLT.DSI_DISPARITY_MAIN],        dsi[TwoQuadCLT.DSI_STRENGTH_MAIN]};
+		String []   titles =   {
+				TwoQuadCLT.DSI_SLICES[TwoQuadCLT.DSI_DISPARITY_MAIN],
+				TwoQuadCLT.DSI_SLICES[TwoQuadCLT.DSI_STRENGTH_MAIN]};
+		double [][] dsi_main = {
+				dsi[TwoQuadCLT.DSI_DISPARITY_MAIN],
+				dsi[TwoQuadCLT.DSI_STRENGTH_MAIN]};
 
-		ImagePlus imp = (new ShowDoubleFloatArrays()).makeArrays(dsi_main, tp.getTilesX(), tp.getTilesY(),  title, titles);
+		ImagePlus imp = (new ShowDoubleFloatArrays()).makeArrays(
+				dsi_main, tp.getTilesX(), tp.getTilesY(),  title, titles);
 		eyesisCorrections.saveAndShow(
 				imp,      // ImagePlus             imp,
 				x3d_path, // String                path,
@@ -1400,7 +1535,8 @@ public class QuadCLTCPU {
 	{
 		String x3d_path = getX3dDirectory();
 		String title = image_name+suffix; // "-DSI_MAIN";
-		ImagePlus imp = (new ShowDoubleFloatArrays()).makeArrays(dsi, tp.getTilesX(), tp.getTilesY(),  title, TwoQuadCLT.DSI_SLICES);
+		ImagePlus imp = (new ShowDoubleFloatArrays()).makeArrays(
+				dsi, tp.getTilesX(), tp.getTilesY(),  title, TwoQuadCLT.DSI_SLICES);
 		eyesisCorrections.saveAndShow(
 				imp,      // ImagePlus             imp,
 				x3d_path, // String                path,
@@ -9584,7 +9720,8 @@ public class QuadCLTCPU {
 					  combo_pass, // CLTPass3d   scan,
 					  "after_multi-tile_disparity_extension");
 		  }
-
+		  // copy second_max from the BG pass to the last one (to be used)
+		  tp.clt_3d_passes.get(tp.clt_3d_passes.size()-1).setSecondMax(tp.clt_3d_passes.get(bg_pass).getSecondMax());
  ///// Refining after all added   - end
 		  Runtime.getRuntime().gc();
 	      System.out.println("preExpandCLTQuad3d(): processing  finished at "+

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

@@ -3912,11 +3912,12 @@ ImageDtt.startAndJoin(threads);
 	{
 		double [] disparity =     scan.getDisparity(use_final?0:1);
 		double [] disparityLMA =  scan.getDisparityLMA();
+		double [] second_max_bg = scan.getSecondMax();
 		if (!use_final) {
 			scan.setStrength(null);
 		}
 		double [] strength =     scan.getStrength();
-		return new double [][] {disparity, strength, disparityLMA}; 
+		return new double [][] {disparity, strength, disparityLMA,second_max_bg}; // second maximal difference between channels in BG 
 	}
 
 	// TODO: update for variable length

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

@@ -80,6 +80,8 @@ public class TwoQuadCLT {
 	public static  int DSI_STRENGTH_MAIN =      6;
 	public static  int DSI_STRENGTH_AUX =       7;
 	public static  int DSI_STRENGTH_RIG =       8;
+	public static  int DSI_SPREAD_MAIN =        9;
+	public static  int DSI_SPREAD_AUX =        10;
 
 	public static String DSI_COMBO_SUFFIX = "-DSI_COMBO";
 	public static String DSI_MAIN_SUFFIX =  "-DSI_MAIN";
@@ -93,7 +95,9 @@ public class TwoQuadCLT {
 				"disparity_x3d",
 			    "strength_main",
 			    "strength_aux",
-			    "strength_rig"};
+			    "strength_rig",
+			    "spread_main",
+			    "spread_aux"};
 
 	public long                                            startTime;     // start of batch processing
 	public long                                            startSetTime;  // start of set processing