Anatomical structure modeling from medical images

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 8 2 ( 2 0 0 6 ) 203–215

journal homepage: www.intl.elsevierhealth.com/journals/cmpb

Anatomical structure modeling from medical images Neculai Archip a,∗ , Robert Rohling b , Vincent Dessenne c , Pierre-Jean Erard d , Lutz Peter Nolte e a

Computational Radiology Laboratory, Harvard Medical School, Brigham and Women’s Hospital, Children’s Hospital, Boston, MA, USA b Department of Electrical and Computer Engineering, University of British, Columbia, Canada c Heraeus GmbH, Switzerland d Computer Science Institute, University of Neuchatel, ˆ Switzerland e Maurice E. Mueller Research Center for Orthopaedic Surgery, Bern, Switzerland

a r t i c l e

i n f o

a b s t r a c t

Article history:

Some clinical applications, such as surgical planning, require volumetric models of anatomi-

Received 29 December 2004

cal structures represented as a set of tetrahedra. A practical method of constructing anatom-

Received in revised form

ical models from medical images is presented. The method starts with a set of contours

21 January 2006

segmented from the medical images by a clinician and produces a model that has high

Accepted 1 April 2006

fidelity with the contours. Unlike most modeling methods, the contours are not restricted to lie on parallel planes. The main steps are a 3D Delaunay tetrahedralization, culling of non-

Keywords:

object tetrahedra, and refinement of the tetrahedral mesh. The result is a high-quality set

3D reconstruction

of tetrahedra whose surface points are guaranteed to match the original contours. The key

Tetrahedral mesh generation

is to use the distance map and bit volume structures that were created along with the con-

Volumetric modeling

tours. The method is demonstrated on computed tomography, MRI and 3D ultrasound data.

Delaunay

Models of 170,000 tetrahedra are constructed on a standard workstation in approximately

Contours

10 s. A comparison with related methods is also provided. © 2006 Published by Elsevier Ireland Ltd.

1.

Introduction

The problem of 3D model reconstruction arises in medical imaging where anatomical structures are scanned by a multislice imaging modality, such as computed tomography or magnetic resonance imaging. Such reconstructions of human organs can be used for surgery planning and simulation, prosthesis milling, radiation therapy planning and volumetric measurements. A wide range of approaches have been taken. The methods can be roughly classified into two groups: surface reconstruction and volumetric reconstruction. Many of these methods use a set of cross-sectional images as the input, and operate either on the original pixel data, or a description of organ boundaries obtained by segmentation. The boundaries

∗

of the structures in the cross-sectional slices appear as one or more non-intersecting planar contours. Methods that simply join contours together encounter problems because of the absence of any concrete a priori knowledge about the imaged structure. For example, branching and correspondence problems often arise [1,2]. Most of the current contour-joining approaches reconstruct the solid connection between two consecutive slices and then concatenate them in series to model the object. Much of the earlier work assumes a single contour without any holes on each slice, thus concentrating primarily on the problem of generating a non-self-intersecting surface under optimization of a chosen criterion. Representative works appear in [3–5]. Various optimization criteria, such as minimum surface area,

Corresponding author. E-mail address: [email protected] (N. Archip).

0169-2607/$ – see front matter © 2006 Published by Elsevier Ireland Ltd. doi:10.1016/j.cmpb.2006.04.009

204

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 8 2 ( 2 0 0 6 ) 203–215

maximum volume and minimum edge length are considered while constructing the surface, but none guarantee a non-selfintersecting surface. In later work [6] this aspect is addressed by exhibiting two contours that cannot be triangulated without self-intersection, unless new vertices are inserted on the contours. Subsequent approaches allow extra points on the contours and concentrate more on the problem of correspondence and branching with multiple contours on slices. Several approaches that deal with multiple contours on two slices are proposed in [1,2,7–11]. The methods proposed by [2,9,11] only work for restricted cases as pointed out in [1]. The algorithm in [7] is reported to work well in practice, based on the substantial experimental evidence presented. However, it requires that the contour boundary be discretized sufficiently finely, and given an arbitrary dataset, it is not clear how to determine a good discretization quantitatively. The method in [1] combines the approaches of [7] and [10] to handle the branching problem. Pioneering work on surface reconstruction from unorganized set of points is presented by Hoppe in [12], followed by [13–15] with methods related to the use of implicit surfaces in 3D reconstruction. These are sophisticated techniques but are usually much slower than contour-joining techniques. Marching cubes [16] is another popular algorithm for generating isosurfaces from volumetric data, with subsequent improvements in computational efficiency [17]. Another class of algorithms is based on the idea of deforming an initial approximation of a shape, under the effect of external forces and internal reactions. These algorithms are usually called deformable models. 3D deformable models have been recently used to perform the surface reconstruction of anatomical structures from medical images. Different representations have been used to fulfill different 3D modeling needs, from deformable 3D lines [18,19] to deformable volumes [20,21]. Deformable models have been introduced by Kass et al. in 2D as explicit deformable contours [22] and generalized to the 3D case by Terzopoulos et al. [23]. Parametric representations such as superquadratics [24,25] and discrete representations [26] have also been proposed for deformable models. Implicit representations have been used due to their ability to handle topology changes [27,28]. Pentland and Schlaroff [29] adopted an approach based on the finite element method and parametric surfaces. They start with a simple solid model and attach virtual springs between each data point and a point on the surface. The equilibrium condition of this dynamic system is the reconstructed shape. They showed how the set of parameters that describe the recovered shape can be used in object recognition. Other physically based approaches are described in [30,31]. Surveys on deformable models can be found in [32,33]. Level set methods were developed by Osher and Sethian [34] and introduced in the medical imaging community by Malladi et al. [35]. Level set methods are thoroughly described in [36]. The main idea of level sets is to embed the deformable model in a higher dimensional space. The surface is represented as the zero level set of a function. The algorithm starts with an initial surface that subsequently evolves under the guidance of a partial differential equation involving the distance field of the initial surface. The main advantage of level set methods is to allow changes of surface topology implic-

itly. The surface may split in several connected components or merge from several components. The level sets have also been recently used to create tetrahedral meshes [37]. Their main drawbacks are the computational cost, and the lack of user interactivity. Advancing front methods start with a boundary discretization and march a “front” inward forming new elements attached to the existing ones [38]. Advancing front techniques conform well to the boundary. This renders them a useful technique when the specific polygonal boundary representation of the geometry must be matched precisely. When the input geometry is not a polygonal boundary, a triangulation of this boundary must first be performed. The quality of this surface triangulation has a large impact on the three-dimensional algorithm’s behaviour. Poorly shaped surface triangles will engender ill-shaped tetrahedra [39]. A central decision in an advancing front algorithm is the placement of an interior point that marches the front further into the interior of the object. Local element control is possible because new nodes are created at the same time that new elements are created. The node and element creation is done as needed based on local procedures. Others have experimented with various metrics and criteria to evaluate the placement of the new node [40,41]. Most of the advancing front techniques have difficulty when fronts merge, which can occur very near important regions with high curvature [42,43]. An important class of algorithms for model reconstruction are those based on Delaunay triangulation. These algorithms can be classified in two sets. First, volume oriented algorithms output the boundary of a set of tetrahedra [44–47]. Second, surface oriented algorithms output a set of explicitly selected triangles [48,49]. Although easy to use, 3D Delaunay tetrahedralization can produce flat sliver tetrahedra when the data points are very irregularly spaced. Shewchuk provides an overview of the 3D Delaunay methods and he discusses how the worst slivers can often be removed [50]. Techniques for sliver removal can be also found in [51,52]. This brief summary shows the need for mesh models and the diversity of choices available to medical researchers. Researchers in surgical simulation, medical device design and biomechanics of movement all require detailed models of the human form. Given the complexity of biological materials, these models are most often obtained from actual subjects via non-invasive scanning technologies [53–55]. This data is usually voxelized or pixelized and must be converted into a more suitable format. In [53] and [55] MR images were manually labeled to create volumetric muscle meshes for accurate muscle length and moment arm computations for gait correction surgery diagnosis. In [54] the segmented visible human dataset was used to create volumetric muscle from which the researchers designed simplified muscle models for computing accurate muscle paths in the upper extremity. Volumetric elements (such as tetrahedra) are also necessary to allow proper interaction tasks, for instance when cutting a virtual model in surgical training systems [56]. In a typical clinical application (such as radiotherapy planning), contours are drawn around anatomical structures on individual slices. Navigation through the examination is done slice-by-slice. Spatial orientation, anatomical knowledge and semi-automatic segmentation tools are often used to aid con-

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 8 2 ( 2 0 0 6 ) 203–215

tour extraction. This approach allows the clinician to fully control the desired model. Time is a clear constraint since the radiologist or the oncologist cannot wait for more than a few seconds for model reconstruction. Moreover, the 3D model of an anatomical structure should have a high degree of fidelity with the contours, be able to handle complex structures, be efficient to compute and be applicable to non-parallel contours, such as those obtained from 3D ultrasound. This paper focuses on the particular problem of creating 3D mesh models of anatomical structures from segmented sets of medical images. Unlike the majority of related methods which focus on CT or MR images (based on parallel slices), we introduce a general-purpose technique that can also work for 3D ultrasound data (non-parallel slices). Moreover, our intention is to introduce a practical method which fits exactly the current clinical protocol. In other words, it should be based on the common method of contour-based segmentation which is still a widely used starting point for clinicians. For example, the above mentioned deformable models require a 3D initialization. If the iteration deformation does not converge to a desired location, it is not clear how to change it, to improve the result. Moreover, most of the techniques are computationally expensive (e.g. the marching cubes decimation). A simple technique that fits the clinical protocol is the one introduced by Boissonnat [8] and refined further in [10]. These are faster approaches based on Delaunay triangulation. However, in these approaches, numerical problems arise, and relatively small datasets sizes are handled. Another limitation is their reliance upon having contours lie on parallel slices. Delaunay triangulation approaches have also been criticized for producing meshes with poor aspect ratio when based on relatively sparse medical images. Yet, modern MR, CT and US scanners can produce finely spaced images so this concern is partially eliminated. Moreover, mesh refinement algorithms have been developed recently to help further. In light of these recent advances, the natural benefits of Delaunay triangulation based approaches deserve more investigation. We follow the same idea as Boissonnat and we will show that with suitable enhancements it can work well for this particular type of medical problem, even for complex structures. The goal is to create an easy-to-use method without requiring iteration, size restrictions, or tunable parameters that often limit the widespread use of other modeling methods. Previous work by the authors looked at producing compactly-represented models using the Delaunay technique and contour smoothing emphasing on efficiency issues [57] and model accuracy in [58]. Results were given for parallel contours with an emphasis on minimizing computation time without a guarantee of fidelity. The previous work also reported inaccuracies due to mesh culling approximations. In this paper, a new algorithm is proposed whereby the construction and refinement of the tetrahedra that fall along the boundaries is improved. Mesh refinement is also performed to make the mesh suitable for finite-element techniques. Results are also shown for both parallel and non-parallel contours. To compare our method with current 3D volumetric reconstruction methods, we consider three well-established methods based on Delaunay triangulation. The first one, Nuages was developed by Boissonnat [8] and Geiger [10]. The input contours must be located on parallel slices. It creates both surface and volumetric models.

205

The medial axis is used in the sculpting process. The second method proposed in [59] under the name of Cocone, is intended to create only surface models from a cloud of points lying on a surface. The third method was developed at INRIA [60]. It creates surface models from a cloud of points. Finally, our

Fig. 1 – Example of Delaunay triangulation obtained from Voronoi diagram: (a) initial Voronoi diagram for 4 points; (b) the Delaunay triangulation for the same points, obtained from Voronoi diagram shown in (a); (c) the empty sphere criterion is respected.

206

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 8 2 ( 2 0 0 6 ) 203–215

Fig. 3 – The points of the tetrahedron which are considered in the sculpting step.

and holes. Then, the cloud of points (the sum of points on all contours) is used to generate a convex hull based on a 3D tetrahedralization. The mesh structure is further culled, such that non-object tetrahedra are removed. The accuracy for some tetrahedra (partially inside and outside the modeled structure) is then improved. Mesh refinement is also considered for finite element methods. Details of the tetrahedralization, culling steps, and mesh refinement are now described.

2.1. Convex hull generation based on 3D Delaunay tetrahedralization Similar to the approaches of Boissonnat [8] and Geiger [10], we employ a 3D Delaunay tetrahedralization to create the convex hull of the cloud of points obtained in the segmentation step. Fig. 2 – (a) Contour from the skull model and (b) distance map. Each color corresponds to a different distance.

method was compared to an implicit surface modeling algorithm based on radial basis function interpolation [14].

2.

Methods

The input is a set of contours obtained from a volumetric image dataset (computed tomography, magnetic resonance and 3D ultrasound), obtained by segmentation. The images are transferred using the DICOM protocol to a standard personal computer. The modeling process is as follows. First, segmentation of the target anatomical structure is performed by a radiologist. The procedure is done slice-by-slice manually, but semi-automatic (such as active contours and region growing) or automatic segmentation methods [61–63] could be also used. With any of these methods, a distance map is calculated for every contour once available. In most of the cases, the object is non-convex, containing ramifications

Fig. 4 – The points A, B, C, P1 , P2 , are points in the final triangulation. P1 P2 is a part of a contour from initial segmentation. It splits the triangle ABC in two parts. P1 and P2 are not part of the initial contour, which means that the Delaunay criterion (empty sphere) is still respected.

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 8 2 ( 2 0 0 6 ) 203–215

207

xi is associated to a cell Vi . A simple example is presented in Fig. 1. It follows, from the duality and from the general position hypotheses, that each 0-dimensional cell of the Voronoi diagram is the center of a sphere circumscribing an m-simplex. An important property that characterizes the Delaunay triangulation in Rm is that the sphere circumscribing an m-simplex does not contain any other point of A in its interior. Moreover, if there exists a sphere passing through k ≤ m + 1 points of A, which does not contain any other point of A in its interior, these k points are the vertices of a (k − l)-simplex of the Delaunay triangulation. A more detailed description of the Delaunay triangulation and Voronoi diagrams can be found in [64,65]. For its speed and robustness the Qhull implementation [66] is preferred and used here.

2.2.

Fig. 5 – The tetrahedron ABCD is separated by a plane (B C D ). AB C D is outside the object of interest, while B C D BCD is completely inside. The tetrahedron AB C D is removed from the mesh model, while B C D BCD is divided into the tetrahedrons DD B C , BB C D and BCDC . All three tetrahedra are kept in the mesh model.

We will briefly give some definitions for the Voronoi diagrams and Delaunay triangulations. A non-void set of points A = x1 , . . ., xn ∈ Rm is in general position if no affine subspace of Rm contains A and there is no sphere Sm−1 through a subset of A with m + k, k > 1 points. The Voronoi diagram for A is a decomposition of Rm in convex m-dimensional cells V1 , . . ., Vn with the following properties: 1. Each Vi contains one single point xi of A. 2. Given x ∈ Rm , x ∈ Vi ⇔ d(x, xi ) ≤ d(x, xj ), for every i = j, where d is the Euclidean distance. From the definition above, it can be shown that intersection of k, (2 ≤ k ≤ m + 1) Voronoi cells is either empty or a cell with dimension m − k + 1 contained in the diagram. A k-simplex is the convex hull of the k + 1 affinely independent points, e.g. a 3-simplex is a tetrahedron, a 2-simplex is a triangle. The Delaunay triangulation is the geometrical dual of the Voronoi diagram, that is, for each p-dimensional cell of the diagram, we can associate, in a natural way, an (m − p)-simplex of the triangulation. For example, in the two-dimensional case, each triangle (2-simplex) is associated to a vertex of the diagram, each edge of the triangulation (a Delaunay edge) is associated to an edge of the diagram, and each vertex

Mesh culling

The triangulation process generates a convex hull of tetrahedral meshes. Tetrahedra that are not inside the object are also contained in the convex hull. Depending on the anatomical structure, up to 40% of the tetrahedra not on the object target must be removed. The decision regarding which tetrahedra to remove is vital for preserving the correct model shape. Boissonnat in Ref. [8] uses a method based on the external medial axis of each polygon. This method is not robust to numerical approximation for complicated contours with a large number of points. Moreover, this method is based on the assumption that the contours are on parallel slices. We propose a method for sculpting the convex hull, using two data structures: bit volume (BV) and distance map (DM). These structures are often used in surface matching methods [67,68], and are adapted for use in our problem. Both structures are created in the segmentation step. The bit volume is a 3D structure with the same size as the original image data. Each bit volume voxel contains a label indicating it is either inside or outside the object. The distance map stores for each voxel the distance to the nearest contour. More precisely, let ˝ be the 3D object model and ˝i = (x, y)|(x, y, zi ) ∈ ˝ a finite set of cross sections. The distance maps at the levels z0 , . . ., zn is defined by


DMi (x, y) =

−dist((x, y)ı˝), if (x, y) ∈ ˝i dist((x, y)ı˝) otherwise

(1)

where ı˝i denotes the boundary of ˝i which is described by the contours and dist denotes the Euclidean distance within the slices. The algorithm is based on chamfer distance proposed initially by Borgefors in [69]. The assumption is that the distance can be computed only from the values at neighboring positions plus a mask constant. To speed up the computation, the multipass mask propagation is replaced by region growing [70]. Both bit volume and distance map structures are created during the segmentation step once a contour is available. For high speed computations the algorithm presented in [71] can be used. However, in the current implementation less than 2 s are typically needed to compute the distance map for a contour. For a clinician performing segmentation slice-by-slice, this delay should be acceptable. An example of a contour (part of the skull) and its distance map calculated are shown in Fig. 2.

208

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 8 2 ( 2 0 0 6 ) 203–215

Fig. 6 – The steps of the modeling method are demonstrated on a CT scan of the skull from 55 slices: (a) point cloud; (b) tetrahedra that are partially inside and outside the skull are highlighted; (c) tetrahedral mesh obtained after culling; (d) refined mesh; (e) surface rendering.

Each tetrahedron of the convex hull is then examined to determine whether it lies inside or outside the target object volume. A bijective function between each voxel of the triangulation tetrahedrons and each voxel from these two structures can be obtained. This provides a way to decide whether a tetrahedron is inside or outside the target object. For each tetrahedron, the mid-edge points and gravity center for each face are considered. In Fig. 3, the considered tetrahedron’s points are illustrated. A point P is said to be inside the anatomical structures represented by the bit volume BV if BV(P) = 1 or DM(P) ≤
, where
is the maximum √ error approximation of the distance map. In practice,
is 2

(because the Euclidean distance is used when computing the distance map). Let P1 , . . ., P10 be the points of the tetrahedron ABCD (displayed in Fig. 3). The tetrahedron ABCD is said to be contained in the anatomical structure given by BV if all Pi are inside the anatomical structure.

2.3.

Mesh refinement

However, some tetrahedra lie partially inside and partially outside the object. This problem is illustrated in 2D in Fig. 4 and in 3D in Fig. 6(b). The points A, B, C, P1 , P2 are points in the triangulation (tetrahedralization in 3D). P1 and P2 are from

Table 1 – Experimental results

CT skull CT pelvis CT spine US fetus leg MRI spleen

Number of points

Number of tetrahedra generated

Average aspect ratio (unrefined a–r)

Time (s)

22900 22300 25038 12960 7125

172500 135840 137484 61500 36693

0.9 (2.0) 1.1 (2.1) 0.8 (1.9) 1.0 (2.2) 0.9 (2.1)

10 9 10 8 8

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 8 2 ( 2 0 0 6 ) 203–215

209

Fig. 7 – (a) Spine mesh model; (b) pelvis mesh model.

the original contour. The line P1 P2 splits the triangle ABC in two parts. P1 and P2 are not points that denned the contour originally, which means that the Delaunay criterion (empty sphere) is still respected. To resolve this problem, these tetrahedra are subdivided into smaller tetrahedra and the external parts are removed. In Fig. 5 the tetrahedron ABCD is separated by a plane B C D . AB C D is outside the object of interest, while B C D BCD is completely inside. The tetrahedron AB C D is removed from the mesh model, while B C D BCD is divided into the tetrahedra DD B C , BB C D and BCDC . All three tetrahedra inside the object are kept in the mesh model. At this point, the volumetric model is complete, but many long narrow tetrahedra may be present. The size and shape of the tetrahedra may not be ideal for some applications, such as finite element analysis [72]. A mesh refinement algorithm is therefore applied to improve the mesh quality. The method based on the circumradius-to-shortest edge ratio criterion is used here [72]. For the purpose of rendering the model, a surface representation is obtained by selecting the set of triangles lying on outer tetrahedra. The pseudocode for the reconstruction method is presented in Algorithm 1 and the modeling steps shown in Fig. 6.

3. Generate 3D Delaunay tetrahedralization for {Pi } 4. For each tetrahedra Ti decide using the distance map and bit volume, whether to add Ti on the model mesh representation M 5. If the Ti not completely inside the object, keep only the portion totally inside and re-triangulate 6. Refine the tetrahedrons based on circumradius-toshortest edge ration criterion 7. For each tetrahedral in M, each external face is added on the surface 8. Return the sets of tetrahedrals and triangles

The methodology presented in Section 2 is applied to both plastic phantoms and real patient image data. CT images are acquired with a Lightspeed scanner (General Electric, Milwaukee, WI) with a slice thickness of 3 mm and an inter-slice distance of 3 mm. The magnetic resonance images are acquired with Gyroscan NT Intera (Philips Medical Systems, Best, NL). The 3D ultrasound data is acquired by a Voluson 730 machine (General Electric). Several volumetric models are constructed to demonstrate the mesh quality and handling of complex structures.

Algorithm 1. 3D volumetric reconstruction Require: A volumetric medical image dataset and a way to segment the target object Ensure: Both surface and volumetric models of the object 1. Segment images to generate a set of contours {Ci } 2. Obtain the points {Pi } from {Ci }

3.

Results

Models of the pelvis and spine from parallel CT slices are shown in Fig. 7. A model of a fetal leg from non-parallel slices from 3D ultrasound is shown in Fig. 8. In all models, the tetrahedra are shown to have good aspect ratio and approximately

210

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 8 2 ( 2 0 0 6 ) 203–215

Fig. 8 – (a) Fetal leg reconstructed from ultrasound data (non-parallel slices); (b) tetrahedral mesh. equal size. Mesh refinement reduced the average aspect ratio (measured as the ratio between the circumradius and the shortest edge length) from 2.05 to 0.95 – see Table 1. A ratio greater than 2.2 is considered poor quality [73]. All points on the contours are guaranteed to lie on the surface by the

Delaunay criterion and the method of culling and refining the tetrahedra. The additional points used to handle subdivided meshes and improve mesh quality do not change the surface topology. The time needed to obtain a model with approximately 170,000 tetrahedra is less than 10 s on a standard workstation (2 GHz PC with 1 GB of RAM). The most time-consuming part of the algorithm is the 3D Delaunay tetrahedralization algorithm which has O(n2 ) complexity, where n is the number of points representing the contours. The culling process and mesh refinement have O(m) complexity (where m is the number of tetrahedra), so these steps are faster. More numerical results are listed in Table 1. Examples of various other reconstruction models are presented in Figs. 9–12. We first validate the surface models using the radial basis function (RBF) based algorithm introduced in [14]. This is a sophisticated commercial product with a high degree of accuracy, but relatively slow speed. The RBF algorithm is applied to the same set of contour points as used in Delaunay triangulation. The RBF approaches use a smooth interpolating function to create an implicit surface model. From the implicit surface a triangular mesh is extracted. This mesh surface is compared to the surface created by the Delaunay based method. The two-sided Hausdorff distance [74] is used to quantitatively measure the difference

Fig. 9 – Skull reconstruction. Contours are obtained by segmenting 55 slices. The contours contain 22,900 points from which 172,500 tetrahedrons are generated. The time needed for model generation is 10 s.

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 8 2 ( 2 0 0 6 ) 203–215

between the surface meshes for each object. Four types of errors are calculated: minimum, maximum, mean and root mean square. The numerical results are presented in Table 2. The mean error is less than 1 mm for all models. The error distribution map for the skull and fetus leg is presented in Fig. 13.

4.

211

Discussion

This paper introduces an automatic method (except for the segmentation phase) for building 3D geometrical models of

Fig. 10 – Spine reconstruction. Contours are obtained by segmenting 256 slices. The contours contain 25,038 points from which 137,484 tetrahedrons are generated. The time needed for model generation is 10 s.

212

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 8 2 ( 2 0 0 6 ) 203–215

Fig. 11 – Pelvis reconstruction. Contours are obtained by segmenting 119 slices. The contours contain 22,300 points from which 135,840 tetrahedrons are generated. The time needed for model generation is 9 s.

anatomical structures. The volumetric meshes are intended for use in applications related to interaction with virtual models, where the interior of an object is needed. The applications include surgical simulation, training, navigation, and planning. To compare our method with other Delaunay-based 3D volumetric reconstruction methods, we consider three wellestablished modeling methods. The first one, Nuages was developed by Boissonnat and Geiger [10]. The authors state that numerical errors may occur

in the medial axis computation for datasets with a large number of points. For a number of points smaller than 3000 the execution time is fast (approximately 2 s on a standard personal computer). Models with more than 3000 points produced numerical errors in the medial axis calculation algorithm. Typically, for objects with complex topologies (e.g. some of anatomical structures shown here) 3000 points are insufficient for an accurate representation, because the points are too sparsely sampled. For this reason, Nuages may not be suit-

Table 2 – Surface model validation between Delaunay based and radial basis function based methods Structure Skull Pelvis Spine Fetus leg Spleen

Max error (mm)

Mean error (mm)

RMS error (mm)

RBF execution time (s)

2.18 1.27 1.63 1.78 5.33

0.14 0.17 0.15 0.22 0.48

0.25 0.30 0.31 0.33 1.07

184 180 186 75 63

Error is calculated using the two-sided Hausdorff distance. Note that the execution time for the RBF models is only for creation of the surface model and not a volumetric model.

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 8 2 ( 2 0 0 6 ) 203–215

Fig. 12 – Spleen reconstructed from MR images. Contours are obtained by segmenting 24 slices. The contours contain 7125 points from which 36,693 tetrahedra are generated. The time needed for mesh generation is 8 s.

able for clinical applications with large models. Yet, for simpler shapes the models’ quality is satisfactory and execution time is fast. The second method proposed in [59] under the name of Cocone, is intended to create only surface models from a cloud of points lying on a surface. The authors state that their method encounters problems if the surface has very small fea-

tures where sampling is usually sparse, or non-smooth [59]. For a dataset with 23,000 points it needs 164 s to create a surface on a standard personal computer. It also created a surface containing several undesired holes and mis-triangulations. The execution time of the Cocone algorithm also depends on the topology of the object, not only on the number of points. The third method was developed at INRIA [60]. It creates surface models from a cloud of points. The authors state that if the sampling is not sufficient, their reconstructed surfaces do not have boundaries edges [60]. For a skull dataset of points with 23,000 points, it needs 22 s to generate a model. Surface irregularities and holes are also present in the generated model, probably due to local undersampling of the contours. Three major advantages arise with the new Delaunay based method: guaranteed fidelity of the models to the original contours, reasonable speed and robustness. It is also noniterative. First, the way the tetrahedra are removed from the convex hull gives us a good criterion for model accuracy. Even though tetrahedra partially inside and outside the structure of interest may be produced, these are detected and resolved. Second, we have shown experimentally that indeed the execution time needed to create both volumetric and surface models from 25,000 points is about 10 s on a personal computer. The execution time depends mainly on the number of points used, and less on the topology of the model. The most time-consuming part is the 3D Delaunay triangulation. The culling process that typically is the most consuming step with conventional approaches, is very fast in our case (O(n), with n the number of tetrahedrons in the convex hull). This is where most of the speed improvement is found. It is fair to say that the time savings is simply moved to the segmentation stage where the distance map is calculated. Yet, this is not noticed by clinician since it is calculated as a slice-by-slice basis with less than 2 s per slice. This amount of time is likely much smaller than the actual segmentation time. Third, our method follows a non-iterative algorithm, which means the clinician is not required to check that the model has converged to a correct solution. The model simply fits the extracted contours.

5.

Fig. 13 – The error distribution map based on the two sided Hausdorff distance, for the fetus leg and for the skull. The left side shows the histogram for the leg dataset.

213

Conclusions

A method for reconstruction of volumetric models of anatomical structures from segmented medical images is presented. Tetrahedral meshes are obtained, based on a 3D Delaunay tetrahedralization algorithm. Both surface (set of triangles) and volumetric models (tetrahedra) are obtained. A method for culling non-object tetrahedra is introduced that uses distance map and bit volume structures. Mesh accuracy improvement is performed. The model obtained is fidel to the anatomical structure contours created with segmentation algorithm. The meshes are also refined for use with finite element methods. CT, MRI and 3D ultrasound images are used to test the introduced method. The method is fast and robust enough to model various anatomical structures with complex shapes. The method uses data in the form that it is currently available in radiology departments and PACS systems. Overall, we consider the method suitable for clinical applications such as the image guided neurosurgery. In this specific application,

214

computer methods and programs in biomedicine

non-rigid registration algorithms based on advanced biomechanical models of the brain are used. A mesh model of the brain is therefore required. Currently we investigate the use this technology by enabling a prospective study during the neurosurgery.

Acknowledgments The authors wish to thank Francesco Termine, Stephanie Cheung and Hanson Kwong for their help on some implementation issues. We also wish to thank Tamal K. Dey for kindly providing us with the Cocone software that we have used to compare our method. We also thank the Prisme Research Group at INRIA for making their algorithms available online. This investigation was supported by NSF ITR 0426558, a research grant from the Whitaker Foundation, a research grant from CIMIT, by NIH grants R21 MH67054, R01 LM007861, P41 RR13218, U41RR019703, Natural Sciences and Engineering Research Council of Canada and the Institute of Robotics and Intelligent Systems.

references

[1] C. Bajaj, K. Lin, E. Coyle, Arbitrary topology shape reconstruction from planar cross sections, Graph. Models Image Process. 58 (1996) 524–543. [2] D. Meyers, S. Skinner, K. Sloan, Surfaces from contours: the correspondence and branching problems, Proc. Graph. Interf. (1991) 246–254. [3] H.N. Christiansen, T.W. Sederberg, Conversion of complex contour line definitions into polygonal element mosaics, Comput. Graph. 12 (1978) 187–192. [4] S. Ganapathy, T.G. Dennehy, A new general triangulation method for planar contours, Comput. Graph. 16 (1982) 69–75. [5] E. Keppel, Approximating complex surfaces by triangulation of contours lines, IBM J. Res. Dev. 19 (1975) 2–11. [6] C. Gitlin, J. O’Rourke, V. Subramanian, On reconstructing polyhedra from parallel slices, Int. J. Comput. Geom. Appl. 6 (1996) 103–122. [7] G. Barequet, M. Sharir, Piecewise-linear interpolation between polygonal slices, Comput. Vis. Image Understand. 63 (1996) 251–272. [8] J.D. Boissonnat, Shape reconstruction from planar cross-section, Comput. Vis. Graph. Image Process. 44 (1988) 1–29. [9] A.B. Ekoule, F.C. Peyrin, C.L. Odet, A triangulation algorithm from arbitrary shaped multiple planar contours, ACM Trans. Graph. 10 (1991) 182–199. [10] B. Geiger, Three-dimensional modeling of human organs and its applications to diagnosis and surgical planning, Technical report, 2105, France, 1993. [11] Y. Shinagawa, T.L. Kunii, The homotopy model: a generalized model for smooth surface generation from cross sectional data, Visual Comput. 7 (1991) 72–86. [12] H. Hoppe, T. DeRose, T. Duchamp, J. McDonald, W. Stuetzle, Surface reconstruction from unorganized points, Comput. Graph. 26 (2) (1992) 71–78. [13] E. Bittar, N. Tsingos, M.-P. Gascuel, Automatic reconstruction of unstructured 3D data: combining medial axis and implicit surfaces, Comput. Graph. Forum 14 (3) (1995) 457–468.

82

( 2 0 0 6 ) 203–215

[14] J.C. Carr, R.K. Beatson, J.B. Cherrie, T.J. Mitchell, W.R. Fright, B.C. McCallum, T.R. Evans, Reconstruction and representation of 3D objects with radial basis functions, in: E. Fiume (Ed.), SIGGRAPH 2001, Computer Graphics Proceedings, ACM Press/ACM SIGGRAPH, 2001, pp. 67–76. [15] Y. Ohtake, A. Belyaev, M. Alexa, G. Turk, H.-P. Seidel, Multi-level partition of unity implicits, in: Computer Graphics (Proceedings SIGGRAPH 2003), 2003, pp. 463–470. [16] W.E. Lorensen, H.E. Cline, Marching cubes: a high resolution 3d surface construction algorithm, Comput. Graph. 21 (1987) 163–169. [17] K.S. Delibasis, G.K. Matsopoulos, N.A. Mouravliansky, K.S. Nikita, A novel and efficient implementation of the marching cubes algorithm, Comput. Med. Imag. Graph. 25 (2001) 343–352. [18] A. Gueziec, N. Ayache, Smoothing and matching of 3-d space curves, in: Proceedings of the Second European Conference on Computer Vision, 1992, pp. 620–629. [19] G. Subsol, J.-P. Thirion, N. Ayache, A scheme for automatically building three-dimensional morphometric atlases: application to skull atlas, Med. Image Anal. 2 (1) (1998) 37–60. [20] G.E. Christensen, R.D. Rabbitt, M.I. Miller, Deformable templates using large deformation kinematics, IEEE Trans. Image Process. 5 (10) (1996) 1435–1447. [21] J.-P. Thirion, Image matching as a diffusion process: an analogy with maxwell’s demons, Med. Image Anal. 2 (3) (1998) 243–260. [22] M. Kass, A. Witkin, D. Terzopoulos, Snakes: active contour models, Int. J. Comput. Vis. (1988) 321–331. [23] D. Terzopoulos, A. Witkin, M. Kass, Constraints on deformable models: recovering 3d shape and nonrigid motion, Artif. Intell. 36 (1) (1988) 91–123. [24] E. Bardinet, L. Cohen, N. Ayache, Tracking and motion analysis of the left ventricle with deformable quadratics, Med. Image Anal. 1 (2) (1996) 129–149. [25] D. Terzopoulos, D. Metaxas, Dynamic 3d models with local and global deformations: deformable superquadratics, IEEE Trans. Pattern Anal. Mach. Intell. 13 (7) (1991) 703– 714. [26] H. Delingette, General object reconstruction based on simplex meshes, Int. J. Comput. Vis. 32 (2) (1999) 111– 146. [27] L. Lorigo, O. Faugeras, W.E.L. Grimson, R. Keriven, R. Kikinis, C.-F. Westin, General object reconstruction based on simplex meshes, Int. J. Comput. Vis. 32 (2) (1999) 111–146. [28] X. Zeng, L. Staib, R. Schultz, J. Duncan, Segmentation and measurement of the cortex from 3d mr images, in: Proceedings of the Medical Imaging Computing and Computer-Assisted Intervention (MICCAI ’98), 1998, pp. 519–530. [29] A. Pentland, S. Schlaroff, Closed-form solutions for physically-based modeling and recognition, IEEE Trans. Pattern Anal. Mach. Intell. 13 (1991) 715–729. [30] C.W. Liao, G. Medioni, Surface approximation of a cloud of 3d points, GMIP 57 (1) (1995) 67–74. [31] R. Poli, G. Coppini, G. Valli, Recovery of 3-d closed surfaces from sparse data, Comput. Vis. Graph. Image Process.: Image Understand. 60 (1) (1994) 1–25. [32] J. Montagnat, H. Delingette, N. Ayache, A review of deformable surfaces: topology, geometry and deformation, Image Vision Comput. 19 (4) (2001) 1023–1040. [33] T. Mclnerney, D. Terzopoulos, Deformable models in medical image analysis: a survey, Med. Image Anal. 1 (2) (1996). [34] S. Osher, J.A. Sethian, Fronts propagating with curvature dependent speed: algorithms based on hamilton-jacobi formulation, J. Comput. Phys. 79 (1988) 12–49.

computer methods and programs in biomedicine

[35] R. Malladi, J.A. Sethian, B.C. Vemuri, Shape modeling with front propagation: a level set approach, IEEE Trans. Pattern Anal. Mach. Intell. 17 (2) (1995) 158–174. [36] J.A. Sethian, Level Set Methods: Evolving Interfaces in Geometry, Fluid Mechanics, Computer Vision and Materials Science, Cambridge University Press, 1996. [37] N. Molino, R. Bridson, J. Teran, R. Fedkiw, A crystalline, red green strategy for meshing highly deformable objects with tetrahedra, in: Proceedings of the 12th International Meshing Roundtable, 2003, pp. 103–114. [38] J. Schoeberl, Netgen – an advancing front 2d/3d mesh generator based on abstract rules, Comput. Vis. Sci. 1 (1997) 41–52. [39] P. Moeller, P. Hansbo, On advancing front mesh generation in three dimensions, Int. J. Num. Meth. Eng. 38 (1995) 3551–3569. [40] S.H. Lo, Volume discretization into tetrahedra – i., verification and orientation of boundary surfaces, Comput. Struct. 39 (1991) 493–500. [41] S.H. Lo, Volume discretization into tetrahedra – ii., 3d triangulation by advancing front approach, Comput. Struct. 39 (1991) 501–511. [42] R. Garimella, M. Sherphard, Boundary layer meshing for viscous flows in complex domain, in: Proceedings of the Seventh International Meshing Roundtable, 1998, pp. 107–118. [43] R. Lohner, J. Cebral, Generation of non-isotropic unstructured grids via directional enrichment, in: Proceedings of the Second Symposium on Trends in Unstructured Mesh Generation, 1999. [44] H. Edelsbrunner, Weighted alpha shapes, Technical Report UIUCDCS-R-92-1760, Sept. Comp. Sci., Univ. Illinois, Urbana, IL, 1992. [45] N. Amenta, M. Bern, D. Eppstein, The crust and the skeleton: combinatorial curve reconstruction, Graph. Models Image Process.: GMIP 60 (2) (1998) 125–135. [46] J.-D. Boissonnat, Geometric structures for three-dimensional shape reconstruction, ACM Trans. Graph. 3 (1984) 266–286. [47] N. Amenta, S. Choi, R. Kolluri, The power crust, unions of balls, and the medial axis transform, Comput. Geom.: Theory Appl. 19 (2001) 127–153. [48] N. Amenta, M.W. Bern, Surface reconstruction by voronoi filtering, in: Proceedings of the Symposium on Computational Geometry, 1998, pp. 39–48. [49] N. Amenta, S. Choi, T.K. Dey, N. Leekha, A simple algorithm for homeomorphic surface reconstruction, Int. J. Comput. Geom. Appl. 12 (1/2) (2002) 125–141. [50] J.R. Shewchuk, Tetrahedral mesh generation by delaunay refinement, in: Proceedings of the Symposium on Computational Geometry, 1998, pp. 86–95. [51] S.-W. Cheng, T.K. Dey, H. Edelsbrunner, M.A. Facello, S.-H. Teng, Sliver exudation, in: Proceedings of the 15th ACM Symposium on Computational Geometry, 1999, pp. 1–13. [52] H. Edelsbrunner, D. Guoy, An experimental study of sliver exudation, Proceedings of 10th International Meshing Roundtable, 2001, pp. 307–316. [53] A. Arnold, S. Salinas, D. Asakawa, S. Delp, Accuracy of muscle moment arms estimated from mri-based musculoskeletal models of the lower extremity, Comput. Aided Surg. 5 (2000) 108–119. [54] B. Garner, M. Pandy, Musculoskeletal model of the upper limb based on the visible human male dataset, Comput. Meth. Biomech. Biomed. Eng. 4 (2001) 93–126.

82

( 2 0 0 6 ) 203–215

215

[55] A. Arnold, S. Salinas, S. Delp, Evaluation of a deformable musculoskeletal model for estimating muscle-tendon lengths during crouch gait, Comput. Aided Surg. 29 (2001) 263–274. [56] D. Bielser, M. Gross, Interactive simulation of surgical cuts, in: Proceedings of the Pacific Graphics, 2000, pp. 116–125. [57] N. Archip, R. Rohling, Fast reconstruction of volumetric models of anatomical structures, in: Proceedings of the IEEE Engineering in Medicine and Biology Society, 2003. [58] N. Archip, R. Rohling, Volumetric anatomical modeling from medical images, in: Proceedings of the IEEE Engineering in Medicine and Biology Society, 2004. [59] N. Amenta, S. Choi, T.K. Dey, N. Leekha, A simple algorithm for homeomorphic surface reconstruction, Int. J. Comput. Geom. Appl. 12 (2002) 125–141. [60] D. Cohen-Steiner, F. Da, A greedy delaunay based surface reconstruction algorithm, Technical report, INRIA, 2002. [61] N. Archip, P.J. Erard, J.-M. Haefliger, J. Germond, Lung metastases detection and visualization on ct images – a knowledge-based method, J. Visual. Comput. Animat. 13 (1) (2002). [62] N. Archip, P.J. Erard, M. Egmont-Petersen, J.-M. Haefliger, J. Germond, A knowledge-based approach for automatic detection of spinal cord in ct images, IEEE Trans. Med. Imag. 21 (12) (2002). [63] N. Archip, P.J. Erard, J.-M. Haefliger, J. Germond, A computer aided diagnostic system for radiotherapy planning, J. Med. Phys. 12 (2002) 252–259. [64] S. Fortune, Voronoi diagrams and delaunay triangulations, in: D.-Z. Du, F. Hwang (Eds.), Computing in Euclidean Geometry, Lecture Notes Series on Computing, vol. 1, World Scientific, 1992. [65] F. Aurenhammer, Voronoi diagrams – a survey of a fundamental geometric data structure, ACM Comput. Surv. 23 (1991) 345–405. [66] C.B. Barber, D.P. Dobkin, H. Huhdanpaa, The quickhull algorithm for convex hulls, ACM Trans. Math. Softw. 22 (1996) 469–483. [67] R. Baechler, H. Bunke, L.P. Nolte, Restricted surface matching – numerical optimization and technical evaluation, Comput. Aided Surg. 6 (2001) 143–152. [68] E. Gautier, R. Baechler, P. Heini, L.P. Nolte, Accuracy of computer-guided screw fixation of the sacroiliac joint, Clin. Orthop. Relat. Res. 393 (2001) 310–317. [69] G. Borgefors, Distance transformations in digital images, Comput. Vis. Graph. Image Process. 34 (1986) 344–371. [70] O. Cuisenaire, B.M. Macq, Applications of the region growing euclidean distance transform: anisotropy and skeletons, in: Proceedings of the ICIP, vol. 1, 1997, pp. 200–203. [71] A. Sud, M.A. Otaduy, D. Manocha, Difi: fast 3d distance field computation using graphics hardware, in: Proceedings of the EUROGRAPHICS, 2004. [72] J.R. Shewchuk, Delaunay refinement algorithms for triangular mesh generation, Comput. Geom.: Theory Appl. 22 (2002) 21–74. [73] S.-W. Cheng, T.K. Dey, E. Ramos, T. Ray, Quality meshing for polyhedra with small angles, in: Proceedings of the 20th Symposium on Computational Geometry, 2004. [74] N. Aspert, D. Santa-Cruz, T. Ebrahimi, Mesh: measuring error between surfaces using the hausdorff distance, in: Proceedings of the IEEE International Conference on Multimedia and Expo (ICME), 2002, pp. 705–708.