Examples with PeakResult gdsc.smlm.results.PeakResult used on opensource projects

Example 1 with PeakResult

use of gdsc.smlm.results.PeakResult 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 2 with PeakResult

use of gdsc.smlm.results.PeakResult in project GDSC-SMLM by aherbert.

the class ResultsImageSampler method getSample.

/**
	 * Gets the sample image. The image is a stack of the samples with an overlay of the localisation positions. The
	 * info property is set with details of the localisations and the image is calibrated.
	 *
	 * @param nNo
	 *            the number of samples with no localisations
	 * @param nLow
	 *            the number of samples with low localisations
	 * @param nHigh
	 *            the number of samples with high localisations
	 * @return the sample image (could be null if no samples were made)
	 */
public ImagePlus getSample(int nNo, int nLow, int nHigh) {
    ImageStack out = new ImageStack(size, size);
    if (!isValid())
        return null;
    list.clearf();
    // empty
    for (int i : Random.sample(nNo, no.length, r)) list.add(ResultsSample.createEmpty(no[i]));
    // low
    for (int i : Random.sample(nLow, lower, r)) list.add(data[i]);
    // high
    for (int i : Random.sample(nHigh, upper, r)) list.add(data[i + lower]);
    if (list.isEmpty())
        return null;
    double nmPerPixel = 1;
    if (results.getCalibration() != null) {
        Calibration calibration = results.getCalibration();
        if (calibration.hasNmPerPixel()) {
            nmPerPixel = calibration.getNmPerPixel();
        }
    }
    // Sort descending by number in the frame
    ResultsSample[] sample = list.toArray(new ResultsSample[list.size()]);
    Arrays.sort(sample, rcc);
    int[] xyz = new int[3];
    Rectangle stackBounds = new Rectangle(stack.getWidth(), stack.getHeight());
    Overlay overlay = new Overlay();
    float[] ox = new float[10], oy = new float[10];
    StringBuilder sb = new StringBuilder();
    if (nmPerPixel == 1)
        sb.append("Sample X Y Z Signal\n");
    else
        sb.append("Sample X(nm) Y(nm) Z(nm) Signal\n");
    for (ResultsSample s : sample) {
        getXYZ(s.index, xyz);
        // Construct the region to extract
        Rectangle target = new Rectangle(xyz[0], xyz[1], size, size);
        target = target.intersection(stackBounds);
        if (target.width == 0 || target.height == 0)
            continue;
        // Extract the frame
        int slice = xyz[2];
        ImageProcessor ip = stack.getProcessor(slice);
        // Cut out the desired pixels (some may be blank if the block overruns the source image)
        ImageProcessor ip2 = ip.createProcessor(size, size);
        for (int y = 0; y < target.height; y++) for (int x = 0, i = y * size, index = (y + target.y) * ip.getWidth() + target.x; x < target.width; x++, i++, index++) {
            ip2.setf(i, ip.getf(index));
        }
        int size = s.size();
        if (size > 0) {
            int position = out.getSize() + 1;
            // Create an ROI with the localisations
            for (int i = 0; i < size; i++) {
                PeakResult p = s.list.get(i);
                ox[i] = p.getXPosition() - xyz[0];
                oy[i] = p.getYPosition() - xyz[1];
                sb.append(position).append(' ');
                sb.append(Utils.rounded(ox[i] * nmPerPixel)).append(' ');
                sb.append(Utils.rounded(oy[i] * nmPerPixel)).append(' ');
                // Z can be stored in the error field
                sb.append(Utils.rounded(p.error * nmPerPixel)).append(' ');
                sb.append(Utils.rounded(p.getSignal())).append('\n');
            }
            PointRoi roi = new PointRoi(ox, oy, size);
            roi.setPosition(position);
            overlay.add(roi);
        }
        out.addSlice(String.format("Frame=%d @ %d,%d px (n=%d)", slice, xyz[0], xyz[1], size), ip2.getPixels());
    }
    if (out.getSize() == 0)
        return null;
    ImagePlus imp = new ImagePlus("Sample", out);
    imp.setOverlay(overlay);
    // Note: Only the info property can be saved to a TIFF file
    imp.setProperty("Info", sb.toString());
    if (nmPerPixel != 1) {
        ij.measure.Calibration cal = new ij.measure.Calibration();
        cal.setUnit("nm");
        cal.pixelHeight = cal.pixelWidth = nmPerPixel;
        imp.setCalibration(cal);
    }
    return imp;
}

Example 3 with PeakResult

use of gdsc.smlm.results.PeakResult in project GDSC-SMLM by aherbert.

the class PSFImagePeakResults method addAll.

/*
	 * (non-Javadoc)
	 * 
	 * @see gdsc.utils.fitting.results.PeakResults#addAll(java.util.Collection)
	 */
public void addAll(Collection<PeakResult> results) {
    if (!imageActive)
        return;
    // TODO - Make this more efficient. It could use worker threads to increase speed.
    int i = 0;
    for (PeakResult result : results) {
        addPeak(result.getFrame(), result.origX, result.origY, result.origValue, result.error, result.noise, result.params, result.paramsStdDev);
        if (++i % 64 == 0) {
            updateImage();
            if (!imageActive)
                return;
        }
    }
    updateImage();
}

Example 4 with PeakResult

use of gdsc.smlm.results.PeakResult in project GDSC-SMLM by aherbert.

the class FIRE method showFrcTimeEvolution.

private void showFrcTimeEvolution(String name, double fireNumber, ThresholdMethod thresholdMethod, double fourierImageScale, int imageSize) {
    IJ.showStatus("Calculating FRC time evolution curve...");
    List<PeakResult> list = results.getResults();
    int nSteps = 10;
    int maxT = list.get(list.size() - 1).getFrame();
    if (maxT == 0)
        maxT = list.size();
    int step = maxT / nSteps;
    TDoubleArrayList x = new TDoubleArrayList();
    TDoubleArrayList y = new TDoubleArrayList();
    double yMin = fireNumber;
    double yMax = fireNumber;
    MemoryPeakResults newResults = new MemoryPeakResults();
    newResults.copySettings(results);
    int i = 0;
    for (int t = step; t <= maxT - step; t += step) {
        while (i < list.size()) {
            if (list.get(i).getFrame() <= t) {
                newResults.add(list.get(i));
                i++;
            } else
                break;
        }
        x.add((double) t);
        FIRE f = this.copy();
        FireResult result = f.calculateFireNumber(fourierMethod, samplingMethod, thresholdMethod, fourierImageScale, imageSize);
        double fire = (result == null) ? 0 : result.fireNumber;
        y.add(fire);
        yMin = FastMath.min(yMin, fire);
        yMax = FastMath.max(yMax, fire);
    }
    // Add the final fire number
    x.add((double) maxT);
    y.add(fireNumber);
    double[] xValues = x.toArray();
    double[] yValues = y.toArray();
    String units = "px";
    if (results.getCalibration() != null) {
        nmPerPixel = results.getNmPerPixel();
        units = "nm";
    }
    String title = name + " FRC Time Evolution";
    Plot2 plot = new Plot2(title, "Frames", "Resolution (" + units + ")", (float[]) null, (float[]) null);
    double range = Math.max(1, yMax - yMin) * 0.05;
    plot.setLimits(xValues[0], xValues[xValues.length - 1], yMin - range, yMax + range);
    plot.setColor(Color.red);
    plot.addPoints(xValues, yValues, Plot.CONNECTED_CIRCLES);
    Utils.display(title, plot);
}

Example 5 with PeakResult

use of gdsc.smlm.results.PeakResult in project GDSC-SMLM by aherbert.

the class FIRE method createImages.

/**
	 * Creates the images to use for the FIRE calculation. This must be called after
	 * {@link #initialise(MemoryPeakResults, MemoryPeakResults)}.
	 *
	 * @param fourierImageScale
	 *            the fourier image scale (set to zero to auto compute)
	 * @param imageSize
	 *            the image size
	 * @param useSignal
	 *            Use the localisation signal to weight the intensity. The default uses a value of 1 per localisation.
	 * @return the fire images
	 */
public FireImages createImages(double fourierImageScale, int imageSize, boolean useSignal) {
    if (results == null)
        return null;
    final SignalProvider signalProvider = (useSignal && (results.getHead().getSignal() > 0)) ? new PeakSignalProvider() : new FixedSignalProvider();
    // Draw images using the existing IJ routines.
    Rectangle bounds = new Rectangle(0, 0, (int) Math.ceil(dataBounds.getWidth()), (int) Math.ceil(dataBounds.getHeight()));
    boolean weighted = true;
    boolean equalised = false;
    double imageScale;
    if (fourierImageScale <= 0) {
        double size = FastMath.max(bounds.width, bounds.height);
        if (size <= 0)
            size = 1;
        imageScale = imageSize / size;
    } else
        imageScale = fourierImageScale;
    IJImagePeakResults image1 = ImagePeakResultsFactory.createPeakResultsImage(ResultsImage.NONE, weighted, equalised, "IP1", bounds, 1, 1, imageScale, 0, ResultsMode.ADD);
    image1.setDisplayImage(false);
    image1.begin();
    IJImagePeakResults image2 = ImagePeakResultsFactory.createPeakResultsImage(ResultsImage.NONE, weighted, equalised, "IP2", bounds, 1, 1, imageScale, 0, ResultsMode.ADD);
    image2.setDisplayImage(false);
    image2.begin();
    float minx = (float) dataBounds.getX();
    float miny = (float) dataBounds.getY();
    if (this.results2 != null) {
        // Two image comparison
        for (PeakResult p : results) {
            float x = p.getXPosition() - minx;
            float y = p.getYPosition() - miny;
            image1.add(x, y, signalProvider.getSignal(p));
        }
        for (PeakResult p : results2) {
            float x = p.getXPosition() - minx;
            float y = p.getYPosition() - miny;
            image2.add(x, y, signalProvider.getSignal(p));
        }
    } else {
        // Block sampling.
        // Ensure we have at least 2 even sized blocks.
        int blockSize = Math.min(results.size() / 2, Math.max(1, FIRE.blockSize));
        int nBlocks = (int) Math.ceil((double) results.size() / blockSize);
        while (nBlocks <= 1 && blockSize > 1) {
            blockSize /= 2;
            nBlocks = (int) Math.ceil((double) results.size() / blockSize);
        }
        if (nBlocks <= 1)
            // This should not happen since the results should contain at least 2 localisations
            return null;
        if (blockSize != FIRE.blockSize)
            IJ.log(TITLE + " Warning: Changed block size to " + blockSize);
        int i = 0;
        int block = 0;
        PeakResult[][] blocks = new PeakResult[nBlocks][blockSize];
        for (PeakResult p : results) {
            if (i == blockSize) {
                block++;
                i = 0;
            }
            blocks[block][i++] = p;
        }
        // Truncate last block
        blocks[block] = Arrays.copyOf(blocks[block], i);
        final int[] indices = Utils.newArray(nBlocks, 0, 1);
        if (randomSplit)
            MathArrays.shuffle(indices);
        for (int index : indices) {
            // Split alternating so just rotate
            IJImagePeakResults image = image1;
            image1 = image2;
            image2 = image;
            for (PeakResult p : blocks[index]) {
                float x = p.getXPosition() - minx;
                float y = p.getYPosition() - miny;
                image.add(x, y, signalProvider.getSignal(p));
            }
        }
    }
    image1.end();
    ImageProcessor ip1 = image1.getImagePlus().getProcessor();
    image2.end();
    ImageProcessor ip2 = image2.getImagePlus().getProcessor();
    if (maxPerBin > 0 && signalProvider instanceof FixedSignalProvider) {
        // We can eliminate over-sampled pixels
        for (int i = ip1.getPixelCount(); i-- > 0; ) {
            if (ip1.getf(i) > maxPerBin)
                ip1.setf(i, maxPerBin);
            if (ip2.getf(i) > maxPerBin)
                ip2.setf(i, maxPerBin);
        }
    }
    return new FireImages(ip1, ip2, nmPerPixel / imageScale);
}