Search in sources :

Example 36 with Statistics

use of uk.ac.sussex.gdsc.core.utils.Statistics in project GDSC-SMLM by aherbert.

the class DiffusionRateTest method run.

@Override
public void run(String arg) {
    SmlmUsageTracker.recordPlugin(this.getClass(), arg);
    pluginSettings = Settings.load();
    pluginSettings.save();
    if (IJ.controlKeyDown()) {
        simpleTest();
        return;
    }
    extraOptions = ImageJUtils.isExtraOptions();
    if (!showDialog()) {
        return;
    }
    lastSimulation.set(null);
    final int totalSteps = (int) Math.ceil(settings.getSeconds() * settings.getStepsPerSecond());
    conversionFactor = 1000000.0 / (settings.getPixelPitch() * settings.getPixelPitch());
    // Diffusion rate is um^2/sec. Convert to pixels per simulation frame.
    final double diffusionRateInPixelsPerSecond = settings.getDiffusionRate() * conversionFactor;
    final double diffusionRateInPixelsPerStep = diffusionRateInPixelsPerSecond / settings.getStepsPerSecond();
    final double precisionInPixels = myPrecision / settings.getPixelPitch();
    final boolean addError = myPrecision != 0;
    ImageJUtils.log(TITLE + " : D = %s um^2/sec, Precision = %s nm", MathUtils.rounded(settings.getDiffusionRate(), 4), MathUtils.rounded(myPrecision, 4));
    ImageJUtils.log("Mean-displacement per dimension = %s nm/sec", MathUtils.rounded(1e3 * ImageModel.getRandomMoveDistance(settings.getDiffusionRate()), 4));
    if (extraOptions) {
        ImageJUtils.log("Step size = %s, precision = %s", MathUtils.rounded(ImageModel.getRandomMoveDistance(diffusionRateInPixelsPerStep)), MathUtils.rounded(precisionInPixels));
    }
    // Convert diffusion co-efficient into the standard deviation for the random walk
    final DiffusionType diffusionType = CreateDataSettingsHelper.getDiffusionType(settings.getDiffusionType());
    final double diffusionSigma = ImageModel.getRandomMoveDistance(diffusionRateInPixelsPerStep);
    ImageJUtils.log("Simulation step-size = %s nm", MathUtils.rounded(settings.getPixelPitch() * diffusionSigma, 4));
    // Move the molecules and get the diffusion rate
    IJ.showStatus("Simulating ...");
    final long start = System.nanoTime();
    final UniformRandomProvider random = UniformRandomProviders.create();
    final Statistics[] stats2D = new Statistics[totalSteps];
    final Statistics[] stats3D = new Statistics[totalSteps];
    final StoredDataStatistics jumpDistances2D = new StoredDataStatistics(totalSteps);
    final StoredDataStatistics jumpDistances3D = new StoredDataStatistics(totalSteps);
    for (int j = 0; j < totalSteps; j++) {
        stats2D[j] = new Statistics();
        stats3D[j] = new Statistics();
    }
    final SphericalDistribution dist = new SphericalDistribution(settings.getConfinementRadius() / settings.getPixelPitch());
    final Statistics asymptote = new Statistics();
    // Save results to memory
    final MemoryPeakResults results = new MemoryPeakResults(totalSteps);
    results.setCalibration(CalibrationHelper.create(settings.getPixelPitch(), 1, 1000.0 / settings.getStepsPerSecond()));
    results.setName(TITLE);
    results.setPsf(PsfHelper.create(PSFType.CUSTOM));
    int peak = 0;
    // Store raw coordinates
    final ArrayList<Point> points = new ArrayList<>(totalSteps);
    final StoredData totalJumpDistances1D = new StoredData(settings.getParticles());
    final StoredData totalJumpDistances2D = new StoredData(settings.getParticles());
    final StoredData totalJumpDistances3D = new StoredData(settings.getParticles());
    final NormalizedGaussianSampler gauss = SamplerUtils.createNormalizedGaussianSampler(random);
    for (int i = 0; i < settings.getParticles(); i++) {
        if (i % 16 == 0) {
            IJ.showProgress(i, settings.getParticles());
            if (ImageJUtils.isInterrupted()) {
                return;
            }
        }
        // Increment the frame so that tracing analysis can distinguish traces
        peak++;
        double[] origin = new double[3];
        final int id = i + 1;
        final MoleculeModel m = new MoleculeModel(id, origin.clone());
        if (addError) {
            origin = addError(origin, precisionInPixels, gauss);
        }
        if (pluginSettings.useConfinement) {
            // Note: When using confinement the average displacement should asymptote
            // at the average distance of a point from the centre of a ball. This is 3r/4.
            // See: http://answers.yahoo.com/question/index?qid=20090131162630AAMTUfM
            // The equivalent in 2D is 2r/3. However although we are plotting 2D distance
            // this is a projection of the 3D position onto the plane and so the particles
            // will not be evenly spread (there will be clustering at centre caused by the
            // poles)
            final double[] axis = (diffusionType == DiffusionType.LINEAR_WALK) ? nextVector(gauss) : null;
            for (int j = 0; j < totalSteps; j++) {
                double[] xyz = m.getCoordinates();
                final double[] originalXyz = xyz.clone();
                for (int n = pluginSettings.confinementAttempts; n-- > 0; ) {
                    if (diffusionType == DiffusionType.GRID_WALK) {
                        m.walk(diffusionSigma, random);
                    } else if (diffusionType == DiffusionType.LINEAR_WALK) {
                        m.slide(diffusionSigma, axis, random);
                    } else {
                        m.move(diffusionSigma, random);
                    }
                    if (!dist.isWithin(m.getCoordinates())) {
                        // Reset position
                        for (int k = 0; k < 3; k++) {
                            xyz[k] = originalXyz[k];
                        }
                    } else {
                        // The move was allowed
                        break;
                    }
                }
                points.add(new Point(id, xyz));
                if (addError) {
                    xyz = addError(xyz, precisionInPixels, gauss);
                }
                peak = record(xyz, id, peak, stats2D[j], stats3D[j], jumpDistances2D, jumpDistances3D, origin, results);
            }
            asymptote.add(distance(m.getCoordinates()));
        } else if (diffusionType == DiffusionType.GRID_WALK) {
            for (int j = 0; j < totalSteps; j++) {
                m.walk(diffusionSigma, random);
                double[] xyz = m.getCoordinates();
                points.add(new Point(id, xyz));
                if (addError) {
                    xyz = addError(xyz, precisionInPixels, gauss);
                }
                peak = record(xyz, id, peak, stats2D[j], stats3D[j], jumpDistances2D, jumpDistances3D, origin, results);
            }
        } else if (diffusionType == DiffusionType.LINEAR_WALK) {
            final double[] axis = nextVector(gauss);
            for (int j = 0; j < totalSteps; j++) {
                m.slide(diffusionSigma, axis, random);
                double[] xyz = m.getCoordinates();
                points.add(new Point(id, xyz));
                if (addError) {
                    xyz = addError(xyz, precisionInPixels, gauss);
                }
                peak = record(xyz, id, peak, stats2D[j], stats3D[j], jumpDistances2D, jumpDistances3D, origin, results);
            }
        } else {
            for (int j = 0; j < totalSteps; j++) {
                m.move(diffusionSigma, random);
                double[] xyz = m.getCoordinates();
                points.add(new Point(id, xyz));
                if (addError) {
                    xyz = addError(xyz, precisionInPixels, gauss);
                }
                peak = record(xyz, id, peak, stats2D[j], stats3D[j], jumpDistances2D, jumpDistances3D, origin, results);
            }
        }
        // Debug: record all the particles so they can be analysed
        // System.out.printf("%f %f %f\n", m.getX(), m.getY(), m.getZ());
        final double[] xyz = m.getCoordinates();
        double d2 = 0;
        totalJumpDistances1D.add(d2 = xyz[0] * xyz[0]);
        totalJumpDistances2D.add(d2 += xyz[1] * xyz[1]);
        totalJumpDistances3D.add(d2 += xyz[2] * xyz[2]);
    }
    final long nanoseconds = System.nanoTime() - start;
    IJ.showProgress(1);
    MemoryPeakResults.addResults(results);
    simulation = new SimulationData(results.getName(), myPrecision);
    // Convert pixels^2/step to um^2/sec
    final double msd2D = (jumpDistances2D.getMean() / conversionFactor) / (results.getCalibrationReader().getExposureTime() / 1000);
    final double msd3D = (jumpDistances3D.getMean() / conversionFactor) / (results.getCalibrationReader().getExposureTime() / 1000);
    ImageJUtils.log("Raw data D=%s um^2/s, Precision = %s nm, N=%d, step=%s s, mean2D=%s um^2, " + "MSD 2D = %s um^2/s, mean3D=%s um^2, MSD 3D = %s um^2/s", MathUtils.rounded(settings.getDiffusionRate()), MathUtils.rounded(myPrecision), jumpDistances2D.getN(), MathUtils.rounded(results.getCalibrationReader().getExposureTime() / 1000), MathUtils.rounded(jumpDistances2D.getMean() / conversionFactor), MathUtils.rounded(msd2D), MathUtils.rounded(jumpDistances3D.getMean() / conversionFactor), MathUtils.rounded(msd3D));
    aggregateIntoFrames(points, addError, precisionInPixels, gauss);
    IJ.showStatus("Analysing results ...");
    if (pluginSettings.showDiffusionExample) {
        showExample(totalSteps, diffusionSigma, random);
    }
    // Plot a graph of mean squared distance
    final double[] xValues = new double[stats2D.length];
    final double[] yValues2D = new double[stats2D.length];
    final double[] yValues3D = new double[stats3D.length];
    final double[] upper2D = new double[stats2D.length];
    final double[] lower2D = new double[stats2D.length];
    final double[] upper3D = new double[stats3D.length];
    final double[] lower3D = new double[stats3D.length];
    final SimpleRegression r2D = new SimpleRegression(false);
    final SimpleRegression r3D = new SimpleRegression(false);
    final int firstN = (pluginSettings.useConfinement) ? pluginSettings.fitN : totalSteps;
    for (int j = 0; j < totalSteps; j++) {
        // Convert steps to seconds
        xValues[j] = (j + 1) / settings.getStepsPerSecond();
        // Convert values in pixels^2 to um^2
        final double mean2D = stats2D[j].getMean() / conversionFactor;
        final double mean3D = stats3D[j].getMean() / conversionFactor;
        final double sd2D = stats2D[j].getStandardDeviation() / conversionFactor;
        final double sd3D = stats3D[j].getStandardDeviation() / conversionFactor;
        yValues2D[j] = mean2D;
        yValues3D[j] = mean3D;
        upper2D[j] = mean2D + sd2D;
        lower2D[j] = mean2D - sd2D;
        upper3D[j] = mean3D + sd3D;
        lower3D[j] = mean3D - sd3D;
        if (j < firstN) {
            r2D.addData(xValues[j], yValues2D[j]);
            r3D.addData(xValues[j], yValues3D[j]);
        }
    }
    // TODO - Fit using the equation for 2D confined diffusion:
    // MSD = 4s^2 + R^2 (1 - 0.99e^(-1.84^2 Dt / R^2)
    // s = localisation precision
    // R = confinement radius
    // D = 2D diffusion coefficient
    // t = time
    final PolynomialFunction fitted2D;
    final PolynomialFunction fitted3D;
    if (r2D.getN() > 0) {
        // Do linear regression to get diffusion rate
        final double[] best2D = new double[] { r2D.getIntercept(), r2D.getSlope() };
        fitted2D = new PolynomialFunction(best2D);
        final double[] best3D = new double[] { r3D.getIntercept(), r3D.getSlope() };
        fitted3D = new PolynomialFunction(best3D);
        // For 2D diffusion: d^2 = 4D
        // where: d^2 = mean-square displacement
        double diffCoeff = best2D[1] / 4.0;
        final String msg = "2D Diffusion rate = " + MathUtils.rounded(diffCoeff, 4) + " um^2 / sec (" + TextUtils.nanosToString(nanoseconds) + ")";
        IJ.showStatus(msg);
        ImageJUtils.log(msg);
        diffCoeff = best3D[1] / 6.0;
        ImageJUtils.log("3D Diffusion rate = " + MathUtils.rounded(diffCoeff, 4) + " um^2 / sec (" + TextUtils.nanosToString(nanoseconds) + ")");
    } else {
        fitted2D = fitted3D = null;
    }
    // Create plots
    plotMsd(totalSteps, xValues, yValues2D, lower2D, upper2D, fitted2D, 2);
    plotMsd(totalSteps, xValues, yValues3D, lower3D, upper3D, fitted3D, 3);
    plotJumpDistances(TITLE, jumpDistances2D, 2, 1);
    plotJumpDistances(TITLE, jumpDistances3D, 3, 1);
    // Show the total jump length for debugging
    // plotJumpDistances(TITLE + " total", totalJumpDistances1D, 1, totalSteps);
    // plotJumpDistances(TITLE + " total", totalJumpDistances2D, 2, totalSteps);
    // plotJumpDistances(TITLE + " total", totalJumpDistances3D, 3, totalSteps);
    windowOrganiser.tile();
    if (pluginSettings.useConfinement) {
        ImageJUtils.log("3D asymptote distance = %s nm (expected %.2f)", MathUtils.rounded(asymptote.getMean() * settings.getPixelPitch(), 4), 3 * settings.getConfinementRadius() / 4);
    }
}
Also used : SphericalDistribution(uk.ac.sussex.gdsc.smlm.model.SphericalDistribution) StoredDataStatistics(uk.ac.sussex.gdsc.core.utils.StoredDataStatistics) ArrayList(java.util.ArrayList) PolynomialFunction(org.apache.commons.math3.analysis.polynomials.PolynomialFunction) StoredDataStatistics(uk.ac.sussex.gdsc.core.utils.StoredDataStatistics) Statistics(uk.ac.sussex.gdsc.core.utils.Statistics) MoleculeModel(uk.ac.sussex.gdsc.smlm.model.MoleculeModel) SimpleRegression(org.apache.commons.math3.stat.regression.SimpleRegression) DiffusionType(uk.ac.sussex.gdsc.smlm.model.DiffusionType) StoredData(uk.ac.sussex.gdsc.core.utils.StoredData) UniformRandomProvider(org.apache.commons.rng.UniformRandomProvider) MemoryPeakResults(uk.ac.sussex.gdsc.smlm.results.MemoryPeakResults) NormalizedGaussianSampler(org.apache.commons.rng.sampling.distribution.NormalizedGaussianSampler)

Example 37 with Statistics

use of uk.ac.sussex.gdsc.core.utils.Statistics in project GDSC-SMLM by aherbert.

the class BenchmarkFit method runFit.

private void runFit() {
    // Initialise the answer.
    answer[Gaussian2DFunction.BACKGROUND] = benchmarkParameters.getBackground();
    answer[Gaussian2DFunction.SIGNAL] = benchmarkParameters.getSignal();
    answer[Gaussian2DFunction.X_POSITION] = benchmarkParameters.x;
    answer[Gaussian2DFunction.Y_POSITION] = benchmarkParameters.y;
    answer[Gaussian2DFunction.Z_POSITION] = benchmarkParameters.z;
    answer[Gaussian2DFunction.X_SD] = benchmarkParameters.sd / benchmarkParameters.pixelPitch;
    answer[Gaussian2DFunction.Y_SD] = benchmarkParameters.sd / benchmarkParameters.pixelPitch;
    // Set up the fit region. Always round down since 0.5 is the centre of the pixel.
    final int x = (int) benchmarkParameters.x;
    final int y = (int) benchmarkParameters.y;
    region = new Rectangle(x - regionSize, y - regionSize, 2 * regionSize + 1, 2 * regionSize + 1);
    if (!new Rectangle(0, 0, imp.getWidth(), imp.getHeight()).contains(region)) {
        // Check if it is incorrect by only 1 pixel
        if (region.width <= imp.getWidth() + 1 && region.height <= imp.getHeight() + 1) {
            ImageJUtils.log("Adjusting region %s to fit within image bounds (%dx%d)", region.toString(), imp.getWidth(), imp.getHeight());
            region = new Rectangle(0, 0, imp.getWidth(), imp.getHeight());
        } else {
            IJ.error(TITLE, "Fit region does not fit within the image");
            return;
        }
    }
    // Adjust the centre & account for 0.5 pixel offset during fitting
    answer[Gaussian2DFunction.X_POSITION] -= (region.x + 0.5);
    answer[Gaussian2DFunction.Y_POSITION] -= (region.y + 0.5);
    // Configure for fitting
    fitConfig.setBackgroundFitting(backgroundFitting);
    fitConfig.setNotSignalFitting(!signalFitting);
    fitConfig.setComputeDeviations(false);
    // Create the camera model
    CameraModel cameraModel = fitConfig.getCameraModel();
    // Crop for speed. Reset origin first so the region is within the model
    cameraModel.setOrigin(0, 0);
    cameraModel = cameraModel.crop(region, false);
    final ImageStack stack = imp.getImageStack();
    final int totalFrames = benchmarkParameters.frames;
    // Create a pool of workers
    final int nThreads = Prefs.getThreads();
    final BlockingQueue<Integer> jobs = new ArrayBlockingQueue<>(nThreads * 2);
    final List<Worker> workers = new LinkedList<>();
    final List<Thread> threads = new LinkedList<>();
    final Ticker ticker = ImageJUtils.createTicker(totalFrames, nThreads, "Fitting frames ...");
    for (int i = 0; i < nThreads; i++) {
        final Worker worker = new Worker(jobs, stack, region, fitConfig, cameraModel, ticker);
        final Thread t = new Thread(worker);
        workers.add(worker);
        threads.add(t);
        t.start();
    }
    // Store all the fitting results
    results = new double[totalFrames * startPoints.length][];
    resultsTime = new long[results.length];
    // Fit the frames
    for (int i = 0; i < totalFrames; i++) {
        // Only fit if there were simulated photons
        if (benchmarkParameters.framePhotons[i] > 0) {
            put(jobs, i);
        }
    }
    // Finish all the worker threads by passing in a null job
    for (int i = 0; i < threads.size(); i++) {
        put(jobs, -1);
    }
    // Wait for all to finish
    for (int i = 0; i < threads.size(); i++) {
        try {
            threads.get(i).join();
        } catch (final InterruptedException ex) {
            Thread.currentThread().interrupt();
            throw new ConcurrentRuntimeException(ex);
        }
    }
    threads.clear();
    if (hasOffsetXy()) {
        ImageJUtils.log(TITLE + ": CoM within start offset = %d / %d (%s%%)", comValid.intValue(), totalFrames, MathUtils.rounded((100.0 * comValid.intValue()) / totalFrames));
    }
    ImageJUtils.finished("Collecting results ...");
    // Collect the results
    Statistics[] stats = null;
    for (int i = 0; i < workers.size(); i++) {
        final Statistics[] next = workers.get(i).stats;
        if (stats == null) {
            stats = next;
            continue;
        }
        for (int j = 0; j < next.length; j++) {
            stats[j].add(next[j]);
        }
    }
    workers.clear();
    Objects.requireNonNull(stats, "No statistics were computed");
    // Show a table of the results
    summariseResults(stats, cameraModel);
    // Optionally show histograms
    if (showHistograms) {
        IJ.showStatus("Calculating histograms ...");
        final WindowOrganiser windowOrganiser = new WindowOrganiser();
        final double[] convert = getConversionFactors();
        final HistogramPlotBuilder builder = new HistogramPlotBuilder(TITLE).setNumberOfBins(histogramBins);
        for (int i = 0; i < NAMES.length; i++) {
            if (displayHistograms[i] && convert[i] != 0) {
                // We will have to convert the values...
                final double[] tmp = ((StoredDataStatistics) stats[i]).getValues();
                for (int j = 0; j < tmp.length; j++) {
                    tmp[j] *= convert[i];
                }
                final StoredDataStatistics tmpStats = StoredDataStatistics.create(tmp);
                builder.setData(tmpStats).setName(NAMES[i]).setPlotLabel(String.format("%s +/- %s", MathUtils.rounded(tmpStats.getMean()), MathUtils.rounded(tmpStats.getStandardDeviation()))).show(windowOrganiser);
            }
        }
        windowOrganiser.tile();
    }
    if (saveRawData) {
        final String dir = ImageJUtils.getDirectory("Data_directory", rawDataDirectory);
        if (dir != null) {
            saveData(stats, dir);
        }
    }
    IJ.showStatus("");
}
Also used : CameraModel(uk.ac.sussex.gdsc.smlm.model.camera.CameraModel) ImageStack(ij.ImageStack) Ticker(uk.ac.sussex.gdsc.core.logging.Ticker) Rectangle(java.awt.Rectangle) StoredDataStatistics(uk.ac.sussex.gdsc.core.utils.StoredDataStatistics) HistogramPlotBuilder(uk.ac.sussex.gdsc.core.ij.HistogramPlot.HistogramPlotBuilder) WindowOrganiser(uk.ac.sussex.gdsc.core.ij.plugin.WindowOrganiser) StoredDataStatistics(uk.ac.sussex.gdsc.core.utils.StoredDataStatistics) Statistics(uk.ac.sussex.gdsc.core.utils.Statistics) LinkedList(java.util.LinkedList) AtomicInteger(java.util.concurrent.atomic.AtomicInteger) ConcurrentRuntimeException(org.apache.commons.lang3.concurrent.ConcurrentRuntimeException) ArrayBlockingQueue(java.util.concurrent.ArrayBlockingQueue)

Example 38 with Statistics

use of uk.ac.sussex.gdsc.core.utils.Statistics in project GDSC-SMLM by aherbert.

the class BenchmarkSmartSpotRanking method summariseResults.

private void summariseResults(TIntObjectHashMap<RankResults> rankResults, CandidateData candidateData) {
    // Summarise the ranking results.
    final StringBuilder sb = new StringBuilder(filterResult.resultPrefix);
    // countPositive and countNegative is the fractional score of the spot candidates
    addCount(sb, (double) candidateData.countPositive + candidateData.countNegative);
    addCount(sb, candidateData.countPositive);
    addCount(sb, candidateData.countNegative);
    addCount(sb, candidateData.fractionPositive);
    addCount(sb, candidateData.fractionNegative);
    final double[] counter1 = new double[2];
    final int[] counter2 = new int[2];
    candidateData.filterCandidates.forEachValue(result -> {
        counter1[0] += result.np;
        counter1[1] += result.nn;
        counter2[0] += result.pos;
        counter2[1] += result.neg;
        return true;
    });
    final int countTp = counter2[0];
    final int countFp = counter2[1];
    // The fraction of positive and negative candidates that were included
    add(sb, (100.0 * countTp) / candidateData.countPositive);
    add(sb, (100.0 * countFp) / candidateData.countNegative);
    // Add counts of the the candidates
    add(sb, countTp + countFp);
    add(sb, countTp);
    add(sb, countFp);
    // Add fractional counts of the the candidates
    double tp = counter1[0];
    double fp = counter1[1];
    add(sb, tp + fp);
    add(sb, tp);
    add(sb, fp);
    // Materialise rankResults
    final int[] frames = new int[rankResults.size()];
    final RankResults[] rankResultsArray = new RankResults[rankResults.size()];
    final int[] counter = new int[1];
    rankResults.forEachEntry((frame, result) -> {
        frames[counter[0]] = frame;
        rankResultsArray[counter[0]] = result;
        counter[0]++;
        return true;
    });
    // Summarise actual and candidate spots per frame
    final Statistics actual = new Statistics();
    final Statistics candidates = new Statistics();
    for (final RankResults rr : rankResultsArray) {
        actual.add(rr.zPosition.length);
        candidates.add(rr.spots.length);
    }
    add(sb, actual.getMean());
    add(sb, actual.getStandardDeviation());
    add(sb, candidates.getMean());
    add(sb, candidates.getStandardDeviation());
    final String resultPrefix = sb.toString();
    // ---
    // TODO: Further explore pre-filtering of spot candidates.
    // Add good label to spot candidates and have the benchmark spot filter respect this before
    // applying the fail count limit.
    // Correlation between intensity and SNR ...
    // SNR is very good at low density
    // SNR fails at high density. The SNR estimate is probably wrong for high intensity spots.
    // Triangle is very good when there are a large number of good spots in a region of the image
    // (e.g. a mask is used).
    // Triangle is poor when there are few good spots in an image.
    // Perhaps we can estimate the density of the spots and choose the correct thresholding method?
    // ---
    // Do a full benchmark through different Spot SNR, image sizes, densities and mask structures
    // and see if there are patterns for a good threshold method.
    // ---
    // Allow using the fitted results from benchmark spot fit. Will it make a difference if we fit
    // the candidates (some will fail if weak).
    // Can this be done by allowing the user to select the input (spot candidates or fitted
    // positions)?
    // Perhaps I need to produce a precision estimate for all simulated spots and then only use
    // those that achieve a certain precision, i.e. are reasonably in focus. Can this be done?
    // Does the image PSF have a width estimate for the entire stack?
    // Perhaps I should filter, fit and then filter all spots using no fail count. These then become
    // the spots to work with for creating a smart fail count filter.
    // ---
    // Pre-compute the results and have optional sort
    final ArrayList<ScoredResult> list = new ArrayList<>(methodNames.length);
    for (int i = 0; i < methodNames.length; i++) {
        tp = 0;
        fp = 0;
        double tn = 0;
        int itp = 0;
        int ifp = 0;
        int itn = 0;
        final Statistics s = new Statistics();
        long time = 0;
        for (final RankResults rr : rankResultsArray) {
            final RankResult r = rr.results.get(i);
            // Some results will not have a threshold
            if (!Float.isInfinite(r.threshold)) {
                s.add(r.threshold);
            }
            time += r.time;
            tp += r.fresult.getTruePositives();
            fp += r.fresult.getFalsePositives();
            tn += r.fresult.getTrueNegatives();
            itp += r.cresult.getTruePositives();
            ifp += r.cresult.getFalsePositives();
            itn += r.cresult.getTrueNegatives();
        }
        sb.setLength(0);
        sb.append(resultPrefix);
        add(sb, methodNames[i]);
        if (methodNames[i].startsWith("SNR")) {
            sb.append('\t');
        } else {
            add(sb, settings.compactBins);
        }
        add(sb, s.getMean());
        add(sb, s.getStandardDeviation());
        add(sb, TextUtils.nanosToString(time));
        // TP are all accepted candidates that can be matched to a spot
        // FP are all accepted candidates that cannot be matched to a spot
        // TN are all accepted candidates that cannot be matched to a spot
        // FN = The number of missed spots
        // Raw counts of match or no-match
        final FractionClassificationResult f1 = new FractionClassificationResult(itp, ifp, itn, simulationParameters.molecules - itp);
        final double s1 = addScores(sb, f1);
        // Fractional scoring
        final FractionClassificationResult f2 = new FractionClassificationResult(tp, fp, tn, simulationParameters.molecules - tp);
        final double s2 = addScores(sb, f2);
        // Store for sorting
        list.add(new ScoredResult(i, (settings.useFractionScores) ? s2 : s1, sb.toString()));
    }
    if (list.isEmpty()) {
        return;
    }
    Collections.sort(list, ScoredResult::compare);
    final TextWindow summaryTable = createTable();
    if (summaryTable.getTextPanel().getLineCount() > 0) {
        summaryTable.append("");
    }
    try (BufferedTextWindow tw = new BufferedTextWindow(summaryTable)) {
        tw.setIncrement(0);
        for (final ScoredResult r : list) {
            tw.append(r.result);
        }
    }
    if (settings.showOverlay) {
        final int bestMethod = list.get(0).index;
        final Overlay o = new Overlay();
        for (int j = 0; j < rankResultsArray.length; j++) {
            final int frame = frames[j];
            final RankResults rr = rankResultsArray[j];
            final RankResult r = rr.results.get(bestMethod);
            final int[] x1 = new int[r.good.length];
            final int[] y1 = new int[r.good.length];
            int c1 = 0;
            final int[] x2 = new int[r.good.length];
            final int[] y2 = new int[r.good.length];
            int c2 = 0;
            final int[] x3 = new int[r.good.length];
            final int[] y3 = new int[r.good.length];
            int c3 = 0;
            final int[] x4 = new int[r.good.length];
            final int[] y4 = new int[r.good.length];
            int c4 = 0;
            for (int i = 0; i < x1.length; i++) {
                if (r.good[i] == TP) {
                    x1[c1] = rr.spots[i].spot.x;
                    y1[c1] = rr.spots[i].spot.y;
                    c1++;
                } else if (r.good[i] == FP) {
                    x2[c2] = rr.spots[i].spot.x;
                    y2[c2] = rr.spots[i].spot.y;
                    c2++;
                } else if (r.good[i] == TN) {
                    x3[c3] = rr.spots[i].spot.x;
                    y3[c3] = rr.spots[i].spot.y;
                    c3++;
                } else if (r.good[i] == FN) {
                    x4[c4] = rr.spots[i].spot.x;
                    y4[c4] = rr.spots[i].spot.y;
                    c4++;
                }
            }
            addToOverlay(o, frame, x1, y1, c1, Color.green);
            addToOverlay(o, frame, x2, y2, c2, Color.red);
            if (IJ.debugMode) {
                // light green
                addToOverlay(o, frame, x3, y3, c3, new Color(153, 255, 153));
            }
            // light red
            addToOverlay(o, frame, x4, y4, c4, new Color(255, 153, 153));
        }
        imp.setOverlay(o);
    }
}
Also used : BufferedTextWindow(uk.ac.sussex.gdsc.core.ij.BufferedTextWindow) Color(java.awt.Color) TDoubleArrayList(gnu.trove.list.array.TDoubleArrayList) ArrayList(java.util.ArrayList) Statistics(uk.ac.sussex.gdsc.core.utils.Statistics) PeakResultPoint(uk.ac.sussex.gdsc.smlm.results.PeakResultPoint) TextWindow(ij.text.TextWindow) BufferedTextWindow(uk.ac.sussex.gdsc.core.ij.BufferedTextWindow) FractionClassificationResult(uk.ac.sussex.gdsc.core.match.FractionClassificationResult) Overlay(ij.gui.Overlay)

Example 39 with Statistics

use of uk.ac.sussex.gdsc.core.utils.Statistics in project GDSC-SMLM by aherbert.

the class BenchmarkSpotFilter method getStats.

private double[][] getStats(ArrayList<BatchResult[]> batchResults) {
    if (settings.selectionMethod < 2) {
        // 1 = Max Jaccard
        return null;
    }
    // 2 = AUC+Max Jaccard so we require the mean and standard deviation of each
    // score to normalise them
    final double[][] stats = new double[2][2];
    for (int index = 0; index < stats.length; index++) {
        final Statistics s = new Statistics();
        for (final BatchResult[] batchResult : batchResults) {
            if (batchResult == null || batchResult.length == 0) {
                continue;
            }
            for (int i = 0; i < batchResult.length; i++) {
                s.add(batchResult[i].getScore(index));
            }
        }
        stats[index][0] = s.getMean();
        stats[index][1] = s.getStandardDeviation();
    }
    return stats;
}
Also used : Statistics(uk.ac.sussex.gdsc.core.utils.Statistics) PeakResultPoint(uk.ac.sussex.gdsc.smlm.results.PeakResultPoint) BasePoint(uk.ac.sussex.gdsc.core.match.BasePoint)

Example 40 with Statistics

use of uk.ac.sussex.gdsc.core.utils.Statistics in project GDSC-SMLM by aherbert.

the class PcPalmAnalysis method computeRadialAverage.

/**
 * Compute the radial average correlation function (gr).
 *
 * @param maxRadius the maximum radius to process (in pixels)
 * @param nmPerPixel covert pixel distances to nm
 * @param correlation auto-correlation
 * @return { distances[], gr[], gr_se[] }
 */
private static double[][] computeRadialAverage(int maxRadius, double nmPerPixel, FloatProcessor correlation) {
    // Perform averaging of the correlation function using integer distance bins
    log("  Computing distance vs correlation curve");
    final int centre = correlation.getHeight() / 2;
    // Count the number of pixels at each distance and sum the correlations
    final Statistics[] gr = new Statistics[maxRadius + 1];
    // Cache distances
    final double[] d2 = new double[maxRadius + 1];
    for (int dy = 0; dy <= maxRadius; dy++) {
        gr[dy] = new Statistics();
        d2[dy] = (double) dy * dy;
    }
    final int[][] distance = new int[maxRadius + 1][maxRadius + 1];
    for (int dy = 0; dy <= maxRadius; dy++) {
        for (int dx = dy; dx <= maxRadius; dx++) {
            distance[dy][dx] = distance[dx][dy] = (int) Math.round(Math.sqrt(d2[dx] + d2[dy]));
        }
    }
    final float[] data = (float[]) correlation.getPixels();
    for (int dy = -maxRadius; dy <= maxRadius; dy++) {
        final int absY = Math.abs(dy);
        int index = (centre + dy) * correlation.getWidth() + centre - maxRadius;
        for (int dx = -maxRadius; dx <= maxRadius; dx++, index++) {
            final int d = distance[absY][Math.abs(dx)];
            if (d > maxRadius || d == 0) {
                continue;
            }
            gr[d].add(data[index]);
        }
    }
    // Create the final data: a curve showing distance (in nm) verses the average correlation
    final double[] x = new double[maxRadius + 1];
    final double[] y = new double[maxRadius + 1];
    final double[] sd = new double[maxRadius + 1];
    for (int i = 0; i < x.length; i++) {
        x[i] = i * nmPerPixel;
        y[i] = gr[i].getMean();
        sd[i] = gr[i].getStandardError();
    }
    y[0] = correlation.getf(centre, centre);
    return new double[][] { x, y, sd };
}
Also used : Statistics(uk.ac.sussex.gdsc.core.utils.Statistics)

Aggregations

Statistics (uk.ac.sussex.gdsc.core.utils.Statistics)46 StoredDataStatistics (uk.ac.sussex.gdsc.core.utils.StoredDataStatistics)16 MemoryPeakResults (uk.ac.sussex.gdsc.smlm.results.MemoryPeakResults)11 Rectangle (java.awt.Rectangle)10 ArrayList (java.util.ArrayList)10 WindowOrganiser (uk.ac.sussex.gdsc.core.ij.plugin.WindowOrganiser)10 LocalList (uk.ac.sussex.gdsc.core.utils.LocalList)10 Plot (ij.gui.Plot)8 PeakResult (uk.ac.sussex.gdsc.smlm.results.PeakResult)8 PeakResultProcedure (uk.ac.sussex.gdsc.smlm.results.procedures.PeakResultProcedure)7 ImagePlus (ij.ImagePlus)6 ExecutorService (java.util.concurrent.ExecutorService)6 Future (java.util.concurrent.Future)6 HistogramPlotBuilder (uk.ac.sussex.gdsc.core.ij.HistogramPlot.HistogramPlotBuilder)6 ExtendedGenericDialog (uk.ac.sussex.gdsc.core.ij.gui.ExtendedGenericDialog)6 Ticker (uk.ac.sussex.gdsc.core.logging.Ticker)6 DistanceUnit (uk.ac.sussex.gdsc.smlm.data.config.UnitProtos.DistanceUnit)6 Trace (uk.ac.sussex.gdsc.smlm.results.Trace)6 TIntArrayList (gnu.trove.list.array.TIntArrayList)5 ImageStack (ij.ImageStack)5