Whole slide image processing with QuPath & ImageJ

QuPath-Detecting regions with the help of ImageJ.groovy

/**
 * Demonstration of how to use ImageJ + QuPath to detect stained regions in
 * a brightfield image.
 *
 * @author Pete Bankhead
 */


import ij.CompositeImage
import ij.IJ
import ij.ImagePlus
import ij.ImageStack
import ij.plugin.filter.ThresholdToSelection
import ij.process.AutoThresholder
import ij.process.ColorProcessor
import ij.process.FloatProcessor
import ij.process.ImageStatistics
import qupath.imagej.color.ColorDeconvolutionIJ
import qupath.imagej.gui.IJExtension
import qupath.imagej.images.servers.ImagePlusServerBuilder
import qupath.imagej.objects.ROIConverterIJ
import qupath.lib.images.ImageData
import qupath.lib.objects.PathAnnotationObject
import qupath.lib.regions.RegionRequest
import qupath.lib.scripting.QPEx


// Just view the image - don't do further processing
// This is useful to help interactively explore potential parameters and processing steps
boolean viewOnly = false

// Define the resolution of the image we want to work with
// A typical 'full-resolution, 40x' pixel size for a whole slide image is around 0.25 microns
// Therefore we will be requesting at a *much* lower resolution here
double requestedPixelSizeMicrons = 20

// The stain that should be thresholded
def stainName = 'DAB'
// Sigma value for a Gaussian filter, used to reduce noise before thresholding
double sigmaMicrons = 20
// Built-in automated threshold method used by ImageJ
def thresholdMethod = AutoThresholder.Method.Otsu

//------------------------------------------------

// Access the relevant QuPath data structures
def imageData = QPEx.getCurrentImageData()
def hierarchy = imageData.getHierarchy()
def server = imageData.getServer()

// For color deconvolution, we need an 8-bit brightfield RGB image, and also stains to be set
// Check for these now, and return if we don't have what we need
def stains = imageData.getColorDeconvolutionStains()
if (!server.isRGB() || !imageData.isBrightfield() || stains == null) {
    println 'An 8-bit RGB brightfield image is required!'
    return
}

// Get the index of the stain we want, based on the specified name
int stainIndex = -1
for (int i = 0; i < 3; i++) {
    // Stains are accessed as 1, 2, 3 (and not 0, 1, 2... sorry...)
    if (stains.getStain(i+1).getName() == stainName) {
        stainIndex = i
        break
    }
}
if (stainIndex < 0) {
    println 'Could not find stain with name ' + stainName + '!'
    return
}

// If we have a selected annotation object with a ROI, use that - otherwise use the entire image
def selectedObject = hierarchy.getSelectionModel().getSelectedObject()
def selectedROI = selectedObject?.isAnnotation() ? selectedObject.getROI() : null

// Convert requestedPixelSizeMicrons into a sensible downsample value
double downsample = requestedPixelSizeMicrons / server.getAveragedPixelSizeMicrons()

// Create a region request, either for the full image or the selected region
def region = selectedROI == null ?
        RegionRequest.createInstance(server.getPath(), downsample, 0, 0, server.getWidth(), server.getHeight()) :
        RegionRequest.createInstance(server.getPath(), downsample, selectedROI)

// Request a PathImage containing an ImagePlus
server = ImagePlusServerBuilder.ensureImagePlusWholeSlideServer(server)
def pathImage = server.readImagePlusRegion(region)
def imp = pathImage.getImage()

// Get the current ImageProcessor - it should be 8-bit RGB, so we can
// ask Groovy to make sure it's an ImageJ ColorProcessor, required for the next step
def ip = imp.getProcessor() as ColorProcessor

// Apply color deconvolution - this gives a list of 3 stain images
def fpDeconvolved = ColorDeconvolutionIJ.colorDeconvolve(ip, stains)

// If we just want to view the images, show them as an ImageJ stack
if (viewOnly) {
    // Ensure the ImageJ user interface is showing
    IJExtension.getImageJInstance()
    // Show the original image
    imp.show()
    // Create an ImageStack for the stains, setting the labels according to stain names
    def stack = new ImageStack(imp.getWidth(), imp.getHeight())
    fpDeconvolved.eachWithIndex { FloatProcessor fpStain, int ind ->
        stack.addSlice(stains.getStain(ind+1).getName(), fpStain)
    }
    // Create a new image for the stains, setting the calibration based on the original image
    // This means we can use 'Send ROI to QuPath' if we want
    def impDeconvolved = new ImagePlus('Color deconvolved ' + imp.getTitle(), stack)
    impDeconvolved.setCalibration(imp.getCalibration().clone())
    // Make the image pseudo-fluorescence, with 3 channels & reset the brightness/contrast for display
    impDeconvolved = new CompositeImage(impDeconvolved, CompositeImage.COMPOSITE)
    impDeconvolved.resetDisplayRanges()
    // Show the image
    impDeconvolved.show()
    return
}

// Extract the stain ImageProcessor
def ipStain = fpDeconvolved[stainIndex].duplicate()

// Convert blur sigma to pixels & apply if > 0
double sigmaPixels = sigmaMicrons / requestedPixelSizeMicrons
if (sigmaPixels > 0)
    ipStain.blurGaussian(sigmaPixels)

// Set the threshold
ipStain.setAutoThreshold(thresholdMethod, true)

// Create a selection; this will use the current threshold
def tts = new ThresholdToSelection()
def roiIJ = tts.convert(ipStain)

// Make some measurements
ipStain.setRoi(roiIJ)
def stats = ImageStatistics.getStatistics(ipStain,
        ImageStatistics.MEAN + ImageStatistics.MIN_MAX + ImageStatistics.AREA, imp.getCalibration())

// Convert ImageJ ROI to a QuPath ROI
// Here, the pathImage comes in handy because it has the calibration info we want
def roi = ROIConverterIJ.convertToPathROI(roiIJ, pathImage)

// Create a QuPath annotation
def annotation = new PathAnnotationObject(roi)

// Add the measurements to the annotation's MeasurementList
// While we're here, add the threshold value as a measurement too
// Warning!  This uses the *smoothed* stain image - it might be better to use the original!
def measurementList = annotation.getMeasurementList()
measurementList.putMeasurement('Threshold (IJ)', ipStain.getMinThreshold())
measurementList.putMeasurement('Area (IJ)', stats.area)
measurementList.putMeasurement('Mean ' + stainName + ' (IJ)', stats.mean)
measurementList.putMeasurement('Min ' + stainName + ' (IJ)', stats.min)
measurementList.putMeasurement('Max ' + stainName + ' (IJ)', stats.max)
// Once we're done adding measurements, it's important to close the list -
// this can allow QuPath to optimize the storage (very important if there are many objects)
measurementList.closeList()

// By default, lock the annotation to make it harder to accidentally edit
annotation.setLocked(true)

// Add the annotation back to the hierarchy
// If we have a parent object, we can tell QuPath to put it below that -
// if we leave it up to QuPath to figure out that relationship, it might put the object
// *beside* rather than within the intended parent, because of potential slight pixel shifts
// in the conversion to/from ImageJ meaning the new annotation is not fully contained inside the parent
if (selectedObject != null)
    imageData.getHierarchy().addPathObjectBelowParent(selectedObject, annotation, false, true)
else
    imageData.getHierarchy().addPathObject(annotation, false)