Fixed old (pre 05.2024) bug that was deleting ERS velocities

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

@@ -123,6 +123,7 @@ import com.elphel.imagej.tileprocessor.ImageDtt;
 import com.elphel.imagej.tileprocessor.MLStats;
 import com.elphel.imagej.tileprocessor.MultisceneLY;
 import com.elphel.imagej.tileprocessor.QuadCLT;
+import com.elphel.imagej.tileprocessor.QuadCLTCPU;
 import com.elphel.imagej.tileprocessor.SymmVector;
 import com.elphel.imagej.tileprocessor.TwoQuadCLT;
 import com.elphel.imagej.tileprocessor.lwoc.LwirWorld;
@@ -855,6 +856,7 @@ public class Eyesis_Correction implements PlugIn, ActionListener {
 			addButton("Generate DATI",  panelLWIRWorld, color_process);
 			addButton("Create mine",  panelLWIRWorld, color_process);
 			addButton("Pattern correlate",  panelLWIRWorld, color_process);
+			addButton("Properties compare",  panelLWIRWorld, color_process);
 			plugInFrame.add(panelLWIRWorld);
 			
 		}
@@ -1593,7 +1595,7 @@ public class Eyesis_Correction implements PlugIn, ActionListener {
 			if (label.equals("Save Clean")) {
 				PROPERTIES = new Properties();
 			}
-			saveProperties(null, CORRECTION_PARAMETERS.resultsDirectory, true, PROPERTIES);
+			saveProperties(null, CORRECTION_PARAMETERS.resultsDirectory, true, PROPERTIES, DEBUG_LEVEL);
 			return;
 			/* ======================================================================== */
 		} else if (label.equals("Save offset")) {
@@ -5541,7 +5543,7 @@ public class Eyesis_Correction implements PlugIn, ActionListener {
 						SYNC_COMMAND.stopRequested, // AtomicInteger stopRequested,
 						MASTER_DEBUG_LEVEL); // int debug_level);
 			}
-			FOOTAGE_ORGANIZE.OrganizeSeries(
+			FootageOrganize.OrganizeSeries(
 					CLT_PARAMETERS,
 					CALIBRATION_ILLUSTRATION);
 		} else if (label.equals("Super batch")) {
@@ -5776,6 +5778,8 @@ public class Eyesis_Correction implements PlugIn, ActionListener {
 				return;
 			}
 			OrthoMap.testPatternCorrelate(imp_sel);
+		} else if (label.equals("Properties compare")) {
+			QuadCLTCPU.compareProperties();			
 		}
 	}
 // 
@@ -9691,14 +9695,17 @@ public class Eyesis_Correction implements PlugIn, ActionListener {
 		saveProperties(path + "_" + IJ.d2s(0.000001 * (System.nanoTime() / 1000), 6).replace('.', '_'), // full path or
 																										// null
 				directory, // use as default directory if path==null
-				useXML, properties);
+				useXML,
+				properties,
+				DEBUG_LEVEL);
 
 	}
 
 	public void saveProperties(String path, // full path or null
-			String directory, // use as default directory if path==null
-			boolean useXML,
-			Properties properties) {
+			String     directory, // use as default directory if path==null
+			boolean    useXML,
+			Properties properties,
+			int        debugLevel) {
 		String[] XMLPatterns = { ".corr-xml", ".xml" };
 		String[] confPatterns = { ".conf" };
 		String[] patterns = useXML ? XMLPatterns : confPatterns;
@@ -9758,7 +9765,7 @@ public class Eyesis_Correction implements PlugIn, ActionListener {
 			// TODO Auto-generated catch block
 			e.printStackTrace();
 		}
-		if (DEBUG_LEVEL > -3)
+		if (debugLevel > -3)
 			System.out.println("Configuration parameters are saved to " + path);
 
 	}

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

@@ -2618,9 +2618,123 @@ public class Correlation2d {
 		}
 		x0 += sx / s0; // relative to top-left
 		y0 += sy / s0;
+		
+		//https://users.cs.utah.edu/~tch/CS4640/resources/A%20geometric%20interpretation%20of%20the%20covariance%20matrix.pdf
+		double cxx = sx2 - sx * sx / s0, cyy= sy2 - sy * sy / s0, cxy = sxy - sx * sy / s0; 
+		Matrix covar = new Matrix(new double[][] {{cxx, cxy},{cxy,cyy}});
+		EigenvalueDecomposition eig = covar.eig(); 
+		double [] eigval = {eig.getD().get(0, 0),eig.getD().get(1, 1)};
+		double [][] eigvec = eig.getV().getArray(); // columns - vectors?
+		int eig_indx = (eigval[0] > eigval[1]) ? 1 : 0;
+		double [] rslt = {
+				x0 - center_xy,
+				y0 - center_xy,
+				mx,
+				eigvec[0][eig_indx],
+				eigvec[1][eig_indx],
+				eigval[eig_indx],
+				eigval[1-eig_indx]};
+		if (debug){
+			System.out.println("getMaxXYCm() -> "+rslt[0]+":"+rslt[1]+" ("+rslt[2]+
+					"), eigv0=["+rslt[3]+","+rslt[4]+"], lambda0="+rslt[5]+", lambda1="+rslt[6]);
+		}
+		if (debug_data != null) {
+			debug_data[0] = data.clone();
+			for (int i = 0; i < data.length; i++) {
+				if (!en_data[i]) {
+					debug_data[0][i] = Double.NaN;		
+				}
+			}
+		}
+		return rslt;
+	}
+	
+	public static double [] getMaxXYCmEig(
+			double []  data, // will be modified if fpn_mask != null;
+			int        data_width,      //  = 2 * transform_size - 1;
+			double     abs_min,
+			double     rel_min,
+			double     min_peak,
+			double     eig_sub_frac, // subtract fraction of threshold {eig_min_abs,eig_min_rel} after selecting by them (0 - select only, will have pedestal)
+			int        refine,       //  re-center window around new maximum. 0 -no refines (single-pass)
+			double     eig_sub_frac1,// subtract during refine (may be 0)
+			double     scale_axes,   // 1.2 scale half-axes of the ellipse: 1.0 <-> sqrt(eigenvalue)
+			double     inc_axes,     // 1.0 add do half-axes
+			boolean [] fpn_mask,
+			boolean    ignore_border,// only if fpn_mask != null - ignore tile if maximum touches fpn_mask
+			double [][] debug_data,  // null or double [1+][]
+			boolean    eig_fast2x2,  // use fast eigenvectors for 2x2 matrices
+			boolean    debug)
+	{
+		int data_height = data.length/data_width;
+		int center_xy = (data_width - 1)/2; //  = transform_size - 1;
+		double x0 = center_xy, y0 = center_xy;
+		int imax= 0;
+		for (int i= 1; i < data.length;i++) {
+			if (Double.isNaN(data[i])) {
+				System.out.println("NaN in getMaxXYCmEig()");
+				return null;
+			}
+			if (data[i] > data[imax]) {
+				imax = i;
+			}
+		}
+		if (data[imax] < min_peak) {
+			return null; // too weak;even before fpn filter
+		}
+		int ix0 = imax % data_width;
+		int iy0 = imax / data_width;
+		x0 = ix0;
+		y0 = iy0;
+		//min_peak
+		// if (fpn_mask != null
+		if (fpn_mask != null) { // modifies data, returns null if hits fpn
+			for (int i = 0; i < fpn_mask.length; i++) if (fpn_mask[i]) {
+				int iy = i / data_width; 
+				int ix = i - iy * data_width; 
+				if (ignore_border) {
+					if(((ix - ix0) <= 1) && ((ix - ix0) >= -1) && ((iy - iy0) <= 1) && ((iy - iy0) >= -1)) {
+						return null; // new double[3];
+					}
+				}
+				int ix1 = 2 * ix0 - ix;
+				if ((ix1 >= 0) && (ix1 < data_width)) {
+					int iy1 = 2 * iy0 - iy;
+					if ((iy1 >= 0) && (iy1 < data_height)) {
+						data[iy1 * data_width + ix1] = 0.0; // zero out symmetrical to fpn mask around integer maximum
+					}
+				}
+			}
+			// update imax
+			imax= 0;
+			for (int i= 1; i < data.length;i++) {
+				if (data[i] > data[imax]) {
+					imax = i;
+				}
+			}
+			if (data[imax] < min_peak) {
+				return null; // too weak after fpn filter
+			}
+			
+		}
+			// create mask of connected to max pixels
+		double mx = data[imax];
+		double min_d = Math.min(abs_min, rel_min*mx);
+		//
+		double sub_pedestal =  min_d * eig_sub_frac;
+		double sub_pedestal1 = min_d * eig_sub_frac1;
+		
+		boolean [] above_threshold = new boolean [data.length];
+		for (int i = 0; i < data.length; i++) {
+			above_threshold[i] = data[i] >= min_d;
+		}
 
-		/*
-		double s0 = 0, sx=0,sy = 0;
+		boolean [] en_data = (new TileNeibs(data_width, data_height)).getConnected(
+				above_threshold, // boolean [] tiles,
+				ix0,             // int        seedX,
+				iy0);            // int        seedY)
+		// find centroid
+		double s0 = 0, sx=0,sy = 0, sx2 = 0, sy2=0, sxy = 0;
 		for (int iy = 0; iy < data_height; iy++) {
 			double y = iy - y0;
 			for (int ix = 0; ix < data_width; ix++) {
@@ -2631,44 +2745,56 @@ public class Correlation2d {
 					s0 += d;
 					sx += d * x;
 					sy += d * y;
-
-				}
-			}
-		}
-		x0 += sx / s0; // relative to top-left
-		y0 += sy / s0;
-		double sx2 = 0, sy2=0, sxy = 0;
-		for (int iy = 0; iy < data_height; iy++) {
-			double y = iy - y0;
-			for (int ix = 0; ix < data_width; ix++) {
-				int indx = iy * data_width + ix;
-				if (en_data[indx]) { // assumes d >0, as it is >= min_d  
-					double x = ix - x0;
-					double d =  data[iy * data_width + ix] - sub_pedestal;
 					sx2 += d * x * x;
 					sy2 += d * y * y;
 					sxy += d * x * y;
 				}
 			}
 		}
-
-		 */
+		x0 += sx / s0; // relative to top-left
+		y0 += sy / s0;
 		
 		//https://users.cs.utah.edu/~tch/CS4640/resources/A%20geometric%20interpretation%20of%20the%20covariance%20matrix.pdf
 		double cxx = sx2 - sx * sx / s0, cyy= sy2 - sy * sy / s0, cxy = sxy - sx * sy / s0; 
-		Matrix covar = new Matrix(new double[][] {{cxx, cxy},{cxy,cyy}});
-		EigenvalueDecomposition eig = covar.eig(); 
-		double [] eigval = {eig.getD().get(0, 0),eig.getD().get(1, 1)};
-		double [][] eigvec = eig.getV().getArray(); // columns - vectors?
-		int eig_indx = (eigval[0] > eigval[1]) ? 1 : 0;
-		double [] rslt = {
-				x0 - center_xy,
-				y0 - center_xy,
-				mx,
-				eigvec[0][eig_indx],
-				eigvec[1][eig_indx],
-				eigval[eig_indx],
-				eigval[1-eig_indx]};
+		double [][] acovar =  {{cxx, cxy},{cxy,cyy}};
+		// TODO: calculate for 2x2 faster?
+		double [] xy = {x0-center_xy, y0-center_xy}; // relative to the center
+		double [] rslt;
+		if (eig_fast2x2) {
+			double [] eig = getEigen2x2(acovar);
+			rslt = new double[] {
+					xy[0], // x0 - center_xy,
+					xy[1], // y0 - center_xy,
+					mx,
+					eig[0],
+					eig[1],
+					eig[2],
+					eig[3]};
+		} else {
+			Matrix covar = new Matrix(acovar);
+			EigenvalueDecomposition eig = covar.eig(); 
+			double [] eigval = {eig.getD().get(0, 0),eig.getD().get(1, 1)};
+			double [][] eigvec = eig.getV().getArray(); // columns - vectors?
+			int eig_indx = (eigval[0] > eigval[1]) ? 1 : 0;
+			double k = (eigvec[0][eig_indx] > 0) ? 1 : -1;
+			rslt = new double[] {
+					xy[0], // x0 - center_xy,
+					xy[1], // y0 - center_xy,
+					mx,
+					k * eigvec[0][eig_indx],
+					k * eigvec[1][eig_indx],
+					eigval[eig_indx],
+					eigval[1-eig_indx]};
+			if (debug){
+				 double [] r = getEigen2x2(acovar);
+				 System.out.println(String.format("%10f %10f", rslt[3], r[0])); 
+				 System.out.println(String.format("%10f %10f", rslt[4], r[1])); 
+				 System.out.println(String.format("%10f %10f", rslt[5], r[2])); 
+				 System.out.println(String.format("%10f %10f", rslt[6], r[3])); 
+			}
+		}
+		
+		
 		if (debug){
 			System.out.println("getMaxXYCm() -> "+rslt[0]+":"+rslt[1]+" ("+rslt[2]+
 					"), eigv0=["+rslt[3]+","+rslt[4]+"], lambda0="+rslt[5]+", lambda1="+rslt[6]);
@@ -2681,8 +2807,89 @@ public class Correlation2d {
 				}
 			}
 		}
+		if (refine > 0) {
+			/*
+			double [] half_axes = {
+					scale_axes*Math.sqrt(eigval[  eig_indx])+inc_axes,
+					scale_axes*Math.sqrt(eigval[1-eig_indx])+inc_axes};
+			double [][] transf  = {
+					{eigvec[0][eig_indx]/half_axes[0],   eigvec[1][eig_indx]/half_axes[0]},
+					{eigvec[1][eig_indx]/half_axes[1],  -eigvec[0][eig_indx]/half_axes[1]}};
+			*/
+			double [] half_axes = {
+					scale_axes*Math.sqrt(rslt[5])+inc_axes,
+					scale_axes*Math.sqrt(rslt[6])+inc_axes};
+			double [][] transf  = {
+					{rslt[3]/half_axes[0],   rslt[4]/half_axes[0]},
+					{rslt[4]/half_axes[1],  -rslt[3]/half_axes[1]}};
+			for (int nref = 1; nref <= refine; nref++) {
+				boolean showdbg =  ((debug_data != null) && (debug_data.length > nref));
+				if (showdbg) {
+					debug_data[nref] = new double [data.length];
+					Arrays.fill(debug_data[nref], Double.NaN);
+				}
+				s0 = 0; sx=0; sy = 0;
+				for (int iy = 0; iy < data_height; iy++) {
+					double y = iy - xy[1] - center_xy;
+					for (int ix = 0; ix < data_width; ix++) {
+						double x = ix - xy[0] - center_xy;
+						double de0 = transf[0][0]*x + transf[0][1]*y;
+						double de1 = transf[1][0]*x + transf[1][1]*y;
+						double er2 = de0*de0+de1*de1;
+						if (er2 < 1.0) {
+							double er = Math.sqrt(er2);
+							double d =  data[iy * data_width + ix] - sub_pedestal1;
+							if (d > 0) { // ignore negative
+								d *=  Math.cos(0.5*Math.PI*er);
+								s0 += d;
+								sx += d * x;
+								sy += d * y;
+								if (showdbg) {
+									debug_data[nref][iy * data_width + ix] = d;
+								}
+							}
+						}
+					}
+				}
+				if (s0 >0) {
+					xy[0] += sx / s0;
+					xy[1] += sy / s0;
+				}
+			}
+			rslt[0] = xy[0];
+			rslt[1] = xy[1];
+		}
+
+		
+		
 		return rslt;
 	}
+	/**
+	 * Calculates eigenvalues and eigenvectors for a 2x2 (covariance) matrix,
+	 * returns normalized eigenvector for the smaller eigenvalue (other is orthogonal)
+	 * followed by 2 eigenvalues (smaller first)  
+	 * @param covar 2x2 covariance matrix as [2][2] array
+	 * @return 4-element array: first eigenvector x, y, lambda0, lambda1 (lambda0 <= lampda1)
+	 */
+	public static double [] getEigen2x2(double [][] A) { //
+		double e = 1E-12;
+		double hapd = (A[0][0]+A[1][1])/2;
+		double hamd = (A[0][0]-A[1][1])/2;
+		double bc = A[1][0]*A[0][1];
+		double d = Math.sqrt(hamd*hamd + bc);
+		if (d/Math.abs(hapd) < e) {
+			return new double [] {1,0,hapd,hapd};
+		}
+		double [] lambda = {hapd-d, hapd+d};
+		double [] v0 = {A[0][0] - lambda[1], A[1][0]};
+		double  k = 1.0/Math.sqrt(v0[0]*v0[0] + v0[1]*v0[1]);
+		if (v0[0] < 0) { // to be compatible with standard?
+			k = -k;
+		}
+		return new double [] {k*v0[0], k*v0[1], lambda[0], lambda[1]};
+	}
+	
+	
 	
 	/**
 	 * Find maximum of the 2d array projected on a specified vector using centroid.

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

@@ -2375,6 +2375,13 @@ public class ImageDtt extends ImageDttCPU {
 			final double              eigen_min_abs,    // 0.05 values below this do not count for covariance
 			final double              eigen_min_rel,    // 0.2  values less than this fraction of tile's max do not count for covariance
 			final double              eig_sub_frac,     // 1.0;   // subtract fraction of threshold {eig_min_abs,eig_min_rel} after selecting by them (0 - select only, will have pedestal)
+			
+			final  int                eig_recenter,     // 2;     // Re-center window around new maximum. 0 -no refines (single-pass)
+			final double              eig_sub_frac1,    // 0.0;   // subtract during refine (may be 0)
+			// Using eigenvectors/values to create ellipse (slightly larger) and use it for a cosine mask (as round before) to refine centers
+			final double              eig_scale_axes,   // 1.2;   // scale half-axes of the ellipse: 1.0 <-> sqrt(eigenvalue)
+			final double              eig_inc_axes,     // 1.0;   // add do half-axes
+			final boolean             eig_fast2x2,      // use fast eigenvectors for 2x2 matrices
 			final double [][]         eigen,            // null or [tilesX*tilesY]{lamb0_x,lamb0_y, lamb0, lamb1} eigenvector0[x,y],lam0,lam1
 			final boolean             eigen_debug,      //
 			final int                 debug_tileX,
@@ -2468,10 +2475,11 @@ public class ImageDtt extends ImageDttCPU {
 
 		// currently execCorr2D_normalize() output has 17 slices for old variant (no neibs) and 18/2 if (use_neibs)
 		final int extra_len =     extra_sum? 1 : 0;
-		final int extra_len_eig = eigen_debug? 4: 0; // all, all_remain, weak, weak_remain 
+		// eig_recenter
+		final int extra_len_eig = eigen_debug? (4 + eig_recenter): 0; // all, all_remain, weak, weak_remain 
 //		final int corrs_len = ((use_partial || use_neibs) ? used_sensors_list.length:1); // without optional extra_len but including GPU sum
 		final int corrs_len = (use_neibs || use_partial) ? used_sensors_list.length:1; // without optional extra_len but including GPU sum
-		final int eigen_indx = (extra_len_eig > 0) ? (corrs_len + extra_len):-1;
+		final int eigen_indx = (extra_len_eig > 0) ? (corrs_len + extra_len):-1; // first of eigen-related data
 		final int indx_sum_pd =      (extra_len > 0) ? corrs_len : -1; 
 		final int indx_sum_td =      use_neibs ? (corrs_len -2): (corrs_len -1); 
 		final int indx_sum_td_neib = use_neibs ? (corrs_len -1): -1; 
@@ -2621,7 +2629,7 @@ public class ImageDtt extends ImageDttCPU {
 						// only two cases:
 						// 1 - TD only (!neib_notd_only)
 						// 2 - TD, then TD-neibs if failed (neib_notd_only)
-						double [][] debug_data = ((dcorr_tiles != null) && (eigen_indx >=0)) ? (new double[1][]):null;
+						double [][] debug_data = ((dcorr_tiles != null) && (eigen_indx >=0)) ? (new double[1+eig_recenter][]):null;
 						double [] stats_mv = Correlation2d.getMaxXYCmEig(
 								corrs[indx_sum_td],    // double []  data, // will be modified if fpn_mask != null;
 								corr_size,             // int        data_width,      //  = 2 * transform_size - 1;
@@ -2629,10 +2637,15 @@ public class ImageDtt extends ImageDttCPU {
 								eigen_min_rel,         // double     rel_min,
 								eig_str_sum,           // double     min_peak,
 								eig_sub_frac,          // double     eig_sub_frac, // subtract fraction of threshold {eig_min_abs,eig_min_rel} after selecting by them (0 - select only, will have pedestal)
+								eig_recenter,          // int        refine,       //  re-center window around new maximum. 0 -no refines (single-pass)
+								eig_sub_frac1,         // double     eig_sub_frac1,// subtract during refine (may be 0)
+								eig_scale_axes,        // double     scale_axes,   // 1.2 scale half-axes of the ellipse: 1.0 <-> sqrt(eigenvalue)
+								eig_inc_axes,          // double     inc_axes,     // 1.0 add do half-axes
 								fpn_mask,              // boolean [] fpn_mask,
 								false,                 // boolean    ignore_border, // only if fpn_mask != null - ignore tile if maximum touches fpn_mask
 								debug_data,            // double [][] debug_data, // null or double [1]
-								false);                // boolean    debug)
+								eig_fast2x2,           // boolean    eig_fast2x2,  // use fast eigenvectors for 2x2 matrices
+								debugTile0);           // boolean    debug)
 						used_td[nTile] = stats_mv != null;
 						if (stats_mv != null) {
 							stats_mv[2] -= eig_str_sum * scale_neibs_td;
@@ -2650,10 +2663,15 @@ public class ImageDtt extends ImageDttCPU {
 									eigen_min_rel,         // double     rel_min,
 									eig_str_neib,          // double     min_peak,
 									eig_sub_frac,          // double     eig_sub_frac, // subtract fraction of threshold {eig_min_abs,eig_min_rel} after selecting by them (0 - select only, will have pedestal)
+									eig_recenter,          // int        refine,       //  re-center window around new maximum. 0 -no refines (single-pass)
+									eig_sub_frac1,         // double     eig_sub_frac1,// subtract during refine (may be 0)
+									eig_scale_axes,        // double     scale_axes,   // 1.2 scale half-axes of the ellipse: 1.0 <-> sqrt(eigenvalue)
+									eig_inc_axes,          // double     inc_axes,     // 1.0 add do half-axes
 									fpn_mask,              // boolean [] fpn_mask,
 									false,                 // boolean    ignore_border, // only if fpn_mask != null - ignore tile if maximum touches fpn_mask
 									debug_data,            // double [][] debug_data, // null or double [1]
-									false);                // boolean    debug)
+									eig_fast2x2,           // boolean    eig_fast2x2,  // use fast eigenvectors for 2x2 matrices
+									debugTile0);           // boolean    debug)
 							weak_tile[nTile] = stats_mv != null;
 							if (stats_mv != null) {
 								stats_mv[2] -= eig_str_neib * scale_neibs_td;
@@ -2671,7 +2689,12 @@ public class ImageDtt extends ImageDttCPU {
 								dcorr_tiles[iCorrTile][eigen_indx+3] = debug_data [0].clone(); // weak REMAIN
 							}
 							//used_td[nTile]
-							
+							int recenter_indx_m1 = eigen_indx+3;
+							if (stats_mv != null) {
+								for (int i = 1; (i < debug_data.length) && (recenter_indx_m1 + i < dcorr_tiles[iCorrTile].length) ; i++) {
+									dcorr_tiles[iCorrTile][recenter_indx_m1+i] = debug_data [i].clone();
+								}
+							}
 						}
 						if (stats_mv != null) {
 							if (eigen != null) {
@@ -2741,6 +2764,7 @@ public class ImageDtt extends ImageDttCPU {
 									int iCorrTile = iCorrTile_index[nTile];
 									dcorr_tiles[iCorrTile][eigen_indx+1] = null;
 									dcorr_tiles[iCorrTile][eigen_indx+3] = null;
+									// keeping extra recenter layers
 								}
 							}
 						}

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

@@ -3355,19 +3355,20 @@ public class Interscene {
 			double      eig_max_sqrt=clt_parameters.imp.eig_max_sqrt;
 			double      eig_min_sqrt=clt_parameters.imp.eig_min_sqrt;
 //			boolean dbg_images = false;
+			double [][] eigen_masked = clt_parameters.imp.eig_xy_lma? null : eigen;
 			String      dbg_prefix = dbg_images? (first_QuadClt.getImageName()+"-"+scene_QuadClt.getImageName()+"-NLMA_"+nlma) : null;
 			intersceneLma.prepareLMA(
-					scene_xyzatr0,         // final double []   scene_xyzatr0,     // camera center in world coordinates (or null to use instance)
-					scene_xyzatr_pull,     // final double []   scene_xyz_pull, // if both are not null, specify target values to pull to
-					ref_xyzatr, // _inv, // ref_xyzatr, 
+					scene_xyzatr0,       // final double []   scene_xyzatr0,     // camera center in world coordinates (or null to use instance)
+					scene_xyzatr_pull,   // final double []   scene_xyz_pull, // if both are not null, specify target values to pull to
+					ref_xyzatr,          // _inv, // ref_xyzatr, 
 					// reference atr, xyz are considered 0.0 not anymore?
 					scene_QuadClt,       // final QuadCLT     scene_QuadClt,
-					first_QuadClt,   // final QuadCLT     reference_QuadClt,
-					param_select_mod, // param_select,        // final boolean[]   param_select,
+					first_QuadClt,       // final QuadCLT     reference_QuadClt,
+					param_select_mod,    // param_select,        // final boolean[]   param_select,
 					param_regweights,    // final double []   param_regweights,
-					eig_max_sqrt, // final double      eig_max_sqrt, //  10;    // for sqrt(lambda) - consider infinity (infinite linear feature, a line)
-					eig_min_sqrt, // final double      eig_min_sqrt, //  1;    // for sqrt(lambda) - consider infinity (infinite linear feature, a line)
-					eigen,               // final double [][] eigen, // [tilesX*tilesY]{lamb0_x,lamb0_y, lamb0, lamb1} eigenvector0[x,y],lam0,lam1
+					eig_max_sqrt,        // final double      eig_max_sqrt, //  10;    // for sqrt(lambda) - consider infinity (infinite linear feature, a line)
+					eig_min_sqrt,        // final double      eig_min_sqrt, //  1;    // for sqrt(lambda) - consider infinity (infinite linear feature, a line)
+					eigen_masked,        // final double [][] eigen, // [tilesX*tilesY]{lamb0_x,lamb0_y, lamb0, lamb1} eigenvector0[x,y],lam0,lam1
 					coord_motion[1],     // final double [][] vector_XYS, // optical flow X,Y, confidence obtained from the correlate2DIterate()
 					coord_motion[0],     // final double [][] centers,    // macrotile centers (in pixels and average disparities
 					(nlma == 0),         // boolean           first_run,
@@ -3546,11 +3547,11 @@ public class Interscene {
 		if ((rms_out != null) && (intersceneLma.getLastRms() != null)) {
 			rms_out[0] = intersceneLma.getLastRms()[0];
 			rms_out[1] = intersceneLma.getLastRms()[1];
-			if (rms_out.length >=2) {
+			if (rms_out.length >=4) {
 				rms_out[2] = intersceneLma.getSumWeights();
 				rms_out[3] = intersceneLma.getNumDefined();
 				
-				if (rms_out.length >=4) {
+				if (rms_out.length >=6) {
 					if (intersceneLma.isEigenNormalized()) {
 						double [] rms_metric = intersceneLma.calcRMS(true);
 						rms_out[4] = rms_metric[0];
@@ -4296,8 +4297,6 @@ public class Interscene {
 		double [][] eigen = use_eigen? (new double[tilesX*tilesY][]) : null; 
 		boolean             eigen_debug = use_eigen && show_2d_correlations;
 		if (use_eigen) {
-			double  eigen_min_abs=  clt_parameters.imp.eig_min_abs ; //eigen_min_abs = 0.05; // 0.05 values below this do not count for covariance
-			double  eigen_min_rel = clt_parameters.imp.eig_min_rel ; //0.2;  // 0.2  values less than this fraction of tile's max do not count for covariance
 			coord_motion = image_dtt.clt_process_tl_interscene(       // convert to pixel domain and process correlations already prepared in fcorr_td and/or fcorr_combo_td
 					clt_parameters.img_dtt,            // final ImageDttParameters  imgdtt_params,   // Now just extra correlation parameters, later will include, most others
 					// only used here to keep extra array element for disparity difference
@@ -4345,9 +4344,15 @@ public class Interscene {
 					scale_neibs_pd,                    // final double              scale_neibs_pd,   // scale threshold for the pixel-domain average maximums  		
 					scale_neibs_td,                    // final double              scale_neibs_td,   // scale threshold for the transform-domain average maximums
 					scale_avg_weight,                  // final double              scale_avg_weight,  // reduce influence of the averaged correlations compared to the single-tile ones
-					eigen_min_abs,                     // final double              eigen_min_abs,    // 0.05 values below this do not count for covariance
-					eigen_min_rel,                     // final double              eigen_min_rel,    // 0.2  values less than this fraction of tile's max do not count for covariance
-					clt_parameters.imp.eig_sub_frac,   //final double              eig_sub_frac,     // 1.0;   // subtract fraction of threshold {eig_min_abs,eig_min_rel} after selecting by them (0 - select only, will have pedestal)					
+					clt_parameters.imp.eig_min_abs,    // final double              eigen_min_abs,    // 0.05 values below this do not count for covariance
+					clt_parameters.imp.eig_min_rel,    // final double              eigen_min_rel,    // 0.2  values less than this fraction of tile's max do not count for covariance
+					clt_parameters.imp.eig_sub_frac,   // final double              eig_sub_frac,     // 1.0;   // subtract fraction of threshold {eig_min_abs,eig_min_rel} after selecting by them (0 - select only, will have pedestal)
+					clt_parameters.imp.eig_recenter,   // final  int                eig_recenter,     // 2;     // Re-center window around new maximum. 0 -no refines (single-pass)
+					clt_parameters.imp.eig_sub_frac1,  // final double              eig_sub_frac1,    // 0.0;   // subtract during refine (may be 0)
+					// Using eigenvectors/values to create ellipse (slightly larger) and use it for a cosine mask (as round before) to refine centers
+					clt_parameters.imp.eig_scale_axes, // final double              eig_scale_axes,   // 1.2;   // scale half-axes of the ellipse: 1.0 <-> sqrt(eigenvalue)
+					clt_parameters.imp.eig_inc_axes,   // final double              eig_inc_axes,     // 1.0;   // add do half-axes
+					clt_parameters.imp.eig_fast2x2,       // final boolean             eig_fast2x2,         // use fast eigenvectots for 2x2 matrices
 					eigen, // null, // eigen,                             // final double [][]         eigen,            // null or [tilesX*tilesY]{lamb0_x,lamb0_y, lamb0, lamb1} eigenvector0[x,y],lam0,lam1
 					eigen_debug,                       // final double [][]         eigen_debug,      // null or [tilesX*tilesY][]
 					clt_parameters.tileX,              // final int                 debug_tileX,

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

@@ -803,6 +803,9 @@ public class IntersceneLma {
 			final double   eig_max_sqrt,
 			final double   eig_min_sqrt,
 			final double [][] eigen){ // [tilesX*tilesY]{lamb0_x,lamb0_y, lamb0, lamb1} eigenvector0[x,y],lam0,lam1
+		if (eigen == null) {
+			return null;
+		}
 		final double [][][] transform = new double[eigen.length][2][2];
 		final Thread[] threads = ImageDtt.newThreadArray(QuadCLT.THREADS_MAX);
 		final AtomicInteger ai = new AtomicInteger(0);
@@ -824,10 +827,6 @@ public class IntersceneLma {
 			};
 		}		      
 		ImageDtt.startAndJoin(threads);
-
-		
-		
-		
 		return transform;
 	}
 	

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

@@ -8767,11 +8767,15 @@ public class OpticalFlow {
 			String [] timestamps = ers_reference.getScenes();
 			for (String ts:timestamps) {
 				double [][] xyzatr = ers_reference.getSceneXYZATR(ts);
+				double [][] xyzatr_dt = ers_reference.getSceneErsXYZATR_dt(ts);
 				for (int i = 0; i < xyzatr[0].length; i++) {
 					xyzatr[0][i] *= scale;
+					xyzatr_dt[0][i] *= scale; // scale linear velocities too 06.25.2024
 				}
 //				xyzatr[0][2] -= zdiff;
-				ers_reference.addScene(ts,xyzatr);
+//				ers_reference.addScene(ts,xyzatr); // should use addScene(String timestamp, double [][] xyzatr, double [][] ers_xyzatr_dt) {
+				ers_reference.addScene(ts,xyzatr,xyzatr_dt); // should use addScene(String timestamp, double [][] xyzatr, double [][] ers_xyzatr_dt) {
+
 			}
 			ref_scene.saveInterProperties( // save properties for interscene processing (extrinsics, ers, ...)
 					null, // String path,             // full name with extension or w/o path to use x3d directory
@@ -9557,7 +9561,8 @@ public class OpticalFlow {
 		String rslt_suffix = "-INTER-INTRA-HISTORIC";
 		rslt_suffix += (clt_parameters.correlate_lma?"-LMA":"-NOLMA");
 		// ref_index from the center is moved to the very end (indx_ref)
-		if (compensate_dsi && !scenes[ref_index].interDsiExists()) {
+		// was bug - used ref_index Maybe that was intentional?
+		if (compensate_dsi && !scenes[indx_ref].interDsiExists()) {
 			System.out.println("Will compensate pose Z as average disparity may change");
 			double z_corr = compensateDSI(
 					clt_parameters,    // CLTParameters        clt_parameters,

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

@@ -61,6 +61,7 @@ import com.elphel.imagej.cameras.ColorProcParameters;
 import com.elphel.imagej.cameras.EyesisCorrectionParameters;
 import com.elphel.imagej.cameras.ThermalColor;
 import com.elphel.imagej.common.DoubleGaussianBlur;
+import com.elphel.imagej.common.GenericJTabbedDialog;
 import com.elphel.imagej.common.PolynomialApproximation;
 import com.elphel.imagej.common.ShowDoubleFloatArrays;
 import com.elphel.imagej.correction.CorrectionColorProc;
@@ -3834,6 +3835,12 @@ public class QuadCLTCPU {
     		boolean silent) {
     	String x3d_path = getX3dDirectory();
     	String file_path = x3d_path + Prefs.getFileSeparator() + image_name + suffix + ".tiff";
+    	if (!(new File(file_path).exists())) {
+    		if (!silent) {
+    			System.out.println ("File does not exist: "+file_path);
+    		}
+    		return -1;
+    	}
     	ImagePlus imp = null;
     	try {
     		imp = new ImagePlus(file_path);
@@ -3987,6 +3994,63 @@ public class QuadCLTCPU {
     	System.out.println("Restored interframe properties from :"+path);
     	return properties;
     }
+    
+    
+    public static void compareProperties() {
+    	String path1 = "/media/elphel/SSD3-4GB/lwir16-proc/berdich3/models/models_1697875868-1697879449-b/1697877436_779067/39-dbg-02-11/1697877436_779067-INTERFRAME.corr-xml";
+    	String path2 = "/media/elphel/SSD3-4GB/lwir16-proc/berdich3/models/models_1697875868-1697879449-b/1697877436_779067/37-dbg-01-11_noeig/1697877436_779067-INTERFRAME.corr-xml";
+    	GenericJTabbedDialog gd = new GenericJTabbedDialog("Set image pair",1200,800);
+    	gd.addStringField ("First  configuration path",     path1, 180, "Full path of the first configuration properties file.");	
+    	gd.addStringField ("Second configuration path",     path2, 180, "Second path of the first configuration properties file.");	
+    	gd.showDialog();
+    	if (gd.wasCanceled()) return;
+    	path1 = gd.getNextString();
+    	path2 = gd.getNextString();
+    	compareProperties(path1,path2);
+    }
+
+    public static void compareProperties(
+    		String path1,
+    		String path2) {
+    	Properties properties1, properties2;
+    	properties1 = loadProperties(
+    			path1,  // String path,
+    			null,   // Properties properties)
+    			false); // boolean silent) 
+    	properties2 = loadProperties(
+    			path2,  // String path,
+    			null,   // Properties properties)
+    			false); // boolean silent)
+    	//	Enumeration<String> en1 = (Enumeration<String>) properties1.propertyNames();
+    	//	Enumeration<String> en2 = (Enumeration<String>) properties2.propertyNames();
+    	System.out.println("path1 = "+path1);
+    	System.out.println("path2 = "+path2);
+    	int num_diff = 0;
+    	for (@SuppressWarnings("unchecked")
+    	Enumeration<String> e = (Enumeration<String>) properties1.propertyNames(); e.hasMoreElements();) {
+    		String key = e.nextElement();
+    		String v1 = properties1.getProperty(key);
+    		String v2 = properties2.getProperty(key);
+    		if (!v1.equals(v2)) {
+    			System.out.println("property1["+key+"] = " + v1+", property2["+key+"] = " +v2);
+    			num_diff++;
+    		}
+    	}
+    	for (@SuppressWarnings("unchecked")
+    	Enumeration<String> e = (Enumeration<String>) properties2.propertyNames(); e.hasMoreElements();) {
+    		String key = e.nextElement();
+    		String v1 = properties1.getProperty(key);
+    		String v2 = properties2.getProperty(key);
+    		if (v1 == null ) { // other already printed
+    			System.out.println("property1["+key+"] = " + v1+", property2["+key+"] = " +v2);
+    			num_diff++;
+    		}
+    	}
+    	System.out.println(num_diff+" differences found");
+    }
+
+
+
 
     public boolean propertiesContainString (
     		String     needle) {