- Email: [email protected]
Deformable registration of brain tumor images via a statistical model of tumor-induced deformation
Medical Image Analysis 10 (2006) 752–763 www.elsevier.com/locate/media
Deformable registration of brain tumor images via a statistical model of tumor-induced deformation Ashraf Mohamed b
a,b,*
, Evangelia I. Zacharaki b, Dinggang Shen
a,b
, Christos Davatzikos
a,b
a Section of Biomedical Image Analysis, University of Pennsylvania School of Medicine, Philadelphia, PA, USA CISST NSF Engineering Research Center, Department of Computer Science, Johns Hopkins University, Baltimore, MD, USA
Received 30 January 2006; received in revised form 4 May 2006; accepted 8 June 2006 Available online 24 July 2006
Abstract An approach to the deformable registration of three-dimensional brain tumor images to a normal brain atlas is presented. The approach involves the integration of three components: a biomechanical model of tumor mass-effect, a statistical approach to estimate the model’s parameters, and a deformable image registration method. Statistical properties of the sought deformation map from the atlas to the image of a tumor patient are first obtained through tumor mass-effect simulations on normal brain images. This map is decomposed into the sum of two components in orthogonal subspaces, one representing inter-individual differences in brain shape, and the other representing tumor-induced deformation. For a new tumor case, a partial observation of the sought deformation map is obtained via deformable image registration and is decomposed into the aforementioned spaces in order to estimate the mass-effect model parameters. Using this estimate, a simulation of tumor mass-effect is performed on the atlas image in order to generate an image that is similar to tumor patient’s image, thereby facilitating the atlas registration process. Results for a real tumor case and a number of simulated tumor cases indicate significant reduction in the registration error due to the presented approach as compared to the direct use of deformable image registration. 2006 Elsevier B.V. All rights reserved. Keywords: Brain image registration; Brain tumor; Finite element model; Atlas registration; Statistical deformation model; Neurosurgical planning
1. Introduction Deformable registration of normal brain images into a common stereotactic space makes possible the construction of statistical atlases that are based on collective morphological, functional, and pathological information (Davatzikos, 1997). Similar atlases constructed from tumor patients’ images can act as tools for optimal planning of therapeutic and neuro-surgical approaches that deal with tumors by statistically linking functional and structural neuroanatomy to variables such as the tumor size, location, and grade to the surgical or treatment approach and outcomes * Corresponding author. Present address: Siemens Corporate Research Inc., Imaging and Visualization, 755 College Road East, Princeton, NJ 08540, USA. Tel.: +1 609 233 8826; fax: +1 832 355 7923. E-mail address: [email protected] (A. Mohamed).
1361-8415/$ - see front matter 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.media.2006.06.005
(Kyriacou et al., 1999; Mohamed et al., 2001; Dawant et al., 2002; Cuadra et al., 2004). A major hurdle preventing the construction of such brain tumor atlases is the unsuitability of currently available deformable registration methods for adapting a tumor-bearing image to the stereotactic space of a normal neuroanatomy atlas image. This is due to the substantial dissimilarity between the two images resulting from topological differences, tissue death and resorption, the confounding effects of edema and tissue infiltration, and severe deformation in the vicinity of the tumor beyond natural anatomical variability. These image dissimilarities typically result in the failure of currently available image registration methods to produce an accurate registration in the vicinity of the tumor. To account for topological differences between the atlas and the patient’s images Dawant et al. (2002) proposed the
A. Mohamed et al. / Medical Image Analysis 10 (2006) 752–763
introduction of a tumor ‘‘seed’’ in the atlas image and relied on image features to drive the registration. Cuadra et al. (2004) extended this idea by adding a radially symmetric model of tumor growth. The lack of a physically realistic model of tumor-induced deformation, as well as the approximate determination of the seed location results in limited accuracy of these approaches for large tumor cases. In an earlier attempt at solving this deformation registration problem, Kyriacou et al. (1999) used a biomechanical model of the deformation caused by tumors to register images of non-infiltrative tumor patients to anatomical atlases. However, this approach was only implemented in 2D and relied on a computationally expensive regression procedure to solve the inverse problem of estimating the tumor location in the atlas. In order to register brain tumor images to a normal anatomical brain atlas, here we present an approach that involves the integration of three components. The first is a biomechanical three-dimensional 3D model for the softtissue deformation caused by the bulk tumor and peritumor edema. This model is implemented using the finite element (FE) method and it is used to generate a number of examples of deformed brain anatomies due to tumors starting from normal brain images. The second component is a statistical model of the desired deformation map. The statistical model approximates this map via the sum of two components which lie in orthogonal subspaces, and which have different statistical properties. The first of these two components represents normal inter-individual differences in brain shape and the second represents the deformation induced by a tumor (tumor mass-effect). For any particular tumor case that should be registered to the atlas, a partial observation of the sought deformation map is obtained via a deformable image registration method, which is the third component of the presented approach. Based on the constructed statistical model of the deformation, this partial observation is used to estimate the corresponding mass-effect model parameters that correspond to the observed deformation. Finally, the sought deformation map is obtained by applying the mass-effect model to the atlas image and the use of deformable image registration to match it to the subject’s image. The rest of this paper is organized as follows. The details of the proposed approach are presented in Section 2. This includes a description of the mass-effect model, statistical training, and statistical estimation approaches. In Section 3, we demonstrate the overall approach on a real and a number of simulated tumor cases, and we show that the registration error decreases significantly with our approach as compared to the direct use of a readily available image registration method. The paper is concluded with a discussion of the results and future work in Section 4. 2. Methods The proposed approach is explained with the aid of Fig. 1. The subject’s brain BSD includes regions TSD (bulk
753
Fig. 1. Illustration of the deformation maps involved in the proposed approach. uf is the map from the atlas to a subject’s tumor-bearing image. Regions TSD and DSD denote the bulk tumor and edema regions in the subject’s images, and TA, DA are the corresponding regions in the atlas. uc is the mapping from the atlas to the subject’s image before tumor masseffect simulation (BS is not known for non-simulated cases), and ud is that obtained through the simulation of tumor mass-effect. Simulating the tumor mass-effect on the atlas results in ua and a deformed atlas image which can then be registered to the deformed subject’s image through ub.
tumor), and possibly DSD (peri-tumor edema). The main goal of the deformable registration problem is to find the invertible transformation uf : BA n T A ! BSD n T SD which maps points with coordinates XA in the atlas image to points with coordinates XSD in the subject image. Another goal is to identify TA, which corresponds to brain tissue that is no longer present in the subject’s image (died or invaded by tumor). As will be demonstrated later in this paper, the use of a readily available brain image registration approach to obtain uf produces an inaccurate warping in and around tumor area (the area denoted by MA on the atlas image in Fig. 1). This is mainly due to the inability of such image matching approaches to account for topological differences between the two images, the effects of edema, tumor infiltration, signal difference between the two images in an around the tumor region, and tumor mass-effect which causes severe deformation atypical of natural inter-subject variability to which these approaches are tuned. If an accurate model of the deformation induced by the tumor is available, it can be used to simulate this deformation in the atlas and to obtain ua, followed by the application of deformable image registration to get ub, and therefore uf = ub ua. A model of the mass-effect caused by tumor growth is described in Section 2.1. Estimates of region TA as well as the other parameters affecting the model’s behavior, such as the extent of peri-tumor edema and the mass-effect of the bulk tumor, are still needed in order to apply this approach. Here, we solve this inverse estimation problem by exploiting the statistical dependency between uf and the mass-effect model parameters.
754
A. Mohamed et al. / Medical Image Analysis 10 (2006) 752–763
Although an approximation of uf obtained through the direct application of deformable image registration is incorrect in and around the tumor, the pattern of this deformation outside that region can guide the estimation of the tumor model parameters. The overall approach for deformable image registration between an atlas and a brain tumor image is summarized in Fig. 2. In Section 2.2, we explain the collection of the statistics on uf = ud uc through tumor mass-effect simulations on images of normal subjects. The use of these statistics to estimate the tumor mass-effect model parameters is explained in Section 2.3. 2.1. Tumor mass-effect model In Mohamed and Davatzikos (2005), we described a 3D FE model of brain tumor mass-effect. This model is initialized with a 3D normal brain image (free of tumor) and it produces an estimate of the deformation due to the masseffect of a simulated tumor. Here, we briefly review the description of this model by assuming that it is applied to the atlas image. As explained later in this paper, this mass-effect model is also applied to other normal images for purposes of statistical model training. Since tumor growth is not a purely mechanical process, but it involves a host of interacting biological, biochemical and mechanical mechanisms at different scales, it is essential to initialize the mass-effect model simulations with a configuration for the brain from which the target configuration (that deformed by the tumor at the desired stage of tumor growth) is reachable by solving a mechanical boundary value problem (BVP).
With the assumption that the tumor mass-effect is due to the bulk tumor and the peri-tumor edema only, regions TA and DA are defined in the undeformed (normal) atlas image. These correspond to the bulk tumor and peritumor edema regions, respectively, in the deformed atlas at the end of the simulation. The bulk tumor is assumed to be composed of proliferative, quiescent, and necrotic tumor cells (Kansal et al., 2000; Wasserman and Acharya, 1996; Greenspan, 1972). The extent of peri-tumor edema for brain tumors is highly variable and it depends on the type of the involved tumor and its grade (Lamszus, 2004). Given the initial shape of the brain, regions TA and DA, the stresses caused by the tumor, and the amount of swelling due to edema, the deformation map ua can be obtained by solving the mechanical BVP. For real tumor cases, values of the parameters which are responsible for producing the deformation in the patient’s brain are not known, but can be estimated via the statistical approach as explained in Section 2.3. To simulate the mass-effect of tumors starting with normal brain images, approximations of these parameters must be used as explained in the following section. 2.1.1. Specifying TA and DA TA corresponds to brain tissue that is no longer present at the target stage of tumor growth (died or infiltrated by tumor cells), while DA corresponds to brain tissue that is swollen by edema. Although regions TA and DA are highly variable for different tumor cases and are not known in general, here for tractability, we assume that these regions are both spherical and concentric with center ct and radii rt
~ f is an initial rough estimate of the map from Fig. 2. A flowchart summarizing the overall steps involved in the described image registration approach. u e is an estimate of the parameters of tumor mass-effect model. A mass-effect the atlas to the patient’s brain obtained through direct image registration. H simulation on the atlas followed by deformable image registration produces the two deformation maps ua and ub, respectively, which are composed to ^ f , the estimated atlas to tumor image deformation map. produce, u
A. Mohamed et al. / Medical Image Analysis 10 (2006) 752–763
and rd, respectively. ct, rt, and rd will now be treated as parameters of the mass-effect model. It is important to note that the assumptions on the shape of regions TA and DA do not restrict our approach to dealing with spherical tumors only since the final simulated tumor shapes depend on the elastic constraints around the tumor (Kyriacou et al., 1999). As illustrated later in this article in Fig. 3, tumors generated by the described masseffect model need not be spherical. Also, it is worth mentioning that for the deformable registration of real brain tumor cases to brain atlases, estimates of regions TA and DA are later deformed through the deformable image registration component of our approach. 2.1.2. Bulk tumor mass-effect and edema To account for the mass-effect of the bulk tumor, we follow the work of Wasserman and Acharya (1996) and assume that the expansive force of the neoplasm is generated by the proliferative tumor cells only. Accordingly, for model simulations, brain tissue in the region TA is removed, and a constant outward pressure P is applied normal to the boundary of Tr. P is a model parameter that
755
determines the mass-effect exerted by the bulk tumor, and therefore, to a large extent, the final tumor size. Depending on the type of tumor and its aggressiveness, TA may be surrounded with a peri-tumor edema region DA. Peri-tumor edema is usually of the vasogenic type, and it causes swelling of white matter (WM) only, mostly perpendicular to the direction in which the fibers run (Kuroiwa et al., 1994; Nagashima et al., 1994). Here, no knowledge of WM fibers’ orientation is assumed and therefore, an isotropic expansive strain e is applied to WM in DA by using analogy to thermal expansion. A thermal conductivity value of zero for brain tissues prevents this expansion from spreading outside DA. Studies that measured a volume expansion of 200–300% in WM due to edema (Kuroiwa et al., 1994; Nagashima et al., 1994) imply that e 2 [0.26,0.44]. For simulations starting with normal brain scans, a value of e = 0.35 is adopted. 2.1.3. Boundary value problem statement Given the time scale of the tumor growth process, which is at least a few days, the deformation of brain tissues may be modeled as a quasi-static process. Additionally, if body
Fig. 3. Illustration of a tumor mass-effect simulation and the associated displacement maps. Upper row (left to right): atlas image, normal subject’s MRI with an introduced small tumor, and resulting image after simulation of tumor mass-effect. Middle row (left to right): displacement map uc, displacement map ud XS, and displacement map uf. Bottom figure: displacement map ud.
756
A. Mohamed et al. / Medical Image Analysis 10 (2006) 752–763
forces, such as gravity, are ignored, the required deformation mapping ua can be found by solving the static equilibrium equation (Marsden and Hughes, 1983): DivðSÞ ¼ 0
ð1Þ
where S is the first Piola–Kirchhoff tensor, which is related to strain via the material constitutive law (ABAQUS Version 6.4, 2003). FE simulations were run to compare the behaviour of the linear elastic material model (Wasserman and Acharya, 1996; Skrinjar et al., 2002; Clatz et al., 2005) and four hyperelastic material models Kyriacou et al. (1999); Mendis et al. (1995); Miller and Chinzei (2002); Prange et al. (2002) used for brain tissues in the literature. In all simulations, nonlinearities caused by large deformations were taken into account. The stability of the tested material models at strain levels encountered during tumor masseffect simulations, the ability of these simulations to converge, and the error in predicting deformations observed in real tumor cases were the factors used to select a material model to use in our mass-effect model. Based on these experiments, which are described in Mohamed (2005), we adopted the hyperelastic model proposed by Miller and Chinzei (2002) while relaxing the perfect incompressibility assumption. Since the cell doubling time of typical brain tumors such as gliomas has been estimated to be between 1 week and 12 months (Swanson et al., 2000), and since the viscous time constants of normal soft tissues are at most of the order of 1 min (Miller and Chinzei, 2002; Prange et al., 2002), viscous effects of the brain tissues may be ignored. Under these assumptions, the strain energy density function of the brain material becomes (ABAQUS Version 6.4, 2003): W ¼
2l a 1 2 ðk þ ka2 þ ka3 3Þ þ ðJ =J th 1Þ ; a2 1 D1
ð2Þ
where ki ¼ J 1=3 ki , ki, i = 1, 2, 3 are the principal material stretches, J ¼ detð~ F Þ is the volume ratio, ~ F is the deformation gradient, Jth = (1 + eth)3 is the thermal volume ratio, and eth is the thermal strain. The constants l, D1 are related to the Young’s modulus at zero strain Eo, and Poisson’s ratio m by Eo l¼ 2ð1 þ mÞ
and
6ð1 2mÞ D1 ¼ : Eo
ð3Þ
The value a = 4.7 determined in Miller and Chinzei (2002) was adopted here. Based on the arguments and the experiments explained in Mohamed and Davatzikos (2005) and Mohamed (2005), we adopted the value of l = 842 Pa suggested by Miller and Chinzei (2002) and a Poisson’s ratio of m = 0.485. A sliding boundary condition is used over the whole brain surface except for locations, where the falx meets the inner surface of the skull which are assumed pinned (Miga et al., 1999a). The sliding boundary condition approximates the interaction between the brain and the inner surface of the skull by preventing motion in the direc-
tion normal to the brain surface and permitting motion in the tangential direction (Miga et al., 1999b). For simulations in this paper, the brain ventricles are assumed to be void with a zero intra-ventricular pressure (Kyriacou et al., 1999). These boundary conditions, the static equilibrium Eq. (1), and the material constitutive relationship (2) formulate a mechanical BVP which is solved via the FE method using ABAQUS (ABAQUS Version 6.4, 2003). FE meshes are generated using the approach described in Mohamed and Davatzikos (2004). 2.2. Statistical model training Let the parameters of the mass-effect model described in Section 2.1 be collectively referred to by H ” (ct, rt, rd, P). The values of these parameters are not known for a real tumor case. The goal of the work described in the section is to create a statistical model for the deformation uf that will be useful in the estimation of H based on the brain tumor image for a patient. Since uc is not related to the tumor model parameters, we collect statistics for ud and uc separately as follows. First, the deformation maps uci ; i ¼ 1; . . . ; ns between the atlas and MRI images of ns normal subjects are obtained using a deformable image registration approach (Shen and Davatzikos, 2002). Simulations of the masseffect of tumor growth are then conducted for each subject i for values Hj, j = 1,. . ., nm covering a range of the model parameters to produce the deformations ud i;j ; i ¼ 1; . . . ; ns ; j ¼ 1; . . . ; nm . A problem preventing the collection of statistics on ud i;j directly is that the domains of these maps are different for different values of i and j. This precludes the point-to-point comparison of these deformation maps. To overcome this problem, for all tumor model simulations, regions T Aj and DAj are defined in the atlas space based on Hj and mapped to each subject’s space via uci ; i ¼ 1; . . . ; ns . Now, for XA 2 BA n T Aj ; i ¼ 1; . . . ; ns ; j ¼ 1; . . . ; nm , we define ufi;j ðXA Þ ud i;j uci ðXA Þ ¼ ud i;j ðuci ðXA ÞÞ
ð4Þ
ufi;j ðXA Þ ufi;j ðXA Þ XA
ð5Þ
uci ðXA Þ uci ðXA Þ XA ud i;j ðXA Þ ufi;j ðXA Þ uci ðXA Þ
ð6Þ ð7Þ
Eqs. (5)–(7) imply that ufi;j ðXA Þ ¼ uci ðXA Þ þ ud i;j ðXA Þ:
ð8Þ
Therefore, we have decomposed ufi;j , the displacement map from the atlas to a tumor image, into the sum of two displacement maps, uci and ud i;j , which are both defined in the atlas space. The displacement uci is due to differences between the shape of the brain of subject i and the atlas. The displacement ud i;j is due to the tumor mass-effect simulation with tumor model parameters Hj on the brain of subject i. For different i = 1, . . ., ns, but the same
A. Mohamed et al. / Medical Image Analysis 10 (2006) 752–763
j = 1, . . ., nm, the domains of ud i;j are exactly the same. In Fig. 3, we show an example of a tumor mass-effect simulation on a normal brain image and we also show the involved displacement maps. We construct discrete versions of the displacement maps uci and ud i;j by sampling their Cartesian components for all voxels in the atlas in BAnMA to yield vectors uci and ud i;j , respectively, for i = 1, . . ., ns, j = 1, . . ., nm. Assuming that uci ; i ¼ 1; . . . ; ns are independent realizations of a Gaussian random vector, principal component analysis (PCA) is applied to these vectors to yield the mean lc and the matrix Vc whose columns are eigenvectors corresponding to the first mc principal components (mc 6 ns 1). Through this assumption, we have approximated the displacement map from the atlas to a normal brain image via a linear subspace in the space of all possible displacements from the atlas. Next, we construct a statistical model of the displacement due to the tumor mass-effect. We compute the component of ud i;j in the subspace orthogonal to the columns of Vc as u0d i;j ¼ Ud i;j Vc VTc Ud i;j :
ð9Þ
We further assume that, for each j, u0d i;j ; i ¼ 1; . . . ; ns are independent realizations of a Gaussian random vector and we perform PCA on these vectors to yield the mean ld j and the matrices Vd j whose columns are eigenvectors corresponding to the first md j principal components associated with eigenvalues kd j ;l ; l ¼ 1; . . . ; md j ðmd j 6 ns 1Þ. Now, we can approximate the discrete displacement map Uf between the atlas and a subject with a simulated tumor with parameters Hj, j = 1,. . .,nm as follows: Uf lc þ Vc a þ ld j þ Vd j bj :
ð10Þ
Therefore, we can now represent the displacement from the atlas to a image of a tumor patient as the sum of two components with different statistical properties that can be learnt from readily available training data. The first component (lc + Vca) represents normal deviation of the shape of the brain of the patient from that of the atlas. The second component ðld j þ Vd j bj Þ represents tumor mass-effect, and is chosen, by construction, to be orthogonal to the space of the first component. This orthogonality, although not strictly essential, makes the estimation of the tumor mass-effect model parameters for a tumor patient easier as will be explained below. The use of PCA implies that the vectors a, and T bj ¼ ½bj;1 ; . . . ; bj;md j each follows a Gaussian distribution with decorrelated components. The probability density function of bj explicitly stated here as: ! md j X b2j;l 1 exp 0:5 fj ðbj Þ ¼ md kd j ;l Qj pffiffiffiffiffiffiffiffiffiffiffiffiffi l¼1 2pkd j ;l l¼1
for j ¼ 1; . . . ; nm :
ð11Þ
757
2.3. Statistical estimation ~ f (between a Given an approximate deformation map u real tumor patient’s images and the atlas) obtained through the direct use of deformable image registration, the goal b of the methods presented here is to obtain an estimate H of the tumor model parameters. The displacement map ~uf defined in a similar manner to Eq. (5) is also discretized over all the atlas voxels in BAnMA and represented by a ~ f . Owing to the orthogonality of Vd to Vc for all vector U j j, we can compute the component of this displacement that is caused by the tumor by projection as: ~d ¼ U ~ f lc Vc ~a; U
ð12Þ
~f VTc ðU
~ d could be where ~a ¼ lc Þ. The likelihood that U generated with tumor model parameters Hj is defined as: Lj fj ð~bj Þ;
ð13Þ
where ~bj ¼ VT ðU ~ d ld;j Þ d;j
for
j ¼ 1; . . . ; nm :
ð14Þ
High Lj implies that the corresponding Hj is likely according to the observed tumor-induced deformation. One possible approach to estimating H would involve computing the likelihoods for all Hj that were used for training and returning the one with the highest Lj. This estimator would however provide only discrete values of H – one of those used during training. To produce an estimator with continuous values of H within the range used in training, we propose to estimate the tumor mass-effect model parameters as: !, ! nm nm X X b ¼ H Lj Hj Lj : ð15Þ j¼1
j¼1
This estimator favors values of Hj used in training and having high likelihood, while allowing the estimated H to take continuous values. 3. Experiments and results Results of applying the approach described above are reported here for a real tumor image and for eight simulated tumor images. The real tumor image is an MRI of patient with a glioma and a large region of peri-tumor edema. The simulated tumor images are obtained by applying the mass-effect model described in Section 2.1 to MR images of normal subjects. All images used are T1weighted MRI. The atlas image dimensions are 256 · 256 · 198 and a voxel size of 1 · 1 · 1 mm. The used and the simulated tumor images are of dimensions 256 · 256 · 124 and voxel size 0.9375 · 0.9375 · 1.5 mm. The FE tumor mass-effect model simulations are the most computationally intensive step of the presented approach. The simulations were carried out using the FE software ABAQUS (ABAQUS Version 6.4, 2003). The average time needed to perform one simulation was
758
A. Mohamed et al. / Medical Image Analysis 10 (2006) 752–763
approximately 35 min on a 900 MHz processor SGI machine. In order to make the statistical training step tractable, we performed tumor simulations on ns = 20 MRI brain images of normal subjects. For each subject, nm = 64 simulations were performed with two values of each of the six model parameters covering the range expected for the real tumor case. The parameter values were rt 2 {3, 5} mm, rd 2 {20, 27} mm, P 2 {2, 5} kPa and corners of a cube in the atlas for the simulated tumor center locations. Therefore, in total, 1280 simulations were performed for statistical training. For the results reported here, all principal components of the displacement Uc were retained and we used md j ¼ 1, for j = 1, . . ., nm. 3.1. Tumor patient image In Figs. 4 and 5, the results of applying the proposed approach to register the image of the real tumor subject to the atlas are illustrated. With the use of deformable registration to directly register the (normal) atlas image to the patient’s MR image, the warping result is inaccurate in the tumor area. Gray matter from the right cingulate region and adjacent cortical CSF in the atlas were stretched to match the tumor and the surrounding edema in the patient’s image. To apply the statistical estimate approach described in Section 2.3, the region MA, where the registra-
tion is inaccurate was manually outlined. The estimated tumor model parameters were ^ct ¼ ð109; 86; 126Þ, ^rt ¼ 3:9 mm, ^rd ¼ 24 mm and Pb ¼ 3:55 kPa. In order to quantitatively assess the improvement in the registration accuracy due to the proposed approach, 21 landmark points were selected around the tumor area in the patient’s images and corresponding points were identified by an expert in the atlas. The point coordinates were mapped through the resulting deformation map with direct deformable registration, and with the approach described above. The location of one landmark point in the atlas image and the patient’s image before and after registration is shown in Fig. 6. The results for all 21 points are presented in Table 1. The maximum error was reduced by 71% by the use of our approach while the mean error was reduced by 57.6%. 3.2. Simulated tumor images Deformable registration experiments were performed for simulated tumor images generated from MRI scans of normal brain images. The brain tumor mass-effect model described in Section 2.1 was used for these simulations. Values of the mass-effect model parameters used were chosen to be in the range of the parameter values used for training the statistical model. The use of simulated
Fig. 4. Three orthogonal 2D images of the atlas and the real tumor patient before and after deformable registration. Left to right: atlas image, subject’s image, atlas image warped to the patient’s space via direct deformable registration between the two images, and the warped atlas image in the patient’s space via the use of the proposed approach.
A. Mohamed et al. / Medical Image Analysis 10 (2006) 752–763
759
Fig. 5. Three orthogonal 2D images of the atlas and the real tumor patient before and after deformable registration. Some labels associated with the atlas image are warped and superimposed on the patient’s images. Left to right: atlas image with all labels, atlas image with five selected labels near the tumor area, patient’s image, patient’s image with superimposed labels warped from the atlas via direct deformable registration between the two images, and the patient’s image with superimposed labels warped from the atlas via the proposed approach. The five selected labels are: the right middle frontal gyrus (green), the right medial frontal gyrus (dark green), the right superior frontal gyrus (cyan), the right cingulate region (magenta), and the left cingulate region (brown).
Fig. 6. Locations of one manually selected landmark point in three orthogonal 2D images of the atlas and the patient before and after deformable registration. The landmark point is marked with the green cross. Left to right: atlas image, patient’s image with the manually selected landmark point, patient’s image with the warped landmark point from the atlas with the direct use of deformable registration, and the patient’s image with the warped landmark point from the atlas with the proposed approach.
760
A. Mohamed et al. / Medical Image Analysis 10 (2006) 752–763
Table 1 Deformable registration error statistics for landmark points in the real tumor (RT) and simulated tumor (ST) cases
RT no model, mm RT with model, mm ST no model, mm ST with model, mm
Minimum
Mean
Maximum
Standard deviation
1.06 0.47 2.54 0.61
8.70 3.69 6.39 3.90
24.87 7.19 10.91 7.79
6.19 1.83 2.62 2.01
For each case, the errors are provided for the direct deformable image registration to the atlas (no model), and the registration using the approach described in this paper (with model). 21 landmark points were used for RT and 25 were used for ST.
mass-effect images provides two benefits. First, since for simulated cases, both uc and ud are available, the mapping ud uc is also available. Therefore, ud uc is treated as ground truth for uf and is used for evaluating the accuracy of the deformable registration, i.e., ub ua. Second, since values of the mass-effect model parameters are known for the simulated cases, this offers a method for evaluating the accuracy of the statistical estimator of the model parameters presented in Section 2.3. First, we report the results for one simulated tumor case with the parameters ct = (106, 86, 128), rt = 4.5 mm, rd = 21 mm and P = 4.5 kPa. Using the statistical approach described above, the estimated values of these parameters were ^ct ¼ ð109; 85; 128Þ, ^rt ¼ 4:1 mm, ^rd ¼ 23 mm and Pb ¼ 3:6 kPa. The subject’s image and that of the atlas before and after registration are shown in Fig. 7. The warping of a few selected labels associated with the atlas to the subject’s image near the tumor area are shown in Fig. 8. Visually, the results indicate an improve-
ment in the image registration accuracy due to the approach described in this paper. To quantitatively evaluate the registration accuracy in a way similar to that used in the real tumor case, 25 points were selected arbitrarily in the area around the simulated tumor, and their corresponding coordinates (found through ud uc) were computed in the atlas image and treated as ground truth. The errors for the direct deformable registration and that obtained by the proposed approach are presented in Table 1 next to the results of the real tumor case. On average, for the simulated case, the maximum error was reduced by 29% using the proposed approach and the corresponding average error was reduced by 39%. In another experiment, seven simulated tumor cases were generated with mass-effect model parameters ct = (109, 84, 129), rt = 4.3 mm, rd = 25 mm and P = 4 kPa starting with brain images of seven normal humans. The estimated model parameters in each case is provided in Table 2. The
Fig. 7. Three orthogonal 2D images of the atlas and the simulated tumor subject before and after deformable registration. Left to right: atlas image, subject’s image before simulation of tumor mass-effect, subject’s image after simulation of tumor mass-effect, atlas image warped to the subject’s space via direct deformable registration between the two images, and the warped atlas image in the subject’s space via the use of the proposed approach.
A. Mohamed et al. / Medical Image Analysis 10 (2006) 752–763
761
Fig. 8. Three orthogonal 2D images of the atlas and the simulated tumor subject before and after deformable registration. Some labels associated with the atlas image are warped and superimposed on the subject’s images. Left to right: atlas image with all labels, atlas image with five selected labels, subject’s image, subject’s image with superimposed labels warped from the atlas via direct deformable registration between the two images, and the subject’s image with superimposed labels warped from the atlas via the proposed approach. The five selected labels are the same as those shown in Fig. 5 for the real tumor case.
Table 2 Estimated values of the tumor mass-effect model parameters for seven simulated tumor images
Subj1 Subj2 Subj3 Subj4 Subj5 Subj6 Subj7
ct
rt, mm
rd, mm
P, Pa
(109, 85, 127) (109, 86, 127) (109, 85, 127) (109, 85, 127) (109, 85, 127) (109, 86, 125) (109, 83, 126)
4.0 4.0 4.0 4.0 4.0 4.0 4.0
23 23 23 23 23 24 24
3456 3432 3463 3481 3467 3544 3476
The actual values of mass-effect model parameters used to generate these images were ct = (109, 84, 129), rt = 4.3 mm, rd = 25 mm and P = 4 kPa.
deformable registration accuracy was evaluated as follows. A spherical region of radius 2.5 cm was selected in the atlas around the center of the simulated tumor but outside the tumor itself. This region, which is the same region where edema is assumed to be present, includes voxels of the atlas, where the direct deformable registration between the atlas the tumor image is likely to be inaccurate. The maximum, minimum, and root-mean-squares (r.m.s.) error in the registration for all voxels in this region are reported in Table 3. Despite differences between the estimated masseffect model parameters and the true values used for simulations, the deformable registration results indicate that on
Table 3 Minimum, maximum, and root-mean-square errors (mm) in the deformable registration for seven simulated tumor images No model
Subj1 Subj2 Subj3 Subj4 Subj5 Subj6 Subj7 Mean
With model
Minimum
Maximum
r.m.s.
Minimum
Maximum
r.m.s.
0.028 0.117 0.138 0.074 0.063 0.090 0.071 0.083
19.70 18.07 22.08 18.39 22.58 22.52 16.37 19.96
7.54 8.26 6.19 5.55 6.87 6.50 4.97 6.55
0.031 0.049 0.043 0.042 0.098 0.080 0.078 0.060
16.33 16.47 14.18 14.71 18.76 9.91 18.22 15.51
3.62 4.43 3.44 4.34 4.08 3.62 4.47 4.00
Results are reported for direct deformable registration between the atlas and the tumor image (No Model) and the approach described in this paper (With Model). Voxels in the atlas over which the error statistics are computed lie within 25 mm from the center of the simulated tumors.
762
A. Mohamed et al. / Medical Image Analysis 10 (2006) 752–763
average the r.m.s. error in the registration is reduced by 38.9% and the maximum is reduced by 22.3%. In one case (Subj6), the maximum error was reduced by 56% with a corresponding 44% reduction in the r.m.s. error.
4. Discussion and future work In this paper, we described an approach for the deformable registration of a brain atlas to brain tumor images. The approach uses a 3D biomechanical FE model of tumor-induced deformation to introduce and simulate the tumor in the atlas followed by the use of a readily available deformable image registration method. In addition to the proposed overall approach for deformable registration for tumor images, the contributions of this work include introducing a statistical approach for solving the inverse problem of determining the mass-effect model parameters. This statistical approach relies on the decomposition of the desired deformation map (between the atlas and the tumorbearing patient’s image) into the sum of two maps in orthogonal subspaces, defined on the same domain, but having different statistical properties. The first deformation map is from the atlas to another normal brain image. The second, is the deformation map from the normal brain image to one that is deformed by the biomechanical model of tumor mass-effect. The statistical properties of both of these deformation maps are learned via PCA from a number of training samples. Owing to the orthogonality of the two components of the modeled deformation map, an initial rough estimate of this map is projected onto the subspace representing tumor-induced deformation and is used to estimate the tumor mass-effect model parameters. The results of applying the proposed approach to a real tumor image and simulated tumor images indicate a significant reduction in the registration error. These experiments should be regarded as a proof-of-concept study. More validation experiments are needed to asses the viability of the proposed approach for a variety of tumor cases of different grades, types and sizes. In particular, more experiments with real tumor patient images are needed, with an assessment of inter-rater and intra-rater variability in the identification of homologous landmark points between the atlas and the patient images. In addition, the sensitivity of the statistical estimator of the mass-effect model parameters to the number of used principal components and the number of training samples also present important directions for future investigations. Understanding the effects of the introduction of the tumor and edema into the atlas image, the biomechanical mass-effect model simulation, and the statistical estimator on the stability of the overall image registration approach needs further analysis. Visual examination of the image of the atlas after warping to the real and simulated patient’s images reveals no significant differences in the registration with and without the proposed approach away from the tumor, in the ipsilateral and contralateral hemi-
spheres. Quantitative analysis of the stability of the overall approach is subject to future work. One of the difficulties associated with the deformable registration of a normal brain atlas image to a brain tumor patient’s image is the topological differences between the two brains. These differences arise from loss of normal tissue from the patient’s brain (due to its death or invasion by tumor) and its replacement with tumor tissue. To compensate for these differences, the approach described in this paper introduces a region of simple spherical geometry into the atlas image. This region approximates tissue that was lost from the patients brain due to tumor-related effects and that was replaced by tumor tissue. The model approximates peri-tumor edema with another spherical region. The two deformations applied to register the atlas to a brain tumor image do not change the topology of these simple spherical regions introduced into the atlas image. While mechanical factors included in the mass-effect simulations can make the final simulated tumor shape nonspherical, in real tumors, other factors that are not mechanical in nature (e.g., tumor infiltration and edema spread), may also play a role in determining the tumor shape and may cause it to deviate from radial symmetry. More accurate mass-effect simulations may require the definition of bulk tumor and edema regions that are not spherical and that may vary from one simulation to another. To overcome this limitation of the current model, two approaches may be investigated in future work. First, models of tumor infiltration and edema spread (e.g., see Clatz et al., 2005; Nagashima et al., 1994), combined with additional image information (e.g., diffusion tensor images), may be used to obtain patient-specific initializations of the mass-effect model. In the second approach, segmentations of the tumor and edema regions in the patient image, if available, may be mapped back to the atlas via the most current estimate of the deformation map uf and used for a new initialization of the mass-effect model. Thus, the bulk tumor and edema regions will not be necessarily spherical, but will be derived from the shapes of the bulk tumor and edema regions in the patient. This suggests an iterative approach, where the initialization of the masseffect model is hopefully improved as better estimates of the sought deformation map are obtained. Another complication associated with the deformable registration of brain tumor images is the significant signal changes associated with edema in MR images. Edema typically causes hypointensity changes in T1-weighted images, which makes it difficult to discern cortical sulci in the affected brain regions. It is therefore not possible to obtain an accurate deformable registration in these regions based on image matching alone. The coordinates of brain structures masked by edema may however, be estimated from the known structures (outside the edema region) through the statistical estimation approaches (Liu et al., 2004). Although the implementation presented above relied on manual outlining of the region MA (where direct deformable registration is inaccurate), automatic identification of
A. Mohamed et al. / Medical Image Analysis 10 (2006) 752–763
this region can be achieved based on the elastic stretching energy or the matching criterion of many deformable image registration algorithms. The developed biomechanical mass-effect model plays an important role in the overall presented approach for image registration. Improvements of this mass-effect model are expected to improve the registration results. For example, improving the accuracy of the boundary conditions by using a stiffer material for the falx compared to the rest of the brain and allowing the brain tissue to slide against it can be expected to produce more accurate simulations of subcortical tissue deformation. Other improvements to the mass-effect model may involve the use of biphasic consolidation theory material model which has been suggested in the literature for modeling edema expansion (Nagashima et al., 1994) and brain shift (Miga et al., 1999a; Miga et al., 1999b). While this constitutive model offers a direct way to simulate dilatation due to edema, the implementation of this material model in ABAQUS requires the use of second-order tetrahedral FE meshes, which are one order higher than the meshes used in the current mass-effect simulations. To avoid the significantly higher computational cost associated with second-order tetrahedral meshes, we chose to use a single phase material model in this work. The comparison of biphasic models to the adopted single phase model with quadratic tetrahedral meshes is an interesting point for future work. Acknowledgements The authors thank Dr. Nick Fox at the University College London, UK, for providing the tumor patient’s images. We also thank Xiaoying Wu at the Section of Biomedical Image Analysis at the University of Pennsylvania for her help in processing the used data. This work was supported in part by the National Science Foundation under Engineering Research Center Grant EEC9731478, and by the National Institutes of Health Grant R01NS42645. References ABAQUS Version 6.4, 2003. User’s Manual. Hibbitt, Karlsson, and Sorensen Inc., USA. Clatz, O., Sermesant, M., Bondiau, P.-Y., Delingette, H., Warfield, S.K., Malandain, G., Ayache, N., 2005. Realistic simulation of the 3d growth of brain tumors in mr images coupling diffusion with mass effect. IEEE Transactions on Medical Imaging 24 (10), 1334–1346. Cuadra, M.B., Pollo, C., Bardera, A., Cuisenaire, O., Villemure, J.-G., Thiran, J.-P., 2004. Atlas-based segmentation of pathological MR brains using a model of lesion growth. IEEE Transactions on Medical Imaging 23 (10), 1301–1314. Davatzikos, C., 1997. Spatial transformation and registration of brain images using elastically deformable models. Comp. Vis. and Image Understanding, Special Issue on Medical Imaging 66 (2), 207–222. Dawant, B.M., Hartmann, S.L., Pan, S., Gadamsetty, S., 2002. Brain atlas deformation in the presence of small and large space-occupying tumors. Computer Aided Surgery 7, 1–10.
763
Greenspan, H.P., 1972. Models for the growth of a solid tumor by diffusion. Studies in Applied Mathematics LI (4), 317–340. Kansal, A.R., Torquato, S., IV, G.R.H., Chiocca, E.A., Deisboeck, T.S., 2000. Simulated brain tumor growth dynamics using a three-dimensional cellular automaton. Journal of Theoretical Biology 203, 367–382. Kuroiwa, T., Ueki, M., Suemasu, H., Taniguchi, I., Okeda, R., 1994. Biomechanical characteristics of brain edema: the difference between vasogenic-type and cytotoxic-type edema. Acta Neurochir. 60 (Suppl.), 158–161. Kyriacou, S.K., Davatzikos, C., Zinreich, S.J., Bryan, R.N., 1999. Nonlinear elastic registration of brain images with tumor pathology using a biomechanical model. IEEE Transactions on Medical Imaging 18 (7), 580–592. Lamszus, K., 2004. Meningioma pathology, genetics, and biology. Journal of Neuropathology and Experimental Neurology 63 (4), 275–286. Liu, T., Shen, D., Davatzikos, C., 2004. Predictive modeling of anatomic structures using canonical correlation analysis. In: Proceedings of the IEEE International Symposium on Biomedical Imaging (ISBI). Arlington, VA, pp. 1279–1282. Marsden, J.E., Hughes, T.J.R., 1983. Mathematical Foundations of Elasticity. Prentice-Hall Inc., Englewood, NJ. Mendis, K.K., Stalnaker, R.L., Advani, S.H., 1995. A constitutive relationship for large deformation finite element modeling of brain tissue. Journal of Biomechanical Engineering 117, 279–285. Miga, M., Paulsen, K., Kennedy, F.E., Hartov, A., Roberts, D., 1999a. Model-updated image-guided neurosurgery using the finite element method: Incorporation of the falx cerebri. In: Medical image computing and computer assisted intervention 1999Lecture Notes in Computer Science, vol. 1679. Springer-Verlag, Cambridge, UK, pp. 900–909. Miga, M.I., Paulsen, K.D., Lemery, J.M., Eisner, S.D., Hartov, A., Kennedy, F.E., Roberts, D.W., 1999b. Model-updated image guidance: initial clinical experience with gravity induced brain deformation. IEEE Transactions on Medical Imaging 18 (10), 866–874. Miller, K., Chinzei, K., 2002. Mechanical properties of brain tissue in tension. Journal of Biomechanics 35, 483–490. Mohamed, A., 2005. Combining statistical and biomechanical models for estimation of anatomical deformations. Ph.D. thesis, Johns Hopkins University. Mohamed, A., Davatzikos, C., 2004. Finite element mesh generation and remeshing from segmented medical images. In: Proceedings of the 2004 IEEE International Symposium on Biomedical Imaging: From Nano To Macro. Arlington, Virginia, USA, pp. 420–423. Mohamed, A.,Davatzikos, C., 2005. Finite element modeling of brain tumor mass-effect from 3D medical images. In: Medical image computing and computer assisted intervention 2005. Lecture Notes in Computer Science. pp. 400–408. Mohamed, A., Kyriacou, S.K., Davatzikos, C., 2001. A statistical approach for estimating brain tumor induced deformation. In: Proceedings of the IEEE Workshop on Mathematical Models in Biomedical Image Analysis, pp. 52–59. Nagashima, T., Tada, Y., Hamano, S., Skakakura, M., Masaoka, K., Tamaki, N., Matsumoto, S., 1994. The finite element an. Acta Neurochir. 60 (Suppl.), 165–167. Prange, M.T., Margulies, S.S., 2002. Regional, directional, and agedependent properties of the brain undergoing large deformation. Journal of Biomechanical Engineering 124, 244–252. Shen, D., Davatzikos, C., 2002. HAMMER: hierarchical attribute matching mechanism for elastic registration. IEEE Transactions on Medical Imaging 21 (11), 1421–1439. Skrinjar, O., Nabavi, A., Duncan, J., 2002. Model-driven brain shift compensation. Medical Image Analysis 6, 361–373. Swanson, K.R., Alvord Jr., E.C., Murray, J.D., 2000. A quantitative model for differential motility of gliomas in grey and white matter. Cell Proliferation 33 (5), 317–329. Wasserman, R., Acharya, R., 1996. A patient-specific in vivo tumor model. Mathematical Biosciences 136, 111–140.
Deformable registration of brain tumor images via a statistical model of tumor-induced deformation Ashraf Mohamed b
a,b,*
, Evangelia I. Zacharaki b, Dinggang Shen
a,b
, Christos Davatzikos
a,b
a Section of Biomedical Image Analysis, University of Pennsylvania School of Medicine, Philadelphia, PA, USA CISST NSF Engineering Research Center, Department of Computer Science, Johns Hopkins University, Baltimore, MD, USA
Received 30 January 2006; received in revised form 4 May 2006; accepted 8 June 2006 Available online 24 July 2006
Abstract An approach to the deformable registration of three-dimensional brain tumor images to a normal brain atlas is presented. The approach involves the integration of three components: a biomechanical model of tumor mass-effect, a statistical approach to estimate the model’s parameters, and a deformable image registration method. Statistical properties of the sought deformation map from the atlas to the image of a tumor patient are first obtained through tumor mass-effect simulations on normal brain images. This map is decomposed into the sum of two components in orthogonal subspaces, one representing inter-individual differences in brain shape, and the other representing tumor-induced deformation. For a new tumor case, a partial observation of the sought deformation map is obtained via deformable image registration and is decomposed into the aforementioned spaces in order to estimate the mass-effect model parameters. Using this estimate, a simulation of tumor mass-effect is performed on the atlas image in order to generate an image that is similar to tumor patient’s image, thereby facilitating the atlas registration process. Results for a real tumor case and a number of simulated tumor cases indicate significant reduction in the registration error due to the presented approach as compared to the direct use of deformable image registration. 2006 Elsevier B.V. All rights reserved. Keywords: Brain image registration; Brain tumor; Finite element model; Atlas registration; Statistical deformation model; Neurosurgical planning
1. Introduction Deformable registration of normal brain images into a common stereotactic space makes possible the construction of statistical atlases that are based on collective morphological, functional, and pathological information (Davatzikos, 1997). Similar atlases constructed from tumor patients’ images can act as tools for optimal planning of therapeutic and neuro-surgical approaches that deal with tumors by statistically linking functional and structural neuroanatomy to variables such as the tumor size, location, and grade to the surgical or treatment approach and outcomes * Corresponding author. Present address: Siemens Corporate Research Inc., Imaging and Visualization, 755 College Road East, Princeton, NJ 08540, USA. Tel.: +1 609 233 8826; fax: +1 832 355 7923. E-mail address: [email protected] (A. Mohamed).
1361-8415/$ - see front matter 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.media.2006.06.005
(Kyriacou et al., 1999; Mohamed et al., 2001; Dawant et al., 2002; Cuadra et al., 2004). A major hurdle preventing the construction of such brain tumor atlases is the unsuitability of currently available deformable registration methods for adapting a tumor-bearing image to the stereotactic space of a normal neuroanatomy atlas image. This is due to the substantial dissimilarity between the two images resulting from topological differences, tissue death and resorption, the confounding effects of edema and tissue infiltration, and severe deformation in the vicinity of the tumor beyond natural anatomical variability. These image dissimilarities typically result in the failure of currently available image registration methods to produce an accurate registration in the vicinity of the tumor. To account for topological differences between the atlas and the patient’s images Dawant et al. (2002) proposed the
A. Mohamed et al. / Medical Image Analysis 10 (2006) 752–763
introduction of a tumor ‘‘seed’’ in the atlas image and relied on image features to drive the registration. Cuadra et al. (2004) extended this idea by adding a radially symmetric model of tumor growth. The lack of a physically realistic model of tumor-induced deformation, as well as the approximate determination of the seed location results in limited accuracy of these approaches for large tumor cases. In an earlier attempt at solving this deformation registration problem, Kyriacou et al. (1999) used a biomechanical model of the deformation caused by tumors to register images of non-infiltrative tumor patients to anatomical atlases. However, this approach was only implemented in 2D and relied on a computationally expensive regression procedure to solve the inverse problem of estimating the tumor location in the atlas. In order to register brain tumor images to a normal anatomical brain atlas, here we present an approach that involves the integration of three components. The first is a biomechanical three-dimensional 3D model for the softtissue deformation caused by the bulk tumor and peritumor edema. This model is implemented using the finite element (FE) method and it is used to generate a number of examples of deformed brain anatomies due to tumors starting from normal brain images. The second component is a statistical model of the desired deformation map. The statistical model approximates this map via the sum of two components which lie in orthogonal subspaces, and which have different statistical properties. The first of these two components represents normal inter-individual differences in brain shape and the second represents the deformation induced by a tumor (tumor mass-effect). For any particular tumor case that should be registered to the atlas, a partial observation of the sought deformation map is obtained via a deformable image registration method, which is the third component of the presented approach. Based on the constructed statistical model of the deformation, this partial observation is used to estimate the corresponding mass-effect model parameters that correspond to the observed deformation. Finally, the sought deformation map is obtained by applying the mass-effect model to the atlas image and the use of deformable image registration to match it to the subject’s image. The rest of this paper is organized as follows. The details of the proposed approach are presented in Section 2. This includes a description of the mass-effect model, statistical training, and statistical estimation approaches. In Section 3, we demonstrate the overall approach on a real and a number of simulated tumor cases, and we show that the registration error decreases significantly with our approach as compared to the direct use of a readily available image registration method. The paper is concluded with a discussion of the results and future work in Section 4. 2. Methods The proposed approach is explained with the aid of Fig. 1. The subject’s brain BSD includes regions TSD (bulk
753
Fig. 1. Illustration of the deformation maps involved in the proposed approach. uf is the map from the atlas to a subject’s tumor-bearing image. Regions TSD and DSD denote the bulk tumor and edema regions in the subject’s images, and TA, DA are the corresponding regions in the atlas. uc is the mapping from the atlas to the subject’s image before tumor masseffect simulation (BS is not known for non-simulated cases), and ud is that obtained through the simulation of tumor mass-effect. Simulating the tumor mass-effect on the atlas results in ua and a deformed atlas image which can then be registered to the deformed subject’s image through ub.
tumor), and possibly DSD (peri-tumor edema). The main goal of the deformable registration problem is to find the invertible transformation uf : BA n T A ! BSD n T SD which maps points with coordinates XA in the atlas image to points with coordinates XSD in the subject image. Another goal is to identify TA, which corresponds to brain tissue that is no longer present in the subject’s image (died or invaded by tumor). As will be demonstrated later in this paper, the use of a readily available brain image registration approach to obtain uf produces an inaccurate warping in and around tumor area (the area denoted by MA on the atlas image in Fig. 1). This is mainly due to the inability of such image matching approaches to account for topological differences between the two images, the effects of edema, tumor infiltration, signal difference between the two images in an around the tumor region, and tumor mass-effect which causes severe deformation atypical of natural inter-subject variability to which these approaches are tuned. If an accurate model of the deformation induced by the tumor is available, it can be used to simulate this deformation in the atlas and to obtain ua, followed by the application of deformable image registration to get ub, and therefore uf = ub ua. A model of the mass-effect caused by tumor growth is described in Section 2.1. Estimates of region TA as well as the other parameters affecting the model’s behavior, such as the extent of peri-tumor edema and the mass-effect of the bulk tumor, are still needed in order to apply this approach. Here, we solve this inverse estimation problem by exploiting the statistical dependency between uf and the mass-effect model parameters.
754
A. Mohamed et al. / Medical Image Analysis 10 (2006) 752–763
Although an approximation of uf obtained through the direct application of deformable image registration is incorrect in and around the tumor, the pattern of this deformation outside that region can guide the estimation of the tumor model parameters. The overall approach for deformable image registration between an atlas and a brain tumor image is summarized in Fig. 2. In Section 2.2, we explain the collection of the statistics on uf = ud uc through tumor mass-effect simulations on images of normal subjects. The use of these statistics to estimate the tumor mass-effect model parameters is explained in Section 2.3. 2.1. Tumor mass-effect model In Mohamed and Davatzikos (2005), we described a 3D FE model of brain tumor mass-effect. This model is initialized with a 3D normal brain image (free of tumor) and it produces an estimate of the deformation due to the masseffect of a simulated tumor. Here, we briefly review the description of this model by assuming that it is applied to the atlas image. As explained later in this paper, this mass-effect model is also applied to other normal images for purposes of statistical model training. Since tumor growth is not a purely mechanical process, but it involves a host of interacting biological, biochemical and mechanical mechanisms at different scales, it is essential to initialize the mass-effect model simulations with a configuration for the brain from which the target configuration (that deformed by the tumor at the desired stage of tumor growth) is reachable by solving a mechanical boundary value problem (BVP).
With the assumption that the tumor mass-effect is due to the bulk tumor and the peri-tumor edema only, regions TA and DA are defined in the undeformed (normal) atlas image. These correspond to the bulk tumor and peritumor edema regions, respectively, in the deformed atlas at the end of the simulation. The bulk tumor is assumed to be composed of proliferative, quiescent, and necrotic tumor cells (Kansal et al., 2000; Wasserman and Acharya, 1996; Greenspan, 1972). The extent of peri-tumor edema for brain tumors is highly variable and it depends on the type of the involved tumor and its grade (Lamszus, 2004). Given the initial shape of the brain, regions TA and DA, the stresses caused by the tumor, and the amount of swelling due to edema, the deformation map ua can be obtained by solving the mechanical BVP. For real tumor cases, values of the parameters which are responsible for producing the deformation in the patient’s brain are not known, but can be estimated via the statistical approach as explained in Section 2.3. To simulate the mass-effect of tumors starting with normal brain images, approximations of these parameters must be used as explained in the following section. 2.1.1. Specifying TA and DA TA corresponds to brain tissue that is no longer present at the target stage of tumor growth (died or infiltrated by tumor cells), while DA corresponds to brain tissue that is swollen by edema. Although regions TA and DA are highly variable for different tumor cases and are not known in general, here for tractability, we assume that these regions are both spherical and concentric with center ct and radii rt
~ f is an initial rough estimate of the map from Fig. 2. A flowchart summarizing the overall steps involved in the described image registration approach. u e is an estimate of the parameters of tumor mass-effect model. A mass-effect the atlas to the patient’s brain obtained through direct image registration. H simulation on the atlas followed by deformable image registration produces the two deformation maps ua and ub, respectively, which are composed to ^ f , the estimated atlas to tumor image deformation map. produce, u
A. Mohamed et al. / Medical Image Analysis 10 (2006) 752–763
and rd, respectively. ct, rt, and rd will now be treated as parameters of the mass-effect model. It is important to note that the assumptions on the shape of regions TA and DA do not restrict our approach to dealing with spherical tumors only since the final simulated tumor shapes depend on the elastic constraints around the tumor (Kyriacou et al., 1999). As illustrated later in this article in Fig. 3, tumors generated by the described masseffect model need not be spherical. Also, it is worth mentioning that for the deformable registration of real brain tumor cases to brain atlases, estimates of regions TA and DA are later deformed through the deformable image registration component of our approach. 2.1.2. Bulk tumor mass-effect and edema To account for the mass-effect of the bulk tumor, we follow the work of Wasserman and Acharya (1996) and assume that the expansive force of the neoplasm is generated by the proliferative tumor cells only. Accordingly, for model simulations, brain tissue in the region TA is removed, and a constant outward pressure P is applied normal to the boundary of Tr. P is a model parameter that
755
determines the mass-effect exerted by the bulk tumor, and therefore, to a large extent, the final tumor size. Depending on the type of tumor and its aggressiveness, TA may be surrounded with a peri-tumor edema region DA. Peri-tumor edema is usually of the vasogenic type, and it causes swelling of white matter (WM) only, mostly perpendicular to the direction in which the fibers run (Kuroiwa et al., 1994; Nagashima et al., 1994). Here, no knowledge of WM fibers’ orientation is assumed and therefore, an isotropic expansive strain e is applied to WM in DA by using analogy to thermal expansion. A thermal conductivity value of zero for brain tissues prevents this expansion from spreading outside DA. Studies that measured a volume expansion of 200–300% in WM due to edema (Kuroiwa et al., 1994; Nagashima et al., 1994) imply that e 2 [0.26,0.44]. For simulations starting with normal brain scans, a value of e = 0.35 is adopted. 2.1.3. Boundary value problem statement Given the time scale of the tumor growth process, which is at least a few days, the deformation of brain tissues may be modeled as a quasi-static process. Additionally, if body
Fig. 3. Illustration of a tumor mass-effect simulation and the associated displacement maps. Upper row (left to right): atlas image, normal subject’s MRI with an introduced small tumor, and resulting image after simulation of tumor mass-effect. Middle row (left to right): displacement map uc, displacement map ud XS, and displacement map uf. Bottom figure: displacement map ud.
756
A. Mohamed et al. / Medical Image Analysis 10 (2006) 752–763
forces, such as gravity, are ignored, the required deformation mapping ua can be found by solving the static equilibrium equation (Marsden and Hughes, 1983): DivðSÞ ¼ 0
ð1Þ
where S is the first Piola–Kirchhoff tensor, which is related to strain via the material constitutive law (ABAQUS Version 6.4, 2003). FE simulations were run to compare the behaviour of the linear elastic material model (Wasserman and Acharya, 1996; Skrinjar et al., 2002; Clatz et al., 2005) and four hyperelastic material models Kyriacou et al. (1999); Mendis et al. (1995); Miller and Chinzei (2002); Prange et al. (2002) used for brain tissues in the literature. In all simulations, nonlinearities caused by large deformations were taken into account. The stability of the tested material models at strain levels encountered during tumor masseffect simulations, the ability of these simulations to converge, and the error in predicting deformations observed in real tumor cases were the factors used to select a material model to use in our mass-effect model. Based on these experiments, which are described in Mohamed (2005), we adopted the hyperelastic model proposed by Miller and Chinzei (2002) while relaxing the perfect incompressibility assumption. Since the cell doubling time of typical brain tumors such as gliomas has been estimated to be between 1 week and 12 months (Swanson et al., 2000), and since the viscous time constants of normal soft tissues are at most of the order of 1 min (Miller and Chinzei, 2002; Prange et al., 2002), viscous effects of the brain tissues may be ignored. Under these assumptions, the strain energy density function of the brain material becomes (ABAQUS Version 6.4, 2003): W ¼
2l a 1 2 ðk þ ka2 þ ka3 3Þ þ ðJ =J th 1Þ ; a2 1 D1
ð2Þ
where ki ¼ J 1=3 ki , ki, i = 1, 2, 3 are the principal material stretches, J ¼ detð~ F Þ is the volume ratio, ~ F is the deformation gradient, Jth = (1 + eth)3 is the thermal volume ratio, and eth is the thermal strain. The constants l, D1 are related to the Young’s modulus at zero strain Eo, and Poisson’s ratio m by Eo l¼ 2ð1 þ mÞ
and
6ð1 2mÞ D1 ¼ : Eo
ð3Þ
The value a = 4.7 determined in Miller and Chinzei (2002) was adopted here. Based on the arguments and the experiments explained in Mohamed and Davatzikos (2005) and Mohamed (2005), we adopted the value of l = 842 Pa suggested by Miller and Chinzei (2002) and a Poisson’s ratio of m = 0.485. A sliding boundary condition is used over the whole brain surface except for locations, where the falx meets the inner surface of the skull which are assumed pinned (Miga et al., 1999a). The sliding boundary condition approximates the interaction between the brain and the inner surface of the skull by preventing motion in the direc-
tion normal to the brain surface and permitting motion in the tangential direction (Miga et al., 1999b). For simulations in this paper, the brain ventricles are assumed to be void with a zero intra-ventricular pressure (Kyriacou et al., 1999). These boundary conditions, the static equilibrium Eq. (1), and the material constitutive relationship (2) formulate a mechanical BVP which is solved via the FE method using ABAQUS (ABAQUS Version 6.4, 2003). FE meshes are generated using the approach described in Mohamed and Davatzikos (2004). 2.2. Statistical model training Let the parameters of the mass-effect model described in Section 2.1 be collectively referred to by H ” (ct, rt, rd, P). The values of these parameters are not known for a real tumor case. The goal of the work described in the section is to create a statistical model for the deformation uf that will be useful in the estimation of H based on the brain tumor image for a patient. Since uc is not related to the tumor model parameters, we collect statistics for ud and uc separately as follows. First, the deformation maps uci ; i ¼ 1; . . . ; ns between the atlas and MRI images of ns normal subjects are obtained using a deformable image registration approach (Shen and Davatzikos, 2002). Simulations of the masseffect of tumor growth are then conducted for each subject i for values Hj, j = 1,. . ., nm covering a range of the model parameters to produce the deformations ud i;j ; i ¼ 1; . . . ; ns ; j ¼ 1; . . . ; nm . A problem preventing the collection of statistics on ud i;j directly is that the domains of these maps are different for different values of i and j. This precludes the point-to-point comparison of these deformation maps. To overcome this problem, for all tumor model simulations, regions T Aj and DAj are defined in the atlas space based on Hj and mapped to each subject’s space via uci ; i ¼ 1; . . . ; ns . Now, for XA 2 BA n T Aj ; i ¼ 1; . . . ; ns ; j ¼ 1; . . . ; nm , we define ufi;j ðXA Þ ud i;j uci ðXA Þ ¼ ud i;j ðuci ðXA ÞÞ
ð4Þ
ufi;j ðXA Þ ufi;j ðXA Þ XA
ð5Þ
uci ðXA Þ uci ðXA Þ XA ud i;j ðXA Þ ufi;j ðXA Þ uci ðXA Þ
ð6Þ ð7Þ
Eqs. (5)–(7) imply that ufi;j ðXA Þ ¼ uci ðXA Þ þ ud i;j ðXA Þ:
ð8Þ
Therefore, we have decomposed ufi;j , the displacement map from the atlas to a tumor image, into the sum of two displacement maps, uci and ud i;j , which are both defined in the atlas space. The displacement uci is due to differences between the shape of the brain of subject i and the atlas. The displacement ud i;j is due to the tumor mass-effect simulation with tumor model parameters Hj on the brain of subject i. For different i = 1, . . ., ns, but the same
A. Mohamed et al. / Medical Image Analysis 10 (2006) 752–763
j = 1, . . ., nm, the domains of ud i;j are exactly the same. In Fig. 3, we show an example of a tumor mass-effect simulation on a normal brain image and we also show the involved displacement maps. We construct discrete versions of the displacement maps uci and ud i;j by sampling their Cartesian components for all voxels in the atlas in BAnMA to yield vectors uci and ud i;j , respectively, for i = 1, . . ., ns, j = 1, . . ., nm. Assuming that uci ; i ¼ 1; . . . ; ns are independent realizations of a Gaussian random vector, principal component analysis (PCA) is applied to these vectors to yield the mean lc and the matrix Vc whose columns are eigenvectors corresponding to the first mc principal components (mc 6 ns 1). Through this assumption, we have approximated the displacement map from the atlas to a normal brain image via a linear subspace in the space of all possible displacements from the atlas. Next, we construct a statistical model of the displacement due to the tumor mass-effect. We compute the component of ud i;j in the subspace orthogonal to the columns of Vc as u0d i;j ¼ Ud i;j Vc VTc Ud i;j :
ð9Þ
We further assume that, for each j, u0d i;j ; i ¼ 1; . . . ; ns are independent realizations of a Gaussian random vector and we perform PCA on these vectors to yield the mean ld j and the matrices Vd j whose columns are eigenvectors corresponding to the first md j principal components associated with eigenvalues kd j ;l ; l ¼ 1; . . . ; md j ðmd j 6 ns 1Þ. Now, we can approximate the discrete displacement map Uf between the atlas and a subject with a simulated tumor with parameters Hj, j = 1,. . .,nm as follows: Uf lc þ Vc a þ ld j þ Vd j bj :
ð10Þ
Therefore, we can now represent the displacement from the atlas to a image of a tumor patient as the sum of two components with different statistical properties that can be learnt from readily available training data. The first component (lc + Vca) represents normal deviation of the shape of the brain of the patient from that of the atlas. The second component ðld j þ Vd j bj Þ represents tumor mass-effect, and is chosen, by construction, to be orthogonal to the space of the first component. This orthogonality, although not strictly essential, makes the estimation of the tumor mass-effect model parameters for a tumor patient easier as will be explained below. The use of PCA implies that the vectors a, and T bj ¼ ½bj;1 ; . . . ; bj;md j each follows a Gaussian distribution with decorrelated components. The probability density function of bj explicitly stated here as: ! md j X b2j;l 1 exp 0:5 fj ðbj Þ ¼ md kd j ;l Qj pffiffiffiffiffiffiffiffiffiffiffiffiffi l¼1 2pkd j ;l l¼1
for j ¼ 1; . . . ; nm :
ð11Þ
757
2.3. Statistical estimation ~ f (between a Given an approximate deformation map u real tumor patient’s images and the atlas) obtained through the direct use of deformable image registration, the goal b of the methods presented here is to obtain an estimate H of the tumor model parameters. The displacement map ~uf defined in a similar manner to Eq. (5) is also discretized over all the atlas voxels in BAnMA and represented by a ~ f . Owing to the orthogonality of Vd to Vc for all vector U j j, we can compute the component of this displacement that is caused by the tumor by projection as: ~d ¼ U ~ f lc Vc ~a; U
ð12Þ
~f VTc ðU
~ d could be where ~a ¼ lc Þ. The likelihood that U generated with tumor model parameters Hj is defined as: Lj fj ð~bj Þ;
ð13Þ
where ~bj ¼ VT ðU ~ d ld;j Þ d;j
for
j ¼ 1; . . . ; nm :
ð14Þ
High Lj implies that the corresponding Hj is likely according to the observed tumor-induced deformation. One possible approach to estimating H would involve computing the likelihoods for all Hj that were used for training and returning the one with the highest Lj. This estimator would however provide only discrete values of H – one of those used during training. To produce an estimator with continuous values of H within the range used in training, we propose to estimate the tumor mass-effect model parameters as: !, ! nm nm X X b ¼ H Lj Hj Lj : ð15Þ j¼1
j¼1
This estimator favors values of Hj used in training and having high likelihood, while allowing the estimated H to take continuous values. 3. Experiments and results Results of applying the approach described above are reported here for a real tumor image and for eight simulated tumor images. The real tumor image is an MRI of patient with a glioma and a large region of peri-tumor edema. The simulated tumor images are obtained by applying the mass-effect model described in Section 2.1 to MR images of normal subjects. All images used are T1weighted MRI. The atlas image dimensions are 256 · 256 · 198 and a voxel size of 1 · 1 · 1 mm. The used and the simulated tumor images are of dimensions 256 · 256 · 124 and voxel size 0.9375 · 0.9375 · 1.5 mm. The FE tumor mass-effect model simulations are the most computationally intensive step of the presented approach. The simulations were carried out using the FE software ABAQUS (ABAQUS Version 6.4, 2003). The average time needed to perform one simulation was
758
A. Mohamed et al. / Medical Image Analysis 10 (2006) 752–763
approximately 35 min on a 900 MHz processor SGI machine. In order to make the statistical training step tractable, we performed tumor simulations on ns = 20 MRI brain images of normal subjects. For each subject, nm = 64 simulations were performed with two values of each of the six model parameters covering the range expected for the real tumor case. The parameter values were rt 2 {3, 5} mm, rd 2 {20, 27} mm, P 2 {2, 5} kPa and corners of a cube in the atlas for the simulated tumor center locations. Therefore, in total, 1280 simulations were performed for statistical training. For the results reported here, all principal components of the displacement Uc were retained and we used md j ¼ 1, for j = 1, . . ., nm. 3.1. Tumor patient image In Figs. 4 and 5, the results of applying the proposed approach to register the image of the real tumor subject to the atlas are illustrated. With the use of deformable registration to directly register the (normal) atlas image to the patient’s MR image, the warping result is inaccurate in the tumor area. Gray matter from the right cingulate region and adjacent cortical CSF in the atlas were stretched to match the tumor and the surrounding edema in the patient’s image. To apply the statistical estimate approach described in Section 2.3, the region MA, where the registra-
tion is inaccurate was manually outlined. The estimated tumor model parameters were ^ct ¼ ð109; 86; 126Þ, ^rt ¼ 3:9 mm, ^rd ¼ 24 mm and Pb ¼ 3:55 kPa. In order to quantitatively assess the improvement in the registration accuracy due to the proposed approach, 21 landmark points were selected around the tumor area in the patient’s images and corresponding points were identified by an expert in the atlas. The point coordinates were mapped through the resulting deformation map with direct deformable registration, and with the approach described above. The location of one landmark point in the atlas image and the patient’s image before and after registration is shown in Fig. 6. The results for all 21 points are presented in Table 1. The maximum error was reduced by 71% by the use of our approach while the mean error was reduced by 57.6%. 3.2. Simulated tumor images Deformable registration experiments were performed for simulated tumor images generated from MRI scans of normal brain images. The brain tumor mass-effect model described in Section 2.1 was used for these simulations. Values of the mass-effect model parameters used were chosen to be in the range of the parameter values used for training the statistical model. The use of simulated
Fig. 4. Three orthogonal 2D images of the atlas and the real tumor patient before and after deformable registration. Left to right: atlas image, subject’s image, atlas image warped to the patient’s space via direct deformable registration between the two images, and the warped atlas image in the patient’s space via the use of the proposed approach.
A. Mohamed et al. / Medical Image Analysis 10 (2006) 752–763
759
Fig. 5. Three orthogonal 2D images of the atlas and the real tumor patient before and after deformable registration. Some labels associated with the atlas image are warped and superimposed on the patient’s images. Left to right: atlas image with all labels, atlas image with five selected labels near the tumor area, patient’s image, patient’s image with superimposed labels warped from the atlas via direct deformable registration between the two images, and the patient’s image with superimposed labels warped from the atlas via the proposed approach. The five selected labels are: the right middle frontal gyrus (green), the right medial frontal gyrus (dark green), the right superior frontal gyrus (cyan), the right cingulate region (magenta), and the left cingulate region (brown).
Fig. 6. Locations of one manually selected landmark point in three orthogonal 2D images of the atlas and the patient before and after deformable registration. The landmark point is marked with the green cross. Left to right: atlas image, patient’s image with the manually selected landmark point, patient’s image with the warped landmark point from the atlas with the direct use of deformable registration, and the patient’s image with the warped landmark point from the atlas with the proposed approach.
760
A. Mohamed et al. / Medical Image Analysis 10 (2006) 752–763
Table 1 Deformable registration error statistics for landmark points in the real tumor (RT) and simulated tumor (ST) cases
RT no model, mm RT with model, mm ST no model, mm ST with model, mm
Minimum
Mean
Maximum
Standard deviation
1.06 0.47 2.54 0.61
8.70 3.69 6.39 3.90
24.87 7.19 10.91 7.79
6.19 1.83 2.62 2.01
For each case, the errors are provided for the direct deformable image registration to the atlas (no model), and the registration using the approach described in this paper (with model). 21 landmark points were used for RT and 25 were used for ST.
mass-effect images provides two benefits. First, since for simulated cases, both uc and ud are available, the mapping ud uc is also available. Therefore, ud uc is treated as ground truth for uf and is used for evaluating the accuracy of the deformable registration, i.e., ub ua. Second, since values of the mass-effect model parameters are known for the simulated cases, this offers a method for evaluating the accuracy of the statistical estimator of the model parameters presented in Section 2.3. First, we report the results for one simulated tumor case with the parameters ct = (106, 86, 128), rt = 4.5 mm, rd = 21 mm and P = 4.5 kPa. Using the statistical approach described above, the estimated values of these parameters were ^ct ¼ ð109; 85; 128Þ, ^rt ¼ 4:1 mm, ^rd ¼ 23 mm and Pb ¼ 3:6 kPa. The subject’s image and that of the atlas before and after registration are shown in Fig. 7. The warping of a few selected labels associated with the atlas to the subject’s image near the tumor area are shown in Fig. 8. Visually, the results indicate an improve-
ment in the image registration accuracy due to the approach described in this paper. To quantitatively evaluate the registration accuracy in a way similar to that used in the real tumor case, 25 points were selected arbitrarily in the area around the simulated tumor, and their corresponding coordinates (found through ud uc) were computed in the atlas image and treated as ground truth. The errors for the direct deformable registration and that obtained by the proposed approach are presented in Table 1 next to the results of the real tumor case. On average, for the simulated case, the maximum error was reduced by 29% using the proposed approach and the corresponding average error was reduced by 39%. In another experiment, seven simulated tumor cases were generated with mass-effect model parameters ct = (109, 84, 129), rt = 4.3 mm, rd = 25 mm and P = 4 kPa starting with brain images of seven normal humans. The estimated model parameters in each case is provided in Table 2. The
Fig. 7. Three orthogonal 2D images of the atlas and the simulated tumor subject before and after deformable registration. Left to right: atlas image, subject’s image before simulation of tumor mass-effect, subject’s image after simulation of tumor mass-effect, atlas image warped to the subject’s space via direct deformable registration between the two images, and the warped atlas image in the subject’s space via the use of the proposed approach.
A. Mohamed et al. / Medical Image Analysis 10 (2006) 752–763
761
Fig. 8. Three orthogonal 2D images of the atlas and the simulated tumor subject before and after deformable registration. Some labels associated with the atlas image are warped and superimposed on the subject’s images. Left to right: atlas image with all labels, atlas image with five selected labels, subject’s image, subject’s image with superimposed labels warped from the atlas via direct deformable registration between the two images, and the subject’s image with superimposed labels warped from the atlas via the proposed approach. The five selected labels are the same as those shown in Fig. 5 for the real tumor case.
Table 2 Estimated values of the tumor mass-effect model parameters for seven simulated tumor images
Subj1 Subj2 Subj3 Subj4 Subj5 Subj6 Subj7
ct
rt, mm
rd, mm
P, Pa
(109, 85, 127) (109, 86, 127) (109, 85, 127) (109, 85, 127) (109, 85, 127) (109, 86, 125) (109, 83, 126)
4.0 4.0 4.0 4.0 4.0 4.0 4.0
23 23 23 23 23 24 24
3456 3432 3463 3481 3467 3544 3476
The actual values of mass-effect model parameters used to generate these images were ct = (109, 84, 129), rt = 4.3 mm, rd = 25 mm and P = 4 kPa.
deformable registration accuracy was evaluated as follows. A spherical region of radius 2.5 cm was selected in the atlas around the center of the simulated tumor but outside the tumor itself. This region, which is the same region where edema is assumed to be present, includes voxels of the atlas, where the direct deformable registration between the atlas the tumor image is likely to be inaccurate. The maximum, minimum, and root-mean-squares (r.m.s.) error in the registration for all voxels in this region are reported in Table 3. Despite differences between the estimated masseffect model parameters and the true values used for simulations, the deformable registration results indicate that on
Table 3 Minimum, maximum, and root-mean-square errors (mm) in the deformable registration for seven simulated tumor images No model
Subj1 Subj2 Subj3 Subj4 Subj5 Subj6 Subj7 Mean
With model
Minimum
Maximum
r.m.s.
Minimum
Maximum
r.m.s.
0.028 0.117 0.138 0.074 0.063 0.090 0.071 0.083
19.70 18.07 22.08 18.39 22.58 22.52 16.37 19.96
7.54 8.26 6.19 5.55 6.87 6.50 4.97 6.55
0.031 0.049 0.043 0.042 0.098 0.080 0.078 0.060
16.33 16.47 14.18 14.71 18.76 9.91 18.22 15.51
3.62 4.43 3.44 4.34 4.08 3.62 4.47 4.00
Results are reported for direct deformable registration between the atlas and the tumor image (No Model) and the approach described in this paper (With Model). Voxels in the atlas over which the error statistics are computed lie within 25 mm from the center of the simulated tumors.
762
A. Mohamed et al. / Medical Image Analysis 10 (2006) 752–763
average the r.m.s. error in the registration is reduced by 38.9% and the maximum is reduced by 22.3%. In one case (Subj6), the maximum error was reduced by 56% with a corresponding 44% reduction in the r.m.s. error.
4. Discussion and future work In this paper, we described an approach for the deformable registration of a brain atlas to brain tumor images. The approach uses a 3D biomechanical FE model of tumor-induced deformation to introduce and simulate the tumor in the atlas followed by the use of a readily available deformable image registration method. In addition to the proposed overall approach for deformable registration for tumor images, the contributions of this work include introducing a statistical approach for solving the inverse problem of determining the mass-effect model parameters. This statistical approach relies on the decomposition of the desired deformation map (between the atlas and the tumorbearing patient’s image) into the sum of two maps in orthogonal subspaces, defined on the same domain, but having different statistical properties. The first deformation map is from the atlas to another normal brain image. The second, is the deformation map from the normal brain image to one that is deformed by the biomechanical model of tumor mass-effect. The statistical properties of both of these deformation maps are learned via PCA from a number of training samples. Owing to the orthogonality of the two components of the modeled deformation map, an initial rough estimate of this map is projected onto the subspace representing tumor-induced deformation and is used to estimate the tumor mass-effect model parameters. The results of applying the proposed approach to a real tumor image and simulated tumor images indicate a significant reduction in the registration error. These experiments should be regarded as a proof-of-concept study. More validation experiments are needed to asses the viability of the proposed approach for a variety of tumor cases of different grades, types and sizes. In particular, more experiments with real tumor patient images are needed, with an assessment of inter-rater and intra-rater variability in the identification of homologous landmark points between the atlas and the patient images. In addition, the sensitivity of the statistical estimator of the mass-effect model parameters to the number of used principal components and the number of training samples also present important directions for future investigations. Understanding the effects of the introduction of the tumor and edema into the atlas image, the biomechanical mass-effect model simulation, and the statistical estimator on the stability of the overall image registration approach needs further analysis. Visual examination of the image of the atlas after warping to the real and simulated patient’s images reveals no significant differences in the registration with and without the proposed approach away from the tumor, in the ipsilateral and contralateral hemi-
spheres. Quantitative analysis of the stability of the overall approach is subject to future work. One of the difficulties associated with the deformable registration of a normal brain atlas image to a brain tumor patient’s image is the topological differences between the two brains. These differences arise from loss of normal tissue from the patient’s brain (due to its death or invasion by tumor) and its replacement with tumor tissue. To compensate for these differences, the approach described in this paper introduces a region of simple spherical geometry into the atlas image. This region approximates tissue that was lost from the patients brain due to tumor-related effects and that was replaced by tumor tissue. The model approximates peri-tumor edema with another spherical region. The two deformations applied to register the atlas to a brain tumor image do not change the topology of these simple spherical regions introduced into the atlas image. While mechanical factors included in the mass-effect simulations can make the final simulated tumor shape nonspherical, in real tumors, other factors that are not mechanical in nature (e.g., tumor infiltration and edema spread), may also play a role in determining the tumor shape and may cause it to deviate from radial symmetry. More accurate mass-effect simulations may require the definition of bulk tumor and edema regions that are not spherical and that may vary from one simulation to another. To overcome this limitation of the current model, two approaches may be investigated in future work. First, models of tumor infiltration and edema spread (e.g., see Clatz et al., 2005; Nagashima et al., 1994), combined with additional image information (e.g., diffusion tensor images), may be used to obtain patient-specific initializations of the mass-effect model. In the second approach, segmentations of the tumor and edema regions in the patient image, if available, may be mapped back to the atlas via the most current estimate of the deformation map uf and used for a new initialization of the mass-effect model. Thus, the bulk tumor and edema regions will not be necessarily spherical, but will be derived from the shapes of the bulk tumor and edema regions in the patient. This suggests an iterative approach, where the initialization of the masseffect model is hopefully improved as better estimates of the sought deformation map are obtained. Another complication associated with the deformable registration of brain tumor images is the significant signal changes associated with edema in MR images. Edema typically causes hypointensity changes in T1-weighted images, which makes it difficult to discern cortical sulci in the affected brain regions. It is therefore not possible to obtain an accurate deformable registration in these regions based on image matching alone. The coordinates of brain structures masked by edema may however, be estimated from the known structures (outside the edema region) through the statistical estimation approaches (Liu et al., 2004). Although the implementation presented above relied on manual outlining of the region MA (where direct deformable registration is inaccurate), automatic identification of
A. Mohamed et al. / Medical Image Analysis 10 (2006) 752–763
this region can be achieved based on the elastic stretching energy or the matching criterion of many deformable image registration algorithms. The developed biomechanical mass-effect model plays an important role in the overall presented approach for image registration. Improvements of this mass-effect model are expected to improve the registration results. For example, improving the accuracy of the boundary conditions by using a stiffer material for the falx compared to the rest of the brain and allowing the brain tissue to slide against it can be expected to produce more accurate simulations of subcortical tissue deformation. Other improvements to the mass-effect model may involve the use of biphasic consolidation theory material model which has been suggested in the literature for modeling edema expansion (Nagashima et al., 1994) and brain shift (Miga et al., 1999a; Miga et al., 1999b). While this constitutive model offers a direct way to simulate dilatation due to edema, the implementation of this material model in ABAQUS requires the use of second-order tetrahedral FE meshes, which are one order higher than the meshes used in the current mass-effect simulations. To avoid the significantly higher computational cost associated with second-order tetrahedral meshes, we chose to use a single phase material model in this work. The comparison of biphasic models to the adopted single phase model with quadratic tetrahedral meshes is an interesting point for future work. Acknowledgements The authors thank Dr. Nick Fox at the University College London, UK, for providing the tumor patient’s images. We also thank Xiaoying Wu at the Section of Biomedical Image Analysis at the University of Pennsylvania for her help in processing the used data. This work was supported in part by the National Science Foundation under Engineering Research Center Grant EEC9731478, and by the National Institutes of Health Grant R01NS42645. References ABAQUS Version 6.4, 2003. User’s Manual. Hibbitt, Karlsson, and Sorensen Inc., USA. Clatz, O., Sermesant, M., Bondiau, P.-Y., Delingette, H., Warfield, S.K., Malandain, G., Ayache, N., 2005. Realistic simulation of the 3d growth of brain tumors in mr images coupling diffusion with mass effect. IEEE Transactions on Medical Imaging 24 (10), 1334–1346. Cuadra, M.B., Pollo, C., Bardera, A., Cuisenaire, O., Villemure, J.-G., Thiran, J.-P., 2004. Atlas-based segmentation of pathological MR brains using a model of lesion growth. IEEE Transactions on Medical Imaging 23 (10), 1301–1314. Davatzikos, C., 1997. Spatial transformation and registration of brain images using elastically deformable models. Comp. Vis. and Image Understanding, Special Issue on Medical Imaging 66 (2), 207–222. Dawant, B.M., Hartmann, S.L., Pan, S., Gadamsetty, S., 2002. Brain atlas deformation in the presence of small and large space-occupying tumors. Computer Aided Surgery 7, 1–10.
763
Greenspan, H.P., 1972. Models for the growth of a solid tumor by diffusion. Studies in Applied Mathematics LI (4), 317–340. Kansal, A.R., Torquato, S., IV, G.R.H., Chiocca, E.A., Deisboeck, T.S., 2000. Simulated brain tumor growth dynamics using a three-dimensional cellular automaton. Journal of Theoretical Biology 203, 367–382. Kuroiwa, T., Ueki, M., Suemasu, H., Taniguchi, I., Okeda, R., 1994. Biomechanical characteristics of brain edema: the difference between vasogenic-type and cytotoxic-type edema. Acta Neurochir. 60 (Suppl.), 158–161. Kyriacou, S.K., Davatzikos, C., Zinreich, S.J., Bryan, R.N., 1999. Nonlinear elastic registration of brain images with tumor pathology using a biomechanical model. IEEE Transactions on Medical Imaging 18 (7), 580–592. Lamszus, K., 2004. Meningioma pathology, genetics, and biology. Journal of Neuropathology and Experimental Neurology 63 (4), 275–286. Liu, T., Shen, D., Davatzikos, C., 2004. Predictive modeling of anatomic structures using canonical correlation analysis. In: Proceedings of the IEEE International Symposium on Biomedical Imaging (ISBI). Arlington, VA, pp. 1279–1282. Marsden, J.E., Hughes, T.J.R., 1983. Mathematical Foundations of Elasticity. Prentice-Hall Inc., Englewood, NJ. Mendis, K.K., Stalnaker, R.L., Advani, S.H., 1995. A constitutive relationship for large deformation finite element modeling of brain tissue. Journal of Biomechanical Engineering 117, 279–285. Miga, M., Paulsen, K., Kennedy, F.E., Hartov, A., Roberts, D., 1999a. Model-updated image-guided neurosurgery using the finite element method: Incorporation of the falx cerebri. In: Medical image computing and computer assisted intervention 1999Lecture Notes in Computer Science, vol. 1679. Springer-Verlag, Cambridge, UK, pp. 900–909. Miga, M.I., Paulsen, K.D., Lemery, J.M., Eisner, S.D., Hartov, A., Kennedy, F.E., Roberts, D.W., 1999b. Model-updated image guidance: initial clinical experience with gravity induced brain deformation. IEEE Transactions on Medical Imaging 18 (10), 866–874. Miller, K., Chinzei, K., 2002. Mechanical properties of brain tissue in tension. Journal of Biomechanics 35, 483–490. Mohamed, A., 2005. Combining statistical and biomechanical models for estimation of anatomical deformations. Ph.D. thesis, Johns Hopkins University. Mohamed, A., Davatzikos, C., 2004. Finite element mesh generation and remeshing from segmented medical images. In: Proceedings of the 2004 IEEE International Symposium on Biomedical Imaging: From Nano To Macro. Arlington, Virginia, USA, pp. 420–423. Mohamed, A.,Davatzikos, C., 2005. Finite element modeling of brain tumor mass-effect from 3D medical images. In: Medical image computing and computer assisted intervention 2005. Lecture Notes in Computer Science. pp. 400–408. Mohamed, A., Kyriacou, S.K., Davatzikos, C., 2001. A statistical approach for estimating brain tumor induced deformation. In: Proceedings of the IEEE Workshop on Mathematical Models in Biomedical Image Analysis, pp. 52–59. Nagashima, T., Tada, Y., Hamano, S., Skakakura, M., Masaoka, K., Tamaki, N., Matsumoto, S., 1994. The finite element an. Acta Neurochir. 60 (Suppl.), 165–167. Prange, M.T., Margulies, S.S., 2002. Regional, directional, and agedependent properties of the brain undergoing large deformation. Journal of Biomechanical Engineering 124, 244–252. Shen, D., Davatzikos, C., 2002. HAMMER: hierarchical attribute matching mechanism for elastic registration. IEEE Transactions on Medical Imaging 21 (11), 1421–1439. Skrinjar, O., Nabavi, A., Duncan, J., 2002. Model-driven brain shift compensation. Medical Image Analysis 6, 361–373. Swanson, K.R., Alvord Jr., E.C., Murray, J.D., 2000. A quantitative model for differential motility of gliomas in grey and white matter. Cell Proliferation 33 (5), 317–329. Wasserman, R., Acharya, R., 1996. A patient-specific in vivo tumor model. Mathematical Biosciences 136, 111–140.