Multistep, automatic and nonrigid image registration method for histology samples acquired using multiple stains

There are many works related to the medical image registration including the use of different transformations, optimization algorithms, similarity metrics (Sotiras et al 2013, Oliveira and Tavares 2014) and recently also methods using deep learning approach (Litjens et al 2017). However, due to the mentioned difficulties, the generally applicable methods have limited usefulness for the registration of histological images. Firstly due to the resulting registration accuracy, and secondly due to a high computation time (Borovec et al 2018, 2020). The ANHIR challenge organizers evaluated several openly available nonrigid registration algorithms like bUnwarpJ (Arganda-Carreras et al 2006), RVSS (Arganda-Carreras et al 2006), NiftyReg (Rueckert 1999), Elastix (Klein et al 2010), ANTs (Avants et al 2008) or DROP (Glocker et al 2011). However, all of them provided rather poor results compared to the dedicated algorithms proposed by the challenge participants (Borovec et al 2020).

Majority of the works about registration of histology images involve the registration between histology and magnetic resonance images. There is much less work done for registration of histological images acquired using multiple stains. However, there are some interesting contributions made prior to the ANHIR challenge. An intriguing idea was presented in Song et al (2014) where the authors introduced a concept about an intensity-based nonrigid registration algorithm driven by an unsupervised content classification enhancing the structural similarity. Another important contribution was presented in Shojaii et al (2011) where authors proposed use of fiducial markers to calculate the target registration error and improve the registration evaluation process. The bUnwarpJ, which is an elastic and consistent image registration algorithm is also noteworthy since a lot people working with histology images use it on a daily basis (Arganda-Carreras et al 2006).

Several other methods were proposed directly for the challenge. The best scoring method proposed by the MEVIS team (Lotz et al 2019), similar to our approach, is a multistep procedure and consists of: (i) exhaustive, initial rotation search, (ii) iterative affine registration based on the normalized gradient field (NGF) similarity metric (Haber and Modersitzki 2006), (iii) B-Splines-based nonrigid registration using the NGF similarity metric. The main advantage of their method is a commercial, well optimized, matrix-free implementation which makes it extremely fast. The method is able to register the images in just a few seconds, even using resolution up to 8000 pixels in the larger dimension. However, since the method is commercial, there is no open implementation or free access available. The method proposed by the UPENN team (Venet et al 2019) also consists of three steps: (i) initial preprocessing involving images resampling and background removal by deconvolution, (ii) a random orientation search and a following affine registration, (iii) use of a general-purpose medical image registration tool named 'Greedy' which is an implementation of a greedy diffeomorphic registration algorithm (Joshi et al 2004, Yushkevich et al 2016). Only one of the best scoring teams decided to use a deep learning-based solution. The TUB team used a modified volume tweening network, described by Zhao et al (2020). The network was first unsupervisedly trained and then fine-tuned using the training landmarks which may be controversial. The description of the methods proposed by TUNI and UA teams can be be found in Kartasalo et al (2018) and Punithakumar et al (2017) respectively.

In this work, we describe in detail the method proposed by the AGH team for the ANHIR challenge and provide access to the method software. The method is fully automatic, robust to the applied dye and the tissue type, and is currently the second-best method proposed for the challenge in terms of the average ranking of the median target registration error, without a statistically significant difference to the top-ranked method (according to the one-sided Wilcoxon signed-rank test with p = 0.01) (Borovec et al 2020).

We propose a method fulfilling all the ANHIR challenge requirements. The method is fully automatic, robust, resistant to staining using multiple dyes, reasonably fast and calculates an accurate dense deformation field. The algorithm works well with predefined parameters, without the necessity to tune them to each image pair separately.

The proposed method consists of several sequential steps. First, both the moving and the fixed images are preprocessed to make the registration process as easy as possible. Then, the algorithm attempts to initially align the images using a feature-based approach. In case of failure, an exhaustive rotation search is performed. As a next step, the initial affine transformation is being refined by a dedicated affine registration resistant to missing data. Finally, a fully nonrigid, multimodal registration is used to calculate the final deformation field. The algorithm flow is presented in algorithm 1. An exemplary visualization of the images after each step is shown in figure 1.

Algorithm 1: The proposed algorithm.
Input : M (moving image), F (fixed image)
Output : u (deformation field)
$\mathbf{M_P}$, $\mathbf{F_P}$ = preprocess the $\mathbf{M}$, $\mathbf{F}$ images (algorithm 2)
u = calculate the feature-based or the exhaustive initial alignment using $\mathbf{M_P}$,
$\mathbf{F_P}$ (algorithm 4)
u = calculate the initial alignment refinement between $\mathbf{M_P}$ and $\mathbf{F_P}$
(algorithm 3)
u = perform the nonrigid registration using $\mathbf{M_P}$, $\mathbf{F_P}$ and u as the initial
deformation (algorithm 5)
return u

Zoom In Zoom Out Reset image size

The preprocessing consists of the following steps. Firstly, the data is converted to the grayscale, smoothed and downsampled to a lower resolution. The downsampling is required because the computation time using the original resolution is extremely high and the improvement of the registration quality is not that significant considering the unacceptable increase of the computational time. The resolution is different for the initial alignment, the exhaustive search, the affine transformation tuning, and the nonrigid registration. The images are resampled in a way that makes the smaller dimension equal to 2048, 512, 1024, 4096 pixels, respectively for the initial alignment, the exhaustive search, the affine transformation tuning, and the nonrigid registration. After the downsampling, the image entropy is calculated for both images and the histogram of the image with a lower entropy is matched to the histogram of the image with higher entropy. This step makes it possible to use a traditional feature-based approach to calculate the similarity or the affine transformation. The step is not required for the following affine transformation tuning or the nonrigid registration. Finally, the images are padded to the same size and the intensity is inverted (so the background is equal to 0). Padding the images to the same shape simplifies operations involving both the images (e.g. the similarity metric calculation) because it limits the necessity to check the boundary conditions and improves the cache memory utilization. The proposed preprocessing is presented in algorithm 2.

Algorithm 2: Preprocessing.
Input : M (moving image), F (fixed image)
Output : $\mathbf{M_P, F_P}$ (preprocessed images)
$\mathbf{M_P}$, $\mathbf{F_P}$ = convert both the $\mathbf{M}$, $\mathbf{F}$ images to the grayscale
$\mathbf{M_P}$, $\mathbf{F_P}$ = resample the gray $\mathbf{M_P}$, $\mathbf{F_P}$ images to a given resolution
$\mathbf{E_M}$, $\mathbf{E_F}$ = calculate the $\mathbf{M_P}$, $\mathbf{F_P}$ images entropy
$\mathbf{M_P}$, $\mathbf{F_P}$ = match the histogram of the image with a lower entropy to the
image with the higher entropy based on the $\mathbf{E_M}$, $\mathbf{E_F}$
$\mathbf{M_P}$, $\mathbf{F_P}$ = pad the $\mathbf{M_P}$, $\mathbf{F_P}$ images to the same shape and invert the intensity
return $\mathbf{M_P}$, $\mathbf{F_P}$

The histological images are often misaligned with a low overlap ratio making the process of the initial alignment essential for the following nonrigid registration. Even though an intensity-based affine or rigid registration is sufficient for many cases, it fails for numerous of them due to large displacements (e.g. due to almost 180 degrees rotation) or missing data resulted from the use of distinct dyes (e.g. 20% of tissue being invisible in the fixed image).

The initial alignment during the proposed algorithm is based on the similarity or the affine transformation. The initial alignment consists of several subsequent steps: (i) the feature-based affine registration, (ii) the exhaustive initial rotation search (in case of a failure of the feature-based registration), and (iii) the iterative, intensity-based affine registration (in all cases). In the first step, three sets of matching features are calculated, using SIFT (Lowe 2004), ORB (Rublee et al 2110), and SURF (Bay et al 2006). The use of three descriptors is motivated by the necessity to make the algorithm fully automatic, without manual method selection or parameter tuning. The use of three feature descriptors is not strictly mandatory but it makes the algorithm more robust, and the computational time is negligible compared to the following registration steps. Then, the RANSAC algorithm is used to calculate the transformation between the filtered keypoints. Finally, the evaluation of the transformation quality is determined using a failure detector (described later in this section). The Dice coefficient is calculated between appropriately thresholded masks, as described in section 2.5. Transformation with the highest coefficient is selected. The use of this procedure successfully aligns more than 70% of the image pairs, the remaining ones are handled by an exhaustive search or skipped because not all pairs require the initial alignment.

The alternative scenario based on the exhaustive search uses the binary version of the moving and fixed images, calculated using the Li thresholding (Li and Tam 1998). Then, the centroid values are calculated for both the binary images and the desired rotation is optimized exhaustively with a fixed angle step equal to 2 degrees. The fixed step size was not tuned because the initial alignment time is negligible compared to the time required for the data loading or the nonrigid registration. The Dice coefficient is used as a similarity measure between the warped moving mask and the fixed mask.

Both the feature-based and the exhaustive initial alignment need further refinement. We noticed that the use of an additional global affine registration resistant to missing data significantly improves subsequent nonrigid registration for many cases, especially ones with significant differences in the appearance. On the other hand, the iterative intensity-based affine registration is unable to address largely misaligned cases (e.g. 180 degrees rotation), therefore it cannot be applied alone and the previous steps are mandatory. The influence of this step alone is shown in the supplementary material (is available online at stacks.iop.org/PMB/66/025006/mmedia). The additional affine registration is based on the method described by Periaswamy and Farid (2006). However, we use it as a global affine registration algorithm. The method is based on an explicit, pyramid-based, iterative optimization of 8 parameters (6 affine and 2 intensity correction parameters). During each iteration the transformation is calculated as a:

$$\begin{equation} \mathbf{T} = (\mathbf{P}\mathbf{P}^\intercal\mathbf{K})^{-1}\mathbf{P}\mathbf{K}^\intercal \end{equation} \tag{ 1 }$$

where $\mathbf{P} = [\nabla \mathbf{G_x}\cdot\mathbf{D_x}, \nabla \mathbf{G_x \cdot D_y}, \nabla \mathbf{G_y \cdot D_x}$, $\nabla \mathbf{G_y \cdot D_y}, \nabla \mathbf{G_x}, \nabla \mathbf{G_y}, -\mathbf{M}, \mathbf{1}]^\intercal$, $\mathbf{K} = \mathbf{M} - \mathbf{F} + \nabla \mathbf{G_x \cdot D_x} + \nabla \mathbf{G_y \cdot D_y}$ and $\nabla \mathbf{G_x}$, $\nabla \mathbf{G_y}$ denotes the image gradient in the x and y direction respectively, $\mathbf{D_x}$, $\mathbf{D_y}$ denotes the image domain in the x and y dimension respectively, $\mathbf{M}$ is the moving image, $\mathbf{F}$ is the fixed image and · denotes the piecewise multiplication. The image gradient is calculated using filters described by Farid and Simoncelli (2004). Then, the calculated transformation is converted to the velocity field and composed with the actual deformation field. The algorithm is run for a predefined, constant number of iterations and pyramid levels. The more detailed description of the initial alignment refinement is presented in algorithm 3 and the full initial alignment algorithm is described in algorithm 4.

Algorithm 3: Global affine registration.
Input : M (moving image), F (fixed image), $\mathbf{u_i}$ (initial deformation field—optionally),
number of resolutions, number of iterations per resolution
Output: u (the refined initial deformation field)
$\mathbf{Ms}$, $\mathbf{Fs}$ = create moving and fixed images for each resolution by resampling
$\mathbf{M}$ and $\mathbf{F}$ respectively
foreach Resolution do
forEach Iteration do
$\mathbf{M_w}$ = warp the $\mathbf{M}$ moving image with the current deformation field u
$\mathbf{T}$ = As in equation (1)
$\mathbf{v}$ = convert $\mathbf{T}$ to a velocity field
u = u ° $\mathbf{v}$
_u = upsample the deformation field u
_return u

Algorithm 4: Initial alignment summary.
Input : M (moving image), F (fixed image)
Output : u (calculated initial deformation field)
$\mathbf{Ft, Kp}$ = calculate SIFT, SURF, ORB keypoints and features using $\mathbf{M}$ and $\mathbf{F}$ images
$\mathbf{Ft, Kp}$ = filter only the possible matches between $\mathbf{Ft}$ and $\mathbf{Kp}$
$\mathbf{T_{SIFT}}, \mathbf{T_{SURF}}, \mathbf{T_{ORB}}$ = estimate the initial transformation for each feature
extraction method using RANSAC and $\mathbf{Ft, Kp}$
$\mathbf{T}$ = choose the best transformation among $\mathbf{T_{SIFT}}, \mathbf{T_{SURF}}, \mathbf{T_{ORB}}$ or detect the
failure
if Failure then
$\mathbf{T}$ = calculate the transformation using the exhaustive search between M
_and F images
u = convert the transformation $\mathbf{T}$ to the deformation field
u = refine the initial transformation u using algorithm 3
return u

After the initial alignment, a nonrigid registration is used to calculate the final, dense deformation field. The proposed method is based on the Demons algorithm (Thirion 1998) but using modality independent neighborhood (MIND) descriptor (Heinrich et al 2012) instead of the traditional unimodal approach. The sum of squared differences (SSD) between the MIND descriptor evaluated for all corresponding pixels can be used as a multimodal similarity metric. The similarity metric performance is superior compared to the mutual information and can be implemented in an efficient, parallel way (Heinrich et al 2012). We use the compositive and symmetric version of the Demons algorithm due to its fastest convergence and an ability to recover more complex deformations (Reaungamornrat et al 2160). However, we decided to not use the diffeomorphic version because we know that the real deformation field is inherently noninvertible due to the method of sample preparation. The fluid and the diffusion smoothing are applied as in the unimodal approach. The velocity field is calculated as follows:

$$\begin{equation} \begin{split} \mathbf{v} & = \frac{\nabla \mathbf{M_{MIND}} \cdot (\mathbf{M_{MIND}} - \mathbf{F_{MIND}})} {\nabla \mathbf{M_{MIND}}^2 + (\mathbf{M_{MIND}} - \mathbf{F_{MIND}})^2} \\ & \quad+ \frac{\nabla \mathbf{F_{MIND}} \cdot (\mathbf{M_{MIND}} - \mathbf{F_{MIND}})} {\nabla \mathbf{F_{MIND}}^2 + (\mathbf{M_{MIND}} - \mathbf{F_{MIND}})^2} \end{split} \end{equation} \tag{ 2 }$$

where $\mathbf{M_{MIND}}$ and $\mathbf{F_{MIND}}$ denote the MIND descriptors of the moving and the fixed image respectively. The MIND Demons algorithm is described in algorithm 5. The proposed nonrigid algorithm is dedicated to relatively small (but complex) deformations. Therefore, the correct initial alignment is crucial to obtain an accurate deformation field. The default parameters can be found in the code repository (Proposed Method Software 2019).

Algorithm 5: MIND demons algorithm.
Input : M (moving image), F (fixed image), $\sigma_\textrm{f}$ (fluid sigma), $\sigma_\textrm{d}$ (diffusion
sigma), $\sigma_\textrm{m}$ (MIND sigma), $\mathbf{r_m}$ (MIND radius), $\mathbf{u_i}$ (initial
deformation field—optionally), number of resolutions (3),
convergence indicator
Output: u (calculated deformation field)
Ms, Fs = create moving and fixed images for each resolution by resampling
M and F images
$\mathbf{F_{MINDs}}$ = calculate and cache the MIND descriptor for fixed images $\mathbf{Fs}$
u = initialize the deformation field using $\mathbf{u_i}$ or identity transform
foreach Resolution do
$\mathbf{G}(\sigma_\textrm{d})$, $\mathbf{G}(\sigma_\textrm{f})$ = initialize convolution kernels (typically Gaussian)
while Not Converged do
$\mathbf{M_w}$ = warp the moving image $\mathbf{M}$ with the current deformation field u
$\mathbf{M_{MIND}}$ = calculate the MIND descriptor for the $\mathbf{M_w}$
$\mathbf{v}$ = As in equation (2)
$\mathbf{v}$ = $\mathbf{G(\sigma_f)}$ $\star$ $\mathbf{v}$
u = u ° $\mathbf{v}$
u = $\mathbf{G}(\sigma_\textrm{d})$ $\star$ u
_u = upsample the deformation field u
_return u

Noteworthy, the original algorithm proposed for the ANHIR challenge was a bit different. The locally affine registration was used together with the MIND-based Demons registration and the better result was chosen automatically by the failure detector. However, after extending the research and improving the performance of the MIND-based Demons we decided to skip the local affine registration which improved the results significantly, taking the second instead of third place in the ANHIR competition (Borovec et al 2020, ANHIR Website 2019). It is important to emphasize that both this article, as well as the ANHIR challenge summary (Borovec et al 2020), present the results for the method described here. For a comparison to the locally affine registration we refer reader to a supplementary material.

One of the most difficult requirement was to make the method fully automatic. It was not allowed to tune parameters individually to a particular case. The methods were supposed to run as-is for all provided image pairs. As a solution to this problem, we proposed failure detection procedures after each registration stage.

For the initial alignment, after several experiments, we decided that the simplest solution is the most successful. Since the intensity histogram easily distinguishes between the image background and the tissue structure, we use adaptive thresholding for both the moving and the fixed image, and then calculate the Dice coefficient between the resulting masks. The adaptive thresholding was performed as in the initial alignment, using the Li method (Li and Tam 1998). If the resulting coefficient is lower than a given, fixed threshold or does not increase compared to the masks of images before the initial alignment, we treat the result as a failure. A similar method but with fixed parameters was used to choose the best output for the exhaustive search. The fixed threshold was chosen experimentally using the desired volume overlap ratio (0.75), without exhaustively tuning this parameter. The results are robust to changes of this parameter, as shown in the supplementary material.

The failure detection for the nonrigid registration or the affine transformation refinement is easier. We experimentally verified that it is enough to use the SSD between the MIND descriptors calculated for the warped moving image and the fixed image. Since the purpose of both the methods is to calculate relatively small deformations, the MIND descriptor can capture beneficial deformations without falling into local minima, as it would be in the initial alignment case.

The proposed algorithm is implemented using Python with standard numerical (NumPy, SciPy) and computer vision libraries (OpenCV), and plain C++ for the MIND-based Demons nonrigid registration. The implementation of the proposed method is made available at a GitHub repository (Proposed Method Software 2019). The MIND-based Demons algorithm is implemented in C++ using the OpenMP library because it is the most expensive part in terms of computational time and a naive, non-parallel Python-based implementation was too slow to take part in the ANHIR challenge. The remaining source code, since the performance is not that essential, is implemented using Python and is calculated using only a single core.

The method parameters were tuned to make the results as good as possible, however, without an exhaustive search over the whole parameter space. We chose the initial parameters using the knowledge about the algorithm behavior and then tuned them using a binary search between reasonable values. Since majority of the parameters describe how particular images, descriptors or deformation fields should be smoothed by the Gaussian filter, it is easy to determine reasonable values due to how Gaussian filtering is being implemented in practice, using the separated 1-D accumulators (the reasonable values for smoothing sigma are in the range [0.5–2.0]). The MIND radius is an integer taking just few reasonable values (from 1 to 3), therefore tuning it is even easier.

Instead of registering the pairs as a source to target (chosen by the challenge organizers) we decided to register targets to sources to avoid the necessity of inverting the calculated deformation field, thus, lowering the target registration error between the target and the warped source landmarks. During all the registration steps, all transformations, both linear and non-linear are composed together, avoiding the addition of vector fields. The moving image is never interpolated more than once to make the interpolation error negligible. The MIND descriptor for the fixed image is calculated once for each resolution and cached to significantly decrease the computation time.

Currently, the whole procedure, including loading the data, all the registration steps and saving the results, takes between 60 and 90 s using i9-9960X, which is the main drawback compared to other methods proposed during the ANHIR challenge (Borovec et al 2020). The algorithm computation time could be lowered significantly by proposing a GPU-based implementation. Since all the algorithm stages are fully parallelizable and there is no need for memory transfers, we forecast that a proper CUDA implementation on a modern GPU (e.g. GeForce RTX 2080 Ti) could result in more than 40 times speedup. This would make registration of a single pair a matter of single seconds resulting in a useful real-time tool for physicians and life scientists. It should be clarified that the processing time reported by the challenge organizers is normalized by the ANHIR evaluation system with respect to the number of cores and the processing power (Borovec et al 2020). It should be considered as processing time relative to other challenge participants, not the real time requires for the analysis, especially considering that all steps of our method, excluding the nonrigid registration, use just a single core.

The dataset consists of 481 image pairs and is split into 230 training and 251 evaluation pairs (Fernandez-Gonzalez et al 2002, Gupta et al 2018, Mikhailov et al 2018, Bueno and Deniz 2019, Borovec et al 2020). For the training pairs, both the source and the target landmarks are widely available, while for the evaluation pairs only the source landmarks can be downloaded and examined. The evaluation is performed server-side, automatically using scripts provided by the ANHIR organizers. The landmarks were annotated by 9 experts and there are on the average 86 landmarks per image. The average error between the same landmarks chosen by two annotators is 0.05% of the image size which can be used as the indicator of the best possible results to achieve by the registration methods, below which the results become indistinguishable (Borovec et al 2018, 2020).

The tissues are provided as RGB images in a .png format, without the histological metadata. The tissues were stained using: (i) clara cell 10 protein, (ii) prosurfactant protein C, (iii) hematoxylin and eosin, (iv) anigen KI-67, (v) platelet endothelial cell adhesion molecule, (vi) human epidermal growth factor receptor 2, (vii) estrogen receptor, (viii) progesterone receptor, (ix) cytokeratin, (x) podocin. The images are divided into 8 classes containing: (i) lung lesion tissue, (ii) lung lobes, (iii) mammary glands, (iv) the colon adenocarcinomas (COADs), (v) mice kidney tissues, (vi) gastric mucosa and adenocarcinomas tissues, (vii) breast tissues, (viii) human kidney tissues (Borovec et al 2020, ANHIR Website 2019). All the tissues were acquired in different acquisition settings making the dataset more diverse. Full description of the dataset can be found in ANHIR Website (2019). An exemplary visualization of images are shown in figure 2. The challenge was performed on medium size images which were approximately 25% of the original image scale. The approximate image size after the resampling varied from 8k to 16k pixels (in one dimension).

Zoom In Zoom Out Reset image size

The target registration error (TRE) was the evaluation metric used for the ANHIR challenge, more concretely a rTRE (dimensionless) defined as:

$$\begin{equation} \textrm{rTRE} = \frac{\textrm{TRE}}{\sqrt{w^2 + h^2}} , \end{equation} \tag{ 3 }$$

where TRE denotes the target registration error, w is the image width and h is the image height. The normalization by the image diagonal made it possible to compare TREs between image pairs with different resolutions. The challenge organizers reported following metrics related to the rTRE: (i) average median rTRE, (ii) average rank of median rTRE (the main evaluation metric), (iii) average max rTRE, (iv) average rank of max rTRE. The median rTRE describes the quality of the final nonrigid registration. The average rTRE shows the initial alignment quality and the resistance to failures. The max rTRE presents the ability to recover large deformations. In this article, we focus on the main ANHIR evaluation metric: the average rank of median rTRE. For a full summary of the challenge results, we refer to Borovec et al (2020).