creating terrain 3d (removing vegetation)

src/main/java/com/elphel/imagej/orthomosaic/OrthoMap.java

@@ -3122,7 +3122,8 @@ public class OrthoMap implements Comparable <OrthoMap>, Serializable{
 
 		double avg_all = getMaskedAverage (
 				elev, // final double [] data,
-				null); // final boolean [] mask)
+				null, // final boolean [] mask)
+				null); // final double []  stdp);
 		double abs_high= avg_all + initial_above;
 		double abs_low=  avg_all - initial_below;
 		boolean [] mask = 	removeAbsoluteLowHigh (
@@ -3236,15 +3237,18 @@ public class OrthoMap implements Comparable <OrthoMap>, Serializable{
 	 * Calculate average of a double array, ignore masked out (!mask[i] and NaNs in the data)
 	 * @param data double data array
 	 * @param mask optional (may be null) boolean mask - ignore data items that have corresponding mask element false
+	 * @param stdp null or double[1] standard deviation
 	 * @return average value of the non-NaN and not disabled by mask elements of data[] 
 	 */
-	private static double getMaskedAverage (
-			final double [] data,
-			final boolean [] mask) {
+	public static double getMaskedAverage (
+			final double []  data,
+			final boolean [] mask,
+			final double []  stdp) {
 		final Thread[] threads = ImageDtt.newThreadArray();
 		final AtomicInteger ai = new AtomicInteger(0);
 		final AtomicInteger ati = new AtomicInteger(0);
 		final double [] avg_arr = new double [threads.length];
+		final double [] std_arr = new double [threads.length];
 		final double [] npix_arr = new double [threads.length];
 		for (int ithread = 0; ithread < threads.length; ithread++) {
 			threads[ithread] = new Thread() {
@@ -3252,18 +3256,25 @@ public class OrthoMap implements Comparable <OrthoMap>, Serializable{
 					int thread_num = ati.getAndIncrement();
 					for (int iPix = ai.getAndIncrement(); iPix < data.length; iPix = ai.getAndIncrement()) if (!Double.isNaN(data[iPix]) && ((mask == null) || (mask[iPix]))){
 						avg_arr[thread_num] += data[iPix];
+						std_arr[thread_num] += data[iPix]*data[iPix];
 						npix_arr[thread_num] += 1.0;
 					}
 				}
 			};
 		}		      
 		ImageDtt.startAndJoin(threads);
-		double avg=0, num=0;
+		double avg=0, num=0,sum2=0;
 		for (int i = 0; i < avg_arr.length; i++) {
 			avg+=avg_arr[i];
 			num+=npix_arr[i];
+			sum2+=std_arr[i];
 		}
-		return avg/num;
+		avg /= num;
+		sum2 /= num;
+		if (stdp != null) {
+			stdp[0] = Math.sqrt(sum2- avg*avg);
+		}
+		return avg;
 	}
 
 	/**
@@ -3310,7 +3321,7 @@ public class OrthoMap implements Comparable <OrthoMap>, Serializable{
 	 * @param threshold_low remove lower elements () and NaNs in the data
 	 * @return updated mask (same instance if not null) that disables (makes false) filtered out data [] elements
 	 */
-	private static boolean [] removeAbsoluteLowHigh (
+	public static boolean [] removeAbsoluteLowHigh (
 			final double [] data,
 			final boolean [] mask,
 			final double threshold_high,
@@ -3386,7 +3397,7 @@ public class OrthoMap implements Comparable <OrthoMap>, Serializable{
 	 * @param num_bins number of the histogram bins
 	 * @return boolean array of the remaining data elements. Input mask_in array (if not null) is modified too.
 	 */
-	private static boolean [] removeRelativeLowHigh (
+	public static boolean [] removeRelativeLowHigh (
 			final double [] data,
 			final boolean [] mask_in,
 			final double abs_high,

src/main/java/com/elphel/imagej/tileprocessor/IntersceneMatchParameters.java

@@ -320,8 +320,13 @@ public class IntersceneMatchParameters {
 	public  double  range_max =       5000.0;
 	// Export 3D model
 	public  boolean export3d =           false; // true;
-	
-	
+	public  boolean export3dterrain =    false; // true;	
+	// Generate CT scan
+	public  boolean export_CT =          false; // true;
+	public  double  ct_min =             3.7; //
+	public  double  ct_max =             4.5; //
+	public  double  ct_step =            0.02;  // pix
+	public  int     ct_expand =          2;    // pix // expand section area by 1 pixel in 8 directions
 	
 	// Debug and visualization
 	public  boolean scene_is_ref_test=              false; // correlate ref-2-ref for testing
@@ -1149,6 +1154,15 @@ min_str_neib_fpn	0.35
 		gd.addMessage  ("3D model generation");
 		gd.addCheckbox ("Generate 3D model",                         this.export3d,
 				"Generate textures and model.");
+		gd.addCheckbox ("Generate 3D model of terrain",              this.export3dterrain,
+				"Generate textures and model, remove vegetation.");
+		
+		gd.addMessage  ("CT scan generation");
+		gd.addCheckbox ("Generate ST scan",                          this.export_CT,	"Generate \"CT scan\" of the vegetation.");
+		gd.addNumericField("Minimal CT scan disparity",              this.ct_min, 5,7,"pix","Start disparity (lowest altitude).");
+		gd.addNumericField("Maximal CT scan disparity",              this.ct_max, 5,7,"pix","End disparity (highest altitude).");
+		gd.addNumericField("CT scan step",                           this.ct_step, 5,7,"pix","Scan step.");
+		gd.addNumericField("CT scan expand section",                 this.ct_expand, 0,3,"tiles","expand section area, 2 - a tile in 8 directions.");
 
 		gd.addMessage  ("Debug and visialization parameters");
 		gd.addCheckbox ("Replace scene with reference scene",        this.scene_is_ref_test,
@@ -2075,18 +2089,25 @@ min_str_neib_fpn	0.35
 		this.range_max =                      gd.getNextNumber();
 		
 		this.export3d =                       gd.getNextBoolean();
+		this.export3dterrain =                gd.getNextBoolean();
 		
-		this.scene_is_ref_test =        gd.getNextBoolean();
-		this.render_ref =               gd.getNextBoolean();
-		this.render_scene =             gd.getNextBoolean();
-		this.toRGB =                    gd.getNextBoolean();
-		this.show_2d_correlations =     gd.getNextBoolean();
-		this.show_motion_vectors =      gd.getNextBoolean();
-		this.debug_level =        (int) gd.getNextNumber();
+		this.export_CT =                      gd.getNextBoolean();
+		this.ct_min =                         gd.getNextNumber();
+		this.ct_max =                         gd.getNextNumber();
+		this.ct_step =                        gd.getNextNumber();
+		this.ct_expand =                (int) gd.getNextNumber();
+
+		this.scene_is_ref_test =              gd.getNextBoolean();
+		this.render_ref =                     gd.getNextBoolean();
+		this.render_scene =                   gd.getNextBoolean();
+		this.toRGB =                          gd.getNextBoolean();
+		this.show_2d_correlations =           gd.getNextBoolean();
+		this.show_motion_vectors =            gd.getNextBoolean();
+		this.debug_level =              (int) gd.getNextNumber();
 
-		this.test_ers =                 gd.getNextBoolean();
-		this.test_ers0 =          (int) gd.getNextNumber();
-		this.test_ers1 =          (int) gd.getNextNumber();
+		this.test_ers =                       gd.getNextBoolean();
+		this.test_ers0 =                (int) gd.getNextNumber();
+		this.test_ers1 =                (int) gd.getNextNumber();
 
 		this.num_bottom =               (int) gd.getNextNumber();
 		this.num_passes =               (int) gd.getNextNumber();
@@ -2684,6 +2705,13 @@ min_str_neib_fpn	0.35
 		properties.setProperty(prefix+"range_min_strength",            this.range_min_strength+"");            // double
 		properties.setProperty(prefix+"range_max",                     this.range_max+"");                     // double
 		properties.setProperty(prefix+"export3d",                      this.export3d+"");                      // boolean
+		properties.setProperty(prefix+"export3dterrain",               this.export3dterrain+"");               // boolean
+		
+		properties.setProperty(prefix+"export_CT",                     this.export_CT+"");                     // boolean
+		properties.setProperty(prefix+"ct_min",                        this.ct_min+"");                        // double
+		properties.setProperty(prefix+"ct_max",                        this.ct_max+"");                        // double
+		properties.setProperty(prefix+"ct_step",                       this.ct_step+"");                       // double
+		properties.setProperty(prefix+"ct_expand",                     this.ct_expand+"");                     // int
 		
 		properties.setProperty(prefix+"scene_is_ref_test",    this.scene_is_ref_test+"");   // boolean
 		properties.setProperty(prefix+"render_ref",           this.render_ref+"");          // boolean
@@ -3256,7 +3284,13 @@ min_str_neib_fpn	0.35
 		if (properties.getProperty(prefix+"range_max")!=null)            this.range_max=Double.parseDouble(properties.getProperty(prefix+"range_max"));
 
 		if (properties.getProperty(prefix+"export3d")!=null)             this.export3d=Boolean.parseBoolean(properties.getProperty(prefix+"export3d"));		
+		if (properties.getProperty(prefix+"export3dterrain")!=null)      this.export3dterrain=Boolean.parseBoolean(properties.getProperty(prefix+"export3dterrain"));		
 		
+		if (properties.getProperty(prefix+"export_CT")!=null)            this.export_CT=Boolean.parseBoolean(properties.getProperty(prefix+"export_CT"));		
+		if (properties.getProperty(prefix+"ct_min")!=null)               this.ct_min=Double.parseDouble(properties.getProperty(prefix+"ct_min"));
+		if (properties.getProperty(prefix+"ct_max")!=null)               this.ct_max=Double.parseDouble(properties.getProperty(prefix+"ct_max"));
+		if (properties.getProperty(prefix+"ct_step")!=null)              this.ct_step=Double.parseDouble(properties.getProperty(prefix+"ct_step"));
+		if (properties.getProperty(prefix+"ct_expand")!=null)            this.ct_expand=Integer.parseInt(properties.getProperty(prefix+"ct_expand"));
 		
 		if (properties.getProperty(prefix+"scene_is_ref_test")!=null)    this.scene_is_ref_test=Boolean.parseBoolean(properties.getProperty(prefix+"scene_is_ref_test"));		
 		if (properties.getProperty(prefix+"render_ref")!=null)           this.render_ref=Boolean.parseBoolean(properties.getProperty(prefix+"render_ref"));		
@@ -3834,8 +3868,13 @@ min_str_neib_fpn	0.35
 		imp.range_min_strength            = this.range_min_strength;
 		imp.range_max                     = this.range_max;
 		imp.export3d                      = this.export3d;
+		imp.export3dterrain               = this.export3dterrain;
 		
-		
+		imp.export_CT                     = this.export_CT;
+		imp.ct_min                        = this.ct_min;
+		imp.ct_max                        = this.ct_max;
+		imp.ct_step                       = this.ct_step;
+		imp.ct_expand                     = this.ct_expand;
 		
 		imp.scene_is_ref_test     = this.scene_is_ref_test;
 		imp.render_ref            = this.render_ref;

src/main/java/com/elphel/imagej/tileprocessor/OpticalFlow.java

@@ -4672,7 +4672,14 @@ public class OpticalFlow {
 		boolean generate_egomotion =         clt_parameters.imp.generate_egomotion; // generate egomotion table (image-based and ims)
 		boolean generate_mapped =            clt_parameters.imp.generate_mapped &&
 				(gen_seq_mono_color[0] || gen_seq_mono_color[1]);   // generate sequences - Tiff and/or video
-		boolean export3d =           clt_parameters.imp.export3d; //  true;
+		boolean export3d =                   clt_parameters.imp.export3d; //  true;
+		boolean export3dterrain =            clt_parameters.imp.export3dterrain; //  true;
+
+		boolean export_CT =                  clt_parameters.imp.export_CT; //  false;
+		double  ct_min =                     clt_parameters.imp.ct_min ;
+		double  ct_max =                     clt_parameters.imp.ct_max ;
+		double  ct_step =                    clt_parameters.imp.ct_step ;
+		int     ct_expand =                  clt_parameters.imp.ct_expand ;
 		
 		boolean [] annotate_mono_color = {clt_parameters.imp.annotate_mono,clt_parameters.imp.annotate_color}; 
 		boolean annotate_transparent_mono = clt_parameters.imp.annotate_transparent_mono;
@@ -6176,9 +6183,60 @@ public class OpticalFlow {
 	        }
 		}
 		
+		if (export_CT) {
+			double [][][] ct_scans = new double [(int) Math.ceil((ct_max-ct_min)/ct_step) + 1][][];
+			boolean ok_ct = TexturedModel.output3d( // quadCLTs have same image name, and everything else
+					clt_parameters,      // CLTParameters                            clt_parameters,
+					colorProcParameters, // ColorProcParameters                      colorProcParameters,
+					rgbParameters,       // EyesisCorrectionParameters.RGBParameters rgbParameters,
+					quadCLTs[ref_index], // quadCLTs.length-1], // quadCLT_main,        // final QuadCLT                            parameter_scene, // to use for rendering parameters in multi-series sequences
+		            // if null - use reference scene 
+					quadCLTs,            // QuadCLT []                               scenes,
+					ref_index,           // final int                                ref_index,
+					combo_dsn_final,     //double [][]                              combo_dsn_final, // null OK, will read file
+					ct_scans,            // final double [][]                        ct_scans,
+					ct_min,              // final double                             ct_min,
+					ct_step,             // final double                             ct_step,
+					ct_expand,           // final int                                ct_expand,
+					false,               // final boolean                            terrain_mode,
+					null, // final double [][]                        hdr_render_size, // { hdr_whs[3],hdr_x0y0[2]}; save/use rendering parameters
+					false, // final boolean                            hdr_render_slave,// use rendering parameters (to match other mode)
+					updateStatus,        // final boolean                            updateStatus,
+					debugLevel); //  + 1);         // final int                                debugLevel)
+			System.out.println ("CT scan-> "+ok_ct);
+			String [] titles = new String[ct_scans.length];
+			for (int i = 0; i < titles.length; i++) {
+				titles[i] = String.format("disparity=%.3f pix",ct_min+ct_step*i);
+			}
+			String suffix = String.format("-CT_SCAN_%.3f_%.3f_%.3f_%d",ct_min,ct_max,ct_step,ct_expand);
+			final int transform_size=           quadCLTs[ref_index].getTileProcessor().getTileSize();
+	 		final int tilesX =                  quadCLTs[ref_index].getTileProcessor().getTilesX();
+			final int tilesY =                  quadCLTs[ref_index].getTileProcessor().getTilesY();
+			int nslices = 0;
+			for (int n = 0; n < ct_scans.length; n++) {
+				nslices = Math.max(nslices,  ct_scans[n].length);
+			}
+			for (int nslice = 0; nslice < nslices; nslice++) {
+				double [][] ct_scan_slice = new double [ct_scans.length][];
+				for (int n = 0; n < ct_scan_slice.length; n++) {
+//					int i = Math.min(nslice, ct_scans[n].length);
+					if (nslice < ct_scans[n].length) {
+						ct_scan_slice[n] = ct_scans[n][nslice];
+					}
+				}
+				quadCLTs[ref_index].saveDoubleArrayInModelDirectory(
+						suffix+"-SLICE"+nslice,   // String      suffix,
+						titles,                   // String []   labels, // or null
+						ct_scan_slice,         // double [][] data,
+						tilesX * transform_size,  // int            width, // int tilesX,
+						tilesY * transform_size); // int            height, // int tilesY,
+			}
+		}
+		
 		// debugging 3D model
 		int [] whs = new int[3];
-		double [] x0y0 = new double[2];		
+		double [] x0y0 = new double[2];
+		double [][] hdr_render_size = new double[2][];
 		if (export3d) { //combo_dsn_final had strength 1.0e-4 where it should not? Reset it?
 			boolean ok_3d = TexturedModel.output3d( // quadCLTs have same image name, and everything else
 					clt_parameters,      // CLTParameters                            clt_parameters,
@@ -6189,9 +6247,38 @@ public class OpticalFlow {
 					quadCLTs,            // QuadCLT []                               scenes,
 					ref_index,           // final int                                ref_index,
 					combo_dsn_final,     //double [][]                              combo_dsn_final, // null OK, will read file
+					null,                // final double [][]                        ct_scans,
+					0,                   // final double                             ct_min,
+					0,                   // final double                             ct_step,
+					0,                   // final int                                ct_expand,
+					false,               // final boolean                            terrain_mode,
+					hdr_render_size,     // final double [][]                        hdr_render_size, // { hdr_whs[3],hdr_x0y0[2]}; save/use rendering parameters
+					false,               // final boolean                            hdr_render_slave,// use rendering parameters (to match other mode)
+					updateStatus,        // final boolean                            updateStatus,
+					debugLevel); //  + 1);         // final int                                debugLevel)
+			System.out.println ("TexturedModel.output3d() -> "+ok_3d+" (full, with vegetation)");
+		}
+		if (export3dterrain) { //combo_dsn_final had strength 1.0e-4 where it should not? Reset it?
+			boolean  hdr_render_slave= true;
+			boolean ok_3d = TexturedModel.output3d( // quadCLTs have same image name, and everything else
+					clt_parameters,      // CLTParameters                            clt_parameters,
+					colorProcParameters, // ColorProcParameters                      colorProcParameters,
+					rgbParameters,       // EyesisCorrectionParameters.RGBParameters rgbParameters,
+					quadCLTs[ref_index], // quadCLTs.length-1], // quadCLT_main,        // final QuadCLT                            parameter_scene, // to use for rendering parameters in multi-series sequences
+		            // if null - use reference scene 
+					quadCLTs,            // QuadCLT []                               scenes,
+					ref_index,           // final int                                ref_index,
+					combo_dsn_final,     //double [][]                              combo_dsn_final, // null OK, will read file
+					null,                // final double [][]                        ct_scans,
+					0,                   // final double                             ct_min,
+					0,                   // final double                             ct_step,
+					0,                   // final int                                ct_expand,
+					true,        // final boolean                            terrain_mode,
+					hdr_render_size,     // final double [][]                        hdr_render_size, // { hdr_whs[3],hdr_x0y0[2]}; save/use rendering parameters
+					hdr_render_slave,    // final boolean                            hdr_render_slave,// use rendering parameters (to match other mode)
 					updateStatus,        // final boolean                            updateStatus,
 					debugLevel); //  + 1);         // final int                                debugLevel)
-			System.out.println ("TexturedModel.output3d() -> "+ok_3d);
+			System.out.println ("TexturedModel.output3d() -> "+ok_3d+" (terrain only, no vegetation)");
 		}
 		
 		

src/main/java/com/elphel/imagej/tileprocessor/QuadCLTCPU.java

@@ -75,6 +75,7 @@ import com.elphel.imagej.ims.Did_ins_2;
 import com.elphel.imagej.ims.Did_pimu;
 import com.elphel.imagej.ims.EventLogger;
 import com.elphel.imagej.ims.Imx5;
+import com.elphel.imagej.orthomosaic.OrthoMap;
 import com.elphel.imagej.readers.ElphelTiffWriter;
 import com.elphel.imagej.readers.ImagejJp4Tiff;
 import com.elphel.imagej.x3d.export.TriMesh;
@@ -1706,6 +1707,85 @@ public class QuadCLTCPU {
 				tilesY);               // int         height)
     }
     
+    public double [] getFlatGround(
+    		double [] disparity,
+			double    rmse_above, // from average
+			double    rmse_below, // from average, // positive value
+			int       num_refine,
+			double    frac_above,
+			double    frac_below,
+			String    debug_title) {
+    	final double [] plane_disparity = new double [disparity.length];
+		final int     num_bins = 1024;
+    	final int tilesX =   tp.getTilesX();	
+    	final int tilesY =   tp.getTilesY();
+    	final double [] center_xy= {0.5*tilesX, 0.5*tilesY};
+    	final double [] stdp = new double[1];
+		double avg_all = OrthoMap.getMaskedAverage (
+				disparity, // final double [] data,
+				null, // final boolean [] mask)
+				stdp); // final double []  stdp);
+		double initial_above = rmse_above*stdp[0];
+		double initial_below = rmse_below*stdp[0];
+		double abs_high= avg_all + initial_below;
+		double abs_low=  avg_all - rmse_below*stdp[0];
+		boolean [] mask = 	OrthoMap.removeAbsoluteLowHigh (
+				disparity,     // final double [] data,
+				null,     // final boolean [] mask,
+				abs_high, // final double threshold_high,
+				abs_low); // final double threshold_low),
+		double [] plane_tiles = {0,0,avg_all,center_xy[0], center_xy[1]}; // 5 elements here
+		for (int ntry = 0; ntry < num_refine; ntry++) {
+			if (ntry > 0) {
+				// remove relatives, start from new mask
+				mask = OrthoMap.removeRelativeLowHigh (
+						disparity,       // final double [] data,
+						null,       // final boolean [] mask_in,
+						initial_above,   //  final double abs_high,
+						-initial_below,    // final double abs_low,
+						frac_above, // final double rhigh,
+						frac_below, // final double rlow,
+						plane_tiles,  // final double []  ground_plane, // tiltx,tilty, offs, x0(pix), y0(pix) or null
+						tilesX,      // final int     width, // only used with ground_plane != null;
+						num_bins);  // final int    num_bins) 				
+				
+			}
+			plane_tiles=  OrthoMap.getPlane(
+					disparity,     // final double []   data,
+					mask ,         // final boolean [] mask,
+					null,          // final double []  weight,
+					tilesX,        // final int        width,
+					center_xy);    // final double []  xy0)
+		}
+    	
+		for (int nTile = 0; nTile < disparity.length; nTile++) {
+			int tileX = nTile % tilesX;
+			int tileY = nTile / tilesX;
+			double x = tileX - plane_tiles[3];
+			double y = tileY - plane_tiles[4];
+			plane_disparity[nTile] = plane_tiles[2]+ x *plane_tiles[0] +y * plane_tiles[1];
+		}
+		if (debug_title != null) {
+			String [] dbg_titles = {"disparity", "plane", "diff"};
+			double [][] dbg_img = new double [dbg_titles.length][disparity.length];
+			dbg_img[0] = disparity;
+			dbg_img[1] = plane_disparity;
+			for (int i = 0; i < disparity.length; i++) {
+				dbg_img[2][i] = disparity[i] - plane_disparity[i];
+			}
+    		ShowDoubleFloatArrays.showArrays(
+    				dbg_img,
+    				tilesX,
+    				tilesY,
+    				true,
+    				debug_title,
+    				dbg_titles);
+		}
+    	return plane_disparity;
+    	
+    }
+    
+    
     
     public double [] getSmoothGround(
     		CLTParameters clt_parameters,

src/main/java/com/elphel/imagej/x3d/export/Render3D.java

@@ -78,12 +78,19 @@ public class Render3D {
 		ground_normal = ground_origin.subtract(
 				ground_x.scalarMultiply(ground_origin.dotProduct(ground_x))).subtract(
 						ground_y.scalarMultiply(ground_origin.dotProduct(ground_y)));
-		ground_normal_unit = ground_normal.normalize(); // unitary vector normal and away from the ground plane
+		Vector3D gnu; 
+		try {
+			gnu = ground_normal.normalize(); // unitary vector normal and away from the ground plane
+		} catch (Exception e) {
+			System.out.println("Zero above ground");
+			gnu = ground_x.crossProduct(ground_y).normalize(); // sets to (0, 0, 1)
+		}
+		ground_normal_unit = gnu; // only used for center projection where ground_origin != 0;
 		above_ground = ground_normal.getNorm();
 		xy_offs = new double[] {ground_x.dotProduct(ground_origin), ground_y.dotProduct(ground_origin)};
 	}
 	
-	public double [] projectToPlaneLinear(Vector3D v3) { // get ground plane pixel coordinate from camera x,y,z
+	public double [] projectToPlaneLinear(Vector3D v3) { // get ground plane pixel coordinate from camera x,y,z - not used
 		double   z =        ground_normal_unit.dotProduct(v3);
 		Vector3D in_plane = v3.scalarMultiply(above_ground/z);
 		double [] xy = new double[3];
@@ -93,7 +100,7 @@ public class Render3D {
 		return xy;
 	}
 	
-	public double [] projectToPlanePixels(Vector3D v3) { // get ground plane pixel coordinate from camera x,y,z
+	public double [] projectToPlanePixels(Vector3D v3) { // get ground plane pixel coordinate from camera x,y,z used in center proj
 		double   z =        ground_normal_unit.dotProduct(v3);
 		Vector3D in_plane = v3.scalarMultiply(above_ground/z);
 		double [] xy = new double[3];