Examples with ImagePSFModel gdsc.smlm.model.ImagePSFModel used on opensource projects

Example 6 with ImagePSFModel

use of gdsc.smlm.model.ImagePSFModel 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[gridSize];
    for (int i = 0; i < grid.length; i++) grid[i] = (double) i / gridSize;
    // Configure fitting region
    final int w = 2 * regionSize + 1;
    centrePixel = w / 2;
    // Check region size using the image PSF
    double newPsfWidth = (double) imp.getWidth() / scale;
    if (Math.ceil(newPsfWidth) > w)
        Utils.log(TITLE + ": Fitted region size (%d) is smaller than the scaled PSF (%.1f)", w, newPsfWidth);
    // Create robust PSF fitting settings
    final double a = psfSettings.nmPerPixel * scale;
    final double sa = PSFCalculator.squarePixelAdjustment(psfSettings.nmPerPixel * (psfSettings.fwhm / Gaussian2DFunction.SD_TO_FWHM_FACTOR), a);
    fitConfig.setInitialPeakStdDev(sa / a);
    fitConfig.setBackgroundFitting(backgroundFitting);
    fitConfig.setNotSignalFitting(false);
    fitConfig.setComputeDeviations(false);
    fitConfig.setDisableSimpleFilter(true);
    // Create the PSF over the desired z-depth
    int depth = (int) Math.round(zDepth / psfSettings.nmPerSlice);
    int startSlice = psfSettings.zCentre - depth;
    int endSlice = psfSettings.zCentre + depth;
    int nSlices = imp.getStackSize();
    startSlice = (startSlice < 1) ? 1 : (startSlice > nSlices) ? nSlices : startSlice;
    endSlice = (endSlice < 1) ? 1 : (endSlice > nSlices) ? nSlices : endSlice;
    ImagePSFModel psf = createImagePSF(startSlice, endSlice);
    int minz = startSlice - psfSettings.zCentre;
    int maxz = endSlice - psfSettings.zCentre;
    final int nZ = maxz - minz + 1;
    final int gridSize2 = grid.length * grid.length;
    total = nZ * gridSize2;
    // Store all the fitting results
    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
    int nThreads = Prefs.getThreads();
    BlockingQueue<Job> jobs = new ArrayBlockingQueue<Job>(nThreads * 2);
    List<Worker> workers = new LinkedList<Worker>();
    List<Thread> threads = new LinkedList<Thread>();
    for (int i = 0; i < nThreads; i++) {
        Worker worker = new Worker(jobs, psf, w, fitConfig);
        Thread t = new Thread(worker);
        workers.add(worker);
        threads.add(t);
        t.start();
    }
    // Fit 
    Utils.showStatus("Fitting ...");
    final int step = Utils.getProgressInterval(total);
    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 (i % step == 0) {
                IJ.showProgress(i, total);
            }
        }
    }
    // If escaped pressed then do not need to stop the workers, just return
    if (Utils.isInterrupted()) {
        IJ.showProgress(1);
        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 (InterruptedException e) {
            e.printStackTrace();
        }
    }
    threads.clear();
    IJ.showProgress(1);
    IJ.showStatus("");
    // Plot the average and SE for the drift curve
    // Plot the recall
    double[] zPosition = new double[nZ];
    double[] avX = new double[nZ];
    double[] seX = new double[nZ];
    double[] avY = new double[nZ];
    double[] seY = new double[nZ];
    double[] recall = new double[nZ];
    for (int z = minz, i = 0; z <= maxz; z++, i++) {
        Statistics statsX = new Statistics();
        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.nmPerSlice;
        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.zCentre) {
            centre = i;
            break;
        }
    }
    if (recall[centre] < recallLimit)
        return;
    int start = centre, end = centre;
    for (int i = centre; i-- > 0; ) {
        if (recall[i] < recallLimit)
            break;
        start = i;
    }
    for (int i = centre; ++i < recall.length; ) {
        if (recall[i] < recallLimit)
            break;
        end = i;
    }
    int iterations = 1;
    LoessInterpolator loess = null;
    if (smoothing > 0)
        loess = new LoessInterpolator(smoothing, iterations);
    double[][] smoothx = displayPlot("Drift X", "X (nm)", zPosition, avX, seX, loess, start, end);
    double[][] smoothy = displayPlot("Drift Y", "Y (nm)", zPosition, avY, seY, loess, start, end);
    displayPlot("Recall", "Recall", zPosition, recall, null, null, start, end);
    WindowOrganiser wo = new WindowOrganiser();
    wo.tileWindows(idList);
    // Ask the user if they would like to store them in the image
    GenericDialog gd = new GenericDialog(TITLE);
    gd.enableYesNoCancel();
    gd.hideCancelButton();
    startSlice = psfSettings.zCentre - (centre - start);
    endSlice = psfSettings.zCentre + (end - centre);
    gd.addMessage(String.format("Save the drift to the PSF?\n \nSlices %d (%s nm) - %d (%s nm) above recall limit", startSlice, Utils.rounded(zPosition[start]), endSlice, Utils.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, positionsToAverage);
    gd.showDialog();
    if (gd.wasOKed()) {
        positionsToAverage = Math.abs((int) gd.getNextNumber());
        ArrayList<PSFOffset> offset = new ArrayList<PSFOffset>();
        final double pitch = psfSettings.nmPerPixel;
        int j = 0, jj = 0;
        for (int i = start, slice = startSlice; i <= end; slice++, i++) {
            j = findCentre(zPosition[i], smoothx, j);
            if (j == -1) {
                Utils.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][j] / pitch;
            double cy = smoothy[1][j] / pitch;
            jj = findOffset(slice, jj);
            if (jj != -1) {
                cx += psfSettings.offset[jj].cx;
                cy += psfSettings.offset[jj].cy;
            }
            offset.add(new PSFOffset(slice, cx, cy));
        }
        addMissingOffsets(startSlice, endSlice, nSlices, offset);
        psfSettings.offset = offset.toArray(new PSFOffset[offset.size()]);
        psfSettings.addNote(TITLE, String.format("Solver=%s, Region=%d", PeakFit.getSolverName(fitConfig), regionSize));
        imp.setProperty("Info", XmlUtils.toXML(psfSettings));
    }
}