Working snapshot, will try to improve matching scenes from different

sequences

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

@@ -793,7 +793,10 @@ min_str_neib_fpn	0.35
 	public double    terr_alpha_init_offs =    0.01; // Initial alpha: opaque/transparent offset from 1.0/0.0 
     
 	public double    terr_alpha_loss =       100.0;
-	public double    terr_alpha_offset =       0.0;
+	public double    terr_alpha_offset =       0.1;
+	public double    terr_alpha_min_veg =      -1; // old version 0.5;  // Minimal vegetation alpha. If (alpha-alpha_offset)/(1-2*alpha_offset) < alpha_min_veg, pull down to lpha_offset
+	
+	
 	public double    terr_alpha_lpf =          2.5;  // pull to average of 4 neighbors
 	public boolean   terr_alpha_piece_linear = true;
 	public double    terr_alpha_scale_avg =    1.0;  // scale average alpha (around 0.5) when pulling to it
@@ -2102,6 +2105,8 @@ min_str_neib_fpn	0.35
 		gd.addNumericField("Alpha loss",           terr_alpha_loss, 5,7,"", 	"Alpha quadratic growing loss for when out of [0,1] range");
 		
 		gd.addNumericField("Alpha offset",         terr_alpha_offset, 5,7,"",   "Start alpha losses above 0.0 and below 1.0 by this value.");
+		gd.addNumericField("Minimal vegetation alpha", terr_alpha_min_veg, 5,7,"","Minimal vegetation alpha. If (alpha-alpha_offset)/(1-2*alpha_offset) < alpha_min_veg, pull down to lpha_offset.");
+		
 		gd.addNumericField("Alpha diffusion",      terr_alpha_lpf, 5,7,"", 	    "Alpha diffusion to 4 ortho neighbors.");
 		gd.addCheckbox    ("Alpha piece-linear",   terr_alpha_piece_linear,     "Piece-linear alpha (_/\u203E, false - 0.0-cosine-1.0.");
 		gd.addNumericField("Scale alpha",          terr_alpha_scale_avg, 5,7,"","Scale target (average of neighbors) alpha before pulling to it (not used now).");
@@ -2818,7 +2823,6 @@ min_str_neib_fpn	0.35
 		terr_fit_alpha =           gd.getNextBoolean();// boolean
 		terr_fit_scenes =          gd.getNextBoolean();// boolean
 		terr_fit_elevations =      gd.getNextBoolean();// boolean
-		
 		for (int i = 0; i < terr_fit_disable.length; i++) {
 			terr_fit_disable[i] =      gd.getNextBoolean();// boolean
 		}
@@ -2847,6 +2851,7 @@ min_str_neib_fpn	0.35
 		
 		terr_alpha_loss =          gd.getNextNumber();// double 
 		terr_alpha_offset =        gd.getNextNumber();// double 
+		terr_alpha_min_veg =       gd.getNextNumber();// double
 		terr_alpha_lpf =           gd.getNextNumber();// double 
 		terr_alpha_piece_linear =  gd.getNextBoolean();// boolean
 		terr_alpha_scale_avg =     gd.getNextNumber();// double 
@@ -3557,6 +3562,7 @@ min_str_neib_fpn	0.35
 		
 		properties.setProperty(prefix+"terr_alpha_loss",             terr_alpha_loss+"");         // double 
 		properties.setProperty(prefix+"terr_alpha_offset",           terr_alpha_offset+"");       // double 
+		properties.setProperty(prefix+"terr_alpha_min_veg",          terr_alpha_min_veg+"");      // double 
 		properties.setProperty(prefix+"terr_alpha_lpf",              terr_alpha_lpf+"");          // double 
 		properties.setProperty(prefix+"terr_alpha_piece_linear",     terr_alpha_piece_linear+""); // boolean
 		properties.setProperty(prefix+"terr_alpha_scale_avg",        terr_alpha_scale_avg+"");    // double 
@@ -4284,6 +4290,7 @@ min_str_neib_fpn	0.35
 		
 		if (properties.getProperty(prefix+"terr_alpha_loss")!=         null) terr_alpha_loss=Double.parseDouble(properties.getProperty(prefix+"terr_alpha_loss"));
 		if (properties.getProperty(prefix+"terr_alpha_offset")!=       null) terr_alpha_offset=Double.parseDouble(properties.getProperty(prefix+"terr_alpha_offset"));
+		if (properties.getProperty(prefix+"terr_alpha_min_veg")!=       null) terr_alpha_min_veg=Double.parseDouble(properties.getProperty(prefix+"terr_alpha_min_veg"));
 		if (properties.getProperty(prefix+"terr_alpha_lpf")!=          null) terr_alpha_lpf=Double.parseDouble(properties.getProperty(prefix+"terr_alpha_lpf"));
 		if (properties.getProperty(prefix+"terr_alpha_piece_linear")!= null) terr_alpha_piece_linear=Boolean.parseBoolean(properties.getProperty(prefix+"terr_alpha_piece_linear"));
 		if (properties.getProperty(prefix+"terr_alpha_scale_avg")!=    null) terr_alpha_scale_avg=Double.parseDouble(properties.getProperty(prefix+"terr_alpha_scale_avg"));
@@ -4976,7 +4983,8 @@ min_str_neib_fpn	0.35
 		imp.terr_alpha_init_offs =               this.terr_alpha_init_offs;
 		
 		imp.terr_alpha_loss =                    this.terr_alpha_loss;                           
-		imp.terr_alpha_offset =                  this.terr_alpha_offset;                           
+		imp.terr_alpha_offset =                  this.terr_alpha_offset;
+		imp.terr_alpha_min_veg =                 this.terr_alpha_min_veg;
 		imp.terr_alpha_lpf =                     this.terr_alpha_lpf;                           
 		imp.terr_alpha_piece_linear =            this.terr_alpha_piece_linear;                           
 		imp.terr_alpha_scale_avg =               this.terr_alpha_scale_avg;                           

src/main/java/com/elphel/imagej/vegetation/VegetationLMA.java

@@ -61,6 +61,10 @@ public class VegetationLMA {
 	public static final int SAMPLES_SIZE =  SAMPLES_TOTAL+1;
 	public static final int SAMPLES_EXTRA = SAMPLES_SCENES; // start of "extra" samples
 	
+	public final static String[] SAMPLES_NAMES = {"SAMPLES_Y", "SAMPLES_Y_HF", "SAMPLES_Y_AVG",
+			"SAMPLES_SCENES","SAMPLES_ALPHA_PULL","SAMPLES_ALPHA_LPF","SAMPLES_TERRAIN_PULL",
+			"SAMPLES_TERRAIN_LPF","SAMPLES_VEGETATION_PULL","SAMPLES_VEGETATION_LPF",
+			"SAMPLES_ELEVATION_PULL","SAMPLES_ELEVATION_LPF","SAMPLES_ELEVATION_ALPHA","SAMPLES_TOTAL"};	
 	public static final int [] SAVE_TYPES = {TVAO_TERRAIN, TVAO_VEGETATION, TVAO_ALPHA, TVAO_ELEVATION, TVAO_SCENE_OFFSET};
 
 	private final Rectangle    full;
@@ -153,6 +157,7 @@ public class VegetationLMA {
 	private boolean []         fits;                      // (+) 
 	private double             alpha_loss =         0;    // (+) not used with cosine alpha
 	private double             alpha_offset =       0;    // (+) if >0,  start losses above 0.0 and below 1.0;
+	private double             alpha_min_veg =      0.5;  // (+) if (alpha-alpha_offset)/(1-2*alpha_offset) < alpha_min_veg, pull down to lpha_offset
 	private double             alpha_lpf =          0;    // (+) 
 	private boolean            alpha_piece_linear = true; // (+) 
 	private double             alpha_scale_avg =    1.0;  // (+) 1.1; // scale average alpha (around 0.5) when pulling to it
@@ -249,7 +254,7 @@ public class VegetationLMA {
 		image_length = full.width * full.height;
 		num_scenes = vegetationModel.terrain_scenes_render.length;
 		vegetation_average =vegetationModel.vegetation_full;      // vegetation_average_render;
-		vegetation_pull =   vegetationModel.vegetation_pull;      // pull LMA to this level
+//		vegetation_pull =   vegetationModel.vegetation_pull;      // pull LMA to this level
 		terrain_average =   vegetationModel.terrain_filled;       // terrain filled under vegetation, use to pull terrain level to it
 		terrain_rendered =  vegetationModel.terrain_scenes_render;  
 		vegetation_offsets =vegetationModel.vegetation_warp;
@@ -274,6 +279,8 @@ public class VegetationLMA {
 		scene_weights = new double [num_scenes];
 		scales_xy = getScalesXY(vegetationModel.scale_dirs);
 		tvao[TVAO_ELEVATION] =  vegetationModel.elevations.clone();
+		// 		vegetation_pull =   vegetationModel.vegetation_pull;      // pull LMA to this level
+		vegetation_pull =tvao[TVAO_VEGETATION].clone();
 		if (debugLevel > -2) {
 			String [] titles = {"TVAO_TERRAIN", "TVAO_VEGETATION", "vegetation_pull", "vegetation_full", "elevations", "TVAO_ALPHA"};
 			double [][] dbg_img = {tvao[TVAO_TERRAIN],tvao[TVAO_VEGETATION],vegetationModel.vegetation_pull,vegetationModel.vegetation_full,vegetationModel.elevations,tvao[TVAO_ALPHA]};
@@ -516,6 +523,7 @@ public class VegetationLMA {
 			final double             reg_weights,    // fraction of the total weight used for regularization
 			final double             alpha_loss,     // quadratic loss when alpha reaches -1.0 or 2.0 
 			final double             alpha_offset,   // quadratic loss when alpha reaches -1.0 or 2.0
+			final double             alpha_min_veg, // 0.5; // if (alpha-alpha_offset)/(1-2*alpha_offset) < alpha_min_veg, pull down to lpha_offset
 			final double             alpha_lpf,      // pull vegetation alpha to average of 4 neighbors
 			final boolean            alpha_piece_linear, // true - piece-linear, false - half-cosine
 			final double             alpha_scale_avg, //  = 1.2; // scale average alpha (around 0.5) when pulling to it
@@ -554,6 +562,7 @@ public class VegetationLMA {
 		this.reg_weights =       reg_weights;    // fraction of the total weight used for regularization
 		this.alpha_loss =        alpha_loss;
 		this.alpha_offset =      alpha_offset;
+		this.alpha_min_veg =     alpha_min_veg;
 		this.alpha_lpf =         alpha_lpf;
 		this.alpha_scale_avg =   alpha_scale_avg;
 		this.alpha_piece_linear= alpha_piece_linear;
@@ -923,6 +932,71 @@ public class VegetationLMA {
 		}
 		return wymfw;
 	}
+
+	public double[][] yMinusFxSummary(
+			final double []   fx,
+			final boolean show) {
+		final double [][] stats = new double [samples_pointers.length][2];
+		final Thread[]      threads =     ImageDtt.newThreadArray();
+		final AtomicInteger ai =          new AtomicInteger(0);
+		final AtomicInteger ati =         new AtomicInteger(0);
+		for (int n = 0; n < SAMPLES_TOTAL; n++) {
+			final int nstat = n;
+			final int nstat_start =  samples_pointers[nstat][0];
+			final int nstat_length = samples_pointers[nstat][1];
+			if (nstat_length > 0) {
+				ai.set(0);
+				ati.set(0);
+				final double [][] stats_threaded = new double [threads.length][2];
+				for (int ithread = 0; ithread < threads.length; ithread++) { // first sum for pairs
+					threads[ithread] = new Thread() {
+						public void run() {
+							int nthread =  ati.getAndIncrement();
+							for (int n = ai.getAndIncrement(); n < nstat_length; n = ai.getAndIncrement()) {
+								int indx = nstat_start + n;
+								double y = (indx < y_vector.length) ? y_vector[indx] : 0;
+								double d = y - fx[indx]; // - fx[n]; // +y_vector[i]
+								double w = weights[indx];
+								double wd = d * w;
+								if (Double.isNaN(wd)) {
+									System.out.println("wd=NaN, n="+n);
+								}
+								stats_threaded[nthread][0] += wd * d;
+								stats_threaded[nthread][1] += w;
+							}
+						}
+					};
+				}		      
+				ImageDtt.startAndJoin(threads);
+				for (int nthread = 0; nthread < stats_threaded.length; nthread++) {
+					stats[nstat][0] += stats_threaded[nthread][0];
+					stats[nstat][1] += stats_threaded[nthread][1];
+				}
+			}
+		}
+		for (int nstat = 0; nstat < SAMPLES_TOTAL; nstat++)  {
+			stats[SAMPLES_TOTAL][0] += stats[nstat][0];
+			stats[SAMPLES_TOTAL][1] += stats[nstat][1];
+		}
+		
+		if (show) {
+			System.out.println(String.format("%2s, %24s, %12s, %12s, %12s, %7s, %7s, %7s", "#", "name", "RMS_avg", "RMS_contrib", "weight", "start","end","length"));
+			for (int nstat = 0; nstat < stats.length; nstat++) if (stats[nstat][1] > 0) {
+				int istart =   samples_pointers[nstat][0];
+				int ilength =  samples_pointers[nstat][1];
+				int iend = istart+ilength-1;
+				double rms_avg =     Math.sqrt(stats[nstat][0]/stats[nstat][1]);
+				double rms_contrib = Math.sqrt(stats[nstat][0] * stats[nstat][1]);
+				double weight = stats[nstat][1];
+				System.out.println(String.format("%2d, %24s, %12.8f, %12.8f, %12.8f, %7d, %7d, %7d",
+						nstat, SAMPLES_NAMES[nstat],
+						rms_avg, rms_contrib, weight,
+						istart, iend, ilength));
+			}			
+		}
+		return stats;
+	}
+	
 	
 	public int runLma( // <0 - failed, >=0 iteration number (1 - immediately)
 			double lambda,           // 0.1
@@ -1020,7 +1094,15 @@ public class VegetationLMA {
 			}
 			System.out.println();
 		}
-
+		if (debug_level > -2) {
+			double [] fx = getFxDerivs(
+					parameters_vector, // double []         vector,
+					null,           // final double [][] jt, // should be null or initialized with [vector.length][]
+					debug_level);      // final int         debug_level)
+			yMinusFxSummary(     
+					fx, // final double []   fx,
+					true); // final boolean show)
+		}
 		return rslt[0]? iter : -1;
 	}
 
@@ -1044,7 +1126,12 @@ public class VegetationLMA {
 					last_rms); // final double []   rms_fp // null or [2]
 			this.initial_rms = this.last_rms.clone();
 			this.good_or_bad_rms = this.last_rms.clone();
-
+			if (debug_level > 0) {
+				yMinusFxSummary(     
+						fx, // final double []   fx,
+						true); // final boolean show)
+			}
+		
 			if (last_ymfx == null) {
 				return null; // need to re-init/restart LMA
 			}
@@ -1175,11 +1262,10 @@ public class VegetationLMA {
 //				vector_XYS, // final double [][] vector_XYS,
 				fx, // final double []   fx,
 				rms); // final double []   rms_fp // null or [2]
-		if (debug_level > 2) {
-			/*
-			dbgYminusFx(this.last_ymfx, "next y-fX");
-			dbgXY(new_vector, "XY-correction");
-			*/
+		if (debug_level > 0) {
+			yMinusFxSummary(     
+					fx, // final double []   fx,
+					true); // final boolean show)
 		}
 
 		if (last_ymfx == null) {
@@ -1223,6 +1309,11 @@ public class VegetationLMA {
 			last_ymfx = getYminusFxWeighted(
 					fx, // final double []   fx,
 					this.last_rms); // final double []   rms_fp // null or [2]
+			if (debug_level > 0) {
+				yMinusFxSummary(     
+						fx, // final double []   fx,
+						true); // final boolean show)
+			}
 			if (last_ymfx == null) {
 				return null; // need to re-init/restart LMA
 			}
@@ -1274,7 +1365,7 @@ public class VegetationLMA {
 		// hifreq_weight
 		String  debug_title = "parameters_vector-x"+woi.x+"-y"+woi.y+"-w"+woi.width+"-h"+woi.height;
 		debug_title += "-hf"+hifreq_weight+"-tc"+terrain_correction;
-		debug_title += "-al"+alpha_loss+"-alo"+alpha_offset+"-alp"+alpha_lpf+"-als"+alpha_scale_avg+(alpha_piece_linear?"-alin":"-acos");
+		debug_title += "-al"+alpha_loss+"-alo"+alpha_offset+"-amv"+alpha_min_veg+"-alp"+alpha_lpf+"-als"+alpha_scale_avg+(alpha_piece_linear?"-alin":"-acos");
 		debug_title += "-ap"+alpha_push+"-apn"+alpha_push_neutral+"-apc"+alpha_push_center+(alpha_en_holes?("-amm"+alpha_mm_hole):"");
 		debug_title += "-tl"+terr_lpf+"-vl"+veget_lpf+"-el"+elevation_lpf+"-tp"+terr_pull0+"-vp"+veget_pull0+"-ep"+elevation_pull0;
 		if (elev_alpha_en) {
@@ -2726,6 +2817,8 @@ public class VegetationLMA {
 		// No need here to add contributors as this should not depend on any other parameters 
 		if (samples_pointers[SAMPLES_ALPHA_PULL][1] > 0) { // fits[TVAO_ALPHA] &&  ((alpha_loss > 0) || (alpha_push > 0))) {
 			int dbg_nx = -76340;
+			final double alpha_threshold = alpha_offset + alpha_min_veg *(1 - 2*alpha_offset);
+			final double alpha_plateau = alpha_threshold * alpha_loss;
 			ai.set(0);
 			for (int ithread = 0; ithread < threads.length; ithread++) {
 				threads[ithread] = new Thread() {
@@ -2740,16 +2833,33 @@ public class VegetationLMA {
 							fX[nx] = 0.0;
 							if (alpha_loss > 0) {
 								double alpha = vector[np];
-								if (alpha < alpha_offset) {
-									d = alpha- alpha_offset;
-								} else if (alpha > (1 - alpha_offset)) {
-									d = alpha - (1.0 - alpha_offset);
-								}
-								if (d != 0) {
-									fX[nx] = d * d * alpha_loss;
-									if (jt != null) {
-										jt[np][nx] = 2 * alpha_loss * d; // d/dalpha[i]		
-									}							
+								if (alpha_min_veg >= 0) {
+									if (alpha < alpha_threshold) {
+										fX[nx] = (alpha - alpha_offset) * alpha_loss;
+										if (jt != null) {
+											jt[np][nx] = alpha_loss; // d/dalpha[i]		
+										}							
+									} else if (alpha > (1 - alpha_offset)) {
+										fX[nx] = (alpha - (1.0 - alpha_offset)) * alpha_loss + alpha_plateau;
+										if (jt != null) {
+											jt[np][nx] = alpha_loss; // d/dalpha[i]		
+										}							
+									} else {
+										fX[nx] = alpha_plateau;
+//										jt[np][nx] = 0.0; // d/dalpha[i]		
+									}
+								} else { // old version, requires alpha_min_veg < 0 
+									if (alpha < alpha_offset) {
+										d = alpha- alpha_offset;
+									} else if (alpha > (1 - alpha_offset)) {
+										d = alpha - (1.0 - alpha_offset);
+									}
+									if (d != 0) {
+										fX[nx] = d * d * alpha_loss;
+										if (jt != null) {
+											jt[np][nx] = 2 * alpha_loss * d; // d/dalpha[i]		
+										}							
+									}
 								}
 							}
 							if (alpha_push > 0) {
@@ -3862,6 +3972,7 @@ public class VegetationLMA {
 		imp.setProperty("FIT_SCENE_OFFSET",          ""+fits[TVAO_SCENE_OFFSET]);
 		imp.setProperty("ALPHA_LOSS",                ""+alpha_loss);
 		imp.setProperty("ALPHA_OFFSET",              ""+alpha_offset);
+		imp.setProperty("ALPHA_MIN_VEG",             ""+alpha_min_veg);
 		imp.setProperty("ALPHA_LPF",                 ""+alpha_lpf);
 		imp.setProperty("ALPHA_PIECE_LINEAR",        ""+alpha_piece_linear);
 		imp.setProperty("ALPHA_SCALE_AVG",           ""+alpha_scale_avg);
@@ -4028,6 +4139,7 @@ public class VegetationLMA {
 				fits[TVAO_SCENE_OFFSET] =   Boolean.parseBoolean((String) imp_pars.getProperty("FIT_SCENE_OFFSET"));
 				alpha_loss =                Double.parseDouble((String)   imp_pars.getProperty("ALPHA_LOSS"));
 				alpha_offset =              Double.parseDouble((String)   imp_pars.getProperty("ALPHA_OFFSET"));
+				if (imp_pars.getProperty("ALPHA_MIN_VEG") != null)  alpha_min_veg = Double.parseDouble((String)  imp_pars.getProperty("ALPHA_MIN_VEG"));
 				alpha_lpf =                 Double.parseDouble((String)   imp_pars.getProperty("ALPHA_LPF"));
 				alpha_piece_linear =        Boolean.parseBoolean((String) imp_pars.getProperty("ALPHA_PIECE_LINEAR"));
 				alpha_scale_avg =           Double.parseDouble((String)   imp_pars.getProperty("ALPHA_SCALE_AVG"));
@@ -5211,7 +5323,6 @@ public class VegetationLMA {
 								System.out.println("y_vector["+(wIndx+terr_start)+"]=NaN");
 							}
 							y_vector[wIndx+terr_start] = swy; // /sw; (fX uses the same weights)
-							
 						}
 					}
 				};
@@ -6380,8 +6491,8 @@ public class VegetationLMA {
 			threads[ithread] = new Thread() {
 				public void run() {
 					for (int indx_decim = ai.getAndIncrement(); indx_decim < data_decimated.length; indx_decim = ai.getAndIncrement()) {
-						int indx = par_map[0][indx_decim];
-						data_full[indx_decim] = (indx >= 0) ? data_decimated[indx] : 0.0;
+						int indx = par_map[1][indx_decim];
+						data_full[indx] = (indx >= 0) ? data_decimated[indx_decim] : 0.0;
 					}
 				}
 			};

src/main/java/com/elphel/imagej/vegetation/VegetationModel.java

@@ -1428,6 +1428,8 @@ public class VegetationModel {
 		double    default_alpha =     clt_parameters.imp.terr_alpha_dflt;          //  0.5; //  0.8;
 		double    alpha_loss =        clt_parameters.imp.terr_alpha_loss;          //100.0; // 10.0; /// 100.0; // 10.0;  // 10000.0; // 1000.0; // 100.; // 10.0;         // quadratic loss when alpha reaches -1.0 or 2.0 
 		double    alpha_offset =      clt_parameters.imp.terr_alpha_offset;        //  0.0; // 0.02; // 0.03;    // if >0,  start losses above 0.0 and below 1.0;
+		double    alpha_min_veg =     clt_parameters.imp.terr_alpha_min_veg;       //  0.5  // Minimal vegetation alpha. If (alpha-alpha_offset)/(1-2*alpha_offset) < alpha_min_veg, pull down to lpha_offset
+		
 		double    alpha_lpf =         clt_parameters.imp.terr_alpha_lpf;           //  2.5; // 5; /// 2; // 5; /// 10; /// 15; // 10.0; // 5.0; // 10.0; // 3; // 10;   // 20; // 6.0; // 3.0; // 2.0;    // 1.5; // 5.0; // 0.5;         // pull to average of 4 neighbors
 		boolean   alpha_piece_linear =clt_parameters.imp.terr_alpha_piece_linear;  // true; // false; // true;
 		double    alpha_scale_avg =   clt_parameters.imp.terr_alpha_scale_avg;     //  1.0; // 1.1; // 0.9; // 2.0; // 1.5; // scale average alpha (around 0.5) when pulling to it
@@ -2007,6 +2009,7 @@ public class VegetationModel {
 					reg_weights,   // final double             reg_weights,        // fraction of the total weight used for regularization
 					alpha_loss,    // final double             alpha_loss,         // quadratic loss when alpha reaches -1.0 or 2.0 
 					alpha_offset,  // final double             alpha_offset,       // quadratic loss when alpha reaches -1.0 or 2.0
+					alpha_min_veg, //final double             alpha_min_veg, // 0.5; // if (alpha-alpha_offset)/(1-2*alpha_offset) < alpha_min_veg, pull down to lpha_offset
 					alpha_lpf,     // final double             alpha_lpf,          // pull to average of 4 neighbors
 					alpha_piece_linear, // final boolean       alpha_piece_linear, // true - piece-linear, false - half-cosine
 					alpha_scale_avg, // final double             alpha_scale_avg, //  = 1.2; // scale average alpha (around 0.5) when pulling to it