Bidirectional Recurrence for Cardiac Motion Tracking with Gaussian Process Latent Coding

Authors: jiewen  yang
yiqun  lin
bin  pu
xiaomeng  li
  frame
Venue: NeurIPS_cc_2024_Conference
Published: 2024-04-30 13:45:34

Abstract: Quantitative analysis of cardiac motion is crucial for assessing cardiac function. This analysis typically uses imaging modalities such as MRI and Echocardiograms that capture detailed image sequences throughout the heartbeat cycle. Previous methods predominantly focused on the analysis of image pairs lacking consideration of the motion dynamics and spatial variability. Consequently, these methods often overlook the long-term relationships and regional motion characteristic of cardiac. To overcome these limitations, we introduce the GPTrack, a novel unsupervised framework crafted to fully explore the temporal and spatial dynamics of cardiac motion. The GPTrack enhances motion tracking by employing the sequential Gaussian Process in the latent space and encoding statistics by spatial information at each time stamp, which robustly promotes temporal consistency and spatial variability of cardiac dynamics. Also, we innovatively aggregate sequential information in a bidirectional recursive manner, mimicking the behavior of diffeomorphic registration to better capture consistent long-term relationships of motions across cardiac regions such as the ventricles and atria. Our GPTrack significantly improves the precision of motion tracking in both 3D and 4D medical images while maintaining computational efficiency. The code is available at: https://github.com/xmed-lab/GPTrack.

Main paper: 
 introduction: 
 Cardiac motion tracking from Cardiac Magnetic Resonance Imaging (MRI) and Echocardiograms is crucial in quantitative cardiac image processing. These imaging techniques provide comprehensive image sequences that cover an entire heartbeat cycle, allowing for detailed analysis of cardiac dynamics. Conventional non-parametric cardiac motion tracking approaches, such as B-splines [1] , Demons algorithms [2] and optical-flow based methods [3, 4, 5] , are commonly utilized due to their flexibility and ability to align detailed structures within images. However, these methods face significant challenges in motion tracking because they lack topology-preserving constraints and temporal coherence. The diffeomorphic registration method [6, 7, 8] , which formulates the registration process as a group of diffeomorphisms in Lagrangian dynamics, is a nice candidate for topology-preserving motion tracking. However, traditional optimization-based diffeomorphic registration methods are computationally intensive and sensitive to noises, hindering their applications in efficient cardiac motion tracking.Current deep learning-based techniques [9, 10, 11, 12, 13, 14, 15, 16] employ these advanced imaging modalities to register image pairs within the same patient, and some [14, 15, 16] adopt the diffeomorphic routine and learn the Lagrangian strain to describe the motion relationship between the reference frame and subsequent frames, aggregating dynamic information into consecutive cardiac motions as Lagrangian displacements. Although the above diffeomorphic methods better model the dynamic and continuous nature of cardiac motion, they still have room for improvement in handling the long-term temporal relationship in videos. For example, the approach [14] requires segmentation annotation to generate dense motion trajectories and calculate Lagrangian strain. Additionally, the methods outlined in [15, 16] are prone to error accumulation as Lagrangian displacements are integrated without regularizing temporal variations. In spatial views, while the global trajectory flow may follow a specific pattern, significant variations exist within each region regarding the phases, amplitudes, and intensities of the motion. For example, Figure 1 shows regions of the Right Atrium (red) and Myocardium (green) performing the opposite trajectories during the heartbeat cycle. Conversely, similar regions across different cases exhibit consistent motions. Hence, ignoring the regional scale may lead to a fragmented understanding of cardiac motion, underscoring the need for more nuanced analytical approaches. Furthermore, as shown in the right of Figure 1 , the deformation is bounded in the space of periodically specific human cardiac motion variation.To leverage discussed temporal and spatial information for cardiac motion tracking, we propose a novel framework named GPTrack. Our GPTrack has several appealing facets: 1) GPTrack employs the Gaussian Process (GP) to formulate the consistent temporal patterns in the latent space of diffeomorphic frameworks, promoting the consistency of cardiac motion; 2) GPTrack utilizes position information in the latent space to encode the statistics of the sequential Gaussian process, by which we model the region-specific motion and obtain a more precise estimation related to cardiac motion; 3) GPTrack leverages the inherent temporal continuity in cardiac motion by aggregating long-term relationships through forward and backward video flows, which mimics the forward-backward manner of classical diffeomorphic registration framework [6] . To evaluate the performance of GPTrack in cardiac motion tracking, we conduct experiments based on 3D Echocardiogram videos [17, 18] and 4D temporal MRI image [19] . Results in Tables 1, 2 and 3, show the GPTrack enhance the accuracy of motion tracking performance in a clear margin, without substantially increasing the computational cost in comparison to other state-of-the-art methods. Our contributions are summarized as follows:1. We propose a novel cardiac motion tracking framework named the GPTrack. This framework employs the Gaussian Process (GP) to promote temporal consistency and regional variability in compact latent space, establishing a robust regularizer to enhance cardiac motion tracking accuracy.The GPTrack framework is designed to capture the long-term relationship of cardiac motion via a bidirectional recursive manner, and its forward-backward manner mimics the workflows of the classical diffeomorphic registration framework. By this approach, our method provides a more accurate and reliable estimation of cardiac motion.Our GPTrack framework achieves state-of-the-art performance on both 3D Echocardiogram videos and 4D temporal MRI datasets, maintaining comparable computational efficiency. The results demonstrate that our method adapts effectively across different medical imaging modalities, proving its utility in different clinical settings. 
cardiac motion tracking via non-parametric registration approach: 
 Extensive works have been proposed to address registration by optimising within the space of displacement vector fields. Models related to elastic matching were proposed by [20, 21] . [22] utilized statistical parametric mapping for improvement. Techniques incorporating free-form deformations with B-splines and Maxwell demons were adapted by [1] and [2] , respectively. The Harmonic phase-based method, which utilizes spectral peaks in the Fourier domain for cardiac motion tracking, was introduced by [23] . This method calculates phase images from the inverse Fourier transforms and is specifically exploited in the analysis of tagged MRI. Popular formulations [6, 8, 24, 7] introduce topology-preserved diffeomorphic transforms. In the realm of diffeomorphic registration, inverse consistent diffeomorphic deformations have been estimated by [6] and [7] . Syn [24] proposed standard symmetric normalization. RDMM [8] considered regional parameterization based on Large Deformation Diffeomorphic Metric Mapping [6] . Despite their remarkable success in computational anatomy studies, these approaches are also time-consuming and susceptible to noise.Recent advancements in medical image registration have increasingly leveraged Deep Learning technologies. Pioneering studies [14, 25, 26, 27, 28] utilize ground truth displacement fields obtained by simulating deformations and deformed images, typically estimated using non-parametric methods. These approaches, however, may be limited by the types of deformations they can effectively model, which can affect both the quality and accuracy of the registration. Unsupervised methods, as discussed by [9, 10, 15, 29, 30, 31, 32] have shown promise by learning deformation through the warping of a fixed image to a moving image using spatial transformation functions [33] . These methods have been extended to include deformable models for single directional deformation field tracking [9, 26, 34] and diffeomorphic models for stationary velocity fields [32, 35, 36] . Further application of diffeomorphic models to cardiac motion tracking has been explored by [13, 15, 16, 30, 31] . These models predict motion fields that are both differentiable and invertible, ensuring one-to-one mappings and topology preservation. Recent studies in denoising diffusion probabilistic models (DDPM), such as [11] and [12] , have achieved considerable success in registration tasks. However, DDPM-based methods face challenges in building temporal connections and demand substantial computational resources. The DL-based optical flow (OF) methods [37, 38, 39, 40] apply widely in nature image motion tracking. However, as illustrated in [7, 15, 16] , due to annotations requirements and photometric constraints, they cannot be adopted in unsupervised cardiac motion tracking in medical image domains (See Section A1 in Appendix for detailed discussion). Diffeomorphic motion tracking and registration techniques are widely used in medical image analysis because they seek topology-preserving mapping between source and target images [6, 8] . Formally, given the source image x 0 and target image x 1 , the diffeomorphic registration aims at a family of differentiable and invertible mappings {ϕ t } t∈[0,1] with the boundary condition ϕ 0 = Id and ϕ 1 (x 0 ) = x 1 , where Id is the identity mapping. The diffeomorphism ϕ t can be parameterized as its derivatives (velocity field) v t as follows:EQUATIONwhere • is the composition operator. For numerical implementation, the associative property of the diffeomorphism group indicates ϕ t1+t2 = ϕ t1 • ϕ t2 , and the integral of Equation 1 can be approximated by ϕ t+δ ≈ ϕ t + v t δ for t ∈ [0, 1) and small enough δ. In this study, we follow the parameter settings of [7, 15, 16, 35, 36] and take δ = 1 2 N , N = 7 to discretize the path of diffeomorphic deformation. The overview pipeline of GPTrack (one layer). The x, h, ḣ and z denote input, forward hidden states, backward hidden states and latent coordinates. Feature ⃗ f t with probabilistic prior on the latent space via Gaussian Process then enters the decoder to predict the motion field ϕ. Subscript t denotes the t-th position in total T moments. elu(•) represent the exponential linear units [42] . 
motion tracking with gaussian process latent coding: 
 Our proposed method adopts the generative variational autoencoder (VAE) framework as the backbone of a diffeomorphic tracking network, as suggested by previous methods [15, 16, 36] . However, as shown in Figure 2 , the first difference between our method and other methods is that ours allows the registration network to aggregate the spatial information temporally, both forward and backward (see Figure 2(a) ). The conventional approaches [9, 10, 15, 16, 24, 31] only conduct between two adjacent frames, which ignore the relationship of long-term dependency of the cardiac motion (see Figure 2 (b)). Secondly, despite the large deformation through the path of diffeomorphism, the space of periodically specific human cardiac motion variation is bounded. Though methods [14, 15, 16] use Lagrangian strain to formulate the continuous dynamic of cardiac motion, however, without considering the motion consistency between two adjacent state spaces, the Lagrangian strain is prone to error accumulation and degrade the tracking performance. To address this problem, we take the simple yet efficient Bayesian approach, which employs the Gaussian Process (GP) to model cardiac motion dynamics in the compact feature space, predicting more consistent motion fields over dynamics parameters. Our proposed GPTrack can also be easily extended to other modalities or motion-tracking tasks, such as 3D Echocardiogram videos and 4D cardiac MRI.In this paper, we follow the research [41] that employs the recursive manner in transformer for sequential data. As shown in Figure 3 left, the GPTrack pipeline comprised the GPTrack layer for feature extraction, the Gaussian Process (GP) layer for modelling the cardiac motion dynamics, and the Decoder for motion field estimation. Given the sequential 4D inputs {x t } T t=1 , x t ∈ R H×W ×D×1 , where H, W, D, T denote the height, width, depth, length of the input. For each x t , we first decompose it to P non-overlapping patches of shape p × p × p, where P = H p × W p × D p and p, H p , W p , D p ∈ Z + . We then embed each patch as a feature with C channels via embedding layers and disentangle patches with the dimension of R P ×C from x t . The GPTrack layer then takes the x, both forward and backward hidden states ⃗ h, ⃗ h ∈ R P ×C as the input, then predicts the motion field ϕ via the decoder after the GP layer. Note that the initial hidden states of forward ⃗ h 0 and backward ⃗ h 0 are set as zero.The GPTrack layer consists of two independent GPTrack cells that respond to forward and backward computation. Similar to the [41, 43] , adapting the hidden state to maintain and aggregate the sequential information allows the input with variable length. Meanwhile, the conventional convolutional neural network or vision in the transformer-based method is limited by the fixed input length. Furthermore, medical images such as Echocardiogram videos usually consist of hundreds of frames that cover multiple heartbeat cycles. Hence, parallel computing all frames requires a large amount of computational consumption, which hinders the application in real scenarios limited by lowcomputational devices. To this end, as shown in Tables 1 and 2 , our GPTrack is able to formulate the variable temporal information while maintaining the comparable computational cost.Using the forward GPTrack cell shown in the right of Figuer 3 as an example, the x t and ⃗ h t-1 followed by the addition of learnable position encoding pos t ∈ R P ×C are respectively normalized by Layer Normalization. The linear self-attention then computes the attentive weight A t ∈ R P ×C of combined x t and ⃗ h t-1 . The above operations can be formulated as follows:EQUATIONwhereW Q , W K , W V ∈ R 2C×Care learnable weights of different projection layers named query, key and value, δ(•) represent the exponential linear units elu(•) [42] , LN(•) and ⊕ are the layer normalization and the concatenation operation, respectively.In order to raise the descendant hidden state ⃗ h t that aggregates information before t + 1 moment. The attentive weight A t then conducts the element-wise addition with ancestral positional encoded ⃗ h t-1 . Additionally, the A t takes the positional encoded x t as the residual connection and applies the addition operation. The Feed Forward Network denoted as FFN(•), is then introduced to output the feature of t-th moment. The formulation can be written as follows:EQUATIONIn the Bidirectional Forward-Backward Recursive Cell, both forward and backward share the same computation processes through the GPTrack cell. The only difference between the two directions is that the forward process starts from the first moment x 0 of input while the backward starts from the last moment x T . Hence, the feature f t in t-th moment is formulated asf t = ⃗ f t + ⃗ f t , f t ∈ R P ×C, where ⃗ f t aggregate the forward information from 0 to t, while the ⃗ f t aggregate the backward information from the last frame to the t-th frame.The primary objective of integrating the Gaussian Process (GP) is to establish a probabilistic prior in the latent space that incorporates prior knowledge. Specifically, it posits that cardiac motion across different individuals within the same region should yield similar motion fields in latent space encodings. Furthermore, as illustrated in Figure 1 , cardiac structures within an individual that are spatially distant or exhibit motions in opposite directions should consistently adhere to the periodic pattern of motion.Initially, we define a covariance (kernel) function for the GP layer as depicted in Figure 3 . We design the prior for the latent space processes to be stationary, mean square continuous, and differentiable in sequential motion fields. This design stems from our expectation that the latent functions should model cardiac motion more prominently than visual features. Consequently, we anticipate the latent space to manifest continuous and relatively smooth behaviour. To this end, we employ the isotropic and stationary Matern kernel (refer to Equation 4) to fulfil the required covariance function structure:EQUATIONwhere ν, σ, l > 0 are the smoothness, magnitude and length scale parameters, K ν is the modified Bessel function, and D(•, •) denotes the distance metric between features of two consecutive motion fields. Our goal is to formulate cardiac motion as robust prior knowledge applicable to unseen data.To address this, we propose a position-related distance measurement.As outlined in section 3.3 and referenced in [44] , we utilize a learnable parameter ṗ os ∈ R P ×C as the spatial positional encoding for each region, which provides relative or absolute positional information about the decomposed patches. To capture the periodic temporal positional information of cardiac motion, distinct from the spatial encoding ṗ os, we apply sine and cosine functions of various frequencies as temporal encoding { p os t } T t=1 , p os t ∈ R P ×C for each moment. The overall positional encoding at moment t is formulated as:EQUATION)The i denotes the i-th position of C channels, and n is the scaling factor. We assign independent GP priors to all values in {z t } T t=1 , z t ∈ R P ×C to disseminate temporal information between frames in the sequence. During this process, we regard the sequential output {f t } T t=1 of GPTrack as noisecorrupted versions of the ideal latent space encodings, formulating the inference as the following GP regression model with noise observations:z t ∼ GP µ(pos t ), κ pos t-1 , pos t , f t = z t + ϵ t , ϵ t ∼ N (0, σ 2 ),where σ 2 is the noise variance of the likelihood model set as the learnable parameter in GPTrack. The above Gaussian process can be treated as the temporal sequence with intrinsic Markov property, and we adopt the methodology of connecting the Gaussian process with state space model [45] to decrease its computational complexity from O(T 3 ) to O(T ), where T is the number of time step. Concretely, the Gaussian process of Equation 6 corresponds to the following linear stochastic differential equation:EQUATIONwith the solution as:EQUATIONwhere z(t) := (z(t), d dt z(t)), w(t) is the zero-mean Gaussian random process, and h := (0, 1) T is used for modelling observation model. The state transition matrix (vector) A and b can be calculated from the covariance function κ of the Gaussian process. For the Matern kernel shown in Equation 4, we take ν = 3 2 , and corresponding state transition matrix A and vector b [46] read as:EQUATIONThen we can discretize Equation 7and get its weakly equivalent state-space model of Equation as:EQUATIONwhere Φ t = exp D(post,pos t-1 )A , n t ∼ N (0, Φ t bQ w (t -1, t)b T Φ T t ), and Q w (t, t -1) is the covariance of w(t). Given the initial value z 0 ∼ N (µ 0 , Σ 0 ) with µ 0 = 0 and Σ 0 = diag( σ 2 2 , 3σ 2 l 2 ), we can sequentially calculate the posterior distribution z t |f 1:t-1 ∼ N (µ t , Σ t ) and z t |f 1:t ∼ N (µ t , Σ t ) using update criterion of Kalman filter for state space model [47] as:EQUATIONwhere k t := Σth hΣth+σ 2 is the optimal Kalman gain at time t. The output of the GP layer in t-th moment thus can be formulated as z GP t = ReLU(k t z t ), where ReLU is the activation function. In the final, the t-th motion field θ t is obtained by decoder from the z GP t .The decoder takes the Gaussian-process coding {z GP t } T t=1 as velocity filed to composite the diffeomorphic motion field ϕ according to the criterion of Section 3.1. Here, we adopt the training loss of [15, 16] , which minimizes the integration of four components summarized as follows: a) Dissimilarities of tracking results between adjacent states from both forward and backward; b) Smoothness of motion fields between adjacent states from both forward and backward; c) Dissimilarities of tracking results between the start state and each state; d) Smoothness of motion fields between the start state and each state. The overall loss function L is formulated as:T -1 t=1 [L kl (x t , x t+1 ) a) +α 1 (L sm (ϕ t:t+1 ) + L sm (ϕ t+1:t ) b) ) + α 2 L nc (x t+1 , x 1 • ϕ 0:t+1 ) c) +α 3 L sm (ϕ 1:t+1 ) d) ],where α 1 , α 2 and α 3 are loss weights, and ϕ t1:t2 is the motion field from state t 1 to t 2 . L kl (x t , x t+1 ) = KL(q(z GP t |x t ; x t+1 )||p(z GP t |x t ; x t+1 ))+KL(q(z GP t |x t+1 ; x t )||p(z GP t |x t+1 ; x t )) is the summation of forward and backward VAE losses with latent coding z GT t , posterior distribution q and conditional distribution p, L nc is the negative normalized local cross-correlation metric, andL sm (ϕ) = ||∇ϕ|| 22 is the ℓ 2 -total variation metric. 
datasets: 
 CardiacUDA [17] . The CardiacUDA dataset collected from two medical centers consists of 314 echocardiogram videos from patients. The video is scanned from the apical four-chamber heart (A4C) view. In this paper, we conduct training and validation in the A4C view that consists of 314 videos with 5 frame annotations in the Left/Right Ventricle and Atrium (LV, LA, RV, RA). For testing, we report our results in 10 videos with full annotation provided by the CardiacUDA.CAMUS [18] . The CAMUS dataset provides pixel-level annotations for the left ventricle, myocardium, and left atrium in the Apical two-chamber view, which consists of 500 echocardiogram videos in total. There are 450 subjects in the training set with 2 frames annotated in the Left Ventricle (LV), Left Atrium (LA) and Myocardium (Myo) in the end-diastole (ED) and end-systole (ES) of the heartbeat cycle. The remaining 50 subjects without any annotation masks belong to the testing set.ACDC [19] . The ACDC dataset consists of 100 4D temporal cardiac MRI cases. All data provide the segmentation annotations corresponding with the Left Ventricle (LV), Left Atrium (LA) and Myocardium (Myo) in the end-diastole (ED) and end-systole (ES) during the heartbeat cycle.Training. We trained the model using the Adam optimizer with betas equal to 0.9 and 0.99. The training batch size of the model was set to 1. We trained for a total of 1000 epochs with an initial learning rate of 5e -4 and decay by a factor of 0.5 in every 50 epochs. During training, for CardiacUDA [17] and CAMUS [18] , we resized each frame to 384 × 384 and then randomly cropped them to 256 × 256. All frames were normalized to [0,1] during training. In temporal augmentation of datasets [17, 18] , we randomly selected 32 frames from an echocardiogram video with a sampling ratio of either 1 or 2. For ACDC [19] , we resampled all scans with a voxel spacing of 1.5 × 1.5 × 3.15mm and cropped them to 128 × 128 × 32, normalized the intensity of all images to [-1, 1]. For spatial data augmentation of all datasets, we randomly applied flipping, rotation and Gaussian blurring.In CardiacUDA, we split the dataset into 8 : 2 for training and validation. During testing, we reported results in 10 fully annotated videos. In the CAMUS [18] dataset, videos without annotation are used for only training, while we randomly split the remaining 450 annotated videos into 300/50/100 for training, validation and testing. In the ACDC [19] , following the [11, 12] , we split the training set in the ratio of 90 and 10 for training and testing. The reproduced methods strictly follow the official code and the description in the paper. For all experiments, We use Intel(R) Xeon(R) Platinum 8375C with 1× RTX3090 for both training and inference. All reproduced methods strictly followed the training settings with their original paper in the same experimental environment.Inference. For CardiacUDA and CAMUS, we resized videos to 384 × 384, cropped to 256 × 256 in central and normalized to [0,1]. We sample 32 frames that cover the segmentation annotation. When the sequence has more than 32 frames, the extra frames will be removed from the sequence, except for the first and the last one. The ACDC dataset remains the same sampling strategy as training in the inference stage, without any argumentation except for normalizing intensity to [-1,1].Evaluation Metrics. For the evaluation of the quality of registered target frames, we follow [12] to use the Peak Signal-to-Noise Ratio (PSNR) and Structural Similarity Index (SSIM) [48] to measure whether the Lagrangian motion field is accurately estimated between the first frame and the following wrapped frames. We also use the Dice [49] score to measure the discrepancy between tracked and ground-truth cardiac segmentation. For CardiacUDA [17] , only the first frame and corresponding segmentation are provided for tracking the following 32 frames, then report the averaged results by the above metrics in these 32 frames. For CAMUS [18] and ACDC [19] , frame and segmentation of the ED stage are used to track the go-after frames, and we report all the metrics in the ES stage. We evaluate diffeomorphic property by computing the percentage of non-positive values of the Jacobian determinant det(J ϕ ) ≤ 0 (%) on the Lagrangian motion field. In order to access the evaluation of comparing the physiological plausibility following the [50, 51] , we also compute the mean absolute difference between the 1 and Jacobian determinant (||J ϕ | -1|) over the tracking areas. For a fair comparison, we evaluate the computational efficiency and report the computational time in seconds (Times), the parameter quantities in millions (Params), and the tera-floating point operations per second (TFlops). We also provide the result evaluated by Hussdorf Distance (HD) in Tables B1, B2and B3 of Appendix Section B. Table 3 : The segmentation performance1 of different cardiac structures in the CAMUS [18] . Result of 3D Echocardiogram Video. In Table 1 , we compare our method with state-of-the-arts 2D registration [9, 10, 12, 15, 31, 52] and non deep learning [6, 8, 24] methods in CardiacUDA dataset. In comparison to DeepTag [15] , our GPTrack-XL reach 82.37% and 12.75 in DICE and HD scores with the best non-positive Jacobian determinant value, which denotes the learned motion field is smooth. The registration quality of our method also achieved the best with 32.03 (3.50 ↑) and 80.04 (3.64 ↑) in PSNR and SSIM compared to the second-best method, respectively. Though GP-Track introduces more learnable parameters and requires a small amount of additional computation, slightly increases inference time and TFlops compared to DeepTag [15] and VM-DIF [9] . GPTrack surpasses all other methods in the registration and verifies the necessity of formulating a strong temporal relationship among frames by using the recursive manner. Tables 1 and 3 show registration results of cardiac structures (LV, RV, LA, RA, Myo). Compared to other methods, our GPTrack can outperform existing baseline methods by a substantial margin. In areas such as the left atrium (LA) and left ventricular (LV), which usually cause larger deformation, our method can also provide better alignment than other approaches. As shown in Figure 4 , the tracking error of our GPTrack and other methods show a significant difference in LA (Labelled by Orange Colour) and LV (Labelled by Green Colour) when compared to the ground truth. The methods [24, 10, 9, 31] present inaccuracy tracking results due to lack of the constraints on the consecutive motion and ignore the long-term temporal information. These results further verify that our specifically designed GPTrack for modelling motion patterns is more suitable for cardiac motion tracking on echocardiography.Result of 4D Temporal Cardiac MRI Dataset. We compare our GPTrack method against with the state-of-the-art deep learning based methods [9, 10, 11, 12, 15] and different Non-rigid approaches [6, 8, 24] . As illustrated in Table 2 , our GPTrack-XL achieves the best average DICE score of 82.65 compared to FSDiffReg [11] with 82.30. In registration quality, our GPTrack-XL reaches the highest scores, 31.52 and 86.19, in PSNR and SSIM, respectively. Moreover, our Jacobian determinant on deformation fields shows numbers comparable to other methods with the diffeomorphic constraint.All results are based on our fast and lightweight model, reducing around 96.93% inference time, 17.2% model parameters and 76.02% computational consumption (TFlops) compared to the secondbest performance. In comparison to the diffusion-based method [12, 11] , which requires enormous computation that hinders real-time inference and is nearly impossible to deploy in real scenarios, the GPTrack preserve light-weight and considerable performance by formulating cardiac motion patterns as the Gaussian process latent coding and bidirectionally understand the cardiac motion. The visualization result in Figure 5 also indicates our GPTrack can achieve better tracking accuracy.The Scale of Model Hyper-parameters. Table 4 shows the settings of the 2D/3D GPTrack-M/L/XL. For the 3D echocardiogram video dataset and 4D cardiac MRI dataset, the GPTrack with different scales has different patch sizes and dimension numbers. Referring to the result performed by Tables 3, 1 and 2, the registration result can be boosted by increasing the layers number and dimension size of GPTracks according to different requirements. The Ablation of Bidirectional Recursive Manner and Gaussian Process. Table 5 illustrates the Bi-directional layer in GPTracks can better formulate both forward and backward motion fields by aggregating the temporal information of the cardiac heartbeat cycle. The Gaussian Process models cardiac motion with strong prior knowledge from data, which makes more accurate predictions of the deformation field. The Figure 6 illustrate the tracking error from 1-st to 32-th frame in CardiacUDA [17] full annotated set. Our GPTrack and the DeepTag [15] both use the Lagrangian strain. However, the accuracy degrades significantly without aggregating the temporal information and GP when the input length increases. As shown in Figure 6 , the tracking error after the 20-th frame becomes larger by using only Lagrangian strain, while introducing GP and bidirectional recursive methods can efficiently eliminate the tracking error by predicting more accurate motion fields. 
conclusion and limitation: 
 In this paper, we proposed a new framework named GPTrack to improve cardiac motion tracking accuracy. GPTrack innovatively aggregates both forward and backward temporal information by using the bidirectional recurrent transformer. Furthermore, we introduce the Gaussian Process to model the variability and predictability of cardiac motion. In experiments, our framework demonstrates the state-of-the-art in 3D echocardiogram and 4D cardiac MRI datasets. The limitation of our framework is that we use positional encoding as the prior knowledge of the cardiac motion, which may degrade the tracking performance in out-of-domain datasets. In our future work, we will build a more robust representation of cardiac motion and further our work across different medical domains. We also provide the illustration of Broader Impacts, please see Section A2 in Appendix. 
 
 Appendix: We provide the experiment reported on Hussdorf Distance (HD) for datasets CardiacUDA [17] , CAMUS [18] and ACDC [19] in Tables B1, B2 and B3, respectively. Furthermore, the consequent tracking results of 3D echocardiogram and 4D Cardiac MRI are presented in Figures B1 and B2 . • For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out of range (e.g. negative error rates). • The answer NA means that there is no societal impact of the work performed.• If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact. • Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations. • The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster. • The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology. • We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution. • For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.Segmentation results are reported in the Dice score (%). Bold, underline denote the best results and the second best performance. The superscript * indicates computational time reported in only CPU implementation. We use the t-test for statistical significance analysis, where the p-value between the two methods is p < 0.05, indicating statistically significant improvements 


Decision: This work proposes an unsupervised framework for cardiac motion tracking by considering the temporal and spatial dynamics of cardiac motion. Reviewers recognise that this work has merits in the innovation of model design, the state-of-the-art results, robust and rigorous methodology, and the clear articulation of findings. At the same time, the reviewers raise issues related to more evaluation and comparison, the impact of the number of frames, deeper analysis and stronger connection, and the statistical significant analysis. The authors have made efforts to provide responses to the raised issues. Two of the three reviewers indicate that the rebuttal has addressed their concerns and give positive final ratings. The third reviewer does not respond to the rebuttal and the rating is at borderline reject. The AC agrees that this work has its novelty and values and that the authors have well responded to the reviews. Meanwhile, it is recommended that this work apply statistical significant test (say, paired samples t test) to the results in Tables 2 and 3 to better assess the improvement. Considering all the factors, AC recommends this paper for acceptance. The authors shall effectively incorporate the key information in the rebuttal to further strengthen this work.

Reviews: review: paper_summary: paper_summary: The authors proposed a latent modeling framework for cardiac motion tracking. They introduced the GPTrack module for image encoding, which considers both forward and backward information flow. A Gaussian process was integrated to describe the motion prior. Extensive experiments were performed on both echocardiography and cardiac MRI datasets, demonstrating good results.
strength_weakness: strength_weakness: 1.It appears that the authors misunderstood the cardiac MRI in Figure 2. The short-axis CINE image shows the myocardium of the left ventricle (green) and right ventricle (red) rather than the left atrium.
2.The authors used PSNR, SSIM, and DICE to evaluate motion tracking performance, which is reasonable. However, they did not evaluate the physics/physiological plausibility. For the myocardium, which is considered an incompressible material, the determinant of the Jacobian |J| should be around 1. I suggest the authors add some evaluation comparing the physiological plausibility [1].
[1] doi: 10.1016/j.media.2022.102682
questions: questions: 1.Are there any results showing the impact of varying the number of frames?
2.The myocardium is the tissue of interest for motion tracking, as the blood pool of the left/right ventricle involves flow filling. Could there be more focus on comparing the myocardium?
limitations: limitations: The authors described the limitations appropriately.
paper_summary: paper_summary: 1. This paper presents GPTrack, an unsupervised framework designed to thoroughly investigate the temporal and spatial dynamics of cardiac movement. 
2. GPTrack refines motion tracking by utilizing sequential Gaussian Processes within the latent space and encoding statistical data with spatial information at each time point, thereby robustly fostering temporal consistency and accommodating spatial variations in cardiac dynamics. 
3. This paper aggregates sequential data in a bidirectional, recursive fashion, emulating the principles of diffeomorphic registration to more effectively capture the enduring motion relationships across cardiac areas, including the ventricles and atria. 
4. The method enhances the accuracy of motion tracking, all while preserving computational efficiency.
strength_weakness: strength_weakness: 1. Line 26, Can optical flow not capture temporal coherence? I doubt that's the case.
2. Line 31, The citation is not formatted correctly.
3. Line 61, "contribuion" repeats the idea of "facets" mentioned earlier, which is redundant.
4. There is no experimental comparison with methods such as denoising diffusion probabilistic models (DDPM).
5. Line 111, How is the relationship between the diffeomorphism \phi_t and the method proposed in the paper demonstrated?
6. The experiment does not include a comparison with the method in reference [36].
7. The paper does not explain, demonstrate, or theoretically analyze long-term dependencies, nor does it clarify the duration of "long-term."
8. Many acronyms are repeated unnecessarily.
9. Line 135, What is meant by "motion consistency," and how do you define and assess its consistency?
10. Line 148, Here, x_{t} represents features; it should not use the same mathematical symbol as the raw data to avoid confusion.
11. Figure 3, in the right subplot, what is meant by "elu(x)"?
12. Line 206, the formula below this line is not numbered. But the author refers to this equation.
13. How is the multidimensional z_{t} output generated? Multidimensional z_{t} involves multitask Gaussian processes, which are not mentioned in this paper.
14. Line 221, sometimes \mu_{t} is a scalar, and sometimes it is a vector.
15. Line 258 and Line 272, there is a conflict. Why is this done? Please explain the motivation.
16. Figure 6, Why does the error decrease and then increase? Moreover, this trend is similar even across different methods.
17. There should be a comparison with optical flow-based methods.
18. Why is FSDiffReg only included in Table B3?
19. The discussion of results is somewhat superficial. A deeper analysis linking the findings back to the research questions and theoretical framework would enrich the paper's contribution to the field.
20. There are inconsistencies in the terminology used throughout the paper. Consistently using the same terms would improve clarity and reduce confusion for the reader.
21. The paper would benefit from a stronger connection between the theoretical framework and the empirical results.
questions: questions: Included in the Strengths and Weaknesses.
limitations: limitations: Included in the Strengths and Weaknesses.
paper_summary: paper_summary: To capture the long-term relationship in cardiac motion, the authors proposed GPTrack, a novel unsupervised framework crafted to fully explore the temporal and spatial dynamics of cardiac motion. They proposed employing the sequential Gaussian Process in the latent space and aggregating sequential information in a bidirectional recursive manner. Through experiments on 3D Echocardiogram videos and 4D temporal MRI datasets, they demonstrate that they achieve state-of-the-art performance.
strength_weakness: strength_weakness: 1. Some details are lacking or not clear in the manuscript, regarding the data processing, how the proposed method is applied to the complete cardiac cycle, and the detailed difference in training and architecture between the proposed method and traditional methods.
2. Lack of comparison to the latest registration methods.
3. No statistical significance analysis is provided.
4. More detailed discussions about the potential application and limitations are needed.
questions: questions: 1.	How is the proposed method applied to the complete cardiac cycle?  What's the difference between the proposed method and the traditional DL-based registration network? Is the comparison to previous pairwise registration approaches a fair comparison?
2.	Does the author use a pair of images or a complete cycle of images as the input? If using the complete cardiac cycle, the incorporation of additional temporal information would significantly increase the computational cost and registration parameters, what is the potential application of the proposed methodology?
3.	The cardiac cycle of different subjects contains different frames, how did the authors process the dataset and train the proposed method?
4.	Why separate the original images into several patches? It would lead to significant unalignment if the motion is large.
5.     Are the improvements significant?
limitations: limitations: The authors only mentioned one sentence of limitation in the conclusion, which needs further discussion.

Comments: title: We are here always ready to answer and address your questions or concerns.
The deadline for the author-reviewer discussion period is now approaching, and we are looking forward to hearing from you.
Your feedback is really important to us
May we kindly know if our detailed explanations and more experimental validations have addressed your concerns?
We have provided a comprehensive response for your question and concerns and we look forward to hearing from you.
Thank you for your review and here are some additional information of our rebuttals.
Thank you for your review and here are some additional information of our rebuttals.
Thank you for your review and here are some additional information of our rebuttals.
comment: **Dear Reviewer Tej6:**

The deadline for the author-reviewer discussion period is almost close. Thank you for your efforts and we are all waiting for your responses. We are here always ready to answer and address your questions or concerns.

Thank you for your time and consideration, and we are looking forward to hearing from you.

Best regards
**Dear Reviewer Tej6:**

The deadline for the author-reviewer discussion period is now approaching. We understand your time is very valuable and we appreciate all your efforts in writing reviews. During the remaining discussion period, we are still waiting for your response and hope our comprehensive response can address your concerns. We will also provide detailed answers to your other questions or concerns.

We all thank you for your time and consideration, and we are looking forward to hearing from you.

Best regards
**Dear Reviewer Tej6:**

We thank you very much that you can review our paper and provide valuable feedback that helps improve our work. As the author-reviewer discussion period is coming to a close, we would like to know whether our responses have addressed your concerns. We know that the discussion may take up your time. However, your feedback is really important to us, and we would like to address your question in our discussion further.

Thank you for your time and consideration, we look forward to hearing from you.

Best regards
We greatly appreciate your efforts and help us improve our work! In our future work, we will further devote ourselves to the cardiac motion tracking task as it is important in cardiac function assessment. Furthermore,  we will follow the previous related research to evaluate the physics/physiological plausibility and make our research and results more reliable in real scenarios. We believe that our work can contribute to follow-up research in medical image analysis.
Thank you for your rebuttal. My concerns have been addressed. It is a nice paper to read. I have increased my rating to accept.
**Dear Reviewer ssxr:**

Thank you again for your valuable comments and your efforts, your kind suggestions help further our paper's quality. We know that there may still be some concerns that need to be addressed and clarified. In our rebuttal, we have provided a comprehensive response to your concerns through detailed explanations and more experimental validations. As the author-reviewer discussion period is about to close, may we kindly know if the responses have addressed your concerns?

Thank you for your time and consideration.

Best regards
**Dear Reviewer Tej6:**

Thank you again for your valuable comments and your efforts. As the author-reviewer discussion period is coming to a close, we know that there exist some concerns that need to be addressed and clarified. In our rebuttal, we have provided a comprehensive response to your concerns through detailed explanations and thorough experimental validations. May we know if the responses have addressed your concerns?

Thank you for your time and consideration.

Best regards
Thank you very much again for your time and efforts, we are delighted that our rebuttal has addressed your concerns! Your suggestions and questions help us better improve our work. We believe that our research can further contribute to cardiac motion tracking. We expect that the employment of the recursive manner as well as the Gaussian Process can bring more insight to the follow-up research related to medical image analysis.
The rebuttal addressed my concerns
Thank you very much for your review and valuable feedback. Your suggestions and comments are important to us and we will address all your concerns in the rebuttal. Our new experiments are included in the ***uploaded PDF file*** to better illustrate your questions. All the ***references and citations*** are included in the content of the paper's rebuttal at the beginning. 

Thank you for your efforts and valuable time!
We appreciate your review and valuable feedback, we will address your concerns in the rebuttal. Our new experiments are included in the ***uploaded PDF file*** to better illustrate your questions. All the ***references and citations*** are included in the content of the paper's rebuttal at the beginning. 

Thank you for your efforts and valuable time!
Thank you for your review and valuable feedback, we will address your concerns in the following. Our new experiments are included in the ***uploaded PDF file*** to better illustrate your questions. All the ***references and citations*** are included in the content of the paper's rebuttal at the beginning.

We thank you very much for all your efforts and valuable time! 

**Q3. A pair of images or a complete cycle of images as the input?**

The input of our model is a sequence of images (frames). However, those frames are not inputted all at once, instead, they are inputted one by one like RNNs and our model is shared for different input frames. The number of frames depends on the task requirements. For example, if we want to estimate the motion of the entire cardiac cycle, the input sequence includes the video frames of the complete cardiac cycle. But if we just want to estimate the motion between two frames, it is ok to only input them in order. In our experiments, motion tracking requires the estimation of the entire cardiac cycle that follows previous works [15, 16] settings. Hence, we unsupervised train network with the entire cardiac cycle as [15, 16], and inference the result between every two adjacent frames in each time step.