Started LMA S/N measurements to compare contrast enhancement

src/main/java/com/elphel/imagej/cameras/ThermalColor.java

@@ -112,7 +112,121 @@ public class ThermalColor {
 				0xfff8c4, 0xfff9c7, 0xfff9ca, 0xfff9cd, 0xfffad1, 0xfffad4, 0xfffbd8, 0xfffcdb,
 				0xfffcdf, 0xfffde2, 0xfffde5, 0xfffde8, 0xfffeeb, 0xfffeee, 0xfffef1, 0xfffef4,
 				0xffffff};
-		int [][] palettes = {white_hot_palette, black_hot_palette, iron_palette};
+		int [] pal01_palette = {0x1f77b4,0xff7f0e,
+				0x2ca02c,
+				0xd62728,
+				0x9467bd,
+				0x8c564b,
+				0xe377c2,
+				0x7f7f7f,
+				0xbcbd22,
+				0x17becf};
+		int [] viridis_palette = {
+				0x46317e,
+				0x365b8c,
+				0x277e8e,
+				0x1fa187,
+				0x49c16d,
+				0x9fd938};
+		
+		int [] viridis_mod_palette = {
+				0x401080,
+				0x46317e,
+				0x365b8c,
+				0x277e8e,
+				0x1fa187,
+				0x49c16d,
+				0x9fd938,
+				0xcc9000,
+				0xff6000
+				};
+		
+		int [] coolwarm_palette = {
+				0x6788ed ,
+				0x99bafe ,
+				0xc8d7ef ,
+				0xedd0c1 ,
+				0xf6a789 ,
+				0xe16852};
+		int [] cubehelix_palette = {
+				0x1a2341 ,
+				0x1b6144 ,
+				0x687a30 ,
+				0xc77a7c ,
+				0xcda2e0 ,
+				0xc6e1f1};
+		int [] hsv_palette = {		
+				0xffd400 ,
+				0x4eff00 ,
+				0xff85 ,
+				0x9dff ,
+				0x3600ff ,
+				0xff00ec};
+		int [] hsv2_palette = {		
+				0xffd400 ,
+				0x4eff00 ,
+				0xff85 ,
+				0x9dff ,
+				0x3600ff ,
+				0xff00ec,
+				0xffd400 ,
+				0x4eff00 ,
+				0xff85 ,
+				0x9dff ,
+				0x3600ff ,
+				0xff00ec};
+		int [] hsv3_palette = {		
+				0xffd400 ,
+				0x4eff00 ,
+				0xff85 ,
+				0x9dff ,
+				0x3600ff ,
+				0xff00ec,
+				0xffd400 ,
+				0x4eff00 ,
+				0xff85 ,
+				0x9dff ,
+				0x3600ff ,
+				0xff00ec ,
+				0xffd400 ,
+				0x4eff00 ,
+				0xff85 ,
+				0x9dff ,
+				0x3600ff ,
+				0xff00ec};
+		
+		int [] spectral_palette = {		
+				0xe1514a ,
+				0xfba55c ,
+				0xfee899 ,
+				0xecf7a2 ,
+				0xa1d9a4 ,
+				0x479fb3};
+
+		int [] spectral_r_palette = {		
+				0x479fb3 ,
+				0xa1d9a4 ,
+				0xecf7a2 ,
+				0xfee899 ,
+				0xfba55c ,
+				0xe1514a};
+
+
+		int [][] palettes = {
+				white_hot_palette,
+				black_hot_palette,
+				iron_palette,
+				pal01_palette,
+				viridis_palette,
+				coolwarm_palette,
+				cubehelix_palette,
+				hsv_palette,
+				hsv2_palette,
+				hsv3_palette,
+				viridis_mod_palette,
+				spectral_palette,
+				spectral_r_palette
+		};
 		if (indx <0) indx = 0;
 		else if (indx >= palettes.length) indx = palettes.length - 1;
 		return palettes[indx];

src/main/java/com/elphel/imagej/correction/Eyesis_Correction.java

@@ -7374,36 +7374,57 @@ private Panel panel1,
 			IJ.showMessage("Error","No images selected");
 			return false;
 		}
-		double disparity0 =         0.75;
-		double disparity_max =     75.0; // pix
+		String sky_mask =            "/home/elphel/lwir16-proc/results-cuda/ERS-debug1/models3/1626032208_613623/v01/1626032208_613623-sky_mask.tiff";
+		double disp_inf =           -0.2; // -1.0; // 0.0;            
+		double disparity0 =         0.75; // 2; // 0.75;
+		double disparity_max =     23; // 75.0; // pix
 		boolean show_dialog2 =     false;
 		boolean linear_disparity = false;
+		int     palette =           4;
 		
-		int legend_width = 0;
-		int legend_gap =   5;
+		
+		int legend_width = 2; // 0;
+		int legend_gap =   2; // 5;
 		GenericJTabbedDialog gd0 = new GenericJTabbedDialog("Ln mode");
-		gd0.addNumericField("disparity0",        disparity0,      5,8,"pix", "Disparity to swict from linear to log. ) - skip log mode");
+		gd0.addStringField ("File path of the sky mask",  sky_mask, 80,"optional file with 0.0 - keep, >0.0 - replace with NaN");
+		gd0.addNumericField("Disparity at infinity", disp_inf,      5,8,"pix", "Disparity at infinity ");
+		gd0.addNumericField("disparity0",        disparity0,      5,8,"pix", "Disparity to switch from linear to log. ) - skip log mode");
 		gd0.addNumericField("Maximal disparity", disparity_max,      5,8,"pix", "Fixed maximal disparity (0 - auto");
 		gd0.addCheckbox("Linear disparity legend", linear_disparity);
 		gd0.addNumericField("Legend width",                         legend_width, 0,3,"", "Optional disparity legend vertical bar width");
 		gd0.addNumericField("Legend gap",                           legend_gap, 0,3,"",   "Optional disparity legend vertical bar gap");
+		gd0.addNumericField("palette",                              palette, 0,3,"", "Palette index");
 		gd0.addCheckbox("Show second dialog", show_dialog2);
 		
 		gd0.showDialog();
 		if (gd0.wasCanceled()) return false;
+		sky_mask =            gd0.getNextString();
+		disp_inf =            gd0.getNextNumber();
 		disparity0=           gd0.getNextNumber();
 		disparity_max=        gd0.getNextNumber();
 		linear_disparity =    gd0.getNextBoolean();
 		legend_width=   (int) gd0.getNextNumber();
 		legend_gap=     (int) gd0.getNextNumber();
+		palette=        (int) gd0.getNextNumber();
 		show_dialog2 =        gd0.getNextBoolean();
 
-		boolean log_mode = disparity0 > 0.0;
-		
-		
+		boolean log_mode = disparity0 > 0.00;
+		float [] fsky_mask = null;
+		if (sky_mask.length() > 0) {
+			ImagePlus imp_sky_mask = new ImagePlus(sky_mask);
+			fsky_mask = (float[]) imp_sky_mask.getProcessor().getPixels();
+		}
+//		float [] pixels=(float []) imp_srcs[srcChannel].getProcessor().getPixels();		
 		int current_slice = imp_sel.getCurrentSlice();
 		ImageStack imageStack = imp_sel.getStack();
-		float [] fpixels = (float[]) imageStack.getPixels(current_slice);
+		float [] fpixels0 = (float[]) imageStack.getPixels(current_slice);
+		float [] fpixels = fpixels0.clone();
+		for (int i = 0; i <fpixels.length; i++) {
+			fpixels[i] -= (float) disp_inf;
+			if ((fsky_mask != null) && (fsky_mask[i] > 0.0f)) {
+				fpixels[i] = Float.NaN;
+			}
+		}
 		int width =  imp_sel.getWidth();
 		int height = imp_sel.getHeight();
 		String title = imp_sel.getShortTitle(); // getTitle();
@@ -7441,11 +7462,10 @@ private Panel panel1,
 				}
 			}
 		}
-		double mn = fpixels[0];
+		double mn = dpixels[0];
 		double mx = mn;
 		double pwr = 1.0;
-		int palette = 2;
-		for (int i = 0; i < fpixels.length; i++) {
+		for (int i = 0; i < dpixels.length; i++) {
 			double d = dpixels[i];
 			if (log_mode) {
 				if (!Double.isNaN(d)) {
@@ -7458,7 +7478,8 @@ private Panel panel1,
 					}
 				}
 			}
-			if (!Double.isNaN(d)) {
+			int px = i % width;
+			if ((px < width0) && !Double.isNaN(d)) { // do not use legend
 				if (!(d <= mx)) mx = d;
 				if (!(d >= mn)) mn = d;
 			}
@@ -7483,13 +7504,11 @@ private Panel panel1,
 			gd.addNumericField("min",                             mn,      5,8,"", "Minimal value to map");
 			gd.addNumericField("max",                             mx,      5,8,"", "Maximal value to map");
 			gd.addNumericField("pwr",                             pwr,     5,8,"", "Exponent power");
-			gd.addNumericField("palette",                         palette, 0,3,"", "Palette index");
 			gd.showDialog();
 			if (gd.wasCanceled()) return false;
 			mn=           gd.getNextNumber();
 			mx=           gd.getNextNumber();
 			pwr=          gd.getNextNumber();
-			palette= (int)gd.getNextNumber();
 		}
 		if (pwr != 1.0) {
 			if (mn < 0) {
@@ -7516,6 +7535,8 @@ private Panel panel1,
 			}
 		}
 		
+		System.out.println("mn="+mn+", mx="+mx);
+//		(new ShowDoubleFloatArrays()) .showArrays(dpixels,  width, height, "test_depth_mn="+mn+"_mx="+mx);
 
 		double [][]  pseudo_pixels = new double [3] [dpixels.length];
 		ThermalColor tc = new ThermalColor(

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

@@ -209,6 +209,12 @@ public class Correlation2d {
 	
 	private static int SUB_SAMPLE = 16;  // subsample source pixel in each direction when generating
 	public static int THREADS_MAX = 100; 
+	public int getPairLength(int pair) {
+		return pair_length[pair];
+	}
+	public int [] getPairLengths() {
+		return pair_length;
+	}
 	
 	// All used pairs (but the diameters) are clockwise (end is clockwise of start)
 	// Orientation calculations are valid for clockwise only
@@ -429,7 +435,7 @@ public class Correlation2d {
      * Add 2D correlation for a pair to the combined correlation tile, applying rotation/scaling
      * @param accum_tile tile for accumulation in line-scan order, same dimension as during generateResample()
      * @param corr_tile correlation tile (currently 15x15)
-     * @param num_pair number of correlation pair for which resamplind data exists
+     * @param num_pair number of correlation pair for which resampling data exists
      * @param weight multiply added tile data by this coefficient before accumulation
      */
     public void accummulatePair(
@@ -472,6 +478,26 @@ public class Correlation2d {
     	return sumw;
     }
     
+    public double accummulatePairs(
+    		double []   accum_tile,
+    		double [][] corr_tiles,
+    		double  []   weights) {
+    	double sumw = 0.0;
+    	for (int num_pair = 0; num_pair < weights.length; num_pair++) {
+    		if ((weights[num_pair]>0.0) && (corr_tiles[num_pair] != null)) {
+    			accummulatePair(
+    					accum_tile,
+    		    		corr_tiles[num_pair],
+    		    		num_pair,
+    		    		weights[num_pair]);
+    			sumw+=weights[num_pair];
+    		}
+    	}
+    	return sumw;
+    }
+    
+    
+    
     public void normalizeAccumulatedPairs(
     		double []   accum_tile,
     		double      sumw) {
@@ -483,6 +509,36 @@ public class Correlation2d {
     	}
     }
     
+    /**
+     * Calculate correlation pair width in the disparity direction to boost weights of the most
+     * informative (for disparity) pairs 
+     * @param corr_tile for selected pair 
+     * @param num_pair selected pair number
+     * @return pair width (in pixels) in the disparity direction after applying rotation and scaling
+     */
+    public double getPairWidth(
+    		double [] corr_tile,
+    		int       num_pair){ // Center line only
+     	if ((resample_indices == null) || (resample_indices[num_pair] == null)) {
+    		throw new IllegalArgumentException ("getPairWidth(): No resample data for num_pair = "+num_pair);
+    	}
+    	double s0 = 0.0, s1 = 0.0, s2 = 0.0;
+		int indx = mcorr_comb_width * (mcorr_comb_offset + ((mcorr_comb_height - 1) / 2));
+    	for (int ix = 0; ix < mcorr_comb_width; ix ++ ) if (resample_indices[num_pair][indx] != null){
+    		for (int j = 0; j < resample_indices[num_pair][indx].length; j++) {
+    			double d = corr_tile[resample_indices[num_pair][indx][j]] * resample_weights[num_pair][indx][j];
+    			if (d < 0) {
+    				d = 0.0; // was get5ting negative (s0*s2- s1*s1) !
+    			}
+    			s0 += d;
+    			s1 += d * ix; 
+    			s2 += d * ix * ix; 
+    		}
+    		indx++;
+    	}
+    	return Math.sqrt(s0*s2- s1*s1)/s0;
+    }
+    
     
     public double [] accumulateInit() {return new double [mcorr_comb_width * mcorr_comb_height]; }
     public int getCombWidth() {return mcorr_comb_width;}

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

@@ -917,7 +917,7 @@ public class ErsCorrection extends GeometryCorrection {
 	 * @param camera_xyz camera lens position during centerline acquisition in world coordinates, null - use instance global
 	 * @param camera_atr camera orientation during centerline acquisition in world frame, null - use instance global
 	 * @param line_err iterate until the line (pY) correction is below this value
-	 * @return {px, py, disparity } (right, down) or null if behind the camera
+	 * @return {px, py, disparity } (right, down) or null if behind the camera of the other camera
 	 */
 	
 	public double [] getImageCoordinatesERS(
@@ -1064,7 +1064,8 @@ public class ErsCorrection extends GeometryCorrection {
 	 * @param line_err iterate until the line (pY) correction is below this value
 	 * @return {px, py, disparity } (right, down)
 	 */
-	public double [] getImageCoordinatesERS( // USED in lwir
+	@Deprecated
+	public double [] getImageCoordinatesERS( // not used?
 			double [] xyzw,
 			boolean correctDistortions, // correct distortion (will need corrected background too !)
 			double    line_err) // threshold error in scan lines (1.0)

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

@@ -16193,7 +16193,11 @@ public class ImageDttCPU {
 								}
 							}
 						}
-						boolean debugTile0 =(tileX == debug_tileX) && (tileY == debug_tileY) && (globalDebugLevel > -3);
+						boolean debugTile0 =(tileX == debug_tileX) && (tileY == debug_tileY) && (globalDebugLevel > -6);
+						if (debugTile0) {
+							System.out.println("clt_process_tl_correlations(): tileX="+tileX+", tileY="+tileY);
+						}
+
 						// TODO: move port coordinates out of color channel loop
 						// double [][] centersXY = tp_tasks[iTile].getDoubleXY(); // isAux());
 						double [][] disp_dist = tp_tasks[iTile].getDoubleDispDist();
@@ -16225,11 +16229,55 @@ public class ImageDttCPU {
 						double [] disp_str = {0.0, 0.0}; // diaprity = 0 will be initial approximation for LMA if no averaging
 						if (combine_corrs) {
 							double [] corr_combo_tile = correlation2d.accumulateInit(); // combine all available pairs
-							double sumw = correlation2d.accummulatePairs(
+							double sumw = 0.0;
+							
+							if (imgdtt_params.mcorr_static_weights || imgdtt_params.mcorr_dynamic_weights) {
+								double [] weights = new double [correlation2d.getNumPairs()];
+								if (imgdtt_params.mcorr_static_weights) {
+									int [] pair_lengths = correlation2d.getPairLengths();
+									for (int npair = 0; npair< weights.length; npair++) if (corrs[npair] != null) {
+										weights[npair] = imgdtt_params.mcorr_weights[pair_lengths[npair] - 1];
+									}
+								} else {
+									for (int npair = 0; npair< weights.length; npair++) if (corrs[npair] != null) {
+										weights[npair] = 1.0;
+									}
+//									Arrays.fill(weights, 1.0);
+								}
+								if (imgdtt_params.mcorr_dynamic_weights) {
+									for (int npair = 0; npair< weights.length; npair++) if (corrs[npair] != null) {
+										double pair_width = correlation2d. getPairWidth(
+												corrs[npair] , //double [] corr_tile,
+									    		npair);        // int       num_pair)
+										if (pair_width < 0.1) {
+											System.out.println("pair_width["+npair+"]="+pair_width);
+										} else {
+										weights[npair] /= Math.pow(pair_width, imgdtt_params.mcorr_weights_power);
+										}
+									}
+									
+								}
+								if (debugTile0) {
+									System.out.println("clt_process_tl_correlations(): per-pair weights:");
+									for (int npair = 0; npair< weights.length; npair++) if (weights[npair] > 0.0) {
+										int [] pair = correlation2d.getPair(npair);
+										int pair_length = correlation2d.getPairLength(npair);
+										System.out.println(String.format("%3d: %2d -> %2d [%2d]: %8.5f",npair,pair[0],pair[1],pair_length,weights[npair]));
+									}
+								}
+								sumw = correlation2d.accummulatePairs(
+										corr_combo_tile,      // double []   accum_tile,
+										corrs,                // double [][] corr_tiles, may be longer than selection, extra will be ignored 
+										weights);             // double []     weights);
+							} else { // old way, same weight for all pairs
+								sumw = correlation2d.accummulatePairs(
 										corr_combo_tile,           // double []   accum_tile,
 										corrs,                     // double [][] corr_tiles, may be longer than selection, extra will be ignored 
 										correlation2d.selectAll(), // boolean []  selection,
 										1.0);             // double      weight);
+							}
+							
+							
 							correlation2d.normalizeAccumulatedPairs(
 									corr_combo_tile,
 									sumw);

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

@@ -116,9 +116,10 @@ public class ImageDttParameters {
 	public boolean mcorr_cons_hor =         true;  // consolidate all horizontal pairs
 	public boolean mcorr_cons_vert =        true;  // consolidate all vertical pairs
 	
-	public double [] mcorr_weights =        {1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 16.0}; // static weights of pairs of length 1...8
+//	public double [] mcorr_weights =        {1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 16.0}; // static weights of pairs of length 1...8
+	public double [] mcorr_weights =        {1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 2.0}; // static weights of pairs of length 1...8
 	public boolean mcorr_static_weights =   true; // when mixing , apply static weights to pairs depending on their lengths
-	public double  mcorr_weights_power =    1.0;  // divide pair by horizontal (disparity) width after rotation/scaling (skip negative when measuring width)
+	public double  mcorr_weights_power =    2.0;  // divide pair by horizontal (disparity) width after rotation/scaling (skip negative when measuring width)
 	public boolean mcorr_dynamic_weights =  true; // Apply weights to pairs dependent on the width in disparity direction
 	
 	/// these are just for testing, actual will be generated specifically for different applications (LY will use closest pairs)  

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

+package com.elphel.imagej.tileprocessor;
+
+import java.util.Arrays;
+
+/**
+ **
+ ** InterIntraLMA - Comparing contrast for different interscene and intrascene
+ ** connfifurations
+ **
+ ** Copyright (C) 2021 Elphel, Inc.
+ **
+ ** -----------------------------------------------------------------------------**
+ **
+ **  InterIntraLMA.java is free software: you can redistribute it and/or modify
+ **  it under the terms of the GNU General Public License as published by
+ **  the Free Software Foundation, either version 3 of the License, or
+ **  (at your option) any later version.
+ **
+ **  This program is distributed in the hope that it will be useful,
+ **  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ **  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ **  GNU General Public License for more details.
+ **
+ **  You should have received a copy of the GNU General Public License
+ **  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ ** -----------------------------------------------------------------------------**
+ **
+ */
+
+
+public class InterIntraLMA {
+	/**
+	 * Estimate noise threshold for each combination of inter/intra and number of
+	 * sensors.
+	 * @param noise_file synthetic noise level per file (either rnd or fpn)
+	 * @param senosor_mode_file sensor mode (0 - 16, 1 - 8, 2 - 4, 3 - 2 sensors) per file
+	 * @param inter_file true - interscene , false - intrascene only
+	 * @param good_file_tile true if tile of this file is considered "good"
+	 * @param min_modes minimal number of modes the tile should be defined for (undefine if less than)
+	 * @return per tile per mode average between highest noise of the good tile and lowest noise
+	 * of the bad tile. Double.NaN if noise threshold can not be determined. null if the tile is
+	 * undefined for all modes   
+	 */
+	public static double [][] getNoiseThreshold(
+			double []    noise_file, //  = new double [noise_files.length];
+			int []       sensor_mode_file,
+			boolean []   inter_file,
+			boolean [][] good_file_tile,
+			int          min_modes)
+	{
+		int num_sensor_modes = 0;
+		int num_tiles = good_file_tile[0].length;
+		for (int i = 0; i < sensor_mode_file.length; i++) {
+			if (sensor_mode_file[i] > num_sensor_modes) {
+				num_sensor_modes = sensor_mode_file[i];
+			}
+		}
+		num_sensor_modes ++;
+		int num_modes = 2 * num_sensor_modes;
+		double [][] rslt = new double [num_tiles][]; // number of tiles  
+		double [][][] noise_interval = new double[num_modes][num_tiles][2]; // [modes][tiles]
+		for (int i = 0; i < num_modes; i++) {
+			for (int j= 0; j < num_tiles; j++) {
+				noise_interval[i][j][0] =Double.NaN;
+				noise_interval[i][j][1] =Double.NaN;
+			}
+		}
+		for (int nf = 0; nf < noise_file.length; nf++) {
+			double noise = noise_file[nf];
+			int mode = sensor_mode_file[nf] + (inter_file[nf] ? 0: num_sensor_modes); 
+			for (int ntile = 0; ntile < num_tiles; ntile++) {
+				if (good_file_tile[nf][ntile]) { // good tile
+					if (!(noise <= noise_interval[mode][ntile][0])){ // including Double.isNaN(noise_interval[mode][ntile][0]
+						noise_interval[mode][ntile][0] = noise;
+					}
+				} else { // bad tile
+					if (!(noise >= noise_interval[mode][ntile][1])){ // including Double.isNaN(noise_interval[mode][ntile][1]
+						noise_interval[mode][ntile][1] = noise;
+					}
+				}
+			}
+		}
+		for (int ntile = 0; ntile < num_tiles; ntile++){ 
+			int num_defined = 0;
+			for (int mode = 0; mode < num_modes; mode++) {
+				if (!Double.isNaN(noise_interval[mode][ntile][0]) && !Double.isNaN(noise_interval[mode][ntile][0])) {
+					num_defined++;
+				}
+			}
+			if (num_defined >= min_modes) {
+				rslt[ntile] = new double [num_modes];
+				for (int mode = 0; mode < num_modes; mode++) {
+					if (!Double.isNaN(noise_interval[mode][ntile][0]) && !Double.isNaN(noise_interval[mode][ntile][0])) {
+						rslt[ntile][mode] = 0.5 * (noise_interval[mode][ntile][0] + noise_interval[mode][ntile][1]);
+					} else {
+						rslt[ntile][mode] = Double.NaN;
+					}
+				}
+			}
+		}
+		return rslt;
+	}
+
+}