working snapshot, will next improve near-horizon targets

src/main/java/com/elphel/imagej/cuas/CuasMotion.java

@@ -2064,7 +2064,7 @@ public class CuasMotion {
 			String              dbg_suffix,
 			int                 dbg_iter_index,     // to generate different image names for different iteration to simplify "save as" images
 			int                 debugLevel) {
-		final int max_index_dbg = 5; // show images for the first iterations only
+		final int max_index_dbg = 3;// 5; //  // show images for the first iterations only
 		final int tileSize = GPUTileProcessor.DTT_SIZE;   // 8 pixels
 		final int tilesX   = cuasMotion.tilesX;
 		final int width    = cuasMotion.gpu_max_width;
@@ -2095,6 +2095,7 @@ public class CuasMotion {
 		final int r1i = (int) Math.ceil(recalc_mv_r1);
 		final int mside = 2 * r1i + 1;
 		final float[] mask_kernel = new float[mside * mside];
+		double sum_w = 0.0;
 		for (int dy = -r1i; dy <= r1i; dy++) {
 			for (int dx = -r1i; dx <= r1i; dx++) {
 				double r = Math.sqrt(dx * dx + dy * dy);
@@ -2103,8 +2104,14 @@ public class CuasMotion {
 				else if (r >= recalc_mv_r1) w = 0.0f;
 				else w = (float)(0.5 * (Math.cos(Math.PI * (r - recalc_mv_r0) / (recalc_mv_r1 - recalc_mv_r0)) + 1.0));
 				mask_kernel[(dy + r1i) * mside + (dx + r1i)] = w;
+				sum_w+=w;
 			}
 		}
+		double scale_w = 4.0*GPUTileProcessor.DTT_SIZE*GPUTileProcessor.DTT_SIZE/sum_w;
+		// maybe scale: cosine integral divided by mask is multiplied by the cosine 
+		for (int i = 0; i < mask_kernel.length; i++) {
+			mask_kernel[i] *= scale_w; 
+		}
 		if ((dbg_nseq >= 0 || dbg_tile >= 0) && (debugLevel >= 0) && (dbg_iter_index<= max_index_dbg)) {
 			ShowDoubleFloatArrays.showArrays(mask_kernel.clone(), mside, mside,
 					"refineMotionVectors-mask-boost"+recalc_mv_boost+"-r0_" + recalc_mv_r0 + "-r1_" + recalc_mv_r1+"-niter"+dbg_iter_index+dbg_suffix);

src/main/java/com/elphel/imagej/cuas/CuasMotionLMA.java

@@ -223,6 +223,75 @@ public class CuasMotionLMA {
 		double y0 = Math.min(Math.max(yc + width/2, 0), width-1);
 		
 
+		int ix0 = (int) Math.round(x0);
+		int iy0 = (int) Math.round(y0);
+		full_vector[INDX_A] = tile_data[ix0+iy0*width];
+		full_vector[INDX_C] = 0;
+		full_vector[INDX_RR0] = Math.PI/(2* r0);
+//		full_vector[INDX_K] = k; // > 1 (~2.0)
+		full_vector[INDX_K] = Math.sqrt(Math.abs(k-1)); // > 1 (~2.0) now k = 1 + k2*k2;
+		full_vector[INDX_X0] = x0;
+		full_vector[INDX_Y0] = y0;
+		weights = new double [width*width];
+		double sw = 0;
+		for (int y = 0; y < width; y++) {
+			int ay = Math.abs(y-iy0);
+			for (int x = 0; x < width; x++) {
+				int ax = Math.abs(x-ix0);
+				double w = window[ay][ax]+pedestal; // window to the nearest integer x,y
+				weights[x + y*width] = w;
+				sw += w;
+			}
+		}
+		for (int i = 0; i < weights.length; i++) {
+			weights[i] /= sw;
+		}
+		int indx = 0;
+		for (int i = 0; i < INDX_LEN; i++) if (param_select[i]){
+			indx++;
+		}
+		rindx = new int [INDX_LEN];
+		pindx = new int [indx];
+		indx= 0;
+		for (int i = 0; i < INDX_LEN; i++) {
+			if (param_select[i]) {
+				pindx[indx] = i;
+				rindx[i] = indx++;
+			} else {
+				rindx[i] = -1;
+			}
+		}
+   		last_jt = new double [pindx.length][];
+		
+		double [] fx = getFxDerivs(
+				getParametersVector(), // double []         vector,
+				last_jt,               // final double [][] jt, // should be null or initialized with [vector.length][]
+				debugLevel);           // final int         debug_level)
+		last_rms =        new double [2];
+		last_ymfx = getYminusFxWeighted(
+				fx, // final double []   fx,
+				last_rms); // final double []   rms_fp // null or [2]
+		
+		initial_rms =     last_rms.clone();
+		good_or_bad_rms = last_rms.clone();
+		return 0;
+	}	
+
+	public int prepareLMA_orig(
+			boolean [] param_select,
+			double []  tile_data,
+			double     xc, // relative to center =width/2
+			double     yc, // relative to center =width/2
+			double     r0,
+			double     k,
+			double     lmax_val, // maximal pixel value near the centroid maximum to be used for comparison with A
+			int        debugLevel) {
+		max_val = lmax_val;
+		y_vector = tile_data;
+		double x0 = Math.min(Math.max(xc + width/2, 0), width-1);
+		double y0 = Math.min(Math.max(yc + width/2, 0), width-1);
+		
+
 		int ix0 = (int) Math.round(x0);
 		int iy0 = (int) Math.round(y0);
 		full_vector[INDX_A] = tile_data[ix0+iy0*width];
@@ -328,6 +397,10 @@ public class CuasMotionLMA {
 	}
 	
 	public double getK() {
+		return 1 + full_vector[INDX_K]*full_vector[INDX_K];
+	}
+
+	public double getK_orig() {
 		return full_vector[INDX_K];
 	}
 	
@@ -795,6 +868,81 @@ public class CuasMotionLMA {
 		double A=    (rindx[INDX_A] >=0)?   vector[rindx[INDX_A]] :   full_vector[INDX_A] ; 
 		double C=    (rindx[INDX_C] >=0)?   vector[rindx[INDX_C]] :   full_vector[INDX_C] ;
 		double RR0 = (rindx[INDX_RR0] >=0)? vector[rindx[INDX_RR0]] : full_vector[INDX_RR0] ;
+//		double K =   (rindx[INDX_K] >=0)?   vector[rindx[INDX_K]] :   full_vector[INDX_K] ; // >1.0
+		double K2 =  (rindx[INDX_K] >=0)?   vector[rindx[INDX_K]] :   full_vector[INDX_K] ; // >1.0
+		double K =   1+ K2*K2; // always >=1
+		double dK_dK2 = 2 * K2;
+		double X0 =  (rindx[INDX_X0] >=0)?  vector[rindx[INDX_X0]] :  full_vector[INDX_X0] ;
+		double Y0 =  (rindx[INDX_Y0] >=0)?  vector[rindx[INDX_Y0]] :  full_vector[INDX_Y0] ;
+		for (int y = 0; y < height; y++) {
+			double dy = y-Y0;
+			double dy2 = dy*dy;
+			for (int x = 0; x < width; x++) {
+				int indx = y*width+x;
+				double dx = x-X0;
+				double dx2 = dx*dx;
+				double r2 = dx2+dy2+R_OFFS;
+				double r = Math.sqrt(r2);
+				if (r*RR0 > Math.PI/2) { // outside circle where function is defined
+					fx[indx] = C;
+					if (jt != null) {
+						double df_dC = 1;
+						if (rindx[INDX_C] >=0) {
+							jt[rindx[INDX_C]][indx] = df_dC;
+						}
+					}					
+				} else {
+					double W1 = Math.cos(r*RR0);
+					double W2 = Math.cos(r*RR0*K); //Math.cos(r*RR0*(1+K2*K2))
+					fx[indx] = A * (W1 * W1 * W2) + C;
+					if (jt != null) {
+						double dr2_dX0 = -2 * dx;
+						double dr2_dY0 = -2 * dy;
+						double dr_dr2 = 0.5 / r; // r !=0
+						double sin_rRR0 = Math.sin(r * RR0);
+						double dW1_dRR0 = -r * sin_rRR0;     // 
+						double dW1_dr = -RR0 *sin_rRR0;
+						double sin_rKRR0 = Math.sin(r * RR0 * K);
+						double dW2_dRR0 =  -K * r * sin_rKRR0;
+//						double dW2_dK =    -r * RR0 * sin_rKRR0;
+						double dW2_dK2 =  (-r * RR0 * sin_rKRR0) * dK_dK2;  
+						double dW2_dr =    -K * RR0 * sin_rKRR0; 
+						double df_dA = W1*W1*W2; // +
+						double df_dC = 1;        // +  
+						double df_dW1 = 2 * A* W1 * W2;
+						double df_dW2 = A * W1 * W1;
+//						double df_dK =  df_dW2 * dW2_dK; 
+						double df_dK2 = df_dW2 * dW2_dK2; 
+						double df_dRR0 = df_dW1*dW1_dRR0+df_dW2*dW2_dRR0;
+						double df_dr = df_dW1 * dW1_dr + df_dW2 * dW2_dr;  
+						double df_dr2 = df_dr * dr_dr2;
+						double df_dX0 = df_dr2 * dr2_dX0;
+						double df_dY0 = df_dr2 * dr2_dY0;
+//						double [] df_dfp = {df_dA,df_dC,df_dRR0, df_dK, df_dX0, df_dY0};
+						double [] df_dfp = {df_dA,df_dC,df_dRR0, df_dK2, df_dX0, df_dY0};
+						for (int i = 0; i < vector.length; i++) {
+							jt[i][indx] = df_dfp[pindx[i]]; 
+						}
+					}
+				}
+			}
+		}
+		return fx;
+	}
+	private double [] getFxDerivs_orig(
+			final double []   vector,
+			final double [][] jt, // should be null or initialized with [vector.length][]
+			final int         debug_level) {
+		int height = y_vector.length/width;
+		final double [] fx = new double [y_vector.length];
+		if (jt != null) {
+			for (int i = 0; i < jt.length; i++) {
+				jt[i] = new double [y_vector.length]; // weights.length];
+			}
+		}
+		double A=    (rindx[INDX_A] >=0)?   vector[rindx[INDX_A]] :   full_vector[INDX_A] ; 
+		double C=    (rindx[INDX_C] >=0)?   vector[rindx[INDX_C]] :   full_vector[INDX_C] ;
+		double RR0 = (rindx[INDX_RR0] >=0)? vector[rindx[INDX_RR0]] : full_vector[INDX_RR0] ;
 		double K =   (rindx[INDX_K] >=0)?   vector[rindx[INDX_K]] :   full_vector[INDX_K] ; // >1.0
 		double X0 =  (rindx[INDX_X0] >=0)?  vector[rindx[INDX_X0]] :  full_vector[INDX_X0] ;
 		double Y0 =  (rindx[INDX_Y0] >=0)?  vector[rindx[INDX_Y0]] :  full_vector[INDX_Y0] ;

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

@@ -822,8 +822,8 @@ min_str_neib_fpn	0.35
     public int      cuas_recalc_mv_num =   2;     // Number of recalculations of the motion vectors before centered targets accumulation by masking far-from target areas
     public double   cuas_recalc_mv_boost = 4.0;   // Scale default number of correlation pairs for motion vectors calculation
     public double   cuas_recalc_mv_corr =  4.0;   // Scale corr_offset for refinement pass (will use (int)Math.round() 
-    public double   cuas_recalc_mv_r0 =    2.0;   // Masking window parameters: for r <= r0, w = 1.0, first (coarse) pass
-    public double   cuas_recalc_mv_r1 =    6.0;   // Masking window parameters: for r >= r1, w = 0.0, r0<r<r1: w = 0.5*(cos(PI*(r-r0)/(r1-r0))+1), first (coarse) pass
+    public double   cuas_recalc_mv_r0 =    4.0;   // Masking window parameters: for r <= r0, w = 1.0, first (coarse) pass
+    public double   cuas_recalc_mv_r1 =    7.0;   // Masking window parameters: for r >= r1, w = 0.0, r0<r<r1: w = 0.5*(cos(PI*(r-r0)/(r1-r0))+1), first (coarse) pass
     public double   cuas_recalc_mv_r0f =   1.5;   // Masking window parameters: for r <= r0, w = 1.0, second (narrow) pass
     public double   cuas_recalc_mv_r1f =   4.0;   // Masking window parameters: for r >= r1, w = 0.0, r0<r<r1: w = 0.5*(cos(PI*(r-r0)/(r1-r0))+1), second (narrow) pass
     public double   cuas_recalc_mv_max2 =  0.2;   // Maximal Vx,Vy corection for the fine pass (only if cuas_recalc_mv_num >= 2)