Making intersceneLma produce same results regardless of threads (decd04ed) · Commits · Elphel / imagej-elphel

pom.xml

+8 −0

Original line number	Original line	Diff line number	Diff line
	@@ -155,6 +155,14 @@
	<artifactId>ant-launcher</artifactId>		<artifactId>ant-launcher</artifactId>
	<version>1.10.10</version>		<version>1.10.10</version>
	</dependency>		</dependency>

			<!-- https://mvnrepository.com/artifact/javax.xml.bind/jaxb-api -->
			<dependency>
			<groupId>javax.xml.bind</groupId>
			<artifactId>jaxb-api</artifactId>
			<version>2.3.1</version>
			</dependency>

	</dependencies>		</dependencies>

	<build>		<build>

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

+9 −0

Original line number	Original line	Diff line number	Diff line
	@@ -250,6 +250,15 @@ public class ErsCorrection extends GeometryCorrection {
	this.ers_wxyz_center_dt = ers_xyz_dt;		this.ers_wxyz_center_dt = ers_xyz_dt;
	this.ers_watr_center_dt = ers_atr_dt;		this.ers_watr_center_dt = ers_atr_dt;
	}		}
			public void setErsDt_test(
			double [] ers_xyz_dt,
			double [] ers_atr_dt) {
			double k = 1.0; // 0.5;
			this.ers_wxyz_center_dt = new double[] {-ers_xyz_dt[0],-ers_xyz_dt[1],-ers_xyz_dt[2]};
			this.ers_watr_center_dt = new double[] {k* ers_atr_dt[0],-kers_atr_dt[1],kers_atr_dt[2]};//ers_atr_dt;
			}


	public void setErsD2t(		public void setErsD2t(
	double [] ers_xyz_d2t,		double [] ers_xyz_d2t,
	double [] ers_atr_d2t) {		double [] ers_atr_d2t) {

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

+225 −37

Original line number	Original line	Diff line number	Diff line
	package com.elphel.imagej.tileprocessor;		package com.elphel.imagej.tileprocessor;

			import java.io.ByteArrayOutputStream;
			import java.io.IOException;
			import java.io.ObjectOutputStream;
			import java.io.Serializable;
			import java.security.MessageDigest;
			import java.security.NoSuchAlgorithmException;
	import java.util.ArrayList;		import java.util.ArrayList;
	import java.util.concurrent.atomic.AtomicInteger;		import java.util.concurrent.atomic.AtomicInteger;
	import java.util.concurrent.atomic.DoubleAdder;		import java.util.concurrent.atomic.DoubleAdder;

			import javax.xml.bind.DatatypeConverter;

	import Jama.Matrix;		import Jama.Matrix;

	public class IntersceneLma {		public class IntersceneLma {
	@@ -27,12 +35,17 @@ public class IntersceneLma {
	private double [][] macrotile_centers = null; // (will be used to pull for regularization)		private double [][] macrotile_centers = null; // (will be used to pull for regularization)
	private double infinity_disparity = 0.1; // treat lower as infinity		private double infinity_disparity = 0.1; // treat lower as infinity
	private int num_samples = 0;		private int num_samples = 0;
			private boolean thread_invariant = true; // Do not use DoubleAdder, provide results not dependent on threads
	public IntersceneLma(		public IntersceneLma(
	OpticalFlow opticalFlow		OpticalFlow opticalFlow,
			boolean thread_invariant
	) {		) {
			this.thread_invariant = thread_invariant;
	this.opticalFlow = opticalFlow;		this.opticalFlow = opticalFlow;
	}		}
			public double[] getLastRms() {
			return last_rms;
			}
	public double [] getSceneXYZ(boolean initial) {		public double [] getSceneXYZ(boolean initial) {
	double [] full_vector = initial? backup_parameters_full: getFullVector(parameters_vector);		double [] full_vector = initial? backup_parameters_full: getFullVector(parameters_vector);
	return new double[] {		return new double[] {
	@@ -170,6 +183,9 @@ public class IntersceneLma {
	setSamplesWeights(vector_XYS); // not regularization yet !		setSamplesWeights(vector_XYS); // not regularization yet !

	last_jt = new double [parameters_vector.length][];		last_jt = new double [parameters_vector.length][];
			if (debug_level > 1) {
			System.out.println("prepareLMA() 1");
			}
	double [] fx = getFxDerivs(		double [] fx = getFxDerivs(
	parameters_vector, // double [] vector,		parameters_vector, // double [] vector,
	last_jt, // final double [][] jt, // should be null or initialized with [vector.length][]		last_jt, // final double [][] jt, // should be null or initialized with [vector.length][]
	@@ -185,6 +201,10 @@ public class IntersceneLma {
	}		}
	normalizeWeights(); // make full weight == 1.0; pure_weight <= 1.0;		normalizeWeights(); // make full weight == 1.0; pure_weight <= 1.0;
	// remeasure fx - now with regularization terms.		// remeasure fx - now with regularization terms.

			if (debug_level > 1) {
			System.out.println("prepareLMA() 2");
			}
	fx = getFxDerivs(		fx = getFxDerivs(
	parameters_vector, // double [] vector,		parameters_vector, // double [] vector,
	last_jt, // final double [][] jt, // should be null or initialized with [vector.length][]		last_jt, // final double [][] jt, // should be null or initialized with [vector.length][]
	@@ -292,6 +312,9 @@ public class IntersceneLma {
	// maybe the following if() branch is not needed - already done in prepareLMA !		// maybe the following if() branch is not needed - already done in prepareLMA !
	if (this.last_rms == null) { //first time, need to calculate all (vector is valid)		if (this.last_rms == null) { //first time, need to calculate all (vector is valid)
	last_rms = new double[2];		last_rms = new double[2];
			if (debug_level > 1) {
			System.out.println("lmaStep(): first step");
			}
	double [] fx = getFxDerivs(		double [] fx = getFxDerivs(
	parameters_vector, // double [] vector,		parameters_vector, // double [] vector,
	last_jt, // final double [][] jt, // should be null or initialized with [vector.length][]		last_jt, // final double [][] jt, // should be null or initialized with [vector.length][]
	@@ -331,6 +354,19 @@ public class IntersceneLma {
	Matrix wjtjlambda = new Matrix(getWJtJlambda(		Matrix wjtjlambda = new Matrix(getWJtJlambda(
	lambda, // *10, // temporary		lambda, // *10, // temporary
	this.last_jt)); // double [][] jt)		this.last_jt)); // double [][] jt)
			if (debug_level > 1) {
			try {
			System.out.println("getFxDerivs(): getChecksum(this.y_vector)="+ getChecksum(this.y_vector));
			System.out.println("getFxDerivs(): getChecksum(this.weights)="+ getChecksum(this.weights));
			System.out.println("getFxDerivs(): getChecksum(this.last_ymfx)="+ getChecksum(this.last_ymfx));
			System.out.println("getFxDerivs(): getChecksum(y_minus_fx_weighted)="+getChecksum(y_minus_fx_weighted));
			System.out.println("getFxDerivs(): getChecksum(wjtjlambda)= "+getChecksum(wjtjlambda));
			} catch (NoSuchAlgorithmException \| IOException e) {
			// TODO Auto-generated catch block
			e.printStackTrace();
			}
			}

	if (debug_level>2) {		if (debug_level>2) {
	System.out.println("JtJ + lambda*diag(JtJ");		System.out.println("JtJ + lambda*diag(JtJ");
	wjtjlambda.print(18, 6);		wjtjlambda.print(18, 6);
	@@ -356,6 +392,17 @@ public class IntersceneLma {
	System.out.println("Jt * (y-fx)");		System.out.println("Jt * (y-fx)");
	jty.print(18, 6);		jty.print(18, 6);
	}		}
			if (debug_level > 1) {
			try {
			System.out.println("getFxDerivs(): getChecksum(jtjl_inv)="+getChecksum(jtjl_inv));
			System.out.println("getFxDerivs(): getChecksum(jty)= "+getChecksum(jty));
			} catch (NoSuchAlgorithmException \| IOException e) {
			// TODO Auto-generated catch block
			e.printStackTrace();
			}
			}



	Matrix mdelta = jtjl_inv.times(jty);		Matrix mdelta = jtjl_inv.times(jty);
	if (debug_level>2) {		if (debug_level>2) {
	@@ -369,6 +416,19 @@ public class IntersceneLma {
	for (int i = 0; i < parameters_vector.length; i++) {		for (int i = 0; i < parameters_vector.length; i++) {
	new_vector[i] += scale * delta[i];		new_vector[i] += scale * delta[i];
	}		}
			if (debug_level > 1) {
			try {
			System.out.println("getFxDerivs(): getChecksum(mdelta)= "+getChecksum(mdelta));
			System.out.println("getFxDerivs(): getChecksum(delta)= "+getChecksum(delta));
			System.out.println("getFxDerivs(): getChecksum(parameters_vector)= "+getChecksum(parameters_vector));
			System.out.println("getFxDerivs(): getChecksum(new_vector)= "+getChecksum(new_vector));
			} catch (NoSuchAlgorithmException \| IOException e) {
			// TODO Auto-generated catch block
			e.printStackTrace();
			}
			}


	double [] fx = getFxDerivs(		double [] fx = getFxDerivs(
	new_vector, // double [] vector,		new_vector, // double [] vector,
	last_jt, // final double [][] jt, // should be null or initialized with [vector.length][]		last_jt, // final double [][] jt, // should be null or initialized with [vector.length][]
	@@ -444,6 +504,27 @@ public class IntersceneLma {

	final Thread[] threads = ImageDtt.newThreadArray(opticalFlow.threadsMax);		final Thread[] threads = ImageDtt.newThreadArray(opticalFlow.threadsMax);
	final AtomicInteger ai = new AtomicInteger(0);		final AtomicInteger ai = new AtomicInteger(0);
			double sum_weights;
			if (thread_invariant) {
			final double [] sw_arr = new double [vector_XYS.length];
			for (int ithread = 0; ithread < threads.length; ithread++) {
			threads[ithread] = new Thread() {
			public void run() {
			for (int iMTile = ai.getAndIncrement(); iMTile < vector_XYS.length; iMTile = ai.getAndIncrement()) if (vector_XYS[iMTile] != null){
			double w = vector_XYS[iMTile][2];
			weights[2 * iMTile] = w;
			// asum_weight.add(w);
			sw_arr[iMTile] = w;
			}
			}
			};
			}
			ImageDtt.startAndJoin(threads);
			sum_weights = 0.0;
			for (double w:sw_arr) {
			sum_weights += w;
			}
			} else {
	final DoubleAdder asum_weight = new DoubleAdder();		final DoubleAdder asum_weight = new DoubleAdder();
	for (int ithread = 0; ithread < threads.length; ithread++) {		for (int ithread = 0; ithread < threads.length; ithread++) {
	threads[ithread] = new Thread() {		threads[ithread] = new Thread() {
	@@ -457,8 +538,12 @@ public class IntersceneLma {
	};		};
	}		}
	ImageDtt.startAndJoin(threads);		ImageDtt.startAndJoin(threads);
			sum_weights = asum_weight.sum();
			}

	ai.set(0);		ai.set(0);
	final double s = 0.5/asum_weight.sum();		// final double s = 0.5/asum_weight.sum();
			final double s = 0.5/sum_weights;
	for (int ithread = 0; ithread < threads.length; ithread++) {		for (int ithread = 0; ithread < threads.length; ithread++) {
	threads[ithread] = new Thread() {		threads[ithread] = new Thread() {
	public void run() {		public void run() {
	@@ -473,7 +558,9 @@ public class IntersceneLma {
	pure_weight = 1.0;		pure_weight = 1.0;
	}		}

	private void normalizeWeights()
			@Deprecated
			private void normalizeWeights_old()
	{		{
	//num_samples		//num_samples
	final Thread[] threads = ImageDtt.newThreadArray(opticalFlow.threadsMax);		final Thread[] threads = ImageDtt.newThreadArray(opticalFlow.threadsMax);
	@@ -514,6 +601,58 @@ public class IntersceneLma {

	}		}

			private void normalizeWeights()
			{
			final Thread[] threads = ImageDtt.newThreadArray(opticalFlow.threadsMax);
			final AtomicInteger ai = new AtomicInteger(0);
			double full_weight, sum_weight_pure;
			if (thread_invariant) {
			sum_weight_pure = 0;
			for (int i = 0; i < num_samples; i++) {
			sum_weight_pure += weights[i];
			}
			full_weight = sum_weight_pure;
			for (int i = 0; i < par_indices.length; i++) {
			int indx = num_samples + i;
			full_weight += weights[indx];
			}
			} else {
			final DoubleAdder asum_weight = new DoubleAdder();
			for (int ithread = 0; ithread < threads.length; ithread++) {
			threads[ithread] = new Thread() {
			public void run() {
			for (int i = ai.getAndIncrement(); i < num_samples; i = ai.getAndIncrement()){
			asum_weight.add(weights[i]);
			}
			}
			};
			}
			ImageDtt.startAndJoin(threads);
			sum_weight_pure = asum_weight.sum();
			for (int i = 0; i < par_indices.length; i++) {
			int indx = num_samples + i;
			asum_weight.add(weights[indx]);
			}
			full_weight = asum_weight.sum();
			}
			pure_weight = sum_weight_pure/full_weight;
			final double s = 1.0/full_weight;
			if (Double.isNaN(s)) {
			System.out.println("normalizeWeights(): s == NaN");
			}
			ai.set(0);
			for (int ithread = 0; ithread < threads.length; ithread++) {
			threads[ithread] = new Thread() {
			public void run() {
			for (int i = ai.getAndIncrement(); i < weights.length; i = ai.getAndIncrement()){
			weights[i] *= s;
			}
			}
			};
			}
			ImageDtt.startAndJoin(threads);
			}



	private double [] getFullVector(double [] vector) {		private double [] getFullVector(double [] vector) {
	@@ -615,6 +754,17 @@ public class IntersceneLma {
	fx [i + 2 * macrotile_centers.length] = vector[i]; // - parameters_initial[i]; // scale will be combined with weights		fx [i + 2 * macrotile_centers.length] = vector[i]; // - parameters_initial[i]; // scale will be combined with weights
	jt[i][i + 2 * macrotile_centers.length] = 1.0; // scale will be combined with weights		jt[i][i + 2 * macrotile_centers.length] = 1.0; // scale will be combined with weights
	}		}
			if (debug_level > 1) {
			try {
			System.out.println ("getFxDerivs(): getChecksum(fx)="+getChecksum(fx));
			if (jt != null) {
			System.out.println("getFxDerivs(): getChecksum(jt)="+getChecksum(jt));
			}
			} catch (NoSuchAlgorithmException \| IOException e) {
			// TODO Auto-generated catch block
			e.printStackTrace();
			}
			}
	return fx;		return fx;
	}		}

	@@ -661,9 +811,30 @@ public class IntersceneLma {
	) {		) {
	final Thread[] threads = ImageDtt.newThreadArray(opticalFlow.threadsMax);		final Thread[] threads = ImageDtt.newThreadArray(opticalFlow.threadsMax);
	final AtomicInteger ai = new AtomicInteger(0);		final AtomicInteger ai = new AtomicInteger(0);
	final DoubleAdder asum_weight = new DoubleAdder();
	final double [] wymfw = new double [fx.length];		final double [] wymfw = new double [fx.length];
			double s_rms;
			if (thread_invariant) {
			final double [] l2_arr = new double [num_samples];
			for (int ithread = 0; ithread < threads.length; ithread++) {
			threads[ithread] = new Thread() {
			public void run() {
			for (int i = ai.getAndIncrement(); i < num_samples; i = ai.getAndIncrement()) {
			double d = y_vector[i] - fx[i];
			double wd = d * weights[i];
			//double l2 = d * wd;
			l2_arr[i] = d * wd;
			wymfw[i] = wd;
			}
			}
			};
			}
			ImageDtt.startAndJoin(threads);
			s_rms = 0.0;
			for (double l2:l2_arr) {
			s_rms += l2;
			}
			} else {
			final DoubleAdder asum_weight = new DoubleAdder();
	for (int ithread = 0; ithread < threads.length; ithread++) {		for (int ithread = 0; ithread < threads.length; ithread++) {
	threads[ithread] = new Thread() {		threads[ithread] = new Thread() {
	public void run() {		public void run() {
	@@ -678,11 +849,11 @@ public class IntersceneLma {
	};		};
	}		}
	ImageDtt.startAndJoin(threads);		ImageDtt.startAndJoin(threads);
	double s_rms = asum_weight.sum();		s_rms = asum_weight.sum();
			}
	double rms_pure = Math.sqrt(s_rms/pure_weight);		double rms_pure = Math.sqrt(s_rms/pure_weight);
	for (int i = 0; i < par_indices.length; i++) {		for (int i = 0; i < par_indices.length; i++) {
	int indx = i + num_samples;		int indx = i + num_samples;
	// double d = parameters_initial[i] - fx[indx]; // fx[indx] == vector[i]
	double d = y_vector[indx] - fx[indx]; // fx[indx] == vector[i]		double d = y_vector[indx] - fx[indx]; // fx[indx] == vector[i]
	double wd = d * weights[indx];		double wd = d * weights[indx];
	s_rms += d * wd;		s_rms += d * wd;
	@@ -696,6 +867,23 @@ public class IntersceneLma {
	return wymfw;		return wymfw;
	}		}

			public static String getChecksum(Serializable object) throws IOException, NoSuchAlgorithmException {
			ByteArrayOutputStream baos = null;
			ObjectOutputStream oos = null;
			try {
			baos = new ByteArrayOutputStream();
			oos = new ObjectOutputStream(baos);
			oos.writeObject(object);
			MessageDigest md = MessageDigest.getInstance("MD5");
			byte[] thedigest = md.digest(baos.toByteArray());
			return DatatypeConverter.printHexBinary(thedigest);
			} finally {
			oos.close();
			baos.close();
			}
			}


	/*		/*
	* par_indices		* par_indices
	double [] rslt = {rms, rms_pure};		double [] rslt = {rms, rms_pure};

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

+32 −2

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

+509 −77

File changed.

Preview size limit exceeded, changes collapsed.

Original line number	Original line	Diff line number	Diff line
	@@ -29,15 +29,19 @@ import java.util.Properties;
	import com.elphel.imagej.common.GenericJTabbedDialog;		import com.elphel.imagej.common.GenericJTabbedDialog;

	public class IntersceneLmaParameters {		public class IntersceneLmaParameters {
			public boolean ilma_thread_invariant = true; // Do not use DoubleAdder, provide results not dependent on threads
	public boolean [] ilma_lma_select = new boolean [ErsCorrection.DP_NUM_PARS]; // first three will not be used		public boolean [] ilma_lma_select = new boolean [ErsCorrection.DP_NUM_PARS]; // first three will not be used
	public double [] ilma_regularization_weights = new double [ErsCorrection.DP_NUM_PARS]; // first three will not be used		public double [] ilma_regularization_weights = new double [ErsCorrection.DP_NUM_PARS]; // first three will not be used
			public boolean ilma_ignore_ers = false; // ignore linear and angular velocities, assume tham zeroes
			public boolean ilma_ers_adj_lin = false; // adjust linear ERS for scene-to-ref (LWIR does not work, check with high-speed)
			public boolean ilma_ers_adj_ang = false; // adjust angular ERS for scene-to-ref (LWIR does not work, check with high-speed)
	public double ilma_lambda = 0.1;		public double ilma_lambda = 0.1;
	public double ilma_lambda_scale_good = 0.5;		public double ilma_lambda_scale_good = 0.5;
	public double ilma_lambda_scale_bad = 8.0;		public double ilma_lambda_scale_bad = 8.0;
	public double ilma_lambda_max = 100;		public double ilma_lambda_max = 100;
	public double ilma_rms_diff = 0.001;		public double ilma_rms_diff = 0.001;
	public int ilma_num_iter = 20;		public int ilma_num_iter = 20;
	public int ilma_num_corr = 10; // maximal number of full correlatiobn+LMA cycles		public int ilma_num_corr = 10; // maximal number of full correlation+LMA cycles
	public int ilma_debug_level = 1;		public int ilma_debug_level = 1;

	public IntersceneLmaParameters() {		public IntersceneLmaParameters() {
	@@ -94,6 +98,8 @@ public class IntersceneLmaParameters {

	public void dialogQuestions(GenericJTabbedDialog gd) {		public void dialogQuestions(GenericJTabbedDialog gd) {
	gd.addMessage("Interframe LMA parameters selection");		gd.addMessage("Interframe LMA parameters selection");
			gd.addCheckbox ("Thread-invariant execution", this.ilma_thread_invariant,
			"Do not use DoubleAdder and provide results not dependent on threads" );
	for (int i = ErsCorrection.DP_DVAZ; i < ErsCorrection.DP_NUM_PARS; i++) {		for (int i = ErsCorrection.DP_DVAZ; i < ErsCorrection.DP_NUM_PARS; i++) {
	gd.addCheckbox (ErsCorrection.DP_DERIV_NAMES[i], this.ilma_lma_select[i],		gd.addCheckbox (ErsCorrection.DP_DERIV_NAMES[i], this.ilma_lma_select[i],
	"Adjust parameter "+ErsCorrection.DP_DERIV_NAMES[i]+" with interscene LMA" );		"Adjust parameter "+ErsCorrection.DP_DERIV_NAMES[i]+" with interscene LMA" );
	@@ -105,6 +111,14 @@ public class IntersceneLmaParameters {
	" will cause error equal to all reprojection ones");		" will cause error equal to all reprojection ones");
	}		}
	gd.addMessage("LMA other parameters");		gd.addMessage("LMA other parameters");
			gd.addCheckbox ("Ignore linear and angular velocities", this.ilma_ignore_ers,
			"Ignore calculated linear and angular velocities when correlating scenes to the reference one" );

			gd.addCheckbox ("Adjust linear ERS for scene-to-ref", this.ilma_ers_adj_lin,
			"Adjust linear velocities during LMA for scene-to-reference matching. So far does not work for LWIR (not stable - effect is samll)" );
			gd.addCheckbox ("Aadjust angular ERS for scene-to-ref", this.ilma_ers_adj_ang,
			"Adjust angular velocities during LMA for scene-to-reference matching. So far does not work for LWIR (not stable - effect is samll)" );

Admin message