This paper demonstrates a robust temporal image registration framework. Large motion presents significant challenges for time-course analysis in dynamic image series. To correct for motion in medical images, registration that produces voxel correspondences is commonly employed. Unfortunately, registration often fails if the motion is substantial. We take advantage of the temporal smoothness of motion to achieve robust alignment of volumetric MRI time series. This work is motivated by a clinical research application in monitoring placental and fetal health, which we introduce below.
I-a Motivating application
In-utero blood oxygenation level dependent (BOLD) magnetic resonance imaging (MRI) is a promising imaging tool for studying functional dynamics of placenta and fetal organs [1, 2, 3, 4]. Changes in fetal and placental oxygenation levels with maternal hyperoxygenation can be used for charactering and detecting placental dysfunction . Temporal MRI data suffers from serious motion artifacts due to maternal respiration, unpredictable fetal movement, and signal non-uniformities , as illustrated in Fig. 1. Different types of motion are presented in these MRI time series. The fetal brains move as rigid objects with a wide range of motion; the placenta moves more locally and deforms non-rigidly. In our study, the image series include about 500 volumes (frames) that are analyzed to study fetal and placental physiology.
I-B Related Work
Temporal registration has been previously investigated in time series of cardiac images [7, 8, 9, 10] and lung images [10, 11, 12, 13]. Both cardiac motion and breathing patterns are somewhat regular and smooth across time, enabling biomechanical modeling [7, 11, 8]. Electrocardiography (ECG) gating or respiratory gating can be incorporated in the cardiac or lung image motion correction, by either enforcing the transformation models to be periodic or penalizing non-periodic transformations in the cost function. These approaches cannot be generalized to images that contain complex non-periodic motion patterns.
Existing registration techniques for dynamic medical imaging data can be categorized into two distinct groups. The first type of methods performs registration based on the differences between each pair of consecutive frames. Each pair is aligned, and the estimated transformations are concatenated to estimate alignment of all images in the series [13, 14]. In application to long image series like our data, this approach leads to substantial accumulated errors after several concatenation steps. These methods only take the moving image and the reference image into account in each pairwise registration step, and ignore essential temporal information from other images in the series. The second type applies groupwise registration to simultaneously align each individual image to a group template – estimated or selected from one of the input images – aiming to yield acceptable registration results for all images in the series [15, 16, 17, 18, 12, 19, 10, 20, 21]
. Taking this approach in our application, registration fails for a great number of images that are substantially different from the template, necessitating outlier detection and rejection. Alternatively, a global objective function based on pairwise differences between each individual image in the series and an implicitly defined template frame can incorporate the temporal information . This leads to a global optimization process, i.e., the algorithm performs pairwise registration of consecutive frames iteratively until the entire series comes into alignment. This type of global optimization is extremely challenging in large image sets like our time series.
The problem of temporal alignment has also been investigated in longitudinal image analysis that aims to characterize temporal variations [23, 24, 20, 25, 26]. The goal is to examine the behavior of a biological system driven by development, aging, or disease processes . The statistical analysis is performed via temporal shape regression on a Riemannian manifold [23, 24, 20, 25] to construct a trajectory that characterizes the shape variation across time. Longitudinal data is often sparsely distributed in time, and temporal smoothness constraint is therefore imposed. Typically, the challenges in longitudinal analysis involve growth or degeneration, which requires modeling tools that capture subtle deformation and are distinct from those in our problem of MRI time series alignment in the presence of substantial motion.
Unlike other applications that can require subjects to remain stationary, fetal and placental imaging is an application with irregular motion that cannot be controlled at the time of acquisition.
I-C Proposed Temporal Registration Method
To overcome the difficulty of correcting large and irregular motion in a long MRI time series, we assume a Markov structure that induces temporal smoothness and offers an efficient way to estimate the transformations by serially solving pairwise registration problems. We derive filtered estimates of the transformations from the hidden Markov model that represents the MRI series. The resulting sequential procedure estimates the transformations of the template to each frame in the series. We provide a flexible framework for temporal dense alignment of image time series.
Unlike the methods that perform pairwise registration of consecutive frames and concatenate the estimated transformations, we use temporal smoothness explicitly at every estimation step. The sequential estimation framework yields an efficient algorithm that makes alignment of long MRI series feasible. Our approach is related to filtering methods in respiratory motion modeling . We do not model the motion explicitly but rather capture it through transformation of organs of interest to template to adapt to complex motion patterns.
We demonstrate and evaluate the proposed method in the context of the in-utero BOLD MRI study by providing registration-based tracking for organs of interest. We employ the rigid motion model for fetal brains and a non-rigid B-Spline motion model for placentae. The experimental results offer significant improvements over pairwise registration in alignment of fetal brains that tend to move substantially, and robust improvements in alignment of placentae.
This paper presents an extended version of the preliminary algorithm presented at the International Conference on Medical Image Computing and Computer-Assisted Intervention (MICCAI) . In this paper, we derive in detail the instantiation of a hidden Markov model in the context of MRI time series. We model the motion as a latent state in the hidden Markov model, derive a temporal registration algorithm using filtered estimates, and discuss the challenges of developing a smoothing algorithm. Here we employ rigid and B-Spline transformation models on the regions of fetal brains and placentae respectively, instead of one single transformation model for the entire volumetric image used in , to adapt to the different motion patterns in the two regions. Finally, extensive experiments investigate the sensitivity of the resulting registration on the model parameters.
This paper is organized as follows. In the next section, we present the generative model for temporal motion and observed images using hidden Markov model, derive the resulting algorithm for filtered estimates, and provide implementation details. In Section III, we introduce the in-utero BOLD MRI data and report the experimental results. In Section IV, we discuss future directions, methodological challenges of developing a backward message passing algorithm for temporal registration, and potential extension for other transformation models. Conclusions follow in Section V.
We employ hidden Markov model (HMM) to represent the temporal changes in the MRI time series. One image in the MRI time series is chosen as a template, whose transformations represent the rest of the images in the series. We treat the transformation as a hidden state, and each image as an observation that depends on the template and the corresponding hidden state. In the following subsections, we present inference in HMMs using filtered estimates [28, 29] in the context of temporal registration, and describe the implementation details.
Ii-a HMM and Filtered Estimates
The temporal dynamics of latent (hidden) states in a HMM exhibit the Markov property, i.e., the current state depends on the history only through the previous state. Thus, inference based on the entire series of past observations is equivalent to integrating the information from the previous state and from the corresponding observation, leading to an efficient sequential estimation algorithm.
In our application, one image is chosen as the template to be treated as a global parameter shared by all observed images. We assume that the template deforms at each time point, and the transformation of the template is the latent state at time , where is the total number of images in the series. The observed image at time is generated by applying the anatomical transformation to the template independently of all other time points. Fig. 2 illustrates this model.
Our aim is to estimate the latent state variable from the observations , where we use and to denote sub-series of the MRI time series and the associated transformations . Formally, the estimation of the latent variable is reduced to maximization of its posterior distribution parameterized by the template . This posterior distribution, also known as the filtered estimate of the state, can be efficiently computed using forward message passing [29, 28], also known as sequential estimation. The posterior distribution is determined by integrating the previous posterior distribution with the temporal dynamics and the likelihood of the current observation:
where Eq. (2) follows from the Markov property:
Applied recursively, the forward message passing produces the posterior distribution for each time point in the number of steps that is linear with . Similarly, backward message passing enables inference of the posterior distribution , which is based on all data, often referred to as smoothing. In this paper, we incorporate temporal smoothness in temporal registration by filtering. We discuss the challenges of developing a smoothing algorithm for temporal registration in Section V.
Ii-B Estimating Temporal Transformations in HMM
Since it is intractable to integrate over all possible transformation fields, we employ a commonly used approach of approximating the posterior distribution as a point estimate that maximizes the probability distribution. In particular, if is the best transformation estimated by the algorithm for time point , the message from the state to state is the optimal transformation :
yielding the (recursive) estimate of the hidden state at time point :
The estimate is then used to compute the optimal transformation at time , and so on until we reach the end of the time series.
Ii-C Model Setup
We model the likelihood term in Eq. (5
), also known as emission probability, as the exponential likelihood on image similarity:
where is a measure of dissimilarity (distance) between images. This probability measures how different the transformed template and the currently observed image are given the template and the transformation .
The state transition probability in Eq. (5) encourages temporal smoothness of transformations and spatial smoothness of :
where is an appropriate norm that encourages transformation to be similar to , and is the regularization term that encourages spatial smoothness and other desired properties. The regularization term is often used in image registration to restrict the estimated motion to a specific transformation group. and are regularization parameters that control the amount of temporal and spatial smoothness.
We manipulate Eq. (5) to obtain
and observe that this optimization problem is reduced to serial pairwise image registration of the template and the observed image .
The algorithm proceeds as follows. Given the transformation estimate of the template to the image , we apply the registration algorithm to and by minimizing the cost function in Eq. (8), while also using as an initialization, resulting in the estimate .
This HMM-based temporal transformation model can be easily generalized to different types of similarity metrics between moving images and fixed images, temporal regularization, and transformation models. The model can be readily augmented to account for signal non-uniformities and a latent template that is estimated jointly with the transformations, similar to prior work in groupwise registration [10, 20, 21, 12]. In the following subsection, we describe our implementation for aligning in-utero BOLD MRI time series.
Ii-D Implementation Details
We implemented our method using the Insight Segmentation and Registration Toolkit (ITK) . Negative cross-correlation  is employed as the metric of image dissimilarity as cross-correlation can handle images with signal non-uniformities that are present in the in-utero BOLD MRI data. Cross-correlation is computed over a 555 voxel volumetric window. We implemented the rigid transformation model to account for the fetal brain motion and the B-Spline transformation model to characterize the placental transformation. We apply rigid transformation model for fetal brain alignment by defining a region-of-interest (ROI) mask in the template that is used to compute the registration cost function. For placentae, the B-Spline transformation model corrects the non-rigid deformation in addition to motion. For the B-Spline model, the grid spacing was set to be 101010 voxels, in the volume of 150150120 voxels of our MRI data.
Our choice of the transformation model offers a natural implementation of spatial smoothness. For the rigid model, the third term in Eq. (8) disappears, and for B-Spline model, it becomes a sum-of-squares of the B-Spline coefficients. The temporal regularization is implemented as the norm of the difference between the current transformation and the previously estimated transformation parameters. Specifically, the six parameters in 3-D rigid transformation model are split into two groups: three parameters representing rotation angles and three parameters representing translation. The weight of the second term in Eq. (8) is then represented by and respectively for the two groups of parameters. and are set to and in the rigid registration. is set to in the B-Spline registration.
Note that the template can be any image in the series, and the algorithm can proceed in either ascending or descending order of the series. In our experiments, we chose the image in the middle of each series as the template . In each series, we register images to the template in a descending order, and register images to the template in an ascending order.
Iii Experimental Results
We evaluated our method on in-utero BOLD MRI time series acquired as part of a placental function study .
Ten pregnant women were consented and scanned on a 3T Skyra Siemens scanner (single-shot GRE-EPI, in-plane resolution, slice thickness, interleaved slice acquisition, TR, TE, FA) using 18-channel body and 12-channel spine receive arrays. Each series contains around 300 volumetric images. This study included three singleton pregnancies, six twin pregnancies, and one triplet pregnancy, between 28 and 37 weeks of gestational age.
To eliminate the effects of slice interleaving, we resampled odd and even slices of each volumetric image onto a common isotropicmm
image grid. We split each interleaved volumetric image into even and odd slices, and linearly interpolated the acquired slices to complete the volume. The number of images in each MRI series was thus doubled to around 600. To enable quantitative evaluation, we manually delineated the fetal brains (total of 18) and placentae (total of 10) in six frames (, , , , , ) in each series.
To evaluate the temporal model, we compare it to a variant of our method that does not assume the temporal structure and instead aligns each image in the series to the reference frame using the same registration algorithm used by our method. This approach is commonly used for fetal MRI series alignment [22, 33]. Algorithmically, this corresponds to setting in Eq. (8) and initializing the registration step with an identity transformation instead of the previously estimated transformation . The other parameter settings in the temporal model are the same with those in the pairwise model.
To quantify the accuracy of the alignment, we transform the manual segmentation labels (fetal brains and placentae) in the template to the five segmented frames () in each series using the estimated transformations. We employ Dice coefficient  to quantify volume overlap between the transferred and the manual segmentations. In our application, the goal is to study the average temporal signal for each ROI, and therefore delineation of an ROI provides an appropriate evaluation target. In addition, we also investigate how the alignment accuracy changes with the weight of the temporal smoothness term in Eq. (8). For B-Spline model, we additionally evaluate the log determinants of Jacobian of the estimated transformation fields.
Fig. 3 provides example results from the study for fetal brains and placentae. We observe that the reference frame is warped accurately by the temporal registration algorithm in the regions of fetal brains and placentae to represent the first frame in the series that is substantially different from the template. The delineations achieved by transferring manual segmentation labels from the reference frame to the coordinate system of the selected frame are in good alignment with the manual segmentation outlines for that selected frame.
Fig. 4 reports volume overlap statistics (Dice coefficient) for the fetal brains in this study. We observe that temporal alignment improves volume overlap in important ROIs over the baseline pairwise registration. We note that temporal alignment offers particularly substantial gains in cases with a lot of motion, i.e., low original volume overlap, in fetal brains. Fig. 4 also shows how alignment accuracy changes with the weights of temporal smoothness regularization in Eq. (8). For the rigid model, the Dice coefficients are more sensitive to , and the registration achieves the best accuracy for . Note that the temporal smoothness of translation is still somewhat maintained by initializing the current registration with the previous transformation estimation, even if . In our tests, the temporal rigid registration yields the best result when and for fetal brains.
Fig. 5 reports volume overlap statistics for the placentae in this study. Our temporal alignment offers consistent improvement for placentae for all cases over pairwise registration to the reference frame. Fig. 5 also reports the sensitivity of the alignment accuracy with respect to the weight of temporal smoothness regularization in Eq. (8). For B-Spline registration on placenta, the algorithm yields the best results when . Finally, Fig. 5
shows the mean and the standard deviation of the log determinants in each case (in total, 10 cases). The mean values of the log determinants in each subject are all betweenand . A log determinant value above indicates volume expansion of the transformation, and a log determinant value below indicates volume compression. All the determinants of Jacabian of our estimated B-Spline transformations for placentae are positive.
Our temporal registration approach improves alignment of fetal brains and placentae in the in-utero MRI times series. The temporal model offers substantial improvement for fetal brain alignment in the presence of large motion, which presents significant challenges for the pairwise registration method. The incorporation of temporal smoothness improves the alignment of the placenta for every subject. In the study, the log determinants of Jacobian of the estimated B-Spline transformations are close to , which indicates no over-expansion or over-compression in the resulting placenta transformation estimates.
One important downstream use of this temporal registration method is modeling of ROI-specific intensity changes in the in-utero MRI time series [32, 22]. Currently the signal dynamics in fetal brains, which tend to have a wide range of motion, include substantial fluctuations, due to spatial signal non-uniformities that affect the estimated signal if the fetus moves substantially in the volumetric image. The cross-correlation measure of image similarity used in our temporal registration captures local image differences, but tracking the intensity changes in corresponding regions across images is still affected by the signal non-uniformities. Future work will obtain robust estimates of the in-utero MRI signal time-courses by augmenting the temporal method with a model of ROI-specific intensity changes.
In addition, the registration of each frame to the template frame can be based on all the images in the time series by including a backward message pass. One of the challenges is to find a template that is close to both the first observation (image) and the last observation (image) in the HMM. Otherwise, estimating the first or the last hidden state (transformation) might fail due to the substantial difference between the moving image and the template image, since there is no temporal information passed to the very first step of the sequential estimation .
Our temporal registration method naturally accepts different types of transformation models besides the rigid model and the B-Spline model. For example, it can be extended to diffeomorphic transformation. The diffeomorphic demons algorithm 
adapts convolutional kernels on the incremental vector field and the updated displacement field at each iteration acts as a regularization in the cost function, which can naturally be implemented with our temporal registration method to incorporate temporal smoothness.
We present a HMM-based registration method to align images in MRI time series. Forward message passing incorporates the temporal model of motion into the estimation procedure. The filtered estimates are therefore based on not only the present image frame and the template, but also on the previous frames in the image series. Our method can be easily extended on other forms of transformation models and image similarity metrics. The experiments on in-utero BOLD MRI time series demonstrate the promise of our approach in this novel and challenging application.
-  V. Schöpf, G. Kasprian, P. Brugger, D. Prayer, Watching the fetal brain at “rest”?, International Journal of Developmental Neuroscience 30 (1) (2012) 11–17.
-  A. Sørensen, D. Peters, C. Simonsen, M. Pedersen, B. Stausbøl-Grøn, O. B. Christiansen, G. Lingman, N. Uldbjerg, Changes in human fetal oxygenation during maternal hyperoxia as estimated by BOLD MRI, Prenatal diagnosis 33 (2) (2013) 141–145.
-  A. Sørensen, D. Peters, E. Fründ, G. Lingman, O. Christiansen, N. Uldbjerg, Changes in human placental oxygenation during maternal hyperoxia estimated by blood oxygen level-dependent magnetic resonance imaging (bold mri), Ultrasound in Obstetrics & Gynecology 42 (3) (2013) 310–314.
-  J. Luo, E. A. Turk, T. Hahn, M. Teulon Gonzalez, B. Gagoski, C. Bibbo, A. Palanisamy, C. Tempany, A. Torrado-Carvajal, N. Malpica, J. Mart nez Gonz lez, J. N. Robinson, J. A. Hernandez-Tamames, E. Adalsteinsson, P. E. Grant, Human placental and fetal response to maternal hyperoxygenation in IUGR pregnancy as measured by BOLD MRI, in: In Proceedings of the 23rd Annual Meeting of ISMRM, Toronto, Ontario, Canada, 2015, International Society of Magnetic Resonance in Medicine (ISMRM), 2015, p. 633.
-  S. Aimot-Macron, L. Salomon, B. Deloison, R. Thiam, C. Cuenod, O. Clement, N. Siauve, In vivo MRI assessment of placental and fetal oxygenation changes in a rat model of growth restriction using blood oxygen level-dependent (bold) magnetic resonance imaging, European radiology 23 (5) (2013) 1335–1342.
-  C. Studholme, Mapping fetal brain development in utero using MRI: the big bang of brain mapping, Annual review of biomedical engineering 13 (2011) 345.
-  H. Sundar, C. Davatzikos, G. Biros, Biomechanically-constrained 4D estimation of myocardial motion, in: Medical Image Computing and Computer-Assisted Intervention–MICCAI 2009, Springer, 2009, pp. 257–265.
-  J. Park, D. Metaxas, A. A. Young, L. Axel, Deformable models with parameter functions for cardiac motion analysis from tagged MRI data, Medical Imaging, IEEE Transactions on 15 (3) (1996) 278–289.
-  R. Chandrashekara, A. Rao, G. I. Sanchez-Ortiz, R. H. Mohiaddin, D. Rueckert, Construction of a statistical model for cardiac motion analysis using nonrigid image registration, in: Information processing in medical imaging, Springer, 2003, pp. 599–610.
-  C. Metz, S. Klein, M. Schaap, T. van Walsum, W. J. Niessen, Nonrigid registration of dynamic medical imaging data using nD+t B-splines and a groupwise optimization approach, Medical image analysis 15 (2) (2011) 238–249.
-  J. R. McClelland, D. J. Hawkes, T. Schaeffter, A. P. King, Respiratory motion models: a review, Medical image analysis 17 (1) (2013) 19–42.
-  E. Rietzel, G. T. Chen, Deformable registration of 4D computed tomography data, Medical physics 33 (11) (2006) 4423–4430.
-  J. M. Reinhardt, K. Ding, K. Cao, G. E. Christensen, E. A. Hoffman, S. V. Bodas, Registration-based estimates of local lung tissue expansion compared to xenon CT measures of specific ventilation, Medical image analysis 12 (6) (2008) 752–763.
-  V. Boldea, G. C. Sharp, S. B. Jiang, D. Sarrut, 4d-ct lung motion estimation with deformable registration: quantification of motion nonlinearity and hysteresis, Medical physics 35 (3) (2008) 1008–1018.
-  S. K. Balci, P. Golland, W. Wells, Non-rigid groupwise registration using b-spline deformation model, Open source and open data for MICCAI (2007) 105–121.
-  K. K. Bhatia, J. V. Hajnal, B. K. Puri, A. D. Edwards, D. Rueckert, Consistent groupwise non-rigid registration for atlas construction, in: Biomedical Imaging: Nano to Macro, 2004. IEEE International Symposium on, IEEE, 2004, pp. 908–911.
E. G. Miller, N. E. Matsakis, P. A. Viola, Learning from one example through shared densities on transforms, in: Computer Vision and Pattern Recognition, 2000. Proceedings. IEEE Conference on, Vol. 1, IEEE, 2000, pp. 464–471.
-  L. Zöllei, E. Learned-Miller, E. Grimson, W. Wells, Efficient population registration of 3d data, in: International Workshop on Computer Vision for Biomedical Image Applications, Springer, 2005, pp. 291–301.
-  S. Marsland, C. J. Twining, C. J. Taylor, A minimum description length objective function for groupwise non-rigid image registration, Image and Vision Computing 26 (3) (2008) 333–346.
-  S. Durrleman, X. Pennec, A. Trouvé, J. Braga, G. Gerig, N. Ayache, Toward a comprehensive framework for the spatiotemporal statistical analysis of longitudinal shape data, International journal of computer vision 103 (1) (2013) 22–59.
-  N. Singh, J. Hinkle, S. Joshi, P. T. Fletcher, Hierarchical geodesic models in diffeomorphisms, International Journal of Computer Vision (2015) 1–23.
-  E. A. Turk, J. Luo, B. Gagoski, J. Pascau, C. Bibbo, J. N. Robinson, P. E. Grant, E. Adalsteinsson, P. Golland, N. Malpica, Spatiotemporal alignment of in utero bold-mri series, Journal of Magnetic Resonance Imaging.
-  B. C. Davis, P. T. Fletcher, E. Bullitt, S. Joshi, Population shape regression from random design data, International journal of computer vision 90 (2) (2010) 255–266.
-  A. R. Khan, M. F. Beg, Representation of time-varying shapes in the large deformation diffeomorphic framework, in: Biomedical Imaging: From Nano to Macro, 2008. ISBI 2008. 5th IEEE International Symposium on, IEEE, 2008, pp. 1521–1524.
-  J. Fishbaugh, S. Durrleman, G. Gerig, Estimation of smooth growth trajectories with controlled acceleration from time series shape data, in: International Conference on Medical Image Computing and Computer-Assisted Intervention, Springer, 2011, pp. 401–408.
-  M. Reuter, B. Fischl, Avoiding asymmetry-induced bias in longitudinal image processing, Neuroimage 57 (1) (2011) 19–21.
-  R. Liao, E. A. Turk, M. Zhang, J. Luo, P. E. Grant, E. Adalsteinsson, P. Golland, Temporal registration in in-utero volumetric mri time series, in: International Conference on Medical Image Computing and Computer-Assisted Intervention, Springer, 2016, pp. 54–62.
L. E. Baum, T. Petrie, Statistical inference for probabilistic functions of finite state markov chains, The annals of mathematical statistics 37 (6) (1966) 1554–1563.
C. M. Bishop, Pattern recognition, Machine Learning.
-  T. S. Yoo, M. J. Ackerman, W. E. Lorensen, W. Schroeder, V. Chalana, S. Aylward, D. Metaxas, R. Whitaker, Engineering and algorithm design for an image processing api: a technical report on itk-the insight toolkit, Studies in health technology and informatics (2002) 586–592.
-  B. B. Avants, C. L. Epstein, M. Grossman, J. C. Gee, Symmetric diffeomorphic image registration with cross-correlation: evaluating automated labeling of elderly and neurodegenerative brain, Medical image analysis 12 (1) (2008) 26–41.
-  J. Luo, E. A. Turk, C. Bibbo, B. Gagoski, D. J. Roberts, M. Vangel, C. M. Tempany-Afdhal, C. Barnewolt, J. Estroff, A. Palanisamy, et al., In vivo quantification of placental insufficiency by bold mri: A human study, Scientific reports 7 (1) (2017) 3713.
-  W. You, I. E. Evangelou, Z. Zun, N. Andescavage, C. Limperopoulos, Robust preprocessing for stimulus-based functional mri of the moving fetus, Journal of Medical Imaging 3 (2) (2016) 026001–026001.
-  L. R. Dice, Measures of the amount of ecologic association between species, Ecology 26 (3) (1945) 297–302.
-  T. Vercauteren, X. Pennec, A. Perchant, N. Ayache, Diffeomorphic demons: Efficient non-parametric image registration, NeuroImage 45 (1) (2009) S61–S72.