1 Introduction
Slicetovolume registration has received increasing attention during the last decades within the medical imaging community. Given a tomographic 2D slice and a 3D volume as input, this challenging problem consists in finding the slice (extracted from the input volume and specified by an arbitrary rigid transformation) that best matches the 2D input image. We stress the fact that we are working with 2D slices (e.g. ultrasound (US)) as opposed to projective 2D images (e.g. Xray images). This is important since both problems are usually refereed as 2D/3D registration, even if they are intrinsically different. In slicetovolume registration, every pixel from the 2D image corresponds to a single voxel in 3D space. However, in a projective 2D image every pixel represents a projection of voxels from a given viewpoint.
One can formulate different versions of slicetovolume registration, depending on several aspects of the problem such as the matching criterion used to determine the similarity between the images, the transformation model we aim at estimating, the optimization strategy used to infer the optimal transformation model (continuous or discrete) and the number of slices given as input. In this work, we propose an iconic method (where the matching criterion is defined as a function of the image intensity values) to infer rigid transformation models (specified using 6DOF). The input consists of a single slice and a single volume, and we formulate it as a discrete optimization problem.
Discrete methods, where the tasks are usually formulated as a discrete labeling problem on a graph, have become a popular strategy to model vision problems [1] (and particularly, biomedical vision problems [2]) thanks to their modularity, efficiency, robustness and theoretical simplicity. This paper presents a graphbased formulation (inspired by the work of [3, 4]) to solve rigid (only) slicetovolume registration using discrete methods. As we will discuss in section 1.2, other works have tackled similar problems. However, to date, no work has shown the potential of discrete methods to deal with rigid slicetovolume registration. Our main contribution is to put a new spin on graphbased registration theory, by demonstrating that discrete methods and graphical models are suitable to estimate rigid transformations mapping slicetovolume. We validate our approach using a dataset of magnetic resonance images (MRI) of the heart, and we compare its performance with a stateoftheart approach based on continuous optimization using simplex method. Moreover, in the spirit of encouraging reproducible research, we make the source code of the application publicly available at the following website: https://gitlab.com/franco.stramana1/slicetovolume.
1.1 Motivation
In the extensive literature of medical image registration, it is possible to identify two main problems which motivated the development of slicetovolume registration methods during the last decades. The first one is the fusion of preoperative highdefinition volumetric images with intraoperative tomographic slices to perform diagnostic and minimally invasive surgeries. In this case, slicetovolume registration is one of the enabling technologies for computeraided image guidance, bringing highresolution preoperative data into the operating room to provide more realistic information about the patient’s anatomy [5]. This technique has been used when dealing with several scenarios such as liver surgery [6], radiofrequency thermal ablation of prostate cancer [7], minimally invasive cardiac procedures [8], among many others.
The second problem is the correction of slicewise motion in volumetric acquisitions. In a variety of situations, inter slice motion may appear when capturing a volumetric image. For example, in case of fetal brain imaging (essential to understand neurodevelopmental disabilities in childhood and infancy) [9], fetus motion generates inconsistencies due to the slice acquisition time. Another case is related to functional magnetic resonance images (fMRI), usually acquired as time series of multislice singleshot echoplanar images (EPI). Patient head motion during the experiments may introduce artifacts on activation signal analysis. Slicetovolume registration can be used to alleviate this problem by registering every slice to an anatomical volumetric reference following the wellknow mapslicetovolume (MSV) method [10].
1.2 Previous Work
Graphbased image registration, where the task is conceived as a discrete labeling problem on a graph, constitutes one of the most efficient and accurate stateoftheart methods for image registration [11]. Even if they have shown to be particularly suitable to estimate deformable nonlinear transformations [12, 13], they were also adapted to the linear case [3]. Most of the publications on the field focus on registering images which are in dimensional correspondence (2D/2D or 3D/3D). In case of projective 2D/3D image registration, only linear transformations were estimated using discrete methods by [4, 3]. More recently, several graphbased approaches to perform deformable slicetovolume registration were introduced in [14, 15, 16]. In these works, rigid transformations were computed as a byproduct of the deformable parameters, leading to unnecessary computational burden (since rigid transformations are by far lower dimensional than deformable ones). To the best of our knowledge, rigid (only) slicetovolume registration has not been formulated within this powerful framework. To date, all the methods focusing on this challenging problem are based on continuous (e.g. simplex [7], gradient descent [9], Powell’s method [17]
, etc.) or heuristic approaches (evolutionary algorithms
[18], simulated annealing [19]), missing the aforementioned advantages offered by discrete optimization. Based on the work of [4], we propose a discrete Markov Random Field (MRF) formulation of this problem, delivering more precise results than the stateoftheart continuous approaches. Moreover, inspired by the work of [20]in the context of vector flow estimation, we discuss how continuous stateoftheart approaches can be used to further refine the rigid transformations obtained through discrete optimization, resulting in more accurate solutions.
2 Rigid SlicetoVolume Registration through Markov Random Fields
We formulate rigid slicetovolume registration as an optimization problem. Given a 2D image , and a 3D image , we aim at recovering the rigid transformation specified by that better aligns both images, by solving:
(1) 
where corresponds to the slice extracted from image (using trilinear interpolation) and specified by the rigid transformation (as explained in Figure 1). is the socalled matching criteria, that quantifies the dissimilarity between the 2D image and the slice . Alternative matching criteria can be adopted depending on the type of images we are registering. For example, in monomodal cases where intensities tend to be linearly correlated in both images, simple functions such as sum of absolute differences (SAD) or sum of squared differences (SSD) may make the job. However, for more complicated cases like multimodal registration (where the relation between intensity values in both images is usually nonlinear), more elaborated functions like mutual information (MI) are applied.
This optimization problem is commonly solved through continuous (gradient or nongradient based) methods, which are considerably sensible to initialization and may be stuck in local minima. As discussed in Section 1.2, in this work we model rigid slicetovolume registration as a discrete labeling problem following the discretization strategy proposed by [3].
2.1 Rigid SlicetoVolume Registration as a Discrete Labeling Problem
Rigid slicetovolume registration, as well as many other problems in computer vision, can be formulated as a discrete labeling problem on a pairwise Markov Random Field (MRF)
[1]. Formally, a discrete pairwise MRF is an undirected graph , where each node represents a discrete variable. Any two variables depend on each other if there is an edge linking the corresponding nodes. The range of values that can be assigned to a discrete variable is determined by the label space . A discrete labeling problem on a pairwise MRF consists on assigning a label to every , such that the following energy is minimized:(2) 
where is a labeling assigning a label to every , are the unary potentials associated to and are the pairwise potentials associated to edges . These functions return scalar values when labels are assigned to variables . Since we pose the optimization as a minimization problem, potentials must associate lower values to labelings representing good solutions, and higher values otherwise.
In the formulation presented in Eq 1, one would like to explore the space of parameters and chose the values giving the best matching. Since we are modeling it as a discrete problem, we must adopt a discretization strategy for the (naturally) continuous space of rigid transformations . In [3], authors proved that it is possible to estimate linear (an particularly, rigid body) transformations by solving a discrete and approximated version of this formulation. Following their proposal, we model rigid slicetovolume registration through a graph , associating every parameter of the rigid transformation to a variable , giving a total of 6 variables (nodes in the graph). is a fully connected pairwise graph where , meaning that all variables (parameters) depend on each other. Note that this pairwise model is clearly an approximation, since the real dependency between the parameters is not pairwise but highorder. However, as stated in [3], similar approximations have shown to be good enough to estimate linear transformations, while making the problem tractable.
In our discrete strategy, every parameter is updated through a discrete variation associated to the label . Given an initial transformation , we explore the space of solutions by sampling discrete variations of , and choosing the one that generates the slice best matching image . Therefore, for a maximum size and a quantization factor , we consider the following variations to the initial estimate of : . The total number of labels results . Note that 0 is always included since we give the possibility of keeping the current parameter estimate. For example, in case that corresponds to , and , the label space of will correspond to .
Ideally, we would like to explore the complete search space around given by all possible combinations of labels. Since we have an exponential number of potential solutions, we adopt a pairwise approximation where only variations for pairs of variables are considered. This variations are encoded in the pairwise terms of the energy defined in Eq 2 as . Here denotes the updated version of , where only and were modified according to the labels and , while the rest of the parameters remained fixed. Unary potentials are not considered since we are only interested in the interaction between variables. Therefore, the discrete version of the optimization problem introduced in Eq 1 becomes:
(3) 
where the optimal labeling represents the final rigid transformation used to extract the solution slice .
2.2 Discrete Optimization
We solve the discrete multilabeling problem from Eq 3
using FastPD. FastPD is a discrete optimization algorithm based on principles from linear programming and primal dual strategies, which at the same time generalizes the well known
expansion [21]. One of the main advantages of FastPD is its modularity/scalability, since it deals with a much wider class of problems than expansion, being an order of magnitude faster while providing the same optimality guarantees when performing metric labeling [22]. Our problem does not fulfill the conditions to be considered a metric labeling problem (we refer the reader to [23] for a complete discussion about metric labeling); however, FastPD has shown promising results for similar formulations [3].FastPD solves a series of maxflow mincut [24] problems on a graph. In that sense, it is similar to expansion which also performs discrete inference on multilabel problems by solving successive binary maxflow mincut problems. The main difference between these approaches is the construction of the graph where maxflow mincut algorithm is applied. expansion constructs the binary problem by restricting the label space, so that the only options for a given variable are to remain in its current assignment, or to take a label (which varies in every iteration). Instead, FastPD constructs these binary problems by performing a Linear Programming Relaxation (LPR) of the integer program that represents the discrete MRF formulation.
2.3 Incremental Approach
Discrete approximations of continuous spaces usually suffer from low accuracy (since it is bounded by the quality of the discretization). Thus, we adopt an incremental approach to explore the space of solutions in a finer way. The idea is to successively solve the problem from Eq 3, using the solution from time as initialization for time , keeping a fixed number of labels but decreasing the maximum sizes in a factor . Moreover, we also adopt a pyramidal approach, where we generate a Gaussian pyramid for both input images and , and we run the complete incremental approach on every level of the pyramid. In that way, we increase the capture range of our method.
2.4 Simplex Refinement Step
Let us advance one of the conclusions of this work, so that we can motivate the last step of our approach. In Section 3.1, we compare the performance of our discrete approach with a continuous method based on simplex [25] algorithm. As we will see, when the initialization is good enough, both continuous and discrete approaches perform well. In fact, in some cases, simplex is delivering more accurate solutions than discrete. However, as we move away from the initialization, discrete optimization gives more and more significant improvements, thanks to its wider capture range. In order to improve the accuracy of our proposal, and inspired by similar conclusions discussed by [20] in the context of vector flow estimation, we refine the results obtained with our approach by optimizing Eq 1 through simplex, using the discrete solution as initialization.
3 Experiments and Results
In this section, we present the results obtained using the proposed method (considering also the refined version), and we compare them with a stateoftheart approach based on continuous optimization trough downhill simplex [25] (also known as NelderMead or amoeba method). Simplex is one of the most popular optimization algorithms used to deal with rigid slicetovolume registration (some examples are [26, 10, 19, 27]). It is a continuous and derivativefree method, which relies on the notion of simplex (which is a special polytope of vertices living in a dimensional space) to explore the space of solutions in a systematic way. We used a dataset composed of MRI images of a beating heart. Given an initial sequence of 3D images of a beating heart (with a resolution of voxels and a voxel size of ), we generated slices which were used for two different experiments.
3.1 Implementation Details and Parameters Setting
We implemented the three versions of the algorithms discussed in this paper (simplex, discrete and refined) mainly using Python and ITK for image manipulation ^{1}^{1}1The source code can be downloaded from
https://gitlab.com/franco.stramana1/slicetovolume. For simplex optimization we used the method implemented in scipy.optimize package, while discrete optimization was performed using a Python wrapped version of the standard C++ implementation of FastPD. In all the experiments, we used a pyramidal approach with 4 Gaussian levels (3D images where not downsampled in z axis because of the low resolution of the images in this direction). The matching criterion adopted in all the experiments was the sum of squared differences, since pixel intensities are equivalent in both 2D and 3D images. The matching criterion based on SSD is simply defined as:
(4) 
For the discrete case, at every pyramid level we decreased the maximum label size for both, rotation (0.02, 0.015, 0.0125, 0.01rad) and translation (7, 6.5, 6, 5mm) parameters. Starting from these maximum sizes, we solved Eq 3 running FastPD several times per level until no improvement is produced or a maximum number of iterations is achieved ([200, 100, 150, 600]), using different label space decreasing factors at every pyramid level (=[0.08, 0.07, 0.05, 0.03]). The total number of labels was fixed to 5 () for all the experiments. For the continuous case (where Eq 1 was optimized using simplex), we used again a 4levels pyramidal approach, with simplex running until convergence in every level. Finally, for the refined experiment, we just ran the simplex experiment initialized with the solution estimated with the discrete method. For every registration case, continuous approach took around 30secs while the discrete version took 9mins, running on a laptop with an Intel i74720HQ and 16GB of RAM.
3.2 Experiments
We performed two different type of experiments, considering individual registration cases as well as image series. For validation, we measured three different indicators: the distance in terms of translation and rotation parameters between the estimated and ground truth transformations, together with the mean of absolute differences (MAD) between the input 2D image and the slice specified by the estimated rigid transformation.
Method  MAD  

Simplex  0,14  0,13  0,12  8,46  9,62  10,69  51,88 
Discrete  0,12  0,08  0,09  5,87  6,72  6,18  42,36 
Refined  0,10  0,07  0,08  5,09  5,96  4,92  36,45 
Individual tests. The first set of experiments measures the accuracy of the three approaches using individual tests, where 100 random slices extracted from the 20 volumes are considered as single images (independently of the series), and registered to the first volume . We run the same experiment for every slice using three different initializations (resulting in 300 registration cases), where ground truth parameters were randomly perturbed in three different ranges ([5, 12), [12, 18), [18, 25) millimeters for translation and [0.1 , 0.2), [0.2, 0.3), [0.3, 0.4) radians for rotation parameters) to guarantee that both, good and bad initializations, are considered for every slice. Quantitative results are reported in Figures 3, 4 and summarized in Table 1. Visual results for qualitative evaluation are reported in Figure 2.
Results in the scatter plot from Figure 3 indicate that, as we go farther away from the initialization (in this case, it is quantified by the MAD between the input 2D image and the slice corresponding to the initialization), discrete and refined methods tend to be more robust. This robustness is clearly reflected by the slope of the trend lines: the refined method presents the trend line with the lower slope, meaning that even for bad initializations it converges to better solutions. The boxplot from Figure 4 and the numerical results from Table 1 confirm that discrete and refined methods perform better not only in terms of MAD, but also with respect to the distance between the rotation/translation estimated and ground truth parameters.
Temporal series test. The idea behind the second experiment is to simulate an image guided surgery (IGS) scenario, where a fixed preoperative volume must be fused with consecutive intraoperative 2D images suffering deformations (in this case, due to heart beating). Given the temporal sequence of 20 volumetric MRI images , we generated a sequence of 20 2D slices to validate our method. It was extracted as in [14]: starting from a random initial translation and rotation , we extracted a 2D slice from the initial volume . Gaussian noise was added to every parameter in order to generate the position used to extract the next slice from the next volume. We used for the rotation and for the translation parameters. All the slices were registered to the first volume . The solution of the registration problem for slice was used as initialization for the slice . The first experiment was initialized randomly perturbing its ground truth transformation with the same noise parameters. As it can be observed in Figure 5, discrete and refined approaches manage to keep a good estimation error while simplex can not. Note that different strategies could be used in real scenarios to obtain good initializations for the first slice. For example, in IGS, physicians could start from a plane showing always the same anatomical structure, or identify landmark correspondences in the first slice and the 3D image useful to estimate an initial transformation.
4 Discussion, Conclusions and Future Works
In this paper we presented a strategy to solve rigid slicetovolume registration as a discrete graph labeling problem, following the discretization strategy introduced by [3]. We validated our proposal using a MRI dataset of a beating heart, where arbitrary 2D slices are fused with a 3D volume. The experimental results showed that our discrete approach produces more accurate and robust estimates for the rigid transformations than a continuous method based on simplex. Moreover, they also reflected that results obtained using such a method can be further refined using a continuous approach like simplex, leading to even more accurate estimations. This is coherent with the conclusions presented by [20] for the case of optical flow estimation.
An interesting discussion about the limitations of our approach, emerges when we observe the results obtained in previous work by [14, 15, 16] using similar images. In these works, both rigid and local deformable parameters are estimated in a oneshot discrete optimization procedure, delivering results which are considerably better than ours, even for the refined approach. Since we are dealing with 2D images which are deformed with respect to the static 3D volume (due to heart beating), estimating both rigid and deformable parameters at the same time seems to be the correct solution since there is a clear mutual dependence between them. However, if we look at the results corresponding to the first slices of the temporal series in Figure 5 (where there is almost no deformation, and even null deformation for the 1st slice), we can see that the quality of the solution is significantly better than in the other cases. In fact, the error is almost 0. It suggests that when the 2D image is not deformed with respect to the input volume, our method is enough to capture slicetovolume mapping. This limitation is somehow inherent to the model we are using: rigid transformations can not deal with local deformations. To improve the results in these cases, we plan to extend our approach to linear transformations where also anisotropic scaling and shearing can be considered. Following the strategy by [3], it will result straightforward.
Finally, a future line of research has to do with applying discrete rigid (or linear) slicetovolume registration to other problems. As discussed in Section 1.1 motion correction for volume reconstruction is another problem requiring mapping slicetovolume. It would be interesting to explore how our approaches performs in this case.
References
 [1] Wang, C., Komodakis, N., Paragios, N.: Markov Random Field modeling, inference & learning in computer vision & image understanding: A survey. CVIU 117(11) (2013) 1610–1627
 [2] Paragios, N., Ferrante, E., Glocker, B., Komodakis, N., Parisot, S., Zacharaki, E.I.: (Hyper)graphical models in biomedical image analysis. Medical Image Analysis (2016)
 [3] Zikic, D., Glocker, B., Kutter, O., Groher, M., Komodakis, N., Kamen, A., Paragios, N., Navab, N.: Linear intensitybased image registration by Markov random fields and discrete optimization. Medical image analysis 14(4) (aug 2010) 550–62
 [4] Zikic, D., Glocker, B., Kutter, O., Groher, M., Komodakis, N., Khamene, A., Paragios, N., Navab, N.: Markov random field optimization for intensitybased 2D3D registration. In: SPIE Medical Imaging, International Society for Optics and Photonics (2010) 762334
 [5] Liao, R., Zhang, L., Sun, Y., Miao, S., Chefd’Hotel, C.: A Review of Recent Advances in Registration Techniques Applied to Minimally Invasive Therapy. IEEE TMM 15(5) (2013)
 [6] Bao, P., Warmath, J., Galloway, R., Herline, A.: Ultrasoundtocomputertomography registration for imageguided laparoscopic liver surgery. Surgical Endoscopy 19 (2005) 424–429
 [7] Fei, B., Duerk, J.L., Boll, D.T., Lewin, J.S., Wilson, D.L.: Slicetovolume registration and its potential application to interventional MRIguided radiofrequency thermal ablation of prostate cancer. IEEE Transactions on Medical Imaging 22(4) (2003) 515–525
 [8] Huang, X., Moore, J., Guiraudon, G., Jones, D.L., Bainbridge, D., Ren, J., Peters, T.M.: Dynamic 2D ultrasound and 3D CT image registration of the beating heart. IEEE TMI 28(8) (2009) 1179–1189
 [9] Rousseau, F., Glenn, O.: A novel approach to high resolution fetal brain MR imaging. Medical Image Computing and ComputerAssisted Intervention (MICCAI) 3749 (2005) 548–555
 [10] Kim, B., Boes, J.L., Bland, P.H., Chenevert, T.L., Meyer, C.R.: Motion correction in fMRI via registration of individual slices into an anatomical volume. Magnetic Resonance in Medicine 41(5) (1999) 964–972
 [11] Sotiras, A., Davatazikos, C., Paragios, N.: Deformable Medical Image Registration: A Survey. IEEE Transactions on Medical Imaging 32 (2013) 1153–1190
 [12] Glocker, B., Sotiras, A.: Deformable Medical Image Registration: Setting the State of the Art with Discrete Methods. Annual review of biomedical engineering 13 (2011) 219–244
 [13] Heinrich, M.P., Jenkinson, M., Brady, M., Schnabel, J.A.: MRFBased deformable registration and ventilation estimation of lung CT. IEEE TMI 32(7) (2013) 1239–1248
 [14] Ferrante, E., Paragios, N.: Nonrigid 2D3D Medical Image Registration Using Markov Random Fields. In: MICCAI. (2013) 163–170
 [15] Ferrante, E., Fecamp, V., Paragios, N.: Implicit planar and inplane deformable mapping in medical images through high order graphs. In: ISBI 2015. (apr 2015) 721–724
 [16] Ferrante, E., Fecamp, V., Paragios, N.: Slicetovolume deformable registration: efficient oneshot consensus between plane selection and inplane deformation. In: IJCARS. Volume 10. (2015) 791–800

[17]
Gholipour, A., Estroff, J.A., Warfield, S.K.:
Robust SuperResolution Volume Reconstruction From Slice Acquisitions: Application to Fetal Brain MRI.
IEEE TMI 29(10) (oct 2010)  [18] Tadayyon, H., Lasso, A., Kaushal, A., Guion, P., Fichtinger, G.: Target motion tracking in MRIguided transrectal robotic prostate biopsy. IEEE TBE 58(11) (2011) 3135–42
 [19] Birkfellner, W., Figl, M., Kettenbach, J., Hummel, J., Homolka, P., Schernthaner, R., Nau, T., Bergmann, H.: Rigid 2D/3D slicetovolume registration and its application on fluoroscopic CT images. Medical Physics 34(1) (2007) 246
 [20] Lempitsky, V., Roth, S., Rother, C.: FusionFlow: Discretecontinuous optimization for optical flow estimation. CVPR (2008) 1–22
 [21] Komodakis, N., Tziritas, G., Paragios, N.: Performance vs computational efficiency for optimizing single and dynamic MRFs: Setting the state of the art with primaldual strategies. Computer Vision and Image Understanding 112(1) (oct 2008) 14–29
 [22] Komodakis, N., Tziritas, G., Paragios, N.: Fast, approximately optimal solutions for single and dynamic MRFs. CVPR (2007)
 [23] Boykov, Y., Veksler, O., Zabih, R.: Fast approximate energy minimization via graph cuts. IEEE Transactions on Pattern Analysis and Machine Intelligence 23(11) (2001) 1222–1239
 [24] Boykov, Y., Kolmogorov, V.: An experimental comparison of mincut/max flow algorithms for energy minimization in vision. IEEE TPAMI 26(9) (sep 2004) 1124–1137
 [25] Nelder, J.A., Mead, R.: A Simplex Method for Function Minimization. The Computer Journal 7(4) (jan 1965) 308–313
 [26] Xu, R., Athavale, P., Nachman, A., Wright, G.A.: Multiscale Registration of RealTime and Prior MRI Data for ImageGuided Cardiac Interventions. IEEE TBE (2014)
 [27] Park, H., Meyer, C.R., Kim, B.: Improved Motion Correction in fMRI by Joint Mapping of Slices into an Anatomical Volume. MICCAI (2004) 745–751
Comments
There are no comments yet.