Examples with Gaussian2DFitter uk.ac.sussex.gdsc.smlm.fitting.Gaussian2DFitter used on opensource projects

Example 1 with Gaussian2DFitter

use of uk.ac.sussex.gdsc.smlm.fitting.Gaussian2DFitter in project GDSC-SMLM by aherbert.

the class GaussianFit method runFinal.

/**
 * Perform fitting using the chosen maxima. Update the overlay if successful.
 *
 * @param ip The input image
 */
private void runFinal(ImageProcessor ip) {
    ip.reset();
    final Rectangle bounds = ip.getRoi();
    // Crop to the ROI
    final float[] data = ImageJImageConverter.getData(ip);
    final int width = bounds.width;
    final int height = bounds.height;
    // Sort the maxima
    float[] smoothData = data;
    if (getSmooth() > 0) {
        // Smoothing destructively modifies the data so create a copy
        smoothData = Arrays.copyOf(data, width * height);
        final BlockMeanFilter filter = new BlockMeanFilter();
        if (settings.smooth <= settings.border) {
            filter.stripedBlockFilterInternal(smoothData, width, height, (float) settings.smooth);
        } else {
            filter.stripedBlockFilter(smoothData, width, height, (float) settings.smooth);
        }
    }
    SortUtils.sortIndices(maxIndices, smoothData, true);
    // Show the candidate peaks
    if (maxIndices.length > 0) {
        final String message = String.format("Identified %d peaks", maxIndices.length);
        if (isLogProgress()) {
            IJ.log(message);
            for (final int index : maxIndices) {
                IJ.log(String.format("  %.2f @ [%d,%d]", data[index], bounds.x + index % width, bounds.y + index / width));
            }
        }
        // Check whether to run if the number of peaks is large
        if (maxIndices.length > 10) {
            final GenericDialog gd = new GenericDialog("Warning");
            gd.addMessage(message + "\nDo you want to fit?");
            gd.showDialog();
            if (gd.wasCanceled()) {
                return;
            }
        }
    } else {
        IJ.log("No maxima identified");
        return;
    }
    results = new ImageJTablePeakResults(settings.showDeviations, imp.getTitle() + " [" + imp.getCurrentSlice() + "]");
    final CalibrationWriter cw = new CalibrationWriter();
    cw.setIntensityUnit(IntensityUnit.COUNT);
    cw.setDistanceUnit(DistanceUnit.PIXEL);
    cw.setAngleUnit(AngleUnit.RADIAN);
    results.setCalibration(cw.getCalibration());
    results.setPsf(PsfProtosHelper.getDefaultPsf(getPsfType()));
    results.setShowFittingData(true);
    results.setAngleUnit(AngleUnit.DEGREE);
    results.begin();
    // Perform the Gaussian fit
    long ellapsed = 0;
    final FloatProcessor renderedImage = settings.showFit ? new FloatProcessor(ip.getWidth(), ip.getHeight()) : null;
    if (!settings.singleFit) {
        if (isLogProgress()) {
            IJ.log("Combined fit");
        }
        // Estimate height from smoothed data
        final double[] estimatedHeights = new double[maxIndices.length];
        for (int i = 0; i < estimatedHeights.length; i++) {
            estimatedHeights[i] = smoothData[maxIndices[i]];
        }
        final FitConfiguration config = new FitConfiguration();
        setupPeakFiltering(config);
        final long time = System.nanoTime();
        final double[] params = fitMultiple(data, width, height, maxIndices, estimatedHeights);
        ellapsed = System.nanoTime() - time;
        if (params != null) {
            // Copy all the valid parameters into a new array
            final double[] validParams = new double[params.length];
            int count = 0;
            int validPeaks = 0;
            validParams[count++] = params[0];
            final double[] initialParams = convertParameters(fitResult.getInitialParameters());
            final double[] paramsDev = convertParameters(fitResult.getParameterDeviations());
            final Rectangle regionBounds = new Rectangle();
            final float[] xpoints = new float[maxIndices.length];
            final float[] ypoints = new float[maxIndices.length];
            int npoints = 0;
            for (int i = 1, n = 0; i < params.length; i += Gaussian2DFunction.PARAMETERS_PER_PEAK, n++) {
                final int y = maxIndices[n] / width;
                final int x = maxIndices[n] % width;
                // Check the peak is a good fit
                if (settings.filterResults && config.validatePeak(n, initialParams, params, paramsDev) != FitStatus.OK) {
                    continue;
                }
                if (settings.showFit) {
                    // Copy the valid parameters before there are adjusted to global bounds
                    validPeaks++;
                    for (int ii = i, j = 0; j < Gaussian2DFunction.PARAMETERS_PER_PEAK; ii++, j++) {
                        validParams[count++] = params[ii];
                    }
                }
                final double[] peakParams = extractParams(params, i);
                final double[] peakParamsDev = extractParams(paramsDev, i);
                addResult(bounds, regionBounds, peakParams, peakParamsDev, npoints, x, y, data[maxIndices[n]]);
                // Add fit result to the overlay - Coords are updated with the region offsets in addResult
                final double xf = peakParams[Gaussian2DFunction.X_POSITION];
                final double yf = peakParams[Gaussian2DFunction.Y_POSITION];
                xpoints[npoints] = (float) xf;
                ypoints[npoints] = (float) yf;
                npoints++;
            }
            setOverlay(npoints, xpoints, ypoints);
            // Draw the fit
            if (validPeaks != 0) {
                addToImage(bounds.x, bounds.y, renderedImage, validParams, validPeaks, width, height);
            }
        } else {
            if (isLogProgress()) {
                IJ.log("Failed to fit " + TextUtils.pleural(maxIndices.length, "peak") + ": " + getReason(fitResult));
            }
            imp.setOverlay(null);
        }
    } else {
        if (isLogProgress()) {
            IJ.log("Individual fit");
        }
        int npoints = 0;
        final float[] xpoints = new float[maxIndices.length];
        final float[] ypoints = new float[maxIndices.length];
        // Extract each peak and fit individually
        final ImageExtractor ie = ImageExtractor.wrap(data, width, height);
        float[] region = null;
        final Gaussian2DFitter gf = createGaussianFitter(settings.filterResults);
        double[] validParams = null;
        final ShortProcessor renderedImageCount = settings.showFit ? new ShortProcessor(ip.getWidth(), ip.getHeight()) : null;
        for (int n = 0; n < maxIndices.length; n++) {
            final int y = maxIndices[n] / width;
            final int x = maxIndices[n] % width;
            final long time = System.nanoTime();
            final Rectangle regionBounds = ie.getBoxRegionBounds(x, y, settings.singleRegionSize);
            region = ie.crop(regionBounds, region);
            final int newIndex = (y - regionBounds.y) * regionBounds.width + x - regionBounds.x;
            if (isLogProgress()) {
                IJ.log("Fitting peak " + (n + 1));
            }
            final double[] peakParams = fitSingle(gf, region, regionBounds.width, regionBounds.height, newIndex, smoothData[maxIndices[n]]);
            ellapsed += System.nanoTime() - time;
            // Output fit result
            if (peakParams != null) {
                if (settings.showFit) {
                    // Copy the valid parameters before there are adjusted to global bounds
                    validParams = peakParams.clone();
                }
                double[] peakParamsDev = null;
                if (settings.showDeviations) {
                    peakParamsDev = convertParameters(fitResult.getParameterDeviations());
                }
                addResult(bounds, regionBounds, peakParams, peakParamsDev, n, x, y, data[maxIndices[n]]);
                // Add fit result to the overlay - Coords are updated with the region offsets in addResult
                final double xf = peakParams[Gaussian2DFunction.X_POSITION];
                final double yf = peakParams[Gaussian2DFunction.Y_POSITION];
                xpoints[npoints] = (float) xf;
                ypoints[npoints] = (float) yf;
                npoints++;
                // Draw the fit
                if (settings.showDeviations) {
                    final int ox = bounds.x + regionBounds.x;
                    final int oy = bounds.y + regionBounds.y;
                    addToImage(ox, oy, renderedImage, validParams, 1, regionBounds.width, regionBounds.height);
                    addCount(ox, oy, renderedImageCount, regionBounds.width, regionBounds.height);
                }
            } else if (isLogProgress()) {
                IJ.log("Failed to fit peak " + (n + 1) + ": " + getReason(fitResult));
            }
        }
        // Update the overlay
        if (npoints > 0) {
            setOverlay(npoints, xpoints, ypoints);
        } else {
            imp.setOverlay(null);
        }
        // Create the mean
        if (settings.showFit) {
            for (int i = renderedImageCount.getPixelCount(); i-- > 0; ) {
                final int count = renderedImageCount.get(i);
                if (count > 1) {
                    renderedImage.setf(i, renderedImage.getf(i) / count);
                }
            }
        }
    }
    results.end();
    if (renderedImage != null) {
        ImageJUtils.display(TITLE, renderedImage);
    }
    if (isLogProgress()) {
        IJ.log("Time = " + (ellapsed / 1000000.0) + "ms");
    }
}

Example 2 with Gaussian2DFitter

use of uk.ac.sussex.gdsc.smlm.fitting.Gaussian2DFitter in project GDSC-SMLM by aherbert.

the class SpotAnalysis method updateCurrentSlice.

private void updateCurrentSlice(int slice) {
    if (slice != currentSlice) {
        currentSlice = slice;
        final double signal = getSignal(slice);
        final double noise = smoothSd[slice - 1];
        currentLabel.setText(String.format("Frame %d: Signal = %s, SNR = %s", slice, MathUtils.rounded(signal, 4), MathUtils.rounded(signal / noise, 3)));
        drawProfiles();
        // Fit the PSF using a Gaussian
        final FitConfiguration fitConfiguration = new FitConfiguration();
        fitConfiguration.setPsf(PsfProtosHelper.defaultOneAxisGaussian2DPSF);
        fitConfiguration.setFixedPsf(true);
        fitConfiguration.setBackgroundFitting(true);
        fitConfiguration.setSignalStrength(0);
        fitConfiguration.setCoordinateShift(rawImp.getWidth() / 4.0f);
        fitConfiguration.setComputeResiduals(false);
        fitConfiguration.setComputeDeviations(false);
        final Gaussian2DFitter gf = new Gaussian2DFitter(fitConfiguration);
        double[] params = new double[1 + Gaussian2DFunction.PARAMETERS_PER_PEAK];
        final double psfWidth = Double.parseDouble(widthTextField.getText());
        params[Gaussian2DFunction.BACKGROUND] = smoothMean[slice - 1];
        params[Gaussian2DFunction.SIGNAL] = (gain * signal);
        params[Gaussian2DFunction.X_POSITION] = rawImp.getWidth() / 2.0f;
        params[Gaussian2DFunction.Y_POSITION] = rawImp.getHeight() / 2.0f;
        params[Gaussian2DFunction.X_SD] = params[Gaussian2DFunction.Y_SD] = psfWidth;
        float[] data = (float[]) rawImp.getImageStack().getProcessor(slice).getPixels();
        FitResult fitResult = gf.fit(SimpleArrayUtils.toDouble(data), rawImp.getWidth(), rawImp.getHeight(), 1, params, new boolean[1]);
        if (fitResult.getStatus() == FitStatus.OK) {
            params = fitResult.getParameters();
            final double spotSignal = params[Gaussian2DFunction.SIGNAL] / gain;
            rawFittedLabel.setText(String.format("Raw fit: Signal = %s, SNR = %s", MathUtils.rounded(spotSignal, 4), MathUtils.rounded(spotSignal / noise, 3)));
            ImageRoiPainter.addRoi(rawImp, slice, new OffsetPointRoi(params[Gaussian2DFunction.X_POSITION], params[Gaussian2DFunction.Y_POSITION]));
        } else {
            rawFittedLabel.setText("");
            rawImp.setOverlay(null);
        }
        // Fit the PSF using a Gaussian
        if (blurImp == null) {
            return;
        }
        params = new double[1 + Gaussian2DFunction.PARAMETERS_PER_PEAK];
        params[Gaussian2DFunction.BACKGROUND] = (float) smoothMean[slice - 1];
        params[Gaussian2DFunction.SIGNAL] = (float) (gain * signal);
        params[Gaussian2DFunction.X_POSITION] = rawImp.getWidth() / 2.0f;
        params[Gaussian2DFunction.Y_POSITION] = rawImp.getHeight() / 2.0f;
        params[Gaussian2DFunction.X_SD] = params[Gaussian2DFunction.Y_SD] = psfWidth;
        data = (float[]) blurImp.getImageStack().getProcessor(slice).getPixels();
        fitResult = gf.fit(SimpleArrayUtils.toDouble(data), rawImp.getWidth(), rawImp.getHeight(), 1, params, new boolean[1]);
        if (fitResult.getStatus() == FitStatus.OK) {
            params = fitResult.getParameters();
            final double spotSignal = params[Gaussian2DFunction.SIGNAL] / gain;
            blurFittedLabel.setText(String.format("Blur fit: Signal = %s, SNR = %s", MathUtils.rounded(spotSignal, 4), MathUtils.rounded(spotSignal / noise, 3)));
            ImageRoiPainter.addRoi(blurImp, slice, new OffsetPointRoi(params[Gaussian2DFunction.X_POSITION], params[Gaussian2DFunction.Y_POSITION]));
        } else {
            blurFittedLabel.setText("");
            blurImp.setOverlay(null);
        }
    }
}

Example 3 with Gaussian2DFitter

use of uk.ac.sussex.gdsc.smlm.fitting.Gaussian2DFitter in project GDSC-SMLM by aherbert.

the class GaussianFit method fit.

/**
 * Fits a single 2D Gaussian to the data. The fit is initialised at the highest value and then
 * optimised.
 *
 * <p>Data must be arranged in yx block order, i.e. height rows of width.
 *
 * <p>The angle parameter is only set if using elliptical Gaussian fitting.
 *
 * <p>Note: The returned fit coordinates should be offset by 0.5 if the input data represents
 * pixels
 *
 * @param data the data
 * @param width the width
 * @param height the height
 * @return Array containing the fitted curve data: Background, Amplitude, PosX, PosY, StdDevX,
 *         StdDevY, Angle. Null if no fit is possible.
 */
@Nullable
public double[] fit(float[] data, int width, int height) {
    if (data == null || data.length != width * height) {
        return null;
    }
    // Get the limits
    float max = Float.MIN_VALUE;
    int maxIndex = -1;
    for (int i = data.length; i-- > 0; ) {
        final float f = data[i];
        if (max < f) {
            max = f;
            maxIndex = i;
        }
    }
    if (maxIndex < 0) {
        return null;
    }
    final Gaussian2DFitter gf = createGaussianFitter(false);
    final FitResult fitResult = gf.fit(SimpleArrayUtils.toDouble(data), width, height, new int[] { maxIndex });
    if (fitResult.getStatus() == FitStatus.OK) {
        chiSquared = fitResult.getError();
        final double[] params = fitResult.getParameters();
        // Check bounds
        final double x = params[Gaussian2DFunction.X_POSITION];
        final double y = params[Gaussian2DFunction.Y_POSITION];
        if (x < 0 || x >= width || y < 0 || y >= height) {
            return null;
        }
        // Re-arrange order for backwards compatibility with old code.
        final double background = params[Gaussian2DFunction.BACKGROUND];
        final double intensity = params[Gaussian2DFunction.SIGNAL];
        final double sx = params[Gaussian2DFunction.X_SD];
        final double sy = params[Gaussian2DFunction.Y_SD];
        final double angle = params[Gaussian2DFunction.ANGLE];
        final double amplitude = Gaussian2DPeakResultHelper.getAmplitude(intensity, sx, sy);
        return new double[] { background, amplitude, x, y, sx, sy, angle };
    }
    return null;
}

Example 4 with Gaussian2DFitter

use of uk.ac.sussex.gdsc.smlm.fitting.Gaussian2DFitter in project GDSC-SMLM by aherbert.

the class GaussianFit method createGaussianFitter.

private Gaussian2DFitter createGaussianFitter(boolean simpleFiltering) {
    final FitConfiguration config = new FitConfiguration();
    config.setFitSolver(FitSolver.LVM_LSE);
    config.setPsf(PsfProtosHelper.getDefaultPsf(getPsfType()));
    config.setMaxIterations(getMaxIterations());
    config.setRelativeThreshold(settings.relativeThreshold);
    config.setAbsoluteThreshold(settings.absoluteThreshold);
    config.setInitialPeakStdDev(getInitialPeakStdDev());
    config.setComputeDeviations(settings.showDeviations);
    // Set-up peak filtering only for single fitting
    config.setDisableSimpleFilter(!simpleFiltering);
    setupPeakFiltering(config);
    if (isLogProgress()) {
        config.setLog(ImageJPluginLoggerHelper.getLogger(getClass()));
    }
    config.setBackgroundFitting(settings.fitBackground);
    return new Gaussian2DFitter(config);
}

Example 5 with Gaussian2DFitter

use of uk.ac.sussex.gdsc.smlm.fitting.Gaussian2DFitter in project GDSC-SMLM by aherbert.

the class GaussianFit method fitMultiple.

/**
 * Fits a 2D Gaussian to the given data. Fits all the specified peaks.
 *
 * <p>Data must be arranged in yx block order, i.e. height rows of width.
 *
 * <p>Note: The fit coordinates should be offset by 0.5 if the input data represents pixels
 *
 * @param data the data
 * @param width the width
 * @param height the height
 * @param maxIndices Indices of the data to fit
 * @param estimatedHeights Estimated heights for the peaks (input from smoothed data)
 * @return Array containing the fitted curve data: The first value is the Background. The
 *         remaining values are Amplitude, PosX, PosY, StdDevX, StdDevY for each fitted peak. If
 *         elliptical fitting is performed the values are Amplitude, Angle, PosX, PosY, StdDevX,
 *         StdDevY for each fitted peak Null if no fit is possible.
 */
@Nullable
private double[] fitMultiple(float[] data, int width, int height, int[] maxIndices, double[] estimatedHeights) {
    if (data == null || data.length != width * height) {
        return null;
    }
    if (maxIndices == null || maxIndices.length == 0) {
        return null;
    }
    final Gaussian2DFitter gf = createGaussianFitter(false);
    this.fitResult = gf.fit(SimpleArrayUtils.toDouble(data), width, height, maxIndices, estimatedHeights);
    if (fitResult.getStatus() == FitStatus.OK) {
        chiSquared = fitResult.getError();
        final double[] params = fitResult.getParameters();
        convertParameters(params);
        return params;
    }
    return null;
}