Search in sources :

Example 31 with Ticker

use of uk.ac.sussex.gdsc.core.logging.Ticker in project GDSC-SMLM by aherbert.

the class DriftCalculator method calculateDriftUsingFrames.

/**
 * Calculate the drift by aligning N consecutive frames with the overall image. Update the current
 * drift parameters.
 *
 * @param blocks the blocks
 * @param blockT the block T
 * @param bounds the bounds
 * @param scale the scale
 * @param dx the dx
 * @param dy the dy
 * @param originalDriftTimePoints the original drift time points
 * @param smoothing the smoothing
 * @param iterations the iterations
 * @return the double
 */
private double calculateDriftUsingFrames(ArrayList<ArrayList<Localisation>> blocks, int[] blockT, Rectangle bounds, float scale, double[] dx, double[] dy, double[] originalDriftTimePoints, double smoothing, int iterations) {
    // Construct images using the current drift
    tracker.status("Constructing images");
    // Build an image for each block of results.
    final ImageProcessor[] images = new ImageProcessor[blocks.size()];
    final List<Future<?>> futures = new LinkedList<>();
    final Ticker ticker = Ticker.createStarted(tracker, images.length, true);
    for (int i = 0; i < images.length; i++) {
        futures.add(executor.submit(new ImageBuilder(blocks.get(i), images, i, bounds, scale, dx, dy, ticker)));
    }
    ConcurrencyUtils.waitForCompletionUnchecked(futures);
    // Build an image with all results.
    final FloatProcessor allIp = new FloatProcessor(images[0].getWidth(), images[0].getHeight());
    for (final ImageProcessor ip : images) {
        allIp.copyBits(ip, 0, 0, Blitter.ADD);
    }
    return calculateDrift(blockT, scale, dx, dy, originalDriftTimePoints, smoothing, iterations, images, allIp, true, ticker);
}
Also used : ImageProcessor(ij.process.ImageProcessor) FloatProcessor(ij.process.FloatProcessor) Ticker(uk.ac.sussex.gdsc.core.logging.Ticker) Future(java.util.concurrent.Future) LinkedList(java.util.LinkedList) Point(java.awt.Point)

Example 32 with Ticker

use of uk.ac.sussex.gdsc.core.logging.Ticker in project GDSC-SMLM by aherbert.

the class EmGainAnalysis method simulateFromPdf.

/**
 * Random sample from the cumulative probability distribution function that is used during
 * fitting.
 *
 * @return The histogram
 */
private int[] simulateFromPdf() {
    final double step = getStepSize(settings.settingPhotons, settings.settingGain, settings.settingNoise);
    final Pdf pdf = pdf(0, step, settings.settingPhotons, settings.settingGain, settings.settingNoise);
    // Debug this
    final double[] g = pdf.probability;
    final double[] x = pdf.x;
    final Plot plot = new Plot(TITLE + " PDF", "ADU", "P");
    plot.addPoints(x, Arrays.copyOf(g, g.length), Plot.LINE);
    ImageJUtils.display(TITLE + " PDF", plot);
    // Get cumulative probability
    double sum = 0;
    for (int i = 0; i < g.length; i++) {
        final double p = g[i];
        g[i] += sum;
        sum += p;
    }
    for (int i = 0; i < g.length; i++) {
        g[i] /= sum;
    }
    // Ensure value of 1 at the end
    g[g.length - 1] = 1;
    // Randomly sample
    final UniformRandomProvider rng = UniformRandomProviders.create();
    final int bias = (int) settings.settingBias;
    final int[] bins = new int[x.length];
    for (int i = 0; i < x.length; i++) {
        bins[i] = bias + (int) x[i];
    }
    final int[] h = new int[bins[bins.length - 1] + 1];
    final int steps = settings.simulationSize;
    final Ticker ticker = ImageJUtils.createTicker(steps, 1);
    for (int n = 0; n < steps; n++) {
        int index = binarySearch(g, rng.nextDouble());
        if (index < 0) {
            index = -(index + 1);
        }
        h[bins[index]]++;
        ticker.tick();
    }
    return h;
}
Also used : Plot(ij.gui.Plot) Ticker(uk.ac.sussex.gdsc.core.logging.Ticker) UniformRandomProvider(org.apache.commons.rng.UniformRandomProvider) Point(java.awt.Point)

Example 33 with Ticker

use of uk.ac.sussex.gdsc.core.logging.Ticker in project GDSC-SMLM by aherbert.

the class PsfDrift method computeDrift.

private void computeDrift() {
    // Create a grid of XY offset positions between 0-1 for PSF insert
    final double[] grid = new double[settings.gridSize];
    for (int i = 0; i < grid.length; i++) {
        grid[i] = (double) i / settings.gridSize;
    }
    // Configure fitting region
    final int w = 2 * settings.regionSize + 1;
    centrePixel = w / 2;
    // Check region size using the image PSF
    final double newPsfWidth = imp.getWidth() / settings.scale;
    if (Math.ceil(newPsfWidth) > w) {
        ImageJUtils.log(TITLE + ": Fitted region size (%d) is smaller than the scaled PSF (%.1f)", w, newPsfWidth);
    }
    // Create robust PSF fitting settings
    final double a = psfSettings.getPixelSize() * settings.scale;
    final double sa = PsfCalculator.squarePixelAdjustment(psfSettings.getPixelSize() * (psfSettings.getFwhm() / Gaussian2DFunction.SD_TO_FWHM_FACTOR), a);
    fitConfig.setInitialPeakStdDev(sa / a);
    fitConfig.setBackgroundFitting(settings.backgroundFitting);
    fitConfig.setNotSignalFitting(false);
    fitConfig.setComputeDeviations(false);
    fitConfig.setDisableSimpleFilter(true);
    // Create the PSF over the desired z-depth
    final int depth = (int) Math.round(settings.zDepth / psfSettings.getPixelDepth());
    int startSlice = psfSettings.getCentreImage() - depth;
    int endSlice = psfSettings.getCentreImage() + depth;
    final int nSlices = imp.getStackSize();
    startSlice = MathUtils.clip(1, nSlices, startSlice);
    endSlice = MathUtils.clip(1, nSlices, endSlice);
    final ImagePsfModel psf = createImagePsf(startSlice, endSlice, settings.scale);
    final int minz = startSlice - psfSettings.getCentreImage();
    final int maxz = endSlice - psfSettings.getCentreImage();
    final int nZ = maxz - minz + 1;
    final int gridSize2 = grid.length * grid.length;
    total = nZ * gridSize2;
    // Store all the fitting results
    final int nStartPoints = getNumberOfStartPoints();
    results = new double[total * nStartPoints][];
    // TODO - Add ability to iterate this, adjusting the current offset in the PSF
    // each iteration
    // Create a pool of workers
    final int threadCount = Prefs.getThreads();
    final Ticker ticker = ImageJUtils.createTicker(total, threadCount, "Fitting...");
    final BlockingQueue<Job> jobs = new ArrayBlockingQueue<>(threadCount * 2);
    final List<Thread> threads = new LinkedList<>();
    for (int i = 0; i < threadCount; i++) {
        final Worker worker = new Worker(jobs, psf, w, fitConfig, ticker);
        final Thread t = new Thread(worker);
        threads.add(t);
        t.start();
    }
    // Fit
    outer: for (int z = minz, i = 0; z <= maxz; z++) {
        for (int x = 0; x < grid.length; x++) {
            for (int y = 0; y < grid.length; y++, i++) {
                if (IJ.escapePressed()) {
                    break outer;
                }
                put(jobs, new Job(z, grid[x], grid[y], i));
            }
        }
    }
    // If escaped pressed then do not need to stop the workers, just return
    if (ImageJUtils.isInterrupted()) {
        ImageJUtils.finished();
        return;
    }
    // Finish all the worker threads by passing in a null job
    for (int i = 0; i < threads.size(); i++) {
        put(jobs, new Job());
    }
    // 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("Unexpected interrupt", ex);
        }
    }
    threads.clear();
    ImageJUtils.finished();
    // Plot the average and SE for the drift curve
    // Plot the recall
    final double[] zPosition = new double[nZ];
    final double[] avX = new double[nZ];
    final double[] seX = new double[nZ];
    final double[] avY = new double[nZ];
    final double[] seY = new double[nZ];
    final double[] recall = new double[nZ];
    for (int z = minz, i = 0; z <= maxz; z++, i++) {
        final Statistics statsX = new Statistics();
        final Statistics statsY = new Statistics();
        for (int s = 0; s < nStartPoints; s++) {
            int resultPosition = i * gridSize2 + s * total;
            final int endResultPosition = resultPosition + gridSize2;
            while (resultPosition < endResultPosition) {
                if (results[resultPosition] != null) {
                    statsX.add(results[resultPosition][0]);
                    statsY.add(results[resultPosition][1]);
                }
                resultPosition++;
            }
        }
        zPosition[i] = z * psfSettings.getPixelDepth();
        avX[i] = statsX.getMean();
        seX[i] = statsX.getStandardError();
        avY[i] = statsY.getMean();
        seY[i] = statsY.getStandardError();
        recall[i] = (double) statsX.getN() / (nStartPoints * gridSize2);
    }
    // Find the range from the z-centre above the recall limit
    int centre = 0;
    for (int slice = startSlice, i = 0; slice <= endSlice; slice++, i++) {
        if (slice == psfSettings.getCentreImage()) {
            centre = i;
            break;
        }
    }
    if (recall[centre] < settings.recallLimit) {
        return;
    }
    int start = centre;
    int end = centre;
    for (int i = centre; i-- > 0; ) {
        if (recall[i] < settings.recallLimit) {
            break;
        }
        start = i;
    }
    for (int i = centre; ++i < recall.length; ) {
        if (recall[i] < settings.recallLimit) {
            break;
        }
        end = i;
    }
    final int iterations = 1;
    LoessInterpolator loess = null;
    if (settings.smoothing > 0) {
        loess = new LoessInterpolator(settings.smoothing, iterations);
    }
    final double[][] smoothx = displayPlot("Drift X", "X (nm)", zPosition, avX, seX, loess, start, end);
    final double[][] smoothy = displayPlot("Drift Y", "Y (nm)", zPosition, avY, seY, loess, start, end);
    displayPlot("Recall", "Recall", zPosition, recall, null, null, start, end);
    windowOrganiser.tile();
    // Ask the user if they would like to store them in the image
    final GenericDialog gd = new GenericDialog(TITLE);
    gd.enableYesNoCancel();
    gd.hideCancelButton();
    startSlice = psfSettings.getCentreImage() - (centre - start);
    endSlice = psfSettings.getCentreImage() + (end - centre);
    ImageJUtils.addMessage(gd, "Save the drift to the PSF?\n \nSlices %d (%s nm) - %d (%s nm) above recall limit", startSlice, MathUtils.rounded(zPosition[start]), endSlice, MathUtils.rounded(zPosition[end]));
    gd.addMessage("Optionally average the end points to set drift outside the limits.\n" + "(Select zero to ignore)");
    gd.addSlider("Number_of_points", 0, 10, settings.positionsToAverage);
    gd.showDialog();
    if (gd.wasOKed()) {
        settings.positionsToAverage = Math.abs((int) gd.getNextNumber());
        final Map<Integer, Offset> oldOffset = psfSettings.getOffsetsMap();
        final boolean useOldOffset = settings.useOffset && !oldOffset.isEmpty();
        final LocalList<double[]> offset = new LocalList<>();
        final double pitch = psfSettings.getPixelSize();
        int index = 0;
        for (int i = start, slice = startSlice; i <= end; slice++, i++) {
            index = findCentre(zPosition[i], smoothx, index);
            if (index == -1) {
                ImageJUtils.log("Failed to find the offset for depth %.2f", zPosition[i]);
                continue;
            }
            // The offset should store the difference to the centre in pixels so divide by the pixel
            // pitch
            double cx = smoothx[1][index] / pitch;
            double cy = smoothy[1][index] / pitch;
            if (useOldOffset) {
                final Offset o = oldOffset.get(slice);
                if (o != null) {
                    cx += o.getCx();
                    cy += o.getCy();
                }
            }
            offset.add(new double[] { slice, cx, cy });
        }
        addMissingOffsets(startSlice, endSlice, nSlices, offset);
        final Offset.Builder offsetBuilder = Offset.newBuilder();
        final ImagePSF.Builder imagePsfBuilder = psfSettings.toBuilder();
        for (final double[] o : offset) {
            final int slice = (int) o[0];
            offsetBuilder.setCx(o[1]);
            offsetBuilder.setCy(o[2]);
            imagePsfBuilder.putOffsets(slice, offsetBuilder.build());
        }
        imagePsfBuilder.putNotes(TITLE, String.format("Solver=%s, Region=%d", PeakFit.getSolverName(fitConfig), settings.regionSize));
        imp.setProperty("Info", ImagePsfHelper.toString(imagePsfBuilder));
    }
}
Also used : LoessInterpolator(org.apache.commons.math3.analysis.interpolation.LoessInterpolator) LocalList(uk.ac.sussex.gdsc.core.utils.LocalList) ArrayBlockingQueue(java.util.concurrent.ArrayBlockingQueue) ImagePSF(uk.ac.sussex.gdsc.smlm.data.config.PSFProtos.ImagePSF) NonBlockingExtendedGenericDialog(uk.ac.sussex.gdsc.core.ij.gui.NonBlockingExtendedGenericDialog) ExtendedGenericDialog(uk.ac.sussex.gdsc.core.ij.gui.ExtendedGenericDialog) GenericDialog(ij.gui.GenericDialog) Ticker(uk.ac.sussex.gdsc.core.logging.Ticker) Statistics(uk.ac.sussex.gdsc.core.utils.Statistics) LinkedList(java.util.LinkedList) Offset(uk.ac.sussex.gdsc.smlm.data.config.PSFProtos.Offset) ConcurrentRuntimeException(org.apache.commons.lang3.concurrent.ConcurrentRuntimeException) ImagePsfModel(uk.ac.sussex.gdsc.smlm.model.ImagePsfModel)

Aggregations

Ticker (uk.ac.sussex.gdsc.core.logging.Ticker)33 LinkedList (java.util.LinkedList)16 Future (java.util.concurrent.Future)12 ImageStack (ij.ImageStack)11 ArrayBlockingQueue (java.util.concurrent.ArrayBlockingQueue)10 LocalList (uk.ac.sussex.gdsc.core.utils.LocalList)10 ConcurrentRuntimeException (org.apache.commons.lang3.concurrent.ConcurrentRuntimeException)9 ExecutorService (java.util.concurrent.ExecutorService)8 AtomicInteger (java.util.concurrent.atomic.AtomicInteger)8 ImagePlus (ij.ImagePlus)7 ArrayList (java.util.ArrayList)6 List (java.util.List)6 Statistics (uk.ac.sussex.gdsc.core.utils.Statistics)6 MemoryPeakResults (uk.ac.sussex.gdsc.smlm.results.MemoryPeakResults)6 TIntObjectHashMap (gnu.trove.map.hash.TIntObjectHashMap)5 ImageProcessor (ij.process.ImageProcessor)5 TextWindow (ij.text.TextWindow)5 Point (java.awt.Point)5 Rectangle (java.awt.Rectangle)5 UniformRandomProvider (org.apache.commons.rng.UniformRandomProvider)5