Inconsistent behavior in QuPath's Stardist implementation

Hi @petebankhead,
First, I want to say thank you for making a Stardist plugin for QuPath- it’s great! But, I noticed some inconsistent behavior in how it deals with region boundaries. This was all done in m11 with the CPU Tensorflow build (thanks @Research_Associate!)

If you detect nuclei only, it will show the complete nuclei, even if part of the nucleus is outside of an annotation border:

If you perform cell expansion, it cuts off both the nuclei and the cell membrane at the annotation boundary:

If you add ignoreCellOverlaps(true), it cuts off the membrane boundary, but not the nuclei:

I’m not sure what the best behavior “should” be, but that 3rd one is wrong. In my particular case with these odd shaped annotations from the pixel classifier, I’d prefer it to include the entire cell, even if some part of the cell is outside of the annotation. But, I can easily see the opposite argument.

Also, I’m seeing the artifact you mentioned in the readthedocs page on Stardist, where there are issues with pieces missing from some shapes:
image
Is there any way to use the nuclei from Stardist and then the original pixel-based cell expansion? This would let us do exactly what we need- grab the full nucleus, even the parts slightly outside the annotation of interest, and then expand without overlap. :slight_smile:

Thanks!

1 Like

Doh, will get you the full 0.2.0 shortly.

Thanks for the feedback, by way of explanation:

  • QuPath doesn’t crop the StarDist nuclei to any annotation in general (since the StarDist output is in the form of vertices, it can predict a vertex beyond the region of interest)
  • If there is cell expansion, QuPath crops the expansion… and then crops the nucleus to make sure it fits inside the cell

These two rules produce weird behavior in the third case… although I suspect that ignoreCellOverlaps(true) isn’t generally a very useful option for anything other than debugging.

The opposite argument is more likely to win, since it is a more frequent request to prevent QuPath from expanding beyond the tissue :slight_smile:

However, probably what is required is another option to either crop or not.

You can adjust the simplify option in the builder; this can have an impact but I have not been able to find a way to eliminate these completely.

'Fraid not… I plan to return to the StarDist code at some point, but I’m not sure when I’ll get to it.

3 Likes

Here’s a pull request that aims to resolve the inconsistent behavior: https://github.com/qupath/qupath/pull/539
It also adds a new constrainToParent(false) option (the default is true, if unspecified).

No progress on the artifacts – how frequently do you see them? I spent quite a lot of time on this at the beginning and managed to reduce the frequency considerably, but I couldn’t eliminate them completely… there is rather a lot of geometry manipulation going on. However, further improvement is difficult because I can replicate the problem only rarely now.

I don’t really see an easy solution, but will look again when I have time. I would prefer to avoid needing to go a pixel-based route if possible, since this opens up a world of painful resolution / rasterization / overlap resolution problems… and ultimately decoupling the cell expansion from the nucleus detection feels like the right way to go, to make it easier to add more cell detection methods in the future.

1 Like

This might seem gimicky, but could you store the nucleus created by stardist, and then replace the “expanded nucleus” with the original after creating the expansion?

I have not looked at the code in detail, so this is purely an off the cuff thought. I am not sure how it would deal with the inner point being cytoplasmic. Hmm. Unless you could remove any ROI points from the cytoplasm that were contained within the nucleus, then rebuild both ROIs.
That might be something we could do in a script though.

If I use a pixel size similar to the image pixel size, it’s maybe ~1% of cells. If I use a larger pixel size, it’s more often than that (~10%). I’m happy to send you an example image that fails frequently :stuck_out_tongue: .

1 Like

That could be handy – or even better if you can find a way to get it to fail consistently with anything I already have.

I started writing something to use the original cell expansion yesterday, before you posted the new update with the constrainToParent() option. I wasn’t able to drop it from my brain just because you solved it already more elegantly

So, in case it helps anyone, here’s a script that takes the nuclei created by Stardist and then applies the original cell expansion. These cells are not constrained to the parent boundary (though it looks like that is now included in the QuPath code). It also avoids the artifact where chunks of cells are occasionally missing. Important notes: it is only written for brightfield, and it assumes you have 1 or more SMALL parent annotations. This doesn’t do any of the tiling that QuPath uses to keep memory issues under control, so it will not work with large annotations.

It’s mostly taken from the QuPath Source Code but greatly simplified. Thanks to @Mike_Nelson for lots of help!


/****
 * Script to use original QuPath cell expansion algorithm on Stardist nuclei
 * This prevents the artifact that sometimes occurs with the cellExpansion parameter in the Stardist script.
 * It also allows cells boundaries to go outside the annotation.
 * Assumes you have a few small annotations in your image that you want to use to detect nuclei.
 *
 * Grabbed from the QuPath source code, but greatly simplified:
 * -Cell boundary created at full resolution. MUST HAVE SMALL ANNOTATIONS because it will load the entire annotation region at once.
 * -Only works for H&E or H-DAB images, not fluorescence.
 *
 * Written by Sara McArdle of the La Jolla Institute, 2020, with help from the incredible Michael Nelson.
 * QuPath 0.2.0
 ****/


/**imports**/
import ij.gui.PolygonRoi
import ij.measure.Calibration
import ij.plugin.filter.EDM
import ij.process.ColorProcessor
import ij.process.FloatPolygon
import ij.process.FloatProcessor
import ij.process.ImageProcessor
import ij.process.ShortProcessor
import qupath.imagej.processing.RoiLabeling
import qupath.imagej.processing.Watershed
import qupath.imagej.tools.IJTools
import qupath.imagej.tools.PixelImageIJ
import qupath.imagej.detect.cells.ObjectMeasurements
import qupath.lib.analysis.images.SimpleImage
import qupath.lib.analysis.stats.RunningStatistics
import qupath.lib.analysis.stats.StatisticsHelper
import qupath.lib.color.ColorDeconvolutionStains
import qupath.lib.color.StainVector
import qupath.lib.gui.dialogs.Dialogs
import qupath.lib.images.PathImage
import qupath.lib.images.servers.PixelCalibration
import qupath.lib.images.servers.ServerTools
import qupath.lib.objects.PathObjectTools
import qupath.lib.regions.ImagePlane
import qupath.lib.regions.RegionRequest
import qupath.lib.roi.ShapeSimplifier

import static qupath.lib.gui.scripting.QPEx.*
import qupath.lib.objects.PathObject
import qupath.lib.objects.PathObjects
import ij.process.ByteProcessor;
import ij.gui.Roi;
import qupath.lib.roi.PolygonROI;

import qupath.tensorflow.stardist.StarDist2D

/**Things to Set**/
def pathModel = "path\\to\\stardist\\model"
double cellExpansionMicrons=3
makeMeasurements=true //only tested with true
smoothBoundaries=true //only tested with true
//need a region slightly larger than actual annotation because the cell boundaries can be outside. There's probably a better way to do this.
int delta = 80 //must be even or ImageJ gets sad. In pixels.
//Don't forget about StarDist settings

/**Adjust Stardist settings to your liking**/
def imageData = getCurrentImageData()
PixelCalibration Qcal=imageData.getServer().getPixelCalibration()
def pixelSize=Qcal.getAveragedPixelSize()

def stardist = StarDist2D.builder(pathModel)
        .threshold(0.15)              // Prediction threshold
        .normalizePercentiles(1, 99) // Percentile normalization
        .pixelSize(pixelSize*2)              // Resolution for detection
        .includeProbability(true)    // Add probability as a measurement (enables later filtering)
        .tileSize(2048)              // Specify width & height of the tile used for prediction
        .build()

// Run detection for the selected objects
def annotations = getAnnotationObjects()
if (annotations.isEmpty()) {
    Dialogs.showErrorMessage("StarDist", "Please select a parent object!")
    return
}
stardist.detectObjects(imageData, annotations)
println 'StarDist Done!'

/**Begin Cell Expansion **/
//https://github.com/qupath/qupath/blob/015765039dd22373ea4a40ba1efb9ac0c1c30e55/qupath-core-processing/src/main/java/qupath/imagej/detect/cells/WatershedCellDetection.java#L844

double cellExpansion = cellExpansionMicrons/pixelSize
imageData=getCurrentImageData()
ImageServer<BufferedImage> server = imageData.getServer()
double downsample = ServerTools.getDownsampleFactor(server, pixelSize);
double downsampleSqrt = Math.sqrt(downsample);
ColorDeconvolutionStains stains = imageData.getColorDeconvolutionStains();

//loop over annotations
def annot=getAnnotationObjects()
annot.each { a ->
    annotROI = a.getROI()


//Need ALL the coordinates
    int width = (int) (annotROI.getBoundsWidth() + delta)
    int height = (int) (annotROI.getBoundsHeight() + delta)
    int left = (int) (annotROI.getBoundsX().toInteger() - delta / 2)
    int top = (int) (annotROI.getBoundsY().toInteger() - delta / 2)
    ImagePlane plane = ImagePlane.getPlane(annotROI)

//black image to put nuclei on top of
    ByteProcessor bp = new ByteProcessor(width, height)
    bp.setValue(255)

//Actual image of the islet
    PathImage pathImage = IJTools.convertToImagePlus(server, RegionRequest.createInstance(server.getPath(), downsample, left, top, width, height));
    ip = pathImage.getImage().getProcessor()
    Calibration cal = pathImage.getImage().getCalibration();

//initializing the color channels
    Map<String, FloatProcessor> channels = new LinkedHashMap<>(); // Map of channels to measure for nuclei only, and their names
    Map<String, FloatProcessor> channelsCell = new LinkedHashMap<>(); // Map of channels to measure for cell/cytoplasm, and their names

//get the deconvolved stians
    if (ip instanceof ColorProcessor && stains != null) {
        FloatProcessor[] fps = IJTools.colorDeconvolve((ColorProcessor) ip, stains);
        for (int i = 0; i < 3; i++) {
            StainVector stain = stains.getStain(i + 1);
            if (!stain.isResidual()) {
                channels.put(stain.getName() + " OD", fps[i]);
                channelsCell.put(stain.getName() + " OD", fps[i]);
            }
        }

        fpDetection = (FloatProcessor) fps[0].duplicate();
        //switch to this to use Optical Density Sum
        //fpDetection = IJTools.convertToOpticalDensitySum((ColorProcessor)ip, stains.getMaxRed(), stains.getMaxGreen(), stains.getMaxBlue());
    }

//turn nuclei objects into a black and white image
    def nucleiObjs=a.getChildObjects().findAll{it.isDetection()}
    def nuclei = nucleiObjs.collect { it.getROI() }
    List<PolygonRoi> roisNuclei = new ArrayList<>()
    nuclei.each {
        //def shift = it.translate(width/2-centerX, height/2-centerY)
        //Roi roiIJ=IJTools.convertToIJRoi(shift,0,0,downsample)
        Roi roiIJ = IJTools.convertToIJRoi(it, cal, downsample)
        roisNuclei.add(roiIJ)
        bp.fill(roiIJ);
    }
//bp is now a black and white image with nuclei.

//create label image of nuclei (1-# nuclei, 0 = background).
    ShortProcessor ipLabels = new ShortProcessor(width, height)
    RoiLabeling.labelROIs(ipLabels, roisNuclei);
//ipLabels is now a label image of the nuclei .

//calculate intensity statistics for nuclei
    Map<String, List<RunningStatistics>> statsMap = new LinkedHashMap<>();
    if (makeMeasurements) {
        SimpleImage imgLabels = new PixelImageIJ(ipLabels);
        for (String key : channels.keySet()) {
            List<RunningStatistics> statsList = StatisticsHelper.createRunningStatisticsList(roisNuclei.size());
            StatisticsHelper.computeRunningStatistics(new PixelImageIJ(channels.get(key)), imgLabels, statsList);
            statsMap.put(key, statsList);
        }
    }

/**Now we are ready to find the cell boundaries**/
//Euclidean Distance Map of distance to nuclei.
    FloatProcessor fpEDM = new EDM().makeFloatEDM(bp, (byte) 255, false)
    fpEDM.multiply(-1);
//fpEDM = 0 at all nuclei pixels, increasingly negative as you get further away

//watershed segmentation on binary nuclei image, up to the distance given by cellExpansion
    ImageProcessor ipLabelsCells = ipLabels.duplicate()
    double cellExpansionThreshold = -cellExpansion;
    Watershed.doWatershed(fpEDM, ipLabelsCells, cellExpansionThreshold, false);
//ipLabelCells now is a label image, like ipLabels, but larger for the cells. Already watershed segmented.

//turn label image into rois (imageJ rois, not yet qupath rois)
    PolygonRoi[] roisCells = RoiLabeling.labelsToFilledROIs(ipLabelsCells, roisNuclei.size());

//calculate intensity measurements on cells
    Map<String, List<RunningStatistics>> statsMapCell = new LinkedHashMap<>();
    if (makeMeasurements) {
        for (String key : channelsCell.keySet()) {
            List<RunningStatistics> statsList = StatisticsHelper.createRunningStatisticsList(roisNuclei.size());
            StatisticsHelper.computeRunningStatistics(new PixelImageIJ(channelsCell.get(key)), new PixelImageIJ(ipLabelsCells), statsList);
            statsMapCell.put(key, statsList);
        }
    }

// Create labelled image for cytoplasm, i.e. remove all nucleus pixels
    for (int i = 0; i < ipLabels.getWidth() * ipLabels.getHeight(); i++) {
        if (ipLabels.getf(i) != 0)
            ipLabelsCells.setf(i, 0f);
    }

//calculate intensity measurements on cytoplasm
    Map<String, List<RunningStatistics>> statsMapCytoplasm = new LinkedHashMap<>();
    if (makeMeasurements) {
        for (String key : channelsCell.keySet()) {
            List<RunningStatistics> statsList = StatisticsHelper.createRunningStatisticsList(roisNuclei.size());
            StatisticsHelper.computeRunningStatistics(new PixelImageIJ(channelsCell.get(key)), new PixelImageIJ(ipLabelsCells), statsList);
            statsMapCytoplasm.put(key, statsList);
        }
    }

//turns the ImageJ rois into QuPath ROIs
    List<PathObject> pathObjects = new ArrayList<>() //initialize
    for (int i = 0; i < roisCells.length; i++) {
        //actually add the measurements to the objects
        PolygonRoi rN = roisNuclei.get(i)
        PolygonRoi r = roisCells[i]
        def nucleusObj = nucleiObjs[i]

        if (r == null)
            continue;

        if (makeMeasurements) {

            measurementList = nucleusObj.getMeasurementList()
            ObjectMeasurements.addShapeStatistics(measurementList, rN, pathImage.getImage().getProcessor(), cal, "Nucleus: ");

            for (String key : channels.keySet()) {
                List<RunningStatistics> statsList = statsMap.get(key);
                RunningStatistics stats = statsList.get(i);
                measurementList.addMeasurement("Nucleus: " + key + " mean", stats.getMean());
                measurementList.addMeasurement("Nucleus: " + key + " sum", stats.getSum());
                measurementList.addMeasurement("Nucleus: " + key + " std dev", stats.getStdDev());
                measurementList.addMeasurement("Nucleus: " + key + " max", stats.getMax());
                measurementList.addMeasurement("Nucleus: " + key + " min", stats.getMin());
                measurementList.addMeasurement("Nucleus: " + key + " range", stats.getRange());
            }
        }

        //smooth boundaries of cell polygons
        if (smoothBoundaries) {
            r = new PolygonRoi(r.getInterpolatedPolygon(1, false), Roi.POLYGON);
            r = smoothPolygonRoi(r);
            r = new PolygonRoi(r.getInterpolatedPolygon(Math.min(2, r.getNCoordinates().toDouble() * 0.1), false), Roi.POLYGON);
        }

        PolygonROI pathROI = IJTools.convertToPolygonROI(r, cal, downsample, plane);
        if (smoothBoundaries)
            pathROI = ShapeSimplifier.simplifyPolygon(pathROI, downsampleSqrt / 2.0);

        if (makeMeasurements) {
            ObjectMeasurements.addShapeStatistics(measurementList, r, fpDetection, pathImage.getImage().getCalibration(), "Cell: ");

            // Add cell measurements
            for (String key : channelsCell.keySet()) {
                if (statsMapCell.containsKey(key)) {
                    RunningStatistics stats = statsMapCell.get(key).get(i);
                    measurementList.addMeasurement("Cell: " + key + " mean", stats.getMean());
                    measurementList.addMeasurement("Cell: " + key + " std dev", stats.getStdDev());
                    measurementList.addMeasurement("Cell: " + key + " max", stats.getMax());
                    measurementList.addMeasurement("Cell: " + key + " min", stats.getMin());
                }
            }

            // Add cytoplasm measurements
            for (String key : channelsCell.keySet()) {
                if (statsMapCytoplasm.containsKey(key)) {
                    RunningStatistics stats = statsMapCytoplasm.get(key).get(i);
                    measurementList.addMeasurement("Cytoplasm: " + key + " mean", stats.getMean());
                    measurementList.addMeasurement("Cytoplasm: " + key + " std dev", stats.getStdDev());
                    measurementList.addMeasurement("Cytoplasm: " + key + " max", stats.getMax());
                    measurementList.addMeasurement("Cytoplasm: " + key + " min", stats.getMin());
                }
            }

            // Add nucleus area ratio, if available
            if (nucleusObj != null && nucleusObj.getROI().isArea()) {
                double nucleusArea = nucleusObj.getROI().getArea();
                double cellArea = pathROI.getArea();
                measurementList.addMeasurement("Nucleus/Cell area ratio", Math.min(nucleusArea / cellArea, 1.0));
            }
        }

        // Create & store the cell object
        PathObject pathObject = PathObjects.createCellObject(pathROI, nucleusObj == null ? null : nucleusObj.getROI(), null, measurementList);
        pathObjects.add(pathObject);
    }

    for (PathObject pathObject : pathObjects)
        pathObject.getMeasurementList().close();

// Sometimes smoothing can cause nuclei of cell boundaries to be removed - in this case,
// filter out the invalid ROIs now
    int sizeBefore = pathObjects.size();
    pathObjects.removeIf(p -> PathObjectTools.getROI(p, false).isEmpty() ||
            PathObjectTools.getROI(p, true).isEmpty());
    int sizeAfter = pathObjects.size();
    if (sizeBefore != sizeAfter) {
        logger.debug("Filtered out {} invalid cells (empty ROIs)", sizeBefore - sizeAfter);
    }

//add the objects to the hierarchy and remove the (now duplicated) nuclei
    addObjects(pathObjects)
    removeObjects(nucleiObjs, true)
}
resolveHierarchy()
print('Cell Expansion Complete')

private static PolygonRoi smoothPolygonRoi(PolygonRoi r) {
    FloatPolygon poly = r.getFloatPolygon();
    FloatPolygon poly2 = new FloatPolygon();
    int nPoints = poly.npoints;
    for (int i = 0; i < nPoints; i += 2) {
        int iMinus = (i + nPoints - 1) % nPoints;
        int iPlus = (i + 1) % nPoints;
        poly2.addPoint((poly.xpoints[iMinus] + poly.xpoints[iPlus] + poly.xpoints[i])/3,
                (poly.ypoints[iMinus] + poly.ypoints[iPlus] + poly.ypoints[i])/3);
    }
//			return new PolygonRoi(poly2, r.getType());
    return new PolygonRoi(poly2, Roi.POLYGON);
}
2 Likes