The use of strain tensor to estimate thoracic tumors deformation Darek Michalski, M. Saiful Huq, Greg Bednarz, and Dwight E. Heron Citation: Medical Physics 41, 073503 (2014); doi: 10.1118/1.4884222 View online: http://dx.doi.org/10.1118/1.4884222 View Table of Contents: http://scitation.aip.org/content/aapm/journal/medphys/41/7?ver=pdfcov Published by the American Association of Physicists in Medicine Articles you may be interested in An initial study on the estimation of time-varying volumetric treatment images and 3D tumor localization from single MV cine EPID images Med. Phys. 41, 081713 (2014); 10.1118/1.4889779 Modeling lung deformation: A combined deformable image registration method with spatially varying Young's modulus estimates Med. Phys. 40, 081902 (2013); 10.1118/1.4812419 Toward in vivo lung's tissue incompressibility characterization for tumor motion modeling in radiation therapy Med. Phys. 40, 051902 (2013); 10.1118/1.4798461 Deformable image registration of heterogeneous human lung incorporating the bronchial tree Med. Phys. 37, 4560 (2010); 10.1118/1.3471020 Reconstruction of a time-averaged midposition CT scan for radiotherapy planning of lung cancer patients using deformable registrationa) Med. Phys. 35, 3998 (2008); 10.1118/1.2966347
The use of strain tensor to estimate thoracic tumors deformation Darek Michalski,a) M. Saiful Huq, Greg Bednarz, and Dwight E. Heron Department of Radiation Oncology, University of Pittsburgh Cancer Institute, Pittsburgh, Pennsylvania 15232
(Received 25 June 2013; revised 2 June 2014; accepted for publication 3 June 2014; published 24 June 2014) Purpose: Respiration-induced kinematics of thoracic tumors suggests a simple analogy with elasticity, where a strain tensor is used to characterize the volume of interests. The application of the biomechanical framework allows for the objective determination of tumor characteristics. Methods: Four-dimensional computed tomography provides the snapshots of the patient’s anatomy at the end of inspiration and expiration. Image registration was used to obtain the displacement vector fields and deformation fields, which allows one for the determination of the strain tensor. Its departure from the identity matrix gauges the departure of the medium from rigidity. The tensorial characteristic of each GTV voxel was determined and averaged. To this end, the standard Euclidean matrix norm as well as the Log-Euclidean norm were employed. Tensorial anisotropy was gauged with the fractional anisotropy measure which is based on the normalized variance of the tensors eigenvalues. Anisotropy was also evaluated with the geodesic distance in the Log-Euclidean framework of a given strain tensor to its closest isotropic counterpart. Results: The averaged strain tensor was determined for each of the 15 retrospectively analyzed thoracic GTVs. The amplitude of GTV motion varied from 0.64 to 4.21 with the average of 1.20 cm. The GTV size ranged from 5.16 to 149.99 cc with the average of 43.19 cc. The tensorial analysis shows that deformation is inconsiderable and that the tensorial anisotropy is small. The Log-Euclidean distance of averaged strain tensors from the identity matrix ranged from 0.06 to 0.31 with the average of 0.19. The Frobenius distance from the identity matrix is similar and ranged from 0.06 to 0.35 with the average of 0.21. Their fractional anisotropy ranged from 0.02 to 0.12 with the average of 0.07. Their geodesic anisotropy ranged from 0.03 to 0.16 with the average of 0.09. These values also indicate insignificant deformation. Conclusions: The tensorial framework allows for direct measurements of tissue deformation. It goes beyond the evaluation of deformation via comparison of shapes. It is an independent and objective determination of tissue properties. This methodology can be used to determine possible changes in lung properties due to radiation therapy and possible toxicities. © 2014 American Association of Physicists in Medicine. [http://dx.doi.org/10.1118/1.4884222] Key words: strain tensor, image fusion, image processing, mechanical properties 1. INTRODUCTION Physiologically induced organ motion in thorax introduces certain degree of uncertainty that must be dealt with during the treatment planning and delivery.1, 2 The kinematics of the lungs affect the tumor, which might move during the breathing cycle even up to 5 cm for supradiaphragmatic masses as we observed in our clinic. Tumor motion can occur in all three anatomical directions and its trajectory might exhibit hysteresis. Tumor motion has profound consequences for the entire clinical workflow of the radiation treatment process. The helical computed tomography (CT) can misrepresent tumor shape, location, and size if its motion magnitude exceeds 0.5 cm. This can be adequately addressed via fourdimensional (4D) CT scans. The knowledge of the tumor motion amplitude allows for the accurate determination of margins about the clinical tumor volume (CTV) to define the integral tumor volume (ITV). This weighs in the decision process of the radiation oncologist whether or not elevated doses can be prescribed, or gated treatment should be applied. In the end, it affects the radiation treatment outcomes, which for lung cancer patients do not exhibit expected efficacy. 073503-1
Med. Phys. 41 (7), July 2014
4D CT scanning methodology allows for direct tumor observation over the entire breathing cycle. The AAPM Task Group 76 report2 points out that tumor deformation is usually not quantified or considered in various aspects of motion management in radiation oncology. It is tacitly assumed that the magnitude of relevant margins would address this. There have been two reports in the literature on the lung tumor changes during breathing cycle utilizing image registration for contour comparisons. Liu et al.3 investigated motion of 166 tumors and concluded that tumor motion does not manifest major deformations to tumor contours. They concluded that rigid-body motion suffices to model tumor motion. Similar in methodology, a study of 30 lung tumors4 reported that the comparison of the GTV contours at the end of exhalation (EE) and inhalation (EI) reveals mainly translational change. We are revisiting the topic of lung tumor deformation applying biomechanical framework. We also use image registration to the EE and EI 4D CT volumes but with totally different methodology of quantitative evaluation. Aforementioned two studies relied on comparison of shapes of tumors. The overlap of polygons that built the tumor contours gauged the tumor deformation. We consider the entire volume of the tumor. This
0094-2405/2014/41(7)/073503/6/$30.00
© 2014 Am. Assoc. Phys. Med.
073503-1
073503-2
Michalski et al.: Strain tensor to estimate thoracic tumors deformation
evaluation is direct and does not resort to comparative methods. Deformation is directly quantified by the strain tensor. We utilize scalar and tensorial metrics to evaluate deformation for the entire tumor volume in contrast to the previous studies. This method goes beyond the heuristic of the comparison, which in practice is equivalent to visual inspection. This strain-based approach can be applied to any region of interest since it generically deals with finding a homologous region of interest in images and directly derives all the required parameters to gauge the deformation. We use the lung tumor cases to illustrate this methodology. The deformation is the change of size and/or shape of a given volume of interest. The kinematics of the thorax during breathing can be characterized by its deformation and motion. The same can be said about thoracic gross tumor volume (GTVs). Applying this analogy one may consider to use mathematical constructs inherent to the description of the deformations, which is a strain tensor. In general, this approach allowed to utilize other 4D imaging modalities to evaluate regional strains in heart and in other organs. The temporal resolution of patient’s imaging provides the same topological view of the anatomy but it is also able to provide the insight into cyclic changes of the anatomy during the scan. In other words, this variability can be seen in terms of regional deformations and can be quantified with a strain tensor. 4DCT scans of thoracic patients enable to view physiologically induced motion of the entire thorax. 2. MATERIALS AND METHODS The deformation of a given tissue element can be determined if its displacement is known. The latter can be obtained from 4DCT scans, which were acquired using a GE Discovery ST scanner (GE Medical Systems, Milwaukee, WI) with a real-time position management (RPM) respiratory gating system (Varian Medical Systems, Palo Alto, CA). Retrospective phase-based sorting yields ten temporally equispaced along the breathing cycle CT phases. The breathing cycle starts at EI. All patients during the scanning were immobilized in a supine position with their arms extended above their heads. Patients were coached using an audio directive (“breathe in”). The frequency of the command was matched to a given patient’s normal respiratory period. For the evaluation of tumor deformation, patients whose tumors exhibited at least 0.5 cm motion were selected. Their resampled isotropic 4D CT phases corresponding to EI (CT0) and EE (CT50) were used. The voxel size was the slice thickness. The latter was used as a reference image scan and was registered to the CT0 scan in the Lagrangian framework. A free form deformation image registration provides the displacement field, u = x − X, where X and x are the locations of the homologous tissue elements in CT50 and CT0, respectively. Equivalently we obtain the mapping, , from the CT50 to the CT0, x = (X). The mapping is characterized by its Jacobian matrix, F = ∇. If the mapping is viewed in terms of the coordinate changes, then from the change of variables theorem, the deformed volume element dv is related to the undeformed one dV with dv = J dV. Medical Physics, Vol. 41, No. 7, July 2014
073503-2
Thus, the Jacobian at a given point, J = det(F ) must be positive to guarantee physical feasibility of the transformation.5 The Jacobian is useful since it provides a first insight into the deformation.6 For example, it was used to capture the tissue deformation variability of the group of patients and incorporated in a mathematical construct of a cost function in image registration.7 Yet, even without volume change (J = 1) the medium might exhibit deformation. The parameter that captures this behavior is a strain tensor. Strain tensors might assume various mathematical forms. For a Lagrangian framework F = RU, where R is a proper orthogonal rotation and U a symmetric positive definite matrix. The latter serves in the definition of the Doyle-Ericksen class of strain tensors with the following form E(n) = 1/n (Un − I), where n can be any real number.5, 8, 9 For n = 2, it is the Green-St. Venant strain tensor, E(2) = 1/2 (FT F − Id ), where Id is an identity matrix. Strain tensors are not affected by the rigid transformation so the only change they capture is the deformation. The departure of the strain tensor from the identity matrix gauges the departure of the medium from rigidity. Similar is the right Cauchy-Green strain tensor, computed in this project, C = FT F. It is the identity matrix for the rigid transformation. For each point in the GTV, we can determine its tensorial characteristic and average them to gauge the GTV behavior. Specific mathematical characteristics of the tensorial matrices are not fully congruent with the Euclidean matrix norm as applied in the traditional elasticity, e.g., averaging the tensors C = 1/N i Ci might result in det(Ci ) < det(C) for all det(Ci ). An ingenious Log-Euclidean metric has been proposed10, 11 that addresses these problems. There are other metrics,12–15 which amend the shortcomings of Euclidean norms but the Log-Euclidean is computationally most efficient. This metric determines the distance, d, between tensors T1 and T2 as follows:10, 11 d(T1 , T2 ) = ||ln(T1 ) − ln(T2 )||, where ||. . . || denotes the Euclidean norm. As it was proposed by Arsigny et al.,11 it is the Frobenius norm applied to logarithms of tensors. The Log-Euclidean metric allows 11 for the computation of the average strain tensor over all N voxels of the GTV as C = exp(1/N i ln(Ci )). The logarithm of the tensor in the Log-Euclidean framework is the inverse of the matrix exponential. The logarithm or the exponential of the tensor, T, is obtained by applying the operation to its eigenvalues. These mathematical constructs are well defined for tensors, which are symmetric definite positive matrices.11, 16 Computation of the averaged strain tensor allows for the examination of its anisotropy. The commonly used gauge is a fractional anisotropy,17 FA, which is a scalar value that aims at describing tensor’s anisotropy in terms of the normalized variance of the tensor’s eigenvalues, λi . It is defined √ as FA = ( 3/2 )[ i (λi − λ)2 ]1/2 /[ i (λi )2 ]1/2 where λ = 1/3 i λi . However, another measure of the anisotropy is the geodesic distance in the Log-Euclidean framework of a given strain tensor to its closest isotropic counterpart.13, 18, 19 The latter is the tensor with the same determinant and identical eigenvalues. This gauge was named geodesic anisotropy and can be expressed as GA = [ i (ln(λi ) − lnλ)2 ]1/2 where lnλ=1/3 i ln(λi ).
Michalski et al.: Strain tensor to estimate thoracic tumors deformation
073503-3
0.4
073503-3
0.4 (a)
(b)
0.3
1
0.2
0.2
0.1
2 11
5
Distance
Distance
0.3
8
3
7 13 12
15
6 14 10 9
0.1
4
0 0
50
100
150
0 0
0.5
1
1.5
2
2.5
3
3.5
4
Motion [cm]
Tumor Volume [cc]
F IG . 1. The plot of the Log-Euclidean (䊏) and Euclidean (◦ ) distance of the averaged GTV strain tensors from Id as a function of the GTV volume [cc] (a) and the plot of Log-Euclidean distance for these tensors as a function of the GTV motion amplitude [cm] (b). For the latter, the points are annotated with the patient number ascending with the GTV magnitude.
3. RESULTS AND DISCUSSION The image registration was robust. The GTV motion ranged from 0.64 to 4.21 cm with its average 1.20 cm. The GTV size ranged from 5.16 to 149.99 cc with its average 43.19. Figure 1(a) shows the Log-Euclidean and the Frobenius distance of the averaged strain tensors to the identity matrix. The distances computed with the Log-Euclidean and the traditional matrix norm are very similar. The Log-Euclidean values range from 0.06 to 0.31 with the average 0.19. The Euclidean values range from 0.06 to 0.35 with the average 0.21. The RMSE between the distances computed with these metrics is 0.02. This plot depicts this distance as a function of the GTV magnitude. The values do not exhibit correlation with the GTV size. The Pearson correlation coefficient is −0.04. Since the Log-Euclidean and the Euclidean distances are similar Fig. 1(b) shows only the former as a function of GTV motion. The points are annotated with a GTV number. The GTVs are ordered along ascending GTV magnitude. There is no correlation with GTV motion amplitude as well. Figure 2 shows GA and FA for all the patients numbered in ascending order of the GTV magnitude. The values of GA range from 0.03 to 0.16 with the average 0.09. The values of FA range Medical Physics, Vol. 41, No. 7, July 2014
from 0.02 to 0.12 with the average at 0.07. For all the cases, GA factor is higher. The RMSE for GA and FA is 0.03. There are no correlations of the magnitudes of GA and FA with the size of the GTVs. The Pearson correlation coefficient is 0.06. There is no correlation with the motion amplitude as well. It is instructive to examine these measures using synthetic data. The distance from the identity matrix can be computed using a synthetic tensor with the following eigenvalues λ1 = exp(x), λ2 = λ3 = exp(−x), where x varies between 0 and 1 (λ1 ∈ [1.0, 2.71], and λ2 = λ3 ∈ [1.0, 0.37]). This corresponds to the variation of the Jacobian from 1.0 to about 0.61. Figure 3(a) shows the plot using the Log-Euclidean as well as the Frobenius norms. The plot in Fig. 3(b) shows the results for the following eigenvalues λ1 = exp(x), λ2 = λ3 = exp(−x/2). This plot shows a scenario for a deforming medium that does not change the volume since its Jacobian remains equal to 1 for x ∈ [0, 1]. The difference between the results obtained with the Log-Euclidean and the Frobenius norm in the latter case is more pronounced. The distance from the identity matrix for the isochoric deformation is smaller than for the volumetric deformation since the
0.15
GA/FA
The registration method employed here is based on an iterative and computationally robust Demons algorithm.20–22 It yields voxel displacement field, u, which in turn provides F = Id + ∇u. The algorithm is fast and accurate as reported in the literature. This image registration method was well-tested and it was used by many research groups for various clinical data.23–31 The implementation uses a hierarchical multiscale approach with a three level Gaussian pyramid.32 In order to prevent unphysical folding, the magnitude of the updates was always smaller than 0.4 of the voxel size.33–36 The gradient of the displacement field was computed using analytical form of the derivatives of cubic B-splines. The voxels of the GTVs were extracted using binary mask derived from the structure dicom files. The image registration program was implemented in C on a Linux box.
0.1
0.05
0
0
5
10
15
Patient Number F IG . 2. The plot of GA (䊏) and FA (◦ ) for all patients. Abscissa refers to patients ordinal number. The patients are numbered in the ascending order of the GTV size.
Michalski et al.: Strain tensor to estimate thoracic tumors deformation
073503-4
2
2 (a)
(b) 1.5
Distance
Distance
1.5
1
0.5
0 0
073503-4
1
0.5
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
0 0
1
0.1
0.2
0.3
0.4
x
0.5
0.6
0.7
0.8
0.9
1
x
F IG . 3. The distance of the synthetic C tensor from Id calculated with the Log-Euclidean (——) and the Frobenius (– – –) norms. The abscissa covers the values of x that parametrize the eigenvalues of the tensor; (a) with eigenvalues λ1 = exp(x), λ2 = λ3 = exp(−x) (volumetric deformation) and (b) with λ1 = exp(x), λ2 = λ3 = exp(−x/2) (distortional deformation).
λ2 for this case varies less. For both types of eigenvalues, the Frobenius distance tends to be higher as x increases but for small values of x both distances are similar. The average distance of the GTV strain tensors from the rigidity for the Log-Euclidean metric was 0.19 and for the Frobenius metric 0.21. Both metrics for the strain tensors also yielded the distances of similar magnitude. One can see that these ordinate values in Fig. 3 correspond to the value of the abscissa of about x = 0.1. In this region, both plots also show similar distance for both metrics. The comparison to the synthetic cases suggests that the deformation is in the proximity of the identity matrix. This indicates that the tumor deformation is inconsiderable. Figure 4(a) and 4(b) show, respectively, corresponding plots of GA and FA as calculated with synthetic eigenvalues as described for Figs. 3(a) and 3(b). The difference between GA and FA is less pronounced for the case when the medium deforms but does not change the volume. This is due to the different values of λ2 ( = λ3 ) for the dilational and the distortional case. The synthetic case shows that GA > FA. However, there is no clear relationship between GA and FA. Their magnitudes are also incongruent, FA ∈ [0, 1] and GA ∈ [0, ∞). As for the values of GA and
FA obtained for the averaged strain tensors they indicate small deviation from an isotropic case. In other words, anisotropy is inconsiderable. Since the strain tensor can be multiplicatively ˜ the distortional comdecomposed as C = [(det(C))1/3 Id ] C, ˜ ponent, C, demonstrates an inconsiderable deformation. This component is responsible only for the deformation and does not change the volume. The examination of thoracic tumors within the biomechanical framework agrees with previous findings3, 4 that GTVs demonstrate rather inconsiderable deformation during respiration. In contrast to these studies, we only considered tumors with motion exceeding 0.5 cm. The thoracic tumors are denser than surrounding tissue, which is unable to deform the GTV during respiration. However, theoretically one cannot exclude the existence of the cases where GTV can deform, e.g., if one end of the GTV is attached to the rib cage. In general, the biomechanical framework allows for the quantitative description of the tissue or anatomical regions of interest. Such regional characterization of volume of interests can be used in the objective evaluation of changes of the anatomy during the treatment or after the treatment. As evidenced37, 38 radiation treatment causes toxicity and change
2
2
(a)
(b) 1.5
GA/FA
GA/FA
1.5
1
0.5
0 0
1
0.5
0.1
0.2
0.3
0.4
0.5
x
0.6
0.7
0.8
0.9
1
0 0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
x
F IG . 4. GA (——) and FA (– – –) gauges of a synthetic tensor anisotropy. The abscissa covers the values of x that parametrize the eigenvalues of the tensor; (a) with eigenvalues λ1 = exp(x), λ2 = λ3 = exp(−x) (volumetric deformation) and (b) with λ1 = exp(x), λ2 = λ3 = exp(−x/2) (distortional deformation). Medical Physics, Vol. 41, No. 7, July 2014
073503-5
Michalski et al.: Strain tensor to estimate thoracic tumors deformation
in pulmonary function. For the latter, the tests can only indicate the overall functional change in respiration. However, lung functionality is not homogenous and its map can be revealed by nuclear imaging or other imaging techniques of pulmonary ventilation to help in developing biomarkers. Radiation-induced pulmonary injury results in anatomical changes in form of the increased lung density39–41 that affects its elasticity and its kinematics. The longitudinal measurement of changes of the strain or its anisotropy can reveal asymptomatic tissue damage and forecast its future abnormality. Data can be used to apply clinical intervention to decrease anticipated delayed organ dysfunction. Interfractional measurements might show regional lung response that could warrant radiation treatment replanning. It might be used for correlation of outcome studies with objective characterization of changes within biomechanical framework. These objective characteristics do not rely on human interpretation. The measured changes might have predictive characteristics for the therapeutic success of the treatment. It can be applied to various anatomical regions. Image registration can provide required mapping for comparison studies. Additional advantage is the regional variability of these measurements, e.g., outcome studies of thoracic patients do not take into the consideration heterogeneous functionality of lungs. The global description of the risk of lungs injury cannot provide the map of the interplay between dose and lungs functional variability. For example, Palma et al.42, 43 studied the changes of lung density after stereotactic body radiotherapy and the correlation with the toxicity. Fiorino et al.44 evaluated parotid shrinkage with Jacobians, their histograms and correlated with dose-volume histograms and toxicity. 4. CONCLUSIONS Respiration induced tumor motion does not result in tumor deformation. The averaged over all the tumor voxels strain tensor demonstrates small deformation. The results were obtained with traditional elasticity using the Frobenius norm for strain tensor matrices as well as with much stricter methodology of tensorial processing with the Log-Euclidean distance. The Cauchy-Green strain tensors are positive definite matrices and the Log-Euclidean metrics provides appropriate computational constructs to obtain their means. The anisotropy of the strain tensor was evaluated with factorial and geodesic anisotropy factors. The values were small indicating small anisotropy and indirectly indicating small deformations as well. a) Author
to whom correspondence should be addressed. Electronic mail: [email protected] 1 K. M. Langen and D. T. Jones, “Organ motion and its management,” Int. J. Radiat. Oncol., Biol., Phys. 50, 265–278 (2001). 2 H. H. Liu, P. Balter, T. Tutt, B. Choi, J. Zhang, C. Wang, M. Chi, D. Luo, T. Pan, S. Hunjan, G. Starkschall, I. Rosen, K. Prado, Z. Liao, J. Chang, J. D. Cox, R. Mohan, and L. Dong, “The management of respiratory motion in radiation oncology report of AAPM Task Group 76,” Med. Phys. 33, 3874–3900 (2006). 3 H. H. Liu, P. Balter, T. Tutt, B. Choi, J. Zhang, C. Wang, M. Chi, D. Luo, T. Pan, S. Hunjan, G. Starkschall, I. Rosenand, K. Prado, Z. Liao, J. Chang, R. Komaki, J. D. Cox, R. Mohan, and L. Dong, “Assessing respirationinduced tumor motion and internal target volume using four-dimensional
Medical Physics, Vol. 41, No. 7, July 2014
073503-5
computed tomography for radiotherapy of lung cancer,” Int. J. Radiat. Oncol., Biol., Phys. 68, 531–540 (2007). 4 J. Wu, P. Lei, R. Shekhar, H. Li, M. Suntharalingam, and W. D. D’Souza, “Do tumors in the lung deform during normal respiration? An image registration investigation,” Int. J. Radiat. Oncol., Biol., Phys. 75, 268–275 (2009). 5 G. Holzapfel, Nonlinear Solid Mechanics: A Continuum Approach for Engineering (Wiley, New York, NY, 2000). 6 C. Davatzikos, M. Vaillant, S. M. Resnick, J. L. Prince, S. Letovsky, and R. N. Bryan, “Morphological analysis of brain structures using spatial normalization,” in Proceedings of the 4th International Conference on Visualization in Biomedical Computing, Hamburg, Germany, 1996, Lecture Notes in Computer Science Vol. 1131, edited by K. H. Höhne and R. Kikinis (Springer-Verlag, London, UK, 1996), pp. 355–60. 7 O. Commowick, R. Stefanescu, P. Fillard, V. Arsigny, N. Ayache, X. Pennec, and G. Malandain, “Incorporating statistical measures of anatomical variability in atlas-to-subject registration for conformal brain radiotherapy,” Image Comput. Comput. Assist. Interv. 8, 927–934 (2005). 8 T. C. Doyle and J. L. EricksenL, “Non-linear elasticity,” Adv. Appl. Mech. 4, 53–115 (1956). 9 B. R. Seth, “Generalized strain measures with applications to physical problems,” in Second-order Effects in Elasticity, Plasticity and Fluid Dynamics, edited by M. Reiner and D. Abir (Pergamon, Oxford, 1964), pp. 162–172. 10 V. Arsigny, P. Fillard, X. Pennec, and N. Ayache, “Geometric means in a novel vector space structure on symmetric positive-definite matrices,” SIAM J. Matrix Anal. Appl. 29, 328–347 (2006). 11 V. Arsigny, P. Fillard, X. Pennec, and N. Ayache, “Log-Euclidean metrics for fast and simple calculus on diffusion tensors,” Magn. Reson. Med. 56, 411–421 (2006). 12 P. T. Fletcher, and S. Joshi, “Principal geodesic analysis on symmetric spaces: Statistics of diffusion tensors,” in ECCV Workshops CVAMIA and MMBIA, Lecture Notes in Computer Science Vol. 3117, edited by M. Sonka, I. A. Kakadiaris, and J. Kybic (Springer-Verlag, Berlin, 2004), pp. 87–98. 13 P. G. Batchelor, M. Moakher, D. Atkinson, F. Calamante, and A. Connelly, “A rigorous framework for diffusion tensor calculus,” Magn. Reson. Med. 53, 221–225 (2005). 14 X. Pennec, P. Fillard, and N. Ayache, “A Riemannian framework for tensor computing,” Int. J. Comput. Vis. 66, 41–66 (2006). 15 C. Lenglet, M. Rousson, R. Deriche, and O. Faugeras, “Statistics on the manifold of multivariate normal distributions: Theory and application to diffusion tensor MRI processing,” J. Math. Imaging Vis. 25, 423–444 (2006). 16 J. Gallier, Geometric Methods and Applications: For Computer Science and Engineering, 1st ed., Text in Applied Mathematics Vol. 38 (SpringerVerlag, New York, NY, 2001). 17 P. J. Basser and C. Pierpaoli, “Microstructural and physiological features of tissues elucidated by quantitative-diffusion-tensor MRI,” J. Magn. Reson. B 111, 209–219 (1996). 18 P. T. Fletcher, “Statistical variability in nonlinear spaces: Application to shape analysis and DT-MRI,” Ph.D. thesis (Department of Computer Science, University of North Carolina, 2004). 19 P. T. Fletcher and S. C. Joshi, “Riemannian geometry for the statistical analysis of diffusion tensor data,” Signal Process. 87, 250–262 (2007). 20 J.-P. Thirion, “Non-rigid matching using demons,” in Proceedings of the Computer Vision and Pattern Recognition, CVPR’96, San Francisco, CA, 1996 (IEEE Computer Society, Los Alamitos, California, 1996), pp. 245– 251. 21 J.-P. Thirion, “Image matching as a diffusion process: An analogy with Maxwell’s demons,” Med. Image Anal. 2, 243–260 (1998). 22 X. Pennec, P. Cachier, and N. Ayache, “Understanding the “Demon’s algorithm”: 3D non-rigid registration by gradient descent,” in Proceedings of the 2nd International Conference on Medical Image Computing and Computer-Assisted Intervention (MICCAI’99), Cambridge, UK, 1999, LNCS Vol. 1679, edited by C. Taylor and A. Colchester (Springer, Cambridge, UK, 1999), pp. 597–605. 23 H. Wang, L. Dong, J. O’Daniel, R. Mohan, A. S. Garden, K. K. Ang, D. A. Kuban, M. Bonnen, J. Y. Chang, and R. Cheung, “Validation of an accelerated ‘demons’ algorithm for deformable image registration in radiation therapy,” Phys. Med. Biol. 50, 2887–2905 (2005). 24 H. Wang, L. Dong, M. F. Lii, A. L. Lee, R. de Crevoisier, R. Mohan, J. D. Cox, D. A. Kuban, and R. R. Cheung, “Implementation and
073503-6
Michalski et al.: Strain tensor to estimate thoracic tumors deformation
validation of a three-dimensional deformable registration algorithm for targeted prostate cancer radiotherapy,” Int. J. Radiat. Oncol., Biol., Phys. 61, 725–735 (2005). 25 V. Boldea, D. Sarrut, and C. Carrie, “Comparison of 3D dense deformable registration methods for breath-hold reproducibility study in radiotherapy,” Proc. SPIE 5747, 222–230 (2005). 26 G. C. Sharp, N. Kandasamy, H. Singh, and M. Folkert, “GPU-based streaming architectures for fast cone-beam CT image reconstruction and demons deformable registration,” Phys. Med. Biol. 52, 5771–5783 (2007). 27 S. S. Samant, J. Xia, P. Muyan-Özçelik, and J. D. Owens, “High performance computing for deformable image registration: Towards a new paradigm in adaptive radiotherapy,” Med. Phys. 35, 3546–3553 (2008). 28 S. Nithiananthan, K. K. Brock, M. J. Daly, H. Chan, J. C. Irish, and J. H. Siewerdsen, “Demons deformable registration for CBCT-guided procedures in the head and neck: Convergence and accuracy,” Med. Phys. 36, 4755–4764 (2009). 29 H. Zhong, J. Kim, and I. J. Chetty, “Analysis of deformable image registration accuracy using computational modeling,” Med. Phys. 37, 970–979 (2010). 30 K. K. Brock, “Results of a multi-institution deformable registration accuracy study (MIDRAS),” Int. J. Radiat. Oncol., Biol., Phys. 576, 583–596 (2010). 31 X. Gu, H. Pan, Y. Liang, R. Castillo, D. Yang, D. Choi, E. Castillo, A. Majumdar, T. Guerrero, and S. B. Jiang, “Implementation and evaluation of various demons deformable image registration algorithms on a GPU,” Phys. Med. Biol. 55, 207–219 (2010). 32 E. H. Adelson, C. H. Anderson, J. R. Bergen, P. J. Burt, and J. M. Ogden, “Pyramid methods in image processing,” RCA Eng. 29, 33–41 (1984). 33 T. F. Cootes, S. Marsland, C. J. Twining, K. Smith, and C. J. Taylor, “Groupwise diffeomorphic non-rigid registration for automatic model building,” in ECCV (4), Lecture Notes in Computer Science Vol. 3024, edited by T. Pajdla and J. Matas (Springer-Verlag, Berlin Heidelberg, 2004), pp. 316–327. 34 R. Stefanescu, X. Pennec, and N. Ayache, “Grid powered nonlinear image registration with locally adaptive regularization,” Med. Image Anal. 8, 325–342 (2004).
Medical Physics, Vol. 41, No. 7, July 2014
35 D.
073503-6
Rueckert, P. Aljabar, R. A. Heckemann, J. V. Hajnal, and A. Hammers, “Diffeomorphic registration using b-splines,” in Proceedings of the MICCAI (2)’06, Copenhagen, Denmark, 2006 (Springer-Verlag, Berlin Heidelberg, 2006), pp. 702–709. 36 D. Yang, H. Li, D. A. Low, J. O. Deasy, and I. E. Naqa, “A fast inverse consistent deformable image registration method based on symmetric optical flow computation,” Phys. Med. Biol. 53, 6143–6165 (2008). 37 M. Fan, L. B. Marks, P. Lind, D. Hollis, R. T. Woel, G. G. Bentel, M. S. Anscher, T. D. Shafman, R. Coleman, R. J. Jaszczak, and M. T. Munley, “Relating radiation-induced regional lung injury to changes in pulmonary function tests,” Int. J. Radiat. Oncol., Biol., Phys. 51, 311–317 (2001). 38 K. D. Jaeger, Y. Seppenwoolde, L. J. Boersma, S. H. Muller, P. Baas, J. S. Belderbos, and J. V. Lebesque, “Pulmonary function following highdose radiotherapy of non-small-cell lung cancer,” Int. J. Radiat. Oncol., Biol., Phys. 55, 1331–1340 (2003). 39 Y. Seppenwoolde, K. D. Jaeger, L. J. Boersma, J. S. Belderbos, and J. V. Lebesque, “Regional differences in lung radiosensitivity after radiotherapy for non-small-cell lung cancer,” Int. J. Radiat. Oncol., Biol., Phys. 60, 748–758 (2004). 40 J. Ma, J. Zhang, S. Zhou, J. L. Hubbs, R. J. Foltz, D. R. Hollis, K. L. Light, T. Z. Wong, C. R. Kelsey, and L. B. Marks, “Regional lung density changes after radiation therapy for tumors in and around thorax,” Int. J. Radiat. Oncol., Biol., Phys. 76, 116–122 (2010). 41 E. C. Phernambucq, D. A. Palma, A. Vincent, E. F. Smit, and S. Senan, “Time and dose-related changes in radiological lung density after concurrent chemoradiotherapy for lung cancer,” Lung Cancer 74, 451–456 (2011). 42 D. A. Palma, J. V. S. D. Koste, W. F. A. R. Verbakel, A. Vincent, and S. Senan, “A new approach to quantifying lung damage after stereotactic body radiation therapy,” Acta Oncol. 50, 509–517 (2011). 43 D. A. Palma, J. V. S. D. Koste, W. F. A. R. Verbakel, A. Vincent, and S. Senan, “Lung density changes after stereotactic radiotherapy: A quantitative analysis in 50 patients,” Int. J. Radiat. Oncol., Biol., Phys. 81, 974– 978 (2011). 44 C. Fiorino, E. Maggiulli, S. Broggi, S. Liberini, G. M. Cattaneo, I. Dell’Oca, E. Faggiano, N. D. Muzio, R. Calandrino, and G. Rizzo, “Introducing the Jacobian-volume-histogram of deforming organs: Application to parotid shrinkage evaluation,” Phys. Med. Biol. 56, 3301–3312 (2011).
The use of strain tensor to estimate thoracic tumors deformation Darek Michalski,a) M. Saiful Huq, Greg Bednarz, and Dwight E. Heron Department of Radiation Oncology, University of Pittsburgh Cancer Institute, Pittsburgh, Pennsylvania 15232
(Received 25 June 2013; revised 2 June 2014; accepted for publication 3 June 2014; published 24 June 2014) Purpose: Respiration-induced kinematics of thoracic tumors suggests a simple analogy with elasticity, where a strain tensor is used to characterize the volume of interests. The application of the biomechanical framework allows for the objective determination of tumor characteristics. Methods: Four-dimensional computed tomography provides the snapshots of the patient’s anatomy at the end of inspiration and expiration. Image registration was used to obtain the displacement vector fields and deformation fields, which allows one for the determination of the strain tensor. Its departure from the identity matrix gauges the departure of the medium from rigidity. The tensorial characteristic of each GTV voxel was determined and averaged. To this end, the standard Euclidean matrix norm as well as the Log-Euclidean norm were employed. Tensorial anisotropy was gauged with the fractional anisotropy measure which is based on the normalized variance of the tensors eigenvalues. Anisotropy was also evaluated with the geodesic distance in the Log-Euclidean framework of a given strain tensor to its closest isotropic counterpart. Results: The averaged strain tensor was determined for each of the 15 retrospectively analyzed thoracic GTVs. The amplitude of GTV motion varied from 0.64 to 4.21 with the average of 1.20 cm. The GTV size ranged from 5.16 to 149.99 cc with the average of 43.19 cc. The tensorial analysis shows that deformation is inconsiderable and that the tensorial anisotropy is small. The Log-Euclidean distance of averaged strain tensors from the identity matrix ranged from 0.06 to 0.31 with the average of 0.19. The Frobenius distance from the identity matrix is similar and ranged from 0.06 to 0.35 with the average of 0.21. Their fractional anisotropy ranged from 0.02 to 0.12 with the average of 0.07. Their geodesic anisotropy ranged from 0.03 to 0.16 with the average of 0.09. These values also indicate insignificant deformation. Conclusions: The tensorial framework allows for direct measurements of tissue deformation. It goes beyond the evaluation of deformation via comparison of shapes. It is an independent and objective determination of tissue properties. This methodology can be used to determine possible changes in lung properties due to radiation therapy and possible toxicities. © 2014 American Association of Physicists in Medicine. [http://dx.doi.org/10.1118/1.4884222] Key words: strain tensor, image fusion, image processing, mechanical properties 1. INTRODUCTION Physiologically induced organ motion in thorax introduces certain degree of uncertainty that must be dealt with during the treatment planning and delivery.1, 2 The kinematics of the lungs affect the tumor, which might move during the breathing cycle even up to 5 cm for supradiaphragmatic masses as we observed in our clinic. Tumor motion can occur in all three anatomical directions and its trajectory might exhibit hysteresis. Tumor motion has profound consequences for the entire clinical workflow of the radiation treatment process. The helical computed tomography (CT) can misrepresent tumor shape, location, and size if its motion magnitude exceeds 0.5 cm. This can be adequately addressed via fourdimensional (4D) CT scans. The knowledge of the tumor motion amplitude allows for the accurate determination of margins about the clinical tumor volume (CTV) to define the integral tumor volume (ITV). This weighs in the decision process of the radiation oncologist whether or not elevated doses can be prescribed, or gated treatment should be applied. In the end, it affects the radiation treatment outcomes, which for lung cancer patients do not exhibit expected efficacy. 073503-1
Med. Phys. 41 (7), July 2014
4D CT scanning methodology allows for direct tumor observation over the entire breathing cycle. The AAPM Task Group 76 report2 points out that tumor deformation is usually not quantified or considered in various aspects of motion management in radiation oncology. It is tacitly assumed that the magnitude of relevant margins would address this. There have been two reports in the literature on the lung tumor changes during breathing cycle utilizing image registration for contour comparisons. Liu et al.3 investigated motion of 166 tumors and concluded that tumor motion does not manifest major deformations to tumor contours. They concluded that rigid-body motion suffices to model tumor motion. Similar in methodology, a study of 30 lung tumors4 reported that the comparison of the GTV contours at the end of exhalation (EE) and inhalation (EI) reveals mainly translational change. We are revisiting the topic of lung tumor deformation applying biomechanical framework. We also use image registration to the EE and EI 4D CT volumes but with totally different methodology of quantitative evaluation. Aforementioned two studies relied on comparison of shapes of tumors. The overlap of polygons that built the tumor contours gauged the tumor deformation. We consider the entire volume of the tumor. This
0094-2405/2014/41(7)/073503/6/$30.00
© 2014 Am. Assoc. Phys. Med.
073503-1
073503-2
Michalski et al.: Strain tensor to estimate thoracic tumors deformation
evaluation is direct and does not resort to comparative methods. Deformation is directly quantified by the strain tensor. We utilize scalar and tensorial metrics to evaluate deformation for the entire tumor volume in contrast to the previous studies. This method goes beyond the heuristic of the comparison, which in practice is equivalent to visual inspection. This strain-based approach can be applied to any region of interest since it generically deals with finding a homologous region of interest in images and directly derives all the required parameters to gauge the deformation. We use the lung tumor cases to illustrate this methodology. The deformation is the change of size and/or shape of a given volume of interest. The kinematics of the thorax during breathing can be characterized by its deformation and motion. The same can be said about thoracic gross tumor volume (GTVs). Applying this analogy one may consider to use mathematical constructs inherent to the description of the deformations, which is a strain tensor. In general, this approach allowed to utilize other 4D imaging modalities to evaluate regional strains in heart and in other organs. The temporal resolution of patient’s imaging provides the same topological view of the anatomy but it is also able to provide the insight into cyclic changes of the anatomy during the scan. In other words, this variability can be seen in terms of regional deformations and can be quantified with a strain tensor. 4DCT scans of thoracic patients enable to view physiologically induced motion of the entire thorax. 2. MATERIALS AND METHODS The deformation of a given tissue element can be determined if its displacement is known. The latter can be obtained from 4DCT scans, which were acquired using a GE Discovery ST scanner (GE Medical Systems, Milwaukee, WI) with a real-time position management (RPM) respiratory gating system (Varian Medical Systems, Palo Alto, CA). Retrospective phase-based sorting yields ten temporally equispaced along the breathing cycle CT phases. The breathing cycle starts at EI. All patients during the scanning were immobilized in a supine position with their arms extended above their heads. Patients were coached using an audio directive (“breathe in”). The frequency of the command was matched to a given patient’s normal respiratory period. For the evaluation of tumor deformation, patients whose tumors exhibited at least 0.5 cm motion were selected. Their resampled isotropic 4D CT phases corresponding to EI (CT0) and EE (CT50) were used. The voxel size was the slice thickness. The latter was used as a reference image scan and was registered to the CT0 scan in the Lagrangian framework. A free form deformation image registration provides the displacement field, u = x − X, where X and x are the locations of the homologous tissue elements in CT50 and CT0, respectively. Equivalently we obtain the mapping, , from the CT50 to the CT0, x = (X). The mapping is characterized by its Jacobian matrix, F = ∇. If the mapping is viewed in terms of the coordinate changes, then from the change of variables theorem, the deformed volume element dv is related to the undeformed one dV with dv = J dV. Medical Physics, Vol. 41, No. 7, July 2014
073503-2
Thus, the Jacobian at a given point, J = det(F ) must be positive to guarantee physical feasibility of the transformation.5 The Jacobian is useful since it provides a first insight into the deformation.6 For example, it was used to capture the tissue deformation variability of the group of patients and incorporated in a mathematical construct of a cost function in image registration.7 Yet, even without volume change (J = 1) the medium might exhibit deformation. The parameter that captures this behavior is a strain tensor. Strain tensors might assume various mathematical forms. For a Lagrangian framework F = RU, where R is a proper orthogonal rotation and U a symmetric positive definite matrix. The latter serves in the definition of the Doyle-Ericksen class of strain tensors with the following form E(n) = 1/n (Un − I), where n can be any real number.5, 8, 9 For n = 2, it is the Green-St. Venant strain tensor, E(2) = 1/2 (FT F − Id ), where Id is an identity matrix. Strain tensors are not affected by the rigid transformation so the only change they capture is the deformation. The departure of the strain tensor from the identity matrix gauges the departure of the medium from rigidity. Similar is the right Cauchy-Green strain tensor, computed in this project, C = FT F. It is the identity matrix for the rigid transformation. For each point in the GTV, we can determine its tensorial characteristic and average them to gauge the GTV behavior. Specific mathematical characteristics of the tensorial matrices are not fully congruent with the Euclidean matrix norm as applied in the traditional elasticity, e.g., averaging the tensors C = 1/N i Ci might result in det(Ci ) < det(C) for all det(Ci ). An ingenious Log-Euclidean metric has been proposed10, 11 that addresses these problems. There are other metrics,12–15 which amend the shortcomings of Euclidean norms but the Log-Euclidean is computationally most efficient. This metric determines the distance, d, between tensors T1 and T2 as follows:10, 11 d(T1 , T2 ) = ||ln(T1 ) − ln(T2 )||, where ||. . . || denotes the Euclidean norm. As it was proposed by Arsigny et al.,11 it is the Frobenius norm applied to logarithms of tensors. The Log-Euclidean metric allows 11 for the computation of the average strain tensor over all N voxels of the GTV as C = exp(1/N i ln(Ci )). The logarithm of the tensor in the Log-Euclidean framework is the inverse of the matrix exponential. The logarithm or the exponential of the tensor, T, is obtained by applying the operation to its eigenvalues. These mathematical constructs are well defined for tensors, which are symmetric definite positive matrices.11, 16 Computation of the averaged strain tensor allows for the examination of its anisotropy. The commonly used gauge is a fractional anisotropy,17 FA, which is a scalar value that aims at describing tensor’s anisotropy in terms of the normalized variance of the tensor’s eigenvalues, λi . It is defined √ as FA = ( 3/2 )[ i (λi − λ)2 ]1/2 /[ i (λi )2 ]1/2 where λ = 1/3 i λi . However, another measure of the anisotropy is the geodesic distance in the Log-Euclidean framework of a given strain tensor to its closest isotropic counterpart.13, 18, 19 The latter is the tensor with the same determinant and identical eigenvalues. This gauge was named geodesic anisotropy and can be expressed as GA = [ i (ln(λi ) − lnλ)2 ]1/2 where lnλ=1/3 i ln(λi ).
Michalski et al.: Strain tensor to estimate thoracic tumors deformation
073503-3
0.4
073503-3
0.4 (a)
(b)
0.3
1
0.2
0.2
0.1
2 11
5
Distance
Distance
0.3
8
3
7 13 12
15
6 14 10 9
0.1
4
0 0
50
100
150
0 0
0.5
1
1.5
2
2.5
3
3.5
4
Motion [cm]
Tumor Volume [cc]
F IG . 1. The plot of the Log-Euclidean (䊏) and Euclidean (◦ ) distance of the averaged GTV strain tensors from Id as a function of the GTV volume [cc] (a) and the plot of Log-Euclidean distance for these tensors as a function of the GTV motion amplitude [cm] (b). For the latter, the points are annotated with the patient number ascending with the GTV magnitude.
3. RESULTS AND DISCUSSION The image registration was robust. The GTV motion ranged from 0.64 to 4.21 cm with its average 1.20 cm. The GTV size ranged from 5.16 to 149.99 cc with its average 43.19. Figure 1(a) shows the Log-Euclidean and the Frobenius distance of the averaged strain tensors to the identity matrix. The distances computed with the Log-Euclidean and the traditional matrix norm are very similar. The Log-Euclidean values range from 0.06 to 0.31 with the average 0.19. The Euclidean values range from 0.06 to 0.35 with the average 0.21. The RMSE between the distances computed with these metrics is 0.02. This plot depicts this distance as a function of the GTV magnitude. The values do not exhibit correlation with the GTV size. The Pearson correlation coefficient is −0.04. Since the Log-Euclidean and the Euclidean distances are similar Fig. 1(b) shows only the former as a function of GTV motion. The points are annotated with a GTV number. The GTVs are ordered along ascending GTV magnitude. There is no correlation with GTV motion amplitude as well. Figure 2 shows GA and FA for all the patients numbered in ascending order of the GTV magnitude. The values of GA range from 0.03 to 0.16 with the average 0.09. The values of FA range Medical Physics, Vol. 41, No. 7, July 2014
from 0.02 to 0.12 with the average at 0.07. For all the cases, GA factor is higher. The RMSE for GA and FA is 0.03. There are no correlations of the magnitudes of GA and FA with the size of the GTVs. The Pearson correlation coefficient is 0.06. There is no correlation with the motion amplitude as well. It is instructive to examine these measures using synthetic data. The distance from the identity matrix can be computed using a synthetic tensor with the following eigenvalues λ1 = exp(x), λ2 = λ3 = exp(−x), where x varies between 0 and 1 (λ1 ∈ [1.0, 2.71], and λ2 = λ3 ∈ [1.0, 0.37]). This corresponds to the variation of the Jacobian from 1.0 to about 0.61. Figure 3(a) shows the plot using the Log-Euclidean as well as the Frobenius norms. The plot in Fig. 3(b) shows the results for the following eigenvalues λ1 = exp(x), λ2 = λ3 = exp(−x/2). This plot shows a scenario for a deforming medium that does not change the volume since its Jacobian remains equal to 1 for x ∈ [0, 1]. The difference between the results obtained with the Log-Euclidean and the Frobenius norm in the latter case is more pronounced. The distance from the identity matrix for the isochoric deformation is smaller than for the volumetric deformation since the
0.15
GA/FA
The registration method employed here is based on an iterative and computationally robust Demons algorithm.20–22 It yields voxel displacement field, u, which in turn provides F = Id + ∇u. The algorithm is fast and accurate as reported in the literature. This image registration method was well-tested and it was used by many research groups for various clinical data.23–31 The implementation uses a hierarchical multiscale approach with a three level Gaussian pyramid.32 In order to prevent unphysical folding, the magnitude of the updates was always smaller than 0.4 of the voxel size.33–36 The gradient of the displacement field was computed using analytical form of the derivatives of cubic B-splines. The voxels of the GTVs were extracted using binary mask derived from the structure dicom files. The image registration program was implemented in C on a Linux box.
0.1
0.05
0
0
5
10
15
Patient Number F IG . 2. The plot of GA (䊏) and FA (◦ ) for all patients. Abscissa refers to patients ordinal number. The patients are numbered in the ascending order of the GTV size.
Michalski et al.: Strain tensor to estimate thoracic tumors deformation
073503-4
2
2 (a)
(b) 1.5
Distance
Distance
1.5
1
0.5
0 0
073503-4
1
0.5
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
0 0
1
0.1
0.2
0.3
0.4
x
0.5
0.6
0.7
0.8
0.9
1
x
F IG . 3. The distance of the synthetic C tensor from Id calculated with the Log-Euclidean (——) and the Frobenius (– – –) norms. The abscissa covers the values of x that parametrize the eigenvalues of the tensor; (a) with eigenvalues λ1 = exp(x), λ2 = λ3 = exp(−x) (volumetric deformation) and (b) with λ1 = exp(x), λ2 = λ3 = exp(−x/2) (distortional deformation).
λ2 for this case varies less. For both types of eigenvalues, the Frobenius distance tends to be higher as x increases but for small values of x both distances are similar. The average distance of the GTV strain tensors from the rigidity for the Log-Euclidean metric was 0.19 and for the Frobenius metric 0.21. Both metrics for the strain tensors also yielded the distances of similar magnitude. One can see that these ordinate values in Fig. 3 correspond to the value of the abscissa of about x = 0.1. In this region, both plots also show similar distance for both metrics. The comparison to the synthetic cases suggests that the deformation is in the proximity of the identity matrix. This indicates that the tumor deformation is inconsiderable. Figure 4(a) and 4(b) show, respectively, corresponding plots of GA and FA as calculated with synthetic eigenvalues as described for Figs. 3(a) and 3(b). The difference between GA and FA is less pronounced for the case when the medium deforms but does not change the volume. This is due to the different values of λ2 ( = λ3 ) for the dilational and the distortional case. The synthetic case shows that GA > FA. However, there is no clear relationship between GA and FA. Their magnitudes are also incongruent, FA ∈ [0, 1] and GA ∈ [0, ∞). As for the values of GA and
FA obtained for the averaged strain tensors they indicate small deviation from an isotropic case. In other words, anisotropy is inconsiderable. Since the strain tensor can be multiplicatively ˜ the distortional comdecomposed as C = [(det(C))1/3 Id ] C, ˜ ponent, C, demonstrates an inconsiderable deformation. This component is responsible only for the deformation and does not change the volume. The examination of thoracic tumors within the biomechanical framework agrees with previous findings3, 4 that GTVs demonstrate rather inconsiderable deformation during respiration. In contrast to these studies, we only considered tumors with motion exceeding 0.5 cm. The thoracic tumors are denser than surrounding tissue, which is unable to deform the GTV during respiration. However, theoretically one cannot exclude the existence of the cases where GTV can deform, e.g., if one end of the GTV is attached to the rib cage. In general, the biomechanical framework allows for the quantitative description of the tissue or anatomical regions of interest. Such regional characterization of volume of interests can be used in the objective evaluation of changes of the anatomy during the treatment or after the treatment. As evidenced37, 38 radiation treatment causes toxicity and change
2
2
(a)
(b) 1.5
GA/FA
GA/FA
1.5
1
0.5
0 0
1
0.5
0.1
0.2
0.3
0.4
0.5
x
0.6
0.7
0.8
0.9
1
0 0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
x
F IG . 4. GA (——) and FA (– – –) gauges of a synthetic tensor anisotropy. The abscissa covers the values of x that parametrize the eigenvalues of the tensor; (a) with eigenvalues λ1 = exp(x), λ2 = λ3 = exp(−x) (volumetric deformation) and (b) with λ1 = exp(x), λ2 = λ3 = exp(−x/2) (distortional deformation). Medical Physics, Vol. 41, No. 7, July 2014
073503-5
Michalski et al.: Strain tensor to estimate thoracic tumors deformation
in pulmonary function. For the latter, the tests can only indicate the overall functional change in respiration. However, lung functionality is not homogenous and its map can be revealed by nuclear imaging or other imaging techniques of pulmonary ventilation to help in developing biomarkers. Radiation-induced pulmonary injury results in anatomical changes in form of the increased lung density39–41 that affects its elasticity and its kinematics. The longitudinal measurement of changes of the strain or its anisotropy can reveal asymptomatic tissue damage and forecast its future abnormality. Data can be used to apply clinical intervention to decrease anticipated delayed organ dysfunction. Interfractional measurements might show regional lung response that could warrant radiation treatment replanning. It might be used for correlation of outcome studies with objective characterization of changes within biomechanical framework. These objective characteristics do not rely on human interpretation. The measured changes might have predictive characteristics for the therapeutic success of the treatment. It can be applied to various anatomical regions. Image registration can provide required mapping for comparison studies. Additional advantage is the regional variability of these measurements, e.g., outcome studies of thoracic patients do not take into the consideration heterogeneous functionality of lungs. The global description of the risk of lungs injury cannot provide the map of the interplay between dose and lungs functional variability. For example, Palma et al.42, 43 studied the changes of lung density after stereotactic body radiotherapy and the correlation with the toxicity. Fiorino et al.44 evaluated parotid shrinkage with Jacobians, their histograms and correlated with dose-volume histograms and toxicity. 4. CONCLUSIONS Respiration induced tumor motion does not result in tumor deformation. The averaged over all the tumor voxels strain tensor demonstrates small deformation. The results were obtained with traditional elasticity using the Frobenius norm for strain tensor matrices as well as with much stricter methodology of tensorial processing with the Log-Euclidean distance. The Cauchy-Green strain tensors are positive definite matrices and the Log-Euclidean metrics provides appropriate computational constructs to obtain their means. The anisotropy of the strain tensor was evaluated with factorial and geodesic anisotropy factors. The values were small indicating small anisotropy and indirectly indicating small deformations as well. a) Author
to whom correspondence should be addressed. Electronic mail: [email protected] 1 K. M. Langen and D. T. Jones, “Organ motion and its management,” Int. J. Radiat. Oncol., Biol., Phys. 50, 265–278 (2001). 2 H. H. Liu, P. Balter, T. Tutt, B. Choi, J. Zhang, C. Wang, M. Chi, D. Luo, T. Pan, S. Hunjan, G. Starkschall, I. Rosen, K. Prado, Z. Liao, J. Chang, J. D. Cox, R. Mohan, and L. Dong, “The management of respiratory motion in radiation oncology report of AAPM Task Group 76,” Med. Phys. 33, 3874–3900 (2006). 3 H. H. Liu, P. Balter, T. Tutt, B. Choi, J. Zhang, C. Wang, M. Chi, D. Luo, T. Pan, S. Hunjan, G. Starkschall, I. Rosenand, K. Prado, Z. Liao, J. Chang, R. Komaki, J. D. Cox, R. Mohan, and L. Dong, “Assessing respirationinduced tumor motion and internal target volume using four-dimensional
Medical Physics, Vol. 41, No. 7, July 2014
073503-5
computed tomography for radiotherapy of lung cancer,” Int. J. Radiat. Oncol., Biol., Phys. 68, 531–540 (2007). 4 J. Wu, P. Lei, R. Shekhar, H. Li, M. Suntharalingam, and W. D. D’Souza, “Do tumors in the lung deform during normal respiration? An image registration investigation,” Int. J. Radiat. Oncol., Biol., Phys. 75, 268–275 (2009). 5 G. Holzapfel, Nonlinear Solid Mechanics: A Continuum Approach for Engineering (Wiley, New York, NY, 2000). 6 C. Davatzikos, M. Vaillant, S. M. Resnick, J. L. Prince, S. Letovsky, and R. N. Bryan, “Morphological analysis of brain structures using spatial normalization,” in Proceedings of the 4th International Conference on Visualization in Biomedical Computing, Hamburg, Germany, 1996, Lecture Notes in Computer Science Vol. 1131, edited by K. H. Höhne and R. Kikinis (Springer-Verlag, London, UK, 1996), pp. 355–60. 7 O. Commowick, R. Stefanescu, P. Fillard, V. Arsigny, N. Ayache, X. Pennec, and G. Malandain, “Incorporating statistical measures of anatomical variability in atlas-to-subject registration for conformal brain radiotherapy,” Image Comput. Comput. Assist. Interv. 8, 927–934 (2005). 8 T. C. Doyle and J. L. EricksenL, “Non-linear elasticity,” Adv. Appl. Mech. 4, 53–115 (1956). 9 B. R. Seth, “Generalized strain measures with applications to physical problems,” in Second-order Effects in Elasticity, Plasticity and Fluid Dynamics, edited by M. Reiner and D. Abir (Pergamon, Oxford, 1964), pp. 162–172. 10 V. Arsigny, P. Fillard, X. Pennec, and N. Ayache, “Geometric means in a novel vector space structure on symmetric positive-definite matrices,” SIAM J. Matrix Anal. Appl. 29, 328–347 (2006). 11 V. Arsigny, P. Fillard, X. Pennec, and N. Ayache, “Log-Euclidean metrics for fast and simple calculus on diffusion tensors,” Magn. Reson. Med. 56, 411–421 (2006). 12 P. T. Fletcher, and S. Joshi, “Principal geodesic analysis on symmetric spaces: Statistics of diffusion tensors,” in ECCV Workshops CVAMIA and MMBIA, Lecture Notes in Computer Science Vol. 3117, edited by M. Sonka, I. A. Kakadiaris, and J. Kybic (Springer-Verlag, Berlin, 2004), pp. 87–98. 13 P. G. Batchelor, M. Moakher, D. Atkinson, F. Calamante, and A. Connelly, “A rigorous framework for diffusion tensor calculus,” Magn. Reson. Med. 53, 221–225 (2005). 14 X. Pennec, P. Fillard, and N. Ayache, “A Riemannian framework for tensor computing,” Int. J. Comput. Vis. 66, 41–66 (2006). 15 C. Lenglet, M. Rousson, R. Deriche, and O. Faugeras, “Statistics on the manifold of multivariate normal distributions: Theory and application to diffusion tensor MRI processing,” J. Math. Imaging Vis. 25, 423–444 (2006). 16 J. Gallier, Geometric Methods and Applications: For Computer Science and Engineering, 1st ed., Text in Applied Mathematics Vol. 38 (SpringerVerlag, New York, NY, 2001). 17 P. J. Basser and C. Pierpaoli, “Microstructural and physiological features of tissues elucidated by quantitative-diffusion-tensor MRI,” J. Magn. Reson. B 111, 209–219 (1996). 18 P. T. Fletcher, “Statistical variability in nonlinear spaces: Application to shape analysis and DT-MRI,” Ph.D. thesis (Department of Computer Science, University of North Carolina, 2004). 19 P. T. Fletcher and S. C. Joshi, “Riemannian geometry for the statistical analysis of diffusion tensor data,” Signal Process. 87, 250–262 (2007). 20 J.-P. Thirion, “Non-rigid matching using demons,” in Proceedings of the Computer Vision and Pattern Recognition, CVPR’96, San Francisco, CA, 1996 (IEEE Computer Society, Los Alamitos, California, 1996), pp. 245– 251. 21 J.-P. Thirion, “Image matching as a diffusion process: An analogy with Maxwell’s demons,” Med. Image Anal. 2, 243–260 (1998). 22 X. Pennec, P. Cachier, and N. Ayache, “Understanding the “Demon’s algorithm”: 3D non-rigid registration by gradient descent,” in Proceedings of the 2nd International Conference on Medical Image Computing and Computer-Assisted Intervention (MICCAI’99), Cambridge, UK, 1999, LNCS Vol. 1679, edited by C. Taylor and A. Colchester (Springer, Cambridge, UK, 1999), pp. 597–605. 23 H. Wang, L. Dong, J. O’Daniel, R. Mohan, A. S. Garden, K. K. Ang, D. A. Kuban, M. Bonnen, J. Y. Chang, and R. Cheung, “Validation of an accelerated ‘demons’ algorithm for deformable image registration in radiation therapy,” Phys. Med. Biol. 50, 2887–2905 (2005). 24 H. Wang, L. Dong, M. F. Lii, A. L. Lee, R. de Crevoisier, R. Mohan, J. D. Cox, D. A. Kuban, and R. R. Cheung, “Implementation and
073503-6
Michalski et al.: Strain tensor to estimate thoracic tumors deformation
validation of a three-dimensional deformable registration algorithm for targeted prostate cancer radiotherapy,” Int. J. Radiat. Oncol., Biol., Phys. 61, 725–735 (2005). 25 V. Boldea, D. Sarrut, and C. Carrie, “Comparison of 3D dense deformable registration methods for breath-hold reproducibility study in radiotherapy,” Proc. SPIE 5747, 222–230 (2005). 26 G. C. Sharp, N. Kandasamy, H. Singh, and M. Folkert, “GPU-based streaming architectures for fast cone-beam CT image reconstruction and demons deformable registration,” Phys. Med. Biol. 52, 5771–5783 (2007). 27 S. S. Samant, J. Xia, P. Muyan-Özçelik, and J. D. Owens, “High performance computing for deformable image registration: Towards a new paradigm in adaptive radiotherapy,” Med. Phys. 35, 3546–3553 (2008). 28 S. Nithiananthan, K. K. Brock, M. J. Daly, H. Chan, J. C. Irish, and J. H. Siewerdsen, “Demons deformable registration for CBCT-guided procedures in the head and neck: Convergence and accuracy,” Med. Phys. 36, 4755–4764 (2009). 29 H. Zhong, J. Kim, and I. J. Chetty, “Analysis of deformable image registration accuracy using computational modeling,” Med. Phys. 37, 970–979 (2010). 30 K. K. Brock, “Results of a multi-institution deformable registration accuracy study (MIDRAS),” Int. J. Radiat. Oncol., Biol., Phys. 576, 583–596 (2010). 31 X. Gu, H. Pan, Y. Liang, R. Castillo, D. Yang, D. Choi, E. Castillo, A. Majumdar, T. Guerrero, and S. B. Jiang, “Implementation and evaluation of various demons deformable image registration algorithms on a GPU,” Phys. Med. Biol. 55, 207–219 (2010). 32 E. H. Adelson, C. H. Anderson, J. R. Bergen, P. J. Burt, and J. M. Ogden, “Pyramid methods in image processing,” RCA Eng. 29, 33–41 (1984). 33 T. F. Cootes, S. Marsland, C. J. Twining, K. Smith, and C. J. Taylor, “Groupwise diffeomorphic non-rigid registration for automatic model building,” in ECCV (4), Lecture Notes in Computer Science Vol. 3024, edited by T. Pajdla and J. Matas (Springer-Verlag, Berlin Heidelberg, 2004), pp. 316–327. 34 R. Stefanescu, X. Pennec, and N. Ayache, “Grid powered nonlinear image registration with locally adaptive regularization,” Med. Image Anal. 8, 325–342 (2004).
Medical Physics, Vol. 41, No. 7, July 2014
35 D.
073503-6
Rueckert, P. Aljabar, R. A. Heckemann, J. V. Hajnal, and A. Hammers, “Diffeomorphic registration using b-splines,” in Proceedings of the MICCAI (2)’06, Copenhagen, Denmark, 2006 (Springer-Verlag, Berlin Heidelberg, 2006), pp. 702–709. 36 D. Yang, H. Li, D. A. Low, J. O. Deasy, and I. E. Naqa, “A fast inverse consistent deformable image registration method based on symmetric optical flow computation,” Phys. Med. Biol. 53, 6143–6165 (2008). 37 M. Fan, L. B. Marks, P. Lind, D. Hollis, R. T. Woel, G. G. Bentel, M. S. Anscher, T. D. Shafman, R. Coleman, R. J. Jaszczak, and M. T. Munley, “Relating radiation-induced regional lung injury to changes in pulmonary function tests,” Int. J. Radiat. Oncol., Biol., Phys. 51, 311–317 (2001). 38 K. D. Jaeger, Y. Seppenwoolde, L. J. Boersma, S. H. Muller, P. Baas, J. S. Belderbos, and J. V. Lebesque, “Pulmonary function following highdose radiotherapy of non-small-cell lung cancer,” Int. J. Radiat. Oncol., Biol., Phys. 55, 1331–1340 (2003). 39 Y. Seppenwoolde, K. D. Jaeger, L. J. Boersma, J. S. Belderbos, and J. V. Lebesque, “Regional differences in lung radiosensitivity after radiotherapy for non-small-cell lung cancer,” Int. J. Radiat. Oncol., Biol., Phys. 60, 748–758 (2004). 40 J. Ma, J. Zhang, S. Zhou, J. L. Hubbs, R. J. Foltz, D. R. Hollis, K. L. Light, T. Z. Wong, C. R. Kelsey, and L. B. Marks, “Regional lung density changes after radiation therapy for tumors in and around thorax,” Int. J. Radiat. Oncol., Biol., Phys. 76, 116–122 (2010). 41 E. C. Phernambucq, D. A. Palma, A. Vincent, E. F. Smit, and S. Senan, “Time and dose-related changes in radiological lung density after concurrent chemoradiotherapy for lung cancer,” Lung Cancer 74, 451–456 (2011). 42 D. A. Palma, J. V. S. D. Koste, W. F. A. R. Verbakel, A. Vincent, and S. Senan, “A new approach to quantifying lung damage after stereotactic body radiation therapy,” Acta Oncol. 50, 509–517 (2011). 43 D. A. Palma, J. V. S. D. Koste, W. F. A. R. Verbakel, A. Vincent, and S. Senan, “Lung density changes after stereotactic radiotherapy: A quantitative analysis in 50 patients,” Int. J. Radiat. Oncol., Biol., Phys. 81, 974– 978 (2011). 44 C. Fiorino, E. Maggiulli, S. Broggi, S. Liberini, G. M. Cattaneo, I. Dell’Oca, E. Faggiano, N. D. Muzio, R. Calandrino, and G. Rizzo, “Introducing the Jacobian-volume-histogram of deforming organs: Application to parotid shrinkage evaluation,” Phys. Med. Biol. 56, 3301–3312 (2011).
