Search in sources :

Example 1 with TurboList

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

the class CMOSAnalysis method simulate.

private void simulate() {
    // Create the offset, variance and gain for each pixel
    int n = size * size;
    float[] pixelOffset = new float[n];
    float[] pixelVariance = new float[n];
    float[] pixelGain = new float[n];
    IJ.showStatus("Creating random per-pixel readout");
    long start = System.currentTimeMillis();
    RandomGenerator rg = new Well19937c();
    PoissonDistribution pd = new PoissonDistribution(rg, offset, PoissonDistribution.DEFAULT_EPSILON, PoissonDistribution.DEFAULT_MAX_ITERATIONS);
    ExponentialDistribution ed = new ExponentialDistribution(rg, variance, ExponentialDistribution.DEFAULT_INVERSE_ABSOLUTE_ACCURACY);
    totalProgress = n;
    stepProgress = Utils.getProgressInterval(totalProgress);
    for (int i = 0; i < n; i++) {
        if (i % n == 0)
            IJ.showProgress(i, n);
        // Q. Should these be clipped to a sensible range?
        pixelOffset[i] = (float) pd.sample();
        pixelVariance[i] = (float) ed.sample();
        pixelGain[i] = (float) (gain + rg.nextGaussian() * gainSD);
    }
    IJ.showProgress(1);
    // Avoid all the file saves from updating the progress bar and status line
    Utils.setShowStatus(false);
    Utils.setShowProgress(false);
    JLabel statusLine = Utils.getStatusLine();
    progressBar = Utils.getProgressBar();
    // Save to the directory as a stack
    ImageStack simulationStack = new ImageStack(size, size);
    simulationStack.addSlice("Offset", pixelOffset);
    simulationStack.addSlice("Variance", pixelVariance);
    simulationStack.addSlice("Gain", pixelGain);
    simulationImp = new ImagePlus("PerPixel", simulationStack);
    // Only the info property is saved to the TIFF file
    simulationImp.setProperty("Info", String.format("Offset=%s; Variance=%s; Gain=%s +/- %s", Utils.rounded(offset), Utils.rounded(variance), Utils.rounded(gain), Utils.rounded(gainSD)));
    IJ.save(simulationImp, new File(directory, "perPixelSimulation.tif").getPath());
    // Create thread pool and workers
    ExecutorService executor = Executors.newFixedThreadPool(getThreads());
    TurboList<Future<?>> futures = new TurboList<Future<?>>(nThreads);
    // Simulate the zero exposure input.
    // Simulate 20 - 200 photon images.
    int[] photons = new int[] { 0, 20, 50, 100, 200 };
    totalProgress = photons.length * frames;
    stepProgress = Utils.getProgressInterval(totalProgress);
    progress = 0;
    progressBar.show(0);
    // For saving stacks
    int blockSize = 10;
    int nPerThread = (int) Math.ceil((double) frames / nThreads);
    // Convert to fit the block size
    nPerThread = (int) Math.ceil((double) nPerThread / blockSize) * blockSize;
    long seed = start;
    for (int p : photons) {
        statusLine.setText("Simulating " + Utils.pleural(p, "photon"));
        // Create the directory
        File out = new File(directory, String.format("photon%03d", p));
        if (!out.exists())
            out.mkdir();
        for (int from = 0; from < frames; ) {
            int to = Math.min(from + nPerThread, frames);
            futures.add(executor.submit(new SimulationWorker(seed++, out.getPath(), simulationStack, from, to, blockSize, p)));
            from = to;
        }
        // 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();
            }
        }
        futures.clear();
    }
    Utils.setShowStatus(true);
    Utils.setShowProgress(true);
    IJ.showProgress(1);
    executor.shutdown();
    Utils.log("Simulation time = " + Utils.timeToString(System.currentTimeMillis() - start));
}
Also used : PoissonDistribution(org.apache.commons.math3.distribution.PoissonDistribution) TurboList(gdsc.core.utils.TurboList) ImageStack(ij.ImageStack) ExponentialDistribution(org.apache.commons.math3.distribution.ExponentialDistribution) JLabel(javax.swing.JLabel) Well19937c(org.apache.commons.math3.random.Well19937c) ImagePlus(ij.ImagePlus) PseudoRandomGenerator(gdsc.core.utils.PseudoRandomGenerator) RandomGenerator(org.apache.commons.math3.random.RandomGenerator) ExecutorService(java.util.concurrent.ExecutorService) Future(java.util.concurrent.Future) File(java.io.File)

Example 2 with TurboList

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

the class ConfigurationTemplate method loadSelectedStandardTemplates.

private void loadSelectedStandardTemplates() {
    final TemplateResource[] templates = listTemplates(false, true);
    if (templates.length == 0)
        return;
    MultiDialog md = new MultiDialog("Select Templates", new MultiDialog.BaseItems() {

        public int size() {
            return templates.length;
        }

        public String getFormattedName(int i) {
            return templates[i].name;
        }
    });
    md.addSelected(selected);
    md.showDialog();
    if (md.wasCanceled())
        return;
    selected = md.getSelectedResults();
    if (selected.isEmpty())
        return;
    // Use list filtering to get the selected templates
    TurboList<TemplateResource> list = new TurboList<TemplateResource>(Arrays.asList(templates));
    list.removeIf(new SimplePredicate<TemplateResource>() {

        public boolean test(TemplateResource t) {
            return !(selected.contains(t.name));
        }
    });
    int count = loadTemplates(list.toArray(new TemplateResource[list.size()]));
    IJ.showMessage("Loaded " + Utils.pleural(count, "standard template"));
}
Also used : TurboList(gdsc.core.utils.TurboList) Point(java.awt.Point)

Example 3 with TurboList

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

the class EJMLLinearSolverTest method runInversionSpeedTest.

private void runInversionSpeedTest(int flags) {
    final Gaussian2DFunction f0 = GaussianFunctionFactory.create2D(1, 10, 10, flags, null);
    int n = f0.size();
    final double[] y = new double[n];
    final TurboList<DenseMatrix64F> aList = new TurboList<DenseMatrix64F>();
    double[] testbackground = new double[] { 0.2, 0.7 };
    double[] testsignal1 = new double[] { 30, 100, 300 };
    double[] testcx1 = new double[] { 4.9, 5.3 };
    double[] testcy1 = new double[] { 4.8, 5.2 };
    double[] testw1 = new double[] { 1.1, 1.2, 1.5 };
    int np = f0.getNumberOfGradients();
    GradientCalculator calc = GradientCalculatorFactory.newCalculator(np);
    final RandomDataGenerator rdg = new RandomDataGenerator(new Well19937c(30051977));
    //double lambda = 10;
    for (double background : testbackground) // Peak 1
    for (double signal1 : testsignal1) for (double cx1 : testcx1) for (double cy1 : testcy1) for (double w1 : testw1) {
        double[] p = new double[] { background, signal1, 0, cx1, cy1, w1, w1 };
        f0.initialise(p);
        f0.forEach(new ValueProcedure() {

            int i = 0;

            public void execute(double value) {
                // Poisson data 
                y[i++] = rdg.nextPoisson(value);
            }
        });
        double[][] alpha = new double[np][np];
        double[] beta = new double[np];
        //double ss = 
        calc.findLinearised(n, y, p, alpha, beta, f0);
        //System.out.printf("SS = %f\n", ss);
        // As per the LVM algorithm
        //for (int i = 0; i < np; i++)
        //	alpha[i][i] *= lambda;
        aList.add(EJMLLinearSolver.toA(alpha));
    }
    DenseMatrix64F[] a = aList.toArray(new DenseMatrix64F[aList.size()]);
    boolean[] ignore = new boolean[a.length];
    double[][] answer = new double[a.length][];
    int runs = 100000 / a.length;
    TimingService ts = new TimingService(runs);
    TurboList<InversionTimingTask> tasks = new TurboList<InversionTimingTask>();
    tasks.add(new PseudoInverseInversionTimingTask(a, ignore, answer));
    tasks.add(new CholeskyInversionTimingTask(a, ignore, answer));
    tasks.add(new LinearInversionTimingTask(a, ignore, answer));
    tasks.add(new CholeskyLDLTInversionTimingTask(a, ignore, answer));
    tasks.add(new DirectInversionInversionTimingTask(a, ignore, answer));
    tasks.add(new DiagonalDirectInversionInversionTimingTask(a, ignore, answer));
    for (InversionTimingTask task : tasks) if (!task.badSolver)
        ts.execute(task);
    ts.repeat();
    ts.report();
}
Also used : ValueProcedure(gdsc.smlm.function.ValueProcedure) TurboList(gdsc.core.utils.TurboList) RandomDataGenerator(org.apache.commons.math3.random.RandomDataGenerator) Well19937c(org.apache.commons.math3.random.Well19937c) DenseMatrix64F(org.ejml.data.DenseMatrix64F) Gaussian2DFunction(gdsc.smlm.function.gaussian.Gaussian2DFunction) GradientCalculator(gdsc.smlm.fitting.nonlinear.gradient.GradientCalculator) TimingService(gdsc.core.test.TimingService)

Example 4 with TurboList

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

the class ErfGaussian2DFunctionTest method functionIsFasterUsingForEach.

// Speed test forEach verses equivalent eval() function calls
@Test
public void functionIsFasterUsingForEach() {
    final ErfGaussian2DFunction f1 = (ErfGaussian2DFunction) this.f1;
    final TurboList<double[]> params = new TurboList<double[]>();
    for (double background : testbackground) // Peak 1
    for (double amplitude1 : testamplitude1) for (double shape1 : testshape1) for (double cx1 : testcx1) for (double cy1 : testcy1) for (double[] w1 : testw1) {
        double[] a = createParameters(background, amplitude1, shape1, cx1, cy1, w1[0], w1[1]);
        params.add(a);
    }
    double[][] x = params.toArray(new double[params.size()][]);
    int runs = 10000 / x.length;
    TimingService ts = new TimingService(runs);
    ts.execute(new FunctionTimingTask(f1, x, 2));
    ts.execute(new FunctionTimingTask(f1, x, 1));
    ts.execute(new FunctionTimingTask(f1, x, 0));
    ts.execute(new ForEachTimingTask(f1, x, 2));
    ts.execute(new ForEachTimingTask(f1, x, 1));
    ts.execute(new ForEachTimingTask(f1, x, 0));
    int size = ts.getSize();
    ts.repeat(size);
    ts.report();
    // in the more used eval(...) functions, so skip for now
    if (!TestSettings.ASSERT_SPEED_TESTS)
        return;
    int n = ts.getSize() - 1;
    Assert.assertTrue("0 order", ts.get(n).getMean() < ts.get(n - 3).getMean());
    n--;
    Assert.assertTrue("1 order", ts.get(n).getMean() < ts.get(n - 3).getMean());
    n--;
    Assert.assertTrue("2 order", ts.get(n).getMean() < ts.get(n - 3).getMean());
}
Also used : TurboList(gdsc.core.utils.TurboList) TimingService(gdsc.core.test.TimingService) Gaussian2DFunctionTest(gdsc.smlm.function.gaussian.Gaussian2DFunctionTest) Test(org.junit.Test)

Example 5 with TurboList

use of gdsc.core.utils.TurboList 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));
}
Also used : RollingArrayMoment(gdsc.core.math.RollingArrayMoment) ArrayMoment(gdsc.core.math.ArrayMoment) RollingArrayMoment(gdsc.core.math.RollingArrayMoment) SimpleArrayMoment(gdsc.core.math.SimpleArrayMoment) ArrayBlockingQueue(java.util.concurrent.ArrayBlockingQueue) SimpleArrayMoment(gdsc.core.math.SimpleArrayMoment) TurboList(gdsc.core.utils.TurboList) ImageStack(ij.ImageStack) SeriesImageSource(gdsc.smlm.ij.SeriesImageSource) JLabel(javax.swing.JLabel) WindowOrganiser(ij.plugin.WindowOrganiser) Statistics(gdsc.core.utils.Statistics) ImagePlus(ij.ImagePlus) ExecutorService(java.util.concurrent.ExecutorService) Future(java.util.concurrent.Future) File(java.io.File)

Aggregations

TurboList (gdsc.core.utils.TurboList)11 TimingService (gdsc.core.test.TimingService)4 Gaussian2DFunction (gdsc.smlm.function.gaussian.Gaussian2DFunction)3 ImagePlus (ij.ImagePlus)3 ExecutorService (java.util.concurrent.ExecutorService)3 Future (java.util.concurrent.Future)3 Statistics (gdsc.core.utils.Statistics)2 GradientCalculator (gdsc.smlm.fitting.nonlinear.gradient.GradientCalculator)2 ValueProcedure (gdsc.smlm.function.ValueProcedure)2 Gaussian2DFunctionTest (gdsc.smlm.function.gaussian.Gaussian2DFunctionTest)2 MemoryPeakResults (gdsc.smlm.results.MemoryPeakResults)2 PeakResult (gdsc.smlm.results.PeakResult)2 ImageStack (ij.ImageStack)2 ExtendedGenericDialog (ij.gui.ExtendedGenericDialog)2 WindowOrganiser (ij.plugin.WindowOrganiser)2 Rectangle (java.awt.Rectangle)2 File (java.io.File)2 JLabel (javax.swing.JLabel)2 Well19937c (org.apache.commons.math3.random.Well19937c)2 ArrayMoment (gdsc.core.math.ArrayMoment)1