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

improving matching from the center to mitigate scale/rotations

parent 5e03168f

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

+1 −1
Original line number Diff line number Diff line
@@ -1299,7 +1299,7 @@ adjusted affines[1] for a pair: 1694564291_293695/1694564778_589341 
				int nSlices = stack_scenes.getSize();

				String [] scene_names = new String [nSlices];

				for (int n = 0; n < scene_names.length; n++) {

					scene_names[n] = stack_scenes.getSliceLabel(n+1).substring(0,17);

					scene_names[n] = stack_scenes.getSliceLabel(n+1).substring(0,17); // wrong image opened-> StringIndexOutOfBoundsException: begin 0, end 17, length 8

				}

				if (all_scenes != null) {

					all_scenes[0] = scene_names;

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

+25 −3
Original line number Diff line number Diff line
@@ -1013,13 +1013,17 @@ public class OrthoMapsCollection implements Serializable{ 
				null, // double [][] jtj,

				Double.NaN, // double rms,

				zoom_lev); // int zoom_lev); 

		double           max_std = 1.5;      // maximal standard deviation to limit center area  

		double           max_std =     0.5;      // maximal standard deviation to limit center area  

		double           min_std_rad = 2.0;  // minimal radius of the central area (if less - fail)



		while ((Math.abs(nx) <= nabs) && (Math.abs(ny) <= nabs)) {

			affine1[0][2] = affines_init[1][0][2] + pix_step * pix_size * nx;

			affine1[1][2] = affines_init[1][1][2] + pix_step * pix_size * ny;

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

				System.arraycopy(affines_init[1][i], 0, affine1[i], 0, affine1[i].length);

			}

			affine1[0][2] += pix_step * pix_size * nx;

			affine1[1][2] += pix_step * pix_size * ny;

			pairwiseOrthoMatch.rms = Double.NaN;

			

			correlateOrthoPair(

					clt_parameters,  // CLTParameters    clt_parameters,

					pairwiseOrthoMatch, //PairwiseOrthoMatch pairwiseOrthoMatch, // will return statistics
@@ -1262,6 +1266,7 @@ public class OrthoMapsCollection implements Serializable{ 
		double    frac_below =    0.1;

		double [][] jtj = null;

		double rms = Double.NaN;

		double last_center_radius = 0;

		for (int ntry = 0; ntry < num_tries; ntry++) {

			if (!batch_mode) {

				if (debugLevel>-3) {
@@ -1500,6 +1505,15 @@ public class OrthoMapsCollection implements Serializable{ 
					min_std_rad,  // double        min_std_rad,  // minimal radius of the central area (if less - fail)



					debugLevel);     // final int     debug_level)

			last_center_radius = orthoPairLMA.getCenterRadius();

			if ((min_std_rad >0) && (last_center_radius < 1)) {

				if (debugLevel>-3) {

					System.out.println("correlateOrthoPair(): center_radius="+

							orthoPairLMA.getCenterRadius()+" < "+ 1);

				}

				return null;

			}

			

			if (num_good_tiles < min_good_tiles) {

				if (debugLevel>-4) {

					System.out.println("correlateOrthoPair(): num_good_tiles="+
@@ -1563,6 +1577,14 @@ public class OrthoMapsCollection implements Serializable{ 
			prev_rms=rms;

		} // for (int ntry = 0; ntry < num_tries; ntry++) {

		

		if ((min_std_rad > 0) && (last_center_radius < min_std_rad)) {

			if (debugLevel>-3) {

				System.out.println("correlateOrthoPair(): center_radius="+

						last_center_radius+" < "+ min_std_rad);

			}

			return null;

		}

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

+194 −3
Original line number Diff line number Diff line
@@ -37,8 +37,8 @@ import Jama.Matrix; 


public class OrthoPairLMA {

	private boolean           thread_invariant= true;

	private int               N =               0;

	private Rectangle         woi =             null;

	private int               N =               0; // number of tiles in WOI:  woi.width * woi.height

	private Rectangle         woi =             null; // in tiles

	private int               width;

	private double []         last_rms =        null; // {rms, rms_pure}, matching this.vector

	private double []         good_or_bad_rms = null; // just for diagnostics, to read last (failed) rms
@@ -53,11 +53,15 @@ public class OrthoPairLMA { 
	private double []         last_ymfx =       null;

	private double [][]       last_jt =         null;

	private boolean           origin_center =   false; // true - origin in overlap center, false - top left corner (as it was)

	private double []         origin =          null; // either {0,0} for top-left or center of the woi

	private double []         origin =          null; // in pixels either {0,0} for top-left or center of the woi

	private double            center_radius =   0;    // center zone radius (in tiles) that has limited standard deviation

	public  int               num_good_tiles = 0;

	public OrthoPairLMA (boolean origin_center) {

		this.origin_center = origin_center;

	}

	public double getCenterRadius() {

		return center_radius;

	}

	public int prepareLMA(

			// will always calculate relative affine, starting with unity

			int           width, // tilesX
@@ -95,6 +99,24 @@ public class OrthoPairLMA { 
		}

		N = woi.width * woi.height;

		parameters_vector = new double [] {1,0,0,1,0,0};

		if (min_std_rad > 0) {

			int         min_tiles = 4;

			getCenterRadius(

					max_std,       // final double      max_std,      // maximal standard deviation to limit center area  

					min_tiles,     // final int         min_tiles,

					vector_XYS,    // final double [][] vector_XYS,

					weights_extra, // final double []   weights_extra, // null or additional weights (such as elevation-based)

					centers);      // final double [][] centers)

//			if (getCenterRadius() < min_std_rad) { // only after last pass of LMA!

//				num_good_tiles = 0;

//				return num_good_tiles;

//			}

			weights_extra = applyRadialWeights( // uses this.center_radius;

					vector_XYS,    // final double [][] vector_XYS,

					weights_extra, // final double []   weights_extra,

					centers);      // final double [][] centers)	

		}

		

		setSamplesWeightsYCenters(

				vector_XYS,

				weights_extra, // null or additional weights (such as elevation-based)
@@ -137,6 +159,175 @@ public class OrthoPairLMA { 
//	lambda, // *10, // temporary

//	this.last_jt)



	private double [] applyRadialWeights( // uses this.center_radius;

			final double [][] vector_XYS,

			final double []   weights_extra,

			final double [][] centers) { // null or additional weights (such as elevation-based)

		return applyRadialWeights(

				center_radius, // final double      radius_tiles,

				vector_XYS,    // final double [][] vector_XYS,

				weights_extra, // final double []   weights_extra,  // null or additional weights (such as elevation-based)

				centers);      // final double [][] centers);

	}

	private double [] applyRadialWeights(

			final double      radius_tiles,

			final double [][] vector_XYS,

			final double []   weights_extra,  // null or additional weights (such as elevation-based)

			final double [][] centers) {

		if (Double.isInfinite(radius_tiles)) {

			return weights_extra; // may be null - should be OK

		}

		final double radius_pix = radius_tiles * GPUTileProcessor.DTT_SIZE;

		final double radius_pix2 = radius_pix*radius_pix;

		final double [] weights = new double [vector_XYS.length];

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

		final AtomicInteger ai = new AtomicInteger(0);

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

			threads[ithread] = new Thread() {

				public void run() {

					for (int iTile = ai.getAndIncrement(); iTile < N; iTile = ai.getAndIncrement()) {

						int tileX = iTile % woi.width + woi.x;

						int tileY = iTile / woi.width + woi.y;

						int aTile = tileY * width + tileX;

						if ((vector_XYS[aTile] != null) && (centers[aTile] != null)) {

							double w = vector_XYS[aTile][2];

							if (weights_extra != null) w *= weights_extra[aTile]; 

							if (Double.isNaN(w)) w = 0;

							double dx = centers[aTile][0] - origin[0];

							double dy = centers[aTile][1] - origin[1];

							double r_pix2 = dx*dx+dy*dy; // radius in pixels, squared

							if (r_pix2 < 4*radius_pix2) {

								double r_pix = Math.sqrt(r_pix2); // radius in pixels

								double wr = 0.5*(1.0 + Math.cos(0.5*Math.PI*r_pix/radius_pix));

								weights[aTile] = w*wr;

							}

						}

					}

				}

			};

		}		      

		ImageDtt.startAndJoin(threads);

		return weights;

	}

	

	

	

	/**

	 * Calculate center radius (in tiles) to use for LMA when initial rotation or

	 * scaling is inaccurate so the peripheral areas do not fit into correlation range.

	 * In that case try to adjust only the central area first, then increase that area

	 * in next iterations. 

	 * @param max_std       maximal standard deviation (average for X and Y) inside center area

	 * @param min_tiles     minimal tiles in the center zone

	 * @param vector_XYS    2D correlation-measured X,Y, and strength

	 * @param weights_extra null or optional additional weights of the samples

	 * @param centers       per tile (matching vector_XYS and weights_extra) tile centers in pixels

	 * @return maximal radius from the center (in pixels) where standard deviation of the inside samples is

	 *         below max_std. Returns Double.POSITIVE_INFINITY if std is not reached

	 */

	private double getCenterRadius(

	        final double      max_std,      // maximal standard deviation to limit center area  

//	        final double      min_std_rad,  // minimal radius of the central area (if less - fail)

	        final int         min_tiles,

			final double [][] vector_XYS,

			final double []   weights_extra, // null or additional weights (such as elevation-based)

			final double [][] centers) {

		final int rad_max = Math.max(woi.width,woi.height)/2 +1; 

		final int rad_length = rad_max +1; 

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

		final AtomicInteger ai = new AtomicInteger(0);

		final AtomicInteger ati = new AtomicInteger(0);

		final double [][] s0_arr =  new double [threads.length][rad_length];

		final double [][] sx_arr =  new double [threads.length][rad_length];

		final double [][] sx2_arr = new double [threads.length][rad_length];

		final double [][] sy_arr =  new double [threads.length][rad_length];

		final double [][] sy2_arr = new double [threads.length][rad_length];

		final int    [][] sn_arr =  new int [threads.length][rad_length];

		

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

			threads[ithread] = new Thread() {

				public void run() {

					int thread_num = ati.getAndIncrement();

					for (int iTile = ai.getAndIncrement(); iTile < N; iTile = ai.getAndIncrement()) {

						int tileX = iTile % woi.width + woi.x;

						int tileY = iTile / woi.width + woi.y;

						int aTile = tileY * width + tileX;

						if ((vector_XYS[aTile] != null) && (centers[aTile] != null)) {

							double w = vector_XYS[aTile][2];

							if (weights_extra != null) w *= weights_extra[aTile]; 

							if (Double.isNaN(w)) w = 0;

							double dx = centers[aTile][0] - origin[0];

							double dy = centers[aTile][1] - origin[1];

							double r_t = Math.sqrt(dx*dx+dy*dy)/GPUTileProcessor.DTT_SIZE; // radius in tiles

							int irt = (int) Math.round(r_t);

							if (irt < rad_length) {

								double fx = vector_XYS[aTile][0];

								double fy = vector_XYS[aTile][1];

								sn_arr [thread_num][irt] += 1;

								s0_arr [thread_num][irt] += w;

								sx_arr [thread_num][irt] += w * fx;

								sx2_arr[thread_num][irt] += w * fx*fx;

								sy_arr [thread_num][irt] += w * fy;

								sy2_arr[thread_num][irt] += w * fy * fy;

							}

						}

					}

				}

			};

		}		      

		ImageDtt.startAndJoin(threads);

		final double [] s0 =  new double [rad_length];

		final double [] sx =  new double [rad_length];

		final double [] sx2 = new double [rad_length];

		final double [] sy =  new double [rad_length];

		final double [] sy2 = new double [rad_length];

		final int []    sn =  new int [rad_length];

		ai.set(0);

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

			threads[ithread] = new Thread() {

				public void run() {

					for (int iRad = ai.getAndIncrement(); iRad < rad_length; iRad = ai.getAndIncrement()) {

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

							s0 [iRad]+=s0_arr [i][iRad];

							sx [iRad]+=sx_arr [i][iRad];

							sx2[iRad]+=sx2_arr[i][iRad];

							sy [iRad]+=sy_arr [i][iRad];

							sy2[iRad]+=sy2_arr[i][iRad];

							sn [iRad]+=sn_arr [i][iRad];

						}

					}

				}

			};

		}		      

		ImageDtt.startAndJoin(threads);

		double sc0=0, scx=0,scx2=0,scy=0,scy2=0,std_prev=0, std=0;

		int    scn=0;

		center_radius = Double.POSITIVE_INFINITY; 

		for (int iRad=0; iRad < rad_length; iRad++) {

			sc0+= s0 [iRad];

			scx+= sx [iRad];

			scx2+=sx2[iRad];

			scy+= sy [iRad];

			scy2+=sy2[iRad];

			scn+= sn [iRad];

			if ((scn > min_tiles) && (sc0 > 0)) {// for one tile gets scrt() of a small negative error

				std = Math.sqrt((scx2*sc0 - scx*scx + scy2*sc0 - scy*scy)/2)/sc0;

			}

			if ((scn >= min_tiles) && (std >= max_std)) {

				if (iRad==0) {

					center_radius = 0;

					break; // probably can not happen

				}

				center_radius = iRad - (std - max_std) / (std - std_prev);

				break;

			}

			std_prev = std;

		}

		return center_radius;

	}

	

	

	

	private void setSamplesWeightsYCenters(

			final double [][] vector_XYS,

			final double []   weights_extra, // null or additional weights (such as elevation-based)