Generating more data for plots in Libreoffice Calc

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

@@ -46,19 +46,187 @@ public class InterIntraLMA {
 	 * 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 [][][] getNoiseThresholdsPartial(
+			double []    noise_file,     //  = new double [noise_files.length];
+			int []       group_indices,  // last points after last file index
+			int []       sensor_mode_file,
+			boolean []   inter_file,
+			boolean [][] good_file_tile,
+			double       min_inter16_noise_level,
+			boolean      apply_min_inter16_to_inter,
+			int          min_modes,
+			boolean      zero_all_bad,     // (should likely be set!) set noise_level to zero if all noise levels result in bad tiles
+			boolean      all_inter,        // tile has to be defined for all inter
+			boolean      need_same_inter,  // = true; // do not use intra sample if same inter is bad for all noise levels
+			boolean need_same_zero,        // do not use samle if it is bad for zero-noise 
+			int          dbg_tile)
+	{
+//		int num_groups = group_indices.length - 1;
+		int num_variants = 1;
+		for (int i = 0; i < group_indices.length - 1; i++) {
+			int ng = group_indices[i+1] - group_indices[i];
+			if (ng > num_variants) {
+				num_variants = ng;
+			}
+		}
+		double [][][] rslt_partial = new double [num_variants][][];
+		for (int nv = 0; nv < num_variants; nv++) {
+			rslt_partial[nv] = getNoiseThreshold(
+					noise_file,              // double []    noise_file,     //  = new double [noise_files.length];
+					group_indices,           // int []       group_indices,  // last points after last file index
+					nv,                      // int          group_index,    // use this index of same noise value variants (or 0 if that does not exist) 
+					sensor_mode_file,        // int []       sensor_mode_file,
+					inter_file,              // boolean []   inter_file,
+					good_file_tile,          // boolean [][] good_file_tile,
+					min_inter16_noise_level, // double       min_inter16_noise_level,
+					apply_min_inter16_to_inter, // boolean      apply_min_inter16_to_inter,
+					min_modes,               // int          min_modes,
+					zero_all_bad,            // boolean      zero_all_bad,     // set noise_level to zero if all noise levels result in bad tiles
+					all_inter,               // boolean      all_inter,        // tile has to be defined for all inter
+					need_same_inter,         // boolean      need_same_inter, // = true; // do not use intra sample if same inter is bad for all noise levels
+					need_same_zero,          // boolean need_same_zero,        // do not use samle if it is bad for zero-noise 
+					dbg_tile);               // int          dbg_tile);
+		}
+
+		return  rslt_partial;
+	}
+	
+	public static double [][][] getNoiseThresholdsPartial(
+			double []    noise_file,     //  = new double [noise_files.length];
+			int []       group_indices,  // last points after last file index
+			int []       sensor_mode_file,
+			boolean []   inter_file,
+			int          outliers,  // may need do modify algorithm to avoid bias - removing same side (densier) outliers 
+			int          min_keep, // remove less outliers if needed to keep this remain 
+			boolean [][][] good_file_tile_range,
+			double       min_inter16_noise_level,
+			boolean      apply_min_inter16_to_inter,
+			int          min_modes,
+			boolean      zero_all_bad,     // (should likely be set!) set noise_level to zero if all noise levels result in bad tiles
+			boolean      all_inter,        // tile has to be defined for all inter
+			boolean      need_same_inter,  // = true; // do not use intra sample if same inter is bad for all noise levels
+			boolean      need_same_zero,   // do not use samle if it is bad for zero-noise 
+			int          dbg_tile)
+	{
+//		int num_groups = group_indices.length - 1;
+		int num_variants = 1;
+		for (int i = 0; i < group_indices.length - 1; i++) {
+			int ng = group_indices[i+1] - group_indices[i];
+			if (ng > num_variants) {
+				num_variants = ng;
+			}
+		}
+		double [][][] rslt_partial = new double [num_variants][][];
+		for (int nv = 0; nv < num_variants; nv++) {
+			rslt_partial[nv] = getNoiseThreshold(
+					noise_file,              // double []    noise_file,     //  = new double [noise_files.length];
+					group_indices,           // int []       group_indices,  // last points after last file index
+					nv,                      // int          group_index,    // use this index of same noise value variants (or 0 if that does not exist) 
+					sensor_mode_file,        // int []       sensor_mode_file,
+					inter_file,              // boolean []   inter_file,
+					outliers,                // int          outliers,  // may need do modify algorithm to avoid bias - removing same side (densier) outliers 
+					min_keep,                 //int          min_keep, // remove less outliers if needed to keep this remain 
+					good_file_tile_range,    // boolean [][][] good_file_tile_range,
+					min_inter16_noise_level, // double       min_inter16_noise_level,
+					apply_min_inter16_to_inter, // boolean      apply_min_inter16_to_inter,
+					min_modes,               // int          min_modes,
+					zero_all_bad,            // boolean      zero_all_bad,     // set noise_level to zero if all noise levels result in bad tiles
+					all_inter,               // boolean      all_inter,        // tile has to be defined for all inter
+					need_same_inter,         //boolean      need_same_inter, // = true; // do not use intra sample if same inter is bad for all noise levels
+					need_same_zero,          // boolean need_same_zero,        // do not use samle if it is bad for zero-noise 
+					dbg_tile);               // int          dbg_tile);
+		}
+
+		return  rslt_partial;
+	}
+	
+	public static double [][] mergeNoiseVariants(
+			double [][][] partial_thresholds,
+			int           num_outliers,
+			int           min_remain){
+		int num_modes = 0;
+		int num_tiles = partial_thresholds[0].length;
+		for (int i = 0; i < num_tiles; i++) {
+			if (partial_thresholds[0][i]!= null) {
+				num_modes = partial_thresholds[0][i].length;
+				break;
+			}
+		}
+		double [][] rslt = new double [num_tiles][]; // number of tiles
+		for (int ntile = 0; ntile < num_tiles; ntile++) {
+			boolean has_data = false;
+			for (int nv = 0; nv < partial_thresholds.length; nv++) {
+				if (partial_thresholds[nv][ntile] != null) {
+					has_data = true;
+					break;
+				}
+			}
+			if (has_data) {
+				rslt[ntile] = new double [num_modes];
+				for (int mode = 0; mode < num_modes; mode ++) {
+					double [] partial_tile = new double [partial_thresholds.length];
+					Arrays.fill(partial_tile, Double.NaN);
+					for (int nv = 0; nv < partial_tile.length; nv++) if (partial_thresholds[nv][ntile] != null){
+						partial_tile[nv] = partial_thresholds[nv][ntile][mode];
+					}
+					int num_defined = 0;
+					double avg = Double.NaN;
+					double s = 0.0;
+					for (int nv = 0; nv < partial_tile.length; nv++) if (!Double.isNaN(partial_tile[nv])) {
+						s+= partial_tile[nv];
+						num_defined++;
+					}
+					if (num_defined > 0) {
+						avg = s / num_defined;
+						for (int num_removed = 0; num_removed < num_outliers; num_removed++) {
+							if (num_defined <= min_remain) {
+								break;
+							}
+							double max_diff2 = 0;
+							int ioutlier = -1;
+							for (int nv = 0; nv < partial_tile.length; nv++) if (!Double.isNaN(partial_tile[nv])) {
+								double diff2 = partial_tile[nv] - s;
+								diff2 *= diff2;
+								if (diff2 > max_diff2) {
+									max_diff2 = diff2;
+									ioutlier = nv;
+								}
+							}							
+							if (ioutlier < 0) {
+								break;
+							}
+							s -=                     partial_tile[ioutlier];
+							partial_tile[ioutlier] = Double.NaN;
+							num_defined--;
+							avg = s / num_defined;
+						}
+					}
+					rslt[ntile][mode] = avg;
+				}
+			}
+		}
+		return rslt;
+	}
+	
 	public static double [][] getNoiseThreshold(
-			double []    noise_file, //  = new double [noise_files.length];
+			double []    noise_file,     //  = new double [noise_files.length];
+			int []       group_indices,  // last points after last file index
+			int          group_index,    // use this index of same noise value variants (or 0 if that does not exist) 
 			int []       sensor_mode_file,
 			boolean []   inter_file,
 			boolean [][] good_file_tile,
 			double       min_inter16_noise_level,
+			boolean      apply_min_inter16_to_inter,
 			int          min_modes,
 			boolean      zero_all_bad,     // set noise_level to zero if all noise levels result in bad tiles
 			boolean      all_inter,        // tile has to be defined for all inter
-			boolean      need_same_inter, // = true; // do not use intra sample if same inter is bad for all noise levels
+			boolean      need_same_inter,  // = true; // do not use intra sample if same inter is bad for all noise levels
+			boolean      need_same_zero,   // do not use samle if it is bad for zero-noise 
 			int          dbg_tile)
 			
 	{
+		// min_inter16_noise_level
+		int num_groups = group_indices.length - 1;
 //		int dbg_tile = 828; // 1222;
 		int num_sensor_modes = 0;
 		int num_tiles = good_file_tile[0].length;
@@ -77,12 +245,17 @@ public class InterIntraLMA {
 				noise_interval[i][j][1] =Double.NaN;
 			}
 		}
-		for (int nf = 0; nf < noise_file.length; nf++) {
+//		for (int nf = 0; nf < noise_file.length; nf++) {
+		for (int ng = 0; ng < num_groups; ng++) {
+			int nf = group_indices[ng] +  group_index;
+			if (nf >= group_indices[ng+1]) {
+				nf = group_indices[ng+1] -1;
+			}
 			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 (ntile == dbg_tile) {
-					System.out.println("ntile = "+ntile+", nf ="+nf);
+					System.out.println("ntile = "+ntile+", nf ="+nf+", mode = "+mode);
 				}
 				if (good_file_tile[nf][ntile]) { // good tile
 					if (!(noise <= noise_interval[mode][ntile][0])){ // including Double.isNaN(noise_interval[mode][ntile][0]
@@ -115,20 +288,15 @@ public class InterIntraLMA {
 			if ((num_defined >= min_modes) && (!all_inter || (num_defined_inter >= 4))) {
 				rslt[ntile] = new double [num_modes];
 				for (int mode = 0; mode < num_modes; mode++) {
-//					if (need_same_inter && (mode >= 4) && Double.isNaN(noise_interval[mode & 3][ntile][0])) {  // no good for same sensors inter
 					if (need_same_inter && Double.isNaN(noise_interval[mode & 3][ntile][0])) { // no good for same sensors inter
 						rslt[ntile][mode] = Double.NaN;
+			        } else if (need_same_zero && (Double.isNaN(noise_interval[mode][ntile][0])
+			        		|| (noise_interval[mode][ntile][1] == 0.0))) { // no good for same sensors
+			        	rslt[ntile][mode] = Double.NaN;
 					} else 	if (!Double.isNaN(noise_interval[mode][ntile][0]) && !Double.isNaN(noise_interval[mode][ntile][1])) {
-						/*
-						rslt[ntile][mode] = 0.5 * (noise_interval[mode][ntile][0] + noise_interval[mode][ntile][1]);
-						if (remove_non_monotonic && (noise_interval[mode][ntile][0] > noise_interval[mode][ntile][1])) {
-							rslt[ntile][mode] = Double.NaN;
-						}
-						*/
 						// use the lowest failed noise level assuming that false positive may happen even for much higher noise level 
 						rslt[ntile][mode] = noise_interval[mode][ntile][1]; // lowest noise for bad
 						
-//					} else if (zero_all_bad && Double.isNaN(noise_interval[mode][ntile][1])) {
 					} else if (zero_all_bad && Double.isNaN(noise_interval[mode][ntile][0])) {
 						rslt[ntile][mode] = 0.0;
 					} else {
@@ -138,7 +306,13 @@ public class InterIntraLMA {
 			}
 			if ((rslt[ntile] != null) && (min_inter16_noise_level >0)){ // filter by to weak inter-16 (mode 0)
 				if (!(rslt[ntile][0] >= min_inter16_noise_level)){
-					rslt[ntile] = null;
+					if (apply_min_inter16_to_inter) {
+						rslt[ntile] = null;
+					} else {
+						for (int mode = 4; mode < rslt[ntile].length; mode++) {
+							rslt[ntile][mode] = Double.NaN;
+						}
+					}
 				}
 			}
 			if (rslt[ntile] != null) {
@@ -172,21 +346,26 @@ public class InterIntraLMA {
 	// trying multi-threshold good_file_tile_range
 	public static double [][] getNoiseThreshold(
 			double []      noise_file, //  = new double [noise_files.length];
+			int []         group_indices,  // last points after last file index
+			int            group_index,    // use this index of same noise value variants (or 0 if that does not exist) 
 			int []         sensor_mode_file,
 			boolean []     inter_file,
 			int            outliers,  // may need do modify algorithm to avoid bias - removing same side (densier) outliers 
 			int            min_keep, // remove less outliers if needed to keep this remain 
 			boolean [][][] good_file_tile_range,
 			double         min_inter16_noise_level,
+			boolean        apply_min_inter16_to_inter,
 			int            min_modes,
 			boolean        zero_all_bad,     // set noise_level to zero if all noise levels result in bad tiles
 			boolean        all_inter,        // tile has to be defined for all inter
 			boolean        need_same_inter, // = true; // do not use intra sample if same inter is bad for all noise levels
+			boolean        need_same_zero,        // do not use samle if it is bad for zero-noise
 			int            dbg_tile)
 			
 	{
 		//int dbg_tile = 828;
-		int num_sensor_modes = 00;
+		int num_groups = group_indices.length - 1;
+		int num_sensor_modes = 0;
 		int num_tiles = good_file_tile_range[0].length;
 		for (int i = 0; i < sensor_mode_file.length; i++) {
 			if (sensor_mode_file[i] > num_sensor_modes) {
@@ -203,7 +382,12 @@ public class InterIntraLMA {
 				lowest_all_bad[i][j] =Double.NaN;
 			}
 		}
-		for (int nf = 0; nf < noise_file.length; nf++) {
+//		for (int nf = 0; nf < noise_file.length; nf++) {
+		for (int ng = 0; ng < num_groups; ng++) {
+			int nf = group_indices[ng] +  group_index;
+			if (nf >= group_indices[ng+1]) {
+				nf = group_indices[ng+1] -1;
+			}
 			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++) {
@@ -285,7 +469,14 @@ public class InterIntraLMA {
 						double [] pre_rslt = new double [noise_intervals[mode][ntile].length];  //null pointer
 						int num_def = 0;
 						for (int stp = 0; stp < pre_rslt.length; stp++) {
-							if (need_same_inter && Double.isNaN(noise_intervals[mode & 3][ntile][stp][0])) { // no good for same sensors inter
+							if (need_same_inter && (
+									(noise_intervals[mode & 3][ntile] == null) ||
+									Double.isNaN(noise_intervals[mode & 3][ntile][stp][0]))) { // no good for same sensors inter
+								pre_rslt[stp] = Double.NaN;
+							} else if (need_same_zero && (
+									(noise_intervals[mode][ntile] == null) || 
+									Double.isNaN(noise_intervals[mode][ntile][stp][0]) ||
+									(noise_intervals[mode][ntile][stp][1] == 0.0))) { // bad for no- noise for same sensors
 								pre_rslt[stp] = Double.NaN;
 							} else 	if (!Double.isNaN(noise_intervals[mode][ntile][stp][0]) && !Double.isNaN(noise_intervals[mode][ntile][stp][1])) {
 								pre_rslt[stp] = noise_intervals[mode][ntile][stp][1]; // lowest noise for bad
@@ -344,16 +535,29 @@ public class InterIntraLMA {
 						} else {
 							rslt[ntile][mode] = Double.NaN;
 						}
-					} else {
-						rslt[ntile][mode] = zero_all_bad ? 0.0 : Double.NaN; // no good in any stp
+					} else { // bad for all stp
+						if (need_same_inter && (noise_intervals[mode & 3][ntile] == null))
+							rslt[ntile][mode] = Double.NaN;
+						else if (need_same_zero && (noise_intervals[mode][ntile] == null)){
+							rslt[ntile][mode] = Double.NaN;
+						} else {
+							rslt[ntile][mode] = zero_all_bad ? 0.0 : Double.NaN; // no good in any stp
+						}
 					}
 				}
 			}
 			if ((rslt[ntile] != null) && (min_inter16_noise_level >0)){ // filter by to weak inter-16 (mode 0)
 				if (!(rslt[ntile][0] >= min_inter16_noise_level)){
-					rslt[ntile] = null;
+					if (apply_min_inter16_to_inter) {
+						rslt[ntile] = null;
+					} else {
+						for (int mode = 4; mode < rslt[ntile].length; mode++) {
+							rslt[ntile][mode] = Double.NaN;
+						}
+					}
 				}
 			}
+			
 			if (rslt[ntile] != null) {
 				boolean all_nan = true;
 				boolean has_nan = false;
@@ -427,6 +631,7 @@ public class InterIntraLMA {
 			double      offset, // for "relative" noise
 			double      n0,    // initial value for N0 0.02
 			int         tilesX, // debug images only
+			double      scale_intra, // scale weight of intra-scene samples ( < 1.0)
 			int         debug_level)
 	{
 		boolean debug_img =  (debug_level > -1);
@@ -473,14 +678,24 @@ public class InterIntraLMA {
 		
 		
 		// create Y, K and weights vectors
+		// scale_intra
+		double sum_w = 0.0;
 		for (int nsample = 0; nsample < num_samples; nsample++) {
 			int tile = tile_index[sample_indx[nsample][0]];
+			int mode = sample_indx[nsample][1];
+			double w = (mode >= 4) ? scale_intra : 1.0;
 			double d = noise_thresh[tile][sample_indx[nsample][1]];
 			K[nsample] = 1.0/(d + offset);
 			Y[nsample] = d * K[nsample];
 			// may be modified, but sum (weights) should be == 1.0;
-			weights[nsample] = 1.0/num_samples;
+			weights[nsample] = w; // 1.0/num_samples;
+			sum_w += w;
 		}
+		double scale_w = 1.0/sum_w;
+		for (int nsample = 0; nsample < num_samples; nsample++) {
+			weights[nsample] *= scale_w;
+		}
+		
 		// initial approximation
 		double N0 =  n0; // offset;
 		double N02 = N0*N0;
@@ -727,8 +942,10 @@ public class InterIntraLMA {
 					if (adjust_N0) {
 						jt[indx++][i] = - K[i];
 					}
-					if (adjust_Gi && (mode > 0)) {
-						jt[indx + mode -1][i] = K[i] *st[itile];
+					if (adjust_Gi) {
+						if (mode > 0) {
+							jt[indx + mode -1][i] = K[i] *st[itile];
+						}
 						indx += gi.length -1;
 					}
 					if (adjust_St) {
@@ -753,8 +970,11 @@ public class InterIntraLMA {
 							jt[indx++][i] = - Amti * n0;
 						}
 						double asg = Amti*st[itile]*gi[mode];
-						if (adjust_Gi && (mode > 0)) {
-							jt[indx + mode -1][i] = asg *st[itile];
+//						if (adjust_Gi && (mode > 0)) {
+						if (adjust_Gi) {
+							if  (mode > 0) {
+								jt[indx + mode -1][i] = asg *st[itile];
+							}
 							indx += gi.length -1;
 						}
 						if (adjust_St) {
@@ -767,6 +987,140 @@ public class InterIntraLMA {
 		return fx;
 	}
 	
+	
+	private boolean debugJt(
+    		double      delta,
+    		double []   vector) {
+    	int num_points = sample_indx.length;
+    	int num_pars = vector.length;
+
+//    	delta = 0.001;
+
+    	double [][] jt =          new double [num_pars][num_points];
+    	double [][] jt_delta =    new double [num_pars][num_points];
+    	double [][] jt_diff =     new double [num_pars][num_points];
+    	double [] max_diff =      new double [num_pars];
+    	boolean [] has_nan =      new boolean [num_pars];
+    	int    [] max_diff_indx = new int [num_pars];
+    	double    worst_diff =  0.0;
+    	int       worst_par  = -1;
+    	boolean   has_any_nan = false;
+    	double [] fx = getFxJt( vector,jt);
+    	if (fx == null) return false;
+    	if (getFxJt(delta, vector,jt_delta) == null) return false;
+    	
+    	for (int npar = 0; npar < jt.length; npar++) {
+    		for (int npoint = 0; npoint < jt[npar].length; npoint++) {
+    			jt_diff[npar][npoint] = jt[npar][npoint] - jt_delta[npar][npoint];
+    			if (Double.isNaN(jt_diff[npar][npoint])) {
+    				has_nan[npar] = true;
+    			} else {
+    				if (Math.abs(jt_diff[npar][npoint]) > max_diff[npar]) {
+    					max_diff[npar] =      Math.abs(jt_diff[npar][npoint]);
+    					max_diff_indx[npar] = npoint;
+    				}
+    			}
+    		}
+    		has_any_nan |= has_nan[npar];
+    		if (max_diff[npar] > worst_diff) {
+    			worst_diff = max_diff[npar];
+    			worst_par = npar;
+    		}
+    	}
+    	System.out.println("Has NaNs = "+has_any_nan);
+    	if (has_any_nan) {
+    		for (int i = 0; i < has_nan.length; i++) {
+    			if (has_nan[i]) {
+    				System.out.print(i+", ");
+    			}
+    		}
+    		System.out.println();
+    	}
+    	System.out.println("Worst diff = "+worst_diff+" for parameter #"+worst_par+", point "+max_diff_indx[worst_par]);
+		System.out.println();
+    	/*
+    	
+    	System.out.println("Test of jt-jt_delta difference,  delta = "+delta+ ":");
+    	System.out.print(String.format("Til P %3s: %10s ", "#", "fx"));
+    	for (int anp = 0; anp< all_pars.length; anp++) if(par_mask[anp]){
+			String parname;
+			if      (anp >= ndisp_index) parname = PAR_NAME_CORRNDISP + (anp - ndisp_index);
+			else if (anp >= ddisp_index) parname = PAR_NAME_CORRDISP +  (anp - ddisp_index);
+			else {
+				int ntile = anp / tile_params;
+				int anpr =  anp % tile_params;
+				if (anpr < G0_INDEX) parname = PAR_NAMES[anpr]+"-"+ntile;
+				else                 parname = PAR_NAME_SCALE +"-"+ntile + ":"+ (anpr - G0_INDEX);
+			}
+        	System.out.print(String.format("| %16s ", parname));
+    	}
+    	System.out.println();
+    	int npair0 = -1;
+    	for (int i = 0; i < num_points; i++) {
+    		if (i < samples.size()) {
+//    			int [] fs = correlation2d.getPair(samples.get(i).pair);
+//        		int npair = used_pairs_map[samples.get(i).tile][samples.get(i).fcam][samples.get(i).scam];
+//        		int npair = used_pairs_map[samples.get(i).tile][fs[0]][fs[1]];
+        		int npair = used_pairs_map[samples.get(i).tile][samples.get(i).pair];
+        		if (npair !=npair0) {
+        			if (npair0 >=0) System.out.println();
+        			npair0 = npair;
+        		}
+    			System.out.print(String.format("%3d %1d %3d: %10.7f ",samples.get(i).tile, npair, i, fx[i]));
+    		} else {
+            	System.out.print(String.format(" -  - %3d: %10.7f ", i, fx[i]));
+    		}
+        	for (int np = 0; np < num_pars; np++) {
+//            	System.out.print(String.format("|%8.5f %8.5f ", jt_delta[np][i], 1000*(jt[np][i] - jt_delta[np][i])));
+            	System.out.print(String.format("|%8.5f %8.5f ", jt_delta[np][i], 1.0 * (jt[np][i] - jt_delta[np][i])));
+            	double adiff = Math.abs(jt[np][i] - jt_delta[np][i]);
+            	if (adiff > max_diff[np]) {
+            		max_diff[np] = adiff;
+            	}
+        	}
+        	System.out.println();
+    	}
+    	double tmd = 0.0;
+    	for (int np = 0; np < num_pars; np++) {
+    		if (max_diff[np] > tmd) tmd= max_diff[np];
+    	}
+//        System.out.print(String.format("      %15s ", "Maximal diff:"));
+        System.out.print(String.format("Max diff.(%10.5f):", tmd));
+
+    	for (int np = 0; np < num_pars; np++) {
+        	System.out.print(String.format("|%8s %8.5f ", "1/1000×",  1000*max_diff[np]));
+    	}
+    	System.out.println();
+    	*/
+    	return true;
+    }
+
+
+	public double [] getFxJt( // not used in lwir
+			double      delta, // for testing derivatives: calculates as delta-F/delta_x
+			double []   vector,
+			double [][] jt) { // should be either [vector.length][samples.size()] or null - then only fx is calculated
+		double [] fx0=getFxJt(vector,null);
+		if (fx0 == null) return null;
+		for (int np = 0; np < vector.length; np++) {
+			double [] vector1 = vector.clone();
+			vector1[np]+= delta;
+			double [] fxp=getFxJt(vector1,null);
+			if (fxp == null) return null;
+			vector1 = vector.clone();
+			vector1[np]-= delta;
+			double [] fxm=getFxJt(vector1,null);
+			if (fxm == null) return null;
+			jt[np] = new double [fxp.length];
+			for (int i = 0; i < fxp.length; i++) {
+				jt[np][i] = (fxp[i] - fxm[i])/delta/2;
+			}
+		}
+		return fx0;
+	}
+	
+	
+	
 	public double [][] getYDbg() {
 		double [][] dbg_Y = new double [gi.length][dbgTilesX*dbgTilesY];
 		for (int mode = 0; mode < dbg_Y.length; mode++) {
@@ -915,7 +1269,11 @@ public class InterIntraLMA {
 							if (i == j) {
 								wjtjl[i][j] += d * lambda;
 								if (d == 0) {
-									System.out.println("Diagonal ZERO for i=j="+i+" absolute tile = "+tile_index[i-8]); // assuming N0, gi[1]...gi[7]
+									if (i >= 8) {
+										System.out.println("Diagonal ZERO for i=j="+i+" absolute tile = "+tile_index[i-8]); // assuming N0, gi[1]...gi[7]
+									}else {
+										System.out.println("Diagonal ZERO for i=j="+i); // assuming N0, gi[1]...gi[7]
+									}
 									wjtjl[i][j] = 1.0; // Jt * (y-fx) will anyway be 0, so any value here should work.
 								}
 							} else {
@@ -1102,13 +1460,11 @@ public class InterIntraLMA {
 						true,
 						"dbg_Fx");
 			}
-			/*
 			if (debug_level > 3) {
 				debugJt(
 						0.000001, // double      delta, // 0.2, //
 						this.vector); // double []   vector);
 			}
-			*/
 		}
 		Matrix y_minus_fx_weighted = new Matrix(last_ymfx, last_ymfx.length);
 

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

@@ -2010,11 +2010,7 @@ public class QuadCLT extends QuadCLTCPU {
 			EyesisCorrectionParameters.RGBParameters        rgbParameters,
 			final int                                       threadsMax,  // maximal number of threads to launch
 			final int                                       debugLevel){
-		String x3d_path= correctionsParameters.selectX3dDirectory( // for x3d and obj
-				correctionsParameters.getModelName(image_name), // quad timestamp. Will be ignored if correctionsParameters.use_x3d_subdirs is false
-				correctionsParameters.x3dModelVersion,
-				true,  // smart,
-				true);  //newAllowed, // save
+		String x3d_path = getX3dDirectory();
 		String file_name = image_name + suffix;
 		String file_path = x3d_path + Prefs.getFileSeparator() + file_name + ".tiff";
 		if ((getGPU() != null) && (getGPU().getQuadCLT() != this)) {
@@ -2224,11 +2220,7 @@ public class QuadCLT extends QuadCLTCPU {
 
 		if (clt_parameters.gen_4_img) { // save 4 JPEG images
 			// Save as individual JPEG images in the model directory
-			String x3d_path= correctionsParameters.selectX3dDirectory(
-					image_name, // quad timestamp. Will be ignored if correctionsParameters.use_x3d_subdirs is false
-					correctionsParameters.x3dModelVersion,
-					true,  // smart,
-					true);  //newAllowed, // save
+			String x3d_path = getX3dDirectory();
 			for (int sub_img = 0; sub_img < imps_RGB.length; sub_img++){
 				EyesisCorrections.saveAndShow(
 						imps_RGB[sub_img],
@@ -2759,11 +2751,7 @@ public class QuadCLT extends QuadCLTCPU {
 
 		if (clt_parameters.gen_4_img) { // save 4 JPEG images
 			// Save as individual JPEG images in the model directory
-			String x3d_path= quadCLT_main.correctionsParameters.selectX3dDirectory(
-					name, // quad timestamp. Will be ignored if correctionsParameters.use_x3d_subdirs is false
-					quadCLT_main.correctionsParameters.x3dModelVersion,
-					true,  // smart,
-					true);  //newAllowed, // save
+			String x3d_path= quadCLT_main.getX3dDirectory();
 			for (int sub_img = 0; sub_img < imps_RGB.length; sub_img++){
 				EyesisCorrections.saveAndShow(
 						imps_RGB[sub_img],

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

@@ -329,6 +329,26 @@ public class QuadCLTCPU {
 		return image_name;
 	}
     
+	public String getX3dDirectory() { // replace direct calculations  
+		String x3d_path = correctionsParameters.selectX3dDirectory( // for x3d and obj
+				correctionsParameters.getModelName(image_name), // quad timestamp. Will be ignored if correctionsParameters.use_x3d_subdirs is false
+				correctionsParameters.x3dModelVersion,
+				true,  // smart,
+				true);  //newAllowed, // save
+		return x3d_path;
+	}
+
+	// maybe will not be needed? TODO: Check
+	public String getX3dDirectory(String name) { // replace direct calculations  
+		String x3d_path = correctionsParameters.selectX3dDirectory( // for x3d and obj
+				name,       // quad timestamp. Will be ignored if correctionsParameters.use_x3d_subdirs is false
+				correctionsParameters.x3dModelVersion,
+				true,  // smart,
+				true);  //newAllowed, // save
+		return x3d_path;
+	}
+	
+	
 	public int restoreDSI(
 			String suffix,
 			boolean silent) // "-DSI_COMBO", "-DSI_MAIN" (DSI_COMBO_SUFFIX, DSI_MAIN_SUFFIX)
@@ -344,11 +364,7 @@ public class QuadCLTCPU {
 			String suffix, // "-DSI_COMBO", "-DSI_MAIN" (DSI_COMBO_SUFFIX, DSI_MAIN_SUFFIX)
 			double [][] dsi,
 			boolean silent) {
-		String x3d_path= correctionsParameters.selectX3dDirectory( // for x3d and obj
-				correctionsParameters.getModelName(image_name), // quad timestamp. Will be ignored if correctionsParameters.use_x3d_subdirs is false
-				correctionsParameters.x3dModelVersion,
-				true,  // smart,
-				true);  //newAllowed, // save
+		String x3d_path = getX3dDirectory();
 		String file_path = x3d_path + Prefs.getFileSeparator() + image_name + suffix + ".tiff";
 		ImagePlus imp = null;
 		try {
@@ -404,12 +420,8 @@ public class QuadCLTCPU {
 
 		}
 		if (!path.contains(Prefs.getFileSeparator())) {
-			  String x3d_path= correctionsParameters.selectX3dDirectory( // for x3d and obj
-					  correctionsParameters.getModelName(image_name), // quad timestamp. Will be ignored if correctionsParameters.use_x3d_subdirs is false
-					  correctionsParameters.x3dModelVersion,
-						  true,  // smart,
-						  true);  //newAllowed, // save
-			  path = x3d_path+Prefs.getFileSeparator()+path;
+			String x3d_path = getX3dDirectory();
+			path = x3d_path+Prefs.getFileSeparator()+path;
 		}
 		Properties	inter_properties = new Properties();
 		String prefix = is_aux?PREFIX_AUX:PREFIX; 
@@ -465,12 +477,8 @@ public class QuadCLTCPU {
 
 		}
 		if (!path.contains(Prefs.getFileSeparator())) {
-			  String x3d_path= correctionsParameters.selectX3dDirectory( // for x3d and obj
-					  correctionsParameters.getModelName(image_name), // quad timestamp. Will be ignored if correctionsParameters.use_x3d_subdirs is false
-					  correctionsParameters.x3dModelVersion,
-						  true,  // smart,
-						  true);  //newAllowed, // save
-			  path = x3d_path+Prefs.getFileSeparator()+path;
+			String x3d_path = getX3dDirectory();
+			path = x3d_path+Prefs.getFileSeparator()+path;
 		}
 		properties = loadProperties(
 				path, // String path,
@@ -827,11 +835,7 @@ public class QuadCLTCPU {
 		final int num_cams = this.image_data.length;
 		final int num_cols = image_data[0].length;
 		final int [] image_wh = geometryCorrection.getSensorWH();
-		String x3d_path= correctionsParameters.selectX3dDirectory( // for x3d and obj
-				correctionsParameters.getModelName(image_name), // quad timestamp. Will be ignored if correctionsParameters.use_x3d_subdirs is false
-				correctionsParameters.x3dModelVersion,
-				true,  // smart,
-				true);  //newAllowed, // save
+		String x3d_path = getX3dDirectory();
 		String noise_suffix = suffix + sigma;
 		String file_name = image_name + noise_suffix;
 		String file_path = x3d_path + Prefs.getFileSeparator() + file_name + ".tiff";
@@ -962,11 +966,7 @@ public class QuadCLTCPU {
 			int         width,
 			int         height)
 	{
-		String x3d_path= correctionsParameters.selectX3dDirectory( // for x3d and obj
-				correctionsParameters.getModelName(image_name), // quad timestamp. Will be ignored if correctionsParameters.use_x3d_subdirs is false
-				correctionsParameters.x3dModelVersion,
-				true,  // smart,
-				true);  //newAllowed, // save
+		String x3d_path = getX3dDirectory();
 		String file_name = image_name + suffix;
 		String file_path = x3d_path + Prefs.getFileSeparator() + file_name + ".tiff";
 		ImageStack imageStack = (new ShowDoubleFloatArrays()).makeStack(data, width, height, labels);
@@ -984,11 +984,7 @@ public class QuadCLTCPU {
 			)
 	{
 //		final int [] image_wh = geometryCorrection.getSensorWH();
-		String x3d_path= correctionsParameters.selectX3dDirectory( // for x3d and obj
-				correctionsParameters.getModelName(image_name), // quad timestamp. Will be ignored if correctionsParameters.use_x3d_subdirs is false
-				correctionsParameters.x3dModelVersion,
-				true,  // smart,
-				true);  //newAllowed, // save
+		String x3d_path = getX3dDirectory();
 		String file_name = image_name + suffix;
 		String file_path = x3d_path + Prefs.getFileSeparator() + file_name + ".tiff";
 		ImagePlus imp = null;
@@ -1030,19 +1026,15 @@ public class QuadCLTCPU {
 			double [][] dsi
 			)
 	{
-		  String x3d_path= correctionsParameters.selectX3dDirectory( // for x3d and obj
-				  correctionsParameters.getModelName(image_name), // quad timestamp. Will be ignored if correctionsParameters.use_x3d_subdirs is false
-				  correctionsParameters.x3dModelVersion,
-					  true,  // smart,
-					  true);  //newAllowed, // save
-		  String title = image_name+TwoQuadCLT.DSI_COMBO_SUFFIX;
-		  ImagePlus imp = (new ShowDoubleFloatArrays()).makeArrays(dsi,tp.getTilesX(), tp.getTilesY(),  title, TwoQuadCLT.DSI_SLICES);
-			eyesisCorrections.saveAndShow(
-					   imp,      // ImagePlus             imp,
-					   x3d_path, // String                path,
-					   false,    // boolean               png,
-					   false,    // boolean               show,
-					   0);       // int                   jpegQuality)
+		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);
+		eyesisCorrections.saveAndShow(
+				imp,      // ImagePlus             imp,
+				x3d_path, // String                path,
+				false,    // boolean               png,
+				false,    // boolean               show,
+				0);       // int                   jpegQuality)
 	}
 
 	public void showDSI(){ showDSI(this.dsi);}
@@ -1056,11 +1048,7 @@ public class QuadCLTCPU {
 	public void saveDSIMain(
 			double [][] dsi) // DSI_SLICES.length
 	{
-		String x3d_path= correctionsParameters.selectX3dDirectory( // for x3d and obj
-				correctionsParameters.getModelName(image_name), // quad timestamp. Will be ignored if correctionsParameters.use_x3d_subdirs is false
-				correctionsParameters.x3dModelVersion,
-				true,  // smart,
-				true);  //newAllowed, // save
+		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]};
@@ -1078,11 +1066,7 @@ public class QuadCLTCPU {
 			String suffix, // "-DSI_MAIN"
 			double [][] dsi) // DSI_SLICES.length
 	{
-		String x3d_path= correctionsParameters.selectX3dDirectory( // for x3d and obj
-				correctionsParameters.getModelName(image_name), // quad timestamp. Will be ignored if correctionsParameters.use_x3d_subdirs is false
-				correctionsParameters.x3dModelVersion,
-				true,  // smart,
-				true);  //newAllowed, // save
+		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);
 		eyesisCorrections.saveAndShow(
@@ -1101,11 +1085,7 @@ public class QuadCLTCPU {
 			QuadCLT quadCLT_aux,
 			double [][] dsi_aux_from_main)
 	{
-		String x3d_path= correctionsParameters.selectX3dDirectory( // for x3d and obj
-				correctionsParameters.getModelName(image_name), // quad timestamp. Will be ignored if correctionsParameters.use_x3d_subdirs is false
-				correctionsParameters.x3dModelVersion,
-				true,  // smart,
-				true);  //newAllowed, // save
+		String x3d_path = getX3dDirectory();
 		String title = quadCLT_aux.image_name+"-DSI_GT-AUX";
 //		String []   titles =   {DSI_SLICES[DSI_DISPARITY_MAIN], DSI_SLICES[DSI_STRENGTH_MAIN]};
 //		double [][] dsi_main = {dsi[DSI_DISPARITY_MAIN],        dsi[DSI_STRENGTH_MAIN]};
@@ -5955,11 +5935,7 @@ public class QuadCLTCPU {
 			  }
 			  if (clt_parameters.gen_4_img) {
 				  // Save as individual JPEG images in the model directory
-				  String x3d_path= correctionsParameters.selectX3dDirectory(
-						  image_name, // quad timestamp. Will be ignored if correctionsParameters.use_x3d_subdirs is false
-						  correctionsParameters.x3dModelVersion,
-						  true,  // smart,
-						  true);  //newAllowed, // save
+					String x3d_path = getX3dDirectory();
 				  for (int sub_img = 0; sub_img < imps_RGB.length; sub_img++){
 					  EyesisCorrections.saveAndShow(
 							  imps_RGB[sub_img],
@@ -11047,12 +11023,7 @@ public class QuadCLTCPU {
 					  tp.clt_3d_passes.get(next_pass-1),   // CLTPass3d   scan,
 					  "after_pass3-"+(next_pass-1)); //String title)
 		  }
- 		  String x3d_path= correctionsParameters.selectX3dDirectory( // for x3d and obj
- 				 correctionsParameters.getModelName(this.image_name), // quad timestamp. Will be ignored if correctionsParameters.use_x3d_subdirs is false
- 				  correctionsParameters.x3dModelVersion,
-				  true,  // smart,
-				  true);  //newAllowed, // save
-
+			String x3d_path = getX3dDirectory();
  			// create x3d file
 		  if (clt_parameters.output_x3d) {
 			  x3dOutput = new X3dOutput(

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

@@ -24,6 +24,7 @@ package com.elphel.imagej.tileprocessor;
 import java.awt.Rectangle;
 import java.io.DataOutputStream;
 import java.io.File;
+import java.io.FileFilter;
 import java.io.FileInputStream;
 import java.io.FileNotFoundException;
 import java.io.FileOutputStream;
@@ -40,6 +41,7 @@ import java.nio.file.StandardCopyOption;
 import java.util.ArrayList;
 import java.util.Arrays;
 import java.util.Collections;
+import java.util.Comparator;
 import java.util.HashMap;
 import java.util.List;
 import java.util.Properties;
@@ -1034,11 +1036,14 @@ public class TwoQuadCLT {
 			}
 			if (clt_parameters.gen_4_img) {
 				// Save as individual JPEG images in the model directory
-				String x3d_path= quadCLT_main.correctionsParameters.selectX3dDirectory(
-						name, // quad timestamp. Will be ignored if correctionsParameters.use_x3d_subdirs is false
-						quadCLT_main.correctionsParameters.x3dModelVersion,
-						true,  // smart,
-						true);  //newAllowed, // save
+//				String x3d_path= quadCLT_main.correctionsParameters.selectX3dDirectory(
+//						name, // quad timestamp. Will be ignored if correctionsParameters.use_x3d_subdirs is false
+//						quadCLT_main.correctionsParameters.x3dModelVersion,
+//						true,  // smart,
+//						true);  //newAllowed, // save
+				String x3d_path = quadCLT_main.getX3dDirectory(name);
+				
+				
 				for (int sub_img = 0; sub_img < imps_RGB.length; sub_img++){
 					EyesisCorrections.saveAndShow(
 							imps_RGB[sub_img],
@@ -1927,11 +1932,15 @@ public class TwoQuadCLT {
 			}
 			if (clt_parameters.gen_4_img) {
 				// Save as individual JPEG images in the model directory
+				/*
 				String x3d_path= quadCLT_main.correctionsParameters.selectX3dDirectory(
 						name, // quad timestamp. Will be ignored if correctionsParameters.use_x3d_subdirs is false
 						quadCLT_main.correctionsParameters.x3dModelVersion,
 						true,  // smart,
 						true);  //newAllowed, // save
+						*/
+				String x3d_path = quadCLT_main.getX3dDirectory(name);
+
 				for (int sub_img = 0; sub_img < imps_RGB.length; sub_img++){
 					EyesisCorrections.saveAndShow(
 							imps_RGB[sub_img],
@@ -8686,6 +8695,7 @@ if (debugLevel > -100) return true; // temporarily !
 				clt_parameters,
 				colorProcParameters, //
 				noise_sigma_level,   // double []            noise_sigma_level,
+				-1,                  // int                  noise_variant, // <0 - no-variants, compatible with old code
 				null,                // QuadCLTCPU           ref_scene, // may be null if scale_fpn <= 0
 				threadsMax,
 				clt_parameters.inp.noise_debug_level); // debugLevel);
@@ -9429,8 +9439,6 @@ if (debugLevel > -100) return true; // temporarily !
 				"-results-rnd_2.5-fpn_0.0-sigma_1.5-offset1.4142-sensors8-inter-nolma",
 				"-results-rnd_2.5-fpn_0.0-sigma_1.5-offset1.4142-sensors8-nointer-nolma",
 				
-				
-				
 				/*
 				"-results-rnd_0.0-fpn_0.0-sigma_1.5-offset1.0-sensors16-inter",
 				"-results-rnd_0.0-fpn_0.0-sigma_1.5-offset1.0-sensors16-nointer",
@@ -9628,11 +9636,53 @@ if (debugLevel > -100) return true; // temporarily !
 				"-results-lev_5.0-sigma_1.5-offset1.0-nointer-mask1"
 				*/
 				};
-		
+		// extend files by noise variants
+		//ref_scene
+
+		ArrayList<String> full_files_list = new ArrayList<String>();
+		String x3d_path = ref_scene.getX3dDirectory()+"";
+		File model_directory = new File(x3d_path);
+		int [] group_indices = new int [noise_files.length + 1];
+		for (int nf = 0; nf < noise_files.length; nf++) {
+			group_indices[nf] = full_files_list.size();
+			String base_suffix = noise_files[nf];
+			final String file_prefix = ref_scene.getImageName()+base_suffix;
+			FileFilter noiseFilefilter = new FileFilter() 
+			{
+				//Override accept method
+				public boolean accept(File file) {
+
+					//if the file extension is .log return true, else false
+					if (file.getName().endsWith(".tiff") && file.getName().startsWith(file_prefix)) {
+						return true;
+					}
+					return false;
+				}
+			};
+			
+			File[] files = model_directory.listFiles(noiseFilefilter);
+			ArrayList<String> flist = new ArrayList<String>();
+			for (File f:files) {
+				String name = f.getName();
+				String suffix = name.substring(ref_scene.getImageName().length(), name.length() - ".tiff".length());
+				flist.add(suffix);
+			}
+
+    		Collections.sort(flist, new Comparator<String>() {
+    		    @Override
+    		    public int compare(String lhs, String rhs) { // ascending
+    		        return rhs.length() > lhs.length()  ? -1 : (rhs.length() < lhs.length()) ? 1 : lhs.compareTo(rhs);
+    		    }
+    		});
+    		full_files_list.addAll(flist);
+		}
+		group_indices[group_indices.length - 1] = full_files_list.size();
+		String [] noise_files_full = full_files_list.toArray(new String[0]);
 		getNoiseStats(
 				clt_parameters,
 				ref_scene, // ordered by increasing timestamps
-				noise_files,
+				noise_files_full, // noise_files,
+				group_indices,
 				debug_level);		
 	}
 
@@ -9739,38 +9789,68 @@ if (debugLevel > -100) return true; // temporarily !
 		return good_tiles;
 	}
 	
-	
+	class DisparityResults{
+		int    []   num_instances;
+		double [][] results;
+	}
+	class NoiseLevel implements Comparable<NoiseLevel>{
+		double rnd;
+		double fpn;
+		NoiseLevel(double rnd, double fpn) {
+			this.rnd = rnd;
+			this.fpn = fpn;
+		}
+		@Override
+		public boolean equals(Object obj) {
+			if ((obj == null) || (getClass() != obj.getClass())){
+				return false;
+			}
+			NoiseLevel other = (NoiseLevel) obj;
+			return (other.rnd == rnd) && (other.fpn == fpn);
+		}
+		@Override
+		public int compareTo(NoiseLevel other) {
+			return (this.fpn > other.fpn)? 1 :((this.fpn < other.fpn) ? -1: ((this.rnd > other.rnd)? 1 : ((this.rnd < other.rnd) ? -1: 0)));
+		}
+        @Override
+        public int hashCode() {
+        	return ((Double)(rnd + 7919*fpn)).hashCode();
+        }
+		
+	}
 	
 	public void getNoiseStats(
 			CLTParameters        clt_parameters,
 			QuadCLT              ref_scene, // ordered by increasing timestamps
 			String []            noise_files,
+			int []               group_indices,
 			int                  debug_level)
 	{
-		int dbg_tile = 829; // 828; // 1222; // 737;
-
-		/*
-		"disp-last",
-		"str_last",
-		"num vlaid" <= 1.0
-		*/
+		int dbg_tile = 194; // 829; // 828; // 1222; // 737;
 		final int tilesX = ref_scene.tp.getTilesX();
-//		final int tilesY = ref_scene.tp.getTilesY();
+		int [] var_map = new int [noise_files.length]; // for each file - group number (matching group_indices).
+		double [] var_weights = new double [noise_files.length];
+		for (int i = 0; i < group_indices.length-1; i++) {
+			int indx_from = group_indices[i];
+//			int indx_to =   (i < (group_indices.length -1))? group_indices[i+1] : noise_files.length ;
+			int indx_to =   group_indices[i+1];
+			double w = 1.0/(indx_to - indx_from);
+			for (int j = indx_from; j < indx_to; j++) {
+				var_map[j] = i;
+				var_weights[j] = w;
+			}
+		}
 
 		double max_diff = 0.01; // 0.001; // 0.01; //  0.04; // 0.01; // last diff >
-		double max_err = 2.0; // 2.5; // 1.5; // 2.0; // 1.0; // 0.5; // pix
-		double max_err1 =0.250; // pix
+		double max_err =  2.0; // 2.5; // 1.5; // 2.0; // 1.0; // 0.5; // pix
+		double max_err1 = 0.250; // pix
 		double min_strength = 0.0; // minimal strength to calculate rmse (ignore weaker)
 		int indx_used =      3;
 		int indx_last =      0;
 		int indx_lma_last =  clt_parameters.correlate_lma? 2 : 0; // if lma was disabled fallback to just disparity
 		int indx_diff_last = 4;
-		
-//		int indx_initial =  3;
-		
 		int indx_strength = 1;
 		double max_disparity = 30.0; // for max_err1
-//		double disp_rel_min = 0.5;
 		
 		double disp_near_rel =     2.5;
 		double disp_max_rel =      0.25;
@@ -9785,31 +9865,35 @@ if (debugLevel > -100) return true; // temporarily !
 		int          min_modes =   4; // 5; // 6; // 5; // 4;//at least half are meaningfull
 		
 		// LMA parameters
-		boolean useLinear =         true;
-		double  noise_offset =      0.05; // 0.1; // 0.03; // 0.10; // 0.03; // 50;
+		boolean useLinear =         false; // true; // false; // true;
+		double  noise_offset =      0.03; // 0.05; // 0.03; // 0.5; // 0.1; // 0.03; //  0.1; // 0.05; // 0.1; // 0.03; // 0.10; // 0.03; // 50;
 		double  n0 =                0.03;
-		boolean adjust_N0 =         true;
+		boolean adjust_N0 =         false; // true;
 		boolean adjust_Gi =         true;
-		boolean adjust_St =         true; // false;
-		
-		double min_inter16_noise_level = 0.10; // 0.3; // tile should have at least this noise level for 1nter16 (mode 0)
+		boolean adjust_St =         true; // false; // true; // false;
 		
+		// change to only apply to intra, inter ahould use weak.
+		double min_inter16_noise_level = 0.5; // 0.3; // 0.1; // 0.3; //  0.1; // 0.3; // tile should have at least this noise level for 1nter16 (mode 0)
+		boolean apply_min_inter16_to_inter = false;
 		boolean zero_all_bad =      true; // false; //  true;    // set noise_level to zero if all noise levels result in bad tiles
 		boolean all_inter =         true;    // tile has to be defined for all inter
 		boolean need_same_inter =   true; // do not use intra sample if same inter is bad for all noise levels  
+		boolean need_same_zero =    false; // true; // do not use sample if it is bad for zero-noise  
 
-		double max_diff_from_ref = 0.20; // 0.06; // 5; // 0.1; // max_err1; // 0.25 pix
+		double max_diff_from_ref = 0.25; // 0.1; // 0.20; // 0.06; // 5; // 0.1; // max_err1; // 0.25 pix
 		boolean use_fpn = false;
 
 		double max_diff_from_ref_range = 0.25*max_diff_from_ref; // trying to stay in linear
 		int    max_diff_from_ref_steps = 21;
 		int    range_outliers =          2;
-		int    range_min_keep =          1; // emove less outliers if needed to keep this remain 
-
-		
+		int    range_min_keep =          1; // remove less outliers if needed to keep this remain 
+		int    num_var_outliers =        4;
+		int    min_var_remain =         10;
+		double scale_intra =             0.2;
+		boolean run_lma =           false;
 		
 		if (use_edges) {
-			disp_max_rel = 100.0;
+			disp_max_rel = 100.00;
 			disp_max_abs = 100.0;
 		}
 		
@@ -9819,7 +9903,6 @@ if (debugLevel > -100) return true; // temporarily !
 		null); // int []      wh);
 		
 		double [][] ref_dsn = ref_scene.readDoubleArrayFromModelDirectory(
-//				"-results-nonoise", // String      suffix,
 				"-results-nonoise-nolma", // String      suffix,
 				0, // int         num_slices, // (0 - all)
 				null); // int []      wh);
@@ -9845,10 +9928,6 @@ if (debugLevel > -100) return true; // temporarily !
 		
 		int num_good_init = 0;
 		for (int i = 0; i < good_tiles_ref.length; i++) {
-//			good_tiles[i] = (ref_dsn[indx_used][i] > 0.999) && (Math.abs(ref_dsn[indx_last_diff][i]) < max_diff);
-//			if (good_tiles[i] && (sky_map != null) && (sky_map[0][i] > 0.0)) {
-//				good_tiles[i] = false;
-//			}
 			if (good_tiles_ref[i]){
 				num_good_init++;
 			}
@@ -9859,48 +9938,17 @@ if (debugLevel > -100) return true; // temporarily !
 		}
 		
 		
-//		double  [] noise_level = new double  [noise_files.length];
-//		boolean [] intra =       new boolean [noise_files.length];
-//		boolean [] inter =       new boolean [noise_files.length];
-		class DisparityResults{
-			double [][] results;
-		}
-		class NoiseLevel implements Comparable<NoiseLevel>{
-			double rnd;
-			double fpn;
-			NoiseLevel(double rnd, double fpn) {
-				this.rnd = rnd;
-				this.fpn = fpn;
-			}
-			@Override
-			public boolean equals(Object obj) {
-				if ((obj == null) || (getClass() != obj.getClass())){
-					return false;
-				}
-				NoiseLevel other = (NoiseLevel) obj;
-				return (other.rnd == rnd) && (other.fpn == fpn);
-			}
-			@Override
-			public int compareTo(NoiseLevel other) {
-				return (this.fpn > other.fpn)? 1 :((this.fpn < other.fpn) ? -1: ((this.rnd > other.rnd)? 1 : ((this.rnd < other.rnd) ? -1: 0)));
-			}
-	        @Override
-	        public int hashCode() {
-	        	return ((Double)(rnd + 7919*fpn)).hashCode();
-	        }
-			
-		}
-		
 //		boolean all_converge = false; // use only tiles that converge for all variants (intra, inter, used sensors)
 //		boolean all_max_err =  false;  // use only tiles that have limited error for all variants (intra, inter, used sensors)
 		boolean [] converged_tiles = good_tiles_ref.clone();
 		boolean [] good_tiles =      good_tiles_ref.clone();
 		boolean [][] good_tiles_mode = {good_tiles_ref.clone(),good_tiles_ref.clone(),good_tiles_ref.clone(),good_tiles_ref.clone()}; // per intra mode
 		double [][][] ref_dsn_mode = new double[4][][];
-		int [] sensor_mode_file = new int [noise_files.length];
+		int [] sensor_mode_file =  new int [noise_files.length];
 		double [] noise_rnd_file = new double [noise_files.length];
 		double [] noise_fpn_file = new double [noise_files.length];
 		boolean [] inter_file =    new boolean [noise_files.length];
+//		int [] var_map = new int [noise_files.length]; // for each file - group number (matching group_indices).
 		
 		// extract data from file names
 		for (int nf = 0; nf < noise_files.length; nf++) {
@@ -9935,10 +9983,9 @@ if (debugLevel > -100) return true; // temporarily !
 		}		
 		
 		
-		if (all_converge || all_max_err || same_num_sensors) { // scan all images
+		if (all_converge || all_max_err || same_num_sensors) { // scan all images, use only interscene with no-noise
 			for (int nf = 0; nf < noise_files.length; nf++) {
 				String fn = noise_files[nf];
-				
 				if (inter_file[nf]) {
 					double [][] noise_dsn = ref_scene.readDoubleArrayFromModelDirectory(
 							fn, // noise_files[nf], // String      suffix,
@@ -10008,7 +10055,6 @@ if (debugLevel > -100) return true; // temporarily !
 		// For each file find boolean good/bad, comparing to zero noise of the same number of sensors, interscene 
 		boolean [][] good_file_tile = new boolean[noise_files.length][]; // [good_tiles.length];
 		boolean [][][] good_file_tile_range = new boolean[noise_files.length][][]; // [good_tiles.length];
-//		double max_err1 =0.25; // pix
 		for (int nf = 0; nf < noise_files.length; nf++) {
 			// common or per number of sensors reference data 
 			String fn = noise_files[nf];
@@ -10018,7 +10064,6 @@ if (debugLevel > -100) return true; // temporarily !
 					fn, // noise_files[nf], // String      suffix,
 					0, // int         num_slices, // (0 - all)
 					null); // int []      wh);
-//			boolean [] good_ref = good_tiles_mode[sensor_mode]; // good tile without noise for this number of sensors
 			good_file_tile[nf] = good_tiles_mode[sensor_mode].clone();
 			good_file_tile_range[nf] = new boolean [good_tiles_mode[sensor_mode].length][];
 			for (int ntile = 0; ntile < good_file_tile[nf].length; ntile++) if (good_file_tile[nf][ntile]) {
@@ -10026,7 +10071,6 @@ if (debugLevel > -100) return true; // temporarily !
 					System.out.println("Finding good tiles: ntile = "+ntile+", nf="+nf+" ("+fn+")");
 					System.out.println("noise_dsn["+indx_last+"]["+ntile+"]="+noise_dsn[indx_last][ntile]+
 							", ref_var["+indx_last+"]["+ntile+"]="+ref_var[indx_last][ntile]);
-					
 				}
 				
 				boolean converged = (Math.abs(noise_dsn[indx_last_diff][ntile]) < max_diff);
@@ -10046,16 +10090,9 @@ if (debugLevel > -100) return true; // temporarily !
 				if (!has_good) {
 					good_file_tile_range[nf][ntile] = null; // all bad
 				}
-				
-				/*
-				boolean good_tiles_this = (Math.abs(noise_dsn[indx_last][ntile] - ref_var[indx_last][ntile]) < max_diff_from_ref);
-				if (!good_tiles_this) {
-					good_file_tile[nf][ntile] = false;
-					continue;
-				}
-				*/
 			}			
 		}
+		
 		{
 			// show number of noise values for each tile, num sensors and intra/inter, discarding tiles that are good/bad for all noise levels
 			double [][] dbg_num_noise_val = new double [good_tiles_mode.length*2][good_tiles.length];
@@ -10073,12 +10110,14 @@ if (debugLevel > -100) return true; // temporarily !
 				}
 
 				if (good_tiles[ntile]) { // do not bother with obviously bad
+					double [] val_good = new double [dbg_num_noise_val.length];
 					int [] num_good =    new int [dbg_num_noise_val.length];
 					boolean [] has_bad = new boolean [dbg_num_noise_val.length];
 					for (int nf = 0; nf < noise_files.length; nf++) {
 						int results_index = sensor_mode_file[nf] + (inter_file[nf]? 0 : 4); // inter; //  ? 0 : (intra? 1 : 2);
 						if (good_file_tile[nf][ntile]) {
 							num_good[results_index]++;
+							val_good[results_index] += var_weights[nf];
 						} else {
 							has_bad[results_index] = true;
 						}
@@ -10094,7 +10133,7 @@ if (debugLevel > -100) return true; // temporarily !
 						continue;
 					}
 					for (int i = 0; i < dbg_num_noise_val.length; i++) if (has_bad[i]) {
-						dbg_num_noise_val[i][ntile] = num_good[i];
+						dbg_num_noise_val[i][ntile] = val_good[i]; // num_good[i];
 					}
 				}
 			}
@@ -10106,7 +10145,7 @@ if (debugLevel > -100) return true; // temporarily !
 					"num_noise_levels",
 					dbg_num_noise_titles);
 
-			for (int i = 00; i < dbg_num_noise_val.length; i++) {
+			for (int i = 0; i < dbg_num_noise_val.length; i++) {
 				Arrays.fill(dbg_num_noise_val[i], Double.NaN);
 			}
 			for (int ntile = 0; ntile < good_tiles.length; ntile++) {
@@ -10116,6 +10155,7 @@ if (debugLevel > -100) return true; // temporarily !
 
 				if (good_tiles[ntile]) { // do not bother with obviously bad
 					int [] num_good =    new int [dbg_num_noise_val.length];
+					double [] val_good = new double [dbg_num_noise_val.length];
 					boolean [] has_bad = new boolean [dbg_num_noise_val.length];
 					for (int nf = 0; nf < noise_files.length; nf++) {
 						int results_index = sensor_mode_file[nf] + (inter_file[nf]? 0 : 4); // inter; //  ? 0 : (intra? 1 : 2);
@@ -10123,6 +10163,7 @@ if (debugLevel > -100) return true; // temporarily !
 							for (int stp = 0; stp < good_file_tile_range[nf][ntile].length; stp++) {
 								if (good_file_tile_range[nf][ntile][stp]) {
 									num_good[results_index]++;
+									val_good[results_index] += var_weights[nf];
 								} else {
 									has_bad[results_index] = true;
 								}
@@ -10143,7 +10184,7 @@ if (debugLevel > -100) return true; // temporarily !
 						continue;
 					}
 					for (int i = 0; i < dbg_num_noise_val.length; i++) if (has_bad[i]) {
-						dbg_num_noise_val[i][ntile] = num_good[i];
+						dbg_num_noise_val[i][ntile] = val_good[i]; // num_good[i];
 					}
 				}
 			}
@@ -10157,18 +10198,28 @@ if (debugLevel > -100) return true; // temporarily !
 			
 			
 			double []    noise_file = use_fpn ? noise_fpn_file : noise_rnd_file;
-			double [][] noise_levels0 = 	InterIntraLMA.getNoiseThreshold(
+			double [][][] noise_levels_partial0 = 	InterIntraLMA.getNoiseThresholdsPartial(
 					noise_file,              // double []    noise_file, //  = new double [noise_files.length];
+					group_indices,           // int []       group_indices,  // last points after last file index
 					sensor_mode_file,        // int []       sensor_mode_file,
 					inter_file,              // boolean []   inter_file,
 					good_file_tile,          // boolean [][] good_file_tile,
 					min_inter16_noise_level, // double       min_inter16_noise_level,
+					apply_min_inter16_to_inter, // boolean      apply_min_inter16_to_inter,
 					min_modes,               // int          min_modes,
 					zero_all_bad,            // boolean zero_all_bad = true;    // set noise_level to zero if all noise levels result in bad tiles
 					all_inter,               // boolean all_inter =    true;    // tile has to be defined for all inter
 					need_same_inter,         // boolean need_same_inter = true; // do not use intra sample if same inter is bad for all noise levels
+					need_same_zero,          // boolean need_same_zero,        // do not use samle if it is bad for zero-noise 
 					dbg_tile);               // int            dbg_tile);
 					
+			double [][] noise_levels0 = InterIntraLMA.mergeNoiseVariants(
+					noise_levels_partial0, // double [][][] partial_thresholds,
+					num_var_outliers,      // int           num_outliers,
+					min_var_remain);       // int           min_remain);
+					
+					
+					
 			double [][] dbg_noise_levels = new double [dbg_num_noise_titles.length][good_tiles.length];
 			for (int i = 0; i < dbg_noise_levels.length; i++) {
 				Arrays.fill(dbg_noise_levels[i], Double.NaN);
@@ -10187,7 +10238,7 @@ if (debugLevel > -100) return true; // temporarily !
 					dbg_num_noise_titles);
 			{
 //				int dbg_tile = 828;
-				for (int mode = 0; mode < 8; mode++) {
+				for (int mode = 00; mode < 8; mode++) {
 					for (int nf = 0; nf < noise_files.length; nf++) if(good_file_tile_range[nf] !=null){ // always
 						int mode_file = sensor_mode_file[nf] + (inter_file[nf]? 0 : 4); // inter; //  ? 0 : (intra? 1 : 2);
 						if (mode_file == mode) {
@@ -10200,27 +10251,32 @@ if (debugLevel > -100) return true; // temporarily !
 								System.out.println(String.format("%1d:%3d %6.4f %s", mode, nf, noise_rnd, s));
 							} else {
 								System.out.println(String.format("%1d:%3d %6.4f", mode, nf, noise_rnd));
-
 							}
 						}
-
 					}
 				}
 			}
 			
-			double [][] noise_levels = 	InterIntraLMA.getNoiseThreshold(
+			double [][][] noise_levels_partial = 	InterIntraLMA.getNoiseThresholdsPartial(
 					noise_file,              // double []    noise_file, //  = new double [noise_files.length];
+					group_indices,           // int []       group_indices,  // last points after last file index
 					sensor_mode_file,        // int []       sensor_mode_file,
 					inter_file,              // boolean []   inter_file,
 					range_outliers,          // int            outliers,  // may need do modify algorithm to avoid bias - removing same side (densier) outliers 
 					range_min_keep,          // int            min_keep, // remove less outliers if needed to keep this remain 
 					good_file_tile_range,    // boolean [][][] good_file_tile_range,
 					min_inter16_noise_level, // double       min_inter16_noise_level,
-					min_modes+0,             // int          min_modes,
+					apply_min_inter16_to_inter, // boolean      apply_min_inter16_to_inter,
+					min_modes,               // int          min_modes,
 					zero_all_bad,            // boolean zero_all_bad = true;    // set noise_level to zero if all noise levels result in bad tiles
 					all_inter,               // boolean all_inter =    true;    // tile has to be defined for all inter
 					need_same_inter,         // boolean need_same_inter = true; // do not use intra sample if same inter is bad for all noise levels
+					need_same_zero,          // boolean need_same_zero,        // do not use samle if it is bad for zero-noise 
 					dbg_tile);               // int            dbg_tile);
+			double [][] noise_levels = InterIntraLMA.mergeNoiseVariants(
+					noise_levels_partial, // double [][][] partial_thresholds,
+					num_var_outliers,      // int           num_outliers,
+					min_var_remain);       // int           min_remain);
 			
 			double [][] dbg_noise_levels_range = new double [dbg_num_noise_titles.length][good_tiles.length];
 			for (int i = 0; i < dbg_noise_levels_range.length; i++) {
@@ -10247,8 +10303,9 @@ if (debugLevel > -100) return true; // temporarily !
 					noise_offset, // double      offset  // for "relative" noise
 					n0 ,          // double       n0,    // initial value for N0
 					tilesX,       // int         tilesX, // debug images only
+					scale_intra, // double      scale_intra, // scale weight of intra-scene samples ( < 1.0)
 					1); // int         debug_level)
-
+			if (run_lma) {
 
 			boolean LMA_OK = interIntraLMA.runLma(
 					adjust_N0, // boolean adjust_N0,
@@ -10257,7 +10314,7 @@ if (debugLevel > -100) return true; // temporarily !
 					0.1, // double lambda,           // 0.1
 					0.5, // double lambda_scale_good,// 0.5
 					8.0, // double lambda_scale_bad, // 8.0
-					100, // double lambda_max,       // 100
+					1000, // double lambda_max,       // 100
 					0.001, // double rms_diff,         // 0.001
 					30, // int    num_iter,         // 20
 					2); // 0); // int    debug_level)
@@ -10279,7 +10336,7 @@ if (debugLevel > -100) return true; // temporarily !
 					tilesX,
 					good_tiles.length/ tilesX,
 					"lma_st_out");
-			
+			}
 
 		}
 
@@ -10352,23 +10409,28 @@ if (debugLevel > -100) return true; // temporarily !
 			if (!results_map.containsKey(noise_level)) {
 				DisparityResults dr = new DisparityResults();
 				dr.results = new double [8][];
+				dr.num_instances = new int[8];
 				results_map.put(noise_level, dr);
 			}
 			// inter 16, inter 8, inter 4, inter 2, intra16, intra 8, intra 4, intra 2 
 			
 			DisparityResults dr = results_map.get(noise_level);
-			dr.results[results_index] = results;
+			if (dr.results[results_index] == null) {
+				dr.results[results_index] = results.clone();
+			} else {
+				for (int i = 0; i < results.length; i++) {
+					dr.results[results_index][i] += results[i];
+				}
+			}
+			dr.num_instances[results_index]++;
 			System.out.println("getNoiseStats(): "+noise_files[nf]+": good_ref= " + num_good_init+
 					", converged= "+num_converged+", good= "+num_good+" good(0.1)= "+num_good1+", num_near="+num_near);
 		}
-//		List<Double> noise_levels_list = new ArrayList<Double>(results_map.keySet());
 		List<NoiseLevel> noise_levels_list = new ArrayList<NoiseLevel>(results_map.keySet());
 		Collections.sort(noise_levels_list);
 		
-//		String [] config_types = {"all_16", "octal", "quad", "binocular"};
 		String [] config_types = {"16", "8", "4", "2"};
 		System.out.println("\n");
-//		System.out.print("noise_level, ");
 		System.out.print("noise_random, noise_fpn, ");
 		
 		for (int it = 0; it < config_types.length; it++) {
@@ -10383,9 +10445,17 @@ if (debugLevel > -100) return true; // temporarily !
 					"("+max_err1+"), intra_rmse_"+config_types[it]+"("+max_err+"),");
 		}
 		
-//		System.out.print("inter("+max_err+"), inter("+max_err1+"), inter_conf("+max_err+"), inter_conf("+max_err1+"), inter_rmse("+max_err+"),");
-//		System.out.print("intra("+max_err+"), intra("+max_err1+"), intra_conf("+max_err+"), intra_conf("+max_err1+"), intra_rmse("+max_err+"),");
-//		System.out.print("binocular("+max_err+"), binocular("+max_err1+"), binocular_conf("+max_err+"), binocular_conf("+max_err1+"), binocular_rmse("+max_err+")");
+		for (NoiseLevel nl:noise_levels_list) {
+			double [][] results = results_map.get(nl).results;
+			int [] num_instances = results_map.get(nl).num_instances;
+			for (int m = 0; m < results.length; m++) {
+				if (num_instances[m] > 1) {
+					for (int i = 0; i < results[m].length; i++){
+						 results[m][i] /= num_instances[m];
+					}
+				}
+			}
+		}
 		System.out.println();
 		for (NoiseLevel nl:noise_levels_list) {
 			System.out.print(nl.rnd+", "+nl.fpn+", ");
@@ -10409,10 +10479,382 @@ if (debugLevel > -100) return true; // temporarily !
 			}
 			System.out.println();
 		}
+		// good (2.0), good(.25) conf(2.0), conf(0.25) rmse(2.0)
+		double [] rslt_err= {max_err1,max_err,max_err1,max_err,max_err}; 
+		int [] used_results = {0,1,4};
+		int [][] sens_to_res = {{7,6,5,4},{3,2,1,0}}; //[inter]{2,4,8,16};
 		System.out.println();
+		double [] noise_lev = new double [noise_levels_list.size()];
+		double [][][][] plots_direct = new double [used_results.length][sens_to_res.length][sens_to_res[0].length][noise_levels_list.size()];
+		int nn = 0;
+		for (NoiseLevel nl:noise_levels_list) {
+			noise_lev[nn] = use_fpn ? nl.fpn : nl.rnd;
+			double [][] results = results_map.get(nl).results;
+			for (int pt = 0; pt < plots_direct.length; pt++) {
+				for (int inter = 0; inter < sens_to_res.length; inter++) {
+					for (int isens = 0; isens < sens_to_res[0].length; isens++) {
+						plots_direct[pt][inter][isens][nn] = results[sens_to_res[inter][isens]][used_results[pt]]; 
+					}
+				}
+			}
+			nn++;
+		}
+		int num_y_steps = 100;
+		// compare RMSE : 4, 8, 16 to binocular
+		double rmse_min =         0.0;
+		double rmse_max =         0.9;
+		double rmse_max_ratio =   8.0;
+		double rmse_max_ratio1 =  4.0;
+		double rmse_max_ratio2 = 15.0;
+		double density_min     =  0.0;
+		double density_max     =  1.0;
+		plotInverted(
+				noise_levels_list, // List<NoiseLevel> noise_levels_list,
+				results_map,       // HashMap <NoiseLevel, DisparityResults> results_map,
+				max_err,           // double max_err,
+				max_err1,          // double max_err1,
+				use_fpn,           // boolean use_fpn,
+				num_y_steps,       // int num_y_steps, //  = 100;
+				rmse_min,          // double rmse_min, //  = 0.0;
+				rmse_max,          // double rmse_max, //  = 0.9;
+				rmse_max_ratio,    // double rmse_max_ratio, // = 6.0;
+				rmse_max_ratio1,   // double rmse_max_ratio1,// = 4.0;
+				rmse_max_ratio2,   // double rmse_max_ratio2 // = 20.0;
+				density_min,       // double density_min,    // == 0
+				density_max);       // double density_max,    // == 1
+
+	
 	}
+
+	public void plotInverted(
+			List<NoiseLevel> noise_levels_list,
+			HashMap <NoiseLevel, DisparityResults> results_map,
+			double max_err,
+			double max_err1,
+			boolean use_fpn,
+			int num_y_steps, //  = 100;
+			double rmse_min, //  = 0.0;
+			double rmse_max, //  = 0.9;
+			double rmse_max_ratio, // = 6.0;
+			double rmse_max_ratio1,// = 4.0;
+			double rmse_max_ratio2,// = 20.0;
+			double density_min,    // == 0
+			double density_max     // == 1
+			) {
+		// good (2.0), good(.25) conf(2.0), conf(0.25) rmse(2.0)
+		double [] rslt_err= {max_err1,max_err,max_err1,max_err,max_err}; 
+		int [] used_results = {0,1,4};
+		int [][] sens_to_res = {{7,6,5,4},{3,2,1,0}}; //[inter]{2,4,8,16};
+		System.out.println();
+		double [] noise_lev = new double [noise_levels_list.size()];
+		double [][][][] plots_direct = new double [used_results.length][sens_to_res.length][sens_to_res[0].length][noise_levels_list.size()];
+		int nn = 0;
+		for (NoiseLevel nl:noise_levels_list) {
+			noise_lev[nn] = use_fpn ? nl.fpn : nl.rnd;
+			double [][] results = results_map.get(nl).results;
+			for (int pt = 0; pt < plots_direct.length; pt++) {
+				for (int inter = 0; inter < sens_to_res.length; inter++) {
+					for (int isens = 0; isens < sens_to_res[0].length; isens++) {
+						plots_direct[pt][inter][isens][nn] = results[sens_to_res[inter][isens]][used_results[pt]]; 
+					}
+				}
+			}
+			nn++;
+		}
+		
+		
+		
+		// compare RMSE : 4, 8, 16 to binocular
+		int pt_rmse = 2;
+		String [] col_titles_rmse = {
+				"rmse("+rslt_err[pt_rmse]+")",
+				"4:2 intra",
+				"8:2 intra",
+				"16:2 intra",
+				"4:2 inter",
+				"8:2 inter",
+				"16:2 inter"};
+		for (int i = 0; i < col_titles_rmse.length; i++) {
+			System.out.print(col_titles_rmse[i]);
+			if (i < (col_titles_rmse.length -1)) {
+				System.out.print(", ");
+			}
+			
+		}
+		System.out.println();
+		double [] rmse_vals = inverseLinTFunc(
+				noise_lev, // double [] x_val,
+				plots_direct[pt_rmse][0][0], // double [] y_val,
+				rmse_min, // double    y_min,
+				rmse_max, // double    y_max,
+				num_y_steps+1)[0]; //int       npoints		) 
+		double [][][] rmse_rslt = new double [sens_to_res.length][sens_to_res[0].length][];
+		for (int inter = 0; inter < sens_to_res.length; inter++) {
+			for (int isens = 0; isens < sens_to_res[0].length; isens++) {
+				rmse_rslt[inter][isens] = inverseLinTFunc(
+						noise_lev, // double [] x_val,
+						plots_direct[pt_rmse][inter][isens], // double [] y_val,
+						rmse_min, // double    y_min,
+						rmse_max, // double    y_max,
+						num_y_steps+1)[1]; //int       npoints		) 
+			}
+		}
+		for (int i = 0; i < rmse_vals.length; i++) {
+			System.out.print(String.format("%8.5f, ", rmse_vals[i]));
+			for (int inter = 0; inter < sens_to_res.length; inter++) {
+				for (int isens = 1; isens < sens_to_res[0].length; isens++) {
+					double ratio = rmse_rslt[inter][isens][i]/rmse_rslt[inter][0][i];
+					if ((ratio > rmse_max_ratio) || Double.isInfinite(ratio)) {
+						ratio = Double.NaN;
+					}
+					if (!Double.isNaN(ratio)) {
+						System.out.print (String.format("%8.5f", ratio));
+					}
+					if (isens < (sens_to_res[0].length - 1)) {
+						System.out.print (", ");
+					}
+				}
+				if (inter < (sens_to_res.length - 1)) {
+					System.out.print (", ");
+				}
+			}
+			System.out.println();
+		}
+		System.out.println("\n");
+		// compare RMSE : 4 to binocular, 8 to 4, 16 to 4
+		String [] col_titles_rmse1 = {
+				"rmse("+rslt_err[pt_rmse]+")",
+				"4:2 intra",
+				"8:4 intra",
+				"16:8 intra",
+				"4:2 inter",
+				"8:4 inter",
+				"16:8 inter"};
+		for (int i = 0; i < col_titles_rmse1.length; i++) {
+			System.out.print(col_titles_rmse1[i]);
+			if (i < (col_titles_rmse1.length -1)) {
+				System.out.print(", ");
+			}
+		}
+		System.out.println();
+		
+		for (int i = 0; i < rmse_vals.length; i++) {
+			System.out.print(String.format("%8.5f, ", rmse_vals[i]));
+			for (int inter = 0; inter < sens_to_res.length; inter++) {
+				for (int isens = 1; isens < sens_to_res[0].length; isens++) {
+					double ratio = rmse_rslt[inter][isens][i]/rmse_rslt[inter][isens - 1][i];
+					if ((ratio > rmse_max_ratio1) || Double.isInfinite(ratio)) {
+						ratio = Double.NaN;
+					}
+					if (!Double.isNaN(ratio)) {
+						System.out.print (String.format("%8.5f", ratio));
+					}
+					if (isens < (sens_to_res[0].length - 1)) {
+						System.out.print (", ");
+					}
+				}
+				if (inter < (sens_to_res.length - 1)) {
+					System.out.print (", ");
+				}
+			}
+			System.out.println();
+		}
+		System.out.println("\n");
+		// compare RMSE : inter to same intra
+		String [] col_titles_rmse2 = {
+				"rmse("+rslt_err[pt_rmse]+")",
+				"inter:intra (2)",
+				"inter:intra (4)",
+				"inter:intra (8)",
+				"inter:intra (16)"};
+		for (int i = 0; i < col_titles_rmse2.length; i++) {
+			System.out.print(col_titles_rmse2[i]);
+			if (i < (col_titles_rmse2.length -1)) {
+				System.out.print(", ");
+			}
+		}
+		System.out.println();
+		
+		for (int i = 0; i < rmse_vals.length; i++) {
+			System.out.print(String.format("%8.5f, ", rmse_vals[i]));
+			for (int isens = 0; isens < sens_to_res[0].length; isens++) {
+				double ratio = rmse_rslt[1][isens][i]/rmse_rslt[0][isens][i];
+				if ((ratio > rmse_max_ratio2) || Double.isInfinite(ratio)) {
+					ratio = Double.NaN;
+				}
+				if (!Double.isNaN(ratio)) {
+					System.out.print (String.format("%8.5f", ratio));
+				}
+				if (isens < (sens_to_res[0].length - 1)) {
+					System.out.print (", ");
+				}
+			}
+			System.out.println();
+		}
+		System.out.println("\n");
+		
+
+		for (int idensity_err = 0; idensity_err < 2; idensity_err++) {
+			
+			// compare density (2.0 and 0.25) : 4, 8, 16 to binocular
+			int pt_density = used_results[idensity_err];
+			String [] col_titles_density = {
+					"density("+rslt_err[pt_density]+")",
+					"4:2 intra",
+					"8:2 intra",
+					"16:2 intra",
+					"4:2 inter",
+					"8:2 inter",
+					"16:2 inter"};
+			for (int i = 0; i < col_titles_density.length; i++) {
+				System.out.print(col_titles_density[i]);
+				if (i < (col_titles_density.length -1)) {
+					System.out.print(", ");
+				}
+				
+			}
+			System.out.println();
+			double [] density_vals = inverseLinTFunc(
+					noise_lev, // double [] x_val,
+					plots_direct[pt_density][0][0], // double [] y_val,
+					rmse_min, // double    y_min,
+					rmse_max, // double    y_max,
+					num_y_steps+1)[0]; //int       npoints		) 
+			double [][][] density_rslt = new double [sens_to_res.length][sens_to_res[0].length][];
+			for (int inter = 0; inter < sens_to_res.length; inter++) {
+				for (int isens = 0; isens < sens_to_res[0].length; isens++) {
+					density_rslt[inter][isens] = inverseLinTFunc(
+							noise_lev, // double [] x_val,
+							plots_direct[pt_density][inter][isens], // double [] y_val,
+							density_min, // double    y_min,
+							density_max, // double    y_max,
+							num_y_steps+1)[1]; //int       npoints		) 
+				}
+			}
+			for (int i = 0; i < density_vals.length; i++) {
+				System.out.print(String.format("%8.5f, ", density_vals[i]));
+				for (int inter = 0; inter < sens_to_res.length; inter++) {
+					for (int isens = 1; isens < sens_to_res[0].length; isens++) {
+						double ratio = density_rslt[inter][isens][i]/density_rslt[inter][0][i];
+						if ((ratio > rmse_max_ratio) || Double.isInfinite(ratio)) {
+							ratio = Double.NaN;
+						}
+						if (!Double.isNaN(ratio)) {
+							System.out.print (String.format("%8.5f", ratio));
+						}
+						if (isens < (sens_to_res[0].length - 1)) {
+							System.out.print (", ");
+						}
+					}
+					if (inter < (sens_to_res.length - 1)) {
+						System.out.print (", ");
+					}
+				}
+				System.out.println();
+			}
+			System.out.println("\n");
+			// compare RMSE : 4 to binocular, 8 to 4, 16 to 4
+			String [] col_titles_density1 = {
+					"density("+rslt_err[pt_density]+")",
+					"4:2 intra",
+					"8:4 intra",
+					"16:8 intra",
+					"4:2 inter",
+					"8:4 inter",
+					"16:8 inter"};
+			for (int i = 0; i < col_titles_density1.length; i++) {
+				System.out.print(col_titles_density1[i]);
+				if (i < (col_titles_density1.length -1)) {
+					System.out.print(", ");
+				}
+			}
+			System.out.println();
+			
+			for (int i = 0; i < density_vals.length; i++) {
+				System.out.print(String.format("%8.5f, ", density_vals[i]));
+				for (int inter = 0; inter < sens_to_res.length; inter++) {
+					for (int isens = 1; isens < sens_to_res[0].length; isens++) {
+						double ratio = density_rslt[inter][isens][i]/density_rslt[inter][isens - 1][i];
+						if ((ratio > rmse_max_ratio1) || Double.isInfinite(ratio)) {
+							ratio = Double.NaN;
+						}
+						if (!Double.isNaN(ratio)) {
+							System.out.print (String.format("%8.5f", ratio));
+						}
+						if (isens < (sens_to_res[0].length - 1)) {
+							System.out.print (", ");
+						}
+					}
+					if (inter < (sens_to_res.length - 1)) {
+						System.out.print (", ");
+					}
+				}
+				System.out.println();
+			}
+			System.out.println("\n");
+			// compare RMSE : inter to same intra
+			String [] col_titles_density2 = {
+					"density("+rslt_err[pt_density]+")",
+					"inter:intra (2)",
+					"inter:intra (4)",
+					"inter:intra (8)",
+					"inter:intra (16)"};
+			for (int i = 0; i < col_titles_density2.length; i++) {
+				System.out.print(col_titles_density2[i]);
+				if (i < (col_titles_density2.length -1)) {
+					System.out.print(", ");
+				}
+			}
+			System.out.println();
+			
+			for (int i = 0; i < density_vals.length; i++) {
+				System.out.print(String.format("%8.5f, ", density_vals[i]));
+				for (int isens = 0; isens < sens_to_res[0].length; isens++) {
+					double ratio = density_rslt[1][isens][i]/density_rslt[0][isens][i];
+					if ((ratio > rmse_max_ratio2) || Double.isInfinite(ratio)) {
+						ratio = Double.NaN;
+					}
+					if (!Double.isNaN(ratio)) {
+						System.out.print (String.format("%8.5f", ratio));
+					}
+					if (isens < (sens_to_res[0].length - 1)) {
+						System.out.print (", ");
+					}
+				}
+				System.out.println();
+			}
+			System.out.println("\n");
+		}
+	}
+		
 	
 	
+	public double [][] inverseLinTFunc(
+			double [] x_val,
+			double [] y_val,
+			double    y_min,
+			double    y_max,
+			int       npoints
+			){
+		double [][] x_y = new double [2][npoints];
+		for (int i = 0; i < npoints; i++) {
+			double y = y_min + i * (y_max - y_min)/(npoints -1);
+			x_y[0][i] = y;
+			x_y[1][i] = Double.NaN;
+			for (int j = 0; j < (x_val.length - 1); j++){
+				if (y_val[j] == y) {
+					x_y[1][i] = x_val[j];
+					break;
+				} else 	if ((y_val[j] - y) * (y_val[j + 1] - y) <= 0) {
+					x_y[1][i] = x_val[j] + (x_val[j+1] - x_val[j]) * (y - y_val[j])/(y_val[j+1]-y_val[j]);
+					break;
+				}
+			}
+			
+		}
+		return x_y;
+	}
+	
 	
 	public void batchLwirRig(
 			QuadCLT                                              quadCLT_main, // tiles should be set
@@ -11531,12 +11973,15 @@ if (debugLevel > -100) return true; // temporarily !
 
 		}
 		if (!path.contains(Prefs.getFileSeparator())) {
+			/*
 			  String x3d_path= quadCLT_main.correctionsParameters.selectX3dDirectory( // for x3d and obj
 					  quadCLT_main.correctionsParameters.getModelName(quadCLT.image_name), // quad timestamp. Will be ignored if correctionsParameters.use_x3d_subdirs is false
 					  quadCLT_main.correctionsParameters.x3dModelVersion,
 						  true,  // smart,
 						  true);  //newAllowed, // save
-			  path = x3d_path+Prefs.getFileSeparator()+path;
+			 */
+			String x3d_path = quadCLT_main.getX3dDirectory();
+			path = x3d_path+Prefs.getFileSeparator()+path;
 		}
 		if (properties == null) {
 			properties = new Properties();