Search in sources :

Example 1 with Cluster

use of org.apache.commons.math3.ml.clustering.Cluster in project GDSC-SMLM by aherbert.

the class PCPALMMolecules method getRunMode.

private boolean getRunMode(boolean resultsAvailable) {
    ExtendedGenericDialog gd = new ExtendedGenericDialog(TITLE);
    gd.addHelp(About.HELP_URL);
    // Build a list of all images with a region ROI
    List<String> titles = new LinkedList<String>();
    if (WindowManager.getWindowCount() > 0) {
        for (int imageID : WindowManager.getIDList()) {
            ImagePlus imp = WindowManager.getImage(imageID);
            if (imp != null && imp.getRoi() != null && imp.getRoi().isArea())
                titles.add(imp.getTitle());
        }
    }
    if (!resultsAvailable) {
        runMode = 3;
        gd.addMessage("Simulate molecules for cluster analysis.\nComputes a binary image from localisation data");
        gd.addNumericField("Molecules", nMolecules, 0);
        gd.addNumericField("Simulation_size (um)", simulationSize, 2);
        gd.addNumericField("Blinking_rate", blinkingRate, 2);
        gd.addChoice("Blinking_distribution", BLINKING_DISTRIBUTION, BLINKING_DISTRIBUTION[blinkingDistribution]);
        gd.addNumericField("Average_precision (nm)", sigmaS, 2);
        gd.addCheckbox("Show_histograms", showHistograms);
        gd.addCheckbox("Distance_analysis", distanceAnalysis);
        gd.addChoice("Cluster_simulation", CLUSTER_SIMULATION, CLUSTER_SIMULATION[clusterSimulation]);
        gd.addNumericField("Cluster_number", clusterNumber, 2);
        gd.addNumericField("Cluster_variation (SD)", clusterNumberSD, 2);
        gd.addNumericField("Cluster_radius", clusterRadius, 2);
        gd.addCheckbox("Show_cluster_mask", showClusterMask);
        Recorder.recordOption("Run_mode", RUN_MODE[runMode]);
    } else {
        gd.addMessage("Prepare molecules for cluster analysis.\nComputes a binary image from raw localisation data");
        ResultsManager.addInput(gd, inputOption, InputSource.MEMORY);
        if (!titles.isEmpty())
            gd.addCheckbox((titles.size() == 1) ? "Use_ROI" : "Choose_ROI", chooseRoi);
        gd.addChoice("Run_mode", RUN_MODE, RUN_MODE[runMode]);
    }
    gd.addMessage("Select options for low resolution image:");
    gd.addSlider("Image_size (px)", 512, 2048, lowResolutionImageSize);
    gd.addSlider("ROI_size (um)", 1.5, 4, roiSizeInUm);
    gd.addMessage("Select options for high resolution image:");
    gd.addCheckbox("Show_high_res_image", showHighResolutionImage);
    gd.addSlider("nm_per_pixel_limit", 0, 20, nmPerPixelLimit);
    gd.addMessage("Optionally remove all analysis results from memory");
    gd.addCheckbox("Clear_results", clearResults);
    gd.showDialog();
    if (gd.wasCanceled())
        return false;
    if (!resultsAvailable) {
        nMolecules = (int) Math.abs(gd.getNextNumber());
        simulationSize = Math.abs(gd.getNextNumber());
        blinkingRate = Math.abs(gd.getNextNumber());
        blinkingDistribution = gd.getNextChoiceIndex();
        sigmaS = Math.abs(gd.getNextNumber());
        showHistograms = gd.getNextBoolean();
        distanceAnalysis = gd.getNextBoolean();
        clusterSimulation = gd.getNextChoiceIndex();
        clusterNumber = Math.abs(gd.getNextNumber());
        clusterNumberSD = Math.abs(gd.getNextNumber());
        clusterRadius = Math.abs(gd.getNextNumber());
        showClusterMask = gd.getNextBoolean();
    } else {
        inputOption = ResultsManager.getInputSource(gd);
        if (!titles.isEmpty())
            chooseRoi = gd.getNextBoolean();
        runMode = gd.getNextChoiceIndex();
    }
    lowResolutionImageSize = (int) gd.getNextNumber();
    roiSizeInUm = gd.getNextNumber();
    showHighResolutionImage = gd.getNextBoolean();
    nmPerPixelLimit = Math.abs(gd.getNextNumber());
    clearResults = gd.getNextBoolean();
    // Check arguments
    try {
        if (!resultsAvailable) {
            Parameters.isAboveZero("Molecules", nMolecules);
            Parameters.isAboveZero("Simulation size", simulationSize);
            Parameters.isEqualOrAbove("Blinking rate", blinkingRate, 1);
            Parameters.isEqualOrAbove("Cluster number", clusterNumber, 1);
        }
        Parameters.isAbove("Image scale", lowResolutionImageSize, 1);
        Parameters.isAboveZero("ROI size", roiSizeInUm);
    } catch (IllegalArgumentException ex) {
        IJ.error(TITLE, ex.getMessage());
        return false;
    }
    if (!titles.isEmpty() && chooseRoi && resultsAvailable) {
        if (titles.size() == 1) {
            roiImage = titles.get(0);
            Recorder.recordOption("Image", roiImage);
        } else {
            String[] items = titles.toArray(new String[titles.size()]);
            gd = new ExtendedGenericDialog(TITLE);
            gd.addMessage("Select the source image for the ROI");
            gd.addChoice("Image", items, roiImage);
            gd.showDialog();
            if (gd.wasCanceled())
                return false;
            roiImage = gd.getNextChoice();
        }
        ImagePlus imp = WindowManager.getImage(roiImage);
        roiBounds = imp.getRoi().getBounds();
        roiImageWidth = imp.getWidth();
        roiImageHeight = imp.getHeight();
    } else {
        roiBounds = null;
    }
    if (!resultsAvailable) {
        if (!getPValue())
            return false;
    }
    if (clearResults) {
        PCPALMAnalysis.results.clear();
        PCPALMFitting.previous_gr = null;
    }
    return true;
}
Also used : ExtendedGenericDialog(ij.gui.ExtendedGenericDialog) ImagePlus(ij.ImagePlus) LinkedList(java.util.LinkedList) WeightedObservedPoint(org.apache.commons.math3.fitting.WeightedObservedPoint) ClusterPoint(gdsc.core.clustering.ClusterPoint)

Example 2 with Cluster

use of org.apache.commons.math3.ml.clustering.Cluster in project GDSC-SMLM by aherbert.

the class BinomialFitter method fitBinomial.

/**
	 * Fit the binomial distribution (n,p) to the input data. Performs fitting assuming a fixed n value and attempts to
	 * optimise p. All n from minN to maxN are evaluated. If maxN is zero then all possible n from minN are evaluated
	 * until the fit is worse.
	 * 
	 * @param data
	 *            The input data (all value must be positive)
	 * @param minN
	 *            The minimum n to evaluate
	 * @param maxN
	 *            The maximum n to evaluate. Set to zero to evaluate all possible values.
	 * @param zeroTruncated
	 *            True if the model should ignore n=0 (zero-truncated binomial)
	 * @return The best fit (n, p)
	 * @throws IllegalArgumentException
	 *             If any of the input data values are negative
	 */
public double[] fitBinomial(int[] data, int minN, int maxN, boolean zeroTruncated) {
    double[] histogram = getHistogram(data, false);
    final double initialSS = Double.POSITIVE_INFINITY;
    double bestSS = initialSS;
    double[] parameters = null;
    int worse = 0;
    int N = (int) histogram.length - 1;
    if (minN < 1)
        minN = 1;
    if (maxN > 0) {
        if (N > maxN) {
            // Limit the number fitted to maximum
            N = maxN;
        } else if (N < maxN) {
            // Expand the histogram to the maximum
            histogram = Arrays.copyOf(histogram, maxN + 1);
            N = maxN;
        }
    }
    if (minN > N)
        minN = N;
    final double mean = getMean(histogram);
    String name = (zeroTruncated) ? "Zero-truncated Binomial distribution" : "Binomial distribution";
    log("Mean cluster size = %s", Utils.rounded(mean));
    log("Fitting cumulative " + name);
    // score several times in succession)
    for (int n = minN; n <= N; n++) {
        PointValuePair solution = fitBinomial(histogram, mean, n, zeroTruncated);
        if (solution == null)
            continue;
        double p = solution.getPointRef()[0];
        log("Fitted %s : N=%d, p=%s. SS=%g", name, n, Utils.rounded(p), solution.getValue());
        if (bestSS > solution.getValue()) {
            bestSS = solution.getValue();
            parameters = new double[] { n, p };
            worse = 0;
        } else if (bestSS != initialSS) {
            if (++worse >= 3)
                break;
        }
    }
    return parameters;
}
Also used : PointValuePair(org.apache.commons.math3.optim.PointValuePair)

Example 3 with Cluster

use of org.apache.commons.math3.ml.clustering.Cluster in project GDSC-SMLM by aherbert.

the class PCPALMMolecules method performDistanceAnalysis.

private void performDistanceAnalysis(double[][] intraHist, int p99) {
    // We want to know the fraction of distances between molecules at the 99th percentile
    // that are intra- rather than inter-molecule.
    // Do single linkage clustering of closest pair at this distance and count the number of 
    // links that are inter and intra.
    // Convert molecules for clustering
    ArrayList<ClusterPoint> points = new ArrayList<ClusterPoint>(molecules.size());
    for (Molecule m : molecules) // Precision was used to store the molecule ID
    points.add(ClusterPoint.newClusterPoint((int) m.precision, m.x, m.y, m.photons));
    ClusteringEngine engine = new ClusteringEngine(Prefs.getThreads(), ClusteringAlgorithm.PARTICLE_SINGLE_LINKAGE, new IJTrackProgress());
    IJ.showStatus("Clustering to check inter-molecule distances");
    engine.setTrackJoins(true);
    ArrayList<Cluster> clusters = engine.findClusters(points, intraHist[0][p99]);
    IJ.showStatus("");
    if (clusters != null) {
        double[] intraIdDistances = engine.getIntraIdDistances();
        double[] interIdDistances = engine.getInterIdDistances();
        int all = interIdDistances.length + intraIdDistances.length;
        log("  * Fraction of inter-molecule particle linkage @ %s nm = %s %%", Utils.rounded(intraHist[0][p99], 4), (all > 0) ? Utils.rounded(100.0 * interIdDistances.length / all, 4) : "0");
        // Show a double cumulative histogram plot
        double[][] intraIdHist = Maths.cumulativeHistogram(intraIdDistances, false);
        double[][] interIdHist = Maths.cumulativeHistogram(interIdDistances, false);
        // Plot
        String title = TITLE + " molecule linkage distance";
        Plot2 plot = new Plot2(title, "Distance", "Frequency", intraIdHist[0], intraIdHist[1]);
        double max = (intraIdHist[1].length > 0) ? intraIdHist[1][intraIdHist[1].length - 1] : 0;
        if (interIdHist[1].length > 0)
            max = FastMath.max(max, interIdHist[1][interIdHist[1].length - 1]);
        plot.setLimits(0, intraIdHist[0][intraIdHist[0].length - 1], 0, max);
        plot.setColor(Color.blue);
        plot.addPoints(interIdHist[0], interIdHist[1], Plot2.LINE);
        plot.setColor(Color.black);
        Utils.display(title, plot);
    } else {
        log("Aborted clustering to check inter-molecule distances");
    }
}
Also used : TDoubleArrayList(gnu.trove.list.array.TDoubleArrayList) ArrayList(java.util.ArrayList) IJTrackProgress(gdsc.core.ij.IJTrackProgress) Cluster(gdsc.core.clustering.Cluster) Plot2(ij.gui.Plot2) WeightedObservedPoint(org.apache.commons.math3.fitting.WeightedObservedPoint) ClusterPoint(gdsc.core.clustering.ClusterPoint) ClusteringEngine(gdsc.core.clustering.ClusteringEngine) ClusterPoint(gdsc.core.clustering.ClusterPoint)

Example 4 with Cluster

use of org.apache.commons.math3.ml.clustering.Cluster 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));
    }
}
Also used : ByteProcessor(ij.process.ByteProcessor) TDoubleArrayList(gnu.trove.list.array.TDoubleArrayList) ArrayList(java.util.ArrayList) MaskDistribution(gdsc.smlm.model.MaskDistribution) Well19937c(org.apache.commons.math3.random.Well19937c) RandomGenerator(org.apache.commons.math3.random.RandomGenerator) MemoryPeakResults(gdsc.smlm.results.MemoryPeakResults) NullSource(gdsc.smlm.results.NullSource) RandomDataGenerator(org.apache.commons.math3.random.RandomDataGenerator) UniformDistribution(gdsc.smlm.model.UniformDistribution) StoredDataStatistics(gdsc.core.utils.StoredDataStatistics) Statistics(gdsc.core.utils.Statistics) StoredDataStatistics(gdsc.core.utils.StoredDataStatistics) DescriptiveStatistics(org.apache.commons.math3.stat.descriptive.DescriptiveStatistics) WeightedObservedPoint(org.apache.commons.math3.fitting.WeightedObservedPoint) ClusterPoint(gdsc.core.clustering.ClusterPoint) StoredData(gdsc.core.utils.StoredData)

Example 5 with Cluster

use of org.apache.commons.math3.ml.clustering.Cluster in project GDSC-SMLM by aherbert.

the class PCPALMClusters method run.

/*
	 * (non-Javadoc)
	 * 
	 * @see ij.plugin.PlugIn#run(java.lang.String)
	 */
public void run(String arg) {
    SMLMUsageTracker.recordPlugin(this.getClass(), arg);
    if (!showDialog())
        return;
    PCPALMMolecules.logSpacer();
    Utils.log(TITLE);
    PCPALMMolecules.logSpacer();
    long start = System.currentTimeMillis();
    HistogramData histogramData;
    if (fileInput) {
        histogramData = loadHistogram(histogramFile);
    } else {
        histogramData = doClustering();
    }
    if (histogramData == null)
        return;
    float[][] hist = histogramData.histogram;
    // Create a histogram of the cluster sizes
    String title = TITLE + " Molecules/cluster";
    String xTitle = "Molecules/cluster";
    String yTitle = "Frequency";
    // Create the data required for fitting and plotting
    float[] xValues = Utils.createHistogramAxis(hist[0]);
    float[] yValues = Utils.createHistogramValues(hist[1]);
    // Plot the histogram
    float yMax = Maths.max(yValues);
    Plot2 plot = new Plot2(title, xTitle, yTitle, xValues, yValues);
    if (xValues.length > 0) {
        double xPadding = 0.05 * (xValues[xValues.length - 1] - xValues[0]);
        plot.setLimits(xValues[0] - xPadding, xValues[xValues.length - 1] + xPadding, 0, yMax * 1.05);
    }
    Utils.display(title, plot);
    HistogramData noiseData = loadNoiseHistogram(histogramData);
    if (noiseData != null) {
        if (subtractNoise(histogramData, noiseData)) {
            // Update the histogram
            title += " (noise subtracted)";
            xValues = Utils.createHistogramAxis(hist[0]);
            yValues = Utils.createHistogramValues(hist[1]);
            yMax = Maths.max(yValues);
            plot = new Plot2(title, xTitle, yTitle, xValues, yValues);
            if (xValues.length > 0) {
                double xPadding = 0.05 * (xValues[xValues.length - 1] - xValues[0]);
                plot.setLimits(xValues[0] - xPadding, xValues[xValues.length - 1] + xPadding, 0, yMax * 1.05);
            }
            Utils.display(title, plot);
            // Automatically save
            if (autoSave) {
                String newFilename = Utils.replaceExtension(histogramData.filename, ".noise.tsv");
                if (saveHistogram(histogramData, newFilename)) {
                    Utils.log("Saved noise-subtracted histogram to " + newFilename);
                }
            }
        }
    }
    // Fit the histogram
    double[] fitParameters = fitBinomial(histogramData);
    if (fitParameters != null) {
        // Add the binomial to the histogram
        int n = (int) fitParameters[0];
        double p = fitParameters[1];
        Utils.log("Optimal fit : N=%d, p=%s", n, Utils.rounded(p));
        BinomialDistribution dist = new BinomialDistribution(n, p);
        // A zero-truncated binomial was fitted.
        // pi is the adjustment factor for the probability density.
        double pi = 1 / (1 - dist.probability(0));
        if (!fileInput) {
            // Calculate the estimated number of clusters from the observed molecules:
            // Actual = (Observed / p-value) / N
            final double actual = (nMolecules / p) / n;
            Utils.log("Estimated number of clusters : (%d / %s) / %d = %s", nMolecules, Utils.rounded(p), n, Utils.rounded(actual));
        }
        double[] x = new double[n + 2];
        double[] y = new double[n + 2];
        // Scale the values to match those on the histogram
        final double normalisingFactor = count * pi;
        for (int i = 0; i <= n; i++) {
            x[i] = i + 0.5;
            y[i] = dist.probability(i) * normalisingFactor;
        }
        x[n + 1] = n + 1.5;
        y[n + 1] = 0;
        // Redraw the plot since the limits may have changed
        plot = new Plot2(title, xTitle, yTitle, xValues, yValues);
        double xPadding = 0.05 * (xValues[xValues.length - 1] - xValues[0]);
        plot.setLimits(xValues[0] - xPadding, xValues[xValues.length - 1] + xPadding, 0, Maths.maxDefault(yMax, y) * 1.05);
        plot.setColor(Color.magenta);
        plot.addPoints(x, y, Plot2.LINE);
        plot.addPoints(x, y, Plot2.CIRCLE);
        plot.setColor(Color.black);
        Utils.display(title, plot);
    }
    double seconds = (System.currentTimeMillis() - start) / 1000.0;
    String msg = TITLE + " complete : " + seconds + "s";
    IJ.showStatus(msg);
    Utils.log(msg);
    return;
}
Also used : Plot2(ij.gui.Plot2) BinomialDistribution(org.apache.commons.math3.distribution.BinomialDistribution) ClusterPoint(gdsc.core.clustering.ClusterPoint)

Aggregations

ArrayList (java.util.ArrayList)7 ClusterPoint (gdsc.core.clustering.ClusterPoint)5 Plot2 (ij.gui.Plot2)3 WeightedObservedPoint (org.apache.commons.math3.fitting.WeightedObservedPoint)3 PointValuePair (org.apache.commons.math3.optim.PointValuePair)3 Material (com.jme3.material.Material)2 ColorRGBA (com.jme3.math.ColorRGBA)2 Geometry (com.jme3.scene.Geometry)2 Mesh (com.jme3.scene.Mesh)2 Node (com.jme3.scene.Node)2 VertexBuffer (com.jme3.scene.VertexBuffer)2 TDoubleArrayList (gnu.trove.list.array.TDoubleArrayList)2 PrintWriter (java.io.PrintWriter)2 Collections (java.util.Collections)2 HashMap (java.util.HashMap)2 HashSet (java.util.HashSet)2 List (java.util.List)2 Map (java.util.Map)2 Set (java.util.Set)2 Collectors (java.util.stream.Collectors)2