Using Cellpose mask in QuPath- cell detection

Hi all,
I want to segment and later on create classification by QuPath.
I checked a few methods, and for my data, it looks like Cellpose doing the best (without specific training).

I used @Mike_Nelson code to adapt the ROIs\ annotation back to QuPath (Thank you!), but I want to use the feature of ‘artificial’ cytoplasm and dilate the nuclei segmentation as in “cell detection” (since a lot of my staining is also in the cytoplasm).

Do you know how can I do it? (I tried to adapt from the StarDist code, but it different)

thank you!

@petebankhead

Mike_Nelson code:

Not sure if this is what you mean about adapting the StarDist code, but @smcardle wrote the following: Inconsistent behavior in QuPath's Stardist implementation - #8 by smcardle
It was intended for StarDist nuclei, but I think it accomplishes the same general idea. I tried to run it quickly on a test sample, but ran into some errors - though it might work for you. It seems I could not get the fpDetection created for my one channel grayscale image (single frame, so no tiling), probably due to it being within:

    if (ip instanceof ColorProcessor && stains != null) {

There might have been some issue for CellPose with casting ShapeRois to PolygonRois.

You may want to remove the StarDist section around line 66, since the script does run StarDist. I have some other meetings so I will take a look at trying to get this to work with my sample image later, but that may not apply to your images, in whatever format they are.

All that said, if I edited out the fpDetection lines and do not force the PolygonRoi conversion, I can get the expansion as shown:
image
Big caveat - no measurements. You may get the original script to run on a brightfield image, if that is what you are analyzing.

Butchered version of script
/****
 * 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 = 0 //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()


/**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.
print roisNuclei.size()
//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
        //print roisNuclei.get(i).getClass().toString()
        //if (roisNuclei.get(i).getClass().toString() != "class ij.gui.ShapeROI"){return}
        /*PolygonRoi*/ rN = roisNuclei.get(i)
        //print "test10"
        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);
}
3 Likes

Come to think of it, if you want measurements, you should be able to add those using this post, and possibly the following posts:

If you have a single channel image like I did, you may need to change the channel name with a setChannelNames() line.

Thank you @Mike_Nelson for your help.

I apologize I didn’t give you enough details,
My images are fluorescent (9 channel) including DAPI, I attached now tile of DAPI channel and the segmentation.

I tried @smcardle code you attached here (and try to change the StarDist parts) it doesn’t work for me, and she also wrote in the description, this code is for H&E and DAB staining and not for fluorescent.

DAPI with cellpose segmentation.tif (3.6 MB) DAPI.tif (3.4 MB)

Not sure what to say, your sample image worked with the edited version of the code for me, though I did need to edit the channel name and add measurements separately as mentioned.

what??? :woozy_face:
Can you send me please the code you used? Maybe I changed it not in the way I should?
Many thanks

The code was posted in the original post in the last text block. It was minimized since it is a rather long script.

I’m sorry I don’t understand where to change the channel name, and without it, the code run and I get the message “Cell expansion complete” without any change. :worried:

Sorry, the only other thing I can think of is if you did not add annotations as indicated in Sara’s script. Without annotations, I would certainly expect no cell expansion.

The function setChannelNames only requires channel names. See here for an example: Multiplexed analysis — QuPath 0.2.3 documentation

I think @Mike_Nelson is correct about needing an annotation. Try putting in a full image annotation (Ctrl + Shift + A or Objects > Annotations > Create Full Image Annotation) before running the script. Please remember that the cells at the boundaries between different images are going to have really odd and unusable measurements and shapes.

To get the measurements in one script, you could try copy-pasting the code from lines 240-275 here after the If statement in Mike’s edit script: if (ip instanceof ColorProcessor && stains != null) { ... } . I don’t have time to test that myself, but I think if you get the correct fpDetection image, the measurements should work from there. Probably?

2 Likes

Finally, understand what was my problem!

The workflow process:

  1. Adapt CP PNG mask with @Mike_Nelson code

Blockquote
QuPath - transferring objects from one image to the other - #3 by Mike_Nelson

  1. The second step, is like @smcardle wrote- to add the ‘full image annotation’

    • my problem was that all the cells\ ROIs were added as a list of annotation so to solve that-
  2. In this step, I organized the Hierarchy and selected all the ROIs\ nuclei and "insert (them) into hierarchy’ (ctrl+shift+I).

  3. Then I run the code you added here, with removing the path and the import of StarDist, and with the change of fpDetection to the IF Ch I want to detect (instead of what was written before) -
    fpDetection = setChannelNames('DAPI');

  4. The next step was the 'measurements ’ part, I just run it in separate code as Mike wrote (the integration of the codes as Sara suggested didn’t work for me, can’t think now what was wrong with that).

  5. And continue to classify objects! :slight_smile:
    (that seems to be working)

Thank you both for all your help :pray:

BTW, I worked with ver 0.2.3

1 Like

I would recommend importing the initial objects, if they are cells, as detections instead of annotations. QuPath will tend to run better, and it will make exporting easier as you will have a clear delineation between detections as cells and annotations as the surrounding whole image annotation.

1 Like

Probably you right.
I add this line:

def detections = getDetectionObjects()

Now all the objects\cells define as “Cell” in both- ‘Hierarchy’ and ‘Annotation’, this is how it looks in the end, very nice :grinning::
@Research_Associate, do you think it’s ok now?

I do not think you would want them to show up in both the annotations tab and the Hierarchy (which indicates duplicates), but only in the Hierarchy. Either way, as long as it is working for you.

For now, it looks ok, but I think you right.
I don’t know how to delete them from only one of the options (Hierarchy or Annotations).

You should not need to delete any annotations, better to not create them to begin with. The original script, from what I can tell, only created detection objects, so it looks like you changed something in order to get annotation objects.