Examples with Statistics gdsc.core.utils.Statistics used on opensource projects

Example 16 with Statistics

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

the class CMOSAnalysis method run.

private void run() {
    long start = System.currentTimeMillis();
    // Avoid all the file saves from updating the progress bar and status line
    Utils.setShowProgress(false);
    Utils.setShowStatus(false);
    JLabel statusLine = Utils.getStatusLine();
    progressBar = Utils.getProgressBar();
    // Create thread pool and workers
    ExecutorService executor = Executors.newFixedThreadPool(getThreads());
    TurboList<Future<?>> futures = new TurboList<Future<?>>(nThreads);
    TurboList<ImageWorker> workers = new TurboList<ImageWorker>(nThreads);
    double[][][] data = new double[subDirs.size()][2][];
    double[] pixelOffset = null, pixelVariance = null;
    Statistics statsOffset = null, statsVariance = null;
    // For each sub-directory compute the mean and variance
    final int nSubDirs = subDirs.size();
    boolean error = false;
    for (int n = 0; n < nSubDirs; n++) {
        SubDir sd = subDirs.getf(n);
        statusLine.setText("Analysing " + sd.name);
        // Open the series
        SeriesImageSource source = new SeriesImageSource(sd.name, sd.path.getPath());
        //source.setLogProgress(true);
        if (!source.open()) {
            error = true;
            IJ.error(TITLE, "Failed to open image series: " + sd.path.getPath());
            break;
        }
        // So the bar remains at 99% when workers have finished
        totalProgress = source.getFrames() + 1;
        stepProgress = Utils.getProgressInterval(totalProgress);
        progress = 0;
        progressBar.show(0);
        ArrayMoment moment = (rollingAlgorithm) ? new RollingArrayMoment() : new SimpleArrayMoment();
        final BlockingQueue<ImageJob> jobs = new ArrayBlockingQueue<ImageJob>(nThreads * 2);
        for (int i = 0; i < nThreads; i++) {
            final ImageWorker worker = new ImageWorker(jobs, moment);
            workers.add(worker);
            futures.add(executor.submit(worker));
        }
        // Process the data
        for (float[] pixels = source.next(); pixels != null; pixels = source.next()) {
            put(jobs, new ImageJob(pixels));
        }
        // Finish all the worker threads by passing in a null job
        for (int i = 0; i < nThreads; i++) {
            put(jobs, new ImageJob(null));
        }
        // Wait for all to finish
        for (int t = futures.size(); t-- > 0; ) {
            try {
                // The future .get() method will block until completed
                futures.get(t).get();
            } catch (Exception e) {
                // This should not happen. 
                e.printStackTrace();
            }
        }
        // Create the final aggregate statistics
        for (ImageWorker w : workers) moment.add(w.moment);
        data[n][0] = moment.getFirstMoment();
        data[n][1] = moment.getVariance();
        // Reset
        futures.clear();
        workers.clear();
        Statistics s = new Statistics(data[n][0]);
        if (n != 0) {
            // Compute mean ADU
            Statistics signal = new Statistics();
            double[] mean = data[n][0];
            for (int i = 0; i < pixelOffset.length; i++) signal.add(mean[i] - pixelOffset[i]);
            Utils.log("%s Mean = %s +/- %s. Signal = %s +/- %s ADU", sd.name, Utils.rounded(s.getMean()), Utils.rounded(s.getStandardDeviation()), Utils.rounded(signal.getMean()), Utils.rounded(signal.getStandardDeviation()));
            // TODO - optionally save
            ImageStack stack = new ImageStack(source.getWidth(), source.getHeight());
            stack.addSlice("Mean", Utils.toFloat(data[n][0]));
            stack.addSlice("Variance", Utils.toFloat(data[n][1]));
            IJ.save(new ImagePlus("PerPixel", stack), new File(directory, "perPixel" + sd.name + ".tif").getPath());
        } else {
            pixelOffset = data[0][0];
            pixelVariance = data[0][1];
            statsOffset = s;
            statsVariance = new Statistics(pixelVariance);
            Utils.log("%s Offset = %s +/- %s. Variance = %s +/- %s", sd.name, Utils.rounded(s.getMean()), Utils.rounded(s.getStandardDeviation()), Utils.rounded(statsVariance.getMean()), Utils.rounded(statsVariance.getStandardDeviation()));
        }
    }
    Utils.setShowStatus(true);
    Utils.setShowProgress(true);
    IJ.showProgress(1);
    executor.shutdown();
    if (error)
        return;
    // Compute the gain
    statusLine.setText("Computing gain");
    double[] pixelGain = new double[pixelOffset.length];
    // Ignore first as this is the 0 exposure image
    for (int i = 0; i < pixelGain.length; i++) {
        // Use equation 2.5 from the Huang et al paper.
        double bibiT = 0;
        double biaiT = 0;
        for (int n = 1; n < nSubDirs; n++) {
            double bi = data[n][0][i] - pixelOffset[i];
            double ai = data[n][1][i] - pixelVariance[i];
            bibiT += bi * bi;
            biaiT += bi * ai;
        }
        pixelGain[i] = biaiT / bibiT;
    }
    Statistics statsGain = new Statistics(pixelGain);
    Utils.log("Gain Mean = %s +/- %s", Utils.rounded(statsGain.getMean()), Utils.rounded(statsGain.getStandardDeviation()));
    // Histogram of offset, variance and gain
    int bins = Utils.getBinsSturges(pixelGain.length);
    WindowOrganiser wo = new WindowOrganiser();
    showHistogram("Offset", pixelOffset, bins, statsOffset, wo);
    showHistogram("Variance", pixelVariance, bins, statsVariance, wo);
    showHistogram("Gain", pixelGain, bins, statsGain, wo);
    wo.tile();
    // Save
    measuredStack = new ImageStack(size, size);
    measuredStack.addSlice("Offset", Utils.toFloat(pixelOffset));
    measuredStack.addSlice("Variance", Utils.toFloat(pixelVariance));
    measuredStack.addSlice("Gain", Utils.toFloat(pixelGain));
    IJ.save(new ImagePlus("PerPixel", measuredStack), new File(directory, "perPixel.tif").getPath());
    // Remove the status from the ij.io.ImageWriter class
    IJ.showStatus("");
    Utils.log("Analysis time = " + Utils.timeToString(System.currentTimeMillis() - start));
}

Example 17 with Statistics

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

the class CreateData method setNoise.

/**
	 * Sets the noise in the results if missing.
	 *
	 * @param results
	 *            the results
	 */
private void setNoise(MemoryPeakResults results, ImagePlus imp) {
    // Loaded results do not have noise
    for (PeakResult r : results.getResults()) if (r.noise != 0)
        return;
    // Compute noise per frame
    ImageStack stack = imp.getImageStack();
    final int width = stack.getWidth();
    final int height = stack.getHeight();
    final IJImageSource source = new IJImageSource(imp);
    final float[] noise = new float[source.getFrames() + 1];
    for (int slice = 1; slice < noise.length; slice++) {
        stack.getPixels(slice);
        float[] data = source.next();
        // Use the trimmed method as there may be a lot of spots in the frame
        noise[slice] = (float) FitWorker.estimateNoise(data, width, height, NoiseEstimator.Method.QUICK_RESIDUALS_LEAST_TRIMMED_OF_SQUARES);
    }
    Statistics stats = new Statistics(Arrays.copyOfRange(noise, 1, noise.length));
    System.out.printf("Noise = %.3f +/- %.3f (%d)\n", stats.getMean(), stats.getStandardDeviation(), stats.getN());
    for (PeakResult p : results.getResults()) {
        if (p.getFrame() < noise.length)
            p.noise = noise[p.getFrame()];
    }
}

Example 18 with Statistics

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

the class TraceDiffusion method summarise.

private void summarise(Trace[] traces, double[] fitMSDResult, int n, double[][] jdParams) {
    IJ.showStatus("Calculating summary ...");
    // Create summary table
    createSummaryTable();
    Statistics[] stats = new Statistics[NAMES.length];
    for (int i = 0; i < stats.length; i++) {
        stats[i] = (settings.showHistograms) ? new StoredDataStatistics() : new Statistics();
    }
    for (Trace trace : traces) {
        stats[T_ON].add(trace.getOnTime() * exposureTime);
        final double signal = trace.getSignal() / results.getGain();
        stats[TOTAL_SIGNAL].add(signal);
        stats[SIGNAL_PER_FRAME].add(signal / trace.size());
    }
    // Add to the summary table
    StringBuilder sb = new StringBuilder(title);
    sb.append('\t').append(createCombinedName());
    sb.append("\t");
    sb.append(Utils.rounded(exposureTime * 1000, 3)).append("\t");
    sb.append(Utils.rounded(settings.distanceThreshold, 3)).append("\t");
    sb.append(Utils.rounded(settings.distanceExclusion, 3)).append("\t");
    sb.append(settings.minimumTraceLength).append("\t");
    sb.append(settings.ignoreEnds).append("\t");
    sb.append(settings.truncate).append("\t");
    sb.append(settings.internalDistances).append("\t");
    sb.append(settings.fitLength).append("\t");
    sb.append(settings.msdCorrection).append("\t");
    sb.append(settings.precisionCorrection).append("\t");
    sb.append(settings.mle).append("\t");
    sb.append(traces.length).append("\t");
    sb.append(Utils.rounded(precision, 4)).append("\t");
    double D = 0, s = 0;
    if (fitMSDResult != null) {
        D = fitMSDResult[0];
        s = fitMSDResult[1];
    }
    sb.append(Utils.rounded(D, 4)).append("\t");
    sb.append(Utils.rounded(s * 1000, 4)).append("\t");
    sb.append(Utils.rounded(settings.jumpDistance * exposureTime)).append("\t");
    sb.append(n).append("\t");
    sb.append(Utils.rounded(beta, 4)).append("\t");
    if (jdParams == null) {
        sb.append("\t\t\t");
    } else {
        sb.append(format(jdParams[0])).append("\t");
        sb.append(format(jdParams[1])).append("\t");
        sb.append(Utils.rounded(ic)).append("\t");
    }
    for (int i = 0; i < stats.length; i++) {
        sb.append(Utils.rounded(stats[i].getMean(), 3)).append("\t");
    }
    if (java.awt.GraphicsEnvironment.isHeadless()) {
        IJ.log(sb.toString());
        return;
    } else {
        summaryTable.append(sb.toString());
    }
    if (settings.showHistograms) {
        IJ.showStatus("Calculating histograms ...");
        for (int i = 0; i < NAMES.length; i++) {
            if (displayHistograms[i]) {
                showHistogram((StoredDataStatistics) stats[i], NAMES[i], alwaysRemoveOutliers[i], ROUNDED[i], false);
            }
        }
    }
    tileNewWindows();
    IJ.showStatus("Finished " + TITLE);
}

Example 19 with Statistics

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

the class SpotAnalysis method extractSpotProfile.

private double[][] extractSpotProfile(ImagePlus imp, Rectangle bounds, ImageStack rawSpot) {
    final int nSlices = imp.getStackSize();
    IJImageSource rawSource = new IJImageSource(imp);
    double[][] profile = new double[2][nSlices];
    for (int n = 0; n < nSlices; n++) {
        IJ.showProgress(n, nSlices);
        float[] data = rawSource.next(bounds);
        rawSpot.setPixels(data, n + 1);
        Statistics stats = new Statistics(data);
        profile[0][n] = stats.getMean() / gain;
        profile[1][n] = stats.getStandardDeviation() / gain;
    }
    return profile;
}

Example 20 with Statistics

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

the class SpotAnalysis method saveSpot.

private void saveSpot() {
    if (onFrames.isEmpty() || !updated)
        return;
    createResultsWindow();
    id++;
    double signal = 0;
    Trace trace = null;
    float psfWidth = Float.parseFloat(widthTextField.getText());
    float cx = areaBounds.x + areaBounds.width / 2.0f;
    float cy = areaBounds.y + areaBounds.height / 2.0f;
    for (Spot s : onFrames) {
        // Get the signal again since the background may have changed.
        final double spotSignal = getSignal(s.frame);
        signal += spotSignal;
        //signal += s.signal;
        float[] params = new float[7];
        params[Gaussian2DFunction.X_POSITION] = cx;
        params[Gaussian2DFunction.Y_POSITION] = cy;
        params[Gaussian2DFunction.SIGNAL] = (float) (spotSignal);
        params[Gaussian2DFunction.X_SD] = params[Gaussian2DFunction.Y_SD] = psfWidth;
        PeakResult result = new PeakResult(s.frame, (int) cx, (int) cy, 0, 0, 0, params, null);
        if (trace == null)
            trace = new Trace(result);
        else
            trace.add(result);
    }
    Statistics tOn = new Statistics(trace.getOnTimes());
    Statistics tOff = new Statistics(trace.getOffTimes());
    resultsWindow.append(String.format("%d\t%.1f\t%.1f\t%s\t%s\t%s\t%d\t%s\t%s\t%s", id, cx, cy, Utils.rounded(signal, 4), Utils.rounded(tOn.getSum() * msPerFrame, 3), Utils.rounded(tOff.getSum() * msPerFrame, 3), trace.getNBlinks() - 1, Utils.rounded(tOn.getMean() * msPerFrame, 3), Utils.rounded(tOff.getMean() * msPerFrame, 3), imp.getTitle()));
    // Save the individual on/off times for use in creating a histogram
    traces.put(id, trace);
    updated = false;
//clearSelectedFrames();
}