Commit c8d3d7d9 authored by Andrey Filippov's avatar Andrey Filippov
Browse files

cleaning up with eigen, adding debug image generation for LMA

parent 477343de

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

+2 −2
Original line number Diff line number Diff line
@@ -2351,6 +2351,7 @@ public class ImageDtt extends ImageDttCPU { 
			final boolean             td_nopd_only,    // only use TD accumulated data if no safe PD is available for the tile.

			final boolean             eig_use_neibs,   // use correlation from 9 tiles with neibs, if single-tile fails

			final boolean             eig_remove_neibs, //remove weak (by-neibs) tiles if they have strong (by-single) neighbor

			final boolean             eig_filt_other,   // apply other before-eigen filters

//			final double              min_str_nofpn,    //  = 0.25;

			final double              eig_str_sum_nofpn,// = 0.8; // 5;

			final double              eig_str_neib_nofpn,// = 0.8; // 5;
@@ -2720,10 +2721,9 @@ public class ImageDtt extends ImageDttCPU { 
			}

			startAndJoin(threads);

		}

		boolean old_filter = false; // true;

		// Reduce weight if differs much from average of 8 neighbors, large disparity, remove too few neibs

		final double scale_num_neib = ((weight_zero_neibs >= 0) && (weight_zero_neibs < 1.0)) ? (weight_zero_neibs * 8/(1.0 - weight_zero_neibs)): 0.0;

		if (old_filter) {

		if (eig_filt_other) {

			ai.set(0);

			for (int ithread = 0; ithread < threads.length; ithread++) {

				threads[ithread] = new Thread() {

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

+172 −33

File changed.

Preview size limit exceeded, changes collapsed.

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

+143 −19
Original line number Diff line number Diff line
@@ -31,11 +31,14 @@ import java.io.Serializable; 
import java.security.MessageDigest;

import java.security.NoSuchAlgorithmException;

import java.util.ArrayList;

import java.util.Arrays;

import java.util.concurrent.atomic.AtomicInteger;

import java.util.concurrent.atomic.DoubleAdder;



import javax.xml.bind.DatatypeConverter;



import com.elphel.imagej.common.ShowDoubleFloatArrays;



import Jama.Matrix;



public class IntersceneLma {
@@ -44,6 +47,7 @@ public class IntersceneLma { 
	private double []         good_or_bad_rms = null; // just for diagnostics, to read last (failed) rms

	private double []         initial_rms =     null; // {rms, rms_pure}, first-calcualted rms

	private double []         y_vector =        null; // sum of fx(initial parameters) and correlation offsets

	private double []         s_vector =        null; // strength component - just for debug images

	private double []         last_ymfx =       null;

	private double [][]       last_jt =         null;

	private double []         weights; // normalized so sum is 1.0 for all - samples and extra regularization
@@ -69,6 +73,8 @@ public class IntersceneLma { 
	private double            disparity_weight = 1.0;  // relative weight of disparity errors

	

	private double [][][]     eig_trans = null;

	private int               tilesX = -1;

	private String            dbg_prefix = null;

	

	public IntersceneLma(

			boolean thread_invariant,
@@ -259,20 +265,21 @@ public class IntersceneLma { 
			// now includes optional Disparity as the last element (for num_components==3) 

			final double [][] vector_XYSDS,// optical flow X,Y, confidence obtained from the correlate2DIterate()

			final double [][] centers,     // macrotile centers (in pixels and average disparities

			

			boolean           first_run,

			String            dbg_prefix, // null or image name prefix

			final int         debug_level) {

		// befolre getFxDerivs

		// before getFxDerivs

		eig_trans = setEigenTransform(

				eig_max_sqrt, // final double   eig_max_sqrt,

				eig_min_sqrt, // final double   eig_min_sqrt,

				eigen); // final double [][] eigen); // [tilesX*tilesY]{lamb0_x,lamb0_y, lamb0, lamb1} eigenvector0[x,y],lam0,lam1



		

		tilesX = reference_QuadClt.getTileProcessor().getTilesX(); // just for debug images

		scenesCLT = new QuadCLT [] {reference_QuadClt, scene_QuadClt};

		par_mask = param_select;

		macrotile_centers = centers;

		num_samples = num_components * centers.length;

		this.dbg_prefix = dbg_prefix;

		s_vector = (this.dbg_prefix != null) ? (new double[centers.length]): null; // original strength

		ErsCorrection ers_ref =   reference_QuadClt.getErsCorrection();

		ErsCorrection ers_scene = scene_QuadClt.getErsCorrection();

		final double []   scene_xyz = (scene_xyzatr0 != null) ? scene_xyzatr0[0] : ers_scene.camera_xyz;
@@ -366,8 +373,11 @@ public class IntersceneLma { 
			if (vector_XYSDS[i] != null){

				y_vector[num_components * i + 0] += vector_XYSDS[i][0];

				y_vector[num_components * i + 1] += vector_XYSDS[i][1];

				if (num_components > 2) {

					y_vector[num_components * i + 2] += vector_XYSDS[i][2];

				if (num_components > 2) { // ****************************** [3] - disparity? Was BUG: [2] !

					y_vector[num_components * i + 2] += vector_XYSDS[i][3]; // vector_XYSDS[i][2]; 

				}

				if (s_vector != null) {

					s_vector[i] = vector_XYSDS[i][2];

				}

			}

		}
@@ -562,6 +572,10 @@ public class IntersceneLma { 
					break; // not used in lwir

				}

			}

			if (dbg_prefix != null) {

				showDebugImage(dbg_prefix+"-"+iter+(rslt[0]?"-GOOD":"-BAD"));

			}



		}

		if (rslt[0]) { // better

			if (iter >= num_iter) { // better, but num tries exceeded
@@ -578,6 +592,9 @@ public class IntersceneLma { 
				if (debug_level > 1) System.out.println("Step "+iter+": Failed to converge");

			}

		}

		if (dbg_prefix != null) {

			showDebugImage(dbg_prefix+"-FINAL");

		}

		boolean show_intermediate = true;

		if (show_intermediate && (debug_level > 0)) {

			System.out.println("LMA: full RMS="+last_rms[0]+" ("+initial_rms[0]+"), pure RMS="+last_rms[1]+" ("+initial_rms[1]+") + lambda="+lambda);
@@ -657,6 +674,9 @@ public class IntersceneLma { 
							true); //boolean graphic)

				*/

			}

			if (dbg_prefix != null) {

				showDebugImage(dbg_prefix+"-INIT");

			}

		}

		Matrix y_minus_fx_weighted = new Matrix(this.last_ymfx, this.last_ymfx.length);
@@ -829,7 +849,7 @@ public class IntersceneLma { 
							w = vector_XYSDS[iMTile][4] * disparity_weight;

							if (Double.isNaN(w) || Double.isNaN(vector_XYSDS[iMTile][3])) {

								w = 0;

								vector_XYSDS[iMTile][3] = 0.0;

								vector_XYSDS[iMTile][3] = 0.0; // disparity

							}

							weights[num_components * iMTile + 2] = w;

							sw_arr[thread_num] += w;
@@ -1036,6 +1056,7 @@ public class IntersceneLma { 
									}

								}

							} else if (mb_mode) {

								if (jt != null) {

									for (int i = 0; i < par_indices.length; i++) {

										jt[i][2 * iMTile + 0] = Double.NaN; // pX

										jt[i][2 * iMTile + 1] = Double.NaN; // ; // pY (disparity is not used)
@@ -1044,6 +1065,7 @@ public class IntersceneLma { 
							}

						}

					}

				}

			};

		}		      

		ImageDtt.startAndJoin(threads);
@@ -1053,8 +1075,10 @@ public class IntersceneLma { 
		// pull to the initial parameter values

		for (int i = 0; i < par_indices.length; i++) {

			fx [i + num_samples] = vector[i]; // - parameters_initial[i]; // scale will be combined with weights

			if (jt != null) { 

				jt[i][i + num_samples] = 1.0; // scale will be combined with weights

			}

		}

///		if (parameters_pull != null){

///			for (int i = 0; i < par_indices.length; i++) {

///				fx [i + num_samples] -= parameters_pull[i]; // - parameters_initial[i]; // scale will be combined with weights
@@ -1102,15 +1126,113 @@ public class IntersceneLma { 
		return wjtjl;

	}

    	

	public void showDebugImage(	String dbg_title ) { // includes "_iteration_number" or "_final"

		if (s_vector == null) {

			return;

		}

		int tilesY =  s_vector.length / tilesX;

		String [] titles = {"dx","dy","dr", "str", "e0", "e1", "e"};

		double [][] dbg_img = new double [titles.length][tilesX*tilesY];

		for (int l = 0; l < dbg_img.length; l++) {

			Arrays.fill(dbg_img[l], Double.NaN);

		}

		double [] fx = getFxDerivs(

				parameters_vector, // double []         vector,

				null,              // final double [][] jt, // should be null or initialized with [vector.length][]

				scenesCLT[1],      // final QuadCLT     scene_QuadClt,

				scenesCLT[0],      // final QuadCLT     reference_QuadClt,

				-1);               // debug_level);      // final int         debug_level)

		double [] ymfxw_m = getYminusFxWeighted(

				fx,     // final double []   fx,

				null,   // final double []   rms_fp // null or [2]

				true);//final boolean     force_metric // when true, ignore transform with eig_trans and use linear dx, dy in pixels

        double [] ymfxw_e = null;

        if (eig_trans != null) {

        	ymfxw_e = getYminusFxWeighted(

    				fx,     // final double []   fx,

    				null,   // final double []   rms_fp // null or [2]

    				false);//final boolean     force_metric // when true, ignore transform with eig_trans and use linear dx, dy in pixels

        }

        for (int nTile = 0; nTile < s_vector.length; nTile++) {

        	int indx = num_components * nTile;

        	double w = weights[num_components * nTile];

        	if ((weights[indx] > 0) && (weights[indx+1] > 0)) {

        		for (int i = 0; i < 2; i++) {

        			ymfxw_m[indx + i] /= weights[indx + i];

        			if (ymfxw_e != null) {

            			ymfxw_e[indx + i] /= weights[indx + i];

        			}

        		}

        		double dx = ymfxw_m[indx + 0]; 

        		double dy = ymfxw_m[indx + 1]; 

        		dbg_img[0][nTile] = dx; 

        		dbg_img[1][nTile] = dy;

        		dbg_img[2][nTile] = Math.sqrt(dx*dx+dy*dy);

        		dbg_img[3][nTile] = s_vector[nTile];

    			if (ymfxw_e != null) {

            		double d0 = ymfxw_e[indx + 0]; 

            		double d1 = ymfxw_e[indx + 1]; 

            		dbg_img[4][nTile] = d0; 

            		dbg_img[5][nTile] = d1;

            		dbg_img[6][nTile] = Math.sqrt(d0*d0+d1*d1);

    			}        	

        	}

        }

		if (ymfxw_e == null) {

			dbg_img[4]=null;

			dbg_img[5]=null;

			dbg_img[6]=null;

		}

		ShowDoubleFloatArrays.showArrays( // out of boundary 15

				dbg_img,

				tilesX,

				tilesY,

				true,

				dbg_title,

				titles);

		

	}

	

	

	public boolean isEigenNormalized() {

		return (eig_trans != null);

	}

	public double [] calcRMS (boolean metric) {

		double [] rms_fp = new double [2];

		double [] fx = getFxDerivs(

				parameters_vector, // double []         vector,

				null,              // final double [][] jt, // should be null or initialized with [vector.length][]

				scenesCLT[1],      // final QuadCLT     scene_QuadClt,

				scenesCLT[0],      // final QuadCLT     reference_QuadClt,

				-1);               // debug_level);      // final int         debug_level)

		last_ymfx = getYminusFxWeighted(

				fx,     // final double []   fx,

				rms_fp, // final double []   rms_fp // null or [2]

				metric);//final boolean     force_metric // when true, ignore transform with eig_trans and use linear dx, dy in pixels

		return rms_fp;

	}

	

	

	private double [] getYminusFxWeighted(

			final double []   fx,

			final double []   rms_fp // null or [2]

			final double []   rms_fp) { // null or [2]

		return getYminusFxWeighted(

				fx,     // final double []   fx,

				rms_fp, // final double []   rms_fp, // null or [2]

				false); // final boolean     force_metric)			



	}



	private double [] getYminusFxWeighted(

			final double []   fx,

			final double []   rms_fp, // null or [2]

			final boolean     force_metric // when true, ignore transform with eig_trans and use linear dx, dy in pixels			

			) {

		double [] ymfxw;

		if (thread_invariant) {

			 ymfxw =  getYminusFxWeightedInvariant(fx,rms_fp); // null or [2]

			 ymfxw =  getYminusFxWeightedInvariant(fx,rms_fp, force_metric); // null or [2]

		} else {

			 ymfxw = getYminusFxWeightedFast     (fx,rms_fp); // null or [2]

			 ymfxw = getYminusFxWeightedFast      (fx,rms_fp, force_metric); // null or [2]

		}

		return ymfxw; 

	}
@@ -1120,14 +1242,15 @@ public class IntersceneLma { 
	

	private double [] getYminusFxWeightedInvariant(

			final double []   fx,

			final double []   rms_fp // null or [2]

			final double []   rms_fp, // null or [2]

			final boolean     force_metric // when true, ignore transform with eig_trans and use linear dx, dy in pixels			

			) {

		final Thread[]      threads =     ImageDtt.newThreadArray(QuadCLT.THREADS_MAX);

		final AtomicInteger ai =          new AtomicInteger(0);

		final double []     wymfw =       new double [fx.length];

		double s_rms; 

		final double [] l2_arr = new double [num_samples]; 

		if (eig_trans != null) {

		if (!force_metric && (eig_trans != null)) {

			for (int ithread = 0; ithread < threads.length; ithread++) {

				threads[ithread] = new Thread() {

					public void run() {
@@ -1213,14 +1336,15 @@ public class IntersceneLma { 
	

	private double [] getYminusFxWeightedFast( // problems. at least with eigen?

			final double []   fx,

			final double []   rms_fp // null or [2]

			final double []   rms_fp, // null or [2]

			final boolean     force_metric // when true, ignore transform with eig_trans and use linear dx, dy in pixels			

			) {

		final Thread[]      threads =     ImageDtt.newThreadArray(QuadCLT.THREADS_MAX);

		final AtomicInteger ai =          new AtomicInteger(0);

		final double []     wymfw =       new double [fx.length];

		final AtomicInteger ati = new AtomicInteger(0);

		final double [] l2_arr = new double [threads.length];

		if (eig_trans != null) {

		if (!force_metric && (eig_trans != null)) {

			for (int ithread = 0; ithread < threads.length; ithread++) {

				threads[ithread] = new Thread() {

					public void run() {

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

+16 −4
Original line number Diff line number Diff line
@@ -451,6 +451,9 @@ public class IntersceneMatchParameters { 
	public double   eig_min_sqrt =                  1.0;   // for sqrt(lambda) - limit minimal sqrt(lambda) - can be sharp for very small max

	public boolean  eig_use_neibs =                true;   // use correlation from 9 tiles with neibs, if single-tile fails

	public boolean  eig_remove_neibs =             true;   // remove weak (by-neibs) tiles if they have strong (by-single) neighbor

	public boolean  eig_filt_other =              false;   // apply other before-eigen filters

	public  double  eig_max_rms =                   2.0;   // eigen-normalized maximal RMS to consider adjustment to be a failure

	



/*

min_str_sum_nofpn	0.22
@@ -1358,8 +1361,10 @@ min_str_neib_fpn 0.35 
				"Use correlation from 9 tiles with neibs, if single-tile fails");

		gd.addCheckbox ("Remove weak by strong neighbors",           this.eig_remove_neibs,

				"Remove weak (by-neibs) tiles if they have strong (by-single) neighbor");

		

		

		gd.addCheckbox ("Apply other filters",                       this.eig_filt_other,

				"Apply other (before-eigen) filters");

		gd.addNumericField("Maximal eigen-normalized RMS",           this.eig_max_rms, 5,7,"",

				"Maximal eigen-normalized RMSE for LMA adjustment. Replaces \"Maximal RMS to fail\" setting below.");

		

		

		gd.addMessage  ("Filtering tiles for interscene matching");
@@ -2048,7 +2053,8 @@ min_str_neib_fpn 0.35 
		this.eig_min_sqrt =             gd.getNextNumber();

		this.eig_use_neibs =            gd.getNextBoolean();

		this.eig_remove_neibs =         gd.getNextBoolean();

		

		this.eig_filt_other =           gd.getNextBoolean();

		this.eig_max_rms =              gd.getNextNumber();

		this.use_combo_reliable =       gd.getNextBoolean();

		this.ref_need_lma =             gd.getNextBoolean();

		this.ref_need_lma_combo =       gd.getNextBoolean();
@@ -2603,6 +2609,8 @@ min_str_neib_fpn 0.35 
		properties.setProperty(prefix+"eig_min_sqrt",         this.eig_min_sqrt+"");        // double

		properties.setProperty(prefix+"eig_use_neibs",        this.eig_use_neibs+"");       // boolean

		properties.setProperty(prefix+"eig_remove_neibs",     this.eig_remove_neibs+"");    // boolean

		properties.setProperty(prefix+"eig_filt_other",       this.eig_filt_other+"");      // boolean

		properties.setProperty(prefix+"eig_max_rms",          this.eig_max_rms+"");         // double

		

		properties.setProperty(prefix+"use_combo_reliable",   this.use_combo_reliable+"");  // boolean

		properties.setProperty(prefix+"ref_need_lma",         this.ref_need_lma+"");        // boolean
@@ -3120,6 +3128,8 @@ min_str_neib_fpn 0.35 
		if (properties.getProperty(prefix+"eig_min_sqrt")!=null)         this.eig_min_sqrt=Double.parseDouble(properties.getProperty(prefix+"eig_min_sqrt"));

        if (properties.getProperty(prefix+"eig_use_neibs")!=null)        this.eig_use_neibs=Boolean.parseBoolean(properties.getProperty(prefix+"eig_use_neibs"));		

		if (properties.getProperty(prefix+"eig_remove_neibs")!=null)     this.eig_remove_neibs=Boolean.parseBoolean(properties.getProperty(prefix+"eig_remove_neibs"));		

		if (properties.getProperty(prefix+"eig_filt_other")!=null)       this.eig_filt_other=Boolean.parseBoolean(properties.getProperty(prefix+"eig_filt_other"));		

		if (properties.getProperty(prefix+"eig_max_rms")!=null)          this.eig_max_rms=Double.parseDouble(properties.getProperty(prefix+"eig_max_rms"));



		if (properties.getProperty(prefix+"use_combo_reliable")!=null)   this.use_combo_reliable=Boolean.parseBoolean(properties.getProperty(prefix+"use_combo_reliable"));

		else if (properties.getProperty(prefix+"use_combo_relaible")!=null)   this.use_combo_reliable=Boolean.parseBoolean(properties.getProperty(prefix+"use_combo_relaible"));
@@ -3645,6 +3655,8 @@ min_str_neib_fpn 0.35 
		imp.eig_min_sqrt =          this.eig_min_sqrt;

		imp.eig_use_neibs =         this.eig_use_neibs;

		imp.eig_remove_neibs =      this.eig_remove_neibs;

		imp.eig_filt_other =        this.eig_filt_other;

		imp.eig_max_rms =           this.eig_max_rms;

		

		imp.use_combo_reliable    = this.use_combo_reliable;

		imp.ref_need_lma          = this.ref_need_lma;