--- a +++ b/Calculate-Transforms.groovy @@ -0,0 +1,465 @@ +/*** + * See https://github.com/MarkZaidi/QuPath-Image-Alignment/blob/main/Calculate-Transforms.groovy for most up-to-date version. + * Please link to github repo when referencing code in forums, as the code will be continually updated. + * + * Script to align 2 or more images in the project with different pixel sizes, using either intensities or annotations. + * Run from any image in a project containing all sets of images that require alignment + * Writes the affine transform to an object inside the Affine subfolder of your project folder. + * Also grabs the detection objects from the template image. Can change this to annotation objects. + * Usage: + * - Load in all sets of images to be aligned. Rename file names such that the only underscore (_) in the image name + * separates the SlideID from stain. Example: N19-1107 30Gy M5_PANEL2.vsi + * - Adjust the inputs specified under "Needed inputs", and run script (can run on any image, iterates over entire project) + * - If script errors due to alignment failing to converge, set 'align_specific' to the SlideID of the image it failed on + * - Set 'skip_image' to 1, rerun script to skip over the error-causing image + * - Set 'skip_image' to 0, and either adjust 'AutoAlignPixelSize' or draw tissue annotations on all stains of images in list + * - run script, verify all moving images contain a transform file located in the 'Affine' folder + * + * Needed inputs: + * - registrationType : Set as "AFFINE" for translations, rotations, scaling, and sheering. Set as "RIGID" for only translations and rotations. + * - refStain : Set to stain name of image to align all subsequent images to + * - wsiExt : file name extension + * - align_specific : see above, set to null for first run through + * - AutoAlignPixelSize : downsample factor when calculating the transform. Greater values result in faster calculation, but may impact quality + * - skip_image see above, value doesn't matter if align_specific is null + * + * + * Script largely adapted from Sara McArdle's callable implementation of QuPath's Interactive Image Alignment, and Yau Mun Lim's method + * of matching reference (static) and overlay (moving) images based on file names. + * + */ + +import qupath.lib.objects.PathCellObject +import qupath.lib.objects.PathDetectionObject +import qupath.lib.objects.PathObject +import qupath.lib.objects.PathObjects +import qupath.lib.objects.PathTileObject +import qupath.lib.objects.classes.PathClassFactory +import qupath.lib.roi.RoiTools +import qupath.lib.roi.interfaces.ROI + +import java.awt.geom.AffineTransform +import javafx.scene.transform.Affine +import qupath.lib.images.servers.ImageServer + +import java.awt.Graphics2D +import java.awt.Transparency +import java.awt.color.ColorSpace +import java.awt.image.BufferedImage + +import org.bytedeco.opencv.global.opencv_core; +import org.bytedeco.opencv.opencv_core.Mat; +import org.bytedeco.opencv.opencv_core.TermCriteria; +import org.bytedeco.opencv.global.opencv_video; +import org.bytedeco.javacpp.indexer.FloatIndexer; +import org.bytedeco.javacpp.indexer.Indexer; + +import qupath.lib.gui.dialogs.Dialogs; +import qupath.lib.images.servers.PixelCalibration; + +import qupath.lib.regions.RegionRequest; +import qupath.opencv.tools.OpenCVTools + +import java.awt.image.ComponentColorModel +import java.awt.image.DataBuffer + +import static qupath.lib.gui.scripting.QPEx.*; + +// Variables to set +////////////////////////////////// +String registrationType="AFFINE" //Specify as "RIGID" or "AFFINE" +String refStain = "PTEN" //stain to use as reference image (all images will be aligned to this) +String wsiExt = ".ndpi" //image name extension +//def align_specific=['N19-1107 30Gy M5']//If auto-align on intensity fails, put the image(s) that it fails on here +def AutoAlignPixelSize = 30 //downsample factor for calculating transform (tform). Does not affect scaling of output image +align_specific=null //When referencing an image, just include the slide name (stuff before _) +skip_image=0 // If 1, skips the images defined by 'align_specific'. If 0, skips all but image(s) in 'align_specific' +//Experimental features +use_single_channel=0 // Use a single channel from each image for alignment (set to channel number to use). Set to 0 to use all channels. +iterations=5 // Number of times to iteratively calculate the transformation +mov_rotation=180 // rotation to apply to ALL moving images before calculating alignment. Strongly recommend ensuring proper orientation before loading into QuPath. +decrement_factor=1.1 // if iterations>1, by what factor to decrease AutoAlignPixelSize (increasing resolution of alignment). Set to 1 to leave AutoAlignPixelSize unchanged across iterations. + +///////////////////////////////// + + +//Lim's code for file name matching +// Get list of all images in project +def projectImageList = getProject().getImageList() + +// Create empty lists +def imageNameList = [] +def slideIDList = [] +def stainList = [] +def missingList = [] + +// Split image file names to desired variables and add to previously created lists +for (entry in projectImageList) { + def name = entry.getImageName() + def (imageName, imageExt) = name.split('\\.') + def (slideID, stain) = imageName.split('_') + imageNameList << imageName + slideIDList << slideID + stainList << stain +} + +// Remove duplicate entries from lists +slideIDList = slideIDList.unique() +stainList = stainList.unique() +print (slideIDList) +print (align_specific) +// Remove specific entries if causing alignment to not converge +if (align_specific != null) + if (skip_image == 1) + slideIDList.removeAll(align_specific) + else + slideIDList.retainAll(align_specific) + + +if (stainList.size() == 1) { + print 'Only one stain detected. Target slides may not be loaded.' + return +} + +// Create Affine folder to put transformation matrix files +path = buildFilePath(PROJECT_BASE_DIR, 'Affine') +mkdirs(path) + +// Process all combinations of slide IDs, tissue blocks, and stains based on reference stain slide onto target slides +for (slide in slideIDList) { + for (stain in stainList) { + if (stain != refStain) { + refFileName = slide + "_" + refStain + wsiExt + targetFileName = slide + "_" + stain + wsiExt + path = buildFilePath(PROJECT_BASE_DIR, 'Affine', targetFileName) + + def refImage = projectImageList.find {it.getImageName() == refFileName} + def targetImage = projectImageList.find {it.getImageName() == targetFileName} + if (refImage == null) { + print 'Reference slide ' + refFileName + ' missing!' + missingList << refFileName + continue + } + if (targetImage == null) { + print 'Target slide ' + targetFileName + ' missing!' + missingList << targetFileName + continue + } + println("Aligning reference " + refFileName + " to target " + targetFileName) + //McArdle's code for image alignment + ImageServer<BufferedImage> serverBase = refImage.readImageData().getServer() + ImageServer<BufferedImage> serverOverlay = targetImage.readImageData().getServer() + def static_img_name = refFileName + def moving_img_name = targetFileName + def project_name = getProject() + def entry_name_static = project_name.getImageList().find { it.getImageName() == static_img_name } + def entry_name_moving = project_name.getImageList().find { it.getImageName() == moving_img_name } + + def serverBaseMark = entry_name_static.readImageData() + def serverOverlayMark = entry_name_moving.readImageData() + Affine affine=[] + + + mov_width=serverOverlayMark.getServer().getWidth() + mov_height=serverOverlayMark.getServer().getHeight() + affine.prependRotation(mov_rotation,mov_width/2,mov_height/2) + + + + for(int i = 0;i<iterations;i++) { + //Perform the alignment. If no annotations present, use intensity. If annotations present, use area + print("autoalignpixelsize:" + AutoAlignPixelSize) + if (serverBaseMark.hierarchy.nObjects() > 0 || serverOverlayMark.hierarchy.nObjects() > 0) + autoAlignPrep(AutoAlignPixelSize, "AREA", serverBaseMark, serverOverlayMark, affine, registrationType, use_single_channel) + else + autoAlignPrep(AutoAlignPixelSize, "notAREA", serverBaseMark, serverOverlayMark, affine, registrationType, use_single_channel) + AutoAlignPixelSize/=decrement_factor + if (AutoAlignPixelSize<1){ + AutoAlignPixelSize=1 + } + } + + def matrix = [] + + + matrix << affine.getMxx() + matrix << affine.getMxy() + matrix << affine.getTx() + matrix << affine.getMyx() + matrix << affine.getMyy() + matrix << affine.getTy() + + new File(path).withObjectOutputStream { + it.writeObject(matrix) + } + } + } +} + +if (missingList.isEmpty() == true) { + print 'Done!' +} else { + missingList = missingList.unique() + print 'Done! Missing slides: ' + missingList +} + + +/*Subfunctions taken from here: +https://github.com/qupath/qupath/blob/a1465014c458d510336993802efb08f440b50cc1/qupath-experimental/src/main/java/qupath/lib/gui/align/ImageAlignmentPane.java + */ + +//creates an image server using the actual images (for intensity-based alignment) or a labeled image server (for annotation-based). +double autoAlignPrep(double requestedPixelSizeMicrons, String alignmentMethod, ImageData<BufferedImage> imageDataBase, ImageData<BufferedImage> imageDataSelected, Affine affine,String registrationType, int use_single_channel) throws IOException { + ImageServer<BufferedImage> serverBase, serverSelected; + + if (alignmentMethod == 'AREA') { + logger.debug("Image alignment using area annotations"); + Map<PathClass, Integer> labels = new LinkedHashMap<>(); + int label = 1; + labels.put(PathClassFactory.getPathClassUnclassified(), label++); + for (def annotation : imageDataBase.getHierarchy().getAnnotationObjects()) { + def pathClass = annotation.getPathClass(); + if (pathClass != null && !labels.containsKey(pathClass)) + labels.put(pathClass, label++); + } + for (def annotation : imageDataSelected.getHierarchy().getAnnotationObjects()) { + def pathClass = annotation.getPathClass(); + if (pathClass != null && !labels.containsKey(pathClass)) + labels.put(pathClass, label++); + } + + double downsampleBase = requestedPixelSizeMicrons / imageDataBase.getServer().getPixelCalibration().getAveragedPixelSize().doubleValue(); + serverBase = new LabeledImageServer.Builder(imageDataBase) + .backgroundLabel(0) + .addLabels(labels) + .downsample(downsampleBase) + .build(); + + double downsampleSelected = requestedPixelSizeMicrons / imageDataSelected.getServer().getPixelCalibration().getAveragedPixelSize().doubleValue(); + serverSelected = new LabeledImageServer.Builder(imageDataSelected) + .backgroundLabel(0) + .addLabels(labels) + .downsample(downsampleSelected) + .build(); + //disable single channel alignment when working with Area annotations, unsure what bugs it can cause + use_single_channel=0 + } else { + // Default - just use intensities + logger.debug("Image alignment using intensities"); + serverBase = imageDataBase.getServer(); + serverSelected = imageDataSelected.getServer(); + } + + scaleFactor=autoAlign(serverBase, serverSelected, registrationType, affine, requestedPixelSizeMicrons,use_single_channel); + return scaleFactor +} + +double autoAlign(ImageServer<BufferedImage> serverBase, ImageServer<BufferedImage> serverOverlay, String regionstrationType, Affine affine, double requestedPixelSizeMicrons, use_single_channel) { + PixelCalibration calBase = serverBase.getPixelCalibration() + double pixelSizeBase = calBase.getAveragedPixelSizeMicrons() + double downsampleBase = 1 + if (!Double.isFinite(pixelSizeBase)) { + // while (serverBase.getWidth() / downsampleBase > 2000) + // downsampleBase++; + // logger.warn("Pixel size is unavailable! Default downsample value of {} will be used", downsampleBase) + pixelSizeBase=50 + downsampleBase = requestedPixelSizeMicrons / pixelSizeBase + } else { + downsampleBase = requestedPixelSizeMicrons / pixelSizeBase + } + + PixelCalibration calOverlay = serverOverlay.getPixelCalibration() + double pixelSizeOverlay = calOverlay.getAveragedPixelSizeMicrons() + double downsampleOverlay = 1 + if (!Double.isFinite(pixelSizeOverlay)) { + // while (serverBase.getWidth() / downsampleOverlay > 2000) + // downsampleOverlay++; + // logger.warn("Pixel size is unavailable! Default downsample value of {} will be used", downsampleOverlay) + pixelSizeOverlay=50 + downsampleOverlay = requestedPixelSizeMicrons / pixelSizeOverlay + } else { + downsampleOverlay = requestedPixelSizeMicrons / pixelSizeOverlay + } + + double scaleFactor=downsampleBase/downsampleOverlay + + BufferedImage imgBase = serverBase.readBufferedImage(RegionRequest.createInstance(serverBase.getPath(), downsampleBase, 0, 0, serverBase.getWidth(), serverBase.getHeight())) + BufferedImage imgOverlay = serverOverlay.readBufferedImage(RegionRequest.createInstance(serverOverlay.getPath(), downsampleOverlay, 0, 0, serverOverlay.getWidth(), serverOverlay.getHeight())) + + //Determine whether to calculate intensity-based alignment using all channels or a single channel + Mat matBase + Mat matOverlay + if (use_single_channel==0) { + //print 'using all channels' + imgBase = ensureGrayScale(imgBase) + imgOverlay = ensureGrayScale(imgOverlay) + matBase = OpenCVTools.imageToMat(imgBase) + matOverlay = OpenCVTools.imageToMat(imgOverlay) + + } else { + + matBase = OpenCVTools.imageToMat(imgBase) + matOverlay = OpenCVTools.imageToMat(imgOverlay) + int channel = use_single_channel-1 + //print ('using channel ' + channel) + matBase = OpenCVTools.splitChannels(matBase)[channel] + matOverlay = OpenCVTools.splitChannels(matOverlay)[channel] + //use this to preview how the channel looks + //OpenCVTools.matToImagePlus('Channel:' + channel.toString(), matBase).show() + } + + + /////pete code block///// + +//// New bit +// int channel = 2 +// matBase = OpenCVTools.splitChannels(matBase)[channel] +// matOverlay = OpenCVTools.splitChannels(matOverlay)[channel] +// ///end pete code block/// + + Mat matTransform = Mat.eye(2, 3, opencv_core.CV_32F).asMat() +// Initialize using existing transform +// affine.setToTransform(mxx, mxy, tx, myx, myy, ty) + try { + FloatIndexer indexer = matTransform.createIndexer() + indexer.put(0, 0, (float)affine.getMxx()) + indexer.put(0, 1, (float)affine.getMxy()) + indexer.put(0, 2, (float)(affine.getTx() / downsampleBase)) + indexer.put(1, 0, (float)affine.getMyx()) + indexer.put(1, 1, (float)affine.getMyy()) + indexer.put(1, 2, (float)(affine.getTy() / downsampleBase)) +// System.err.println(indexer) + } catch (Exception e) { + logger.error("Error closing indexer", e) + } + + TermCriteria termCrit = new TermCriteria(TermCriteria.COUNT, 100, 0.0001) + + try { + int motion + switch (regionstrationType) { + case "AFFINE": + motion = opencv_video.MOTION_AFFINE + break + case "RIGID": + motion = opencv_video.MOTION_EUCLIDEAN + break + default: + logger.warn("Unknown registraton type {} - will use {}", regionstrationType, RegistrationType.AFFINE) + motion = opencv_video.MOTION_AFFINE + break + } + double result = opencv_video.findTransformECC(matBase, matOverlay, matTransform, motion, termCrit, null) + logger.info("Transformation result: {}", result) + } catch (Exception e) { + Dialogs.showErrorNotification("Estimate transform", "Unable to estimated transform - result did not converge") + logger.error("Unable to estimate transform", e) + return + } + +// To use the following function, images need to be the same size +// def matTransform = opencv_video.estimateRigidTransform(matBase, matOverlay, false); + Indexer indexer = matTransform.createIndexer() + affine.setToTransform( + indexer.getDouble(0, 0), + indexer.getDouble(0, 1), + indexer.getDouble(0, 2) * downsampleBase, + indexer.getDouble(1, 0), + indexer.getDouble(1, 1), + indexer.getDouble(1, 2) * downsampleBase + ) + indexer.release() + + matBase.release() + matOverlay.release() + matTransform.release() + + return scaleFactor +} + +//to gather detection objects instead of annotation, change line ~250 to def pathObjects = otherHierarchy.getDetectionObjects() +def GatherObjects(boolean deleteExisting, boolean createInverse, File f){ + f.withObjectInputStream { + matrix = it.readObject() + + // Get the project & the requested image name + def project = getProject() + def entry = project.getImageList().find {it.getImageName()+".aff" == f.getName()} + if (entry == null) { + print 'Could not find image with name ' + f.getName() + return + } + + def otherHierarchy = entry.readHierarchy() + def pathObjects = otherHierarchy.getDetectionObjects() //OR getAnnotationObjects() + + // Define the transformation matrix + def transform = new AffineTransform( + matrix[0], matrix[3], matrix[1], + matrix[4], matrix[2], matrix[5] + ) + if (createInverse) + transform = transform.createInverse() + + if (deleteExisting) + clearAllObjects() + + def newObjects = [] + for (pathObject in pathObjects) { + newObjects << transformObject(pathObject, transform) + } + addObjects(newObjects) + } +} + +//other subfunctions + +PathObject transformObject(PathObject pathObject, AffineTransform transform) { + // Create a new object with the converted ROI + def roi = pathObject.getROI() + def roi2 = transformROI(roi, transform) + def newObject = null + if (pathObject instanceof PathCellObject) { + def nucleusROI = pathObject.getNucleusROI() + if (nucleusROI == null) + newObject = PathObjects.createCellObject(roi2, pathObject.getPathClass(), pathObject.getMeasurementList()) + else + newObject = PathObjects.createCellObject(roi2, transformROI(nucleusROI, transform), pathObject.getPathClass(), pathObject.getMeasurementList()) + } else if (pathObject instanceof PathTileObject) { + newObject = PathObjects.createTileObject(roi2, pathObject.getPathClass(), pathObject.getMeasurementList()) + } else if (pathObject instanceof PathDetectionObject) { + newObject = PathObjects.createDetectionObject(roi2, pathObject.getPathClass(), pathObject.getMeasurementList()) + newObject.setName(pathObject.getName()) + } else { + newObject = PathObjects.createAnnotationObject(roi2, pathObject.getPathClass(), pathObject.getMeasurementList()) + newObject.setName(pathObject.getName()) + } + // Handle child objects + if (pathObject.hasChildren()) { + newObject.addPathObjects(pathObject.getChildObjects().collect({transformObject(it, transform)})) + } + return newObject +} + +ROI transformROI(ROI roi, AffineTransform transform) { + def shape = RoiTools.getShape(roi) // Should be able to use roi.getShape() - but there's currently a bug in it for rectangles/ellipses! + shape2 = transform.createTransformedShape(shape) + return RoiTools.getShapeROI(shape2, roi.getImagePlane(), 0.5) +} + +static BufferedImage ensureGrayScale(BufferedImage img) { + if (img.getType() == BufferedImage.TYPE_BYTE_GRAY) + return img + if (img.getType() == BufferedImage.TYPE_BYTE_INDEXED) { + ColorSpace cs = ColorSpace.getInstance(ColorSpace.CS_GRAY) + def colorModel = new ComponentColorModel(cs, 8 as int[], false, true, + Transparency.OPAQUE, + DataBuffer.TYPE_BYTE) + return new BufferedImage(colorModel, img.getRaster(), false, null) + } + BufferedImage imgGray = new BufferedImage(img.getWidth(), img.getHeight(), BufferedImage.TYPE_BYTE_GRAY) + Graphics2D g2d = imgGray.createGraphics() + g2d.drawImage(img, 0, 0, null) + g2d.dispose() + return imgGray +}