Examples with SimpleArrayMoment gdsc.core.math.SimpleArrayMoment used on opensource projects

Example 1 with SimpleArrayMoment

use of gdsc.core.math.SimpleArrayMoment 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 2 with SimpleArrayMoment

use of gdsc.core.math.SimpleArrayMoment in project GDSC-SMLM by aherbert.

the class BaseFunctionSolverTest method canFitSingleGaussian.

void canFitSingleGaussian(FunctionSolver solver, boolean applyBounds, NoiseModel noiseModel) {
    // Allow reporting the fit deviations
    boolean report = false;
    double[] crlb = null;
    SimpleArrayMoment m = null;
    double[] noise = getNoise(noiseModel);
    if (solver.isWeighted())
        solver.setWeights(getWeights(noiseModel));
    randomGenerator.setSeed(seed);
    for (double s : signal) {
        double[] expected = createParams(1, s, 0, 0, 1);
        double[] lower = createParams(0, s * 0.5, -0.2, -0.2, 0.8);
        double[] upper = createParams(3, s * 2, 0.2, 0.2, 1.2);
        if (applyBounds)
            solver.setBounds(lower, upper);
        if (report) {
            // Compute the CRLB for a Poisson process
            PoissonGradientProcedure gp = PoissonGradientProcedureFactory.create((Gradient1Function) ((BaseFunctionSolver) solver).getGradientFunction());
            gp.computeFisherInformation(expected);
            FisherInformationMatrix f = new FisherInformationMatrix(gp.getLinear(), gp.n);
            crlb = f.crlbSqrt();
            // Compute the deviations
            m = new SimpleArrayMoment();
        }
        double[] data = drawGaussian(expected, noise, noiseModel);
        for (double db : base) for (double dx : shift) for (double dy : shift) for (double dsx : factor) {
            double[] p = createParams(db, s, dx, dy, dsx);
            double[] fp = fitGaussian(solver, data, p, expected);
            for (int i = 0; i < expected.length; i++) {
                if (fp[i] < lower[i])
                    Assert.assertTrue(String.format("Fit Failed: [%d] %.2f < %.2f: %s != %s", i, fp[i], lower[i], Arrays.toString(fp), Arrays.toString(expected)), false);
                if (fp[i] > upper[i])
                    Assert.assertTrue(String.format("Fit Failed: [%d] %.2f > %.2f: %s != %s", i, fp[i], upper[i], Arrays.toString(fp), Arrays.toString(expected)), false);
                if (report)
                    fp[i] = expected[i] - fp[i];
            }
            // Store the deviations
            if (report)
                m.add(fp);
        }
        // Report
        if (report)
            System.out.printf("%s %s %f : CRLB = %s, Devaitions = %s\n", solver.getClass().getSimpleName(), noiseModel, s, Arrays.toString(crlb), Arrays.toString(m.getStandardDeviation()));
    }
}