Simulation of pseudo-CT images based on deformable image registration of ultrasound images

Purpose: Imaging of patient anatomy during treatment is a necessity for position verification and for adaptive radiotherapy based on daily dose recalculation. Ultrasound (US) image guided radiotherapy systems are currently available to collect US images at the simulation stage (US sim ), coregistered with the simulation computed tomography (CT)


INTRODUCTION
Image guidance has become an essential part of radiotherapy (RT) treatment to allow for safe delivery of radiation doses.
Image guided RT (IGRT) is often performed for several or all treatment fractions to position the patient correctly. Beyond the aim of image guidance, the availability of daily imaging also allows for the possibility of adaptive RT (ART). 1,2 The goal of ART is to improve RT treatment by systematically monitoring dose discrepancies and incorporating them to reoptimize the treatment plan. Normally only the planning computed tomography (CT) image, acquired at simulation stage, is available for the dose calculation, but both interfraction and intrafraction patient anatomy motion and changes (like tumor shrinkage, nodal volume changes, and weight loss) may alter the dose distribution. [3][4][5][6] In ART, the anatomy from the planning CT is updated by the anatomy from the daily imaging, acquired during the IGRT workflow to monitor dose distribution and if necessary adapt the treatment plan.
CT scanners are usually not available in the treatment room. Instead, cone-beam computed tomography (CBCT) can be used for dose calculations either directly [7][8][9][10] or indirectly with deformable image registration (DIR) 11,12 even though they offer a lower image quality when compared to CT scanners. In some studies, using the CBCT directly for dose calculations, the inaccuracies in the Hounsfield units (HUs) are large enough to result in clinically relevant dose errors. [13][14][15] In this paper, a workflow is introduced to produce pseudo-CT images based on deformable registration of ultrasound (US) volumes. A 3D US IGRT system can acquire volumetric, high-contrast soft-tissue images noninvasively on a daily basis without using ionizing radiation (Fig. 1). Subsequently, deformable registration of these volumes can reveal changes in tissue distribution that occurred over time.
Relatively few papers on US to US deformable registration can be found in the literature and as far as we could find, there are presently no papers involving deformable registration of pelvic or abdominal US volumes in RT. In other medical fields, however, some publications are available. For example, Shekhar et al. 16 proposed a nonrigid method based on mutual information to register cardiac US images in different phases throughout the complete cardiac cycle.
A similar workflow as proposed in this study was presented for brain surgery applications by Pennec et al. 17 In this study, preoperative magnetic resonance (MR) images and US images were acquired. Subsequently, intraoperative US images were used to create pseudo-MR images of the brain. This resulted in acceptable representations of the brain anatomy during surgery.
As these results were promising, we used a similar approach to create pseudo-CT (CT ps ) images. We hypothesize that a pseudo-CT image can be created based on CT sim using a deformation field calculated between US sim and US tx . We expect that the CT ps so created gives a better representation of the patient's anatomy during treatment delivery than the planning CT sim .

2.A. The concept
In the proposed workflow (Fig. 2) for CT ps image creation, DIR has to be performed to calculate a deformation field between US sim and US tx . Subsequently, this deformation field has to be applied to CT sim which results in the creation of CT ps .

2.B. Patient scans
Clinical examples with multiple coregistered US-CT combinations at the simulation stage (instead of the treatment stage) were used to validate the concept. In this study, three prostate cancer patients from a previous study 18 were used. Due to clinical reasons, these patients underwent additional US and CT imaging next to US sim and CT sim acquisitions. In the normal clinical workflow, these extra CT and US images are not acquired. The extra CT scans were used as ground truth (CT gt ) scans to which the derived CT ps scans can be compared in this proof of concept study. In Table I, the method used to calculate and evaluate the result from the deformations is described.
The coregistered CT-US images were acquired at two time points for patients 1 and 2 (three and one weeks apart, respectively). Acquisitions for patient 3 were made for five time points where the first two were two weeks apart and the following three time points were one week apart.
F. 1. Workflow of acquisition of CT sim , US sim , and US tx images (Clarity US system; Elekta) (adapted from Elekta with their permission).
F. 2. DIR is computed between the two US images (US sim and US tx ) and then applied to CT sim ; a new pseudo-CT ps is obtained. The question is whether this CT ps is indeed representative for the patient anatomy during treatment.
All coregistered US-CT combinations were acquired in the CT-room with the patient's external skin markers positioned along the room lasers. The 3D US scans (Clarity system; Elekta, Stockholm, Sweden, voxels: 1 × 1 mm 2 × 3 mm slice thickness; US probe type C5-2/60, center frequency 3.5 MHz, Sonix Series; Ultrasonix Medical Corporation, Richmond, BC, Canada) were performed transabdominally immediately before or after the CT scan. The number of voxels of the US images varied between [512, 512, 90] and [512, 512, 131]. For each patient, the images were resampled to match the dimensions of the first acquired US volume (US sim ).
The CT scans were acquired using a SOMATOM Sensation Open (Syngo CT 2006A, Siemens, Germany; voxels: 1 × 1 mm 2 × 3 mm slice thickness). Both scans were performed in the same supine patient position, stabilized with knee fix and foot support (Combifix, Civco Medical Solutions, Kalona, IA, USA), resulting in a correct automatic fusion of the US and CT images. 19 In all US images, the prostate was delineated. All CT images had delineations of the body contour, prostate, seminal vesiculae (SV, except for patient 3), anus, rectum, and bladder (except for patient 1).

2.C. Deformation
For each US-CT combination (as detailed in Table I method from ElastiX; Utrecht, The Netherlands). 20,21 Prior to the deformation field calculation, all volumes were resampled to the same image dimensions per patient. In addition, segmentation of the CT sim images resulted in a binary mask of the bones and the region of interest (ROI) was defined as the overlapping parts of the US images (ROI: US sim ∩ US tx ). All these preprocessing steps were performed in the  (MathWorks, Inc., Natick, MA) software. During the acquisition of the different US-CT combinations, the patients were in the same position with the body markers aligned to the lasers. For this reason, no rigid transformation was performed prior to the deformable registration, in particular, to prevent erroneous full body shifts based on internal shifts of the prostate. 22 As mentioned before, the deformable registration was performed using the ElastiX software. This software package requires three inputs: fixed image (US tx ), moving image (US sim ), and a parameter file. The parameter file contains all the parameters that determine the characteristics of the registration. In Sec. A of the supplementary material, 23 an example of such a parameter file is detailed.
In this study, the deformable registration was performed either on the overlapping parts of the US images or on binary masks of the delineated prostate volumes only. In total, five different parameter sets (parameters A-E in Table II) were defined for this purpose using the file in Sec. A of the supplementary material 23 as a basis.
The deformation field calculations were based on the overlapping parts of the US images, but were propagated further through the image (Fig. 3). Also bones were sometimes present in these overlapping parts. As bones are in principle rigid structures, they are not expected to undergo deformations. Therefore, the binary bone mask defined during preprocessing was input in the rigidity penalty 23 of ElastiX to prevent bones from deforming.

2.D. Evaluation of the deformation
The created CT ps and the deformed CT delineations were then compared to the ground truth, i.e., the corresponding CT gt and its delineations. The contours were evaluated using the dice similarity coefficient [DSC = (2|X ∩Y |)/(|X | + |Y |)]. A DSC ratio of 1 indicates complete overlap, while 0 indicates no overlap. T II. Five different parameter sets (A-E) were used during the deformable registration. This registration could be based on the whole US volume or on the binary mask of the delineated prostate volume only (reported in the columns: fixed image and moving image). In addition, both the metric and iterations were varied among the different sets. In addition, the prostate contours were also evaluated using both the maximum absolute distance (MAX) and the mean absolute distance (MAD). 24 The MAX defines the largest difference between two contours, e.g., prostate contour A and prostate contour B. For each point a on prostate contour A, the minimal distance to all points on prostate contour B was calculated. The same was repeated for each point b on prostate contour B with respect to prostate contour A. This resulted in a set of minimal distances and the maximum of this set is referred to as MAX. Calculating the mean of this set gave the MAD.
The CT sim and CT ps images were compared to CT gt using a gamma (γ) index evaluation. 25,26 The γ index is commonly used for dose evaluations. Prior to the index calculation, two acceptance criteria need to be set: voxel-by-voxel numerical dose difference and distance-to-agreement (DTA: distance between a voxel on one volume and the nearest voxel in the other volume that has the same dose). The resulting index gives information on a voxel scale, while taking the voxels in the vicinity into account as well.
In this case, not only dose was evaluated with the γ index but also HU (γ CT ). The γ values were calculated using an inhouse developed method 27,28 using  and ++. The used method allows the sign of the γ value to indicate whether an overdose (γ > 0) or underdose (γ < 0) is found for each voxel. 28 In this case, because we evaluate HU, a γ > 0 means that the HU is relatively higher than the reference and γ < 0 means that the HU is relatively lower. A value |γ| > 1 in a voxel indicates that the voxel fails to meet the acceptance criteria; in this case, a 50 HU voxel intensity difference and a 3 mm distance-to-agreement. (The 50 HU is a conservative measure based on that for typical radiotherapy beams; to produce a 1% error in dosimetry would require errors of over 8% in bone electron density 29 and hence HU. The 3 mm distanceto-agreement is a commonly used criterion in dosimetry. 26 ) The percentages of the volume with a |γ CT | > 1 within the contours "intersection body contours," "prostate," "anus and rectum," and "bladder" were reported. The percentages of gamma failure and DSC evaluations are reported using the contours of the CT gt , except for the intersection body contours which is the overlapping part of the body contours of both CT sim and CT gt .

2.E. Dose calculation and evaluation
Dose distributions were obtained by recalculating the original treatment plans (five-beam IMRT plans; XiO CMS 4.51, Elekta, Stockholm, Sweden) designed on the planning CT sim , on the CT sim , CT ps , and CT gt scans. For this, an inhouse developed software was used, based on Monte Carlo F. 3. Example of overlap between CT (gray) and US (color) (a) and between two US images (b) of patient 1. US-based DIR can only be performed on the area where both CT and US information (of both US sim and US tx ) is available. In this example, only the prostate and its surrounding tissue, e.g., a part of the bladder, are present in both US images. In (c), only the overlapping area of both US images (yellow contour) contains information where the deformation field (2D representation with red arrows) is based on. The field propagates further beyond this border (see color version online). simulation using the XVMC code. 30,31 Dose distributions on the CT sim and CT ps images were compared to the dose on CT gt using a γ evaluation (γ Dose ), 25 with acceptance criteria of 3% dose difference and 3 mm distance-to-agreement. Again the percentage of the volume with a |γ Dose | > 1 within the contours intersection body contours, prostate, anus and rectum, and bladder was reported.

RESULTS
In most cases, deformation did improve the results according to all evaluation methods, although these improvements were in some cases very small or even negligible. Only for patient 1, there was a large improvement (more than 10% decrease in the volume with |γ Dose | > 1) in the dose of the prostate when the intensity based normalized-correlation metric with 100 iterations (parameter set C) was used (Table III).
In Fig. 4, an example is given for patient 1 using parameter set C. In the second column, the overlap of the prostate and anus and rectum contours is shown. DSC increased by 0.3 when the deformations were used. The third and fourth columns show the γ CT and γ Dose values. In the overlapping body contours, the percentage of γ CT failure decreased by 1.7% in volume. For the prostate and anus and rectum contours, there was a γ CT failure decrease of 9% and 8.4%, respectively. For the dose, the volume percentage of γ dose failure decreased by 11.2% in volume for the prostate. Yet the percentage of γ dose failure decreased by only 0.6% and 0.0% for the overlapping body contours and anus and rectum contours, respectively.
All available results for patient 1 are summarized in Fig. 5. Figure 5(A) shows that the DSC improved for all parameter sets. For prostate, the best results were obtained with parameter set E; for anus and rectum, set C performed best. Both the MAD and the MAX where smaller compared to the reference situation [ Fig. 5

(B)]. Figures 5(C) and 5(D)
detail results on gamma failure, respectively, based on CT values and dose. In case of CT based evaluations, the best results were achieved using parameter set B for prostate and anus and rectum and using parameter set D for the body contours. For the dose based evaluations, parameter set C gave the best results in all cases. The analyses were repeated for all available patient data and the overview of the results is detailed in Fig. B of the supplementary material. 23 Evaluation of all patient cases [  23 ] shows that the DSC of the prostate increased the most for the two contour based parameter sets (D and E). (Parameter set E with 300 iterations did not succeed in the deformation of patient 2 because there was a too small overlapping volume. Therefore not enough voxels could be mapped and the registration failed to find a solution.) Only for patient 3a, none of the parameter sets gave an improvement for any of the contours. Overall, the maximum changes in DSC for the intensity based normalizedcorrelation parameter sets were a decrease of −0.5 or an improvement of +0.3. For the contour based parameter sets, these were −0.3 and +0.4. T III. Five evaluation methods were used to evaluate the delineated prostate contours. The first and second columns detail the patient and the used evaluation method. Both gamma index values show the volume percentage of gamma failure, [γ CT(50 HU,3 mm) > 1] and [γ Dose(3%,3 mm) > 1], respectively. In the third column, the reference situation (comparison between CT sim and CT gt ) can be found. In the final five columns, the results for each of the parameter sets (A-E) are detailed. The bold numbers indicate which parameter sets resulted in the same result or in an improvement with respect to the reference. 10.3 9.9 9.5 8.7 6.9 6.5 Note: DSC, dice similarity coefficient; MAD, mean absolute distance; MAX, maximum absolute distance; and γ CT , γ CT(50 HU,3 mm)>1 ; γ Dose , γ Dose(3%,3 mm)>1 .
For the changes in CT HU values, the percentage of the volume with a |γ CT(50 HU,3 mm) | > 1 for prostate is shown in Table III and for the other contours, in the supplementary material [ Table B and Figs. B(C,G,K,O,S,W) 23 ]. A maximum improvement was seen of 20.5% (14.6% for contour based) and the poorest results gave an increase of 3.2% (2.2% for contour based) in the volume with |γ CT(50 HU,3 mm) | > 1.
Looking at the prostate results as shown in Table III, in case an improvement was achieved, the contour parameter set (D, 100 iterations) seemed to give an improvement in most cases, yet it was not always the best one. The results for the other contours (body, anus and rectum, and bladder) that can be found in Table B in the supplementary material 23 confirm this as well.
F. 4. Results for patient 1 (parameter set C). In the first column, the CT sim and CT ps (in pink) are compared to CT gt (in green). In the second column, the contours of prostate (P) and anus/rectum ( A/R) are compared. When the images are grayscale (column 1) or white (column 2), there is overlap between the compared images. The third and fourth columns show the γ CT (column 3) and γ Dose (column 4). In green, the γ values are between −1 and 1. In red and blue are the voxels in which the γ failed to meet the criteria of (50 HU, 3 mm) for the CT values and (3%, 3 mm) for the dose. For column 4, the areas where there is an underdosages compared to CT gt (γ Dose < −1) are shown in blue. In red, there is an overdosage compared to CT gt (γ Dose > 1) (see color version online).

DISCUSSION
We have evaluated the impact of applying US-derived tissue deformations to approximate CT images to the real anatomical organ position of prostate patients during radiation therapy. As noted before, a similar workflow was presented by Pennec et al. 17 for brain surgery applications. However, in that study, pseudo-MR images of the brain were created. To our knowledge, this is the first time a similar method is used for RT applications.
In this study, patients 1 and 3d would have benefited most from the deformations (>3% volume decrease for the volume with a |γ Dose | > 1). In addition, the difference in dose between CT sim and CT gt was there also the largest (>10% volume with a |γ Dose | > 1). For the other patient cases, the improvements were not clinically relevant.
Ideally, one should be able to evaluate beforehand which patients would benefit from applying the deformations. The only metric that is available prior to DIR and could be suitable is the DSC of the prostate contours on US sim and US tx . A statistical evaluation was performed to find a possible correlation between these DSCs and the effect on the dose deposition on the prostate (|γ Dose | > 1). Unfortunately such a correlation was not found, possibly due to the limited number of patients. However, there seems to be a trend that the patients with the largest geometric changes benefit most from deformations, but a future study with a larger image database will be necessary to validate the predictive power of this DSC parameter to get a clearer indication when it is worthwhile to perform DIR.
Besides a larger database to perform statistics, such a database could be used to find an optimal metric and parameter set for the DIR. For this proof-of-principle study, two deformation metrics were used and only the number of iterations varied. Optimization of the metrics and parameter set may improve the results. In the current study, the results of the evaluation methods were not always in agreement. Even between the CT and dose values, there were some differences due to the cumulative effect of the dose along the beam path. The differences between change in γ CT and γ Dose are caused by the fact that the dose in the organs is not only dependent on the local HU but also on the HU along the beam path. The best evaluation method is dependent on the purpose; the evaluation of the best parameter set should therefore always be assessed with the correct evaluation method. In case of ART, this could be γ Dose(3%,3 mm) .
A limitation of an US-based deformation field is that the volume of the CT on which one can directly calculate the deformation field is limited to the volume of the US data available (Fig. 3). The deformation field propagates further, but this is not based on image data and is therefore maybe less reliable. For patient 2, a small overlap of US volumes resulted in a failure in parameter set E. Standardization of scanning, so that at least the complete prostate is visible and the US volume overlap is maximal, and US images with larger fields of view may improve the results. Transperineal scanning with a larger image sector or perhaps even fusion of multiple US scans from different directions can extend the field of view.
However the US image will never completely overlap the CT image, therefore part of the deformation field will still be based on only an extrapolated deformation field. For an ideal exact extrapolation, it may be crucial to take into account the mechanical properties of tissues and organs, such as skin, bones, and bladder, which are positioned outside of the overlapping US images. In this work, some deformation field propagation outside of the overlapping US volumes is already inherently taken into account, due to the use of the so-called multiresolution approach during the deformable registration. In this approach, the registration starts with images that have a lower complexity. For example, images that were smoothed and possibly down sampled. During the registration, a B-spline control point grid is overlaid on the fixed image. This grid is always rectangular. Control points that are outside of the region of interest (overlapping parts of the US volumes) are in principle not affected. However, due to the multiresolution approach, the control point spacing is larger at lower resolutions than at higher resolutions. For this reason, a larger area around the region of interest is affected at lower resolutions, which typically produces deformations outside of the region of interest.
Another reason why it is important to have standardization of the US scanning is that, just like with the IGRT usage of the US images, it is important to have reproducible US images. In particular, the probe pressure 18,32 and speedof-sound aberration 33 along the imaging beam should be comparable. One cannot distinguish between the US imaging dependent changes caused by nonstandardized procedures or a real anatomy changes. Therefore it is best to prevent them or correct [34][35][36][37][38] for them before the DIR procedure. For our specific cases, preliminary inspection revealed that these corrections were not necessary.
Validation of the DIR methods in general is also still necessary to reliably perform DIR for ART. Different deformation algorithms lead to different results, therefore more research is necessary.

CONCLUSIONS
It was possible to generate a pseudo-CT ps with the use of DIR based on US imaging which was more representative of CT gt than CT sim . For the patients with the smaller prostate change over time, the procedure did not improve the dose calculations much. The largest improvements were seen for patients with the largest anatomical changes. More research with a larger image database is necessary to find an optimal deformation metric and parameter set. With a larger database, it might be possible to find a predictive measure and criteria to decide whether DIR is worthwhile for individual patients.