Examples with Fht uk.ac.sussex.gdsc.core.ij.process.Fht used on opensource projects

Example 1 with Fht

use of uk.ac.sussex.gdsc.core.ij.process.Fht in project GDSC-SMLM by aherbert.

the class JTransformsTest method jtransforms2DDhtIsFasterThanFht2.

@SpeedTag
@SeededTest
void jtransforms2DDhtIsFasterThanFht2(RandomSeed seed) {
    Assumptions.assumeTrue(TestSettings.allow(TestComplexity.MEDIUM));
    // Test the forward DHT of data. and reverse transform or the pre-computed correlation.
    final int size = 256;
    final int w = size / 4;
    final UniformRandomProvider r = RngUtils.create(seed.getSeed());
    // Blob in the centre
    FloatProcessor fp = createProcessor(size, size / 2 - w / 2, size / 2 - w / 2, w, w, null);
    final Fht fht2 = new Fht((float[]) fp.getPixels(), size, false);
    fht2.transform();
    fht2.initialiseFastMultiply();
    // Random blobs, original and correlated
    final int N = 40;
    final float[][] data = new float[N * 2][];
    final int lower = w;
    final int upper = size - w;
    final int range = upper - lower;
    for (int i = 0, j = 0; i < N; i++) {
        final int x = lower + r.nextInt(range);
        final int y = lower + r.nextInt(range);
        fp = createProcessor(size, x, y, w, w, r);
        final float[] pixels = (float[]) fp.getPixels();
        data[j++] = pixels.clone();
        final Fht fht1 = new Fht(pixels, size, false);
        fht1.copyTables(fht2);
        fht2.transform();
        final float[] pixels2 = new float[pixels.length];
        fht2.conjugateMultiply(fht2, pixels2);
        data[j++] = pixels2;
    }
    // CommonUtils.setThreadsBeginN_1D_FFT_2Threads(Long.MAX_VALUE);
    // CommonUtils.setThreadsBeginN_1D_FFT_4Threads(Long.MAX_VALUE);
    CommonUtils.setThreadsBeginN_2D(Long.MAX_VALUE);
    final TimingService ts = new TimingService();
    ts.execute(new ImageJFhtSpeedTask(size, data));
    ts.execute(new ImageJFht2SpeedTask(size, data));
    ts.execute(new JTransformsDhtSpeedTask(size, data));
    ts.repeat();
    if (logger.isLoggable(Level.INFO)) {
        logger.info(ts.getReport());
    }
    // Assertions.assertTrue(ts.get(-1).getMean() < ts.get(-2).getMean());
    final TimingResult slow = ts.get(-2);
    final TimingResult fast = ts.get(-1);
    logger.log(TestLogUtils.getTimingRecord(slow, fast));
}

Example 2 with Fht

use of uk.ac.sussex.gdsc.core.ij.process.Fht in project GDSC-SMLM by aherbert.

the class DriftCalculator method calculateUsingImageStack.

/**
 * Calculates drift using images from a reference stack aligned to the overall z-projection.
 *
 * @param stack the stack
 * @param limits the limits
 * @return the drift { dx[], dy[] }
 */
@Nullable
private double[][] calculateUsingImageStack(ImageStack stack, int[] limits) {
    // Update the limits using the stack size
    final int upperT = settings.startFrame + settings.frameSpacing * (stack.getSize() - 1);
    limits[1] = Math.max(limits[1], upperT);
    // TODO - Truncate the stack if there are far too many frames for the localisation limits
    tracker.status("Constructing images");
    executor = Executors.newFixedThreadPool(Prefs.getThreads());
    // Built an image and FHT image for each slice
    final ImageProcessor[] images = new ImageProcessor[stack.getSize()];
    final Fht[] fhtImages = new Fht[stack.getSize()];
    final List<Future<?>> futures = new LinkedList<>();
    final Ticker ticker = Ticker.createStarted(tracker, images.length, true);
    final int imagesPerThread = getImagesPerThread(images);
    final AlignImagesFft aligner = new AlignImagesFft();
    final FloatProcessor referenceIp = stack.getProcessor(1).toFloat(0, null);
    // We do not care about the window method because this processor will not
    // actually be used for alignment, it is a reference for the FHT size
    aligner.initialiseReference(referenceIp, WindowMethod.NONE, false);
    for (int i = 0; i < images.length; i += imagesPerThread) {
        futures.add(executor.submit(new ImageFhtInitialiser(stack, images, aligner, fhtImages, i, i + imagesPerThread, ticker)));
    }
    ConcurrencyUtils.waitForCompletionUnchecked(futures);
    tracker.progress(1);
    if (tracker.isEnded()) {
        return null;
    }
    final double[] dx = new double[limits[1] + 1];
    final double[] dy = new double[dx.length];
    final double[] originalDriftTimePoints = new double[dx.length];
    final int[] blockT = new int[stack.getSize()];
    for (int i = 0, t = settings.startFrame; i < stack.getSize(); i++, t += settings.frameSpacing) {
        originalDriftTimePoints[t] = 1;
        blockT[i] = t;
    }
    final double smoothing = updateSmoothingParameter(originalDriftTimePoints);
    lastdx = null;
    // For the first iteration calculate drift to the first image in the stack
    // (since the average projection may have a large drift blurring the image)
    double change = calculateDriftUsingImageStack(referenceIp, images, fhtImages, blockT, dx, dy, originalDriftTimePoints, smoothing, settings.iterations);
    if (Double.isNaN(change) || tracker.isEnded()) {
        return null;
    }
    plotDrift(limits, dx, dy);
    ImageJUtils.log("Drift Calculator : Initial drift " + MathUtils.rounded(change));
    for (int i = 1; i <= settings.maxIterations; i++) {
        change = calculateDriftUsingImageStack(null, images, fhtImages, blockT, dx, dy, originalDriftTimePoints, smoothing, settings.iterations);
        if (Double.isNaN(change)) {
            return null;
        }
        plotDrift(limits, dx, dy);
        if (converged(i, change, getTotalDrift(dx, dy, originalDriftTimePoints))) {
            break;
        }
    }
    if (tracker.isEnded()) {
        return null;
    }
    plotDrift(limits, dx, dy);
    return new double[][] { dx, dy };
}

Example 3 with Fht

use of uk.ac.sussex.gdsc.core.ij.process.Fht 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) {
    final FloatProcessor taperedDataImage = getSquareTaperedImage(ip);
    final Fht fht = new Fht(taperedDataImage);
    fht.transform();
    return fht.getComplexTransformProcessors();
}

Example 4 with Fht

use of uk.ac.sussex.gdsc.core.ij.process.Fht 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
    final TrackProgress progess = getTrackProgress();
    progess.incrementProgress(0);
    progess.status("Calculating complex FFT images...");
    // Pad images to the same size if different
    final int maxWidth = Math.max(ip1.getWidth(), ip2.getWidth());
    final int maxHeight = Math.max(ip1.getHeight(), ip2.getHeight());
    if (Math.max(maxWidth, maxHeight) > MAX_SIZE) {
        progess.status("Error calculating FRC curve...");
        progess.incrementProgress(1);
        return null;
    }
    final int fieldOfView = Math.max(maxWidth, maxHeight);
    ip1 = pad(ip1, maxWidth, maxHeight);
    ip2 = pad(ip2, maxWidth, maxHeight);
    // The mean of each image after applying the taper
    double mean1;
    double mean2;
    // Real and imaginary components
    float[] re1 = null;
    float[] im1;
    float[] re2;
    float[] im2;
    // Do the first image
    ip1 = getSquareTaperedImage(ip1);
    mean1 = taperedImageMean;
    final int size = ip1.getWidth();
    if (fourierMethod == FourierMethod.JTRANSFORMS && JTransformsLoader.JTRANSFORMS_AVAILABLE) {
        // Speed up by reusing the FFT object which performs pre-computation
        final float[] data = new float[size * size * 2];
        final FloatFFT_2D fft = new FloatFFT_2D(size, size);
        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
        final float[] f1 = (float[]) ip1.getPixels();
        final Fht fht1 = new Fht(f1, ip1.getWidth(), false);
        fht1.transform();
        FloatProcessor[] fft = fht1.getComplexTransformProcessors();
        re1 = (float[]) fft[0].getPixels();
        im1 = (float[]) fft[1].getPixels();
        progess.incrementProgress(THIRD);
        ip2 = getSquareTaperedImage(ip2);
        mean2 = taperedImageMean;
        final float[] f2 = (float[]) ip2.getPixels();
        final Fht fht2 = new Fht(f2, ip2.getWidth(), false);
        fht2.copyTables(fht1);
        fft = fht2.getComplexTransformProcessors();
        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
    final float[] conjMult = new float[re1.length];
    final float[] absFft1 = new float[re1.length];
    final 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;
    }
    final 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;
    final FrcCurveResult[] results = new FrcCurveResult[max];
    if (samplingMethod == SamplingMethod.INTERPOLATED_CIRCLE) {
        // Set the results for the centre pixel
        final int cx = size * centre + centre;
        results[0] = new FrcCurveResult(0, 1, conjMult[cx], absFft1[cx], absFft2[cx]);
        final 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) {
                final double cosA = Math.cos(angle);
                final double x = centre + radius * cosA;
                final double sinA = getSine(angle, cosA);
                final double y = centre + radius * sinA;
                final 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
        final double[][] sum = RadialStatisticsUtils.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 5 with Fht

use of uk.ac.sussex.gdsc.core.ij.process.Fht in project GDSC-SMLM by aherbert.

the class FhtFilter method filterInternal.

/**
 * Compute the filter.
 *
 * <p>Note: the input data is destructively modified
 *
 * @param data the data
 * @param maxx the maxx
 * @param maxy the maxy
 * @param border the border
 */
private void filterInternal(float[] data, final int maxx, final int maxy, int border) {
    initialiseKernel(maxx, maxy);
    final Fht dataFht = createFht(data, maxx, maxy, border);
    final int maxN = kernelFht.getWidth();
    final Fht result = compute(dataFht);
    // Do the transform using JTransforms as it is faster
    dht.inverse(result.getData(), true);
    result.swapQuadrants();
    if (maxx < maxN || maxy < maxN) {
        final int x = getInsert(maxN, maxx);
        final int y = getInsert(maxN, maxy);
        for (int to = 0, from = y * maxN + x; to < data.length; to += maxx, from += maxN) {
            System.arraycopy(tmp, from, data, to, maxx);
        }
    } else {
        System.arraycopy(tmp, 0, data, 0, tmp.length);
    }
}