Patient and CT imaging
A total of 19 patients with tumors of the prostate receiving scanned carbon-ion beam treatment at our hospital were enrolled. The pelvic region was immobilized with a urethane resin cushion (Moldcare®, Alcare, Tokyo, Japan) and low-temperature thermoplastic shells (Shell Fitter®, Kuraray Co., Ltd., Osaka, Japan). No water restriction was imposed, and the rectum was emptied by the patient’s effort or a laxative or enema. Two or three fiducial markers were implanted into the prostate.
Weekly serial CT scans (repeat CTs: Weeks 1–4) were acquired under breath-hold conditions at exhalation using a 320-detector CT (Aquilion One Vision®, Canon Medical Systems, Otawara, Japan) in the simulation room. CT imaging conditions were based on our clinical protocols, with a tube voltage of 120 kV and section thickness of 2.0 mm. X-ray tube current was adjusted by automatic exposure control. An initial CT (CT0) was acquired for treatment planning purposes. After the initial treatment planning scan, weekly CT scans (repeat CTs: Weeks 1–4) were acquired. Weekly CT scans were not performed on the same weekday for respective patients. The elapsed times between the planning CT and weekly CT scans are summarized in Table 3. The study involved human participants and was approved by the Institutional Review Board of our institution (N21-007) and performed in accordance with the Declaration of Helsinki. All patients provided informed consent to participate in this study.
Volume of interest (VOI) definition
Prostate, seminal vesicles, rectum, and bladder were manually delineated by certified radiation oncologists on the planning CT (CT0). The CTV included the prostate and seminal vesicles (Fig. 8A). The initial PTV (PTV1) was defined by adding a 10-mm margin to the anterior and lateral aspects and a 5-mm margin to the posterior aspect of the CTV.
Schematic drawing of the VOI. (A) The CTV (yellow line) included the prostate and seminal vesicles. Initial PTV (PTV1) was defined by adding a 10-mm margin to the anterior, lateral, and superior-inferior sides and a 5-mm margin to the posterior side of the CTV. The second PTV (PTV2) was designed as follows: A horizontal line (cut line) was set 5 mm posterior to the CTV posterior aspect on CT image sections. The part of the PTV1 region posterior to the cut line was removed (colored white in (A)). (B) The rectal ROI was expanded by 4 mm from each lateral border (colored white in (B)) and this volume was removed (PTV2′). (C) PTV2 was calculated by expanding the prostate region by 10 mm lateral expansion, 2 mm ant/post expansion, forming an approximate semicircle with max diameter at the anterior portion of the rectum. The rectum itself is cropped from the PTV2 and (D) deleted 2 mm in inferior and superior sides. This image is created by Affinity Designer version2 (Affinity, Nottingham, England, https://affinity.serif.com/en-us/designer/).
The second PTV (PTV2) was designed as follows: A horizontal line (cut line) was set 5 mm from the posterior aspect of the CTV on CT images. It is the region above the horizontal line that is colored light pink (Fig. 8A). The rectum was expanded 4 mm on laterally (both sides) and the two 4 mm-thick areas were removed (PTV2′) (Fig. 8B). Finally, PTV2 was calculated by expanding the prostate region by 10 mm lateral expansion, 2 mm ant/post expansion, forming an approximate semicircle with max diameter at the anterior portion of the rectum. The rectum itself is cropped from the PTV2 (Fig. 8C). PTV2 was calculated by reducing rectum VOI from PTV2’ and deleted 2 mm in inferior and superior side (Fig. 8D). The PTV design included many steps, however, it can be created in a systematic way and used in our clinical treatment to reduce rectum dose. A ring-shaped VOI between 2 and 6 mm away from PTV1 was created (Ring-PTV1)37.
A 2-mm reduced rectum VOI was calculated to dose distribution (Rectum-2 mm).
All VOIs on the weekly CT images were then automatically calculated from the VOIs on the planning CT using DIR employing B-spline technique with our in-house image registration software39. The DIR was performed on Windows 10 environment and is installed on a workstation (Intel® Core i7 CPU@3.7 GHz, 64 GB physical memory) and a single GPU on a board (NVIDIA GTX1080Ti ©, NVIDIA Corporation, Santa Clara, CA, USA), which is equipped with 3584 compute unified device architecture (CUDA) cores, and has 11 GB of memory. Then oncologist and medical physicist checked all VOIs on the weekly CT and modified them, if necessary. Especially, we modified prostate and seminal vesicles in several cases.
To distinguish VOIs on the plan and the treatment CTs, we use “plan-” and “treatment-” such as “plan-CTV”, “plan-PTV” and “treatment-CTV”.
Image registration
We describe the three types of image registration algorithms, which were based on the rigid image registration. And we compared to register the treatment CT to the planning CT below.
All image registrations were programmed using software developed using the C++ program language (Microsoft Visual Studio 2012®, Microsoft, Redmond WA, USA). It works under a Windows 10 environment and is installed on a workstation (Intel® Core i7 CPU3.7 GHz, 64 GB physical memory) and a single GPU on a board (NVIDIA TitanV©, NVIDIA Corporation, Santa Clara, CA, USA), which is equipped with 5120 CUDA cores, and has 12 GB of memory.
Intensity-based image registration
In general, image registration is performed by aligning one image to the other such that the differences between the two images are minimized. The most common intensity-based image registration minimizes the sum of the squared difference of plan and treatment CT images40,41,42,43. This metric requires minimal computation power and is intuitively understandable for registering two images.
The alignment value was calculated in translations and rotations by
$$\Delta V=\mathrm{arg}\underset{\Delta V}{\mathrm{min}}E\left(\Delta V\right),$$
(1)
where \(\Delta V\) is an alignment vector that denotes CT position and rotation in the treatment room.
\(E\left(\Delta V\right)\) was defined as
$$E\left(\Delta V\right)={\Sigma }_{i\in\Omega }{\left|{T}_{i}\left({V}_{p}+\Delta V\right)-{P}_{i}\left({V}_{p}\right)\right|}^{2},$$
(2)
where \({T}_{i}\left({V}_{p}\right)\) and \({P}_{i}\left({V}_{p}\right)\) were treatment and planning CT images respectively, which were set at the at the position \({V}_{p}\) determined at the planning; i is a 3D vector which denotes a point in the room; \(\Omega\) is a set of calculation points, which was assigned within the body.
Target-based image registration
In the first step, target-based image registration aligns the treatment-CTV to the plan-CTV. The translation elements of \(\Delta V\) were simply calculated by
$$\Delta V={{\varvec{g}}}_{p}-{{\varvec{g}}}_{t},$$
(3)
where \({{\varvec{g}}}_{p}\) is the center of mass of the plan-CTV.\({{\varvec{g}}}_{p}\) is calculated by
$${g}_{p}=\frac{1}{N\left({\Omega }_{p,{V}_{p}}\right)}{\Sigma }_{{\Omega }_{p,{V}_{p}}}{\varvec{i}},$$
(4)
where i is the 3D position in the room; \(V\) is a six-dimensional vector that denotes CT position and rotation in the treatment room; \({\Omega }_{p,{V}_{p}}\) is the set of 3D calculation points inside the plan-CTV when the planning CT position is \({V}_{p}\); \(N\left({\Omega }_{p,{V}_{p}}\right)\) derives the number of voxels within the \({\Omega }_{p,{V}_{p}}\). \({{\varvec{g}}}_{t}\) is the center of mass of the CTV on the treatment CT, calculated by
$${g}_{t}=\frac{1}{N\left({\Omega }_{t,{V}_{t}}\right)}{\Sigma }_{{\Omega }_{t,{V}_{t}}}i,$$
(5)
where \({\Omega }_{t,{V}_{t}}\) is the set of 3D calculation points inside the treatment-CTV when the treatment CT position is \({V}_{t}\).
Note that the rotation value of \(\Delta V\) calculated by target-based image registration is 0, which means that treatment CT moved only by translation.
WEPL-based image registration
Since WEPL is related to the charged particle beam stopping position, WEPL analysis is effective for measuring the range of variations affecting the penetration of charged particle beams. Using WEPL analysis for each beam, WEPL-based image registration registered treatment CT to the plan CT by minimizing the WEPL differences over the CTV and OAR. This metric was applied to our image registration algorithm. WEPL is calculated by the physical length (d) to the radiological water equivalent pathlength (d’) using the following equation:
$$WEPL=\sum_{i}\left(\Delta {d}_{i}\right)\cdot {\rho }_{i},$$
(6)
where i is an iterator over the points for the ray path of the treatment beam; ∆di is the pathlength in tissue with the relative stopping power ρi. The value ρi is calculated from the transformed CT image pixel values by a conversion table using polybinary calibration methods44.
To minimize the range variations, each pixel value of the CT images is changed into WEPL values and the alignment values are calculated in the same manner as the intensity-based image registration algorithm.
Depending on the patient position at the time of treatment, the center of mass of the treatment CTV may not be within the irradiation field defined by the plan-PTV before the patient setup procedure. When the treatment CTV is not within the beam field as defined by the planning PTV, the optimizer might not find the global minimum. To avoid this problem, we designed a registration error function \(E\left(\Delta V\right)\) involving a penalty function which returns the volume of the CTV’s center of mass extending from the plan-PTV (Fig. 9). The function was defined as the following;
$$E\left(\Delta V\right)={\Sigma }_{i\in\Omega }{\left|{T}_{i }\left(V+\Delta V\right)-{P}_{i}\left({V}_{p}\right)\right|}^{2}+f\left(V+\Delta V\right)+\lambda {\left|\Delta V\right|}^{2},$$
(7)
where \({T}_{i}\left(V\right) \mathrm{and} {P}_{i}\left({V}_{p}\right)\) are WEPL values transformed from treatment and planning CT images for each voxel by using Eq. (6), whose position in the simulation room is represented by \(V\) and \({V}_{p}\), which is the position determined on the planning CT; \(\Delta V\) is an alignment vector that denotes CT position and rotation in the treatment room.; i is an iterator over the set points in the room; \(\Omega\) is a set of calculation points, which is assigned:
$$\Omega ={\Omega }_{p,{V}_{p}}\cup {\Omega }_{t,{V}_{t}}\cup {\widetilde{\Omega }}_{p,{V}_{p}}\cup {\widetilde{\Omega }}_{t,{V}_{t},}$$
(8)
where \({\widetilde{\Omega }}_{p,{V}_{p}}\) is the set of calculation points in the OAR (e.g., rectum) on the planning CT when the planning CT position is \({V}_{p}\), and \({\widetilde{\Omega }}_{t,{V}_{t}}\) is the set of calculation points in the OAR on the treatment CT when the treatment CT position is \({V}_{t}\).
Schematic drawing of the penalty function: \(f\left(V+\Delta V\right)\) is the CTV on the treatment CT extending from the PTV on the planning CT. \({f}_{r}\left(V+\Delta V\right)\) is the OAR on the treatment CT overlapped with the plan PTV. This image is created by Affinity Designer version2 (Affinity, Nottingham, England, https://affinity.serif.com/en-us/designer/).
In Eq. (7), the first term on the right hand was the WEPL difference and \(f(V)\) is the penalty function, which returns the sum of the number of the treatment-CTV voxels outside the plan-PTV and the number of the voxels in the OAR on the treatment CT over the plan-PTV as below;
$$f\left(V\right)={\eta }_{1}N\left(\neg {\mathrm{\Omega {^{\prime}}}}_{p,{V}_{p}}\cap {\Omega }_{t,{V}_{t}}\right)+{\eta }_{2}N\left({\mathrm{\Omega {^{\prime}}}}_{p,{V}_{p}}\cap \widetilde{\Omega }{ }_{t,{V}_{t}}\right),$$
(9)
where \({\Omega {^{\prime}}}_{p,{V}_{p}}\) is the set of calculation points in the PTV on the planning CT when the planning CT position is \({V}_{p}\) and the hyper-parameters \({\eta }_{1}\) and \({\eta }_{2}\) were set at 1000 and 1000, respectively.
In Eq. (7) the third term is a regularizer to avoid the calculation of large alignment values \(\Delta V\) and the hyper-parameter \(\lambda\) was set at 100,000. We minimized the registration error function by using the Lucas-Kanade algorithm45 and obtained:
$$\Delta V={H}^{-1}\left({\Sigma }_{i\in\Omega }{\nabla }_{i}\left(V\right)\left({T}_{i}\left(V\right)-{P}_{i}\left({V}_{p}\right)\right)-\frac{1}{2}{\nabla }_{f}\left(V\right)\right).$$
(10)
In this expression, \({\nabla }_{i}\left(V\right)={\left(\frac{\partial {T}_{i}\left(V\right)}{\partial {r}_{x}},\frac{\partial {T}_{i}\left(V\right)}{\partial {r}_{y}},\frac{\partial {T}_{i}\left(V\right)}{\partial {r}_{z}},\frac{\partial {T}_{i}\left(V\right)}{\partial {t}_{x}},\frac{\partial {T}_{i}\left(V\right)}{\partial {t}_{y}},\frac{\partial {T}_{i}\left(V\right)}{\partial {t}_{z}}\right)}^{T},\) is the gradient of image \({T}_{i}\left(V\right)\) computed in the 6-degrees of freedom coordinate; \({\nabla }_{f}\left(V\right)\) is the gradient of function \(f\left(V\right)\) in the 6-degrees of freedom coordinate; \(H={\Sigma }_{i\in\Omega }{\nabla }_{i}\left(V\right){\nabla }_{i}^{T}\left(V\right)+\lambda I\) is the Hessian matrix, where \(I\) is the unit matrix.
To minimize Eq. (7) the Lucas-Kanade algorithm initially sets value \(V\) as \({V}_{p}\) and then iteratively solves for increments to the parameter \(\Delta V\) by calculating Eq. (11) and then the parameter \(V\) is updated:
$$V\leftarrow V+\Delta V.$$
(11)
Since we added the penalty function to minimize excessive dose for the OAR, WEL difference from WEPL-based image registration might not be always lower than that from the intensity-based image registration, target-based image registration. However, WEPL-based image registration could provide better registration position balanced respective conditions.
Treatment planning
Initial dose calculation
Initial dose distributions for carbon ion pencil-beam scanning therapy were calculated using the planning CT data, and the details of treatment planning were previously reported46. A prescribed dose of 51.6 Gy (RBE) was administered to the PTV1 with D50%. The number of treatment fractions used was 12. Two different beam angles (90° and 270°) and single field uniform dose were selected. The treatment planning parameters (beam spot position and beam spot weight) were optimized to satisfy the dose constraints using the RBE-weighted absorbed dose47 to satisfy the following metrics: PTV2-D95% and CTV-D95% > 95% of the prescribed dose, PTV1 minimum dose was > 43.8 Gy (RBE), and rectal dose was less than the reference dose volume histogram, which was prevented rectal complications derived from our prostate treatment with carbon-ion beam treatment48,49. Treatment planning optimization objectives are summarized in Table 4. All dose calculations were calculated with a pencil beam algorithm50 and performed on a workstation (RayStation 6.99; RaySearch Laboratories AB, Stockholm, Sweden).
Weekly dose calculation
Weekly dose distributions were calculated by applying the treatment planning parameters (beam angle, beam weighted spots in respective energy layers) to the treatment CTs (CT1–CT4), which were registered three times using the three different image registration techniques (see “WEPL variation with the image registration”). Since the intensity-based image registration and target-based image registration provided a single registered position for respective beam angles, doses from the two beams were summed for the respective beam angles on the same registered CT data (Fig. 10A,B).
(A) Intensity-based image registration designed to minimize the difference between the planning CT (CT0) and the treatment CT (CTn). Dose distributions for beam 1 and beam 2 were calculated on the registered treatment CT. (B) Target-based image registration of the treatment CTV isocenter to the planning CTV. Dose distributions were calculated on the registered treatment CT. (C) Treatment CTs registered with the WEPL-based image registration for Beam1 and Beam 2. Dose distribution for each beam was calculated by using the treatment CT registered by the WEPL-based method. This image is created by Affinity Designer version2 (Affinity, Nottingham, England, https://affinity.serif.com/en-us/designer/).
While the WEPL-based image registration gave different registered positions for respective beam angles, the treatment planning parameters on the workstation were defined for each beam angle, and it is impossible to change beam angles with the treatment planning parameters. Dose distributions for the respective beam angles were calculated on the treatment CT data with different registered positions (CTn for beam1 and CTn for beam2, n is integer between 1 and 4) over the four treatments (Fig. 10C). Then to calculate summed dose distribution, the dose distributions for the respective beam angles were registered to the same CT coordinates. Due to the limitations of DIR accuracy, we did not calculate time-accumulated dose distributions throughout the treatment course (CT1–CT4)51,52.
Evaluation
CTV displacements
The treatment CT was registered to the planning CT by the rigid bone matching. And then the center of mass of the CTV on the treatment CT (\({{\varvec{g}}}_{t}\)) was subtracted from that on the planning CT (\({{\varvec{g}}}_{c}\)), and it was expressed by the Euclidian distance as follows:
$$\left|\Delta {V}_{c}\right|=\left|{{\varvec{g}}}_{c}-{{\varvec{g}}}_{t}\right|.$$
(12)
Volume variations
Volume variations were calculated using prostate and bladder VOIs on the plan and the treatment CTs (CTn minus CT0).
WEPL variation with image registration
The accuracy of intensity-based image registration and target-based image registration was generally assessed visually. Since WEPL based image registration did not always produce registrations where anatomical structures matched, visual evaluation is not appropriate. Instead, WEPL variations were compared in between the plan and treatment CTs with respective CT image registration algorithms.
Computation time
Before calculating three types of CT image registrations, VOIs on the planning CT should be propagated to the weekly CTs by using DIR. Therefore, the computation time for DIR and the image registration were evaluated.
Dose assessment
Dose assessments for respective image registrations were performed using D95% to the CTV, and the percentage of the VOI irradiated with > n Gy (RBE) (Vn) for the rectum (V40, V30, and V20). These metrics were evaluated on the planning CT and the treatment CTs (CT1–CT4). All metrics in the respective plans were compared using the Wilcoxon signed-rank test (MATLAB R2021a®, Mathworks, Natick MA, USA). All p values were two-sided and those < 0.05 were regarded as statistically significant. The relationship between CTV displacement and dose metrics were evaluated using Pearson correlation coefficients.
Ethics statement
The study involved human participants and was approved by the Institutional Review Board of National Institutes for Quantum Science and Technology (N21-007) and performed in accordance with the Declaration of Helsinki. All patients provided informed consent to participate in this study.