Cleaning up photometric equalization

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

@@ -176,7 +176,6 @@ public class IntersceneMatchParameters {
 	public  double  centroid_radius =               4.0;   // 
 	public  int     n_recenter =                    2;     // when cosine window, re-center window this many times
 	// filtering motion vectors
-	public  double  min_ref_str =                   0.22;  // For orientations: use only tiles of the reference scene DSI_MAIN is stronger  
 	// TD accumulation of the inter-scene correlations  demonstrated artifacts (horizontally offset by 8 pixels
 	// false maximum that is sharper than the real one. Still not understood - maybe float precision related. 
 	public  double  td_weight =                     0.5;   // mix correlations accumulated in TD with 
@@ -194,7 +193,18 @@ public class IntersceneMatchParameters {
 	public  double  half_disparity =                5.0;   // Reduce weight twice for this disparity
 	public  double  half_avg_diff =                 0.2;   // when L2 of x,y difference from average of neibs - reduce twice
 
+	// Photometric calibration (move elsewhere)?
+	public boolean  photo_en =                      false; // perform photogrammetric calibration to equalize pixel values
+	public int      photo_num_full =                1;     // Number of full recalibrations with re-processing of the images  
+	public int      photo_num_refines =             3;     // Calibrate, remove outliers, recalibrate, ... 
+	public double   photo_min_strength =            0.0;   // maybe add to filter out weak tiles
+	public double   photo_max_diff =                40.0;  // To filter mismatches. Normal (adjusted) have RMSE ~9
+	public boolean  photo_debug =                   false; // Generate images and text
+	
+	
+	
 	// Detect initial match
+	public  double  min_ref_str =                   0.22;  // For orientations: use only tiles of the reference scene DSI_MAIN is stronger  
 	public  int     pix_step =                      4;     // Azimuth/tilt search step in pixels
 	public  int     search_rad =                   10;     // Search radius in steps
 	public  double  maybe_sum =                     1.0;   // minimal sum of strengths (will search for the best)
@@ -593,9 +603,9 @@ public class IntersceneMatchParameters {
 		gd.addNumericField("Refine centroids",                       this.n_recenter, 0,5,"",
 				"Repeat centroids after moving the window center to the new centroid location this many times (0 - calculate once)");
 
-		gd.addMessage  ("Filter tiles to use for scene poses");
-		gd.addNumericField("DSI_MAIN minimal strength",              this.min_ref_str, 5,7,"",
-				"Match only tiles where DSI_MAIN is stronger than that (and has LMA).");
+	//	gd.addMessage  ("Filter tiles to use for scene poses");
+	//	gd.addNumericField("DSI_MAIN minimal strength",              this.min_ref_str, 5,7,"",
+	//			"Match only tiles where DSI_MAIN is stronger than that (and has LMA).");
 		
 		
 		gd.addMessage  ("Mixing TD and PD accumulation of 2d correlations");
@@ -625,7 +635,25 @@ public class IntersceneMatchParameters {
 		gd.addNumericField("Difference from neighbors average ",     this.half_avg_diff, 5,7,"",
 				"Reduce twice for high difference from neighbors average.");
 
+		gd.addMessage  ("Photometric calibration (move elsewhere?)");
+		gd.addCheckbox ("Enable photometric calibration",            this.photo_en,
+				"Equalize per- sensor gains and offsets. Requires disparity map. Save to reference scene and with current scene (to .corr-zml).");
+		gd.addNumericField("Full photometric (re)calibrations",      this.photo_num_full, 0,3,"pix",
+				"Full recalibratrions include re-importing raw images with updated offsets/gains");
+		gd.addNumericField("Refines",                                this.photo_num_refines, 0,3,"pix",
+				"Calculate calibration, remove outliers (e.g. FG/BG) and repeat");
+		gd.addNumericField("Minimal DSI strength",                   this.photo_min_strength, 5,7,"",
+				"Do not use weak tiles.");
+		gd.addNumericField("Maximal channel mismatch",               this.photo_max_diff, 5,7,"",
+				"Detect (and remove outliers). Adjusted images have RMSE ~9 counts.");
+		gd.addCheckbox ("Debug pphotometric calibration",            this.photo_debug,
+				"Generate debug images an text output.");
+		
+		
 		gd.addMessage  ("Initial search for the inter-scene match");
+		gd.addNumericField("DSI_MAIN minimal strength",              this.min_ref_str, 5,7,"",
+					"Match only tiles where DSI_MAIN is stronger than that (and has LMA).");
+		
 		gd.addNumericField("Azimuth/tilt step",                      this.pix_step, 0,3,"pix",
 				"Search in a spiral starting with no-shift with this step between probes, in approximate pixels");
 		gd.addNumericField("Search spiral radius",                   this.search_rad, 0,3,"steps",
@@ -931,7 +959,7 @@ public class IntersceneMatchParameters {
 		this.centroid_radius =          gd.getNextNumber();
 		this.n_recenter =         (int) gd.getNextNumber();
 
-		this.min_ref_str =              gd.getNextNumber();
+//		this.min_ref_str =              gd.getNextNumber();
 
 		this.td_weight =                gd.getNextNumber();
 		this.pd_weight =                gd.getNextNumber();
@@ -946,6 +974,15 @@ public class IntersceneMatchParameters {
 		this.weight_zero_neibs =        gd.getNextNumber();
 		this.half_disparity =           gd.getNextNumber();
 		this.half_avg_diff =            gd.getNextNumber();
+		
+		this.photo_en =                 gd.getNextBoolean();
+		this.photo_num_full =     (int) gd.getNextNumber();
+		this.photo_num_refines =  (int) gd.getNextNumber();
+		this.photo_min_strength =       gd.getNextNumber();
+		this.photo_max_diff =           gd.getNextNumber();
+		this.photo_debug =              gd.getNextBoolean();
+		
+		this.min_ref_str =              gd.getNextNumber();
 		this.pix_step =           (int) gd.getNextNumber();
 		this.search_rad =         (int) gd.getNextNumber();
 		this.maybe_sum =                gd.getNextNumber();
@@ -1157,6 +1194,14 @@ public class IntersceneMatchParameters {
 		properties.setProperty(prefix+"weight_zero_neibs",    this.weight_zero_neibs+"");   // double
 		properties.setProperty(prefix+"half_disparity",       this.half_disparity+"");      // double
 		properties.setProperty(prefix+"half_avg_diff",        this.half_avg_diff+"");       // double
+		
+		properties.setProperty(prefix+"photo_en",             this.photo_en+"");            // boolean
+		properties.setProperty(prefix+"photo_num_full",       this.photo_num_full+"");      // int
+		properties.setProperty(prefix+"photo_num_refines",    this.photo_num_refines+"");   // int
+		properties.setProperty(prefix+"photo_min_strength",   this.photo_min_strength+"");  // double
+		properties.setProperty(prefix+"photo_max_diff",       this.photo_max_diff+"");      // double
+		properties.setProperty(prefix+"photo_debug",          this.photo_debug+"");         // boolean
+		
 		properties.setProperty(prefix+"pix_step",             this.pix_step+"");            // int
 		properties.setProperty(prefix+"search_rad",           this.search_rad+"");          // int
 		properties.setProperty(prefix+"maybe_sum",            this.maybe_sum+"");           // double
@@ -1394,6 +1439,14 @@ public class IntersceneMatchParameters {
 		if (properties.getProperty(prefix+"weight_zero_neibs")!=null)    this.weight_zero_neibs=Double.parseDouble(properties.getProperty(prefix+"weight_zero_neibs"));
 		if (properties.getProperty(prefix+"half_disparity")!=null)       this.half_disparity=Double.parseDouble(properties.getProperty(prefix+"half_disparity"));
 		if (properties.getProperty(prefix+"half_avg_diff")!=null)        this.half_avg_diff=Double.parseDouble(properties.getProperty(prefix+"half_avg_diff"));
+		
+		if (properties.getProperty(prefix+"photo_en")!=null)             this.photo_en=Boolean.parseBoolean(properties.getProperty(prefix+"photo_en"));		
+		if (properties.getProperty(prefix+"photo_num_full")!=null)       this.photo_num_full=Integer.parseInt(properties.getProperty(prefix+"photo_num_full"));
+		if (properties.getProperty(prefix+"photo_num_refines")!=null)    this.photo_num_refines=Integer.parseInt(properties.getProperty(prefix+"photo_num_refines"));
+		if (properties.getProperty(prefix+"photo_min_strength")!=null)   this.photo_min_strength=Double.parseDouble(properties.getProperty(prefix+"photo_min_strength"));
+		if (properties.getProperty(prefix+"photo_max_diff")!=null)       this.photo_max_diff=Double.parseDouble(properties.getProperty(prefix+"photo_max_diff"));
+		if (properties.getProperty(prefix+"photo_debug")!=null)          this.photo_debug=Boolean.parseBoolean(properties.getProperty(prefix+"photo_debug"));		
+		
 		if (properties.getProperty(prefix+"pix_step")!=null)             this.pix_step=Integer.parseInt(properties.getProperty(prefix+"pix_step"));
 		if (properties.getProperty(prefix+"search_rad")!=null)           this.search_rad=Integer.parseInt(properties.getProperty(prefix+"search_rad"));
 		if (properties.getProperty(prefix+"maybe_sum")!=null)            this.maybe_sum=Double.parseDouble(properties.getProperty(prefix+"maybe_sum"));
@@ -1590,6 +1643,14 @@ public class IntersceneMatchParameters {
 		imp.weight_zero_neibs     = this.weight_zero_neibs;
 		imp.half_disparity        = this.half_disparity;
 		imp.half_avg_diff         = this.half_avg_diff;
+		
+		imp.photo_en              = this.photo_en;
+		imp.photo_num_full        = this.photo_num_full;
+		imp.photo_num_refines     = this.photo_num_refines;
+		imp.photo_min_strength    = this.photo_min_strength;
+		imp.photo_max_diff        = this.photo_max_diff;
+		imp.photo_debug           = this.photo_debug;
+		
 		imp.pix_step =              this.pix_step;
 		imp.search_rad =            this.search_rad;
 		imp.maybe_sum =             this.maybe_sum;

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

@@ -4045,7 +4045,14 @@ public class OpticalFlow {
     	boolean export_images =              clt_parameters.imp.export_images;
     	boolean export_dsi_image =           clt_parameters.imp.export_ranges;
     	boolean export_ml_files =            clt_parameters.imp.export_ml_files;
-    	boolean calibrate_photometric =      true;
+
+    	//boolean calibrate_photometric =      true;
+    	boolean  photo_en =                  clt_parameters.imp.photo_en; //          false; // perform photogrammetric calibration to equalize pixel values
+    	int      photo_num_full =            clt_parameters.imp.photo_num_full; //     1;    // Number of full recalibrations with re-processing of the images  
+    	int      photo_num_refines =         clt_parameters.imp.photo_num_refines; //  3;    // Calibrate, remove outliers, recalibrate, ... 
+    	double   photo_min_strength =        clt_parameters.imp.photo_min_strength; // 0.0;  // maybe add to filter out weak tiles
+    	double   photo_max_diff =            clt_parameters.imp.photo_max_diff; //    40.0;  // To filter mismatches. Normal (adjusted) have RMSE ~9
+    	boolean  photo_debug =               clt_parameters.imp.photo_debug; //       false; // Generate images and text
     	
     	boolean show_dsi_image =             clt_parameters.imp.show_ranges && !batch_mode;
     	boolean show_images =                clt_parameters.imp.show_images && !batch_mode;
@@ -4507,7 +4514,7 @@ public class OpticalFlow {
 			}
 		}
 
-		if (calibrate_photometric) {
+		if (photo_en && !reuse_video) {
 			if (combo_dsn_final == null) { // always re-read?
 				combo_dsn_final =quadCLTs[ref_index].readDoubleArrayFromModelDirectory( // always re-read?
 						"-INTER-INTRA-LMA", // String      suffix,
@@ -4521,45 +4528,46 @@ public class OpticalFlow {
 							combo_dsn_final,     // double [][]    combo_dsn_final, // dls,
 							quadCLTs[ref_index], // QuadCLT        scene,
 							debugLevel); // int            debugLevel);// > 0
-			QuadCLT.calibratePhotometric(
-					clt_parameters,           // CLTParameters     clt_parameters,
-					quadCLTs[ref_index],      // final QuadCLT     ref_scene, // now - may be null - for testing if scene is rotated ref
-					0.0,                      // final double      min_strength,
-					100.0,                    // final double      max_diff,    // 30.0
-					2,                        // final int         num_refines, // 2
-					combo_dsn_final_filtered, // final double [][] combo_dsn_final,     // double [][]    combo_dsn_final, // dls,			
-					threadsMax,               // int               threadsMax,
-					true);                    //final boolean     debug)
-			// copy offsets to the current to be saved with other properties. Is it correct/needed?
-			quadCLTs[ref_index].saveInterProperties( // save properties for interscene processing (extrinsics, ers, ...)
-					null, // String path,             // full name with extension or w/o path to use x3d directory
-					debugLevel+1);
-			quadCLT_main.setLwirOffsets(quadCLTs[ref_index].getLwirOffsets());
-			quadCLT_main.setLwirScales(quadCLTs[ref_index].getLwirScales());
-// Re-read reference and other scenes using new offsets			
-			quadCLTs[ref_index] = (QuadCLT) quadCLT_main.spawnQuadCLT( // restores dsi from "DSI-MAIN"
-					set_channels[ref_index].set_name,
-					clt_parameters,
-					colorProcParameters, //
-					threadsMax,
-					debugLevel);
-			quadCLTs[ref_index].setDSRBG(
-					clt_parameters, // CLTParameters  clt_parameters,
-					threadsMax,     // int            threadsMax,  // maximal number of threads to launch
-					updateStatus,   // boolean        updateStatus,
-					debugLevel);    // int            debugLevel)
-			ers_reference = quadCLTs[ref_index].getErsCorrection();
-			if (!reuse_video) { // reuse_video only uses reference scene
-				for (int scene_index =  ref_index - 1; scene_index >= earliest_scene ; scene_index--) {
-					//						quadCLTs[scene_index] = (QuadCLT) quadCLT_main.spawnNoModelQuadCLT( // restores image data
-					// to include ref scene photometric calibration
-					quadCLTs[scene_index] = quadCLTs[ref_index].spawnNoModelQuadCLT( // restores image data
-							set_channels[scene_index].set_name,
-							clt_parameters,
-							colorProcParameters, //
-							threadsMax,
-							debugLevel-2);
-				}
+			for (int nrecalib = 0; nrecalib < photo_num_full; nrecalib++) {
+				QuadCLT.calibratePhotometric(
+						clt_parameters,           // CLTParameters     clt_parameters,
+						quadCLTs[ref_index],      // final QuadCLT     ref_scene, // now - may be null - for testing if scene is rotated ref
+						photo_min_strength,       // final double      min_strength,
+						photo_max_diff,           // final double      max_diff,    // 30.0
+						photo_num_refines,        // final int         num_refines, // 2
+						combo_dsn_final_filtered, // final double [][] combo_dsn_final,     // double [][]    combo_dsn_final, // dls,			
+						threadsMax,               // int               threadsMax,
+						photo_debug);             //final boolean     debug)
+				// copy offsets to the current to be saved with other properties. Is it correct/needed?
+				quadCLTs[ref_index].saveInterProperties( // save properties for interscene processing (extrinsics, ers, ...)
+						null, // String path,             // full name with extension or w/o path to use x3d directory
+						debugLevel+1);
+				quadCLT_main.setLwirOffsets(quadCLTs[ref_index].getLwirOffsets());
+				quadCLT_main.setLwirScales(quadCLTs[ref_index].getLwirScales());
+				// Re-read reference and other scenes using new offsets			
+				quadCLTs[ref_index] = (QuadCLT) quadCLT_main.spawnQuadCLT( // restores dsi from "DSI-MAIN"
+						set_channels[ref_index].set_name,
+						clt_parameters,
+						colorProcParameters, //
+						threadsMax,
+						debugLevel);
+				quadCLTs[ref_index].setDSRBG(
+						clt_parameters, // CLTParameters  clt_parameters,
+						threadsMax,     // int            threadsMax,  // maximal number of threads to launch
+						updateStatus,   // boolean        updateStatus,
+						debugLevel);    // int            debugLevel)
+				ers_reference = quadCLTs[ref_index].getErsCorrection();
+			}
+			// propagate to other scenes in this sequence
+			for (int scene_index =  ref_index - 1; scene_index >= earliest_scene ; scene_index--) {
+				//						quadCLTs[scene_index] = (QuadCLT) quadCLT_main.spawnNoModelQuadCLT( // restores image data
+				// to include ref scene photometric calibration
+				quadCLTs[scene_index] = quadCLTs[ref_index].spawnNoModelQuadCLT( // restores image data
+						set_channels[scene_index].set_name,
+						clt_parameters,
+						colorProcParameters, //
+						threadsMax,
+						debugLevel-2);
 			}
 		}
 		

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

@@ -152,6 +152,7 @@ public class QuadCLTCPU {
     boolean                                                is_aux = false;
     double  []                                             lwir_offsets = null; // per image subtracted values
     double []                                              lwir_scales = null;  // per image scales
+    @Deprecated
     double                                                 lwir_offset =  Double.NaN; // average of lwir_offsets[]
     // hot and cold are calculated during autoranging (when generating 4 images for restored (added lwir_offset)
     // absolute temperatures to be used instead of colorProcParameters lwir_low and lwir_high if autoranging
@@ -1478,7 +1479,15 @@ public class QuadCLTCPU {
     public boolean isLwir()       {return geometryCorrection.isLwir();} // clt_kernels // USED in lwir
     public boolean isMonochrome() {return geometryCorrection.isMonochrome();} // clt_kernels // USED in lwir
     
-    public double  getLwirOffset() {return lwir_offset;} // USED in lwir
+    public double  getLwirOffset() {
+    	double s = 0;
+    	if (lwir_offsets == null) {
+    		throw new IllegalArgumentException ("Trying to use non-existing lwir_offsets!");
+    	}
+    	for (double offset:lwir_offsets) s+=offset;
+    	lwir_offset = s/lwir_offsets.length;
+    	return lwir_offset;
+    }
     
     public boolean isLwirCalibrated() {
     	return lwir_offsets != null;
@@ -5256,6 +5265,7 @@ public class QuadCLTCPU {
 								  set_name, // just for debug messages == setNames.get(nSet)
 								  threadsMax,
 								  debugLevel));
+				  // will actually is now calculated by setLwirOffsets()
 				  int num_avg = 0;
 				  this.lwir_offset = 0.0;
 				  for (int srcChannel=0; srcChannel < channelFiles.length; srcChannel++){