Further author information: (Send correspondence to M.G.)
M.G.: Email<EMAIL_ADDRESS>
# An End-to-end Deep Learning Approach for Landmark Detection and Matching in
Medical Images

Monika Grewal  Life Sciences & Health Research Group, Centrum Wiskunde &
Informatica, P.O. Box 94079, 1090 GB Amsterdam, The Netherlands  Timo M. Deist
Life Sciences & Health Research Group, Centrum Wiskunde & Informatica, P.O.
Box 94079, 1090 GB Amsterdam, The Netherlands  Jan Wiersma  Department of
Radiation Oncology, Amsterdam UMC, University of Amsterdam, P.O. Box 22660,
1100 DD Amsterdam, The Netherlands  Peter A. N. Bosman  Life Sciences & Health
Research Group, Centrum Wiskunde & Informatica, P.O. Box 94079, 1090 GB
Amsterdam, The Netherlands  Faculty of Electrical Engineering, Mathematics and
Computer Science, Delft University of Technology, P.O. Box 5, 2600 AA Delft,
The Netherlands  Tanja Alderliesten

###### Abstract

Anatomical landmark correspondences in medical images can provide additional
guidance information for the alignment of two images, which, in turn, is
crucial for many medical applications. However, manual landmark annotation is
labor-intensive. Therefore, we propose an end-to-end deep learning approach to
automatically detect landmark correspondences in pairs of two-dimensional (2D)
images. Our approach consists of a Siamese neural network, which is trained to
identify salient locations in images as landmarks and predict matching
probabilities for landmark pairs from two different images. We trained our
approach on 2D transverse slices from 168 lower abdominal Computed Tomography
(CT) scans. We tested the approach on 22,206 pairs of 2D slices with varying
levels of intensity, affine, and elastic transformations. The proposed
approach finds an average of 639, 466, and 370 landmark matches per image pair
for intensity, affine, and elastic transformations, respectively, with spatial
matching errors of at most 1 mm. Further, more than 99% of the landmark pairs
are within a spatial matching error of 2 mm, 4 mm, and 8 mm for image pairs
with intensity, affine, and elastic transformations, respectively. To
investigate the utility of our developed approach in a clinical setting, we
also tested our approach on pairs of transverse slices selected from follow-up
CT scans of three patients. Visual inspection of the results revealed landmark
matches in both bony anatomical regions as well as in soft tissues lacking
prominent intensity gradients.

###### keywords:

end-to-end, landmark detection, CT, deep learning, deformable image
registration

##  1 INTRODUCTION

Deformable Image Registration (DIR) can be extremely valuable in work-flows
related to image-guided diagnostics and treatment planning. However, DIR in
medical imaging can be challenging due to large anatomical variations between
images. This is particularly the case in the lower abdomen, where internal
structures can undergo large deformations between two scans of a patient due
to physical conditions like presence of gas pockets and bladder filling. Such
scenarios are particularly challenging for intensity based registration, as
there are many local optima to overcome. Landmark correspondences between
images can provide additional guidance information to the DIR methods[1, 2]
and increase the probabilty of finding the right transformation by adding
landmark matches as an additional constraint or objective in the optimization.
Since the manual annotation of anatomical landmarks is labor-intensive and
requires expertise, developing methods for finding landmark correspondences
automatically has great potential benefits.

The existing methods[3, 4, 5, 6, 7] for obtaining landmark correspondences in
medical images are based on large and time-consuming pipelines that involve
identifying landmark locations followed by matching local feature
descriptors[8] within a restricted neighborhood. These methods rely upon
multiple pre- and post-processing steps, multi-resolution search, and manual
checking to achieve robustness; each step adding more heuristics and empirical
hyperparameters to an already complex pipeline. Further, existing methods for
landmark detection that restrict the definition of landmarks to certain
intensity gradient patterns specific to the underlying data set or anatomical
region may not be easily adaptable to other contexts [9]. Generalizing the
definition of landmarks and reducing the number of heuristics would allow for
faster adaptation of automated methods for different clinical settings. In
addition, faster execution times for landmark detection and matching could
benefit their clinical application.

Recently, deep Convolutional Neural Networks (CNNs) have shown promising
results for classification and segmentation tasks in medical imaging due to
their capability of learning discriminant feature descriptors from raw images
[10, 11, 12]. There exist a few deep learning approaches for finding landmarks
in medical images [13, 14]. However, in these approaches a neural network is
trained in a supervised manner to learn a small number of manually annotated
landmarks. It is to be noted that a high density of landmark correspondences
is desirable to effectively provide additional guidance to the DIR methods. In
a supervised setting, it means annotating thousands of landmarks per CT scan,
which is intractable in terms of required manual efforts. On the other hand,
many deep learning approaches have been developed for automatically finding
object landmarks in natural images [15, 16, 17, 18] that do not require manual
annotations. Some of these approaches focus on discovering a limited number of
landmarks in an image dataset. Whereas, others either fine-tune a pre-trained
network or make use of incremental training in a self-supervised fashion.

Our proposed approach is based on the above-mentioned approaches developed for
natural images and tailored to meet the specific requirements relating to the
medical images. We propose a two-headed Siamese neural network that based on a
pair of images simultaneously predicts the landmarks and their feature
descriptors corresponding to each image. These are then sent to another module
to predict their matching probabilities. We train the neural network from
scratch and gradients are back-propagated from end-to-end. To the best of our
knowledge, this is first endeavour to develop an end-to-end deep learning
approach for finding landmark correspondences in medical images. Our approach
has the following distinct advantages compared to existing methods for finding
landmark correspondences:

  * ·

Our approach is end-to-end deep learning based; therefore, the need for data
pre- and post-processing during inference is avoided. In addition, the
proposed approach is faster at run-time and has fewer hyperparameters than
traditional approaches.

  * ·

We do not impose any prior on the definition of a landmark in an image.
Instead, we train the network in a way that the landmarks represent salient
regions in the image that can be found repeatedly despite potential intensity
variations, and deformations.

  * ·

The proposed approach does not require manual annotations for training and
learns from data in a self-supervised manner.

  * ·

Our approach improves over the existing approaches for natural images by
avoiding the need for pre-training, or incremental fine-tuning of the neural
network.

##  2 MATERIALS AND METHODS

###  2.1 Data

In total 222 lower abdominal Computed Tomography (CT) scans of female patients
acquired for radiation treatment planning purposes were retrospectively
included: 168 scans (24,923 two-dimensional (2D) slices) were used for
training and 54 scans (7,402 2D slices) were used for testing. For a separate
set of three patients, one original scan along with a follow-up CT scan was
included. The scans of these three patients were used for testing the approach
in a clinical setting. All CT scans had an in-plane resolution from 0.91 mm
$\times$ 0.91 mm to 1.31 mm $\times$ 1.31 mm. All the 2D slices were resampled
to 1 mm $\times$ 1 mm in-plane resolution.

###  2.2 Approach

Figure 1:  Schematic representation of our approach. The weights are shared
between two branches of the Siamese neural network. The transformation is
required only during training for calculating the ground truths. Abbreviations
of the data input and output at various stages follow the description in the
text.

In Figure 1, the different modules of our approach are illustrated along with
the data flow between them. Our approach comprises a Siamese architecture
consisting of CNN branches with shared weights. The outputs of the CNN
branches are sent to a module named Sampling Layer followed by another module
named Feature Descriptor Matching Module. The network takes two images $I_{1}$
and $I_{2}$ as inputs and predicts $K_{1}$ and $K_{2}$ landmarks in $I_{1}$
and $I_{2}$, respectively. In addition, the network predicts matching
probabilities ($\hat{c}_{i,j}$) for each landmark $i\in\\{1,2,...,K_{1}\\}$ in
$I_{1}$ to a landmark $j\in\\{1,2,...,K_{2}\\}$ in $I_{2}$. In the following
paragraphs, a description of each module is provided.

####  2.2.1 CNN branches

The CNN branches of the Siamese neural network have shared weights and consist
of an encoder-decoder type network similar to the U-Net[10] architecture. The
only difference from the original implementation is that the number of
convolutional filters in each layer is reduced by a factor of four to avoid
overfitting. The implemented architecture contains 16, 32, 64, 128, and 256
convolutional filters in successive downsampling blocks respectively. The CNN
branches give two outputs for each input image: a landmark probability map,
and feature descriptors. The landmark probability map is computed at the end
of the upsampling path after applying the sigmoid non-linearity and the
feature descriptors are computed by concatenation of feature maps from the
last two downsampling blocks. The feature maps from different downsampling
blocks intrinsically allow for feature matching at multiple resolutions and
abstraction levels.

####  2.2.2 Sampling Layer

The sampling layer is a parameter-free module of the network. It performs the
following tasks:

  1. 1.

It samples $K_{1}$ and $K_{2}$ landmark locations in $I_{1}$ and $I_{2}$,
respectively, which correspond to the highest probability score locations in
the predicted landmark probability maps.

  2. 2.

It extracts predicted landmark probabilities $\hat{p}^{I_{1}}_{i}$, and
$\hat{p}^{I_{2}}_{j}$ corresponding to $K_{1}$ and $K_{2}$ locations in
landmark probability maps of image $I_{1}$ and $I_{2}$.

  3. 3.

It extracts feature descriptors $f^{I_{1}}_{i}$ and $f^{I_{2}}_{j}$
corresponding to the sampled landmark locations in $I_{1}$ and $I_{2}$,
respectively, and creates feature descriptor pairs
$(f^{I_{1}}_{i},f^{I_{2}}_{j})$ for each $i\in\\{1,2,...,K_{1}\\}$ and
$j\in\\{1,2,...,K_{2}\\}$.

  4. 4.

During training, it generates the ground truths for landmark probabilities and
feature descriptor matching probabilities on-the-fly as mentioned in Georgakis
et al [17]. Briefly, the sampled landmark locations of $I_{2}$ are projected
onto $I_{1}$ based on the known transformation between the images. A landmark
location $i$ in $I_{1}$ is decided to be matching to a landmark location $j$
in $I_{2}$ if the Euclidean distance between $i$ and the projection of $j$ on
image $I_{1}$ is less than a predefined pixel threshold ($thresh_{pixels}$).

####  2.2.3 Feature Descriptor Matching Module

All the feature descriptor pairs $(f^{I_{1}}_{i},f^{I_{2}}_{j})$ are fed to
the feature descriptor matching module. The feature descriptor matching module
consists of a single fully connected layer that predicts the matching
probability for each feature descriptor pair.

###  2.3 Training

Training image pairs were generated on-the-fly by sampling a reference image
randomly and generating the target image by transforming the reference image
with a known transformation (randomly simulated brightness or contrast jitter,
rotation, scaling, shearing, or elastic transformation). During training, the
ground truths for landmark probabilities and feature descriptor matching
probabilities are generated in the sampling layer as described above. We
trained the network by minimizing a multi-task loss defined as follows:

$Loss={LandmarkProbabilityLoss}_{I_{1}}+{LandmarkProbabilityLoss}_{I_{2}}+{DescriptorMatchingLoss}$
(1)

The $LandmarkProbabilityLoss_{I_{n}}$ for the probabilities of landmarks in
image $I_{n},n\in\\{1,2\\}$ is defined as:

${LandmarkProbabilityLoss}_{I_{n}}=\frac{1}{K_{n}}\sum_{i=1}^{K_{n}}\left((1-\hat{p}^{I_{n}}_{i})+CrossEntropy(\hat{p}^{I_{n}}_{i},p^{I_{n}}_{i})\right)$
(2)

where $CrossEntropy$ is the cross entropy loss between predicted landmark
probabilities $\hat{p}^{I_{n}}_{i}$ and ground truths $p^{I_{n}}_{i}$. The
term $(1-\hat{p}^{I_{n}}_{i})$ in (2) encourages high probability scores at
all the sampled landmark locations, whereas the cross entropy loss term forces
low probability scores at the landmark locations that do not have a
correspondence in the other image. As a consequence, the network is forced to
predict high landmark probabilities only at the salient locations that have
correspondence in the other image as well.

Hinge loss is widely used for learning discriminant landmark descriptors
between matching and non-matching landmark pairs. We observed that a positive
margin for the matching pairs in the hinge loss encourages the network to
focus on hard positive examples (i.e., non-trivial landmark matches).
Therefore, we defined $DescriptorMatchingLoss$ (equation 3) as a linear
combination of hinge loss with a positive margin $m_{pos}$ on the L2-norm of
feature descriptor pairs and cross entropy loss on matching probabilities
predicted by the feature descriptor matching module.

$\begin{split}DescriptorMatchingLoss&=\sum_{i=1,j=1}^{K_{1},K_{2}}\left(\frac{c_{i,j}max(0,||f^{I_{1}}_{i}-f^{I_{2}}_{j}||^{2}-m_{pos})}{K_{pos}}\right.\\\
&+\frac{(1-c_{i,j})max(0,m_{neg}-||f^{I_{1}}_{i}-f^{I_{2}}_{j}||^{2})}{K_{neg}}\\\
&+\left.\frac{WeightedCrossEntropy(\hat{c}_{i,j},c_{i,j})}{(K_{pos}+K_{neg})}\right)\end{split}$
(3)

where $\hat{c}_{i,j}$, and $c_{i,j}$ are the predicted and the ground truth
matching probabilities, respectively, for the feature descriptor pair
$(f^{I_{1}}_{i},f^{I_{2}}_{j})$; $K_{pos}$ and $K_{neg}$ are the number of
matching (positive class) and non-matching (negative class) feature descriptor
pairs; $m_{pos}$ and $m_{neg}$ are the margins for the L2-norm of matching and
non-matching feature descriptor pairs. $WeightedCrossEntropy$ is the binary
cross entropy loss where the loss corresponding to positive class is weighted
by the frequency of negative examples and vice versa. The gradients are back-
propagated from end-to-end as indicated by the dashed arrows in Figure 1.

###  2.4 Constraining Landmark Locations

A naive implementation of the approach may find all the landmarks clustered in
a single anatomical region, which is not desirable. Therefore, to learn
landmarks in all anatomical regions during training, we sample the landmarks
on a coarse grid in the sampling layer, i.e., in each $8\times 8$ pixel
section of the grid, only one landmark location with the maximum landmark
probability is sampled.

Another challenge in the CT scan imaging data comes from a large number of
pixels belonging to the background. Traditionally, the image is cropped to the
center to avoid prediction of landmarks in the background or on the patient
table. However, this strategy requires an additional pre-processing step
during inference. To avoid this, we computed a valid mask for each image,
which contained the value 1 at the location of body pixels and 0 elsewhere.
The valid mask was generated by image binarization using intensity
thresholding and removing small connected components in the binarized image.
The network is trained to predict high landmark probabilities as well as
feature descriptor matching probabilities only in the matching locations that
correspond to a value of 1 in the valid mask. This allows the network to learn
a content-based prior on the landmark locations and avoids the need for image
pre-processing during inference.

###  2.5 Inference

During inference, only the locations in $I_{1}$ and $I_{2}$ with landmark
probabilities above a threshold ($thresh_{landmark}$) are considered. Further,
landmark pairs from different images are only matched if their matching is
inverse consistent. Suppose, locations $i\in\\{1,..,K_{1}\\}$ in $I_{1}$ and
locations $j\in\\{1,..,K_{2}\\}$ in $I_{2}$ have landmark probabilities above
$thresh_{landmark}$. A pair $(i^{\ast},j^{\ast})$ is considered matching if
there is no other pair $(i^{\ast},j^{\prime})$ where
$j^{\prime}\in\\{1,..,K_{2}\\}$ or $(i^{\prime},j^{\ast})$ where
$i^{\prime}\in\\{1,..,K_{1}\\}$ with higher descriptor matching probabilities
or lower L2-norms for their feature descriptor pairs
$(f_{i^{\ast}}^{I_{1}},f_{j^{\prime}}^{I_{2}})$ or
$(f_{i^{\prime}}^{I_{1}},f_{j^{\ast}}^{I_{2}})$.

###  2.6 Implementation Details

We implemented our approach using PyTorch[19]. We trained the network for 50
epochs using the Adam[20] optimizer with learning rate $10^{-3}$ and a weight
decay of $10^{-4}$. The training was done with a batchsize of 4 and took 28
GPU (NVIDIA GeForce RTX 2080 Ti) hours. To allow for batching, a constant $K$
(set to 400) landmarks were sampled from all the images. The threshold for
Euclidean distance while generating the ground truth ($thresh_{pixels}$) was 2
pixels. The margin for the L2-norm of matching feature descriptors ($m_{pos}$)
was set to 0.1 and the margin for the L2-norm of non-matching pairs
($m_{neg}$) was set to 1. During inference, $thresh_{landmark}$ = 0.5 was
used.

The empirical values for the hyperparameters were decided based on experience
in the preliminary experiments. For example, the number for landmarks to be
sampled during training ($K$) was decided such that the entire image was
covered with sufficient landmark density, which was inspected visually.
Similarly, the decision for $thresh_{pixels}$ was motivated by the fact that a
threshold less than 2 pixels did not yield any matching landmarks in the first
few iterations of the training and hence the network could not be trained. We
initially trained the network with default values of $m_{pos}$, and $m_{neg}$
($m_{pos}=0$, and $m_{neg}=1$). However, we noticed on the validation set that
all the predicted landmark pairs were clustered in regions of no deformation.
To avoid this behaviour, we trained the network with $m_{pos}=0.1$ and
$m_{pos}=0.2$ so that the gradients were not affected by the hinge loss
corresponding to easy landmark matches. The final results are reported
corresponding to the run with $m_{pos}=0.1$ as it had a better trade off
between number of landmarks per image pair and difficulty of landmark
locations. The value of $thresh_{landmark}$ was chosen to give the best trade
off between the number of landmarks per image pair and the spatial matching
error on the validation set.

##  3 Experiments

###  3.1 Baseline

Scale Invariant Feature Transform (SIFT[21]) based keypoint detectors and
feature descriptors are prevalent approaches used in both natural image
analysis as well as in medical image analysis [6]. Therefore, we used the
OpenCV[22] implementation of SIFT as the baseline approach for comparison. We
used two matching strategies for SIFT: a) brute-force matching with inverse
consistency (similar to our approach, we refer to this approach as SIFT-
InverseConsistency), b) brute-force matching with ratio test (as described in
the original paper[21], we refer to this approach as SIFT-RatioTest). Default
values provided in the OpenCV implementation were used for all other
hyperparameters.

###  3.2 Datasets

The performance is evaluated on two test sets. First, for quantitative
evaluation, we transformed all 7,402 testing images from 54 CT scans with
three different types of transformations corresponding to intensity (jitter in
pixel intensities = $\pm 20\%$ maximum intensity), affine (pixel displacement:
median = 29 mm, Inter Quartile Range (IQR) = 14 mm - 51 mm), and elastic
transformations (pixel displacement: median = 12 mm, IQR = 9 mm - 15 mm),
respectively. Elastic transformations were generated by deforming the original
image according to a deformation vector field representing randomly-generated
2D Gaussian deformations. The extent of transformations was decided such that
the intensity variations and the displacement of pixels represented the
typical variations in thoracic and abdominal CT scan images [23, 24]. This
resulted in three sets of 7,402 2D image pairs (total 22,206 pairs).

Second, to test the generalizability of our approach in a clinical setting,
image pairs were taken from two CT scans of the same patient but acquired on
different days. The two scans were aligned with each other using affine
registration in the SimpleITK [25] package. This process was repeated for
three patients.

###  3.3 Evaluation

For quantitative evaluation, we projected the predicted landmarks in the
target images to the reference images and calculated the Euclidean distance to
their corresponding matches in the reference images. We report the cumulative
distribution of landmark pairs with respect to the Euclidean distance between
them.

The performance of our approach on clinical data was assessed visually. We
show the predicted results on four transverse slices belonging to different
anatomical regions. To visually trace the predicted correspondences of
landmarks, the colors of the landmarks in both the images vary according to
their location in the original CT slice. Similarly colored dots between slices
from original and follow-up image represent matched landmarks.

##  4 RESULTS

Table 1:  Description of predicted landmark matches. Median number of landmark matches per image pair with Inter Quartile Range (IQR) in parentheses are provided together with the spatial matching error. The entries in bold represent the best value among all approaches. Transformations |  Intensity |  Affine |  Elastic
---|---|---|---
No. of landmarks |  Proposed Approach |  639 (547 - 729) |  466 (391 - 555) |  370 (293 - 452)
SIFT - InverseConsistency |  711 (594 - 862) |  610 (509 - 749) |  542 (450 - 670)
SIFT - RatioTest |  698 (578 - 849) |  520 (426 - 663) |  418 (330 - 541)
Spatial matching error (mm) |  Proposed Approach |  0.0 (0.0 - 0.0) |  1.0 (0.0 - 1.4) |  1.0 (1.0 - 1.4)
SIFT - InverseConsistency |  1.0 (1.0 - 1.4) |  1.0 (1.0 - 1.4) |  1.0 (1.0 - 2.0)
SIFT - RatioTest |  1.0 (1.0 - 1.4) |  1.0 (1.0 - 1.4) |  1.0 (1.0 - 1.4)

The inference time of our approach per 2D image pair is within 10 seconds on a
modern CPU without any parallelization. On the GPU the inference time is
$\sim$20 milliseconds. The model predicted on average 639 (IQR = 547 - 729),
466 (IQR = 391 - 555), and 370 (IQR = 293 - 452) landmark matches per image
pair for intensity, affine, and elastic transformations, respectively.

###  4.1 Simulated Transformations

Table 1 describes the number of landmark matches per image pair and the
spatial matching error for both our approach and the two variants of SIFT.
Though our approach finds less landmarks per image as compared to the two
variants of SIFT, the predicted landmarks have smaller spatial matching error
than the SIFT variants. Further, Figure 2 shows the cumulative distribution of
landmark pairs with respect to the Euclidean distance between them. All the
approaches are able to find more than 90% of landmark matches within 2 mm
error for intensity transformations. Predicting landmark correspondences under
affine and elastic transformations is considerably more difficult; this can
also be seen in the worse performance of all approaches. However, our approach
is still able to find more than 99% of landmark matches within a spatial
matching error of 4 mm and 8 mm, respectively for affine and elastic
transformations. However, a noticeable percentage (about 2% for affine
transformations and 3% for elastic transformations) of landmarks detected by
SIFT-RatioTest are wrongly matched with landmarks from far apart regions (more
than 64 mm). It should be noted that if landmark matches with such high
inaccuracies are used for providing guidance to a registration method, it may
have a deteriorating effect on the registration if the optimizer is not
sufficiently regularized.

Figure 2:  Cumulative distribution of landmarks. The cross-hairs in (b) and
(c) correspond to the percentile of landmarks in SIFT-RatioTest at 64 mm.

For visual comparison, the landmark correspondences in pairs of original and
elastic transformed images are shown in Figure 3 (rows a-b) for our approach
as well as for SIFT. As can be seen, the cases of mismatch in predictions from
our approach (i.e., the number of landmarks in transformed slices not
following the color gradient in the original slice) are rather scarce in
comparison to the baseline approaches. Another interesting point to note is
the difference in the landmark locations from our approach and the two
baseline approaches. Since SIFT is designed to predict landmarks at locations
of local extrema, the landmark matches are concentrated on the edges in the
images. Our approach, however, predicts matches in soft tissue regions as
well. Further inspection reveals that our approach predicts a considerable
number of landmark matches even in the deformed regions in contrast to the
baseline approaches. The capability to establish landmark correspondences in
the soft tissues and deformed regions is important because DIR methods can
especially benefit from guidance information in these regions.

###  4.2 Clinical Transformations

Rows c-f in Figure 3 show landmark correspondences in pairs of transverse
slices corresponding to the lower abdominal region in the original and follow-
up CT for our approach as well as for SIFT. As can be seen, the original and
follow-up slices have large differences in local appearance of structures
owing to contrast agent, bladder filling, presence or absence of gas pockets,
which was not part of the training procedure. It is notable that the model is
able to find considerable landmark matches in image pairs despite these
changes in local appearance. Moreover, the spatial matching error of landmarks
seems similar to that of images with simulated transformations, in contrast to
the baseline approach SIFT-InverseConsistency. Further, SIFT-RatioTest
predicts fewer mismatched landmarks compared to SIFT-InverseConsistency, but
this is achieved at the cost of a large decrease in the number of landmark
matches per image pair.

Figure 3:  Landmark correspondences for pairs of different transverse slices
in abdominal CT scans. The landmark correspondences predicted by our approach
are shown in comparison with two variants of SIFT. Rows (a-b) show predictions
on pairs of original (left) and elastic transformed (right) slices. Rows (c-f)
show transverse slices taken from different anatomical regions. The slices in
the original CT (left) are matched with a similar slice from a follow-up CT
scan (right) by affine registration.

##  5 DISCUSSION AND CONCLUSIONS

With a motivation to provide additional guidance information for DIR methods
of medical images, we developed an end-to-end deep learning approach for the
detection and matching of landmarks in an image pair. To the best of our
knowledge, this is the first approach that simultaneously learns landmark
locations as well as the feature descriptors for establishing landmarks
correspondences in medical imaging. While the final version of this manuscript
was being prepared, we came across one research on retinal images [26], whose
approach for landmark detection using UNet architecture in a semi-supervised
manner is partly similar to ours. However, our approach not only learns the
landmark locations, but also the feature descriptors and the feature matching
such that the entire pipeline for finding landmark correspondences can be
replaced by a neural network. Therefore, our approach can be seen as an
essential extension to the mentioned approach.

Our proposed approach does not require any expert annotation or prior
knowledge regarding the appearance of landmarks in the learning process.
Instead, it learns landmarks based on their distinctiveness in feature space
despite local transformations. Such a definition of landmarks is generic so as
to be applicable in any type of image and sufficient for the underlying
application of establishing correspondences between image pairs. Further, in
contrast to the traditional unsupervised approaches for landmark detection in
medical imaging, the proposed approach does not require any pre- or post-
processing steps, and has fewer hyperparameters.

The main challenge for intensity based DIR methods is to overcome local optima
caused by multiple low contrast regions in the image, which result in image
folding and unrealistic transformations in the registered image. It can be
speculated that the availability of landmark correspondences in the low
contrast image regions may prove to be beneficial for DIR methods. Moreover, a
uniform coverage of entire image is desirable for improved performance. Upon
visual inspection of the landmarks predicted by our approach, we observed that
our approach not only finds landmark correspondences in bony anatomical
regions but also in soft tissue regions lacking intensity gradients. Moreover,
a considerable density of landmarks (approximately 400 landmarks per image
pair) was observed despite the presence of intensity, affine, or elastic
transformations. Based on these observations, we are optimistic about the
potential added value of our approach to the DIR methods.

We validated our approach on images with simulated intensity, affine, and
elastic transformations. The quantitative results show low spatial matching
error of the landmarks predicted by our approach. Additionally, the results on
clinical data demonstrate the generalization capability of our approach. We
compared the performance of our approach with the two variants of widely used
SIFT keypoint detection approach. Our approach not only outperforms the SIFT
based approach in terms of matching error under simulated transformations, but
also finds more accurate matches in the clinical data. As such the results
look quite promising. However, the current approach is developed for 2D images
i.e., it overlooks the possibility of the out-of-plane correspondences in two
CT scans, which is quite likely especially in lower abdominal regions. The
extension of the approach to 3D is, therefore, imperative so as to speculate
into its benefits in providing additional guidance information to the DIR
methods.

##  6 ACKNOWLEDGEMENTS

The research is part of the research programme, Open Technology Programme with
project number 15586, which is financed by the Dutch Research Council (NWO),
Elekta, and Xomnia. Further, the work is co-funded by the public-private
partnership allowance for top consortia for knowledge and innovation (TKIs)
from the Ministry of Economic Affairs.

## References

  * [1] Alderliesten, T., Bosman, P. A. N., and Bel, A., “Getting the most out of additional guidance information in deformable image registration by leveraging multi-objective optimization,” in [Medical Imaging 2015: Image Processing ], Proc. SPIE 9413, 94131R, International Society for Optics and Photonics (2015).
  * [2] Han, D., Gao, Y., Wu, G., Yap, P.-T., and Shen, D., “Robust anatomical landmark detection with application to MR brain image registration,” Comput. Med. Imaging Graph. 46, 277–290 (2015).
  * [3] Yang, D., Zhang, M., Chang, X., Fu, Y., Liu, S., Li, H. H., Mutic, S., and Duan, Y., “A method to detect landmark pairs accurately between intra-patient volumetric medical images,” Med. Phys. 44(11), 5859–5872 (2017).
  * [4] Werner, R., Duscha, C., Schmidt-Richberg, A., Ehrhardt, J., and Handels, H., “Assessing accuracy of non-linear registration in 4D image data using automatically detected landmark correspondences,” in [Medical Imaging 2013: Image Processing ], Proc. SPIE 8669, 86690Z, International Society for Optics and Photonics (2013).
  * [5] Rühaak, J., Polzin, T., Heldmann, S., Simpson, I. J. A., Handels, H., Modersitzki, J., and Heinrich, M. P., “Estimation of large motion in lung CT by integrating regularized keypoint correspondences into dense deformable registration,” IEEE Trans. Med. Imaging 36(8), 1746–1757 (2017).
  * [6] Ghassabi, Z., Shanbehzadeh, J., Sedaghat, A., and Fatemizadeh, E., “An efficient approach for robust multimodal retinal image registration based on UR-SIFT features and PIIFD descriptors,” EURASIP J. Image Video Process. 2013(1), 25 (2013).
  * [7] Chen, J., Tian, J., Lee, N., Zheng, J., Smith, R. T., and Laine, A. F., “A partial intensity invariant feature descriptor for multimodal retinal image registration,” IEEE Trans. Biomed. Eng. 57(7), 1707–1718 (2010).
  * [8] Guo, Y., Bennamoun, M., Sohel, F., Lu, M., Wan, J., and Kwok, N. M., “A comprehensive performance evaluation of 3D local feature descriptors,” Int. J. Comput. Vis. 116(1), 66–89 (2016).
  * [9] Hervella, Á. S., Rouco, J., Novo, J., and Ortega, M., “Multimodal registration of retinal images using domain-specific landmarks and vessel enhancement,” Procedia Comput. Sci. 126, 97–104 (2018).
  * [10] Ronneberger, O., Fischer, P., and Brox, T., “U-Net: Convolutional networks for biomedical image segmentation,” in [International Conference on Medical Image Computing and Computer-Assisted Intervention ], 234–241, Springer (2015).
  * [11] Gulshan, V., Peng, L., Coram, M., Stumpe, M. C., Wu, D., Narayanaswamy, A., Venugopalan, S., Widner, K., Madams, T., Cuadros, J., Kim, R., Raman, R., Nelson, P. C., Mega, J. L., and Webster, D. R., “Development and validation of a deep learning algorithm for detection of diabetic retinopathy in retinal fundus photographs,” JAMA Netw. Open 316(22), 2402–2410 (2016).
  * [12] Esteva, A., Kuprel, B., Novoa, R. A., Ko, J., Swetter, S. M., Blau, H. M., and Thrun, S., “Dermatologist-level classification of skin cancer with deep neural networks,” Nature 542(7639), 115 (2017).
  * [13] Bier, B., Unberath, M., Zaech, J.-N., Fotouhi, J., Armand, M., Osgood, G., Navab, N., and Maier, A., “X-ray-transform invariant anatomical landmark detection for pelvic trauma surgery,” in [International Conference on Medical Image Computing and Computer-Assisted Intervention ], 55–63, Springer (2018).
  * [14] Tuysuzoglu, A., Tan, J., Eissa, K., Kiraly, A. P., Diallo, M., and Kamen, A., “Deep adversarial context-aware landmark detection for ultrasound imaging,” in [International Conference on Medical Image Computing and Computer-Assisted Intervention ], 151–158, Springer (2018).
  * [15] Thewlis, J., Bilen, H., and Vedaldi, A., “Unsupervised learning of object landmarks by factorized spatial embeddings,” in [The IEEE International Conference on Computer Vision ], 5916–5925 (2017).
  * [16] Zhang, Y., Guo, Y., Jin, Y., Luo, Y., He, Z., and Lee, H., “Unsupervised discovery of object landmarks as structural representations,” in [Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition ], 2694–2703 (2018).
  * [17] Georgakis, G., Karanam, S., Wu, Z., Ernst, J., and Košecká, J., “End-to-end learning of keypoint detector and descriptor for pose invariant 3D matching,” in [Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition ], 1965–1973 (2018).
  * [18] DeTone, D., Malisiewicz, T., and Rabinovich, A., “Superpoint: Self-supervised interest point detection and description,” in [Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition Workshops ], 224–236 (2018).
  * [19] Paszke, A., Gross, S., Chintala, S., Chanan, G., Yang, E., DeVito, Z., Lin, Z., Desmaison, A., Antiga, L., and Lerer, A., “Automatic differentiation in PyTorch,” in [Advances in Neural Information Processing Systems-W ], (2017).
  * [20] Kingma, D. P. and Ba, J., “Adam: A method for stochastic optimization,” arXiv preprint arXiv:1412.6980 (2014).
  * [21] Lowe, D. G., “Distinctive image features from scale-invariant keypoints,” Int. J. Comput. Vis. 60(2), 91–110 (2004).
  * [22] Bradski, G., “The OpenCV library,” Dr. Dobb’s J. 25, 120–125 (2000).
  * [23] Vásquez Osorio, E. M., Hoogeman, M. S., Méndez Romero, A., Wielopolski, P., Zolnay, A., and Heijmen, B. J. M., “Accurate CT/MR vessel-guided nonrigid registration of largely deformed livers,” Med. Phys. 39(5), 2463–2477 (2012).
  * [24] Polzin, T., Rühaak, J., Werner, R., Strehlow, J., Heldmann, S., Handels, H., and Modersitzki, J., “Combining automatic landmark detection and variational methods for lung CT registration,” in [Fifth International Workshop on Pulmonary Image Analysis ], 85–96 (2013).
  * [25] Lowekamp, B. C., Chen, D. T., Ibáñez, L., and Blezek, D., “The design of SimpleITK,” Frontiers in Neuroinformatics 7, 45 (2013).
  * [26] Truong, P., Apostolopoulos, S., Mosinska, A., Stucky, S., Ciller, C., and Zanet, S. D., “Glampoints: Greedily learned accurate match points,” in [Proceedings of the IEEE International Conference on Computer Vision ], 10732–10741 (2019).