Making a memory efficient script to get the inner and outer margins of an annotation in qupath

Question

I'd like to make a script in QuPath that accepts an overall tissue annotation and a "tumor" annotation within that annotation and finds 6 different regions inside and around the cancer tissue:

1. Non-tumor, non-margin tissue
2. 60-120 micron region around the tumor
3. 0-60 micron region around the tumor
4. 0-60 micron region within the tumor
5. 60-120 micron region within the tumor
6. Remaining tumor tissue

So far I've made the Groovy script below, but it's quite memory and time inefficient.

import org.locationtech.jts.geom.Geometry
import qupath.lib.common.GeneralTools
import qupath.lib.objects.PathObject
import qupath.lib.objects.PathObjects
import qupath.lib.roi.GeometryTools
import qupath.lib.roi.ROIs

import static qupath.lib.gui.scripting.QPEx.*

//-----
// Some things you might want to change

// How much to expand each region
double expandMargin60 = 60.0
double expandMargin120 = 120.0

// Define the colors
def colorInnerMargin60 = getColorRGB(0, 150, 200)
def colorInnerMargin120 = getColorRGB(150, 0, 200)
def colorInnerMargin612 = getColorRGB(0, 0, 200)
def colorOuterMargin60 = getColorRGB(0, 200, 0)
def colorOuterMargin120 = getColorRGB(0, 200, 150)
def colorOuterMargin612 = getColorRGB(150, 200, 0)
def colorCore = getColorRGB(200, 0, 0)
def colorRemainder = getColorRGB(200, 200, 200)
def colorNonTumor = getColorRGB(200, 0, 200)

// Choose whether to lock the annotations or not (it's generally a good idea to avoid accidentally moving them)
def lockAnnotations = true

//-----

// Extract the main info we need
def imageData = getCurrentImageData()
def hierarchy = imageData.getHierarchy()
def server = imageData.getServer()

// We need the pixel size
def cal = server.getPixelCalibration()
if (!cal.hasPixelSizeMicrons()) {
  print 'We need the pixel size information here!'
  return
}
if (!GeneralTools.almostTheSame(cal.getPixelWidthMicrons(), cal.getPixelHeightMicrons(), 0.0001)) {
  print 'Warning! The pixel width & height are different; the average of both will be used'
}

// Get annotation & detections
def annotations = getAnnotationObjects()
def selected = getSelectedObject()
if (selected == null || !selected.isAnnotation()) {
  print 'Please select an annotation object!'
  return
}

// We need one selected annotation as a starting point; if we have other annotations, they will constrain the output
annotations.remove(selected)

// Extract the ROI & plane
def roiOriginal = selected.getROI()
def plane = roiOriginal.getImagePlane()

// If we have at most one other annotation, it represents the tissue
Geometry areaTissue
PathObject tissueAnnotation
if (annotations.isEmpty()) {
  areaTissue = ROIs.createRectangleROI(0, 0, server.getWidth(), server.getHeight(), plane).getGeometry()
} else if (annotations.size() == 1) {
  tissueAnnotation = annotations.get(0)
  areaTissue = tissueAnnotation.getROI().getGeometry()
} else {
  print 'Sorry, this script only support one selected annotation for the tumor region, and at most one other annotation to constrain the expansion'
  return
}
// Getting area of original annotation
def areaTumor = roiOriginal.getGeometry()

// Calculate how much to expand
double expandPixels60 = expandMargin60 / cal.getAveragedPixelSizeMicrons()
double expandPixels120 = expandMargin120 / cal.getAveragedPixelSizeMicrons()

// Making classes for Hierarchy Resolution
def nonTumorClass = getPathClass('Non_tumor')
def nonMarginClass = getPathClass('Non_margin')
def outer120Class = getPathClass('Outer120')
def outer60Class = getPathClass('Outer60')
def out60120Class = getPathClass('Outer60-120')
def inner120Class = getPathClass('Inner120')
def inner60Class = getPathClass('Inner60')
def inner60120Class = getPathClass('Inner60-120')
def coreClass = getPathClass('Core')

//----------------------OUTER AREA------------------------

// Status update
print("Calculating outer area")

// Get the stroma-side 120 microns
def geomOuter120 = areaTumor.buffer(expandPixels120)
geomOuter120 = geomOuter120.difference(areaTumor)
geomOuter120 = geomOuter120.intersection(areaTissue)
geomOuter120 = GeometryTools.ensurePolygonal(geomOuter120)
def roiOuter120 = GeometryTools.geometryToROI(geomOuter120, plane)
def annotationOuter120 = PathObjects.createAnnotationObject(roiOuter120)
annotationOuter120.setName("Outer_margin_120")
annotationOuter120.setPathClass(outer120Class)
annotationOuter120.setColorRGB(colorOuterMargin120)

// Get the stroma-side 60 microns
def geomOuter60 = areaTumor.buffer(expandPixels60)
geomOuter60 = geomOuter60.difference(areaTumor)
geomOuter60 = geomOuter60.intersection(areaTissue)
geomOuter60 = GeometryTools.ensurePolygonal(geomOuter60)
def roiOuter60 = GeometryTools.geometryToROI(geomOuter60, plane)
def annotationOuter60 = PathObjects.createAnnotationObject(roiOuter60)
annotationOuter60.setName("Outer_margin_60")
annotationOuter60.setPathClass(outer60Class)
annotationOuter60.setColorRGB(colorOuterMargin60)

// Get the stroma-side 60-120 micron region
//def geomOuter612 = geomOuter60.buffer(expandPixels60) // This seems to work better in general
def geomOuter612 = geomOuter120 // For some reason, this one works better with hierarchy resolution
geomOuter612 = geomOuter612.difference(annotationOuter60.getROI().getGeometry())
geomOuter612 = geomOuter612.intersection(annotationOuter120.getROI().getGeometry())
geomOuter612 = GeometryTools.ensurePolygonal(geomOuter612)
def roiOuter612 = GeometryTools.geometryToROI(geomOuter612, plane)
def annotationOuter612 = PathObjects.createAnnotationObject(roiOuter612)
annotationOuter612.setName("Outer_margin_60-120")
annotationOuter612.setPathClass(out60120Class)
annotationOuter612.setColorRGB(colorOuterMargin612)

// Get non-tumor tissue
def geomNonTumor = areaTissue
geomNonTumor = geomNonTumor.difference(areaTumor)
geomNonTumor = geomNonTumor.intersection(areaTissue)
def roiNonTumor = GeometryTools.geometryToROI(geomNonTumor, plane)
def annotationNonTumor = PathObjects.createAnnotationObject(roiNonTumor)
annotationNonTumor.setName("Non_tumor")
annotationNonTumor.setPathClass(nonTumorClass)
annotationNonTumor.setColorRGB(colorNonTumor)

// Get non-tumor, non-margin tissue
def geomRemainder = areaTissue
geomRemainder = geomRemainder.difference(areaTumor)
geomRemainder = geomRemainder.difference(annotationOuter120.getROI().getGeometry())
geomRemainder = geomRemainder.intersection(areaTissue)
def roiRemainder = GeometryTools.geometryToROI(geomRemainder, plane)
def annotationRemainder = PathObjects.createAnnotationObject(roiRemainder)
annotationRemainder.setName("Non_margin")
annotationRemainder.setPathClass(nonMarginClass)
annotationRemainder.setColorRGB(colorRemainder)

//----------------------INNER AREA------------------------

// Status update
print("Calculating inner area")

// Get the central areas
def geomCentral60 = areaTumor.buffer(-expandPixels60)
geomCentral60 = geomCentral60.intersection(areaTumor)
def geomCentral120 = areaTumor.buffer(-expandPixels120)
geomCentral120 = geomCentral120.intersection(areaTumor)

// Get the inner 120 microns
def geomInner120 = areaTumor
geomInner120 = geomInner120.difference(geomCentral120)
geomInner120 = geomInner120.intersection(areaTissue)
geom = GeometryTools.homogenizeGeometryCollection(geomInner120)
def roiInner120 = GeometryTools.geometryToROI(geomInner120, plane)
def annotationInner120 = PathObjects.createAnnotationObject(roiInner120)
annotationInner120.setPathClass(inner120Class)
annotationInner120.setName("Inner_margin_120")
annotationInner120.setColorRGB(colorInnerMargin120)

// Get the inner 60 microns
def geomInner60 = areaTumor
geomInner60 = geomInner60.difference(geomCentral60)
geomInner60 = geomInner60.intersection(areaTissue)
def roiInner60 = GeometryTools.geometryToROI(geomInner60, plane)
def annotationInner60 = PathObjects.createAnnotationObject(roiInner60)
annotationInner60.setPathClass(inner60Class)
annotationInner60.setName("Inner_margin_60")
annotationInner60.setColorRGB(colorInnerMargin60)

// Get the inner 60-120 micron region
def geomInner612 = geomInner120
geomInner612 = geomInner612.difference(annotationInner60.getROI().getGeometry())
geomInner612 = geomInner612.intersection(annotationInner120.getROI().getGeometry())
def roiInner612 = GeometryTools.geometryToROI(geomInner612, plane)
def annotationInner612 = PathObjects.createAnnotationObject(roiInner612)
annotationInner612.setPathClass(inner60120Class)
annotationInner612.setName("Inner_margin_60-120")
annotationInner612.setColorRGB(colorInnerMargin612)

// Getting remainder of epithelium core
def geomCore = areaTumor
geomCore = geomCore.difference(annotationInner120.getROI().getGeometry())
def roiCore = GeometryTools.geometryToROI(geomCore, plane)
def annotationCore = PathObjects.createAnnotationObject(roiCore)
annotationCore.setPathClass(coreClass)
annotationCore.setName("Core")
annotationCore.setColorRGB(colorCore)

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

// Status update
print("Adding annotations")

// Add the annotations
hierarchy.getSelectionModel().clearSelection()
// hierarchy.removeObject(selected, true)
//def annotationsToAdd = [annotationOuter120, annotationOuter60, annotationOuter612, annotationInner60, annotationInner120, annotationInner612, annotationCore];
def annotationsToAdd = [annotationOuter120, annotationOuter60, annotationOuter612, annotationInner60, annotationInner120, annotationInner612, annotationCore, annotationRemainder, annotationNonTumor];
annotationsToAdd.each {it.setLocked(lockAnnotations)}
if (tissueAnnotation == null) {
  hierarchy.addChildObjects(annotationsToAdd)
} else {
  tissueAnnotation.addChildObjects(annotationsToAdd)
  hierarchy.fireHierarchyChangedEvent(this, tissueAnnotation)
  if (lockAnnotations)
    tissueAnnotation.setLocked(true)
}

// Status update
print("Resolving hierarchy")

// Resolving hierarchy
resolveHierarchy()

// Print
tissueAnnotation

// Manually setting hierarchy
// Tumor and Non-tumor to Original Tissue annotation
parentAnnotation = tissueAnnotation
print(parentAnnotation)
childAnnotation = getAnnotationObjects().findAll{it.getPathClass() == getPathClass("Non_tumor")}[0]
print(childAnnotation)
getCurrentHierarchy().addPathObjectBelowParent(parentAnnotation, childAnnotation, true)
childAnnotation = getAnnotationObjects().findAll{it.getPathClass() == getPathClass("Tumor")}[0]
print(childAnnotation)
getCurrentHierarchy().addPathObjectBelowParent(parentAnnotation, childAnnotation, true)

parentAnnotation = getAnnotationObjects().findAll{it.getPathClass() == getPathClass("Non_tumor")}[0]
print(parentAnnotation)
childAnnotation = getAnnotationObjects().findAll{it.getPathClass() == getPathClass("Non_margin")}[0]
print(childAnnotation)
getCurrentHierarchy().addPathObjectBelowParent(parentAnnotation, childAnnotation, true)
childAnnotation = getAnnotationObjects().findAll{it.getPathClass() == getPathClass("Outer120")}[0]
print(childAnnotation)
getCurrentHierarchy().addPathObjectBelowParent(parentAnnotation, childAnnotation, true)

parentAnnotation = getAnnotationObjects().findAll{it.getPathClass() == getPathClass("Outer120")}[0]
print(parentAnnotation)
childAnnotation = getAnnotationObjects().findAll{it.getPathClass() == getPathClass("Outer60")}[0]
print(childAnnotation)
getCurrentHierarchy().addPathObjectBelowParent(parentAnnotation, childAnnotation, true)
childAnnotation = getAnnotationObjects().findAll{it.getPathClass() == getPathClass("Outer60-120")}[0]
print(childAnnotation)
getCurrentHierarchy().addPathObjectBelowParent(parentAnnotation, childAnnotation, true)


// Status update
print("Finished")

I've also found these two scripts that seem to be more memory efficient, but I'm struggling to find a way to combine them. The first one gets the outer margin and the second gets the inner margin.

Script 1 (outer margin).

/**
0.2.0
 * Script to help with annotating tumor regions, expanding the tumor margin.
 * SEE THREAD HERE FOR DESCRIPTION ON USE: https://forum.image.sc/t/reduce-annotations/24305/12
 * Here, each of the margin regions is approximately 60 microns in width.
 **************************************************************************
 *When starting this script, have one "Tissue" and one "Tumor" annotation.*
 **************************************************************************
 * @author Pete Bankhead
 * @mangled by Svidro because reasons
 * @editor Mike Nelson
 */
classifier = "Tissue"
//createAnnotationsFromPixelClassifier(classifier, 1000000.0, 0.0)

//-----
// Some things you might want to change

// How much to expand each region
double expandMarginMicrons = 60.0
// How many times you want to chop into your annotation. Edit color script around line 115 if you go over 5
int howManyTimes = 2
// Define the colors
// Inner layers are given scripted colors, but gretaer than 6 or 7 layers may require adjustments
def colorOuterMargin = getColorRGB(0, 200, 0)



// Extract the main info we need
def imageData = getCurrentImageData()
def hierarchy = imageData.getHierarchy()
def server = imageData.getServer()
// We need the pixel size
def cal = server.getPixelCalibration()
if (!cal.hasPixelSizeMicrons()) {
  print 'We need the pixel size information here!'
  return
}


// Choose whether to lock the annotations or not (it's generally a good idea to avoid accidentally moving them)
def lockAnnotations = true

//-----
//Setup - Merge all Tumor objects into one, they can be split later. Get Geometries for each object
selectObjectsByClassification("Tumor")
mergeSelectedAnnotations()
double expandPixels = expandMarginMicrons / cal.getAveragedPixelSizeMicrons()
initialTumorObject = getAnnotationObjects().find{it.getPathClass() == getPathClass("Tumor")}
def tumorGeom = getAnnotationObjects().find{it.getPathClass() == getPathClass("Tumor")}.getROI().getGeometry()
def plane = ImagePlane.getDefaultPlane()
def tissueGeom = getAnnotationObjects().find{it.getPathClass() == getPathClass("Other")}.getROI().getGeometry()

//Clean up the Tumor geometry
cleanTumorGeom = tissueGeom.intersection(tumorGeom)
tumorROIClean = GeometryTools.geometryToROI(cleanTumorGeom, plane)


cleanTumor = PathObjects.createAnnotationObject(tumorROIClean, getPathClass("Tumor"))
cleanTumor.setName("CleanTumor")

//Create a list of objects we need to add back in at the end, keep adding to it as we go proceed
annotationsToAdd = []
annotationsToAdd << cleanTumor
/*
addObject(cleanTumor)*/

for (i=0; i<howManyTimes;i++){
    currentArea = annotationsToAdd[annotationsToAdd.size()-1].getROI().getGeometry()
    println(currentArea)
    //Expand from the current area, starting with the tumor
    areaExpansion = currentArea.buffer(expandPixels)
    //Clip off anything outside of the tissue
    areaExpansion = areaExpansion.intersection(tissueGeom)
    //Remove anything that intersects with the tumor
    areaExpansion = areaExpansion.difference(cleanTumorGeom)
    //If we have already expanded once, include the prevous geometry in the exclusion
    if(i>=1){
        for (k=1; k<=i;k++){
            remove = annotationsToAdd[annotationsToAdd.size()-k].getROI().getGeometry()
            areaExpansion = areaExpansion.difference(remove)
            
        }
    }   

    roiExpansion = GeometryTools.geometryToROI(areaExpansion, plane)
    j = i+1
    int nameValue = j*expandMarginMicrons
    annotationExpansion = PathObjects.createAnnotationObject(roiExpansion, getPathClass(nameValue.toString()))

    annotationExpansion.setName("Margin "+nameValue+" microns")
    annotationExpansion.setColorRGB(getColorRGB(20*i, 40*i, 200-30*i))
    annotationsToAdd << annotationExpansion

}
remainingTissueGeom = tissueGeom.difference(cleanTumorGeom)
annotationsToAdd.each{
    remainingTissueGeom = remainingTissueGeom.difference(it.getROI().getGeometry())
}
remainingTissueROI = GeometryTools.geometryToROI(remainingTissueGeom, plane)
remainingTissue = PathObjects.createAnnotationObject(remainingTissueROI)
remainingTissue.setName("Other Tissue")
remainingTissue.setPathClass(getPathClass("Other Tissue"))
addObject(remainingTissue)

// Add the annotations

addObjects(annotationsToAdd)
removeObject(initialTumorObject, true)

fireHierarchyUpdate()
println("Done!")

import org.locationtech.jts.geom.Geometry
import qupath.lib.common.GeneralTools
import qupath.lib.objects.PathObject
import qupath.lib.objects.PathObjects
import qupath.lib.roi.GeometryTools
import qupath.lib.roi.ROIs

import java.awt.Rectangle
import java.awt.geom.Area

Script 2 (inner margin).

/**
0.2.0
 * Script to help with annotating tumor regions, chopping increasing chunks into the tumor.
 * SEE THREAD HERE FOR DESCRIPTION ON USE: https://forum.image.sc/t/reduce-annotations/24305/12?u=research_associate
 * Here, each of the margin regions is approximately 100 microns in width.
Recommended BEFORE running the script if your tumor is on the tissue border:
    1. Do have a simple tissue detection made to limit the borders.
    2. Make an inverse area, delete the original tissue, and use CTRL+SHIFT BRUSH TOOL mostly to create the border of your tumor touching the simple tissue detection edge. You can use the wand tool for the interior of the tumor, but any little pixel that is missed by the wand tool defining the whitespace border of the tumor will result in an expansion. Brush tool them all away.
    3. Invert the inverted tissue detection to get the original tissue back, and then delete the inverted annotation. You now have a tumor annotation that is right up against the tissue border.
    4. Select the tumor and run the script
 * @author Pete Bankhead
 * @mangled by Svidro because reasons
 */

//-----
// Some things you might want to change

// How much to expand each region
double expandMarginMicrons = 20.0
// How many times you want to chop into your annotation. Edit color script around line 115 if you go over 5
int howManyTimes = 2
// Define the colors
// Inner layers are given scripted colors, but gretaer than 6 or 7 layers may require adjustments
def colorOuterMargin = getColorRGB(0, 200, 0)


import org.locationtech.jts.geom.Geometry
import qupath.lib.common.GeneralTools
import qupath.lib.objects.PathObject
import qupath.lib.objects.PathObjects
import qupath.lib.roi.GeometryTools
import qupath.lib.roi.ROIs

import java.awt.Rectangle
import java.awt.geom.Area


// Extract the main info we need
def imageData = getCurrentImageData()
def hierarchy = imageData.getHierarchy()
def server = imageData.getServer()
// We need the pixel size
def cal = server.getPixelCalibration()
if (!cal.hasPixelSizeMicrons()) {
  print 'We need the pixel size information here!'
  return
}


// Choose whether to lock the annotations or not (it's generally a good idea to avoid accidentally moving them)
def lockAnnotations = true

//-----




if (!GeneralTools.almostTheSame(cal.getPixelWidthMicrons(), cal.getPixelHeightMicrons(), 0.0001)) {
    print 'Warning! The pixel width & height are different; the average of both will be used'
}

// Get annotation & detections
def annotations = getAnnotationObjects()
def selected = getSelectedObject()
if (selected == null || !selected.isAnnotation()) {
    print 'Please select an annotation object!'
    return
}

// We need one selected annotation as a starting point; if we have other annotations, they will constrain the output
annotations.remove(selected)

// If we have at most one other annotation, it represents the tissue
Geometry areaTissue
PathObject tissueAnnotation
// Calculate how much to expand
double expandPixels = expandMarginMicrons / cal.getAveragedPixelSizeMicrons()
def roiOriginal = selected.getROI()
def plane = roiOriginal.getImagePlane()
def areaTumor = roiOriginal.getGeometry()

if (annotations.isEmpty()) {
  areaTissue = ROIs.createRectangleROI(0, 0, server.getWidth(), server.getHeight(), plane).getGeometry()
} else if (annotations.size() == 1) {
  tissueAnnotation = annotations.get(0)
  areaTissue = tissueAnnotation.getROI().getGeometry()
} else {
  print 'Sorry, this script only support one selected annotation for the tumor region, and at most one other annotation to constrain the expansion'
  return
}
println("Working, give it some time")


// Get the outer margin area
def geomOuter = areaTumor.buffer(expandPixels)

geomOuter = geomOuter.difference(areaTumor)
geomOuter = geomOuter.intersection(areaTissue)
def roiOuter = GeometryTools.geometryToROI(geomOuter, plane)
def annotationOuter = PathObjects.createAnnotationObject(roiOuter)
annotationOuter.setName("Outer margin")
annotationOuter.setColorRGB(colorOuterMargin)

innerAnnotations = []
innerAnnotations << annotationOuter
//innerAnnotations << selected

for (i=0; i<howManyTimes;i++){


    //select the current expansion, which the first time is outside of the tumor, then expand it and intersect it
    currentArea = innerAnnotations[innerAnnotations.size()-1].getROI().getGeometry()
    println(currentArea)
    /*
    if (getQuPath().getBuildString().split()[1]<"0.2.0-m2"){
        areaExpansion = PathROIToolsAwt.shapeMorphology(currentArea, expandPixels)
    }else {areaExpansion = PathROIToolsAwt.getArea(PathROIToolsAwt.roiMorphology(innerAnnotations[innerAnnotations.size()-1].getROI(), expandPixels))}
    */
    areaExpansion = currentArea.buffer(expandPixels)
    
    areaExpansion = areaExpansion.intersection(areaTumor)
    //println(areaExpansion)
    areaExpansion = areaExpansion.intersection(areaTissue)
    //println(areaExpansion)
    //remove outer areas previously defined as other innerAnnotations
    if(i>=1){
        for (k=1; k<=i;k++){
            remove = innerAnnotations[innerAnnotations.size()-k].getROI().getGeometry()
            areaExpansion = areaExpansion.difference(remove)
            
        }
    }
    roiExpansion = GeometryTools.geometryToROI(areaExpansion, plane)
    j = i+1
    annotationExpansion = PathObjects.createAnnotationObject(roiExpansion)
    int nameValue = j*expandMarginMicrons
    annotationExpansion.setName("Inner margin "+nameValue+" microns")
    annotationExpansion.setColorRGB(getColorRGB(20*i, 40*i, 200-30*i))
    innerAnnotations << annotationExpansion
    //println("innerannotations size "+innerAnnotations.size())
}
//add one last inner annotation that contains the rest of the tumor
core = areaTumor

for (i=1; i<=howManyTimes;i++){
    core = core.difference(innerAnnotations[i].getROI().getGeometry())
}
coreROI = GeometryTools.geometryToROI(core, plane)
coreAnno = PathObjects.createAnnotationObject(coreROI)
coreAnno.setName("Remaining Tumor")
innerAnnotations << coreAnno



// Add the annotations
hierarchy.getSelectionModel().clearSelection()
//hierarchy.removeObject(selected, true)
def annotationsToAdd = innerAnnotations;
annotationsToAdd.each {it.setLocked(lockAnnotations)}
if (tissueAnnotation == null) {
    hierarchy.addPathObjects(annotationsToAdd, false)
} else {
    tissueAnnotation.addPathObjects(annotationsToAdd)
    hierarchy.fireHierarchyChangedEvent(this, tissueAnnotation)
    if (lockAnnotations)
        tissueAnnotation.setLocked(true)
}
println("Done! Wheeeee!")

Basic question; can anyone help me make the first script more time and memory efficient or combine the bottom two scipts to get the inner and outer margins of the tumor?

Answer 1

Just realized I never put out a solution to this question. With a great deal of help from Mike Nelson and Pete Bankhead on the Image.sc forum, I have a working script. Thanks for everyone's input!

Validated in QuPath version 0.4.4
 * Script to help with annotating tumor regions, chopping increasing chunks into the tumor.
 * SEE THREAD HERE FOR DESCRIPTION ON USE: https://forum.image.sc/t/reduce-annotations/24305/12
 **************************************************************************
 *When starting this script, have one "Tissue" and one or more "Tumor" annotations.*
 **************************************************************************
 * @author Pete Bankhead
 * @editor Mike Nelson
 * @editor Fredrick Gootkind
 * Readability improved with ChatGPT4 Code Interpreter
 */

//ADD CODE TO DETECT TISSUE AND OR TUMOR ANNOTATIONS HERE. They must have the classes of "Tumor" and "Tissue"
// Imports baby
import qupath.lib.gui.commands.Commands
import qupath.lib.roi.RoiTools.CombineOp

// Clearing stuff
print("Clearing annotations")
selectObjectsByClassification("Tumor", "Stroma", "Necrosis", "Epithelium")
clearSelectedObjects(false)
clearSelectedObjects()
selectObjectsByClassification(null)
mergeSelectedAnnotations()
selectAnnotations()

// Simple tissue detection
print("Performing simple tissue detection")
createAnnotationsFromPixelClassifier("simple_tissue_detection", 500.0, 500.0, "DELETE_EXISTING", "INCLUDE_IGNORED", "SELECT_NEW")
selectObjectsByClassification(null)
clearSelectedObjects(true)
clearSelectedObjects()
selectObjectsByClassification("Region*")

// Make new annotations
print("Making new tissue annotations")
//createAnnotationsFromPixelClassifier("tumor_stroma_epi_nsclc_v1", 500.0, 1000.0)
createAnnotationsFromPixelClassifier("lf014", 500.0, 1000.0)

//Get single object of each class
print("Subtracting necrosis")
annotation1 = getAnnotationObjects().findAll{p -> p.getPathClass() == getPathClass("Necrosis")}[0]
annotation2 = getAnnotationObjects().findAll{p -> p.getPathClass() == getPathClass("Region*")}[0]

//Select one object, then the other, then subtract the two
getCurrentHierarchy().getSelectionModel().setSelectedObject(annotation1, true);
getCurrentHierarchy().getSelectionModel().setSelectedObject(annotation2, true);
Commands.combineSelectedAnnotations(getCurrentImageData(), CombineOp.SUBTRACT);
fireHierarchyUpdate()
print "Subtractions complete"

// Reclassifying
selectObjectsByClassification("Epithelium", "Stroma")
clearSelectedObjects(false)
clearSelectedObjects()

// Configuration parameters

// How much to expand each region
double expandMarginMicrons = 60
// Number of expansions and erosions. Edit color script around line 115 if you go over 5
int expandCount = 2  //Out beyond the tumor
int erodeCount = 2   //Into the tumor
// Define the colors
// Inner layers are given scripted colors, but gretaer than 6 or 7 layers may require adjustments
def colorOuterMargin = getColorRGB(0, 200, 0)

//Added to to prevent errors, see https://forum.image.sc/t/reduce-annotations-invasion-assay-tumor-margins/24305/26
PrecisionModel PM = new PrecisionModel(PrecisionModel.FIXED)

// Extract the main info we need
def imageData = getCurrentImageData()
def hierarchy = imageData.getHierarchy()
def server = imageData.getServer()
// We need the pixel size
def cal = server.getPixelCalibration()
if (!cal.hasPixelSizeMicrons()) {
  print 'We need the pixel size information here!'
  return
}


// Choose whether to lock the annotations or not (it's generally a good idea to avoid accidentally moving them)
def lockAnnotations = true

//-----
//Setup - Merge all Tumor objects into one, they can be split later. Get Geometries for each object
print("Selecting and merging tumor")
selectObjectsByClassification("Tumor")
mergeSelectedAnnotations()
double expandPixels = expandMarginMicrons / cal.getAveragedPixelSizeMicrons()
initialTumorObject = getAnnotationObjects().find{it.getPathClass() == getPathClass("Tumor")}
def tumorGeom = getAnnotationObjects().find{it.getPathClass() == getPathClass("Tumor")}.getROI().getGeometry()
def plane = ImagePlane.getDefaultPlane()
def tissueGeom = getAnnotationObjects().find{it.getPathClass() == getPathClass("Region*")}.getROI().getGeometry()

//Clean up the Tumor geometry
print("Cleaning up tumor geometry")
cleanTumorGeom = tissueGeom.intersection(tumorGeom)
cleanTumorGeom = GeometryTools.ensurePolygonal(cleanTumorGeom)
tumorROIClean = GeometryTools.geometryToROI(cleanTumorGeom, plane)


cleanTumor = PathObjects.createAnnotationObject(tumorROIClean, getPathClass("Tumor"))
cleanTumor.setName("CleanTumor")

annotationsToAdd = [cleanTumor]

// Run the expansion function for a set number of times (expandCount)
(0..<expandCount).each { i ->
    annotationsToAdd = expandAndAddAnnotation(i, annotationsToAdd, expandPixels, tissueGeom, cleanTumorGeom, PM, plane, expandMarginMicrons)
}

// Process the remaining tissue geometry
print("Processing remaining geometry")
remainingTissueGeom = tissueGeom.difference(cleanTumorGeom)
annotationsToAdd.each {
    remainingTissueGeom = remainingTissueGeom.difference(it.getROI().getGeometry())
}
remainingTissueGeom = GeometryTools.ensurePolygonal(remainingTissueGeom)
remainingTissueROI = GeometryTools.geometryToROI(remainingTissueGeom, plane)
remainingTissue = PathObjects.createAnnotationObject(remainingTissueROI)
remainingTissue.setName("Other Tissue")
remainingTissue.setPathClass(getPathClass("Other Tissue"))
addObject(remainingTissue)

// Add the annotations
print("Adding outer annotations")
addObjects(annotationsToAdd)
removeObject(initialTumorObject, true)


innerAnnotations = [annotationsToAdd[1]]

// Erode the annotations and add the inner ones
innerAnnotations = erodeAndAddInnerAnnotations(erodeCount, innerAnnotations, expandPixels, tumorGeom, tissueGeom, plane, expandMarginMicrons, PM)

// Optionally remove the original tumor annotation
removeObjects(getAnnotationObjects().findAll { it.getPathClass() == getPathClass("Tumor") }, true)
//Finally, add the inner tumor margins
print("Adding inner annotations")
addObjects(innerAnnotations)
resetSelection()

print 'Expansion and dilation complete, see the image.sc forum if you have any issues. \n https://forum.image.sc/t/resolving-image-hierarchy-for-margin-annotation-in-qupath/84127/19?u=mike_nelson'


/**
 * Expands the current annotation area and adds the new annotation to the annotations list.
 *
 * @param iteration The current iteration count.
 * @param annotationsToAdd The list to which new annotations are added.
 * @param expandPixels The number of pixels by which to expand.
 * @param tissueGeom The geometry of the tissue.
 * @param cleanTumorGeom The clean geometry of the tumor.
 * @param precisionModel The precision model.
 * @param plane The image plane.
 * @param expandMarginMicrons The micron count for expansion.
 * @return Updated list of annotations.
 */
def expandAndAddAnnotation(iteration, annotationsToAdd, expandPixels, tissueGeom, cleanTumorGeom, precisionModel, plane, expandMarginMicrons) {
    print("Expanding annotations")
    // Get the last annotation's geometry
    def currentAnnotationGeometry = annotationsToAdd[-1].getROI().getGeometry()

    // Expand the current geometry by the defined pixel amount
    def expandedGeometry = currentAnnotationGeometry.buffer(expandPixels)

    // Restrict the expansion to within the tissue boundary
    expandedGeometry = expandedGeometry.intersection(tissueGeom)

    // Create a geometry that will be used to exclude areas from the expansion
    def exclusionGeometry = cleanTumorGeom
    
    // If there have been previous expansions, exclude those areas as well
    if (iteration >= 1) {
        (1..iteration).each { idx ->
            exclusionGeometry = exclusionGeometry.union(annotationsToAdd[annotationsToAdd.size() - idx].getROI().getGeometry())
        }
    }
    // Subtract any overlapping regions from the expanded geometry
    expandedGeometry = expandedGeometry.difference(exclusionGeometry)

    // Refine the geometry to ensure it's defined with a consistent level of detail
    expandedGeometry = GeometryPrecisionReducer.reduce(expandedGeometry, precisionModel)
    expandedGeometry = GeometryTools.ensurePolygonal(expandedGeometry)
    // Convert the refined geometry to a Region of Interest (ROI)
    def expandedROI = GeometryTools.geometryToROI(expandedGeometry, plane)

    // Create and configure a new annotation based on the expanded ROI
    def annotationNameValue = (iteration + 1) * expandMarginMicrons
    def newAnnotation = PathObjects.createAnnotationObject(expandedROI, getPathClass(annotationNameValue.toString()))
    newAnnotation.setName("Outer Margin ${annotationNameValue} microns")
    newAnnotation.setColor(getColorRGB(200 - 30 * iteration, 40 * iteration, 20 * iteration))

    return annotationsToAdd << newAnnotation
}


/**
 * Erodes the inner annotations iteratively and adds them to the list.
 *
 * @param erodeCount The number of times to erode.
 * @param innerAnnotations The list to which new inner annotations are added.
 * @param expandPixels The number of pixels by which to expand.
 * @param tumorGeom The geometry of the tumor.
 * @param tissueGeom The geometry of the tissue.
 * @param plane The image plane.
 * @param expandMarginMicrons The micron count for expansion.
 * @return Updated list of inner annotations.
 */
def erodeAndAddInnerAnnotations(erodeCount, innerAnnotations, expandPixels, tumorGeom, tissueGeom, plane, expandMarginMicrons, precisionModel) {
    print("Eroding annotations")
    // Iterate over the desired number of erosions
    for (i = 0; i < erodeCount; i++) {
        // Get the last annotation's geometry
        currentAnnotationGeometry = innerAnnotations[-1].getROI().getGeometry()
        
        // Erode (expand inwardly) the current geometry
        erodedGeometry = currentAnnotationGeometry.buffer(expandPixels)

        // Intersect the eroded geometry with the tumor and tissue boundaries to ensure it remains within those bounds
        erodedGeometry = erodedGeometry.intersection(tumorGeom).intersection(tissueGeom)

        // Exclude any previously defined inner annotations from the eroded geometry
        if (i >= 1) {
            (1..i).each { idx ->
                def removeGeom = innerAnnotations[innerAnnotations.size() - idx].getROI().getGeometry()
                removeGeom = GeometryTools.ensurePolygonal(removeGeom)
                erodedGeometry = erodedGeometry.difference(removeGeom)
            }
        }
        erodedGeometry = GeometryTools.ensurePolygonal(erodedGeometry)
        erodedGeometry= GeometryPrecisionReducer.reduce(erodedGeometry, precisionModel)
        // Convert the eroded geometry to a Region of Interest (ROI)
        erodedGeometry = GeometryTools.ensurePolygonal(erodedGeometry)
        erodedROI = GeometryTools.geometryToROI(erodedGeometry, plane)
        
        // Create and configure a new inner annotation based on the eroded ROI
        annotationNameValue = (i + 1) * expandMarginMicrons
        newInnerAnnotation = PathObjects.createAnnotationObject(erodedROI, getPathClass("-${annotationNameValue}"))
        newInnerAnnotation.setName("Inner margin ${annotationNameValue} microns")
        newInnerAnnotation.setColorRGB(getColorRGB(20 * i, 40 * i, 200 - 30 * i))

        innerAnnotations << newInnerAnnotation
    }

    // Compute the core tumor region after all erosions
    def coreGeometry = GeometryTools.ensurePolygonal(tumorGeom)
    (1..erodeCount).each { i ->
        
        coreGeometry= GeometryPrecisionReducer.reduce(coreGeometry, precisionModel)
        coreGeometry = coreGeometry.difference(innerAnnotations[i].getROI().getGeometry())
        coreGeometry = GeometryTools.ensurePolygonal(coreGeometry)
    }
    
    def coreROI = GeometryTools.geometryToROI(coreGeometry, plane)
    def coreAnnotation = PathObjects.createAnnotationObject(coreROI)
    coreAnnotation.setName("Remaining Tumor")

    return innerAnnotations << coreAnnotation
}

// Detect cells
print("Beginning cell detection")
selectObjectsByClassification("Other Tissue", "Tumor Core", "-120.0", "-60.0", "120.0", "60.0")
runPlugin('qupath.imagej.detect.cells.WatershedCellDetection', '{"detectionImage": "DAPI",  "requestedPixelSizeMicrons": 0.5,  "backgroundRadiusMicrons": 8.0,  "medianRadiusMicrons": 0.0,  "sigmaMicrons": 0.5,  "minAreaMicrons": 10.0,  "maxAreaMicrons": 150.0,  "threshold": 10.0,  "watershedPostProcess": false,  "cellExpansionMicrons": 2.0,  "includeNuclei": false,  "smoothBoundaries": true,  "makeMeasurements": true}')

// Classify cells
print("Beginning cell classification")
runObjectClassifier("cell_classifier_v2")

import org.locationtech.jts.precision.GeometryPrecisionReducer
import org.locationtech.jts.geom.PrecisionModel```