Z Jiang

J.-F Witz

P Lecomte-Grosbras

J Dequidt

C Duriez

M Cosson

S Cotin

M Brieu

B-spline Based Multi-organ Detection in Magnetic Resonance Imaging

Keywords: B-spline, MR images, pelvic organs detection, virtual image correlation

In the context of the female pelvic medicine, non-invasive Magnetic Resonance Imaging (MRI) is widely used for the diagnosis of pelvic floor disorders. Nowadays in the clinical routine, diagnoses rely largely on human interpretation of medical images, on the experience of physicians, with sometimes subjective interpretations. Hence, image correlation methods would be an alternative way to assist physicians to provide more objective analyses with standard procedures and parametrization for patient-specific cases. Moreover, the main symptoms of pelvic system pathologies are abnormal mobilities. The FEM (Finite Element Model) simulation is a powerful tool for understanding such mobilities. Both the patient-specific simulation and the image analysis require accurate and smooth geometries of the pelvic organs. This paper introduces a new method that can be classified as a model-to-image correlation approach. The method performs fast semi-automatic detection of the bladder, vagina and rectum from MR images for geometries reconstruction and further study of the mobilities. The approach consists of fitting a B-spline model to the organ shapes in real images via a generated virtual image. We provided e cient, adaptive and consistent segmentation on a dataset of 19 patient images (healthy and pathological).

Introduction

Female pelvic disorder such as the organ prolapse is a common problem that occurs especially as a woman ages: 20% to 30% of women of all ages combined su↵er from severe degree of prolapse [START_REF] Samuelsson | Signs of genital prolapse in a Swedish population of women 20 to 59 years of age and possible related factors[END_REF] and more than 60% of women over 60 years of age are a↵ected by this pathology [START_REF] Swift | The distribution of pelvic organ support in a population of female subjects seen for routine gynecologic health care[END_REF]. These problems are related to the mobility of female pelvic system where a physician's diagnosis is usually based on his meaningful analysis of medical images (e.g. dynamic MRI of a patient breathing or pushing). The subjectivity related to human perception or medical experience can not be avoided and can cause variability in the diagnoses. The advantages of involving numerical models relies not only in visualization, but also in quantitative analysis on patient-specific medical images, which can make the diagnoses more objective. Besides that, recent works [START_REF] Venugopala Rao | Experiments and finite element modelling for the study of prolapse in the pelvic floor system[END_REF][START_REF] Vallet | Simulation of pelvic mobility: topology optimization of ligamentous system[END_REF][START_REF] Lecomte-Grosbras | Towards a better understanding of pelvic system disorders using numerical simulation[END_REF] have attested the feasibility and have provided promising results of finite elements model simulations for evaluating the pelvic organs mobility. In order to perform patient-specific FEM simulation or image analysis, the reconstruction of an accurate and smooth geometry of the pelvic organs is required, and it should be e cient. In this paper, we propose a semi-automatic method for organ shapes detection, which is an important preliminary step for further studies and modelling. Not yet can this method provide direct diagnoses, but it is an e cient assistance with common procedure for analysing medical images of patients, especially the pathological cases. Existing work [START_REF] Bay | Geometric modeling of pelvic organs[END_REF] has introduced an algorithm creating thick surfaces of hollow pelvic organs, using periodic B-splines and o↵sets. This modelling was a step between segmentation and physical modelling. The reconstruction of geometry was based on a cloud of points resulted from presegmentation. This modelling was a step between manual segmentation, needed by the method, and physical modelling. More generally, the surface modelling of patient-specific organ geometries has become an increasingly important issue for computer-assisted simulation. This objective is usually achieved in two steps: segmentation and geometric modelling. Regarding the first step, a multitude of algorithms have been proposed to convert original image data to segmented ones (labelled voxels). From Level Set Methods [START_REF] Malladi | Shape modeling with front propagation: A level set approach[END_REF][START_REF] Sethian | Level Set Methods and Fast Marching Methods: Evolving Interfaces in Computational Geometry[END_REF] to Graph Cuts [START_REF] Paragios | Handbook of Mathematical Models in Computer Vision[END_REF], algorithms have achieved promising results for segmentation and grouping. Then for the latter, Marching Cubes [START_REF] Lorensen | Marching Cubes: A High Resolution 3D Surface Construction Algorithm[END_REF] algorithm has become the standard approach to generate surface meshes from a scalar field. Despite of their success, these methods are all pixel-based approach, which do not provide directly smooth and accurate contour. However these aspects are required for further processing such as extraction of displacement along contour and 3D model reconstruction. Another modelling method, "model-to-image correlation",is carried out in two steps:

1. Initialization of a fine model in terms of geometry and topology 2. Deformation of the model for fitting the original image by means of cost function minimization Such idea as matching deformable model to images was introduced by Kass et al. [START_REF] Kass | Snakes: Active contour models[END_REF]. In this approach an "Active contour (snakes)" model is driven by internal and external "Forces" to fit the patterns in images. This idea was improved and developed by a host of techniques [START_REF] Mcinerney | Deformable models in medical image analysis: a survey[END_REF][START_REF] Davatzikos | Spatial transformation and registration of brain images using elastically deformable models[END_REF]. In [START_REF] Blake | Active Contours[END_REF], a deformable model was described as a variation of degrees of freedom and coupled with probabilistic modelling. Segmentation results were hence based on some prior knowledge. Brigger et al. [START_REF] Brigger | B-Spline Snakes: A Flexible Tool for Parametric Contour Detection[END_REF] introduced powerful B-Spline functions for building the snake model, which required less degrees of freedom and produced higher continuity. Li et al. [START_REF] Li | Active Contour External Force Using Vector Field Convolution for Image Segmentation[END_REF] proposed an external force based on vector field convolution for the active contours to deal with the image noises and the shape concavities. With the classical active contour, the external force is merely guided by the pixel intensity of the image. As a consequence, the algorithm does not produce desired results if the points situated in a blur background or far from the optimal positions. Generally, one needs to perform additional processing on the original image to promote the right deformation of snakes (as an example applying gradient operator to enhance the edges). Moreover, these algorithms require a good initialization of the models, which is not simple to create. Chan et al. [START_REF] Chan | Active contours without edges[END_REF] proposed a model to deal with this problem. An implicit curve model can detect contours without gradient edgedetector and the initial curve can be placed anywhere in the image. However, we are introducing a new B-spline like method which is more consistent with our numerical approach. Practically in these ill-posed problem, the choice of the number of the degrees of freedom is not trivial. More degrees of freedom enable to create a finer geometry descriptor; however, a finer regularization of these variables is also required. Such model-based approach is widely used for medical image processing [START_REF] Smith | Segmentation and Tracking of Cytoskeletal Filaments Using Open Active Contours[END_REF][START_REF] Muralidhar | Snakules: A Model-Based Active Contour Algorithm for the Annotation of Spicules on Mammography[END_REF][START_REF] Berendsen | Free-form image registration regularized by a statistical shape model: application to organ segmentation in cervical MR[END_REF] and it is well suited in our case, whereas in this paper, we introduce a new model-to-image approach. The idea is to propose a method that combines the method of "Active Contour" with the "Virtual Image" approach, the latter was firstly introduced by [START_REF] Semin | Analytical shape determination of fiberlike objects with Virtual Image Correlation[END_REF][START_REF] Réthoré | Curve and boundaries measurement using B-splines and virtual images[END_REF]. Our approach relies on several contributions which are: (1) A fast and semi-automatic organ geometries modelling is proposed, which is helpful for our further study of the mobility and modelling. (2) A virtual image is generated from the mathematical model, in order to avoid the local optimum problem and the pre-processing. Based on the virtual image, we use a cost function to find the best correlation between the virtual image and the real one. Finally, the model-to-image correlation can be written as an energy minimization problem, consistent with a typical image registration framework [START_REF] Hill | Medical image registration[END_REF][START_REF] Ibanez | The ITK software guide[END_REF][START_REF] Klein | Elastix: a toolbox for intensity-based medical image registration[END_REF]. Concerning the medical imaging, the correlation method is denoted as image registration. The objective is to find the spatial correspondence of two medical images, one of which is the virtual image generated from the model in our case. Similar with the principle in [START_REF] Ibanez | The ITK software guide[END_REF][START_REF] Klein | Elastix: a toolbox for intensity-based medical image registration[END_REF], we fit the model to image spatially by measuring the similarity between the virtual image and the medical image. (3) Besides the aspect of the image preprocessing, our algorithm consists of a more general framework including the active contour [START_REF] Kass | Snakes: Active contour models[END_REF][START_REF] Brigger | B-Spline Snakes: A Flexible Tool for Parametric Contour Detection[END_REF]. For example, the cost function can be replaced by the one of the active contour. The adaptive refinements and decimations of control points can be included thanks to the flexibility of B-spline model, for the purpose of adjusting the number of the degrees of freedom during optimization. (4) A problem in multi-organ detection is to deal with the unwanted inter-penetration between two di↵erent structures which is not favourable for the mechanical simulation. We use a signed distance map to compute an energy term quadratically increasing when organs approached towards each other. A general view of our approach is presented in the MRI and Correlation Method Section. The geometric model is described in the B-spline Model Section. The generation of virtual image is detailed in the Virtual Image Section. In the Optimization Section, we present how to couple the two aspects in a common energy minimization framework. The Results Section reports our experimental results on a data set of 19 patient images. Finally, we give a conclusion and perspectives of our work in the last section.

Medical Image

Model

MRI and Correlation Method

In the point of view of computation, the proposed approach can be formulated as an optimization procedure. Four key parts are involved in such optimization (see Figure 1):

input data (3D static and 2D dynamic MR images of the patient) -a mathematical model with variables to be optimized (B-spline Model Section) -a cost function which links the model to the input data (ex., integration on the domain of the model) (Cost Function Formulation Section)

an optimizer which finds the optimal values of the parameters to minimize the cost function (Optimization Section) The input experimental image data consist of T2-weighted static and dynamic MRI of the female pelvic system. The use of these images of volunteers for our research was approved by the National Ethics Committee (Comité d'Ethique de la Recherche en Obstétrique et Gynécologie CEROG 2012-GYN-06-01-R1). The static MR images are obtained in the three planes (sagittal, axial and coronal), which provide information on the 3D anatomical organization. The dynamic MR images (or Cine MRI) refer to a temporal sequence of 2D images, which are widely used for diagnosing pelvic sagging pathologies [START_REF] Pannu | Dynamic MR Imaging of Pelvic Organ Prolapse: Spectrum of Abnormalities 1[END_REF]. For our study, these 2D dynamic images were obtained in the same midline sagittal plane of the patient during pushing, for a given frequency (about 1.2 fps). The choice of this middle plane keeps a standard procedure in the clinical routine because it's the most representative of the organ shapes and motions. In order to increase the image contrast of pelvic organs, a gel had been injected into the vaginal cavity and the rectum while patient drank liquid to highlight the bladder. In this way, pelvic organs were easily observable during MRI acquisition. We were interested in 2D segmentation of the three organs (bladder, vagina and rectum), for which both 2D sagittal static and dynamic MR images were used (spatial resolution from 0.47 to 1.17 mm/pixel, see Figure 2). Our registration procedure is a multi-scale optimization in two steps: an a ne transformation of the model for the coarse registration and a B-spline deformation for the finer registration. For the first step, our model is defined by the 6 degrees of freedom (DOF) of the a ne transformation applied on the initial contour (a circle-like curve whose center and radius are predefined by the user). For the second step, as we use a B-spline model for the description of geometry, the control points are considered as DOF of the mathematical model. The optimization procedure updates the positions of control points to find the minimum of the cost function. This part ensures the flexibility and quality of our model. Moreover, we use an adaptive B-spline model to adjust the number of control points during the optimization. This is hoped to be a powerful improvement: integrating the insertion and removal of control points with respect to the maximum error along the contour. Due to the feature of multi-level framework, the registration is more stable to reach the desired minimum and with a faster convergence. However, a major problem occurs when the contours of organs intersect each other, which leads to divergence of the algorithm. This will also create unwanted artifacts for the simulation. We avoid the problem by adding a collision energy term that prevents the inter-penetrations.

B-spline Model

The optimization process can be described as iterative updates of its degrees of freedom (DOF). As we deform the geometric model, we optimize the DOFs that describe the geometrical shape. The geometric model is considered as a mapping from DOFs to the geometry. In this paper, we introduce two mappings: (I) an a ne transformation for the coarse registration and (II) a cubic B-spline as the shape descriptor for the deformable registration. The control points are the researched DOFs (see [START_REF] Piegl | The NURBS Book (2nd Ed.). Monographs in visual communication[END_REF] for definitions). Thus similar with the typical image registration, our model-to-image registration method is also based on a multi-scale framework. Concretely the geometries of the relevant organs are represented in an analytical manner. In order to generate a virtual image in the neighbourhood of the contour, the DOF vector V is mapped into a set of discrete points in the neighbourhood of the contour, which presents a narrow band containing the contour, denoted by ⌦ in this paper. These points form a bended grid over ⌦ as the domain of integration (details presented in the Virtual Image Section). A position on the curve is defined using the curvilinear coordinate u based on the parametric curve. For modelling the geometries of the three organs, we used 2D B-spline of 3-degree that is a classical geometric modelling tool. In our case, each organ is modelled by a closed B-spline. A uniformly spaced knot vector U = {u0, . . . , um} is used to define the basis functions of p-degree, denoted as Ni,p : u 2 [0, 1] ! [0, 1]. The Ni,p(u) is the ith p-degree piecewise polynomial which is a linear combination of two (p 1)-degree basis functions, defined as:

Ni,p(u) = u ui ui+p ui Ni,p 1(u) + ui+p+1 u ui+p+1 ui+1 Ni+1,p 1(u). (1)
Hence the basis function is computed recursively as the degree p increases. Each interval [ui, ui+1] defines a span of the B-spline curve. The multiplicity of a knot determinates the continuity of the curve on this point. Geometrically each knot ui in the vector U is mapped into a breakpoint of the curve. Generally the knot vector is setted from 0 to 1 and with the multiplicity p + 1 at the first and the last knot to form a clamped curve. Figure 3 shows an example of the basis functions of degree 0, 1, 2 and 3 defined on a uniform knot vector with 16 knots. Due to the multiplicity of the knots u0, u1, u2 and u3, the first three basis functions of low degree are identically zero. Thus with no loss of generality, N3 and N4 (defined on spans of distinct knots) are used for illustration.

Each organ is presented by a parametric B-spline curve

M : u 2 [0, 1] ! M(u) 2 R 2 .
Each position on the curve can be calculated by the value of u:

M(u, [P0, P1, • • • P l ]) = l X i=0 Ni,p(u)Pi (2)
where [P0, P1, • • • P l ] are the l + 1 control points. However, the first and last control point of each curve are attached to form a closed curve: M(0) = M(1) = P0 = P l . Hence our model is analytical and of C 2 continuity except at the end point. Under this configuration, the l control points are supposed to be the degrees of freedom in the registration algorithm. Another advantage of the B-spline is its locality property which means each control point is highly independent. The change of one's position merely influences a certain spans of the contour (see Figure 4). In addition, for the a ne transformation step, we have the relevant mapping MA, 6 DOFs of which are evaluated in this coarse step: [T11, T12, T21, T22, Tx, Ty]. Suppose M 0 is the mapping function corresponding to the initial B-spline curve, and C the center, then the mapping function of the a ne transformation becomes:

MA(u, [T11, T12, T21, T22, Tx, Ty]) =  T11 T12 T21 T22 (M 0 (u) C) +  Tx Ty + C. (3)
4 Cost Function Formulation

Virtual Image

In order to define an appropriate cost function, we introduce a virtual image to form the first energy term. A so-called virtual image is generated from the model for calculation of a cost function which finds the best correlation with the real image. This virtual image is created in the narrow band ⌦ along the contour with a certain width (see Figure 6(b,c)). The value of gray levels in this narrow band image is chosen to be similar to the ones in the real image (gray level MR image). This idea was first introduced and proved by Semin et al. [START_REF] Semin | Analytical shape determination of fiberlike objects with Virtual Image Correlation[END_REF] to detect fibre-like objects in digital images. This approach was developed then by Réthoré et al. [START_REF] Réthoré | Curve and boundaries measurement using B-splines and virtual images[END_REF]. However both the width value and the gray levels are not necessary to be exactly close to the real image. Let us consider the one-dimension case as an example. To detect the o↵set of a discrete 1-D sinusoidal signal, we can define a "virtual" smooth and analytical sinusoid of which the o↵set is optimized to find the best correlation between these two signals. The optimal value of o↵set can be obtained with a range of values chosen close to the magnitude of the discrete signal, thus the parametrization of the analytical signal is not required to be precisely equal with the discrete signal. The analytical sinusoid presents an approximate intensity profile depicting how the discrete signal evolves in one direction. Similarly, in 2D manner considering a discrete image of an object, we create such a virtual image by associating an intensity profile at each point on the object border with a certain width to simulate the variation of gray level from foreground to background (ex. decreasing gradually from white to black). We modelled the object with a constant intensity profile which implies an implicit assumption that the variation of gray levels in the direction perpendicular to the object border propagates equally in the direction of the tangent. There may be concerns for certain cases where image data do not provide the homogeneous property for an organ; however, in practice, a single organ has a similar intensity in our images. The intensity profile of the virtual image is defined by a levelset function L : v 2 R ! R where v is the distance from a position X to the contour (see Equation 4). In [START_REF] Réthoré | Curve and boundaries measurement using B-splines and virtual images[END_REF], the choice of the levelset intensity was discussed, a bell-shaped function was used to approach the intensity profile in the real image. The choice of sinusoidal function is kept in this paper for its simplicity and the consistency is proved (see Figure 5). In the application, the band width and the o↵set of the levelset function are to be adjusted for a better performance.

L(v) = 8 > > > > > < > > > > > : Amax if v < 0, Amin + Amax Amin 2 ⇤ (1 + cos( ⇡v ⌧ )) if 0  v  ⌧, Amin if v > ⌧. (4)
In Equation ( 4), ⌧ is the estimated width of narrow band ⌦ and Amax, Amin the gray levels inside and outside organ. The value of Amax and Amin are chosen by clicking on the foreground and background in the real image. The virtual image Iv is then computed on the discrete points cloud over the narrow band of the organ border. The discrete points form actually a bended grid, a position X(u, v) of which is defined in the curvilinear coordinate system based on tangent and normal direction of the contour. The virtual image value Iv(X) is computed in the analytical way (see Equation 5) so that the grid spacing size can even be much smaller than the pixel spacing in real image. Unlike the virtual image, the real one has a pixel-wise square grid, we should then project X and evaluate real image value Ir(X) using a bi-cubic interpolation (see Figure 6(a)). Also for these reasons, the virtual image is considered "virtual" and provides subpixel information compared to the real one. Furthermore, gray level values are merely evaluated in its definition domain ⌦ instead of the whole domain of the real image (see Figure 6(b)).

Iv(X) = Iv(u, v) = L(v) ( 5 )
The positions X can be evaluated by the analytical formulation of the parametric geometric model (see Equation 6). Suppose a parametric curve :

M : u 2 R ! R 2 . X(u, v) = M(u) + v ! N(u) ( 6 )
Thus, the computation of the virtual image is more natural and direct with a parametric model M but can be applied to any type of model. The cost function related to this virtual image for each organ can be computed with:

Eimage([P0, P1, • • • P l ]) = Z X2⌦(Band) [Ir(X) Iv(X)] 2 d⌦. (7)
Thanks to the virtual image, our method is exempt from image preprocessing and based merely on reliable original data (no need of preprocessing the images). The method requires less computation and in a mathematical point of view, the virtual image can be considered as a smoothing of gradient field which avoid local optimum in this non-convex problem.

Collision Term

To avoid the inter-penetration of the analytical functions defining the contours of organs during deformable registration, we introduce a collision term to penalize the movement of DOFs (control points). The second energy term is computed by using signed distance maps D : X 2 R 2 ! D(X) 2 R. Each contour retains its own distance map describing how far a position in the image domain is from the contour: D bl for the bladder, Dva for the vagina and Dre for the rectum. The collision energy term is defined by a quadratic function with regard to the distance D,

F : D 2 R ! R + : F (D) = 8 < : Amax if D < , Amax 4 2 (D ) 2 if  D  , 0 if D > . ( 8
)
where defines the width of the possible intersection zone associated with the contour, so that it penalises the possible intersections when one contour happen to approach the narrow band of another closely. As shown in Figure 2, the vagina is at the middle of the three organs, so the bladder and the rectum are never in contact. We use the following equations to compute the collision term for each contour: As the virtual image, F is used as a levelset function for computing the energy associated distance map. The map for each organ is redefined once at each iteration. Thus for example, the value of F (Dva(X)) for a given position X is obtained directly by projecting X onto the map of the vagina. We do not recompute the distance Dva(X) for each discrete point at the same iteration. Finally these two energy terms are summed to form a global energy term as the cost function.

E Bl = Z X2⌦ Bl (Band) F (Dva(X)) d⌦ Bl ( 9
)
EV a = Z X2⌦ V a (Band) [F (D bl (X)) + F (Dre(X))] d⌦V a (
E global = Eimage + ↵E collision ( 12
)
Eimage is to be associated with each organ. E collision is to be replaced by E Bl , EV a and ERe respectively. The coe cient ↵ is chosen empirically during tests, the same value is used for the images of the same spatial resolution.

Optimization

In our registration procedure, optimization is the step connecting the cost function (energy) and the model (Figure 1). It updates iteratively the DOFs of the model. For such a non-linear optimisation problem, the algorithm is based on gradient. We opt for a gradient descent algorithm in this paper, which is suitable for optimizing the position and the number of the control points, and for controlling the convergence with regards to the adaptive refinement step. However, the Newton-Raphson method can also be computed similarly, in Figure 1, the Optimization component can be flexibly replaced by another optimizer. The searched direction of the DOFs depends on the derivatives of the cost function. Let Ti be one of the DOFs in Equation 3, the derivative is computed as follows:

@Eimage @Ti = Z X2⌦(Band) 2[Ir(X) Iv(X)] • @Iv @X |{z} A • @X @Ti |{z} B d⌦ (13)
As mentioned in the B-spline Model Section, our model is optimized via a multi-scale registration procedure. At the initialization step, user defines a point at the inside of each organ as the center of B-spline curve. We have tested a new initialization, while asking the user to put the point approximatively at the zone where the portion of organ is larger. We have observed a better increase of the initial contours at the first steps. Concerning the number of the control points, thanks to the adaptive refinement, that we consider as one of our major contributions, we think it may not have a real importance because the points are to be added and deleted automatically after. And we give 10 by default. Then a circle-like B-spline with 10 control points is created for each organ, which is a basic and simple one-click initialization (see Figure 7(a)).

For the coarse registration, the parameters of a ne transformation (Equation 3) are optimized (search direction computed by Equation 13). In Figure 7(b), one can observe that after the a ne transformation, the three circle-like shapes are transformed to ellipse-like ones which give an approximate initialization of the detection procedure. Following the first step, a finer registration is applied using deformations based on B-spline model. Thus the positions of control points are DOFs to be updated to fit the organ contours by means of minimizing a cost function:

@Eimage @Pi = Z X2⌦(Band) 2[Ir(X) Iv(X)] • @Iv @X |{z} A • @X @Pi | {z } B d⌦ (14)
The partial derivative A is the gradient of virtual image which can be computed by the derivative of the levelset function (Equation 4). The partial derivative B is the mapping function (Equation 2) which means B-spline interpolation function in this case. Thanks to the power of analytical formulation, the required derivatives are simple to compute. In addition, the chain rule of derivation shows the independence of geometric model and the image energy computed by the di↵erence between the virtual image and the real one. Hence one can use the virtual image approach with di↵erent geometric representations (or mapping function). The collision term E collision can be computed in a similar way. Due to the nature of gradient descent, after certain iterations, it shows a phenomenon of local minimum (the energy does not evolve). We use this phenomenon as the criterion for adding control points. Figure 8 illustrates the cost function with respect to the iterations. A test of local minimum is performed automatically for every 10 iterations. Let Emax be the maximum value of the cost function during this 10 iterations, Emin the minimum value and E0 the initial value:

kEmax Emink E0 < ✏. (15)
If Equation 15 is satisfied, the optimization is supposed to reach a local minimum (low variation of the cost function), then the curve is refined where the cost function needs to be reduced. The refinement step is treated automatically by inserting and removing control points based on the evaluation of cost function along the contour (details can be found in the next two paragraphs).

Once the curve is refined, the optimization passes the local stability and continues to find the best correlation as shown in Figure 8. This test and the adaptive refinement are performed repeatedly until the refinement can no more a↵ect the cost function. In order to show this point, in Figure 8, we set the optimization to 800 iterations and it can be observed that after certain iteration, the insertion and removal (simplification) of the points become periodic, however the cost function is finally invariant. Practically for the stop condition test, we can elevate the number of iterations from 10 to 40 by using the criteria of Equation 15. At each iteration, the variance of the cost function during the last 40 iterations is tested. As any further iteration or refinement does not provide any change, the algorithm reaches its convergence. Finally in Figure 8, the optimal positions at iteration 404 are found. Moreover, by integrating insertion and removal of control points, our method can refine adaptively the geometric model with respect to images. The major advantage is that our B-splines fit properly the forms even in case of local aberrant curvature changes and avoid useless control points for description of the curve. This feature makes our method capable of a balance between model complexity and simplicity of optimization.

The insert of control points does not change the B-spline curve by using its geometric algorithms. Every time when the gradient descent algorithm reaches its local minimum detected by a threshold of variation of the cost function, we refine the curve by subdividing the knot interval where appears the maximum value of cost function. On the other hand, when two control points approach each other within a distance threshold, one is deleted for eliminating inter-cross of curve. In Figure 7(d,e,f), we illustrate the distribution of cost function along the contour at the three stages (initialization, after a ne transformation and the final position). The values of the cost function decrease significantly during the optimization. Moreover, the distance from a position on the contour to the organ border can be indicated by evaluating the metric map along contour. By zooming inside the contour, higher values occur when it is far from the real image border, which can be the indication for refinement.

Results

We validate our method of segmentation on a database of 19 patient images. The types of these patients include the witness, the prolapse and the endometriosis (a pathology of the pelvic system related to the hypo-mobility). Each patient data set consisted of a T2-weighted MRI of the pelvic system We used for this 2D segmentation test the image in the midline sagittal plane of each patient which is representative and useful for further study of the pelvic mobility. By the end of each optimization process, a Medical Doctor (MD) manually corrected the shape of organs by modifying the positions of control points so that the final curves and the contours of organs, to the MD opinion, were well fitted. These final curves were used as the reference in our validation.

In Figure 10, we illustrated the detection and the validation for each type. Then numerically we compared the curves before manual corrections with the final ones by evaluating error with several metrics: Dice Coe cient (DC), Hausdor↵ Distance (HD) (see [START_REF] Kolawoleo | Comparison and Evaluation of Segmentation Techniques for Subcortical Structures in Brain MRI[END_REF] for the definition) and Average Curve Distance (ACD). Figure 9 shows the computation of these metrics by using the right contour (green) and an arbitrary contour (red, and manually created to show a larger error). In this paper, suppose that F is the set of edge points on the semi-automatically segmented region and G is one of the manually segmented region. The ACD values represent the average (L 1 norm) of distances on millimeter between the edge points on the manually segmented organ and the automatically detected organ. The minimum distance for the i th point in F to the set F is d F G i , thus the ACD values can be computed as in Equation 16. The number NF and NG were chosen to make sure that approximately an edge point was used by pixel.

ACD F G = 1 NF N F X i=1 d F G i , ACD GF = 1 NG N G X i=1 d GF i ACD = (ACD F G + ACD GF )/2 ( 1 6 )
Similarly the HD is considered as a L 1 norm. In addition, for the DC values, a higher Dice value assesses a better segmentation quality. On the 19 patients data, we got a similarity of over 90% for bladder and 80-90% for vagina and rectum especially since the ends of these two organs are di cult to define. In Figure 11, we illustrated some exceptions of the detection due to the unwanted e↵ect presented in the images. Generally the bladder is well highlighted under acquisition, while as the injected gel can exit the cavity with the patient pushing, blurs and incertitude will occur in the image for the vagina and the rectum (see Figure 11(b)). Because of the complexity of the shape, crease will occur in the 2D plane of the rectum (see Figure 11(a)) and intrinsic material generates spots in the rectum in the image (see Figure 11(c)). All these unwanted e↵ects would make the algorithm limited in certain cases. In these cases, the algorithm detected contours. However, the contours are not the anatomical contours of the organs. The quantitative results of our tests are shown in Table 1. An example of registration (Figure 8) was done on a quad-core desktop computer running at 3.10 GHz, with 16G memory. The computational time for one image is 53.8 seconds without counting the time consumed on updating the graphical user interface (GUI) bu↵er and reading/writing files and images (which took about 2.5 minutes). In this example, 4800 (400*12) discrete points were generated in the narrow band of each contour and it took about 400 iterations.

Conclusion and Perspectives

The correlation between geometric models and images remains a major challenge in the applications of medical imaging and computer-assisted diagnoses. In this paper, we present a B-spline based model-to-image registration approach for segmenting pelvic organs using MR images. The fitting procedure was guided by the minimization of an energy. The quality of fitting was ensured by a virtual image correlation approach and adaptive refinement. Thanks to the optimization Patient [START_REF] Berendsen | Free-form image registration regularized by a statistical shape model: application to organ segmentation in cervical MR[END_REF] Endometriosis ( 8) Prolapse ( 8 process the B-spline curves used to model the organs are registered on each image, without any intersection. The models provide a C 2 continuity (smooth property) for the contours. This method was developed and tested for segmentation of organs in 2D images.

The first perspective of this work is to assess the approach with more data and particularly with pathological images. We plan to lead a small clinical test to evaluate the interest of this new image processing technique in the diagnostic routines. Another possible use of this technique is tracking the motion of pelvic organs in dynamic MR images. The quantification of the displacements and the localization of the pathological zones related to the suspension devices between the organs could thus be possible, which could help in early diagnostic of pelvic organ prolapse. Finally, despite of the fact that the performance of the method is shown merely on 2D situation, it will be interesting to extend this algorithm for 3D reconstruction of multi-organ based on 3D images under the same framework. This would be more valuable for the purpose of patient-specific medical simulation.

Fig. 1 .

 1 Fig.1. The overview of our method (Compared to the image correlation, the deformable image is replaced by a virtual image generated repeatedly from the model. The cost function measures the similarity between the two images in our method.)

Fig. 2 .

 2 Fig. 2. Sagittal dynamic MR image (spatial resolution = 1.17mm/pixel, 256*256 pixels)

Fig. 3 .

 3 Fig. 3. The B-spline basis functions

Fig. 4 .

 4 Fig. 4. Cubic B-spline curve -Change of the point 3 influences 4 spans of the curve

Fig. 5 .

 5 Fig. 5. Choice of the levelset function with di↵erent configuration

FFig. 6 .

 6 Fig. 6. (a) Illustration of the virtual image computing (b) Virtual image generated on narrow bands ⌦ of the three organs (c) Virtual image superimposed on the real one.

Fig. 7 .

 7 Fig. 7. Contour detection with the cost functions. (a) Three small circles initialized by user with 10 control points for each (b) Positions after an a ne transformation (c) Final contours with optimal number of control points. (d,e,f) Illustration of the cost function metrics Eimage along the organ contours: [Ir(u, v) Iv(u, v)] 2 , with u evolves in the tangent direction of the contour and v in the direction perpendicular to the contour (see Virtual Image Section).

Fig. 8 .

 8 Fig. 8. Illustration of the convergence of the gradient descent method.

Fig. 9 .

 9 Fig. 9. Measurement of the detection error between a correct contour (green) and an arbitrary contour (red): (a) Hausdor↵ distance = length of s (b) Average curve distance = average length of s (c) Dice coe cient = s/(s + t).

Fig. 10 .Fig. 11 .

 1011 Fig. 10. Illustration of results: green curve shows the manually segmented organ as the reference, red curve shows the result of semi-automatic algorithm, yellow line shows the Hausdor↵ Distance. (a) Witness patient (1.17mm/px) (b) Patient with endometriosis (0.49mm/px) (c) Patient with prolapse (1.17mm/px).

Table 1 .

 1 Summary table for results of the error measures between manual and automatic segmentation -The three sub-tables present results of each organ in di↵erent pathological cases.

	) Witness (3) Average