Examples with PSFSettings gdsc.smlm.ij.settings.PSFSettings used on opensource projects

Example 1 with PSFSettings

use of gdsc.smlm.ij.settings.PSFSettings in project GDSC-SMLM by aherbert.

the class CreateData method getImageHWHM.

/**
	 * Get the PSF half-width at half-maxima from the Image PSF
	 * 
	 * @return
	 */
private double getImageHWHM() {
    ImagePlus imp = WindowManager.getImage(settings.psfImageName);
    if (imp == null) {
        IJ.error(TITLE, "Unable to create the PSF model from image: " + settings.psfImageName);
        return -1;
    }
    Object o = XmlUtils.fromXML(imp.getProperty("Info").toString());
    if (!(o != null && o instanceof PSFSettings)) {
        IJ.error(TITLE, "Unknown PSF settings for image: " + imp.getTitle());
        return -1;
    }
    PSFSettings psfSettings = (PSFSettings) o;
    if (psfSettings.fwhm <= 0) {
        IJ.error(TITLE, "Unknown PSF FWHM setting for image: " + imp.getTitle());
        return -1;
    }
    if (psfSettings.nmPerPixel <= 0) {
        IJ.error(TITLE, "Unknown PSF nm/pixel setting for image: " + imp.getTitle());
        return -1;
    }
    // output image
    return 0.5 * psfSettings.fwhm * psfSettings.nmPerPixel / settings.pixelPitch;
}

Example 2 with PSFSettings

use of gdsc.smlm.ij.settings.PSFSettings in project GDSC-SMLM by aherbert.

the class PSFCreator method run.

/*
	 * (non-Javadoc)
	 * 
	 * @see ij.plugin.filter.PlugInFilter#run(ij.process.ImageProcessor)
	 */
public void run(ImageProcessor ip) {
    loadConfiguration();
    BasePoint[] spots = getSpots();
    if (spots.length == 0) {
        IJ.error(TITLE, "No spots without neighbours within " + (boxRadius * 2) + "px");
        return;
    }
    ImageStack stack = getImageStack();
    final int width = imp.getWidth();
    final int height = imp.getHeight();
    final int currentSlice = imp.getSlice();
    // Adjust settings for a single maxima
    config.setIncludeNeighbours(false);
    fitConfig.setDuplicateDistance(0);
    ArrayList<double[]> centres = new ArrayList<double[]>(spots.length);
    int iterations = 1;
    LoessInterpolator loess = new LoessInterpolator(smoothing, iterations);
    // TODO - The fitting routine may not produce many points. In this instance the LOESS interpolator
    // fails to smooth the data very well. A higher bandwidth helps this but perhaps 
    // try a different smoothing method.
    // For each spot
    Utils.log(TITLE + ": " + imp.getTitle());
    Utils.log("Finding spot locations...");
    Utils.log("  %d spot%s without neighbours within %dpx", spots.length, ((spots.length == 1) ? "" : "s"), (boxRadius * 2));
    StoredDataStatistics averageSd = new StoredDataStatistics();
    StoredDataStatistics averageA = new StoredDataStatistics();
    Statistics averageRange = new Statistics();
    MemoryPeakResults allResults = new MemoryPeakResults();
    allResults.setName(TITLE);
    allResults.setBounds(new Rectangle(0, 0, width, height));
    MemoryPeakResults.addResults(allResults);
    for (int n = 1; n <= spots.length; n++) {
        BasePoint spot = spots[n - 1];
        final int x = (int) spot.getX();
        final int y = (int) spot.getY();
        MemoryPeakResults results = fitSpot(stack, width, height, x, y);
        allResults.addAllf(results.getResults());
        if (results.size() < 5) {
            Utils.log("  Spot %d: Not enough fit results %d", n, results.size());
            continue;
        }
        // Get the results for the spot centre and width
        double[] z = new double[results.size()];
        double[] xCoord = new double[z.length];
        double[] yCoord = new double[z.length];
        double[] sd = new double[z.length];
        double[] a = new double[z.length];
        int i = 0;
        for (PeakResult peak : results.getResults()) {
            z[i] = peak.getFrame();
            xCoord[i] = peak.getXPosition() - x;
            yCoord[i] = peak.getYPosition() - y;
            sd[i] = FastMath.max(peak.getXSD(), peak.getYSD());
            a[i] = peak.getAmplitude();
            i++;
        }
        // Smooth the amplitude plot
        double[] smoothA = loess.smooth(z, a);
        // Find the maximum amplitude
        int maximumIndex = findMaximumIndex(smoothA);
        // Find the range at a fraction of the max. This is smoothed to find the X/Y centre
        int start = 0, stop = smoothA.length - 1;
        double limit = smoothA[maximumIndex] * amplitudeFraction;
        for (int j = 0; j < smoothA.length; j++) {
            if (smoothA[j] > limit) {
                start = j;
                break;
            }
        }
        for (int j = smoothA.length; j-- > 0; ) {
            if (smoothA[j] > limit) {
                stop = j;
                break;
            }
        }
        averageRange.add(stop - start + 1);
        // Extract xy centre coords and smooth
        double[] smoothX = new double[stop - start + 1];
        double[] smoothY = new double[smoothX.length];
        double[] smoothSd = new double[smoothX.length];
        double[] newZ = new double[smoothX.length];
        for (int j = start, k = 0; j <= stop; j++, k++) {
            smoothX[k] = xCoord[j];
            smoothY[k] = yCoord[j];
            smoothSd[k] = sd[j];
            newZ[k] = z[j];
        }
        smoothX = loess.smooth(newZ, smoothX);
        smoothY = loess.smooth(newZ, smoothY);
        smoothSd = loess.smooth(newZ, smoothSd);
        // Since the amplitude is not very consistent move from this peak to the 
        // lowest width which is the in-focus spot.
        maximumIndex = findMinimumIndex(smoothSd, maximumIndex - start);
        // Find the centre at the amplitude peak
        double cx = smoothX[maximumIndex] + x;
        double cy = smoothY[maximumIndex] + y;
        int cz = (int) newZ[maximumIndex];
        double csd = smoothSd[maximumIndex];
        double ca = smoothA[maximumIndex + start];
        // The average should weight the SD using the signal for each spot
        averageSd.add(smoothSd[maximumIndex]);
        averageA.add(ca);
        if (ignoreSpot(n, z, a, smoothA, xCoord, yCoord, sd, newZ, smoothX, smoothY, smoothSd, cx, cy, cz, csd)) {
            Utils.log("  Spot %d was ignored", n);
            continue;
        }
        // Store result - it may have been moved interactively
        maximumIndex += this.slice - cz;
        cz = (int) newZ[maximumIndex];
        csd = smoothSd[maximumIndex];
        ca = smoothA[maximumIndex + start];
        Utils.log("  Spot %d => x=%.2f, y=%.2f, z=%d, sd=%.2f, A=%.2f\n", n, cx, cy, cz, csd, ca);
        centres.add(new double[] { cx, cy, cz, csd, n });
    }
    if (interactiveMode) {
        imp.setSlice(currentSlice);
        imp.setOverlay(null);
        // Hide the amplitude and spot plots
        Utils.hide(TITLE_AMPLITUDE);
        Utils.hide(TITLE_PSF_PARAMETERS);
    }
    if (centres.isEmpty()) {
        String msg = "No suitable spots could be identified centres";
        Utils.log(msg);
        IJ.error(TITLE, msg);
        return;
    }
    // Find the limits of the z-centre
    int minz = (int) centres.get(0)[2];
    int maxz = minz;
    for (double[] centre : centres) {
        if (minz > centre[2])
            minz = (int) centre[2];
        else if (maxz < centre[2])
            maxz = (int) centre[2];
    }
    IJ.showStatus("Creating PSF image");
    // Create a stack that can hold all the data.
    ImageStack psf = createStack(stack, minz, maxz, magnification);
    // For each spot
    Statistics stats = new Statistics();
    boolean ok = true;
    for (int i = 0; ok && i < centres.size(); i++) {
        double progress = (double) i / centres.size();
        final double increment = 1.0 / (stack.getSize() * centres.size());
        IJ.showProgress(progress);
        double[] centre = centres.get(i);
        // Extract the spot
        float[][] spot = new float[stack.getSize()][];
        Rectangle regionBounds = null;
        for (int slice = 1; slice <= stack.getSize(); slice++) {
            ImageExtractor ie = new ImageExtractor((float[]) stack.getPixels(slice), width, height);
            if (regionBounds == null)
                regionBounds = ie.getBoxRegionBounds((int) centre[0], (int) centre[1], boxRadius);
            spot[slice - 1] = ie.crop(regionBounds);
        }
        int n = (int) centre[4];
        final float b = getBackground(n, spot);
        if (!subtractBackgroundAndWindow(spot, b, regionBounds.width, regionBounds.height, centre, loess)) {
            Utils.log("  Spot %d was ignored", n);
            continue;
        }
        stats.add(b);
        // Adjust the centre using the crop
        centre[0] -= regionBounds.x;
        centre[1] -= regionBounds.y;
        // This takes a long time so this should track progress
        ok = addToPSF(maxz, magnification, psf, centre, spot, regionBounds, progress, increment, centreEachSlice);
    }
    if (interactiveMode) {
        Utils.hide(TITLE_INTENSITY);
    }
    IJ.showProgress(1);
    if (threadPool != null) {
        threadPool.shutdownNow();
        threadPool = null;
    }
    if (!ok || stats.getN() == 0)
        return;
    final double avSd = getAverage(averageSd, averageA, 2);
    Utils.log("  Average background = %.2f, Av. SD = %s px", stats.getMean(), Utils.rounded(avSd, 4));
    normalise(psf, maxz, avSd * magnification, false);
    IJ.showProgress(1);
    psfImp = Utils.display("PSF", psf);
    psfImp.setSlice(maxz);
    psfImp.resetDisplayRange();
    psfImp.updateAndDraw();
    double[][] fitCom = new double[2][psf.getSize()];
    Arrays.fill(fitCom[0], Double.NaN);
    Arrays.fill(fitCom[1], Double.NaN);
    double fittedSd = fitPSF(psf, loess, maxz, averageRange.getMean(), fitCom);
    // Compute the drift in the PSF:
    // - Use fitted centre if available; otherwise find CoM for each frame
    // - express relative to the average centre
    double[][] com = calculateCentreOfMass(psf, fitCom, nmPerPixel / magnification);
    double[] slice = Utils.newArray(psf.getSize(), 1, 1.0);
    String title = TITLE + " CoM Drift";
    Plot2 plot = new Plot2(title, "Slice", "Drift (nm)");
    plot.addLabel(0, 0, "Red = X; Blue = Y");
    //double[] limitsX = Maths.limits(com[0]);
    //double[] limitsY = Maths.limits(com[1]);
    double[] limitsX = getLimits(com[0]);
    double[] limitsY = getLimits(com[1]);
    plot.setLimits(1, psf.getSize(), Math.min(limitsX[0], limitsY[0]), Math.max(limitsX[1], limitsY[1]));
    plot.setColor(Color.red);
    plot.addPoints(slice, com[0], Plot.DOT);
    plot.addPoints(slice, loess.smooth(slice, com[0]), Plot.LINE);
    plot.setColor(Color.blue);
    plot.addPoints(slice, com[1], Plot.DOT);
    plot.addPoints(slice, loess.smooth(slice, com[1]), Plot.LINE);
    Utils.display(title, plot);
    // TODO - Redraw the PSF with drift correction applied. 
    // This means that the final image should have no drift.
    // This is relevant when combining PSF images. It doesn't matter too much for simulations 
    // unless the drift is large.
    // Add Image properties containing the PSF details
    final double fwhm = getFWHM(psf, maxz);
    psfImp.setProperty("Info", XmlUtils.toXML(new PSFSettings(maxz, nmPerPixel / magnification, nmPerSlice, stats.getN(), fwhm, createNote())));
    Utils.log("%s : z-centre = %d, nm/Pixel = %s, nm/Slice = %s, %d images, PSF SD = %s nm, FWHM = %s px\n", psfImp.getTitle(), maxz, Utils.rounded(nmPerPixel / magnification, 3), Utils.rounded(nmPerSlice, 3), stats.getN(), Utils.rounded(fittedSd * nmPerPixel, 4), Utils.rounded(fwhm));
    createInteractivePlots(psf, maxz, nmPerPixel / magnification, fittedSd * nmPerPixel);
    IJ.showStatus("");
}

Example 3 with PSFSettings

use of gdsc.smlm.ij.settings.PSFSettings in project GDSC-SMLM by aherbert.

the class PSFDrift method createImageList.

public static List<String> createImageList() {
    List<String> titles = new LinkedList<String>();
    int[] ids = WindowManager.getIDList();
    if (ids != null) {
        for (int id : ids) {
            ImagePlus imp = WindowManager.getImage(id);
            if (imp != null) {
                // Image must be greyscale
                if (imp.getType() == ImagePlus.GRAY8 || imp.getType() == ImagePlus.GRAY16 || imp.getType() == ImagePlus.GRAY32) {
                    // Image must be square and a stack of a single channel
                    if (imp.getWidth() == imp.getHeight() && imp.getNChannels() == 1) {
                        // Check if these are PSF images created by the SMLM plugins
                        PSFSettings psfSettings = getPSFSettings(imp);
                        if (psfSettings != null) {
                            if (psfSettings.zCentre <= 0) {
                                Utils.log(TITLE + ": Unknown PSF z-centre setting for image: " + imp.getTitle());
                                continue;
                            }
                            if (psfSettings.nmPerPixel <= 0) {
                                Utils.log(TITLE + ": Unknown PSF nm/pixel setting for image: " + imp.getTitle());
                                continue;
                            }
                            if (psfSettings.nmPerSlice <= 0) {
                                Utils.log(TITLE + ": Unknown PSF nm/slice setting for image: " + imp.getTitle());
                                continue;
                            }
                            if (psfSettings.fwhm <= 0) {
                                Utils.log(TITLE + ": Unknown PSF FWHM setting for image: " + imp.getTitle());
                                continue;
                            }
                            titles.add(imp.getTitle());
                        }
                    }
                }
            }
        }
    }
    return titles;
}

Example 4 with PSFSettings

use of gdsc.smlm.ij.settings.PSFSettings in project GDSC-SMLM by aherbert.

the class PSFCombiner method combineImages.

private void combineImages() {
    double nmPerPixel = getNmPerPixel();
    if (nmPerPixel <= 0)
        return;
    double nmPerSlice = getNmPerSlice();
    if (nmPerPixel <= 0)
        return;
    // Find the lowest start point
    int min = 0;
    for (PSF psf : input) {
        if (min > psf.start)
            min = psf.start;
    }
    // Shift all stacks and find the dimensions
    final int shift = -min;
    int max = 0, size = 0;
    int totalImages = 0;
    for (PSF psf : input) {
        psf.start += shift;
        totalImages += psf.psfSettings.nImages;
        if (max < psf.getEnd())
            max = psf.getEnd();
        if (size < psf.getSize())
            size = psf.getSize();
    }
    // Create a stack to hold all the images
    ImageStack stack = new ImageStack(size, size, max);
    for (int n = 1; n <= max; n++) stack.setPixels(new float[size * size], n);
    // Insert all the PSFs
    IJ.showStatus("Creating combined image ...");
    int imageNo = 0;
    double fraction = 1.0 / input.size();
    for (PSF psf : input) {
        double progress = imageNo * fraction;
        ImageStack psfStack = psf.psfStack;
        final int offsetXY = (psf.getSize() - size) / 2;
        final int offsetZ = psf.start;
        final int w = psf.getSize();
        final double weight = (1.0 * psf.psfSettings.nImages) / totalImages;
        final double increment = fraction / psfStack.getSize();
        for (int n = 1; n <= psfStack.getSize(); n++) {
            IJ.showProgress(progress += increment);
            // Get the data and adjust using the weight
            float[] psfData = ImageConverter.getData(psfStack.getProcessor(n));
            for (int i = 0; i < psfData.length; i++) psfData[i] *= weight;
            // Insert into the combined PSF
            ImageProcessor ip = stack.getProcessor(n + offsetZ);
            ip.copyBits(new FloatProcessor(w, w, psfData, null), offsetXY, offsetXY, Blitter.ADD);
        }
        imageNo++;
    }
    // IJ.showStatus("Normalising ...");
    // PSFCreator.normalise(stack, 1 + shift);
    IJ.showProgress(1);
    IJ.showStatus("");
    ImagePlus imp = Utils.display("Combined PSF", stack);
    imp.setSlice(1 + shift);
    imp.resetDisplayRange();
    imp.updateAndDraw();
    final double fwhm = getFWHM();
    imp.setProperty("Info", XmlUtils.toXML(new PSFSettings(imp.getSlice(), nmPerPixel, nmPerSlice, totalImages, fwhm)));
    Utils.log("%s : z-centre = %d, nm/Pixel = %s, nm/Slice = %s, %d images, FWHM = %s\n", imp.getTitle(), imp.getSlice(), Utils.rounded(nmPerPixel), Utils.rounded(nmPerSlice), totalImages, Utils.rounded(fwhm));
}

Example 5 with PSFSettings

use of gdsc.smlm.ij.settings.PSFSettings in project GDSC-SMLM by aherbert.

the class CreateData method createImagePSF.

/**
	 * Create a PSF model from the image that contains all the z-slices needed to draw the given localisations
	 * 
	 * @param localisationSets
	 * @return
	 */
private ImagePSFModel createImagePSF(List<LocalisationModelSet> localisationSets) {
    ImagePlus imp = WindowManager.getImage(settings.psfImageName);
    if (imp == null) {
        IJ.error(TITLE, "Unable to create the PSF model from image: " + settings.psfImageName);
        return null;
    }
    try {
        Object o = XmlUtils.fromXML(imp.getProperty("Info").toString());
        if (!(o != null && o instanceof PSFSettings))
            throw new RuntimeException("Unknown PSF settings for image: " + imp.getTitle());
        PSFSettings psfSettings = (PSFSettings) o;
        // Check all the settings have values
        if (psfSettings.nmPerPixel <= 0)
            throw new RuntimeException("Missing nmPerPixel calibration settings for image: " + imp.getTitle());
        if (psfSettings.nmPerSlice <= 0)
            throw new RuntimeException("Missing nmPerSlice calibration settings for image: " + imp.getTitle());
        if (psfSettings.zCentre <= 0)
            throw new RuntimeException("Missing zCentre calibration settings for image: " + imp.getTitle());
        if (psfSettings.fwhm <= 0)
            throw new RuntimeException("Missing FWHM calibration settings for image: " + imp.getTitle());
        // To save memory construct the Image PSF using only the slices that are within 
        // the depth of field of the simulation
        double minZ = Double.POSITIVE_INFINITY, maxZ = Double.NEGATIVE_INFINITY;
        for (LocalisationModelSet l : localisationSets) {
            for (LocalisationModel m : l.getLocalisations()) {
                final double z = m.getZ();
                if (minZ > z)
                    minZ = z;
                if (maxZ < z)
                    maxZ = z;
            }
        }
        int nSlices = imp.getStackSize();
        // z-centre should be an index and not the ImageJ slice number so subtract 1
        int zCentre = psfSettings.zCentre - 1;
        // Calculate the start/end slices to cover the depth of field
        // This logic must match the ImagePSFModel.
        final double unitsPerSlice = psfSettings.nmPerSlice / settings.pixelPitch;
        // We assume the PSF was imaged axially with increasing z-stage position (moving the stage 
        // closer to the objective). Thus higher z-coordinate are for higher slice numbers.
        int lower = (int) Math.round(minZ / unitsPerSlice) + zCentre;
        int upper = (int) Math.round(maxZ / unitsPerSlice) + zCentre;
        upper = (upper < 0) ? 0 : (upper >= nSlices) ? nSlices - 1 : upper;
        lower = (lower < 0) ? 0 : (lower >= nSlices) ? nSlices - 1 : lower;
        // Image PSF requires the z-centre for normalisation
        if (!(lower <= zCentre && upper >= zCentre)) {
            // Ensure we include the zCentre
            lower = Math.min(lower, zCentre);
            upper = Math.max(upper, zCentre);
        }
        final double noiseFraction = 1e-3;
        final ImagePSFModel model = new ImagePSFModel(extractImageStack(imp, lower, upper), zCentre - lower, psfSettings.nmPerPixel / settings.pixelPitch, unitsPerSlice, psfSettings.fwhm, noiseFraction);
        // Add the calibrated centres
        if (psfSettings.offset != null) {
            final int sliceOffset = lower + 1;
            for (PSFOffset offset : psfSettings.offset) {
                model.setRelativeCentre(offset.slice - sliceOffset, offset.cx, offset.cy);
            }
        }
        // Initialise the HWHM table so that it can be cloned
        model.initialiseHWHM();
        return model;
    } catch (Exception e) {
        IJ.error(TITLE, "Unable to create the image PSF model:\n" + e.getMessage());
        return null;
    }
}