Fixing small bugs (52ace80b) · Commits · Elphel / imagej-elphel

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

+197 −97

Original line number	Diff line number	Diff line
		@@ -283,7 +283,8 @@ public class VegetationLMA {
		private double [] last_ymfx = null;
		private double [][] last_jt = null;
		public double terrain_offset = Double.NaN; // average offset from the initial

		public boolean debug_source_vegetation = false;
		public int [][][] debug_source_veget = null;


		/**
		@@ -2490,7 +2491,7 @@ public class VegetationLMA {
		if (need_derivs) {
		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]) {
		for (int wvindex = 0; wvindex < woi_veg_length; wvindex++) if (valid_vegetation[wvindex]) { //
		int wvx = wvindex % woi_veg.width; // relative to woi_veg
		int wvy = wvindex / woi_veg.width; // relative to woi_veg
		int x = wvx + woi_veg.x; // relative to full
		@@ -2510,7 +2511,7 @@ public class VegetationLMA {
		}

		double [] scales = scales_xy[nScene][fnpix];
		double radius = elev_radius[fnpix];
		double radius = elev_radius[fnpix]; // *1.1;
		double px = x + scales[0] * elevation; // ELEV-SIGN
		double py = y + scales[1] * elevation; // ELEV-SIGN
		double px0 = px - radius, py0 = py - radius;
		@@ -3371,6 +3372,7 @@ public class VegetationLMA {
		terr_elev_radius, // final double radius,
		(jt != null), // final boolean calc_derivatives,
		debug_level); // final int debugLevel)
		terr_elev_woi4 = null; // do not try common terrain elevation
		} else {
		terr_elpix_woi4 = null;
		if (use_terr_elev) { //(getIndxTerrElev() >=0) && !fits_disable[TVAO_TERR_ELEV]) {
		@@ -3583,7 +3585,7 @@ public class VegetationLMA {
		} // if (woi_veg != null) {

		if (terr_elev_woi4 != null) { // common terrain elevation
		for (int wtindex = 0; wtindex < woi_terr_length; wtindex++) if (terr_elev_woi4[nScene][wtindex] != null) {
		for (int wtindex = 0; wtindex < woi_terr_length; wtindex++) if ((terr_elev_woi4[nScene] != null) &&(terr_elev_woi4[nScene][wtindex] != null)) {
		int wtx = wtindex % woi_terr.width; // relative to woi_terr
		int wty = wtindex / woi_terr.width; // relative to woi_terr
		int x = wtx + woi_terr.x; // relative to full
		@@ -3642,7 +3644,7 @@ public class VegetationLMA {
		} // if (terr_elev_woi4 != null) {

		if (terr_elpix_woi4 != null) { // per-pixel terrain elevation. Make mutually-exclusive with terr_elev_woi4
		for (int wtindex = 0; wtindex < woi_terr_length; wtindex++) if (terr_elpix_woi4[nScene][wtindex] != null) {
		for (int wtindex = 0; wtindex < woi_terr_length; wtindex++) if ((terr_elpix_woi4[nScene] != null) && (terr_elpix_woi4[nScene][wtindex] != null)) {
		int wtx = wtindex % woi_terr.width; // relative to woi_terr
		int wty = wtindex / woi_terr.width; // relative to woi_terr
		int x = wtx + woi_terr.x; // relative to full
		@@ -8491,12 +8493,31 @@ public class VegetationLMA {
		ImageDtt.startAndJoin(threads);
		Rectangle dbg_rectangle = new Rectangle(80,61,3,3);
		int dbg_scene = -81;
		// final boolean debug_source_veg = false; // this.debug_source_vegetation
		final AtomicInteger ati = new AtomicInteger(0);
		final ArrayList<ArrayList<ArrayList<HashSet<Integer>>>> aaahi = this.debug_source_vegetation ? (new ArrayList<ArrayList<ArrayList<HashSet<Integer>>>>()) : null;
		if (this.debug_source_vegetation) {
		for (int ithread = 0; ithread < threads.length; ithread++) {
		aaahi.add(new ArrayList<ArrayList<HashSet<Integer>>>(num_scenes));
		}
		}
		if (debug_source_vegetation) {
		ai.set(0);
		for (int ithread = 0; ithread < threads.length; ithread++) { // first sum for pairs
		threads[ithread] = new Thread() {
		public void run() {
		int nthread = ati.getAndIncrement();
		ArrayList<ArrayList<HashSet<Integer>>> aahi = null;
		if (debug_source_vegetation) {
		aahi = aaahi.get(nthread);
		for (int nscene = 0; nscene < num_scenes; nscene++) {
		ArrayList<HashSet<Integer>> ahi = new ArrayList<HashSet<Integer>>(woi_ext_len);
		for (int i = 0; i < woi_length; i++) {
		ahi.add(new HashSet<Integer>());
		}
		aahi.add(ahi);
		}
		}
		for (int indx = ai.getAndIncrement(); indx < woi_ext_len; indx = ai.getAndIncrement()) {
		int y = indx / woi_ext.width + woi_ext.y;
		int x = indx % woi_ext.width + woi_ext.x;
		@@ -8537,6 +8558,9 @@ public class VegetationLMA {
		if (valid_pixels[wpix]) {
		has_veg_thread[nthread][nscene][wpix] = true;
		has_valid_dependent = true;
		if (debug_source_vegetation) {
		aahi.get(nscene).get(wpix).add(npix);
		}
		}
		}
		}
		@@ -8554,7 +8578,72 @@ public class VegetationLMA {
		};
		}
		ImageDtt.startAndJoin(threads);

		} else {
		// final HashSet<Integer> [][][] dbg_src_veg = debug_source_veg ? (new HashSet<Integer>[threads.length][num_scenes][woi_length]): null;
		ai.set(0);
		for (int ithread = 0; ithread < threads.length; ithread++) { // first sum for pairs
		threads[ithread] = new Thread() {
		public void run() {
		int nthread = ati.getAndIncrement();
		// ArrayList<ArrayList<HashSet<Integer>>> aahi = debug_source_veg ?
		for (int indx = ai.getAndIncrement(); indx < woi_ext_len; indx = ai.getAndIncrement()) {
		int y = indx / woi_ext.width + woi_ext.y;
		int x = indx % woi_ext.width + woi_ext.x;
		int npix = x + y * full.width;
		double radius = elev_radius[npix];
		double radius2 = radius*radius;
		int iradius = (int) Math.ceil (radius);
		double elevation = tvao[TVAO_ELEVATION][npix];
		double vegetation = tvao[TVAO_VEGETATION][npix];
		if (!Double.isNaN(elevation) && !Double.isNaN(vegetation)) {
		for (int nscene = 0; nscene < num_scenes; nscene++) if (valid_scenes[nscene]){
		if ((nscene==dbg_scene) && dbg_rectangle.contains(x,y)) {
		System.out.println("dbg_scene="+dbg_scene+" x="+x+", y="+y);
		}
		// for (int nscene = 0; nscene < 62; nscene++) if (valid_scenes[nscene]){
		double [] scales = scales_xy[nscene][npix];
		if (scales != null) {
		double dx = scales[0] * elevation;
		double dy = scales[1] * elevation;
		int ax0 = (int) Math.floor(x + dx - radius);
		int ax1 = (int) Math.ceil (x + dx + radius);
		int ay0 = (int) Math.floor(y + dy - radius);
		int ay1 = (int) Math.ceil (y + dy + radius);
		Rectangle ar = new Rectangle(ax0, ay0, ax1-ax0 +1, ay1-ay0 +1);
		if (woi.intersects(ar)) {
		boolean has_valid_dependent = false;
		for (int ay = ay0; ay <= ay1; ay++) {
		double ddy = ay - (y + dy);
		double ddy2 = ddy*ddy;
		for (int ax = ax0; ax <= ax1; ax++) {
		if (woi.contains(ax,ay)) {
		double ddx = ax - (x + dx);
		// double ddr2 = ddy2 + ddx*ddx;
		if ((ddy2 <= radius2) && (ddx*ddx < radius2)) { // if (ddr2 <= radius2) {
		int wx = ax - woi.x;
		int wy = ay - woi.y;
		int wpix =wx + wy*woi.width;
		if (valid_pixels[wpix]) {
		has_veg_thread[nthread][nscene][wpix] = true;
		has_valid_dependent = true;
		}
		}
		}
		}
		}
		if (has_valid_dependent) {
		veg_scene_used_prefilter[indx] = true;
		}
		}
		}
		}
		}
		}
		}
		};
		}
		ImageDtt.startAndJoin(threads);
		}

		// combine from threads
		final boolean [][] has_veg = new boolean [num_scenes][woi_length];
		@@ -8577,6 +8666,35 @@ public class VegetationLMA {
		};
		}
		ImageDtt.startAndJoin(threads);
		if (debug_source_vegetation) {
		debug_source_veget = new int [num_scenes][][];
		ai.set(0);
		for (int ithread = 0; ithread < threads.length; ithread++) {
		threads[ithread] = new Thread() {
		public void run() {
		for (int nScene = ai.getAndIncrement(); nScene < num_scenes; nScene = ai.getAndIncrement()) if (valid_scenes[nScene]){
		debug_source_veget[nScene] = new int [woi_length][];
		for (int wpix = 0; wpix < woi_length; wpix++) {
		HashSet<Integer> hs = new HashSet<Integer>();
		for (int nthread =0; nthread < has_veg_thread.length; nthread++) {
		// ArrayList<HashSet<Integer>> ahi = aaahi.get(nthread).get(nScene); // ArrayList<ArrayList<HashSet<Integer>>> aahi = null;
		ArrayList<HashSet<Integer>> ahi = aaahi.get(nthread).get(nScene); // ArrayList<ArrayList<HashSet<Integer>>> aahi = null;
		if (!ahi.isEmpty()) {
		hs.addAll(ahi.get(wpix));
		if (!hs.isEmpty()) {
		// debug_source_veget[nScene][wpix] = new int[woi_length];
		debug_source_veget[nScene][wpix] = hs.stream().mapToInt(Integer::intValue).toArray();
		}
		}
		}
		}
		}
		}
		};
		}
		ImageDtt.startAndJoin(threads);
		}

		final int first_scene = afirst_scene.get();
		final int last_scene = alast_scene.get();
		// now similar to orange brunch make all scenes have the same pixels that have vegetation - remove "unstable" (border) pixels
		@@ -8628,14 +8746,12 @@ public class VegetationLMA {
		};
		}
		ImageDtt.startAndJoin(threads);
		// int num_removed = 0;
		if (anum_offsets.get() < anum_offsets_middle.get()) {
		System.out.println("getValidScenesPixels(): Had to remove some scenes that do not have offset data for pixels that have it in most scenes.");
		System.out.println("Number of pixels with offset data in most scenes: "+anum_offsets_middle.get()+", in all scenes: "+anum_offsets.get()+".");
		for (int nscene = 0; nscene < num_scenes; nscene++) if (valid_scenes[nscene]) {
		if (bad_offset_scenes[nscene].get()) {
		valid_scenes[nscene] = false;
		// num_removed++;
		System.out.println("Removing scene "+nscene+".");
		}
		}
		@@ -8654,6 +8770,9 @@ public class VegetationLMA {
		int y = indx / woi_ext.width + woi_ext.y;
		int x = indx % woi_ext.width + woi_ext.x;
		int npix = x + y * full.width;

		double radius = elev_radius[npix];
		double radius2 = radius*radius;
		double elevation = tvao[TVAO_ELEVATION][npix];
		double vegetation = tvao[TVAO_VEGETATION][npix];
		if (!Double.isNaN(elevation) && !Double.isNaN(vegetation)) {
		@@ -8662,17 +8781,23 @@ public class VegetationLMA {
		if (scales != null) {
		double dx = scales[0] * elevation;
		double dy = scales[1] * elevation;
		int px05 = (int) Math.floor(x + dx); // ELEV-SIGN
		int py05 = (int) Math.floor(y + dy); // ELEV-SIGN
		if (woi_plus1.contains(px05,py05)) { // up to 1 pixel to the nearest
		int ax0 = (int) Math.floor(x + dx - radius);
		int ax1 = (int) Math.ceil (x + dx + radius);
		int ay0 = (int) Math.floor(y + dy - radius);
		int ay1 = (int) Math.ceil (y + dy + radius);
		Rectangle ar = new Rectangle(ax0, ay0, ax1-ax0 +1, ay1-ay0 +1);
		if (woi.intersects(ar)) {
		// valid_pixels
		boolean has_valid_dependent = false;
		for (int yc = 0; yc <2; yc++) {
		int wy = py05+yc-woi.y;
		if ((wy >= 0) && (wy < woi.height)) {
		for (int xc = 0; xc < 2; xc++) {
		int wx = px05+xc-woi.x;
		if ((wx >=0) && (wx < woi.width)) {
		for (int ay = ay0; ay <= ay1; ay++) {
		double ddy = ay - (y + dy);
		double ddy2 = ddy*ddy;
		for (int ax = ax0; ax <= ax1; ax++) {
		if (woi.contains(ax,ay)) {
		double ddx = ax - (x + dx);
		if ((ddy2 <= radius2) && (ddx*ddx < radius2)) { // if (ddr2 <= radius2) {
		int wx = ax - woi.x;
		int wy = ay - woi.y;
		int wpix =wx + wy*woi.width;
		if (valid_pixels[wpix] && has_vegetation[wpix]) {
		has_valid_dependent = true;
		@@ -8701,10 +8826,6 @@ public class VegetationLMA {
		ImageDtt.startAndJoin(threads);

		Rectangle woi_veg_new = new Rectangle(amin_x.get(), amin_y.get(), amax_x.get()-amin_x.get()+1, amax_y.get()-amin_y.get()+1);
		// woi_veg_new.x=amin_x.get();
		// woi_veg_new.y=amin_y.get();
		// woi_veg_new.width= amax_x.get()-amin_x.get()+1;
		// woi_veg_new.height= amax_y.get()-amin_y.get()+1;
		boolean continue_no_veg = true;
		boolean terrain_only = (woi_veg_new.width <= 0) \|\| (woi_veg_new.width <= 0);
		if (terrain_only) {
		@@ -8779,7 +8900,6 @@ public class VegetationLMA {
		final boolean [] terr_scene_used_prefilter = new boolean [woi_terr_ext_len]; // which of the source terrain pixel may be used
		final boolean [] terr_scene_used = new boolean [woi_terr_ext_len]; // which of the source terrain pixel may be used
		final boolean [][][] has_terr_thread = new boolean [threads.length][num_scenes][woi_length];
		/// final boolean [][] has_terr_scene = new boolean [num_scenes][woi_length]; // which image pixels in woi have source terrain
		ai.set(0);
		ati.set(0);
		for (int ithread = 0; ithread < threads.length; ithread++) { // first sum for pairs
		@@ -8791,8 +8911,6 @@ public class VegetationLMA {
		int x = wtindx % woi_terr_ext.width + woi_terr_ext.x;
		int findx = x + y * full.width;
		double radius = elev_radius[findx]; // reuse same as for vegetation =1.5
		/// double radius2 = radius*radius;
		/// int iradius = (int) Math.ceil (radius);
		double terrain = tvao[TVAO_TERRAIN][findx];
		if (!Double.isNaN(terrain)) {
		for (int nscene = 0; nscene < num_scenes; nscene++) if (valid_scenes[nscene]){
		@@ -8845,24 +8963,6 @@ public class VegetationLMA {
		};
		}
		ImageDtt.startAndJoin(threads);
		/*
		// combine from threads - which image pixels in woi may have source terrain
		ai.set(0);
		for (int ithread = 0; ithread < threads.length; ithread++) {
		threads[ithread] = new Thread() {
		public void run() {
		for (int nScene = ai.getAndIncrement(); nScene < num_scenes; nScene = ai.getAndIncrement()) if (valid_scenes[nScene]){
		for (int wpix = 0; wpix < woi_length; wpix++) {
		for (int nthread =0; nthread < has_veg_thread.length; nthread++) {
		has_terr_scene[nScene][wpix] \|=has_terr_thread[nthread][nScene][wpix];
		}
		}
		}
		}
		};
		}
		ImageDtt.startAndJoin(threads);
		*/
		boolean [] has_terrain = new boolean [woi_length]; // which image pixels in woi have source terrain
		ai.set(0);
		for (int ithread = 0; ithread < threads.length; ithread++) {

Admin message