Equalizing intensities

src/main/java/com/elphel/imagej/common/PolynomialApproximation.java

 package com.elphel.imagej.common;
+import java.util.concurrent.atomic.AtomicInteger;
+
+import com.elphel.imagej.tileprocessor.ImageDtt;
+import com.elphel.imagej.tileprocessor.QuadCLT;
+
 import Jama.LUDecomposition;
 import Jama.Matrix;
 
 public class PolynomialApproximation {
 	public int debugLevel=1;
-// TODO Move other methods here
+	// TODO Move other methods here
 	public PolynomialApproximation(){}
 	public PolynomialApproximation(int debugLevel){
 		this.debugLevel=debugLevel;
@@ -53,7 +58,7 @@ public class PolynomialApproximation {
 			System.out.println("polynomialApproximation1d() B:");
 			B.print(10, 5);
 		}
-//		while (!(new LUDecomposition(M)).isNonsingular() && (N1>0)){
+		//		while (!(new LUDecomposition(M)).isNonsingular() && (N1>0)){
 		while (!(new LUDecomposition(M)).isNonsingular() && (N1>=0)){ // make N=0 legal ?
 			aM=new double [N1][N1];
 			aB=new double [N1][1];
@@ -81,7 +86,7 @@ public class PolynomialApproximation {
 		for (int i=0;i<=N;i++) result[i]=(i<=N1)?aR[i][0]:0.0;
 		return result;
 	}
-/**
+	/**
 	 * Linear approximates each of 3 functions of 3 variables and finds where they are all zero
 	 * @param data: for each sample (1-st index):
 	 *               0 - {x,y,z}
@@ -312,7 +317,7 @@ public class PolynomialApproximation {
 				1.0E-15,  // threshold ratio of matrix determinant to norm for quadratic approximation (det too low - fail)
 				this.debugLevel);
 	}
-/*
+	/*
 	   public double [][] quadraticApproximation( // no use
 			   double [][][] data,
 			   boolean forceLinear,  // use linear approximation
@@ -328,7 +333,7 @@ public class PolynomialApproximation {
 				   1.0E-15,  // threshold ratio of matrix determinant to norm for quadratic approximation (det too low - fail)
 				   this.debugLevel);
 	   }
-*/
+	 */
 	public double [][] quadraticApproximation(
 			double [][][] data,
 			boolean forceLinear,  // use linear approximation
@@ -528,7 +533,7 @@ public class PolynomialApproximation {
 		if (debugLevel>3) System.out.println(">>> n="+n+" det_lin="+mLin.det()+" norm_lin="+normMatix(mAarrayL));
 		double nmL=normMatix(mAarrayL);
 		if ((nmL==0.0) || (Math.abs(mLin.det())/nmL<thresholdLin)){
-// return average value for each channel
+			// return average value for each channel
 			if (S00==0.0) return null; // not even average
 			double [][] ABCDEF=new double[zDim][3];
 			for (i=0;i<zDim;i++) {
@@ -540,7 +545,7 @@ public class PolynomialApproximation {
 		}
 		double []zAarrayL=new double [3];
 		double [][] ABCDEF=new double[zDim][];
-//			   double [] zAarrayL={SZ10,SZ01,SZ00};
+		//			   double [] zAarrayL={SZ10,SZ01,SZ00};
 		for (i=0;i<zDim;i++) {
 			zAarrayL[0]=SZ10[i];
 			zAarrayL[1]=SZ01[i];
@@ -575,7 +580,7 @@ public class PolynomialApproximation {
 					", Math.abs(M.det())/normMatix(mAarrayQ)="+(Math.abs(mQuad.det())/normMatix(mAarrayQ))); //did not happen
 			return ABCDEF; // not enough data for the quadratic approximation, return linear
 		}
-//			   double [] zAarrayQ={SZ20,SZ02,SZ11,SZ10,SZ01,SZ00};
+		//			   double [] zAarrayQ={SZ20,SZ02,SZ11,SZ10,SZ01,SZ00};
 		double [] zAarrayQ=new double [6];
 		for (i=0;i<zDim;i++) {
 			zAarrayQ[0]=SZ20[i];
@@ -589,7 +594,7 @@ public class PolynomialApproximation {
 		}
 		return ABCDEF;
 	}
-//			calculate "volume" made of the matrix row-vectors, placed orthogonally
+	//			calculate "volume" made of the matrix row-vectors, placed orthogonally
 	// to be compared to determinant
 	public double normMatix(double [][] a) {
 		double d,norm=1.0;
@@ -601,6 +606,161 @@ public class PolynomialApproximation {
 		return norm;
 	}
 
+	public static double [] getYXRegression(
+			final double []  data_x,
+			final double []  data_y,
+			final boolean [] mask) {
+		final Thread[] threads = ImageDtt.newThreadArray(QuadCLT.THREADS_MAX);
+		final AtomicInteger ai = new AtomicInteger(0);
+		final double [] as0 =  new double[threads.length];
+		final double [] asx =  new double[threads.length];
+		final double [] asx2 = new double[threads.length];
+		final double [] asy =  new double[threads.length];
+		final double [] asxy=  new double[threads.length];
+		final AtomicInteger ati = new AtomicInteger(0);
+		ai.set(0);
+		for (int ithread = 0; ithread < threads.length; ithread++) {
+			threads[ithread] = new Thread() {
+				public void run() {
+					int thread_num = ati.getAndIncrement();
+					for (int ipix = ai.getAndIncrement(); ipix < data_x.length; ipix = ai.getAndIncrement()) if ((mask == null) || mask[ipix]){
+						double x = data_x[ipix];
+						double y = data_y[ipix];
+						if (!Double.isNaN(x) && !Double.isNaN(y)) {
+							as0 [thread_num] += 1;
+							asx [thread_num] += x;
+							asx2[thread_num] += x * x;
+							asy [thread_num] += y;
+							asxy[thread_num] += x * y;
+						}
+					} // for (int ipix
+				}
+			};
+		}		      
+		ImageDtt.startAndJoin(threads);
+		double s0 = 0.0;
+		double sx = 0.0;
+		double sx2 = 0.0;
+		double sy = 0.0;
+		double sxy= 0.0;
+		for (int i = 0; i < threads.length; i++) {
+			s0+=  as0[i];
+			sx+=  asx[i];
+			sx2+= asx2[i];
+			sy+=  asy[i];
+			sxy+= asxy[i];
+		}
+		double dnm = s0 * sx2 - sx*sx;
+		double a = (sxy * s0 - sy * sx) / dnm;
+		double b = (sy * sx2 - sxy * sx) / dnm;
+		return new double[] {a,b};
+	}
 	
-//RuntimeException
+	/**
+	 * Get best fit ax+b symmetrical for X and Y, same weights
+	 * https://en.wikipedia.org/wiki/Deming_regression#Orthogonal_regression
+	 * @param data_x
+	 * @param data_y
+	 * @param mask
+	 * @return
+	 */
+	
+	public static double [] getOrthoRegression( // symmetrical for X and Y, same eror weights
+			final double []  data_x,
+			final double []  data_y,
+			final boolean [] mask_in) {
+		final boolean [] mask = new boolean [data_x.length];
+		final Thread[] threads = ImageDtt.newThreadArray(QuadCLT.THREADS_MAX);
+		final AtomicInteger ai = new AtomicInteger(0);
+		final double [] as0 =  new double[threads.length];
+		final double [] asx =  new double[threads.length];
+		final double [] asy =  new double[threads.length];
+		final AtomicInteger ati = new AtomicInteger(0);
+		ai.set(0);
+		// Find centroid first
+		for (int ithread = 0; ithread < threads.length; ithread++) {
+			threads[ithread] = new Thread() {
+				public void run() {
+					int thread_num = ati.getAndIncrement();
+					for (int ipix = ai.getAndIncrement(); ipix < data_x.length; ipix = ai.getAndIncrement()) if ((mask_in == null) || mask_in[ipix]){
+						
+						double x = data_x[ipix];
+						double y = data_y[ipix];
+						if (!Double.isNaN(x) && !Double.isNaN(y)) {
+							as0 [thread_num] += 1;
+							asx [thread_num] += x;
+							asy [thread_num] += y;
+							mask[ipix] = true;
+						}
+					} // for (int ipix
+				}
+			};
+		}		      
+		ImageDtt.startAndJoin(threads);
+		double s0 = 0.0;
+		double sx = 0.0;
+		double sy = 0.0;
+		for (int i = 0; i < threads.length; i++) {
+			s0+=  as0[i];
+			sx+=  asx[i];
+			sy+=  asy[i];
+		}
+		double [] z_mean = {sx/s0, sy/s0}; // complex
+		final double [] as_re =  new double[threads.length];
+		final double [] as_im =  new double[threads.length];
+		ai.set(0);
+		ati.set(0);
+		for (int ithread = 0; ithread < threads.length; ithread++) {
+			threads[ithread] = new Thread() {
+				public void run() {
+					int thread_num = ati.getAndIncrement();
+					for (int ipix = ai.getAndIncrement(); ipix < data_x.length; ipix = ai.getAndIncrement()) if ((mask == null) || mask[ipix]){
+						double x = data_x[ipix]-z_mean[0];
+						double y = data_y[ipix]-z_mean[1];
+						as_re[thread_num] += x*x - y*y;
+						as_im[thread_num] += 2*x*y;
+					} // for (int ipix
+				}
+			};
+		}		      
+		ImageDtt.startAndJoin(threads);
+		double s_re = 0.0;
+		double s_im = 0.0;
+		for (int i = 0; i < threads.length; i++) {
+			s_re+=  as_re[i];
+			s_im+=  as_im[i];
+		}
+		// https://en.wikipedia.org/wiki/Square_root#Algebraic_formula
+		// sqrt (s_re+i*s_im)
+		double sqrt_re = Math.sqrt(0.5 * (Math.sqrt(s_re*s_re + s_im*s_im) +s_re));
+		double sqrt_im = ((s_im > 0)? 1 : -1) *  Math.sqrt(0.5 * (Math.sqrt(s_re*s_re + s_im*s_im) - s_re));
+		double a = sqrt_im/sqrt_re;
+		double b = z_mean[1]-  a * z_mean[0];
+		return new double[] {a,b};
+	}
+	
+	
+	
+	
+	
+	
+	public static void applyRegression(
+			final double [] data, // clone by caller
+			final double [] regression) {
+		final Thread[] threads = ImageDtt.newThreadArray(QuadCLT.THREADS_MAX);
+		final AtomicInteger ai = new AtomicInteger(0);
+		for (int ithread = 0; ithread < threads.length; ithread++) {
+			threads[ithread] = new Thread() {
+				public void run() {
+					for (int ipix = ai.getAndIncrement(); ipix < data.length; ipix = ai.getAndIncrement()){
+						data[ipix] = regression[0] * data[ipix] + regression[1];
+					}
+				}
+			};
+		}		      
+		ImageDtt.startAndJoin(threads);
+	}
+	public static double [] invertRegression(double [] regression) {
+		return new double [] {1.0/regression[0], -regression[1]/regression[0]};
+	}
 }

src/main/java/com/elphel/imagej/orthomosaic/ComboMatch.java

@@ -171,7 +171,7 @@ public class ComboMatch {
 		boolean pattern_match =       true; // false;
 		
 		boolean bounds_to_indices =   true;
-		int     temp_mode =           1;
+		int     temp_mode =           0;
 		boolean restore_temp =        true;
 		double  frac_remove  =        clt_parameters.imp.pmtch_frac_remove; //   0.15;
 		double  metric_error =        clt_parameters.imp.pmtch_metric_err; // 0.05; // 0.02;//  2 cm
@@ -1030,10 +1030,13 @@ public class ComboMatch {
         				}
         				if (render_match) {
         					String title=String.format("multi_%03d-%03d_%s-%s_zoom%d_%d",gpu_pair[0],gpu_pair[1],gpu_spair[0],gpu_spair[1],min_zoom_lev,zoom_lev);
+        					// Avoid renderMulti() - it duplicates code renderMultiDouble()
+        					int eq_mode = 2; // calculate
         					ImagePlus imp_img_pair = 	maps_collection.renderMulti (
         							//_zoom<integer> is needed for opening with "Extract Objects" command
         							title,             // String      title,
-        							OrthoMapsCollection.MODE_IMAGE,  // int           mode,    // 0 - regular image, 1 - altitudes, 2 - black/white mask       // boolean     use_alt,
+//        							OrthoMapsCollection.MODE_IMAGE,  // int           mode,    // 0 - regular image, 1 - altitudes, 2 - black/white mask       // boolean     use_alt,
+        							eq_mode,           //int           eq_mode, // 0 - ignore equalization, 1 - use stored equalization, 2 - calculate equalization 
         							gpu_pair,          // int []        indices, // null or which indices to use (normally just 2 for pairwise comparison)
         							bounds_to_indices, // boolean       bounds_to_indices,
         							temp_mode,         // int           temp_mode, // 0 - do nothing, 1 - equalize average,2 - try to correct

src/main/java/com/elphel/imagej/orthomosaic/OrthoMap.java

@@ -115,6 +115,7 @@ public class OrthoMap implements Comparable <OrthoMap>, Serializable{
 	public transient double                     agl = Double.NaN;
 	public transient int                        num_scenes = -1;;   // number of scenes that made up this image
 	public transient double                     sfm_gain = Double.NaN;     // maximal SfM gain of this map
+	public transient double []                  equalize = {1,0}; // rectified value = equalize[0]*source_value+equalize[1]
 	private void writeObject(ObjectOutputStream oos) throws IOException {
 		oos.defaultWriteObject();
 		oos.writeObject(path);
@@ -137,6 +138,7 @@ public class OrthoMap implements Comparable <OrthoMap>, Serializable{
 		oos.writeObject(agl);
 		oos.writeObject(num_scenes);
 		oos.writeObject(sfm_gain);
+		oos.writeObject(equalize);
 	}
 	
 	private void readObject(ObjectInputStream ois) throws ClassNotFoundException, IOException {
@@ -160,6 +162,8 @@ public class OrthoMap implements Comparable <OrthoMap>, Serializable{
 		agl = (double)  ois.readObject();
 		num_scenes = (int) ois.readObject();
 		sfm_gain = (double) ois.readObject();
+//		equalize = new double[] {1,0};
+		equalize = (double []) ois.readObject();
 		images = new HashMap <Integer, FloatImageData>(); // field images was not saved
 		averageImagePixel = Double.NaN; // average image pixel value (to combine with raw)
 		
@@ -167,6 +171,18 @@ public class OrthoMap implements Comparable <OrthoMap>, Serializable{
 //		pairwise_matches = new HashMap<Double, PairwiseOrthoMatch>();
 	}
 
+	double getEqualized(double d) {
+		return d * equalize[0] + equalize[1];
+	}
+	
+	double [] getEqualize() {
+		return equalize;
+	}
+	
+	void setEqualize(double [] equalize) {
+		this.equalize = equalize;
+	}
+	
     @Override
     public int compareTo(OrthoMap otherPlayer) {
         return Double.compare(ts, otherPlayer.ts);
@@ -1717,9 +1733,9 @@ public class OrthoMap implements Comparable <OrthoMap>, Serializable{
 			};
 		}		      
 		ImageDtt.startAndJoin(threads);
-		double [] ab = getDatiRegression(
-				temp, // final double []  temp,
-				dati, // final double []  dati,
+		double [] ab = PolynomialApproximation.getYXRegression(
+				temp, // final double []  data_x,
+				dati, // final double []  data_y,
 				flat); // final boolean [] mask);
 		double a = ab[0];
 		double b = ab[1];
@@ -1774,56 +1790,6 @@ public class OrthoMap implements Comparable <OrthoMap>, Serializable{
 		return flat;
 	}
 	
-	private static double [] getDatiRegression(
-			final double []  temp,
-			final double []  dati,
-			final boolean [] mask) {
-		final Thread[] threads = ImageDtt.newThreadArray(QuadCLT.THREADS_MAX);
-		final AtomicInteger ai = new AtomicInteger(0);
-		final double [] as0 =  new double[threads.length];
-		final double [] asx =  new double[threads.length];
-		final double [] asx2 = new double[threads.length];
-		final double [] asy =  new double[threads.length];
-		final double [] asxy=  new double[threads.length];
-		final AtomicInteger ati = new AtomicInteger(0);
-		ai.set(0);
-		for (int ithread = 0; ithread < threads.length; ithread++) {
-			threads[ithread] = new Thread() {
-				public void run() {
-					int thread_num = ati.getAndIncrement();
-					for (int ipix = ai.getAndIncrement(); ipix < temp.length; ipix = ai.getAndIncrement()) if (mask[ipix]){
-						double x = temp[ipix];
-						double y = dati[ipix];
-						if (!Double.isNaN(x+y)) {
-							as0 [thread_num] += 1;
-							asx [thread_num] += x;
-							asx2[thread_num] += x * x;
-							asy [thread_num] += y;
-							asxy[thread_num] += x * y;
-						}
-					} // for (int ipix
-				}
-			};
-		}		      
-		ImageDtt.startAndJoin(threads);
-		double s0 = 0.0;
-		double sx = 0.0;
-		double sx2 = 0.0;
-		double sy = 0.0;
-		double sxy= 0.0;
-		for (int i = 0; i < threads.length; i++) {
-			s0+=  as0[i];
-			sx+=  asx[i];
-			sx2+= asx2[i];
-			sy+=  asy[i];
-			sxy+= asxy[i];
-		}
-		double dnm = s0 * sx2 - sx*sx;
-		double a = (sxy * s0 - sy * sx) / dnm;
-		double b = (sy * sx2 - sxy * sx) / dnm;
-		return new double[] {a,b};
-	}
-	
 	
 	private static int [][] getCirclePoints(
 			double radius) {

src/main/java/com/elphel/imagej/orthomosaic/OrthoMultiLMA.java

@@ -41,6 +41,7 @@ import ij.gui.PointRoi;
 
 public class OrthoMultiLMA {
 	final   boolean           move_only;
+	final   boolean           corr_avg; // correct average skew, tilt, scale
 	private double []         last_rms =        null; // {rms, rms_pure}, matching this.vector
 	private double []         good_or_bad_rms = null; // just for diagnostics, to read last (failed) rms
 	private double []         initial_rms =     null; // {rms, rms_pure}, first-calcualted rms
@@ -61,7 +62,11 @@ public class OrthoMultiLMA {
 	private double [][]       offsets =         null; // scene offsets (rd)
 	private int               num_scenes =      0;
 	private int               num_pairs =       0;
-	public OrthoMultiLMA (boolean move_only) {
+	public OrthoMultiLMA (
+			boolean corr_avg,
+			boolean move_only
+			) {
+		this.corr_avg = corr_avg && !move_only;
 		this.move_only = move_only;
 	}
 	public static boolean testMultiLMA(
@@ -122,6 +127,8 @@ public class OrthoMultiLMA {
 				indices,    // int []        indices, // null or which indices to use (normally just 2 for pairwise comparison)
 				true,       // boolean       bounds_to_indices,
 				affines_a, // affines_0d, // affines,    // double [][][] affines, // null or [indices.length][2][3]
+				null,              // double [][]   equalize,
+				true,              // boolean       ignore_equalize,
 				null,       // warp,       // FineXYCorr    warp,,
 				zoom_level, // int         zoom_level,
 				wh,         // int []      wh,
@@ -147,6 +154,8 @@ public class OrthoMultiLMA {
 				indices,    // int []        indices, // null or which indices to use (normally just 2 for pairwise comparison)
 				true,       // boolean       bounds_to_indices,
 				affines_sym,// affines,    // double [][][] affines, // null or [indices.length][2][3]
+				null,              // double [][]   equalize,
+				true,              // boolean       ignore_equalize,
 				null,       // warp,       // FineXYCorr    warp,,
 				zoom_level, // int         zoom_level,
 				wh,         // int []      wh,
@@ -213,19 +222,26 @@ public class OrthoMultiLMA {
 		int                 debugLevel = 1;
 		boolean             move_only = false;
 		double []           val_coord = null; // 1 - valid, 0 - invalid, minimize coordinates errors
-		double              position_pull = 0.001;
         boolean             use_inv = false;
         double              overlap_pow = 2.0; // match weight as overlap fraction to this power
+        
+        double              skew_pull =     1.0;
+        double              tilt_pull =     1.0;
+        double              scale_pull =    0.1; // .0;
+		double              position_pull = 0.0001;
+		boolean             corr_avg= (skew_pull > 0) || (tilt_pull > 0) || (scale_pull > 0);
 		int [] indices = maps_collection.getScenesSelection(
 				null, // boolean select_all,
 				" to build a map"); // String purpose)
-		OrthoMultiLMA oml = new OrthoMultiLMA(move_only);
+		OrthoMultiLMA oml = new OrthoMultiLMA(
+				corr_avg,
+				move_only);
 		double lambda =            0.1;
 		double lambda_scale_good = 0.5;
 		double lambda_scale_bad =  8.0;
 		double lambda_max =      100;
 		double rms_diff =          0.000001;
-		int    num_iter =         20;
+		int    num_iter =        100; //  50;
 		boolean last_run =       false;
 		oml.prepareLMA (
 				clt_parameters,   // CLTParameters       clt_parameters,
@@ -233,6 +249,9 @@ public class OrthoMultiLMA {
 				indices,          // int []              indices,
 				val_coord,        // double []           val_coord, // 1 - valid, 0 - invalid, minimize coordinates errors
 				position_pull,    // double              position_pull,
+				skew_pull,        // double              skew_pull,
+				tilt_pull,        // double              tilt_pull,
+				scale_pull,       // double              scale_pull,
 				use_inv,          // boolean             use_inv,
 				overlap_pow,      // double              overlap_pow, // match weight as overlap fraction to this power
 				debugLevel);      // int                 debugLevel)
@@ -256,7 +275,7 @@ public class OrthoMultiLMA {
         //Get and apply affines
     	//oml.updateAffines(maps_collection);
         double [][][] affines = oml.getAffines();
-        
+        double [] fx = oml.getFx();
         // test
         /*
 		PairwiseOrthoMatch match = maps_collection.ortho_maps[indices[0]].getMatch(maps_collection.ortho_maps[indices[1]].getName());
@@ -284,6 +303,8 @@ public class OrthoMultiLMA {
 				indices,    // int []        indices, // null or which indices to use (normally just 2 for pairwise comparison)
 				true, // boolean       bounds_to_indices,
 				affines, // null, // affines,    // double [][][] affines, // null or [indices.length][2][3]
+				null,              // double [][]   equalize,
+				true,              // boolean       ignore_equalize,
 				null, // warp,       // FineXYCorr    warp,,
 				zoom_level, // int         zoom_level,
 				wh,         // int []      wh,
@@ -307,6 +328,7 @@ public class OrthoMultiLMA {
     	roi.setOptions("label");
     	imp_multi.setRoi(roi);
 		imp_multi.show();
+		/*
 		double [][] affine = 		getAffineAndDerivatives(
 				move_only,      //boolean       move_only,
 				affines[0],     // double [][]   affine00,
@@ -317,7 +339,7 @@ public class OrthoMultiLMA {
 				affines[0],
 				affines[1],
 				oml.offsets[0]);
-		
+		*/
 		
 		
 		/*
@@ -336,6 +358,16 @@ public class OrthoMultiLMA {
         return 0;
 	}
 	
+	public double [] getParameters() {
+		return parameters_vector;
+	}
+	
+	public double [] getFx() {
+		return getFxDerivs(
+				parameters_vector, // double []         vector,
+				null,              // final double [][] jt, // should be null or initialized with [vector.length][]
+				0);               // final int         debug_level)
+	}
 	
 	public int prepareLMA (
 			CLTParameters       clt_parameters,
@@ -343,6 +375,9 @@ public class OrthoMultiLMA {
 			int []              indices,
 			double []           val_coord, // 1 - valid, 0 - invalid, minimize coordinates errors
 			double              position_pull,
+	        double              skew_pull,
+	        double              tilt_pull,
+	        double              scale_pull,
             boolean             use_inv,
             double              overlap_pow, // match weight as overlap fraction to this power
 			int                 debugLevel) {
@@ -383,7 +418,6 @@ public class OrthoMultiLMA {
 		for (int np = 0; np < num_pairs; np++) {
 			Point pair = pairs_list.get(np);
 			pairs[np] = new int[] {pair.x, pair.y};
-//			double [] enuOffset = maps_collection.ortho_maps[indices[pair.x]].enuOffsetTo(maps_collection.ortho_maps[indices[pair.y]]);
 			double [] enuOffset = maps_collection.ortho_maps[indices[pair.y]].enuOffsetTo(maps_collection.ortho_maps[indices[pair.x]]);
 			offsets[np] = new double[] {enuOffset[0], -enuOffset[1]}; // {right,down} of the image 
 		}
@@ -396,9 +430,8 @@ public class OrthoMultiLMA {
 		int n62 = move_only ? 2 : 6;
 		int n13 = move_only ? 1 : 3;
 		
-		
-		y_vector =          new double [n62*num_pairs + 2 * num_scenes]; // last 2 * num_scenes will stay 0
-		weights =           new double [n62*num_pairs + 2 * num_scenes];
+		y_vector =          new double [n62*num_pairs + 2 * num_scenes + (corr_avg ? 3:0)]; // last 2 * num_scenes + 3 will stay 0
+		weights =           new double [n62*num_pairs + 2 * num_scenes + (corr_avg ? 3:0)];
 		parameters_vector = new double [n62*num_scenes]; // maybe will need move-only mode?
 		for (int n = 0; n < num_scenes; n++) {
 			double [][] affine = maps_collection.ortho_maps[indices[n]].getAffine();
@@ -442,6 +475,11 @@ public class OrthoMultiLMA {
 			weights[wi++] = w;
 			sw+= 2*w;
 		}
+		if (corr_avg) {
+			weights[wi++] = skew_pull;
+			weights[wi++] = tilt_pull;
+			weights[wi++] = scale_pull;
+		}
 		pure_weight = swp/sw;
 		double s = 1.0/sw;
 		for (int i = 0; i <weights.length;i++) {
@@ -1166,6 +1204,49 @@ public class OrthoMultiLMA {
 			};
 		}		      
 		ImageDtt.startAndJoin(threads);
+		ai.set(0);
+		if (corr_avg) { //  && ((skew_pull > 0) || (tilt_pull > 0) || (scale_pull > 0))) {
+			final int pull_indices = weights.length-3;
+			for (int nScene = 0; nScene < num_scenes; nScene++){
+				int indx = 6 * nScene;
+				double a00 = vector[indx+0],a01=vector[indx+1],a10=vector[indx+3],a11=vector[indx+4];
+				double aa1 = 0.5*(a00*a00 + a10*a10);
+				double aa2 = 0.5*(a01*a01 + a11*a11);
+				fx[pull_indices + 0] += a00*a01+a10*a11;
+				fx[pull_indices + 1] += aa1-aa2;
+				fx[pull_indices + 2] += aa1+aa2 - 1.0;
+			}
+			
+			if (jt!= null) {
+				for (int ithread = 0; ithread < threads.length; ithread++) {
+					threads[ithread] = new Thread() {
+						public void run() {
+							for (int nScene = ai.getAndIncrement(); nScene < num_scenes; nScene = ai.getAndIncrement()){
+								int indx = 6 * nScene;
+								double a00=vector[indx+0],a01=vector[indx+1],a10=vector[indx+3],a11=vector[indx+4];
+								jt[indx + 0][pull_indices + 0] = a01;
+								jt[indx + 1][pull_indices + 0] = a00;
+								jt[indx + 3][pull_indices + 0] = a11;
+								jt[indx + 4][pull_indices + 0] = a10;
+								
+								jt[indx + 0][pull_indices + 1] = a00;
+								jt[indx + 1][pull_indices + 1] =-a01;
+								jt[indx + 3][pull_indices + 1] = a10;
+								jt[indx + 4][pull_indices + 1] =-a11;
+								
+								jt[indx + 0][pull_indices + 2] = a00;
+								jt[indx + 1][pull_indices + 2] = a01;
+								jt[indx + 3][pull_indices + 2] = a10;
+								jt[indx + 4][pull_indices + 2] = a11;
+							}
+						}
+					};
+				}		      
+				ImageDtt.startAndJoin(threads);
+				ai.set(0);
+			}
+		}
+		
 		return fx;
 	}	
 

src/main/java/com/elphel/imagej/orthomosaic/PairwiseOrthoMatch.java

@@ -15,6 +15,8 @@ public class PairwiseOrthoMatch implements Serializable {
 	public double                rms = Double.NaN;
 	public transient int []      nxy = null; // not saved, just to communicate for logging
 	public transient double      overlap = 0.0;
+	public transient double []   equalize1to0 = {1,0}; // value1 = equalize2to1[0]*value2+equalize2to1[1]
+
 	public PairwiseOrthoMatch() {
 		
 	}
@@ -30,6 +32,11 @@ public class PairwiseOrthoMatch implements Serializable {
 		this.rms =     rms;
 		this.overlap = overlap;
 	}
+	
+	public void setEqualize2to1 (double [] equalize2to1) {
+		this.equalize1to0 = equalize2to1;
+	}
+	
 	public PairwiseOrthoMatch clone() {
 		double [][] affine = {this.affine[0].clone(),this.affine[1].clone()};
 		double [][] jtj = new double [this.jtj.length][];
@@ -45,6 +52,7 @@ public class PairwiseOrthoMatch implements Serializable {
 		if (nxy != null) {
 			pom.nxy = nxy.clone();
 		}
+		pom.equalize1to0 = this.equalize1to0.clone();
 		return pom;
 	}
 	/**
@@ -102,6 +110,7 @@ public class PairwiseOrthoMatch implements Serializable {
 				rd[0] * affine[1][0]+       rd[1]*(affine[1][1]-1.0)};
 		affine[0][2] += corr[0];
 		affine[1][2] += corr[1];
+		inverted_match.setEqualize2to1(new double [] {1/equalize1to0[0], -equalize1to0[1]/equalize1to0[0]}); 
 		return inverted_match;	
 	}
 	
@@ -168,10 +177,16 @@ public class PairwiseOrthoMatch implements Serializable {
 		affine = new double[][] {
 				{A.get(0,0),A.get(0,1), B.get(0,0)},
 				{A.get(1,0),A.get(1,1), B.get(1,0)}};
-//				jtj = null;
-//				rms = Double.NaN; // double rms,
-//				zoom_lev = 0; // int zoom_lev)
 	}
+	
+	public void combineEqualize(
+			double [] equalize0,
+			double [] equalize1	) {
+		setEqualize2to1(new double[] {
+				equalize1[0]/equalize0[0],
+				equalize1[1]-equalize1[0]/equalize0[0]*equalize0[1]});
+	}
+	
 	public static double [][] combineAffines( 
 			double [][] affine0,
 			double [][] affine, // differential
@@ -203,6 +218,8 @@ public class PairwiseOrthoMatch implements Serializable {
 			}
 		}
 		oos.writeObject(overlap);
+		oos.writeObject(equalize1to0);
+		
 	}
 	private void readObject(ObjectInputStream ois) throws ClassNotFoundException, IOException {
 		ois.defaultReadObject();
@@ -216,6 +233,8 @@ public class PairwiseOrthoMatch implements Serializable {
 			}
 		}
 		overlap = (Double) ois.readObject();
+//		equalize1to0 = new double[] {1,0};
+		equalize1to0 = (double[]) ois.readObject();
 	}
 	//private void readObjectNoData() throws ObjectStreamException; // used to modify default values
 }

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

@@ -1356,4 +1356,5 @@ public class TileNeibs{
 		}
 		return filled;
 	}
+	
 }