Examples with PeakResultProcedure uk.ac.sussex.gdsc.smlm.results.procedures.PeakResultProcedure used on opensource projects

Example 6 with PeakResultProcedure

use of uk.ac.sussex.gdsc.smlm.results.procedures.PeakResultProcedure in project GDSC-SMLM by aherbert.

the class TraceExporter method export.

private void export(MemoryPeakResults results) {
    // Copy to allow manipulation
    results = results.copy();
    // Strip results with no trace Id
    results.removeIf(result -> result.getId() <= 0);
    // Sort by ID then time
    results.sort(IdFramePeakResultComparator.INSTANCE);
    // Split traces with big jumps and long tracks into smaller tracks
    results = splitTraces(results);
    results = splitLongTraces(results);
    // Count each ID and remove short traces
    int id = 0;
    int count = 0;
    int tracks = 0;
    int maxLength = 0;
    final TIntHashSet remove = new TIntHashSet();
    for (int i = 0, size = results.size(); i < size; i++) {
        final PeakResult result = results.get(i);
        if (result.getId() != id) {
            if (count < settings.minLength) {
                remove.add(id);
            } else {
                tracks++;
                maxLength = Math.max(maxLength, count);
            }
            count = 0;
            id = result.getId();
        }
        count += getLength(result);
    }
    // Final ID
    if (count < settings.minLength) {
        remove.add(id);
    } else {
        tracks++;
        maxLength = Math.max(maxLength, count);
    }
    if (!remove.isEmpty()) {
        results.removeIf(result -> remove.contains(result.getId()));
        results.sort(IdFramePeakResultComparator.INSTANCE);
    }
    if (settings.wobble > 0) {
        // Just leave any exceptions to trickle up and kill the plugin
        final TypeConverter<DistanceUnit> c = results.getDistanceConverter(DistanceUnit.NM);
        final double w = c.convertBack(settings.wobble);
        final UniformRandomProvider rng = UniformRandomProviders.create();
        final NormalizedGaussianSampler gauss = SamplerUtils.createNormalizedGaussianSampler(rng);
        final boolean is3D = results.is3D();
        results.forEach((PeakResultProcedure) peakResult -> {
            peakResult.setXPosition((float) (peakResult.getXPosition() + w * gauss.sample()));
            peakResult.setYPosition((float) (peakResult.getYPosition() + w * gauss.sample()));
            if (is3D) {
                peakResult.setZPosition((float) (peakResult.getZPosition() + w * gauss.sample()));
            }
        });
    }
    if (settings.format == 2) {
        exportVbSpt(results);
    } else if (settings.format == 1) {
        exportAnaDda(results);
    } else {
        exportSpotOn(results);
    }
    ImageJUtils.log("Exported %s: %s in %s", results.getName(), TextUtils.pleural(results.size(), "localisation"), TextUtils.pleural(tracks, "track"));
    if (settings.showTraceLengths) {
        // We store and index (count-1)
        final int[] h = new int[maxLength];
        id = 0;
        for (int i = 0, size = results.size(); i < size; i++) {
            final PeakResult result = results.get(i);
            if (result.getId() != id) {
                h[count - 1]++;
                count = 0;
                id = result.getId();
            }
            count += getLength(result);
        }
        h[count - 1]++;
        final String title = TITLE + ": " + results.getName();
        final Plot plot = new Plot(title, "Length", "Frequency");
        plot.addPoints(SimpleArrayUtils.newArray(h.length, 1, 1.0f), SimpleArrayUtils.toFloat(h), Plot.BAR);
        plot.setLimits(SimpleArrayUtils.findIndex(h, i -> i != 0), maxLength + 1, 0, Double.NaN);
        ImageJUtils.display(title, plot);
    }
}

Example 7 with PeakResultProcedure

use of uk.ac.sussex.gdsc.smlm.results.procedures.PeakResultProcedure in project GDSC-SMLM by aherbert.

the class PsfCreator method runUsingFitting.

private void runUsingFitting() {
    if (!showFittingDialog()) {
        return;
    }
    if (!loadConfiguration()) {
        return;
    }
    final BasePoint[] spots = getSpots(0, true);
    if (spots.length == 0) {
        IJ.error(TITLE, "No spots without neighbours within " + (boxRadius * 2) + "px");
        return;
    }
    final ImageStack stack = getImageStack();
    final int width = imp.getWidth();
    final int height = imp.getHeight();
    final int currentSlice = imp.getSlice();
    // Adjust settings for a single maxima
    config.setIncludeNeighbours(false);
    final ArrayList<double[]> centres = new ArrayList<>(spots.length);
    final int iterations = 1;
    final LoessInterpolator loess = new LoessInterpolator(settings.getSmoothing(), iterations);
    // TODO - The fitting routine may not produce many points. In this instance the LOESS
    // interpolator
    // fails to smooth the data very well. A higher bandwidth helps this but perhaps
    // try a different smoothing method.
    // For each spot
    ImageJUtils.log(TITLE + ": " + imp.getTitle());
    ImageJUtils.log("Finding spot locations...");
    ImageJUtils.log("  %d spot%s without neighbours within %dpx", spots.length, ((spots.length == 1) ? "" : "s"), (boxRadius * 2));
    final StoredDataStatistics averageSd = new StoredDataStatistics();
    final StoredDataStatistics averageA = new StoredDataStatistics();
    final Statistics averageRange = new Statistics();
    final MemoryPeakResults allResults = new MemoryPeakResults();
    allResults.setCalibration(fitConfig.getCalibration());
    allResults.setPsf(fitConfig.getPsf());
    allResults.setName(TITLE);
    allResults.setBounds(new Rectangle(0, 0, width, height));
    MemoryPeakResults.addResults(allResults);
    for (int n = 1; n <= spots.length; n++) {
        final BasePoint spot = spots[n - 1];
        final int x = (int) spot.getX();
        final int y = (int) spot.getY();
        final MemoryPeakResults results = fitSpot(stack, width, height, x, y);
        allResults.add(results);
        if (results.size() < 5) {
            ImageJUtils.log("  Spot %d: Not enough fit results %d", n, results.size());
            continue;
        }
        // Get the results for the spot centre and width
        final double[] z = new double[results.size()];
        final double[] xCoord = new double[z.length];
        final double[] yCoord = new double[z.length];
        final double[] sd;
        final double[] a;
        final Counter counter = new Counter();
        // We have fit the results so they will be in the preferred units
        results.forEach(new PeakResultProcedure() {

            @Override
            public void execute(PeakResult peak) {
                final int i = counter.getAndIncrement();
                z[i] = peak.getFrame();
                xCoord[i] = peak.getXPosition() - x;
                yCoord[i] = peak.getYPosition() - y;
            }
        });
        final WidthResultProcedure wp = new WidthResultProcedure(results, DistanceUnit.PIXEL);
        wp.getW();
        sd = SimpleArrayUtils.toDouble(wp.wx);
        final HeightResultProcedure hp = new HeightResultProcedure(results, IntensityUnit.COUNT);
        hp.getH();
        a = SimpleArrayUtils.toDouble(hp.heights);
        // Smooth the amplitude plot
        final double[] smoothA = loess.smooth(z, a);
        // Find the maximum amplitude
        int maximumIndex = findMaximumIndex(smoothA);
        // Find the range at a fraction of the max. This is smoothed to find the X/Y centre
        int start = 0;
        int stop = smoothA.length - 1;
        final double limit = smoothA[maximumIndex] * settings.getAmplitudeFraction();
        for (int j = 0; j < smoothA.length; j++) {
            if (smoothA[j] > limit) {
                start = j;
                break;
            }
        }
        for (int j = smoothA.length; j-- > 0; ) {
            if (smoothA[j] > limit) {
                stop = j;
                break;
            }
        }
        averageRange.add(stop - start + 1);
        // 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];
        final double[] newZ = new double[smoothX.length];
        for (int j = start, k = 0; j <= stop; j++, k++) {
            smoothX[k] = xCoord[j];
            smoothY[k] = yCoord[j];
            smoothSd[k] = sd[j];
            newZ[k] = z[j];
        }
        smoothX = loess.smooth(newZ, smoothX);
        smoothY = loess.smooth(newZ, smoothY);
        smoothSd = loess.smooth(newZ, smoothSd);
        // Since the amplitude is not very consistent move from this peak to the
        // lowest width which is the in-focus spot.
        maximumIndex = findMinimumIndex(smoothSd, maximumIndex - start);
        // Find the centre at the amplitude peak
        final double cx = smoothX[maximumIndex] + x;
        final double cy = smoothY[maximumIndex] + y;
        int cz = (int) newZ[maximumIndex];
        double csd = smoothSd[maximumIndex];
        double ca = smoothA[maximumIndex + start];
        // The average should weight the SD using the signal for each spot
        averageSd.add(smoothSd[maximumIndex]);
        averageA.add(ca);
        if (ignoreSpot(n, z, a, smoothA, xCoord, yCoord, sd, newZ, smoothX, smoothY, smoothSd, cx, cy, cz, csd)) {
            ImageJUtils.log("  Spot %d was ignored", n);
            continue;
        }
        // Store result - it may have been moved interactively
        maximumIndex += this.slice - cz;
        cz = (int) newZ[maximumIndex];
        csd = smoothSd[maximumIndex];
        ca = smoothA[maximumIndex + start];
        ImageJUtils.log("  Spot %d => x=%.2f, y=%.2f, z=%d, sd=%.2f, A=%.2f", n, cx, cy, cz, csd, ca);
        centres.add(new double[] { cx, cy, cz, csd, n });
    }
    if (settings.getInteractiveMode()) {
        imp.setSlice(currentSlice);
        imp.setOverlay(null);
        // Hide the amplitude and spot plots
        ImageJUtils.hide(TITLE_AMPLITUDE);
        ImageJUtils.hide(TITLE_PSF_PARAMETERS);
    }
    if (centres.isEmpty()) {
        final String msg = "No suitable spots could be identified";
        ImageJUtils.log(msg);
        IJ.error(TITLE, msg);
        return;
    }
    // Find the limits of the z-centre
    int minz = (int) centres.get(0)[2];
    int maxz = minz;
    for (final double[] centre : centres) {
        if (minz > centre[2]) {
            minz = (int) centre[2];
        } else if (maxz < centre[2]) {
            maxz = (int) centre[2];
        }
    }
    IJ.showStatus("Creating PSF image");
    // Create a stack that can hold all the data.
    final ImageStack psf = createStack(stack, minz, maxz, settings.getMagnification());
    // For each spot
    final Statistics stats = new Statistics();
    boolean ok = true;
    for (int i = 0; ok && i < centres.size(); i++) {
        final double increment = 1.0 / (stack.getSize() * centres.size());
        setProgress((double) i / centres.size());
        final double[] centre = centres.get(i);
        // Extract the spot
        final float[][] spot = new float[stack.getSize()][];
        Rectangle regionBounds = null;
        for (int slice = 1; slice <= stack.getSize(); slice++) {
            final ImageExtractor ie = ImageExtractor.wrap((float[]) stack.getPixels(slice), width, height);
            if (regionBounds == null) {
                regionBounds = ie.getBoxRegionBounds((int) centre[0], (int) centre[1], boxRadius);
            }
            spot[slice - 1] = ie.crop(regionBounds);
        }
        if (regionBounds == null) {
            // Empty stack
            continue;
        }
        final int n = (int) centre[4];
        final float b = getBackground(n, spot);
        if (!subtractBackgroundAndWindow(spot, b, regionBounds.width, regionBounds.height, centre, loess)) {
            ImageJUtils.log("  Spot %d was ignored", n);
            continue;
        }
        stats.add(b);
        // Adjust the centre using the crop
        centre[0] -= regionBounds.x;
        centre[1] -= regionBounds.y;
        // This takes a long time so this should track progress
        ok = addToPsf(maxz, settings.getMagnification(), psf, centre, spot, regionBounds, increment, settings.getCentreEachSlice());
    }
    if (settings.getInteractiveMode()) {
        ImageJUtils.hide(TITLE_INTENSITY);
    }
    IJ.showProgress(1);
    if (!ok || stats.getN() == 0) {
        return;
    }
    final double avSd = getAverage(averageSd, averageA, 2);
    ImageJUtils.log("  Average background = %.2f, Av. SD = %s px", stats.getMean(), MathUtils.rounded(avSd, 4));
    normalise(psf, maxz, avSd * settings.getMagnification(), false);
    IJ.showProgress(1);
    psfImp = ImageJUtils.display(TITLE_PSF, psf);
    psfImp.setSlice(maxz);
    psfImp.resetDisplayRange();
    psfImp.updateAndDraw();
    final double[][] fitCom = new double[2][psf.getSize()];
    Arrays.fill(fitCom[0], Double.NaN);
    Arrays.fill(fitCom[1], Double.NaN);
    final double fittedSd = fitPsf(psf, loess, maxz, averageRange.getMean(), fitCom);
    // Compute the drift in the PSF:
    // - Use fitted centre if available; otherwise find CoM for each frame
    // - express relative to the average centre
    final double[][] com = calculateCentreOfMass(psf, fitCom, nmPerPixel / settings.getMagnification());
    final double[] slice = SimpleArrayUtils.newArray(psf.getSize(), 1, 1.0);
    final String title = TITLE + " CoM Drift";
    final Plot plot = new Plot(title, "Slice", "Drift (nm)");
    plot.addLabel(0, 0, "Red = X; Blue = Y");
    // double[] limitsX = Maths.limits(com[0]);
    // double[] limitsY = Maths.limits(com[1]);
    final double[] limitsX = getLimits(com[0]);
    final double[] limitsY = getLimits(com[1]);
    plot.setLimits(1, psf.getSize(), Math.min(limitsX[0], limitsY[0]), Math.max(limitsX[1], limitsY[1]));
    plot.setColor(Color.red);
    plot.addPoints(slice, com[0], Plot.DOT);
    plot.addPoints(slice, loess.smooth(slice, com[0]), Plot.LINE);
    plot.setColor(Color.blue);
    plot.addPoints(slice, com[1], Plot.DOT);
    plot.addPoints(slice, loess.smooth(slice, com[1]), Plot.LINE);
    ImageJUtils.display(title, plot);
    // TODO - Redraw the PSF with drift correction applied.
    // This means that the final image should have no drift.
    // This is relevant when combining PSF images. It doesn't matter too much for simulations
    // unless the drift is large.
    // Add Image properties containing the PSF details
    final double fwhm = getFwhm(psf, maxz);
    psfImp.setProperty("Info", ImagePsfHelper.toString(ImagePsfHelper.create(maxz, nmPerPixel / settings.getMagnification(), settings.getNmPerSlice(), stats.getN(), fwhm, createNote())));
    ImageJUtils.log("%s : z-centre = %d, nm/Pixel = %s, nm/Slice = %s, %d images, " + "PSF SD = %s nm, FWHM = %s px\n", psfImp.getTitle(), maxz, MathUtils.rounded(nmPerPixel / settings.getMagnification(), 3), MathUtils.rounded(settings.getNmPerSlice(), 3), stats.getN(), MathUtils.rounded(fittedSd * nmPerPixel, 4), MathUtils.rounded(fwhm));
    createInteractivePlots(psf, maxz, nmPerPixel / settings.getMagnification(), fittedSd * nmPerPixel);
    IJ.showStatus("");
}

Example 8 with PeakResultProcedure

use of uk.ac.sussex.gdsc.smlm.results.procedures.PeakResultProcedure in project GDSC-SMLM by aherbert.

the class Fire method createImages.

/**
 * Creates the images to use for the FIRE calculation. This must be called after
 * {@link #initialise(MemoryPeakResults, MemoryPeakResults)}.
 *
 * @param fourierImageScale the fourier image scale (set to zero to auto compute)
 * @param imageSize the image size
 * @param useSignal Use the localisation signal to weight the intensity. The default uses a value
 *        of 1 per localisation.
 * @return the fire images
 */
public FireImages createImages(double fourierImageScale, int imageSize, boolean useSignal) {
    if (results == null) {
        return null;
    }
    final SignalProvider signalProvider = (useSignal && (results.hasIntensity())) ? new PeakSignalProvider() : new FixedSignalProvider();
    // Draw images using the existing IJ routines.
    final Rectangle bounds = new Rectangle((int) Math.ceil(dataBounds.getWidth()), (int) Math.ceil(dataBounds.getHeight()));
    final ResultsImageSettings.Builder builder = ResultsImageSettings.newBuilder().setImageType(ResultsImageType.DRAW_NONE).setWeighted(true).setEqualised(false).setImageMode(ResultsImageMode.IMAGE_ADD);
    if (fourierImageScale > 0) {
        builder.setImageSizeMode(ResultsImageSizeMode.SCALED);
        builder.setScale(fourierImageScale);
    } else {
        builder.setImageSizeMode(ResultsImageSizeMode.IMAGE_SIZE);
        builder.setImageSize(imageSize);
    }
    ImageJImagePeakResults image1 = createPeakResultsImage(bounds, builder, "IP1");
    ImageJImagePeakResults image2 = createPeakResultsImage(bounds, builder, "IP2");
    final float minx = (float) dataBounds.getX();
    final float miny = (float) dataBounds.getY();
    if (this.results2 != null) {
        // Two image comparison
        final ImageJImagePeakResults i1 = image1;
        results.forEach((PeakResultProcedure) result -> {
            final float x = result.getXPosition() - minx;
            final float y = result.getYPosition() - miny;
            i1.add(x, y, signalProvider.getSignal(result));
        });
        final ImageJImagePeakResults i2 = image2;
        results2.forEach((PeakResultProcedure) result -> {
            final float x = result.getXPosition() - minx;
            final float y = result.getYPosition() - miny;
            i2.add(x, y, signalProvider.getSignal(result));
        });
    } else {
        // Block sampling.
        // Ensure we have at least 2 even sized blocks.
        int blockSize = Math.min(results.size() / 2, Math.max(1, settings.blockSize));
        int nblocks = (int) Math.ceil((double) results.size() / blockSize);
        while (nblocks <= 1 && blockSize > 1) {
            blockSize /= 2;
            nblocks = (int) Math.ceil((double) results.size() / blockSize);
        }
        if (nblocks <= 1) {
            // This should not happen since the results should contain at least 2 localisations
            return null;
        }
        if (blockSize != settings.blockSize) {
            IJ.log(pluginTitle + " Warning: Changed block size to " + blockSize);
        }
        final Counter i = new Counter();
        final Counter block = new Counter();
        final int finalBlockSize = blockSize;
        final PeakResult[][] blocks = new PeakResult[nblocks][blockSize];
        results.forEach((PeakResultProcedure) result -> {
            if (i.getCount() == finalBlockSize) {
                block.increment();
                i.reset();
            }
            blocks[block.getCount()][i.getAndIncrement()] = result;
        });
        // Truncate last block
        blocks[block.getCount()] = Arrays.copyOf(blocks[block.getCount()], i.getCount());
        final int[] indices = SimpleArrayUtils.natural(block.getCount() + 1);
        if (settings.randomSplit) {
            MathArrays.shuffle(indices);
        }
        for (final int index : indices) {
            // Split alternating so just rotate
            final ImageJImagePeakResults image = image1;
            image1 = image2;
            image2 = image;
            for (final PeakResult p : blocks[index]) {
                final float x = p.getXPosition() - minx;
                final float y = p.getYPosition() - miny;
                image.add(x, y, signalProvider.getSignal(p));
            }
        }
    }
    image1.end();
    final ImageProcessor ip1 = image1.getImagePlus().getProcessor();
    image2.end();
    final ImageProcessor ip2 = image2.getImagePlus().getProcessor();
    if (settings.maxPerBin > 0 && signalProvider instanceof FixedSignalProvider) {
        // We can eliminate over-sampled pixels
        for (int i = ip1.getPixelCount(); i-- > 0; ) {
            if (ip1.getf(i) > settings.maxPerBin) {
                ip1.setf(i, settings.maxPerBin);
            }
            if (ip2.getf(i) > settings.maxPerBin) {
                ip2.setf(i, settings.maxPerBin);
            }
        }
    }
    return new FireImages(ip1, ip2, nmPerUnit / image1.getScale());
}

Example 9 with PeakResultProcedure

use of uk.ac.sussex.gdsc.smlm.results.procedures.PeakResultProcedure in project GDSC-SMLM by aherbert.

the class FilterMolecules method run.

@Override
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
    MemoryPeakResults results = ResultsManager.loadInputResults(settings.inputOption, false, null, null);
    if (MemoryPeakResults.isEmpty(results)) {
        IJ.error(TITLE, "No results could be loaded");
        return;
    }
    // Allow reordering when filtering
    results = results.copy();
    if (settings.removeSingles) {
        results.removeIf(p -> p.getId() <= 0);
        final TIntIntHashMap count = new TIntIntHashMap(results.size());
        results.forEach((PeakResultProcedure) r -> count.adjustOrPutValue(r.getId(), 1, 1));
        results.removeIf(p -> count.get(p.getId()) == 1);
        if (results.isEmpty()) {
            IJ.error(TITLE, "No results after filtering singles");
            return;
        }
    }
    switch(settings.filterMode) {
        case D:
            new FilterDiffusionCoefficient().run(results);
            break;
        default:
            IJ.error(TITLE, "Unknown filter mode: " + settings.filterMode);
    }
}

Example 10 with PeakResultProcedure

use of uk.ac.sussex.gdsc.smlm.results.procedures.PeakResultProcedure in project GDSC-SMLM by aherbert.

the class Fire method calculatePrecisionHistogram.

/**
 * Calculate a histogram of the precision. The precision can be either stored in the results or
 * calculated using the Mortensen formula. If the precision method for Q estimation is not fixed
 * then the histogram is fitted with a Gaussian to create an initial estimate.
 *
 * @return The precision histogram
 */
private PrecisionHistogram calculatePrecisionHistogram() {
    final boolean logFitParameters = false;
    final String title = results.getName() + " Precision Histogram";
    // Check if the results has the precision already or if it can be computed.
    final boolean canUseStored = canUseStoredPrecision(results);
    final boolean canCalculatePrecision = canCalculatePrecision(results);
    // Set the method to compute a histogram. Default to the user selected option.
    PrecisionMethod method = null;
    if ((canUseStored && precisionMethod == PrecisionMethod.STORED) || (canCalculatePrecision && precisionMethod == PrecisionMethod.CALCULATE)) {
        method = precisionMethod;
    }
    if (method == null) {
        // We only have two choices so if one is available then select it.
        if (canUseStored) {
            method = PrecisionMethod.STORED;
        } else if (canCalculatePrecision) {
            method = PrecisionMethod.CALCULATE;
        }
        // If the user selected a method not available then log a warning
        if (method != null && precisionMethod != PrecisionMethod.FIXED) {
            IJ.log(String.format("%s : Selected precision method '%s' not available, switching to '%s'", pluginTitle, precisionMethod, method.getName()));
        }
        if (method == null) {
            // This does not matter if the user has provide a fixed input.
            if (precisionMethod == PrecisionMethod.FIXED) {
                final PrecisionHistogram histogram = new PrecisionHistogram(title);
                histogram.mean = settings.mean;
                histogram.sigma = settings.sigma;
                return histogram;
            }
            // No precision
            return null;
        }
    }
    // We get here if we can compute precision.
    // Build the histogram
    StoredDataStatistics precision = new StoredDataStatistics(results.size());
    if (method == PrecisionMethod.STORED) {
        final StoredDataStatistics p = precision;
        results.forEach((PeakResultProcedure) result -> p.add(result.getPrecision()));
    } else {
        precision.add(pp.precisions);
    }
    double yMin = Double.NEGATIVE_INFINITY;
    double yMax = 0;
    // Set the min and max y-values using 1.5 x IQR
    final DescriptiveStatistics stats = precision.getStatistics();
    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) {
        final int n = 5;
        yMin = Math.max(stats.getMin(), stats.getMean() - n * stats.getStandardDeviation());
        yMax = Math.min(stats.getMax(), stats.getMean() + n * stats.getStandardDeviation());
        if (logFitParameters) {
            ImageJUtils.log("  Data range: %f - %f. Plotting mean +/- %dxSD: %f - %f", stats.getMin(), stats.getMax(), n, yMin, yMax);
        }
    }
    // Get the data within the range
    final double[] data = precision.getValues();
    precision = new StoredDataStatistics(data.length);
    for (final double d : data) {
        if (d < yMin || d > yMax) {
            continue;
        }
        precision.add(d);
    }
    final int histogramBins = HistogramPlot.getBins(precision, HistogramPlot.BinMethod.SCOTT);
    final float[][] hist = HistogramPlot.calcHistogram(precision.getFloatValues(), yMin, yMax, histogramBins);
    final PrecisionHistogram histogram = new PrecisionHistogram(hist, precision, title);
    if (precisionMethod == PrecisionMethod.FIXED) {
        histogram.mean = settings.mean;
        histogram.sigma = settings.sigma;
        return histogram;
    }
    // Fitting of the histogram to produce the initial estimate
    // 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);
    if (logFitParameters) {
        ImageJUtils.log("  Initial Statistics: %f +/- %f", stats2[0], stats2[1]);
    }
    histogram.mean = stats2[0];
    histogram.sigma = stats2[1];
    // Standard Gaussian fit
    final double[] parameters = fitGaussian(x, y);
    if (parameters == null) {
        ImageJUtils.log("  Failed to fit initial Gaussian");
        return histogram;
    }
    final double newMean = parameters[1];
    final double error = Math.abs(stats2[0] - newMean) / stats2[1];
    if (error > 3) {
        ImageJUtils.log("  Failed to fit Gaussian: %f standard deviations from histogram mean", error);
        return histogram;
    }
    if (newMean < yMin || newMean > yMax) {
        ImageJUtils.log("  Failed to fit Gaussian: %f outside data range %f - %f", newMean, yMin, yMax);
        return histogram;
    }
    if (logFitParameters) {
        ImageJUtils.log("  Initial Gaussian: %f @ %f +/- %f", parameters[0], parameters[1], parameters[2]);
    }
    histogram.mean = parameters[1];
    histogram.sigma = parameters[2];
    return histogram;
}