Creating kernel and matching images from different altitude (1c62dfa5) · Commits · Elphel / imagej-elphel

src/main/java/com/elphel/imagej/common/DoubleFHT.java

+75 −2

Original line number	Original line	Diff line number	Diff line
	@@ -346,7 +346,7 @@ public class DoubleFHT {
	return phaseCorrelate(first, phaseCoeff, filter,fht_save);		return phaseCorrelate(first, phaseCoeff, filter,fht_save);
	}		}

	public double [] phaseCorrelate (		public double [] phaseCorrelate ( // never
	double [] first,		double [] first,
	double phaseCoeff,		double phaseCoeff,
	double [] filter) { // high/low pass filtering		double [] filter) { // high/low pass filtering
	@@ -561,7 +561,7 @@ public class DoubleFHT {
	if (debug) ShowDoubleFloatArrays.showArrays(this.translateFHT, "translateFHT-"+IJ.d2s(dx,3)+":"+IJ.d2s(dy,3));		if (debug) ShowDoubleFloatArrays.showArrays(this.translateFHT, "translateFHT-"+IJ.d2s(dx,3)+":"+IJ.d2s(dy,3));
	}		}

	private boolean updateMaxN(double [] data){		public boolean updateMaxN(double [] data){ // was private
	if (data==null) return false; // do nothing		if (data==null) return false; // do nothing
	if (!powerOf2Size(data)) {		if (!powerOf2Size(data)) {
	String msg="Image is not power of 2 size";		String msg="Image is not power of 2 size";
	@@ -2146,6 +2146,41 @@ public class DoubleFHT {
	return result;		return result;
	}		}

			public double [] divide(double [] h1, double [] h2, double fat_zero) {
			int rowMod, colMod;
			double mag, h2e, h2o;
			double fz2=fat_zero*fat_zero;
			double[] result = new double[maxN*maxN];
			for (int r=0; r<maxN; r++) {
			rowMod = (maxN - r) % maxN;
			for (int c=0; c<maxN; c++) {
			colMod = (maxN - c) % maxN;
			mag =h2[rmaxN+c] h2[rmaxN+c] + h2[rowModmaxN+colMod] * h2[rowMod*maxN+colMod]+fz2;
			if (mag<1e-20)
			mag = 1e-20;
			h2e = (h2[rmaxN+c] + h2[rowModmaxN+colMod]);
			h2o = (h2[rmaxN+c] - h2[rowModmaxN+colMod]);
			double tmp = (h1[rmaxN+c] h2e - h1[rowModmaxN+colMod] h2o);
			result[r*maxN+c] = tmp/mag;
			}
			}
			return result;
			}

			public double [] setReal (double [] amp) { // only first half used
			double[] result = new double[maxN*maxN];
			int rowMod, colMod;
			for (int r=0; r<maxN/2; r++) {
			rowMod = (maxN - r) % maxN;
			for (int c=0; c<maxN; c++) {
			colMod = (maxN - c) % maxN;
			result[rmaxN+c] = amp[rmaxN+c];
			result[rowModmaxN+colMod] = amp[rmaxN+c];
			}
			}
			return result;
			}


	public double [] calculateAmplitude(double [] fht) {		public double [] calculateAmplitude(double [] fht) {
	int size=(int) Math.sqrt(fht.length);		int size=(int) Math.sqrt(fht.length);
	@@ -2156,6 +2191,31 @@ public class DoubleFHT {
	swapQuadrants(amp);		swapQuadrants(amp);
	return amp;		return amp;
	}		}
			public double [] calculateAmplitudeNoSwap(double [] fht) {
			int size=(int) Math.sqrt(fht.length);
			double[] amp = new double[size*size];
			for (int row=0; row<size; row++) {
			amplitude(row, size, fht, amp);
			}
			return amp;
			}
			public static double [] calculatePhaseNoSwap(double [] fht) {
			int size=(int) Math.sqrt(fht.length);
			double[] phs = new double[size*size];
			for (int row=0; row<size; row++) {
			phase(row, size, fht, phs);
			}
			return phs;
			}
			public double [] calculatePhase(double [] fht) {
			int size=(int) Math.sqrt(fht.length);
			double[] phs = new double[size*size];
			for (int row=0; row<size; row++) {
			phase(row, size, fht, phs);
			}
			swapQuadrants(phs);
			return phs;
			}

	public double [] calculateAmplitudeHalf(double [] fht) {		public double [] calculateAmplitudeHalf(double [] fht) {
	int size=(int) Math.sqrt(fht.length);		int size=(int) Math.sqrt(fht.length);
	@@ -2186,6 +2246,19 @@ public class DoubleFHT {
	}		}
	}		}

			static void phase(int row, int size, double[] fht, double[] phase) {
			int base = row*size;
			int l;
			for (int c=0; c<size; c++) {
			l = ((size-row)%size) * size + (size-c)%size;
			double re=0.5*(fht[base+c]+fht[l]);
			double im=0.5*(fht[base+c]-fht[l]);
			phase[base+c] = Math.atan2(im,re);;
			}
			}



	/* Squared amplitude of one row from 2D Hartley Transform. */		/* Squared amplitude of one row from 2D Hartley Transform. */
	void amplitude2(int row, int size, double[] fht, double[] amplitude) {		void amplitude2(int row, int size, double[] fht, double[] amplitude) {
	int base = row*size;		int base = row*size;

src/main/java/com/elphel/imagej/correction/Eyesis_Correction.java

+16 −1

Original line number	Original line	Diff line number	Diff line
	@@ -849,6 +849,7 @@ public class Eyesis_Correction implements PlugIn, ActionListener {
	addButton("Read Tiff", panelLWIRWorld, color_process);		addButton("Read Tiff", panelLWIRWorld, color_process);
	addButton("Set pair GPS", panelLWIRWorld, color_process);		addButton("Set pair GPS", panelLWIRWorld, color_process);
	addButton("Test video", panelLWIRWorld, color_process);		addButton("Test video", panelLWIRWorld, color_process);
			addButton("Deconvolve Slices", panelLWIRWorld, color_process);
	plugInFrame.add(panelLWIRWorld);		plugInFrame.add(panelLWIRWorld);

	}		}
	@@ -5733,8 +5734,22 @@ public class Eyesis_Correction implements PlugIn, ActionListener {
	return;		return;
	}		}
	OrthoMap.testVideo(imp_sel);		OrthoMap.testVideo(imp_sel);
			} else if (label.equals("Deconvolve Slices")) {
			ImagePlus imp_sel = WindowManager.getCurrentImage();
			if (imp_sel == null) {
			IJ.showMessage("Error", "No images selected");
			return;
			}
			OrthoMap.testDeconvolveSlices(
			imp_sel, // ImagePlus imp,
			512, // int [] slices,
			new int [] {1,2}, // int [] slices, deconvolve slice 2 with slice1
			4, // int kernel_radius,
			true, // boolean hor_sym,
			true, // boolean vert_sym,
			true, // false, // boolean all_sym,
			DEBUG_LEVEL); // int debugLevel) { // >0
	}		}

	}		}

	public boolean debugInitOneScene() {		public boolean debugInitOneScene() {

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

+426 −0

Original line number	Original line	Diff line number	Diff line
	@@ -2,6 +2,7 @@ package com.elphel.imagej.orthomosaic;

	import java.awt.Color;		import java.awt.Color;
	import java.awt.Font;		import java.awt.Font;
			import java.awt.Rectangle;
	import java.io.File;		import java.io.File;
	import java.io.IOException;		import java.io.IOException;
	import java.io.ObjectInputStream;		import java.io.ObjectInputStream;
	@@ -22,6 +23,7 @@ import java.util.Properties;
	import java.util.concurrent.atomic.AtomicInteger;		import java.util.concurrent.atomic.AtomicInteger;

	import com.elphel.imagej.cameras.CLTParameters;		import com.elphel.imagej.cameras.CLTParameters;
			import com.elphel.imagej.common.DoubleFHT;
	import com.elphel.imagej.common.GenericJTabbedDialog;		import com.elphel.imagej.common.GenericJTabbedDialog;
	import com.elphel.imagej.common.PolynomialApproximation;		import com.elphel.imagej.common.PolynomialApproximation;
	import com.elphel.imagej.common.ShowDoubleFloatArrays;		import com.elphel.imagej.common.ShowDoubleFloatArrays;
	@@ -39,6 +41,7 @@ import ij.ImagePlus;
	import ij.ImageStack;		import ij.ImageStack;
	import ij.Prefs;		import ij.Prefs;
	import ij.gui.PointRoi;		import ij.gui.PointRoi;
			import ij.gui.Roi;
	import ij.plugin.filter.AVI_Writer;		import ij.plugin.filter.AVI_Writer;
	import ij.plugin.filter.GaussianBlur;		import ij.plugin.filter.GaussianBlur;
	import ij.process.ColorProcessor;		import ij.process.ColorProcessor;
	@@ -759,6 +762,429 @@ public class OrthoMap implements Comparable <OrthoMap>, Serializable{
	ImageDtt.startAndJoin(threads);		ImageDtt.startAndJoin(threads);
	return padded_gpu;		return padded_gpu;
	}		}
			/**
			* Extracts ROI, rounds to fft_size and deconvolves 2-nd slice (lower resolution) with
			* the first one (higher resolution)
			* @param imp
			*/
			public static void testDeconvolveSlices(
			ImagePlus imp,
			final int fft_size,
			int [] slices,
			int kernel_radius,
			boolean hor_sym,
			boolean vert_sym,
			boolean all_sym,
			int debugLevel) { // >0

			ImageStack stack = imp.getStack();
			final int width = stack.getWidth();
			final int height = stack.getHeight();
			Roi roi= imp.getRoi();
			boolean good_ROI = false;
			Rectangle rroi = null;
			if (roi != null) {
			good_ROI = roi.getType()== Roi.RECTANGLE;
			}
			if (good_ROI) {
			rroi=roi.getBounds();
			if ((rroi.width != fft_size) \|\| (rroi.height != fft_size)) {
			good_ROI=false;
			}
			}
			if (slices == null) {
			slices = new int [] {1,2};
			}
			// String kernel_path = "/media/elphel/NVME/lwir16-proc/ortho_videos/kernel_50_75.tiff";
			String kernel_path = "/media/elphel/NVME/lwir16-proc/ortho_videos/kernel_50_100.tiff";

			double [] kernel = null;
			if (!good_ROI) {
			ImagePlus imp_kernel = new ImagePlus(kernel_path);
			if (imp_kernel.getWidth() == 0) {
			System.out.println("testDeconvolveSlices(): precomputed kernel "+kernel_path+
			" is not found, and for calculation");
			System.out.println("it needs rectangular selection of "+fft_size+"x"+fft_size);
			return;
			} else {
			float [] kernel_pixels = (float[]) imp_kernel.getProcessor().getPixels();
			kernel = new double[kernel_pixels.length];
			for (int i = 0; i < kernel.length; i++) {
			kernel[i] = kernel_pixels[i];
			}
			}
			}
			final Thread[] threads = ImageDtt.newThreadArray(QuadCLT.THREADS_MAX);
			final AtomicInteger ai = new AtomicInteger(0);
			if (kernel==null) {
			roi= imp.getRoi(); // retry roi
			rroi=roi.getBounds();
			if (rroi.width != fft_size) {
			rroi.width = fft_size;
			}
			if (rroi.height != fft_size) {
			rroi.height = fft_size;
			}
			final double [][] dpixels = new double [2][fft_size*fft_size];
			// ImageStack stack = imp.getStack();
			// final int width = stack.getWidth();
			// final int height = stack.getHeight();
			final int tl = rroi.y*width+rroi.x; // top left corner index
			for (int n = 0; n < 2; n++) {
			final int fn = n;
			final float [] fpixels = (float[]) stack.getPixels(slices[fn]);
			ai.set(0);
			for (int ithread = 0; ithread < threads.length; ithread++) {
			threads[ithread] = new Thread() {
			public void run() {
			for (int ipix = ai.getAndIncrement(); ipix < dpixels[fn].length; ipix = ai.getAndIncrement()) {
			int ix = ipix % fft_size;
			int iy = ipix / fft_size;
			dpixels[fn][ipix] = fpixels[tl+iy*width+ix];
			}
			}
			};
			}
			ImageDtt.startAndJoin(threads);
			}
			double fat_zero = 1000; // 300 - too small
			double [] deconvolved= deconvolvePair(
			dpixels,
			fat_zero,
			debugLevel);
			if (debugLevel>0) {
			double [][] test_img = new double [][] {dpixels[0], dpixels[1], deconvolved}; // dpixels modified - DC, window
			String [] test_titles= {"high_res","low_res","deconvolved"};
			ShowDoubleFloatArrays.showArrays(
			test_img,
			fft_size,
			fft_size,
			true,
			removeKnownExtension(imp.getTitle())+"deconvolved",
			test_titles);
			}

			kernel = extractKernel(
			deconvolved, // double [] data,
			kernel_radius, // int kernel_radius,
			hor_sym, // boolean hor_sym,
			vert_sym, // boolean vert_sym,
			all_sym, // boolean all_sym,
			debugLevel); // int debugLevel)
			}
			final double [][] full_img = new double [4][width*height];
			final int [] out_slices = {0,2,1,3}; // {high-res, low-res, convoled_high_res,data}
			for (int n = 0; n < 2; n++) {
			final int fn = out_slices[n];

			final float [] fpixels = (float[]) stack.getPixels(slices[n]);
			ai.set(0);
			for (int ithread = 0; ithread < threads.length; ithread++) {
			threads[ithread] = new Thread() {
			public void run() {
			for (int ipix = ai.getAndIncrement(); ipix < full_img[fn].length; ipix = ai.getAndIncrement()) {
			full_img[fn][ipix] = fpixels[ipix];
			}
			}
			};
			}
			ImageDtt.startAndJoin(threads);
			}


			full_img[out_slices[2]] = convolveWithKernel(
			full_img[out_slices[0]], // final double [] data,
			kernel, // final double [] kernel,
			width); // final int width)
			// calculate DATI
			ai.set(0);
			for (int ithread = 0; ithread < threads.length; ithread++) {
			threads[ithread] = new Thread() {
			public void run() {
			for (int ipix = ai.getAndIncrement(); ipix < full_img[0].length; ipix = ai.getAndIncrement()) {
			full_img[out_slices[3]][ipix] = full_img[out_slices[1]][ipix]-full_img[out_slices[2]][ipix];
			}
			}
			};
			}
			ImageDtt.startAndJoin(threads);
			if (debugLevel>0) {
			String [] test_titles= {"high_res","convolved", "low_res", "DATI"};
			ShowDoubleFloatArrays.showArrays(
			full_img,
			width,
			height,
			true,
			removeKnownExtension(imp.getTitle())+"-DATI",
			test_titles);
			}
			System.out.println("testDeconvolveSlices() Done");
			}

			public static double [] convolveWithKernel(
			final double [] data,
			final double [] kernel,
			final int width) {
			final int height = data.length / width;
			final double [] convolved = new double [data.length];
			Arrays.fill(convolved, Double.NaN);
			final int kernel_size = (int) Math.sqrt(kernel.length); // supposed to be odd
			final int kernel_radius = (kernel_size-1)/2;

			final Thread[] threads = ImageDtt.newThreadArray(QuadCLT.THREADS_MAX);
			final AtomicInteger ai = new AtomicInteger(0);
			final int kernel_center= kernel_radius * (kernel_size + 1);
			for (int ithread = 0; ithread < threads.length; ithread++) {
			threads[ithread] = new Thread() {
			public void run() {
			for (int ipix = ai.getAndIncrement(); ipix < data.length; ipix = ai.getAndIncrement()) if (!Double.isNaN(data[ipix])){
			int ix = ipix % width;
			int iy = ipix / width;
			int dx_min=-kernel_radius,dx_max=kernel_radius, dy_min=-kernel_radius,dy_max=kernel_radius;
			if (ix < kernel_radius) dx_max = ix;
			if (iy < kernel_radius) dy_max = iy;
			if (ix >= (width - kernel_radius)) dx_min = ix - width + 1;
			if (iy >= (height - kernel_radius)) dy_min = iy - height + 1;
			double swd = 0, sw = 0;
			for (int dy = dy_min; dy <= dy_max; dy++) {
			for (int dx = dx_min; dx <= dx_max; dx++) {
			int src_dindex = ipix - dy * width - dx;
			double d = data[src_dindex];
			if (!Double.isNaN(d)) {
			int src_kindex = kernel_center + dy * kernel_size + dx;
			double k = kernel [src_kindex];
			sw += k;
			swd += d * k;
			}
			convolved[ipix] = swd/sw;
			}
			}
			}
			}
			};
			}
			ImageDtt.startAndJoin(threads);
			return convolved;
			}

			public static double [] extractKernel(
			double [] data,
			int kernel_radius,
			boolean hor_sym,
			boolean vert_sym,
			boolean all_sym,
			int debugLevel) {
			hor_sym \|= all_sym;
			vert_sym \|= all_sym;
			double [][] dbg_img = (debugLevel > 0)? new double [4][]: null;
			int kernel_size = 2*kernel_radius + 1;
			int fft_size = (int) Math.sqrt(data.length);
			double [] kernel = new double [kernel_size*kernel_size];
			int data_tl = (fft_size/2 - kernel_radius) * ( fft_size + 1);
			for (int row = 0; row < kernel_size; row++) {
			System.arraycopy(
			data,
			data_tl + row * fft_size,
			kernel,
			row* kernel_size,
			kernel_size);
			}
			double sum = 0;
			for (int i = 0; i < kernel.length; i++) {
			sum += kernel[i];
			}
			double s = 1.0/sum;
			for (int i = 0; i < kernel.length; i++) {
			kernel[i] *= s;
			}
			if (dbg_img != null) dbg_img[0] = kernel.clone();
			if (hor_sym) {
			for (int row = 0; row < kernel_size; row++) {
			for (int col = 1; col <= kernel_radius; col++) {
			int indx = row*kernel_size + kernel_radius; // center of the line
			double d = 0.5*(kernel[indx-col] + kernel[indx+col]);
			kernel[indx-col] = d;
			kernel[indx+col] = d;
			}
			}
			}
			if (dbg_img != null) dbg_img[1] = kernel.clone();
			if (vert_sym) {
			for (int row = 1; row <= kernel_radius; row++) {
			int indx0= (kernel_radius-row)*kernel_size;
			int indx1= (kernel_radius+row)*kernel_size;
			for (int col = 0; col < kernel_size; col++) {
			double d = 0.5*(kernel[indx0+col] + kernel[indx1+col]);
			kernel[indx0+col] = d;
			kernel[indx1+col] = d;
			}
			}
			}
			if (dbg_img != null) dbg_img[2] = kernel.clone();
			if (all_sym) {
			for (int row = 0; row < (kernel_size-1); row++) {
			for (int col = row+1; col < kernel_size; col++) {
			int indx0 = row * kernel_size + col;
			int indx1 = col * kernel_size + row;
			double d = 0.5*(kernel[indx0] + kernel[indx1]);
			kernel[indx0] = d;
			kernel[indx1] = d;
			}
			}
			}
			if (dbg_img != null) dbg_img[3]= kernel.clone();
			if (dbg_img != null) {
			String [] dbg_titles= {"orig","hor","vert","all"};
			ShowDoubleFloatArrays.showArrays(
			dbg_img,
			kernel_size,
			kernel_size,
			true,
			"kernel",
			dbg_titles);
			}
			return kernel;
			}
			public static double [] deconvolvePair(
			double [][] data,
			double fat_zero, // 1000
			int debugLevel) {
			boolean zero_phase=true;
			int fft_size = (int) Math.sqrt(data[0].length);
			DoubleFHT doubleFHT = new DoubleFHT();
			doubleFHT.updateMaxN(data[1]);
			double [] w1d = doubleFHT.getHamming1d(fft_size);
			for (int n= 0; n < data.length; n++) {
			removeDC(data[n]);
			multiplySquareBy1D(data[n], w1d);
			removeDC(data[n]); // once more?
			}
			double [][] data_bkp = null;
			if (debugLevel > 0) {
			data_bkp = new double[][] {data[0].clone(),data[1].clone()};
			}
			doubleFHT.swapQuadrants(data[0]);
			doubleFHT.swapQuadrants(data[1]);
			if (!doubleFHT.transform(data[1],false)) return null; // direct FHT
			if (!doubleFHT.transform(data[0],false)) return null; // direct FHT {
			if (debugLevel > 0) {
			double [] amp = doubleFHT.calculateAmplitude(data[0]);
			ShowDoubleFloatArrays.showArrays(
			amp,
			fft_size,
			fft_size,
			"amp");
			}
			double [] deconv = doubleFHT.divide(data[1],data[0], fft_size);
			if (debugLevel > 0) {
			ShowDoubleFloatArrays.showArrays(
			deconv,
			fft_size,
			fft_size,
			"deconv");
			}
			// if (zero_phase) {
			double [] div_amp = doubleFHT.calculateAmplitudeNoSwap(deconv);
			double [] div_phase = doubleFHT.calculatePhaseNoSwap(deconv);
			double [][] amp_phase = {div_amp, div_phase};
			ShowDoubleFloatArrays.showArrays(
			amp_phase,
			fft_size,
			fft_size,
			true,
			"deconvolved-amp-phase",
			new String[] {"amp","phase"});
			/*
			if (debugLevel > 0) {
			ShowDoubleFloatArrays.showArrays(
			div_amp,
			fft_size,
			fft_size,
			"div_amp");
			}
			deconv = doubleFHT.setReal(div_amp);
			if (debugLevel > 0) {
			ShowDoubleFloatArrays.showArrays(
			deconv,
			fft_size,
			fft_size,
			"deconv_real");
			}
			*/

			// }

			// try with zero phase


			//if (filter!=null) multiplyByReal(first, filter); // add filter if needed
			doubleFHT.transform(deconv,true) ; // inverse transform
			doubleFHT.swapQuadrants(deconv);
			if (debugLevel > 0) {
			data[0]= data_bkp[0];
			data[1]= data_bkp[1];
			}
			return deconv;
			}


			public static void removeDC(
			final double [] data) {
			final Thread[] threads = ImageDtt.newThreadArray(QuadCLT.THREADS_MAX);
			final AtomicInteger ai = new AtomicInteger(0);
			final AtomicInteger ati = new AtomicInteger(0);
			final double [] avg_arr = new double [threads.length];
			final double [] npix_arr = new double [threads.length];
			for (int ithread = 0; ithread < threads.length; ithread++) {
			threads[ithread] = new Thread() {
			public void run() {
			int thread_num = ati.getAndIncrement();
			for (int iPix = ai.getAndIncrement(); iPix < data.length; iPix = ai.getAndIncrement()) if (!Double.isNaN(data[iPix])){
			avg_arr[thread_num] += data[iPix];
			npix_arr[thread_num] += 1.0;
			}
			}
			};
			}
			ImageDtt.startAndJoin(threads);
			double avg=0, num=0;
			for (int i = 0; i < avg_arr.length; i++) {
			avg+=avg_arr[i];
			num+=npix_arr[i];
			}
			final double favg = avg/num;
			ai.set(0);
			for (int ithread = 0; ithread < threads.length; ithread++) {
			threads[ithread] = new Thread() {
			public void run() {
			for (int iPix = ai.getAndIncrement(); iPix < data.length; iPix = ai.getAndIncrement()){
			data[iPix] -= favg;
			}
			}
			};
			}
			ImageDtt.startAndJoin(threads);
			}

			public static void multiplySquareBy1D(
			final double [] data,
			final double[] wnd1d) {
			final int width = wnd1d.length;
			final Thread[] threads = ImageDtt.newThreadArray(QuadCLT.THREADS_MAX);
			final AtomicInteger ai = new AtomicInteger(0);
			for (int ithread = 0; ithread < threads.length; ithread++) {
			threads[ithread] = new Thread() {
			public void run() {
			for (int iPix = ai.getAndIncrement(); iPix < data.length; iPix = ai.getAndIncrement()) if (!Double.isNaN(data[iPix])){
			int x = iPix % width;
			int y = iPix / width;
			data[iPix] = wnd1d[x] wnd1d[y];
			}
			}
			};
			}
			ImageDtt.startAndJoin(threads);
			}


	public static void testVideo(ImagePlus imp) {		public static void testVideo(ImagePlus imp) {

Admin message