Informatica 39 (2015) 261-269 261 Cervix Cancer Spatial Modelling for Brachytherapy Applicator Analysis Peter Rogelj University of Primorska, Faculty of Mathematics, Natural Sciences and Information Technologies Glagoljaška 8, 6000 Koper, Slovenia E-mail: peter.rogelj@upr.si Muhamed Barakovic University of Verona, Department of Biotechnology Ca' Vignal 1, Strada Le Grazie 15, 37134 Verona, Italia E-mail: muhamed.barakovic@studenti.univr.it Keywords: cervix cancer, spatial distribution, PCA, BT applicators Received: June 18, 2014 SStandard applicators for cervix cancer brachytherapy (BT) do not always enable a sufficient radiation dose coverage of the target structure (HR-CTV). The aim of this study was to develop methodology for building models oftheBT target from a cohort of cervix cancer patients, which would enable BT applicator testing. In this paper we propose two model types, a spatial distribution model and a principal component model. Each of them can be built from data of several patients that includes medical images of arbitrary resolution and modality supplemented with delineations of HR-CTV structure, reconstructed applicator structure and eventual organs at risk (OAR) structures. The spatial distribution model is a static model providing probability distribution of the target in the applicator coordinate system, and as such provides information of the target region that applicators must be able to cover. The principal component model provides information of the target spatial variability described by only a few parameters. It can be used to predict specific extreme situations in the scope of sufficient applicator radiation dose coverage in the target structure as well as radiation dose avoidance in OARs. The results are generated 3D images that can be imported into existent BT planning systems for further BT applicator analysis and eventual improvements. Povzetek: Razvita sta dva modela za izboljšanje brahoterapije. 1 Introduction Applicators for cervix cancer brachytherapy (BT) enable cancer treatment that in comparison with external beam radiotherapy (EBRT) provides better radiation coverage of the high risk clinical target volume (HR-CTV) and better avoidance of organs at risk [1]. During the last decade remarkable progress has been made in radiotherapy, including cervix cancer BT [2]. Standard BT applicators for cervix cancer, as shown in Fig. 1, however still do not always enable a sufficient radiation dose coverage of the target, especially in cases of locally advanced cervical cancer. Improvements are searched in the direction of incorporating additional application needles. A development of new applicators that would enable better target dose coverage requires knowledge of cervix cancer spatial distribution and variation. Furthermore, as the applicators should be able to avoid organs at risk, the information of their variability should also be taken into account. In this work we aimed to develop methodology to obtain this information statistically using available data of past and present cervix cancer patients. The information required from each patient includes BT planning medical image, delineated HR-CTV structure, reconstructed BT applicator structure, and organs at risk (OAR) structures. HR-CTV and OAR structures are in each 3D image delineated on each image slice, wherever the specific structure is present and, thus, available as a set of closed planar contours. BT applicators are reconstructed such that an applicator model is positioned in the 3D image inside the BT planning system. The applicator models consist of a ring, applicator tandem and needles, which are reconstructed independently. The actual position of the applicator is evident from the position of the applicator ring. Because the applicator must be positioned directly to the cervix and because the purpose of the models is applicator analysis, the spatial distribution and variation must be defined in the applicator coordinate system. The significance of tumor distribution depend on the tumor type. It can help in development of tumor treatment and biopsy strategies and techniques [3, 4]. In the case of the cervical cancer it is also important due to BT applicator design. Representation of 3D structures by sets of closed planar contours is not convenient for further spatial analysis. Other representations can be used instead, e.g., by tensors [5], Gaussian random spheres [6], signed distance maps [7] and others. However, due to eventual high complexity of BT target structures, we have selected the most common 262 Informatica 39 (2015) 261-269 P. Rogelj et al. Figure 1: An image of a standard BT cervix cancer applicator with indicated parts: tandem, ring and optional needles. representation of structures by binary images. In the following sections our approach to model the spatial configuration of cervix cancer is described first. We describe the proposed methods of constructing the spatial distribution model and the principal component model. Then we show some test results; for the spatial distribution model based on real patient data, while the principal component model is illustrated using synthetic data. We conclude with discussion that includes the analysis of provided benefits and limitations. 2 Methods Our approach to build spatial models of cervix cancer consists of the following processing steps that are described below: data input, applicator coordinate system definition, structure processing, modelling and data export. 2.1 Data input The input data for building the models consists of patient medical data sets that comprise all the required information of each patient, i.e., a 3D medical image, delineated HR-CTV structure, OAR structures and a reconstructed applicator structure. This data is typically provided in the DI-COM file format, which can be imported using DICOM libraries, e.g., the GDCM library [8], or by Matlab using Image processing toolbox. Medical images are needed only to obtain the image configuration, i.e., transformation of image coordinate system according to the patient coordinate system and image slice positions, which are required to correctly interpret the structures. Structures are given in a form of structure sets that include all the structure data required for BT treatment. The target structure (HR-CTV), OAR structures and the applicator ring structure can be identified from all the structures according to their names that must be known in advance for each individual data set; the structure naming is not standardized. The target and OAR Figure 2: The applicator coordinate system is defined according to the applicator ring structure contour, with origin in the applicator ring center defined by the last point of the contour A(N), xy plane in the ring plane with x axis pointing towards the contour starting point A(1) and z axis in the direction of the tandem. structures are defined as sets of closed planar contours, i.e., each contour is positioned on one slice of the corresponding 3D image. The applicator ring structures are described with a single open nonplanar contour. All contours are defined in the patient coordinate system. 2.2 Applicator coordinate system definition Because the spatial configuration of cervix cancer needs to be defined according to the applicator perspective, an applicator coordinate system needs to be defined. The applicator reconstruction [9, 10] is performed on radiotherapy planning systems by importing predefined geometry structures. The applicator consists of tandem, ring and eventual additional needles, see Figure 1, which are all reconstructed independently. The ring structure, when inserted, tightly fits to the cervix anatomy, and provides a good base for defining the applicator coordinate system. Different applicator types may have different ring diameter, may be described with different number of contour points, however in practice the point ordering is always the same. For the illustration see Fig. 2. We propose that the applicator coordinate system is defined with origin in the ring center (the last point of the contour), xy plane in the ring plane, x axis in the direction towards ring contour starting point and z axis in the direction of the tandem. The transformation that defines the applicator coordinate system according to the patient coordinate system can be for each applicator type computed from its ring contour coordinates. Actually, only three noncollinear contour points are required to compute the applicator coordinate system transformation TA, i.e., the last point, A(N) that is positioned at the ring center as the coordinate system origin, A(1) that defines the applicator coordinate system x axis, and any other point in the applicator contour circumference, A(M) for defining the xy plane. The procedure is Cervix Cancer Spatial Modelling for... Informatica 39 (2015)261-269 263 the following: Oa = A(N) (1) V1 = A(1) - Oa (2) V2 = A(M) - Oa (3) Vz Vi x V2 II V x V2II (4) Vy Vz x Vi II Vz x Vi I (5) Vx Vy X Vz II Vy X Vz I (6) T A = Vx(x) V*(y) V*(Z) 0 Vy (x) Vy (y) Vy (z) 0 Vz (x) Oa (x) Vz (y) OA(y) Vz (z) Oa(z) 0 1 (7) where, V and O are three dimensional vectors with components x, y, and z, such that OA represents applicator coordinate system origin while Vx, Vy, and Vz are applicator coordinate system axes. Vector (cross) products assured coordinate axis perpendicularity, defining a Cartesian coordinate system. The obtained transformation matrix TA is needed for transforming BT structures to the applicator coordinate system in which the cervix cancer needs to be modelled. 2.3 Structure processing For each image the corresponding BT target structure and OAR structures must be mapped into the applicator coordinate system. These structures are created by drawing contours on individual image slices and are provided as point sequences in the patient coordinate system. Such vector definition of structures is difficult to process statistically in coordinate system that is not parallel to the coordinate system of the originating image. Our solution is to present the structures in bitmap instead of vector format and process them as 3D (binary) images with voxel values 1 representing regions inside structures and 0 representing the surrounding. The approach is illustrated in Figure 3. The binary images cover the same region as the original medical image, except that, they may have different resolution in x and y image direction to control discretization error and data size. The resolution in z image direction must remain unchanged in order to preserve location of slices on which contours are defined. The process of converting certain structure into a binary image starts with mapping the structure to the coordinate system of the original image. The transformation T/, form patient to image coordinate system can be obtained from image meta information, i.e., from DICOM tags Image Position Patient (0020,0032) and Image Orientation Patient (0020,0037). Thus, each contour C of structure S gets defined in its image coordinate system as C/: All points of the the same contour gets equal image coordinate z/ that is equivalent to the position of the image slice on which the contour was defined. The obtained structure can, as such, be drawn to the binary image, contour by contour. The process initiates by initializing all voxel values of the binary structure image to 0, followed by drawing the contours by checking each slice voxel if it is positioned inside of a polygon of contour points. Voxels inside the polygon gets negated to correctly interpret even complex structure shapes, e.g., shapes that include holes. Binary structure images enable further data integration towards spatial cervix cancer models. To integrate the structure binary images of all patients into a spatial model, they all need to be mapped into the common applicator coordinate system (A), because patient coordinate systems (P) and image coordinate systems (I) are specific for each study/patient. Transformation between the coordinate systems are illustrated in Fig. 4. The data defined in image coordinate system (I) can be transformed to the applicator coordinate system (A) through the patient coordinate system (P) using transformation T/A: T IA T-1 Ti . (9) Structure binary images do not differ only according to their coordinate systems, but also according to image size and voxel size (resolution). For further analysis they need to be unified. The target region of interest and required precision define configuration of the resulting model (image) size and voxel size. All binary images must be resampled into this common spatial configuration. We recommend resampling by linear interpolation in reverse direction such that intensity corresponding to each voxel in the model configuration is interpolated from voxel intensities in the binary image. Note, that linear interpolation transforms a binary image into an image with real voxel values in interval [0,1]. The result of structure processing is, therefore, a set of structure images SA in a common applicator coordinate system and with common size and voxel size, i.e., each structure of each patient results in one structure image with common applicator (ring) position: Sa = mterp(T-A S), (10) Where S denotes a structure binary image in the coordi-mate system of the original image, interp a linear image interpolation, and SA a structure image in the applicator coordinate system. 2.4 Spatial distribution model The purpose of the spatial distribution model is to provide an overview of BT target spatial extent. It is given in a form of a spatial distribution image D, i.e., an image of the region of interest whose voxel values denote probability of voxel being inside the BT target region. It is obtained from images of the HR-CTV structure by averaging: Ci T-1C. (8) D =1 p n ^ p= 1 Sa,p,hr-ctv , (11) 264 Informatica 39 (2015) 261-269 P. Rogelj et al. T .-1 structure drawing Figure 3: Illustration of structure processing: first, contours provided in the patient coordinate system (left) are transformed to the image coordinate system (center). Then, contours are drawn on image slices in 2D, which resulted in a 3D binary image of the structure (right). Figure 4: Illustration of patient (P), image (I) and applicator (A) coordinate systems and their transformations: T= T-1T/. where Sa,p,hr-ctv represents HR-CTV structure data of p-th patient resampled into an applicator coordinate system, and P is a total number of patients included in the analysis. 2.5 Principal component model The principal component model provides information of the BT target spatial variability expressed by only a small number of parameters. The general idea is to be able to reconstruct any target configuration, i.e., position and extent of HR-CTV as well as OAR structures, by correctly setting the model parameters. As such the principal component model can be used to predict various target configurations, e.g., extreme situations in the scope of sufficient applicator radiation dose coverage in the target structure as well as radiation dose avoidance in OARs. Such situations may be crucial for testing real applicator efficiency. The principal component model tends to extract a minimal set of orthogonal components of spatial variations in the region of interest using the principal component analysis (PCA). PCA projects the data into a lower dimensional linear space such that the variance of the projected data is maximized, or equivalently, it is the linear projection that minimizes the mean squared distance between the data points and their projections. PCA provides a full set of components that enable perfect data reconstruction, however, it also orders the components according to their importance, i.e., according to their contribution to the data description. It turns out that majority of the components have low importance and only a small error is made when only a few most important components are used. In this case the important components can be computed more efficiently using singular value docomposition (SVD) [11]. Our input data for the PCA analysis of the BT target are the HR-CTV structure images in the applicator coordinate system SA,p,HR-CTV. The data of each image is reordered into a row vector and joined for all the patients into a matrix XPxL, with L being the number of pixels in the image. Then the mean vector X is computed and subtracted from each data row to obtain the matrix X0 representing the zero-mean data variation. Here, the mean vector X corresponds to the reordered data of the spatial distribution model D. SVD decomposes X0 into three matrices; matrix V with orthogonal columns that represent principal components, diagonal matrix S with singular values that represent importance of the components, and matrix U providing component weights for reconstructing the input data: X0 = USVT (12) The efficient SVD implementations, e.g., Matlab svds function, enable computation of only a given number of principal components R, and as such provide approximate solutions: X0 P xi ~ UP xRSRxR VLxR (13) The obtained matrices S and V represent a principal component model of the HR-CTV structure, such that HR-CTV structure of any patient can be represented with the R components, i.e., the columns of V, with weights: U' = X0VS-1 (14) where X'0 = X' - X represents deviation of the data from the average. Similarly, BT target data can also be simulated Cervix Cancer Spatial Modelling for... Informatica 39 (2015)261-269 265 by manual setting of component weights in U, following equation (13) and adding the mean vector X. Component weights form a low dimensional linear space with a certain region around the origin that corresponds to realistic data variation. The limits of this realistic subspace can be estimated by analyzing large amount of data, i.e., large number of patients. Values at the border of the realistic subspace can be used in U to simulate specific extreme situations suitable for BT applicator analysis. Testing of BT applicators on their ability to radiate the HR-CTV regions may be biased, as applicators should also be able to avoid radiation of OAR structures. It is important that HR-CTV and OAR structures cannot overlap. This property can be used to simultaneously model both structure types , i.e., HR-CTV as well as OAR structures, without increasing the amount of data in the PCA analysis. For this purpose the input vector X must be constructed from all the structure images, such that positive values represent target regions and negative values represent OAR regions: X = Xhr-ctv — Xoar (15) Here, XHR-CTV and XOAR are constructed from structure images with reordering into row vectors as described earlier, using HR-CTV and available OAR structures. Typically, OAR structures include bladder, rectum and sigmoid colon. The simulated or reconstructed data that results from the principal component model, as well as principal components themselves, can be reordered back into 3D images. Due to interpolations and approximations the reconstructed structures are not presented only with values 1,-1 and 0 for target structures, OAR structures and surrounding respectively. Consequently, we recommend completing the reconstruction procedure with thresholding using thresholds -0.5 and 0.5. 2.6 Data export The resulting images, i.e., an image of the spatial distribution model and simulated or reconstructed BT target configurations can be used in BT planning systems, e.g., Brachyvision®, for further analysis. BT planning systems include functionalities that enable radiation simulations using different radiation plans and can be used to test the efficiency of different applicators. To enable this procedures the images shall be exported to DICOM image format, which can be done using DICOM libraries, e.g. GDCM library [8], or Matlab image processing toolbox. 3 Results We have tested the proposed methods on real and simulated data. First a spatial distribution model has been created from real data of 264 consecutive cervix cancer patients. Due to relatively large number of patients, the obtained estimate of spatial cervix cancer distribution was named a Figure 5: Illustration of the cervix cancer spatial distribution representing a virtual patient, the central coronal slice. virtual patient (VP). Imported to the BT planning system isosurfaces that connect voxels with the same values were created and labeled as percentage of encompassed voxels. VPn was defined as VP subvolume, encompassed by the n% isosurface, see the illustration in Figure 5. The obtained VP data was used for analysis and development of BT applicators for cervix cancer [12]. It was found out, that standard tandem and ring (T&R) applicator enables adequate treatment of VP60 subvolume, additional needles parallel to tandem extend adequate treatment to VP95 and additional oblique needles, inserted at points, angles and depths extend adequate treatment to VP99 subvolume. The principal component model was not built for this dataset, such that applicators were tested only for their general capability of radiating the HR-CTV, not considering the capability of avoiding the radiation of OAR structures. The principal component model was tested using a simulated dataset that we have created for this purpose. Note that the simulated structure images presented here do not realistically simulate the BT target configuration, however enable illustration of the concept and testing of its suitability for creating a realistic model. The simulated data was generated using four random parameters where three of them were used to simulate variability of the HR-CTV structure and the additional one for the variability of one OAR structure. The HR-CTV structure was simulated as an ellipsoid with the three parameters representing the semi-axes lengths while its center was always in the applicator coordinate system origin. The OAR region was simulated as a sphere with the given parameter representing its radius, while its center was defined such that the distance between the edges of BT target and the OAR structure was constant. For the illustration see Figure 6. A principal component model was generated from a 266 Informatica 39 (2015) 261-269 P. Rogelj et al. Figure 6: Illustration of the simulated dataset configuration. The HR-CTV was simulated with an ellipsoid and one additional OAR structure as a sphere with a constant distance d from the HR-CTV. dataset of 400 simulated 3D images with 100 x 100 x 50 voxels. The simulation parameters were selected randomly in the following ranges: a e [40, 73], b e [35, 59], c G [35,49], r e [15,20] and d =5. The computation of the principal component model was restricted to 11 principal components. The mean image X and the components are illustrated in Figure 7. The singular values that represent the distribution of the dataset's energy among the principal components indicate that the component energy gradually decreases with the component number, see Figure 8. However, although not all of the energy was considered, the reconstructed images did not differ considerably from the images from the training set as shown in Figure 9, where a randomly selected input structure image is compared with its reconstructed approximations obtained using three and eleven principal components. We can notice minor differences even when reconstructing from three components only. If we observe the component weights (the values of matrix U), we can see that they are spread over a limited PCA subspace, see Figure 10, which corresponds to valid structure images. According to the shape of the subspace, we can conclude that component weights of valid images are not fully independent, although the components are orthogonal. By selecting weights manually, additional structure images can be simulated. If the selected weights are from the subspace of valid structure images, the simulated images follow the concepts of the input dataset, else the results may include major deviations as demonstrated in Figure 11. The possibility to simulate structure images and have control over its validity offers good opportunity to generate specific synthetic images of the BT target region that represent extreme situations for BT applicator testing. In that case the principal component model should be created from real patient data and the test cases selected at the border of the populated PCA subregion. The realistic model has not been created, yet, however we are looking forward to create it in collaboration with medical institutions that maintain large databases of their cervix cancer patients. 4 Discussion and conclusion Cancer spatial distributions must be considered whenever cancer treatment tools and procedures are being developed. Unfortunately, statistical analysis of spatial distributions related to specific organs is in general tedious due to difficulties defining the reference coordinate systems because of their complex shapes and their high variability. In the specific case of cervix cancer the organ geometry enables unambiguous coordinate system definition that agrees with the applicator ring structure. Analysis of other cancer types would require definition of analysis coordinate system according to organ geometry and data integration performed by image registration with a reference or atlas [13]. Similarly, image registration has already been used for analyzing interfraction variation of high dose regions of OARs [14], and could be extended to intersubject analysis of cancer distributions. The spatial distribution model provides useful information about target region that needs to be radiated, and has already been used for development of novel applicator types [12]. However, this model does not consider OARs and difficulties of restricting radiation dose in these structures. If a distribution model would be made for OARs as well, it would most probably overlap with the cancer distribution model due to closeness of some OAR structures to HR-CTV and due to anatomical variability. Better applicator testing must, therefore, take into account the BT target variability, e.g., by testing on diverse specific target configurations, which can correspond to real patients or obtained by modelling. The proposed principal component model has advantages over using the real patients' data, because of the established control over the specificity of the cases, a possibility to simulate the non-existent cases and deper-sonalization. The limitation of the principal model is in its high computational cost. Computation of all the PCA components would require enormous amount of memory, only the V matrix would have the size of 500k x 500k elements (assuming 2 x 2 x 2 mm voxel size), which in float data format requires 1TB of memory. Using the SVD approach with computation of the most important components only, drastically reduces the memory requirements; in our simulated case matrix V occupied only 22MB. Such reduction of components is possible due to final thresholding, which is applicable due to binary nature of the structures. When the computational cost remains a problem, high efficient PCA solutions [15] or alternative structure representations could be used. A principal component model of real cervix cancer has not been made, yet. A large number of patient datasets is required and in contrast to the spatial distribution model the OAR structures must be included. The preparation of such data is tedious due to non-standardized structure naming. Cervix Cancer Spatial Modelling for... Informatica 39 (2015)261-269 267 X component 1 component 2 component 3 Figure 7: Components of the simulated dataset (central slices only). The mean image X is presented in a scale from -1 (black) to +1 (white) and components with a scale from -0.01 (black) to +0.01 (white). • • Figure 8: Singular values corresponding to the first 11 components; singular values represent the distribution of the dataset's energy. Figure 9: The central slice of an input structure image (top) and its reconstruction using 3 and 11 components (bottom left and right respectively). 1400 800 300 400 200 10 component 268 Informatica 39 (2015) 261-269 P. Rogelj et al. 0.1 -0.080.060.04-| 002- IP 0£ -0.02-0.04-0.06-0.08- " . a : . 0.15 -0.1 Figure 10: Weights for the first three components of the simulated structure images. The small dots correspond to images from the input dataset, squares and large dots represent selected values inside and outside the populated subspace for further simulations. However the benefits of such dataset are not only in the support of applicator development, but also in outcomes of further statistical analysis that could support clinical process, e.g., structure delineation or radiation planing, as well as making of clinical decisions. To conclude, it may be widely accepted that reducing dose at organs of risk is difficult without reducing dose at large tumors [16], we believe that applicator improvements based on spatial modelling could provide better alternatives. References [1] R. Potter, "Image-guided brachytherapy sets benchmarks in advanced radiotherapy." Radiother Oncol, vol. 91, no. 2, pp. 141-146, May 2009. [2] A. H. Sadozye and N. Reed, "A review of recent developments in image-guided radiation therapy in cervix cancer." Curr Oncol Rep, vol. 14, no. 6, pp. 519-526, Dec 2012. [3] M. B. Opell, J. Zeng, J. J. Bauer, R. R. Connelly, W. Zhang, I. A. Sesterhenn, S. K. Mun, J. W. Moul, and J. H. Lynch, "Investigating the distribution of prostate cancer using three-dimensional computer simulation." Prostate Cancer Prostatic Dis, vol. 5, no. 3, pp. 204-208, 2002. [4] Y. Ou, D. Shen, J. Zeng, L. Sun, J. Moul, and C. Da-vatzikos, "Sampling the spatial patterns of cancer: optimized biopsy procedures for estimating prostate cancer volume and gleason score." Med Image Anal, vol. 13, no. 4, pp. 609-620, Aug 2009. Figure 11: Central slices of simulated structure images using selected component weights from the subspace of valid structure images (left) and from other parts of the PCA space (right). [5] R. Xu and Y.-W. Chen, "Appearance models for medical volumes with few samples by generalized 3d-pca," in Neural Information Processing, ser. Lecture Notes in Computer Science, M. Ishikawa, K. Doya, H. Miyamoto, and T. Yamakawa, Eds. Springer Berlin Heidelberg, 2008, vol. 4984, pp. 821-830. [6] R. C. Conceicao, M. O'Halloran, E. Jones, and G. M., "Investigation of classifiers for early-stage breast cancer based on radar target signatures," Progress In Electromagnetics Research, vol. 105, pp. 295-311, 2010. [7] K. M. Pohl, S. K. Warfield, R. Kikinis, W. L. Grim-son, and W. M. Wells, "Coupling statistical segmentation and pca shape modeling," in Medical Image Computing and Computer-Assisted Intervention MICCAI2004, ser. Lecture Notes in Computer Science, C. Barillot, D. Haynor, and P. Hellier, Eds. Springer Berlin Heidelberg, 2004, vol. 3216, pp. 151159. [8] "GDCM library home page (version 1.x)." [Online]. Available: http://www.creatis.insa-lyon.fr/software/ public/Gdcm/ [9] D. Berger, J. Dimopoulos, R. Potter, and C. Kirisits, "Direct reconstruction of the vienna applicator on mr d a e f b d c Component 2 Component 1 Cervix Cancer Spatial Modelling for... Informatica 39 (2015)261-269 269 images." Radiother Oncol, vol. 93, no. 2, pp. 347351, Nov 2009. [10] S. Haack, S. K. Nielsen, J. C. Lindegaard, J. Gelineck, and K. Tanderup, "Applicator reconstruction in mri 3d image-based dose planning of brachytherapy for cervical cancer." Radiother Oncol, vol. 91, no. 2, pp. 187-193, May 2009. [11] M. E. Wall, A. Rechtsteiner, and L. M. Rocha, "Singular Value Decomposition and Principal Component Analysis," ArXiv Physics e-prints, Aug. 2002. [12] P. Petric, R. Hudej, P. Rogelj, J. Lindegaard, K. Tanderup, C. Kirisits, D. Berger, J. C. A. Di-mopoulos, and R. Potter, "Frequency-distribution mapping of hr ctv in cervix cancer : possibilities and limitations of existent and prototype applicators." Radiother. oncol., vol. 96, suppl. 1, p. S70, 2010. [13] K. Diaz, B. Castaneda, M. L. Montero, J. Yao, J. Joseph, D. Rubens, and K. J. Parker, "Analysis of the spatial distribution of prostate cancer obtained from histopathological images," in Proc. SPIE 8676, Medical Imaging 2013: Digital Pathology, 86760V, March 2013. [14] J. Swamidas, U. Mahantshetty, D. Deshpande, and S. Shrivastava, "Inter-fraction variation of high dose regions of oars in mr image based cervix brachyther-apy using rigid registration," Medical Physics, vol. 39, no. 6, p. 3802, JUN 2012. [15] V. Zipunnikov, B. Caffo, D. M. Yousem, C. Da-vatzikos, B. S. Schwartz, and C. Crainiceanu, "Multilevel functional principal component analysis for high-dimensional data," Journal of Computational and Graphical Statistics, vol. 20, no. 4, pp. 852-873, 2011. [16] R. Kim, A. F. Dragovic, and S. Shen, "Spatial distribution of hot spots for organs at risk with respect to the applicator: 3-d image-guided treatment planning of brachytherapy for cervical cancer," International Journal of Radiation Oncology, Biology, Physics, vol. 81, no. 2, S, p. S466, 2011.