before per-pixel terrain elevation (fa638d8f) · Commits · Elphel / imagej-elphel

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

+49 −54

Original line number	Diff line number	Diff line
		@@ -107,10 +107,9 @@ public class VegetationLMA {
		private int [][][] elev_woi4; // [scene][woi_veg][~4] indices (in woi) of 4 (or more/less) neighbors of projections (negatives)
		private double [][][] elev_weights4; // [scene][woi_veg][~4] weights of influence for 4 neighbors
		private double [][] elev_sum_weights; // [scene][woi] sum weights from up to 4 elevation pixels (even more with overlap)
		// derivatives, only defined when elevations are adjusted NOT USED???
		private double [][][] elev_dweights4; // [scene][woi_veg][4] derivatives of weights of influence for 4 neighbors with respect to elevations
		// derivatives, only defined when elevations are adjusted NOT USED yes , not used (with derivatives verified by delta)
		/// private double [][][] elev_dweights4; // [scene][woi_veg][4] derivatives of weights of influence for 4 neighbors with respect to elevations
		private double [][][] elev_dsum_weights; // [scene][elev_par][woi] derivative of elev_sum_weights by each of the elevation parameter
		// private int [][][] elev_dindx4; // [scene][woi_veg][4] derivatives of weights of influence for 4 neighbors with respect to elevations


		// elevation-dependent parameters, calculated once if terrain elevation is not adjusted or each time if they are
		@@ -2262,7 +2261,7 @@ public class VegetationLMA {
		elev_weights4 = new double [num_scenes][][];
		elev_sum_weights = new double [num_scenes][];
		if (need_derivs) {
		elev_dweights4 = new double [num_scenes][][];
		// elev_dweights4 = new double [num_scenes][][];
		elev_dsum_weights = new double [num_scenes][][];
		}
		//elev_radius[num_pixels];
		@@ -2298,8 +2297,8 @@ public class VegetationLMA {
		elev_weights4[nScene] = new double [woi_veg_length][];
		elev_sum_weights[nScene] = new double [woi_length];
		if (need_derivs) {
		elev_dweights4[nScene] = new double [woi_veg_length][];
		elev_dsum_weights[nScene] = new double [num_pars[TVAO_ELEVATION]][woi_length];
		/// elev_dweights4[nScene] = new double [woi_veg_length][]; // not used ?
		elev_dsum_weights[nScene] = new double [num_pars[TVAO_ELEVATION]][woi_length]; // used
		}
		for (int wvindex = 0; wvindex < woi_veg_length; wvindex++) if (valid_vegetation[wvindex]) {
		int wvx = wvindex % woi_veg.width; // relative to woi_veg
		@@ -2358,9 +2357,9 @@ public class VegetationLMA {
		elev_woi4[nScene][wvindex] = new int[dia2]; // meter2];
		Arrays.fill(elev_woi4[nScene][wvindex], -1);
		elev_weights4[nScene][wvindex] = new double[dia2]; // meter2];meter2];
		if (need_derivs) {
		elev_dweights4[nScene][wvindex] = new double[dia2]; // meter2];meter2];
		}
		/// if (need_derivs) {
		/// elev_dweights4[nScene][wvindex] = new double[dia2]; // meter2];meter2];
		/// }
		}
		int dx = ipx - ipx0;
		int dindx = dx + dy*dia_x; // meter;
		@@ -2376,8 +2375,9 @@ public class VegetationLMA {
		elev_sum_weights[nScene][wpindex] += ww;
		if (need_derivs) {
		double dww = wnd_y * dwnd_x[dx] + dwnd_y*wnd_x[dx];
		elev_dweights4[nScene][wvindex][dindx] = dww;
		/// elev_dweights4[nScene][wvindex][dindx] = dww;
		int par_indx_elev = par_index[TVAO_ELEVATION][fnpix]; // ipx+ipy*full.width];
		// As par_indx_elev <-> fnpix <->wvindx, there is only one term in a sum below
		elev_dsum_weights[nScene][par_indx_elev-ind_pars[TVAO_ELEVATION]][wpindex] += dww;
		contrib_thread.get(wpindex).add(par_indx_elev); // full parameter index
		}
		@@ -3046,7 +3046,6 @@ public class VegetationLMA {
		final int dbg_cols=2,dbg_rows=2;
		final int dbg_ntpar = -44;
		final double alpha_threshold = alpha_0offset + alpha_min_veg *(1 - alpha_offset - alpha_0offset);
		// final double alpha_max_terrain = 0.75;
		final double alpha_threshold_pull = alpha_0offset + alpha_max_terrain *(1 - alpha_offset - alpha_0offset);

		for (int ithread = 0; ithread < threads.length; ithread++) {
		@@ -3055,9 +3054,7 @@ public class VegetationLMA {
		// vegetation projections to terrain/y_vector/fX pixels
		double [] alpha_woi = new double [woi_length];
		double [] veget_woi = new double [woi_length];
		// double [] terrain_woi = offset_terrain ? (new double [woi_length]) : null;
		double [] terrain_woi = new double [woi_length]; // needs to be divided by terrain_woi_weights
		// double [] terrain_weight_woi = offset_terrain ? (new double [woi_length]) : null;
		double [][] dalpha_woi = null;
		double [][] dveget_woi = null;
		double [][] dterrain_woi_dterrain = null; // derivatives of terrain in image [woi] per terrain [woi_terr]
		@@ -3147,15 +3144,15 @@ public class VegetationLMA {
		} else {
		alpha = tvao[TVAO_ALPHA][npix];
		}
		int [] terrain_woi_ind = elev_woi4[nScene][wvindex];
		if (terrain_woi_ind != null) {
		double [] terrain_woi_weights = elev_weights4[nScene][wvindex];// same vegetation pixel, per terrain neighbor
		for (int i = 0; i < terrain_woi_ind.length; i++) {
		int windx = terrain_woi_ind[i];
		int [] result_woi_ind = elev_woi4[nScene][wvindex]; // list of output pixels influenced by each vegetation pixel
		if (result_woi_ind != null) {
		double [] result_woi_weights = elev_weights4[nScene][wvindex];// same vegetation pixel, per result neighbor
		for (int i = 0; i < result_woi_ind.length; i++) {
		int windx = result_woi_ind[i];
		if (windx >= 0) {
		double w = terrain_woi_weights[i];
		veget_woi[windx] += w * vegetation; // accumulate vegetation for the same terrain output pixel
		alpha_woi[windx] += w * alpha; // accumulate alpha for the same terrain output pixel
		double w = result_woi_weights[i];
		veget_woi[windx] += w * vegetation; // accumulate vegetation for the same result pixel
		alpha_woi[windx] += w * alpha; // accumulate alpha for the same result pixel
		if (need_derivs) {
		if (nvpar >= 0) {
		contrib_thread.get(windx).add(nvpar); // full parameter index - this result pixel depends on this vegetation parameter
		@@ -3302,12 +3299,6 @@ public class VegetationLMA {
		contrib_thread.get(windx).add(ind_pars[TVAO_TERR_ELEV]); // full parameter index - this result pixel depends on this terrain parameter
		jt[ind_pars[TVAO_TERR_ELEV]][y_indx] += dterrain_woi_dtelev[windx] * (1 - k);// sum_a; // d/dterrain
		}
		/*
		if (ntpar >= 0) {
		jt[ntpar][y_indx] = 1 - k;// sum_a; // d/dterrain
		}
		*/


		if (vegetation_parameters \|\| alpha_parameters) {
		HashSet<Integer> pars_set = contrib_thread.get(windx); // may have unrelated (from other scenes) but still a small number
		@@ -3371,26 +3362,41 @@ public class VegetationLMA {
		/*
		* avg_veget, avg_alpha - weighted-averaged (divided by sw) from elevated contributors
		* result pixel value d combining vegetation and terrain:
		* d = terrain * (1.0 - k) + avg_veget * k;
		* ( 1) d = terrain * (1.0 - k) + avg_veget * k;
		* derivative of the result pixel d with respect to iterated elevation parameter epar:
		* dd_depar = k * davg_veget_depar +dk_depar * (avg_veget-terrain)
		* ( 2) dd_depar = k * davg_veget_depar +dk_depar * (avg_veget-terrain)
		* derivative of k (vegetation contribution derived from alpha_avg) with respect
		* to the iterated elevation parameter epar:
		* dk_depar = dk_dalpha_avg * dalpha_avg_depar
		* dk_dalpha_avg = 1.0 \| 0.5 * PI * sin(avg_alpha * PI) (depending on mode)
		* dalpha_avg_depar =dsw_dnepar/sw * (alpha_src - avg_alpha) ??
		* dalpha_avg_depar= (d_swa_d_elev * sw - swadsw_delev)/ (sw sw) =
		* (d_swa_d_elev_i - (swa/sw) * dsw_delev) / sw =
		* (d_swa_d_elev_i - avg_alpha * dsw_delev) / sw =
		* (elev_dweights4[nScene][wvindex][dindx] - avg_alpha * elev_dweights4[nScene][windx][dindx])/ sw
		* ( 3) dk_depar = dk_dalpha_avg * dalpha_avg_depar
		* ( 4) dk_dalpha_avg = 1.0 \| 0.5 * PI * sin(avg_alpha * PI) (depending on mode)
		* ( 5) avg_veget = sum(veget[epar]*w[epar])/sum(w[epar]), were [epar] means "corresponding to"
		* ( 6) davg_veget_depar = d_(sum(veget[epar]*w[epar])/sum(w(epar)))_depar
		* ( 7) davg_veget_depar = (d_(sum(veget[epar]w[epar])_depar sum(w[epar]) - sum(veget[epar]w[epar]) d_(sum(w[epar]))_depar)/
		* (sum(w[epar]) * sum(w[epar]))
		* ( 8) sw = sum(w(epar)
		* ( 9) dsw_dnepar = d_(sum(w[epar]))_depar has only one term, so it is
		* (10) dsw_dnepar = d_w[epar]_depar
		* In d_(sum(veget[epar]*w[epar])_depar there is only one term (others do not depend on epar)
		* so
		* (11) d_(sum(veget[epar]w[epar])_depar = d_(veget[epar]w[epar])_depar
		* veget[epar] is just veget_src and does not depend on epar:
		* (12) d_(sum(veget[epar]w[epar])_depar = veget_src d_w[epar]_depar
		* using (11) that dsw_dnepar has just one term:
		* (13) d_(sum(veget[epar]w[epar])_depar = veget_src dsw_dnepar
		* and (7) becomes
		* (14) davg_veget_depar = (veget_srcdsw_dneparsw -sum(veget[epar]w[epar])dsw_dnepar)/(sw * sw)
		* using (5)
		* (15) davg_veget_depar = (veget_srcdsw_dneparsw - (avg_alphasw)dsw_dnepar)/(sw * sw)
		* (16) davg_veget_depar = (veget_srcdsw_dnepar - avg_alphadsw_dnepar)/sw
		* (17) davg_veget_depar = dsw_dnepar(veget_src - avg_alpha)/sw
		* similarly if alpha >=0 and alpha <= 1.0:
		* (17) davg_alpha_depar = dsw_dnepar(alpha_src - avg_alpha)/sw
		*/
		double davg_veget_depar= dsw_dnepar/sw * (veget_src - avg_veget);
		double dd_depar = k * davg_veget_depar;
		if ((avg_alpha >= 0) && (avg_alpha <= 1.0)) {
		double dk_dalpha_avg = alpha_piece_linear ? 1.0 : (0.5 * Math.PI * Math.sin(avg_alpha*Math.PI));
		double dalpha_avg_depar = dsw_dnepar/sw * (alpha_src - avg_alpha);
		// double dalpha_avg_depar1 = (elev_dweights4[nScene][windx][dindx] - avg_alpha * elev_dsum_weights[nScene][windx][dindx])/ sw;

		double dk_depar = dk_dalpha_avg * dalpha_avg_depar;
		dd_depar += dk_depar * (avg_veget-terrain);
		}
		@@ -3423,16 +3429,6 @@ public class VegetationLMA {
		jt[ind_pars[TVAO_TERR_ELEV]][y_indx] += dterrain_woi_dtelev[windx];// sum_a; // d/dterrain
		}
		}
		/*
		if (need_derivs && (ntpar >= 0)) {
		jt[ntpar][y_indx] = 1.0;
		}
		*/
		/*
		if (extra_data != null) {
		extra_data[0][y_indx] = terrain + scene_offset;
		}
		*/
		if (extra_data != null) {
		double k = 0.0;
		if (extra_data.length > 0) {
		@@ -3451,12 +3447,6 @@ public class VegetationLMA {

		}
		fX[y_indx] = d + scene_offset;
		/*
		// already added
		if (need_derivs && (ntpar >= 0)) {
		contrib_thread.get(windx).add(ntpar); // full parameter index
		}
		*/
		if (need_derivs && (nspar >= 0)) {
		jt[nspar][y_indx] = 1.0;
		contrib_thread.get(windx).add(nspar); // full parameter index
		@@ -5867,6 +5857,7 @@ public class VegetationLMA {
		int debugLevel,
		boolean show_img) {
		double [][] jt = new double [vector.length][];
		boolean dbg_reset_before_elevation = false;
		getFxDerivs(
		vector, // final double [] vector,
		jt, // final double [][] jt, // should be null or initialized with [vector.length][]
		@@ -5878,6 +5869,9 @@ public class VegetationLMA {
		double [][] jt_diff = new double [jt_delta.length][jt_delta[0].length];
		double max_err=0;
		for (int n = 0; n < jt_diff.length; n++) {
		if (dbg_reset_before_elevation && (n == ind_pars[TVAO_ELEVATION])) {
		max_err = 0;
		}
		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) \|\| (Math.abs(jt_diff[n][w]) > .001)) {
		@@ -5913,6 +5907,7 @@ public class VegetationLMA {
		System.out.println("debugDerivs(): "+s1+s3+s2);
		}
		max_err = Math.max(max_err, Math.abs(jt_diff[n][w]));

		}
		// max_err = Math.max(max_err, Math.abs(jt_diff[n][w]));
		}

Admin message