Examples with MaximumFinder ij.plugin.filter.MaximumFinder used on opensource projects

Example 1 with MaximumFinder

use of ij.plugin.filter.MaximumFinder in project qupath by qupath.

the class SubcellularDetection method processObject.

/**
 * Initial version of subcellular detection processing.
 *
 * @param pathObject
 * @param params
 * @param imageWrapper
 * @return
 * @throws InterruptedException
 * @throws IOException
 */
static boolean processObject(final PathObject pathObject, final ParameterList params, final ImageWrapper imageWrapper) throws InterruptedException, IOException {
    // Get the base classification for the object as it currently stands
    PathClass baseClass = PathClassTools.getNonIntensityAncestorClass(pathObject.getPathClass());
    // Variable to hold estimated spot count
    double estimatedSpots;
    // We assume that after this processing, any previous sub-cellular objects should be removed
    pathObject.clearPathObjects();
    // Ensure we have no existing subcellular detection measurements - if we do, remove them
    String[] existingMeasurements = pathObject.getMeasurementList().getMeasurementNames().stream().filter(n -> n.startsWith("Subcellular:")).toArray(n -> new String[n]);
    if (existingMeasurements.length > 0) {
        pathObject.getMeasurementList().removeMeasurements(existingMeasurements);
        pathObject.getMeasurementList().close();
    }
    // // If we're part of a TMA core, request the whole core...
    // if (pathObject.getParent() instanceof TMACoreObject && pathObject.getParent().hasROI()) {
    // regionStore.getImage(server, RegionRequest.createInstance(server.getPath(), 1, pathObject.getParent().getROI()), 25, true);
    // }
    ROI pathROI = pathObject.getROI();
    if (pathROI == null || pathROI.isEmpty())
        return false;
    // double downsample = 0.5;
    double downsample = 1;
    // Determine spot size
    ImageServer<BufferedImage> server = imageWrapper.getServer();
    PixelCalibration cal = server.getPixelCalibration();
    double minSpotArea, maxSpotArea, singleSpotArea;
    double pixelWidth, pixelHeight;
    if (cal.hasPixelSizeMicrons()) {
        double spotSizeMicrons = params.getDoubleParameterValue("spotSizeMicrons");
        double minSpotSizeMicrons = params.getDoubleParameterValue("minSpotSizeMicrons");
        double maxSpotSizeMicrons = params.getDoubleParameterValue("maxSpotSizeMicrons");
        pixelWidth = cal.getPixelWidthMicrons() * downsample;
        pixelHeight = cal.getPixelHeightMicrons() * downsample;
        singleSpotArea = spotSizeMicrons / (pixelWidth * pixelHeight);
        minSpotArea = minSpotSizeMicrons / (pixelWidth * pixelHeight);
        maxSpotArea = maxSpotSizeMicrons / (pixelWidth * pixelHeight);
    } else {
        singleSpotArea = params.getDoubleParameterValue("spotSizePixels");
        minSpotArea = params.getDoubleParameterValue("minSpotSizePixels");
        maxSpotArea = params.getDoubleParameterValue("maxSpotSizePixels");
        pixelWidth = downsample;
        pixelHeight = downsample;
    }
    boolean includeClusters = Boolean.TRUE.equals(params.getBooleanParameterValue("includeClusters"));
    boolean doSmoothing = Boolean.TRUE.equals(params.getBooleanParameterValue("doSmoothing"));
    boolean splitByIntensity = Boolean.TRUE.equals(params.getBooleanParameterValue("splitByIntensity"));
    boolean splitByShape = Boolean.TRUE.equals(params.getBooleanParameterValue("splitByShape"));
    // Get region to request - give a pixel as border
    int xStart = (int) Math.max(0, pathROI.getBoundsX() - 1);
    int yStart = (int) Math.max(0, pathROI.getBoundsY() - 1);
    int width = (int) Math.min(server.getWidth() - 1, pathROI.getBoundsX() + pathROI.getBoundsWidth() + 1.5) - xStart;
    int height = (int) Math.min(server.getHeight() - 1, pathROI.getBoundsY() + pathROI.getBoundsHeight() + 1.5) - yStart;
    if (width <= 0 || height <= 0) {
        logger.error("Negative ROI size for {}", pathROI);
        pathObject.setPathClass(baseClass);
        return false;
    }
    int z = pathROI.getZ();
    int t = pathROI.getT();
    // Don't associate with channel
    int c = -1;
    RegionRequest region = RegionRequest.createInstance(server.getPath(), 1.0, xStart, yStart, width, height, z, t);
    // Mask to indicate pixels within the cell
    byte[] cellMask = null;
    for (String channelName : imageWrapper.getChannelNames(true, true)) {
        double detectionThreshold = params.getDoubleParameterValue("detection[" + channelName + "]");
        if (Double.isNaN(detectionThreshold) || detectionThreshold < 0)
            continue;
        // // TODO: Consider whether to use channel numbers for non-brightfield images
        // if (!imageWrapper.imageData.isBrightfield())
        // c++;
        SimpleImage img = imageWrapper.getRegion(region, channelName);
        // Get an ImageJ-friendly calibration for ROI conversion
        Calibration calIJ = new Calibration();
        calIJ.xOrigin = -xStart / downsample;
        calIJ.yOrigin = -yStart / downsample;
        // Create a cell mask
        if (cellMask == null) {
            BufferedImage imgMask = new BufferedImage(img.getWidth(), img.getHeight(), BufferedImage.TYPE_BYTE_GRAY);
            Graphics2D g2d = imgMask.createGraphics();
            if (downsample != 1)
                g2d.scale(1.0 / downsample, 1.0 / downsample);
            g2d.translate(-xStart, -yStart);
            Shape shape = RoiTools.getShape(pathROI);
            g2d.setColor(Color.WHITE);
            g2d.fill(shape);
            g2d.dispose();
            cellMask = (byte[]) ((DataBufferByte) imgMask.getRaster().getDataBuffer()).getData(0);
        }
        // Get a buffer containing the image pixels
        int w = img.getWidth();
        int h = img.getHeight();
        // Identify (& try to separate) spots
        // Mask out non-cell areas as we go
        FloatProcessor fpDetection = new FloatProcessor(w, h);
        if (doSmoothing) {
            for (int i = 0; i < w * h; i++) fpDetection.setf(i, img.getValue(i % w, i / w));
            fpDetection.smooth();
            for (int i = 0; i < w * h; i++) {
                if (cellMask[i] == (byte) 0)
                    fpDetection.setf(i, 0f);
            }
        } else {
            for (int i = 0; i < w * h; i++) {
                if (cellMask[i] == (byte) 0)
                    fpDetection.setf(i, 0f);
                else
                    fpDetection.setf(i, img.getValue(i % w, i / w));
            }
        }
        ByteProcessor bpSpots;
        if (splitByIntensity)
            bpSpots = new MaximumFinder().findMaxima(fpDetection, detectionThreshold / 10.0, detectionThreshold, MaximumFinder.SEGMENTED, false, false);
        else
            bpSpots = SimpleThresholding.thresholdAboveEquals(fpDetection, (float) detectionThreshold);
        if (splitByShape) {
            new EDM().toWatershed(bpSpots);
        }
        // Loop through spot ROIs & make a decision
        bpSpots.setThreshold(1, ImageProcessor.NO_THRESHOLD, ImageProcessor.NO_LUT_UPDATE);
        List<PolygonRoi> possibleSpotRois = RoiLabeling.getFilledPolygonROIs(bpSpots, Wand.FOUR_CONNECTED);
        List<PathObject> spotObjects = new ArrayList<>();
        List<PathObject> clusterObjects = new ArrayList<>();
        estimatedSpots = 0;
        for (PolygonRoi spotRoi : possibleSpotRois) {
            fpDetection.setRoi(spotRoi);
            ImageStatistics stats = fpDetection.getStatistics();
            // In v0.2
            // ImagePlane plane = ImagePlane.getPlaneWithChannel(spotRoi.getCPosition(), spotRoi.getZPosition(), spotRoi.getTPosition());
            // In v0.3
            ImagePlane plane = ImagePlane.getPlaneWithChannel(c, z, t);
            PathObject spotOrCluster = null;
            if (stats.pixelCount >= minSpotArea && stats.pixelCount <= maxSpotArea) {
                ROI roi = IJTools.convertToROI(spotRoi, calIJ, downsample, plane);
                // cluster = new SubcellularObject(roi, 1);
                spotOrCluster = createSubcellularObject(roi, 1);
                estimatedSpots += 1;
            } else if (includeClusters && stats.pixelCount >= minSpotArea) {
                // Add a cluster
                ROI roi = IJTools.convertToROI(spotRoi, calIJ, downsample, plane);
                double nSpots = stats.pixelCount / singleSpotArea;
                estimatedSpots += nSpots;
                // cluster = new SubcellularObject(roi, nSpots);
                spotOrCluster = createSubcellularObject(roi, nSpots);
            }
            if (spotOrCluster != null) {
                boolean isCluster = spotOrCluster.getMeasurementList().getMeasurementValue("Num spots") > 1;
                int rgb = imageWrapper.getChannelColor(channelName);
                rgb = isCluster ? ColorTools.makeScaledRGB(rgb, 0.5) : ColorTools.makeScaledRGB(rgb, 1.5);
                PathClass pathClass = PathClassFactory.getDerivedPathClass(spotOrCluster.getPathClass(), channelName + " object", rgb);
                spotOrCluster.setPathClass(pathClass);
                spotOrCluster.getMeasurementList().putMeasurement("Subcellular cluster: " + channelName + ": Area", stats.pixelCount * pixelWidth * pixelHeight);
                spotOrCluster.getMeasurementList().putMeasurement("Subcellular cluster: " + channelName + ": Mean channel intensity", stats.mean);
                // cluster.getMeasurementList().putMeasurement("Subcellular cluster: " + channelName +  ": Max channel intensity", stats.max);
                spotOrCluster.getMeasurementList().close();
                if (isCluster)
                    clusterObjects.add(spotOrCluster);
                else
                    spotObjects.add(spotOrCluster);
            }
        }
        // Add measurements
        MeasurementList measurementList = pathObject.getMeasurementList();
        measurementList.putMeasurement("Subcellular: " + channelName + ": Num spots estimated", estimatedSpots);
        measurementList.putMeasurement("Subcellular: " + channelName + ": Num single spots", spotObjects.size());
        measurementList.putMeasurement("Subcellular: " + channelName + ": Num clusters", clusterObjects.size());
        // Add spots
        pathObject.addPathObjects(spotObjects);
        pathObject.addPathObjects(clusterObjects);
    }
    return true;
}

Example 2 with MaximumFinder

use of ij.plugin.filter.MaximumFinder in project qupath by qupath.

the class TMADearrayer method refineGridCoordinatesByShifting.

private static void refineGridCoordinatesByShifting(ByteProcessor bp, Polygon polyGrid, int nHorizontal, double coreDiameterPx) {
    // Create a binary image containing regions to consider - exclude very small & very large regions
    // Use watershed segmentation to split up round structures
    new EDM().toWatershed(bp);
    double estimatedArea = Math.PI * coreDiameterPx * coreDiameterPx * 0.25;
    FloatProcessor fpDensity = identifyGoodCores(bp, estimatedArea * 0.1, estimatedArea * 2.0, 0, false, null).convertToFloatProcessor();
    fpDensity.max(1);
    // new ImagePlus("Density-Orig", fpDensity.duplicate()).show();
    // Loop through the grid & see examples where the estimated point already falls within a likely core
    double minDiameter = coreDiameterPx * 0.7;
    Wand wand = new Wand(bp);
    boolean[] confirmed = new boolean[polyGrid.npoints];
    for (int i = 0; i < polyGrid.npoints; i++) {
        int x = polyGrid.xpoints[i];
        int y = polyGrid.ypoints[i];
        // If grid location already falls within a convincingly large ROI, use its bounding box centre
        // IJ.log(new Point(x, y).toString());
        boolean inside = fpDensity.getf(x, y) > 0;
        if (inside) {
            wand.autoOutline(x, y, 0.0, Wand.FOUR_CONNECTED);
            Roi roi = roiFromWand(wand);
            Rectangle bounds = roi.getBounds();
            if (bounds.width >= minDiameter && bounds.height >= minDiameter) {
                polyGrid.xpoints[i] = bounds.x + bounds.width / 2;
                polyGrid.ypoints[i] = bounds.y + bounds.height / 2;
                // Introduce a large density penalty to prevent missing cores from encroaching on this area
                // (Scale the ROI slightly first to make this more effective)
                roi = RoiScaler.scale(roi, 1, 1, true);
                fpDensity.setValue(-10000);
                fpDensity.fill(roi);
                // Mark ROI as confirmed
                confirmed[i] = true;
            }
        }
    }
    // Create (effectively) a Voronoi image based on the current centroid estimate
    ByteProcessor bpTest = new ByteProcessor(bp.getWidth(), bp.getHeight());
    for (int i = 0; i < polyGrid.npoints; i++) {
        if (!confirmed[i]) {
            int x = polyGrid.xpoints[i];
            int y = polyGrid.ypoints[i];
            bpTest.setf(x, y, 255);
        }
    }
    // Need to try/catch this, as might return null if thread is interrupted
    try {
        FloatProcessor fpEDM = new EDM().makeFloatEDM(bpTest, (byte) 255, false);
        // fpEDM.copyBits(bpTest, 0, 0, Blitter.ADD);
        for (int i = 0; i < fpEDM.getWidth() * fpEDM.getHeight(); i++) {
            fpEDM.setf(i, bpTest.getf(i) - fpEDM.getf(i));
        }
        ByteProcessor bpRegions = new MaximumFinder().findMaxima(fpEDM, 255, ImageProcessor.NO_THRESHOLD, MaximumFinder.SEGMENTED, false, false);
        // Penalise boundaries - when refining uncertain cores, we don't want to move into the territory of another core
        byte[] pxRegions = (byte[]) bpRegions.getPixels();
        float[] pxDensity = (float[]) fpDensity.getPixels();
        for (int i = 0; i < pxRegions.length; i++) {
            if (pxRegions[i] == 0)
                pxDensity[i] = -10000;
        }
    } catch (Exception e) {
        return;
    }
    // Apply a mean filter to determine local unassigned densities
    // new ImagePlus("Density_before", fpDensity.duplicate()).show();
    // long start = System.currentTimeMillis();
    // Note: this is another bottleneck... filter size can be large
    new RankFilters().rank(fpDensity, coreDiameterPx * 0.5, RankFilters.MEAN);
    // Find local maxima within each unassigned core region, with a preference towards the maximum closest to the original estimate
    for (int i = 0; i < polyGrid.npoints; i++) {
        if (!confirmed[i]) {
            int x = polyGrid.xpoints[i];
            int y = polyGrid.ypoints[i];
            OvalRoi roiRegion = new OvalRoi(x - coreDiameterPx * 0.5, y - coreDiameterPx * 0.5, coreDiameterPx, coreDiameterPx);
            Point maxPoint = findClosestMaximumInROI(fpDensity, roiRegion, new Point(x, y));
            if (maxPoint != null) {
                polyGrid.xpoints[i] = maxPoint.x;
                polyGrid.ypoints[i] = maxPoint.y;
                confirmed[i] = true;
            }
        }
    }
}