Implementing scene high frequency decay calculation/reporting (b70fdc39) · Commits · Elphel / imagej-elphel

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

+23 −3

Original line number	Diff line number	Diff line
		@@ -2424,7 +2424,27 @@ public class DoubleFHT {
		return amp;
		}

		public double[] calculateAmplitudeNoSwap(double[] fht) {
		public double [] getFreqAmplitude(
		double [] data) {
		updateMaxN(data);
		swapQuadrants(data);
		if (!transform(data, false))
		return null; // direct FHT
		double [] amp = calculateAmplitude(data);
		return amp;
		}
		public double [] getFreqAmplitude2(
		double [] data) {
		updateMaxN(data);
		swapQuadrants(data);
		if (!transform(data, false))
		return null; // direct FHT
		double [] amp = calculateAmplitude2(data);
		return amp;
		}


		public static double[] calculateAmplitudeNoSwap(double[] fht) {
		int size = (int) Math.sqrt(fht.length);
		double[] amp = new double[size * size];
		for (int row = 0; row < size; row++) {
		@@ -2473,7 +2493,7 @@ public class DoubleFHT {
		}

		/* Amplitude of one row from 2D Hartley Transform. */
		void amplitude(int row, int size, double[] fht, double[] amplitude) {
		static void amplitude(int row, int size, double[] fht, double[] amplitude) {
		int base = row * size;
		int l;
		for (int c = 0; c < size; c++) {
		@@ -2495,7 +2515,7 @@ public class DoubleFHT {
		}

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

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

+341 −0

Original line number	Diff line number	Diff line
		@@ -3531,6 +3531,347 @@ public class OrthoMap implements Comparable <OrthoMap>, Serializable{
		return null;
		}

		/**
		* Trying to estimate image OTF to modify correlation results. Some images are better, some - worse
		* (blurred because of elevation errors)? For all image or parts of it? So on some images all
		* real objects (and false ones) get higher correlation, on some - all get lower. So some compensation
		* on image quality may help to discriminate
		* @param data image to process (may have NaNs)
		* @param width image width
		* @param size FFT size (now 128)
		* @param center_period period (in pixels) corresponding to the frequency to measure OTF derivative
		* @param range_period relative frequency range to average: low band from center/range_period to center,
		* high band - from center to center*range_period
		* @param wh if not null, should be int[2] - will return {tilesX,tilesY} for the result
		* @param debugLevel
		* @return per tile: null or a pair of high_frequency_response/low_frequency_response (around center)
		* for horizontal and vertical directions
		*/
		public static double [][] getHiFreq(
		final double [] data,
		final int width,
		final int size, // power of 2, such as 64
		final double center_period,// center frequency is size/center_period
		final double range_period, // ~1.5 - from center/range to center*range
		final int [] wh, // result size
		final int debugLevel){
		final int dbg_x = -2668;
		final int dbg_y = 256;

		final int height = data.length/width;
		final int tilesX = (int) Math.ceil(width/(size/2)) + 1;
		final int tilesY = (int) Math.ceil(height/(size/2)) + 1;
		if (wh != null) {
		wh[0] = tilesX;
		wh[1] = tilesY;
		}
		final double [][] hi_feq = new double [tilesX*tilesY][];
		int center_freq = (int) Math.round(size/center_period);
		int low_freq = (int) Math.round(size/center_period/range_period);
		int high_freq = (int) Math.round(size/center_period*range_period);
		final int [][][] ranges = { // {first, second},{low freq, high freq}, {start, end}
		{ {size/2 - center_freq + 1, size/2 - low_freq},
		{size/2 - high_freq, size/2 - center_freq - 1}
		},
		{{size/2 + low_freq, size/2 + center_freq - 1},
		{size/2 + center_freq + 1, size/2 + high_freq}}};
		final double [] range_npix = {
		size*(center_freq-low_freq),
		size*(high_freq-center_freq)};
		final double [] window = new double [size*size];
		final double [] wnd1d = new double[size/2];
		for (int i = 0; i < size/2; i++) {
		wnd1d[i] = Math.sin((i+0.5)*Math.PI/size);
		}
		for (int i = 0; i < size/2; i++) {
		for (int j = 0; j < size/2; j++) {
		double w = wnd1d[i]*wnd1d[j];
		int k = i*size + j;
		window[k] = w;
		window[window.length-1-k] = w;
		k = (i+ 1)*size - 1 -j;
		window[k] = w;
		window[window.length-1-k] = w;
		}
		}

		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() {
		double [] dtile = new double [size*size];
		TileNeibs tn = new TileNeibs(size,size);
		DoubleFHT doubleFHT = new DoubleFHT();
		for (int nTile = ai.getAndIncrement(); nTile <hi_feq.length; nTile = ai.getAndIncrement()){
		int tileX = nTile % tilesX;
		int tileY = nTile / tilesX;
		int px0 = (size/2) * tileX; // absolute in the original/result image
		int py0 = (size/2) * tileY;
		int x0 = Math.max(0, -px0);
		int y0 = Math.max(0, -py0);
		int x1 = Math.min(size, width- px0);
		int y1 = Math.min(size, height-py0);
		boolean dbg_tile = (Math.abs((px0 + size/2) - dbg_x) < size/4) && (Math.abs((py0 + size/2) - dbg_y) < size/4);
		if (dbg_tile) {
		System.out.println("getHiFreq(): tileX="+tileX+", tileY="+tileY);
		System.out.println("getHiFreq(): px0="+px0+", py0="+py0);
		}
		int lwidth=x1-x0;
		boolean has_NaN = false;
		if ((x0>0) \|\| (y0>0) \|\| (x1 < size) \|\| (y1 < size)) {
		Arrays.fill(dtile,Double.NaN);
		}
		for (int y = y0; y < y1; y++) {
		System.arraycopy(
		data,
		(py0 + y)*width+(px0+x0),
		dtile,
		y*size+x0,
		lwidth);
		}
		for (int i = 0; i < dtile.length; i++) {
		if (Double.isNaN(dtile[i])) {
		has_NaN = true;
		break;
		}
		}
		if (has_NaN) {
		fillNaNs(dtile, tn, 3);
		}
		for (int i = 0; i < dtile.length; i++) {
		dtile[i] *= window[i];
		}
		if (dbg_tile) {
		String [] rslt_titles= {"windowed"};
		ShowDoubleFloatArrays.showArrays(
		new double[][] {dtile},
		size,
		size,
		true,
		"windowed_data_tx"+tileX+"_ty"+tileY,
		rslt_titles);
		}
		double [] amp = doubleFHT.getFreqAmplitude(dtile);
		if (dbg_tile) {
		String [] rslt_titles= {"amplitude"};
		ShowDoubleFloatArrays.showArrays(
		new double[][] {amp},
		size,
		size,
		true,
		"amplitude_tx"+tileX+"_ty"+tileY,
		rslt_titles);
		}
		double [][] lo_hi_avg = new double[2][2]; // {x,y}{low,high}
		for (int hl = 0; hl < 2; hl++) { // 0 - low, 1 - high
		for (int sf = 0; sf < 2; sf++) { // 0 - fist, 1 second range
		for (int i = ranges[sf][hl][0]; i <=ranges[sf][hl][1]; i++) {
		for (int j = 0; j < size; j++) {
		lo_hi_avg[0][hl] += amp[j*size + i];
		lo_hi_avg[1][hl] += amp[i*size + j];
		}
		}
		}
		}
		hi_feq[nTile] = new double[2];
		for (int yx = 0; yx < 2; yx++) { // 0 - y, 1 - x
		for (int hl = 0; hl < 2; hl++) { // 0 - low, 1 - high
		lo_hi_avg[yx][hl] /= range_npix[hl];
		}
		hi_feq[nTile][yx] = lo_hi_avg[yx][1]/lo_hi_avg[yx][0];
		}
		}
		}
		};
		}
		ImageDtt.startAndJoin(threads);
		return hi_feq;
		}


		/**
		* Similar to above, but calculates for 2 rings
		* @param data
		* @param width
		* @param size
		* @param center_period
		* @param range_period
		* @param blank_xy discard data long x and y axes (probably remaining row/column noise?)
		* @param wh
		* @param debugLevel
		* @return
		*/
		public static double [][] getHiFreqCirc(
		final double [] data,
		final int width,
		final int size, // power of 2, such as 64
		final double center_period,// center frequency is size/center_period
		final double range_period, // ~1.5 - from center/range to center*range
		final int blank_xy, //
		final int [] wh, // result size
		final int debugLevel){
		final int dbg_x = 1144; // -2668;
		final int dbg_y = 199; // 256;

		final int height = data.length/width;
		final int tilesX = (int) Math.ceil(width/(size/2)) + 1;
		final int tilesY = (int) Math.ceil(height/(size/2)) + 1;
		if (wh != null) {
		wh[0] = tilesX;
		wh[1] = tilesY;
		}
		final double [][] hi_feq = new double [tilesX*tilesY][];
		int center_freq = (int) Math.round(size/center_period);
		int low_freq = (int) Math.round(size/center_period/range_period);
		int high_freq = (int) Math.round(size/center_period*range_period);
		final double blank2 = (blank_xy-1) * (blank_xy-1) + 0.5;
		final double [][] masks= new double [2][size*size];
		for (int y = 0; y<size; y++) {
		double y2 = (y-size/2);
		y2*=y2;
		if ((blank_xy == 0) \|\| (y2 > blank2)) {
		for (int x = 0; x < size; x++) {
		double x2 = (x-size/2);
		x2*=x2;
		if ((blank_xy == 0) \|\| (x2 > blank2)) {
		double r = Math.sqrt(x2+y2);
		if ((r >= low_freq) && (r <= high_freq)) {
		int indx = y*size + x;
		if (r < center_freq) {
		masks[0][indx] = Math.sin(Math.PI *(center_freq - r)/(center_freq-low_freq));
		} else {
		masks[1][indx] = Math.sin(Math.PI *(r - center_freq)/(high_freq -center_freq));
		}
		}
		}
		}
		}
		}
		if ((dbg_x >=0) && (dbg_y >=0)) {
		String [] rslt_titles= {"low_mask","high_mask"};
		ShowDoubleFloatArrays.showArrays(
		masks,
		size,
		size,
		true,
		"masks",
		rslt_titles);
		}
		for (int n = 0; n < 2; n++) {
		double s=0.0;
		for (int i = 0; i < masks[n].length; i++) {
		s+=masks[n][i];
		}
		s = 1/s;
		for (int i = 0; i < masks[n].length; i++) {
		masks[n][i]*=s;
		}
		}
		final double [] window = new double [size*size];
		final double [] wnd1d = new double[size/2];
		for (int i = 0; i < size/2; i++) {
		wnd1d[i] = Math.sin((i+0.5)*Math.PI/size);
		}
		for (int i = 0; i < size/2; i++) {
		for (int j = 0; j < size/2; j++) {
		double w = wnd1d[i]*wnd1d[j];
		int k = i*size + j;
		window[k] = w;
		window[window.length-1-k] = w;
		k = (i+ 1)*size - 1 -j;
		window[k] = w;
		window[window.length-1-k] = w;
		}
		}

		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() {
		double [] dtile = new double [size*size];
		TileNeibs tn = new TileNeibs(size,size);
		DoubleFHT doubleFHT = new DoubleFHT();
		for (int nTile = ai.getAndIncrement(); nTile <hi_feq.length; nTile = ai.getAndIncrement()){
		int tileX = nTile % tilesX;
		int tileY = nTile / tilesX;
		int px0 = (size/2) * tileX; // absolute in the original/result image
		int py0 = (size/2) * tileY;
		int x0 = Math.max(0, -px0);
		int y0 = Math.max(0, -py0);
		int x1 = Math.min(size, width- px0);
		int y1 = Math.min(size, height-py0);
		boolean dbg_tile = (Math.abs((px0 + size/2) - dbg_x) < size/4) && (Math.abs((py0 + size/2) - dbg_y) < size/4);
		if (dbg_tile) {
		System.out.println("getHiFreq(): tileX="+tileX+", tileY="+tileY);
		System.out.println("getHiFreq(): px0="+px0+", py0="+py0);
		}
		int lwidth=x1-x0;
		boolean has_NaN = false;
		if ((x0>0) \|\| (y0>0) \|\| (x1 < size) \|\| (y1 < size)) {
		Arrays.fill(dtile,Double.NaN);
		}
		for (int y = y0; y < y1; y++) {
		System.arraycopy(
		data,
		(py0 + y)*width+(px0+x0),
		dtile,
		y*size+x0,
		lwidth);
		}
		for (int i = 0; i < dtile.length; i++) {
		if (Double.isNaN(dtile[i])) {
		has_NaN = true;
		break;
		}
		}
		if (has_NaN) {
		continue;
		// fillNaNs(dtile, tn, 3);
		}
		for (int i = 0; i < dtile.length; i++) {
		dtile[i] *= window[i];
		}
		if (dbg_tile) {
		String [] rslt_titles= {"windowed"};
		ShowDoubleFloatArrays.showArrays(
		new double[][] {dtile},
		size,
		size,
		true,
		"windowed_data_tx"+tileX+"_ty"+tileY,
		rslt_titles);
		}
		double [] amp2 = doubleFHT.getFreqAmplitude2(dtile);
		if (dbg_tile) {
		String [] rslt_titles= {"amplitude2","low_mask","high_mask"};
		ShowDoubleFloatArrays.showArrays(
		new double[][] {amp2,masks[0],masks[1]},
		size,
		size,
		true,
		"amplitude2_tx"+tileX+"_ty"+tileY,
		rslt_titles);
		}
		hi_feq[nTile] = new double[masks.length];
		for (int n = 0; n < masks.length; n++) {
		for (int i = 0; i < masks[n].length; i++) {
		hi_feq[nTile][n] += masks[n][i] * amp2[i];
		}
		hi_feq[nTile][n] = Math.sqrt(hi_feq[nTile][n]);
		}
		}
		}
		};
		}
		ImageDtt.startAndJoin(threads);
		return hi_feq;
		}




		public static double [] correlateWithPattern(
		final double [] data,
		final int width,

Admin message