Examples with FloatFFT_2D org.jtransforms.fft.FloatFFT_2D used on opensource projects

Example 1 with FloatFFT_2D

use of org.jtransforms.fft.FloatFFT_2D in project GDSC-SMLM by aherbert.

the class FRC method calculateFrcCurve.

/**
	 * Calculate the Fourier Ring Correlation curve for two images.
	 *
	 * @param ip1
	 *            The first image
	 * @param ip2
	 *            The second image
	 * @param nmPerPixel
	 *            the nm per pixel for the super-resolution images
	 * @return the FRC curve
	 */
public FRCCurve calculateFrcCurve(ImageProcessor ip1, ImageProcessor ip2, double nmPerPixel) {
    // Allow a progress tracker to be input
    TrackProgress progess = getTrackProgress();
    progess.incrementProgress(0);
    progess.status("Calculating complex FFT images...");
    // Pad images to the same size if different
    final int maxWidth = FastMath.max(ip1.getWidth(), ip2.getWidth());
    final int maxHeight = FastMath.max(ip1.getHeight(), ip2.getHeight());
    if (FastMath.max(maxWidth, maxHeight) > MAX_SIZE) {
        progess.status("Error calculating FRC curve...");
        progess.incrementProgress(1);
        return null;
    }
    final int fieldOfView = FastMath.max(maxWidth, maxHeight);
    ip1 = pad(ip1, maxWidth, maxHeight);
    ip2 = pad(ip2, maxWidth, maxHeight);
    // The mean of each image after applying the taper
    double mean1, mean2;
    // Real and imaginary components
    float[] re1 = null, im1, re2, im2;
    // Do the first image
    ip1 = getSquareTaperedImage(ip1);
    mean1 = taperedImageMean;
    final int size = ip1.getWidth();
    if (JTRANSFORMS && fourierMethod == FourierMethod.JTRANSFORMS) {
        // Speed up by reusing the FFT object which performs pre-computation
        float[] data = new float[size * size * 2];
        FloatFFT_2D fft = new FloatFFT_2D(size, size);
        // For quadrant swap
        FHT2 fht = new FHT2();
        float[] pixels = (float[]) ip1.getPixels();
        System.arraycopy(pixels, 0, data, 0, pixels.length);
        fft.realForwardFull(data);
        // Get the data
        re1 = pixels;
        im1 = new float[pixels.length];
        for (int i = 0, j = 0; i < data.length; j++) {
            re1[j] = data[i++];
            im1[j] = data[i++];
        }
        fht.swapQuadrants(new FloatProcessor(size, size, re1));
        fht.swapQuadrants(new FloatProcessor(size, size, im1));
        progess.incrementProgress(THIRD);
        ip2 = getSquareTaperedImage(ip2);
        mean2 = taperedImageMean;
        pixels = (float[]) ip2.getPixels();
        System.arraycopy(pixels, 0, data, 0, pixels.length);
        for (int i = pixels.length; i < data.length; i++) data[i] = 0;
        fft.realForwardFull(data);
        // Get the data
        re2 = pixels;
        im2 = new float[pixels.length];
        for (int i = 0, j = 0; i < data.length; j++) {
            re2[j] = data[i++];
            im2[j] = data[i++];
        }
        fht.swapQuadrants(new FloatProcessor(size, size, re2));
        fht.swapQuadrants(new FloatProcessor(size, size, im2));
        progess.incrementProgress(THIRD);
    } else {
        // Simple implementation. This is left for testing.
        //FloatProcessor[] fft = getComplexFFT(ip1);
        //mean1 = taperedImageMean;
        //re1 = (float[]) fft[0].getPixels();
        //im1 = (float[]) fft[1].getPixels();
        //progess.incrementProgress(THIRD);
        //
        //fft = getComplexFFT(ip2);
        //mean2 = taperedImageMean;
        //re2 = (float[]) fft[0].getPixels();
        //im2 = (float[]) fft[1].getPixels();
        //progess.incrementProgress(THIRD);
        // Speed up by reusing the FHT object which performs pre-computation
        FHT2 fht = new FHT2();
        fht.setShowProgress(false);
        float[] f1 = (float[]) ip1.getPixels();
        fht.rc2DFHT(f1, false, size);
        FHT2 fht1 = new FHT2(ip1, true);
        FloatProcessor[] fft = getProcessors(fht1.getComplexTransform2());
        re1 = (float[]) fft[0].getPixels();
        im1 = (float[]) fft[1].getPixels();
        progess.incrementProgress(THIRD);
        ip2 = getSquareTaperedImage(ip2);
        mean2 = taperedImageMean;
        float[] f2 = (float[]) ip2.getPixels();
        fht.rc2DFHT(f2, false, size);
        FHT2 fht2 = new FHT2(ip2, true);
        fft = getProcessors(fht2.getComplexTransform2());
        re2 = (float[]) fft[0].getPixels();
        im2 = (float[]) fft[1].getPixels();
        progess.incrementProgress(THIRD);
    }
    progess.status("Preparing FRC curve calculation...");
    final int centre = size / 2;
    // In-line for speed 
    float[] conjMult = new float[re1.length];
    float[] absFFT1 = new float[re1.length];
    float[] absFFT2 = new float[re1.length];
    // Normalise the FFT to the field of view, i.e. normalise by 1/sqrt(N) for each dimension
    final double norm = 1.0 / fieldOfView;
    for (int i = 0; i < re1.length; i++) {
        re1[i] *= norm;
        im1[i] *= norm;
        re2[i] *= norm;
        im2[i] *= norm;
    }
    boolean basic = false;
    if (basic) {
        compute(conjMult, absFFT1, absFFT2, re1, im1, re2, im2);
    } else {
        computeMirroredFast(size, conjMult, absFFT1, absFFT2, re1, im1, re2, im2);
    }
    progess.status("Calculating FRC curve...");
    final int max = centre - 1;
    FRCCurveResult[] results = new FRCCurveResult[max];
    if (samplingMethod == SamplingMethod.INTERPOLATED_CIRCLE) {
        // Set the results for the centre pixel
        int cx = size * centre + centre;
        results[0] = new FRCCurveResult(0, 1, conjMult[cx], absFFT1[cx], absFFT2[cx]);
        float[][] images = new float[][] { conjMult, absFFT1, absFFT2 };
        for (int radius = 1; radius < max; radius++) {
            // Inline the calculation for speed
            double sum0 = 0;
            double sum1 = 0;
            double sum2 = 0;
            // Note: The image has 2-fold radial symmetry. So we only need to sample
            // angles from 0-pi. To sample the perimeter at pixel intervals we need
            // pi*r samples. So the angle step is max_angle / samples == pi / (pi*r) == 1 / r.
            // The number of samples is increased using the sampling factor.
            final double angleStep = 1 / (perimeterSamplingFactor * radius);
            double angle = 0;
            int numSum = 0;
            while (angle < Math.PI) {
                double cosA = FastMath.cos(angle);
                double x = centre + radius * cosA;
                //double sinA = FastMath.sin(angle);
                double sinA = getSine(angle, cosA);
                double y = centre + radius * sinA;
                double[] values = getInterpolatedValues(x, y, images, size);
                sum0 += values[0];
                sum1 += values[1];
                sum2 += values[2];
                numSum++;
                angle += angleStep;
            }
            results[radius] = new FRCCurveResult(radius, numSum, sum0, sum1, sum2);
        }
    } else {
        // Compute the radial sum as per the DIP image Matlab toolbox
        double[][] sum = RadialStatistics.radialSumAndCount(size, conjMult, absFFT1, absFFT2);
        for (int radius = 0; radius < max; radius++) {
            results[radius] = new FRCCurveResult(radius, (int) sum[3][radius], sum[0][radius], sum[1][radius], sum[2][radius]);
        }
    }
    progess.incrementProgress(LAST_THIRD);
    progess.status("Finished calculating FRC curve...");
    return new FRCCurve(nmPerPixel, fieldOfView, mean1, mean2, results);
}

Example 2 with FloatFFT_2D

use of org.jtransforms.fft.FloatFFT_2D in project GDSC-SMLM by aherbert.

the class PCPALMAnalysis method computeAutoCorrelationFFT.

/**
	 * Compute the auto-correlation using the JTransforms FFT library
	 * 
	 * @param ip
	 * @return
	 */
private FloatProcessor computeAutoCorrelationFFT(ImageProcessor ip) {
    FloatProcessor paddedIp = pad(ip);
    if (paddedIp == null)
        return null;
    final int size = paddedIp.getWidth();
    boolean doubleFFT = false;
    float[] pixels = (float[]) paddedIp.getPixels();
    float[] correlation = new float[size * size];
    if (doubleFFT) {
        DoubleFFT_2D fft = new DoubleFFT_2D(size, size);
        double[] data = new double[size * size * 2];
        for (int i = 0; i < pixels.length; i++) data[i] = pixels[i];
        fft.realForwardFull(data);
        // Re-use data
        for (int i = 0, j = 0; i < data.length; i += 2, j++) {
            data[j] = data[i] * data[i] + data[i + 1] * data[i + 1];
        }
        // Zero fill
        for (int j = correlation.length; j < data.length; j++) data[j] = 0;
        // Re-use the pre-computed object
        //fft = new DoubleFFT_2D(size, size);
        fft.realInverseFull(data, true);
        // Get the real part of the data
        for (int i = 0, j = 0; i < data.length; i += 2, j++) {
            correlation[j] = (float) data[i];
        }
    } else {
        FloatFFT_2D fft = new FloatFFT_2D(size, size);
        float[] data = new float[size * size * 2];
        System.arraycopy(pixels, 0, data, 0, pixels.length);
        fft.realForwardFull(data);
        // Re-use data
        for (int i = 0, j = 0; i < data.length; i += 2, j++) {
            data[j] = data[i] * data[i] + data[i + 1] * data[i + 1];
        }
        // Zero fill
        for (int j = correlation.length; j < data.length; j++) data[j] = 0;
        // Re-use the pre-computed object
        //fft = new FloatFFT_2D(size, size);
        fft.realInverseFull(data, true);
        // Get the real part of the data
        for (int i = 0, j = 0; i < data.length; i += 2, j++) {
            correlation[j] = data[i];
        }
    }
    // Swap quadrants
    FloatProcessor fp = new FloatProcessor(size, size, correlation, null);
    new FHT2().swapQuadrants(fp);
    return fp;
}