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

Example 81 with PeakResult

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

the class DiffusionRateTest method aggregateIntoFrames.

private void aggregateIntoFrames(ArrayList<Point> points, boolean addError, double precisionInPixels, RandomGenerator[] random) {
    if (myAggregateSteps < 1)
        return;
    MemoryPeakResults results = new MemoryPeakResults(points.size() / myAggregateSteps);
    Calibration cal = new Calibration(settings.pixelPitch, 1, myAggregateSteps * 1000.0 / settings.stepsPerSecond);
    results.setCalibration(cal);
    results.setName(TITLE + " Aggregated");
    MemoryPeakResults.addResults(results);
    lastSimulatedDataset[1] = results.getName();
    int id = 0;
    int peak = 1;
    int n = 0;
    double cx = 0, cy = 0;
    // Get the mean square distance
    double sum = 0;
    int count = 0;
    PeakResult last = null;
    for (Point result : points) {
        final boolean newId = result.id != id;
        if (n >= myAggregateSteps || newId) {
            if (n != 0) {
                final float[] params = new float[7];
                double[] xyz = new double[] { cx / n, cy / n };
                if (addError)
                    xyz = addError(xyz, precisionInPixels, random);
                params[Gaussian2DFunction.X_POSITION] = (float) xyz[0];
                params[Gaussian2DFunction.Y_POSITION] = (float) xyz[1];
                params[Gaussian2DFunction.SIGNAL] = n;
                params[Gaussian2DFunction.X_SD] = params[Gaussian2DFunction.Y_SD] = 1;
                final float noise = 0.1f;
                PeakResult r = new ExtendedPeakResult(peak, (int) params[Gaussian2DFunction.X_POSITION], (int) params[Gaussian2DFunction.Y_POSITION], n, 0, noise, params, null, peak, id);
                results.add(r);
                if (last != null) {
                    sum += last.distance2(r);
                    count++;
                }
                last = r;
                n = 0;
                cx = cy = 0;
                peak++;
            }
            if (newId) {
                // Increment the frame so that tracing analysis can distinguish traces
                peak++;
                last = null;
                id = result.id;
            }
        }
        n++;
        cx += result.x;
        cy += result.y;
    }
    // Final peak
    if (n != 0) {
        final float[] params = new float[7];
        double[] xyz = new double[] { cx / n, cy / n };
        if (addError)
            xyz = addError(xyz, precisionInPixels, random);
        params[Gaussian2DFunction.X_POSITION] = (float) xyz[0];
        params[Gaussian2DFunction.Y_POSITION] = (float) xyz[1];
        params[Gaussian2DFunction.SIGNAL] = n;
        params[Gaussian2DFunction.X_SD] = params[Gaussian2DFunction.Y_SD] = 1;
        final float noise = 0.1f;
        PeakResult r = new ExtendedPeakResult(peak, (int) params[Gaussian2DFunction.X_POSITION], (int) params[Gaussian2DFunction.Y_POSITION], n, 0, noise, params, null, peak, id);
        results.add(r);
        if (last != null) {
            sum += last.distance2(r);
            count++;
        }
    }
    // MSD in pixels^2 / frame
    double msd = sum / count;
    // Convert to um^2/second
    Utils.log("Aggregated data D=%s um^2/s, Precision=%s nm, N=%d, step=%s s, mean=%s um^2, MSD = %s um^2/s", Utils.rounded(settings.diffusionRate), Utils.rounded(myPrecision), count, Utils.rounded(results.getCalibration().getExposureTime() / 1000), Utils.rounded(msd / conversionFactor), Utils.rounded((msd / conversionFactor) / (results.getCalibration().getExposureTime() / 1000)));
    msdAnalysis(points);
}

Example 82 with PeakResult

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

the class PSFCreator method fitPSF.

/**
	 * Fit the new PSF image and show a graph of the amplitude/width
	 * 
	 * @param psf
	 * @param loess
	 * @param averageRange
	 * @param fitCom
	 * @return The width of the PSF in the z-centre
	 */
private double fitPSF(ImageStack psf, LoessInterpolator loess, int cz, double averageRange, double[][] fitCom) {
    IJ.showStatus("Fitting final PSF");
    // is not appropriate for a normalised PSF. 
    if (fitConfig.getFitSolver() == FitSolver.MLE) {
        Utils.log("  Maximum Likelihood Estimation (MLE) is not appropriate for final PSF fitting.");
        Utils.log("  Switching to Least Square Estimation");
        fitConfig.setFitSolver(FitSolver.LVM);
        if (interactiveMode) {
            GlobalSettings settings = new GlobalSettings();
            settings.setFitEngineConfiguration(config);
            PeakFit.configureFitSolver(settings, null, false, false);
        }
    }
    // Update the box radius since this is used in the fitSpot method.
    boxRadius = psf.getWidth() / 2;
    int x = boxRadius, y = boxRadius;
    FitConfiguration fitConfig = config.getFitConfiguration();
    final double shift = fitConfig.getCoordinateShiftFactor();
    fitConfig.setInitialPeakStdDev0(fitConfig.getInitialPeakStdDev0() * magnification);
    fitConfig.setInitialPeakStdDev1(fitConfig.getInitialPeakStdDev1() * magnification);
    // Need to be updated after the widths have been set
    fitConfig.setCoordinateShiftFactor(shift);
    fitConfig.setBackgroundFitting(false);
    // Since the PSF will be normalised
    fitConfig.setMinPhotons(0);
    //fitConfig.setLog(new IJLogger());
    MemoryPeakResults results = fitSpot(psf, psf.getWidth(), psf.getHeight(), x, y);
    if (results.size() < 5) {
        Utils.log("  Final PSF: Not enough fit results %d", results.size());
        return 0;
    }
    // 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;
    // Set limits for the fit
    final float maxWidth = (float) (FastMath.max(fitConfig.getInitialPeakStdDev0(), fitConfig.getInitialPeakStdDev1()) * magnification * 4);
    // PSF is normalised to 1  
    final float maxSignal = 2;
    for (PeakResult peak : results.getResults()) {
        // Remove bad fits where the width/signal is above the expected
        final float w = FastMath.max(peak.getXSD(), peak.getYSD());
        if (peak.getSignal() > maxSignal || w > maxWidth)
            continue;
        z[i] = peak.getFrame();
        fitCom[0][peak.getFrame() - 1] = xCoord[i] = peak.getXPosition() - x;
        fitCom[1][peak.getFrame() - 1] = yCoord[i] = peak.getYPosition() - y;
        sd[i] = w;
        a[i] = peak.getAmplitude();
        i++;
    }
    // Truncate
    z = Arrays.copyOf(z, i);
    xCoord = Arrays.copyOf(xCoord, i);
    yCoord = Arrays.copyOf(yCoord, i);
    sd = Arrays.copyOf(sd, i);
    a = Arrays.copyOf(a, i);
    // Extract the average smoothed range from the individual fits
    int r = (int) Math.ceil(averageRange / 2);
    int start = 0, stop = z.length - 1;
    for (int j = 0; j < z.length; j++) {
        if (z[j] > cz - r) {
            start = j;
            break;
        }
    }
    for (int j = z.length; j-- > 0; ) {
        if (z[j] < cz + r) {
            stop = j;
            break;
        }
    }
    // 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[] smoothA = new double[smoothX.length];
    double[] newZ = new double[smoothX.length];
    int smoothCzIndex = 0;
    for (int j = start, k = 0; j <= stop; j++, k++) {
        smoothX[k] = xCoord[j];
        smoothY[k] = yCoord[j];
        smoothSd[k] = sd[j];
        smoothA[k] = a[j];
        newZ[k] = z[j];
        if (newZ[k] == cz)
            smoothCzIndex = k;
    }
    smoothX = loess.smooth(newZ, smoothX);
    smoothY = loess.smooth(newZ, smoothY);
    smoothSd = loess.smooth(newZ, smoothSd);
    smoothA = loess.smooth(newZ, smoothA);
    // Update the widths and positions using the magnification
    final double scale = 1.0 / magnification;
    for (int j = 0; j < xCoord.length; j++) {
        xCoord[j] *= scale;
        yCoord[j] *= scale;
        sd[j] *= scale;
    }
    for (int j = 0; j < smoothX.length; j++) {
        smoothX[j] *= scale;
        smoothY[j] *= scale;
        smoothSd[j] *= scale;
    }
    showPlots(z, a, newZ, smoothA, xCoord, yCoord, sd, newZ, smoothX, smoothY, smoothSd, cz);
    // Store the data for replotting
    this.z = z;
    this.a = a;
    this.smoothAz = newZ;
    this.smoothA = smoothA;
    this.xCoord = xCoord;
    this.yCoord = yCoord;
    this.sd = sd;
    this.newZ = newZ;
    this.smoothX = smoothX;
    this.smoothY = smoothY;
    this.smoothSd = smoothSd;
    //maximumIndex = findMinimumIndex(smoothSd, maximumIndex - start);
    return smoothSd[smoothCzIndex];
}

Example 83 with PeakResult

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

the class LoadLocalisations method run.

/*
	 * (non-Javadoc)
	 * 
	 * @see ij.plugin.PlugIn#run(java.lang.String)
	 */
public void run(String arg) {
    SMLMUsageTracker.recordPlugin(this.getClass(), arg);
    GlobalSettings globalSettings = SettingsManager.loadSettings();
    CreateDataSettings settings = globalSettings.getCreateDataSettings();
    String[] path = Utils.decodePath(settings.localisationsFilename);
    OpenDialog chooser = new OpenDialog("Localisations_File", path[0], path[1]);
    if (chooser.getFileName() == null)
        return;
    settings.localisationsFilename = chooser.getDirectory() + chooser.getFileName();
    SettingsManager.saveSettings(globalSettings);
    LocalisationList localisations = loadLocalisations(settings.localisationsFilename);
    if (localisations == null)
        // Cancelled
        return;
    if (localisations.isEmpty()) {
        IJ.error(TITLE, "No localisations could be loaded");
        return;
    }
    MemoryPeakResults results = localisations.toPeakResults();
    // Ask the user what depth to use to create the in-memory results
    if (!getZDepth(results))
        return;
    if (myLimitZ) {
        MemoryPeakResults results2 = new MemoryPeakResults(results.size());
        results.setName(name);
        results.copySettings(results);
        for (PeakResult peak : results.getResults()) {
            if (peak.error < minz || peak.error > maxz)
                continue;
            results2.add(peak);
        }
        results = results2;
    }
    // Create the in-memory results
    if (results.size() > 0) {
        MemoryPeakResults.addResults(results);
    }
    IJ.showStatus(String.format("Loaded %d localisations", results.size()));
    if (myLimitZ)
        Utils.log("Loaded %d localisations, z between %.2f - %.2f", results.size(), minz, maxz);
    else
        Utils.log("Loaded %d localisations", results.size());
}

Example 84 with PeakResult

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

the class SpotInspector method run.

/*
	 * (non-Javadoc)
	 * 
	 * @see ij.plugin.PlugIn#run(java.lang.String)
	 */
public void run(String arg) {
    SMLMUsageTracker.recordPlugin(this.getClass(), arg);
    if (MemoryPeakResults.isMemoryEmpty()) {
        IJ.error(TITLE, "No localisations in memory");
        return;
    }
    if (!showDialog())
        return;
    // Load the results
    results = ResultsManager.loadInputResults(inputOption, false);
    if (results == null || results.size() == 0) {
        IJ.error(TITLE, "No results could be loaded");
        IJ.showStatus("");
        return;
    }
    // Check if the original image is open
    ImageSource source = results.getSource();
    if (source == null) {
        IJ.error(TITLE, "Unknown original source image");
        return;
    }
    source = source.getOriginal();
    if (!source.open()) {
        IJ.error(TITLE, "Cannot open original source image: " + source.toString());
        return;
    }
    final float stdDevMax = getStandardDeviation(results);
    if (stdDevMax < 0) {
        // TODO - Add dialog to get the initial peak width
        IJ.error(TITLE, "Fitting configuration (for initial peak width) is not available");
        return;
    }
    // Rank spots
    rankedResults = new ArrayList<PeakResultRank>(results.size());
    final double a = results.getNmPerPixel();
    final double gain = results.getGain();
    final boolean emCCD = results.isEMCCD();
    for (PeakResult r : results.getResults()) {
        float[] score = getScore(r, a, gain, emCCD, stdDevMax);
        rankedResults.add(new PeakResultRank(r, score[0], score[1]));
    }
    Collections.sort(rankedResults);
    // Prepare results table. Get bias if necessary
    if (showCalibratedValues) {
        // Get a bias if required
        Calibration calibration = results.getCalibration();
        if (calibration.getBias() == 0) {
            ExtendedGenericDialog gd = new ExtendedGenericDialog(TITLE);
            gd.addMessage("Calibrated results requires a camera bias");
            gd.addNumericField("Camera_bias (ADUs)", calibration.getBias(), 2);
            gd.showDialog();
            if (!gd.wasCanceled()) {
                calibration.setBias(Math.abs(gd.getNextNumber()));
            }
        }
    }
    IJTablePeakResults table = new IJTablePeakResults(false, results.getName(), true);
    table.copySettings(results);
    table.setTableTitle(TITLE);
    table.setAddCounter(true);
    table.setShowCalibratedValues(showCalibratedValues);
    table.begin();
    // Add a mouse listener to jump to the frame for the clicked line
    textPanel = table.getResultsWindow().getTextPanel();
    // We must ignore old instances of this class from the mouse listeners
    id = ++currentId;
    textPanel.addMouseListener(this);
    // Add results to the table
    int n = 0;
    for (PeakResultRank rank : rankedResults) {
        rank.rank = n++;
        PeakResult r = rank.peakResult;
        table.add(r.getFrame(), r.origX, r.origY, r.origValue, r.error, r.noise, r.params, r.paramsStdDev);
    }
    table.end();
    if (plotScore || plotHistogram) {
        // Get values for the plots
        float[] xValues = null, yValues = null;
        double yMin, yMax;
        int spotNumber = 0;
        xValues = new float[rankedResults.size()];
        yValues = new float[xValues.length];
        for (PeakResultRank rank : rankedResults) {
            xValues[spotNumber] = spotNumber + 1;
            yValues[spotNumber++] = recoverScore(rank.score);
        }
        // Set the min and max y-values using 1.5 x IQR 
        DescriptiveStatistics stats = new DescriptiveStatistics();
        for (float v : yValues) stats.addValue(v);
        if (removeOutliers) {
            double lower = stats.getPercentile(25);
            double upper = stats.getPercentile(75);
            double iqr = upper - lower;
            yMin = FastMath.max(lower - iqr, stats.getMin());
            yMax = FastMath.min(upper + iqr, stats.getMax());
            IJ.log(String.format("Data range: %f - %f. Plotting 1.5x IQR: %f - %f", stats.getMin(), stats.getMax(), yMin, yMax));
        } else {
            yMin = stats.getMin();
            yMax = stats.getMax();
            IJ.log(String.format("Data range: %f - %f", yMin, yMax));
        }
        plotScore(xValues, yValues, yMin, yMax);
        plotHistogram(yValues, yMin, yMax);
    }
    // Extract spots into a stack
    final int w = source.getWidth();
    final int h = source.getHeight();
    final int size = 2 * radius + 1;
    ImageStack spots = new ImageStack(size, size, rankedResults.size());
    // To assist the extraction of data from the image source, process them in time order to allow 
    // frame caching. Then set the appropriate slice in the result stack
    Collections.sort(rankedResults, new Comparator<PeakResultRank>() {

        public int compare(PeakResultRank o1, PeakResultRank o2) {
            if (o1.peakResult.getFrame() < o2.peakResult.getFrame())
                return -1;
            if (o1.peakResult.getFrame() > o2.peakResult.getFrame())
                return 1;
            return 0;
        }
    });
    for (PeakResultRank rank : rankedResults) {
        PeakResult r = rank.peakResult;
        // Extract image
        // Note that the coordinates are relative to the middle of the pixel (0.5 offset)
        // so do not round but simply convert to int
        final int x = (int) (r.params[Gaussian2DFunction.X_POSITION]);
        final int y = (int) (r.params[Gaussian2DFunction.Y_POSITION]);
        // Extract a region but crop to the image bounds
        int minX = x - radius;
        int minY = y - radius;
        int maxX = FastMath.min(x + radius + 1, w);
        int maxY = FastMath.min(y + radius + 1, h);
        int padX = 0, padY = 0;
        if (minX < 0) {
            padX = -minX;
            minX = 0;
        }
        if (minY < 0) {
            padY = -minY;
            minY = 0;
        }
        int sizeX = maxX - minX;
        int sizeY = maxY - minY;
        float[] data = source.get(r.getFrame(), new Rectangle(minX, minY, sizeX, sizeY));
        // Prevent errors with missing data
        if (data == null)
            data = new float[sizeX * sizeY];
        ImageProcessor spotIp = new FloatProcessor(sizeX, sizeY, data, null);
        // Pad if necessary, i.e. the crop is too small for the stack
        if (padX > 0 || padY > 0 || sizeX < size || sizeY < size) {
            ImageProcessor spotIp2 = spotIp.createProcessor(size, size);
            spotIp2.insert(spotIp, padX, padY);
            spotIp = spotIp2;
        }
        int slice = rank.rank + 1;
        spots.setPixels(spotIp.getPixels(), slice);
        spots.setSliceLabel(Utils.rounded(rank.originalScore), slice);
    }
    source.close();
    ImagePlus imp = Utils.display(TITLE, spots);
    imp.setRoi((PointRoi) null);
    // Make bigger		
    for (int i = 10; i-- > 0; ) imp.getWindow().getCanvas().zoomIn(imp.getWidth() / 2, imp.getHeight() / 2);
}

Example 85 with PeakResult

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

the class SpotInspector method mouseClicked.

public void mouseClicked(MouseEvent e) {
    if (id != currentId)
        return;
    // Show the result that was double clicked in the result table
    if (e.getClickCount() > 1) {
        int rank = textPanel.getSelectionStart() + 1;
        // Show the spot that was double clicked
        ImagePlus imp = WindowManager.getImage(TITLE);
        if (imp != null && rank > 0 && rank <= imp.getStackSize()) {
            imp.setSlice(rank);
            if (imp.getWindow() != null)
                imp.getWindow().toFront();
            // Add an ROI to to the image containing all the spots in that part of the frame
            PeakResult r = rankedResults.get(rank - 1).peakResult;
            final int x = (int) (r.params[Gaussian2DFunction.X_POSITION]);
            final int y = (int) (r.params[Gaussian2DFunction.Y_POSITION]);
            // Find bounds
            int minX = x - radius;
            int minY = y - radius;
            int maxX = x + radius + 1;
            int maxY = y + radius + 1;
            // Create ROIs
            ArrayList<float[]> spots = new ArrayList<float[]>();
            for (PeakResult peak : results.getResults()) {
                if (peak.getXPosition() > minX && peak.getXPosition() < maxX && peak.getYPosition() > minY && peak.getYPosition() < maxY) {
                    // Use only unique points
                    final float xPosition = peak.getXPosition() - minX;
                    final float yPosition = peak.getYPosition() - minY;
                    if (contains(spots, xPosition, yPosition))
                        continue;
                    spots.add(new float[] { xPosition, yPosition });
                }
            }
            int points = spots.size();
            float[] ox = new float[points];
            float[] oy = new float[points];
            for (int i = 0; i < points; i++) {
                ox[i] = spots.get(i)[0];
                oy[i] = spots.get(i)[1];
            }
            imp.setRoi(new PointRoi(ox, oy, points));
        }
    }
}