A Laplacian Surface Deformation and Optimization Based 3D Registration Algorithm for Image Guided Prostate Radiotherapy ()
1. Introduction
This Prostate cancer is the most commonly diagnosed cancer in men in the United States and is the second leading cause of cancer related deaths in males [1]. The external beam radiotherapy (EBRT) is one of the primary treatment modalities for prostate cancer and can be used to deliver conformal and precise radiation doses to the target volume of prostate while sparing adjacent critical normal organs. Intensity modulated radiation therapy (IMRT), one of the clinical radiation delivery technologies, is being increasingly used in the definitive treatment of prostate cancer compared to the three dimensional conformal radiotherapy (3DCRT). The beam intensity is modulated and optimized with sophisticated algorithms not only to produce more conformal dose coverage to the target volume but also to further reduce radiation doses to the surrounding critical structures. However, the organ motions and patient positioning deviations during treatment impose challenges in the EBRT and may compromise the efficacy of these advanced radiation delivery technologies [2-5].
Image guided radiotherapy (IGRT) has been developed to improve the detection of target deviations relative to the target position in an approved radiotherapy treatment plan [6]. When using IGRT, verification images may be acquired at the treatment and compared to the planning images, and the deviations can be detected and corrected. The combination of advanced imaging tools with the state of the art radiation delivery technologies may potentially improve both cancer control and quality of life for cancer patients. The success of this combination hinges on the following prerequisite: precise and fast delineation and registration of the target volume between the planning and the verification images. This is because not only planned dose needs to be delivered to the correct target volume but also the process of target localization at treatment cannot be prolonged to an unacceptable level to compromise patient comfort and treatment time.
Various imaging based rigid registration methods have been developed to automate the process of determining patient positioning [7-12]. These methods have been aimed at registering 3D CT data to 2D portal images. The significant downfall of these methods is that only rigid pelvic bony landmarks are used while the location of prostate itself is not considered. In one study [13], deformable soft tissues were evaluated using an automated CT-CT image registration method; however, that work only considered translations alone. Rigid registration can be an acceptable approach to a number of medical applications. However, the detection of target deviations in radiation treatment of prostate cancers is an area where the rigid registration alone is most likely not to be successful because the dense registration and correspondence problems are involved and need to be solved.
Various full 3D deformable registration methods have been developed to register images in the prostate region, such as biomechanical models, conformal mapping, thin plate spline, diffusion based models, and free form deformation [14-26].
Biomechanical models can be time consuming since the mechanical information of the tissue is required to perform the deformation [14-16]. Conformal mapping method is applied to map bladder deformations and difficulties can be encountered to model the interfaces between structures [17]. Thin plate spline (TPS) approach formulates landmark based image registration as a nonrigid point matching system [18-20]. The main concern related to this approach is the automatic finding of landmark correspondences. In diffusion based “demons” deformation method, images are processed to be a set of isointensity contours and the source image diffuses into the target image [21,22]. The free form deformation registration method is to deform the shape of an object by manipulating a regular control lattice overlaid on its volumetric embedding space [23-26].
In one of our previous studies [26], a global-to-local shape registration model was developed to formulate shape based non-rigid registration in a hierarchical manner: first, the mutual information criterion supporting the transformation models was optimized to perform global registration; then the information of the segmented volume was implicitly used by fitting deformable superquadric mesh models and the Euclidian distances of the corresponding nodal points from the global registration process were minimized in the local registration process. However, in the first step of global deformation, maximizing mutual information can be time consuming and make the approach impractical in a clinical treatment environment. In the second step of this previous approach, the use of deformable superquadric mesh models of local deformation can be limited and become ineffective in cases where deformation becomes excessive or the target shapes are irregular beyond certain limits.
To overcome the limitations of the previous method, we developed a new global-to-local deformable registration method which incorporated both computationally efficient landmark points based global registration method and a new computationally efficient local deformable mesh model based local registration algorithm.
2. Materials and Methods
In the global registration, the corresponding landmark points were first acquired by applying the anatomy-constrained robust active shape models (ACRASM) segmentation method [26]. Then, the corresponding landmark points based global registration method was applied. The speed of this global registration method was 50 - 70 times faster than that of volume based mutual information method.
In the local registration, the Laplacian surface deformation (LSD) [27] method was employed for local deformation and Laplacian surface optimization (LSO) [28] method was employed for refining mesh model to reach an optimal registration solution.
The flow chart of the segmentation and registration process is shown in Figure 1.
2.1. Patient Data
Eight different patients undergoing external beam prostate
Figure 1. The flow chart of the segmentation and registration method of the prostate.
radiotherapy were randomly selected for the study and evaluations. Each of the patients had one set of planning CT images and one set of CBCT images. The planning CT (pCT) images were acquired with an AcQSim CT scanner (Philips Medical Systems, Andover, MA). The number of slices was different for different patients, ranging from 85 slices to 121 slices. Image slice was 512 × 512 pixels with a voxel dimension of 0.94 × 0.94 × 3 mm3. The treatment CBCT images were acquired by an On-Board Imager (OBI, Varian Medical Systems). The number of slices was also different for different patients, ranging from 55 slices to 61 slices. Each slice had 512 × 512 pixels with a voxel dimension of 0.3 × 0.3 × 2.5 mm3.
The target volumes in pCT images were contoured for treatment planning purposes on an Eclipse treatment planning system (Varian Medical Systems, CA) by an experienced radiation oncologist specializing in prostate cancer management. The target volumes were also manually delineated by the same radiation oncologist on the CBCT images. These manually delineated target volumes from the CBCT images were used as the benchmarks for comparison and evaluation. All the contours were rviewed and confirmed by another radiation oncologist.
2.2. Image Preprocessing
2.2.1. Mage Resampling and Image Contrast Enhancement
Since the pCT and CBCT image slice had different voxel dimensions and different field-of-view size, all the image datasets were resampled in both pCT and CBCT to a uniform voxel dimension 0.3 × 0.3 × 0.3 mm3. After the resampling process, the pCT and CBCT will have same voxel dimensions and similar field-of-view size.
The image intensity was first normalized to [0,1] and the histogram equalization method was applied to enhance the contrast of images in the pCT and CBCT before the image segmentation and deformable registration computation. This image contrast enhancement process helped to improve the accuracies of both the segmentation and the deformable registration.
2.2.2. Rigid Alignment
As a second preprocessing step of the segmentation and deformable registration algorithm, were all automatically aligned in a common reference frame, using a rigid transformation based on the pelvic bony structures. These bony structures normally show less variability compared to soft tissues, and can be automatically identified.
All uniformly resampled image datasets will be first cropped into the volume of interest shown in Figure 2 after rigidly aligned in a common reference frame. The anatomy-constrained robust active shape models (ACRASM) based image segmentation method will be applied on the first cropped volume [26]. The first cropped volumes will be second cropped into the approximate same dimensional volume to largely facilitate the image registration process. The second cropped volumes were generated around the region of interest (include the prostate, seminal vesicles, posterior part of the bladder and anterior part of the rectum). The location and size of the target volume (prostate) and normal organs (bladder and rectum) were used to crop the image volume. For most patients, the cropped volumes will be the same dimensional volume approximately of 150 × 150 × 100 pixels to largely facilitate the image registration process.
(a)
(b)
Figure 2. (a) The robust active shape model includes prostate (S1), posterior half internal anal sphincter (S2), posterior half external anal sphincter (S3), the outline borders between the prostate and rectum (S4), in which the anterior border is defined as the anterior part of prostate, the posterior border is defined as the posterior part of iliococoygeus muscle, and the lateral border is defined as the iliococoygeus portion of the levator ani muscle. The four red points on the model S1 is the landmark points. (b) Enhanced planning CT images unfolded from a single 3D planning CT prostate volume and the image segmentation based the anatomy-constrained robust active shape models (ACRASM).
2.3. 3D Automatic Segmentation of Treatment
Planning Images and CBCT Images
The detailed description of the method of segmentation was presented in a previous study [26]. A brief description is presented as follows.
Enhanced planning CT images unfolded from a single 3D planning CT prostate volume were segmented using the anatomy-constrained robust active shape models (ACRASM) [26]. The ACRASM shows in Figure 2(a) includes prostate (model S1), posterior half internal anal sphincter (model S2), posterior half external anal sphincter (model S3), and the outline borders between the prostate and rectum (model S4). The upper, bottom, left and right four points are landmark points of model S1 and the rest points of model S1 are regular points in each image of the 3D planning CT prostate volume. The points in model S2, S3 and S4 are not used in this paper. The landmark points of model S1 in the 3D planning prostate volume will be used in global registration as well as the local registration. Figure 2(b) illustrates the image segmentation of the prostate based the anatomy-constrained robust active shape models and shows model S1 to S4 as well as the four landmark points of model S1.
Enhanced CBCT images unfolded from a single 3D CBCT prostate volume were also segmented using ACRASM method, the same segmentation method by which the enhanced planning CT images were segmented.
The upper, bottom, left and right four points were also selected as landmark points of model S1 and the rest points of model S1 were regular points in each image of the 3D CBCT prostate volume. The landmark points of model S1 in the 3D CBCT prostate volume would be used in global registration as well as the local registration.
2.4. Nonrigid Registration
A global-to-local landmark points based deformable registration framework was utilized. We denote the 3D mesh in the planning CT images as the source shape and the 3D mesh in the CBCT treatment images as the target shape.
2.4.1. Global Registration
Global registration was developed by applying Procrustes analysis to approximately find similarity transformation, such as the translation, scaling and rotation matrix, based on the corresponding landmark points from the planning CT images and CBCT images. Let P = {p1, p2, 
, pn} and Q = {q1, q2, 
, qn} be two sets of corresponding landmark points from the planning CT images and CBCT images respectively. The goal is to find an optimal similarity transformation based on the two sets of corresponding landmark points by minimizing the following the least square function:
(1)
where, T, S, R are translation matrix, scaling matrix and rotation matrix respectively.
After solving this nine dimensional similarity search problem, the entire model can be transformed using these transformation matrices and roughly registered with the other one. After the global deformation and registration, non-rigid local deformable registration using Laplacian surface deformation optimization was employed to achieve optimal accuracy.
2.4.2. Local Deformable Mesh Model Based Registration
2.4.2.1. Local deformation using Laplacian Surface Deformation (LSD)
Laplacian surface is defined on differential coordinates. It represents each vertex point of a mesh as the difference between the point and its neighborhoods. The inputs of LSD are two sets of the control points in the planning CT and CBCT images and the initial mesh of the source shape. The mesh vertices are the vertex points which are divided into two parts. One part is control points or landmark points and the other part is regular points. In the Laplacian surface deformation process, the control points in the source images are moved to target images directly and the deformation of the regular points are calculated by LSD. Note that after global deformation, the displacements of control points are restricted in a local range. The output of LSD is the deformed mesh of the source shape.
Let the initial source mesh MS be described by a pair (VS, ES), where VS ={v1, v2, 
, vn} describes the geometric positions of the vertex points in R3 and ES describes the connectivity. The neighborhood ring of a vertex point i is the set of adjacent vertex points 
and the degree di of this vertex point is the number of elements in NS,i. Instead of using absolute coordinates VS, the mesh geometry is described as a set of differentials Δ =
, 
is the coordinate i, which will be represented by the difference between vi and the average of its neighbors:
(2)
The deformed mesh can be generated by minimizing the quadratic energy function:
(3)
where, 
' is Laplacian coordinate after deformation and vi' is Cartesian coordinate after deformation, C is the set of control points, Ti is transformation matrix for each vertex point i. The first half is to keep the similarity of the current shape and the shape generated in the previous time step and the second half try to move control points to target position. We chosed n equals 420 to present the prostate mesh model and 84 points were selected as the control points which was explained in Section 2.3.
Ti needs to be well constrained to avoid a membrane solution, which shows losing all geometric details. Because of this Ti should include rotation, isotropic scales and translations, we find anisotropic scales should not be allowed, as they will allow removing the normal component from Laplacian coordinates. Specifically, Ti is defined as:
(4)
This matrix is a good linear approximation for rotations with small angles. Furthermore, it only allows isotropic scales by enforcing the same values for diagonal elements.
The quadratic energy function can be minimized iteratively by finding s, h and t for Ti and applying the transformation on each vertex coordinates. We use the changing of residuals as the convergence criterion. When its value is smaller than a threshold (e.g., 1E - 6), we assume that the system converges and this minimization problem is solved. Note that the transformation Ti is an approximation of the isotropic scaling and rotations when the rotation angle is small. In this paper, the major rotation has been handled in the global deformation so that the local rotation fits the small angle assumption well.
2.4.2.2. Local Deformation Using Laplacian Surface Optimization (LSO)
LSO is used to improve triangle quality of a surface mesh. The inputs are landmark points and the initial surface mesh and the output is an optimized surface mesh. We will use the same notation as in Section 2.5.2.1.
Assume V is the matrix representation of VS. The transformation between vertex coordinates V and Laplacian coordinates Δ can be described in matrix algebra. Let N be the mesh adjacency matrix and D = diag(d1, ···, dn) be the degree matrix. Then, 
, where 
for the uniform weights.
Using a subset 
of m landmark points, a mesh can be reconstructed from connectivity information alone. Here the selection of A is the same as the control points selected in II.E.2.1. Positions of the reconstructed object 
can be solved separately by minimizing the quadratic energy:
(5)
where Vp' is the vertex coordinates of the reconstructed object, vap are landmark points. The first half is Laplacian constrains to smooth the object by minimizing the difference, and the second half is positional constrains to keep landmark points unchanged. In practice, with m landmark points, the 
overdetermine linear system: 
(6)
is solved in the least squares method using the method of normal equations
. The first n rows of 
are the Laplacian constrains, corresponding to the first part
, while the next m rows are the positional constrains, corresponding to the second part 
.
Iap is the index matrix of vap, which maps each v'ap to vap. The reconstructed shape is generally smooth, with the possible exception of small areas around landmark points. The minimization procedure moves each vertex point to the centroid of its ring, since the uniform Laplacian L is used, resulting in good inner fairness. The main computation cost of this algorithm is big matrix multiplication and inverse. Since A is sparse matrix, ATA is sparse symmetric definite matrix. The conjugate gradient algorithm can be employed to solve the system.
3. Validation Studies
The quantitative validation of the proposed deformable registration algorithm for prostate cancer image volumes was designed to include three components: distancebased estimators, volume-based estimators and registration efficiency. For all three evaluations, the proposed method was compared against the corresponding benchmarks as well as the registration results obtained using our previous deformable mesh model based image Registration method.
The Distance-Based Estimators and the Volume-Based Estimators
The distance-based estimators include the mean distance and the root mean square error (RMSE). The mean distance, the root-mean-square error (RMSE) and the max distance are calculated the corresponding prostate points between the treatment CBCT images and the registered planning images. The corresponding prostate points include two parts of points, corresponding landmark points and corresponding non-landmark points. The corresponding landmark points was acquired for the global registration step in the Sections 2.3. and 2.4. The corresponding non-landmark points which are defined as the closest point to the original mesh point.
The volume-based estimators include the volumetric true positive (TP), volumetric false negative (FN), volumetric the false positive (FP) and Dice similarity coefficient
between the results of proposed deformable registration method and benchmarks. We used this similarity measure since it accurately reflected the overall volumetric overlapping between binary objects and also included the information about the major causes of the differences, such as overestimate or underestimate.
4. Results
4.1. Qualitative Results and Quantitative Results of Proposed Registration Method
Figure 3 displays the qualitative results of the proposed method applied to three of the datasets. The prostate benchmarks from the planning CT and the treatment CBCT are shown in Figures 3 (a) and (b), respectively. The prostates in column (a) were correspondingly registered to column (b) using the proposed method. Figure 3 (c) shows the registered prostates (in blue) superimposed on the corresponding prostate benchmarks from the treatment CBCT images. Figure 4 shows the axial view of the enlarged image results of Figure 3. Figure 5 shows the spatial error (red lines) between themesh point of the registered model (blue) from pCT images and the closest point of the segmented prostate model from CBCT images.
The quantitative results are summarized in Table 1 and Table 2. Table 1 shows the values of the mean distance (ranged from 0.43 mm to 2.23 mm), the rootmean-square error (RMSE) (ranged from 0.51 mm to 2.46 mm) and the max distance (ranged from 2.48 mm to 5.32 mm) of the corresponding prostate points between the treatment CBCT images and the registered planning images using the proposed deformable registration method. Table 1 also shows the values of the mean distance (ranged from 0.38 mm to 2.20 mm), the root-mean-square error (RMSE) (ranged from 0.45 mm to 2.36 mm) and the max distance (ranged from 1.59 mm to 5.98 mm) of the corresponding prostate points between the treatment CBCT images and the registered planning images using the previous deformable mesh model based image registration method [26]. These results demonstrate that the developed method produced reasonably good agreement between the registered prostates and the corresponding benchmarks using the distance-based estimators.
(a) (b) (c)
Figure 3. Deformable registration results of Laplacian surface deformation and optimization: (a) segmented prostate model (green) from planning CT images; (b) segmented prostate model (red) from CBCT images; and (c) registered model (blue) superimposed on (b).
Figure 4. The enlarged image results of the axial view of of Figure 3. The registered model (blue) from pCT images superimposed on the segmented prostate model (red) from CBCT images using Laplacian surface deformation and optimization deformable registration.
Figure 5. The spatial error (red lines) between the mesh point of the registered model (blue) in the pCT images and the corresponding closest mesh point of the segmented prostate model of CBCT images.
Table 2 shows the values of the volumetric false positive volumetric true positive (TP ranged from 84.5% to 93.7%) and Dice similarity coefficient (ranged from 85.5% to 93.2%) between the results of proposed deformable registration method and benchmarks. Table 2 also shows the values of the volumetric false positive (volumetric FP ranged from 7.53% to 14.56%), volumetric false negative (FN ranged from 5.03% to 14.8%), volumetric true positive (TP ranged from 85.2% to 94.97%) and Dice similarity coefficient (ranged from 86.0% to 93.8%) between the results of the previous deformable mesh model based registration method and benchmarks. These results further confirm that the developed method produced reasonably good agreement between the registered prostates and the corresponding benchmarks.