Search in sources :

Example 1 with MaskDistribution

use of gdsc.smlm.model.MaskDistribution in project GDSC-SMLM by aherbert.

the class FilterResults method filterResults.

/**
	 * Apply the filters to the data
	 */
private void filterResults() {
    checkLimits();
    MemoryPeakResults newResults = new MemoryPeakResults();
    newResults.copySettings(results);
    newResults.setName(results.getName() + " Filtered");
    // Initialise the mask
    ByteProcessor mask = getMask(filterSettings.maskTitle);
    MaskDistribution maskFilter = null;
    final float centreX = results.getBounds().width / 2.0f;
    final float centreY = results.getBounds().height / 2.0f;
    if (mask != null) {
        double scaleX = (double) results.getBounds().width / mask.getWidth();
        double scaleY = (double) results.getBounds().height / mask.getHeight();
        maskFilter = new MaskDistribution((byte[]) mask.getPixels(), mask.getWidth(), mask.getHeight(), 0, scaleX, scaleY);
    }
    int i = 0;
    final int size = results.size();
    final double maxVariance = filterSettings.maxPrecision * filterSettings.maxPrecision;
    for (PeakResult result : results.getResults()) {
        if (i % 64 == 0)
            IJ.showProgress(i, size);
        if (getDrift(result) > filterSettings.maxDrift)
            continue;
        if (result.getSignal() < filterSettings.minSignal)
            continue;
        if (getSNR(result) < filterSettings.minSNR)
            continue;
        if (getVariance(result) > maxVariance)
            continue;
        final float width = getWidth(result);
        if (width < filterSettings.minWidth || width > filterSettings.maxWidth)
            continue;
        if (maskFilter != null) {
            // Check the coordinates are inside the mask
            double[] xy = new double[] { result.getXPosition() - centreX, result.getYPosition() - centreY };
            if (!maskFilter.isWithinXY(xy))
                continue;
        }
        // Passed all filters. Add to the results
        newResults.add(result);
    }
    IJ.showProgress(1);
    IJ.showStatus(newResults.size() + " Filtered localisations");
    MemoryPeakResults.addResults(newResults);
}
Also used : ByteProcessor(ij.process.ByteProcessor) MemoryPeakResults(gdsc.smlm.results.MemoryPeakResults) MaskDistribution(gdsc.smlm.model.MaskDistribution) PeakResult(gdsc.smlm.results.PeakResult)

Example 2 with MaskDistribution

use of gdsc.smlm.model.MaskDistribution in project GDSC-SMLM by aherbert.

the class CreateData method createMaskDistribution.

private SpatialDistribution createMaskDistribution(ImagePlus imp, double sliceDepth, boolean updateArea) {
    // Calculate the scale of the mask
    final int w = imp.getWidth();
    final int h = imp.getHeight();
    final double scaleX = (double) settings.size / w;
    final double scaleY = (double) settings.size / h;
    // Use an image for the distribution
    if (imp.getStackSize() > 1) {
        ImageStack stack = imp.getImageStack();
        List<int[]> masks = new ArrayList<int[]>(stack.getSize());
        int[] maxMask = new int[w * h];
        for (int slice = 1; slice <= stack.getSize(); slice++) {
            int[] mask = extractMask(stack.getProcessor(slice));
            if (updateArea) {
                for (int i = 0; i < mask.length; i++) if (mask[i] != 0)
                    maxMask[i] = 1;
            }
            masks.add(mask);
        }
        if (updateArea)
            updateArea(maxMask, w, h);
        if (sliceDepth == 0)
            // Auto configure to the full depth of the simulation
            sliceDepth = settings.depth / masks.size();
        return new MaskDistribution3D(masks, w, h, sliceDepth / settings.pixelPitch, scaleX, scaleY, createRandomGenerator());
    } else {
        int[] mask = extractMask(imp.getProcessor());
        if (updateArea)
            updateArea(mask, w, h);
        return new MaskDistribution(mask, w, h, settings.depth / settings.pixelPitch, scaleX, scaleY, createRandomGenerator());
    }
}
Also used : ImageStack(ij.ImageStack) ArrayList(java.util.ArrayList) MaskDistribution(gdsc.smlm.model.MaskDistribution) MaskDistribution3D(gdsc.smlm.model.MaskDistribution3D)

Example 3 with MaskDistribution

use of gdsc.smlm.model.MaskDistribution 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 4 with MaskDistribution

use of gdsc.smlm.model.MaskDistribution in project GDSC-SMLM by aherbert.

the class PCPALMAnalysis method cropToRoi.

/**
	 * Extract all the PC-PALM molecules that are within the image ROI region. The coordinates bounds are
	 * converted using relative scaling to the limits of the PC-PALM molecules. If a non-rectangular ROI is used then a
	 * mask is extracted and used for the crop. If no image is provided then the full set of molecules is returned.
	 * <p>
	 * Set the area property to the region covered by the molecules.
	 * 
	 * @param imp
	 * @return
	 */
ArrayList<Molecule> cropToRoi(ImagePlus imp) {
    croppedArea = PCPALMMolecules.area;
    if (PCPALMMolecules.molecules == null || PCPALMMolecules.molecules.isEmpty()) {
        return PCPALMMolecules.molecules;
    }
    double pcw = PCPALMMolecules.maxx - PCPALMMolecules.minx;
    double pch = PCPALMMolecules.maxy - PCPALMMolecules.miny;
    Roi roi = (imp == null) ? null : imp.getRoi();
    if (roi == null || !roi.isArea()) {
        log("Roi = %s nm x %s nm = %s um^2", Utils.rounded(pcw, 3), Utils.rounded(pch, 3), Utils.rounded(croppedArea, 3));
        minx = PCPALMMolecules.minx;
        maxx = PCPALMMolecules.maxx;
        miny = PCPALMMolecules.miny;
        maxy = PCPALMMolecules.maxy;
        return PCPALMMolecules.molecules;
    }
    int w = imp.getWidth();
    int h = imp.getHeight();
    Rectangle bounds = roi.getBounds();
    // Construct relative coordinates
    minx = PCPALMMolecules.minx + pcw * ((double) bounds.x / w);
    miny = PCPALMMolecules.miny + pch * ((double) bounds.y / h);
    maxx = PCPALMMolecules.minx + pcw * ((double) (bounds.x + bounds.width) / w);
    maxy = PCPALMMolecules.miny + pch * ((double) (bounds.y + bounds.height) / h);
    final double roix = maxx - minx;
    final double roiy = maxy - miny;
    ArrayList<Molecule> newMolecules = new ArrayList<Molecule>(PCPALMMolecules.molecules.size() / 2);
    // Support non-rectangular ROIs
    if (roi.getMask() != null) {
        MaskDistribution dist = createMaskDistribution(roi.getMask(), roix, roiy);
        double fraction = (double) dist.getSize() / (roi.getMask().getWidth() * roi.getMask().getHeight());
        log("Roi is a mask of %d pixels", dist.getSize());
        croppedArea = fraction * roix * roiy / 1e6;
        log("Roi area %s x %s nm x %s nm = %s um^2", Utils.rounded(fraction), Utils.rounded(roix, 3), Utils.rounded(roiy, 3), Utils.rounded(croppedArea, 3));
        double[] xyz = new double[3];
        // The mask is 0,0 in the centre so offset by this as well
        double xoffset = minx + roix / 2;
        double yoffset = miny + roiy / 2;
        for (Molecule m : PCPALMMolecules.molecules) {
            xyz[0] = m.x - xoffset;
            xyz[1] = m.y - yoffset;
            if (dist.isWithinXY(xyz))
                newMolecules.add(m);
        }
    } else {
        croppedArea = roix * roiy / 1e6;
        log("Roi = %s nm x %s nm = %s um^2", Utils.rounded(roix, 3), Utils.rounded(roiy, 3), Utils.rounded(croppedArea, 3));
        for (Molecule m : PCPALMMolecules.molecules) {
            if (m.x < minx || m.x >= maxx || m.y < miny || m.y >= maxy)
                continue;
            newMolecules.add(m);
        }
    }
    return newMolecules;
}
Also used : Rectangle(java.awt.Rectangle) ArrayList(java.util.ArrayList) MaskDistribution(gdsc.smlm.model.MaskDistribution) Roi(ij.gui.Roi)

Aggregations

MaskDistribution (gdsc.smlm.model.MaskDistribution)4 ArrayList (java.util.ArrayList)3 MemoryPeakResults (gdsc.smlm.results.MemoryPeakResults)2 ByteProcessor (ij.process.ByteProcessor)2 ClusterPoint (gdsc.core.clustering.ClusterPoint)1 Statistics (gdsc.core.utils.Statistics)1 StoredData (gdsc.core.utils.StoredData)1 StoredDataStatistics (gdsc.core.utils.StoredDataStatistics)1 MaskDistribution3D (gdsc.smlm.model.MaskDistribution3D)1 UniformDistribution (gdsc.smlm.model.UniformDistribution)1 NullSource (gdsc.smlm.results.NullSource)1 PeakResult (gdsc.smlm.results.PeakResult)1 TDoubleArrayList (gnu.trove.list.array.TDoubleArrayList)1 ImageStack (ij.ImageStack)1 Roi (ij.gui.Roi)1 Rectangle (java.awt.Rectangle)1 WeightedObservedPoint (org.apache.commons.math3.fitting.WeightedObservedPoint)1 RandomDataGenerator (org.apache.commons.math3.random.RandomDataGenerator)1 RandomGenerator (org.apache.commons.math3.random.RandomGenerator)1 Well19937c (org.apache.commons.math3.random.Well19937c)1