Examples with ByteProcessor ij.process.ByteProcessor used on opensource projects

Example 86 with ByteProcessor

use of ij.process.ByteProcessor in project TrakEM2 by trakem2.

the class ImageArrayConverter method ImageToFloatArray2D.

public static FloatArray2D ImageToFloatArray2D(ImageProcessor ip) {
    FloatArray2D image;
    Object pixelArray = ip.getPixels();
    int count = 0;
    if (ip instanceof ByteProcessor) {
        image = new FloatArray2D(ip.getWidth(), ip.getHeight());
        byte[] pixels = (byte[]) pixelArray;
        for (int y = 0; y < ip.getHeight(); y++) for (int x = 0; x < ip.getWidth(); x++) image.data[count] = pixels[count++] & 0xff;
    } else if (ip instanceof ShortProcessor) {
        image = new FloatArray2D(ip.getWidth(), ip.getHeight());
        short[] pixels = (short[]) pixelArray;
        for (int y = 0; y < ip.getHeight(); y++) for (int x = 0; x < ip.getWidth(); x++) image.data[count] = pixels[count++] & 0xffff;
    } else if (ip instanceof FloatProcessor) {
        image = new FloatArray2D(ip.getWidth(), ip.getHeight());
        float[] pixels = (float[]) pixelArray;
        for (int y = 0; y < ip.getHeight(); y++) for (int x = 0; x < ip.getWidth(); x++) image.data[count] = pixels[count++];
    } else // RGB
    {
        image = new FloatArray2D(ip.getWidth(), ip.getHeight());
        int[] pixels = (int[]) pixelArray;
    // still unknown how to do...
    /*
            for (int y = 0; y < ip.getHeight(); y++)
                for (int x = 0; x < ip.getWidth(); x++)
                        image[x][y] = pixels[count++];// & 0xffffff;
            */
    }
    return image;
}

Example 87 with ByteProcessor

use of ij.process.ByteProcessor in project GDSC-SMLM by aherbert.

the class SplitResults method splitResults.

private void splitResults(MemoryPeakResults results, ImageProcessor ip) {
    IJ.showStatus("Splitting " + TextUtils.pleural(results.size(), "result"));
    // Create an object mask
    final ObjectAnalyzer objectAnalyzer = new ObjectAnalyzer(ip, false);
    final int maxx = ip.getWidth();
    final int maxy = ip.getHeight();
    final Rectangle bounds = results.getBounds();
    final double ox = bounds.getX();
    final double oy = bounds.getY();
    final double scaleX = bounds.getWidth() / maxx;
    final double scaleY = bounds.getHeight() / maxy;
    // Create a results set for each object
    final int maxObject = objectAnalyzer.getMaxObject();
    final MemoryPeakResults[] resultsSet = new MemoryPeakResults[maxObject + 1];
    for (int object = 0; object <= maxObject; object++) {
        final MemoryPeakResults newResults = new MemoryPeakResults();
        newResults.copySettings(results);
        newResults.setName(results.getName() + " " + object);
        resultsSet[object] = newResults;
    }
    final int[] mask = objectAnalyzer.getObjectMask();
    if (settings.showObjectMask) {
        final ImageProcessor objectIp = (maxObject <= 255) ? new ByteProcessor(maxx, maxy) : new ShortProcessor(maxx, maxy);
        for (int i = 0; i < mask.length; i++) {
            objectIp.set(i, mask[i]);
        }
        final ImagePlus imp = ImageJUtils.display(settings.objectMask + " Objects", objectIp);
        imp.setDisplayRange(0, maxObject);
        imp.updateAndDraw();
    }
    // Process the results mapping them to their objects
    final Counter i = new Counter();
    final int size = results.size();
    final int step = ImageJUtils.getProgressInterval(size);
    results.forEach(DistanceUnit.PIXEL, (XyrResultProcedure) (xx, yy, result) -> {
        if (i.incrementAndGet() % step == 0) {
            IJ.showProgress(i.getCount(), size);
        }
        // Map to the mask objects
        final int object;
        final int x = (int) ((xx - ox) / scaleX);
        final int y = (int) ((yy - oy) / scaleY);
        if (x < 0 || x >= maxx || y < 0 || y >= maxy) {
            object = 0;
        } else {
            final int index = y * maxx + x;
            // is within the bounds of the image processor?
            if (index < 0 || index >= mask.length) {
                object = 0;
            } else {
                object = mask[index];
            }
        }
        resultsSet[object].add(result);
    });
    IJ.showProgress(1);
    // Add the new results sets to memory
    i.reset();
    for (int object = (settings.nonMaskDataset) ? 0 : 1; object <= maxObject; object++) {
        if (resultsSet[object].isNotEmpty()) {
            MemoryPeakResults.addResults(resultsSet[object]);
            i.increment();
        }
    }
    IJ.showStatus("Split " + TextUtils.pleural(results.size(), "result") + " into " + TextUtils.pleural(i.getCount(), "set"));
}

Example 88 with ByteProcessor

use of ij.process.ByteProcessor in project GDSC-SMLM by aherbert.

the class DrawClusters method run.

@Override
public void run(String arg) {
    SmlmUsageTracker.recordPlugin(this.getClass(), arg);
    if (MemoryPeakResults.isMemoryEmpty()) {
        IJ.error(TITLE, "No localisations in memory");
        return;
    }
    if (!showDialog()) {
        return;
    }
    // Load the results
    final MemoryPeakResults results = ResultsManager.loadInputResults(settings.inputOption, false, DistanceUnit.PIXEL);
    if (MemoryPeakResults.isEmpty(results)) {
        IJ.error(TITLE, "No results could be loaded");
        return;
    }
    // Get the traces
    final Trace[] traces = TraceManager.convert(results);
    if (traces == null || traces.length == 0) {
        IJ.error(TITLE, "No traces could be loaded");
        return;
    }
    // Filter traces to a min size
    int maxFrame = 0;
    int count = 0;
    final int myMaxSize = (settings.maxSize < settings.minSize) ? Integer.MAX_VALUE : settings.maxSize;
    final boolean myDrawLines = myMaxSize > 1 && settings.drawLines;
    for (int i = 0; i < traces.length; i++) {
        if (settings.expandToSingles) {
            traces[i].expandToSingles();
        }
        if (traces[i].size() >= settings.minSize && traces[i].size() <= myMaxSize) {
            traces[count++] = traces[i];
            traces[i].sort();
            if (maxFrame < traces[i].getTail().getFrame()) {
                maxFrame = traces[i].getTail().getFrame();
            }
        }
    }
    if (count == 0) {
        IJ.error(TITLE, "No traces achieved the size limits");
        return;
    }
    final String msg = String.format(TITLE + ": %d / %s (%s)", count, TextUtils.pleural(traces.length, "trace"), TextUtils.pleural(results.size(), "localisation"));
    IJ.showStatus(msg);
    final Rectangle bounds = results.getBounds(true);
    ImagePlus imp = WindowManager.getImage(settings.title);
    boolean isUseStackPosition = settings.useStackPosition;
    if (imp == null) {
        // Create a default image using 100 pixels as the longest edge
        final double maxD = (bounds.width > bounds.height) ? bounds.width : bounds.height;
        int width;
        int height;
        if (maxD == 0) {
            // Note that imageSize can be zero (for auto sizing)
            width = height = (settings.imageSize == 0) ? 20 : settings.imageSize;
        } else if (settings.imageSize == 0) {
            // Note that imageSize can be zero (for auto sizing)
            width = bounds.width;
            height = bounds.height;
        } else {
            width = (int) (settings.imageSize * bounds.width / maxD);
            height = (int) (settings.imageSize * bounds.height / maxD);
        }
        final ByteProcessor bp = new ByteProcessor(width, height);
        if (isUseStackPosition) {
            final ImageStack stack = new ImageStack(width, height, maxFrame);
            for (int i = 1; i <= maxFrame; i++) {
                // Do not clone as the image is empty
                stack.setPixels(bp.getPixels(), i);
            }
            imp = ImageJUtils.display(TITLE, stack);
        } else {
            imp = ImageJUtils.display(TITLE, bp);
        }
        // Enlarge
        final ImageWindow iw = imp.getWindow();
        for (int i = 9; i-- > 0 && iw.getWidth() < 500 && iw.getHeight() < 500; ) {
            iw.getCanvas().zoomIn(imp.getWidth() / 2, imp.getHeight() / 2);
        }
    // Check if the image has enough frames for all the traces
    } else if (maxFrame > imp.getNFrames()) {
        isUseStackPosition = false;
    }
    final float xScale = (float) (imp.getWidth() / bounds.getWidth());
    final float yScale = (float) (imp.getHeight() / bounds.getHeight());
    // Create ROIs and store data to sort them
    final Roi[] rois = new Roi[count];
    final int[][] frames = (isUseStackPosition) ? new int[count][] : null;
    final int[] indices = SimpleArrayUtils.natural(count);
    final double[] values = new double[count];
    for (int i = 0; i < count; i++) {
        final Trace trace = traces[i];
        final int npoints = trace.size();
        final float[] xpoints = new float[npoints];
        final float[] ypoints = new float[npoints];
        int ii = 0;
        if (frames != null) {
            frames[i] = new int[npoints];
        }
        for (int k = 0; k < trace.size(); k++) {
            final PeakResult result = trace.get(k);
            xpoints[ii] = (result.getXPosition() - bounds.x) * xScale;
            ypoints[ii] = (result.getYPosition() - bounds.y) * yScale;
            if (frames != null) {
                frames[i][ii] = result.getFrame();
            }
            ii++;
        }
        Roi roi;
        if (myDrawLines) {
            roi = new PolygonRoi(xpoints, ypoints, npoints, Roi.POLYLINE);
            if (settings.splineFit) {
                ((PolygonRoi) roi).fitSpline();
            }
        } else {
            roi = new OffsetPointRoi(xpoints, ypoints, npoints);
            ((PointRoi) roi).setShowLabels(false);
        }
        rois[i] = roi;
        switch(settings.sort) {
            case // Sort by ID
            1:
                values[i] = traces[i].getId();
                break;
            case // Sort by time
            2:
                values[i] = traces[i].getHead().getFrame();
                break;
            case // Sort by size descending
            3:
                values[i] = -traces[i].size();
                break;
            case // Sort by length descending
            4:
                values[i] = -roi.getLength();
                break;
            case // Mean Square Displacement
            5:
                values[i] = -traces[i].getMsd();
                break;
            case // Mean / Frame
            6:
                values[i] = -traces[i].getMeanDistance();
                break;
            // No sort
            case 0:
            default:
                break;
        }
    }
    if (settings.sort > 0) {
        SortUtils.sortIndices(indices, values, true);
    }
    // Draw the traces as ROIs on an overlay
    final Overlay o = new Overlay();
    final LUT lut = LutHelper.createLut(settings.lut);
    final double scale = 256.0 / count;
    if (frames != null) {
        // Add the tracks on the frames containing the results
        final boolean isHyperStack = imp.isDisplayedHyperStack();
        for (int i = 0; i < count; i++) {
            final int index = indices[i];
            final Color c = LutHelper.getColour(lut, (int) (i * scale));
            final PolygonRoi roi = (PolygonRoi) rois[index];
            roi.setFillColor(c);
            roi.setStrokeColor(c);
            // roi.setStrokeWidth(settings.lineWidth);
            roi.updateWideLine(settings.lineWidth);
            final FloatPolygon fp = roi.getNonSplineFloatPolygon();
            // For each frame in the track, add the ROI track and a point ROI for the current position
            for (int j = 0; j < frames[index].length; j++) {
                addToOverlay(o, (Roi) roi.clone(), isHyperStack, frames[index][j]);
                final PointRoi pointRoi = new OffsetPointRoi(fp.xpoints[j], fp.ypoints[j]);
                pointRoi.setPointType(3);
                pointRoi.setFillColor(c);
                pointRoi.setStrokeColor(Color.black);
                addToOverlay(o, pointRoi, isHyperStack, frames[index][j]);
            }
        }
    } else {
        // Add the tracks as a single overlay
        for (int i = 0; i < count; i++) {
            final Roi roi = rois[indices[i]];
            roi.setStrokeColor(new Color(lut.getRGB((int) (i * scale))));
            // roi.setStrokeWidth(settings.lineWidth);
            roi.updateWideLine(settings.lineWidth);
            o.add(roi);
        }
    }
    imp.setOverlay(o);
    IJ.showStatus(msg);
}

Example 89 with ByteProcessor

use of ij.process.ByteProcessor 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 (settings.blinkingDistribution == 3) {
        log("  - Clusters = %d", settings.numberOfMolecules);
        log("  - Simulation size = %s um", MathUtils.rounded(settings.simulationSize, 4));
        log("  - Molecules/cluster = %s", MathUtils.rounded(settings.blinkingRate, 4));
        log("  - Blinking distribution = %s", Settings.BLINKING_DISTRIBUTION[settings.blinkingDistribution]);
        log("  - p-Value = %s", MathUtils.rounded(settings.pvalue, 4));
    } else {
        log("  - Molecules = %d", settings.numberOfMolecules);
        log("  - Simulation size = %s um", MathUtils.rounded(settings.simulationSize, 4));
        log("  - Blinking rate = %s", MathUtils.rounded(settings.blinkingRate, 4));
        log("  - Blinking distribution = %s", Settings.BLINKING_DISTRIBUTION[settings.blinkingDistribution]);
    }
    log("  - Average precision = %s nm", MathUtils.rounded(settings.sigmaS, 4));
    log("  - Clusters simulation = " + Settings.CLUSTER_SIMULATION[settings.clusterSimulation]);
    if (settings.clusterSimulation > 0) {
        log("  - Cluster number = %s +/- %s", MathUtils.rounded(settings.clusterNumber, 4), MathUtils.rounded(settings.clusterNumberStdDev, 4));
        log("  - Cluster radius = %s nm", MathUtils.rounded(settings.clusterRadius, 4));
    }
    final double nmPerPixel = 100;
    final double width = settings.simulationSize * 1000.0;
    final UniformRandomProvider rng = UniformRandomProviders.create();
    final UniformDistribution dist = new UniformDistribution(null, new double[] { width, width, 0 }, rng.nextInt());
    final NormalizedGaussianSampler gauss = SamplerUtils.createNormalizedGaussianSampler(rng);
    settings.molecules = new ArrayList<>(settings.numberOfMolecules);
    // Create some dummy results since the calibration is required for later analysis
    settings.results = new MemoryPeakResults(PsfHelper.create(PSFType.CUSTOM));
    settings.results.setCalibration(CalibrationHelper.create(nmPerPixel, 1, 100));
    settings.results.setSource(new NullSource("Molecule Simulation"));
    settings.results.begin();
    int count = 0;
    // Generate a sequence of coordinates
    final ArrayList<double[]> xyz = new ArrayList<>((int) (settings.numberOfMolecules * 1.1));
    final Statistics statsRadius = new Statistics();
    final Statistics statsSize = new Statistics();
    final String maskTitle = TITLE + " Cluster Mask";
    ByteProcessor bp = null;
    double maskScale = 0;
    if (settings.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 = settings.lowResolutionImageSize;
        int[] mask = null;
        // scale is in nm/pixel
        maskScale = width / maskSize;
        final ArrayList<double[]> clusterCentres = new ArrayList<>();
        int totalSteps = 1 + (int) Math.ceil(settings.numberOfMolecules / settings.clusterNumber);
        if (settings.clusterSimulation == 2 || settings.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, rng);
            double[] centre;
            IJ.showStatus("Computing clusters mask");
            final int roiRadius = (int) Math.round((settings.clusterRadius * 2) / maskScale);
            if (settings.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 * MathUtils.pow2(settings.clusterRadius / maskScale)));
            }
            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
                final double cx = centre[0] / maskScale;
                final double cy = centre[1] / maskScale;
                fillMask(mask, maskSize, (int) cx, (int) cy, roiRadius, 0);
                try {
                    maskDistribution = new MaskDistribution(mask, maskSize, maskSize, 0, maskScale, maskScale, rng);
                } catch (final IllegalArgumentException ex) {
                    // 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;
                }
            }
            ImageJUtils.finished();
        } else {
            // Pick centres randomly from the distribution
            while (clusterCentres.size() < totalSteps) {
                clusterCentres.add(dist.next());
            }
        }
        final double scaledRadius = settings.clusterRadius / maskScale;
        if (settings.showClusterMask || settings.clusterSimulation == 3) {
            // Show the mask for the clusters
            if (mask == null) {
                mask = new int[maskSize * maskSize];
            } else {
                Arrays.fill(mask, 0);
            }
            final int roiRadius = (int) Math.round(scaledRadius);
            for (final double[] c : clusterCentres) {
                final double cx = c[0] / maskScale;
                final double cy = c[1] / maskScale;
                fillMask(mask, maskSize, (int) cx, (int) cy, roiRadius, 1);
            }
            if (settings.clusterSimulation == 3) {
                // We have the mask. Now pick points at random from the mask.
                final MaskDistribution maskDistribution = new MaskDistribution(mask, maskSize, maskSize, 0, maskScale, maskScale, rng);
                // Allocate each molecule position to a parent circle so defining clusters.
                final int[][] clusters = new int[clusterCentres.size()][];
                final int[] clusterSize = new int[clusters.length];
                for (int i = 0; i < settings.numberOfMolecules; i++) {
                    final 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++) {
                        final 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
                    final double[] com = new double[2];
                    for (int n = 0; n < size; n++) {
                        final 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++) {
                        final double dx = xyz.get(clusters[j][n])[0] - com[0];
                        final double dy = xyz.get(clusters[j][n])[1] - com[1];
                        statsRadius.add(Math.sqrt(dx * dx + dy * dy));
                    }
                }
            }
            if (settings.showClusterMask) {
                bp = new ByteProcessor(maskSize, maskSize);
                for (int i = 0; i < mask.length; i++) {
                    if (mask[i] != 0) {
                        bp.set(i, 128);
                    }
                }
                ImageJUtils.display(maskTitle, bp);
            }
        }
        // Use the simulated cluster centres to create clusters of the desired size
        if (settings.clusterSimulation == 1 || settings.clusterSimulation == 2) {
            for (final double[] clusterCentre : clusterCentres) {
                final int clusterN = (int) Math.round((settings.clusterNumberStdDev > 0) ? settings.clusterNumber + gauss.sample() * settings.clusterNumberStdDev : settings.clusterNumber);
                if (clusterN < 1) {
                    continue;
                }
                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.
                    final double[] com = new double[3];
                    int size = 0;
                    while (size < clusterN) {
                        // Generate a random point within a circle uniformly
                        // http://stackoverflow.com/questions/5837572/generate-a-random-point-within-a-circle-uniformly
                        final double t = 2.0 * Math.PI * rng.nextDouble();
                        final double u = rng.nextDouble() + rng.nextDouble();
                        final double r = settings.clusterRadius * ((u > 1) ? 2 - u : u);
                        final double x = r * Math.cos(t);
                        final double y = r * Math.sin(t);
                        final double[] xy = new double[] { clusterCentre[0] + x, clusterCentre[1] + y };
                        xyz.add(xy);
                        for (int k = 0; k < 2; k++) {
                            com[k] += xy[k];
                        }
                        size++;
                    }
                    // 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(size);
                    if (size == 1) {
                        statsRadius.add(0);
                    } else {
                        for (int k = 0; k < 2; k++) {
                            com[k] /= size;
                        }
                        while (size > 0) {
                            final double dx = xyz.get(xyz.size() - size)[0] - com[0];
                            final double dy = xyz.get(xyz.size() - size)[1] - com[1];
                            statsRadius.add(Math.sqrt(dx * dx + dy * dy));
                            size--;
                        }
                    }
                }
            }
        }
    } else {
        // Random distribution
        for (int i = 0; i < settings.numberOfMolecules; 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 = settings.sigmaS / Math.sqrt(2);
    // Show optional histograms
    StoredDataStatistics intraDistances = null;
    StoredData blinks = null;
    if (settings.showHistograms) {
        final int capacity = (int) (xyz.size() * settings.blinkingRate);
        intraDistances = new StoredDataStatistics(capacity);
        blinks = new StoredData(capacity);
    }
    final Statistics statsSigma = new Statistics();
    for (int i = 0; i < xyz.size(); i++) {
        int occurrences = getBlinks(rng, settings.blinkingRate);
        if (blinks != null) {
            blinks.add(occurrences);
        }
        final int size = settings.molecules.size();
        // Get coordinates in nm
        final double[] moleculeXyz = xyz.get(i);
        if (bp != null && occurrences > 0) {
            bp.putPixel((int) Math.round(moleculeXyz[0] / maskScale), (int) Math.round(moleculeXyz[1] / maskScale), 255);
        }
        while (occurrences-- > 0) {
            final double[] localisationXy = Arrays.copyOf(moleculeXyz, 2);
            // Add random precision
            if (sigma1D > 0) {
                final double dx = gauss.sample() * sigma1D;
                final double dy = gauss.sample() * 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];
            settings.molecules.add(new Molecule(x, y, i, 1));
            // Store in pixels
            final float xx = (float) (x / nmPerPixel);
            final float yy = (float) (y / nmPerPixel);
            final float[] params = PeakResult.createParams(0, 0, xx, yy, 0);
            settings.results.add(i + 1, (int) xx, (int) yy, 0, 0, 0, 0, params, null);
        }
        if (settings.molecules.size() > size) {
            count++;
            if (intraDistances != null) {
                final int newCount = settings.molecules.size() - size;
                if (newCount == 1) {
                    // No intra-molecule distances
                    continue;
                }
                // Get the distance matrix between these molecules
                final double[][] matrix = new double[newCount][newCount];
                for (int ii = size, x = 0; ii < settings.molecules.size(); ii++, x++) {
                    for (int jj = size + 1, y = 1; jj < settings.molecules.size(); jj++, y++) {
                        final double d2 = settings.molecules.get(ii).distance2(settings.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));
            }
        }
    }
    settings.results.end();
    if (bp != null) {
        final ImagePlus imp = ImageJUtils.display(maskTitle, bp);
        final Calibration cal = imp.getCalibration();
        cal.setUnit("nm");
        cal.pixelWidth = cal.pixelHeight = maskScale;
    }
    log("Simulation results");
    log("  * Molecules = %d (%d activated)", xyz.size(), count);
    log("  * Blinking rate = %s", MathUtils.rounded((double) settings.molecules.size() / xyz.size(), 4));
    log("  * Precision (Mean-displacement) = %s nm", (statsSigma.getN() > 0) ? MathUtils.rounded(Math.sqrt(statsSigma.getMean()), 4) : "0");
    if (intraDistances != null) {
        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);
            final double[][] intraHist = plot(intraDistances, "Intra-molecule particle linkage distance", false);
            // Determine 95th and 99th percentile
            // Will not be null as we requested a non-integer histogram.
            int p99 = intraHist[0].length - 1;
            final double limit1 = 0.99 * intraHist[1][p99];
            final 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)", MathUtils.rounded(intraDistances.getMean(), 4), MathUtils.rounded(intraHist[0][p95], 4), MathUtils.rounded(intraHist[0][p99], 4), MathUtils.rounded(intraHist[0][intraHist[0].length - 1], 4));
            if (settings.distanceAnalysis) {
                performDistanceAnalysis(intraHist, p99);
            }
        }
    }
    if (settings.clusterSimulation > 0) {
        log("  * Cluster number = %s +/- %s", MathUtils.rounded(statsSize.getMean(), 4), MathUtils.rounded(statsSize.getStandardDeviation(), 4));
        log("  * Cluster radius = %s +/- %s nm (mean distance to centre-of-mass)", MathUtils.rounded(statsRadius.getMean(), 4), MathUtils.rounded(statsRadius.getStandardDeviation(), 4));
    }
}

Example 90 with ByteProcessor

use of ij.process.ByteProcessor in project GDSC-SMLM by aherbert.

the class PsfCreator method plotSignalAtSpecifiedSd.

/**
 * Show a plot of the amount of signal within N x SD for each z position. This indicates how much
 * the PSF has spread from the original Gaussian shape.
 *
 * @param psf The PSF
 * @param fittedSd The width of the PSF (in pixels)
 * @param factor The factor to use
 * @param slice The slice used to create the label
 */
private void plotSignalAtSpecifiedSd(ImageStack psf, double fittedSd, double factor, int slice) {
    if (signalZ == null) {
        // Get the bounds
        final int radius = (int) Math.round(fittedSd * factor);
        final int min = Math.max(0, psf.getWidth() / 2 - radius);
        final int max = Math.min(psf.getWidth() - 1, psf.getWidth() / 2 + radius);
        // Create a circle mask of the PSF projection
        final ByteProcessor circle = new ByteProcessor(max - min + 1, max - min + 1);
        circle.setColor(255);
        circle.fillOval(0, 0, circle.getWidth(), circle.getHeight());
        final byte[] mask = (byte[]) circle.getPixels();
        // Sum the pixels within the mask for each slice
        signalZ = new double[psf.getSize()];
        signal = new double[psf.getSize()];
        for (int i = 0; i < psf.getSize(); i++) {
            double sum = 0;
            final float[] data = (float[]) psf.getProcessor(i + 1).getPixels();
            for (int y = min, ii = 0; y <= max; y++) {
                int index = y * psf.getWidth() + min;
                for (int x = min; x <= max; x++, ii++, index++) {
                    if (mask[ii] != 0 && data[index] > 0) {
                        sum += data[index];
                    }
                }
            }
            double total = 0;
            for (final float f : data) {
                if (f > 0) {
                    total += f;
                }
            }
            signalZ[i] = i + 1;
            signal[i] = 100 * sum / total;
        }
        signalTitle = String.format("%% PSF signal at %s x SD", MathUtils.rounded(factor, 3));
        signalLimits = MathUtils.limits(signal);
    }
    // Plot the sum
    final boolean alignWindows = (WindowManager.getFrame(signalTitle) == null);
    final double total = signal[slice - 1];
    final Plot plot = new Plot(signalTitle, "z", "Signal");
    plot.addPoints(signalZ, signal, Plot.LINE);
    plot.addLabel(0, 0, String.format("Total = %s. z = %s nm", MathUtils.rounded(total), MathUtils.rounded((slice - zCentre) * settings.getNmPerSlice())));
    plot.setColor(Color.green);
    plot.drawLine(slice, signalLimits[0], slice, signalLimits[1]);
    plot.setColor(Color.blue);
    final PlotWindow plotWindow = ImageJUtils.display(signalTitle, plot);
    if (alignWindows && plotWindow != null) {
        final PlotWindow otherWindow = getPlot(TITLE_AMPLITUDE);
        if (otherWindow != null) {
            // Put the two plots tiled together so both are visible
            final Point l = plotWindow.getLocation();
            l.x = otherWindow.getLocation().x + otherWindow.getWidth();
            l.y = otherWindow.getLocation().y;
            plotWindow.setLocation(l);
        }
    }
}