Examples with SkewNormalFunction uk.ac.sussex.gdsc.smlm.function.SkewNormalFunction used on opensource projects

Example 1 with SkewNormalFunction

use of uk.ac.sussex.gdsc.smlm.function.SkewNormalFunction 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 2 with SkewNormalFunction

use of uk.ac.sussex.gdsc.smlm.function.SkewNormalFunction 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;
}