next working snapshot

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

@@ -100,7 +100,7 @@ public class VegetationLMA {
 	
 //[2][image_length] - alpha, [3][num_scenes] - per-scene offset (mostly for vignetting. MAYBE add per-scene pair of tilts
 	private int [][]          par_index;       // indices -  [0..4][full_pixel] - same as for tvao, value - parameter index
-	private int [][]          par_rindex;      // parameter -> pair of type (0..3) and full window pixel index
+	private int [][]          par_rindex;      // parameter -> pair of type (0..4) and full window pixel index
 	private int               num_pars_terrain; 
 	private int               num_pars_vegetation; 
 	private int               num_pars_alpha;            //== num_pars_vegetation
@@ -186,7 +186,9 @@ public class VegetationLMA {
 	private boolean            alpha_en_holes   =   true; // Search for small semi-transparent holes, disable diffusion of local alpha minimums
 	private double             alpha_diff_hole =    0.01; // minimal alpha difference between min and max neighbor to be considered a hole
 
-	
+	private double             w_alpha_neib =       0.7; // weight of each of the neighbors relative to the center when calculating average alpha
+	private double             holes_pwr =          2.0; // Raise alpha (0..1) to this power when replacing vegetation cold spots with holes
+
 	private boolean            from_file =       false;
 
 	
@@ -2535,8 +2537,6 @@ 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 int ind_y_alpha_loss = samples_pointers[SAMPLES_ALPHA_PULL][0]; // ind_next;
-//			ind_next += num_pars_alpha;
 			ai.set(0);
 			for (int ithread = 0; ithread < threads.length; ithread++) {
 				threads[ithread] = new Thread() {
@@ -2563,41 +2563,7 @@ public class VegetationLMA {
 									}							
 								}
 							}
-							double avg = 0;
-							int nn = 0;
-							// removing all dependence on neighbors!
-							/*
-							double neib_min = Double.POSITIVE_INFINITY, neib_max = Double.NEGATIVE_INFINITY;
-							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;
-									if (d < neib_min) neib_min = d; 
-									if (d > neib_max) neib_max = d; 
-									nn++;
-								} else if (di < -1) {
-									d = tvao[TVAO_ALPHA][-di - 2]; 
-									avg+=d;
-									if (d < neib_min) neib_min = d; 
-									if (d > neib_max) neib_max = d; 
-									nn++;
-								}
-							}
-							avg /= nn; // average
-							*/
 							if (alpha_push > 0) {
-								/*
-								 * p = alpha_push
-								 * n = alpha_push_neutral
-								 * a0 =        p/(2 * n)
-								 * a1 =        p/(2 * (1 - n))
-								 * f1(x) =     p/(2 * n)* x * (1 - x/(2*n))
-								 * f2(x) =     p * (x - 2*n + 1) * (1 - x) / (4 * (1-n)*(1-n) 
-								 * df1(x)/dx = p/(2 * n^2 )*( n - x)
-								 * df2(x)/dx = p* (n-x)/(2*(1-n)^2)
-								 */
 								// average including center:
 ///								double sum_w = (nn + alpha_push_center);
 ///								double avg_c = (avg * nn + vector[np] * alpha_push_center) / sum_w;
@@ -2622,14 +2588,6 @@ public class VegetationLMA {
 											dd = alpha_push / (sum_w * 2 * one_minus_n * one_minus_n)  * (alpha_push_neutral  - avg_c); 
 										}
 										jt[np][nx] += dd * alpha_push_center; // 2 * alpha_loss * d; // d/dalpha[i]	
-										/*
-										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] += dd; // 2 * alpha_loss * d; // d/dalpha[i]
-											}
-										}
-										*/
 									}								
 								}
 							}
@@ -2834,10 +2792,8 @@ public class VegetationLMA {
 			ImageDtt.startAndJoin(threads);
 		}
 
-		
+/*
 		if (samples_pointers[SAMPLES_VEGETATION_LPF][1] > 0) { // fits[TVAO_VEGETATION] && (veget_lpf >= 0)) { // should be positive for pull0 and terr_pull_cold (difference between vegetation and terrain)
-//			final int ind_y_veget = samples_pointers[SAMPLES_VEGETATION_LPF][0]; // ind_next;
-			//			ind_next += num_pars_vegetation;
 			ai.set(0);
 			for (int ithread = 0; ithread < threads.length; ithread++) {
 				threads[ithread] = new Thread() {
@@ -2880,6 +2836,135 @@ public class VegetationLMA {
 			}
 			ImageDtt.startAndJoin(threads);
 		}
+	
+ */
+
+
+		if (samples_pointers[SAMPLES_VEGETATION_LPF][1] > 0) { // should be positive for pull0 and terr_pull_cold (difference between vegetation and terrain)
+			if  (num_pars_alpha <= 0) {
+				System.out.println("getFxDerivs(): BUG: Asjusting SAMPLES_VEGETATION_LPF requires alpha parameters aslo to be fitted");
+			} else {
+				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 = n + samples_pointers[SAMPLES_VEGETATION_LPF][2];
+								int nx = n + samples_pointers[SAMPLES_VEGETATION_LPF][0];
+								int findx = par_rindex[np][1];//
+								int npa = par_index[TVAO_ALPHA][findx]; // parameter for same-pixel alpha
+								int na = npa - ind_pars_alpha;
+								if (npa < 0) {
+									System.out.println("getFxDerivs(): BUG: no alpha for vegetation, n="+n+", np="+np+", nx ="+nx+
+											" findx="+findx+", npa="+npa);
+									continue;
+								}
+								double d = 0;
+								if (veget_lpf > 0) {
+									double veget_avg4 = 0;
+									int veget_nn4 = 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
+											veget_avg4+=d;
+											veget_nn4++;
+										} else if ((di < -1) && lpf_fixed[TVAO_VEGETATION]) {
+											d = tvao[TVAO_VEGETATION][-di - 2]; 
+											veget_avg4+=d;
+											veget_nn4++;
+										}
+									}
+									if (veget_nn4 <= 0) {
+										System.out.println("getFxDerivs(): BUG: veget_nn4="+veget_nn4+", n="+n+", np="+np+", nx ="+nx+
+												" findx="+findx+", npa="+npa);
+										continue;
+									}
+									veget_avg4 /= veget_nn4; // average
+									double veget_lap = vector[np] - veget_avg4;
+									if (veget_lap >=0) { // only for cold spots
+										continue;
+									}
+									//								fX[nx] += veget_lpf * (vector[np] - veget_avg4); //   + veget_pull0 * vector[np];
+									//								double w_alpha_neib = 0.7;
+									//								double holes_pwr =    2.0;
+									d = vector[npa];
+									double alpha_avg5; 
+									if (alpha_piece_linear) {
+										alpha_avg5 = (d < 0)? 0: ((d > 1) ? 1.0: d);	
+									} else {
+										alpha_avg5 = (d < 0)? 0: ((d > 1) ? 1.0: 0.5 * (1.0 - Math.cos(d*Math.PI)));	
+									}
+//									double alpha_avg5 = vector[npa]; // center alpha, weight = 1.0;
+									double alpha_w5 = 1.0;
+//									double [] derivs_neibs = new double [ alpha_neibs[n].length]; 
+									for (int i = 0; i < alpha_neibs[n].length; i++) { // now 4, may be increased
+										int di =  alpha_neibs[na][i];
+										d=0;
+										if (di >= 0) {
+											d = vector[di];  // d - full parameter index
+										} else if ((di < -1) && lpf_fixed[TVAO_ALPHA]) {
+											d = tvao[TVAO_ALPHA][-di - 2]; 
+										}
+										double alpha;
+
+										if (alpha_piece_linear) {
+											alpha = (d < 0)? 0: ((d > 1) ? 1.0: d);	
+										} else {
+											alpha = (d < 0)? 0: ((d > 1) ? 1.0: 0.5 * (1.0 - Math.cos(d*Math.PI)));	
+										}
+										alpha_avg5 += alpha * w_alpha_neib;
+										alpha_w5 += w_alpha_neib;
+										// TODO: calculate derivatives
+									}
+									// limit to 0..1 and apply cos
+									alpha_avg5 /= alpha_w5; // average
+									double alpha_pwr =  Math.pow(alpha_avg5, holes_pwr);
+									double val = veget_lap * alpha_pwr; 
+									fX[nx] += veget_lpf * val;
+
+									if (jt != null) {
+										jt[np][nx] += veget_lpf * alpha_pwr;
+										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* alpha_pwr/veget_nn4;
+											}
+										}
+										/* dval_dai = dval_dalpha_pwr * dalpha_pwr_dalpha * dalpha_dfai * dfai_dai; */
+										double dval_dalpha_pwr = veget_lpf * veget_lap;
+										double dalpha_pwr_dalpha = holes_pwr* Math.pow(alpha_avg5, holes_pwr - 1);
+										double dval_dalpha = dval_dalpha_pwr * dalpha_pwr_dalpha;
+										double dalpha_dai = 1.0 / alpha_w5; // dval_dalpha / alpha_w5;
+										if ((vector[npa] > 0) && (vector[npa] < 1.0)) {
+											if (!alpha_piece_linear) {
+												dalpha_dai *= 0.5 * Math.PI * Math.sin(vector[npa]*Math.PI);
+											}
+											jt[npa][nx] += dval_dalpha * dalpha_dai;
+										}
+										for (int i = 0; i < alpha_neibs[na].length; i++) { // now 4, may be increased
+											int di =  alpha_neibs[na][i];
+											if (di >= 0) {
+												if ((vector[di] > 0) && (vector[di] < 1.0)) {
+//													dalpha_dai = dval_dalpha * w_alpha_neib / alpha_w5;
+													dalpha_dai = w_alpha_neib / alpha_w5;
+													if (!alpha_piece_linear) {
+														dalpha_dai *= 0.5 * Math.PI * Math.sin(vector[di]*Math.PI);
+													}
+													jt[di][nx] += dval_dalpha * dalpha_dai;
+												}
+											}
+										}
+									}
+								}
+							}
+						}
+					};
+				}
+				ImageDtt.startAndJoin(threads);
+			}
+		}
 		
 		if (samples_pointers[SAMPLES_ELEVATION_PULL][1]>0) { // fits[TVAO_ELEVATION] && (elevation_lpf >= 0)) { // should be positive for pull0 and terr_pull_cold (difference between vegetation and terrain)
 			ai.set(0);
@@ -4280,7 +4365,7 @@ public class VegetationLMA {
 				int [] ti = par_rindex[n];
 				if (Math.abs(jt_diff[n][w]) > max_err) {
 					System.out.println("debugDerivs(): n="+n+" ("+ti[0]+":"+ti[1]+"), w = "+w +", diff="+jt_diff[n][w]+" jt="+jt[n][w]+" jt_delta="+jt_delta[n][w]);
-				} else if (Math.abs(jt_diff[n][w]) > .1) {
+				} else if (Math.abs(jt_diff[n][w]) > .001) {
 					System.out.println("debugDerivs(): n="+n+" ("+ti[0]+":"+ti[1]+"), w = "+w +", diff="+jt_diff[n][w]+" jt="+jt[n][w]+" jt_delta="+jt_delta[n][w]);
 				}
 				max_err = Math.max(max_err, Math.abs(jt_diff[n][w]));