package com.elphel.imagej.common; /** ** ** DoubleFHT - Determine simulation pattern parameters to match ** the acquired image ** ** Copyright (C) 2010-2011 Elphel, Inc. ** ** -----------------------------------------------------------------------------** ** ** MatchSimulatedPattern.java is free software: you can redistribute it and/or modify ** it under the terms of the GNU General Public License as published by ** the Free Software Foundation, either version 3 of the License, or ** (at your option) any later version. ** ** This program is distributed in the hope that it will be useful, ** but WITHOUT ANY WARRANTY; without even the implied warranty of ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the ** GNU General Public License for more details. ** ** You should have received a copy of the GNU General Public License ** along with this program. If not, see <http://www.gnu.org/licenses/>. ** -----------------------------------------------------------------------------** ** */ import java.util.ArrayList; import java.util.Arrays; import java.util.List; import ij.IJ; import ij.process.FHT; //import ij.plugin.FFT; //import ij.plugin.ContrastEnhancer; //import java.awt.image.ColorModel; //import ij.plugin.*; public class DoubleFHT { // private boolean isFrequencyDomain; private int maxSize = 20; private double[] freqMask = null; // frequency mask cache, will be reused if all parameters are the same // private int freqMaskN=0; // private double freqHighPass=0; // private double freqLowPass=0; private double[] freqPass = null; // {freqHighPass,freqLowPass}; private int maxN = -1; // undefined private int ln2 = -1; // undefined private double[] C; // TODO: make arrays for each reasonable FHT size (or use cache) private double[] S; private int[] bitrev; private double[] tempArr; public boolean showProgress = false; private double[][][] CS_cache = new double[maxSize][][]; // CS cache for different sizes, private double[][] freqMask_cache = new double[maxSize][]; // frequency mask cache, will be reused if all parameters // are the same private double[][] freqPass_cache = new double[maxSize][]; // cache for low/high frequencies private int[][] bitrev_cache = new int[maxSize][]; private double[][] hamming1d = new double[maxSize][]; // will cache 1d Hamming window; private double[] hamming1dMin = new double[maxSize]; // will cache 1d Hamming window; private double[][] gaussian1d = new double[maxSize][]; // will cache 1d Gaussian window; private double[] gaussian1dWidths = new double[maxSize]; // will cache 1d Gaussian window; private double[] translateFHT = null; // cached array for sub-pixel translate private double[] translateFHTdXY = null; private int translateFHTN = 0; public boolean debug = false; public DoubleFHT() { this.C = null; this.S = null; for (int i = 0; i < CS_cache.length; i++) { this.CS_cache[i] = null; this.freqMask_cache[i] = null; this.freqPass_cache[i] = null; this.bitrev_cache[i] = null; this.hamming1d[i] = null; this.gaussian1d[i] = null; this.gaussian1dWidths[i] = -1.0; } } public boolean powerOf2Size(double[] data) { return powerOf2Size(data.length); } public boolean powerOf2Size(int len) { int i = 4; while (i < len) i *= 4; return i == len; } public void swapQuadrants(double[] data) { if (data == null) return; int size = (int) Math.sqrt(data.length); int hsize = size / 2; int shift03 = (size + 1) * hsize; int shift12 = (size - 1) * hsize; int i, j, index; double d; for (i = 0; i < hsize; i++) for (j = 0; j < hsize; j++) { index = i * size + j; d = data[index]; data[index] = data[index + shift03]; data[index + shift03] = d; index += hsize; d = data[index]; data[index] = data[index + shift12]; data[index + shift12] = d; } } /* * Performs a forward transform, converting this image into the frequency * domain. The image contained in data must be square and its width must be a * power of 2. */ public boolean transform(double[] data) { return transform(data, false); } /* * Performs an inverse transform, converting this image into the space domain. * The image contained in data must be square and its width must be a power of * 2. */ public boolean inverseTransform(double[] data) { return transform(data, true); } /* * Returns an inverse transform of this image, which is assumed to be in the * frequency domain. */ public boolean transform(double[] data, boolean inverse) { // IJ.log("transform: "+maxN+" "+inverse); updateMaxN(data); rc2DFHT(data, inverse, this.maxN); return true; } public double[] createFrequencyFilter(double highPass, double lowPass) { return createFrequencyFilter(null, highPass, lowPass); } public double[] createFrequencyFilter(double[] data, // int n, double highPass, double lowPass) { if (data != null) updateMaxN(data); // int n; if ((this.freqMask != null) && (this.freqPass[0] == highPass) && (this.freqPass[1] == lowPass)) return this.freqMask; this.freqMask = new double[(this.maxN + 1) * this.maxN / 2 + 1]; double[] lo = new double[this.maxN]; double[] hi = new double[this.maxN]; int i, j; double kHi = (highPass > 0) ? (0.5 / highPass / highPass) : 0; double kLo = (lowPass > 0) ? (1.0 / lowPass / lowPass / this.maxN / this.maxN) : 0; double mx2; // System.out.println("createFrequencyFilter: highPass ="+highPass+ " lowPass= "+lowPass+ " kHi= "+kHi+ " kLo= "+kLo); for (i = 0; i < this.maxN; i++) { mx2 = -i * i; hi[i] = (kHi > 0) ? Math.exp(mx2 * kHi) : 0.0; lo[i] = (kLo > 0) ? Math.exp(mx2 * kLo) : 1.0; } // more precise version should add data from 4 corners, below is the simple version that may have discontinuities in the middle for (int index = 0; index < this.freqMask.length; index++) { i = index / this.maxN; j = index % this.maxN; if (j > this.maxN / 2) j = this.maxN - j; this.freqMask[index] = (1.0 - hi[i] * hi[j]) * lo[i] * lo[j]; } this.freqPass = new double[2]; this.freqPass[0] = highPass; this.freqPass[1] = lowPass; return this.freqMask; } /** * Shift the data by phase modification in the frequency domain * * @param data square 2**N array in line scan order. Data is modified in-place * @param dx amount of the horizontal shift (left) in pixels * @param dy amount of the vertical shift (down) in pixels * @return modified (shifted) data TODO: add shift+up-sample combination */ public double[] shift(double[] data, double dx, double dy) { return shift(data, 1, dx, dy); } /** * Upsample input array by padding in the frequency domain * * @param first input square power of 2 array * @param scale upsample scale (power of 2) * @return upsampled array, scan */ public double[] upsample(double[] first, int scale) { return shift(first, scale, 0.0, 0.0); } /** * Upsample data and shift it by phase modification in the frequency domain * * @param data square 2**N array in line scan order. Data is modified in-place * @param scale upsample scale * @param dx amount of the horizontal shift (left) in pixels (before scaling) * @param dy amount of the vertical shift (down) in pixels (before scaling) * @return modified (shifted) data */ public double[] shift(double[] data, int scale, double dx, double dy) { if ((scale == 0) && (dx == 0.0) && (dy == 0.0)) return data; updateMaxN(data); swapQuadrants(data); if (!transform(data, false)) return null; // direct FHT double[] result; // scale first if (scale > 1) { int halfN = this.maxN / 2; int shift = this.maxN * (scale - 1); int scaledN = this.maxN * scale; result = new double[data.length * scale * scale]; for (int i = 0; i < result.length; i++) result[i] = 0.0; double scale2 = scale * scale; for (int i = 0; i < data.length; i++) { int iy = i / this.maxN; int ix = i % this.maxN; if (ix > halfN) ix += shift; if (iy > halfN) iy += shift; result[scaledN * iy + ix] = scale2 * data[i]; } updateMaxN(result); } else { result = data; } // now shift if ((dx != 0.0) || (dy != 0.0)) { dx *= scale; dy *= scale; double sX = 2 * Math.PI * dx / this.maxN; double sY = 2 * Math.PI * dy / this.maxN; int halfN = this.maxN / 2; double[] cosDX = new double[this.maxN]; double[] sinDX = new double[this.maxN]; double[] cosDY = new double[halfN + 1]; double[] sinDY = new double[halfN + 1]; for (int i = 0; i <= halfN; i++) { // need less? cosDX[i] = Math.cos(sX * i); sinDX[i] = Math.sin(sX * i); cosDY[i] = Math.cos(sY * i); sinDY[i] = Math.sin(sY * i); } for (int i = 1; i < halfN; i++) { // need less? cosDX[this.maxN - i] = cosDX[i]; sinDX[this.maxN - i] = -sinDX[i]; } for (int row = 0; row <= halfN; row++) { int rowMod = (this.maxN - row) % this.maxN; int maxCol = (row < halfN) ? (this.maxN - 1) : halfN; for (int col = 0; col <= maxCol; col++) { int colMod = (this.maxN - col) % this.maxN; int index = row * this.maxN + col; int indexMod = rowMod * this.maxN + colMod; double re = 0.5 * (data[index] + data[indexMod]); double im = 0.5 * (data[index] - data[indexMod]); if ((col == halfN) || (row == halfN)) im = 0; double cosDelta = cosDX[col] * cosDY[row] - sinDX[col] * sinDY[row]; // cos(deltaX)*cos(deltaY)-sin(deltaX)*sin(deltaY) double sinDelta = sinDX[col] * cosDY[row] + cosDX[col] * sinDY[row]; // sin(deltaX)*cos(deltaY)+cos(deltaX)*sin(deltaY) double reMod = re * cosDelta - im * sinDelta; double imMod = re * sinDelta + im * cosDelta; data[index] = reMod + imMod; data[indexMod] = reMod - imMod; } } } if (!transform(result, true)) return null; // inverse FHT swapQuadrants(result); return result; } public double [] getFrequencyFilter ( double [] data, // just to find out size double highPassSigma, double lowPassSigma ) { double [] filter= null; if ((highPassSigma> 0.0) || (lowPassSigma> 0.0)) { updateMaxN(data); createFrequencyFilter(highPassSigma, lowPassSigma); // for repetitive calls will reuse mask filter =this.freqMask; } return filter; } /** * Invert kernel for deconvolutions. Kernel is assumed to fade near edges * @param data source kernel, square, power of 2 sides. Data is destroyed * @param fat_zero add to denominator to prevent noise amplification * @param highPassSigma high-pass sigma to create frequency filter or 0 * @param lowPassSigma low-pass sigma to create frequency filter or 0 * @return inverted kernel */ public double [] invert ( double [] data, double fat_zero, double highPassSigma, double lowPassSigma ) { double [] filter= null; if ((highPassSigma> 0.0) || (lowPassSigma> 0.0)) { updateMaxN(data); createFrequencyFilter(highPassSigma, lowPassSigma); // for repetitive calls will reuse mask filter =this.freqMask; } return invert (data, fat_zero, filter); } /** * Invert kernel for deconvolutions. Kernel is assumed to fade near edges * @param data source kernel, square, power of 2 sides. Data is destroyed * @param fat_zero add to denominator to prevent noise amplification * @param filter low-pass/high-pass filter or null * @return inverted kernel */ public double [] invert ( double [] data, // destroyed double fat_zero, double [] filter ) { updateMaxN(data); int size = (int) Math.sqrt(data.length); double nominator_FD_value = size*size; // fill FD of the nominator. what is the best value? double [] nominator_FD = new double[data.length]; Arrays.fill(nominator_FD, nominator_FD_value); swapQuadrants(data); if (!transform(data,false)) return null; // direct FHT double [] inverted = divide(nominator_FD,data, fat_zero); if (filter != null) { multiplyByReal(inverted, filter); } transform(inverted, true); // inverse transform swapQuadrants(inverted); return inverted; } /* calculates correlation, destroys original arrays */ public double[] correlate(double[] first, double[] second, double highPassSigma, double lowPassSigma, double phaseCoeff) { updateMaxN(first); createFrequencyFilter(highPassSigma, lowPassSigma); // for repetitive calls will reuse mask if (this.freqMask.length < first.length / 2) { System.out.println("Error: first.length=" + first.length + " second.length=" + second.length + " this.maxN=" + this.maxN + " this.freqMask.length=" + this.freqMask.length); } return correlate(first, second, phaseCoeff, this.freqMask); } public double[] correlate(double[] first, double[] second, double phaseCoeff) { return correlate(first, second, phaseCoeff, null); } public double[] correlate(double[] first, double[] second, double phaseCoeff, double[] filter) { if (phaseCoeff <= 0.0) return correlate(first, second, filter); else return phaseCorrelate(first, second, phaseCoeff, filter); } // asymmetrical - will divide by squared second amplitude (pattern to match) public double[] phaseCorrelate( // old double[] first, double[] second, double phaseCoeff, double[] filter) { // high/low pass filtering if (first.length != second.length) { IJ.showMessage("Error", "Correlation arrays should be the same size"); return null; } // System.out.println("phaseCorrelate, phaseCoeff="+phaseCoeff); swapQuadrants(first); swapQuadrants(second); if (!transform(first, false)) return null; // direct FHT if (!transform(second, false)) return null; // direct FHT first = phaseMultiply(first, second, phaseCoeff); // correlation, not convolution if (filter != null) multiplyByReal(first, filter); transform(first, true); // inverse transform swapQuadrants(first); return first; } public double[] phaseCorrelate(double[] first, double phaseCoeff, double high_pass, double low_pass) { // high/low // pass // filtering return phaseCorrelate(first, phaseCoeff, high_pass, low_pass, null); } public double[] phaseCorrelate(double[] first, double phaseCoeff, double high_pass, double low_pass, // high/low double[] fht_save) { // null-OK updateMaxN(first); double[] filter = null; if ((high_pass > 0) || (low_pass > 0)) { filter = createFrequencyFilter(high_pass, low_pass); } return phaseCorrelate(first, phaseCoeff, filter, fht_save); } public double[] phaseCorrelate( // never double[] first, double phaseCoeff, double[] filter) { // high/low pass filtering return phaseCorrelate(first, phaseCoeff, filter, null); } public double[] phaseCorrelate(double[] first, double phaseCoeff, double[] filter, // high/low pass filtering double[] fht_save) { // null-OK updateMaxN(first); swapQuadrants(first); if (!transform(first, false)) return null; // direct FHT if (fht_save != null) { System.arraycopy(first, 0, fht_save, 0, first.length); } first = phaseMultiplyNorm(first, first, phaseCoeff); // correlation, not convolution if (filter != null) multiplyByReal(first, filter); transform(first, true); // inverse transform swapQuadrants(first); return first; } public double[] phaseCorrelate(double[] first, double[] second, double phaseCoeff, double highPassSigma, double lowPassSigma, double[] first_save, double[] second_save) { // null-OK updateMaxN(first); double[] filter = null; if ((highPassSigma > 0) || (lowPassSigma > 0)) { filter = createFrequencyFilter(highPassSigma, lowPassSigma); } return phaseCorrelate(first, second, phaseCoeff, filter, // high/low pass filtering first_save, second_save); } public double[] phaseCorrelate( // new double[] first, double[] second, double phaseCoeff, double[] filter, // high/low pass filtering double[] first_save, double[] second_save) { // null-OK if (first.length != second.length) { IJ.showMessage("Error", "Correlation arrays should be the same size"); return null; } updateMaxN(first); swapQuadrants(first); swapQuadrants(second); if (!transform(first, false)) return null; // direct FHT if (!transform(second, false)) return null; // direct FHT if (first_save != null) System.arraycopy(first, 0, first_save, 0, first.length); if (second_save != null) System.arraycopy(second, 0, second_save, 0, second.length); first = phaseMultiplyNorm(first, second, phaseCoeff); // correlation, not convolution if (filter != null) multiplyByReal(first, filter); transform(first, true); // inverse transform swapQuadrants(first); return first; } /** * Transform pattern once, and then use FD of the pattern for each overlapping * image tile * * @param second pattern to be transform * @return transformed pattern ( the original is also modified) */ public double[] transformPattern( // new double[] second) { // null-OK updateMaxN(second); swapQuadrants(second); if (!transform(second, false)) return null; // direct FHT return second; } /** * Phase correlate with pattern, that is transformed separately with * transformPattern() * * @param first square array to be phase-correlated, modified but not with * result! * @param secondFD pattern already transformed to FD * @param phaseCoeff 0 - regular correlation, 1.0 - pure phase correlation * @param filter frequency filter or null * @param first_save null or array to accommodate FD of the first array before * updating it with correlation * @return phase correlation, first still contains original transformed! */ public double[] phaseCorrelatePattern( // new double[] first, double[] secondFD, double phaseCoeff, double[] filter, // high/low pass filtering double[] first_save) { // null-OK if (first.length != secondFD.length) { IJ.showMessage("Error", "Correlation arrays should be the same size"); return null; } updateMaxN(first); swapQuadrants(first); if (!transform(first, false)) return null; // direct FHT if (first_save != null) System.arraycopy(first, 0, first_save, 0, first.length); first = phaseMultiplyNorm(first, secondFD, phaseCoeff); // correlation, not convolution if (filter != null) multiplyByReal(first, filter); transform(first, true); // inverse transform swapQuadrants(first); return first; } /** * Convolve with pattern, that is transformed separately with transformPattern() * * @param first square array to be phase-correlated, modified but not with * result! * @param secondFD pattern already transformed to FD * @param filter frequency filter or null * @param first_save null or array to accommodate FD of the first array before * updating it with correlation * @return convolution result, first still contains original transformed! */ public double[] convolvePattern( // new double[] first, double[] secondFD, double[] filter, // high/low pass filtering double[] first_save) { // null-OK if (first.length != secondFD.length) { IJ.showMessage("Error", "Correlation arrays should be the same size"); return null; } updateMaxN(first); swapQuadrants(first); if (!transform(first, false)) return null; // direct FHT if (first_save != null) System.arraycopy(first, 0, first_save, 0, first.length); first = multiply(first, secondFD, false); // convolution, not correlation if (filter != null) multiplyByReal(first, filter); transform(first, true); // inverse transform swapQuadrants(first); return first; } // public double[] applyFreqFilter(double[] first, double[] filter) { if (filter == null) filter = this.freqMask; if (this.freqMask.length < first.length / 2) { System.out.println("Error in applyFreqFilter(): first.length=" + first.length + " this.maxN=" + this.maxN + " filter.length=" + filter.length); } multiplyByReal(first, filter); return first; } // modify - instead of a sigma - correlationHighPassSigma (0.0 - regular correlation, 1.0 - phase correlation) public double[] correlate(double[] first, double[] second) { return correlate(first, second, null); } public double[] correlate(double[] first, double[] second, double[] filter) { // high/low pass filtering // System.out.println("correlate"); if (first.length != second.length) { IJ.showMessage("Error", "Correlation arrays should be the same size"); return null; } swapQuadrants(first); swapQuadrants(second); if (!transform(first, false)) return null; // direct FHT if (!transform(second, false)) return null; // direct FHT first = multiply(first, second, true); // correlation, not convolution if (filter != null) multiplyByReal(first, filter); transform(first, true); // inverse transform swapQuadrants(first); return first; } public double[] convolve( double[] first, double[] second) { return convolve(first, second, null); } public double[] convolve( double[] first, double[] second, double[] filter) { // high/low pass filtering // System.out.println("correlate"); if (first.length != second.length) { IJ.showMessage("Error", "Correlation arrays should be the same size"); return null; } swapQuadrants(first); swapQuadrants(second); if (!transform(first, false)) return null; // direct FHT if (!transform(second, false)) return null; // direct FHT first = multiply(first, second, false); // convolution, not correlation if (filter != null) multiplyByReal(first, filter); transform(first, true); // inverse transform swapQuadrants(first); return first; } /** * Translate real array (zero in the center) by a sub-pixel vector (preferably * in +/- 0.5 range for each coordinate) by multiplying spectrum by the linear * phase gradient. Caches array, so multiple arrays can be shifted by the same * vector faster * * @param data square array to be translated by a sub-pixel value. in the end it * holds FHT (before multiplication)! * @param dx translation in x-direction * @param dy translation in y direction * @return modified data array, translated by specified vector (in-place) */ public double[] translateSubPixel(double[] data, double dx, double dy) { if (debug) ShowDoubleFloatArrays.showArrays(data, "source-" + IJ.d2s(dx, 3) + ":" + IJ.d2s(dy, 3)); swapQuadrants(data); if (!transform(data, false)) return null; // direct FHT if ((this.maxN != translateFHTN) || (dx != translateFHTdXY[0]) || (dy != translateFHTdXY[1])) { calcTranslateFHT(dx, dy); } if (debug) ShowDoubleFloatArrays.showArrays(data, "fht-source-" + IJ.d2s(dx, 3) + ":" + IJ.d2s(dy, 3)); data = multiply(data, this.translateFHT, false); // convolution, not correlation if (debug) ShowDoubleFloatArrays.showArrays(data, "fht-multiplied-" + IJ.d2s(dx, 3) + ":" + IJ.d2s(dy, 3)); transform(data, true); // inverse transform swapQuadrants(data); if (debug) ShowDoubleFloatArrays.showArrays(data, "result-" + IJ.d2s(dx, 3) + ":" + IJ.d2s(dy, 3)); return data; } /** * Calculating (and caching) FHT array for subpixel translation * * @param dx translation in x-direction * @param dy translation in y direction */ private void calcTranslateFHT(double dx, double dy) { double k = 2.0 * Math.PI / this.maxN; double kx = k * dx; double ky = k * dy; this.translateFHT = new double[this.maxN * this.maxN]; double sqrt2 = Math.sqrt(2.0); double pi4 = Math.PI / 4; int hSize = this.maxN / 2; double a, b, p; int rowMod, colMod; for (int r = 0; r < hSize; r++) { rowMod = (maxN - r) % maxN; for (int c = 0; c < hSize; c++) { colMod = (maxN - c) % maxN; p = kx * c + ky * r; a = sqrt2 * Math.sin(pi4 + p); b = sqrt2 * Math.sin(pi4 - p); this.translateFHT[r * this.maxN + c] = a; this.translateFHT[rowMod * this.maxN + colMod] = b; if (c == 0) { this.translateFHT[r * this.maxN + c + hSize] = a; this.translateFHT[rowMod * this.maxN + colMod + hSize] = b; } else { p = -kx * c + ky * r; a = sqrt2 * Math.sin(pi4 + p); b = sqrt2 * Math.sin(pi4 - p); this.translateFHT[r * this.maxN + colMod] = a; this.translateFHT[rowMod * this.maxN + c] = b; } } } // update cache this.translateFHTdXY = new double[2]; this.translateFHTdXY[0] = dx; this.translateFHTdXY[1] = dy; this.translateFHTN = this.maxN; if (debug) ShowDoubleFloatArrays.showArrays(this.translateFHT, "translateFHT-" + IJ.d2s(dx, 3) + ":" + IJ.d2s(dy, 3)); } public boolean updateMaxN(double[] data) { // was private if (data == null) return false; // do nothing if (!powerOf2Size(data)) { String msg = "Image is not power of 2 size"; IJ.showMessage("Error", msg); throw new IllegalArgumentException(msg); } int n = (int) Math.sqrt(data.length); boolean differentSize = (n != this.maxN); this.maxN = n; if (differentSize) { //System.out.println("Changing FFT size: old ln2="+this.ln2+" new maxN="+this.maxN); // save length old mask and mask parameters if (this.ln2 >= 0) { this.freqMask_cache[this.ln2] = this.freqMask; this.freqPass_cache[this.ln2] = this.freqPass; } // int ln2; for (this.ln2 = 0; maxN > (1 << this.ln2); this.ln2++) ; if ((this.ln2 >= CS_cache.length) || (CS_cache[this.ln2] == null)) { if (this.ln2 >= CS_cache.length) { String msg = "Too large image, increase this.maxSize (it is now " + this.maxSize + ", wanted " + this.ln2 + ")"; IJ.showMessage("Error", msg); throw new IllegalArgumentException(msg); } this.freqMask_cache[this.ln2] = null; this.freqPass_cache[this.ln2] = new double[2]; this.freqPass_cache[this.ln2][0] = 0.0; this.freqPass_cache[this.ln2][1] = 0.0; makeSinCosTables(this.maxN); this.CS_cache[this.ln2] = new double[2][]; this.CS_cache[this.ln2][0] = this.C; this.CS_cache[this.ln2][1] = this.S; makeBitReverseTable(this.maxN); this.bitrev_cache[this.ln2] = this.bitrev; } else { this.C = this.CS_cache[this.ln2][0]; this.S = this.CS_cache[this.ln2][1]; this.bitrev = this.bitrev_cache[this.ln2]; } this.freqMask = this.freqMask_cache[this.ln2]; // for the new one there will be just null this.freqPass = this.freqPass_cache[this.ln2]; this.tempArr = new double[maxN]; } return differentSize; } public double[] getHamming1d() { if (this.maxN >= 0) return getHamming1d(this.maxN); return null; } /** * * @param n power of 2 FHT size. * @return one-dimensional array with Hamming window for low-pass frequency * filtering (maximum in the center) */ public double[] getHamming1d(int n) { return getHamming1d(n, 0.08); // normal Hamming window } public double[] getHamming1d(int n, double min) { // min==0.08 - normal Hamming, 0.0 - pure shifted cosine int ln2; for (ln2 = 0; n > (1 << ln2); ln2++) ; if (n != (1 << ln2)) { String msg = "Not a power of 2 :" + n; IJ.showMessage("Error", msg); throw new IllegalArgumentException(msg); } if (ln2 >= this.hamming1d.length) { String msg = "Too large array length, increase this.maxSize (it is now " + this.maxSize + ", wanted " + ln2 + ")"; IJ.showMessage("Error", msg); throw new IllegalArgumentException(msg); } if ((this.hamming1d[ln2] == null) || (this.hamming1dMin[ln2] != min)) { this.hamming1d[ln2] = new double[n]; this.hamming1dMin[ln2] = min; double C054 = 0.5 * (1 + min); double C046 = 0.5 * (1 - min); for (int i = 0; i <= n / 2; i++) this.hamming1d[ln2][i] = (C054 - C046 * Math.cos((i * 2.0 * Math.PI) / n)); for (int i = 1; i <= n / 2; i++) this.hamming1d[ln2][n - i] = this.hamming1d[ln2][i]; } return this.hamming1d[ln2]; } /** * Creates one-dimensional array for low-pass filtering in the frequency domain, * by multiplying by Gaussian Zero is in the center, same as for Hamming above * * @param lowPass relative frequency. with lowPass==1.0, Gaussian sigma will be * sqrt(2)*Nyquist frequency (==2/n) * @param n - number of FFT points (2^2) * @return one dimensional array, with 1.0 at 0, and minimum in the middle */ public double[] getGaussian1d(double lowPass) { return getGaussian1d(lowPass, this.maxN); } public double[] getGaussian1d(double lowPass, int n) { int ln2; for (ln2 = 0; n > (1 << ln2); ln2++) ; if (n != (1 << ln2)) { String msg = "Not a power of 2 :" + n; IJ.showMessage("Error", msg); throw new IllegalArgumentException(msg); } if (ln2 >= this.gaussian1d.length) { String msg = "Too large array length, increase this.maxSize (it is now " + this.maxSize + ", wanted " + ln2 + ")"; IJ.showMessage("Error", msg); throw new IllegalArgumentException(msg); } if (this.gaussian1d[ln2] == null) { this.gaussian1d[ln2] = new double[n]; this.gaussian1dWidths[ln2] = -1.0; } if (this.gaussian1dWidths[ln2] != lowPass) { this.gaussian1dWidths[ln2] = lowPass; double kLo = (lowPass > 0) ? (1.0 / lowPass / lowPass / n / n) : 0; double mx2; for (int i = 0; i <= n / 2; i++) { mx2 = -i * i; this.gaussian1d[ln2][n / 2 - i] = (kLo > 0.0) ? Math.exp(mx2 * kLo) : 1.0; } for (int i = 1; i <= n / 2; i++) this.gaussian1d[ln2][n - i] = this.gaussian1d[ln2][i]; } return this.gaussian1d[ln2]; } private void makeSinCosTables(int maxN) { int ln2; for (ln2 = 0; maxN > (1 << ln2); ln2++) ; int n = maxN / 4; this.C = new double[n]; this.S = new double[n]; double theta = 0.0; double dTheta = 2.0 * Math.PI / maxN; for (int i = 0; i < n; i++) { this.C[i] = Math.cos(theta); this.S[i] = Math.sin(theta); theta += dTheta; } } private void makeBitReverseTable(int maxN) { bitrev = new int[maxN]; int nLog2 = log2(maxN); for (int i = 0; i < maxN; i++) bitrev[i] = bitRevX(i, nLog2); } /* Performs a 2D FHT (Fast Hartley Transform). */ public void rc2DFHT(double[] x, boolean inverse, int maxN) { // IJ.write("FFT: rc2DFHT (row-column Fast Hartley Transform)"); for (int row = 0; row < maxN; row++) dfht3(x, row * maxN, inverse, maxN); progress(0.4); transposeR(x, maxN); progress(0.5); for (int row = 0; row < maxN; row++) dfht3(x, row * maxN, inverse, maxN); progress(0.7); transposeR(x, maxN); progress(0.8); int mRow, mCol; double A, B, C, D, E; for (int row = 0; row <= maxN / 2; row++) { // Now calculate actual Hartley transform for (int col = 0; col <= maxN / 2; col++) { mRow = (maxN - row) % maxN; mCol = (maxN - col) % maxN; A = x[row * maxN + col]; // see Bracewell, 'Fast 2D Hartley Transf.' IEEE Procs. 9/86 B = x[mRow * maxN + col]; C = x[row * maxN + mCol]; D = x[mRow * maxN + mCol]; E = ((A + D) - (B + C)) / 2; x[row * maxN + col] = A - E; x[mRow * maxN + col] = B + E; x[row * maxN + mCol] = C + E; x[mRow * maxN + mCol] = D - E; } } progress(0.95); } void progress(double percent) { if (showProgress) IJ.showProgress(percent); } /* Performs an optimized 1D FHT. */ public void dfht3(double[] x, int base, boolean inverse, int maxN) { int i, stage, gpNum, gpSize, numGps, Nlog2; int bfNum, numBfs; int Ad0, Ad1, Ad2, Ad3, Ad4, CSAd; double rt1, rt2, rt3, rt4; Nlog2 = log2(maxN); BitRevRArr(x, base, Nlog2, maxN); // bitReverse the input array gpSize = 2; // first & second stages - do radix 4 butterflies once thru numGps = maxN / 4; for (gpNum = 0; gpNum < numGps; gpNum++) { Ad1 = gpNum * 4; Ad2 = Ad1 + 1; Ad3 = Ad1 + gpSize; Ad4 = Ad2 + gpSize; rt1 = x[base + Ad1] + x[base + Ad2]; // a + b rt2 = x[base + Ad1] - x[base + Ad2]; // a - b rt3 = x[base + Ad3] + x[base + Ad4]; // c + d rt4 = x[base + Ad3] - x[base + Ad4]; // c - d x[base + Ad1] = rt1 + rt3; // a + b + (c + d) x[base + Ad2] = rt2 + rt4; // a - b + (c - d) x[base + Ad3] = rt1 - rt3; // a + b - (c + d) x[base + Ad4] = rt2 - rt4; // a - b - (c - d) } if (Nlog2 > 2) { // third + stages computed here gpSize = 4; numBfs = 2; numGps = numGps / 2; // IJ.write("FFT: dfht3 "+Nlog2+" "+numGps+" "+numBfs); for (stage = 2; stage < Nlog2; stage++) { for (gpNum = 0; gpNum < numGps; gpNum++) { Ad0 = gpNum * gpSize * 2; Ad1 = Ad0; // 1st butterfly is different from others - no mults needed Ad2 = Ad1 + gpSize; Ad3 = Ad1 + gpSize / 2; Ad4 = Ad3 + gpSize; rt1 = x[base + Ad1]; x[base + Ad1] = x[base + Ad1] + x[base + Ad2]; x[base + Ad2] = rt1 - x[base + Ad2]; rt1 = x[base + Ad3]; x[base + Ad3] = x[base + Ad3] + x[base + Ad4]; x[base + Ad4] = rt1 - x[base + Ad4]; for (bfNum = 1; bfNum < numBfs; bfNum++) { // subsequent BF's dealt with together Ad1 = bfNum + Ad0; Ad2 = Ad1 + gpSize; Ad3 = gpSize - bfNum + Ad0; Ad4 = Ad3 + gpSize; CSAd = bfNum * numGps; if (((base + Ad2) >= x.length) || ((base + Ad4) >= x.length) || (CSAd >= this.C.length) || (CSAd >= this.S.length)) { System.out.println("dfht3 Error: ln2=" + this.ln2 + " maxN=" + this.maxN + " x.length=" + x.length + " this.C.length=" + this.C.length + " this.S.length=" + this.S.length + " base=" + base + " Ad2=" + Ad2 + " Ad4=" + Ad4 + " CSAd=" + CSAd); } rt1 = x[base + Ad2] * C[CSAd] + x[base + Ad4] * S[CSAd]; rt2 = x[base + Ad4] * C[CSAd] - x[base + Ad2] * S[CSAd]; x[base + Ad2] = x[base + Ad1] - rt1; x[base + Ad1] = x[base + Ad1] + rt1; x[base + Ad4] = x[base + Ad3] + rt2; x[base + Ad3] = x[base + Ad3] - rt2; } /* end bfNum loop */ } /* end gpNum loop */ gpSize *= 2; numBfs *= 2; numGps = numGps / 2; } /* end for all stages */ } /* end if Nlog2 > 2 */ if (inverse) { for (i = 0; i < maxN; i++) x[base + i] = x[base + i] / maxN; } } void transposeR(double[] x, int maxN) { int r, c; double rTemp; for (r = 0; r < maxN; r++) { for (c = r; c < maxN; c++) { if (r != c) { rTemp = x[r * maxN + c]; x[r * maxN + c] = x[c * maxN + r]; x[c * maxN + r] = rTemp; } } } } int log2(int x) { int count = 15; while (!btst(x, count)) count--; return count; } private boolean btst(int x, int bit) { // int mask = 1; return ((x & (1 << bit)) != 0); } void BitRevRArr(double[] x, int base, int bitlen, int maxN) { for (int i = 0; i < maxN; i++) tempArr[i] = x[base + bitrev[i]]; for (int i = 0; i < maxN; i++) x[base + i] = tempArr[i]; } private int bitRevX(int x, int bitlen) { int temp = 0; for (int i = 0; i <= bitlen; i++) if ((x & (1 << i)) != 0) temp |= (1 << (bitlen - i - 1)); return temp & 0x0000ffff; } public int[] nonZeroAmpIndices(double[] ampMask) { int numNonZero = 0; for (int i = 0; i < ampMask.length; i++) if (ampMask[i] > 0.0) numNonZero++; int[] result = new int[numNonZero]; numNonZero = 0; for (int i = 0; i < ampMask.length; i++) if (ampMask[i] > 0.0) result[numNonZero++] = i; return result; } /** * Calculate feature matching quality * * @param h FHT array * @param directionAngle direction of the linear phase gradient * @param distance Distance to the linear feature in the direction * directionAngle (determines the phase gradient * absolute value) * @param distance Distance to the linear feature in the direction * directionAngle (determines the phase gradient * absolute value) * @param phaseIntegrationWidth Width of the band corresponding to the linear * feature * @param amHPhase array with amplitudes and (modified) phases for * each index or null. Used to mask out some of the * points * @param phaseStepCorr Zero-based linear phase values plus these values * approximate real phases * @param bandMask if not null, multiply amplitudes by this mask * @return pair of correlation and self-RMS (to calculate relative feature * strength) */ public double[] linearFeatureStrength(double[] h, double directionAngle, double distance, double phaseIntegrationWidth, // zero - no limit int[] nonZeroIndices, double[][] amHPhase, double[] phaseStepCorr, double[] bandMask) { double halfWidth = (phaseIntegrationWidth > 0.0) ? phaseIntegrationWidth / 2 : this.maxN; double sin = Math.sin(directionAngle); double cos = Math.cos(directionAngle); int halfN = this.maxN / 2; double diffPhase = distance * 2 * Math.PI / this.maxN; double corrSum = 0.0; double sumSquares = 0.0; int numSamples = 0; for (int numPoint = 0; numPoint < nonZeroIndices.length; numPoint++) if ((amHPhase == null) || (amHPhase[numPoint][0] > 0.0)) { int index = nonZeroIndices[numPoint]; int x = (index + halfN) % this.maxN - halfN; int y = index / this.maxN; double t = cos * x + sin * y; if (Math.abs(cos * y - sin * x) <= halfWidth) { int indexMod = ((this.maxN - (index / this.maxN)) % this.maxN) * this.maxN + ((this.maxN - (index % this.maxN)) % this.maxN); double re = 0.5 * (h[index] + h[indexMod]); double im = 0.5 * (h[index] - h[indexMod]); double amp2 = re * re + im * im; if (bandMask != null) amp2 *= bandMask[numPoint]; // sumSquares+=((y!=0)?2:1)*amp2; sumSquares += amp2; double phase = phaseStepCorr[numPoint] + diffPhase * t; double cosPhase = Math.cos(phase); double sinPhase = Math.sin(phase); corrSum += re * cosPhase + im * sinPhase; // debugprint numSamples++; if (this.debug) System.out.println(" Used x=" + x + " y=" + y + " cos*y+sin*x=" + (cos * y + sin * x) + " cos=" + cos + " sin=" + sin + " t=" + t); } else { if (this.debug) System.out.println("Dropped x=" + x + " y=" + y + " cos*y+sin*x=" + (cos * y + sin * x) + " cos=" + cos + " sin=" + sin + " t=" + t); } } double[] result = { corrSum, Math.sqrt(numSamples * sumSquares) }; // second is n*RMS return result; } /** * Reconstruct FHT of the linear feature using the masked amplitudes of the * original and linear phase approximation * * @param h FHT array * @param directionAngle direction of the linear phase gradient * @param distance Distance to the linear feature in the direction * directionAngle (determines the phase gradient * absolute value) * @param phaseIntegrationWidth Width of the band corresponding to the linear * feature * @param nonZeroIndices array of non-zero indices in the FHT array * @param amHPhase array with amplitudes and (modified) phases for * each index or null. Used to mask out some of the * points * @param phaseStepCorr Zero-based linear phase values plus these values * approximate real phases * @param bandMask if not null, multiply amplitudes by this mask * @return FHT Reconstructed array transformed to space domain with quadrants * swapped; */ public double[] reconstructLinearFeature(double[] h, double directionAngle, double distance, double phaseIntegrationWidth, // zero - no limit int[] nonZeroIndices, double[][] amHPhase, double[] phaseStepCorr, double[] bandMask) { double halfWidth = (phaseIntegrationWidth > 0.0) ? phaseIntegrationWidth / 2 : this.maxN; double sin = Math.sin(directionAngle); double cos = Math.cos(directionAngle); int halfN = this.maxN / 2; double diffPhase = distance * 2 * Math.PI / this.maxN; double[] reconstructedFHT = new double[h.length]; for (int i = 0; i < reconstructedFHT.length; i++) reconstructedFHT[i] = 0.0; for (int numPoint = 0; numPoint < nonZeroIndices.length; numPoint++) if ((amHPhase == null) || (amHPhase[numPoint][0] > 0.0)) { int index = nonZeroIndices[numPoint]; int x = (index + halfN) % this.maxN - halfN; int y = index / this.maxN; if (Math.abs(cos * y - sin * x) <= halfWidth) { int indexMod = ((this.maxN - (index / this.maxN)) % this.maxN) * this.maxN + ((this.maxN - (index % this.maxN)) % this.maxN); double re = 0.5 * (h[index] + h[indexMod]); double im = 0.5 * (h[index] - h[indexMod]); double amp = Math.sqrt(re * re + im * im); if (bandMask != null) amp *= bandMask[numPoint]; double t = cos * x + sin * y; double phase = phaseStepCorr[numPoint] + diffPhase * t; double reconstructedRe = amp * Math.cos(phase); double reconstructedIm = amp * Math.sin(phase); reconstructedFHT[index] = reconstructedRe + reconstructedIm; reconstructedFHT[indexMod] = reconstructedRe - reconstructedIm; // if (this.debug) System.out.println("reconstructLinearFeature(): x="+x+" y="+y +" index="+index+ " indexMod="+indexMod+ " re="+re+" im="+im +" amp="+amp+ // " t="+t+ " phase="+ phase+" reconstructedRe="+reconstructedRe+" reconstructedIm="+reconstructedIm); } } transform(reconstructedFHT, true); swapQuadrants(reconstructedFHT); return reconstructedFHT; // space domain } /** * Adds phase shift to stitch symmetrical around 0 areas (there is zero at 0) * * @param directionAngle perpendicular to the linear feature (direction of high * amplitude in frequency domain) * @param zeroPhase extrapolated phase at 0 for positive direction at * angle (half step) * @param distance Estimated distance to line (only used with * zeroBinHalfSize) * @param phaseTolerance - if >0.0, will zero amplitude if the phase is too off * @param amHPhase [numLixel]{amplitude, phase} - phase will be modified! * May be null - will work faster * @param zeroBinHalfSize - if abs(projection) is less than zeroBinHalfSize, use * the best (of 2 options) phase * @param nonZeroIndices * @return applied phase correction, added to the modified phase will result * actual phase */ public double[] compensatePhaseStep(double directionAngle, double zeroPhase, double distance, double phaseTolerance, double[][] amHPhase, int[] nonZeroIndices, double zeroBinHalfSize) { double sin = Math.sin(directionAngle); double cos = Math.cos(directionAngle); int halfN = this.maxN / 2; double diffPhase = distance * 2 * Math.PI / this.maxN; double[] phaseCorr = new double[nonZeroIndices.length]; for (int numPoint = 0; numPoint < nonZeroIndices.length; numPoint++) { int index = nonZeroIndices[numPoint]; int x = (index + halfN) % this.maxN - halfN; int y = index / this.maxN; double t = cos * x + sin * y; double dPhase = (t > 0) ? zeroPhase : (-zeroPhase); if (amHPhase != null) { if ((y != 0) && (Math.abs(t) < zeroBinHalfSize)) { double flattenedPhase = amHPhase[numPoint][1] - diffPhase * t; double diffPlus = flattenedPhase - zeroPhase; double diffMinus = flattenedPhase + zeroPhase; diffPlus -= (2 * Math.PI * Math.round(diffPlus / (2 * Math.PI))); diffMinus -= (2 * Math.PI * Math.round(diffMinus / (2 * Math.PI))); dPhase = (Math.abs(diffPlus) < Math.abs(diffMinus)) ? zeroPhase : (-zeroPhase); /* * System.out.println("numPoint="+numPoint+" x="+x+" y="+y+ * " phase(was)="+amHPhase[numPoint][1]+ " flattenedPhase="+flattenedPhase+ * " diffPlus="+diffPlus+" diffMinus="+diffMinus+ " dPhase="+dPhase); */ } amHPhase[numPoint][1] -= dPhase; amHPhase[numPoint][1] -= (2 * Math.PI * Math.round(amHPhase[numPoint][1] / (2 * Math.PI))); double phaseError = amHPhase[numPoint][1] - diffPhase * t; phaseError -= (2 * Math.PI * Math.round(phaseError / (2 * Math.PI))); if ((phaseTolerance > 0.0) && (Math.abs(phaseError) > phaseTolerance)) { amHPhase[numPoint][0] = 0.0; // bad phase - make zero; } } phaseCorr[numPoint] = dPhase; } return phaseCorr; } public int[] getBandIndices(boolean keepDC, double ribbonWidth, double angle) { int halfN = this.maxN / 2; double halfWidth = ribbonWidth / 2; int length = halfN * this.maxN; int[] tmpIndices = new int[length]; double sin = Math.sin(angle); double cos = Math.cos(angle); int pointNumber = 0; for (int index = 0; index < length; index++) { int x = (index + halfN) % this.maxN - halfN; int y = index / this.maxN; double t = cos * x + sin * y; double d = sin * x - cos * y; int it = (int) Math.round(Math.abs(t)); if ((keepDC || (y > 0) || (x > 0)) && (Math.abs(d) < halfWidth) && (it < halfN)) { tmpIndices[pointNumber++] = index; } } int[] indices = new int[pointNumber]; for (int i = 0; i < pointNumber; i++) indices[i] = tmpIndices[i]; return indices; } public double[] getRibbonMask(boolean keepDC, int[] bandIndices, double ribbonWidth, double resultHighPass0, double angle) { double resultHighPass = keepDC ? 0.0 : resultHighPass0; int halfN = this.maxN / 2; double halfWidth = ribbonWidth / 2; double sin = Math.sin(angle); double cos = Math.cos(angle); if (sin < 0) { sin = -sin; cos = -cos; } double[] ribbonMask = new double[bandIndices.length]; for (int iIndex = 0; iIndex < bandIndices.length; iIndex++) { int index = bandIndices[iIndex]; int x = (index + halfN) % this.maxN - halfN; int y = index / this.maxN; double t = Math.abs(cos * x + sin * y); double d = sin * x - cos * y; // int it=(int) Math.round(Math.abs(t)); int it = (int) Math.round(t); if ((keepDC || (y > 0) || (x > 0)) && (Math.abs(d) < halfWidth) && (it < halfN)) { double w = (0.54 + 0.46 * Math.cos(Math.PI * d / halfWidth)); // (masked) amplitude; if (t < resultHighPass) w *= 0.54 + 0.46 * Math.cos(Math.PI * (resultHighPass - t) / resultHighPass); else if (t > (halfN - halfWidth)) w *= 0.54 + 0.46 * Math.cos(Math.PI * (t - (halfN - halfWidth)) / halfWidth); ribbonMask[iIndex] = w; // if (this.debug) System.out.println("getRibbonMask(): x="+x+" y="+y +" t="+t+" d="+d +" w="+w+ " w0="+ (0.54+0.46*Math.cos(Math.PI*d/halfWidth))); } else { ribbonMask[iIndex] = 0.0; System.out.println("getRibbonMask(): should not get here"); } } return ribbonMask; } /** * Detect linear features in the frequency domain sample * * @param ribbonWidth Width of the band for processing phase (data inside * will be multiplied by the Hamming filter * perpendicular to the direction of the band * @param fht frequency-domain FHT array of the sample * @param minRMSFrac Process only frequencies with amplitude above this * ratio of the RMS of all frequencies. (0 - do not use * this threshold) * @param minAbs Process only frequencies with (absolute) amplitude * above this value (0 - do not use this threshold) * @param maxPhaseMismatch Maximal mismatch between the diagonal pair of 2 pixel * pairs to consider pixels valid (distance between * crossing diagonals) * @param dispertionCost - >0 - use phase dispersion when calculating quality * of phase match, 0 - only weights. * @param filter Bitmask enabling different phase half-steps (+1 - 0, * +2 - pi/2, +4 - pi, +8 - 3pi/2. Value zero allows * arbitrary step step 0 corresponds to thin white line, * pi - thin black line, +pi/2 - edge black-to-white in * the direction of directionAngle, 3pi/2 - * white-to-black in the direction of directionAngle * @return array of the following values: {bestAngle, angle of the perpendicular * towards the linear feature: 0 in the direction of pixel x (right), * pi/2 - in the direction of the pixel Y (down) distance, - distance * from the sample center to the linear feature in the direction * bnestAngle (in pixels) phaseAtZero, - extrapolated phase at zero * frequency (DC) , see description of the parameter "filter" above * maxStrength}; strength of the feature (combines multiple factors) */ public double[] calcPhaseApproximation(double ribbonWidth, double[] fht, double minRMSFrac, double minAbs, double maxPhaseMismatch, double dispertionCost, int filter) { createPolarMasksIndices(ribbonWidth); // all but first time will return immediately. int halfN = this.maxN / 2; double dispersionFatZero = (dispertionCost > 0) ? (Math.PI / halfN / dispertionCost) : 0.0; double[][] dPHdXdY = new double[5][this.maxN * this.maxN / 2]; // amp, phase, phaseWeight, dPhase/dX, dPhase/dY for (int n = 0; n < dPHdXdY.length; n++) for (int i = 0; i < dPHdXdY[0].length; i++) dPHdXdY[n][i] = 0.0; double sum2 = 0.0; for (int index = 0; index < dPHdXdY[0].length; index++) { int indexMod = ((this.maxN - (index / this.maxN)) % this.maxN) * this.maxN + ((this.maxN - (index % this.maxN)) % this.maxN); double re = 0.5 * (fht[index] + fht[indexMod]); double im = 0.5 * (fht[index] - fht[indexMod]); // int x=(index+halfN)%this.maxN-halfN; // int y=index/this.maxN; double amp2 = re * re + im * im; sum2 += amp2; dPHdXdY[0][index] = Math.sqrt(amp2); dPHdXdY[1][index] = Math.atan2(im, re); } double rms = Math.sqrt(sum2 / dPHdXdY[0].length); double thresholdAmplitude = 0.0; if (this.debug) System.out.println("calc_dPHdXdY): rms=" + rms); if (minRMSFrac > 0) thresholdAmplitude = minRMSFrac * rms; if ((minAbs > 0) && ((thresholdAmplitude == 0) || (thresholdAmplitude > minAbs))) thresholdAmplitude = minAbs; for (int index = 0; index < dPHdXdY[0].length; index++) if (dPHdXdY[0][index] < thresholdAmplitude) dPHdXdY[0][index] = 0.0; double oneThird = 1.0 / 3; for (int y = 0; y < (halfN - 1); y++) { for (int x = -halfN + 1; x < (halfN - 1); x++) { // int index00=y*this.maxN+(x+halfN); int index00 = y * this.maxN + (x + this.maxN) % this.maxN; int indexX0 = y * this.maxN + (x + this.maxN + 1) % this.maxN; int index0Y = index00 + this.maxN; int indexXY = (y + 1) * this.maxN + (x + this.maxN + 1) % this.maxN; if ((dPHdXdY[0][index00] > 0) && (dPHdXdY[0][indexX0] > 0) && (dPHdXdY[0][index0Y] > 0) && (dPHdXdY[0][indexXY] > 0)) { double dP00X0 = dPHdXdY[1][indexX0] - dPHdXdY[1][index00]; double dP000Y = dPHdXdY[1][index0Y] - dPHdXdY[1][index00]; double dPX0XY = dPHdXdY[1][indexXY] - dPHdXdY[1][indexX0]; double dP0YXY = dPHdXdY[1][indexXY] - dPHdXdY[1][index0Y]; dP00X0 -= (2 * Math.PI) * Math.round(dP00X0 / (2 * Math.PI)); dP000Y -= (2 * Math.PI) * Math.round(dP000Y / (2 * Math.PI)); dPX0XY -= (2 * Math.PI) * Math.round(dPX0XY / (2 * Math.PI)); dP0YXY -= (2 * Math.PI) * Math.round(dP0YXY / (2 * Math.PI)); if (this.debug) { System.out.println(" index00=" + index00 + " indexX0=" + indexX0 + " index0Y=" + index0Y + " indexXY=" + indexXY); System.out.println(" dPHdXdY[1][index00]=" + dPHdXdY[1][index00] + " dPHdXdY[1][indexX0]=" + dPHdXdY[1][indexX0] + " dPHdXdY[1][index0Y]=" + dPHdXdY[1][index0Y] + " dPHdXdY[1][indexXY]=" + dPHdXdY[1][indexXY]); System.out.println( " dP00X0=" + dP00X0 + " dP000Y=" + dP000Y + " dPX0XY=" + dPX0XY + " dP0YXY=" + dP0YXY); } double curv = Math.abs((dPX0XY - dP000Y) / 2); if ((maxPhaseMismatch == 0) || (curv < maxPhaseMismatch)) { double w = Math.pow(dPHdXdY[0][index00] * dPHdXdY[0][indexX0] * dPHdXdY[0][index0Y], oneThird); if (this.debug) System.out.print(">>> x=" + x + " y=" + y + " curv=" + curv); dPHdXdY[2][index00] += w; dPHdXdY[3][index00] += w * dP00X0; dPHdXdY[4][index00] += w * dP000Y; if (this.debug) System.out.print(" index00=" + index00 + " w=" + w); w = Math.pow(dPHdXdY[0][index00] * dPHdXdY[0][indexX0] * dPHdXdY[0][indexXY], oneThird); dPHdXdY[2][indexX0] += w; dPHdXdY[3][indexX0] += w * dP00X0; dPHdXdY[4][indexX0] += w * dPX0XY; if (this.debug) System.out.print(" indexX0=" + indexX0 + " w=" + w); w = Math.pow(dPHdXdY[0][index00] * dPHdXdY[0][index0Y] * dPHdXdY[0][indexXY], oneThird); dPHdXdY[2][index0Y] += w; dPHdXdY[3][index0Y] += w * dP0YXY; // dP000Y; dPHdXdY[4][index0Y] += w * dP000Y; // dP0YXY; if (this.debug) System.out.print(" index0Y=" + index0Y + " w=" + w); w = Math.pow(dPHdXdY[0][indexX0] * dPHdXdY[0][index0Y] * dPHdXdY[0][indexXY], oneThird); dPHdXdY[2][indexXY] += w; dPHdXdY[3][indexXY] += w * dP0YXY; // dPX0XY; dPHdXdY[4][indexXY] += w * dPX0XY; // dP0YXY; if (this.debug) System.out.print(" indexXY=" + indexXY + " w=" + w); } else { if (this.debug) System.out.println(">-- x=" + x + " y=" + y + " curv=" + curv + " w=" + Math .pow(dPHdXdY[0][index00] * dPHdXdY[0][indexX0] * dPHdXdY[0][index0Y], oneThird)); } } } } for (int i = 1; i < halfN; i++) { int j = this.maxN - i; if ((dPHdXdY[2][i] + dPHdXdY[2][j]) > 0.0) { dPHdXdY[3][i] = (dPHdXdY[3][i] * dPHdXdY[2][i] + dPHdXdY[3][j] * dPHdXdY[2][j]) / (dPHdXdY[2][i] + dPHdXdY[2][j]); dPHdXdY[4][i] = (dPHdXdY[4][i] * dPHdXdY[2][i] + dPHdXdY[4][j] * dPHdXdY[2][j]) / (dPHdXdY[2][i] + dPHdXdY[2][j]); dPHdXdY[2][i] = dPHdXdY[2][i] + dPHdXdY[2][j]; dPHdXdY[3][j] = dPHdXdY[3][i]; dPHdXdY[4][j] = dPHdXdY[4][i]; dPHdXdY[2][j] = dPHdXdY[2][i]; } } for (int index = 0; index < dPHdXdY[0].length; index++) if (dPHdXdY[2][index] > 0) { dPHdXdY[3][index] /= dPHdXdY[2][index]; dPHdXdY[4][index] /= dPHdXdY[2][index]; } // generate show image if (this.debug) { double[][] debugdPHdXdY = new double[5][this.maxN * this.maxN]; for (int i = 0; i < debugdPHdXdY[0].length; i++) { int x = (i % this.maxN) - halfN; int y = (i / this.maxN) - halfN; boolean other = (y < 0); if (y < 0) { y = -y; x = -x; } if ((y < halfN) && (x > -halfN) && (x < halfN)) { int j = (x + this.maxN) % this.maxN + this.maxN * y; for (int n = 0; n < debugdPHdXdY.length; n++) { debugdPHdXdY[n][i] = (other && ((n == 1))) ? (-dPHdXdY[n][j]) : dPHdXdY[n][j]; } } } String[] titles = { "amp", "phase", "weight", "d/dx", "d/dy" }; ShowDoubleFloatArrays.showArrays(debugdPHdXdY, this.maxN, this.maxN, true, "dbgDiffPhases", titles); } int numAngles = this.polarPhaseMasks.length; double[] sinAngle = new double[numAngles]; double[] cosAngle = new double[numAngles]; double[] S0 = new double[numAngles]; double[] SF = new double[numAngles]; double[] SF2 = new double[numAngles]; double[] strength = new double[numAngles]; double[] mean = new double[numAngles]; for (int iAngle = 0; iAngle < numAngles; iAngle++) { double angle = iAngle * Math.PI / numAngles; sinAngle[iAngle] = Math.sin(angle); cosAngle[iAngle] = Math.cos(angle); S0[iAngle] = 0.0; SF[iAngle] = 0.0; SF2[iAngle] = 0.0; strength[iAngle] = 0.0; mean[iAngle] = 0.0; } for (int index = 0; index < this.polarPhaseIndices.length; index++) if (dPHdXdY[2][index] > 0.0) { if (this.debug) { int x = ((index + halfN) % this.maxN) - halfN; int y = (index / this.maxN); System.out.println("\n+++ calc_dPHdXdY()x/y/ |" + x + "|" + y); } for (int iiAngle = 0; iiAngle < this.polarPhaseIndices[index].length; iiAngle++) { int iAngle = this.polarPhaseIndices[index][iiAngle]; double w = this.polarPhaseMasks[iAngle][index] * dPHdXdY[2][index]; S0[iAngle] += w; if (this.debug) { System.out.println(" calc_dPHdXdY() iAngle/mask/weight/w: |" + "|" + iAngle + "|" + this.polarPhaseMasks[iAngle][index] + "|" + dPHdXdY[2][index] + "|" + w + "|" + S0[iAngle]); } double F = dPHdXdY[3][index] * cosAngle[iAngle] + dPHdXdY[4][index] * sinAngle[iAngle]; SF[iAngle] += w * F; SF2[iAngle] += w * F * F; } } int iAngleMax = -1; double maxStrength = 0.0; for (int iAngle = 0; iAngle < numAngles; iAngle++) { double disp = 0; if (S0[iAngle] > 0.0) { mean[iAngle] = SF[iAngle] / S0[iAngle]; disp = Math.sqrt(S0[iAngle] * SF2[iAngle] - SF[iAngle] * SF[iAngle]) / S0[iAngle]; // TODO: reduce dispersion influence if it is small strength[iAngle] = (dispersionFatZero > 0.0) ? (dispersionFatZero * S0[iAngle] / (dispersionFatZero + disp)) : S0[iAngle]; } if (this.debug) System.out.println("calc_dPHdXdY(): |" + iAngle + "|" + S0[iAngle] + "|" + SF[iAngle] + "|" + SF2[iAngle] + "|" + mean[iAngle] + "|" + disp + "|" + strength[iAngle]); if (strength[iAngle] > maxStrength) { maxStrength = strength[iAngle]; iAngleMax = iAngle; } } if (iAngleMax < 0) { if (this.debug) System.out.println("calcPhaseApproximation(): Could not find maximal strength direction"); return null; } double b = 0.5 * (strength[(iAngleMax + 1) % numAngles] - strength[(iAngleMax + numAngles - 1) % numAngles]); double a = 0.5 * (strength[(iAngleMax + 1) % numAngles] + strength[(iAngleMax + numAngles - 1) % numAngles]) - strength[iAngleMax]; // negative // double x=-0.5*b/a; double dAngle = (-0.5 * b / a); double bestAngle = iAngleMax + dAngle; maxStrength += a * dAngle * dAngle + b * dAngle; if (this.debug) System.out.println( "calcPhaseApproximation(): |" + iAngleMax + "|" + dAngle + "|" + bestAngle + "|" + maxStrength); int secondAngle = (iAngleMax + ((dAngle < 0) ? -1 : 1) + numAngles) % numAngles; double secondMean = mean[secondAngle]; if ((dAngle * (secondAngle - iAngleMax)) < 0) secondMean = -secondMean; // roll over double bestMean = mean[iAngleMax] * (1.0 - Math.abs(dAngle)) + secondMean * Math.abs(dAngle); if (bestMean < 0) { bestAngle += numAngles; bestMean = -bestMean; } bestAngle *= Math.PI / numAngles; bestAngle -= (2 * Math.PI) * Math.round(bestAngle / (2 * Math.PI)); double distance = bestMean * this.maxN / (2 * Math.PI); if (this.debug) System.out.println( "calcPhaseApproximation(): bestAngle=" + bestAngle + " (" + IJ.d2s(180 * bestAngle / Math.PI, 2) + "), distance=" + distance + " (bestMean=" + bestMean + "), maxStrength=" + maxStrength); // bin weights, Re and Im for the points with non-zero amplitude and within the // band in the direction angle along the direction angle double[][] wReIm = new double[3][halfN]; // weight/re/im double sin = Math.sin(bestAngle); double cos = Math.cos(bestAngle); double halfWidth = this.polarRibbonWidth / 2; boolean conj = (sin < 0); if (conj) { sin = -sin; cos = -cos; } for (int index = 0; index < dPHdXdY[0].length; index++) { int indexMod = ((this.maxN - (index / this.maxN)) % this.maxN) * this.maxN + ((this.maxN - (index % this.maxN)) % this.maxN); int x = (index + halfN) % this.maxN - halfN; int y = index / this.maxN; double t = cos * x + sin * y; double d = sin * x - cos * y; int it = (int) Math.round(Math.abs(t)); if (((y > 0) || (x > 0)) && (Math.abs(d) < halfWidth) && (it < halfN)) { double re = 0.5 * (fht[index] + fht[indexMod]); double im = 0.5 * (fht[index] - fht[indexMod]); if (conj) im = -im; if (t < 0) { im = -im; t = -t; } double w = dPHdXdY[0][index] * (0.54 + 0.46 * Math.cos(Math.PI * d / halfWidth)); // (masked) amplitude; if (t < halfWidth) w *= 0.54 + 0.46 * Math.cos(Math.PI * (halfWidth - t) / halfWidth); else if (t > (halfN - halfWidth)) w *= 0.54 + 0.46 * Math.cos(Math.PI * (t - (halfN - halfWidth)) / halfWidth); wReIm[0][it] += w; wReIm[1][it] += w * re; wReIm[2][it] += w * im; } } if (this.debug) { System.out.println("calcPhaseApproximation(): binnig w/re/im along selected direction"); for (int it = 0; it < wReIm[0].length; it++) { System.out.println("calcPhaseApproximation(): |" + it + "|" + wReIm[0][it] + "|" + wReIm[1][it] + "|" + wReIm[2][it] + "|" + Math.atan2(wReIm[2][it], wReIm[1][it])); } } double[] approximatedZeroPhaseDist = approximateLinearPhase(wReIm, filter); if (approximatedZeroPhaseDist == null) return null; if (this.debug) { System.out.println("calcPhaseApproximation(): initial distance=" + distance + " updated distance=" + approximatedZeroPhaseDist[1] + " phase at zero=" + approximatedZeroPhaseDist[0] + " (" + IJ.d2s(approximatedZeroPhaseDist[0] * 180.0 / Math.PI, 1) + ")"); for (int it = 0; it < wReIm[0].length; it++) { System.out.println("calcPhaseApproximation(): |" + it + "|" + wReIm[0][it] + "|" + wReIm[1][it] + "|" + wReIm[2][it] + "|" + Math.atan2(wReIm[2][it], wReIm[1][it])); } } distance = approximatedZeroPhaseDist[1]; double phaseAtZero = approximatedZeroPhaseDist[0]; double[] result = { bestAngle, distance, phaseAtZero, maxStrength }; return result; } private double polarRibbonWidth = -1; private double[][] polarPhaseMasks = null; private int[][] polarPhaseIndices = null; public void createPolarMasksIndices(double ribbonWidth) { // <0 - invalidate if (this.polarRibbonWidth == ribbonWidth) return; // already set; this.polarRibbonWidth = ribbonWidth; if (ribbonWidth < 0.0) { this.polarPhaseMasks = null; this.polarPhaseIndices = null; return; } int halfN = this.maxN / 2; int length = halfN * this.maxN; double halfWidth = ribbonWidth / 2; int numAngles = (int) Math.round(Math.PI * halfN); this.polarPhaseMasks = new double[numAngles][length]; this.polarPhaseIndices = new int[length][]; for (int i = 0; i < this.polarPhaseIndices.length; i++) { this.polarPhaseIndices[i] = null; for (int iAngle = 0; iAngle < numAngles; iAngle++) this.polarPhaseMasks[iAngle][i] = 0.0; } for (int iAngle = 0; iAngle < numAngles; iAngle++) { double angle = iAngle * Math.PI / numAngles; double sin = Math.sin(angle); double cos = Math.cos(angle); double sumWeights = 0.0; for (int index = 0; index < length; index++) { int x = (index + halfN) % this.maxN - halfN; int y = index / this.maxN; double t = cos * x + sin * y; double d = sin * x - cos * y; if (((y > 0) || (x > 0)) && (Math.abs(d) < halfWidth) && (Math.abs(t) < halfN)) { this.polarPhaseMasks[iAngle][index] = 0.54 + 0.46 * Math.cos(d / halfWidth * Math.PI); // normalize if (Math.abs(t) > (halfN - halfWidth)) this.polarPhaseMasks[iAngle][index] *= 0.54 + 0.46 * Math.cos((Math.abs(t) - (halfN - halfWidth)) / halfWidth * Math.PI); sumWeights += this.polarPhaseMasks[iAngle][index]; } } sumWeights = 1.0 / sumWeights; for (int index = 0; index < length; index++) this.polarPhaseMasks[iAngle][index] *= sumWeights; } for (int index = 0; index < length; index++) { int numUsedAngles = 0; for (int iAngle = 0; iAngle < numAngles; iAngle++) if (this.polarPhaseMasks[iAngle][index] > 0.0) numUsedAngles++; this.polarPhaseIndices[index] = new int[numUsedAngles]; // may be 0 int iiAngle = 0; for (int iAngle = 0; iAngle < numAngles; iAngle++) if (this.polarPhaseMasks[iAngle][index] > 0.0) this.polarPhaseIndices[index][iiAngle++] = iAngle; } if (this.debug) { double[][] debugPolarPhaseMasks = new double[numAngles + 1][this.maxN * this.maxN]; for (int i = 0; i < debugPolarPhaseMasks[0].length; i++) { int x = (i % this.maxN) - halfN; int y = (i / this.maxN) - halfN; if (y < 0) { y = -y; x = -x; } if ((y < halfN) && (x > -halfN) && (x < halfN)) { int j = (x + this.maxN) % this.maxN + this.maxN * y; for (int n = 0; n < numAngles; n++) { debugPolarPhaseMasks[n][i] = this.polarPhaseMasks[n][j]; } debugPolarPhaseMasks[numAngles][i] = this.polarPhaseIndices[j].length; } } String[] titles = new String[numAngles + 1]; for (int n = 0; n < numAngles; n++) titles[n] = IJ.d2s(180.0 * n / numAngles, 1); titles[numAngles] = "usage number"; ShowDoubleFloatArrays.showArrays(debugPolarPhaseMasks, this.maxN, this.maxN, true, "dbgDiffPhases", titles); } } /** * Project and accumulate Re, Im on the direction perpendicular to the linear * feature * * @param directionAngle direction of the perpendicular line * @param phaseIntegrationWidth Disregard pixels outside of the band centered * along directionAngle * @param h FHT array * @param nonZeroIndices List of the above-threshold indices in the FHT * array * @return arrays of {weight, re, Im} along the line with the step equal to 1 * pixel */ public double[][] binReIm(double directionAngle, double phaseIntegrationWidth, // zero - no limit double[] h, int[] nonZeroIndices) { int halfN = this.maxN / 2; double[][] wReIm = new double[3][halfN]; // weight/re/im for (int n = 0; n < wReIm.length; n++) for (int i = 0; i < wReIm[0].length; i++) wReIm[n][i] = 0.0; double sin = Math.sin(directionAngle); double cos = Math.cos(directionAngle); double amp, re, im; double halfWidth = (phaseIntegrationWidth > 0.0) ? phaseIntegrationWidth / 2 : this.maxN; for (int numPoint = 0; numPoint < nonZeroIndices.length; numPoint++) { int index = nonZeroIndices[numPoint]; int indexMod = ((this.maxN - (index / this.maxN)) % this.maxN) * this.maxN + ((this.maxN - (index % this.maxN)) % this.maxN); re = 0.5 * (h[index] + h[indexMod]); im = 0.5 * (h[index] - h[indexMod]); int x = (index + halfN) % this.maxN - halfN; int y = index / this.maxN; amp = Math.sqrt(re * re + im * im); if (y == 0) amp *= 0.5; double t = cos * x + sin * y; if (Math.abs(cos * y - sin * x) <= halfWidth) { int itl = (int) Math.floor(t); int ith = itl + 1; double wl = amp * (t - itl); double wh = amp * (itl + 1 - t); boolean conj = (itl < 0); if (conj) { itl = -itl; ith = -ith; } if (itl < wReIm[0].length) { wReIm[0][itl] += wl; wReIm[1][itl] += wl * re; if (conj) wReIm[2][itl] -= wl * im; else wReIm[2][itl] += wl * im; } if (ith < wReIm[0].length) { wReIm[0][ith] += wh; wReIm[1][ith] += wh * re; if (conj) wReIm[2][ith] -= wh * im; else wReIm[2][ith] += wh * im; } } else { System.out.println( "binReIm(): dropped point x=" + x + " y=" + y + " as it is too far from the center line"); } } for (int i = 0; i < wReIm[0].length; i++) { if (wReIm[0][i] > 0.0) { wReIm[1][i] /= wReIm[0][i]; wReIm[2][i] /= wReIm[0][i]; } } return wReIm; } /** * Linear approximate phase, optionally filtering types of features (black line, * white line black-to-white edge, white-to-black edge * * @param wReIm - array of weights, Re and Im along the selected direction * @param filter Bitmask enabling different phase half-steps (+1 - 0, +2 - pi/2, * +4 - pi, +8 - 3pi/2. Value zero allows arbitrary step step 0 * corresponds to thin white line, pi - thin black line, +pi/2 - * edge black-to-white in the direction of directionAngle, 3pi/2 - * white-to-black in the direction of directionAngle * @return pair of phase at 0 and slope expressed in distance (in pixels) from * the center to the linear feature along the selected direction (may be * negative) */ public double[] approximateLinearPhase(double[][] wReIm, int filter) { double[] phase = new double[wReIm[0].length]; for (int i = 0; i < phase.length; i++) if (wReIm[0][i] > 0.0) phase[i] = Math.atan2(wReIm[2][i], wReIm[1][i]); double diffPhase = 0.0; double sumWeights = 0.0; for (int i = 0; i < (phase.length - 1); i++) { double w = wReIm[0][i] * wReIm[0][i + 1]; if (w > 0.0) { double diff = phase[i + 1] - phase[i]; if (diff > Math.PI) while (diff > Math.PI) diff -= 2 * Math.PI; else if (diff < -Math.PI) while (diff < -Math.PI) diff += 2 * Math.PI; sumWeights += w; diffPhase += diff * w; } } if (sumWeights == 0.0) return null; diffPhase /= sumWeights; double[] binPhases = { 0.0, Math.PI / 2, Math.PI, -Math.PI / 2 }; double[] zeroPhaseBins = { 0.0, 0.0, 0.0, 0.0 }; double zeroPhase = Double.NaN; switch (filter) { case 0: double re = 0.0, im = 0.0; for (int n = 0; n < phase.length; n++) if (wReIm[0][n] > 0.0) { double phZero = phase[n] - diffPhase * n; re += wReIm[0][n] * Math.cos(phZero); im += wReIm[0][n] * Math.sin(phZero); } zeroPhase = Math.atan2(im, re); if (Double.isNaN(zeroPhase)) zeroPhase = 0.0; // ?? break; case 1: zeroPhase = binPhases[0]; break; case 2: zeroPhase = binPhases[1]; break; case 4: zeroPhase = binPhases[2]; break; case 8: zeroPhase = binPhases[3]; break; default: // double [] zeroPhaseBins={0.0,0.0,0.0,0.0}; boolean[] binEnabled = { (filter & 1) != 0, (filter & 2) != 0, (filter & 4) != 0, (filter & 8) != 0 }; for (int n = 0; n < phase.length; n++) if (wReIm[0][n] > 0.0) { double phZero = phase[n] - diffPhase * n; phZero -= 2 * Math.PI * Math.floor(phZero / (2 * Math.PI)); // now in the range 0..2*PI double bestBinDif = 2 * Math.PI; int iBestBin = 0; for (int i = 0; i < binEnabled.length; i++) if (binEnabled[i]) { double d = phZero - binPhases[i]; d -= 2 * Math.PI * Math.round(d / (2 * Math.PI)); d = Math.abs(d); if (d < bestBinDif) { bestBinDif = d; iBestBin = i; } } zeroPhaseBins[iBestBin] += wReIm[0][n]; } double maxBinValue = 0; int iMaxBin = -1; for (int i = 0; i < zeroPhaseBins.length; i++) if (zeroPhaseBins[i] > maxBinValue) { maxBinValue = zeroPhaseBins[i]; iMaxBin = i; } zeroPhase = iMaxBin * 2 * Math.PI / zeroPhaseBins.length; if (zeroPhase > Math.PI) zeroPhase -= 2 * Math.PI; } double[] fullPhase = phase.clone(); for (int i = 0; i < fullPhase.length; i++) if (wReIm[0][i] > 0.0) { double diff = fullPhase[i] - (zeroPhase + diffPhase * i); fullPhase[i] -= 2 * Math.PI * Math.round(diff / (2 * Math.PI)); } if (this.debug) { for (int i = 0; i < zeroPhaseBins.length; i++) System.out.println("Phase " + i + " bin: " + zeroPhaseBins[i]); System.out.println("Initial estimated distance to line " + (this.maxN * diffPhase / (2 * Math.PI)) + ", zeroPhase=" + zeroPhase); for (int i = 0; i < phase.length; i++) if (wReIm[0][i] > 0.0) { System.out.println(i + "|" + IJ.d2s(phase[i], 3) + "|" + IJ.d2s(fullPhase[i], 3) + "|" + IJ.d2s((zeroPhase + diffPhase * i), 3)); } } // re-adjust slope (distance) after zero phase is selected double SFX = 0.0, SX2 = 0.0; for (int i = 0; i < fullPhase.length; i++) if (wReIm[0][i] > 0.0) { SFX += wReIm[0][i] * i * (fullPhase[i] - zeroPhase); SX2 += wReIm[0][i] * i * i; } diffPhase = SFX / SX2; double distance = this.maxN * diffPhase / (2 * Math.PI); double[] result = { zeroPhase, distance }; if (this.debug) { System.out.println("Recalculated distance to line " + distance + " (slope=" + diffPhase + " rad/pix)"); } return result; } public double[][] calcIndAmReIm(double[] h, int[] nonZeroIndices) { double[][] ampReIm = new double[3][nonZeroIndices.length]; for (int numPoint = 0; numPoint < nonZeroIndices.length; numPoint++) { int index = nonZeroIndices[numPoint]; int indexMod = ((this.maxN - (index / this.maxN)) % this.maxN) * this.maxN + ((this.maxN - (index % this.maxN)) % this.maxN); ampReIm[1][numPoint] = 0.5 * (h[index] + h[indexMod]); ampReIm[2][numPoint] = 0.5 * (h[index] - h[indexMod]); ampReIm[0][numPoint] = Math .sqrt(ampReIm[1][index] * ampReIm[1][index] + ampReIm[2][index] * ampReIm[2][index]); } return ampReIm; } public double[][] calcIndAmHPhase(double[] h, int[] nonZeroIndices) { double re, im; double[][] amHPhase = new double[nonZeroIndices.length][2]; for (int numPoint = 0; numPoint < nonZeroIndices.length; numPoint++) { int index = nonZeroIndices[numPoint]; int indexMod = ((this.maxN - (index / this.maxN)) % this.maxN) * this.maxN + ((this.maxN - (index % this.maxN)) % this.maxN); re = 0.5 * (h[index] + h[indexMod]); im = 0.5 * (h[index] - h[indexMod]); amHPhase[numPoint][0] = Math.sqrt(re * re + im * im); if ((index / this.maxN) > 0) amHPhase[numPoint][0] *= 2; // double weight amHPhase[numPoint][1] = Math.atan2(im, re); /* * System.out.println("numPoint="+numPoint+" index="+nonZeroIndices[numPoint]+ * " x="+(nonZeroIndices[numPoint]%this.maxN)+" y="+(nonZeroIndices[numPoint]/ * this.maxN)+ " re="+re+" im="+im+" phase="+amHPhase[numPoint][1]); */ } return amHPhase; } public int updateFullCycles(int[] fullCycles, // may be modified double px, double py, int[] nonZeroIndices, double[][] amHPhase) { int numChanges = 0; double e = 1E-8; double[] linearPhase = tiltedPhase(px, py, nonZeroIndices); for (int numPoint = 0; numPoint < nonZeroIndices.length; numPoint++) { // double delta=amHPhase[numPoint][1]+fullCycles[numPoint]*2*Math.PI -linearPhase[numPoint]; double delta = linearPhase[numPoint] - (amHPhase[numPoint][1] + fullCycles[numPoint] * 2 * Math.PI); if (Math.abs(delta) > (Math.PI + e)) { System.out.print("numPoint=" + numPoint + " index=" + nonZeroIndices[numPoint] + " x=" + (nonZeroIndices[numPoint] % this.maxN) + " y=" + (nonZeroIndices[numPoint] / this.maxN) + " fullCycles[" + numPoint + "]: " + fullCycles[numPoint] + " -> "); fullCycles[numPoint] += (int) Math.round(delta / (2 * Math.PI)); System.out.println(fullCycles[numPoint]); numChanges++; } } return numChanges; } public double calcWeightedPhaseRMS(int[] fullCycles, // may be modified double px, double py, int[] nonZeroIndices, double[][] ampPhase) { double[] linearPhase = tiltedPhase(px, py, nonZeroIndices); double S0 = 0.0, S2 = 0.0; for (int numPoint = 0; numPoint < ampPhase.length; numPoint++) { double delta = ampPhase[numPoint][1] + fullCycles[numPoint] * 2 * Math.PI - linearPhase[numPoint]; S0 += ampPhase[numPoint][0]; S2 += ampPhase[numPoint][0] * delta * delta; /* * System.out.println("nonZeroIndices["+numPoint+"]="+nonZeroIndices[numPoint]+ * " x="+(nonZeroIndices[numPoint]%this.maxN)+" y="+(nonZeroIndices[numPoint]/ * this.maxN)+ * " ampPhase["+numPoint+"][0]="+ampPhase[numPoint][0]+" ampPhase["+numPoint+ * "][1]="+ampPhase[numPoint][1]+ " delta="+delta+" S0="+S0+" S2="+S2); */ } double weightedRMS = Math.sqrt(S2 / S0); return weightedRMS; } public double[] calcPhaseShift(int[] fullCycles, // may be modified int[] nonZeroIndices, double[][] ampPhase) { // double [] linearPhase= tiltedPhase (px, py, nonZeroIndices); // minimize weighted RMS with a plane through (0,0) int halfN = this.maxN / 2; double SFX = 0.0, SFY = 0.0, SX2 = 0.0, SY2 = 0.0, SXY = 0.0; for (int numPoint = 0; numPoint < ampPhase.length; numPoint++) { int index = nonZeroIndices[numPoint]; int x = (index + halfN) % this.maxN - halfN; int y = index / this.maxN; double fXY = ampPhase[numPoint][1] + Math.PI * 2 * fullCycles[numPoint]; double w = ampPhase[numPoint][0]; SFX += w * x * fXY; SFY += w * y * fXY; SX2 += w * x * x; SY2 += w * y * y; SXY += w * x * y; } double denominator = SXY * SXY - SX2 * SY2; double[] phaseTilt = { (SFY * SXY - SFX * SY2) / denominator, (SFX * SXY - SFY * SX2) / denominator }; return phaseTilt; } public double calcPhaseShift(double px, double py, int[] fullCycles, // may be modified int[] nonZeroIndices, double[][] ampPhase) { // double [] linearPhase= tiltedPhase (px, py, nonZeroIndices); // minimize weighted RMS with a plane through (0,0) int halfN = this.maxN / 2; double SFX = 0.0, SFY = 0.0, SX2 = 0.0, SY2 = 0.0, SXY = 0.0; for (int numPoint = 0; numPoint < ampPhase.length; numPoint++) { int index = nonZeroIndices[numPoint]; int x = (index + halfN) % this.maxN - halfN; int y = index / this.maxN; double fXY = ampPhase[numPoint][1] + Math.PI * 2 * fullCycles[numPoint]; double w = ampPhase[numPoint][0]; SFX += w * x * fXY; SFY += w * y * fXY; SX2 += w * x * x; SY2 += w * y * y; SXY += w * x * y; } return (py * SFY + px * SFX - px * py * SXY) / (px * px * SX2 + py * py * SY2); } public double[] tiltedPhase(double px, double py, int[] nonZeroIndices) { int halfN = this.maxN / 2; double[] phase = new double[nonZeroIndices.length]; for (int numPoint = 0; numPoint < nonZeroIndices.length; numPoint++) { int index = nonZeroIndices[numPoint]; int x = (index + halfN) % this.maxN - halfN; int y = index / this.maxN; phase[numPoint] = px * x + py * y; } return phase; } /* * public double [][] calcAmReIm (double [] h, double [] ampMask){ int rowMod, * colMod; int halfN=this.maxN/2; double [][] ampReIm=new * double[3][this.maxN*halfN]; for (int r =0; r<this.maxN/2; r++) { rowMod = * (this.maxN - r) % this.maxN; for (int c=0; c<this.maxN; c++){ int index= r * * this.maxN + c; if ((ampMask==null) || (ampMask[index]>0.0)){ colMod = * (this.maxN - c) % this.maxN; int indexMod=rowMod * this.maxN + colMod; * ampReIm[1][index]=0.5*(h[index]+h[indexMod]); * ampReIm[2][index]=0.5*(h[index]-h[indexMod]); * ampReIm[0][index]=Math.sqrt(ampReIm[1][index]*ampReIm[1][index]+ampReIm[2][ * index]*ampReIm[2][index]); } else { ampReIm[1][index]=0.0; * ampReIm[2][index]=0.0; ampReIm[0][index]=0.0; } } } return ampReIm; } */ public double[] filterAmplitude(double[] h, boolean removeIslands, double threshold) { // relative to RMS value int rowMod, colMod; int halfN = this.maxN / 2; double[] amplitude = new double[this.maxN * halfN]; double max = 0; int iMax = -1; double sum2 = 0.0; double amp2; for (int r = 0; r < this.maxN / 2; r++) { rowMod = (this.maxN - r) % this.maxN; for (int c = 0; c < this.maxN; c++) { int index = r * this.maxN + c; if (c == this.maxN / 2) { amplitude[index] = 0.0; } else { colMod = (this.maxN - c) % this.maxN; int indexMod = rowMod * this.maxN + colMod; if (r > 0) { amp2 = h[index] * h[index] + h[indexMod] * h[indexMod]; // count twice (for the second // symmetrical half) amplitude[index] = 2.0 * amp2; // so square root will be twice } else { amp2 = 0.5 * (h[index] * h[index] + h[indexMod] * h[indexMod]); amplitude[index] = amp2; } sum2 += amp2; if (amp2 > max) { max = amp2; iMax = index; } } } } double minAmp2 = threshold * sum2 / (this.maxN - 1) / (this.maxN - 1); for (int i = 0; i < amplitude.length; i++) { if (amplitude[i] < minAmp2) amplitude[i] = 0.0; else amplitude[i] = Math.sqrt(amplitude[i]); } int numDefined = 0; if (this.debug) { for (int i = 0; i < amplitude.length; i++) if (amplitude[i] > 0.0) numDefined++; System.out.println( "quadraticAmplitude(): number of defined cells=" + numDefined + " ( of " + amplitude.length + ")"); } if (removeIslands) { boolean[] confirmed = new boolean[amplitude.length]; for (int i = 0; i < confirmed.length; i++) confirmed[i] = false; int[][] dirs = { { 1, 0 }, { 1, 1 }, { 0, 1 }, { -1, 1 }, { -1, 0 }, { -1, -1 }, { 0, -1 }, { 1, -1 } }; List<Integer> ampList = new ArrayList<Integer>(confirmed.length); int x = (iMax + halfN) % this.maxN - halfN; int y = iMax / this.maxN; Integer Index = iMax; ampList.add(Index); confirmed[Index] = true; while (ampList.size() > 0) { int index = ampList.remove(0); x = (index + halfN) % this.maxN - halfN; y = index / this.maxN; for (int iDir = 0; iDir < dirs.length; iDir++) { int x1 = x + dirs[iDir][0]; int y1 = y + dirs[iDir][1]; if (y1 < 0) { y1 = -y1; x1 = -x1; } if ((x1 > -halfN) && (x1 < halfN) && (y1 >= 0) && (y1 < halfN)) { Index = (x1 + this.maxN) % this.maxN + this.maxN * y1; if (!confirmed[Index] && (amplitude[Index] > 0.0)) { confirmed[Index] = true; ampList.add(Index); } } } } for (int i = 0; i < amplitude.length; i++) { if (!confirmed[i]) amplitude[i] = 0.0; } if (this.debug) { numDefined = 0; for (int i = 0; i < amplitude.length; i++) if (amplitude[i] > 0.0) numDefined++; System.out.println("quadraticAmplitude(): number of remaining cells=" + numDefined + " ( of " + amplitude.length + ")"); } } /* * if (this.debug){ ShowDoubleFloatArrays.showArrays( amplitude, this.maxN, * halfN, "fltAmp"); double [] debugAmplitude=new double [this.maxN*this.maxN]; * for (int i=0;i<debugAmplitude.length;i++){ debugAmplitude[i]=0.0; int x= * (i%this.maxN)-halfN; int y= (i/this.maxN)-halfN; if (y<0) { y=-y; x=-x; } if * ((y<halfN) && (x>-halfN) && ( x<halfN)){ int j= (x+this.maxN)%this.maxN + * this.maxN*y; debugAmplitude[i]=amplitude[j]; if (y>0) debugAmplitude[i]*=0.5; * } * * } ShowDoubleFloatArrays.showArrays( debugAmplitude, this.maxN, this.maxN, * "fltAmp"); } */ return amplitude; } public double[] quadraticAmplitude(double[] amplitude, boolean removeIslands, double threshold) { // relative to RMS // value int halfN = this.maxN / 2; double[] result = new double[3]; // CX2, CY2, CXY for (int i = 0; i < result.length; i++) result[i] = 0.0; double s = 0.0; double amp; for (int r = 0; r < this.maxN / 2; r++) { for (int c = 0; c < this.maxN; c++) if (c != this.maxN / 2) { int index = r * this.maxN + c; amp = amplitude[index]; s += amp; int x = ((c + halfN) % this.maxN) - halfN; result[0] += amp * x * x; result[1] += amp * r * r; result[2] += amp * x * r; } } s *= halfN * halfN; // normalize so result will be related to pixels, not to the FFT size for (int i = 0; i < result.length; i++) result[i] /= s; return result; } /* * public double [] quadraticAmplitude( double [] h){ int rowMod, colMod; * double[] result = new double[3]; // CX2, CY2, CXY for (int * i=0;i<result.length;i++) result[i]=0.0; double s=0.0; double amp; int * halfN=this.maxN/2; for (int r =0; r<this.maxN/2; r++) { rowMod = (this.maxN - * r) % this.maxN; for (int c=0; c<this.maxN; c++) if (c!=this.maxN/2){ colMod = * (this.maxN - c) % this.maxN; int index= r * this.maxN + c; int * indexMod=rowMod * this.maxN + colMod; * amp=Math.sqrt(0.5*(h[index]*h[index]+h[indexMod]*h[indexMod])); if (r>0) * amp*=2; // count twice (for the second symmetrical half) s+=amp; int * x=((c+halfN)%this.maxN)-halfN; result[0]+=amp*x*x; result[1]+=amp*r*r; * result[2]+=amp*x*r; } } s*=halfN*halfN; // normalize so result will be * related to pixels, not to the FFT size for (int i=0;i<result.length;i++) * result[i]/=s; return result; } */ public double[] ellipseAmplitude(double[] quadratic) { double[] result = new double[3]; // angle (clockwise from X, large half-diameter, small/large ratio double alpha = 0.5 * Math.atan2(2 * quadratic[2], quadratic[0] - quadratic[1]); // 1/2*atan(2*Sxy/(sx2-sy2)) double cos = Math.cos(alpha), sin = Math.sin(alpha); double I1 = sin * sin * quadratic[0] + cos * cos * quadratic[1] - 2 * cos * sin * quadratic[2]; double I2 = cos * cos * quadratic[0] + sin * sin * quadratic[1] + 2 * cos * sin * quadratic[2]; boolean rot = (I1 > I2); result[0] = alpha + (rot ? Math.PI : 0); result[1] = Math.sqrt(rot ? I1 : I2); result[2] = Math.sqrt(rot ? (I2 / I1) : (I1 / I2)); return result; } public double[] multiply(double[] h1, double[] h2, boolean conjugate) { int rowMod, colMod; double h2e, h2o; double[] product = 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; h2e = (h2[r * maxN + c] + h2[rowMod * maxN + colMod]) / 2; h2o = (h2[r * maxN + c] - h2[rowMod * maxN + colMod]) / 2; if (conjugate) product[r * maxN + c] = h1[r * maxN + c] * h2e - h1[rowMod * maxN + colMod] * h2o; else product[r * maxN + c] = h1[r * maxN + c] * h2e + h1[rowMod * maxN + colMod] * h2o; } } return product; } public double[] phaseMultiply(double[] h1, double[] h2, double phaseCoeff) { int rowMod, colMod; double h2e, h2o, d; double[] product = 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; h2e = (h2[r * maxN + c] + h2[rowMod * maxN + colMod]) / 2; h2o = (h2[r * maxN + c] - h2[rowMod * maxN + colMod]) / 2; d = phaseCoeff * (h2e * h2e + h2o * h2o) + (1.0 - phaseCoeff); product[r * maxN + c] = (h1[r * maxN + c] * h2e - h1[rowMod * maxN + colMod] * h2o) / d; } } return product; } public double[] phaseMultiplyNorm(double[] h1, double[] h2, double phaseCoeff) { int size = maxN; int size2 = size >> 1; int size21 = size * size2; int rowMod, colMod, base, baseMod; double h2e, h2o; double[] product = new double[size * size]; for (int r = 0; r < size; r++) { rowMod = (size - r) % size; for (int c = 0; c < size; c++) { colMod = (size - c) % size; h2e = (h2[r * size + c] + h2[rowMod * size + colMod]) / 2; h2o = (h2[r * size + c] - h2[rowMod * size + colMod]) / 2; product[r * size + c] = (h1[r * size + c] * h2e - h1[rowMod * size + colMod] * h2o); } } double[] amplitude = calculateAmplitudeHalf(product); double avg_ampl = 0.0; for (int row = 1; row < size2; row++) { base = row * size; for (int col = 0; col < size; col++) { avg_ampl += amplitude[col + base]; } } avg_ampl *= 2; for (int col = 0; col < size; col++) { avg_ampl += amplitude[col] + amplitude[col + size21]; } avg_ampl /= size * size; double aoffs = (1.0 - phaseCoeff) * avg_ampl; double ampl; for (int row = 0; row <= size2; row++) { base = row * size; rowMod = (size - row) % size; baseMod = rowMod * size; for (int col = 0; col < size; col++) { ampl = phaseCoeff * amplitude[col + base] + aoffs; product[col + base] /= ampl; if ((row > 0) && (row < size2)) { colMod = (size - col) % size; product[colMod + baseMod] /= ampl; } } } return product; } // Multiply by real array (i.e. filtering in frequency domain). Array m length should be at least maxN*maxN/2+1 public void multiplyByReal(double[] h, double[] m) { int rowMod, colMod; int index = 0, indexM; for (int r = 0; r < maxN; r++) { rowMod = (maxN - r) % maxN; for (int c = 0; c < maxN; c++) { colMod = (maxN - c) % maxN; indexM = rowMod * maxN + colMod; if (indexM > index) indexM = index; h[index] *= m[indexM]; index++; } } } /* nothing to prevent division by small values */ /* * Returns the image resulting from the point by point Hartley division of this * image by the specified image. Both images are assumed to be in the frequency * domain. Division in the frequency domain is equivalent to deconvolution in * the space domain. */ public double[] divide(double[] h1, double[] h2) { int rowMod, colMod; double mag, h2e, h2o; 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[r * maxN + c] * h2[r * maxN + c] + h2[rowMod * maxN + colMod] * h2[rowMod * maxN + colMod]; if (mag < 1e-20) mag = 1e-20; h2e = (h2[r * maxN + c] + h2[rowMod * maxN + colMod]); h2o = (h2[r * maxN + c] - h2[rowMod * maxN + colMod]); double tmp = (h1[r * maxN + c] * h2e - h1[rowMod * maxN + colMod] * h2o); result[r * maxN + c] = tmp / mag; } } 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[r * maxN + c] * h2[r * maxN + c] + h2[rowMod * maxN + colMod] * h2[rowMod * maxN + colMod] + fz2; if (mag < 1e-20) mag = 1e-20; h2e = (h2[r * maxN + c] + h2[rowMod * maxN + colMod]); h2o = (h2[r * maxN + c] - h2[rowMod * maxN + colMod]); double tmp = (h1[r * maxN + c] * h2e - h1[rowMod * maxN + 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[r * maxN + c] = amp[r * maxN + c]; result[rowMod * maxN + colMod] = amp[r * maxN + c]; } } return result; } public double[] calculateAmplitude(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); } swapQuadrants(amp); return amp; } 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++) { 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) { int size = (int) Math.sqrt(fht.length); int size2 = (size >> 1) + 1; double[] amp = new double[size * size2]; for (int row = 0; row < size2; row++) { amplitude(row, size, fht, amp); } return amp; } public double[] calculateAmplitude2(double[] fht) { int size = (int) Math.sqrt(fht.length); double[] amp = new double[size * size]; for (int row = 0; row < size; row++) { amplitude2(row, size, fht, amp); } swapQuadrants(amp); return amp; } /* Amplitude of one row from 2D Hartley Transform. */ static void amplitude(int row, int size, double[] fht, double[] amplitude) { int base = row * size; int l; for (int c = 0; c < size; c++) { l = ((size - row) % size) * size + (size - c) % size; amplitude[base + c] = Math.sqrt(fht[base + c] * fht[base + c] + fht[l] * fht[l]); } } 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. */ static void amplitude2(int row, int size, double[] fht, double[] amplitude) { int base = row * size; int l; for (int c = 0; c < size; c++) { l = ((size - row) % size) * size + (size - c) % size; amplitude[base + c] = fht[base + c] * fht[base + c] + fht[l] * fht[l]; } } // Other FHT-related methods moved here /* * converts FHT results (frequency space) to complex numbers of * [fftsize/2+1][fftsize] */ public double[][][] FHT2FFTHalf(FHT fht, int fftsize) { float[] fht_pixels = (float[]) fht.getPixels(); double[][][] fftHalf = new double[(fftsize >> 1) + 1][fftsize][2]; int row1, row2, col1, col2; for (row1 = 0; row1 <= (fftsize >> 1); row1++) { row2 = (fftsize - row1) % fftsize; for (col1 = 0; col1 < fftsize; col1++) { col2 = (fftsize - col1) % fftsize; fftHalf[row1][col1][0] = 0.5 * (fht_pixels[row1 * fftsize + col1] + fht_pixels[row2 * fftsize + col2]); fftHalf[row1][col1][1] = 0.5 * (fht_pixels[row2 * fftsize + col2] - fht_pixels[row1 * fftsize + col1]); } } return fftHalf; } public double[][][] FHT2FFTHalf(double[] fht_pixels, int fftsize) { double[][][] fftHalf = new double[(fftsize >> 1) + 1][fftsize][2]; int row1, row2, col1, col2; for (row1 = 0; row1 <= (fftsize >> 1); row1++) { row2 = (fftsize - row1) % fftsize; for (col1 = 0; col1 < fftsize; col1++) { col2 = (fftsize - col1) % fftsize; fftHalf[row1][col1][0] = 0.5 * (fht_pixels[row1 * fftsize + col1] + fht_pixels[row2 * fftsize + col2]); fftHalf[row1][col1][1] = 0.5 * (fht_pixels[row2 * fftsize + col2] - fht_pixels[row1 * fftsize + col1]); } } return fftHalf; } /* * converts FFT arrays of complex numbers of [fftsize/2+1][fftsize] to FHT * arrays */ public float[] floatFFTHalf2FHT(double[][][] fft, int fftsize) { float[] fht_pixels = new float[fftsize * fftsize]; int row1, row2, col1, col2; for (row1 = 0; row1 <= (fftsize >> 1); row1++) { row2 = (fftsize - row1) % fftsize; for (col1 = 0; col1 < fftsize; col1++) { col2 = (fftsize - col1) % fftsize; fht_pixels[row1 * fftsize + col1] = (float) (fft[row1][col1][0] - fft[row1][col1][1]); fht_pixels[row2 * fftsize + col2] = (float) (fft[row1][col1][0] + fft[row1][col1][1]); } } return fht_pixels; } public double[] FFTHalf2FHT(double[][][] fft, int fftsize) { double[] fht_pixels = new double[fftsize * fftsize]; int row1, row2, col1, col2; for (row1 = 0; row1 <= (fftsize >> 1); row1++) { row2 = (fftsize - row1) % fftsize; for (col1 = 0; col1 < fftsize; col1++) { col2 = (fftsize - col1) % fftsize; fht_pixels[row1 * fftsize + col1] = fft[row1][col1][0] - fft[row1][col1][1]; fht_pixels[row2 * fftsize + col2] = fft[row1][col1][0] + fft[row1][col1][1]; } } return fht_pixels; } // Uses just first half and one line public double[] FFTHalf2FHT(double[][] fft, int fftsize) { double[] fht_pixels = new double[fftsize * fftsize]; int row1, row2, col1, col2; for (row1 = 0; row1 <= (fftsize >> 1); row1++) { row2 = (fftsize - row1) % fftsize; for (col1 = 0; col1 < fftsize; col1++) { col2 = (fftsize - col1) % fftsize; fht_pixels[row1 * fftsize + col1] = fft[0][row1 * fftsize + col1] - fft[1][row1 * fftsize + col1]; fht_pixels[row2 * fftsize + col2] = fft[0][row1 * fftsize + col1] + fft[1][row1 * fftsize + col1]; } } return fht_pixels; } /* Amplitude/phase related methods */ public double[] interpolateFHT(double[] fht0, // first FHT array double[] fht1, // second FHT array double ratio) { // array of interpolation points - 0.0 - fht0, 1.0 - fht1 double[] points = { ratio }; double[][] results = interpolateFHT(fht0, // first FHT array fht1, // second FHT array points, // array of interpolation points - 0.0 - fht0, 1.0 - fht1 true); // do not clone 0.0 and 1.0 (ends) return results[0]; } /* * returns array of interpolated FHTs between fht0 and fht1, endpoints if * present (0.0, 1.0) are referenced, not cloned */ public double[][] interpolateFHT(double[] fht0, // first FHT array double[] fht1, // second FHT array double[] points) { // array of interpolation points - 0.0 - fht0, 1.0 - fht1 return interpolateFHT(fht0, // first FHT array fht1, // second FHT array points, // array of interpolation points - 0.0 - fht0, 1.0 - fht1 false); // do not clone 0.0 and 1.0 (ends) } public double[][] interpolateFHT(double[] fht0, // first FHT array double[] fht1, // second FHT array double[] points, // array of interpolation points - 0.0 - fht0, 1.0 - fht1 boolean cloneTrivial) { double[] fht_div = divide(fht1, fht0); int size = (int) Math.sqrt(fht0.length); int hsize = size / 2; double[][][] aphase = new double[hsize + 1][size][2]; double[][][] amp01 = new double[hsize + 1][size][2]; /* squared amplitudes of fht0 and fht1 */ double[][] phase = new double[hsize + 1][size]; /* +/-pi phase of the first array */ double[][][] fft0 = FHT2FFTHalf(fht0, size); double[][][] fft1 = FHT2FFTHalf(fht1, size); double[][][] fft_div = FHT2FFTHalf(fht_div, size); int i, j, k; double a, c, p; /* use mul for amplitudes, div - for phases */ for (i = 0; i <= hsize; i++) for (j = 0; j < size; j++) { amp01[i][j][0] = fft0[i][j][0] * fft0[i][j][0] + fft0[i][j][1] * fft0[i][j][1]; amp01[i][j][1] = fft1[i][j][0] * fft1[i][j][0] + fft1[i][j][1] * fft1[i][j][1]; if (amp01[i][j][0] > 0.0) phase[i][j] = Math.atan2(fft0[i][j][1], fft0[i][j][0]); else phase[i][j] = 0.0; aphase[i][j][0] = amp01[i][j][0] * amp01[i][j][1]; // product of squared amplitudes is OK for phase // restorations, we just need to know where // amplitudes are higher a = fft_div[i][j][0] * fft_div[i][j][0] + fft_div[i][j][1] * fft_div[i][j][1]; if (a > 0.0) aphase[i][j][1] = Math.atan2(fft_div[i][j][1], fft_div[i][j][0]); else aphase[i][j][1] = 0.0; } aphase[0][0][1] = 0.0; /* calculate full phases */ fullPhase(aphase); double[][] result = new double[points.length][]; for (k = 0; k < result.length; k++) { if (points[k] == 0.0) result[k] = cloneTrivial ? fht0.clone() : fht0; else if (points[k] == 1.0) result[k] = cloneTrivial ? fht1.clone() : fht1; else { /* interpolate */ c = points[k]; for (i = 0; i <= hsize; i++) for (j = 0; j < size; j++) { if ((amp01[i][j][0] == 0.0) || (amp01[i][j][1] == 0.0)) a = 0.0; /* * Extrapolation is defined here only in the direction of decreasing of the * spectral amplitudes (outside, to the wider PSF), so additional limit to * prevent division of small values */ /* * Seems to work, possible improvements: 1-filter spectrum in high-freq areas. 2 * - use farther inner points for farther approximation */ else if ((c < 0.0) && (amp01[i][j][0] > amp01[i][j][1])) a = Math.sqrt(amp01[i][j][0]); else if ((c > 1.0) && (amp01[i][j][0] < amp01[i][j][1])) a = Math.sqrt(amp01[i][j][1]); else a = Math.pow(amp01[i][j][0], 0.5 * (1.0 - c)) * Math.pow(amp01[i][j][1], 0.5 * c); p = phase[i][j] + c * aphase[i][j][1]; fft0[i][j][0] = a * Math.cos(p); fft0[i][j][1] = a * Math.sin(p); } result[k] = FFTHalf2FHT(fft0, size); } } return result; } /** * * @param fht FHT data * @param fullCentered - if false - returns minimal (size*(size/2-1) arrays, if * true - fills the full square arrays and swaps quadrants * @return 2-d array, first index: 0 - amplitude, 1 - phase. second index range * depends on fullCentered */ public double[][] fht2AmpHase(double[] fht, boolean fullCentered) { int size = (int) Math.sqrt(fht.length); int hsize = size / 2; double[][][] fft = FHT2FFTHalf(fht, size); double[][] aphase = new double[2][size * (fullCentered ? size : (hsize + 1))]; int index = 0; double a; for (int i = 0; i <= hsize; i++) for (int j = 0; j < size; j++) { a = Math.sqrt(fft[i][j][0] * fft[i][j][0] + fft[i][j][1] * fft[i][j][1]); aphase[0][index] = a; aphase[1][index++] = (a > 0.0) ? Math.atan2(fft[i][j][1], fft[i][j][0]) : 0.0; } if (fullCentered) { for (int i = 1; i < hsize; i++) for (int j = 0; j < size; j++) { int j1 = (size - j) % size; aphase[0][(size - i) * size + j] = aphase[0][i * size + j1]; aphase[1][(size - i) * size + j] = -aphase[1][i * size + j1]; } swapQuadrants(aphase[0]); swapQuadrants(aphase[1]); } return aphase; } public double[][] fht2ReIm(double[] fht, boolean fullCentered) { int size = (int) Math.sqrt(fht.length); int hsize = size / 2; double[][][] fft = FHT2FFTHalf(fht, size); double[][] reIm = new double[2][size * (fullCentered ? size : (hsize + 1))]; int index = 0; for (int i = 0; i <= hsize; i++) for (int j = 0; j < size; j++) { reIm[0][index] = fft[i][j][0]; reIm[1][index++] = fft[i][j][1]; } if (fullCentered) { for (int i = 1; i < hsize; i++) for (int j = 0; j < size; j++) { int j1 = (size - j) % size; reIm[0][(size - i) * size + j] = reIm[0][i * size + j1]; reIm[1][(size - i) * size + j] = -reIm[1][i * size + j1]; } swapQuadrants(reIm[0]); swapQuadrants(reIm[1]); } return reIm; } /* * replace +/-pi phase with the full phase, using amplitude to guide grouth of * the covered area, so amplitude and phase does not need to be a pair from the * same FFT array */ public void fullPhase(double[][][] aphase) { int size = aphase[0].length; int hsize = aphase.length - 1; boolean[][] map = new boolean[hsize + 1][size]; int i, j; aphase[0][0][1] = 0.0; int ix, iy, ix1, iy1, ix1n, iy1n; List<Integer> pixelList = new ArrayList<Integer>(100); Integer Index; int[][] dirs = { { -1, 0 }, { -1, -1 }, { 0, -1 }, { 1, -1 }, { 1, 0 }, { 1, 1 }, { 0, 1 }, { -1, 1 } }; ix = 0; iy = 0; int clusterSize = 0; boolean noNew = true; int maxX = 0; int maxY = 0; int oldX = 0; int oldY = 0; boolean oldConj = false; int listIndex; Index = iy * size + ix; pixelList.clear(); pixelList.add(Index); clusterSize++; map[iy][ix] = true; noNew = true; double phase, oldPhase, fullCyclesPhase; double maxValue = -1.0; while (pixelList.size() > 0) { /* Find maximal new neighbor */ maxValue = -1.0; listIndex = 0; while (listIndex < pixelList.size()) { Index = pixelList.get(listIndex); iy = Index / size; ix = Index % size; noNew = true; for (j = 0; j < 8; j++) { ix1 = (ix + dirs[j][0] + size) % size; iy1 = (iy + dirs[j][1] + size) % size; if ((iy1 > hsize) || (((iy1 == 0) || (iy1 == hsize)) && (ix1 > hsize))) { ix1n = (size - ix1) % size; iy1n = (size - iy1) % size; } else { /* But phase will be opposite sign */ ix1n = ix1; iy1n = iy1; } if (!map[iy1n][ix1n]) { noNew = false; if (aphase[iy1n][ix1n][0] > maxValue) { maxValue = aphase[iy1n][ix1n][0]; maxX = ix1n; maxY = iy1n; // if (DEBUG_LEVEL>4) System.out.println(" amplPhase(): iy="+iy+ " ix="+ix+" maxY="+maxY+" maxX="+maxX); } } } if (noNew) pixelList.remove(listIndex); // remove current list element else listIndex++; // increase list index } if (pixelList.size() == 0) break; /* * To calculate the phase - find already processed neighbor with the highest * amplitude */ maxValue = -1.0; for (j = 0; j < 8; j++) { ix1 = (maxX + dirs[j][0] + size) % size; iy1 = (maxY + dirs[j][1] + size) % size; if ((iy1 > hsize) || (((iy1 == 0) || (iy1 == hsize)) && (ix1 > hsize))) { ix1n = (size - ix1) % size; iy1n = (size - iy1) % size; } else { /* But phase will be opposite sign */ ix1n = ix1; iy1n = iy1; } if (map[iy1n][ix1n]) { if (aphase[iy1n][ix1n][0] > maxValue) { maxValue = aphase[iy1n][ix1n][0]; oldX = ix1n; oldY = iy1n; oldConj = (iy1 != iy1n) || (ix1 != ix1n); // point on the other half (conjugate) } } } /* Calculate the phase from the closest neighbor */ oldPhase = (oldConj ? -1 : 1) * aphase[oldY][oldX][1]; fullCyclesPhase = 2 * Math.PI * Math.floor(oldPhase / (2 * Math.PI) + 0.5); oldPhase -= fullCyclesPhase; // +/- pi phase = aphase[maxY][maxX][1]; if ((phase - oldPhase) > Math.PI) fullCyclesPhase -= 2 * Math.PI; else if ((oldPhase - phase) > Math.PI) fullCyclesPhase += 2 * Math.PI; aphase[maxY][maxX][1] += fullCyclesPhase; /* * if (DEBUG_LEVEL>3) { * System.out.println(" amplPhase():Old:["+oldConj+"] ("+pixelList.size()+") " * +oldX+":"+oldY+" "+IJ.d2s(aphase[oldY][oldX][0],2)+":"+ * IJ.d2s(aphase[oldY][oldX][1],2)+ * " New:"+maxX+":"+maxY+" "+IJ.d2s(aphase[maxY][maxX][0],2)+":"+ * IJ.d2s(aphase[maxY][maxX][1],2)+ * " Diff="+IJ.d2s((aphase[maxY][maxX][1]-(oldConj?-1:1)*aphase[oldY][oldX][1]), * 2)); } */ /* Add this new point to the list */ Index = maxY * size + maxX; pixelList.add(Index); clusterSize++; map[maxY][maxX] = true; } // end of while (pixelList.size()>0) /* Fix remaining phases for y=0 and y=hsize */ for (i = 1; i < hsize; i++) { aphase[0][size - i][1] = -aphase[0][i][1]; aphase[hsize][size - i][1] = -aphase[hsize][i][1]; } } }