improved edge filtering, working version

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

@@ -6,6 +6,7 @@ import java.io.FileFilter;
 import java.util.Arrays;
 import java.util.Collections;
 import java.util.Comparator;
+import java.util.concurrent.atomic.AtomicBoolean;
 import java.util.concurrent.atomic.AtomicInteger;
 
 import com.elphel.imagej.common.DoubleGaussianBlur;
@@ -99,6 +100,7 @@ public class VegetationLMA {
 	
 //	private int               num_used_scenes; // -> num_par_scenes
 	private boolean []        used_scenes;    // encode unused as NaN
+	private boolean [][]      valid_woi = null; // transitional parameters - [0] - valid terrain and y_vector (per woi), [1] - valid vegetation, [2] - terrain has vegetation
 //	private int []            used_scenes_indices;
 	public final double []    scene_weights;
 	
@@ -466,20 +468,46 @@ public class VegetationLMA {
 		int     min_scenes_used = min_scenes * 4;
 		
 		boolean all_terrain = true; // use all terrain tiles even not used with vegetation
+		boolean [][] valid_woi;
+		int [] valid_stats = new int[4];
+		int            max_sacrifice = 5;
+		boolean [][] valid_scenes_pixels = getValidScenesPixels (
+				this.woi,                 // final Rectangle      woi,
+				this.full,                // final Rectangle      full,
+				this.terrain_rendered,    // final double [][]    terrain_rendered,    // terrain rendered for scenes (has nulls)
+				this.vegetation_offsets,  // final double [][][]  vegetation_offsets // [num_scenes][pixle]{dx,dy} differential offsets of vegetation to terrain grid
+				max_sacrifice,            // final int            max_sacrifice,       // may sacrifice first/last scenes to have the same valid_offset
+				valid_stats);             // final int []         stats // number of valid scenes, valid pixels if int[2]
+		if (debugLevel > -2) {
+			System.out.println("Got "+valid_stats[0]+" valid scenes (of "+this.terrain_rendered.length+
+					"), "+valid_stats[1]+" valid WOI pixels (of "+(this.woi.width*this.woi.height)+
+					"), of them "+valid_stats[3]+" have vegetation influence, valid vegetation pixels "+
+					valid_stats[2] + " (of "+(this.woi.width*this.woi.height)+").");
+		}
+		if (valid_scenes_pixels == null) {
+			System.out.println("Too few valid scenes/pixels.");
+			return -1;
+		}
+		this.used_scenes = valid_scenes_pixels[0];
 		
-		boolean [][] valid_woi = filterValidWoi (
+		if (clean_y) {
+			valid_woi = new boolean [][] {valid_scenes_pixels[1],valid_scenes_pixels[2],valid_scenes_pixels[3]};
+		} else { // will not work correctly
+			boolean [][] valid_terr_veg = filterValidWoi (
 				min_scenes_uses, // int     min_scenes_uses,
 				min_scenes_used, // int     min_scenes_used,
 				min_valid_pixels, // int     min_valid_pixels,
 				debugLevel > 3); // boolean debug_img){ // 4x?
-		if (valid_woi == null) {
+			if (valid_terr_veg == null) {
 				System.out.println("Too few valid pixels (< min_valid_pixels="+min_valid_pixels+").");
 				return -1;
 			}
-		
+			valid_woi = new boolean [][] {valid_terr_veg[0],valid_scenes_pixels[1],null}; // incorrect, should always be clean_y
+		}
+		this.valid_woi = valid_woi; // now [3]: valid terrain/y, valid vegetation, terrain has vegetation
 		setupParametersIndices(
 				valid_woi, // needs to know number of used scenes // all but scenes
-				all_terrain); // boolean all_terrain
+				!clean_y && all_terrain); // boolean all_terrain disable all_terrain when clean_y
 		
 		alpha_neibs = getNeighbors(
 				TVAO_ALPHA,      // final int tvao,        // TVAO_VEGETATION_ALPHA
@@ -3254,7 +3282,7 @@ public class VegetationLMA {
 		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()) {
+					for (int nScene = ai.getAndIncrement(); nScene < num_scenes; nScene = ai.getAndIncrement()) if (!clean_y || (used_scenes[nScene])){
 						for (int drow = 0; drow < woi.height; drow++) {
 							int row = woi.y + drow;
 							for (int dcol = 0; dcol < woi.width; dcol++) {
@@ -3263,7 +3291,9 @@ public class VegetationLMA {
 								int indx = row*full.width + col;
 								// now par_index[TVAO_TERRAIN][indx] may be -1 
 //								if (!Double.isNaN(terrain_rendered[nScene][indx]) && (par_index[TVAO_TERRAIN][indx] >= 0)) { //
-								if (!Double.isNaN(terrain_rendered[nScene][indx])) { //
+								boolean good_pixel = clean_y ? valid_woi[0][windx] : (!Double.isNaN(terrain_rendered[nScene][indx]));
+//								if (!Double.isNaN(terrain_rendered[nScene][indx])) { //
+								if (good_pixel) { //
 									boolean has_data = false;
 									int [][] data_record = new int[DATA_SOURCE_ITEMS][];
 									data_record[DATA_SOURCE_HEAD] = new int [DATA_SOURCE_HEAD_SIZE];
@@ -3278,7 +3308,7 @@ public class VegetationLMA {
 									data_record[DATA_SOURCE_HEAD][DATA_SOURCE_HEAD_SINDEX] = nScene;
 									// try building record, if no data - do not add it
 									double [] veg_xy = vegetation_offsets[nScene][indx];
-									if (veg_xy != null) {
+									if ((veg_xy != null) && (!clean_y || valid_woi[2][windx])) {
 										veg_xy = veg_xy.clone();
 										if (diff_offsets) {
 											veg_xy[0] += col;
@@ -3328,33 +3358,38 @@ public class VegetationLMA {
 			};
 		}		      
 		ImageDtt.startAndJoin(threads);
+
+		// should be no condensing in clean_y mode as it was pre-filtered
 		// now condense removing nulls in data_src (cornw may still have nulls)
 		int max_samples = 0;
 		for (int nscene = 0; nscene < num_scenes; nscene++) {
 			max_samples = Math.max(max_samples, scene_samples[nscene]);
 		}
 		System.out.println("Maximal samples in a scene = "+max_samples);
+
 		int allowed_less = 0;
 		int need_samples = max_samples-allowed_less;
 		int num_samples = 0;
-//		int dbg_early_sample = 0; // 10;
+		//		int dbg_early_sample = 0; // 10;
+		if (!clean_y) {
 			used_scenes = new boolean [num_scenes];
+		}
 		final int [] used_scene_indices = new int [num_scenes];
 		int num_used_scenes = 0;
 		for (int nscene = 0; nscene < num_scenes; nscene++) {
 			int num_prev = num_samples;
 			int defined = scene_samples[nscene];
-			boolean enough = defined >= need_samples;
+			boolean enough = clean_y ?  used_scenes[nscene] : (defined >= need_samples);
 			if (!enough) {
 				System.out.println ("Not enough samples in scene "+nscene+" ("+defined+" < "+need_samples+"), removing it.");
 			}
 			used_scene_indices[nscene] = num_used_scenes;
-//			if ((defined >= min_samples_scene) && (nscene >= dbg_early_sample)) {
+			//			if ((defined >= min_samples_scene) && (nscene >= dbg_early_sample)) {
 			if (enough) {
 				used_scenes[nscene] = true;
 				num_samples+= scene_samples[nscene];
 				scene_samples[nscene] = num_prev; // start index
-				if (keepScene(nscene)) {
+				if (clean_y || keepScene(nscene)) {
 					num_used_scenes++;
 				} else {
 					used_scene_indices[nscene] = -1;
@@ -3688,10 +3723,13 @@ public class VegetationLMA {
 					int windx =dcol + drow * woi.width;
 					if (valid_woi[1][windx]) {
 						int indx = row*full.width + col;
-						if (!Double.isNaN(vegetation_average[indx])) {
+//						if (!Double.isNaN(vegetation_average[indx])) { // should always be non-NaN
+						if (!Double.isNaN(tvao[TVAO_VEGETATION][indx])) { // should always be non-NaN
 							par_rindex[par_num][0] = TVAO_VEGETATION;
 							par_rindex[par_num][1] = indx;
 							par_index[TVAO_VEGETATION][indx] = par_num++;
+						} else if (clean_y){
+							System.out.println("Double.isNaN(vegetation_average["+indx+"] 1 - should not happen");
 						}
 					}
 				}
@@ -3707,10 +3745,12 @@ public class VegetationLMA {
 					int windx =dcol + drow * woi.width;
 					if (valid_woi[1][windx]) {
 						int indx = row*full.width + col;
-						if (!Double.isNaN(vegetation_average[indx])) {
+						if (!Double.isNaN(tvao[TVAO_ALPHA][indx])) { // vegetation_average[indx])) {
 							par_rindex[par_num][0] = TVAO_ALPHA;
 							par_rindex[par_num][1] = indx;
 							par_index[TVAO_ALPHA][indx] = par_num++;
+						} else if (clean_y){
+							System.out.println("Double.isNaN(vegetation_average["+indx+"] 2 - should not happen");
 						}
 					}
 				}
@@ -3726,10 +3766,12 @@ public class VegetationLMA {
 					int windx =dcol + drow * woi.width;
 					if (valid_woi[1][windx]) {
 						int indx = row*full.width + col;
-						if (!Double.isNaN(vegetation_average[indx])) {
+						if (!Double.isNaN(tvao[TVAO_ELEVATION][indx])) { //vegetation_average[indx])) {
 							par_rindex[par_num][0] = TVAO_ELEVATION;
 							par_rindex[par_num][1] = indx;
 							par_index[TVAO_ELEVATION][indx] = par_num++;
+						} else if (clean_y){
+							System.out.println("Double.isNaN(vegetation_average["+indx+"] 3 - should not happen");
 						}
 					}
 				}
@@ -3753,7 +3795,7 @@ public class VegetationLMA {
 		int par_num = ind_pars_scenes;
 		int [][] par_rindex = this.par_rindex;
 		for (int i = 0; i < par_index[TVAO_SCENE_OFFSET].length; i++) {
-			if (used_scenes[i] && keepScene(i)) {
+			if (used_scenes[i] && (clean_y || keepScene(i))) {
 				par_rindex[par_num][0] = TVAO_SCENE_OFFSET;
 				par_rindex[par_num][1] = i;
 				par_index[TVAO_SCENE_OFFSET][i] = par_num++;
@@ -3913,108 +3955,342 @@ public class VegetationLMA {
 	/**
 	 * Get valid scenes and woi pixels. Valid pixel is inside woi and does not depend on any
 	 * vegetation data outside of woi for all valid scenes. Valid scene, in turn, has all image
-	 * data for all of the valid pixels
+	 * data for all of the valid pixels.
+	 * 
+	 * The same pixel may have vegetation offset in some scenes, but not in the others (on the edge),
+	 * This causes problems in later LMA fitting, so all the remaining valid_pixels should either
+	 * have or do not have such links. This requires additional boolean array for terrain/y_vector
+	 * pixels that have vegetation offset data . 
+	 * 
+	 * Assuming vegetation image is defined everywhere where vegetation_offsets may point.
 	 * @param woi  WOI to be fully covered by the rendered data, usually this.woi
 	 * @param full full image WOI (0,0,640,512);
 	 * @param terrain_rendered  per-scene rendered images (for terrain elevation)
 	 * @param vegetation_offsets per-scene, per terrain pixel - x,y offsets to the corresponding vegetation pixel
-	 * @return 2d array - {valid scenes, valid pixels}
+	 * @return 2d array - {valid scenes, valid pixels, valid vegetation pixels, valid offset (vegetation influence)}
 	 */
 	public static boolean [][] getValidScenesPixels (
 			final Rectangle      woi,
 			final Rectangle      full,
 			final double [][]    terrain_rendered,    // terrain rendered for scenes (has nulls)
-			final double [][][]  vegetation_offsets // [num_scenes][pixle]{dx,dy} differential offsets of vegetation to terrain grid
+			final double [][][]  vegetation_offsets,  // [num_scenes][pixle]{dx,dy} differential offsets of vegetation to terrain grid
+			final int            max_sacrifice,       // may sacrifice first/last scenes to have the same valid_offset
+			final int []         stats // number of valid scenes, valid pixels if int[4] 
 			) {
 		final int num_scenes = terrain_rendered.length;
+		final int woi_length = woi.width*woi.height;
 		final boolean [] valid_scenes = new boolean [num_scenes]; 
-		final boolean [] valid_pixels = new boolean [woi.width*woi.height]; 
+//		final boolean [] valid_pixels = new boolean [woi_length];
+//		final boolean [] valid_vegetation = new boolean [woi_length];
+		final boolean [][] valid = new boolean [num_scenes][woi_length];
+		// The same pixel may have vegetation offset in some scenes, but not in the others (on the edge),
+		// This causes problems in later LMA fitting, so all the remaining valid_pixels
+		final int [] num_scenes_used = new int [woi_length];
 		final Thread[] threads = ImageDtt.newThreadArray(QuadCLT.THREADS_MAX);
 		final AtomicInteger ai = new AtomicInteger(0);
-		// first find pixels that are valid in at least one scene
+		final AtomicInteger amax_scenes =  new AtomicInteger(0);
 		for (int ithread = 0; ithread < threads.length; ithread++) {
 			threads[ithread] = new Thread() {
 				public void run() {
-					for (int wPix = ai.getAndIncrement(); wPix < valid_pixels.length; wPix = ai.getAndIncrement()) {
+					for (int wPix = ai.getAndIncrement(); wPix < woi_length; wPix = ai.getAndIncrement()) {
 						int x = wPix % woi.width + woi.x;
 						int y = wPix / woi.width + woi.y;
 						int npix = x + y*full.width;
-						check_pixel: {
 						for (int nscene = 0; nscene < num_scenes; nscene++) if (!Double.isNaN(terrain_rendered[nscene][npix])){
 							if (vegetation_offsets[nscene][npix] == null) {
-									valid_pixels[wPix] = true;
-									break check_pixel; // no vegetation, valid
-								}
+								valid[nscene][wPix] = true;
+								num_scenes_used[wPix]++;
+							} else {
 								int vx0 = (int) Math.floor(x + vegetation_offsets[nscene][npix][0]);
 								int vy0 = (int) Math.floor(y + vegetation_offsets[nscene][npix][1]);
 								Rectangle veg = new Rectangle(vx0,vy0,2,2);
-								if (!full.contains(veg)) {
-									continue; // to next scene. Maybe this scene will be removed later and do not disqualify this pixel.
+								if (woi.contains(veg)) { // assuming vegetation data is everywhere defined
+									valid[nscene][wPix] = true;
+									num_scenes_used[wPix]++;
+								}
+							}
+						}
+						amax_scenes.getAndAccumulate(num_scenes_used[wPix], Math::max);
+					}
+				}
+			};
+		}		      
+		ImageDtt.startAndJoin(threads);
+		// set valid pixels to the intersection of 
+		final int [][] start_end_scene = new int [woi_length][2];
+		ai.set(0);
+		for (int ithread = 0; ithread < threads.length; ithread++) {
+			threads[ithread] = new Thread() {
+				public void run() {
+//					int [] ladder = new int[num_scenes+1];
+					for (int wPix = ai.getAndIncrement(); wPix < woi_length; wPix = ai.getAndIncrement()) {
+						int x = wPix % woi.width + woi.x;
+						int y = wPix / woi.width + woi.y;
+						int npix = x + y*full.width;
+						int max_i = 0;
+						int max_l = 0;
+						int len =  0;
+						for (int i = 0; i < num_scenes; i++) {
+							if (valid[i][wPix]) {
+								len++;
+							} else {
+								if (len > max_l) {
+									max_l = len;
+									max_i = i - len; 
+								}
+								len = 0;
 							}
-								int vindx0 = vx0 + vy0 * full.width;
-								int [] vcorners = {vindx0, vindx0+1, vindx0+full.width, vindx0+full.width+1};
-								check_corners: {
-									for (int vindx:vcorners) {
-										if (Double.isNaN(terrain_rendered[nscene][vindx])) {
-											break check_corners; // no data for vegetation
 						}
+						if (len > max_l) {
+							max_l = len;
+							max_i = num_scenes - len; 
 						}
-									valid_pixels[wPix] = true;
-									break check_pixel; // no vegetation, valid
+						start_end_scene[wPix][0] = max_i;
+						start_end_scene[wPix][1] = max_i + max_l;
 					}
 				}
+			};
+		}		      
+		ImageDtt.startAndJoin(threads);
+		
+		final int [][] scene_pixel_product = new int [num_scenes][2]; // {max area, end}
+		final AtomicInteger amax_area =  new AtomicInteger(0);
+		ai.set(0);
+		for (int ithread = 0; ithread < threads.length; ithread++) {
+			threads[ithread] = new Thread() {
+				public void run() {
+					boolean [] selection = new boolean [woi_length];
+					for (int nScene = ai.getAndIncrement(); nScene < num_scenes; nScene = ai.getAndIncrement()) {
+						Arrays.fill(selection, true);
+						int max_area = 0;
+						int last_scene = 0;
+						for (int nscene1 = nScene; nscene1 < num_scenes; nscene1++) {
+							int num_pix = 0;
+							for (int wpix = 0; wpix < woi_length; wpix++) {
+								selection[wpix] &= valid[nscene1][wpix];
+								if (selection[wpix]) {
+									num_pix++;
+								}
+							}
+							int area = (nscene1 - nScene + 1) * num_pix;
+							if (area > max_area) {
+								max_area=area;
+								last_scene = nscene1;
+							}
+							if (area == 0) {
+								break; // no more intersections
+							}
+						}
+						scene_pixel_product[nScene][0] = max_area;
+						scene_pixel_product[nScene][1] = last_scene + 1; // inclusive-> exclusive
+						amax_area.getAndAccumulate(scene_pixel_product[nScene][0], Math::max);
+					}
+				}
+			};
+		}		      
+		ImageDtt.startAndJoin(threads);
+		
+		
+		int max_area = amax_area.get();
+		int scene0 = -1, scene1 = -1;
+		for (int nscene = 0; nscene < num_scenes; nscene++) {
+			if (scene_pixel_product[nscene][0] == max_area) {
+				scene0 = nscene;
+				scene1 = scene_pixel_product[nscene][1];
+				break;
+			}
+		}
+		System.out.println ("scene0="+scene0+", scene1="+scene1);
+		// AND all valid_pixels for the scene range
+		final AtomicInteger anum_pix = new AtomicInteger(0);
+		final boolean [] valid_pixels = new boolean [woi_length];
+		ai.set(0);
+		final int fscene0=scene0, fscene1 = scene1;
+		for (int ithread = 0; ithread < threads.length; ithread++) {
+			threads[ithread] = new Thread() {
+				public void run() {
+					for (int nPix = ai.getAndIncrement(); nPix < woi_length; nPix = ai.getAndIncrement()) {
+						valid_pixels[nPix] = true;
+						for (int nscene = fscene0; nscene < fscene1; nscene++) {
+							valid_pixels[nPix] &= valid[nscene][nPix];
+						}
+						if (valid_pixels[nPix]) {
+							anum_pix.getAndIncrement();
 						}
 					}
 				}
 			};
 		}		      
 		ImageDtt.startAndJoin(threads);
-		// now remove scenes that have at least one of the valid pixels missing
+		// Now mark selected and try to add more scenes if they contain all valid pixels
+		final AtomicInteger anum_scenes = new AtomicInteger(0);
 		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()) {
-						double []   terrain = terrain_rendered[nScene];
-						double [][] offsets = vegetation_offsets[nScene];
-						check_scene: {
-							for (int wy = 0; wy < woi.height; wy++) {
-								int y = wy+woi.y;
-								for (int wx = 0; wx < woi.width; wx++) {
-									int wpix = wx + wy * woi.width;
-									if (valid_pixels[wpix]) {
-										int x = wx+woi.x;
-										//valid_pixels[npix]
-										int indx = x+ y*full.width;
-										if (Double.isNaN(terrain[indx])) {
-											break check_scene;
-										}
-										if (offsets[indx] != null) { // null is OK - terrain only
-											int vx0 = (int) Math.floor(x + offsets[indx][0]);
-											int vy0 = (int) Math.floor(y + offsets[indx][1]);
+						if ((nScene >= fscene0) && (nScene < fscene1)) {
+							valid_scenes[nScene] = true;
+							anum_scenes.getAndIncrement();
+						} else {
+							check_all_pixels: {
+								for (int npix = 0; npix < woi_length; npix++) {
+									if (!valid[nScene][npix]) {
+										break check_all_pixels;
+									}
+								}
+								valid_scenes[nScene] = true; // new added scene
+								anum_scenes.getAndIncrement();
+								System.out.println("getValidScenesPixels(): added scene "+nScene+" that has all valid pixels");
+							}
+						}
+					}
+				}
+			};
+		}		      
+		ImageDtt.startAndJoin(threads);
+
+		if (stats != null) {
+			stats[0] = anum_scenes.get();
+			stats[1] = anum_pix.get();
+			stats[2] = 0;
+			stats[3] = 0;
+		}
+		if ((anum_pix.get() == 0) || (anum_scenes.get()==0)) {
+			return null;
+		}
+		
+		final AtomicInteger ati = new AtomicInteger(0);
+		final AtomicInteger afirst_scene = new AtomicInteger(num_scenes);
+		final AtomicInteger alast_scene =  new AtomicInteger(0);
+		final boolean [][]  valid_vegetation_threaded = new boolean [threads.length][woi_length];
+		ai.set(0);
+		for (int ithread = 0; ithread < threads.length; ithread++) {
+			threads[ithread] = new Thread() {
+				public void run() {
+					int nthread = ati.getAndIncrement();
+					for (int nScene = ai.getAndIncrement(); nScene < num_scenes; nScene = ai.getAndIncrement()) if (valid_scenes[nScene]){
+						afirst_scene.getAndAccumulate(nScene, Math::min);
+						alast_scene.getAndAccumulate(nScene, Math::max);
+						for (int wpix = 0; wpix < woi_length; wpix++) if (valid_pixels[wpix]) {
+							int x = wpix % woi.width + woi.x;
+							int y = wpix / woi.width + woi.y;
+							int npix = x + y*full.width;
+							if (vegetation_offsets[nScene][npix] != null) {
+								int vx0 = (int) Math.floor(x + vegetation_offsets[nScene][npix][0]);
+								int vy0 = (int) Math.floor(y + vegetation_offsets[nScene][npix][1]);
 								Rectangle veg = new Rectangle(vx0,vy0,2,2);
-											if (!full.contains(veg)) {
-												break check_scene; // required data does not fit in the image (woi too close to the edge?)
+								if (!woi.contains(veg)) { // assuming vegetation data is everywhere defined
+									throw new IllegalArgumentException(
+											"getValidScenesPixels() bug: vegetation quad does not fit: scene="+nScene+", x="+x+", y="+y+
+											", vx0="+vx0+", vy0="+vy0+", woi=("+woi.x+", "+woi.y+", "+woi.width+", woi.height");
 								}
-											int vindx0 = vx0 + vy0 * full.width;
-											int [] vcorners = {vindx0, vindx0+1, vindx0+full.width, vindx0+full.width+1};
+								int wvindx0 = (vx0-woi.x) + (vy0-woi.y) * woi.width;
+								int [] vcorners = {wvindx0, wvindx0+1, wvindx0+woi.width, wvindx0+woi.width+1};
 								for (int vindx:vcorners) {
-												if (Double.isNaN(terrain[vindx])) {
-													break check_scene; // no data for vegetation
+									valid_vegetation_threaded[nthread][vindx] =true;
 								}
 							}
 						}
 					}
 				}
+			};
+		}		      
+		ImageDtt.startAndJoin(threads);
+		final int first_scene = afirst_scene.get();
+		final int last_scene =  alast_scene.get();
+
+// combine
+		final boolean [] valid_vegetation = new boolean [woi_length];
+		final AtomicInteger anum_veget = new AtomicInteger(0);
+		ai.set(0);
+		for (int ithread = 0; ithread < threads.length; ithread++) {
+			threads[ithread] = new Thread() {
+				public void run() {
+					for (int nPix = ai.getAndIncrement(); nPix < woi_length; nPix = ai.getAndIncrement()) {
+						check_valid: {
+							for (int nthread = 0; nthread < valid_vegetation_threaded.length; nthread++) {
+								if  (valid_vegetation_threaded[nthread][nPix]) {
+									valid_vegetation[nPix] = true;
+									anum_veget.getAndIncrement();
+									break check_valid;
+								}
+
 							}
-							valid_scenes[nScene] = true;
 						}
 					}
 				}
 			};
 		}		      
 		ImageDtt.startAndJoin(threads);
-		return new boolean[][] {valid_scenes, valid_pixels};
+		
+		
+		final AtomicInteger anum_offsets = new AtomicInteger(0);
+		final boolean [] valid_offset =    new boolean [woi_length];
+		final AtomicInteger anum_offsets_middle = new AtomicInteger(0);
+		final AtomicBoolean [] bad_offset_scenes = new AtomicBoolean[num_scenes];
+		for (int nscene = first_scene; nscene <= last_scene; nscene++) {
+			bad_offset_scenes[nscene] = new AtomicBoolean(false);
+		}
+		ai.set(0);
+		for (int ithread = 0; ithread < threads.length; ithread++) {
+			threads[ithread] = new Thread() {
+				public void run() {
+					for (int wPix = ai.getAndIncrement(); wPix < woi_length; wPix = ai.getAndIncrement()) if (valid_pixels[wPix]){
+						int x = wPix % woi.width + woi.x;
+						int y = wPix / woi.width + woi.y;
+						int npix = x + y*full.width;
+						boolean has_offset_middle = true;
+						for (int nscene = first_scene + max_sacrifice; nscene <= last_scene-max_sacrifice; nscene++) if (valid_scenes[nscene]){
+							if (vegetation_offsets[nscene][npix] == null) {
+								has_offset_middle = false;
+								break;
+							}
+						}
+						if (has_offset_middle) {
+							boolean has_offset = true;
+							for (int nscene = first_scene; nscene < first_scene + max_sacrifice; nscene++) if (valid_scenes[nscene]){
+								if (vegetation_offsets[nscene][npix] == null) {
+									has_offset = false;
+									bad_offset_scenes[nscene].set(true);
+								}
+							}
+							for (int nscene = last_scene-max_sacrifice + 1; nscene <= last_scene; nscene++) if (valid_scenes[nscene]){
+								if (vegetation_offsets[nscene][npix] == null) {
+									has_offset = false;
+									bad_offset_scenes[nscene].set(true);
+								}
+							}
+							valid_offset[wPix] = true;
+							anum_offsets_middle.getAndIncrement();
+							if (has_offset) {
+								anum_offsets.getAndIncrement();
+							}
+						}
+					}
+				}
+			};
+		}		      
+		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+".");
+				}
+			}
+		}
+		if (stats != null) {
+			stats[0] -= num_removed;
+			stats[2] = anum_veget.get();
+			stats[3] = anum_offsets.get();
+		}
+		return new boolean[][] {valid_scenes, valid_pixels, valid_vegetation, valid_offset};
 	}