FHT multiply(FHT fht, boolean conjugate) { int rowMod, cMod, colMod; double h2e, h2o; float[] h1 = (float[])getPixels(); float[] h2 = (float[])fht.getPixels(); float[] tmp = new float[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) tmp[r * maxN + c] = (float)(h1[r * maxN + c] * h2e - h1[rowMod * maxN + colMod] * h2o); else tmp[r * maxN + c] = (float)(h1[r * maxN + c] * h2e + h1[rowMod * maxN + colMod] * h2o); } } FHT fht2 = new FHT(new FloatProcessor(maxN, maxN, tmp, null)); fht2.isFrequencyDomain = true; return fht2; }
void transform(boolean inverse) { //IJ.log("transform: "+maxN+" "+inverse); if (!powerOf2Size()) throw new IllegalArgumentException("Image not power of 2 size or not square: "+width+"x"+height); maxN = width; if (S==null) initializeTables(maxN); float[] fht = (float[])getPixels(); rc2DFHT(fht, inverse, maxN); isFrequencyDomain = !inverse; }
void transform(boolean inverse) { //IJ.log("transform: "+maxN+" "+inverse); if (!powerOf2Size()) throw new IllegalArgumentException("Image not power of 2 size or not square: "+width+"x"+height); maxN = width; if (S==null) initializeTables(maxN); float[] fht = (float[])getPixels(); rc2DFHT(fht, inverse, maxN); isFrequencyDomain = !inverse; }
/** 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 FHT divide(FHT fht) { int rowMod, cMod, colMod; double mag, h2e, h2o; float[] h1 = (float[])getPixels(); float[] h2 = (float[])fht.getPixels(); float[] out = new float[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); out[r*maxN+c] = (float)(tmp/mag); } } FHT fht2 = new FHT(new FloatProcessor(maxN, maxN, out, null)); fht2.isFrequencyDomain = true; return fht2; }
FHT multiply(FHT fht, boolean conjugate) { int rowMod, cMod, colMod; double h2e, h2o; float[] h1 = (float[])getPixels(); float[] h2 = (float[])fht.getPixels(); float[] tmp = new float[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) tmp[r * maxN + c] = (float)(h1[r * maxN + c] * h2e - h1[rowMod * maxN + colMod] * h2o); else tmp[r * maxN + c] = (float)(h1[r * maxN + c] * h2e + h1[rowMod * maxN + colMod] * h2o); } } FHT fht2 = new FHT(new FloatProcessor(maxN, maxN, tmp, null)); fht2.isFrequencyDomain = true; return fht2; }
/** 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 FHT divide(FHT fht) { int rowMod, cMod, colMod; double mag, h2e, h2o; float[] h1 = (float[])getPixels(); float[] h2 = (float[])fht.getPixels(); float[] out = new float[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); out[r*maxN+c] = (float)(tmp/mag); } } FHT fht2 = new FHT(new FloatProcessor(maxN, maxN, out, null)); fht2.isFrequencyDomain = true; return fht2; }
/** Complex to Complex Inverse Fourier Transform * Author: Joachim Wesner */ void c2c2DFFT(float[] rein, float[] imin, int maxN, float[] reout, float[] imout) { FHT fht = new FHT(new FloatProcessor(maxN,maxN)); float[] fhtpixels = (float[])fht.getPixels(); // Real part of inverse transform for (int iy = 0; iy < maxN; iy++) cplxFHT(iy, maxN, rein, imin, false, fhtpixels); fht.inverseTransform(); // Save intermediate result, so we can do a "in-place" transform float[] hlp = new float[maxN*maxN]; System.arraycopy(fhtpixels, 0, hlp, 0, maxN*maxN); // Imaginary part of inverse transform for (int iy = 0; iy < maxN; iy++) cplxFHT(iy, maxN, rein, imin, true, fhtpixels); fht.inverseTransform(); System.arraycopy(hlp, 0, reout, 0, maxN*maxN); System.arraycopy(fhtpixels, 0, imout, 0, maxN*maxN); }
/** Complex to Complex Inverse Fourier Transform * Author: Joachim Wesner */ void c2c2DFFT(float[] rein, float[] imin, int maxN, float[] reout, float[] imout) { FHT fht = new FHT(new FloatProcessor(maxN,maxN)); float[] fhtpixels = (float[])fht.getPixels(); // Real part of inverse transform for (int iy = 0; iy < maxN; iy++) cplxFHT(iy, maxN, rein, imin, false, fhtpixels); fht.inverseTransform(); // Save intermediate result, so we can do a "in-place" transform float[] hlp = new float[maxN*maxN]; System.arraycopy(fhtpixels, 0, hlp, 0, maxN*maxN); // Imaginary part of inverse transform for (int iy = 0; iy < maxN; iy++) cplxFHT(iy, maxN, rein, imin, true, fhtpixels); fht.inverseTransform(); System.arraycopy(hlp, 0, reout, 0, maxN*maxN); System.arraycopy(fhtpixels, 0, imout, 0, maxN*maxN); }
static Point2D.Double findMaximaWithSubpixelPrecision(Point2D.Double maximumCoords, int roiSize, FHT crossCorrelationImage) { double[] subImageData = new double[roiSize * roiSize]; float[] pixels = (float[]) crossCorrelationImage.getPixels(); int roiX = (int) maximumCoords.x - (roiSize - 1) / 2; int roiY = (int) maximumCoords.y - (roiSize - 1) / 2; if(isCloseToBorder((int) maximumCoords.x, (int) maximumCoords.y, (roiSize - 1) / 2, crossCorrelationImage)) { return maximumCoords; } for(int ys = roiY; ys < roiY + roiSize; ys++) { int offset1 = (ys - roiY) * roiSize; int offset2 = ys * crossCorrelationImage.getWidth() + roiX; for(int xs = 0; xs < roiSize; xs++) { subImageData[offset1++] = pixels[offset2++]; } } SubImage subImage = new SubImage(roiSize, roiSize, null, null, subImageData, 0, 0); RadialSymmetryFitter radialSymmetryFitter = new RadialSymmetryFitter(); Molecule psf = radialSymmetryFitter.fit(subImage); return new Point2D.Double((int) maximumCoords.x + psf.getX(), (int) maximumCoords.y + psf.getY()); }
void customFilter(FHT fht) { int size = fht.getWidth(); showStatus("Filtering"); fht.swapQuadrants(filter); float[] fhtPixels = (float[])fht.getPixels(); boolean isFloat = filter.getBitDepth()==32; for (int i=0; i<fhtPixels.length; i++) { if (isFloat) fhtPixels[i] = fhtPixels[i]*filter.getf(i); else fhtPixels[i] = (float)(fhtPixels[i]*(filter.get(i)/255.0)); } fht.swapQuadrants(filter); }
void customFilter(FHT fht) { int size = fht.getWidth(); showStatus("Filtering"); fht.swapQuadrants(filter); float[] fhtPixels = (float[])fht.getPixels(); boolean isFloat = filter.getBitDepth()==32; for (int i=0; i<fhtPixels.length; i++) { if (isFloat) fhtPixels[i] = fhtPixels[i]*filter.getf(i); else fhtPixels[i] = (float)(fhtPixels[i]*(filter.get(i)/255.0)); } fht.swapQuadrants(filter); }
/** Converts this FHT to a complex Fourier transform and returns it as a two slice stack. * Author: Joachim Wesner */ public ImageStack getComplexTransform() { if (!isFrequencyDomain) throw new IllegalArgumentException("Frequency domain image required"); float[] fht = (float[])getPixels(); float[] re = new float[maxN*maxN]; float[] im = new float[maxN*maxN]; for (int i=0; i<maxN; i++) { FHTreal(i, maxN, fht, re); FHTimag(i, maxN, fht, im); } swapQuadrants(new FloatProcessor(maxN, maxN, re, null)); swapQuadrants(new FloatProcessor(maxN, maxN, im, null)); ImageStack stack = new ImageStack(maxN, maxN); stack.addSlice("Real", re); stack.addSlice("Imaginary", im); return stack; }
/** Converts this FHT to a complex Fourier transform and returns it as a two slice stack. * Author: Joachim Wesner */ public ImageStack getComplexTransform() { if (!isFrequencyDomain) throw new IllegalArgumentException("Frequency domain image required"); float[] fht = (float[])getPixels(); float[] re = new float[maxN*maxN]; float[] im = new float[maxN*maxN]; for (int i=0; i<maxN; i++) { FHTreal(i, maxN, fht, re); FHTimag(i, maxN, fht, im); } swapQuadrants(new FloatProcessor(maxN, maxN, re, null)); swapQuadrants(new FloatProcessor(maxN, maxN, im, null)); ImageStack stack = new ImageStack(maxN, maxN); stack.addSlice("Real", re); stack.addSlice("Imaginary", im); return stack; }
float[] fps = new float[maxN*maxN]; byte[] ps = new byte[maxN*maxN]; float[] fht = (float[])getPixels();
float[] fps = new float[maxN*maxN]; byte[] ps = new byte[maxN*maxN]; float[] fht = (float[])getPixels();
/** * Multiplies a Fourier domain image by a filter * @param imp A frequency domain image, which is modified. * @param filter The filter, 32-bits (0-1) or 8-bits (0-255) * @see #forward * @see #inverse * @see #filter */ public static void multiply(ImagePlus imp, ImageProcessor filter) { Object obj = imp.getProperty("FHT"); FHT fht = obj!=null&&(obj instanceof FHT)?(FHT)obj:null; if (fht==null) return; int size = fht.getWidth(); boolean isFloat = filter.getBitDepth()==32; if (!isFloat) filter = filter.convertToByte(true); filter = filter.resize(size, size); fht.swapQuadrants(filter); float[] fhtPixels = (float[])fht.getPixels(); for (int i=0; i<fhtPixels.length; i++) { if (isFloat) fhtPixels[i] = fhtPixels[i]*filter.getf(i); else fhtPixels[i] = (float)(fhtPixels[i]*(filter.get(i)/255.0)); } fht.swapQuadrants(filter); imp.setProcessor(null, fht.getPowerSpectrum()); }
/** * Multiplies a Fourier domain image by a filter * @param imp A frequency domain image, which is modified. * @param filter The filter, 32-bits (0-1) or 8-bits (0-255) * @see #forward * @see #inverse * @see #filter */ public static void multiply(ImagePlus imp, ImageProcessor filter) { Object obj = imp.getProperty("FHT"); FHT fht = obj!=null&&(obj instanceof FHT)?(FHT)obj:null; if (fht==null) return; int size = fht.getWidth(); boolean isFloat = filter.getBitDepth()==32; if (!isFloat) filter = filter.convertToByte(true); filter = filter.resize(size, size); fht.swapQuadrants(filter); float[] fhtPixels = (float[])fht.getPixels(); for (int i=0; i<fhtPixels.length; i++) { if (isFloat) fhtPixels[i] = fhtPixels[i]*filter.getf(i); else fhtPixels[i] = (float)(fhtPixels[i]*(filter.get(i)/255.0)); } fht.swapQuadrants(filter); imp.setProcessor(null, fht.getPowerSpectrum()); }
void doMasking(FHT ip) { if (stackSize>1) return; float[] fht = (float[])ip.getPixels(); ImageProcessor mask = imp.getProcessor(); mask = mask.convertToByte(false);
void doMasking(FHT ip) { if (stackSize>1) return; float[] fht = (float[])ip.getPixels(); ImageProcessor mask = imp.getProcessor(); mask = mask.convertToByte(false);