Examples with ImageExtractor uk.ac.sussex.gdsc.core.utils.ImageExtractor used on opensource projects

Example 1 with ImageExtractor

use of uk.ac.sussex.gdsc.core.utils.ImageExtractor in project GDSC-SMLM by aherbert.

the class ResultsImageSampler method createResultSamples.

/**
 * Creates the result samples. Do this by storing the coordinates at the region index.
 */
private void createResultSamples() {
    final TLongObjectHashMap<ResultsSample> map = new TLongObjectHashMap<>(results.size());
    ResultsSample next = ResultsSample.create(-1);
    // For SNR computation
    Object[] pixelArray = stack.getImageArray();
    int width = stack.getWidth();
    int height = stack.getHeight();
    float[] buffer = null;
    // Use a null float[] as this is not used for the getBoxRegionBounds method
    ImageExtractor ie = ImageExtractor.wrap(null, width, height);
    for (final PeakResult p : results.toArray()) {
        // Avoid invalid slices
        if (p.getFrame() < 1 || p.getFrame() > stack.getSize()) {
            continue;
        }
        // Avoid low SNR results. Get the SNR using a 3x3 region around the spot.
        final Rectangle bounds = ie.getBoxRegionBounds((int) p.getXPosition(), (int) p.getYPosition(), 1);
        buffer = ImageJImageConverter.getData(pixelArray[p.getFrame()], width, height, bounds, buffer);
        final Statistics stats = new Statistics();
        stats.add(buffer, 0, bounds.width * bounds.height);
        // SNR will be NaN if the region has no size
        final double snr = stats.getMean() / stats.getStandardDeviation();
        if (snr > 3) {
            final long index = getIndex(p.getXPosition(), p.getYPosition(), p.getFrame());
            ResultsSample current = map.putIfAbsent(index, next);
            if (current == null) {
                // If the return value is null then this is a new insertion.
                // Set the current value as the one we just added and create the next insertion object.
                current = next;
                current.index = index;
                next = ResultsSample.create(-1);
            }
            current.add(p);
        }
    }
    // Create an array of all the sample entries.
    // This is used to sample regions by density.
    data = map.values(new ResultsSample[map.size()]);
}

Example 2 with ImageExtractor

use of uk.ac.sussex.gdsc.core.utils.ImageExtractor in project GDSC-SMLM by aherbert.

the class ImageConverterTest method canGetCropData.

@SeededTest
void canGetCropData(RandomSeed seed) {
    final ImageConverterTestData data = (ImageConverterTestData) dataCache.computeIfAbsent(seed, ImageConverterTest::createData);
    final byte[] bdata = data.bdata;
    final short[] sdata = data.sdata;
    final float[] fdata = data.fdata;
    final UniformRandomProvider rand = RngUtils.create(seed.getSeed());
    final ImageExtractor ie = ImageExtractor.wrap(fdata, w, h);
    for (int i = 0; i < 10; i++) {
        final Rectangle bounds = ie.getBoxRegionBounds(10 + rand.nextInt(w - 20), 10 + rand.nextInt(h - 20), 5 + rand.nextInt(5));
        final FloatProcessor ip = new FloatProcessor(w, h, fdata.clone());
        ip.setRoi(bounds);
        final float[] fe = (float[]) (ip.crop().getPixels());
        Assertions.assertArrayEquals(fe, ImageJImageConverter.getData(bdata, w, h, bounds, null));
        Assertions.assertArrayEquals(fe, ImageJImageConverter.getData(sdata, w, h, bounds, null));
        Assertions.assertArrayEquals(fe, ImageJImageConverter.getData(fdata, w, h, bounds, null));
        // Check the double format
        final double[] de = SimpleArrayUtils.toDouble(fe);
        Assertions.assertArrayEquals(de, ImageJImageConverter.getDoubleData(bdata, w, h, bounds, null));
        Assertions.assertArrayEquals(de, ImageJImageConverter.getDoubleData(sdata, w, h, bounds, null));
        Assertions.assertArrayEquals(de, ImageJImageConverter.getDoubleData(fdata, w, h, bounds, null));
    }
}

Example 3 with ImageExtractor

use of uk.ac.sussex.gdsc.core.utils.ImageExtractor in project GDSC-SMLM by aherbert.

the class PsfCreator method findSpotOverlap.

/**
 * Find all the ROI points that have a box region overlapping with any other spot.
 *
 * @param roiPoints the roi points
 * @param excluded the excluded
 * @return the overlap array
 */
private boolean[] findSpotOverlap(BasePoint[] roiPoints, boolean[] excluded) {
    final int n = roiPoints.length;
    final boolean[] bad = new boolean[n];
    if (n == 1) {
        return bad;
    }
    // Check overlap of box regions
    final int w = imp.getWidth();
    final int h = imp.getHeight();
    final ImageExtractor ie = ImageExtractor.wrap(null, w, h);
    final Rectangle[] regions = new Rectangle[n];
    // Check size if not fitting
    final int size = (settings.getMode() != MODE_FITTING) ? 2 * boxRadius + 1 : Integer.MAX_VALUE;
    for (int i = 0; i < n; i++) {
        if (excluded != null && excluded[i]) {
            continue;
        }
        final Rectangle r = ie.getBoxRegionBounds(roiPoints[i].getXint(), roiPoints[i].getYint(), boxRadius);
        regions[i] = r;
        if (r.width < size || r.height < size) {
            ImageJUtils.log("Warning: Spot %d region extends beyond the image, border pixels will be duplicated", i + 1);
        }
    }
    // Add support for 3D overlap analysis. Only do this if a zRadius has been specified.
    final boolean is3D = (settings.getMode() == MODE_ALIGNMENT && zRadius > 0);
    for (int i = 0; i < n; i++) {
        if (excluded != null && excluded[i]) {
            continue;
        }
        if (bad[i]) {
            continue;
        }
        // Check intersect with others
        for (int j = i; ++j < n; ) {
            if (excluded != null && excluded[j]) {
                continue;
            }
            boolean overlap = regions[i].intersects(regions[j]);
            if (overlap && is3D) {
                // Reset (assume non-overlapping)
                overlap = false;
                // Check for 3D overlap:
                // iiiiiiiiiiiiiiiii
                // jjjjjjjjjjjjjjjjjjj
                final int mini = roiPoints[i].getZint() - zRadius;
                final int maxi = roiPoints[i].getZint() + zRadius;
                final int minj = roiPoints[j].getZint() - zRadius;
                final int maxj = roiPoints[j].getZint() + zRadius;
                if (mini <= minj) {
                    overlap = (maxi >= minj);
                } else {
                    overlap = (maxj >= mini);
                }
            }
            if (overlap) {
                ImageJUtils.log("Warning: Spot %d region overlaps with spot %d, ignoring both", i + 1, j + 1);
                bad[i] = bad[j] = true;
                break;
            }
        }
    }
    return bad;
}

Example 4 with ImageExtractor

use of uk.ac.sussex.gdsc.core.utils.ImageExtractor in project GDSC-SMLM by aherbert.

the class PsfCreator method runUsingFitting.

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

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

Example 5 with ImageExtractor

use of uk.ac.sussex.gdsc.core.utils.ImageExtractor in project GDSC-SMLM by aherbert.

the class PsfCreator method relocateCentres.

/**
 * Extract the stack for each centre and try and guess the z-centre based on the type of PSF.
 * Relocate the XY centre using the centre-of-mass around the pixels close to the z-centre slice.
 *
 * @param image the image
 * @param centres the centres
 * @return the new centres
 */
private BasePoint[] relocateCentres(float[][] image, BasePoint[] centres) {
    final int w = imp.getWidth();
    final int h = imp.getHeight();
    // Just for getting the bounds
    final ImageExtractor ie = ImageExtractor.wrap(image[0], w, h);
    // This can be reused as a buffer
    final float[][] psf = new float[image.length][];
    for (int i = 0; i < centres.length; i++) {
        // Extract stack
        final int x = centres[i].getXint();
        final int y = centres[i].getYint();
        final Rectangle bounds = ie.getBoxRegionBounds(x, y, boxRadius);
        for (int z = 0; z < image.length; z++) {
            psf[z] = ImageJImageConverter.getData(image[z], w, h, bounds, psf[z]);
        }
        if (settings.getInteractiveMode()) {
            final Overlay o = new Overlay();
            final Roi roi = new Roi(bounds);
            o.add(roi);
            imp.setOverlay(o);
            imp.updateAndDraw();
        }
        // Create a PSF and select the z-centre
        zSelector.setPsf(new ExtractedPsf(psf, bounds.width, bounds.height));
        zSelector.analyse();
        zSelector.guessZCentre();
        if (settings.getInteractiveMode()) {
            // Ask user for z-centre confirmation
            final double dz = zSelector.run("Confirm PSF Z-centre", true, false, false, Integer.toString(i + 1));
            if (dz == -1) {
                resetImp();
                return null;
            }
            if (dz == -2) {
                // Exclude this PSF
                centres[i] = null;
                continue;
            }
        }
        zCentre = zSelector.getCentreSlice() - 1;
        // Subtract background
        final float adjust = -zSelector.background;
        for (int z = 0; z < image.length; z++) {
            SimpleArrayUtils.add(psf[z], adjust);
        }
        // Update centre
        final double[] com = getCentreOfMassXy(psf, bounds.width, bounds.height, zCentre, settings.getComWindow(), getComXyBorder(bounds.width, bounds.height));
        float dx = (float) (com[0] + bounds.x - centres[i].getX());
        float dy = (float) (com[1] + bounds.y - centres[i].getY());
        float dz = (float) (zSelector.zCentre - centres[i].getZ());
        BasePoint newCentre = centres[i].shift(dx, dy, dz);
        if (settings.getSubPixelPrecision() > 0) {
            newCentre = new BasePoint((float) MathUtils.round(newCentre.getX(), settings.getSubPixelPrecision()), (float) MathUtils.round(newCentre.getY(), settings.getSubPixelPrecision()), (float) MathUtils.round(newCentre.getZ(), settings.getSubPixelPrecision()));
            dx = newCentre.getX() - centres[i].getX();
            dy = newCentre.getY() - centres[i].getY();
            dz = newCentre.getZ() - centres[i].getZ();
        }
        ImageJUtils.log("Centre %d : %s,%s,%s updated by %s,%s,%s", i + 1, rounder.toString(centres[i].getX()), rounder.toString(centres[i].getY()), rounder.toString(centres[i].getZ()), rounder.toString(dx), rounder.toString(dy), rounder.toString(dz));
        centres[i] = newCentre;
    }
    if (settings.getInteractiveMode()) {
        imp.setOverlay(null);
        // Check if any centres were excluded
        int size = 0;
        for (int i = 0; i < centres.length; i++) {
            if (centres[i] == null) {
                continue;
            }
            centres[size++] = centres[i];
        }
        if (size == 0) {
            resetImp();
            IJ.error(TITLE, "No remaining PSF centres");
            return null;
        }
        if (size < centres.length) {
            centres = Arrays.copyOf(centres, size);
        }
    }
    return centres;
}