Unsupervised Landmark Detection Based Spatiotemporal Motion Estimation for 4D Dynamic Medical Images

Image registration is a basic tool in motion estimation, which can be mainly divided into two categories, i.e., non-learning-based methods and learning-based methods. Classic non-learning-based approaches mainly include diffeomorphic Demons [14], Free-Form Deformation (FFD) [27], etc. These methods derive the deformation field by iteratively maximizing the image similarity under certain smoothness regularization. Both methods are widely used in medical image registration, and are often regarded as the comparing baselines in related researches. However, the hyper-parameters of them are not easy to tune, while tremendous optimization costs high in computation time.

Inspired by the great success of deep learning, several learning-based studies have been proposed in recent years. Many of them leverage the supervised learning strategy, which requires the ground-truth (deformation field) derived from conventional registration tools to train the model. For instance, Yang et al. [28] utilized a U-Net architecture [29] to predict the deformation field for brain MR images. Similarly, Rohe et al. [30] also leveraged U-Net to generate the motion field for cardiac volumes. Krebs et al. [31] present a supervised agent learing based non-rigid registration method to learn the deformation field in low-dimension representation. However, the supervision can be cumbersome to acquire, which restricts the deformation types to be learned.

Recently, unsupervised learning solutions for registration have been proposed to eliminate high reliance on the ground-truth. Li et al. [32] and Vos et al. [33] introduced fully convolutional networks to estimate the deformation field for brain volumetric images and 4D cardiac images, respectively. Guha et al. [34] developed the popular VoxelMorph method for brain image registration. Zhang et al. [35] proposed an inverse-consistency loss based on a similar network of VoxelMorph. These unsupervised registration methods utilize a spatial transform layer with an image similarity loss to train the deep model. While the existing methods usually derive the dense deformation field over the entire image domain, the performance is prone to degrade for the situations such as large deformation.

Another category of registration is based on sparse landmarks [36, 37]. Landmarks are salient points of geometrical and anatomical uniqueness. Basically, the landmark-based registration methods are mainly applicable to estimate rigid or affine transformation. With the raising number of landmarks, they can also be used to describe more complex (e.g., non-rigid) deformation [38]. Note, the identification of landmarks is particularly important for these methods.

The landmarks can be manually or automatically detected. Most existing landmark-based registration methods rely on manually labeled landmarks. One may directly derive the sparse displacement from the manually labelled landmarks, and further obtain the dense motion field with the help of interpolation such as thin plate spline (TPS) [39]. However, the manual labeling of landmarks is inevitably time-consuming and prone to intra- and inter-observer uncertainty, which limits its real application. For automatic keypoint detection methods, they generally capture the keypoints through the image local features [40]. For instance, Yang et al. [41] presented a CNN feature-based motion tracking method for temporal images. Hosseini et al. [42] utilize feature points to estimation the cardiac motion problem. Li et al. [43] introduced a Multiview-based Parameter Free framework, where pixel-wise motion estimation has been used to derive feature points for identifying the number groups of people in a crowd video. In another study, Liu et al. [44] presented a point trajectories generation method by exploring the motion boundaries for video processing. However, the automatic landmark detection usually comes with complex point matching [45], such as iterative closest point [46], to track the movement of corresponding landmarks and then to complete the following image registration. Although it does not require manual annotation, the errors caused by automatic landmark detection and subsequent point matching limit the accuracy and robustness in motion estimation, especially for medical images with large motion.

Landmark detection is widely used in computer vision, e.g., for human face/pose recognition [47, 48, 49]. Initially, regression-based CNN methods utilize a fully convolutional layer to estimate the distinct positions. For instance, Yu et al. [50] presented a deformable landmark regression network that utilized a TPS layer to achieve the cross-talk between point-based deformation and images. Lv et al. [51] introduced a regression CNN with global stage and local stage to improve the performance on facial landmark detection.

Different from the regression approaches, localizing the landmark by predicting its heatmap is an alternative way. These methods usually use a kernel function to convert the coordinate information into the form of heatmap. Nie et al. [52] proposed a hierarchical contextual refinement network for landmarks to efficiently predict the human pose. Siarohin et al. [53] recently presented a heatmap-based landmark detection method to transfer the animation from video to still images. However, compared to exploring discriminate landmarks in 2D space, it is admittedly much more challenging for 3D landmark detection concerning the exponential growth in data complexity. Besides, compared with 2D facial images, 3D/4D medical images usually lack enough visual saliency, as the anatomical landmarks are severely interfered by noisy background.

The landmark detection serves as a core step in our DSD motion estimation framework. To perform the landmark detection from dense volumetric image, we propose an unsupervised learning-based landmark detection network. A standard multi-scale encoder-decoder architecture is used as the backbone, as shown in Fig. 3. To obtain representative landmarks [54], we firstly introduce a self-supervised pre-training strategy to help the encoder of the landmark detection network to focus on the potential motion-affected regions. Then, for obtaining accurate and consistent landmarks, multiple losses are introduced to help train the network to perform landmark localization. Note that none of the annotation information is used to supervise any training step.

In this section, the estimation of the dense motion field $\varphi_{t_{m}\leftarrow t_{n}}$ from the detected landmarks ($L_{t_{m}},L_{t_{n}}$) and their sparse correspondences are introduced. We firstly compute the sparse displacement $V_{t_{m}k\leftarrow t_{n}k},\{k=1,...,K\}$ from the detected landmarks $L_{t_{m}},L_{t_{n}}$. Then, we utilize a dense motion reconstruction network $\mathcal{M}$ (standard U-Net architecture) to estimate the dense motion field from the displacement of the detected landmarks.

The sparse displacement $V_{t_{m}k\leftarrow t_{n}k},\{k=1,...,K\}$ links two temporally-consecutive landmarks $L_{t_{n}k}$ and $L_{t_{m}k}$ in $I_{t_{n}}$ and $I_{t_{m}}$, respectively. As the displacement $V_{t_{m}k\leftarrow t_{n}k}$ is spatially sparse and cannot directly input to our motion network $\mathcal{M}$, we first filling the entire $\widetilde{\varphi}_{t_{m}k\leftarrow t_{n}k}$ grid with the value of $V_{t_{m}k\leftarrow t_{n}k}$ and acquire a dense volumetric deformation field $\widetilde{\varphi}_{t_{m}k\leftarrow t_{n}k}$. We then use $\widetilde{\varphi}_{t_{m}k\leftarrow t_{n}k}$ to translate the image $I_{t_{n}}$ via the warp layer $\mathcal{W}$, which results in $K$ different deformed images $\widetilde{I}_{t_{m}k\leftarrow t_{n}k}$ in total. Thus, the deformed images $\widetilde{I}_{t_{m}k\leftarrow t_{n}k}$ are input to $\mathcal{M}$.

However, the deformed $\widetilde{I}_{t_{m}k\leftarrow t_{n}k}$ only contains displacement information in images yet the landmark positions are missing. We thus further inject the corresponding landmarks’ heatmaps into $\mathcal{F}$. We specifically derive the variation between the corresponding heatmaps as

Finally, we concatenate the deformed images $\widetilde{I}_{t_{m}k\leftarrow t_{n}k}$ and the variation of heatmaps $\mathcal{F}_{t_{m}k\leftarrow t_{n}k}$ along the channel dimensionality, and feed them into the motion reconstruction network $\mathcal{M}$ to fuse the $\widetilde{\varphi}_{t_{m}k\leftarrow t_{n}k}$ into final dense motion field $\varphi_{t_{m}\leftarrow t_{n}}$ as

The training of $\mathcal{M}$ is supervised by a voxel-wise similarity loss calculated from the intensity image and its corresponding heatmaps as

In overall, by integrating the landmark detection network and the dense motion reconstruction network, our final training loss $\mathcal{L}$ is defined as

The weights $\lambda_{1}=1,\lambda_{2}=0.5,\lambda_{3}=20,\lambda_{4}=100,\lambda_{d}=500$ have been set empirically. As shown in Fig. 2, we thus achieve the dense-sparse-dense transferring and estimate the motion field between the input images. In the final, the reconstructed motion field can be employed as an initial transformation field and further refined via existing registration methods [60]. The refinement aims to boost the quality in estimating the motion field in detail, while the initial transformation produced by DSD has already addressed large deformation. Thus, the refinement only needs a small amount of iterative computation, which is highly efficient.

We demonstrate our method with three datasets: 4D Cardiac CT (4D-C-CT), 4D Lung CT (4D-L-CT) [61], and 4D MR cardiac ACDC [62].

We conduct comprehensive experiments to validate our proposed method - DSD (without subsequent refinement) and DSD + Demons (with refinement by Demons). First, we compare our proposed method with three different categories of the state-of-the-art motion estimation methods, including conventional non-rigid image registration method, unsupervised CNN-based registration, and feature point based registration. We implement all the compared methods following their default parameter settings in the same environment, including the same training set, and same testing set:

1. 1.

   Conventional non-rigid image registration methods - Demons [14] and Free-Form Deformation (FFD) [27]. Demons utilizes the intensity similarity of the images to derive the driving force, which deforms the moving image to the fixed image. FFD method requires to initialize a group of control points, and optimizes the control points’ displacement to interpolate the deformation field, which warps the moving image to the fixed image.

2. 2.

   Unsupervised deep CNN-based method – VoxelMorph [34]. The deep CNN network is trained via the loss from the similarity of image appearance and the smoothness of the deformation field in an unsupervised manner. In the inference stage, it can directly yield the entire deformation field given the moving and the fixed images.

3. 3.

   Feature point based method - FeaturePoints [42]. Image local feature points are firstly extracted via the popular scale-invariant feature transform (SIFT) descriptors, and then utilized to register the two images.

Next, for unsupervised landmark detection, we compare it with the state-of-the-art first-order motion model (FOMM) [53]³³3https://aliaksandrsiarohin.github.io/first-order-model-website/ in the ablation study. Note that, since FOMM is originally presented for 2D images, we have adapted FOMM to 3D and regarded 3D-FOMM as the baseline for comparison.

For all comparisons, we implement them using their default parameter settings. For Demons and FFD, we implement them via SimpleITK library. For evaluation, we use the popular metrics including Hausdorff Distance (HD), Dice Similarity Coefficient (DSC), Average Surface Distance (ASD), and Relative volume difference (RVD). Lower HD, ASD, RVD and higher DSC indicate better performance.

We have applied our DSD framework to both 4D-C-CT and 4D-L-CT datasets for the evaluation of its performance. Furthermore, we have employed the motion field produced by deep-learning-based methods as an initialization, and then refine it by Demons, following the strategy in [60]. Note that, the FFD is not selected to refine our motion field due to it does not support the initial deformation setting.

We also conduct an ablation study for our unsupervised landmark detection network, which is a critical component in our framework. To effectively discuss the impact of our proposed losses on the accuracy of motion estimation and landmark detection, we analyze our proposed method with three groups of different combinations of loss functions following the scheme in [64], i.e., Only-Keeping, Removing and Accumulating, as shown in Table IV shown. In the ‘Only-Keeping’ subgroup, we only keep one selected loss function to combine with the baseline model. For the second ‘Removing’ subgroup, we remove the selected loss function and employ all other proposed loss functions. In the last ‘Accumulating’ subgroup, we increasingly add the loss terms to the baseline model in the order of $\mathcal{L}_{1}$, $\mathcal{L}_{2}$ and $\mathcal{L}_{3}$. Note that, the baseline model (3D-FOMM) uses $\mathcal{L}_{4}$ and $\mathcal{L}_{d}$. We thus only investigate $\mathcal{L}_{1}$, $\mathcal{L}_{2}$ and $\mathcal{L}_{3}$ in this ablation experiment.

The motion estimation results in the ablation study are listed in Table IV for 4D-C-CT. Compared to the baseline model 3D-FOMM, our $\mathcal{L}_{1}$ loss greatly improves the effectiveness of landmark detection and also the accuracy in motion estimation. In contrast, using $\mathcal{L}_{2}$ and $\mathcal{L}_{3}$ alone does not significantly improve the model performance compared to the baseline, as shown in the ’Only-Keeping’ subgroup in Table IV. However, only relying on $\mathcal{L}_{1}$ and $\mathcal{L}_{2}$ losses is not able to handle the paired images with large motion. The simultaneous use of $\mathcal{L}_{1}$ and $\mathcal{L}_{3}$ can further significantly improve the accuracy of motion estimation. We attribute this to that our specific designed topological loss $\mathcal{L}_{3}$ effectively establish the conversation between the detected landmarks and the target organ’s anatomy. Therefore, the detected landmarks accurately reflect the dynamic changes of the target organ which can be further used to guide the motion estimation. Compared with $\mathcal{L}_{1}$ and $\mathcal{L}_{3}$, $\mathcal{L}_{2}$ does not greatly improve the accuracy of the model. We attribute this to that its effect mainly manifested in improving the description of landmark representation, especially without $\mathcal{L}_{3}$, as shown in Fig. 11. Meanwhile, $\mathcal{L}_{d}$ and subsequent added $\mathcal{L}_{3}$ will weaken its effect on the model performance because they also can guide model to pay attention to the motion-affected region.

Fig. 11 presents the landmark distributions with different loss settings. Adding the driven force, $\mathcal{L}_{1}$, greatly improves the space exploration of the landmark detection as compared with $BK$. Although the loss $\mathcal{L}_{1}$ can effectively spread the landmarks in 3D space, they are evenly distributed which cannot match the shape of the target organ (see Fig. 11(b)). On this basis, the temporal coherence motion loss, $\mathcal{L}_{2}$, further guides the model to explore the regions that may provide more effective guidance to the motion estimation. Finally, to make the detected landmarks conform to the anatomical topology, we add the reference topological loss, $\mathcal{L}_{3}$, which achieves the best quantitative and qualitative performance in all settings under comparison as shown in Table IV, Fig. 11(d), and Fig. 11(e). Our detected landmarks exhibit strong distribution consistency with the reference landmarks which are extracted from the manual segmentation annotations. This also reflects that the detected landmarks by our unsupervised landmark detection network can conform to the anatomical topology of the target organ.

In addition to verifying the proposed losses, we also evaluate the effects of different numbers of landmarks and report the results in Table V. We use 12, 16, 20 and 24 landmarks to represent the target organ, respectively. Note that, to effectively evaluate the impact of the numbers of landmarks on motion estimation, we leverage $BK+\mathcal{L}_{f}+\mathcal{L}_{b}$ and omit $\mathcal{L}_{o}$ to constrain the landmark detection and motion reconstruction. As expected, the more landmarks result in higher accuracy of the estimated motion fields (12-landmarks: DSC=0.722, HD=11.51; 16-landmarks: DSC=0.772, HD=9.26; 20-landmarks: DSC=0.814, HD=8.25; 24-landmarks: DSC=0.842, HD=7.19). The 24-landmark set achieves the best results in our experiments. We have used 24 landmarks in our study and achieved the state-of-the-art results. Nevertheless, we suggest that a higher number will further improve the performance, albeit incrementally.

We present a novel unsupervised topology extraction guided motion estimation method for 4D dynamic medical images. Our main findings are that our DSD method: 1) effectively captures the sparse topological anatomy of the volumetric image through unsupervised landmark detection with multiple losses, especially the landmark topological loss, which aligns the detected landmarks to conform to the anatomical prior; 2) is able to better reconstruct the dense deformation field, with motion tracking guided by the detected landmarks.

In our framework, landmark detection can be regarded as a process of sparse encoding of dense image representation, while dense motion reconstruction is a process of decoding from sparse representation. As shown in Fig. 11, we have found that even if the landmarks are not directly related to the shape of the target organ, the motion field between the paired input images can be relatively accurate (see Fig.11(b)). We attribute this to the network training where the conversation between the image coding (landmark detection) and the motion decoding (dense motion reconstruction) is established, even if the result of image coding (detected landmarks) is not morphologically meaningful. However, this compulsory learning limits the accuracy of the motion estimation. In contrast, our unsupervised landmark detection network can effectively capture the sparse landmarks to represent the obvious anatomy of target organ. This ensures that the motion field estimated from the detected landmarks is more reasonable.

In the comparison with the refinement scheme VoxelMorph + Demons and directly derive the motion field by Demons (iteratively derivation), our DSD + Demons achieves a better performance across all metrics on both of 4D-C-CT dataset and 4D-L-CT dataset (Table I and Table II). This is attributed to the fact that these methods, VoxelMorph and Demons, estimate the deformation only by calculating the similarity of local appearance while ignoring the global structural changes. This may lead to their inability to accurately estimate the moving direction of the target organ. Thus, this is prone to incorrect estimation and makes the subsequent optimization unable to further improve the estimation (see Fig. 8 and Fig. 9). In contrast, the initial motion field estimated from our DSD can effectively handle the large deformation and keep the correct anatomical shape under the guidance of the whole anatomical topology change. Then, our initial estimation can be further iteratively optimized to improve the detailed estimation accuracy. In addition, our ablation study further exemplifies that our unsupervised 3D landmark detection network has high consistency and anatomical significance in temporal landmark detection even when the target object shape changes greatly (see Fig. 11).

In conclusion, we introduce a novel learning-based method to estimate the motion for the dynamic medical images. We present an unsupervised 3D landmark detection method. Particularly, our proposed topological constraint can inject prior knowledge of human anatomies into landmark detection. This strategy can make the detected landmarks conform to the anatomical prior of the organ. In addition, our proposed landmark detection network is adaptable to manual landmarks, implying that the landmark detection network can be trained leveraging the manually annotated landmarks in supervised or weakly-supervised manner. We note that, based on our experiment, our proposed self-supervised pre-training strategy has a positive contribution to the model’s convergence, yet no significant impact upon the accuracy of motion estimation.

Several limitations in our work should be discussed. First, the location of each landmark is represented by an independent image-wise heatmap. This leads to a huge requirement in GPU memory usage when the number of landmarks is increased. Next, although our method can improve the accuracy of large motion estimation, it is not sensitive to detailed structure changes, thus requiring further refinement from conventional methods. In the future, we suggest improving the method by reducing the GPU memory. Therefore, better estimation results can be potentially obtained by increasing the number of landmarks.