|Home | About | Journals | Submit | Contact Us | Français|
Although a variety of diffeomorphic deformable registration methods exist in the literature, application of these methods in the presence of space-occupying lesions is not straightforward. The motivation of this work is spatial normalization of MR images from patients with brain tumors in a common stereotaxic space, aiming to pool data from different patients into a common space in order to perform group analyses. Additionally, transfer of structural and functional information from neuroanatomical brain atlases into the individual patient's space can be achieved via the inverse mapping, for the purpose of segmenting brains and facilitating surgical or radiotherapy treatment planning. A method that estimates the brain tissue loss and replacement by tumor is applied for achieving equivalent image content between an atlas and a patient's scan, based on a biomechanical model of tumor growth. Automated estimation of the parameters modeling brain tissue loss and displacement is performed via optimization of an objective function reflecting feature-based similarity and elastic stretching energy, which is optimized in parallel via APPSPACK (Asynchronous Parallel Pattern Search). The results of the method, applied to 21 brain tumor patients, indicate that the registration accuracy is relatively high in areas around the tumor, as well as in the healthy portion of the brain. Also, the calculated deformation in the vicinity of the tumor is shown to correlate highly with expert-defined visual scores indicating the tumor mass effect, thereby potentially leading to an objective approach to quantification of mass effect, which is commonly used in diagnosis.
Population-based statistical image analysis methods have been utilized in a variety of studies of normal brain development and aging, as well as of brain diseases such as Alzheimer's or mild cognitive impairment , but they haven't been applied in studies of brain cancer. For example, studying the tumor origin and location relative to brain structures, especially to white matter pathways, for different types of tumors could potentially have predictive value in terms of tumor progression . Moreover studying imaging patterns in brain tumor patients can potentially help identify tissue that has been significantly infiltrated and should be treated more aggressively. Finally, evaluating the relationship between spatial distributions of radiation dose and clinical progression/outcome can potentially offer insights into better optimizing radiation treatment. These potential studies require the integration of a variety of patient data, such as conventional MRI, perfusion, and DTI of a large number of patients, into the same space, using deformable registration methods.
Deformable registration can also be useful in assisting neurosurgical treatment planning. Specifically, atlases with segmented structures of interest or with integrated information about anatomical and functional variability can be mapped into the patient's image space, and be further utilized in determining treatment plans that minimize the risk for significant functional impairment.
While the problem of co-registering brain images of healthy subjects has been addressed in the literature, the spatial normalization of images affected by tumor pathology is still a very challenging problem that has motivated our work. The application of registration methods designed to register generally normal neuroanatomies  can lead to poor registration when applied to brain images with tumors, due to large deformations and lack of clear definition of anatomical detail in a patient's images owing to edema and tumor infiltration. Specifically, in the images with tumor, the fundamental assumption of topological equivalence between the atlas and the patient's image, which is almost ubiquitous in deformable registration methods, is violated due to (i) the anatomical changes caused by tissue death and tumor emergence and (ii) the large distortions caused by the mass effect of a growing tumor, which are not in agreement with the usual assumption of smoothness of the deformation fields.
In order to address the first violation point, most methods first create topologically equivalent images by either removing the tumor from the patient , or by placing a small tumor seed in the atlas. Other methods just ignore the regions within and around the tumor during image matching since they are regarded as non-informative or unreliable  or use feature points that exist in each image to establish the correspondence relationship . As pointed in , seeding the atlas with a tumor mass is an essential step, since it ensures the continuity of the transformation in the tumor area and preserves from any irregularities that could appear in this region due to the lack of equivalent image content between atlas and patient images. The atlas seeding represents the biological phenomenon of brain tissue loss and replacement by tumor and renders the registration process non-diffeomorphic. Although modeling of brain tissue loss hasn't been emphasized before, it is an important step towards ensuring high accuracy in very small distances around the tumor.
Regarding the second violation point, some methods simulate the tumor-induced deformation, in order to resolve the geometric discrepancies from the physiologic process of tumor growth prior to registration. They either use advanced tumor simulation models of mass-effect and invasion without further accounting for the inter-subject differences , or simplified radial growth models  refined by a non-rigid deformation based on optical flow . Other methods, in order to allow large deformations around the tumor during registration, control locally the amount of regularization  instead of simulating tumor growth. A more analytical survey and classification of those methods according to the type of registration and main limitations is presented by Cuadra et al. .
It is reasonable to expect that incorporating a model of deformation induced by the tumor is desirable and likely to lead to better registration accuracy. In particular, biomechanical models of tumor-induced deformation based on elasticity properties of brain tissue and structures (e.g. some structures, like the falx, are more rigid than others) can guide the registration more successfully compared to oversimplified growth models, such as radial expansion, especially due to the fact that brain tumor images often lack distinct features around the tumor which would guide the registration process with success. Like all models, tumor growth models that simulate tissue loss and displacement utilize a set of parameters. Methods that use the minimum set of parameters (e.g. only the location of a single voxel seed)  simulate pure displacement of the brain structures with zero tissue loss and can therefore be applied for extracerebral lesion (such as meningioma) growth, but are not appropriate for gliomas and brain metastases. The simulation of brain tissue loss requires a larger number of parameters in order to characterize the size and shape of the seed. Also, the number of parameters increases as more advanced biophysical models, that reflect the effects of peritumoral edema and tumor infiltration, are incorporated . In this study we apply a modeling framework that simulates tumor emergence and tumor growth, and also simplistically differentiates between tumor mass effect and tumor infiltration.
In particular, the study presented in this paper is based on the ORBIT framework  and includes (i) estimation of tumor model parameters (for tumor emergence and growth), (ii) simulation of tumor-induced deformation and (iii) calculation of a dense deformation field that maps the (deformed) atlas with simulated tumor to the patient's image. The registration component is based on the assumption that there is equivalent image content between the atlas with simulated tumor and the patient's image, and the deformation between those images is smooth, similar to normal-to-normal image registration. Previous work of our group in this area focused on the development of statistical models for simulating tumor growth by training PCA models across subjects  or within the same subject . The statistical approach was chosen to reduce the high computational cost of the finite element based biomechanical models for tumor growth simulation  leaving the burden of simulations to off-line training. Statistical models, however, are not very accurate and also are limited by the parameters used during training. For example, training a model for irregularly shaped seeds would require an inhibitive large number of training cases. Recently, Hogea et al. proposed in  a biomechanical model developed in an Eulerian formulation and solved using regular grids, which is significantly faster than common finite element models. Thus, here we employ this model as constraints for an objective function in a model-based registration framework that attempts to maximize the similarity between atlas and patient's images. Also, in comparison to , here we focus on increasing the speed of the estimation of the tumor model parameters by optimizing the objective function with the parallel optimization method, APPSPACK (Asynchronous Parallel Pattern Search).
Furthermore, we show a clinical application of this registration framework, namely how to quantitatively characterize the tumor mass effect. The tumor mass effect has been used as a descriptor for classifying gliomas according to their clinical grade  or as an independent predictor of survival  and is therefore an important factor in the characterization of brain tumors. In this study, we investigate how well the estimated parameters (tumor seed and deformation field) help predicting the mass effect.
The basic components of the proposed method, i.e., the tumor emergence and growth simulation, the registration method for images of tumor pathology, and the estimation criterion for the tumor model parameters, are described in Section 2. Results of the method are presented in Section 3, and involve the registration assessment of 21 brain tumor cases (gliomas and metastases), prediction of the tumor origin, and quantification of the mass effect induced by each tumor. The mass effect is estimated based on the calculated deformation fields and is compared against expert human readings.
Fig. 1 illustrates the whole framework for registration of a normal atlas with a brain image of tumor pathology (patient's image), which involves three components: (i) simulation of tumor growth in the atlas image for a set of tumor model parameters, (ii) deformable registration between the atlas with simulated tumor and the patient's image, and (iii) registration assessment, i.e. evaluation of the optimality criterion. By iterating between the three components, the best tumor model parameters are determined as the ones optimizing the objective function (optimality criterion). Brain tissue loss and tumor emergence is simulated by seeding the atlas. The approach we use to simulate the tumor growth is based on a biomechanical model employing incremental linear elasticity . The atlas with simulated tumor is subsequently registered with the patient's image using a deformable registration method, as in the ORBIT registration framework  which follows a deformation strategy that is robust to the confounding factors caused by the presence of tumor.
Next, the overall framework is described more analytically. The first task is to remove skull from the brain  and segment the atlas and individual MR images into white matter (WM), gray matter (GM) and cerebrospinal fluid (CSF), e.g. using FAST (FMRIB's Automated Segmentation Tool) . After tissue segmentation, different labels are assigned to ventricular CSF and cortical CSF by using a modified version of HAMMER . The tumor is manually delineated by an expert in the patient's original image. We should note that we didn't use a separate label for the segmentation of edema for the results presented in this paper, although the algorithm provides this option. Subsequently, the segmented images are co-registered as illustrated in Fig. 2. The patient's image is first registered globally with the tumor-free template (indicated as normal atlas A0 in Fig. 2) by applying an affine transformation . The transformed image is denoted as affine transformed subject in Fig. 2. The tumor model parameters, θ, are estimated by optimizing a function that reflects feature-based similarity and elastic stretching energy. The optimization is implemented in a coarse to fine resolution scheme. For each resolution level, the optimization is performed in parallel via APPSPACK . An initial estimate of the deformation field at each level is obtained by elastically registering the affine transformed subject to the normal atlas in the whole image domain initiated from the upsampled deformation field of the previous resolution level. Since the deformation field displays almost negligible changes in the regions far away from the tumor during the iterative process of optimizing θ, the optimization is performed only in a subdomain of the subject space in order to considerably speed up the implementation.
The three components of this framework (tumor growth simulation, deformable registration and estimation of model parameters) are described with more details in sections 2.1, 2.2 and 2.3, respectively.
Since brain tumor images often exhibit large tumors, the biomechanical simulator needs to be robust to large deformations and also computationally efficient, particularly for registration purposes. In  we have evaluated the impact that two biomechanical simulators have on the accuracy of deformable registration in order to determine whether potential gain in accuracy of the computationally more expensive simulator warrants the additional computational load it imposes. The first simulator , used in , is based on a finite element model of nonlinear elasticity  and unstructured meshes using the commercial software package ABAQUS. The main draw-backs of this simulator are that: (i) unstructured meshes deteriorate significantly in the presence of large deformations induced by a growing brain tumor, thus frequent remeshing may be needed , and (ii) the method is computationally slow, since construction of efficient solvers for the resulting algebraic system of equations is difficult. The second approach  was proposed to bypass these inherent difficulties associated with the previous simulator. An incremental pressure, linear elasticity, model was developed in an Eulerian formulation and solved using regular grids. This approach circumvents the need for mesh generation and remeshing. Thus, large deformations can be captured and efficient solvers can be employed. The comparison of the two simulators  on a limited number of subjects showed that, although the simulator in  was more accurate in simulating tumor growth, after registration (which furthermore improves the structure's displacement) the gain in accuracy was insignificant compared to the additional computational load imposed by it. Therefore, although for modeling tumor growth in a single subject over time we would prefer a more accurate nonlinear model, such as the simulator in , for atlas registration methods, such as in this study, which require iterative tumor simulations for all possible parameters, we prefer the computationally inexpensive simulation framework in  that we shall refer to as Piecewise Linear Eulerian (PLE).
More specifically, in the PLE simulator, the brain is approximated as an inhomogeneous isotropic linear elastic medium, with different material properties in white matter, gray matter and ventricles. The ventricles are treated as a soft compressible elastic material and zero displacements are imposed at the skull. The target domain (brain) is embedded in a larger computational cubic domain (box), with material properties and distributed forces chosen so that the imposed boundary conditions on the true boundary (here consisting of the brain surface and the tumor boundary, respectively) are approximated. An Eulerian formulation is employed to capture large deformations, with a level-set based approach for evolving fronts. The problem is solved using a regular grid discretization with a fast matrix-free multigrid solver for the resulting algebraic system of equations. The methodology is described in detail in . The following material properties were used for white matter, gray matter and ventricles, in order to produce the results of this paper: Young's modulus (stiffness) Ewhite = Egray = 2100 Pa, Eventricles = 500 Pa, and Poisson's ratio (compressibility) νwhite = νgray = 0.45, νventricles = 0.1.
For the simulations in this paper, we assume that the parameters pertaining to the biomechanical model, such as material properties of brain and ventricles, are similar across patients and we don't include them in the optimization process. The remaining parameters, that need to be estimated, relate to the patient-specific tumor characteristics. Specifically, the tumor seed is created by first uniformly eroding the tumor in the patient's image and then (optionally) creating the convex hull of it . A more complex model could potentially use additional parameters to describe the shape of the tumor seed, at the expense of having to estimate these parameters along with the rest. The size of the eroded tumor and the location in the atlas are the 4 parameters pertaining to the tumor seed. The tumor seed is placed in the atlas and the tumor-induced deformation is calculated by the act of outward forces normal to the tumor boundary. In order to simplistically distinguish between tumor expansion (causing brain tissue displacement) and tumor infiltration (not causing displacement) we introduce a 5th parameter, the mass effect factor. This factor is a parameter that allows the simulation of tumor growth to terminate before the tumor in the atlas has reached the size of the tumor in the subject, motivated by the fact that part of the tumor might be infiltrating and therefore causing minimal displacement. In other words, the mass effect factor reflects the amount of deformation that is applied by a growing tumor to nearby brain tissue, and does not necessarily represent the total amount of cancerous tissue in the brain. Since for many tumor types it is very difficult to distinguish the boundary between tumor bulk and tumor infiltration/edema, this parameter is estimated from the observed pattern of deformation.
The deformable registration method is described in  and also summarized here for completeness. The elastic deformation field is calculated in a multiresolution scheme according to the hierarchical approximation of an energy function, which consists of the similarity matching criterion defined in the template space, a constraint on the inverse matching, and smoothness constraints on the displacement field, following the general framework of the HAMMER algorithm . The similarity criterion is designed based on the similarity of attribute vectors, which are defined for each voxel in the image in order to capture the anatomical context (including healthy and malignant tissue) around it. Specifically, the attribute vector, a(x) = [a1 (x) a2 (x) a3 (x) a4 (x) a5 (x)], reflects edge type (a1), tissue type (a2), and geometric moment invariants from all tissue types, respectively. a1 and a2 are scalars taking discrete labels, whereas a3 is a 1×K vector comprising the geometric moment invariants of each tissue and is used to capture shape information, as described more analytically in . In this application, only zero-order regular moments are used. The number of tissue types depends on the segmentation method applied to labeling brain tissue. In this study we used K = 5 tissue types (white matter, gray matter, ventricular CSF, cortical CSF, and tumor). Besides the attributes that capture brain structure information, the attribute vector captures also the geometric location relative to the brain tumor. Specifically, the signed distance from the tumor boundary (a4) and the angular location with respect to the tumor center (a5) are also used as attributes. All attributes are normalized in the range [0, 1] by linear scaling.
The elastic deformation field that spatially warps the template to the patient's image is calculated by maximizing a similarity criterion reflecting the distance of attributes, constrained by a smoothness term. The similarity of two voxels x and y is defined as the weighted summation of two terms: (i) a similarity criterion matching the brain structures, SimB, which is based on the distance of attributes a1, a2, a3 and (ii) a similarity criterion matching the tumor geometry, SimT, which measures how well the tumors are matched based on the attributes a4 and a5. The two similarity terms are non-linear functions of the distance of attributes expressed by the L1-norm, di (x, y) = |ai (x) − ai (y)|, i = 1,‥,5, as shown below:
The constant c1 and c2 determine how much the attribute vector dissimilarity is penalized. We selected c1 = c2 = 10, which cause a rapid decrease in the similarity value with increasing distance of attributes a4 or a5. The combined similarity function is given below:
where w(x, y) is a weighting factor which decreases with the distance of x and y from each tumor respectively:
If at least one of the two images is normal (without tumor), the distance from tumor boundary becomes infinite, w becomes zero (for the region outside the tumor) and the similarity criterion matches only the brain structures. The use of spatially adapted weights ensures that the identification of corresponding points is driven mainly by one of the two matching criteria, whereas the spatially smooth decrease of w makes the total similarity, Sim, smooth. The constant c3 determines the distance in which the two similarity criteria, SimB and SimT, get equal weight (w(x, y)=0.5). We chose this distance to be 8 mm from the tumor boundary and therefore set c3 = 4 (we scaled it linearly in the range [0, 1] using the same normalization factor as used for a4).
In order to determine correspondences, distinctive points (landmarks) are selected by the algorithm on the tumor boundaries and on the healthy structures following the hierarchical mechanism described in . However, if two images are not topologically equivalent, there is no diffeomorphic deformation field that can map the one to the other. Hence, if the tumor seed corresponds to the tissue that actually died, then a diffeomorphic field exists that can warp the one image to the other, which can be calculated by maximizing the proposed similarity criterion, whereas if the seed placed in the atlas is not correct, such a field does not exist. In the latter case, maximization of SimT (matching tumor boundaries) will cause minimization of SimB, and therefore the total similarity can never achieve a high value. If the total similarity drops below a threshold, no landmark correspondences will be found and the deformation will be primarily driven by the deformations of the neighboring structures. This mechanism is similar to warping the cortex which also presents fundamental morphological differences between subjects (e.g. “one versus two” sulci).
Moreover, the displacement of voxels is not determined directly by the similarity forces, but is subject to constraints as in the case of normal brain images, i.e. the deformation is performed in small steps following local and global transformations (to control the amount of elasticity) and smoothing (to minimize the Laplacian of the deformation). Therefore, although there is no explicit step that guarantees diffeomorphism, the applied deformation mechanism prevents deformation field intersections, even for topologically non-equivalent images. A registration example of two images that are not topologically equivalent is shown in Fig. 3, which illustrates that complete overlap of structures cannot be achieved in that case. Other details on the deformation mechanism and the mechanism for improving the robustness of the method when tumors are present can be found in .
Upon tumor parameters estimation and tumor growth simulation, the registration is performed by relaxing all forces that prioritize the matching of the tumor boundaries. The reason is that the final registration should not be affected by not accurately estimated tumor seeds and by the residual variability in the tumor vicinity which is primarily due to fundamental differences in the growth process between a real and a simulated tumor (e.g. presence of edema and tumor infiltration). Forces that prioritize the matching of tumors come from: (i) the selection of landmark points on the tumor boundaries and (ii) the similarity criterion matching the tumors (SimT). Therefore after estimating the tumor model parameters and simulating tumor growth the algorithm does not select landmark points on the tumor boundaries and also sets the weighting factor w(x, y) to zero. According to this mechanism, the smooth continuity of the deformation field in the tumor neighborhood is preserved.
If a lot of edema is present in the subject's image, an additional label can be used to indicate this area of low confidence in order to avoid operations in this area (similar to masking out). This step however is not very crucial, because the deformation mechanism relaxes the matching forces when the degree of similarity is low, as is the case in tissue displaying edema. The deformation is then driven primarily by the deformations of the neighboring landmark points.
We described the forward model for tumor growth. However, given a patient's image, the tumor model parameters θ n that best represent the specific tumor are unknown and must be estimated. We achieve this goal by solving a constrained nonlinear optimization problem:
where f: n → is an empirical nonlinear function designed to assess the registration accuracy and validity for a given θ, as described next, and L and U are lower and upper bounds on θ, respectively. In this application, we focus only on the model parameters that are patient-specific, such as the origin of the tumor (3D Cartesian coordinates), the amount of brain tissue that died due to the tumor appearance, and the mass effect factor. The bounds L enforce physically valid tumors; also, they constrain the solution to be close to the initial estimate for the solution. For example, the search for the optimal tumor center is bounded within 12mm from the center of mass of tumor in the rigidly registered subject's image. Non-linear constraints on the tumor center describing the complex brain domain - which excludes the ventricles - are not explicitly handled; instead, a penalty function is imposed to inhibit invalid tumor simulations (e.g. outside the brain). The upper bound on the mass effect factor (amount of expansion) is such that the simulated tumor does not exceed in size the tumor in the subject's image.
The objective function used for optimizing the tumor model parameters is based on the hypothesis that the optimal tumor parameters minimize the discrepancies between the co-registered images and also produce realistic deformation maps when trying to match the atlas with simulated tumor with the patient's image. The objective function, f, reflects the success of registration and the validity of the calculated deformation field and is expressed as feature-based similarity and elastic stretching energy, respectively. Specifically, it is defined as the combination of (i) the overlap of the co-registered segmented images (E1), (ii) the feature-based similarity as in equation (1) (E2) , and (iii) the Laplacian of the inter-subject deformation field (E3) which is a measure of smoothness. Overall, the objective function is given by
E1 and E2 are calculated from the warped template image and the subject image, whereas E3 is based on the deformation field values. Specifically, if h(x) defines the mapping from the tumor-deformed template domain T into the domain of the tumor-deformed subject S, h : ΩT → ΩS, then the template image IT at x should look similar to the subject image IS at h(x). For methods that are using segmented images, such as our approach, IS and IT are actually indicator functions and the degree of non-overlap can be defined as
E2 is the attribute vectors' dissimilarity between the template and subject's image
where Sim is defined in (1).
E3 reflects smoothness of the deformation field and is measured by the magnitude of the vector Laplacian (2) of the displacement field u(x) = h(x) − x, as shown below
The constants ck are used to assign different weights on different measures. They were determined by following three steps and satisfying at each step. First, they were set to act as normalization constants scaling each measure Ek, k = 1,‥,3 in the same range. The calculated values were then adjusted to equally penalize the dissimilarity term, which consists of two measures (E1 and E2) with the smoothness term, which consists of a single measure (E3), i.e. c3 received double weight. Finally, experiments were performed with perturbations around the calculated values for a representative range of tumor types. The experiments showed that for a few cases, usually having large tumors, the resulting deformation fields (for the optimized tumor model parameters) were not smooth enough close to the tumor. Therefore the constant c3 was increased more in order to promote smooth deformations for most tumor cases.
wk (x) is used to assign different weights according to the voxel's location x and is selected to decrease with the distance from the tumor boundary for all three measures, and to increase on voxels lying on edges particularly for the image-related measures, i.e., E1 and E2. The assignment of higher weights on voxels lying on edges is due to the distinctiveness of the features of those voxels. More details on the selection of ck and definition of wk (x) can be found in . Ω is the domain over which f is calculated, and it is defined in the subject's image outside the tumor and within a small distance from tumor boundary, where the effects of mis-registration are expected to be more prominent. The distance from the tumor boundary defining Ω is decreased after each resolution level since the estimate of the tumor model parameters is improved. The part of the image, which has no tissue label due to low confidence in tissue segmentation (e.g. peritumoral edema), is excluded from Ω.
The tumor parameters θ are optimized by solving equation (2). We chose a derivative-free method to obtain a solution to equation (2), because the objective function contains discontinuities and the approximations of the derivatives with finite differencing may be unreliable. Pattern Search methods are suitable for such problems . They are directional methods that make use of a finite number of directions with appropriate descent properties. Since the majority of the computational cost is the function evaluation, we apply an Asynchronous Parallel Pattern Search method, called APPSPACK , which takes advantage of parallel platforms. APPSPACK targets simulation-based optimization problems characterized by a small number of parameters. The “asynchronous” algorithm dynamically initiates actions in response to messages, rather than routinely cycling through a fixed set of steps, allowing to effectively balancing the computational load across all available processors. This property is important for our application, since the function evaluation cost varies significantly according to the parameters applied (simulations of small tumors or simulations starting with large initial seeds terminate faster).
Due to the complicated form of our objective function, it is not guaranteed that a global optimum exists. We set an upper limit of 120 iterations in the optimization process. Usually a minimum was reached and the optimization terminated before reaching this limit. However in a few cases, usually with tumors very close to boundaries (e.g. between the ventricles), the defined upper limit was not large enough to achieve convergence, so the best value until this moment was selected as optimum.
One application of inter-subject registration of brain tumor images may involve identifying differences in the tumor growth process among populations of subjects. For example, the tumor mass effect has been used as a descriptor for classifying gliomas according to their clinical grade  or as an independent predictor of survival . There is therefore the need for an objective method that quantitatively characterizes the mass effect. In this study, we investigated how well the estimated tumor seeds and the calculated transformations help predicting the mass effect.
In order to obtain an indicator of mass effect based on the proposed registration method, we evaluated how much the calculated deformation fields deviated from the range observed in a normal population. Specifically we calculated the Jacobian determinant of the deformation fields that spatially warp the atlas (used in this study) to the brain tumor images, as well as to a population of healthy subjects (60 individuals with age less than 68 years). We smoothed all Jacobian images with a Gaussian filter (with standard deviation = 1), in order to reduce the noise and small misregistration effects, and calculated the voxel-wise standard score, in order to normalize the amount of deformation at each brain location, since the tumors appear in different locations in the brain for different subjects. The standard score Zi for each voxel i is Zi = (Ji-μi)/σi, where Ji is the smoothed Jacobian at voxel i, and μ and σ are the mean and standard deviation of the normal population at voxel i. The sum of the standard scores over the tumor seed defined in the common atlas space is a quantity that represents the distance between the total deformation and the population mean and is used as an indicator of mass effect.
In order to assess the results, the mass effect estimated by our method was compared against expert ratings. Specifically, similarly to the ratings definition provided in , two independent expert neuroradiologists provided their scores by examining the images and quantifying the tumor mass effect using the following scale:
|0 – Absent||: adjacent tissue is not compressed, the structure is preserved|
|1 – Light||: shape change, thickening or thinning of adjacent structure|
|2 – Middle||: displacements, adjacent tissue and/or ventricles are slightly shifted|
|3 – Heavy||: large displacements, middle line structure is shifted to the opposite direction|
Under the assumption that both measurements of mass effect (through our modeling framework or by expert ratings) are susceptible to errors, high correlation between the two measurements is an indicator that the measurements are consistent and therefore acts as a means of validation of our registration/estimation framework. Moreover, it should be noted that our method for quantitatively characterizing the mass effect allows studying the tumor growth process among populations of subjects based on reproducible and rater independent techniques.
We applied the proposed framework for registration of 21 brain MRI datasets including tumors of different types, grade and sizes. The brain masses were histologically diagnosed and graded based on World Health Organization (WHO) criteria. Specifically the data sets included 5 metastases (MET), 4 oligodendrogliomas (OLIGO) (1 of grade II and 3 of grade III), 2 astrocytomas (ASTRO) (1 of grade II and 3 of grade III), 1 ganglioglioma (GANGLIO), 1 oligoastrocytoma grade III, and 8 glioblastomas (GBM) (grade IV). All images were registered with a normal brain image serving as a template, with image size 256×256×198 voxels and voxel size 1×1×1 mm3. One registration example of a diffuse and infiltrative tumor is shown in Fig. 1 as part of the illustration of the proposed framework. Fig. 4 shows the registration of another brain tumor case and specifically it illustrates the three main steps involved after estimation of the optimal tumor model parameters: (a) seeding the atlas, (b) simulating tumor growth in the atlas and (c) registering the deformed atlas to the subject.
Fig. 5 illustrates atlas-based segmentation results obtained by warping the normal atlas (shown on the first row) onto three brain images with tumor. Three labeled regions on the atlas, i.e., thalamus, caudate nuclei, and ventricles (shown with lila, green and pink, respectively), are mapped after registration onto the patient's images, in order to visualize the registration performance on anatomical structures deformed by the tumor. The results show that the segmented structures are consistent with the patient's anatomy.
Besides providing the qualitative results based on visual assessment, we calculated a quantitative rater-independent measure, such as the surface distance of the ventricles (VN-dist) between the co-registered images. We calculated VN-dist as the mean Euclidean distance of the ventricular boundaries in both directions, from the patient's image to the warped atlas image and reversely. We selected the ventricles for validation, because they are structures with distinct boundaries providing accurate (automatic) segmentation. Table 1 shows VN-dist for 21 patient's images calculated over the whole ventricular boundary (top rows), as well as only on the part that lies closer to the tumor, in order to emphasize possible limitations of the method due to the presence of tumor (bottom rows). These results show that the distance of ventricles is larger in the tumor vicinity than over the total ventricular boundary. Specifically, the error is at voxel accuracy for the ventricular part that is far from the tumor and at the order of the diagonal voxel distance for the region close to the tumor, for all cases except for one. This worst case, which is highlighted in Table 1, has an error of 2.758 mm and 3.662 mm respectively and is shown in Fig. 6. It can be seen that this is a very difficult case, from image registration perspective, since the anatomy is highly obscured by tumor infiltration and edema, and the ventricles cannot be clearly segmented.
Moreover, the current method is compared with ORBIT for the 10 patient images presented in . The change of VN-dist when using the current method is shown as a percentage rate in Table 2 for the whole tumor boundary as well as for the tumor vicinity. Positive numbers indicate decrease of VN-dist, whereas negative numbers indicate increase. The numbering of the patients corresponds to the numbering in . It can be noticed that for the majority of cases the error was reduced with the current method by up to 55%. The overall error reduction might be due to a more accurate estimate of the tumor model parameters or due to a better simulation of the tumor-induced deformation providing a better (possibly smoother) initialization to the registration algorithm. In some cases, the error in the tumor vicinity increased, whereas the total error almost stayed constant (patient 6) or decreased (patient 3). First, it should be noted that the tumor in those two patients was located far from the ventricles and therefore the error in the tumor vicinity, evaluated on a very small part of ventricular boundaries, is not very representative. Moreover, the value of VN-dist for patient 3 was especially small (only 0.578 mm which is the smallest value across all patients) and therefore the 20% change (= 0.116 mm) is not considered significant.
The parameters estimated with the current method are compared with the ones estimated in . The difference (Δ) between the corresponding parameters (θ) is shown at the bottom part of Table 2. The median difference (over the 10 patients) in tumor seed location is 11.7 mm, in seed size 2.15 mm and in estimated final tumor size (due to mass effect) 1.15 mm. The difference in the estimated parameters possibly indicates the presence of local minima in the objective function. By using different tumor growth simulation approaches and optimization methods the optimization process converged to different solutions.
In order to further evaluate the registration accuracy and compare with our previous framework, we also calculated the volume of overlap between the subject and the warped template for all three brain tissue classes (WM, GM and VN) using the DICE metric . For the same 10 patient shown in Table 2 the average and standard deviation for the DICE score (without removing areas with confounding effects, like edema) was 0.744±0.025, 0.668±0.025 and 0.735±0.101 for WM, GM and VN respectively, when using the current framework, versus 0.714±0.026, 0.637±0.026,0.703±0.106 when using the framework in . The volume of overlap has increased for all three brain tissue classes. We don't measure the volume of overlap for the tumors because we treat the tumors as infiltrating structures with only approximate boundaries.
The construction of diffeomorphic registration fields is important in order to preserve topology (one-to-one mapping) and allow inversion of the mapping. Ashburner presents in  a review on diffeomorphic registration approaches and provides detailed definitions. In order to study the properties of the deformation field constructed by our method, we calculated the determinant of the Jacobian of the total deformation field, which maps the normal atlas to the subject with tumor (forward field), and also for the inverse field, which normalizes the images of tumor pathology into the atlas space. Negative Jacobian determinants in the deformation component simulating tumor growth are mostly attributed to the selection of a linear elasticity model for the brain in conjunction with not adequately small pressure increments. Decreasing the pressure increments or incorporating a non-linear elastic model might preserve the one-to-one mapping property, but also increases the computational cost. Since the parameters were not adjusted for each case separately (e.g. same pressure increments were applied for all cases), 2 (out of 21) forward deformation fields have small negative Jacobian determinants. Deformation fields with negative Jacobian determinants could be reconstructed subsequently by applying the method described in  to recover topology preserving mappings. The results presented here are produced without the reconstruction step.
The inverse of the spatial transformation is calculated by an iterative algorithm using linear interpolation and is applied for spatial normalization of the brain tumor images. Fig. 7 shows examples of brain images with tumor warped in the normal atlas space. This warping causes relaxation of the mass effect and removal of the inter-subject differences facilitating the detection of two tumorous regions: (i) the initial seed (as estimated by our method) which indicates the location of initial tumor appearance and brain tissue loss, and (ii) the surrounding region that is infiltrated by tumor or edema. The spatially normalized images look realistic without artificial deformations due to tumor shrinkage. The estimated tumor origin is displayed for all subjects in Table 3. However it should be noted that these results are mostly qualitative since the main focus of the method is to improve registration accuracy and not to accurately estimate the complex biological process of tumor growth.
It is interesting to investigate whether the estimated tumor seed location and the atlas-to-subject deformation field can be used to predict the mass effect of the respective tumor. We used as an indicator of mass effect the sum, over all voxels in the estimated tumor seed region, of the normalized (using statistics from a healthy population) Jacobian determinant, as explained in section 2.4. We calculated these scores for the 21 subjects in this study and compared them against visual ratings provided by 2 expert neuroradiologists. The visual ratings of mass effect, used as gold standard and averaged over the two raters, are shown in Table 3. None of the cases was rated with 0 (absent mass effect) and the values never fluctuated by more than 1 scale between the two raters. The correlation between our scores estimating tumor mass effect and the visual ratings is 0.763 and 0.618 for the two raters respectively and 0.744 for the average ratings, whereas the correlation between the 2 visual ratings is 0.679.
We also investigated if the brain size is correlated with the mass effect. Specifically, in order to remove the effects due to differences in brain size we normalized our measurements by dividing with the average Jacobian determinant inside the brain. This additional step didn't help increasing the correlation with the visual ratings. This could be due to the noise present in the data or the noise introduced by the processing steps, but it also could indicate that the mass effect is not correlated with the brain size.
Error! Reference source not found. shows on the top row the normal population mean and standard deviation of the Jacobian determinant and on the bottom row, left, one example of normalized Jacobian determinant calculated from the deformation field that maps the atlas to a brain image with tumor. It is easy to notice that the deformation in the tumor seed region is significantly larger than the one observed in the normal population. The average image of the Jacobian determinants of the normal population has a standard deviation over all voxels equal to 1.042 with a mean value of 1.0 (when size differences are removed) and a maximum value of 2.703.
In this study we present a method for brain atlas registration in the presence of tumors. Such a method (i) makes possible the pooling of data from different patients into a common space, thereby enabling the performance of group analyses, (ii) allows the mapping of atlases with segmented structures of interest into the patient's image space for optimization of the surgical or radiotherapy treatment planning, (iii) estimates the brain tissue loss and replacement by tumor, as well as the mass effect of the tumor, thereby helps identifying differences in the tumor growth process among populations of subjects. The method is based on the idea of decoupling the total deformation into the tumor-induced deformation and the deformation due to inter-subject differences, similarly to . The current framework adds to previous work of our group in this area, in the following respects:
The current framework uses a piecewise linear elasticity model and regular grids (PLE simulator), versus a finite element model of nonlinear elasticity and unstructured meshes in  and a local-PCA based model in . PCA-based tumor growth simulations are extremely fast, since they are achieved via linear combination of principal components (deformations), thereby leaving the burden of simulations to off-line training using costly biomechanical models. However PCA-based simulations are not very accurate, since the expansion coefficients are not known in advance and can only be approximated. Moreover, PCA-based models are not flexible, since they can only reproduce deformations in the range of parameters used during training. On the other hand, nonlinear biomechanical simulators  are flexible and more accurate, but also computationally very expensive. As shown in , the use of the robust and computationally efficient PLE simulator did not significantly affect the final registration accuracy in comparison to nonlinear biomechanical simulators. Therefore the PLE simulator seems to provide the best trade off between accuracy, flexibility and computational cost and was therefore chosen here. It should also be noted that the PLE simulator, is a stand-alone program which does not require the use of commercial packages, such as ABAQUS.
It applies the APPSPACK optimization package for parallel optimization (using mpi) of 5 tumor-related parameters, versus the Downhill Simplex method in  for serial optimization of 4 parameters, and a statistical approach in , not based on optimization, applied for a different sets of parameters (edema was also defined in ). This allows us to search over a large range of parameters in a computationally efficient way.
It applies irregular shaped seeds, versus spherical seeds in . The use of irregular seeds allows the creation of an atlas with tumor that is more similar to the subject's image.
It maximizes a similarity criterion that matches both the brain structures and the tumor geometry using locally adapted weights, as in , versus the regular HAMMER registration algorithm developed for normal brains in . Since the tumor in atlas and patient's image might be in non corresponding anatomical locations (e.g. if the estimate of the tumor location is inaccurate), it is not guaranteed that a diffeomorphic deformation field, that maps one image to the other, will exist. Therefore, the registration method here, and in , applies different deformation strategy close and far from the tumor, in order to maintain robust registration of the healthy part of the brain.
The implementation of different registration methods and coupling of them with a tumor growth simulation model is also possible. We chose a registration method based on labeled (segmented) images because (i) methods that are based on gray level information (e.g. maximization of mutual information) will be affected by the choice of intensity profile for the simulated tumor; modeling intensity changes due to tumor evolution is an open and challenging problem, and (ii) the tumor growth model applied is a pressure-based biomechanical one, which preserves labels. Extensions of the tumor simulation methodology to include modeling of tumor infiltration  might require a different or modified registration method.
The warped images in atlas space or subject's space, shown in figures 4 to to7,7, indicate that the proposed method performs well in both areas around the tumor and in the healthy portion of the brain. The distance of the ventricular boundaries between warped atlas and patient's images (Table 1) over the 21 patients was 1.097±0.459 mm. When evaluated only in the tumor vicinity, where the largest errors are expected, the ventricular distance was 1.459±0.752 mm. These results show that the average registration error of the ventricles is at the order of diagonal voxel distance. Also, the results in Table 2, as well as the DICE metric measuring overlap of WM, GM and CSF, indicate that the registration accuracy overall has improved as compared to ORBIT . More validation is required in order to assess the accuracy on other structures of interest.
The computational time is considerably decreased in comparison to our previous work , due to the use of parallel optimization for the parameters of the tumor growth model. The execution time for optimizing the parameters in the subdomain, only in the middle resolution level and using 20 parallel processors (CPUs) in a Linux cluster consisting of 102 CPUs in total (Intel Xeon CPU 3.4 GHz and Intel Woodcrest Core2Duo CPU 1.6 GHz) with 4 or 8 GB RAM, ranged from 1.1h to 4h. The cost depends on the load of the processors, the closeness of the initial estimate to the optimum and the size of the tumor (tumor growth simulation is more expensive for larger tumors). The whole framework which includes all steps (implemented serially in a single processor) before and after optimization required ~2 additional hours. These steps include the calculation of the global alignment, an initial estimate of the deformation field in the whole image domain before tumor parameters optimization and the calculation of the deformation field in the whole image domain after optimization (refinement in full resolution).
The tumor model parameters are important for achieving topological equivalence between the atlas and patient's image, and any assumptions on those limit the possible simulation outcomes. However, it seems that the estimation of the tumor parameters is more important in the case of spatially normalizing brain tumor images into a common (atlas) space, than in the case of atlas-based segmentation (when a normal atlas is mapped into the patient's image space). Small variations in the tumor seed have small effect in the registration of the atlas to the subject space, although the image differences are prominent in the original, undeformed atlas space. Therefore, methods that have been originally developed with the purpose of atlas-based segmentation often apply simplistic models of tumor growth (e.g. radial expansion) , ignore the information contained within the lesion , or do not perform any algorithmic adjustments due to the presence of lesions . On the other hand, for the purpose of performing population-based statistics, the flexibility of the tumor growth simulation and the ability of estimating the model parameters describing the tumor development, becomes essential. Although the parameters estimation requires additional computational cost, we believe that the better understanding of the tumor development, such as knowledge of the origin of the tumor or the amount of tissue death, is of high clinical importance.
As an application example we used the spatial transformations to quantify the mass effect of the growing tumors. The mass effect has been used as a feature in classification of brain tumors into different types. The automated calculation of the mass effect helps avoid descriptive criteria that are rater-dependent and require prior knowledge (the help of experts). We based our calculations on the deviation of the deformations from the range observed in a normal population. The correlation between the estimated mass effect and the expert ratings was higher than the inter-rater variability.
A current limitation of our approach is that it is based on the prior tissue segmentation, which poses considerable difficulties in practice, especially in the region around the tumor that often displays edema and infiltration. The proposed framework is mostly suitable for tumors with distinct tumor boundaries, which are not very challenging from segmentation perspective. However, it is important to note that the current implementation is robust to inaccuracies in tumor segmentation because the simulated tumor is not forced to expand until it reaches the manually segmented tumor in the patient's image. The amount of tumor expansion is rather determined by optimizing the defined optimality criterion. Therefore we do not expect the manual segmentation of the tumor boundary to match with the simulated tumor. One future extension of the proposed registration method could be the transition from hard tissue segmentation into a fuzzy or probabilistic segmentation framework, which is more appropriate for the inherently diffuse and infiltrative brain tumors cases, since in those cases the tumorous area can only be characterized through probabilistic tissue abnormality maps. On-going work in our laboratory  investigates pattern classification methods that use multi-acquisition imaging profiles, including T1, T1-enhanced, FLAIR, DTI, and aims to achieve a more accurate tissue classification, thereby assisting in the registration process.
The authors would like to thank Drs. Sumei Wang and Naomi Morita in the Department of Radiology, University of Pennsylvania for providing the visual ratings of mass effect and Drs. Elias R. Melhem and Ron Wolf in the Department of Radiology, University of Pennsylvania, for providing the patient's data. We also thank Dr. R. Nick Bryan for discussions regarding the estimation of mass effect.
A portion of this work was presented at the SPIE Medical Imaging 2008: Image Processing, San Diego, California, 2008.
Grant support: NIH Grant RO1- NS042645.
Publisher's Disclaimer: This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.