working on motion vectors strength equalization

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

@@ -2802,6 +2802,19 @@ public class ImageDtt extends ImageDttCPU {
 							}
 							if (half_disparity > 0.0) { // scale by disparity
 								mv[nTile][2] *= half_disparity/(half_disparity + pxd[nTile][2]);
+								//							} else { // temporary, testing, make equalization?
+								//								double disp0 = 5.0;
+								//								if (pxd[nTile][2] > disp0) {
+								//									mv[nTile][2] *= pxd[nTile][2]/disp0;
+								//								}
+								/*
+							} else {
+								double disp0 = 5.0;
+								if (pxd[nTile][2] > disp0) {
+									mv[nTile][2] *= 10; // pxd[nTile][2]/disp0;
+								}
+								 */
+
 							}
 							if ((half_avg_diff > 0.0) &&(num_neibs > 0)) {
 								double dx = mv[nTile][0] - sx / num_neibs;

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

@@ -185,8 +185,8 @@ public class IntersceneMatchParameters {
 	public  double  min_str_fpn =                   0.2; // 0.25; // minimal correlation strength for all but TD-accumulated layer
 	public  double  min_str_sum_fpn =               0.5; // 0.8;  // minimal correlation strength for TD-accumulated layer
 
-	public  double  min_str =                       0.18;  // tiles w/o FPN: minimal correlation strength for all but TD-accumulated layer
-	public  double  min_str_sum =                   0.33; // tiles w/o FPN: minimal correlation strength for TD-accumulated layer
+	public  double  min_str =                       0.12; //18;  // tiles w/o FPN: minimal correlation strength for all but TD-accumulated layer
+	public  double  min_str_sum =                   0.2;  // 0.33; // tiles w/o FPN: minimal correlation strength for TD-accumulated layer
 
 	public  int     min_neibs =                     2;	   // minimal number of strong neighbors (> min_str)
 	public  double  weight_zero_neibs =             0.2;   // Reduce weight for no-neib (1.0 for all 8)
@@ -194,7 +194,7 @@ public class IntersceneMatchParameters {
 	public  double  half_avg_diff =                 0.2;   // when L2 of x,y difference from average of neibs - reduce twice
 	
 	// Detect initial match
-	public  double  min_ref_str =                   0.22;  // For orientations: use only tiles of the reference scene DSI_MAIN is stronger  
+	public  double  min_ref_str =                   0.15;  // 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)

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

@@ -4382,7 +4382,7 @@ public class OpticalFlow {
 						debugLevel-2);
 			} // split cycles to remove output clutter
 			int debug_scene = -15;
-			boolean debug2 = false; // true;
+			boolean debug2 = !batch_mode; // false; // true;
 			boolean [] reliable_ref = null;
 			if (min_ref_str > 0.0) {
 				reliable_ref = quadCLTs[ref_index].getReliableTiles( // will be null if does not exist.
@@ -4431,9 +4431,26 @@ public class OpticalFlow {
 					}
 					scenes_xyzatr[scene_index] = new double [][] {new double[3], use_atr};
 				} else { // assume linear motion
+					int num_avg = 1; // 3;
+					double scale_xyz = 0.0; // 0.5;
+					int na = num_avg;
+					if ((scene_index + 1 + na) > ref_index) {
+						na = ref_index - (scene_index + 1); 
+					}
+					
 					double [][] last_diff = ErsCorrection.combineXYZATR(
-							scenes_xyzatr[scene_index+1],
-							ErsCorrection.invertXYZATR(scenes_xyzatr[scene_index+2]));
+							scenes_xyzatr[scene_index + 1],
+							ErsCorrection.invertXYZATR(scenes_xyzatr[scene_index+1 + na]));
+					for (int i = 0; i < 3; i++) {
+						last_diff[0][i] /= na;
+						last_diff[1][i] /= na;
+					}
+					for (int i = 0; i < 3; i++) {
+						last_diff[0][i] *= scale_xyz;
+					}					
+//					last_diff[0][0] = 0.0;
+//					last_diff[0][1] = 0.0;
+//					last_diff[0][2] = 0.0;
 					scenes_xyzatr[scene_index] = ErsCorrection.combineXYZATR(
 							scenes_xyzatr[scene_index+1],
 							last_diff);
@@ -4606,7 +4623,7 @@ public class OpticalFlow {
 							quadCLTs[ref_index], // QuadCLT        scene,
 							debugLevel); // int            debugLevel);// > 0
 			for (int nrecalib = 0; nrecalib < photo_num_full; nrecalib++) {
-				QuadCLT.calibratePhotometric(
+				QuadCLT.calibratePhotometric2(
 						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,
@@ -4620,7 +4637,8 @@ public class OpticalFlow {
 						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());
+				quadCLT_main.setLwirScales (quadCLTs[ref_index].getLwirScales ());
+				quadCLT_main.setLwirScales2(quadCLTs[ref_index].getLwirScales2());
 				// Re-read reference and other scenes using new offsets	
 				quadCLTs[ref_index].saveQuadClt(); // to re-load new set of Bayer images to the GPU (do nothing for CPU) and Geometry
 				quadCLTs[ref_index] = (QuadCLT) quadCLT_main.spawnQuadCLT( // restores dsi from "DSI-MAIN"
@@ -12478,6 +12496,242 @@ public double[][] correlateIntersceneDebug( // only uses GPU and quad
 		return coord_motion;
 	}
 	
+	/**
+	 * Equalize weights of the motion vectors to boost that of important buy weak one.
+	 * Process overlapping (by half, using shifted cosine weight function) supertiles
+	 * independetly and if it qualifies, increase its tiles weights equlaizing (with
+	 * certain limitations) total supertile weights.
+	 * @param coord_motion       [2][tilesX*tilesY][3] input/output arrays.[0][tile][] is
+	 *                           pXpYD triplet (
+	 * @param stride_hor         half of a supertile width
+	 * @param stride_vert        half of a supertile height
+	 * @param min_stile_weight   minimal total weight of the tiles in a supertile (lower
+	 *                           will not be modified)
+	 * @param min_stile_number   minimal number of defined tiles in a supertile
+	 * @param min_stile_fraction minimal total tile strength compared to the average one
+	 * @param min_disparity      minimal disparity of tiles to consider (after filtering)
+	 * @param max_disparity      maximal disparity of tiles to consider (after filtering)
+	 * @param weight_add         add to each tile after scaling (if total weight < average)
+	 * @param weight_scale       scale each tile (if total weight < average)
+	 * If total new weight of a supertile exceeds average - scale each tile to match. If
+	 * lower - keep as is. Only after this step remove tiles (replace with original weight)
+	 * that are discarded by the disparity filter. Multiply result by the the window and
+	 * accumulate (in 4 passes to prevent contentions for the same destination array  
+	 */
+	public void equalizeMotionVectorsWeights(
+			final double [][][] coord_motion,
+			final int           tilesX,
+			final int           stride_hor,
+			final int           stride_vert,
+			final double        min_stile_weight,
+			final int           min_stile_number,
+			final double        min_stile_fraction,
+			final double        min_disparity,
+			final double        max_disparity,
+			final double        weight_add,
+			final double        weight_scale)
+	{
+		final int tiles = coord_motion[0].length;
+		final int tilesY = tiles/tilesX;
+		final int stile_width =  2 * stride_hor;
+		final int stile_height = 2 * stride_vert;
+		double [] whor =  new double [stride_hor];
+		double [] wvert = new double [stride_vert];
+		for (int i = 0; i < stride_hor;  i++) whor[i] =  0.5 *(1.0 - Math.cos(Math.PI*(i+0.5)/stride_hor)); 
+		for (int i = 0; i < stride_vert; i++) wvert[i] = 0.5 *(1.0 - Math.cos(Math.PI*(i+0.5)/stride_vert));
+		final int stilesX = (tilesX - 1)/stride_hor;  // 9
+		final int stilesY = (tilesY - 1)/stride_vert; // 7
+		final int stiles = stilesX * stilesY;
+		int indx = 0;
+		final double[] wind = new double[stile_height*stile_width];
+		for (int iy = 0; iy < stile_height;iy++) {
+			int iy1 = (iy >= stride_vert) ? (stile_height - 1 - iy) : iy;
+			for (int ix = 0; ix < stile_width; ix++) {
+				int ix1 = (ix >= stride_hor) ? (stile_width - 1 - ix) : ix;
+				wind[indx++] = wvert[iy1] * whor[ix1];
+			}
+		}
+		final double [] stile_weight = new double [stiles];
+		final int [] stile_number = new int [stiles];
+		final Thread[] threads = ImageDtt.newThreadArray(threadsMax);
+		final AtomicInteger ai = new AtomicInteger(0);
+		final int dbg_stile =  -1;
+		for (int ithread = 0; ithread < threads.length; ithread++) {
+			threads[ithread] = new Thread() {
+				public void run() {
+					for (int sTile = ai.getAndIncrement(); sTile < stiles; sTile = ai.getAndIncrement()) {
+						if (sTile == dbg_stile) {
+							System.out.println("normalizeMotionVectorsWeights (): dbg_stile = "+dbg_stile);
+						}
+						int sTileX = sTile % stilesX;
+						int sTileY = sTile / stilesX;
+						stile_weight[sTile] = 0.0;
+						int num_tiles = 0; // for partial stiles
+						for (int iy = 0; iy < stile_height;iy++) {
+							int tileY = sTileY * stride_vert + iy;
+							if (tileY >= tilesY) continue;
+							for (int ix = 0; ix < stile_width; ix++) {
+								int tileX = sTileX * stride_hor + ix;
+								if (tileX >= tilesX) continue;
+								num_tiles++; // for partial stiles
+								int tile = tileX + tilesX*tileY;
+								if ((coord_motion[1][tile] != null) && !Double.isNaN(coord_motion[0][tile][2])) {
+									stile_weight[sTile] += coord_motion[1][tile][2];
+									stile_number[sTile]++;
+								}
+							}
+						}
+						stile_weight[sTile] *= 1.0 * stile_height * stile_height / num_tiles; // increase for partiaL stiles
+						stile_number[sTile] *= stile_height*stile_height;
+						stile_number[sTile] /= num_tiles;
+					}
+				}
+			};
+		}		      
+		ImageDtt.startAndJoin(threads);
+		/*
+		double avg_tile_str = 0.0;
+		int num_defined = 0;
+		for (int tile = 0; tile < coord_motion[0].length; tile++) {
+			if ((coord_motion[1][tile] != null) && !Double.isNaN(coord_motion[0][tile][2])) {
+				num_defined++;
+				avg_tile_str +=coord_motion[1][tile][2];
+			}
+		}
+		*/
+		double sum_stile_str = 0.0;
+		int num_stiles_defined = 0;
+		for (int sTile = 0; sTile < stiles; sTile++) {
+			if (stile_weight[sTile] > 0) {
+				num_stiles_defined++;
+				sum_stile_str+=stile_weight[sTile];
+			}
+		}
+		if (num_stiles_defined <=0) {
+			System.out.println("normalizeMotionVectorsWeights(): no defined supertiles, bailing out");
+			return;
+		}
+		final boolean [] dbg_mod = new boolean [stiles];
+		final double [] tile_new_strength = new double [tiles]; // accumulate new strengths here
+		final double avg_stile_str = sum_stile_str / num_stiles_defined;
+		final double min_combo_weight = Math.max(min_stile_weight, min_stile_fraction * avg_stile_str);
+		for (int offsy = 0; offsy < 2; offsy++) { // avoiding overlap
+			final int foffsy = offsy;
+			final int sty2 = (stilesY + 1 - offsy) / 2;
+			for (int offsx = 0; offsx < 2; offsx++) {
+				final int foffsx = offsx;
+				final int stx2 = (stilesX + 1 - offsx) / 2;
+				final int st2 = sty2*stx2;
+				ai.set(0);
+				for (int ithread = 0; ithread < threads.length; ithread++) {
+					threads[ithread] = new Thread() {
+						public void run() {
+							double [] mod_weights =  new double [stile_height*stile_width];
+							for (int indx = ai.getAndIncrement(); indx < st2; indx = ai.getAndIncrement()) {
+								int stileY = (indx/stx2)*2+foffsy;
+								int stileX = (indx%stx2)*2+foffsx;
+								int sTile = stileX + (stilesX * stileY);
+								if (sTile == dbg_stile) {
+									System.out.println("normalizeMotionVectorsWeights (): dbg_stile = "+dbg_stile);
+								}
+								boolean keep_old = false;
+								// check this tile qualifies:
+								//stile_number
+								if (stile_number[sTile] < min_stile_number) {
+									keep_old = true;
+								}
+								if (stile_weight[sTile] < min_combo_weight) {
+									keep_old = true;
+								}
+								if (stile_weight[sTile] >= avg_stile_str) { // already strong enough
+									keep_old = true;
+								}
+								Arrays.fill(mod_weights,0);
+								double sum_weights = 0.0;
+								int num_tiles = 0; // for partial stiles
+								for (int iy = 0; iy < stile_height;iy++) {
+									int tileY = stileY * stride_vert + iy;
+									if (tileY >= tilesY) continue;
+									for (int ix = 0; ix < stile_width; ix++) {
+										int tileX = stileX * stride_hor + ix;
+										if (tileX >= tilesX) continue;
+										int tile = tileX + tilesX*tileY;
+										num_tiles ++;
+										if ((coord_motion[1][tile] != null) && !Double.isNaN(coord_motion[0][tile][2])) {
+											int ltile = ix + iy * stile_width;
+											mod_weights[ltile] = keep_old ?
+													coord_motion[1][tile][2] :
+														(weight_add + coord_motion[1][tile][2] ) * weight_scale;
+											sum_weights += mod_weights[ltile];
+										}
+									}
+								}
+								sum_weights *= 1.0 * stile_height * stile_height / num_tiles  ; // increase for partial tiles 
+								if (!keep_old) {
+									if (sum_weights > avg_stile_str) { // scale back
+										double s = avg_stile_str/sum_weights;
+										for (int ltile = 0; ltile < mod_weights.length; ltile++) {
+											mod_weights[ltile] *= s;
+										}
+									}
+									for (int iy = 0; iy < stile_height;iy++) {
+										int tileY = stileY * stride_vert + iy;
+										if (tileY >= tilesY) continue;
+										for (int ix = 0; ix < stile_width; ix++) {
+											int tileX = stileX * stride_hor + ix;
+											if (tileX >= tilesX) continue;
+											int tile = tileX + tilesX*tileY;
+											int ltile = ix + iy * stile_width;
+											// remove out of range disparity
+											if (coord_motion[0][tile] != null) {
+												double disp = coord_motion[0][tile][2]; // shopuld be non-null
+												if ((disp < min_disparity) || (disp > max_disparity)) { // min/max = NaN - OK
+													mod_weights[ltile] = coord_motion[1][tile][2]; // use original value
+												}
+											}
+										}
+									}
+									dbg_mod[sTile] = true;
+								}
+								// multiply by window and accumulate
+								for (int iy = 0; iy < stile_height;iy++) {
+									int tileY = stileY * stride_vert + iy;
+									if (tileY >= tilesY) continue;
+									for (int ix = 0; ix < stile_width; ix++) {
+										int tileX = stileX * stride_hor + ix;
+										if (tileX >= tilesX) continue;
+										int tile = tileX + tilesX*tileY;
+										int ltile = ix + iy * stile_width;
+										// wind
+										tile_new_strength[tile] += wind[ltile] * mod_weights[ltile];
+									}
+								}
+							}
+						}
+					};
+				}		      
+				ImageDtt.startAndJoin(threads);
+			}
+		}
+		ai.set(0);
+		for (int ithread = 0; ithread < threads.length; ithread++) {
+			threads[ithread] = new Thread() {
+				public void run() {
+					for (int nTile = ai.getAndIncrement(); nTile < tiles; nTile = ai.getAndIncrement()) {
+						if (tile_new_strength[nTile] > 0.0) { // assuming ((coord_motion[1][nTile] != null) && !Double.isNaN(coord_motion[0][nTile][2]))
+							if (coord_motion[1][nTile] == null) {
+								System.out.println("coord_motion[1]["+nTile+"] == null, tileX="+(nTile%tilesX)+", tileY="+(nTile/tilesX));
+							} else {
+								coord_motion[1][nTile][2] = tile_new_strength[nTile]; // replace modified
+							}
+						}
+					}
+				}
+			};
+		}		      
+		ImageDtt.startAndJoin(threads);
+		return;
+	}
 	
 	public static boolean [] getMovementMask(
 			CLTParameters clt_parameters,
@@ -13123,6 +13377,7 @@ public double[][] correlateIntersceneDebug( // only uses GPU and quad
 		double[] camera_atr0 = camera_atr.clone();
 		double [][][] coord_motion = null;
 		boolean show_corr_fpn = debug_level > -1; //  -3; *********** Change to debug FPN correleation *** 
+		boolean       run_equalize = false;
 		
 		int nlma = 00;
 		for (; nlma < clt_parameters.imp.max_cycles; nlma ++) {
@@ -13148,6 +13403,75 @@ public double[][] correlateIntersceneDebug( // only uses GPU and quad
 	        	System.out.println("adjustPairsLMAInterscene() returned null");
 	        	return null;
 	        }
+			int           eq_stride_hor =         8;
+			int           eq_stride_vert =        8;
+			double        eq_min_stile_weight =   1.0;
+			int           eq_min_stile_number =   10;
+			double        eq_min_stile_fraction = 0.05;
+			double        eq_min_disparity =      5;
+			double        eq_max_disparity =    100;
+			double        eq_weight_add =         0.1;
+			double        eq_weight_scale =      10;
+			
+			if (run_equalize && near_important) {
+				TileProcessor tp = reference_QuadClt.getTileProcessor();
+				int tilesX = tp.getTilesX();
+				int tilesY = tp.getTilesY();
+				// backup coord_motion[1][][2] // strength
+				double [] strength_backup = new double [coord_motion[1].length];
+				for (int i = 0; i < strength_backup.length; i++) if (coord_motion[1][i] != null) {
+					strength_backup[i] = coord_motion[1][i][2];
+				}
+
+				while (run_equalize) {
+					// restore
+					for (int i = 0; i < strength_backup.length; i++) if (coord_motion[1][i] != null) {
+						coord_motion[1][i][2] = strength_backup[i];
+					}
+					equalizeMotionVectorsWeights(
+							coord_motion, // final double [][][] coord_motion,
+							tilesX, // final int           tilesX,
+							eq_stride_hor, // final int           stride_hor,
+							eq_stride_vert, // final int           stride_vert,
+							eq_min_stile_weight, // final double        min_stile_weight,
+							eq_min_stile_number, // final int           min_stile_number,
+							eq_min_stile_fraction, // final double        min_stile_fraction,
+							eq_min_disparity, // final double        min_disparity,
+							eq_max_disparity, // final double        max_disparity,
+							eq_weight_add, // final double        weight_add,
+							eq_weight_scale); // final double        weight_scale)
+					String [] mvTitles = {"dx", "dy","conf", "conf0", "pX", "pY","Disp","defined"}; // ,"blurX","blurY", "blur"};
+					double [][] dbg_img = new double [mvTitles.length][tilesX*tilesY];
+					for (int l = 0; l < dbg_img.length; l++) {
+						Arrays.fill(dbg_img[l], Double.NaN);
+					}
+					for (int nTile = 0; nTile <  coord_motion[0].length; nTile++) {
+						if (coord_motion[0][nTile] != null) {
+							for (int i = 0; i <3; i++) {
+								dbg_img[4+i][nTile] = coord_motion[0][nTile][i];
+							}
+						}
+						dbg_img[3] = strength_backup;
+						if (coord_motion[1][nTile] != null) {
+							for (int i = 0; i <3; i++) {
+								dbg_img[0+i][nTile] = coord_motion[1][nTile][i];
+							}
+						}
+						dbg_img[7][nTile] = ((coord_motion[0][nTile] != null)?1:0)+((coord_motion[0][nTile] != null)?2:0);
+					}
+					(new ShowDoubleFloatArrays()).showArrays( // out of boundary 15
+							dbg_img,
+							tilesX,
+							tilesY,
+							true,
+							scene_QuadClt.getImageName()+"-"+reference_QuadClt.getImageName()+"-coord_motion-eq",
+							mvTitles);
+				}
+
+			}
+
+	        
+	        
 			intersceneLma.prepareLMA(
 					camera_xyz0,         // final double []   scene_xyz0,     // camera center in world coordinates (or null to use instance)
 					camera_atr0,         // final double []   scene_atr0,     // camera orientation relative to world frame (or null to use instance)
@@ -13255,6 +13579,12 @@ public double[][] correlateIntersceneDebug( // only uses GPU and quad
 					true,
 					scene_QuadClt.getImageName()+"-"+reference_QuadClt.getImageName()+"-CORR-FPN",
 					fpn_dbg_titles);
+			
+			
+			
+			
+			
+			
 		}
 		
 		

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

@@ -45,6 +45,7 @@ import com.elphel.imagej.cameras.CLTParameters;
 import com.elphel.imagej.cameras.ColorProcParameters;
 import com.elphel.imagej.cameras.EyesisCorrectionParameters;
 import com.elphel.imagej.common.DoubleGaussianBlur;
+import com.elphel.imagej.common.PolynomialApproximation;
 import com.elphel.imagej.common.ShowDoubleFloatArrays;
 import com.elphel.imagej.correction.CorrectionColorProc;
 import com.elphel.imagej.correction.EyesisCorrections;
@@ -2426,20 +2427,234 @@ public class QuadCLT extends QuadCLTCPU {
 			}
 			double [] offsets_old = ref_scene.getLwirOffsets();
 			double [] scales_old =  ref_scene.getLwirScales();
+			double [] scales2_old = ref_scene.getLwirScales2();
 			double [] offsets_new = new double[num_sens];
 			double [] scales_new =  new double[num_sens];
+			double [] scales2_new = new double[num_sens];
 			for (int n = 0; n < num_sens; n++) {
 				scales_new[n] = scales_old[n]/scales[n];
-//				offsets_new[n] = offsets_old[n]  - offsets[n] / scales_old[n];
 				offsets_new[n] = offsets_old[n]  + offsets[n] / scales_old[n];
 			}
 			System.out.println("calibratePhotometric() Updated calibration:");
 			for (int n = 0; n < num_sens; n++) {
-				System.out.println(String.format("%2d: %8.4f %8.6f", n,offsets_new[n],scales_new[n]));
+				System.out.println(String.format("%2d: %8.4f %8.6f %8.6f", n,offsets_new[n],scales_new[n], scales2_new[n]));
 			}
 			System.out.println();
-			ref_scene.setLwirOffsets(offsets_new);
+			ref_scene.setLwirOffsets (offsets_new);
 			ref_scene.setLwirScales  (scales_new);
+			ref_scene.setLwirScales2 (scales2_new);
+			return true;
+	}
+	public static boolean calibratePhotometric2( // quadratic
+			CLTParameters     clt_parameters,
+			final QuadCLT     ref_scene, // now - may be null - for testing if scene is rotated ref
+			final double      min_strength,
+			final double      max_diff,    // 30.0
+			final int         num_refines, // 2
+			final double [][] combo_dsn_final,     // double [][]    combo_dsn_final, // dls,			
+			int               threadsMax,
+			final boolean     debug)
+	{
+		
+		// filter disparity by LMA only, same bg and fg
+			double []         disparity_ref= combo_dsn_final[OpticalFlow.COMBO_DSN_INDX_DISP].clone();
+			for (int i = 00; i < disparity_ref.length; i++) {
+				if (combo_dsn_final[OpticalFlow.COMBO_DSN_INDX_STRENGTH][i] < min_strength) {
+					disparity_ref[i] = Double.NaN;
+				} else if (combo_dsn_final[OpticalFlow.COMBO_DSN_INDX_DISP_BG_ALL][i] !=
+						combo_dsn_final[OpticalFlow.COMBO_DSN_INDX_DISP][i]) {
+					disparity_ref[i] = Double.NaN;
+				}
+			}
+			
+			ImagePlus img_ref = renderGPUFromDSI(
+					-1,                 // final int         sensor_mask,
+					false,              // final boolean     merge_channels,
+					null,               // final Rectangle   full_woi_in,      // show larger than sensor WOI in tiles (or null)
+					 clt_parameters,    // CLTParameters     clt_parameters,
+					 disparity_ref,     // double []         disparity_ref,
+					 OpticalFlow.ZERO3, // final double []   scene_xyz, // camera center in world coordinates
+					 OpticalFlow.ZERO3, // final double []   scene_atr, // camera orientation relative to world frame
+					 ref_scene,         // final QuadCLT     scene,
+					 ref_scene,         // final QuadCLT     ref_scene, // now - may be null - for testing if scene is rotated ref
+					false,              // final boolean     toRGB,
+					true,               // final boolean     show_nan,
+					"PHOTOMETRIC",      // String            suffix,
+					threadsMax,         // int               threadsMax,
+					-2);        // final int         debugLevel);
+			img_ref.show();
+			ImageStack imageStack = img_ref.getStack();
+			int num_sens=imageStack.getSize();
+			float [] fpixels;
+			double [][] dpixels = new double [num_sens][];
+			for (int n = 0; n < num_sens; n++) {
+				fpixels = (float[]) imageStack.getPixels(n + 1);
+				dpixels[n] = new double [fpixels.length];
+				for (int i = 0; i < fpixels.length; i++) {
+					dpixels[n][i] = fpixels[i];
+				}
+			}
+			boolean quadratic = true;
+//			double [][] dpix_orig = new double[num_sens][];
+//			for (int n = 0; n < num_sens; n++) {
+//				dpix_orig[n] =  dpixels[n].clone();
+//			}					
+			int width =  img_ref.getWidth();
+			int height = img_ref.getHeight();
+			int len = width* height;
+			double [] avg_pix = new double [len];
+			double [] offsets = new double[dpixels.length];
+			double [] scales = new double[dpixels.length];
+			double    s0 =  0.0;
+			double    sx=   0.0;
+			double    sx2 = 0.0;
+			double [] sy = new double[num_sens];
+			double [] sxy = new double[num_sens];
+			boolean [] good_pix = new boolean[len];
+			double [][] pa_coeff = new double[num_sens][];
+			Arrays.fill(good_pix, true);
+			for (int nref = 0; nref < num_refines; nref++) {
+				Arrays.fill(avg_pix, 0.0);
+				int num_good = 0;
+				for (int i = 0; i < len; i++) {
+					for (int n = 0; n < num_sens; n++) {
+						avg_pix[i]+=dpixels[n][i];
+						good_pix[i] &= !Double.isNaN(dpixels[n][i]);
+					}
+					avg_pix[i] /= dpixels.length;
+					if (good_pix[i]) {
+						num_good++;
+					}
+				}
+				double [][] pa_data = new double [num_good][2];  
+				for (int nsens = 0; nsens < num_sens; nsens++) {
+					int indx = 0;
+					for (int i = 0; i < len; i++) if (good_pix[i]){
+						pa_data[indx][0] = dpixels[nsens][i]; 
+						pa_data[indx][1] = avg_pix[i];
+						indx++;
+					}
+					// quadratic
+					pa_coeff[nsens] =(new PolynomialApproximation(0)).polynomialApproximation1d(pa_data, quadratic ? 2 : 1);
+				}
+				if (debug) {
+					System.out.println("calibratePhotometric() nref="+nref);
+					for (int n = 0; n < num_sens; n++) {
+						System.out.println(String.format("%2d: %8.4f %8.6f %8.6f", n,pa_coeff[n][0],pa_coeff[n][1], 1e6*pa_coeff[n][2]));
+					}
+					System.out.println();
+					double [][] diffs = new double [num_sens][len];
+					for (int n = 0; n < num_sens; n++) {
+						Arrays.fill(diffs[n], Double.NaN);
+					}
+					for (int i = 0; i < len; i++)  if (good_pix[i]) { // if (!Double.isNaN(avg_pix[i])){
+						for (int n = 0; n < num_sens; n++) {
+//							diffs[n][i] = dpixels[n][i] - (scales[n] * avg_pix[i] + offsets[n]);
+							diffs[n][i] = -(avg_pix[i] - pa_coeff[n][2] * dpixels[n][i] * dpixels[n][i] - pa_coeff[n][1] * dpixels[n][i] - pa_coeff[n][0]); 
+						}
+					}
+					(new ShowDoubleFloatArrays()).showArrays( // out of boundary 15
+							diffs,
+							width,
+							height,
+							true,
+							"photometric-quad-err"+nref);
+				}				
+				
+				
+				Arrays.fill(offsets,0.0);
+				Arrays.fill(scales,0.0);
+				s0 =  0.0;
+				sx=   0.0;
+				sx2 = 0.0;
+				Arrays.fill(sy,0.0);
+				Arrays.fill(sxy,0.0);
+				for (int i = 0; i < len; i++)  if (good_pix[i]) { // !Double.isNaN(avg_pix[i])){
+					s0 +=  1.0;
+					sx +=  avg_pix[i];
+					sx2 += avg_pix[i] * avg_pix[i];
+					for (int n = 0; n < num_sens; n++) {
+						sy[n] +=  dpixels[n][i];
+						sxy[n] += dpixels[n][i] * avg_pix[i];
+					}
+				}
+				for (int n = 0; n < num_sens; n++) {
+					double d = s0 * sx2 + sx * sy[n];
+					scales[n] =  (sxy[n] * s0 + sy[n] * sy[n]) / d;
+					offsets[n] = (sy[n] * sx2 -sxy[n]*sx) / d;
+				}
+				
+				
+				if (debug) {
+					System.out.println("calibratePhotometric() nref="+nref);
+					for (int n = 0; n < num_sens; n++) {
+						System.out.println(String.format("%2d: %8.4f %8.6f", n,offsets[n],scales[n]));
+					}
+					System.out.println();
+					double [][] diffs = new double [num_sens][len];
+					for (int n = 0; n < num_sens; n++) {
+						Arrays.fill(diffs[n], Double.NaN);
+					}
+					for (int i = 0; i < len; i++)  if (good_pix[i]) { // if (!Double.isNaN(avg_pix[i])){
+						for (int n = 0; n < num_sens; n++) {
+							diffs[n][i] = dpixels[n][i] - (scales[n] * avg_pix[i] + offsets[n]);
+						}
+					}
+					(new ShowDoubleFloatArrays()).showArrays( // out of boundary 15
+							diffs,
+							width,
+							height,
+							true,
+							"photometric-err"+nref);
+
+					// dpix_orig[n] =  dpixels[n].clone();
+					double [][] corrected = new double [num_sens][len];
+					for (int n = 0; n < num_sens; n++) {
+						Arrays.fill(corrected[n], Double.NaN);
+						for (int i = 0; i < len; i++)  if (!Double.isNaN(dpixels[n][i])){
+							corrected[n][i] = (dpixels[n][i] - offsets[n])/scales[n];
+						}
+					}
+
+					(new ShowDoubleFloatArrays()).showArrays( // out of boundary 15
+							corrected,
+							width,
+							height,
+							true,
+							"photometric-corr"+nref);
+				}
+				//max_diff
+				if (nref < (num_refines-1)) {
+					for (int i = 0; i < len; i++)  if (good_pix[i]) { // !Double.isNaN(avg_pix[i])){
+						for (int n = 0; n < num_sens; n++) {
+							double diff = Math.abs(dpixels[n][i] - (scales[n] * avg_pix[i] + offsets[n]));
+							if (diff > max_diff) {
+//								dpixels[n][i] = Double.NaN;
+								good_pix[i] = false;
+								break;
+							}
+						}
+					}
+				}
+			}
+			double [] offsets_old = ref_scene.getLwirOffsets();
+			double [] scales_old =  ref_scene.getLwirScales();
+			double [] scales2_old = ref_scene.getLwirScales2();
+			double [] offsets_new = new double[num_sens];
+			double [] scales_new =  new double[num_sens];
+			double [] scales2_new = new double[num_sens];
+			for (int n = 0; n < num_sens; n++) {
+				scales_new[n] = scales_old[n]/scales[n];
+				offsets_new[n] = offsets_old[n]  + offsets[n] / scales_old[n];
+			}
+			System.out.println("calibratePhotometric() Updated calibration:");
+			for (int n = 0; n < num_sens; n++) {
+				System.out.println(String.format("%2d: %8.4f %8.6f %8.6f", n,offsets_new[n],scales_new[n], scales2_new[n]));
+			}
+			System.out.println();
+			ref_scene.setLwirOffsets (offsets_new);
+			ref_scene.setLwirScales  (scales_new);
+			ref_scene.setLwirScales2 (scales2_new);
 			return true;
 	}
 	

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

@@ -153,7 +153,7 @@ public class QuadCLTCPU {
     boolean                                                is_aux = false;
     double  []                                             lwir_offsets = null; // per image subtracted values
     double []                                              lwir_scales =  null;  // per image scales
-//double []                                              lwir_scales2 = null;  // per image quadratic scales
+    double []                                              lwir_scales2 = null;  // per image quadratic 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)
@@ -266,6 +266,7 @@ public class QuadCLTCPU {
 //		this.is_mono =                 qParent.is_mono;
 		this.lwir_offsets =           ErsCorrection.clone1d(qParent.lwir_offsets);
 		this.lwir_scales =            ErsCorrection.clone1d(qParent.lwir_scales);
+		this.lwir_scales2 =           ErsCorrection.clone1d(qParent.lwir_scales2);
 		this.lwir_offset =            qParent.lwir_offset;
 		this.lwir_cold_hot =          ErsCorrection.clone1d(qParent.lwir_cold_hot);
 		this.ds_from_main =           ErsCorrection.clone2d(qParent.ds_from_main);
@@ -1798,6 +1799,12 @@ public class QuadCLTCPU {
     	}
     	return lwir_scales;
     }
+    public double [] getLwirScales2() {
+    	if (lwir_scales2 == null) {
+    		lwir_scales2 = new double [getNumSensors()]; // all 0;
+    	}
+    	return lwir_scales2;
+    }
     public void setLwirOffsets(double [] offsets) {
     	lwir_offsets = offsets; // will need to update properties!
     	if (offsets != null) {
@@ -1812,6 +1819,9 @@ public class QuadCLTCPU {
     public void setLwirScales(double [] scales) {
     	lwir_scales = scales; // will need to update properties!
     }
+    public void setLwirScales2(double [] scales2) {
+    	lwir_scales2 = scales2; // will need to update properties!
+    }
 
     public double [] getColdHot() { // USED in lwir
     	return lwir_cold_hot;
@@ -2011,6 +2021,10 @@ public class QuadCLTCPU {
 			properties.setProperty(prefix+"lwir_scales",
 					IntersceneMatchParameters.doublesToString(this.lwir_scales));
 		}
+		if (this.lwir_scales2 != null) {
+			properties.setProperty(prefix+"lwir_scales2",
+					IntersceneMatchParameters.doublesToString(this.lwir_scales2));
+		}
 
 	}
 
@@ -2053,6 +2067,12 @@ public class QuadCLTCPU {
 			properties.setProperty(this_prefix+"lwir_scales",
 					IntersceneMatchParameters.doublesToString(this.lwir_scales));
 		}
+		if (other_properties.getProperty(other_prefix+"lwir_scales2")!=null) {
+			this.lwir_scales2= IntersceneMatchParameters.StringToDoubles(
+					other_properties.getProperty(other_prefix+"lwir_scales2"),0);
+			properties.setProperty(this_prefix+"lwir_scales2",
+					IntersceneMatchParameters.doublesToString(this.lwir_scales2));
+		}
 		
 		
 		/*
@@ -2173,6 +2193,10 @@ public class QuadCLTCPU {
 			this.lwir_scales= IntersceneMatchParameters.StringToDoubles(
 					properties.getProperty(prefix+"lwir_scales"),0);
 		}
+		if (properties.getProperty(prefix+"lwir_scales2")!=null) {
+			this.lwir_scales2= IntersceneMatchParameters.StringToDoubles(
+					properties.getProperty(prefix+"lwir_scales2"),0);
+		}
 			
 		
 		
@@ -5572,11 +5596,13 @@ public class QuadCLTCPU {
 			  }
 			  double [] offsets =  getLwirOffsets();
 			  double [] scales =   getLwirScales();
+			  double [] scales2 =  getLwirScales2();
 			  channelLwirApplyEqualize( // now apply (was part of channelLwirEqualize() )
 					  channelFiles, // int [] channelFiles,
 					  imp_srcs,     // ImagePlus [] imp_srcs,
 					  offsets,      // double []    offsets,
 					  scales,       // double []    scales,
+					  scales2,      // double []    scales2,
 					  threadsMax,
 					  debugLevel);
 		  }
@@ -6162,6 +6188,7 @@ public class QuadCLTCPU {
 			  ImagePlus [] imp_srcs,
 			  double []    offsets,
 			  double []    scales,
+			  double []    scales2,
 			  final int    threadsMax,
 			  int          debugLevel){
 		  final Thread[] threads = ImageDtt.newThreadArray(threadsMax);
@@ -6174,9 +6201,11 @@ public class QuadCLTCPU {
 						  if (nFile >=0) {
 							  float fd = (float)offsets[srcChannel];
 							  float fscale =  (float) scales[srcChannel];
+							  float fscale2 = (float) scales2[srcChannel];
 							  float [] pixels = (float []) imp_srcs[srcChannel].getProcessor().getPixels();
 							  for (int i = 0; i < pixels.length; i++) {
-								  pixels[i] = (pixels[i] - fd) * fscale;
+								  float poffs = (pixels[i] - fd);
+								  pixels[i] = poffs * fscale + poffs*poffs*fscale2;
 							  }
 						  }
 					  }