| | Segmentation of kidneys using a new active shape model generation technique based on non-rigid image registrationReceived 12 October 2007; received in revised form 16 July 2008; accepted 3 October 2008. Abstract Active shape models (ASMs) are widely used for applications in the field of image segmentation. Building an ASM requires to determine point correspondences for input training data, which usually results in a set of landmarks distributed according to the statistical variations. State-of-the-art methods solve this problem by minimizing the description length of all landmarks using a parametric mapping of the target shape (e.g. a sphere). In case of models composed of multiple sub-parts or highly non-convex shapes, these techniques feature substantial drawbacks. This article proposes a novel technique for solving the crucial correspondence problem using non-rigid image registration. Unlike existing approaches the new method yields more detailed ASMs and does not require explicit or parametric formulations of the problem. Compared to other methods, the already built ASM can be updated with additional prior knowledge in a very efficient manner. For this work, a training set of 3-D kidney pairs has been manually segmented from 41 CT images of different patients and forms the basis for a clinical evaluation. The novel registration based approach is compared to an already established algorithm that uses a minimum description length (MDL) formulation. The presented results indicate that the use of non-rigid image registration to solve the point correspondence problem leads to improved ASMs and more accurate segmentation results. The sensitivity could be increased by approximately 10%. Experiments to analyze the dependency on the user initialization also show a higher sensitivity of 5–15%. The mean squared error of the segmentation results and the ground truth manually classified data could also be reduced by 20–34% with respect to varying numbers of training samples. 1. Introduction  The segmentation of medical images is an important preprocessing step for further medical analyses or diagnoses. Many different methods are used for the classification of structures or organs of interest, ranging from active shape models (ASMs) to region growing or level set segmentations. Due to their inherent statistical regularization using prior knowledge, ASMs are very robust to leaking problems if adjacent image structures are not clearly delineated from each other. The core of ASMs is the statistical knowledge gained from the variations in the shapes that are extracted from the training input. Determining the correspondences between the shapes is therefore a crucial aspect for the generation of the entire model. Errors made during this phase directly lead to wrong statistical values for the shape variation that has to be known as exactly as possible to achieve a good regularization. A state-of-the-art solution consists of the optimization of the minimum description length (MDL) measure between a set of points (landmarks) placed on a parametric surface onto which all training shapes are mapped. Although this technique allows to process surface mesh representations of the training data with potentially different numbers of vertices, finding a solution for the problem is highly complex and – even on newest hardware – may be very time consuming. In fact, the procedure has to be rerun all over again if an existing model is updated with additional learning data. The mapping into the parametric space for the MDL optimization may also pose a problem for some applications, e.g. if a suitable transformation cannot be found or if the objects highly deviate from the parametric target space. From another perspective, the correspondence problem between the training shapes may also be regarded as an image registration task. Depending on the degrees of freedom of the spatial transform, the registration algorithm optimizes the match between corresponding structures within multiple images. The main contributions of this article consist of introducing a non-rigid registration step into the ASM generation phase to solve the point correspondence problem, formulating a suitable distance measure between the training shapes and a comparison of the proposed approach with an already established MDL method. This article demonstrates that the proposed approach yields more detailed models and results of high sensitivity that outperform a state-of-the-art MDL based technique. 2. Related work The clinical evaluation presented in this article focuses on the application of the generated ASMs to kidney segmentation from CT images. The application is primarily used as a means for the comparison of different algorithms to solve the point correspondence problem, however, there is a clinical need for such a segmentation system. Nephrologists are actually interested in some properties of the kidneys, e.g. the size, volume or perfusion. In literature, several approaches towards kidney segmentations can be found that make use of adaptive region-growing, knowledge-based or deformable models. The first class of algorithms is, for instance, described in Pohle and Tönnies [1]. They developed a region-growing algorithm that optimizes a homogeneity criterion automatically from characteristics of the area to be segmented. Kobashi and Shapiro [2] present a knowledge-based recognition system that utilizes the information about properties of the CT images and the anatomy. The very popular deformable model approach has previously been applied to kidney imaging by Tsagaan and Shimizu [3], [4]. Their proposed method is based on a non-uniform rational B-spline surface representation and prior knowledge about the shape of the organ, e.g. mean and variation, which is then incorporated into the objective function for the model fitting process as an additional energy term. Compared to rather heuristic approaches like region growing or threshold based segmentations that may easily leak into adjacent image structures, ASMs have been introduced by Cootes et al. [5] to achieve a regularization using prior knowledge about the statistical variation of the target shape. ASMs may therefore yield robust results also on images with low signal to noise ratio or at blurry organ boundaries. The incorporation of information about statistical organ characteristics additionally addresses the problem of segmenting structures from different patients. As introductory mentioned, the successful application of ASMs heavily depends on the quality of the solution for the point correspondence problem of the training data, which is the major challenge in the model generation phase. Difficulties arise from explicit representations of the surfaces as meshes with potentially varying triangulations or numbers of vertices. Davies et al. [6] published a description of an automatic method for the construction of optimal 3-D statistical shape models. The authors propose a technique to deduce correct point correspondences between different training shapes by estimating a parameterization of each shape surface that is optimal according to the description length measure [7]. Unfortunately, this technique does not establish dense and uniform correspondences across the set of training shapes and may result in poor modes of variation. Heimann et al. [8], [9] address this drawback by extending a parametric re-meshing technique that ensures uniformly distributed landmark positions across the training data. A MDL criterion is optimized to distribute point correspondences within the parameter space to optimally describe the statistical variations. Taking the another direction and realizing the idea of image registration to solve the correspondence problem is challenging. The training data is usually acquired from different patients, which requires a non-rigid spatial transform with a large number of degrees of freedom. In order to incorporate inter-patient shape variations, the applied registration algorithm has to estimate a non-rigid transform for each correspondence problem. Suitable non-rigid registration approaches can mainly be divided into parametric and non-parametric techniques. Parametric approaches incorporate an inherent regularization by the choice of the parametric model, whereas non-parametric methods have to be constrained by additional regularization terms. Comprehensive descriptions about this topic can be found, for instance, in the works of Modersitzki [10], Hermosillo et al. [11] or Clarenz et al. [12]. Among various non-rigid transformations, dense deformation fields provide the largest flexibility regarding the spatial transform. They allow to specify a translation vector for each image element. In the following formulation this leads to a highly ill-posed optimization problem that needs further regularization. Fischer and Modersitzki [13] proposed a curvature regularization for the usage of their non-rigid technique within the field of medical imaging. In the following sections, this non-rigid curvature regularized image registration is incorporated into the ASM generation to compute the shape variations between the meshes of the training set. The major steps of the proposed segmentation system for the evaluation are depicted in Fig. 1. Note that this article focuses on a solution for the point correspondence problem during the model generation phase, which is accentuated in the figure. 3. Methods The beginning of this section provides an overview of the main principles of ASMs followed by a brief description of a well established MDL criterion based approach. The core part of this section deals with a novel registration based method for estimating the point correspondences. 3.1. Basics of statistical shape models The point distribution model by Cootes et al. [5] forms the basis of an ASM that is built from a set of N training shapes. Each shape is represented by n sampled surface points. For 3-D type problems, the components of their position vectors  , k =1,…, n are used to create a representation for a single shape x as follows: first all the x-components of the n position vectors pk are given followed by all y-components and z-components, i.e.  . Denoting the shape representations xi, i=1, …, N, for the training set as column vectors, one can obtain a landmark configuration matrix L=[x1x2…xN]. Applying a principal component analysis (PCA) to L yields the principal modes of variation of the training data, i.e. the eigenvectors ej. The mutually orthogonal eigenvectors ej are sorted in descending order of their respective eigenvalues λj. If the number of eigenvectors is restricted to the vectors that belong to the T largest eigenvalues, a linear combination of these T principal modes of variation with the mean shape  spans the subset of shapes composed of the given modes of variation. A new shape x* contained within this subspace can therefore be expressed by where  is the weighting factor for the corresponding ith variation. A very important precondition to solve the PCA for the matrix L is the knowledge about the corresponding points of the training shapes. These correspondences determine the order of the point coordinates within the shape representations xi and also their dimensionality. It is therefore not possible to select points that occur only in a subset of the training shapes, which directly leads to gaps in L. If the information about the location of these points is not given by construction, computing the boundary conditions for the correspondences is a highly non-trivial task. Nonetheless, it is crucial for the quality of the resulting ASM. 3.2. Description length This section briefly describes an MDL criterion based method for creating 3-D statistical shape models, which was first introduced by Davies et al. [7]. This MDL technique leads to a reformulation of the point correspondence problem to find an optimal mapping of each training shape onto a sphere. This mapping is manipulated such that the description length of all points becomes a minimum. Thodberg [14] provides a simplified and more efficient version of the MDL approach that is based on a cost function F of the description length for the generated model being defined as follows: with where λm denotes an eigenvalue of L and m is the number of used modes of variation. λcut specifies a free parameter that represents the expected noise in the training data. Heimann et al. [8], [9] made use of Thodberg's cost function F and extended the commonly known MDL approach by a procedure that modifies the landmark positions locally while trying to preserve already established correspondences. All of the mentioned MDL methods apply a mesh parameterization of the object to be segmented. They implicitly assume that the target object can be topologically mapped onto a sphere. An overview of various mesh parameterizations can, for instance, be found in [15]. The entire parameterization approach has previously been described by Gu et al. [16] and adapted for the context of medical image processing by Heimann et al. [8]. It is worth mentioning that the cost function F is defined within the domain of the landmark points. Hence, the accuracy of the optimization correlates with the number of landmark points. Likewise, the more landmarks are used for the representation, the more demanding the computation of the correspondences will be. 3.3. Registration of shape models There exist several equivalent ways of representing the shape of an object. A shape may, for instance, be defined as a set of points that form a surface mesh, as it was discussed in the previous section. Or it can be described by a bi-valued function Φ in a discrete image domain Ω that is specified by spatial sampling properties. Let Ωx denote the image domain of the shape x, with Ωx⊂Ω. The corresponding discrete representation of x within this image domain is defined as a spatial region X∈Ωx with an appropriate resolution and size such that x∈Ωx. As X contains the discretization of x, Φ may be denoted as where Hence, X may be regarded as a discrete binary image of the shape x. Let y denote a different shape and Y its corresponding representation in the image domain. The point correspondence problem between x and y can now be formulated as the problem of finding a suitable spatial transform between X and Y that maps corresponding structures onto each other. That means for every point pj∈ Ωx a corresponding point pk has to be found with pk∈ Ωy. This is the classical image registration problem. The non-rigid registration used in this work yields a dense deformation field  that provides a solution for the following point correspondence problem: The transformation is represented as a spatial function u(p) that maps between the image domains of the shapes x and the x. It is worth noticing that this formulation of the point correspondence problem within the image domain does not require an explicit representation of an extracted shape. It can be applied to both implicit and explicit segmentation results, whereas the previously mentioned methods rely on an explicit surface representation. This is an advantageous property, since it allows to relate shapes that are composed of several closed sub-shapes (e.g. as a result of an implicit level set segmentation). Commonly, the training data is usually acquired by a supervised segmentation of the images and stored either as an image or a surface mesh. Both representations may be transformed into each other, however, some attention has to be paid to the discretization of the image domain in order to avoid a loss of substantial structural information due to undersampling. This article focuses on a regularized, non-parametric, non-rigid registration formulation based on the work of Fischer and Modersitzki [10] who proposed a framework which we will now briefly summarize. The registration problem can be stated as the search for a non-parametric mapping from one image domain into another, usually referred to as reference R and template space T with their corresponding image domains ΩR and ΩT. We will refer to the deformation field between R and T simply as u. An objective function has to be formulated with respect to u which accounts for the similarity between R(p) and T(p−u(p)) with p∈ΩR. In an intensity based energy formulation, this similarity is expressed by a distance measure  being explained later on, that is minimal if the mapping yields the best result between the two images: In general, this problem is ill-posed and the solution may not be unique or even continuous. An additional regularization term, called the smoother S, is introduced to address this drawback. With an appropriate regularization it is now possible to penalize transformations that do not seem to be suitable for the given application. Hence, the overall registration problem is to find a spatial mapping u that minimizes the joint functional  : where the curvature term (8) weighted by the scalar α∈ IR determines the amount of regularization applied to the deformation, i.e. it controls the smoothness of the resulting deformation field. The larger the value for α, the more rigid the registration will be. The transform u defines a dense deformation field that assigns a translation vector to each element of the reference image. The nonlinear optimization of the objective functional (7) is driven by a suitable distance measure. The calculus of variations is then applied to solve the minimization problem. The functional  yields a global energy and a minimum is obtained when small changes in the solution u do not increase the energy. This means that the Gâteaux derivative  vanishes for all perturbations  . The weak form of the global energy  is directly related to the Gâteaux derivative and leads to the Euler–Lagrange equation: with the matrix  being the discretized partial derivative operator resulting from  and f the derivative of the distance measure, also called force. A solution of this semi-linear partial differential equation fulfills the necessary condition for the minimization of  . Computing a solution u for (7) usually involves three steps: (a) an initial placement that results in an overlap between the image domains, (b) a rigid registration for rotational and translational parts and (c) the non-rigid registration to solve the point correspondence problem. In practice, there exist many possibilities for an initial placement, for example the alignment using the center points of the bounding boxes of the shapes or their centers of gravity, just to name a few. Step (b) is primarily required if the initial alignment has to be substantially improved in order to get a good starting point for the numerical optimization in (c). Besides choosing the smoother and a numerical optimization scheme, a crucial step of the shape registration approach is the selection of a suitable distance measure. Regarding the point correspondence problem, the distance measure has to operate on binary images that represent the shapes in the image domain. Points on the surface of the shapes have to be registered correctly and intrinsic properties (e.g. the curvature) have to be retained between corresponding surface regions as closely as possible. As the image domains are of the same modality and intensity, the straight forward sum of squared differences (SSD) distance measure can be used. Its definition for a specific point in the image domain is provided in the following equation: A drawback of the SSD in the context of finding shape correspondences is that it does not account for any surface properties. In order to account for this information in addition, we therefore propose a novel extended form of the pure SSD distance measure. It incorporates the surface property of the mean curvature of the shape at a specific point as used by Sethian [17]: The novel distance measure for the shape registration includes both the similarity between the discrete shape representations (10) and their surface curvature (12), abbreviated by CSSD: where κR( p) denotes the curvature of the reference image at position p∈ ΩR and κT( p− u( p) the curvature of the template image at the mapped position, respectively. As depicted in Fig. 2, the curvature extension results in a deformation field that maps between regions with corresponding surface properties. In contrast, optimizing the SSD distance alone results in the most efficient deformation field with respect to the regularization energy. For example, a standard SSD approach might just smooth out a bulge while the CSSD would try to match it with a corresponding bulge first. It is worth noting for implementation purposes that the curvature is calculated on the original (i.e. undeformed) template image and interpolated using the current deformation field within each iteration. Computing the curvature on the already deformed template image is incompatible with the objective of retaining surface properties. In practice, solving Eq. (11) can be very sensitive to noise due to the second order derivatives involved. Well known techniques for low pass filtering, e.g. with a Gaussian kernel, help to alleviate these problems. The derivatives should be computed by analytically derived versions of the corresponding filtering kernel in order to reduce noise. However, the kernel width has to be chosen with care in order to retain the details of the shape. The registration scheme for the point correspondence problem is depicted in Fig. 3. X1 is used as the reference, which is successively registered with the remaining N−1 shapes images of the training set. The resulting deformation field ui then specifies the point correspondences between the shapes X1 and Xi with i=2, …, N. The registration problems have to be formulated so that the shapes Xi are mapped onto X1. The shape X1 has to be chosen carefully and ideally it should be close to the mean of the training samples. As affine transformations are contained in the kernel of the curvature regularization [13], initial rigid mis-alignments up to a certain degree should be of no concern for the registration. A state-of-the-art multi resolution optimization strategy increases the global character of the entire deformation. The influence of the reference shape on the resulting ASM is therefore reduced. For the conducted experiments of this article, choosing a reference shape just by looking at the smoothness of the surface (mean curvature) worked fine for the registration of the entire training set. Finally, the estimated solution for the point correspondence problem is incorporated into the ASM model generation by creating a registered matrix L. In Section 3.1 we briefly described how this matrix is composed of the shape representations xi, i=1, …, N. However, this implies that the correspondence problem is already incorporated into these point distributions on the surfaces. This requirement can now be addressed by applying the deformation fields to the point distribution of the reference shape x1. If only the image domain representation x1 is known, its explicit surface representation can be created using a suitable mesh extraction technique (e.g. the Marching Cubes algorithm [18]). The point corresponding surface representations of the remaining set of training samples can now be achieved by transforming x1 according to the non-rigid deformation fields: withThe operator  denotes the interpolated application of the deformation field to the points of the surface representation x1. As the point correspondences have been established for all training samples using the novel registration approach, the PCA on the registered L can now be applied to compute the principal modes of variation. 3.4. Gray-level appearance model Once the correspondence problem is properly solved the statistical shape model can be generated incorporating the variations seen within the training data. The applied segmentation system (see Fig. 1) incorporates the gray-level appearance model introduced by Cootes et al. [19] to establish a link between the shape variations and the intensities within the images. In this context, a hierarchical multi-resolution approach is applied to improve the efficiency and robustness of the search for the boundary of the shape within the image space. Hence, it is necessary to examine the gray values along the normal direction at each registered model point for all training samples and each level of the resolution hierarchy. The resulting vectors of gray values with intensity samples along the positive and negative normal direction are normalized over all training samples to form a mean gray-level profile for each model point. A gray-level appearance model is thus given by the intensity distributions along the surface normals at the corresponding points for all training shapes. One challenge during the ASM segmentation is the spatial adaption of the model points to minimize a distance criterion between the learned and the intensities observed at their current locations. 3.5. Image search and segmentation The gray-level appearance model is placed into the image domain of the target object that is to be segmented by a user provided seed point. Usually, the model does not coincide with the target object due to variations between the patients and also the placements of the seed. Therefore, a displacement for each model point has to be estimated, that moves the surface of the model towards the boundary of the target shape. This movement is calculated for each point given its trained gray-level profile and a new vector of gray values computed for the current spatial location in the image. The similarity between the jth element of the gray-level appearance model dj and its current observation oj is assumed to be maximal if the model surface is located exactly at the boundary of the shape. The Mahalanobis distance is used to measure this similarity: where ∑ j is the covariance matrix of the jth model profile. It can be regarded as a distance measure that has to be minimized for each object location. The discrete displacement Δ j of the jth model point for each iteration step has to minimize the Mahalanobis distance (16): As the observed gray-level profiles are computed only for a specific sampling range, the problem usually does not have a closed form solution. A nonlinear optimization scheme is applied in order to adapt the model gradually to the shape of the new object. Within each iteration of the nonlinear optimization, the jth element of the gray-level appearance model is displaced by Δj along its normal direction to minimize (16). A corresponding displacement is calculated for each of the n points of the ASM. The nonlinear displacement of the model points minimizes (17), however, it does not incorporate the prior knowledge about the statistical variation within the training shapes. As a result, the displaced model has to be mapped back into the space that is spanned by the modes of variation (see also Eq. (1)). The linear part of this iteration step estimates the pose and shape parameters of the model in a least-squares manner with respect to the training set variations, i.e. determining the weighting factors for the eigenvectors. The computational details of the least-squares update scheme can be found in Cootes and Taylor [20]. Various numerical schemes may be applied to solve the nonlinear problem, for instance the first order gradient descent algorithm or higher order techniques like the quasi-Newton methods. 4. Results Three different algorithms to solve the point correspondence are compared to each other based on an evaluation using 3-D medical images. Each algorithm has been embedded into the described ASM framework. The medical data consists of 3-D abdominal CT images of kidneys from 41 different patients of mixed gender and age. The images have been acquired using two different Siemens CT scanners (Sensation 10 and Sensation 16) with a resolution in x/y/z from 0.6/0.6/5 to 0.75/0.75/5 (in mm) and provided in Dicom format. The volume sizes for the experiments range between 512×512×120 and 512×512×300. In order to evaluate differences in the point correspondence algorithms, all 41 kidney pairs have been manually segmented and approved by a nephrologist. In all experiments, the manual segmentation results are used as the gold standard. The entire set of labeled segmentation data is divided into two disjoint parts: one set of varying size between 10 and 20 is used for the training of the ASM and the remaining one for testing. The novel idea of using a non-rigid registration based algorithm to solve the point correspondence problem (see 3.3) is compared to an established MDL approach (see Section 3.2) from Heimann et al. [8]. For the non-rigid registration, two distance measures are applied: a straight-forward SSD criterion-based implementation using (10) and the novel surface curvature extended SSD (13). In the following text, we will refer to the resulting ASMs as MDL, SSD and CSSD. All ASMs have been created on the same training sets and evaluated on the same test data with equal initialization parameters. Table 1 provides a brief description of the properties for the MDL approach for one side of the kidney pairs. For the other two methods the surface mesh for x1 has been extracted from X1 with 2000 vertices using the Marching Cubes algorithm. The resulting ASM therefore consists of 2000 vertices, as well.  | Mean radius in voxels | 22 | | | Number of samples | 41 | | | Sample complexity for the MDL (# vertices) | 2000–3000 | | | Model complexity for the MDL (# landmarks) | 2562 | | | | |
4.1. Methods of evaluation The evaluation of the ASMs for the segmentation of kidneys is based on the comparison with the gold standard. The differences to the optimal segmentations are measured with respect to the generalization of the model and the error to the gold standard. The generalization assessment yields information about the ability of the model to adapt to a new shape that deviates from the training data. This depends on both the number of used eigenvectors and the diversity of variations within the training samples. A set of models with different numbers of training samples is constructed with all three approaches. A series of leave-one-out cross validations on the testing data is used to measure the distance to the gold standard using two measures: the mean squared error (MSE) and the sensitivity. The MSE is used to estimate the error between the segmentation result of the ASM and the corresponding gold standard. This measure depends on the resolution that is used for the discretization. Therefore, the same image parameters have been applied that were used for the manual segmentation to achieve the gold standard (i.e. the same image resolution, size, position and orientation). In the following formula, this is simply denoted by the domain Ω: Here, X and Y are the image domain representations of two segmented shapes x1 and y. X is always a gold standard segmentation and the second input is the result from the corresponding model. The sensitivity is used as the second evaluation criterion, which is defined by where TP is the number of true positives and FN the number of false negatives. In our case, TP is given by the number of voxels that are segmented consistently as kidney tissue both by the specific ASM and within the gold standard. FN is the number of voxels that have been falsely classified as background. The sensitivity is used to measure the conditional probability that the given shape model has classified kidney structure that actually belongs to the kidney according to the gold standard. Therefore, this measurement states how well the shape model generalizes to new kidneys that are not contained within the training set. In literature, the specificity is often used as an additional criterion for quantifying segmentation results. In practice, however, there is a problem with this particular measure, due to a normalization issue: if the background of a segmented shape is arbitrarily enlarged, the specificity values can be driven close to 1.0. This is due to the fact that background voxels are usually recognized as lying outside the structure anyway. For this reason, the specificity is not used as a criterion in our evaluation. 4.2. Experimental results This section presents the segmentation results for the generated ASMs. All experiments have been performed on a Pentium 4, 2.8GHz with 2GB of main memory. A single registration for both the SSD and CSSD approach took approximately 8min on the given hardware. The complete registration for the largest training set with 20 kidneys took approximately 170min, compared to 16–20h using the MDL approach. The runtimes refer to a non-optimized implementation of the algorithms and are only used as a rough indication of the efficiency of the ASM generation algorithms. As several tasks may be performed in parallel, further improvements can be achieved by utilizing multi-core processor architectures or graphics hardware. Since the presented segmentation system (see Fig. 1) has been initialized by a seed point, the experimental phase was divided into two parts. The ASMs have firstly been tested for sensitivity according to varying seed point locations on one test sample volume. Secondly, the seed points were placed ideally into the center of gravity of the corresponding test images in order to analyze the models with respect to a larger test set and different numbers of training shapes. The training and test samples, as well as the parameters for the ASM, were the same for all three models throughout the corresponding experiment. As mentioned in Section 3.4, a multi-resolution technique has been applied to increase the attraction range for the optimization and the efficiency. The numerical convergence criterion for each level of detail is based on the variation of the MSE of the shape between two subsequent iterations. In general, the image search algorithm converged in less than 30s, where at most 70 iterations in each level were allowed. Analyzing the results, 30 iterations for one level were generally sufficient to achieve numerical convergence. 4.2.1. Sensitivity to seed point variations As the ASM segmentation relies on a manually specified seed point, a series of tests has been applied to evaluate the robustness of this initialization. The proposed ASMs are evaluated for various seed positions with respect to the MSE (18) and the sensitivity (19). The modes of variation have been chosen to cover 99.9% of the training set. 13 different positions for the seed point have been chosen to reflect typical user inputs. From the test set, one left and right kidney were chosen for this experiment. Results for the left and right kidney ASMs generated from 20 training samples are presented in Table 2, Table 3. As can be seen, the registration approaches yield better results compared to the MDL method, with a slight improvement with the curvature extension. According to these experiments, the models generated by the registration approach are less sensitive to variations in the seed points. For the right kidney model, the sensitivity results were additionally compared between a learning set size of 10 and 20. The values for the sensitivity presented in Table 3 between the different size of the training set are also illustrated in Fig. 4. The generalization ability of the models clearly benefits from a larger training set. Regarding the relative increase in performance due to a larger set of training shapes, the registration approaches seem to exploit the learning data more efficiently than the MDL model. Fig. 5 shows an example from the test set used to generate the results for varying seed locations. It illustrates the progress on one slice through the volume during the segmentation of a right kidney using the CSSD model trained with 20 samples. The figures show the initial placement of the mean shape, intermediate results after 10 and 20 iterations and the finally converged segmentation. Fig. 6 presents a 3-D visualization of the corresponding segmentation result. | | | | Initialization sensitivity—left kidney | |
|---|
| N=20 | MDL | SSD | CSSD | |
|---|
| S.E. | 0.67±0.12 | 0.94±0.02 | 0.95±0.01 | | | MSE | 100,490±33,599 | 22,991±3545 | 22,835±3046 | | | | |
| | | | Initialization sensitivity—right kidney | |
|---|
| N=10 | MDL | SSD | CSSD | |
|---|
| S.E. | 0.69±0.18 | 0.74±0.13 | 0.74±0.12 | | | MSE | 79,330±44,930 | 61,064±33,263 | 60,084±33,043 | | | | |
| N=20 | MDL | SSD | CSSD | |
|---|
| S.E. | 0.77±0.15 | 0.91±0.01 | 0.91±0.01 | | | MSE | 54,680±37,344 | 13,866±2159 | 13,501±1098 | | | | |
4.2.2. Cross validation segmentation results In the second part of the experimental results, the 41 kidney pairs have been divided into two disjoint groups using a cross validation strategy. Each method for the model generation has been used to generate ASMs from 7, 10, 15 and 20 different training samples. The remaining samples have been used as the test set. Considering the MDL, SSD and CSSD approaches, this results in a total of 12 models. Compared to the previous tests where the initial placement was varied, the seed point is now automatically placed into the center of gravity of the gold standard segmentation of the corresponding test sample. This initializes the image search and segmentation phase (see also Fig. 1) in an equal manner for all compared models. Fig. 7, Fig. 8 show graphs of the sensitivity and the MSE for models from different numbers of training samples. The corresponding numerical values for each data point are given in the Table 4. It is clear from these results that the ability of an ASM in general to adapt to a new shape depends on the number of training samples used for its creation. The graphs show that both registration approaches deliver better results for models trained with 10 or more samples. In the case of seven training samples, the SSD distance measure performs better than its competitors. However, for 10 or more models, the values for the CSSD approach lead to better results. The graph in Fig. 8 illustrates that the MSE made with the models created using the registration approaches is smaller for all numbers of training samples. Considering ASMs trained with 10, 15 or 20 training samples, the MSE values decrease almost linearly for SSD and CSSD, whereas it seems that the MDL mesh is not capable of further improvements with more learning data. One reason for this effect may become clearer by comparing the mean shapes between the CSSD and the MDL models created from 20 samples. Fig. 9 illustrates that the CSSD model provides more morphological detail on concave parts of the surface compared to the MDL approach. Although the surface of the MDL model features the same number of vertices, the distribution of points on the surface of the CSSD model leads to a better representation of the surface properties of the kidneys. The higher detail of the registration based meshes and the resulting accuracy in concave areas is visualized in Fig. 10, which shows a comparison between CSSD and MDL for a typical segmentation result from the test set. Here, the displayed surface outlines the corresponding gold standard segmentation. The colors indicate the distance between the ASM segmentation to the gold standard given in mm. The registration approach features a higher accuracy at detailed surface structures. The resolution of the MDL mesh is much coarser at concave surface areas, for example at areas around the renal hilus, which leads to larger errors. | | | | # | MDL | SSD | CSSD | |
|---|
| Sensitivity | | | 7 | 0.8 | 0.81 | 0.78 | | | 10 | 0.88 | 0.9 | 0.92 | | | 15 | 0.89 | 0.92 | 0.92 | | | 20 | 0.88 | 0.91 | 0.91 | | | | | | MSE | | | 7 | 40,037 | 26,366 | 31,930 | | | 10 | 18,892 | 15,883 | 15,162 | | | 15 | 16,955 | 13,425 | 13,357 | | | 20 | 19,820 | 12,635 | 13,043 | | | | |
To get an impression of the statistical information contained within the training samples, Fig. 11 illustrates the variations along the three most significant principal modes. The model has been created from 20 training samples using the CSSD approach. 5. Discussion A novel approach has been presented to solve the point correspondence problem between the training samples for an ASM segmentation by the application of non-rigid image registration techniques. For comparison, an already established method using an MDL criterion for this task has been applied. Results have been presented for the automatic segmentation of kidneys from CT images. Compared to an already established MDL method, the registration based solution for the correspondence problem not only allows to include more surface detail within the model itself, but it also provides the ability to easily update the model with new training samples. In comparison, the objective function for the MDL formulation contains dependencies to all model points, which leads to a high computation effort if new training data has to be incorporated. This effort is drastically reduced if the registration approach is applied for the learning. The conducted experiments show that the accuracy of the models generated with both registration approaches lead to slightly more accurate segmentation results compared to the MDL mesh. The overall performance of ASMs for the task of kidney segmentation was very good. The interaction with the algorithms is rather simple and reduces to just placing a seed point as the desired start position for segmentation. This leads to a high acceptance among physicians using ASM based segmentation applications e.g., for perfusion analyses, size or volume measurements or to acquire statistical information. Even if the contrast of the image at the organ boundary is very low, which may happen if the kidney and the liver are touching each other, the image search algorithm regularized by the ASM is robust enough to handle these cases and does not leak into adjacent tissue. Moreover, the image search algorithm is based on a multi-level resolution approach that delivers a large attraction range even if the seed point is misplaced. The execution time of a full segmentation run currently takes less than 1.5min, which may still fit into the time-frame of a physician's examination without causing unacceptable delay. Nonetheless, the accuracy of the specific ASM depends on the statistical information of the training data. The point correspondence problem therefore has to be solved accurately. Parametric models to describe the shapes or mappings into parametric spaces may sometimes be too restrictive. In case of the kidneys, the concave parts pose the biggest problem for a parametric mapping onto a sphere in order to solve the point correspondences using an MDL formulation. Despite these structural problems, the advantage of the MDL method is the opportunity to optimize the location of the points in order to cover the most statistical variation. This is not yet incorporated into the registration based approach, however, the resulting deformation fields in the image domain already provide the necessary information to estimate such a distribution. 6. Conclusion The novel approach to solve the point correspondence problem using a non-rigid, curvature-based image registration provides an attractive alternative to MDL based techniques. It leads to models with higher segmentation accuracy and a considerably faster model generation. In practice, the need for updating an established model with new training data may easily be met, as the proposed approach only requires one additional registration per added sample. Compared to the relatively time consuming MDL methods, this allows to update the models right within a clinical workflow. In the experiments with the given MDL software, the registration approach led to models with an increased segmentation accuracy. The generalization quality was also shown in comparisons to gold standard segmentations using a cross validation comparison. Acknowledgements The authors would like to thank Dr. med. R. Zeltner (Medical Clinic IV—Nephrology and Hypertensiology, FAU) for his medical advice and for helping to manually segment the kidneys. The authors are also thankful to HipGraphics for providing the volume rendering software InSpace and Tobias Heimann (German Cancer Research Center Heidelberg) for supplying the software from his segmentation workshop [21] to solve the point correspondence problem based on the MDL. Our special thanks go to Prof. E. Angelopoulou (Chair of Pattern Recognition, FAU) for extensive proof-readings and helpful suggestions. The authors gratefully acknowledge funding of the Erlangen Graduate School in Advanced Optical Technologies (SAOT) by the German National Science Foundation (DFG) in the framework of the excellence initiative. Reference [1]. [1]Pohle R, Tönnies KD. A new approach for model-based adaptive region growing in medical image analysis. In: Skarbek W editors. CAIP, vol. 2124 of lecture notes in computer science. Springer; 2001;p. 238–246. [2]. [2]Kobashi L, Shapiro M. Knowledge-based organ identification from CT images. In: Loew M editors. Proceedings of SPIE on medical imaging. vol. 1652:1992;p. 544–554. [3]. [3]Tsagaan B, Shimizu A, Kobatake H, Kunihisa M, Hanzawa Y. Segmentation of kidney by using a deformable model. In: ICIP: international conference on image processing, Part 3. 2001;p. 1059–1062. [4]. [4]Tsagaan B, Shimizu A, Kobatake H, Miyakawa K. An automated segmentation method of kidney using statistical information. In: MICCAI ‘02: proceedings of the 5th international conference on medical image computing and computer-assisted intervention—Part I. 2002;p. 556–563. [5]. [5]Cootes TF, Taylor CJ, Cooper DH, Graham J. Active shape models—their training and application. Computer Vision and Image Understanding. 1995;61(1):38–59. [6]. [6]Davies R, Twining C, Cootes T, Waterton J, Taylor C. 3D statistical shape models using direct optimisation of description length. In: European conference on computer vision, vol. 3; 2002. p. 3–20. [7]. [7]Davies R, Twining C, Cootes T, Waterton J, Taylor C. A minimum description length approach to statistical shape modelling. IEEE Transactions on Medical Imaging. 2002;21(5):525–537. MEDLINE |
CrossRef
[8]. [8]Heimann T, Wolf I, Williams T, Meinzer H. 3D active shape models using gradient descent optimization of description length. In: Christensen GE, Sonka M editor. IPMI, vol. 3565 of lecture notes in computer science. Springer; 2005;p. 566–577. [9]. [9]Heimann T, Wolf I, Williams T, Meinzer H. Optimal landmark distributions for statistical shape model construction. In: Reinhardt J, Pluim J editor. Proceedings of SPIE on medical imaging. vol. 6144:2006;p. 518–528. [10]. [10]Modersitzki J. Numerical methods for image registration. Oxford: Oxford University Press; 2004;. [11]. [11]Hermosillo G, Chefd’Hotel C, Faugeras O. Variational methods for multimodal image matching. International Journal of Computer Vision. 2002;50(3):329–343. [12]. [12]Clarenz U, Droske M, Henn S, Rumpf M, Witsch K. Computational methods for nonlinear image registration. In: Scherzer O editors. Mathematical method for registration and applications to medical imaging, mathematics in industry. vol. 10:2006;p. 81–101. [13]. [13]Fischer B, Modersitzki J. A unified approach to fast image registration and a new curvature based registration technique. Linear Algebra and its Applications. 2004;380:107–124. [14]. [14]Thodberg HH. Minimum description length shape and appearance models. In: Taylor CJ, Noble JA editor. IPMI, vol. 2732 of lecture notes in computer science. Springer; 2003;p. 51–62. [15]. [15]Floater MS, Hormann K. Surface parameterization: a tutorial and survey. In: Dodgson NA, Floater MS, Sabin MA editor. Advances in multiresolution for geometric modelling. Springer-Verlag; 2005;p. 157–186. [16]. [16]Gu X, Wang Y, Chan T, Thompson P, Yau S. Genus zero surface conformal mapping and its application to brain surface mapping. IEEE Transactions on Medical Imaging. 2004;23(8):949–958. MEDLINE |
CrossRef
[17]. [17]Sethian JA. Level set methods and fast marching methods. Cambridge University Press; 1999;. [18]. [18]William EL, Harvey EC. Marching cubes: a high resolution 3D surface construction algorithm. Computer Graphics. 1987;21(4):163–169. [19]. [19]Cootes TF, Hill A, Taylor CJ, Haslam J. The use of active shape models for locating structures in medical images. Image and Vision Computing. 1994;12(6):355–366. [20]. [20]Cootes T, Taylor C. Statistical models of appearance for computer vision. Tech. Rep., Imaging Science and Biomedical Engineering, University of Manchester; March 2004. [21]. [21]Heimann T, Oguz I, Wolf I, Styner M, Meinzer H. Implementing the automatic generation of 3D statistical shape models with ITK. In: IJ-MICCAI open science workshop; 2006. Martin Spiegel was born in Ingolstadt, Germany, in 1981. Between 2002 and 2007 he studied computer science at the Friedrich-Alexander University Erlangen-Nuremberg (FAU) with main subjects in pattern recognition and medical image processing. Currently, he is working on his Ph.D. thesis in collaboration with Siemens Healthcare Sector, the department of neuroradiology Erlangen and the Chair of Pattern Recognition of Prof. Dr.-Ing. J. Hornegger concerning vessel segmentation and blood flow simulations. His research interests focus on model-based segmentation techniques, mesh generation methods and blood flow simulations in the field of neuroradiology. Dieter A. Hahn was born in Kronach, Germany, in 1979. He received his diploma degree in computer science at the FAU in 2005. Since 2005 he is working as a Ph.D. student for Prof. Dr.-Ing. J. Hornegger at the Chair of Pattern Recognition at the Department of Computer Science of the FAU. His research interests besides general medical image processing include in particular image registration and segmentation. He currently works in collaboration with the Clinics of Nuclear Medicine (FAU) on multi modal image processing techniques for CT, MR and SPECT images. Volker Daum was born in Bavaria, Germany in 1978. In 2004 he received his diploma degree in computer science at the FAU. From 2004 to 2006 he worked at the Fraunhofer Institute IIS in Erlangen-Tennenlohe. Since 2006 he is pursuing his doctoral degree in Computer Science at the Chair of Pattern Recognition (FAU) as a Max-Planck fellow. His research interests in the field of medical imaging include rigid and nonrigid registration, as well as variational techniques for image processing. Jakob Wasza was born in Munich, Germany, in 1982. Since 2003 he studies computer science at the FAU with major interest in pattern recognition and medical image processing. His student thesis was focusing on the application of registration for shape model generation. Joachim Hornegger graduated in Theoretical Computer Science/Mathematics (1992) and received his Ph.D. degree in Applied Computer Science (1996) at the University of Erlangen-Nuremberg (Germany). His Ph.D. thesis was on statistical learning, recognition and pose estimation of 3D objects. Joachim was a visiting scholar and lecturer at Stanford University (Stanford, CA, USA) in the academic year 1997/98, and a visiting professor at Stanford's Radiological Science Lab (RSL) in winter 2007/2008. In 1998 he joined Siemens Medical Solutions Inc. where he was working on 3D angiography. In parallel to his responsibilities in industry he was a lecturer at the Universities of Erlangen (1998–1999), Eichstaett–Ingolstadt (2000), and Mannheim (2000–2003). In 2003 Joachim became Professor of Medical Imaging Processing at the FAU and since 2005 he is a chaired professor heading the Chair of Pattern Recognition. His main research topics are currently pattern recognition methods in medicine and sports. a Erlangen Graduate School in Advanced Optical Technologies (SAOT), Erlangen, Germany b Department of Computer Science, Chair of Pattern Recognition, Friedrich-Alexander University Erlangen-Nuremberg, Germany  Corresponding author. Tel.: +49 9131 8527874; fax: +49 9131 303811.
PII: S0895-6111(08)00100-6 doi:10.1016/j.compmedimag.2008.10.002 © 2008 Elsevier Ltd. All rights reserved. | |
|
FULL TEXT00100-6/journalimage?img=PIIS0895611108001006.gr1.lrg.gif&fig=fig1&kwhquery=&issn=0895-6111&ishighres=false&allhighres=false&free=yes','journalimage');)
00100-6/journalimage?img=PIIS0895611108001006.gr2.lrg.gif&fig=fig2&kwhquery=&issn=0895-6111&ishighres=false&allhighres=false&free=yes','journalimage');)
00100-6/journalimage?img=PIIS0895611108001006.gr3.lrg.gif&fig=fig3&kwhquery=&issn=0895-6111&ishighres=false&allhighres=false&free=yes','journalimage');)
00100-6/journalimage?img=PIIS0895611108001006.gr4.lrg.gif&fig=fig4&kwhquery=&issn=0895-6111&ishighres=false&allhighres=false&free=yes','journalimage');)
00100-6/journalimage?img=PIIS0895611108001006.gr5.lrg.gif&fig=fig5&kwhquery=&issn=0895-6111&ishighres=false&allhighres=false&free=yes','journalimage');)
00100-6/journalimage?img=PIIS0895611108001006.gr6.lrg.gif&fig=fig6&kwhquery=&issn=0895-6111&ishighres=false&allhighres=false&free=yes','journalimage');)
00100-6/journalimage?img=PIIS0895611108001006.gr7.lrg.gif&fig=fig7&kwhquery=&issn=0895-6111&ishighres=false&allhighres=false&free=yes','journalimage');)
00100-6/journalimage?img=PIIS0895611108001006.gr8.lrg.gif&fig=fig8&kwhquery=&issn=0895-6111&ishighres=false&allhighres=false&free=yes','journalimage');)
00100-6/journalimage?img=PIIS0895611108001006.gr9.lrg.gif&fig=fig9&kwhquery=&issn=0895-6111&ishighres=false&allhighres=false&free=yes','journalimage');)
00100-6/journalimage?img=PIIS0895611108001006.gr10.lrg.gif&fig=fig10&kwhquery=&issn=0895-6111&ishighres=false&allhighres=false&free=yes','journalimage');)
00100-6/journalimage?img=PIIS0895611108001006.gr11.lrg.gif&fig=fig11&kwhquery=&issn=0895-6111&ishighres=false&allhighres=false&free=yes','journalimage');)