Examples with MemoryPeakResults gdsc.smlm.results.MemoryPeakResults used on opensource projects

Example 36 with MemoryPeakResults

use of gdsc.smlm.results.MemoryPeakResults in project GDSC-SMLM by aherbert.

the class BenchmarkSpotFilter method addSpotsToMemory.

/**
	 * Add all the true-positives to memory as a new results set
	 * 
	 * @param filterResults
	 */
void addSpotsToMemory(TIntObjectHashMap<FilterResult> filterResults) {
    final MemoryPeakResults results = new MemoryPeakResults();
    results.setName(TITLE + " TP " + id++);
    filterResults.forEachEntry(new TIntObjectProcedure<FilterResult>() {

        public boolean execute(int peak, FilterResult filterResult) {
            for (ScoredSpot spot : filterResult.spots) {
                if (spot.match) {
                    final float[] params = new float[] { 0, spot.getIntensity(), 0, spot.spot.x, spot.spot.y, 0, 0 };
                    results.addf(peak, spot.spot.x, spot.spot.y, spot.getIntensity(), 0d, 0f, params, null);
                }
            }
            return true;
        }
    });
    MemoryPeakResults.addResults(results);
}

Example 37 with MemoryPeakResults

use of gdsc.smlm.results.MemoryPeakResults in project GDSC-SMLM by aherbert.

the class BenchmarkFilterAnalysis method createResults.

/**
	 * Create peak results.
	 *
	 * @param filterResults
	 *            The results from running the filter (or null)
	 * @param filter
	 *            the filter
	 */
private MemoryPeakResults createResults(PreprocessedPeakResult[] filterResults, DirectFilter filter, boolean withBorder) {
    if (filterResults == null) {
        final MultiPathFilter multiPathFilter = createMPF(filter, minimalFilter);
        //multiPathFilter.setDebugFile("/tmp/filter.txt");
        filterResults = filterResults(multiPathFilter);
    }
    MemoryPeakResults results = new MemoryPeakResults();
    results.copySettings(this.results);
    results.setName(TITLE);
    if (withBorder) {
        // To produce the same results as the PeakFit plugin we must implement the border
        // functionality used in the FitWorker. This respects the border of the spot filter.
        FitEngineConfiguration config = new FitEngineConfiguration(new FitConfiguration());
        updateAllConfiguration(config);
        MaximaSpotFilter spotFilter = config.createSpotFilter(true);
        final int border = spotFilter.getBorder();
        int[] bounds = getBounds();
        final int borderLimitX = bounds[0] - border;
        final int borderLimitY = bounds[1] - border;
        for (PreprocessedPeakResult spot : filterResults) {
            if (spot.getX() > border && spot.getX() < borderLimitX && spot.getY() > border && spot.getY() < borderLimitY) {
                double[] p = spot.toGaussian2DParameters();
                float[] params = new float[p.length];
                for (int j = 0; j < p.length; j++) params[j] = (float) p[j];
                int frame = spot.getFrame();
                int origX = (int) p[Gaussian2DFunction.X_POSITION];
                int origY = (int) p[Gaussian2DFunction.Y_POSITION];
                results.addf(frame, origX, origY, 0, 0, spot.getNoise(), params, null);
            }
        }
    } else {
        for (PreprocessedPeakResult spot : filterResults) {
            double[] p = spot.toGaussian2DParameters();
            float[] params = new float[p.length];
            for (int j = 0; j < p.length; j++) params[j] = (float) p[j];
            int frame = spot.getFrame();
            int origX = (int) p[Gaussian2DFunction.X_POSITION];
            int origY = (int) p[Gaussian2DFunction.Y_POSITION];
            results.addf(frame, origX, origY, 0, 0, spot.getNoise(), params, null);
        }
    }
    return results;
}

Example 38 with MemoryPeakResults

use of gdsc.smlm.results.MemoryPeakResults 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 (blinkingDistribution == 3) {
        log("  - Clusters = %d", nMolecules);
        log("  - Simulation size = %s um", Utils.rounded(simulationSize, 4));
        log("  - Molecules/cluster = %s", Utils.rounded(blinkingRate, 4));
        log("  - Blinking distribution = %s", BLINKING_DISTRIBUTION[blinkingDistribution]);
        log("  - p-Value = %s", Utils.rounded(p, 4));
    } else {
        log("  - Molecules = %d", nMolecules);
        log("  - Simulation size = %s um", Utils.rounded(simulationSize, 4));
        log("  - Blinking rate = %s", Utils.rounded(blinkingRate, 4));
        log("  - Blinking distribution = %s", BLINKING_DISTRIBUTION[blinkingDistribution]);
    }
    log("  - Average precision = %s nm", Utils.rounded(sigmaS, 4));
    log("  - Clusters simulation = " + CLUSTER_SIMULATION[clusterSimulation]);
    if (clusterSimulation > 0) {
        log("  - Cluster number = %s +/- %s", Utils.rounded(clusterNumber, 4), Utils.rounded(clusterNumberSD, 4));
        log("  - Cluster radius = %s nm", Utils.rounded(clusterRadius, 4));
    }
    final double nmPerPixel = 100;
    double width = simulationSize * 1000.0;
    // Allow a border of 3 x sigma for +/- precision
    //if (blinkingRate > 1)
    width -= 3 * sigmaS;
    RandomGenerator randomGenerator = new Well19937c(System.currentTimeMillis() + System.identityHashCode(this));
    RandomDataGenerator dataGenerator = new RandomDataGenerator(randomGenerator);
    UniformDistribution dist = new UniformDistribution(null, new double[] { width, width, 0 }, randomGenerator.nextInt());
    molecules = new ArrayList<Molecule>(nMolecules);
    // Create some dummy results since the calibration is required for later analysis
    results = new MemoryPeakResults();
    results.setCalibration(new gdsc.smlm.results.Calibration(nmPerPixel, 1, 100));
    results.setSource(new NullSource("Molecule Simulation"));
    results.begin();
    int count = 0;
    // Generate a sequence of coordinates
    ArrayList<double[]> xyz = new ArrayList<double[]>((int) (nMolecules * 1.1));
    Statistics statsRadius = new Statistics();
    Statistics statsSize = new Statistics();
    String maskTitle = TITLE + " Cluster Mask";
    ByteProcessor bp = null;
    double maskScale = 0;
    if (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 = lowResolutionImageSize;
        int[] mask = null;
        // scale is in nm/pixel
        maskScale = width / maskSize;
        ArrayList<double[]> clusterCentres = new ArrayList<double[]>();
        int totalSteps = 1 + (int) Math.ceil(nMolecules / clusterNumber);
        if (clusterSimulation == 2 || 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, randomGenerator);
            double[] centre;
            IJ.showStatus("Computing clusters mask");
            int roiRadius = (int) Math.round((clusterRadius * 2) / maskScale);
            if (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 * Math.pow(clusterRadius / maskScale, 2)));
            }
            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
                double cx = centre[0] / maskScale;
                double cy = centre[1] / maskScale;
                fillMask(mask, maskSize, (int) cx, (int) cy, roiRadius, 0);
                //Utils.display("Mask", new ColorProcessor(maskSize, maskSize, mask));
                try {
                    maskDistribution = new MaskDistribution(mask, maskSize, maskSize, 0, maskScale, maskScale, randomGenerator);
                } catch (IllegalArgumentException e) {
                    // 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;
                }
            }
            IJ.showProgress(1);
            IJ.showStatus("");
        } else {
            // Pick centres randomly from the distribution 
            while (clusterCentres.size() < totalSteps) clusterCentres.add(dist.next());
        }
        if (showClusterMask || clusterSimulation == 3) {
            // Show the mask for the clusters
            if (mask == null)
                mask = new int[maskSize * maskSize];
            else
                Arrays.fill(mask, 0);
            int roiRadius = (int) Math.round((clusterRadius) / maskScale);
            for (double[] c : clusterCentres) {
                double cx = c[0] / maskScale;
                double cy = c[1] / maskScale;
                fillMask(mask, maskSize, (int) cx, (int) cy, roiRadius, 1);
            }
            if (clusterSimulation == 3) {
                // We have the mask. Now pick points at random from the mask.
                MaskDistribution maskDistribution = new MaskDistribution(mask, maskSize, maskSize, 0, maskScale, maskScale, randomGenerator);
                // Allocate each molecule position to a parent circle so defining clusters.
                int[][] clusters = new int[clusterCentres.size()][];
                int[] clusterSize = new int[clusters.length];
                for (int i = 0; i < nMolecules; i++) {
                    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++) {
                        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
                    double[] com = new double[2];
                    for (int n = 0; n < size; n++) {
                        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++) {
                        double dx = xyz.get(clusters[j][n])[0] - com[0];
                        double dy = xyz.get(clusters[j][n])[1] - com[1];
                        statsRadius.add(Math.sqrt(dx * dx + dy * dy));
                    }
                }
            }
            if (showClusterMask) {
                bp = new ByteProcessor(maskSize, maskSize);
                for (int i = 0; i < mask.length; i++) if (mask[i] != 0)
                    bp.set(i, 128);
                Utils.display(maskTitle, bp);
            }
        }
        // Use the simulated cluster centres to create clusters of the desired size
        if (clusterSimulation == 1 || clusterSimulation == 2) {
            for (double[] clusterCentre : clusterCentres) {
                int clusterN = (int) Math.round((clusterNumberSD > 0) ? dataGenerator.nextGaussian(clusterNumber, clusterNumberSD) : clusterNumber);
                if (clusterN < 1)
                    continue;
                //double[] clusterCentre = dist.next();
                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.
                    double[] com = new double[3];
                    int j = 0;
                    while (j < clusterN) {
                        // Generate a random point within a circle uniformly
                        // http://stackoverflow.com/questions/5837572/generate-a-random-point-within-a-circle-uniformly
                        double t = 2.0 * Math.PI * randomGenerator.nextDouble();
                        double u = randomGenerator.nextDouble() + randomGenerator.nextDouble();
                        double r = clusterRadius * ((u > 1) ? 2 - u : u);
                        double x = r * Math.cos(t);
                        double y = r * Math.sin(t);
                        double[] xy = new double[] { clusterCentre[0] + x, clusterCentre[1] + y };
                        xyz.add(xy);
                        for (int k = 0; k < 2; k++) com[k] += xy[k];
                        j++;
                    }
                    // 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(j);
                    if (j == 1) {
                        statsRadius.add(0);
                    } else {
                        for (int k = 0; k < 2; k++) com[k] /= j;
                        while (j > 0) {
                            double dx = xyz.get(xyz.size() - j)[0] - com[0];
                            double dy = xyz.get(xyz.size() - j)[1] - com[1];
                            statsRadius.add(Math.sqrt(dx * dx + dy * dy));
                            j--;
                        }
                    }
                }
            }
        }
    } else {
        // Random distribution
        for (int i = 0; i < nMolecules; 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 = sigmaS / Math.sqrt(2);
    // Show optional histograms
    StoredDataStatistics intraDistances = null;
    StoredData blinks = null;
    if (showHistograms) {
        int capacity = (int) (xyz.size() * blinkingRate);
        intraDistances = new StoredDataStatistics(capacity);
        blinks = new StoredData(capacity);
    }
    Statistics statsSigma = new Statistics();
    for (int i = 0; i < xyz.size(); i++) {
        int nOccurrences = getBlinks(dataGenerator, blinkingRate);
        if (showHistograms)
            blinks.add(nOccurrences);
        final int size = molecules.size();
        // Get coordinates in nm
        final double[] moleculeXyz = xyz.get(i);
        if (bp != null && nOccurrences > 0) {
            bp.putPixel((int) Math.round(moleculeXyz[0] / maskScale), (int) Math.round(moleculeXyz[1] / maskScale), 255);
        }
        while (nOccurrences-- > 0) {
            final double[] localisationXy = Arrays.copyOf(moleculeXyz, 2);
            // Add random precision
            if (sigma1D > 0) {
                final double dx = dataGenerator.nextGaussian(0, sigma1D);
                final double dy = dataGenerator.nextGaussian(0, 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];
            molecules.add(new Molecule(x, y, i, 1));
            // Store in pixels
            float[] params = new float[7];
            params[Gaussian2DFunction.X_POSITION] = (float) (x / nmPerPixel);
            params[Gaussian2DFunction.Y_POSITION] = (float) (y / nmPerPixel);
            results.addf(i + 1, (int) x, (int) y, 0, 0, 0, params, null);
        }
        if (molecules.size() > size) {
            count++;
            if (showHistograms) {
                int newCount = molecules.size() - size;
                if (newCount == 1) {
                    //intraDistances.add(0);
                    continue;
                }
                // Get the distance matrix between these molecules
                double[][] matrix = new double[newCount][newCount];
                for (int ii = size, x = 0; ii < molecules.size(); ii++, x++) {
                    for (int jj = size + 1, y = 1; jj < molecules.size(); jj++, y++) {
                        final double d2 = molecules.get(ii).distance2(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));
            }
        }
    }
    results.end();
    if (bp != null)
        Utils.display(maskTitle, bp);
    // Used for debugging
    //System.out.printf("  * Molecules = %d (%d activated)\n", xyz.size(), count);
    //if (clusterSimulation > 0)
    //	System.out.printf("  * Cluster number = %s +/- %s. Radius = %s +/- %s\n",
    //			Utils.rounded(statsSize.getMean(), 4), Utils.rounded(statsSize.getStandardDeviation(), 4),
    //			Utils.rounded(statsRadius.getMean(), 4), Utils.rounded(statsRadius.getStandardDeviation(), 4));
    log("Simulation results");
    log("  * Molecules = %d (%d activated)", xyz.size(), count);
    log("  * Blinking rate = %s", Utils.rounded((double) molecules.size() / xyz.size(), 4));
    log("  * Precision (Mean-displacement) = %s nm", (statsSigma.getN() > 0) ? Utils.rounded(Math.sqrt(statsSigma.getMean()), 4) : "0");
    if (showHistograms) {
        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);
            double[][] intraHist = plot(intraDistances, "Intra-molecule particle linkage distance", false);
            // Determine 95th and 99th percentile
            int p99 = intraHist[0].length - 1;
            double limit1 = 0.99 * intraHist[1][p99];
            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)", Utils.rounded(intraDistances.getMean(), 4), Utils.rounded(intraHist[0][p95], 4), Utils.rounded(intraHist[0][p99], 4), Utils.rounded(intraHist[0][intraHist[0].length - 1], 4));
            if (distanceAnalysis) {
                performDistanceAnalysis(intraHist, p99);
            }
        }
    }
    if (clusterSimulation > 0) {
        log("  * Cluster number = %s +/- %s", Utils.rounded(statsSize.getMean(), 4), Utils.rounded(statsSize.getStandardDeviation(), 4));
        log("  * Cluster radius = %s +/- %s nm (mean distance to centre-of-mass)", Utils.rounded(statsRadius.getMean(), 4), Utils.rounded(statsRadius.getStandardDeviation(), 4));
    }
}

Example 39 with MemoryPeakResults

use of gdsc.smlm.results.MemoryPeakResults in project GDSC-SMLM by aherbert.

the class PCPALMClusters method doClustering.

/**
	 * Extract the results from the PCPALM molecules using the area ROI and then do clustering to obtain the histogram
	 * of molecules per cluster.
	 * 
	 * @return
	 */
private HistogramData doClustering() {
    // Perform clustering analysis to generate the histogram of cluster sizes
    PCPALMAnalysis analysis = new PCPALMAnalysis();
    ArrayList<Molecule> molecules = analysis.cropToRoi(WindowManager.getCurrentImage());
    if (molecules.size() < 2) {
        error("No results within the crop region");
        return null;
    }
    Utils.log("Using %d molecules (Density = %s um^-2) @ %s nm", molecules.size(), Utils.rounded(molecules.size() / analysis.croppedArea), Utils.rounded(distance));
    long s1 = System.nanoTime();
    ClusteringEngine engine = new ClusteringEngine(1, clusteringAlgorithm, new IJTrackProgress());
    if (multiThread)
        engine.setThreadCount(Prefs.getThreads());
    engine.setTracker(new IJTrackProgress());
    IJ.showStatus("Clustering ...");
    ArrayList<Cluster> clusters = engine.findClusters(convertToPoint(molecules), distance);
    IJ.showStatus("");
    if (clusters == null) {
        Utils.log("Aborted");
        return null;
    }
    nMolecules = molecules.size();
    Utils.log("Finished : %d total clusters (%s ms)", clusters.size(), Utils.rounded((System.nanoTime() - s1) / 1e6));
    // Save cluster centroids to a results set in memory. Then they can be plotted.
    MemoryPeakResults results = new MemoryPeakResults(clusters.size());
    results.setName(TITLE);
    // Set an arbitrary calibration so that the lifetime of the results is stored in the exposure time
    // The results will be handled as a single mega-frame containing all localisation. 
    results.setCalibration(new Calibration(100, 1, PCPALMMolecules.seconds * 1000));
    // Make the standard deviation such that the Gaussian volume will be 95% at the distance threshold
    final float sd = (float) (distance / 1.959964);
    int id = 0;
    for (Cluster c : clusters) {
        results.add(new ExtendedPeakResult((float) c.x, (float) c.y, sd, c.n, ++id));
    }
    MemoryPeakResults.addResults(results);
    // Get the data for fitting
    float[] values = new float[clusters.size()];
    for (int i = 0; i < values.length; i++) values[i] = clusters.get(i).n;
    float yMax = (int) Math.ceil(Maths.max(values));
    int nBins = (int) (yMax + 1);
    float[][] hist = Utils.calcHistogram(values, 0, yMax, nBins);
    HistogramData histogramData = (calibrateHistogram) ? new HistogramData(hist, frames, area, units) : new HistogramData(hist);
    saveHistogram(histogramData);
    return histogramData;
}

Example 40 with MemoryPeakResults

use of gdsc.smlm.results.MemoryPeakResults in project GDSC-SMLM by aherbert.

the class TraceMolecules method fitTraces.

private void fitTraces(MemoryPeakResults results, Trace[] traces) {
    // Check if the original image is open and the fit configuration can be extracted
    ImageSource source = results.getSource();
    if (source == null)
        return;
    if (!source.open())
        return;
    FitEngineConfiguration config = (FitEngineConfiguration) XmlUtils.fromXML(results.getConfiguration());
    if (config == null)
        return;
    // Show a dialog asking if the traces should be refit
    ExtendedGenericDialog gd = new ExtendedGenericDialog(TITLE);
    gd.addMessage("Do you want to fit the traces as a single peak using a combined image?");
    gd.addCheckbox("Fit_closest_to_centroid", !fitOnlyCentroid);
    gd.addSlider("Distance_threshold", 0.01, 3, distanceThreshold);
    gd.addSlider("Expansion_factor", 1, 4.5, expansionFactor);
    // Allow fitting settings to be adjusted
    FitConfiguration fitConfig = config.getFitConfiguration();
    gd.addMessage("--- Gaussian fitting ---");
    String[] filterTypes = SettingsManager.getNames((Object[]) DataFilterType.values());
    gd.addChoice("Spot_filter_type", filterTypes, filterTypes[config.getDataFilterType().ordinal()]);
    String[] filterNames = SettingsManager.getNames((Object[]) DataFilter.values());
    gd.addChoice("Spot_filter", filterNames, filterNames[config.getDataFilter(0).ordinal()]);
    gd.addSlider("Smoothing", 0, 2.5, config.getSmooth(0));
    gd.addSlider("Search_width", 0.5, 2.5, config.getSearch());
    gd.addSlider("Border", 0.5, 2.5, config.getBorder());
    gd.addSlider("Fitting_width", 2, 4.5, config.getFitting());
    String[] solverNames = SettingsManager.getNames((Object[]) FitSolver.values());
    gd.addChoice("Fit_solver", solverNames, solverNames[fitConfig.getFitSolver().ordinal()]);
    String[] functionNames = SettingsManager.getNames((Object[]) FitFunction.values());
    gd.addChoice("Fit_function", functionNames, functionNames[fitConfig.getFitFunction().ordinal()]);
    String[] criteriaNames = SettingsManager.getNames((Object[]) FitCriteria.values());
    gd.addChoice("Fit_criteria", criteriaNames, criteriaNames[fitConfig.getFitCriteria().ordinal()]);
    gd.addNumericField("Significant_digits", fitConfig.getSignificantDigits(), 0);
    gd.addNumericField("Coord_delta", fitConfig.getDelta(), 4);
    gd.addNumericField("Lambda", fitConfig.getLambda(), 4);
    gd.addNumericField("Max_iterations", fitConfig.getMaxIterations(), 0);
    gd.addNumericField("Fail_limit", config.getFailuresLimit(), 0);
    gd.addCheckbox("Include_neighbours", config.isIncludeNeighbours());
    gd.addSlider("Neighbour_height", 0.01, 1, config.getNeighbourHeightThreshold());
    gd.addSlider("Residuals_threshold", 0.01, 1, config.getResidualsThreshold());
    //gd.addSlider("Duplicate_distance", 0, 1.5, fitConfig.getDuplicateDistance());
    gd.addMessage("--- Peak filtering ---\nDiscard fits that shift; are too low; or expand/contract");
    gd.addCheckbox("Smart_filter", fitConfig.isSmartFilter());
    gd.addCheckbox("Disable_simple_filter", fitConfig.isDisableSimpleFilter());
    gd.addSlider("Shift_factor", 0.01, 2, fitConfig.getCoordinateShiftFactor());
    gd.addNumericField("Signal_strength", fitConfig.getSignalStrength(), 2);
    gd.addNumericField("Min_photons", fitConfig.getMinPhotons(), 0);
    gd.addSlider("Min_width_factor", 0, 0.99, fitConfig.getMinWidthFactor());
    gd.addSlider("Width_factor", 1.01, 5, fitConfig.getWidthFactor());
    gd.addNumericField("Precision", fitConfig.getPrecisionThreshold(), 2);
    gd.addCheckbox("Debug_failures", debugFailures);
    gd.showDialog();
    if (!gd.wasOKed()) {
        source.close();
        return;
    }
    // Get parameters for the fit
    fitOnlyCentroid = !gd.getNextBoolean();
    distanceThreshold = (float) gd.getNextNumber();
    expansionFactor = (float) gd.getNextNumber();
    config.setDataFilterType(gd.getNextChoiceIndex());
    config.setDataFilter(gd.getNextChoiceIndex(), Math.abs(gd.getNextNumber()), 0);
    config.setSearch(gd.getNextNumber());
    config.setBorder(gd.getNextNumber());
    config.setFitting(gd.getNextNumber());
    fitConfig.setFitSolver(gd.getNextChoiceIndex());
    fitConfig.setFitFunction(gd.getNextChoiceIndex());
    fitConfig.setFitCriteria(gd.getNextChoiceIndex());
    fitConfig.setSignificantDigits((int) gd.getNextNumber());
    fitConfig.setDelta(gd.getNextNumber());
    fitConfig.setLambda(gd.getNextNumber());
    fitConfig.setMaxIterations((int) gd.getNextNumber());
    config.setFailuresLimit((int) gd.getNextNumber());
    config.setIncludeNeighbours(gd.getNextBoolean());
    config.setNeighbourHeightThreshold(gd.getNextNumber());
    config.setResidualsThreshold(gd.getNextNumber());
    fitConfig.setSmartFilter(gd.getNextBoolean());
    fitConfig.setDisableSimpleFilter(gd.getNextBoolean());
    fitConfig.setCoordinateShiftFactor(gd.getNextNumber());
    fitConfig.setSignalStrength(gd.getNextNumber());
    fitConfig.setMinPhotons(gd.getNextNumber());
    fitConfig.setMinWidthFactor(gd.getNextNumber());
    fitConfig.setWidthFactor(gd.getNextNumber());
    fitConfig.setPrecisionThreshold(gd.getNextNumber());
    // Check arguments
    try {
        Parameters.isAboveZero("Distance threshold", distanceThreshold);
        Parameters.isAbove("Expansion factor", expansionFactor, 1);
        Parameters.isAboveZero("Search_width", config.getSearch());
        Parameters.isAboveZero("Fitting_width", config.getFitting());
        Parameters.isAboveZero("Significant digits", fitConfig.getSignificantDigits());
        Parameters.isAboveZero("Delta", fitConfig.getDelta());
        Parameters.isAboveZero("Lambda", fitConfig.getLambda());
        Parameters.isAboveZero("Max iterations", fitConfig.getMaxIterations());
        Parameters.isPositive("Failures limit", config.getFailuresLimit());
        Parameters.isPositive("Neighbour height threshold", config.getNeighbourHeightThreshold());
        Parameters.isPositive("Residuals threshold", config.getResidualsThreshold());
        Parameters.isPositive("Coordinate Shift factor", fitConfig.getCoordinateShiftFactor());
        Parameters.isPositive("Signal strength", fitConfig.getSignalStrength());
        Parameters.isPositive("Min photons", fitConfig.getMinPhotons());
        Parameters.isPositive("Min width factor", fitConfig.getMinWidthFactor());
        Parameters.isPositive("Width factor", fitConfig.getWidthFactor());
        Parameters.isPositive("Precision threshold", fitConfig.getPrecisionThreshold());
    } catch (IllegalArgumentException e) {
        IJ.error(TITLE, e.getMessage());
        source.close();
        return;
    }
    debugFailures = gd.getNextBoolean();
    if (!PeakFit.configureSmartFilter(globalSettings, filename))
        return;
    if (!PeakFit.configureDataFilter(globalSettings, filename, false))
        return;
    if (!PeakFit.configureFitSolver(globalSettings, filename, false))
        return;
    // Adjust settings for a single maxima
    config.setIncludeNeighbours(false);
    fitConfig.setDuplicateDistance(0);
    // Create a fit engine
    MemoryPeakResults refitResults = new MemoryPeakResults();
    refitResults.copySettings(results);
    refitResults.setName(results.getName() + " Trace Fit");
    refitResults.setSortAfterEnd(true);
    refitResults.begin();
    // No border since we know where the peaks are and we must not miss them due to truncated searching 
    FitEngine engine = new FitEngine(config, refitResults, Prefs.getThreads(), FitQueue.BLOCKING);
    // Either : Only fit the centroid
    // or     : Extract a bigger region, allowing all fits to run as normal and then 
    //          find the correct spot using Euclidian distance.
    // Set up the limits
    final double stdDev = FastMath.max(fitConfig.getInitialPeakStdDev0(), fitConfig.getInitialPeakStdDev1());
    float fitWidth = (float) (stdDev * config.getFitting() * ((fitOnlyCentroid) ? 1 : expansionFactor));
    IJ.showStatus("Refitting traces ...");
    List<JobItem> jobItems = new ArrayList<JobItem>(traces.length);
    int singles = 0;
    int fitted = 0;
    for (int n = 0; n < traces.length; n++) {
        Trace trace = traces[n];
        if (n % 32 == 0)
            IJ.showProgress(n, traces.length);
        // Skip traces with one peak
        if (trace.size() == 1) {
            singles++;
            // Use the synchronized method to avoid thread clashes with the FitEngine
            refitResults.addSync(trace.getHead());
            continue;
        }
        Rectangle bounds = new Rectangle();
        double[] combinedNoise = new double[1];
        float[] data = buildCombinedImage(source, trace, fitWidth, bounds, combinedNoise, false);
        if (data == null)
            continue;
        // Fit the combined image
        FitParameters params = new FitParameters();
        params.noise = (float) combinedNoise[0];
        float[] centre = trace.getCentroid();
        if (fitOnlyCentroid) {
            int newX = (int) Math.round(centre[0]) - bounds.x;
            int newY = (int) Math.round(centre[1]) - bounds.y;
            params.maxIndices = new int[] { newY * bounds.width + newX };
        } else {
            params.filter = new ArrayList<float[]>();
            params.filter.add(new float[] { centre[0] - bounds.x, centre[1] - bounds.y });
            params.distanceThreshold = distanceThreshold;
        }
        // This is not needed since the bounds are passed using the FitJob
        //params.setOffset(new float[] { bounds.x, bounds.y });
        int startT = trace.getHead().getFrame();
        params.endT = trace.getTail().getFrame();
        ParameterisedFitJob job = new ParameterisedFitJob(n, params, startT, data, bounds);
        jobItems.add(new JobItem(job, trace, centre));
        engine.run(job);
        fitted++;
    }
    engine.end(false);
    IJ.showStatus("");
    IJ.showProgress(1);
    // Check the success ...
    FitStatus[] values = FitStatus.values();
    int[] statusCount = new int[values.length + 1];
    ArrayList<String> names = new ArrayList<String>(Arrays.asList(SettingsManager.getNames((Object[]) values)));
    names.add(String.format("No maxima within %.2f of centroid", distanceThreshold));
    int separated = 0;
    int success = 0;
    final int debugLimit = 3;
    for (JobItem jobItem : jobItems) {
        int id = jobItem.getId();
        ParameterisedFitJob job = jobItem.job;
        Trace trace = jobItem.trace;
        int[] indices = job.getIndices();
        FitResult fitResult = null;
        int status;
        if (indices.length < 1) {
            status = values.length;
        } else if (indices.length > 1) {
            // Choose the first OK result. This is all that matters for the success reporting
            for (int n = 0; n < indices.length; n++) {
                if (job.getFitResult(n).getStatus() == FitStatus.OK) {
                    fitResult = job.getFitResult(n);
                    break;
                }
            }
            // Otherwise use the closest failure. 
            if (fitResult == null) {
                final float[] centre = traces[id].getCentroid();
                double minD = Double.POSITIVE_INFINITY;
                for (int n = 0; n < indices.length; n++) {
                    // Since the fit has failed we use the initial parameters
                    final double[] params = job.getFitResult(n).getInitialParameters();
                    final double dx = params[Gaussian2DFunction.X_POSITION] - centre[0];
                    final double dy = params[Gaussian2DFunction.Y_POSITION] - centre[1];
                    final double d = dx * dx + dy * dy;
                    if (minD > d) {
                        minD = d;
                        fitResult = job.getFitResult(n);
                    }
                }
            }
            status = fitResult.getStatus().ordinal();
        } else {
            fitResult = job.getFitResult(0);
            status = fitResult.getStatus().ordinal();
        }
        // All jobs have only one peak
        statusCount[status]++;
        // Debug why any fits failed
        if (fitResult == null || fitResult.getStatus() != FitStatus.OK) {
            refitResults.addAll(trace.getPoints());
            separated += trace.size();
            if (debugFailures) {
                FitStatus s = (fitResult == null) ? FitStatus.UNKNOWN : fitResult.getStatus();
                // Only display the first n per category to limit the number of images
                double[] noise = new double[1];
                if (statusCount[status] <= debugLimit) {
                    Rectangle bounds = new Rectangle();
                    buildCombinedImage(source, trace, fitWidth, bounds, noise, true);
                    float[] centre = trace.getCentroid();
                    Utils.display(String.format("Trace %d (n=%d) : x=%f,y=%f", id, trace.size(), centre[0], centre[1]), slices);
                    switch(s) {
                        case INSUFFICIENT_PRECISION:
                            float precision = (Float) fitResult.getStatusData();
                            IJ.log(String.format("Trace %d (n=%d) : %s = %f", id, trace.size(), names.get(status), precision));
                            break;
                        case INSUFFICIENT_SIGNAL:
                            if (noise[0] == 0)
                                noise[0] = getCombinedNoise(trace);
                            float snr = (Float) fitResult.getStatusData();
                            IJ.log(String.format("Trace %d (n=%d) : %s = %f (noise=%.2f)", id, trace.size(), names.get(status), snr, noise[0]));
                            break;
                        case COORDINATES_MOVED:
                        case OUTSIDE_FIT_REGION:
                        case WIDTH_DIVERGED:
                            float[] shift = (float[]) fitResult.getStatusData();
                            IJ.log(String.format("Trace %d (n=%d) : %s = %.3f,%.3f", id, trace.size(), names.get(status), shift[0], shift[1]));
                            break;
                        default:
                            IJ.log(String.format("Trace %d (n=%d) : %s", id, trace.size(), names.get(status)));
                            break;
                    }
                }
            }
        } else {
            success++;
            if (debugFailures) {
                // Only display the first n per category to limit the number of images
                double[] noise = new double[1];
                if (statusCount[status] <= debugLimit) {
                    Rectangle bounds = new Rectangle();
                    buildCombinedImage(source, trace, fitWidth, bounds, noise, true);
                    float[] centre = trace.getCentroid();
                    Utils.display(String.format("Trace %d (n=%d) : x=%f,y=%f", id, trace.size(), centre[0], centre[1]), slices);
                }
            }
        }
    }
    IJ.log(String.format("Trace fitting : %d singles : %d / %d fitted : %d separated", singles, success, fitted, separated));
    if (separated > 0) {
        IJ.log("Reasons for fit failure :");
        // Start at i=1 to skip FitStatus.OK
        for (int i = 1; i < statusCount.length; i++) {
            if (statusCount[i] != 0)
                IJ.log("  " + names.get(i) + " = " + statusCount[i]);
        }
    }
    refitResults.end();
    MemoryPeakResults.addResults(refitResults);
    source.close();
}