Examples with FHT2 ij.process.FHT2 used on opensource projects

Example 1 with FHT2

use of ij.process.FHT2 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 FHT2

use of ij.process.FHT2 in project GDSC-SMLM by aherbert.

the class PCPALMAnalysis method padToFHT2.

/**
	 * Pads the image to the next power of two and transforms into the frequency domain
	 * 
	 * @param ip
	 * @return An FHT2 image in the frequency domain
	 */
private FHT2 padToFHT2(ImageProcessor ip) {
    FloatProcessor im2 = pad(ip);
    if (im2 == null)
        return null;
    FHT2 FHT2 = new FHT2(im2);
    FHT2.setShowProgress(false);
    FHT2.transform();
    return FHT2;
}

Example 3 with FHT2

use of ij.process.FHT2 in project GDSC-SMLM by aherbert.

the class FRC method getComplexFFT.

/**
	 * Convert an image into a Fourier image with real and imaginary parts
	 * 
	 * @param ip
	 *            The image
	 * @return the real and imaginary parts
	 */
public FloatProcessor[] getComplexFFT(ImageProcessor ip) {
    FloatProcessor taperedDataImage = getSquareTaperedImage(ip);
    FHT2 fht = new FHT2(taperedDataImage);
    fht.setShowProgress(false);
    fht.transform();
    ImageStack stack1 = fht.getComplexTransform();
    return getProcessors(stack1);
}

Example 4 with FHT2

use of ij.process.FHT2 in project GDSC-SMLM by aherbert.

the class PCPALMAnalysis method computeAutoCorrelationCurveFHT.

/**
	 * Compute the auto-correlation curve using FHT (ImageJ built-in). Computes the correlation
	 * image and then samples the image at radii up to the specified length to get the average
	 * correlation at a given radius.
	 * 
	 * @param im
	 * @param w
	 * @param maxRadius
	 * @param nmPerPixel
	 * @param density
	 * @return { distances[], gr[], gr_se[] }
	 */
private double[][] computeAutoCorrelationCurveFHT(ImageProcessor im, ImageProcessor w, int maxRadius, double nmPerPixel, double density) {
    log("Creating Hartley transforms");
    FHT2 fht2Im = padToFHT2(im);
    FHT2 fht2W = padToFHT2(w);
    if (fht2Im == null || fht2W == null) {
        error("Unable to perform Hartley transform");
        return null;
    }
    log("Performing correlation");
    FloatProcessor corrIm = computeAutoCorrelationFHT(fht2Im);
    FloatProcessor corrW = computeAutoCorrelationFHT(fht2W);
    IJ.showProgress(1);
    final int centre = corrIm.getHeight() / 2;
    Rectangle crop = new Rectangle(centre - maxRadius, centre - maxRadius, maxRadius * 2, maxRadius * 2);
    if (showCorrelationImages) {
        displayCorrelation(corrIm, "Image correlation", crop);
        displayCorrelation(corrW, "Window correlation", crop);
    }
    log("Normalising correlation");
    FloatProcessor correlation = normaliseCorrelation(corrIm, corrW, density);
    if (showCorrelationImages)
        displayCorrelation(correlation, "Normalised correlation", crop);
    return computeRadialAverage(maxRadius, nmPerPixel, correlation);
}

Example 5 with FHT2

use of ij.process.FHT2 in project GDSC-SMLM by aherbert.

the class PCPALMAnalysis method computeAutoCorrelationFHT.

/**
	 * @param fftIm
	 *            in frequency domain
	 * @return
	 */
private FloatProcessor computeAutoCorrelationFHT(FHT2 fftIm) {
    FHT2 FHT2 = fftIm.conjugateMultiply(fftIm);
    FHT2.inverseTransform();
    FHT2.swapQuadrants();
    return FHT2;
}