Next working snapshot (3f011b38) · Commits · Elphel / imagej-elphel

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

+392 −62

Original line number	Diff line number	Diff line
		@@ -21,8 +21,6 @@ public class VegetationLMA {
		public static final int TVAO_SCENE_OFFSET = 3;

		public final Rectangle full;
		// public final int image_width;
		// public final int image_height;
		public final int image_length;
		public final int num_scenes;
		public Rectangle woi = null;
		@@ -32,14 +30,10 @@ public class VegetationLMA {
		public double [][][] vegetation_offsets; // [num_scenes][pixle]{dx,dy} differential offsets of vegetation to terrain grid
		public boolean diff_offsets;
		// next 3 arrays are full [image_width*image_height], but will be modified only inside woi
		// public double [] terrain; //
		// public double [] vegetation; //
		// public double [] alpha; //
		public double [][] tvao; // [0][image_length] - terrain, [1][image_length] - vegetation,
		//[2][image_length] - alpha, [3][num_scenes] - per-scene offset (mostly for vignetting. MAYBE add per-scene pair of tilts
		public int [][] par_index; // indices - [0..3][full_pixel] - same as for tvao, value - parameter index
		public int [][] par_rindex; // parameter -> pair of type (0..3) and full window pixel index
		// private double default_alpha;
		private int num_pars_terrain;
		private int num_pars_vegetation;
		private int num_pars_vegetation_alpha; //== num_pars_vegetation
		@@ -61,9 +55,15 @@ public class VegetationLMA {
		public double alpha_loss;
		public double alpha_offset = 0; // if >0, start losses above 0.0 and below 1.0;
		public double alpha_lpf = 0;
		// data used to calculate lpf pull of the alpha pixel to average of four neighbors
		public double terr_lpf = 0;
		public double veget_lpf = 0;
		// data used to calculate lpf pull of the alpha pixel to average of four neighbors. Below similar (weaker pull) for terrain and vegetation
		// to smooth areas where there is no data from available images.
		public int [][] alpha_neibs; // corresponds to parameters for alpha (num_pars_vegetation_alpha), each has 4 ortho neibs, -1 - border, >= 0
		// - parameter index, <=-2 -2-full_pixel_index
		public int [][] terr_neibs; // corresponds to parameters for terrain, for smoothing undefined terrain, similar to alpha_neibs
		public int [][] veget_neibs; // corresponds to parameters for vegetation, for smoothing undefined vegetation, similar to alpha_neibs


		public double delta= 1e-5; // 7;

		@@ -123,7 +123,9 @@ 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_lpf, // pull to average of 4 neighbors
		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 boost_parallax, // increase weight of scene with maximal parallax relative to the reference scene
		final String parameters_read_path,
		final int debugLevel) {
		@@ -131,6 +133,8 @@ public class VegetationLMA {
		this.alpha_loss = alpha_loss;
		this.alpha_offset = alpha_offset;
		this.alpha_lpf = alpha_lpf;
		this.terr_lpf = terr_lpf;
		this.veget_lpf = veget_lpf;
		final double [] scene_weights = setupSceneWeights(
		boost_parallax); // double boost_parallax)
		int min_scenes_uses = min_scenes;
		@@ -140,8 +144,19 @@ public class VegetationLMA {
		min_scenes_used, // int min_scenes_used,
		debugLevel > 3); // boolean debug_img){ // 4x?
		setupParametersIndices(valid_woi);
		alpha_neibs = getAlphaNeighbors();

		// alpha_neibs = getAlphaNeighbors();
		alpha_neibs = getNeighbors(
		TVAO_VEGETATION_ALPHA, // final int tvao, // TVAO_VEGETATION_ALPHA
		ind_pars_vegetation_alpha, // final int ind_samples, // ind_pars_vegetation_alpha
		num_pars_vegetation_alpha); // final int num_samples); // num_pars_vegetation_alpha
		terr_neibs = getNeighbors(
		TVAO_TERRAIN, // final int tvao, // TVAO_VEGETATION_ALPHA
		ind_pars_terrain, // final int ind_samples, // ind_pars_vegetation_alpha
		num_pars_terrain); // final int num_samples); // num_pars_vegetation_alpha
		veget_neibs = getNeighbors(
		TVAO_VEGETATION, // final int tvao, // TVAO_VEGETATION_ALPHA
		ind_pars_vegetation, // final int ind_samples, // ind_pars_vegetation_alpha
		num_pars_vegetation); // final int num_samples); // num_pars_vegetation_alpha

		if (!keep_parameters) {
		setupParametersVector (default_alpha); // areas where both terrain and vegetation are available
		@@ -162,6 +177,10 @@ public class VegetationLMA {
		return weights.length;
		}

		public double [] getParametersVector() {
		return parameters_vector;
		}

		private double [] setupSceneWeights(
		double boost_parallax) {
		// find max parallax
		@@ -490,6 +509,7 @@ public class VegetationLMA {
		if (debug_level > -2) { //
		String save_dir = debug_path;
		String debug_title = "parameters_vector-x"+woi.x+"-y"+woi.y+"-w"+woi.width+"-h"+woi.height;
		debug_title +="-al"+alpha_loss+"-alo"+alpha_offset+"-alp"+alpha_lpf+"-tl"+terr_lpf+"-vl"+veget_lpf;
		boolean save_all = (debug_image != null);
		boolean show_this = debug_level > 3;
		boolean show_all = debug_level > 4;
		@@ -633,6 +653,173 @@ public class VegetationLMA {
		return initial_rms;
		}

		public void updateFromParameters(
		double [] vector) {

		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() {
		for (int pindx = ai.getAndIncrement(); pindx < vector.length; pindx = ai.getAndIncrement()) {
		int [] indices = par_rindex[pindx];
		tvao[indices[0]][indices[1]] = vector[pindx];
		}
		}
		};
		}
		ImageDtt.startAndJoin(threads);
		}

		public void showResults(
		String title,
		double [] vector,
		double threshold_terrain,
		double min_max_terrain, //0.1
		double min_terrain, //0.001
		double min_vegetation) { // 0.5
		double [][] tvas = getTVASWoi(
		vector); // final double [] vector)
		double [][] terrain_weights_max = getTerrainWeights(
		threshold_terrain, // final double alpha_threshold, // discard images with too low transparency
		vector); // final double [] vector)
		double [] terrain_masked = tvas[0].clone();
		double [] vegetation_masked = tvas[0].clone();
		for (int i = 0; i < terrain_masked.length; i++) {
		if (!(terrain_weights_max[0][i] >= min_terrain)) terrain_masked[i] = Double.NaN;
		if (!(terrain_weights_max[1][i] >= min_max_terrain)) terrain_masked[i] = Double.NaN;
		if (!(tvas[2][i] >= min_vegetation)) vegetation_masked[i] = Double.NaN;
		}
		double [][] tvww = {terrain_masked, vegetation_masked, tvas[0], tvas[1], terrain_weights_max[0], terrain_weights_max[1],tvas[2]}; // no scene offsets
		String [] titles = {"terr>"+min_terrain,"veget>"+min_vegetation,"terr","veget","terr_str","terr_max","alpha"};

		ShowDoubleFloatArrays.showArrays(
		tvww,
		woi.width,
		woi.height,
		true,
		title,
		titles);
		return;
		}

		public double [][] getTVASWoi(
		final double [] vector){
		final int woi_length = woi.width*woi.height;
		double [][] tvas = new double [4][];
		for (int i = 0; i < 3; i++) tvas[i] = new double [woi_length];
		tvas[3] = new double [num_scenes];
		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() {
		for (int windx = ai.getAndIncrement(); windx < woi_length; windx = ai.getAndIncrement()) {
		int x = windx % woi.width + woi.x;
		int y = windx / woi.width + woi.y;
		int indx = x + y * full.width;
		for (int i = 0; i < 3; i++) {
		tvas[i][windx] = tvao[i][indx]; // non-adjusted parameters
		}
		}
		}
		};
		}
		ImageDtt.startAndJoin(threads);


		ai.set(0);
		for (int ithread = 0; ithread < threads.length; ithread++) {
		threads[ithread] = new Thread() {
		public void run() {
		for (int pindx = ai.getAndIncrement(); pindx < vector.length; pindx = ai.getAndIncrement()) {
		int [] indices = par_rindex[pindx];
		if (indices[0] <3) {
		int wx = indices[1] % full.width - woi.x;
		int wy = indices[1] / full.width - woi.y;
		int windx = wx + wy * woi.width;
		tvas[indices[0]][windx] = vector[pindx]; // all adjusted parameters, including offsets
		} else {
		tvas[indices[0]][ indices[1]] = vector[pindx]; // all adjusted parameters, including offsets
		}
		}
		}
		};
		}
		ImageDtt.startAndJoin(threads);
		return tvas;
		}


		public double [][] getTerrainWeights(
		final double alpha_threshold, // discard images with too low transparency
		final double [] vector)
		{
		final int woi_length = woi.width*woi.height;
		final double [][] terrain_scene_weights = new double [num_scenes][woi_length];
		final double [][] terrain_scene_max = new double [num_scenes][woi_length];
		final double [][] full_scene_weights = new double [num_scenes][woi_length];
		final double [][] terrain_weights = new double [2][woi_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 [] 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];
		int x = indx % full.width - woi.x;
		int y = indx / full.width - woi.y;
		int windx = x + y * woi.width;
		double [] cw = corners_weights[n];
		int [] indx_alpha = data_source[n][2];
		double sum_a =0;
		if (cw != null) {
		for (int i = 0; i < 4; i++ ) {
		int ia=indx_alpha[i];
		alpha[i] = cw[i] * ((ia >= 0) ? vector[ia]: tvao[TVAO_VEGETATION_ALPHA][-1-ia]);
		sum_a += alpha[i];
		}
		double tw = Math.max(0, 1.0 - sum_a);
		terrain_scene_max[nscene][windx] = tw;
		if (tw >= alpha_threshold) {
		terrain_scene_weights[nscene][windx] = weights[n] * tw;
		}
		full_scene_weights[nscene][windx] = weights[n];
		}
		}
		}
		};
		}
		ImageDtt.startAndJoin(threads);
		// sum for all scenes, divide (undefined - also 0)
		Arrays.fill(terrain_weights[0], Double.NaN);
		Arrays.fill(terrain_weights[1], Double.NaN);
		ai.set(0);
		for (int ithread = 0; ithread < threads.length; ithread++) {
		threads[ithread] = new Thread() {
		public void run() {
		for (int windx = ai.getAndIncrement(); windx < woi_length; windx = ai.getAndIncrement()) {
		double sum_w = 0, sum_wt = 0, max_t = 0;
		for (int nscene = 0; nscene < num_scenes;nscene++) {
		sum_wt += terrain_scene_weights[nscene][windx];
		sum_w += full_scene_weights[nscene][windx];
		max_t = Math.max(max_t, terrain_scene_max[nscene][windx]);
		}
		if (sum_w > 0) {
		terrain_weights[0][windx] = sum_wt/sum_w;
		terrain_weights[1][windx] = max_t;
		}
		}
		}
		};
		}
		ImageDtt.startAndJoin(threads);

		// regularization weights and derivatives
		return terrain_weights;
		}


		private double [] getFxDerivs(
		@@ -703,14 +890,18 @@ public class VegetationLMA {
		}
		ImageDtt.startAndJoin(threads);
		// regularization weights and derivatives
		ai.set(0);

		int ind_next = y_vector.length;
		if (alpha_lpf >= 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 + y_vector.length; // x - index
		int nx = n + ind_y_alpha; // y_vector.length; // x - index
		double d = 0;
		fX[nx] = 0.0;
		double alpha = vector[np];
		@@ -759,7 +950,99 @@ public class VegetationLMA {
		};
		}
		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);
		if (jt != null) {
		jt[np][nx] += terr_lpf;
		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);
		if (jt != null) {
		jt[np][nx] += veget_lpf;
		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;
		}

		@@ -977,8 +1260,14 @@ public class VegetationLMA {
		private void setupWeights( // after setupParametersIndices
		final double [] scene_weights,
		final double reg_weights) {
		int extra_samples = num_pars_vegetation_alpha; // in the future may be more regularization
		double reg_sample_weight = reg_weights/num_pars_vegetation_alpha; // weight of each regularization sample
		// int extra_samples = num_pars_vegetation_alpha; // in the future may be more regularization
		int extra_samples = 0;
		// using >=0 no use 0 as NOP but reserve space, <0 - do not reserve space
		if (alpha_lpf >= 0) extra_samples+= num_pars_vegetation_alpha;
		if (terr_lpf >= 0) extra_samples+= num_pars_terrain;
		if (veget_lpf >= 0) extra_samples+= num_pars_vegetation;

		double reg_sample_weight = reg_weights/extra_samples; // weight of each regularization sample
		final double [] sw = (scene_weights != null) ? scene_weights : new double [num_scenes];
		if (scene_weights == null) {
		Arrays.fill (sw, 1.0);
		@@ -1034,6 +1323,8 @@ public class VegetationLMA {
		ImageDtt.startAndJoin(threads);
		return;
		}

		@SuppressWarnings("unused")
		private int [][] getAlphaNeighbors(){
		final int [][] alpha_neibs = new int[num_pars_vegetation_alpha][4];
		final Thread[] threads = ImageDtt.newThreadArray(QuadCLT.THREADS_MAX);
		@@ -1066,6 +1357,45 @@ public class VegetationLMA {
		return alpha_neibs;
		}


		private int [][] getNeighbors(
		final int tvao, // TVAO_VEGETATION_ALPHA
		final int ind_samples, // ind_pars_vegetation_alpha
		final int num_samples // num_pars_vegetation_alpha
		){
		final int [][] alpha_neibs = new int[num_samples][4];
		final Thread[] threads = ImageDtt.newThreadArray(QuadCLT.THREADS_MAX);
		final AtomicInteger ai = new AtomicInteger(0);

		for (int ithread = 0; ithread < threads.length; ithread++) {
		threads[ithread] = new Thread() {
		TileNeibs tn = new TileNeibs(full.width, full.height);
		public void run() {
		for (int i = ai.getAndIncrement(); i < num_samples; i = ai.getAndIncrement()) {
		int pindx = ind_samples + i; // full parameter index
		int indx = par_rindex[pindx][1]; // full pixel index
		for (int dir4 = 0; dir4 < 4; dir4++) { // ortho only
		int dir = 2 * dir4;
		int nindx = tn.getNeibIndex(indx, dir);
		if (nindx < 0) {
		alpha_neibs[i][dir4] = -1;
		} else if (par_index[tvao][nindx] >= 0) {
		alpha_neibs[i][dir4] = par_index[tvao][nindx];
		} else { // pull to fixed alpha value
		alpha_neibs[i][dir4] = -2 -indx; // OFFSET by 2 !
		}
		}
		}
		}
		};
		}
		ImageDtt.startAndJoin(threads);
		//ind_pars_vegetation_alpha
		return alpha_neibs;
		}



		private void setupDataSource() {
		final int woi_length = woi.width*woi.height;
		final int [][][][] data_src = new int [num_scenes][woi_length][][];
		@@ -1252,7 +1582,7 @@ public class VegetationLMA {
		// int num_pars_pixel = par_num;
		ind_pars_scenes = par_num;
		for (int i = 0; i < par_index[TVAO_SCENE_OFFSET].length; i++) {
		par_rindex[par_num][0] = 1;
		par_rindex[par_num][0] = 3;
		par_rindex[par_num][1] = i;
		par_index[TVAO_SCENE_OFFSET][i] = par_num++;
		}

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

+28 −5

Original line number	Diff line number	Diff line
		@@ -698,14 +698,26 @@ 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 = 1000.0; // 100.; // 10.0; // quadratic loss when alpha reaches -1.0 or 2.0
		double alpha_offset = 0.0; // if >0, start losses above 0.0 and below 1.0;
		double alpha_lpf = 2.0; // 1.5; // 5.0; // 0.5; // pull to average of 4 neighbors
		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 = 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 boost_parallax = 5;
		boolean exit_loop = debugLevel < 1000;
		boolean next_run = false;
		boolean read_pars = true;
		String parameters_path = "/media/elphel/SSD3-4GB/lwir16-proc/berdich3/debug/vegetation/lma/parameters_vector_data.tiff";
		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;

		// 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";


		do {
		String par_path = (read_pars && !next_run)? parameters_path : null;
		int num_samples = vegetationLMA.prepareLMA(
		@@ -717,6 +729,8 @@ public class VegetationModel {
		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_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)
		boost_parallax, // final double boost_parallax, // increase weight of scene with maximal parallax relative to the reference scene
		par_path, // final String parameters_read_path,
		debugLevel); // final int debugLevel);
		@@ -724,6 +738,15 @@ public class VegetationModel {
		vegetationLMA.showYfX(
		null, // double [] vector,
		"reconstructed_model"); // String title)
		vegetationLMA.showResults(
		"terr_split", // String title,
		vegetationLMA.getParametersVector(), // double [] vector,
		threshold_terrain, // double threshold_terrain,
		min_max_terrain, // double min_max_terrain, //0.1
		min_terrain, //double min_terrain, //0.001
		min_vegetation); // double min_vegetation) { // 0.5


		}
		next_run = true;
		double lambda = 5.0; // 0.1;

Admin message