minor bug fixes, verified jacobian

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

@@ -52,11 +52,17 @@ public class VegetationLMA {
 	private double []         y_vector = null;
 	private double            weight_pure =     0;
 	private double []         weights; // normalized so sum is 1.0 for all - samples and extra regularization
-	public double             alpha_loss;
+	public double             alpha_loss =      0; // not used with cosine alpha
 	public double             alpha_offset =    0; // if >0,  start losses above 0.0 and below 1.0;
 	public double             alpha_lpf =       0;
 	public double             terr_lpf =        0;
 	public double             veget_lpf =       0;
+	
+	// when unsharp mask is applied , pulling to 0 (when alpha is 0 (for vegetation) or 1.0 (for terrain) makes sense
+	public double             terr_pull0 =      0;
+	public double             veget_pull0 =     0;
+	
+	
 	public double             boost_parallax =  1;
 	public double             um_sigma =        0;       // just use in debug image names
 	public double             um_weight =       0;
@@ -130,6 +136,8 @@ public class VegetationLMA {
 			final double             alpha_lpf,      // pull vegetation alpha to average of 4 neighbors
 			final double             terr_lpf,       // pull terrain to average of 4 neighbors (very small)
 			final double             veget_lpf,      // pull vegetation to average of 4 neighbors (very small - maybe not needed)
+			final double             terr_pull0,     // pull terrain to zero (makes sense with UM
+			final double             veget_pull0,    // pull vegetation to zero (makes sense with UM
 			final double             boost_parallax, // increase weight of scene with maximal parallax relative to the reference scene
 			final double             um_sigma,       // just use in debug image names
 			final double             um_weight,
@@ -831,11 +839,240 @@ public class VegetationLMA {
 		ImageDtt.startAndJoin(threads);
 		
 		// regularization weights and derivatives
-    return terrain_weights;
+		return terrain_weights;
     }
 	
 	
 	private double [] getFxDerivs(
+			final double []   vector,
+			final double [][] jt, // should be null or initialized with [vector.length][]
+			final int         debug_level)	{
+		// using 0.5*(1-cos(alpha/2pi) instead of alpha. alpha < 0 -> 0, alpha > 1 -> 1. Depends on other terms for stability 
+		double [] fX = new double [weights.length]; // num_pairs + vector.length];
+		if (jt != null) {
+			for (int i = 0; i < jt.length; i++) {
+				jt[i] = new double [weights.length]; // weights.length];
+			}
+		}
+		final Thread[] threads = ImageDtt.newThreadArray();
+		final AtomicInteger ai = new AtomicInteger(0);
+		for (int ithread = 0; ithread < threads.length; ithread++) {
+			threads[ithread] = new Thread() {
+				public void run() {
+					double [] vegetation = new double [4]; 
+					double [] alpha =      new double [4];
+					for (int n = ai.getAndIncrement(); n < y_vector.length; n = ai.getAndIncrement()) {
+//						int nscene =   data_source[n][0][0];
+//						int indx =     data_source[n][0][1];
+						double terrain = vector[data_source[n][0][2]];
+						double [] cw = corners_weights[n]; 
+						double d; 
+						int [] indx_vegetation = data_source[n][1]; 
+						int [] indx_alpha =      data_source[n][2];
+						double sum_v =0, sum_a =0;
+						if (cw != null) {
+							for (int i = 0; i < 4; i++ ) {
+								int iv = indx_vegetation[i], ia=indx_alpha[i];
+								vegetation[i] = cw[i] * ((iv >= 0) ? vector[iv]: tvao[TVAO_VEGETATION][-1-iv]);
+								alpha[i] =      cw[i] * ((ia >= 0) ? vector[ia]: tvao[TVAO_VEGETATION_ALPHA][-1-ia]);
+								sum_v += vegetation[i];
+								sum_a += alpha[i];
+							}
+							double k = (sum_a < 0)? 0: ( (sum_a > 1) ? 1.0: 0.5 * (1.0 - Math.cos(sum_a*Math.PI)));
+//							d = terrain * (1.0 - sum_a) + sum_v * sum_a;
+							d = terrain * (1.0 - k) + sum_v * k;
+							
+							if (jt != null) {
+								jt[data_source[n][0][2]][n] = 1 - k;// sum_a; // d/dterrain
+								for (int i = 0; i < 4; i++ ) {
+									if (indx_vegetation[i] >= 0) {
+										jt[data_source[n][1][i]][n] = cw[i] * k; // sum_a; // d/dvegetation[i]
+									}
+									if ((indx_alpha[i] >= 0) && (sum_a > 0) && (sum_a < 1.0)) {
+//										jt[data_source[n][2][i]][n] = cw[i] * (sum_v - terrain); // d/dalpha[i]
+										
+										jt[data_source[n][2][i]][n] = cw[i] * (sum_v - terrain) *0.5*Math.PI *Math.sin(sum_a*Math.PI); // d/dalpha[i]
+									}
+								}								
+							}
+
+						} else {
+							d = terrain; 
+							if (jt != null) {
+								jt[data_source[n][0][2]][n] = 1;  // d/dterrain
+							}
+						}
+						double scene_offs = vector[data_source[n][0][3]];
+						fX[n] = d + scene_offs;
+						if (jt != null) {
+							jt[data_source[n][0][3]][n] = 1;
+						}
+					}
+				}
+			};
+		}
+		ImageDtt.startAndJoin(threads);
+		// regularization weights and derivatives
+		
+		int ind_next = y_vector.length;
+		if ((alpha_lpf >= 0) || (alpha_loss > 0)) {
+			final int ind_y_alpha = ind_next;
+			ind_next += num_pars_vegetation_alpha;
+			ai.set(0);
+			for (int ithread = 0; ithread < threads.length; ithread++) {
+				threads[ithread] = new Thread() {
+					public void run() {
+						for (int n = ai.getAndIncrement(); n < num_pars_vegetation_alpha; n = ai.getAndIncrement()) {
+							int np = ind_pars_vegetation_alpha + n; // index of the alpha parameter
+							int nx = n + ind_y_alpha; //  y_vector.length; // x - index
+							double d = 0;
+							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]		
+									}							
+								}
+							}
+							// add cost for difference between this alpha and average of 4 neighbors (when they exist
+							// applies to alpha before cosine, so it will pull borders even when alpha<0 or alpha > 1 (zero derivatives)
+							if (alpha_lpf > 0) { // should always be > 0 to provide stability for out-of-range alpha 
+								double avg = 0;
+								int nn = 0;
+								for (int i = 0; i < alpha_neibs[n].length; i++) { // now 4, may be increased
+									int di =  alpha_neibs[n][i];
+									d=0;
+									if (di >= 0) {
+										d = vector[di];  // d - full parameter index
+										avg+=d;
+										nn++;
+									} else if (di < -1) {
+										d = tvao[TVAO_VEGETATION_ALPHA][-di - 2]; 
+										avg+=d;
+										nn++;
+									}
+								}
+								avg /= nn; // average
+								fX[nx] += alpha_lpf * (vector[np] - avg);
+								if (jt != null) {
+									jt[np][nx] += alpha_lpf;
+									for (int i = 0; i < alpha_neibs[n].length; i++) { // now 4, may be increased
+										int di =  alpha_neibs[n][i];
+										if (di >= 0) {
+											jt[di][nx] -= alpha_lpf/nn;
+										}
+									}								
+								}							
+							}
+						}
+					}
+				};
+			}
+			ImageDtt.startAndJoin(threads);
+		} // if (alpha_lpf >= 0) {
+		
+		if (terr_lpf >= 0) {
+			final int ind_y_terr = ind_next;
+			ind_next += num_pars_terrain;
+			ai.set(0);
+			for (int ithread = 0; ithread < threads.length; ithread++) {
+				threads[ithread] = new Thread() {
+					public void run() {
+						for (int n = ai.getAndIncrement(); n < num_pars_terrain; n = ai.getAndIncrement()) {
+							int np = ind_pars_terrain + n; // index of the alpha parameter
+							int nx = n + ind_y_terr; //  y_vector.length; // x - index
+							double d = 0;
+							if (terr_lpf > 0) {
+								double avg = 0;
+								int nn = 0;
+								for (int i = 0; i < terr_neibs[n].length; i++) { // now 4, may be increased
+									int di =  terr_neibs[n][i];
+									d=0;
+									if (di >= 0) {
+										d = vector[di];  // d - full parameter index
+										avg+=d;
+										nn++;
+									} else if (di < -1) {
+										d = tvao[TVAO_TERRAIN][-di - 2]; 
+										avg+=d;
+										nn++;
+									}
+								}
+								avg /= nn; // average
+								fX[nx] += terr_lpf * (vector[np] - avg) + terr_pull0 * vector[np];
+								if (jt != null) {
+									jt[np][nx] += terr_lpf + terr_pull0;
+									for (int i = 0; i < terr_neibs[n].length; i++) { // now 4, may be increased
+										int di =  terr_neibs[n][i];
+										if (di >= 0) {
+											jt[di][nx] -= terr_lpf/nn;
+										}
+									}								
+								}							
+							}
+						}
+					}
+				};
+			}
+			ImageDtt.startAndJoin(threads);
+		}
+		if (veget_lpf >= 0) {
+			final int ind_y_veget = ind_next;
+			ind_next += num_pars_vegetation;
+			ai.set(0);
+			
+			for (int ithread = 0; ithread < threads.length; ithread++) {
+				threads[ithread] = new Thread() {
+					public void run() {
+						for (int n = ai.getAndIncrement(); n < num_pars_vegetation; n = ai.getAndIncrement()) {
+							int np = ind_pars_vegetation + n; // index of the alpha parameter
+							int nx = n + ind_y_veget; //  y_vector.length; // x - index
+							double d = 0;
+							if (veget_lpf > 0) {
+								double avg = 0;
+								int nn = 0;
+								for (int i = 0; i < veget_neibs[n].length; i++) { // now 4, may be increased
+									int di =  veget_neibs[n][i];
+									d=0;
+									if (di >= 0) {
+										d = vector[di];  // d - full parameter index
+										avg+=d;
+										nn++;
+									} else if (di < -1) {
+										d = tvao[TVAO_VEGETATION][-di - 2]; 
+										avg+=d;
+										nn++;
+									}
+								}
+								avg /= nn; // average
+								fX[nx] += veget_lpf * (vector[np] - avg)  + veget_pull0 * vector[np];
+								if (jt != null) {
+									jt[np][nx] += veget_lpf + veget_pull0;
+									for (int i = 0; i < veget_neibs[n].length; i++) { // now 4, may be increased
+										int di =  veget_neibs[n][i];
+										if (di >= 0) {
+											jt[di][nx] -= veget_lpf/nn;
+										}
+									}								
+								}							
+							}
+						}
+					}
+				};
+			}
+			ImageDtt.startAndJoin(threads);
+		}
+		return fX;		
+	}
+	
+	private double [] getFxDerivs_precos(
 			final double []   vector,
 			final double [][] jt, // should be null or initialized with [vector.length][]
 			final int         debug_level)
@@ -1058,6 +1295,8 @@ public class VegetationLMA {
 		}
 		return fX;		
 	}
+
+	
 	
 	private double [][] getFxDerivsDelta(
 			double []         vector,
@@ -1255,17 +1494,20 @@ public class VegetationLMA {
 		for (int n = 0; n < jt_diff.length; n++) {
 			for (int w = 0; w < jt_diff[0].length; w++) {
 				jt_diff[n][w] = jt[n][w] - jt_delta[n][w];
+				if (Math.abs(jt_diff[n][w]) > max_err) {
+					System.out.println("debugDerivs(): n="+n+", w = "+w +", diff="+jt_diff[n][w]);
+				}
 				max_err = Math.max(max_err, Math.abs(jt_diff[n][w]));
 			}
 		}
 		System.out.println("delta = "+delta+", max_err = "+max_err);
 		// think of visualization
-		/*
+		/**/
 		if (show_img) {
 			String [] frame_titles = {"jt","jt_delta", "jt_diff"}; 
 			double [][][] debug_img = new double [frame_titles.length][][];
 		}
-		*/
+		/**/
 		return max_err;
 	}
 	

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

@@ -681,28 +681,31 @@ public class VegetationModel {
 		int       min_scenes =   10;
 		double    default_alpha = 0.8;
 		double    reg_weights =   0.25;        // fraction of the total weight used for regularization
-		double    alpha_loss =    10000.0; // 1000.0; // 100.; // 10.0;         // quadratic loss when alpha reaches -1.0 or 2.0 
-		double    alpha_offset =  0.02; // 0.03;    // if >0,  start losses above 0.0 and below 1.0;
-		double    alpha_lpf =     10; // 20; // 6.0; // 3.0; // 2.0;    // 1.5; // 5.0; // 0.5;         // pull to average of 4 neighbors
+		double    alpha_loss =   10.0;  // 10000.0; // 1000.0; // 100.; // 10.0;         // quadratic loss when alpha reaches -1.0 or 2.0 
+		double    alpha_offset =  0.0; // 0.02; // 0.03;    // if >0,  start losses above 0.0 and below 1.0;
+		double    alpha_lpf =     10;   // 20; // 6.0; // 3.0; // 2.0;    // 1.5; // 5.0; // 0.5;         // pull to average of 4 neighbors
 		double    terr_lpf =      0.1;    // pull terrain to average of 4 neighbors (very small)
 		double    veget_lpf =     0.1;    // pull vegetation to average of 4 neighbors (very small - maybe not needed)
-
+		double    terr_pull0 =    0.1;    //pull terrain to zero (makes sense with UM
+		double    veget_pull0 =   0.1;    // pull vegetation to zero (makes sense with UM
 		double    boost_parallax = 1.0; //  5;
-		boolean exit_loop = debugLevel < 1000;
-		boolean next_run = false;
-		boolean read_pars = false; // true; //  false; // true;
-		double      threshold_terrain = 0.05;
-		double      min_max_terrain= 0.1;
-		double      min_terrain =    0.001;
-		double      min_vegetation = 0.5;				
+		boolean   next_run =      false;
+		boolean   read_pars =     true; //   false; // true;
+		double    threshold_terrain = 0.05;
+		double    min_max_terrain= 0.1;
+		double    min_terrain =    0.001;
+		double    min_vegetation = 0.5;				
+		boolean   um_en =     true;
+		double    um_sigma =  1.0; 
+		double    um_weight = 0.8;
+
+		boolean   exit_loop = debugLevel < 1000;
 
 		//		String  parameters_path = "/media/elphel/SSD3-4GB/lwir16-proc/berdich3/debug/vegetation/lma/parameters_vector_data_1000-0.03.tiff";
 		//		String  parameters_path = "/media/elphel/SSD3-4GB/lwir16-proc/berdich3/debug/vegetation/essential/parameters_vector_data_1000-0.03.tiff";
-		String  parameters_path = "/media/elphel/SSD3-4GB/lwir16-proc/berdich3/debug/vegetation/essential/parameters_vector_data_x143-y317-w35-h35-live_10000_0.02_3.0.tiff";
-
-		final boolean     um_en =     true;
-		final double      um_sigma =  1.0; 
-		final double      um_weight = 0.8;
+//		String  parameters_path = "/media/elphel/SSD3-4GB/lwir16-proc/berdich3/debug/vegetation/essential/parameters_vector_data_x143-y317-w35-h35-live_10000_0.02_3.0.tiff";
+//		String  parameters_path = "/media/elphel/SSD3-4GB/lwir16-proc/berdich3/debug/vegetation/essential/parameters_vector_data_x143-y317-w35-h35-al0.0-alo0.02-alp10.0-tl0.1-vl0.1-bp1.0-um1.0_0.8.tiff";
+		String  parameters_path = "/media/elphel/SSD3-4GB/lwir16-proc/berdich3/debug/vegetation/essential/parameters_vector_data_x143-y317-w35-h35-al10.0-alo0.0-alp10.0-tl0.1-vl0.1-bp1.0-um1.0_0.8-live.tiff";
 		if (um_en) {
 			double [][] um_data = new double [terrain_rendered.length+2][];
 			System.arraycopy(terrain_rendered, 0, um_data, 0, terrain_rendered.length);
@@ -786,6 +789,8 @@ public class VegetationModel {
 					alpha_lpf,     // final double             alpha_lpf,           // pull to average of 4 neighbors
 					terr_lpf,      // final double             terr_lpf,       // pull terrain to average of 4 neighbors (very small)
 					veget_lpf,     // final double             veget_lpf,      // pull vegetation to average of 4 neighbors (very small - maybe not needed)
+					terr_pull0,     // final double             terr_pull0,     // pull terrain to zero (makes sense with UM
+					veget_pull0,    // final double             veget_pull0,    // pull vegetation to zero (makes sense with UM
 					boost_parallax,// final double             boost_parallax, // increase weight of scene with maximal parallax relative to the reference scene
 					um_sigma,      // final double             um_sigma,       // just use in debug image names
 					(um_en? um_weight: 0.0),     // final double             um_weight,