Calling Image Alignment function from a script

The problem is most likely due to scanning artifacts on the slides I tested the script on. Many of them have parts of the tissue touching or cut off at the edge of the scanned image.

The transformation improved slightly when I tried to use pixel size 5 in the alignment, which suggests that createInverse is set correctly.

I’m working on clinical samples, so I’m hesitant to share the slides.

But if you do have a set of serial sections named in the format of slideID_tissue_stain, e.g.:

A_brain_H&E.ndpi
A_brain_GFAP.ndpi
A_brain_IBA1.ndpi
B_brain_H&E.ndpi
B_brain_GFAP.ndpi
B_brain_IBA1.ndpi
C_brain_H&E.ndpi
C_brain_GFAP.ndpi
C_brain_IBA1.ndpi

the scripts should work for you to test!

I’m glad it’s (mostly) working!

How do the transformations look when you use the original plugin?

Here are my test alignments at:

Pixel size 20

Pixel size 5

Reference slide is H&E. It aligns well for the middle image, but not so much the right image. It did improve with pixel size 5 vs 20, but I wish it could be a bit better.

That looks pretty good given the damage/differences to the tissue between images. To get better you would need deformable image registration, and even that often makes things worse in the case of situations like the middle image where areas of tissue are “missing” (folded back). as the deformable alignment tries to stretch existing tissue into that space to fill it up.

While it doesn’t work with the automatic alignment, you might want to manually try Point annotations to register particular landmarks that are most important to you.
image

@Mike_Nelson I tested the alignment on more images and I agree with you that on pixel size 5, they are pretty good! I can accept if the alignment does not conform around missing tissue areas, I understand that it is not easy to implement.

It also seems that if the tissue touches or is clipped at the edges, the alignment gets thrown off. I guess the solution here is to rescan the slides and ensure that the whole tissue area is scanned.

1 Like

New version!

The script does 3 new things:

  1. Allows for annotation-based alignment within a script
  2. Handles aligning images with very different pixel sizes.
  3. Automatically grabs the detection objects in the second image and transfers them to the “base” image.

Basic workflow: have 2 images in a project of the same slide acquired at different resolutions. Create a tissue annotation for each. Create detection objects in image1. Change otherFile in the script to be the name of image1. Go to image2. Run this script. See where the image1 detections appear in image 2.

For annotation-based alignment, I copied the QuPath source code with a little bit of modification for it to run without the interactive window.

I also modified the QuPath code for dealing with the different pixel sizes. It now uses different downsample amounts for the 2 images (to be the same pixel size), returns the ratio of the two as a value, and scales the final transform by that ratio. Note: the alignment will fail if one of your images does not have a pixel size in microns.

Thanks to @Mike_Nelson for his help, especially with the function to transfer objects between images!

/***
 * Script to align 2 images in the project with different pixel sizes, using either intensities or annotations.
 * Run from base image and include the name of the template image in otherFile.
 * 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.
 *
 * Needed inputs: name of template image, what values to use for alignment (intensity or annotation), what type of registration (RIGID or AFFINE), what pixel sizes to use for iterative registration refinement. 
 *
 * Written by Sara McArdle of the La Jolla Institute, 2020, with lots of help from Mike Nelson.
 *
 */
import static qupath.lib.gui.scripting.QPEx.*

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 qupath.lib.gui.dialogs.Dialogs
import qupath.lib.images.servers.ImageServer
import qupath.lib.images.servers.PixelCalibration
import qupath.lib.regions.RegionRequest
import qupath.opencv.tools.OpenCVTools

import java.awt.Graphics2D
import java.awt.Transparency
import java.awt.color.ColorSpace
import java.awt.geom.AffineTransform
import java.awt.image.ComponentColorModel
import java.awt.image.DataBuffer
import java.awt.image.BufferedImage
import javafx.scene.transform.Affine

import org.bytedeco.javacpp.indexer.FloatIndexer
import org.bytedeco.opencv.opencv_core.Mat
import org.bytedeco.opencv.opencv_core.TermCriteria
import org.bytedeco.opencv.global.opencv_core
import org.bytedeco.opencv.global.opencv_video
import org.bytedeco.javacpp.indexer.Indexer

//set these things
String alignmentType="AREA"  //"AREA" will use annotations. Any other string will use intensities.
otherFile='Zinc_raw.tiff' //Name of template image

//collect basic information
baseImageName = getProjectEntry().getImageName()
def imageDataBase=getCurrentImageData()
def imageDataSelected=project.getImageList().find{it.getImageName()==otherFile}.readImageData()
double otherPixelSize=imageDataSelected.getServer().getPixelCalibration().getAveragedPixelSizeMicrons()


//perform affine transforms. Start with a large pixel size and iteratively refine the transform. NOTE: the final pixel size in the last autoAlignPrep SHOULD NOT be smaller than the pixel size of your least resolved image
Affine affine= []
autoAlignPrep(200.0,alignmentType,imageDataBase,imageDataSelected,affine,"RIGID")
autoAlignPrep(100.0,alignmentType,imageDataBase,imageDataSelected,affine,"RIGID")
def scaleFactor=autoAlignPrep(otherPixelSize,alignmentType,imageDataBase,imageDataSelected,affine,"AFFINE") //may want to change to RIGID, depending on your application.

//deal with the differing downsample amounts
affine.prependScale(1/scaleFactor,1/scaleFactor)

//save the transform as an object in the project folder
def matrix = []
matrix << affine.getMxx()
matrix << affine.getMxy()
matrix << affine.getTx()
matrix << affine.getMyx()
matrix << affine.getMyy()
matrix << affine.getTy()

affinepath = buildFilePath(PROJECT_BASE_DIR, 'Affine '+baseImageName)
mkdirs(affinepath)
path = buildFilePath(PROJECT_BASE_DIR, 'Affine '+baseImageName,  otherFile+".aff")

new File(path).withObjectOutputStream {
    it.writeObject(matrix)
}


//gather all annotation objects in the otherFile
new File(affinepath).eachFile{ f->
    GatherObjects(false, true, f)
}


/*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) 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();


    } else {
        // Default - just use intensities
        logger.debug("Image alignment using intensities");
        serverBase = imageDataBase.getServer();
        serverSelected = imageDataSelected.getServer();
    }

    scaleFactor=autoAlign(serverBase, serverSelected, registrationType, affine, requestedPixelSizeMicrons);
    return scaleFactor
}

double autoAlign(ImageServer<BufferedImage> serverBase, ImageServer<BufferedImage> serverOverlay, String regionstrationType, Affine affine, double requestedPixelSizeMicrons) {
    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()))

    imgBase = ensureGrayScale(imgBase)
    imgOverlay = ensureGrayScale(imgOverlay)

    Mat matBase = OpenCVTools.imageToMat(imgBase)
    Mat matOverlay = OpenCVTools.imageToMat(imgOverlay)

    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
}
4 Likes

Hi Sara, sorry to revive an old post.

I just got round to trying the alignment script in a project and I was wondering what the result from the code below actually reports.

I found some details on the OpenCV documentation for findTransformECC() but I don’t quite understand how a single number can represent the “geometric transform (warp) between two images”.

Can the number reported in result be interpreted as how good a fit the transformation/alignment is? E.g. if result = 0.951762537, would that mean the alignment is ~95.17% accurate?

From my observations, larger numbers from result do seem to have a better alignment compared to smaller ones.

EDIT: I checked and saw that this code is present in the QuPath source code, so I assume it’s been implemented by @petebankhead ; by any chance would you know what result reports?

Many thanks,
Yau

Hi @ym.lim I’m not sure what it means myself… the OpenCV docs don’t seem to explain, although it is used in an example here: OpenCV: samples/cpp/image_alignment.cpp

It seems -1 is interpreted as a special value.

My guess is that it may be the ‘Enhanced Correlation Coefficient’ referenced here:
https://docs.opencv.org/4.5.0/dc/d6b/group__video__track.html#gae94247c0014ff6497a3f85e60eab0a66
But I don’t know for sure – you might need to dig into the OpenCV source code for a definitive answer.

The transform is stored in matTransform, which is what QuPath really needs.

If it is indeed the Enhanced Correlation Coefficient (ECC), then I found some information about it in OpenCV’s Image Alignment tutorial.

It says:

The ECC image alignment algorithm introduced in OpenCV 3 is based on a 2008 paper titled Parametric Image Alignment using Enhanced Correlation Coefficient Maximization by Georgios D. Evangelidis and Emmanouil Z. Psarakis. They propose using a new similarity measure called Enhanced Correlation Coefficient (ECC) for estimating the parameters of the motion model.

According to the emphasis I made, I think it is a measure of how similar the reference and target images are, and it probably makes sense that the more similar the images are, the better the alignment performance. Anyone more familiar please correct me if I am wrong.

@petebankhead I read through that tutorial and was wondering if there is any particular reason you did not implement the homography motion model?

2 Likes

Java has an AffineTransform class and Java Topology Suite has an AffineTransformation – QuPath can use these to efficiently support image rendering and ROI transforms. I expect that a lot of additional code would be required to support the homography motion model, and I didn’t think it was justified.

3 Likes

Hi @smcardle and @Mike_Nelson, I have a quick question. I have dozens of tumor sections, each section was stained separately with a marker (IHC). The whole slide images are either in .ndpi or .tif. Now, I would like to align all the images of each section and get a single file with each marker in a separate channel. Do you think your script can be applied to this task? Thanks a lot

@kaizen89 Have a look into

Maybe it will help during your investigation and maybe - if there are only dozens and not hundreds of sections - using the ImageCombiner will be more easy than writing a script.

4 Likes

This can definitely be done, but the resulting image size would be huge (I ended up creating downsampled versions of multiple DAB stains for visualization, but they were not used for analysis). There are other posts about converting the resulting images to 8bit in order to make the file sizes somewhat more reasonable, but if your whole slide images have JPEG encoding, and many do, the final quasi-IF images will be enormous. @phaub’s solution is nice if you can get it to work for your project, but right off the bat I would be worried since you said sequential slices.

Those almost never go well - I have had scientists start wanting to look for double (or more) positive cells and the like, which does not work with sequential slices! It does not even work well with restains on the same slice without deformable image alignment (sometimes simply due to stitching variation in tiles). The two posts here summarize that nicely.

The scripts above scripts were more intended to bypass most of these issues passing objects back and forth between the images, and can be run using the automated image alignment - resulting in scripts that were entirely automated (across dozens of projects). Specifically, objects over fairly large areas compared to the errors in the slide alignment.

Depending on the volume of images you have to analyze for your project, you may want a more or less user-interactive experience.

3 Likes

My images are actually from the same slice. But each slice has been sequentially stained with one marker. The idea is now to align those images to overlay those markers.

1 Like

Ah, good, that’s what we did. Still could not be used for double/triple positive cells due to variations in stitching (and the tissue contracts during the wash and restain phases at times), but at least it was close enough for general metrics within broad areas. Maybe your scanner will be better.

And transferring the objects between images was “cleaner” than creating a new image with all of the channels, since the pixel alignment would never quite match up and sacrifices in image quality would need to be made. Resolving two images that are shifted about half a pixel and rotated a few degrees will result in a bit of mismatch.

1 Like

If you think that’s interesting, you could try this workflow (Registration of whole slide images at cellular resolution with Fiji and QuPath - #33 by AntoineRG, shameless self promotion :grimacing:). It’s not doing exactly what you want because the images are not transformed. Instead, you save, for each registration, a transformation file (affine or spline) which can be used to transfer annotations or detections from one slide to any other registered one (exactly like @Mike_Nelson’s suggestion) . By NOT computing the transformed image, you can gain a lot of time.

If you really need to, there’s a way to re-save all the transformed images, sampled identically to your “reference” image, but it’s only doable within Fiji (they can be re-imported in QuPath, but the full computation will take time).

2 Likes

It sounds like he only wants the transformations, so that part sounds like a great fit!

I only looked briefly at the protocol and my main question is how manual labor intensive it might be. If our group (as an example, not related to @kaizen89) was to go back to the affine transforms that were calculated and redid the analysis using non-affine registration, how automated would the process be across 50 whole slide samples, each of which has 9 re-stains?

Also, @smcardle and I had just been talking about your method, so I am glad you chimed in :slight_smile:

1 Like

I only looked briefly at the protocol and my main question is how manual labor intensive it might be.

It really depends on the sample. I demo on purpose a pretty annoying situation: one sample is acquired over a larger area than the other one, with a different camera ( → different pixel size), with an offset in the initial positioning, and the stitching was not good for the DAB image, while it was ok for the fluorescent image.

In this situation, you’ll most probably need some manual work, but the fact that we have a UI at our disposal to perform manual correction if needed, is great! (the result of the automated registration is used as BigWarp’s input).

Now if your sample is initially approximately aligned and covers more or less the same physical area, making a script which performs spline registrations in batch should be easy.

However two things are still missing for a good automated workflow:

  • there’s no ‘cost’ function which tells you how good or bad is your automated registration, (and in which region). To check if the results are ok, you still need to rely on “manual” investigation
  • you need to place landmarks in ‘information rich’ places. For automation, if you place a grid of landmarks with a regular spacing, you need to find a criteria to get rid of the landmarks located at low information places. It should be doable to implement by thresholding the image and excluding landmarks outside the mask, but it’s not there.

For your specific example, making an attempt at full automation should be pretty easy (no results guaranteed!) with a few scripts. “Inputing” the previous affine transform result in the workflow should be easy and will help.

Finally there is this Bio-Formats limitation… if your files cannot be read with bio-formats, then a conversion is required, and that makes this workflow suddenly much less interesting.

1 Like

If I only have the transformations, I will be able to make an overlay of all the stains of one slide? At this time, it does not matter if the process can not be automated.

I do not know - but it does not look like it. You might be trading accuracy of the alignment for the ability to view it all together easily. But who knows, maybe @phaub and @NicoKiaru together could come up with something that would let someone view multiple images at the same time using the mesh.

Remember that creating the combined image will necessarily result in worse quality images vs the way both @NicoKiaru and our pipelines were designed to work. Pixels themselves are not really that malleable.

Open to other points of view, though.

1 Like