Examples with EllipticalGaussian2DFunction gdsc.smlm.function.gaussian.EllipticalGaussian2DFunction used on opensource projects

Example 1 with EllipticalGaussian2DFunction

use of gdsc.smlm.function.gaussian.EllipticalGaussian2DFunction in project GDSC-SMLM by aherbert.

the class EllipticalGaussian2DFunctionTest method init.

protected void init() {
    flags = GaussianFunctionFactory.FIT_SIMPLE_ELLIPTICAL;
    f1 = new EllipticalGaussian2DFunction(1, maxx, maxy);
    f2 = new EllipticalGaussian2DFunction(2, maxx, maxy);
}

Example 2 with EllipticalGaussian2DFunction

use of gdsc.smlm.function.gaussian.EllipticalGaussian2DFunction in project GDSC-SMLM by aherbert.

the class GradientCalculatorSpeedTest method mleCalculatorComputesLogLikelihoodRatio.

@Test
public void mleCalculatorComputesLogLikelihoodRatio() {
    EllipticalGaussian2DFunction func = new EllipticalGaussian2DFunction(1, blockWidth, blockWidth);
    int n = blockWidth * blockWidth;
    double[] a = new double[7];
    rdg = new RandomDataGenerator(new Well19937c(30051977));
    for (int run = 5; run-- > 0; ) {
        a[0] = random(Background);
        a[1] = random(Amplitude);
        a[2] = random(Angle);
        a[3] = random(Xpos);
        a[4] = random(Ypos);
        a[5] = random(Xwidth);
        a[6] = random(Ywidth);
        // Simulate Poisson process
        func.initialise(a);
        double[] x = Utils.newArray(n, 0, 1.0);
        double[] u = new double[x.length];
        for (int i = 0; i < n; i++) {
            u[i] = func.eval(i);
            if (u[i] > 0)
                x[i] = rdg.nextPoisson(u[i]);
        }
        GradientCalculator calc = GradientCalculatorFactory.newCalculator(func.getNumberOfGradients(), true);
        double[][] alpha = new double[7][7];
        double[] beta = new double[7];
        double llr = PoissonCalculator.logLikelihoodRatio(u, x);
        double llr2 = calc.findLinearised(n, x, a, alpha, beta, func);
        //System.out.printf("llr=%f, llr2=%f\n", llr, llr2);
        Assert.assertEquals("Log-likelihood ratio", llr, llr2, llr * 1e-10);
    }
}

Example 3 with EllipticalGaussian2DFunction

use of gdsc.smlm.function.gaussian.EllipticalGaussian2DFunction in project GDSC-SMLM by aherbert.

the class GradientCalculatorSpeedTest method doubleCreateGaussianData.

/**
	 * Create random elliptical Gaussian data an returns the data plus an estimate of the parameters.
	 * Only the chosen parameters are randomised and returned for a maximum of (background, amplitude, angle, xpos,
	 * ypos, xwidth, ywidth }
	 *
	 * @param npeaks
	 *            the npeaks
	 * @param params
	 *            set on output
	 * @param randomiseParams
	 *            Set to true to randomise the params
	 * @return the double[]
	 */
private double[] doubleCreateGaussianData(int npeaks, double[] params, boolean randomiseParams) {
    int n = blockWidth * blockWidth;
    // Generate a 2D Gaussian
    EllipticalGaussian2DFunction func = new EllipticalGaussian2DFunction(npeaks, blockWidth, blockWidth);
    params[0] = random(Background);
    for (int i = 0, j = 1; i < npeaks; i++, j += 6) {
        params[j] = random(Amplitude);
        params[j + 1] = random(Angle);
        params[j + 2] = random(Xpos);
        params[j + 3] = random(Ypos);
        params[j + 4] = random(Xwidth);
        params[j + 5] = random(Ywidth);
    }
    double[] dy_da = new double[params.length];
    double[] y = new double[n];
    func.initialise(params);
    for (int i = 0; i < y.length; i++) {
        // Add random Poisson noise
        y[i] = rdg.nextPoisson(func.eval(i, dy_da));
    }
    if (randomiseParams) {
        params[0] = random(params[0]);
        for (int i = 0, j = 1; i < npeaks; i++, j += 6) {
            params[j] = random(params[j]);
            params[j + 1] = random(params[j + 1]);
            params[j + 2] = random(params[j + 2]);
            params[j + 3] = random(params[j + 3]);
            params[j + 4] = random(params[j + 4]);
            //params[j + 4];
            params[j + 5] = random(params[j + 5]);
        }
    }
    return y;
}

Example 4 with EllipticalGaussian2DFunction

use of gdsc.smlm.function.gaussian.EllipticalGaussian2DFunction in project GDSC-SMLM by aherbert.

the class SpeedTest method doubleCreateGaussianData.

/**
	 * Create random elliptical Gaussian data an returns the data plus an estimate of the parameters.
	 * Only the chosen parameters are randomised and returned for a maximum of (background, amplitude, angle, xpos,
	 * ypos, xwidth, ywidth }
	 * 
	 * @param params
	 *            set on output
	 * @return
	 */
private double[] doubleCreateGaussianData(int npeaks, double[] params) {
    int n = blockWidth * blockWidth;
    // Generate a 2D Gaussian
    EllipticalGaussian2DFunction func = new EllipticalGaussian2DFunction(npeaks, blockWidth, blockWidth);
    params[0] = Background + rand.nextFloat() * 5f;
    for (int i = 0, j = 1; i < npeaks; i++, j += 6) {
        params[j] = Amplitude + rand.nextFloat() * 5f;
        //(double) (Math.PI / 4.0); // Angle
        params[j + 1] = 0f;
        params[j + 2] = Xpos + rand.nextFloat() * 2f;
        params[j + 3] = Ypos + rand.nextFloat() * 2f;
        params[j + 4] = Xwidth + rand.nextFloat() * 2f;
        params[j + 5] = params[j + 4];
    }
    double[] dy_da = new double[params.length];
    double[] y = new double[n];
    func.initialise(params);
    for (int i = 0; i < y.length; i++) {
        // Add random noise
        y[i] = func.eval(i, dy_da) + ((rand.nextFloat() < 0.5f) ? -rand.nextFloat() * 5f : rand.nextFloat() * 5f);
    }
    // Randomise only the necessary parameters (i.e. not angle and X & Y widths should be the same)
    params[0] += ((rand.nextFloat() < 0.5f) ? -rand.nextFloat() : rand.nextFloat());
    for (int i = 0, j = 1; i < npeaks; i++, j += 6) {
        params[j + 1] += ((rand.nextFloat() < 0.5f) ? -rand.nextFloat() : rand.nextFloat());
        params[j + 3] += ((rand.nextFloat() < 0.5f) ? -rand.nextFloat() : rand.nextFloat());
        params[j + 4] += ((rand.nextFloat() < 0.5f) ? -rand.nextFloat() : rand.nextFloat());
        params[j + 5] = params[j + 4];
    }
    return y;
}

Example 5 with EllipticalGaussian2DFunction

use of gdsc.smlm.function.gaussian.EllipticalGaussian2DFunction 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();
    Rectangle bounds = ip.getRoi();
    // Crop to the ROI
    float[] data = ImageConverter.getData(ip);
    int width = bounds.width;
    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);
        AverageFilter filter = new AverageFilter();
        //filter.blockAverage(smoothData, width, height, smooth);
        if (smooth <= border)
            filter.stripedBlockAverageInternal(smoothData, width, height, (float) smooth);
        else
            filter.stripedBlockAverage(smoothData, width, height, (float) smooth);
    }
    Sort.sort(maxIndices, smoothData);
    // Show the candidate peaks
    if (maxIndices.length > 0) {
        String message = String.format("Identified %d peaks", maxIndices.length);
        if (isLogProgress()) {
            IJ.log(message);
            for (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) {
            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 IJTablePeakResults(showDeviations, imp.getTitle() + " [" + imp.getCurrentSlice() + "]");
    results.begin();
    // Perform the Gaussian fit
    long ellapsed = 0;
    if (!singleFit) {
        if (isLogProgress())
            IJ.log("Combined fit");
        // Estimate height from smoothed data
        double[] estimatedHeights = new double[maxIndices.length];
        for (int i = 0; i < estimatedHeights.length; i++) estimatedHeights[i] = smoothData[maxIndices[i]];
        FitConfiguration config = new FitConfiguration();
        setupPeakFiltering(config);
        long time = System.nanoTime();
        double[] params = fitMultiple(data, width, height, maxIndices, estimatedHeights);
        ellapsed = System.nanoTime() - time;
        if (params != null) {
            // Copy all the valid parameters into a new array
            double[] validParams = new double[params.length];
            int c = 0;
            int validPeaks = 0;
            validParams[c++] = params[0];
            double[] initialParams = convertParameters(fitResult.getInitialParameters());
            double[] paramsDev = convertParameters(fitResult.getParameterStdDev());
            Rectangle regionBounds = new Rectangle();
            int[] xpoints = new int[maxIndices.length];
            int[] ypoints = new int[maxIndices.length];
            int nMaxima = 0;
            for (int i = 1, n = 0; i < params.length; i += 6, n++) {
                int y = maxIndices[n] / width;
                int x = maxIndices[n] % width;
                // Check the peak is a good fit
                if (filterResults && config.validatePeak(n, initialParams, params) != FitStatus.OK)
                    continue;
                if (showFit) {
                    // Copy the valid parameters
                    validPeaks++;
                    for (int ii = i, j = 0; j < 6; ii++, j++) validParams[c++] = params[ii];
                }
                double[] peakParams = extractParams(params, i);
                double[] peakParamsDev = extractParams(paramsDev, i);
                addResult(bounds, regionBounds, data, peakParams, peakParamsDev, nMaxima, x, y, data[maxIndices[n]]);
                // Add fit result to the overlay - Coords are updated with the region offsets in addResult
                double xf = peakParams[3];
                double yf = peakParams[4];
                xpoints[nMaxima] = (int) (xf + 0.5);
                ypoints[nMaxima] = (int) (yf + 0.5);
                nMaxima++;
            }
            setOverlay(nMaxima, xpoints, ypoints);
            // Draw the fit
            if (showFit && validPeaks != 0) {
                double[] pixels = new double[data.length];
                EllipticalGaussian2DFunction f = new EllipticalGaussian2DFunction(validPeaks, width, height);
                invertParameters(validParams);
                f.initialise(validParams);
                for (int x = 0; x < pixels.length; x++) pixels[x] = f.eval(x);
                FloatProcessor fp = new FloatProcessor(width, height, pixels);
                // Insert into a full size image
                FloatProcessor fp2 = new FloatProcessor(ip.getWidth(), ip.getHeight());
                fp2.insert(fp, bounds.x, bounds.y);
                Utils.display(TITLE, fp2);
            }
        } else {
            if (isLogProgress()) {
                IJ.log("Failed to fit " + Utils.pleural(maxIndices.length, "peak") + getReason(fitResult));
            }
            imp.setOverlay(null);
        }
    } else {
        if (isLogProgress())
            IJ.log("Individual fit");
        int nMaxima = 0;
        int[] xpoints = new int[maxIndices.length];
        int[] ypoints = new int[maxIndices.length];
        // Extract each peak and fit individually
        ImageExtractor ie = new ImageExtractor(data, width, height);
        float[] region = null;
        Gaussian2DFitter gf = createGaussianFitter(filterResults);
        for (int n = 0; n < maxIndices.length; n++) {
            int y = maxIndices[n] / width;
            int x = maxIndices[n] % width;
            long time = System.nanoTime();
            Rectangle regionBounds = ie.getBoxRegionBounds(x, y, singleRegionSize);
            region = ie.crop(regionBounds, region);
            int newIndex = (y - regionBounds.y) * regionBounds.width + x - regionBounds.x;
            if (isLogProgress()) {
                IJ.log("Fitting peak " + (n + 1));
            }
            double[] peakParams = fitSingle(gf, region, regionBounds.width, regionBounds.height, newIndex, smoothData[maxIndices[n]]);
            ellapsed += System.nanoTime() - time;
            // Output fit result
            if (peakParams != null) {
                double[] peakParamsDev = null;
                if (showDeviations) {
                    peakParamsDev = convertParameters(fitResult.getParameterStdDev());
                }
                addResult(bounds, regionBounds, data, peakParams, peakParamsDev, n, x, y, data[maxIndices[n]]);
                // Add fit result to the overlay - Coords are updated with the region offsets in addResult
                double xf = peakParams[3];
                double yf = peakParams[4];
                xpoints[nMaxima] = (int) (xf + 0.5);
                ypoints[nMaxima] = (int) (yf + 0.5);
                nMaxima++;
            } else {
                if (isLogProgress()) {
                    IJ.log("Failed to fit peak " + (n + 1) + getReason(fitResult));
                }
            }
        }
        // Update the overlay
        if (nMaxima > 0)
            setOverlay(nMaxima, xpoints, ypoints);
        else
            imp.setOverlay(null);
    }
    results.end();
    if (isLogProgress())
        IJ.log("Time = " + (ellapsed / 1000000.0) + "ms");
}