/** * Convert an image into a Fourier image with real and imaginary parts * * @param ip - image * @return the real and imaginary parts */ public FloatProcessor[] getComplexFFT(ImageProcessor ip) { FloatProcessor taperedDataImage = getSquareTaperedImage(ip); FHT fht = new FHT(taperedDataImage); fht.setShowProgress( false ); fht.transform(); FloatProcessor[] ret = new FloatProcessor[2]; ImageStack stack1 = fht.getComplexTransform(); ret[0] = ((FloatProcessor) stack1.getProcessor(1)); ret[1] = ((FloatProcessor) stack1.getProcessor(2)); return ret; }
private static FHT computeCrossCorrelationImage(FHT image1FFT, FHT image2FFT) { FHT crossCorrelationImage = image1FFT.conjugateMultiply(image2FFT); crossCorrelationImage.setShowProgress(false); crossCorrelationImage.inverseTransform(); crossCorrelationImage.swapQuadrants(); return crossCorrelationImage; }
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; }
fht1 = (ImageProcessor)imp1.getProperty("FHT"); if (fht1!=null) h1 = new FHT(fht1); else { IJ.showStatus("Converting to float"); ImageProcessor ip1 = imp1.getProcessor(); h1 = new FHT(ip1); h2 = new FHT(fht2); else { ImageProcessor ip2 = imp2.getProcessor(); if (imp2!=imp1) h2 = new FHT(ip2); if (!h1.powerOf2Size()) { IJ.error("FFT Math", "Images must be a power of 2 size (256x256, 512x512, etc.)"); return; h1.transform(); h2 = new FHT(h1.duplicate()); else { IJ.showStatus("Transform image2"); h2.transform(); case CONJUGATE_MULTIPLY: IJ.showStatus("Complex conjugate multiply"); result = h1.conjugateMultiply(h2); break;
FHT fht = new FHT(ip2); fht.setShowProgress(false); fht.transform(); IJ.showProgress(9,20); fht.inverseTransform(); IJ.showProgress(19,20); fht.setRoi(fitRect); ip2 = fht.crop(); if (doScaling) { ImagePlus imp2 = new ImagePlus(imp.getTitle()+"-filtered", ip2);
private static FHT createPaddedFFTImage(FloatProcessor nextImage, int paddedSize) { FHT imageFFT = new FHT(Padding.PADDING_ZERO.padToBiggerSquare(nextImage, paddedSize)); imageFFT.setShowProgress(false); imageFFT.transform(); return imageFFT; }
void doMasking(FHT ip) { if (stackSize>1) return; float[] fht = (float[])ip.getPixels(); ImageProcessor mask = imp.getProcessor(); mask = mask.convertToByte(false); if (mask.getWidth()!=ip.getWidth() || mask.getHeight()!=ip.getHeight()) return; ImageStatistics stats = ImageStatistics.getStatistics(mask, MIN_MAX, null); ip.swapQuadrants(mask); byte[] maskPixels = (byte[])mask.getPixels(); for (int i=0; i<fht.length; i++) {
void doInverseTransform(FHT fht) { fht = fht.getCopy(); doMasking(fht); showStatus("Inverse transform"); fht.inverseTransform(); if (fht.quadrantSwapNeeded) fht.swapQuadrants(); fht.resetMinAndMax(); ImageProcessor ip2 = fht; if (fht.originalWidth>0) { fht.setRoi(0, 0, fht.originalWidth, fht.originalHeight); ip2 = fht.crop();
void swapQuadrants(ImageStack stack) { FHT fht = new FHT(new FloatProcessor(1, 1)); for (int i=1; i<=stack.getSize(); i++) fht.swapQuadrants(stack.getProcessor(i)); }
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; }
/** 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); }
/** * 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()); }
float[] fps = new float[maxN*maxN]; byte[] ps = new byte[maxN*maxN]; float[] fht = (float[])getPixels(); FHTps(row, maxN, fht, fps); base = row * maxN; for (int col=0; col<maxN; col++) { swapQuadrants(ip); if (FFT.displayRawPS) { ImageProcessor ip2 = new FloatProcessor(maxN, maxN, fps, null); swapQuadrants(ip2); new ImagePlus("PS of "+FFT.fileName, ip2).show(); ImageStack ct = getComplexTransform(); ImagePlus imp2 = new ImagePlus("Complex of "+FFT.fileName, ct); (new ContrastEnhancer()).stretchHistogram(imp2, 0.1);
void doFHTInverseTransform() { FHT fht = new FHT(imp.getProcessor().duplicate()); fht.inverseTransform(); fht.resetMinAndMax(); String name = WindowManager.getUniqueName(imp.getTitle().substring(7)); new ImagePlus(name, fht).show(); }
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; }
/** Returns a clone of this FHT. */ public FHT getCopy() { ImageProcessor ip = super.duplicate(); FHT fht = new FHT(ip); fht.isFrequencyDomain = isFrequencyDomain; fht.quadrantSwapNeeded = quadrantSwapNeeded; fht.rgb = rgb; fht.originalWidth = originalWidth; fht.originalHeight = originalHeight; fht.originalBitDepth = originalBitDepth; fht.originalColorModel = originalColorModel; return fht; }
void doInverseTransform(FHT fht, ImageProcessor ip) { showStatus("Inverse transform"); fht.inverseTransform(); //if (fht.quadrantSwapNeeded) // fht.swapQuadrants(); fht.resetMinAndMax(); ImageProcessor ip2 = fht; fht.setRoi(rect.x, rect.y, rect.width, rect.height); ip2 = fht.crop(); int bitDepth = fht.originalBitDepth>0?fht.originalBitDepth:imp.getBitDepth(); switch (bitDepth) { case 8: ip2 = ip2.convertToByte(true); break; case 16: ip2 = ip2.convertToShort(true); break; case 24: showStatus("Setting brightness"); fht.rgb.setBrightness((FloatProcessor)ip2); ip2 = fht.rgb; fht.rgb = null; break; case 32: break; } ip.insert(ip2, 0, 0); }
public void run(ImageProcessor ip) { slice++; if (done) return; FHT fht = newFHT(ip); if (slice==1) { filter = getFilter(fht.getWidth()); if (filter==null) { done = true; return; } } ((FHT)fht).transform(); customFilter(fht); doInverseTransform(fht, ip); if (slice==1) ip.resetMinAndMax(); if (slice==stackSize) { new ContrastEnhancer().stretchHistogram(imp, 0.0); imp.updateAndDraw(); } IJ.showProgress(1.0); if (Recorder.record && slice==1) Recorder.recordCall("FFT.filter(imp,filter); //see Help/Examples/JavaScript/FFT Filter"); }