Search in sources :

Example 1 with TrackProgress

use of uk.ac.sussex.gdsc.core.logging.TrackProgress in project GDSC-SMLM by aherbert.

the class DarkTimeAnalysis method analyse.

private void analyse(MemoryPeakResults results) {
    // Find min and max time frames
    results.sort();
    final int min = results.getFirstFrame();
    final int max = results.getLastFrame();
    // Trace results:
    // TODO - The search distance could have units to avoid assuming the results are in pixels
    final double d = settings.searchDistance / results.getCalibrationReader().getNmPerPixel();
    int range = max - min + 1;
    if (settings.maxDarkTime > 0) {
        range = Math.max(1, (int) Math.round(settings.maxDarkTime * 1000 / msPerFrame));
    }
    final TrackProgress tracker = SimpleImageJTrackProgress.getInstance();
    tracker.status("Analysing ...");
    tracker.log("Analysing (d=%s nm (%s px) t=%s s (%d frames)) ...", MathUtils.rounded(settings.searchDistance), MathUtils.rounded(d), MathUtils.rounded(range * msPerFrame / 1000.0), range);
    Trace[] traces;
    if (settings.method == 0) {
        final TraceManager tm = new TraceManager(results);
        tm.setTracker(tracker);
        tm.traceMolecules(d, range);
        traces = tm.getTraces();
    } else {
        final ClusteringEngine engine = new ClusteringEngine(Prefs.getThreads(), algorithms[settings.method - 1], tracker);
        final List<Cluster> clusters = engine.findClusters(TraceMolecules.convertToClusterPoints(results), d, range);
        traces = TraceMolecules.convertToTraces(results, clusters);
    }
    tracker.status("Computing histogram ...");
    // Build dark-time histogram
    final int[] times = new int[range];
    final StoredData stats = new StoredData();
    for (final Trace trace : traces) {
        if (trace.getBlinks() > 1) {
            for (final int t : trace.getOffTimes()) {
                times[t]++;
            }
            stats.add(trace.getOffTimes());
        }
    }
    plotDarkTimeHistogram(stats);
    // Cumulative histogram
    for (int i = 1; i < times.length; i++) {
        times[i] += times[i - 1];
    }
    final int total = times[times.length - 1];
    // Plot dark-time up to 100%
    double[] x = new double[range];
    double[] y = new double[range];
    int truncate = 0;
    for (int i = 0; i < x.length; i++) {
        x[i] = i * msPerFrame;
        if (times[i] == total) {
            // Final value at 100%
            y[i] = 100.0;
            truncate = i + 1;
            break;
        }
        y[i] = (100.0 * times[i]) / total;
    }
    if (truncate > 0) {
        x = Arrays.copyOf(x, truncate);
        y = Arrays.copyOf(y, truncate);
    }
    final String title = "Cumulative Dark-time";
    final Plot plot = new Plot(title, "Time (ms)", "Percentile");
    plot.addPoints(x, y, Plot.LINE);
    ImageJUtils.display(title, plot);
    // Report percentile
    for (int i = 0; i < y.length; i++) {
        if (y[i] >= settings.percentile) {
            ImageJUtils.log("Dark-time Percentile %.1f @ %s ms = %s s", settings.percentile, MathUtils.rounded(x[i]), MathUtils.rounded(x[i] / 1000));
            break;
        }
    }
    tracker.status("");
}
Also used : Plot(ij.gui.Plot) SimpleImageJTrackProgress(uk.ac.sussex.gdsc.core.ij.SimpleImageJTrackProgress) TrackProgress(uk.ac.sussex.gdsc.core.logging.TrackProgress) Cluster(uk.ac.sussex.gdsc.core.clustering.Cluster) TraceManager(uk.ac.sussex.gdsc.smlm.results.TraceManager) Trace(uk.ac.sussex.gdsc.smlm.results.Trace) StoredData(uk.ac.sussex.gdsc.core.utils.StoredData) ClusteringEngine(uk.ac.sussex.gdsc.core.clustering.ClusteringEngine)

Example 2 with TrackProgress

use of uk.ac.sussex.gdsc.core.logging.TrackProgress 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);
}
Also used : FloatProcessor(ij.process.FloatProcessor) TrackProgress(uk.ac.sussex.gdsc.core.logging.TrackProgress) NullTrackProgress(uk.ac.sussex.gdsc.core.logging.NullTrackProgress) FloatFFT_2D(org.jtransforms.fft.FloatFFT_2D) Fht(uk.ac.sussex.gdsc.core.ij.process.Fht)

Aggregations

TrackProgress (uk.ac.sussex.gdsc.core.logging.TrackProgress)2 Plot (ij.gui.Plot)1 FloatProcessor (ij.process.FloatProcessor)1 FloatFFT_2D (org.jtransforms.fft.FloatFFT_2D)1 Cluster (uk.ac.sussex.gdsc.core.clustering.Cluster)1 ClusteringEngine (uk.ac.sussex.gdsc.core.clustering.ClusteringEngine)1 SimpleImageJTrackProgress (uk.ac.sussex.gdsc.core.ij.SimpleImageJTrackProgress)1 Fht (uk.ac.sussex.gdsc.core.ij.process.Fht)1 NullTrackProgress (uk.ac.sussex.gdsc.core.logging.NullTrackProgress)1 StoredData (uk.ac.sussex.gdsc.core.utils.StoredData)1 Trace (uk.ac.sussex.gdsc.smlm.results.Trace)1 TraceManager (uk.ac.sussex.gdsc.smlm.results.TraceManager)1