Adeline Paiement

Majid Mirmehdi

Senior Member, IEEE Xianghua Xie
email: x.xie@swansea.ac.uk.

Mark C K Hamilton
email: mark.hamilton@uhbristol.nhs.uk

Registration and Modeling from Spaced and Misaligned Image Volumes

Keywords: Modeling methodologies, registration, segmentation, shape interpolation, level set methods. !

come

INTRODUCTION

W E present an approach for object modeling from sin- gle patient's 3D or 3D+T medical volume made up of several misaligned sequences of spaced 2D images. Sparse datasets are commonly produced to decrease the acquisition time, for example when imaging infants or requiring breath holds. Misalignments may occur between sequences produced at the same visit due to patient movements between their independent acquisitions. If, in addition, the subjects are required to repeatedly hold their breath, as is usual in cardiac MRI (a major application of our method), different amounts of air may remain in the lungs, inducing different positions of the heart in the thorax 1 . The combined effects of patient movement and inconsistent breath holds lead to rigid misalignments between the image sequences, as illustrated in Fig. 1. Object modeling from such data is typically used to e.g. establish a diagnosis or plan surgery, and it requires solving the three closely intertwined tasks of registration, segmentation, and interpolation.

We identify three cases of sequences which may be misaligned due to patient movement between successive acquisitions. Stacks of slices (where "slice" denotes the location of a 2D image in a 3D volume) may be acquired with different scanner parameters in order to highlight different properties of the object, as illustrated in Fig. 2b. They may also be acquired with different slice orientations in order to obtain satisfactory views of different aspects of the object, as in Fig. 2c. In these two cases, each of the many sequences that make up the data comprises a set of slices, acquired at the same time, and spans the whole volume. We refer to such sequences as "spatial sequences". "Temporal sequences" may also be encountered, where one sequence is defined as all the time-frames of a single slice, as in Fig. 2a. In the case of cardiac MRIs, the acquisitions of the different temporal sequences (i.e. different slice positions) are synchronized by an electrocardiogram. Therefore, the sequences can be assumed to be temporally aligned with each other, i.e. frame t in all sequences represents the same phase t in the cardiac cycle. In addition, all the images of a sequence, regardless of the type of sequence, are acquired in a single measure (and breath hold if needed), and the patient should be immobile during this short acquisition time. Thus, we can reasonably assume that the images of the same sequence are spatially aligned. However, spatial misalignments may still exist between these independently acquired sequences. Therefore, the data needs to be registered prior to any segmentation or interpolation. This registration is challenging, because very diverse numbers, positions, and orientations of the slices may be encountered in order to account for the variety of imaged objects and of patient morphologies -see examples in Fig. 3 -and the registration method must be able to deal with arbitrary spatial configurations. In addition, the images of such datasets contain different parts of the object of interest and may have only a very limited number of intersections when their sequences have different orientations. In such cases, intensity information at slice intersections does not capture the geometry of the object, and it is therefore not reliable for driving registration, as will be demonstrated in Section 4. We argue that 3D segmentation information would provide the shape information that is needed for a robust registration. So, while 3D segmentation requires registered images, registration may also benefit from 3D segmentation.

The segmentation process itself is difficult, as it seeks to recover the 3D or 4D shape of the object from spaced 2D images. Independent segmentations of the images are not suitable since these segmentations may not be consistent with each other. Thus the 2D images are better segmented simultaneously in 3D or 4D. However, they commonly have relatively large spacings in order to shorten the acquisition and reduce patient discomfort, resulting in incomplete datasets as shown in Fig. 3. Such missing data do not offer support for 3D or 4D segmentation, and thus interpolation is necessary, either on the 2D images or on the shape of the object being segmented, after corrections of any misalignment. We showed in [START_REF] Paiement | Integrated segmentation and interpolation of sparse data[END_REF], [START_REF]Simultaneous level set interpolation and segmentation of short and long-axis MRI[END_REF] that combining 2D segmentation and shape interpolation produces more robust 3D segmentations than the sequential interpolation and segmentation approach. Clearly, the success of each of these three stages (registration, segmentation, and interpolation) depends on the robustness of the others.

We propose to model objects from spaced and misaligned volumes by combining registration, segmentation, and shape interpolation in a level set framework. The novel contributions of this work are: 1) a new registration method for internal alignment of volumes made up of spaced 2D slices having arbitrary relative positions and orientations and containing different parts of the object of interest. This method exploits segmentation results rather than pixel intensity, to better deal with largely spaced images that have too few intersections between them to allow computing reliable similarity measures. It also accounts for the global geometry of the object of interest, hence being more robust to large misalignments and local minima. 2) the full integration of the three stages of registration, segmentation, and shape interpolation into a novel and general level set framework that allows handling spaced and misaligned data with an increased robustness for all three processes. To the best of our knowledge, this is the first time that the three processes are integrated in a level set framework.

Next in Section 2, we review previous works on the analysis of spaced and misaligned image data, with particular reference to integrated methods. Then we introduce our method in Section 3. Section 4 presents validation tests and comparative results, and Section 6 concludes the paper.

PREVIOUS WORKS

Sequential approach -Traditionally, registration, segmentation, and interpolation, are performed sequentially. Usually, the different sequences of a dataset are registered first. Notably, registration based on similarity measures such as the Mutual Information (MI) and Normalized Mutual Information (NMI) [START_REF] Studholme | An overlap invariant entropy measure of 3D medical image alignment[END_REF] is very popular. Lotjonen et al. [START_REF] Lotjonen | Correction of movement artifacts from 4-D cardiac short-and long-axis MR data[END_REF] used NMI to register by translation two orthogonal sets of largely spaced temporal sequences (i.e. 2D slices) in cardiac MRIs. This method performs in two iterative steps: first, one slice is chosen randomly, and second, it is shifted in the direction which increases its similarity measure with all of its intersecting slices. It yields satisfactory results provided that the number of intersections between the different slices is high enough to compute meaningful similarity measures. Similarly, Slomka et al. [START_REF] Slomka | Patient motion correction for multiplanar, multi-breath-hold cardiac cine MR imaging[END_REF] and Elen et al. [START_REF] Elen | 3D breath hold related motion correction of multi-view 2D functional MRI[END_REF] registered, by in-plane translations and out-of-plane rotations, and by 3D translations, respectively, cardiac MRIs that comprise two orthogonal long-axis (LA) slices and a stack of short-axis (SA) slices, by optimizing the match of intensity profiles along the intersections of the various slices. Kim et al. [START_REF] Kim | Intersection based motion correction of multisclice MRI for 3-D in utero fetal brain image formation[END_REF], [START_REF] Kim | Bias field inconsistency correction of motion-scattered multislice MRI for improved 3D image reconstruction[END_REF] also optimized intensity profile matchings to rigidly align fetal brain MRI slices that form three orthogonal stacks. Woo et al. [START_REF] Woo | Reconstruction of high-resolution tongue volumes from MRI[END_REF] performed MI based affine registration followed by deformation on three volumes formed by orthogonal stacks of slices. Li and Denney [START_REF] Li | Registration of short and long-axis images in cine cardiac MRI[END_REF] registered intersecting and largely spaced 2D slices using the segmentation contours obtained from their preliminary independent segmentations, with the aforementioned drawback that the independent automatic 2D segmentations may not be consistent with each other.

Following registration, modeling an object from sets of spaced images requires interpolating the gaps in the data. Li et al. [START_REF] Li | 3D prostate boundary reconstruction from 2D TRUS images[END_REF] and Pihuit et al. [START_REF] Pihuit | Reconstruction implicite de surfaces 3D à partir de régions 2D dans des plans parallèles[END_REF], amongst others, proposed to analyze sets of parallel slices by first segmenting the slices independently, and then interpolating a surface from the 2D contours. More commonly, the image slices are first interpolated in order to reconstruct a 3D volume, e.g. as in [START_REF] Kim | Intersection based motion correction of multisclice MRI for 3-D in utero fetal brain image formation[END_REF], [START_REF] Kim | Bias field inconsistency correction of motion-scattered multislice MRI for improved 3D image reconstruction[END_REF], [START_REF] Woo | Reconstruction of high-resolution tongue volumes from MRI[END_REF], [START_REF] Yuan | Fusion of multi-planar images for improved three-dimensional object reconstruction[END_REF], [START_REF] Cordero-Grande | A Markov Random Field approach for topology-preserving registration: application to object-based tomographic image interpolation[END_REF]. Kim et al. [START_REF] Kim | Intersection based motion correction of multisclice MRI for 3-D in utero fetal brain image formation[END_REF] addressed the problem of super-resolution from anisotropic voxel data by fusing slice information using a point spread function weighted averaging method that selects information based on estimated local blurring. The problem of the interpolation between images having different orientations is ill-posed and, most other image interpolation works apply only to one or two stacks of parallel slices. In addition, while methods such as [START_REF] Kim | Intersection based motion correction of multisclice MRI for 3-D in utero fetal brain image formation[END_REF] work under the assumption that image slices have identical contrasts, this assumption does not hold in general cases [START_REF]Simultaneous level set interpolation and segmentation of short and long-axis MRI[END_REF], [START_REF] Kim | Bias field inconsistency correction of motion-scattered multislice MRI for improved 3D image reconstruction[END_REF], [START_REF] Woo | Reconstruction of high-resolution tongue volumes from MRI[END_REF]. It has been demonstrated in [START_REF]Simultaneous level set interpolation and segmentation of short and long-axis MRI[END_REF] and [START_REF] Woo | Reconstruction of high-resolution tongue volumes from MRI[END_REF] that the interpolation of intensities may be very sensitive to differences of gain and contrast in the images, thus producing interpolation artefacts which may bias the subsequent segmentation. Thus, [START_REF] Kim | Bias field inconsistency correction of motion-scattered multislice MRI for improved 3D image reconstruction[END_REF], [START_REF] Woo | Reconstruction of high-resolution tongue volumes from MRI[END_REF] introduced an extra step of intensity normalization before the volume interpolation.

This sequential approach lacks robustness, since a failure at an earlier stage is cascaded to all subsequent stages. Thus, a few attempts have been made to combine two or three of these stages in a single framework.

Integrated registration and segmentation -Priorknowledge based segmentation usually requires the registration of a shape model to the image being segmented. Guyader and Vese [START_REF] Guyader | A combined segmentation and registration framework with a nonlinear elasticity smoother[END_REF], and Rousson, Paragios, and coworkers [START_REF] Rousson | Shape priors for level set representations[END_REF], [START_REF] Paragios | Matching distance functions: A shape-to-area variational approach for global-to-local registration[END_REF] achieved this in a level set framework.

Registration and co-segmentation of overlapping images have been combined, mostly using statistical frameworks, to update the two stages simultaneously, e.g. [START_REF] Xiaohua | Simultaneous segmentation and registration for medical image[END_REF] and [START_REF] Gong | Texture-based simultaneous registration and segmentation of breast DCE-MRI[END_REF] amongst others, or iteratively in turn, e.g. [START_REF] Sun | Integrated registration of dynamic renal perfusion MR images[END_REF]. Yezzi et al. [START_REF] Yezzi | A variational framework for integrating segmentation and registration through active contours[END_REF] presented a level set based method to register and cosegment two overlapping images or dense volumes simultaneously by minimizing an energy containing a region-based segmentation term and an alignment term. This method does not apply to our case of spaced images which contain different parts of an object of interest that may not have the same 2D geometry and topology. However, its integration strategy inspired our proposed integrated registration method for spaced data, as explained in Section 3.1. Prevost et al. [START_REF] Prevost | Joint co-segmentation and registration of 3D ultrasound images[END_REF] developed a similar method for co-segmenting the kidney from contrast-enhanced and 3D ultrasounds using additional prior-knowledge constraints on the recovered shape. Unal et al. [START_REF] Unal | Coupled PDEs for non-rigid registration and segmentation[END_REF] proposed a variation of the work in [START_REF] Yezzi | A variational framework for integrating segmentation and registration through active contours[END_REF] for co-segmentation and regularized non-rigid registra-tion of two volumes.

Integrated segmentation and shape interpolation -Methods which register and deform a model rely on the continuity of the model in the gaps between images. This can be viewed as simultaneous segmentation and shape interpolation. Amongst others, [START_REF] Van Assen | A 3-D active shape model driven by fuzzy inference: Application to cardiac CT and MR[END_REF] used Active Shape Models (ASMs), [START_REF] Zambal | A fully automatic system for segmentation and analysis of the left and right ventricles of the heart using a bi-temporal two-component model[END_REF] combined 2D Active Appearance Models (AAMs) and a 3D shape model, and [START_REF] Lotjonen | Evaluation of cardiac biventricular segmentation from multiaxis MRI data: A multicenter study[END_REF] used a triangulated surface model. These methods require a training phase and lack flexibility.

In our previous works in [START_REF]Simultaneous level set interpolation and segmentation of short and long-axis MRI[END_REF] and [START_REF] Paiement | Integrated segmentation and interpolation of sparse data[END_REF], we proposed to integrate segmentation and shape interpolation in a level set framework in order to benefit from the accuracy and flexibility of level set methods. The level set segmenting surface was interpolated in the gaps in the data using a Mean Curvature Flow (MCF) in [START_REF]Simultaneous level set interpolation and segmentation of short and long-axis MRI[END_REF] and Radial Basis Functions (RBFs) in [START_REF] Paiement | Integrated segmentation and interpolation of sparse data[END_REF]. The shape interpolation strategy of [START_REF] Paiement | Integrated segmentation and interpolation of sparse data[END_REF] will be exploited in our proposed, fully integrated, framework.

When dealing with misaligned data, these integrated segmentation and shape interpolation methods require an independent registration to be performed beforehand. Such semi-sequential approaches with a separate registration stage will be shown to lack robustness, and a full integration of all three stages would be desirable. To the best of our knowledge, this has only been attempted by one other work to date which we describe next.

Integrated registration, segmentation, and shape interpolation -Zambal et al. [START_REF] Zambal | Improving segmentation of the left ventricle using a two-component statistical model[END_REF] registered, segmented, and interpolated 3D cardiac MRI data made up of spaced parallel SA slices using a two component deformable model, relying on the continuity of this model to interpolate in the gaps. The model components were updated iteratively: 2D AAMs segmented the individual slices, and a global shape model linked the 2D AAMs together. Registration was performed at every iteration by aligning the AAMs on the global shape model. We also propose to perform registration, segmentation, and shape interpolation simultaneously, but unlike Zambal et al. [START_REF] Zambal | Improving segmentation of the left ventricle using a two-component statistical model[END_REF] who used parametric models and prior knowledge, we combine the three stages in a level set framework, offering greater flexibility and accuracy, while requiring no prior knowledge.

PROPOSED METHOD

The registration method of [START_REF] Yezzi | A variational framework for integrating segmentation and registration through active contours[END_REF] presented in Section 2 applies to the co-segmentation of dense and overlapping data only. However, we will show in Section 3.1 that this approach has physical grounds which we will use first in Section 3.2 to propose a variation that is more robust to local minima, and then in Section 3.3 to derive a new and more general registration method suited for the registration of slices with arbitrary spatial configurations and containing different parts of an object. This new registration method will be integrated in the segmentation and shape interpolation framework of [START_REF] Paiement | Integrated segmentation and interpolation of sparse data[END_REF] in Section 3.4.

Justification for the Registration Approach

As explained in Section 1, this work addresses misalignments that may be approximated by rigid motions. Let us examine in this case the action of the registration of [START_REF] Yezzi | A variational framework for integrating segmentation and registration through active contours[END_REF]:

       ∂C ∂t = f N + |g | f • g N -κN (1a) ∂g i ∂t = C f • g ∂g ∂gi , (g ) -1 |g | T N (1b)
where C is the level set contour, N its normal, κ its curvature, and f and f may be seen as its speed in the target and registered image spaces respectively. g is the mapping from the registered image to the target image, with g i its i th parameter, and |g | is the determinant of g the Jacobian matrix of g with respect to the spatial position. Equation (1a) makes the contour evolve under the sum of the contributions of the two images. When aligning by translation T only, (1b) becomes

∂T i ∂t = C u i , f • g N , (2)
where T i is the i th component of the alignment vector T, and u i is the unit vector oriented in the i th direction. f • g (x) N (x) is the displacement imposed on the contour at point x by the segmentation of the registered image. So (2) makes T evolve by the total segmentation displacement integrated on the contour. When aligning by rotation around one axis only, for example the z-axis, with the rotation matrix being

R θ =   cos θ -sin θ 0 sin θ cos θ 0 0 0 1   , (3)
then, using trigonometric identities, (1b) becomes

∂θ ∂t = C R θ x, f • g R θ N = C R π 2 +θ x, f • g R θ N = C x, f • g R -π 2 N . (4)
x, f • g (x) R -π 2 N (x) may be seen as the torque of a force applied to the contour at point x, directed tangentially, and with an intensity of f • g (x). The sum of these torques along the contour C makes it rotate towards a better alignment with the registered image. Thus, the action of these two types of alignment (translation and rotation) can be understood as having physical grounds and it would make sense to apply them to more general cases than [START_REF] Yezzi | A variational framework for integrating segmentation and registration through active contours[END_REF], such as non-overlapping images in a spaced volume.

Domain Extension of Image Force

In the original registration method of [START_REF] Yezzi | A variational framework for integrating segmentation and registration through active contours[END_REF] expressed in (1b), only the points located on the level set contour are used to update the image registration by summing the pulls and torques of their "segmentation force" identified in Section 3.1. Since registration relies exclusively on these segmentation information, this implies that, before the contour converges on the edges of the object, the global geometry of the object is unknown 2 , as illustrated in Fig. 4b. This makes 2. After convergence, the contour has the shape of the object, thus the segmentation information does contain information on the global geometry of the object. the registration more likely to get trapped in local minima, especially in cases of high initial misalignments.

We propose to address this issue by exploiting the "segmentation force" for registration in areas that provide a complementary description of the geometry of the object. Thus, we use the speed f • g both on the contour and in the domain

Ω diff Ω diff = x ∈ I| f • g (x) • φ (x) < 0 . (5)
Ω diff contains the relevant speeds f • g that describe the global geometry of the object in the registered image as well as the current geometry of C, as seen in Fig. 4c. Equation (1b) then becomes

∂g i ∂t = C∪Ω diff f • g ∂g ∂g i , (g ) -1 |g | T N . (6)
As will be shown in Section 4, this use of additional segmentation information in Ω diff , that describe the geometry of the object of interest even before it is matched by C, reduces the number of local minima which may bias the registration. We refer to this variant of the registration method as the "global variant", as opposed to the original "local variant" which minimizes [START_REF] Yezzi | A variational framework for integrating segmentation and registration through active contours[END_REF]'s co-segmentation energy and uses local information of the contour only. Note that "global" here refers to the domain rather than to a "global solution" of the optimization problem in [START_REF] Yezzi | A variational framework for integrating segmentation and registration through active contours[END_REF]. Instead, this global variant may be seen as the minimization of an energy that penalizes a poor match between the shapes of the contour C and of the object in the image, similar to the energy used in [START_REF] Rousson | Shape priors for level set representations[END_REF], [START_REF] Paragios | Matching distance functions: A shape-to-area variational approach for global-to-local registration[END_REF] for rigidly aligning shapes represented by two distance functions:

E (g) = Ωreg (φ (x) -φ I (g (x))) 2 , (7)
where φ and φ I would be in our case the distance functions of contour C and of the shape in the image, respectively, and Ω reg = Ω diff . Indeed, by performing a similar examination of (7) as in Section 3.1, we find that

∂T i ∂t = Ωreg D(x) u i , ∇φ I (g(x)) ≈ Ωreg D(x) u i , -N I (8)
in the case of alignment by translation, and

∂θ ∂t = Ωreg D(x) x, R -π 2 R -θ ∇φ I (g(x)) ≈ Ωreg D(x) x, -R -π 2 N I (9)
for aligning by rotation, with D(x) = 2 (φ(x) -φ I (g(x)))

and N I the normal to the object contour in the reference space. Since the shape of the object of interest is unknown before segmentation is complete, its normal N I cannot be computed in our case. In our formulation, this is not necessary and registration relies on the normal to the segmenting contour N instead. Nevertheless, the link of our global variant of registration with registration-based energy minimization demonstrates its theoretical foundation.

The contour normal N may be computed following the usual formula for level sets 3 : N = -∇φ. However, as N is now used potentially in the whole volume rather than only on contour C, we must make sure that ∇φ is directed pointing inwards of the contour everywhere. We enforce this by adding the classic level set normalization term sign (φ) (1 -∇φ ) to the level set evolution function.

When registering cluttered images using the global variant, one must ensure that any background objects similar to the object of interest do not attract the contour and bias the registration. We enforce this by simply setting to zero any areas of positive speeds f • g which are not connected to the interior of the level set contour, so that only the object being segmented can drive registration. The local variant does not require this as it ignores the erroneous segmentation of spurious false object pixels in the background. This makes it more accurate, as was shown by our experiments.

Therefore, we recommend to use a combination of the two variants of registration in order to benefit from their two complementary advantages. The global variant should be used first to avoid local minima when there is significant misalignments. Then, when the level set contour matches the shape of the object sufficiently well, thus indicating a low risk of local minimum for registration, the local variant should be used to achieve best local fitting. We use the following shape matching condition between the contour C and the object of interest:

Ω ( f •g)>0 ∩ Ω φ>0 > |Ω diff | , (10)
where |.| denotes cardinality, i.e. the intersection between the shape of the contour C and the shape of the object is greater than the union of their relative complements.

Proposed Registration Method

Based on the analysis of the registration method in [START_REF] Yezzi | A variational framework for integrating segmentation and registration through active contours[END_REF] that we performed in Section 3.1, we propose to align internally, by translation and rotation, a 3D or 4D volume made up of any number of spaced 2D slices that have arbitrary relative positions and orientations. At the same time, we segment an object in the volume and interpolate its shape where no image information is available. For each sequence n of M n 2D images, we define a mapping function 3. We use here the convention φ(x) > 0 inside the contour.

g n (x) = R n x + T n (where R = R χ R ξ R ζ is the 3D rotation matrix)
from the common 3D space to the local 3D space associated to the sequence. The update of g n moves all M n images of sequence n in unison since, as explained in the introduction, all the images of the same sequence are considered as internally registered. The evolution of the level set function φ and of the mapping functions g n is given by our generalization of Yezzi et al. [START_REF] Yezzi | A variational framework for integrating segmentation and registration through active contours[END_REF]'s evolution equations for the case of registration by translation and rotation, i.e.

f = n f n • g n , (11)
                           ∂T n i ∂t = 1 |Ω n reg | Ω n reg (f n • g n ) u i , R n N ∂ζ n ∂t = 1 |Ω n reg | Ω n reg (f n • g n ) R χ n R ξ n R π 2 +ζ n x + T n , R n N ∂ξ n ∂t = 1 |Ω n reg | Ω n reg (f n • g n ) R χ n R π 2 +ξ n R ζ n x + T n , R n N ∂χ n ∂t = 1 |Ω n reg | Ω n reg (f n • g n ) R π 2 +χ n R ξ n R ζ n x + T n , R n N (12)
with Ω n reg = Ω n ∩ C for the local variant of the registration method, and Ω n reg = (Ω n ∩ C) ∪ Ω n diff for its global variant. Ω n is the domain of the images of sequence n within the common 3D space Ω. Adapting the definition of Ω diff in [START_REF] Slomka | Patient motion correction for multiplanar, multi-breath-hold cardiac cine MR imaging[END_REF],

Ω n diff is defined as Ω n diff = {x ∈ Ω n | (f n • g n ) (x) • φ (x) < 0} . (13)
Note that Yezzi et al. [START_REF] Yezzi | A variational framework for integrating segmentation and registration through active contours[END_REF] with rigid registration is a special case of [START_REF] Pihuit | Reconstruction implicite de surfaces 3D à partir de régions 2D dans des plans parallèles[END_REF], with the local variant of Ω n reg , two images or volumes n 1 and n 2 , and Ω n1 = Ω n2 .

In [START_REF] Yezzi | A variational framework for integrating segmentation and registration through active contours[END_REF], the amount of displacement applied to the registered image is proportional to its remaining misalignment with the target image. Indeed, in (2) and ( 4) the amount of misalignment between the registered image and the contour C is an implicit approximation of the misalignment between the two images in [START_REF] Yezzi | A variational framework for integrating segmentation and registration through active contours[END_REF]'s case of co-segmentation through registration based regularization. Although this proportionality of displacement with misalignment may be seen as a way to improve convergence speed, it takes a much more fundamental role in the registration process in our case of aligning multiple (more than two) non-overlapping images. It inactivates or activates registration for images that are well-aligned or misaligned respectively. We noticed in our experiments that this avoids the former to be sent away, jeopardizing the whole registration and segmentation process. It also avoids the latter to stay immobile, giving the segmenting contour C the opportunity to reach its segmentation optimum in most of the image although its position is not consistent with the rest of the dataset. This would block the image to its misaligned position and also endanger the whole process. Thus, moving the sequences by amounts proportional to their remaining misalignment is essential to ensure registration behaves as expected.

However, in our case, the implicit approximation of misalignment of [START_REF] Yezzi | A variational framework for integrating segmentation and registration through active contours[END_REF] does not hold because our 2D images contain different parts of an object, and because very few points of the segmenting contour in the images of sequence n intersect other sequences. Thus, in our first experiments which relied on the implicit approximation of misalignment, registration suffered from the issues described previously and could not be performed. Therefore, we need to compute an explicit estimate of the misalignment λ n of sequence n. We propose to do so by exploiting the contour's speeds f n and f m (m = n) that contain information on the shape of the object of interest in both sequences. We examine their dissimilarity at the intersections between Ω n reg and Ω m reg . When using the local variant of the registration method, λ n = λ local n is estimated as

λ local n = Ω (f n •g n )•(f m •g m )<0 ∩ C |Ω n ∩ Ω m ∩ C| , (14)
and for the global variant, the estimation of

λ n = λ global n is λ global n = Ω (f n •g n )•(f m •g m )<0 ∩ (C ∪ Ω n diff ∪ Ω m diff ) Ω n ∩ Ω m ∩ C ∪ Ω n diff ∪ Ω m diff . (15)
These explicit estimations of misalignment λ n are then used to modulate the amount of registration, represented as the amounts of translation (in pixels) and rotation (in degrees). This modulation is achieved through the update of the norms of ∂T n ∂t and ∂θ n ∂t (with θ n = ζ n , ξ n , χ n T being the vector of rotation angles of R n ). Although ∂θ n ∂t only approximates the amount of rotation of R n , its computation is simpler than finding the rotation axis and angle of rotation of R n , and the scheme we present next in ( 17) and ( 19) was found to satisfactorily match the amount of rotation of a sequence n to its degree of misalignment.

We use piecewise linear schemes for the computation of the amounts of shift and rotation in order to respect the principle of proportionality of [START_REF] Yezzi | A variational framework for integrating segmentation and registration through active contours[END_REF]. We empirically established the following threshold values that provide a satisfactory behavior of our registration method for all our applications:

shift n =      1 λ n ≥ 0.6 linear 0.1→1 0.3 ≤ λ n < 0.6 min linear 0→0.1, ∂T n ∂t λ n < 0.3 , ( 16
)
rotation n =      r 1 λ n ≥ 0.6 linear r 2 →r 1 0.2 ≤ λ n < 0.6 min linear 0→r 2 , ∂θ n ∂t λ n < 0.2 , ( 17
)
with r 1 and r 2 representing rotations of 1 • and 0.1 • respectively around each axis:

r 1 = √ d • 1 • 2 and r 2 = √ d • 0.1 • 2 (d being
the number of degrees of freedom for rotation, i.e. the number of dimensions of the data). The modified shift and rotation vectors are then

∂T n ∂t = ∂T n ∂t ∂T n ∂t • shift n , ( 18
)
∂θ n ∂t = ∂θ n ∂t ∂θ n ∂t • rotation n . ( 19
)
It is desirable to limit the amount of deformation that the level set surface must undergo at each iteration in order to adapt to the displacements of the images. We achieve this by restricting translations and rotations to their minimal useful parts, by subtracting an offset to the update of each component g n i of g n . This offset is the minimum displacement ∂gi ∂t min , in absolute values, amongst all sequences, i.e.

offset i = ∂g i ∂t min , (20)
with ∂gi ∂t min being such that

∂g i ∂t min = min n ∂g n i ∂t . (21)
This offset subtraction has the additional advantage of avoiding the modeled object to drift out of the boundaries of the volume. Thus, it is not necessary to arbitrarily choose a sequence to remain fixed in the volume. Finally, T n and R n are updated as:

   T n t+1 = T n t + dt ∂T n ∂t -offset T θ n t+1 = θ n t + dt ∂θ n ∂t -offset θ , (22)
with offset T = offset Tx , offset Ty , offset Tz T and

offset θ = offset ζ , offset ξ , offset χ T .

Proposed Integrated Registration, Segmentation, and Interpolation

The integration of this registration method with the segmentation and shape interpolation framework of [START_REF] Paiement | Integrated segmentation and interpolation of sparse data[END_REF] is performed by providing the speed f of the contour C in [START_REF] Li | 3D prostate boundary reconstruction from 2D TRUS images[END_REF] to the level set evolution equations ( 23) and ( 24) of [START_REF] Paiement | Integrated segmentation and interpolation of sparse data[END_REF], in order to compute the evolution of the level set function φ both on the domains of the images and in the gaps between them:

dα i dt = C f (x) ψ i (x) ≈ Ω δ (φ (x)) f (x) ψ (x i -x) ≈ ((δ (φ) • f ) ψ) (x i ) , (23)
∂φ ∂t (x) = N i=1 dα i dt ψ(x -x i ) = dα dt ψ (x). ( 24
)
δ ε is an approximation of the Dirac function δ with width , and denotes convolution. ψ i , x i , and α i are the N RBFs, their centers, and their coefficients, introduced in [START_REF] Paiement | Integrated segmentation and interpolation of sparse data[END_REF] to interpolate φ:

φ(x) = N i=1 α i ψ i (x) = N i=1 α i ψ(x -x i ) . ( 25
)
As in [START_REF] Paiement | Integrated segmentation and interpolation of sparse data[END_REF] we use inverse multiquadric RBFs of the form

ψ(x) = x 2 + γ 2 -β 2
, with β being the number of dimensions of the data, to ensure the decrease of the RBF is not too steep, and γ being a variable which must be set by the user according to the data in order to fine tune the flatness of the RBF.

As seen in [START_REF] Li | 3D prostate boundary reconstruction from 2D TRUS images[END_REF], f combines all the speeds f n that are computed for the independent sequences n. They may be produced by any segmentation method of choice. The term κN in (1a) was dropped in [START_REF] Li | 3D prostate boundary reconstruction from 2D TRUS images[END_REF] because the level set framework of [START_REF] Paiement | Integrated segmentation and interpolation of sparse data[END_REF] does not require contour smoothing due to the intrinsic smoothing of the convolutions by the RBF.

As stated in Section 3.2, when using the global variant of the registration method, the classic level set normalization term needs to be used in order to ensure that N can be computed easily everywhere. This term does not depend on the availability of data and thus is defined everywhere and does not require interpolation. In addition, its smoothing by convolutions with the RBF is not desirable and may degrade its performance. So we add it to [START_REF] Van Assen | A 3-D active shape model driven by fuzzy inference: Application to cardiac CT and MR[END_REF] to update φ directly rather than incorporating it into f in [START_REF] Unal | Coupled PDEs for non-rigid registration and segmentation[END_REF]:

∂φ ∂t (x) = dα dt ψ (x) + sign (φ(x)) (1 -∇φ(x) ) . ( 26
)
We use the same timestep dt to update C and all g n . As in [START_REF] Paiement | Integrated segmentation and interpolation of sparse data[END_REF], it is set so that C moves by 1 pixel per iteration at the beginning of the process, and it is reduced progressively when oscillations of C are detected in order to allow a finer convergence of the registration. Convergence of the algorithm is reached at the same time as convergence of the level set segmentation, since any additional motion of an image sequence would force the segmentation to adapt to the new configuration of the data. The automatic detection of the convergence of a level set contour is not in the scope of this work, and we use the simplest scheme with fixed number of iterations in our implementation.

The overall algorithm for the proposed framework is outlined in Algorithm 1 and our software implementation is available on our project web page 4 .

EXPERIMENTAL RESULTS AND ANALYSIS

To assess and evaluate the proposed methodology, we first focus on the accuracy and robustness of the proposed registration approach, and then we evaluate the 3D modeling of the complete framework.

Registration methods that are designed to register two target and reference images containing the same object e.g. [START_REF] Guyader | A combined segmentation and registration framework with a nonlinear elasticity smoother[END_REF], [START_REF] Rousson | Shape priors for level set representations[END_REF], [START_REF] Paragios | Matching distance functions: A shape-to-area variational approach for global-to-local registration[END_REF], [START_REF] Xiaohua | Simultaneous segmentation and registration for medical image[END_REF], [START_REF] Gong | Texture-based simultaneous registration and segmentation of breast DCE-MRI[END_REF], [START_REF] Sun | Integrated registration of dynamic renal perfusion MR images[END_REF], [START_REF] Yezzi | A variational framework for integrating segmentation and registration through active contours[END_REF], [START_REF] Prevost | Joint co-segmentation and registration of 3D ultrasound images[END_REF], [START_REF] Unal | Coupled PDEs for non-rigid registration and segmentation[END_REF], are not suitable for comparison as they cannot align 2D images that contain different parts of a 3D volume and of an object of interest. The methods presented in [START_REF] Slomka | Patient motion correction for multiplanar, multi-breath-hold cardiac cine MR imaging[END_REF], [START_REF] Elen | 3D breath hold related motion correction of multi-view 2D functional MRI[END_REF] are limited to temporal sequences forming two orthogonal stacks of slices. The methods in [START_REF] Kim | Intersection based motion correction of multisclice MRI for 3-D in utero fetal brain image formation[END_REF], [START_REF] Kim | Bias field inconsistency correction of motion-scattered multislice MRI for improved 3D image reconstruction[END_REF], [START_REF] Woo | Reconstruction of high-resolution tongue volumes from MRI[END_REF] only apply to three orthogonal stacks of slices. We compare against the algorithm of [4] that was described in Section 2 and that is popular for aligning cardiac MRIs. Although this method was originally designed for two orthogonal stacks of temporal sequences and was applied to translation only, the same algorithm can be applied to more general cases. It is therefore suitable for comparison, and it is the only method we know of that can align slices that contain different parts of a 3D volume and object, and that have arbitrary spatial configurations. We refer to our extended use of [START_REF] Lotjonen | Correction of movement artifacts from 4-D cardiac short-and long-axis MR data[END_REF] as "SR [START_REF] Lotjonen | Correction of movement artifacts from 4-D cardiac short-and long-axis MR data[END_REF]" ("Sequential Registration [START_REF] Lotjonen | Correction of movement artifacts from 4-D cardiac short-and long-axis MR data[END_REF]"). The method is based on similarity measures and in [START_REF] Lotjonen | Correction of movement artifacts from 4-D cardiac short-and long-axis MR data[END_REF] it is implemented with NMI that is commonly used to register multi-modal images. We retain this similarity measure in our experiments, and, in Experiment 1 with uni-modal data, we also compare against for all image sequence n do 12:

Compute ∂T n ∂t and ∂θ n ∂t using [START_REF] Pihuit | Reconstruction implicite de surfaces 3D à partir de régions 2D dans des plans parallèles[END_REF] 13:

if global variant = true then 14:

Estimate misalignment λ n using [START_REF] Guyader | A combined segmentation and registration framework with a nonlinear elasticity smoother[END_REF] 15:

if condition (10) = true then Compute shift n and rotation n using ( 16) and [START_REF] Paragios | Matching distance functions: A shape-to-area variational approach for global-to-local registration[END_REF] 22:

Compute ∂T n ∂t and ∂θ n ∂t using ( 18) and [START_REF] Gong | Texture-based simultaneous registration and segmentation of breast DCE-MRI[END_REF] 23:

end for

24:

for all degree of freedom i do 25:

Compute offset i using ( 20) and ( 21)

26:
end for

27:

for all image sequence n do 28:

Update T n and R n using [START_REF] Prevost | Joint co-segmentation and registration of 3D ultrasound images[END_REF] 29:

end for 30: until convergence [START_REF] Lotjonen | Correction of movement artifacts from 4-D cardiac short-and long-axis MR data[END_REF] implemented with two other similarity measures that only apply to uni-modal data, namely Sum of Squared intensity Differences (SSD) and Normalized Cross-Correlation (NCC). These added tests aim at demonstrating that the lack of robustness that we find in [START_REF] Lotjonen | Correction of movement artifacts from 4-D cardiac short-and long-axis MR data[END_REF] is not due to the choice of similarity measure, but rather to the method relying on a local computation of a similarity measure.

We do not evaluate the quality of the segmentation, since this quality is highly dependent on the segmentation algorithm, which must be chosen by the user depending on the nature of the data. The choice of an optimum segmentation algorithm is not within the scope of this work.

The shape interpolation contribution of our framework was already evaluated in [START_REF] Paiement | Integrated segmentation and interpolation of sparse data[END_REF]. Indeed, the combined segmentation and interpolation of [START_REF] Paiement | Integrated segmentation and interpolation of sparse data[END_REF] is now used in our fully integrated framework. In [START_REF] Paiement | Integrated segmentation and interpolation of sparse data[END_REF], we compared this integrated method against a state-of-the-art sequential image interpolation [START_REF] Cordero-Grande | A Markov Random Field approach for topology-preserving registration: application to object-based tomographic image interpolation[END_REF] followed by 3D level set segmentation. This revealed that the combined approach, by taking into account the global geometry of the objects, produces more robust and accurate results than the sequential one.

We do not repeat this demonstration in this work, and instead we focus on the synergistic benefit of integrating registration into the framework of [START_REF] Paiement | Integrated segmentation and interpolation of sparse data[END_REF]. Therefore, when assessing the performance of the complete framework, inte- grating all stages of registration, segmentation, and interpolation, we compare against the semi-sequential approach which successively applies the registration SR [START_REF] Lotjonen | Correction of movement artifacts from 4-D cardiac short-and long-axis MR data[END_REF] based on a similarity measure, and the integrated segmentation and interpolation method of [1] 5 . We refer to this framework as "SR+ISI", where ISI is for "Integrated Segmentation and Interpolation". Based on the comparison results in [START_REF] Paiement | Integrated segmentation and interpolation of sparse data[END_REF], we can infer that the fully sequential method SR [START_REF] Lotjonen | Correction of movement artifacts from 4-D cardiac short-and long-axis MR data[END_REF] followed successively by image interpolation and 3D segmentation would produce less robust and accurate results than SR+ISI. Thus, the comparison against SR+ISI is fairer for the sequential approach.

Data and Evaluation Protocol

In our experiments, we distinguish between two cases:

• Spatial sequences, as in Figs. 2b and2c, that are made up of a stack of slices that move together during registration. In the rest of this article, we will denote such registration as "stack-wise" registration.

• Temporal sequences made up of a number of timeframes of a single slice position and orientation, as in Fig. 2a. In this case, the different slices which span the volume are registered with each other. We denote such registration as "slice-wise" registration.

The registration stage alone and the complete integrated framework are both evaluated quantitatively on artificial data in Section 4.2 and on artificial data made up of real images in Section 4.3. Section 4.4 presents qualitative results on real datasets. The datasets used in these experiments are presented in Table 1.

The artificial data used in Section 4.2 was created to allow the assessment of our registration method and of our integrated framework independently of the quality of segmentation. Thus, it is used to establish the limits of performance of our method in the case of a perfect segmentation. Indeed, as long as the segmentation is correct, the 5. Note, we did not compare with the fully integrated framework presented in [START_REF] Zambal | Improving segmentation of the left ventricle using a two-component statistical model[END_REF] since the authors informed us that their code and prior (AAM and 3D shape) models were no longer available. registration quality depends only on the geometry of the data. As explained in Section 4, the quality and robustness of segmentation depends strongly on the chosen segmentation algorithm and this choice is not within the scope of this work. We avoid the problem of finding a perfect segmentation algorithm by producing artificial data that are very easily segmented, i.e. images consisting of piecewise constant regions. We generate such images from a model of a left ventricle (LV) cavity of a heart and from a model of both ventricles of a neonatal brain that were created by segmenting, using the method of [START_REF] Paiement | Integrated segmentation and interpolation of sparse data[END_REF], two real MRI datasets which were visually inspected for no misalignments and requiring only a minimal amount of interpolation. Our artificial datasets were extracted from these two volumes using the positions and orientations of slices of 14 heart MRIs and 9 brain MRIs respectively. Therefore, although our artificial images are made up of piecewise constant regions, they contain the shape of real organs and have realistic positions and orientations. Since cardiac and brain MRIs usually comprise temporal and spatial sequences respectively, the artificial cardiac and brain datasets are used to evaluate the slicewise and stack-wise registration respectively. Only one timeframe was generated for the artificial cardiac data, which is the worst case for slice-wise registration as it restricts the amount of available information.

Our quantitative evaluation of Section 4.3 is on artificial data made up of real images, in order to provide quantitative evaluation on more plausible images and segmentation quality while still having registration and interpolation groundtruths. We create datasets with artificial gaps and misalignments by extracting spaced axial and sagittal slices from full-volume real CT-scans of hip and kidney which have no misalignments. Thus, the resulting spaced datasets are made up of real images but have artificial gaps.

In Sections 4.2 and 4.3, artificial initial misalignments are induced by initializing the registration vectors T n and θ n of the individual sequences to random values. A perfect registration should move all the sequences of a dataset to the same position, which we call "reference sequence position" (RSP). Since we do not use a target sequence, the RSP may ±(4.5e-3, 7.2e-4, 1.5e-4) ±(2e-5, 3.8e-5, 2.6e-5) ±(2e-3, 2.8e-3, 7e-4) take any value in the 6D space of translations and rotations.

T n i ini = rand [-20, 20] θ n i ini = rand [-20
In case of misalignments, we define the RSP as the mean position of the sequences. We then compute the registration error as the average and standard deviation of the distances of the sequences to the RSP. The complete framework is evaluated by computing the Jaccard coefficients of the shapes reconstructed from the initial misalignments. To do this, the groundtruth shape is first placed at the RSP. In Section 4.4, we use real datasets for which no registration and interpolation groundtruth is available, thus we assess the performance of the compared methods qualitatively. The datasets were obtained from hospital patients in real and routine clinical examinations. They include cardiac MRIs that comprise temporal sequences, and neonatal brain and hip MRIs that consist of spatial sequences.

We use the same segmentation algorithm for the proposed and sequential methods. As in [START_REF] Paiement | Integrated segmentation and interpolation of sparse data[END_REF], we use the piecewise constant model from [START_REF] Chan | Active contours without edges[END_REF] for our artificial images, and a method based on a piecewise model and a Parzen window estimator of pixel intensity for our real images. The level set surface is initialized as a small sphere, and its initial position initializes the segmentation models. For artificial data, since the parameters of the segmentation model are known constant intensities, the initial sphere is arbitrarily placed at the center of the volume. For real data, it is placed by the user by clicking at the center of the object on a single central image. Note that the initial sphere does not need to be positioned inside the object for all slices (and it is usually not, especially in case of large misalignments) as long as it allows a suitable initialization of the segmentation parameters. The heart datasets are processed with the scheme in two steps that was introduced in [START_REF] Paiement | Integrated segmentation and interpolation of sparse data[END_REF]. Other datasets are processed with the simple scheme. We refer the reader to [START_REF] Paiement | Integrated segmentation and interpolation of sparse data[END_REF] for details on both schemes. The schemes and values of RBF parameter γ are reported in Table 1.

Quantitative Evaluation on Artificial Data

Registration accuracy

Stack-wise registration -We perform two tests of stackwise registration on the 9 artificial brain datasets, with two different ranges of initial misalignments. In the first test, the sequences are initially randomly shifted by up to 10 pixels in all directions (T n i ini = rand [-10, 10]) and randomly rotated by up to 10 • along all axes (θ n i ini = rand [-10 • , 10 • ], with θ n i the i th component of θ n ). In the second test, we have [START_REF] Sun | Integrated registration of dynamic renal perfusion MR images[END_REF][START_REF] Sun | Integrated registration of dynamic renal perfusion MR images[END_REF] and In the first test, the SR [START_REF] Lotjonen | Correction of movement artifacts from 4-D cardiac short-and long-axis MR data[END_REF] method with NMI is slightly more accurate than the proposed registration method, but both methods achieve sub-pixel and sub-degree accuracy. We believe this small difference may be explained by the fact that in this "not so sparse" scenario (see Fig. 5) the many pixel intensity information that SR [START_REF] Lotjonen | Correction of movement artifacts from 4-D cardiac short-and long-axis MR data[END_REF] uses both inside and outside of the object may constrain it better than the 3D shape information used by our method. In the second test with more considerable initial misalignments, the overall performance of the SR [START_REF] Lotjonen | Correction of movement artifacts from 4-D cardiac short-and long-axis MR data[END_REF] method with all tested similarity measures is lower than the proposed registration method. This is due to the SR [START_REF] Lotjonen | Correction of movement artifacts from 4-D cardiac short-and long-axis MR data[END_REF] method (regardless of the similarity measure used) being sometimes trapped in local minima, even in the absence of any background object. This is illustrated in Fig. 5 where two different ventricles of the brain are confused as one single object and superimposed. Conversely, the proposed registration method exploits information on the global shape of the objects of interest, and therefore is more robust to such local minima and maintains a good accuracy.

T n i ini = rand [-
θ n i ini = rand [-20 • , 20 • ].
Slice-wise registration -We test slice-wise registration on the artificial cardiac datasets with alignment by translation in the SA plane first, and then by 3D translations and rotations, in both cases with a small then large range of initial misalignments. Table 3 presents the registration errors of these 4 tests, averaged on all 14 datasets.

The first two tests of SA plane translations are motivated by the fact that real cardiac MRIs acquired during a breath hold at end-expiration tend to suffer mainly from shifts in the SA plane, and have negligible shifts along the LA direction and rotations. For the first test with a small range of initial misalignments (T n xy ini = rand [-5, 5]), ±(0.17, 0.13, 0.17 while all methods achieve sub-pixel accuracy, the proposed registration method obtains results that are almost twice as precise as the best results of the SR [START_REF] Lotjonen | Correction of movement artifacts from 4-D cardiac short-and long-axis MR data[END_REF] method. The second test (second row of Table 3) is a quite extreme case, with shifts of up to 10 pixels, equivalent to around 20 mm given the image resolution. However, the performance of the proposed registration method remains as good as in the first test, while the accuracy of the SR [START_REF] Lotjonen | Correction of movement artifacts from 4-D cardiac short-and long-axis MR data[END_REF] method is strongly degraded for all 3 similarity measures. A similar behavior is observed for the tests with full 3D translations and rotations, although the overall accuracy of all methods is reduced in this more difficult registration case. For a small range of initial misalignments (third row of Table 3), the two methods achieve a similar accuracy, with the proposed registration method being slightly more accurate than the best results of the SR [START_REF] Lotjonen | Correction of movement artifacts from 4-D cardiac short-and long-axis MR data[END_REF] method. In the fourth row of Table 3 with a larger range of misalignments, the proposed registration method maintains a good accuracy while the average errors of the SR [START_REF] Lotjonen | Correction of movement artifacts from 4-D cardiac short-and long-axis MR data[END_REF] method are higher.

T n xy ini = rand [-5, 5] T n z ini = rand [-3, 3] θ n i ini = rand [-5 • , 5 • ] (1.
Note that the use of three different similarity measures (NMI, SSD, and NCC) did not significantly change the results of SR [START_REF] Lotjonen | Correction of movement artifacts from 4-D cardiac short-and long-axis MR data[END_REF] in all 4 tests. We believe that the differences in the results between SR [START_REF] Lotjonen | Correction of movement artifacts from 4-D cardiac short-and long-axis MR data[END_REF] and the proposed method on the artificial cardiac datasets is due to the fact that SR [START_REF] Lotjonen | Correction of movement artifacts from 4-D cardiac short-and long-axis MR data[END_REF] uses the intensity at the intersections of images to compute a similarity measure, and therefore requires a minimum number of intersection points in order to obtain meaningful similarity measures. Indeed, in [START_REF] Lotjonen | Correction of movement artifacts from 4-D cardiac short-and long-axis MR data[END_REF] the datasets were made up of two stacks of 6 and 8 parallel SA and LA slices respectively, both spanning the whole LV of the heart, thus the numbers of intersections in our datasets are significantly lower. On the contrary, the proposed registration method uses global 3D segmentation results rather than local pixel intensity, and therefore it can handle datasets where the LA slices are less numerous (e.g. 2 or 3) and cross the imaged object at various angles and locations.

We note in the last two rows of Table 3 and on Fig. 8 that the rotation accuracy of our method seems, visually and subjectively, less satisfactory than for its translationalthough always within a 3 • error. This may be explained by considering the nearly cylindrical shape of the modeled object which leaves more degrees of freedom to the rotation of the slices for recovering an approximated cylinder shape.

A difficulty which sometimes arises when an image intersects very few other images, is that a position ambiguity can happen, as illustrated in Fig. 6a, where two different positions can be viewed as equally correct. Such underconstrained situations happened several times during our slice-wise registration tests, and a few datasets had one or two of their LA slices caught in a local minimum at an ambiguous position, especially in the third and fourth tests. This resulted in a wider distribution of the slices around their mean position in these datasets, which explains the large average standard deviations presented in the last two rows of Table 3. Most of the time, this did not prevent the contour to recover the proper shape of the object since both ambiguous positions were acceptable in a topological point of view, as will be illustrated in Section 4.2.2.

The flatness of the RBF has been found to have an influence on the registration results. When the RBF is too sharp, the contour can sometimes fold around a slice, as in Fig. 6b, trapping the slice in its current misalignment. This problem can be solved by using a flatter RBF, as illustrated by Table 4 which presents average registration errors when aligning the artificial heart datasets in the SA plane using two different RBF flatness parameters (note that contrary to the previous tests, only one value for γ is used per processing). We see that the flat RBF (with γ = 1) avoids contour folding traps, and thus yields better results than the sharp RBF (with γ = 0.5). In the 4 previous tests, the segmentation and interpolation stages of the integrated framework required the use of a sharp and a flat RBF successively, therefore the second processing step with a flat RBF corrected any possible folding of the contour.

Complete framework accuracy

Table 6 shows the Jaccard coefficients averaged on all datasets for each test and both methods.

Stack-wise registration -In the first test with small initial misalignments, both the SR+ISI method, with NCC and NMI, and the proposed method recover the shape of the brain model satisfactorily, and yield similar mean Jaccard coefficients (first row of Table 6). SR+ISI with SSD performs less well due to its poorer registration. For larger initial misalignments (second row of Table 6), the SR [START_REF] Lotjonen | Correction of movement artifacts from 4-D cardiac short-and long-axis MR data[END_REF] method sometimes falls into local minima, regardless of the similarity measure, due to its use of local shape information only, as was illustrated in Fig. 5 and discussed in Section 4.2.1. In such cases, the resulting reconstructed shapes, shown in the second column of Fig. 7, suffer from the contradictory information provided by the misaligned sequences, thus significantly decreasing the mean Jaccard coefficients. On the contrary, the use of global shape information by the proposed registration method makes it more robust to such local minima and allows the framework to recover satisfactory shapes (third column of Fig. 7) and to maintain the same accuracy.

Slice-wise registration -The proposed framework recovers better shapes than SR+ISI in all four tests, due to its more accurate registration (4 last rows of Table 6). In particular, as seen in the last 4 rows of Table 6, it is not impaired by larger initial misalignments in tests 2 and 4 and yields mean Jaccard coefficients similar to those of tests 1 and 3 with smaller misalignments. In addition, in tests 3 and 4, the more significant number of degrees of freedom (3 translations and 3 rotations) makes the registration problem less constrained than in tests 1 and 2, and causes some slices to have several acceptable positions, as was discussed in Section 4.2.1 and illustrated in Fig. 6a. However, the proposed framework recovers quite accurate shapes even in such cases, as shown on the right of Fig. 8. Indeed, it manages to place the slices in positions which make them contribute to the overall shape of the object while limiting any conflicts with other slices. The SR+ISI method however, attempts to align the slices without taking into account the global shape of the object, and therefore obtains poor registration and shape recovery in such cases, as illustrated in the middle columns of Fig. 8.

Quantitative Evaluation on Artificial Datasets Made up of Real Images

In this experiment we show that our proposed framework is not limited to images consisting of piecewise constant regions, but can also deal with real images from a variety of modalities that may exhibit large variations in appearance and various levels of noise.

Registration accuracy

Stack-wise registration -We perform two tests, with a small and larger range of initial misalignments respectively. In the first test, only a small ROI around the object of interest is used by both methods in order to avoid possible local minima caused by background objects and thus test the accuracy of both methods. In the second test, the entire images are used in order to evaluate the robustness to local minima due to cluttered backgrounds. Table 5 presents the average registration errors sorted by slice spacings.

The results show the same trend as the corresponding experiment on artificial data. For small initial misalignments and slice spacings (first 5 rows of Table 5), both methods generally achieve sub-pixel and sub-degree accuracy. However, for larger spacings or misalignments (last 5 rows of Table 5), the SR [START_REF] Lotjonen | Correction of movement artifacts from 4-D cardiac short-and long-axis MR data[END_REF] method falls into local minima, especially when the image contains several objects of similar intensities, as illustrated in Fig. 9a. Therefore, its average accuracy in the second test is significantly reduced. The proposed registration method avoids such issues due to its global approach and yields satisfactory registration for all our datasets except one, which contains an object with a rather symmetrical shape (Fig. 9b). This symmetry causes an ambiguity on the pose of the sequences, similar to that discussed in Section 4.2 and illustrated in Fig. 6a. Thus, the method results in a registration error for this dataset of (1.76, 11.75, 2.16) pixels and (9.25 • , 0.23 • , 16.23 • ), but the recovered shape is acceptable, as shown in Fig. 9e. This case increases the overall error of the proposed registration method in test 2 to (0.32, 1.67, 0.58) pixels and (1.3 • , 0.2 • , 2.2 • ). If we discount it, the mean error becomes (0.08, 0.34, 0.49) pixels and (0.41 • , 0.23 • , 0.35 • ), thus the method can yield a similar accuracy than in the first test with smaller initial misalignments if no local minimum occurs due to shape symmetry.

Slice-wise registration -We initially shift all the slices of a dataset independently in the (x,y)-directions, by a random amount bounded by the slice spacing of the dataset: T n xy ini = rand [-SP, SP ] with SP the slice spacing, and T n z ini = 0 and θ n i ini = 0 • . Registration results are presented in Table 7. This test is different to the slice-wise registration tests on ar-

T n i ini = rand [-30, 30], θ n i ini = rand [-20 • , 20 • ]
5 pixels (0.001, 0.004, 0.12) (0.02 tificial data, since in Section 4.2 the datasets were made up of one stack plus a varying small number of roughly orthogonal slices having arbitrary positions and orientations, whilst here the datasets are composed of two orthogonal stacks of parallel slices. This spatial configuration contains more intersections between the slices and is the configuration for which the method in [START_REF] Lotjonen | Correction of movement artifacts from 4-D cardiac short-and long-axis MR data[END_REF] was originally designed. Therefore, both SR [START_REF] Lotjonen | Correction of movement artifacts from 4-D cardiac short-and long-axis MR data[END_REF] and the proposed registration method perform similarly well for slice spacings and initial misalignments of up to 15 pixels, achieving sub-pixel accuracy. In the last test, with 20 pixels spacing and maximum misalignments, the SR [START_REF] Lotjonen | Correction of movement artifacts from 4-D cardiac short-and long-axis MR data[END_REF] method encounters local minima problems when objects of the background are confused with the object of interest. At the same time, the number of intersection points becomes less suitable to compute similarity measures. As a result, its average registration error increases dramatically. The proposed registration method does not suffer from these problems and maintains a sub-pixel accuracy.

Complete framework

The Jaccard coefficients for both methods, averaged on all datasets for each test, are presented in Table 8. When the registration yielded by the SR+ISI method and the proposed framework are similarly accurate, as in the first stack-wise registration test, both frameworks reconstruct satisfactory shapes with similar Jaccard coefficients. However, when the SR [START_REF] Lotjonen | Correction of movement artifacts from 4-D cardiac short-and long-axis MR data[END_REF] method falls into local minima because it confuses background objects with the object of interest, as in Fig. 9a, the two sequences which make the dataset provide incoherent information to the segmenting level set surface, lowering the Jaccard coefficient to 0.12 in the case of Fig. 9d. The proposed framework does not suffer from this problem due to its use of global shape information. Indeed, as explained in Section 4.3.1, the only local minimum its registration encountered was caused by the symmetry of the object. However, the registration method takes into account the global shape of the object and attempts to reduce the discrepancy in the segmentation information that the images provide. Therefore, in such cases, the proposed framework still recovers a plausible shape, as illustrated in Fig. 9e where it obtained a Jaccard coefficient of 0.83. Thus, the mean Jaccard coefficient for the second stack-wise registration test remains as high as for smaller initial misalignments.

Similarly, in the slice-wise registration test, the SR [4] method encountered local minima problems for the largest initial misalignments, and the reconstructions of the full framework suffer from the resulting badly aligned sequences. On the contrary, the greater robustness of the proposed registration method to local minima offered a better support to the shape reconstruction of the full framework, and provides a better average Jaccard measure.

Qualitative Evaluation on Real Data

The brain and hip datasets have only very small initial misalignments of less than 5 pixels and a few degrees. The registration and 3D reconstruction results were inspected visually, and we subjectively assessed that both methods perform similarly well on these datasets.

The range of initial misalignments of the temporal sequences of the cardiac MRIs is more important, with shifts of 1 to 10 pixels. We found that registration by full translations and rotations produced a similarly good quality of alignment on these data than registration by translation in the SA plane only. This is supposedly due to the sequences being acquired at end-expiration, which is a position where the heart undergoes negligible or no rotation and shift along its long-axis. Therefore, we show in Fig. 10 results obtained using translation in the SA plane only.

The datasets in Fig. 10 suffered from severe initial misalignments of 5 to 10 pixels. The SR+ISI method performs less satisfactorily in such cases and yields even poorer results than with the artificial heart datasets. We explain this by the higher complexity of the data whose background contained several objects which can confuse the SR [START_REF] Lotjonen | Correction of movement artifacts from 4-D cardiac short-and long-axis MR data[END_REF] registration method. On the contrary, the proposed registration method registers the datasets very well in all cases and the complete framework recovers plausible shapes of the LV cavity of the heart.

DISCUSSION

We now present a cohesive review of the strengths and limitations of our registration method and its proposed implementation. The main advantage of this method against the co-segmentation via registration based regularization method of Yezzi et al. [START_REF] Yezzi | A variational framework for integrating segmentation and registration through active contours[END_REF] is that it can handle sparse volumes made up of any number of images located in arbitrary planes which contain different parts of the object being segmented. The only requirement on image location is that the images have to intersect each other, although one image does not need to intersect all other images.

A challenge of this more general application field is that the implicit estimation of remaining misalignment, that holds in Yezzi et al.'s special case of images containing the same parts of the object, cannot be used in the general case. One can intuitively see that it makes sense to register images by amounts corresponding to their respective misalignment, so that images that are strongly misaligned get strongly moved towards a better alignment, while ones that are well aligned remain undisturbed. Our experiments showed that our general formulation can compute the direction of alignment as well as Yezzi et al.'s method. However, the remaining misalignment, and therefore amount of registration to be applied to the images, cannot be computed reliably in cases where the images contain different parts of the object and have very few intersection points. Therefore, we proposed a scheme in Section 3.3 to (a) estimate this remaining misalignment (Eqs. ( 14) and ( 15)) and (b) account for it in the amount of registration applied to the individual images (Eqs. ( 16), ( 17), ( 18) and ( 19)).

There are two main limitations to the scheme we propose in ( 16), ( 17), ( 18) and ( 19) for setting the amount of registration. Firstly, the amount of registration is approximated by two separate treatments of translation and rotation. In some of our tests, we even performed registration by translation only, sometimes adding rotation after a number of iterations. This approach is quite appropriate to our case of cardiac MRIs acquired during breath-holds, since the initial amounts of misalignment by translation and by rotation are quite unbalanced, with misalignment by translation being predominant. Our method might benefit from a joint treatment of translation and rotation for a simpler and more general application. Furthermore, although we did not encounter gradient descent issues in our tests, a coupled treatment of translation and rotation, such as by quasi-Newton approaches [START_REF] Nocedal | Numerical Optimization[END_REF], has been suggested as a way to ensure that no such issue arise in other cases.

Secondly, the proposed update scheme in ( 16) and ( 17) is very simple and potentially non-optimum, and its parameters were determined empirically. We simply set the amount of registration to be proportional to the estimated misalignment in the general regime, with a saturation at 1 pixel or 1 degree per iteration to avoid risks of divergence of the gradient descent scheme. We also use either a slower linear scheme or an unaltered amount of registration (whichever is smaller) close to perfect alignment in order to allow for finer registration. The rational for potentially using the unaltered amount of registration in that case is that when images are becoming well aligned, towards the end of the full registration and segmentation process, the 3D segmenting contour starts matching the global shape of the object reasonably well, making the implicit estimation of misalignment valid again. The threshold values on the estimated misalignment were determined empirically. It was found that the method is not very sensitive to these values, and they have not been optimised. This piecewise linear scheme may be replaced by a more elaborate and/or parameter-free one with little impact on the principles of our proposed registration method.

Finally, the optimisation scheme used here is purposefully simple in order to keep the focus on our proposed registration method. We use a simple gradient descent scheme, together with the aforementioned independent processing of translations and rotations, and piece-wise linear scheme for setting the amount of registration. A more elaborate numerical optimisation method could be used to implement our registration strategy and may provide better performance, notably in term of convergence speed.

CONCLUSION AND FUTURE WORK

We have presented an integrated registration, segmentation, and shape interpolation framework to model objects from 3D and 4D volumes made up of spaced and misaligned slices having arbitrary relative positions. We first proposed a new methodology for rigid registration of images that contain different parts of an object. This novel method combines a global approach for better handling of local minima, and a local approach for a better final accuracy. Second, a new level set framework was also proposed, that integrates this registration method and our previous segmentation and shape interpolation method for handling spaced data. It can exploit segmentation results for better robustness of both registration and shape interpolation to large gaps in the data and to limited number of intersections between the images. The framework was validated on artificial data and tested on real MRI and CT scans. In particular, registration is particularly well suited for handling both local minima and fewer image intersections, and it obtained more robust and accurate results than SR [START_REF] Lotjonen | Correction of movement artifacts from 4-D cardiac short-and long-axis MR data[END_REF] in these cases. It yielded similar results to SR [START_REF] Lotjonen | Correction of movement artifacts from 4-D cardiac short-and long-axis MR data[END_REF] in easier scenarios. The complete framework performed significantly better than the sequential approach of registration followed by segmentation and shape interpolation.

As in [START_REF] Yezzi | A variational framework for integrating segmentation and registration through active contours[END_REF], our local variant of registration could be used for non-rigid registration if derived from the non-rigid formulation in [START_REF] Yezzi | A variational framework for integrating segmentation and registration through active contours[END_REF]. However, as in [START_REF] Yezzi | A variational framework for integrating segmentation and registration through active contours[END_REF] also 6 , a regularization method should be chosen with care to avoid arbitrary deformations, and overall our registration method is better suited for rigid registration.

Both registration and shape interpolation rely heavily on segmentation results and this is a limitation of our method. In particular, if a background object is wrongly segmented in an image, then this image cannot be registered properly. Therefore it is important to use reliable segmentation algorithms. The proposed framework can make use of any segmentation scheme, and the choice should be made by the user based on the characteristics of the data. Future work includes investigating robust segmentation algorithms for specific modeling problems.

Fig. 1 .

 1 Fig. 1. Example of misaligned dataset -one short-axis (horizontal) and two long-axis (vertical) intersecting slices in a cardiac MRI. (a) global view of the three slices, and (b) zoom on the heart region. Only one timeframe is shown. Color denotes pixel intensity. Note the misalignment between the horizontal slice and the other two vertical slices, which is due to a shift of the heart between their respective acquisitions.

Fig. 2 .Fig. 3 .

 23 Fig. 2. Illustration of the three types of sequences that may make up a dataset. Temporal sequence: (a) three time-frames of a slice in a cine cardiac MRI. Spatial sequences: (b) central slices of two sequences acquired with different scanner parameters in an MRI of the pelvis, (c) central slices of two sequences acquired with different slice orientations in a brain MRI -top: axial, bottom: sagittal.

Fig. 4 .

 4 Fig. 4. Speeds used for registration in the local and global variants of the registration method. (a) image and segmenting contour in yellow. (b) and (c) segmentation speeds f • g, computed using the piecewise constant model from [28], on the contour C in (b) (local variant) and on C ∪ Ω diff in (c) (global variant) -red: negative; green: positive. Note that the global variant (c) exploits the extra information about the geometry of the object provided by f • g in Ω diff .

Fig. 5 .

 5 Fig.5. Registration of artificial brain datasets that have local minima. SR[START_REF] Lotjonen | Correction of movement artifacts from 4-D cardiac short-and long-axis MR data[END_REF] is implemented with NMI, and implementations with SSD and NCC suffered from the same local minima issues. Colors correspond to different sequences, and may not be all visible in case of near perfect alignment. For visualization purposes, only pixels containing the object are displayed. Due to space constraints, the groundtruth registration is not shown as it has no visible differences with the results of the proposed method.

Fig. 6 .

 6 Fig. 6. Two possible difficulties in slice-wise registration: (a) ambiguity in the position of a slice due to insufficient constraints, and (b) contour folding around individual slices due to a too sharp RBF.

Fig. 7 .

 7 Fig. 7. Shape reconstructed from the two datasets of Fig. 5. Top: dataset 1, view from the top. Bottom: dataset 2, view from a front-right position.

Fig. 8 .

 8 Fig. 8. Shape reconstruction from a dataset containing slices with ambiguous poses. Left: slices after registration. Note that only one LA (vertical) slice is available, hence the numerous ambiguities in the pose of the SA slices. Right: reconstructed 3D shape.

Fig. 9 .

 9 Fig. 9. Registration and shape reconstruction from a kidney CT-scan made up of two spatial sequences with large (15 pixels) gaps, severe initial misalignment, and local minima due to cluttered background. Top: registration result -for visualization purposes, only one central slice of each spatial sequence is displayed. (a) SR [4] and (b) proposed registration methods. Bottom: reconstructed shape -(c) original shape reconstructed from the volume with groundtruth alignment, (d) SR+ISI, (e) proposed framework.

Fig. 10 .

 10 Fig. 10. Registration, segmentation, and shape interpolation of the LV cavity of the heart from two real 4D cardiac MRIs. Bumps in (a) and (c) are due to misaligned slices. The bump at the bottom of shape (b) is not due to a misaligned slice but is a leak of the segmentation surface which could be avoided by choosing a more robust segmentation algorithm.

  Algorithm 1 -Integrated Registration, Segmentation, and Shape Interpolation for Spaced and Misaligned Data 1: Initialize a contour C 2: Initialize all g n to the identity transformation 3: Initialize a boolean global variant to true if local minima are expected, false otherwise 4: repeat

	5: 6:	Compute dα dt using (11) and (23) if global variant = true then
	7:	Update C using (26)
	8:	else
	9:	Update C using (24)
	10:	end if
	11:

4. See www.cs.bris.ac.uk/ ∼ paiement/IReSISD.html.

TABLE 1

 1 Datasets used in the experiments

	Dataset	Orientation	Nb. slices	Pixel size (mm)	Slice spacing	Time-frames	Nb. datasets	ISI scheme
	Artificial cardiac MRI	SA LA	4 to 19 1 to 12	1.77 to 2.08	5 mm to 20 mm -	1	14	2 steps: γ = 0.7 then 3
	Artificial brain MRI	Axial (2 sets) Sagittal (1 or 2 sets)	11 to 18 13 to 20	0.71 to 0.94 0.37 to 0.90	3.3 mm to 6.6 mm 3.3 mm to 4.8 mm	1	9	Simple: γ = 2
	Kidney CT	Axial Sagittal	29, 15, 10 and 8 103, 52, 35 and 26	0.78	5, 10, 15 and 20 pixels	1	4	Simple: γ = 3, 4, 5 and 6
	Hip bone CT	Axial Sagittal	27, 14, 9 and 7 103, 52, 35 and 26	0.49	5, 10, 15 and 20 pixels	1	4	Simple: γ = 3, 5, 7 and 9
	SA+LA cardiac MRI	SA LA	4 to 19 1 to 12	1.77 to 2.08	5 mm to 20 mm -	25	14	2 steps: γ = 1 then 4
	Radial cardiac MRI	LA SA	12 1	1.77	∼ 15 degrees	25	1	2 steps: γ = 1 then 4
	Brain MRI	Axial (2 sets) Sagittal (1 or 2 sets)	11 to 18 13 to 20	0.71 to 0.94 0.37 to 0.90	3.3 mm to 6.6 mm 3.3 mm to 4.8 mm	1	9	Simple: γ = 1.5
		Axial	19	0.74	6 mm
	Hip MRI	Coronal (2 sets)	19	1.64	4.8 mm	1	1	Simple: γ = 1.5
		Diagonal (2 sets)	15	0.86	3.3 mm

TABLE 2

 2 Stack-wise registration accuracy on artificial data: mean and std of registration error (in pixels and degrees) averaged over all 9 datasets (0.21 • , 0.23 • , 1.69 • ) (0.40 • , 0.24 • , 0.10 • ) (0.05 • , 0.08 • , 0.05 • ) (0.27 • , 0.35 • , 0.22 • ) ±(6.7e-4, 6.0e-4, 2.5e-4)

	Initial misalignment	SR [4] with SSD	SR [4] with NCC	SR [4] with NMI	Proposed

T n i ini = rand [-10, 10] θ n i ini = rand [-10 • , 10 • ]

(0.50, 1.07, 0.36) (0.06, 0.13, 0.13) (0.02, 0.05, 0.08) (0.05, 0.08, 0.09) ±(1.8e-3, 1.8e-2,6.7e-2 ) ±(4.1e-3, 3.3e-2, 7.9e-3) ±(1.2e-4, 4.3e-3, 2.3e-3) ±(3e-3, 6.7e-3, 2.5e-3)

  , 4.76 • , 5.78 • ) (3.34 • , 2.27 • , 1.14 • ) (2.15 • , 4.70 • , 2.14 • ) (0.34 • , 0.40 • , 0.22 • ) ±(0.88, 0.22, 0.54)

	(8.09, 8.68, 4.37)	(0.63, 2.26, 1.2)	(0.74, 2.61, 1.99)	(0.05, 0.08, 0.07)
	±(132, 57.5, 15)	±(0.3, 5.39, 1.04)	±(0.33, 6.81, 2.11)	±(2e-3, 4.6e-3, 3.4e-3)
	(7.76

• , 20 • ] • ±(0.13, 0.08, 0.02) ±(0.05, 0.23, 19.47) ±(2.7e-3, 3.9e-3, 6.9e-4)

  Table 2 shows registration errors for all degrees of freedom in pixels and degrees, averaged on all 9 datasets.

TABLE 3

 3 Slice-wise registration accuracy on artificial data: mean and std of registration error (in pixels and degrees) averaged over all 14 datasets

	Initial misalignment	SR [4] with SSD	SR [4] with NCC	SR [4] with NMI	Proposed
	T n xy ini = rand [-5, 5] T n xy ini = rand [-10, 10]	(0.72, 0.98) ±(0.71, 2.59) (1.92, 2.29) ±(8.67, 7.75)	(0.62, 0.64) ±(0.65, 0.72) (1.27, 1.71) ±(4.73, 7.79)	(0.64, 0.77) ±(0.65, 0.87) (1.87, 1.99) ±(3.72, 6.11)	(0.36, 0.39) ±(0.42, 0.2) (0.36, 0.36) ±(0.44, 0.18)

  68, 1.72, 1.72) (1.4, 1.33, 2.09) (1.29, 1.51, 2.02) (1.05, 0.78, 0.88) ±(3.8, 4.12, 1.71) ±(1.37, 1.82, 3.12) ±(1.36, 2.16, 2.54) ±(1.12, 0.68, 0.56) (3.48 • , 3.33 • , 3.71 • ) (3.31 • , 3.44 • , 4.05 • ) (3.23 • , 2.96 • , 3.2 • ) (3.1 • , 2.97 • , 2.65 • ) ±(0.11, 0.1, 0.13) , 3.57 • , 4.11 • ) (3.95 • , 3.97 • , 4.44 • ) (3.04 • , 2.95 • , 2.95 • ) (3.12 • , 2.52 • , 2.79 • ) ±(0.14, 0.11, 0.22)

	(3.19, 2.91, 2.82)	(2.56, 1.89, 2.56)	(2.91, 2.72, 2.08)	(1.2, 1.15, 1.13)
	±(11.9, 11, 7.12)	±(5.05, 3.48, 3.88)	±(4.86, 4.45, 1.95)	±(3.65, 1.86, 0.78)
	(3.87

±(0.11, 0.11, 0.16) ±(0.09, 0.1, 0.09) ±(0.09, 0.09, 0.07)

T n xy ini = rand [-10, 10] T n z ini = rand [-4, 4] θ n i ini = rand [-5 • , 5 • ] •

TABLE 5

 5 Stack-wise registration accuracy on real data: mean error (in pixels and degrees) averaged over all datasets for each test and slice spacing

	Initial misalignment	Slice spacing	SR [4]	Proposed
		5 pixels	(0, 0, 0.11) (0.01 • , 0.001 • , 0.001 • )	(0.01, 0.01, 0.12) (0.02 • , 0.02 • , 0.02 • )
	T n i ini = rand [-10, 10], θ n i ini = rand [-10 • , 10 • ]	10 pixels 15 pixels 20 pixels	(0, 0, 0.11) (0.002 • , 0.001 • , 0.001 • ) (0.001, 0.001, 0.11) (0.01 • , 0.003 • , 0.001 • ) (0.11, 0.56, 0.42) (3.33 • , 0.47 • , 1.95 • )	(0.01, 0.06, 0.13) (0.06 • , 0.09 • , 0.1 • ) (0.09, 0.05, 0.08) (0.11 • , 0.16 • , 0.1 • ) (0.08, 0.28, 0.25) (1.39 • , 0.26 • , 1.09 • )
		Average	(0.03, 0.14, 0.19) (0.84 • , 0.12 • , 0.49 • )	(0.05, 0.1, 0.14) (0.4 • , 0.13 • , 0.33 • )

TABLE 7

 7 Slice-wise registration accuracy on real data: mean and std of registration error (in pixels) averaged over all datasets for each test

	Slice
	spacing	SR [4]	Proposed
	& shift max
	5 pixels	(0.34, 0.54)

±(0.22, 0.69) (0.12, 0.2) ±(0.01, 0.03) 10 pixels (0.34, 0.92) ±(0.14, 3.58) (0.17, 0.35) ±(0.02, 0.06) 15 pixels (0.27, 0.27) ±(0.02, 0.02) (0.26, 0.77) ±(0.13, 0.71) 20 pixels (4.22, 3.45) ±(27.84, 27.25) (0.48, 0.84) ±(0.19, 0.56) Average (1.29, 1.29) ±(7.05, 7.89) (0.26, 0.54) ±(0.09, 0.34)

TABLE 8

 8 Jaccard coefficients yielded on real data

	Test	SR+ISI	Proposed
	Stack-wise 1	0.94±0.05	0.94±0.04
	Stack-wise 2	0.72±0.35	0.93±0.05
	Slice-wise	0.88±0.13	0.93±0.04

ACKNOWLEDGMENTS

The CT-scans datasets used in the experiments were from the Laboratory of Human Anatomy and Embryology, University of Brussels (ULB), Belgium.

6. see Section 3.4 of [21]