Examples with Gaussian org.apache.commons.math3.analysis.function.Gaussian used on opensource projects

Example 56 with Gaussian

use of org.apache.commons.math3.analysis.function.Gaussian in project GDSC-SMLM by aherbert.

the class PcPalmMolecules method addToPlot.

/**
 * Add the skewed gaussian to the histogram plot.
 *
 * @param plot the plot
 * @param x the x
 * @param parameters Gaussian parameters
 * @param shape the shape
 */
private static void addToPlot(Plot plot, float[] x, double[] parameters, int shape) {
    final SkewNormalFunction sn = new SkewNormalFunction(parameters);
    final float[] y = new float[x.length];
    for (int i = 0; i < x.length; i++) {
        y[i] = (float) sn.evaluate(x[i]);
    }
    plot.addPoints(x, y, shape);
}

Example 57 with Gaussian

use of org.apache.commons.math3.analysis.function.Gaussian in project GDSC-SMLM by aherbert.

the class PcPalmMolecules method runSimulation.

private void runSimulation(boolean resultsAvailable) {
    if (resultsAvailable && !showSimulationDialog()) {
        return;
    }
    startLog();
    log("Simulation parameters");
    if (settings.blinkingDistribution == 3) {
        log("  - Clusters = %d", settings.numberOfMolecules);
        log("  - Simulation size = %s um", MathUtils.rounded(settings.simulationSize, 4));
        log("  - Molecules/cluster = %s", MathUtils.rounded(settings.blinkingRate, 4));
        log("  - Blinking distribution = %s", Settings.BLINKING_DISTRIBUTION[settings.blinkingDistribution]);
        log("  - p-Value = %s", MathUtils.rounded(settings.pvalue, 4));
    } else {
        log("  - Molecules = %d", settings.numberOfMolecules);
        log("  - Simulation size = %s um", MathUtils.rounded(settings.simulationSize, 4));
        log("  - Blinking rate = %s", MathUtils.rounded(settings.blinkingRate, 4));
        log("  - Blinking distribution = %s", Settings.BLINKING_DISTRIBUTION[settings.blinkingDistribution]);
    }
    log("  - Average precision = %s nm", MathUtils.rounded(settings.sigmaS, 4));
    log("  - Clusters simulation = " + Settings.CLUSTER_SIMULATION[settings.clusterSimulation]);
    if (settings.clusterSimulation > 0) {
        log("  - Cluster number = %s +/- %s", MathUtils.rounded(settings.clusterNumber, 4), MathUtils.rounded(settings.clusterNumberStdDev, 4));
        log("  - Cluster radius = %s nm", MathUtils.rounded(settings.clusterRadius, 4));
    }
    final double nmPerPixel = 100;
    final double width = settings.simulationSize * 1000.0;
    final UniformRandomProvider rng = UniformRandomProviders.create();
    final UniformDistribution dist = new UniformDistribution(null, new double[] { width, width, 0 }, rng.nextInt());
    final NormalizedGaussianSampler gauss = SamplerUtils.createNormalizedGaussianSampler(rng);
    settings.molecules = new ArrayList<>(settings.numberOfMolecules);
    // Create some dummy results since the calibration is required for later analysis
    settings.results = new MemoryPeakResults(PsfHelper.create(PSFType.CUSTOM));
    settings.results.setCalibration(CalibrationHelper.create(nmPerPixel, 1, 100));
    settings.results.setSource(new NullSource("Molecule Simulation"));
    settings.results.begin();
    int count = 0;
    // Generate a sequence of coordinates
    final ArrayList<double[]> xyz = new ArrayList<>((int) (settings.numberOfMolecules * 1.1));
    final Statistics statsRadius = new Statistics();
    final Statistics statsSize = new Statistics();
    final String maskTitle = TITLE + " Cluster Mask";
    ByteProcessor bp = null;
    double maskScale = 0;
    if (settings.clusterSimulation > 0) {
        // Simulate clusters.
        // Note: In the Veatch et al. paper (Plos 1, e31457) correlation functions are built using
        // circles with small radii of 4-8 Arbitrary Units (AU) or large radii of 10-30 AU. A
        // fluctuations model is created at T = 1.075 Tc. It is not clear exactly how the particles
        // are distributed.
        // It may be that a mask is created first using the model. The particles are placed on the
        // mask using a specified density. This simulation produces a figure to show either a damped
        // cosine function (circles) or an exponential (fluctuations). The number of particles in
        // each circle may be randomly determined just by density. The figure does not discuss the
        // derivation of the cluster size statistic.
        // 
        // If this plugin simulation is run with a uniform distribution and blinking rate of 1 then
        // the damped cosine function is reproduced. The curve crosses g(r)=1 at a value equivalent
        // to the average distance to the centre-of-mass of each drawn cluster, not the input cluster
        // radius parameter (which is a hard upper limit on the distance to centre).
        final int maskSize = settings.lowResolutionImageSize;
        int[] mask = null;
        // scale is in nm/pixel
        maskScale = width / maskSize;
        final ArrayList<double[]> clusterCentres = new ArrayList<>();
        int totalSteps = 1 + (int) Math.ceil(settings.numberOfMolecules / settings.clusterNumber);
        if (settings.clusterSimulation == 2 || settings.clusterSimulation == 3) {
            // Clusters are non-overlapping circles
            // Ensure the circles do not overlap by using an exclusion mask that accumulates
            // out-of-bounds pixels by drawing the last cluster (plus some border) on an image. When no
            // more pixels are available then stop generating molecules.
            // This is done by cumulatively filling a mask and using the MaskDistribution to select
            // a new point. This may be slow but it works.
            // TODO - Allow clusters of different sizes...
            mask = new int[maskSize * maskSize];
            Arrays.fill(mask, 255);
            MaskDistribution maskDistribution = new MaskDistribution(mask, maskSize, maskSize, 0, maskScale, maskScale, rng);
            double[] centre;
            IJ.showStatus("Computing clusters mask");
            final int roiRadius = (int) Math.round((settings.clusterRadius * 2) / maskScale);
            if (settings.clusterSimulation == 3) {
                // Generate a mask of circles then sample from that.
                // If we want to fill the mask completely then adjust the total steps to be the number of
                // circles that can fit inside the mask.
                totalSteps = (int) (maskSize * maskSize / (Math.PI * MathUtils.pow2(settings.clusterRadius / maskScale)));
            }
            while ((centre = maskDistribution.next()) != null && clusterCentres.size() < totalSteps) {
                IJ.showProgress(clusterCentres.size(), totalSteps);
                // The mask returns the coordinates with the centre of the image at 0,0
                centre[0] += width / 2;
                centre[1] += width / 2;
                clusterCentres.add(centre);
                // Fill in the mask around the centre to exclude any more circles that could overlap
                final double cx = centre[0] / maskScale;
                final double cy = centre[1] / maskScale;
                fillMask(mask, maskSize, (int) cx, (int) cy, roiRadius, 0);
                try {
                    maskDistribution = new MaskDistribution(mask, maskSize, maskSize, 0, maskScale, maskScale, rng);
                } catch (final IllegalArgumentException ex) {
                    // This can happen when there are no more non-zero pixels
                    log("WARNING: No more room for clusters on the mask area (created %d of estimated %d)", clusterCentres.size(), totalSteps);
                    break;
                }
            }
            ImageJUtils.finished();
        } else {
            // Pick centres randomly from the distribution
            while (clusterCentres.size() < totalSteps) {
                clusterCentres.add(dist.next());
            }
        }
        final double scaledRadius = settings.clusterRadius / maskScale;
        if (settings.showClusterMask || settings.clusterSimulation == 3) {
            // Show the mask for the clusters
            if (mask == null) {
                mask = new int[maskSize * maskSize];
            } else {
                Arrays.fill(mask, 0);
            }
            final int roiRadius = (int) Math.round(scaledRadius);
            for (final double[] c : clusterCentres) {
                final double cx = c[0] / maskScale;
                final double cy = c[1] / maskScale;
                fillMask(mask, maskSize, (int) cx, (int) cy, roiRadius, 1);
            }
            if (settings.clusterSimulation == 3) {
                // We have the mask. Now pick points at random from the mask.
                final MaskDistribution maskDistribution = new MaskDistribution(mask, maskSize, maskSize, 0, maskScale, maskScale, rng);
                // Allocate each molecule position to a parent circle so defining clusters.
                final int[][] clusters = new int[clusterCentres.size()][];
                final int[] clusterSize = new int[clusters.length];
                for (int i = 0; i < settings.numberOfMolecules; i++) {
                    final double[] centre = maskDistribution.next();
                    // The mask returns the coordinates with the centre of the image at 0,0
                    centre[0] += width / 2;
                    centre[1] += width / 2;
                    xyz.add(centre);
                    // Output statistics on cluster size and number.
                    // TODO - Finding the closest cluster could be done better than an all-vs-all comparison
                    double max = distance2(centre, clusterCentres.get(0));
                    int cluster = 0;
                    for (int j = 1; j < clusterCentres.size(); j++) {
                        final double d2 = distance2(centre, clusterCentres.get(j));
                        if (d2 < max) {
                            max = d2;
                            cluster = j;
                        }
                    }
                    // Assign point i to cluster
                    centre[2] = cluster;
                    if (clusterSize[cluster] == 0) {
                        clusters[cluster] = new int[10];
                    }
                    if (clusters[cluster].length <= clusterSize[cluster]) {
                        clusters[cluster] = Arrays.copyOf(clusters[cluster], (int) (clusters[cluster].length * 1.5));
                    }
                    clusters[cluster][clusterSize[cluster]++] = i;
                }
                // Generate real cluster size statistics
                for (int j = 0; j < clusterSize.length; j++) {
                    final int size = clusterSize[j];
                    if (size == 0) {
                        continue;
                    }
                    statsSize.add(size);
                    if (size == 1) {
                        statsRadius.add(0);
                        continue;
                    }
                    // Find centre of cluster and add the distance to each point
                    final double[] com = new double[2];
                    for (int n = 0; n < size; n++) {
                        final double[] xy = xyz.get(clusters[j][n]);
                        for (int k = 0; k < 2; k++) {
                            com[k] += xy[k];
                        }
                    }
                    for (int k = 0; k < 2; k++) {
                        com[k] /= size;
                    }
                    for (int n = 0; n < size; n++) {
                        final double dx = xyz.get(clusters[j][n])[0] - com[0];
                        final double dy = xyz.get(clusters[j][n])[1] - com[1];
                        statsRadius.add(Math.sqrt(dx * dx + dy * dy));
                    }
                }
            }
            if (settings.showClusterMask) {
                bp = new ByteProcessor(maskSize, maskSize);
                for (int i = 0; i < mask.length; i++) {
                    if (mask[i] != 0) {
                        bp.set(i, 128);
                    }
                }
                ImageJUtils.display(maskTitle, bp);
            }
        }
        // Use the simulated cluster centres to create clusters of the desired size
        if (settings.clusterSimulation == 1 || settings.clusterSimulation == 2) {
            for (final double[] clusterCentre : clusterCentres) {
                final int clusterN = (int) Math.round((settings.clusterNumberStdDev > 0) ? settings.clusterNumber + gauss.sample() * settings.clusterNumberStdDev : settings.clusterNumber);
                if (clusterN < 1) {
                    continue;
                }
                if (clusterN == 1) {
                    // No need for a cluster around a point
                    xyz.add(clusterCentre);
                    statsRadius.add(0);
                    statsSize.add(1);
                } else {
                    // Generate N random points within a circle of the chosen cluster radius.
                    // Locate the centre-of-mass and the average distance to the centre.
                    final double[] com = new double[3];
                    int size = 0;
                    while (size < clusterN) {
                        // Generate a random point within a circle uniformly
                        // http://stackoverflow.com/questions/5837572/generate-a-random-point-within-a-circle-uniformly
                        final double t = 2.0 * Math.PI * rng.nextDouble();
                        final double u = rng.nextDouble() + rng.nextDouble();
                        final double r = settings.clusterRadius * ((u > 1) ? 2 - u : u);
                        final double x = r * Math.cos(t);
                        final double y = r * Math.sin(t);
                        final double[] xy = new double[] { clusterCentre[0] + x, clusterCentre[1] + y };
                        xyz.add(xy);
                        for (int k = 0; k < 2; k++) {
                            com[k] += xy[k];
                        }
                        size++;
                    }
                    // Add the distance of the points from the centre of the cluster.
                    // Note this does not account for the movement due to precision.
                    statsSize.add(size);
                    if (size == 1) {
                        statsRadius.add(0);
                    } else {
                        for (int k = 0; k < 2; k++) {
                            com[k] /= size;
                        }
                        while (size > 0) {
                            final double dx = xyz.get(xyz.size() - size)[0] - com[0];
                            final double dy = xyz.get(xyz.size() - size)[1] - com[1];
                            statsRadius.add(Math.sqrt(dx * dx + dy * dy));
                            size--;
                        }
                    }
                }
            }
        }
    } else {
        // Random distribution
        for (int i = 0; i < settings.numberOfMolecules; i++) {
            xyz.add(dist.next());
        }
    }
    // The Gaussian sigma should be applied so the overall distance from the centre
    // ( sqrt(x^2+y^2) ) has a standard deviation of sigmaS?
    final double sigma1D = settings.sigmaS / Math.sqrt(2);
    // Show optional histograms
    StoredDataStatistics intraDistances = null;
    StoredData blinks = null;
    if (settings.showHistograms) {
        final int capacity = (int) (xyz.size() * settings.blinkingRate);
        intraDistances = new StoredDataStatistics(capacity);
        blinks = new StoredData(capacity);
    }
    final Statistics statsSigma = new Statistics();
    for (int i = 0; i < xyz.size(); i++) {
        int occurrences = getBlinks(rng, settings.blinkingRate);
        if (blinks != null) {
            blinks.add(occurrences);
        }
        final int size = settings.molecules.size();
        // Get coordinates in nm
        final double[] moleculeXyz = xyz.get(i);
        if (bp != null && occurrences > 0) {
            bp.putPixel((int) Math.round(moleculeXyz[0] / maskScale), (int) Math.round(moleculeXyz[1] / maskScale), 255);
        }
        while (occurrences-- > 0) {
            final double[] localisationXy = Arrays.copyOf(moleculeXyz, 2);
            // Add random precision
            if (sigma1D > 0) {
                final double dx = gauss.sample() * sigma1D;
                final double dy = gauss.sample() * sigma1D;
                localisationXy[0] += dx;
                localisationXy[1] += dy;
                if (!dist.isWithinXy(localisationXy)) {
                    continue;
                }
                // Calculate mean-squared displacement
                statsSigma.add(dx * dx + dy * dy);
            }
            final double x = localisationXy[0];
            final double y = localisationXy[1];
            settings.molecules.add(new Molecule(x, y, i, 1));
            // Store in pixels
            final float xx = (float) (x / nmPerPixel);
            final float yy = (float) (y / nmPerPixel);
            final float[] params = PeakResult.createParams(0, 0, xx, yy, 0);
            settings.results.add(i + 1, (int) xx, (int) yy, 0, 0, 0, 0, params, null);
        }
        if (settings.molecules.size() > size) {
            count++;
            if (intraDistances != null) {
                final int newCount = settings.molecules.size() - size;
                if (newCount == 1) {
                    // No intra-molecule distances
                    continue;
                }
                // Get the distance matrix between these molecules
                final double[][] matrix = new double[newCount][newCount];
                for (int ii = size, x = 0; ii < settings.molecules.size(); ii++, x++) {
                    for (int jj = size + 1, y = 1; jj < settings.molecules.size(); jj++, y++) {
                        final double d2 = settings.molecules.get(ii).distance2(settings.molecules.get(jj));
                        matrix[x][y] = matrix[y][x] = d2;
                    }
                }
                // Get the maximum distance for particle linkage clustering of this molecule
                double max = 0;
                for (int x = 0; x < newCount; x++) {
                    // Compare to all-other molecules and get the minimum distance
                    // needed to join at least one
                    double linkDistance = Double.POSITIVE_INFINITY;
                    for (int y = 0; y < newCount; y++) {
                        if (x == y) {
                            continue;
                        }
                        if (matrix[x][y] < linkDistance) {
                            linkDistance = matrix[x][y];
                        }
                    }
                    // Check if this is larger
                    if (max < linkDistance) {
                        max = linkDistance;
                    }
                }
                intraDistances.add(Math.sqrt(max));
            }
        }
    }
    settings.results.end();
    if (bp != null) {
        final ImagePlus imp = ImageJUtils.display(maskTitle, bp);
        final Calibration cal = imp.getCalibration();
        cal.setUnit("nm");
        cal.pixelWidth = cal.pixelHeight = maskScale;
    }
    log("Simulation results");
    log("  * Molecules = %d (%d activated)", xyz.size(), count);
    log("  * Blinking rate = %s", MathUtils.rounded((double) settings.molecules.size() / xyz.size(), 4));
    log("  * Precision (Mean-displacement) = %s nm", (statsSigma.getN() > 0) ? MathUtils.rounded(Math.sqrt(statsSigma.getMean()), 4) : "0");
    if (intraDistances != null) {
        if (intraDistances.getN() == 0) {
            log("  * Mean Intra-Molecule particle linkage distance = 0 nm");
            log("  * Fraction of inter-molecule particle linkage @ 0 nm = 0 %%");
        } else {
            plot(blinks, "Blinks/Molecule", true);
            final double[][] intraHist = plot(intraDistances, "Intra-molecule particle linkage distance", false);
            // Determine 95th and 99th percentile
            // Will not be null as we requested a non-integer histogram.
            int p99 = intraHist[0].length - 1;
            final double limit1 = 0.99 * intraHist[1][p99];
            final double limit2 = 0.95 * intraHist[1][p99];
            while (intraHist[1][p99] > limit1 && p99 > 0) {
                p99--;
            }
            int p95 = p99;
            while (intraHist[1][p95] > limit2 && p95 > 0) {
                p95--;
            }
            log("  * Mean Intra-Molecule particle linkage distance = %s nm" + " (95%% = %s, 99%% = %s, 100%% = %s)", MathUtils.rounded(intraDistances.getMean(), 4), MathUtils.rounded(intraHist[0][p95], 4), MathUtils.rounded(intraHist[0][p99], 4), MathUtils.rounded(intraHist[0][intraHist[0].length - 1], 4));
            if (settings.distanceAnalysis) {
                performDistanceAnalysis(intraHist, p99);
            }
        }
    }
    if (settings.clusterSimulation > 0) {
        log("  * Cluster number = %s +/- %s", MathUtils.rounded(statsSize.getMean(), 4), MathUtils.rounded(statsSize.getStandardDeviation(), 4));
        log("  * Cluster radius = %s +/- %s nm (mean distance to centre-of-mass)", MathUtils.rounded(statsRadius.getMean(), 4), MathUtils.rounded(statsRadius.getStandardDeviation(), 4));
    }
}

Example 58 with Gaussian

use of org.apache.commons.math3.analysis.function.Gaussian in project GDSC-SMLM by aherbert.

the class PcPalmMolecules method calculateAveragePrecision.

/**
 * Calculate the average precision by fitting a skewed Gaussian to the histogram of the precision
 * distribution.
 *
 * <p>A simple mean and SD of the histogram is computed. If the mean of the Skewed Gaussian does
 * not fit within 3 SDs of the simple mean then the simple mean is returned.
 *
 * @param molecules the molecules
 * @param title the plot title (null if no plot should be displayed)
 * @param histogramBins the histogram bins
 * @param logFitParameters Record the fit parameters to the ImageJ log
 * @param removeOutliers The distribution is created using all values within 1.5x the
 *        inter-quartile range (IQR) of the data
 * @return The average precision
 */
public double calculateAveragePrecision(List<Molecule> molecules, String title, int histogramBins, boolean logFitParameters, boolean removeOutliers) {
    // Plot histogram of the precision
    final float[] data = new float[molecules.size()];
    final DescriptiveStatistics stats = new DescriptiveStatistics();
    double yMin = Double.NEGATIVE_INFINITY;
    double yMax = 0;
    for (int i = 0; i < data.length; i++) {
        data[i] = (float) molecules.get(i).precision;
        stats.addValue(data[i]);
    }
    // Set the min and max y-values using 1.5 x IQR
    if (removeOutliers) {
        final double lower = stats.getPercentile(25);
        final double upper = stats.getPercentile(75);
        if (Double.isNaN(lower) || Double.isNaN(upper)) {
            if (logFitParameters) {
                ImageJUtils.log("Error computing IQR: %f - %f", lower, upper);
            }
        } else {
            final double iqr = upper - lower;
            yMin = Math.max(lower - iqr, stats.getMin());
            yMax = Math.min(upper + iqr, stats.getMax());
            if (logFitParameters) {
                ImageJUtils.log("  Data range: %f - %f. Plotting 1.5x IQR: %f - %f", stats.getMin(), stats.getMax(), yMin, yMax);
            }
        }
    }
    if (yMin == Double.NEGATIVE_INFINITY) {
        yMin = stats.getMin();
        yMax = stats.getMax();
        if (logFitParameters) {
            ImageJUtils.log("  Data range: %f - %f", yMin, yMax);
        }
    }
    int bins;
    if (histogramBins <= 0) {
        bins = (int) Math.ceil((stats.getMax() - stats.getMin()) / HistogramPlot.getBinWidthScottsRule(stats.getStandardDeviation(), (int) stats.getN()));
    } else {
        bins = histogramBins;
    }
    final float[][] hist = HistogramPlot.calcHistogram(data, yMin, yMax, bins);
    Plot plot = null;
    if (title != null) {
        plot = new Plot(title, "Precision", "Frequency");
        final float[] xValues = hist[0];
        final float[] yValues = hist[1];
        plot.addPoints(xValues, yValues, Plot.BAR);
        ImageJUtils.display(title, plot);
    }
    // Extract non-zero data
    float[] x = Arrays.copyOf(hist[0], hist[0].length);
    float[] y = Arrays.copyOf(hist[1], hist[1].length);
    int count = 0;
    for (int i = 0; i < y.length; i++) {
        if (y[i] > 0) {
            x[count] = x[i];
            y[count] = y[i];
            count++;
        }
    }
    x = Arrays.copyOf(x, count);
    y = Arrays.copyOf(y, count);
    // Sense check to fitted data. Get mean and SD of histogram
    final double[] stats2 = HistogramPlot.getHistogramStatistics(x, y);
    double mean = stats2[0];
    if (logFitParameters) {
        log("  Initial Statistics: %f +/- %f", stats2[0], stats2[1]);
    }
    // Standard Gaussian fit
    final double[] parameters = fitGaussian(x, y);
    if (parameters == null) {
        log("  Failed to fit initial Gaussian");
        return mean;
    }
    double newMean = parameters[1];
    double error = Math.abs(stats2[0] - newMean) / stats2[1];
    if (error > 3) {
        log("  Failed to fit Gaussian: %f standard deviations from histogram mean", error);
        return mean;
    }
    if (newMean < yMin || newMean > yMax) {
        log("  Failed to fit Gaussian: %f outside data range %f - %f", newMean, yMin, yMax);
        return mean;
    }
    mean = newMean;
    if (logFitParameters) {
        log("  Initial Gaussian: %f @ %f +/- %f", parameters[0], parameters[1], parameters[2]);
    }
    final double[] initialSolution = new double[] { parameters[0], parameters[1], parameters[2], -1 };
    // Fit to a skewed Gaussian (or appropriate function)
    final double[] skewParameters = fitSkewGaussian(x, y, initialSolution);
    if (skewParameters == null) {
        log("  Failed to fit Skewed Gaussian");
        return mean;
    }
    final SkewNormalFunction sn = new SkewNormalFunction(skewParameters);
    if (logFitParameters) {
        log("  Skewed Gaussian: %f @ %f +/- %f (a = %f) => %f +/- %f", skewParameters[0], skewParameters[1], skewParameters[2], skewParameters[3], sn.getMean(), Math.sqrt(sn.getVariance()));
    }
    newMean = sn.getMean();
    error = Math.abs(stats2[0] - newMean) / stats2[1];
    if (error > 3) {
        log("  Failed to fit Skewed Gaussian: %f standard deviations from histogram mean", error);
        return mean;
    }
    if (newMean < yMin || newMean > yMax) {
        log("  Failed to fit Skewed Gaussian: %f outside data range %f - %f", newMean, yMin, yMax);
        return mean;
    }
    // Use original histogram x-axis to maintain all the bins
    if (plot != null) {
        plot.setColor(Color.red);
        addToPlot(plot, hist[0], skewParameters, Plot.LINE);
        plot.setColor(Color.black);
        ImageJUtils.display(title, plot);
    }
    // Return the average precision from the fitted curve
    return newMean;
}

Example 59 with Gaussian

use of org.apache.commons.math3.analysis.function.Gaussian in project GDSC-SMLM by aherbert.

the class PsfDrift method showHwhm.

private void showHwhm() {
    // Build a list of suitable images
    final List<String> titles = createImageList(false);
    if (titles.isEmpty()) {
        IJ.error(TITLE, "No suitable PSF images");
        return;
    }
    final GenericDialog gd = new GenericDialog(TITLE);
    gd.addMessage("Approximate the volume of the PSF as a Gaussian and\n" + "compute the equivalent Gaussian width.");
    settings = Settings.load();
    gd.addChoice("PSF", titles.toArray(new String[0]), settings.title);
    gd.addCheckbox("Use_offset", settings.useOffset);
    gd.addSlider("Smoothing", 0, 0.5, settings.smoothing);
    gd.addHelp(HelpUrls.getUrl("psf-hwhm"));
    gd.showDialog();
    if (gd.wasCanceled()) {
        return;
    }
    settings.title = gd.getNextChoice();
    settings.useOffset = gd.getNextBoolean();
    settings.smoothing = gd.getNextNumber();
    settings.save();
    imp = WindowManager.getImage(settings.title);
    if (imp == null) {
        IJ.error(TITLE, "No PSF image for image: " + settings.title);
        return;
    }
    psfSettings = getPsfSettings(imp);
    if (psfSettings == null) {
        IJ.error(TITLE, "No PSF settings for image: " + settings.title);
        return;
    }
    final int size = imp.getStackSize();
    final ImagePsfModel psf = createImagePsf(1, size, 1);
    final double[] w0 = psf.getAllHwhm0();
    final double[] w1 = psf.getAllHwhm1();
    // Get current centre
    final int centre = psfSettings.getCentreImage();
    // Extract valid values (some can be NaN)
    double[] sw0 = new double[w0.length];
    double[] sw1 = new double[w1.length];
    final TDoubleArrayList s0 = new TDoubleArrayList(w0.length);
    final TDoubleArrayList s1 = new TDoubleArrayList(w0.length);
    int c0 = 0;
    int c1 = 0;
    for (int i = 0; i < w0.length; i++) {
        if (Double.isFinite(w0[i])) {
            s0.add(i + 1);
            sw0[c0++] = w0[i];
        }
        if (Double.isFinite(w1[i])) {
            s1.add(i + 1);
            sw1[c1++] = w1[i];
        }
    }
    if (c0 == 0 && c1 == 0) {
        IJ.error(TITLE, "No computed HWHM for image: " + settings.title);
        return;
    }
    double[] slice0 = s0.toArray();
    sw0 = Arrays.copyOf(sw0, c0);
    double[] slice1 = s1.toArray();
    sw1 = Arrays.copyOf(sw1, c1);
    // Smooth
    if (settings.smoothing > 0) {
        final LoessInterpolator loess = new LoessInterpolator(settings.smoothing, 1);
        sw0 = loess.smooth(slice0, sw0);
        sw1 = loess.smooth(slice1, sw1);
    }
    final TDoubleArrayList minWx = new TDoubleArrayList();
    final TDoubleArrayList minWy = new TDoubleArrayList();
    for (int i = 0; i < w0.length; i++) {
        double weight = 0;
        if (Double.isFinite(w0[i])) {
            if (Double.isFinite(w1[i])) {
                weight = w0[i] * w1[i];
            } else {
                weight = w0[i] * w0[i];
            }
        } else if (Double.isFinite(w1[i])) {
            weight = w1[i] * w1[i];
        }
        if (weight != 0) {
            minWx.add(i + 1);
            minWy.add(Math.sqrt(weight));
        }
    }
    // Smooth the combined line
    final double[] cx = minWx.toArray();
    double[] cy = minWy.toArray();
    if (settings.smoothing > 0) {
        final LoessInterpolator loess = new LoessInterpolator(settings.smoothing, 1);
        cy = loess.smooth(cx, cy);
    }
    final int newCentre = SimpleArrayUtils.findMinIndex(cy);
    // Convert to FWHM
    final double fwhm = psfSettings.getFwhm();
    // Widths are in pixels
    final String title = TITLE + " HWHM";
    final Plot plot = new Plot(title, "Slice", "HWHM (px)");
    double[] limits = MathUtils.limits(sw0);
    limits = MathUtils.limits(limits, sw1);
    final double maxY = limits[1] * 1.05;
    plot.setLimits(1, size, 0, maxY);
    plot.setColor(Color.red);
    plot.addPoints(slice0, sw0, Plot.LINE);
    plot.setColor(Color.blue);
    plot.addPoints(slice1, sw1, Plot.LINE);
    plot.setColor(Color.magenta);
    plot.addPoints(cx, cy, Plot.LINE);
    plot.setColor(Color.black);
    plot.addLabel(0, 0, "X=red; Y=blue, Combined=Magenta");
    final PlotWindow pw = ImageJUtils.display(title, plot);
    // Show a non-blocking dialog to allow the centre to be updated ...
    // Add a label and dynamically update when the centre is moved.
    final NonBlockingExtendedGenericDialog gd2 = new NonBlockingExtendedGenericDialog(TITLE);
    final double scale = psfSettings.getPixelSize();
    // @formatter:off
    ImageJUtils.addMessage(gd2, "Update the PSF information?\n \n" + "Current z-centre = %d, FHWM = %s px (%s nm)\n", centre, MathUtils.rounded(fwhm), MathUtils.rounded(fwhm * scale));
    // @formatter:on
    gd2.addSlider("z-centre", cx[0], cx[cx.length - 1], newCentre);
    final TextField tf = gd2.getLastTextField();
    gd2.addMessage("");
    gd2.addAndGetButton("Reset", event -> tf.setText(Integer.toString(newCentre)));
    final Label label = gd2.getLastLabel();
    gd2.addCheckbox("Update_centre", settings.updateCentre);
    gd2.addCheckbox("Update_HWHM", settings.updateHwhm);
    gd2.enableYesNoCancel();
    gd2.hideCancelButton();
    final UpdateDialogListener dl = new UpdateDialogListener(cx, cy, maxY, newCentre, scale, pw, label);
    gd2.addDialogListener(dl);
    gd.addHelp(HelpUrls.getUrl("psf-hwhm"));
    gd2.showDialog();
    if (gd2.wasOKed() && (settings.updateCentre || settings.updateHwhm)) {
        final ImagePSF.Builder b = psfSettings.toBuilder();
        if (settings.updateCentre) {
            b.setCentreImage(dl.centre);
        }
        if (settings.updateHwhm) {
            b.setFwhm(dl.getFwhm());
        }
        imp.setProperty("Info", ImagePsfHelper.toString(b));
    }
}

Example 60 with Gaussian

use of org.apache.commons.math3.analysis.function.Gaussian in project GDSC-SMLM by aherbert.

the class PoissonGammaGaussianFisherInformation method findUpperLimit.

/**
 * Find the upper limit of the integrated function A^2/P. P is the Poisson-Gamma convolution, A is
 * the partial gradient.
 *
 * <p>When both A and P are convolved with a Gaussian kernel, the integral of this function - 1 is
 * the Fisher information.
 *
 * <p>This method is used to determine the upper limit of the function using a binary search.
 *
 * @param theta the Poisson mean
 * @param max the max of the function (returned from {@link #findMaximum(double, double)})
 * @param rel Relative threshold
 * @return [upper,upper value,evaluations]
 */
public double[] findUpperLimit(final double theta, double[] max, double rel) {
    if (rel < MIN_RELATIVE_TOLERANCE) {
        throw new IllegalArgumentException("Relative tolerance too small: " + rel);
    }
    final UnivariateFunction f = new UnivariateFunction() {

        double[] dgDp = new double[1];

        @Override
        public double value(double x) {
            final double G = PoissonGammaFunction.unscaledPoissonGammaPartial(x, theta, gain, dgDp);
            return getF(G, dgDp[0]);
        }
    };
    int eval = 0;
    // Increase from the max until the tolerance is achieved.
    // Use the mean to get a rough initial step size
    final double mean = theta * gain;
    double step = Math.max(mean, 1);
    double upper = max[0];
    double upperValue = max[1];
    final double threshold = upperValue * rel;
    double lower = upper;
    while (upperValue > threshold) {
        lower = upper;
        upper += step;
        step *= 2;
        eval++;
        upperValue = f.value(upper);
    }
    // Binary search the bracket between lower and upper
    while (lower + 1 < upper) {
        final double mid = (lower + upper) * 0.5;
        eval++;
        final double midg = f.value(mid);
        if (midg > threshold) {
            lower = mid;
        } else {
            upper = mid;
            upperValue = midg;
        }
    }
    return new double[] { upper, upperValue, eval };
}