Examples with Statistics uk.ac.sussex.gdsc.core.utils.Statistics used on opensource projects

Example 21 with Statistics

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

the class CubicSplineFunctionTest method functionComputesTargetGradient1.

private void functionComputesTargetGradient1(int targetParameter) {
    final int gradientIndex = findGradientIndex(f1, targetParameter);
    final Statistics s = new Statistics();
    final StandardValueProcedure p1a = new StandardValueProcedure();
    final StandardValueProcedure p1b = new StandardValueProcedure();
    final StandardGradient1Procedure p2 = new StandardGradient1Procedure();
    for (final double background : testbackground) {
        // Peak 1
        for (final double signal1 : testsignal1) {
            for (final double cx1 : testcx1) {
                for (final double cy1 : testcy1) {
                    for (final double cz1 : testcz1) {
                        final double[] a = createParameters(background, signal1, cx1, cy1, cz1);
                        // System.out.println(java.util.Arrays.toString(a));
                        // Evaluate all gradients
                        p2.getValues(f1, a);
                        // Numerically solve gradient.
                        // Calculate the step size h to be an exact numerical representation
                        final double xx = a[targetParameter];
                        // Get h to minimise roundoff error
                        final double h = Precision.representableDelta(xx, stepH);
                        // Evaluate at (x+h) and (x-h)
                        a[targetParameter] = xx + h;
                        p1a.getValues(f1, a);
                        a[targetParameter] = xx - h;
                        p1b.getValues(f1, a);
                        // Only test close to the XY centre
                        for (final int x : testx) {
                            for (final int y : testy) {
                                final int i = y * maxx + x;
                                final double high = p1a.values[i];
                                final double low = p1b.values[i];
                                final double gradient = (high - low) / (2 * h);
                                final double dyda = p2.gradients[i][gradientIndex];
                                final double error = DoubleEquality.relativeError(gradient, dyda);
                                s.add(error);
                                if ((gradient * dyda) < 0) {
                                    Assertions.fail(String.format("%s sign != %s", gradient, dyda));
                                }
                                // gradient, gradientIndex, dyda, error);
                                if (!eq.almostEqualRelativeOrAbsolute(gradient, dyda)) {
                                    Assertions.fail(String.format("%s != %s", gradient, dyda));
                                }
                            }
                        }
                    }
                }
            }
        }
    }
    logger.info(() -> {
        return String.format("functionComputesTargetGradient1 %s %s (error %s +/- %s)", f1.getClass().getSimpleName(), CubicSplineFunction.getName(targetParameter), MathUtils.rounded(s.getMean()), MathUtils.rounded(s.getStandardDeviation()));
    });
}

Example 22 with Statistics

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

the class CubicSplineFunctionTest method functionComputesTargetGradient2With2Peaks.

private void functionComputesTargetGradient2With2Peaks(int targetParameter) {
    final int gradientIndex = findGradientIndex(f2, targetParameter);
    final Statistics s = new Statistics();
    final StandardGradient1Procedure p1a = new StandardGradient1Procedure();
    final StandardGradient1Procedure p1b = new StandardGradient1Procedure();
    final StandardGradient2Procedure p2 = new StandardGradient2Procedure();
    for (final double background : testbackground) {
        // Peak 1
        for (final double signal1 : testsignal1) {
            for (final double cx1 : testcx1) {
                for (final double cy1 : testcy1) {
                    for (final double cz1 : testcz1) {
                        // Peak 2
                        for (final double signal2 : testsignal2) {
                            for (final double cx2 : testcx2) {
                                for (final double cy2 : testcy2) {
                                    for (final double cz2 : testcz2) {
                                        final double[] a = createParameters(background, signal1, cx1, cy1, cz1, signal2, cx2, cy2, cz2);
                                        // System.out.println(java.util.Arrays.toString(a));
                                        f2.initialise2(a);
                                        final boolean test = !f2.isNodeBoundary(gradientIndex);
                                        // Comment out when printing errors
                                        if (!test) {
                                            continue;
                                        }
                                        // Evaluate all gradients
                                        p2.getValues(f2, a);
                                        // Numerically solve gradient.
                                        // Calculate the step size h to be an exact numerical representation
                                        final double xx = a[targetParameter];
                                        // Get h to minimise roundoff error
                                        final double h = Precision.representableDelta(xx, stepH);
                                        // Evaluate at (x+h) and (x-h)
                                        a[targetParameter] = xx + h;
                                        p1a.getValues(f2, a);
                                        a[targetParameter] = xx - h;
                                        p1b.getValues(f2, a);
                                        // Only test close to the XY centre
                                        for (final int x : testx) {
                                            for (final int y : testy) {
                                                final int i = y * maxx + x;
                                                final double high = p1a.gradients[i][gradientIndex];
                                                final double low = p1b.gradients[i][gradientIndex];
                                                final double gradient = (high - low) / (2 * h);
                                                final double d2yda2 = p2.gradients2[i][gradientIndex];
                                                final double error = DoubleEquality.relativeError(gradient, d2yda2);
                                                // y, gradient, gradientIndex, d2yda2, error);
                                                if (test) {
                                                    s.add(error);
                                                    if ((gradient * d2yda2) < 0) {
                                                        Assertions.fail(String.format("%s sign != %s", gradient, d2yda2));
                                                    }
                                                    // x, y, gradient, gradientIndex, d2yda2, error);
                                                    if (!eq.almostEqualRelativeOrAbsolute(gradient, d2yda2)) {
                                                        Assertions.fail(String.format("%s != %s", gradient, d2yda2));
                                                    }
                                                }
                                            }
                                        }
                                    }
                                }
                            }
                        }
                    }
                }
            }
        }
    }
    logger.info(() -> {
        return String.format("functionComputesTargetGradient2With2Peaks %s %s (error %s +/- %s)", f1.getClass().getSimpleName(), CubicSplineFunction.getName(targetParameter), MathUtils.rounded(s.getMean()), MathUtils.rounded(s.getStandardDeviation()));
    });
}

Example 23 with Statistics

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

the class Gaussian2DFunctionTest method functionComputesTargetGradientWith2Peaks.

private void functionComputesTargetGradientWith2Peaks(int targetParameter) {
    final int gradientIndex = findGradientIndex(f2, targetParameter);
    final double[] dyda = new double[f2.gradientIndices().length];
    final double[] dyda2 = new double[dyda.length];
    double[] params;
    final Gaussian2DFunction f2a = GaussianFunctionFactory.create2D(2, maxx, maxy, flags, zModel);
    final Gaussian2DFunction f2b = GaussianFunctionFactory.create2D(2, maxx, maxy, flags, zModel);
    final Statistics s = new Statistics();
    for (final double background : testbackground) {
        // Peak 1
        for (final double signal1 : testsignal1) {
            for (final double cx1 : testcx1) {
                for (final double cy1 : testcy1) {
                    for (final double cz1 : testcz1) {
                        for (final double[] w1 : testw1) {
                            for (final double angle1 : testangle1) {
                                // Peak 2
                                for (final double signal2 : testsignal2) {
                                    for (final double cx2 : testcx2) {
                                        for (final double cy2 : testcy2) {
                                            for (final double cz2 : testcz2) {
                                                for (final double[] w2 : testw2) {
                                                    for (final double angle2 : testangle2) {
                                                        params = createParameters(background, signal1, cx1, cy1, cz1, w1[0], w1[1], angle1, signal2, cx2, cy2, cz2, w2[0], w2[1], angle2);
                                                        f2.initialise(params);
                                                        // Numerically solve gradient.
                                                        // Calculate the step size h to be an exact numerical
                                                        // representation
                                                        final double xx = params[targetParameter];
                                                        // Get h to minimise roundoff error
                                                        final double h = Precision.representableDelta(xx, stepH);
                                                        // Evaluate at (x+h) and (x-h)
                                                        params[targetParameter] = xx + h;
                                                        f2a.initialise(params.clone());
                                                        params[targetParameter] = xx - h;
                                                        f2b.initialise(params.clone());
                                                        for (final int x : testx) {
                                                            for (final int y : testy) {
                                                                final int i = y * maxx + x;
                                                                f2.eval(i, dyda);
                                                                final double value2 = f2a.eval(i, dyda2);
                                                                final double value3 = f2b.eval(i, dyda2);
                                                                final double gradient = (value2 - value3) / (2 * h);
                                                                final double error = DoubleEquality.relativeError(gradient, dyda2[gradientIndex]);
                                                                s.add(error);
                                                                if ((gradient * dyda2[gradientIndex]) < 0) {
                                                                    Assertions.fail(String.format("%s sign != %s", gradient, dyda2[gradientIndex]));
                                                                }
                                                                // dyda[gradientIndex]);
                                                                if (!eq.almostEqualRelativeOrAbsolute(gradient, dyda[gradientIndex])) {
                                                                    Assertions.fail(String.format("%s != %s", gradient, dyda[gradientIndex]));
                                                                }
                                                            }
                                                        }
                                                    }
                                                }
                                            }
                                        }
                                    }
                                }
                            }
                        }
                    }
                }
            }
        }
    }
    logger.info(() -> {
        return String.format("functionComputesTargetGradientWith2Peaks %s [%d] %s (error %s +/- %s)", f2.getClass().getSimpleName(), Gaussian2DFunction.getPeak(targetParameter), Gaussian2DFunction.getName(targetParameter), MathUtils.rounded(s.getMean()), MathUtils.rounded(s.getStandardDeviation()));
    });
}

Example 24 with Statistics

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

the class CmosAnalysis method runAnalysis.

private void runAnalysis() {
    final long start = System.currentTimeMillis();
    // Create thread pool and workers. The system is likely to be IO limited
    // so reduce the computation threads to allow the reading thread in the
    // SeriesImageSource to run.
    // If the images are small enough to fit into memory then 3 threads are used,
    // otherwise it is 1.
    final int nThreads = Math.max(1, getThreads() - 3);
    final ExecutorService executor = Executors.newFixedThreadPool(nThreads);
    final LocalList<Future<?>> futures = new LocalList<>(nThreads);
    final LocalList<ImageWorker> workers = new LocalList<>(nThreads);
    final double[][] data = new double[subDirs.size() * 2][];
    double[] pixelOffset = null;
    double[] pixelVariance = null;
    Statistics statsOffset = null;
    Statistics statsVariance = null;
    // For each sub-directory compute the mean and variance
    final int nSubDirs = subDirs.size();
    boolean error = false;
    int width = 0;
    int height = 0;
    for (int n = 0; n < nSubDirs; n++) {
        ImageJUtils.showSlowProgress(n, nSubDirs);
        final SubDir sd = subDirs.unsafeGet(n);
        ImageJUtils.showStatus(() -> "Analysing " + sd.name);
        final StopWatch sw = StopWatch.createStarted();
        // Option to reuse data
        final File file = new File(settings.directory, "perPixel" + sd.name + ".tif");
        boolean found = false;
        if (settings.reuseProcessedData && file.exists()) {
            final Opener opener = new Opener();
            opener.setSilentMode(true);
            final ImagePlus imp = opener.openImage(file.getPath());
            if (imp != null && imp.getStackSize() == 2 && imp.getBitDepth() == 32) {
                if (n == 0) {
                    width = imp.getWidth();
                    height = imp.getHeight();
                } else if (width != imp.getWidth() || height != imp.getHeight()) {
                    error = true;
                    IJ.error(TITLE, "Image width/height mismatch in image series: " + file.getPath() + String.format("\n \nExpected %dx%d, Found %dx%d", width, height, imp.getWidth(), imp.getHeight()));
                    break;
                }
                final ImageStack stack = imp.getImageStack();
                data[2 * n] = SimpleArrayUtils.toDouble((float[]) stack.getPixels(1));
                data[2 * n + 1] = SimpleArrayUtils.toDouble((float[]) stack.getPixels(2));
                found = true;
            }
        }
        if (!found) {
            // Open the series
            final SeriesImageSource source = new SeriesImageSource(sd.name, sd.path.getPath());
            if (!source.open()) {
                error = true;
                IJ.error(TITLE, "Failed to open image series: " + sd.path.getPath());
                break;
            }
            if (n == 0) {
                width = source.getWidth();
                height = source.getHeight();
            } else if (width != source.getWidth() || height != source.getHeight()) {
                error = true;
                IJ.error(TITLE, "Image width/height mismatch in image series: " + sd.path.getPath() + String.format("\n \nExpected %dx%d, Found %dx%d", width, height, source.getWidth(), source.getHeight()));
                break;
            }
            // So the bar remains at 99% when workers have finished use frames + 1
            final Ticker ticker = ImageJUtils.createTicker(source.getFrames() + 1L, nThreads);
            // Open the first frame to get the bit depth.
            // Assume the first pixels are not empty as the source is open.
            Object pixels = source.nextRaw();
            final int bitDepth = ImageJUtils.getBitDepth(pixels);
            ArrayMoment moment;
            if (settings.rollingAlgorithm) {
                moment = new RollingArrayMoment();
            // We assume 16-bit camera at the maximum
            } else if (bitDepth <= 16 && IntegerArrayMoment.isValid(IntegerType.UNSIGNED_16, source.getFrames())) {
                moment = new IntegerArrayMoment();
            } else {
                moment = new SimpleArrayMoment();
            }
            final BlockingQueue<Object> jobs = new ArrayBlockingQueue<>(nThreads * 2);
            for (int i = 0; i < nThreads; i++) {
                final ImageWorker worker = new ImageWorker(ticker, jobs, moment);
                workers.add(worker);
                futures.add(executor.submit(worker));
            }
            // Process the raw pixel data
            long lastTime = 0;
            while (pixels != null) {
                final long time = System.currentTimeMillis();
                if (time - lastTime > 150) {
                    if (ImageJUtils.isInterrupted()) {
                        error = true;
                        break;
                    }
                    lastTime = time;
                    IJ.showStatus("Analysing " + sd.name + " Frame " + source.getStartFrameNumber());
                }
                put(jobs, pixels);
                pixels = source.nextRaw();
            }
            source.close();
            if (error) {
                // Kill the workers
                workers.stream().forEach(worker -> worker.finished = true);
                // Clear the queue
                jobs.clear();
                // Signal any waiting workers
                workers.stream().forEach(worker -> jobs.add(ImageWorker.STOP_SIGNAL));
                // Cancel by interruption. We set the finished flag so the ImageWorker should
                // ignore the interrupt.
                futures.stream().forEach(future -> future.cancel(true));
                break;
            }
            // Finish all the worker threads cleanly
            workers.stream().forEach(worker -> jobs.add(ImageWorker.STOP_SIGNAL));
            // Wait for all to finish
            ConcurrencyUtils.waitForCompletionUnchecked(futures);
            // Create the final aggregate statistics
            for (final ImageWorker w : workers) {
                moment.add(w.moment);
            }
            data[2 * n] = moment.getMean();
            data[2 * n + 1] = moment.getVariance();
            // Get the processing speed.
            sw.stop();
            // ticker holds the number of number of frames processed
            final double bits = (double) bitDepth * source.getFrames() * source.getWidth() * source.getHeight();
            final double bps = bits / sw.getTime(TimeUnit.SECONDS);
            final SiPrefix prefix = SiPrefix.getSiPrefix(bps);
            ImageJUtils.log("Processed %d frames. Time = %s. Rate = %s %sbits/s", moment.getN(), sw.toString(), MathUtils.rounded(prefix.convert(bps)), prefix.getPrefix());
            // Reset
            futures.clear();
            workers.clear();
            final ImageStack stack = new ImageStack(width, height);
            stack.addSlice("Mean", SimpleArrayUtils.toFloat(data[2 * n]));
            stack.addSlice("Variance", SimpleArrayUtils.toFloat(data[2 * n + 1]));
            IJ.save(new ImagePlus("PerPixel", stack), file.getPath());
        }
        final Statistics s = Statistics.create(data[2 * n]);
        if (pixelOffset != null) {
            // Compute mean ADU
            final Statistics signal = new Statistics();
            final double[] mean = data[2 * n];
            for (int i = 0; i < pixelOffset.length; i++) {
                signal.add(mean[i] - pixelOffset[i]);
            }
            ImageJUtils.log("%s Mean = %s +/- %s. Signal = %s +/- %s ADU", sd.name, MathUtils.rounded(s.getMean()), MathUtils.rounded(s.getStandardDeviation()), MathUtils.rounded(signal.getMean()), MathUtils.rounded(signal.getStandardDeviation()));
        } else {
            // Set the offset assuming the first sub-directory is the bias image
            pixelOffset = data[0];
            pixelVariance = data[1];
            statsOffset = s;
            statsVariance = Statistics.create(pixelVariance);
            ImageJUtils.log("%s Offset = %s +/- %s. Variance = %s +/- %s", sd.name, MathUtils.rounded(s.getMean()), MathUtils.rounded(s.getStandardDeviation()), MathUtils.rounded(statsVariance.getMean()), MathUtils.rounded(statsVariance.getStandardDeviation()));
        }
        IJ.showProgress(1);
    }
    ImageJUtils.clearSlowProgress();
    if (error) {
        executor.shutdownNow();
        IJ.showStatus(TITLE + " cancelled");
        return;
    }
    executor.shutdown();
    if (pixelOffset == null || pixelVariance == null) {
        IJ.showStatus(TITLE + " error: no bias image");
        return;
    }
    // Compute the gain
    ImageJUtils.showStatus("Computing gain");
    final double[] pixelGain = new double[pixelOffset.length];
    final double[] bibiT = new double[pixelGain.length];
    final double[] biaiT = new double[pixelGain.length];
    // Ignore first as this is the 0 exposure image
    for (int n = 1; n < nSubDirs; n++) {
        // Use equation 2.5 from the Huang et al paper.
        final double[] b = data[2 * n];
        final double[] a = data[2 * n + 1];
        for (int i = 0; i < pixelGain.length; i++) {
            final double bi = b[i] - pixelOffset[i];
            final double ai = a[i] - pixelVariance[i];
            bibiT[i] += bi * bi;
            biaiT[i] += bi * ai;
        }
    }
    for (int i = 0; i < pixelGain.length; i++) {
        pixelGain[i] = biaiT[i] / bibiT[i];
    }
    final Statistics statsGain = Statistics.create(pixelGain);
    ImageJUtils.log("Gain Mean = %s +/- %s", MathUtils.rounded(statsGain.getMean()), MathUtils.rounded(statsGain.getStandardDeviation()));
    // Histogram of offset, variance and gain
    final int bins = 2 * HistogramPlot.getBinsSturgesRule(pixelGain.length);
    final WindowOrganiser wo = new WindowOrganiser();
    showHistogram("Offset (ADU)", pixelOffset, bins, statsOffset, wo);
    showHistogram("Variance (ADU^2)", pixelVariance, bins, statsVariance, wo);
    showHistogram("Gain (ADU/e)", pixelGain, bins, statsGain, wo);
    wo.tile();
    // Save
    final float[] bias = SimpleArrayUtils.toFloat(pixelOffset);
    final float[] variance = SimpleArrayUtils.toFloat(pixelVariance);
    final float[] gain = SimpleArrayUtils.toFloat(pixelGain);
    measuredStack = new ImageStack(width, height);
    measuredStack.addSlice("Offset", bias);
    measuredStack.addSlice("Variance", variance);
    measuredStack.addSlice("Gain", gain);
    final ExtendedGenericDialog egd = new ExtendedGenericDialog(TITLE);
    egd.addMessage("Save the sCMOS camera model?");
    if (settings.modelDirectory == null) {
        settings.modelDirectory = settings.directory;
        settings.modelName = "sCMOS Camera";
    }
    egd.addStringField("Model_name", settings.modelName, 30);
    egd.addDirectoryField("Model_directory", settings.modelDirectory);
    egd.showDialog();
    if (!egd.wasCanceled()) {
        settings.modelName = egd.getNextString();
        settings.modelDirectory = egd.getNextString();
        saveCameraModel(width, height, bias, gain, variance);
    }
    // Remove the status from the ij.io.ImageWriter class
    IJ.showStatus("");
    ImageJUtils.log("Analysis time = " + TextUtils.millisToString(System.currentTimeMillis() - start));
}

Example 25 with Statistics

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

the class BackgroundEstimator method plot.

private void plot(WindowOrganiser wo, double[] xvalues, double[] data1, double[] data2, double[] data3, String title, String title1, String title2, String title3) {
    // Get limits
    final double[] xlimits = MathUtils.limits(xvalues);
    double[] ylimits = MathUtils.limits(data1);
    ylimits = MathUtils.limits(ylimits, data2);
    if (data3 != null) {
        ylimits = MathUtils.limits(ylimits, data3);
    }
    title = imp.getTitle() + " " + title;
    final Plot plot = new Plot(title, "Slice", title);
    double range = ylimits[1] - ylimits[0];
    if (range == 0) {
        range = 1;
    }
    plot.setLimits(xlimits[0], xlimits[1], ylimits[0] - 0.05 * range, ylimits[1] + 0.05 * range);
    plot.setColor(Color.blue);
    plot.addPoints(xvalues, data1, Plot.LINE);
    plot.draw();
    Statistics stats = Statistics.create(data1);
    @SuppressWarnings("resource") final Formatter label = new Formatter().format("%s (Blue) = %s +/- %s", title1, MathUtils.rounded(stats.getMean()), MathUtils.rounded(stats.getStandardDeviation()));
    plot.setColor(Color.red);
    plot.addPoints(xvalues, data2, Plot.LINE);
    stats = Statistics.create(data2);
    label.format(", %s (Red) = %s +/- %s", title2, MathUtils.rounded(stats.getMean()), MathUtils.rounded(stats.getStandardDeviation()));
    if (data3 != null) {
        plot.setColor(Color.green);
        plot.addPoints(xvalues, data3, Plot.LINE);
        stats = Statistics.create(data3);
        label.format(", %s (Green) = %s +/- %s", title3, MathUtils.rounded(stats.getMean()), MathUtils.rounded(stats.getStandardDeviation()));
    }
    plot.setColor(Color.black);
    plot.addLabel(0, 0, label.toString());
    ImageJUtils.display(title, plot, wo);
}