WO2015081079A1 - Software for using magnetic resonance images to generate a synthetic computed tomography image - Google Patents

Abstract

A computing device may assign voxels of at least one magnetic resonance image to tissue region types from a set of region types; and optimize voxel intensity value weights for each of the region types to minimize a difference between (i) actual voxel intensity values of the planning computed tomography image and (ii) synthetic computed tomography voxel values determined by weighting the voxels of the region types of the at least one magnetic resonance image according to the voxel intensity value weights of the corresponding region types. The device may also assign voxels of magnetic resonance images of a treatment area of a patient to tissue region types from a set of region types; and generate a synthetic computed tomography image of the treatment area of the patient according to the magnetic resonance images and optimized voxel intensity value weights for each of the tissue region types.

Description

SOFTWARE FOR USING MAGNETIC RESONANCE IMAGES TO GENERATE A

SYNTHETIC COMPUTED TOMOGRAPHY IMAGE

CROSS-REFERENCE TO RELATED APPLICATIONS

[0001] This application claims the benefit of U.S. provisional application Serial No.

61/909,138 filed November 26, 2013, the disclosure of which is hereby incorporated in its entirety by reference herein.

TECHNICAL FIELD

[0002] The present disclosure relates generally to the field of medical imaging systems and, more specifically, to generation of synthetic computed tomography (CT) images using magnetic resonance (MR) image data.

BACKGROUND

[0003] Tomography refers to a technique for capturing a two-dimensional slice or cross- sectional image of an object, through the use of radiation or any kind of penetrating wave. The word tomography is likely derived from the Greek words, tomos (slice, section) and grapho or graphein (to write). Computed Tomography (CT) refers to a medical imaging technique that uses X-rays and computer processors to collect and display a series of two-dimensional images or tomograms of an object. Additional computer processing and mathematics can be used to generate a three- dimensional image. Magnetic Resonance (MR) imaging refers to a tomographic technique that uses a powerful magnetic field and a radio frequency transmitter to capture detailed images of organs and structures inside the body.

BRIEF DESCRIPTION OF THE DRAWINGS

[0004] FIGURE 1 illustrates an exemplary image construction workflow;

[0005] FIGURE 2A illustrates an example procedure for weight optimization and optimizing a number of patients used during weight optimization; [0006] FIGURE 2B illustrates an example relationship between full field-of-view MAE and number of patients used during optimization;

[0007] FIGURE 2C illustrates an example relationship between dosages for planning CTs and synthetic CTs as a function of number of patients used during optimization;

[0008] FIGURE 3A, 3B and 3C illustrates an exemplary synthetic CT image at various stages during the image construction workflow;

[0009] FIGURES 4A and 4B illustrate an exemplary comparison of dose distributions for a planning-CT based treatment plan (FIGURE 4A) and a synthetic CT-based plan of an exemplary patient (FIGURE 4B);

[0010] FIGURE 5 illustrates an exemplary representative dose-volume histogram including prostate and bony structures;

[0011] FIGURE 6 illustrates an exemplary system for creation of synthetic CT images;

[0012] FIGURE 7 illustrates an exemplary validation workflow performed using the exemplary system; and

[0013] FIGURE 8 illustrates an exemplary clinical workflow performed using the exemplary system.

SUMMARY

[0014] In at least one embodiment, a computing device may acquire, of a treatment area of a patient, a planning computed tomography image and at least one magnetic resonance image; assign voxels of the at least one magnetic resonance image to tissue region types from a set of region types; and optimize voxel intensity value weights for each of the region types to minimize a difference between (i) actual voxel intensity values of the planning computed tomography image and (ii) synthetic computed tomography voxel values determined by weighting the voxels of the region types of the at least one magnetic resonance image according to the voxel intensity value weights of the corresponding region types. [0015] In another embodiment, a computing device may assign voxels of magnetic resonance images of a treatment area of a patient to tissue region types from a set of region types; and generate a synthetic computed tomography image of the treatment area of the patient according to the magnetic resonance images and optimized voxel intensity value weights for each of the tissue region types, the optimized voxel intensity value weights computed according to a minimization of a difference between (i) actual voxel intensity values of training computed tomography images and (ii) magnetic resonance voxel intensity values of training magnetic resonance images weighted according to the optimized voxel intensity value weights.

[0016] In another embodiment, a method includes assigning voxels of magnetic resonance images of a treatment area of a patient to tissue region types from a set of region types; and generating a synthetic computed tomography image of the treatment area of the patient according to the magnetic resonance images and optimized voxel intensity value weights for each of the tissue region types, the optimized voxel intensity value weights computed according to a minimization of a difference between (i) actual voxel intensity values of training computed tomography images and (ii) magnetic resonance voxel intensity values of training magnetic resonance images weighted according to the optimized voxel intensity value weights.

DETAILED DESCRIPTION

[0017] As required, detailed embodiments of the present invention are disclosed herein; however, it is to be understood that the disclosed embodiments are merely exemplary of the invention that may be embodied in various and alternative forms. The Figures are not necessarily to scale; some features may be exaggerated or minimized to show details of particular components. Therefore, specific structural and functional details disclosed herein are not to be interpreted as limiting, but merely as a representative basis for teaching one skilled in the art to variously employ the present invention.

[0018] A system and computer-implemented method is described utilizing an approach for generating synthetic CT images from MR data in order to create accurate electron density maps for dose calculation in treatment planning systems. By way of explanation and not limitation, the disclosure discusses nine selected retrospective patients that had been treated using a CT-only workflow while undergoing additional MR imaging. As discussed in detail below, synthetic CT images generated by way of the approach may be used for creating treatment plans comparable to those based on planning CT images.

[0019] The disclosure describes a system and computer-implemented method for generating a synthetic CT image from a set of common clinical MR acquisitions that may be used in treatment planning systems for accurate dose calculations. The disclosure discusses data captured from nine retrospective patients with prostate cancer, where the exemplary patients underwent the full CT- based radiotherapy workflow as well as additional MR imaging prior to treatment.

[0020] The system may use the imaging from a subset of the exemplary patients as a training set to optimize a set of weights, and may generate synthetic CT images for patients according to the optimized weights. The system may further analyze CT value differences between the original and synthetic CT images, and the synthetic CT images were then imported into a treatment planning system. The system may generate or receive a recontouring of key structures on the MR or synthetic CT images, and may copy the original planning CT-based treatment plans that had been used for the patients during radiotherapy onto the synthetic CT images. The system may recalculate the dose based on the synthetic CT images, and compare the dosimetric results for the prostate in the original and recalculated plans.

[0021] Using the exemplary set of patients, the system generated Synthetic CT images whose average whole -body mean absolute error (MAE) was approximately 73 Houns field Units (HU). The average difference in HU between the planning CT and synthetic CT was less than 1% for the prostate, bladder, and rectum and less than 5% for the left and right femora. Moreover, the averages for the D99 as well as the mean and max dose to the prostate of the synthetic CT-based dose calculations were all within 1% of the planning CT-based dose calculations.

[0022] Accordingly, synthetic CT images created using the described systems and computer- implemented method may be used in treatment planning systems for accurate dose calculations. Thus, the synthetic CT images may be used to support the effort to move to an MR-only workflow for radiotherapy planning. [0023] Computed tomography (CT) serves as a fundamental imaging modality used in radiation treatment planning simulation for 3D conformal radiation therapy, but there has been recent interest in removing CT from radiotherapy planning altogether and instead using a magnetic resonance (MR)-only planning workflow^1"5. One advantage of MR imaging is its excellent soft tissue contrast, which contrasts with the poor soft tissue contrast provided by CT images that often renders it difficult, if not impossible, to accurately delineate certain key structures. Another advantage of MR imaging is that MR avoids the radiation dose to the patient inherent in receiving a CT scan. While MR imaging can be used in combination with CT for contouring of soft tissue structures^6"9, the process introduces systematic errors during MR-to-CT registration. These errors may be especially prominent in prostate and gynecological patients¹⁰' ^u. Additionally, the added expense of having multiple simulation modalities is often prohibitive.

[0024] For an MR-only workflow to be feasible, there are certain issues that should be

12

addressed . One of these issues may be the geometric uncertainty of MR imaging caused by gradient nonlinearity, inhomogeneities in the main magnetic field, and magnetic susceptibilities in the patient. A method to correct for gradient nonlinearities is to use a spherical harmonic expansion to model the magnetic fields generated by the reference coils¹³'¹⁴ which has been shown to be able to reduce the errors caused by gradient nonlinearities to subvoxel distances. Manufacturer MR scanner development may also be able to reduce the errors caused by main field inhomogeneities in modern scanners to one millimeter or less after calibration¹⁵'¹⁶. While distortion caused by patient-specific magnetic susceptibility can be significant, modern MR sequences such as 3D turbo spin echo sequences may provide acceptable tolerance for radiotherapy treatment purposes¹¹.

[0025] Another issue with MR-only workflow, and an aspect addressed by the systems and methods of the disclosure, is that modern treatment planning systems rely on electron density values

17 for accurate dose calculations, especially in the presence of large contrast inhomogeneities in the patient. Whereas CT intensity values in Hounsfield units (HU) may be converted to electron density values using a look-up table, there is no equivalent relationship between MR intensity and electron density values. Therefore, a method must be developed for generating an electron density map from a set of MR images. [0026] Several methods have been introduced for generating a synthetic CT image that can be used as a basis for treatment planning from a set of MR acquisitions. One approach has been to either manually¹¹'¹²'¹⁸ or automatically¹⁹ segment an MR image into a number of regions and to then assign bulk CT values to all those voxels lying within each of the regions. Another approach is based on gathering a large patient cohort that had undergone both CT and MR imaging and using them to

20 generate reference CT and MR atlases that will be deformably registered to each other. The reference MR atlas would next be deformably registered to a new MR image, and those deformations would then be applied to the reference CT atlas to derive the "pseudo CT" image for the new patient. A disadvantage for this type of method is that it suffers from added uncertainty when the patient has non-standard geometry. A recently introduced approach uses a Gaussian

21 22

regression model to derive CT numbers directly from MR intensity values ' , but this method depends on the use of specialized ultra short echo time (UTE) sequences (TE < 0.1s).

[0027] The disclosed systems and methods utilize common clinical MR acquisitions to generate a synthetic CT image using a voxel-based approach. This approach neither requires the use of specialized sequences, such as UTE, that are not widely available, nor necessitates using a large number of sequences for which the patient would have to remain on the table for a prolonged period of time. Using the disclosed approach, synthetic CT images may be generated for a set of prostate cancer patients that have undergone MR imaging. The disclosed approach may further perform a dosimetric evaluation after applying the original treatment plans to the new synthetic CT images. While the disclosure presents prostate image processing in several examples, it should be noted that the disclosed approach is also applicable to other targets and critical structures needing enhanced tissue differentiation and boundary delineation, of which lungs and brain may be additional examples.

[0028] Regarding methods and materials, and by way of explanation and not limitation, the disclosed approach is discussed with respect to certain exemplary acquired data. For example, the analysis may be performed retrospectively on data from an exemplary set of nine prostate cancer patients enrolled in an IRB-approved study to evaluate the performance of MR simulation. In addition to CT simulation, the analysis may be performed on MR data from the prostate patients including MR imaging of the whole pelvis using Tl -weighted, T2 -weighted, and balanced turbo field echo (bTFE) sequences. In another example, synthetic CTs may be generated for brain patients, although since bTFE sequences may be unavailable for brain scans, a fluid attenuated inversion recovery scan used for brain patients may be substituted for the bTFE sequences.

[0029] Acquisition of CT images may be performed using a Philips Brilliance Big Bore.

Exemplary settings for CT acquisition may include a tube voltage of 140 kVp and a tube current- time product of 500 mAs. As one possibility, a slice thickness of 3 mm may be used, where each slice has dimensions of 512 x 512 pixels with isotropic spatial resolution of 1.2754 mm. To ensure reproducibility in patient positioning for CT acquisition and treatment, a flat table-top insert (Civco, Orange City, IA) may be utilized. For instance, patients undergoing scans may be immobilized using leg bands, a ring to hold on the chest, and a shaped foam pad for leg immobilization.

[0030] Acquisition of MR images may be performed on a one Tesla (IT) open-platform MRI scanner (Panorama, Philips Healthcare). As an exemplary possibility, three sequences may be used during scanning for each of the patients: a Tl -weighted fast field echo (FFE) sequence, a high resolution T2-weighted turbo spin echo (TSE) sequence, and a bTFE sequence using spectral presaturation with inversion recovery (SPIR). These sequences may be used as standard of care in clinics (e.g., for prostate patients) and are compatible with the disclosed approach. In addition, an inverse intensity volumetric image may be created from the Tl sequence for use in the approach, (e.g., by subtracting weighted intensity values from an intensity value below which lay 95% of the total area under the curve on the histogram of the Tl -weighted image). An exemplary slice thickness for the MR sequences may be 2.5 mm. As one possibility, both the Tl and T2 images may be collected with slice dimensions of 640 x 640 pixels and isotropic spatial resolution of 0.6568 mm, while the bTFE slice dimensions may be 432 x 432 pixels with isotropic spatial resolution of 0.9491 mm. To ensure patient positioning remains consistent with CT-SIM and treatment, the same immobilization devices and a flat tabletop insert may be utilized used during MR scanning as during CT scanning.

[0031] With respect to image modification, to perform the image processing, the MR images may be registered to the CT image. This may be accomplished for example, by using a mutual information algorithm such as the one provided in the MATLAB® computer software distributed by Mathworks of Natick, MA. In cases where the images are of different resolutions and dimensions, the registered images may be normalized and rebinned to have uniform dimensions and spatial resolution. In an example, the registered images may be interpolated onto a common image grid to ensure uniform dimensions and resolution, and features (e.g., peaks) in the MR intensity value histograms may be normalized to preselected reference values to minimize intensity variation occurring between scans. Continuing with the exemplary image properties discussed above, as one possibility the Tl -weighted, T2-weighted and CT images may be linearly interpolated onto a grid matching the dimensions and resolution of the smaller bTFE image.

[0032] Because of the difficulty in distinguishing between bone and air in MR images, bones in the Tl -weighted MR image may be manually contoured using a treatment planning system, such as the Eclipse® system distributed by Varian Medical Systems of Palo Alto, CA. Both the bone contours and the MR images may then be exported into DICOM files that may, in turn, be imported into MATLAB®. Those transformations applied to the Tl -weighted image during the previously mentioned registration and interpolation steps may also be applied to the bone contours. In another example, a segmentation routine may utilize an ultra-short echo time (UTE) sequence, alone or in combination with the other available images, to identify bone tissue segments.

[0033] In order to reduce the computational expense of the method, a technique for using a

21 mask to ignore non-body voxels may be utilized, such as the technique given in Johansson as one possibility. For example, the Tl -weighted image may be converted into a binary image using Otsu thresholding. The image may then be filtered using a uniform 27-element cubic filter. The resulting image may then be eroded using a spherical structuring element and then thresholded again to create a binary image. A connected component analysis may be performed on this binary image, and the largest connected component may be assumed to be the body. The mask may accordingly consist of only voxels within the body. Using the mask, only these voxels within the may be used for optimizing the weights, generating the synthetic CT image, and quantitatively evaluating the synthetic CT images.

[0034] With respect to generation of a synthetic CT in the exemplary approach, a synthetic

CT intensity value, E_u for each voxel may be defined as a weighted sum of the acquired MR images (e.g., the three Tl, T2, and bTFE images; the three Tl, T2, and fluid attenuated inversion recovery images, etc.) and one derived MR image (inverseTl): E_i =∑w_k, M_k . (1)

k where M_k,i is the i voxel value from MR image k, and w_k is the weight assigned to voxel i from MR image k.

[0035] An objective for generating synthetic CT images may be to minimize the difference between the actual CT voxel intensity values, Ij, and the synthetic CT voxel values calculated using Eq. 1. For example, one approach may optimize the weight matrix W composed of elements Wk,u However, finding the optimum weight for each voxel in the MR images may be both time- consuming and sensitive to differences in patient geometry. Instead, the exemplary approach may assign the voxels to one of a set of region types, such as the five region types of: air, bone, soft tissue, muscle, and water. Using the aforementioned set of region types as an example mapping, only five sets of weights corresponding to each of the five body regions would need to be optimized. In another example, air voxels may be assigned a bulk value, whereas other voxels (e.g., bone, soft tissue, muscle, water) may be calculated using a weighted sum of the acquired and derived MR image. Thus, the formula for calculating the voxel intensity value in the synthetic CT images given in Eq. 1 may be modified to become:

^^■ =∑¾ (2)

k where r(i) is the region assigned to voxel i.

[0036] A set of criteria for assigning a voxel to a region may be developed from the acquired

MR images (Tl , T2, bTFE) and bone contours. Only those voxels assigned to a given region would be used in that region's weight optimization. In an example, bone classification may be performed using manually drawn contours, and very low intensity voxels not falling within the bone contours may be classified as air. As one possibility, a truth table may be constructed based on voxel intensity levels in the acquired MR images may be utilized to determine which region to assign the remaining voxels.

[0037] An alternate approach may optimize either the selection of MR sequences that would be used for the method or the MR criteria set that would be used to assign voxels to a region. In such an alternate approach, a bulk density value would then be assigned to all voxels within each region as in Rank¹⁹. However, to obtain the necessary criteria for dividing the patient into so many regions may require the use of several extra MR sequences, many of which may not be available (e.g., unavailable in sequences used as the standard of care). The additional imaging may also require more time than a prostate (or other) patient could reasonably be expected to sit through during one session.

[0038] With respect to image evaluation, the approach may optimize the image weights using a training set consisting of a subset of the patients used in the training (as one example, four patients out of the exemplary nine). For instance, to generate a synthetic CT for a given patient, a subset of the remaining patients may be used to optimize the weights following the procedure in Figure 2A. In an example, the weights may be initialized with random values, and used to generate a synthetic CT for the first patient in the optimization set. The weights may then be optimized by minimizing the Euclidean distance of calculated voxel value differences between the synthetic CT and the original CT image over a number of iterations (e.g., one hundred in an example). The resulting weights may then be used as initial weights for the next patient in the set. This process may be repeated until all members of the set have been used. The fewer the number of patients used during optimization, the greater the likelihood that the resulting weights may be tailored mainly to those patients rather than to the full patient population. Therefore, it may be important to determine the number of included patients that would optimize weights for the full population.

[0039] After performing the weight optimization, the MR images of a new patient may be loaded, and a synthetic CT image may be generated from these images using Eq. 2. The resulting synthetic CT image may then be rigidly registered to the original CT image using a mean square algorithm. After registration, a mean absolute error (MAE) metric may be used to evaluate the synthetic CT images. As one possibility, Hounsfield unit (HU) value differences between the original CT image and the synthetic CT image may be analyzed via MAE. The formula for MAE is given by: yj CT(i) - synCT{i

MAE (3)

N where N is the number of image voxels. The workflow for this procedure is displayed in Figure 1.

[0040] As an example, for each number of patients in the optimization set, weights may be randomly initialized five times, optimized, and used to generate five different synthetic CTs for each patient, resulting in the example error bars seen in Figure 2B. In the illustrated example, MAE was calculated for each synthetic CT. To evaluate dosimetric effects resulting from increasing the patient number used for optimization, treatment plans were created for a set of synthetic CTs. Each synthetic CT was generated using weights optimized with a different number of patients. As shown in Figure 2C, the prostate dose volume histograms show closer agreement with original dose as number of patients used during optimization increases.

[0041] With respect to treatment planning, the exemplary patients may be planned using intensity modulated radiation therapy (IMRT) with nine fields, and the anatomical structures used for the planning may be contoured on the CT simulation images. As one possibility, the planning may be performed using Eclipse®. For instance, the planned number of fractions may be varied between patients, using a prescribed dose of 1.8 Gy/fraction for each of the patients.

[0042] The generated synthetic CT images may be imported into Eclipse® to be used for planning. As the MR and CT scanning may be performed at different times and in different locations, there may often be large differences in the anatomy of key structures such as the bladder and rectum. Therefore, select structures may be contoured on the Tl -weighted MR images and then propagated onto the synthetic CT images. In the example of scans related to a prostate region, the selected structures may be the external body structure, left and right femora, bladder, rectum, and prostate. Because the geometry of the synthetic CT image may be generated to be substantially identical to that of the Tl -weighted image, no registration to the Tl -weighted image may be required and, therefore, no systematic errors may be introduced during the contour propagation step. Alternatively, the aforementioned selected structures may be contoured directly on the synthetic CT images, but that may not be ideal for a patient, as the soft tissue contrast of the MR images may be sacrificed to an extent during the generation of a synthetic CT image. The Dice similarity coefficient (DSC) may be used to compare the shape of the prostate contours produced using the planning CT and synthetic CT images. [0043] In order to compare treatment planning performance of the synthetic CT images to the planning CT images, each synthetic CT image volume may be rigidly registered to the planning CT volume. Following registration, the original treatment plan may be copied onto the synthetic CT image so that all beam angles, leaf positions, and monitor units used in the original plan may remain constant for retrospective calculations in the synthetic CT-based treatment plans. The beam angles, leaf positions, and monitor units may be recalculated to the isocenter based on the rigid registration of the synthetic and planning CT images. Dose volume histograms (DVH) of the prostate may then be compared in terms of the D99, max dose, and mean dose to the prostate. In addition, to reduce the dosimetric changes resulting from the differences in the shapes of the external body contours between the planning and synthetic CTs, the union of these two contours may be used as a common external body structure for both images. Voxels falling within the new structure, but located outside the intersection of the two original contours, may be assigned the bulk value of water.

[0044] As some exemplary results of generation of the synthetic CT, Figure 3 displays an image at various stages during the image construction workflow. Original MR images are illustrated in Figure 3A, followed by an intermediate CT image in Figure 3B with voxels assigned to one of the five regions. The Figure 3C illustrates a final synthetic CT image for a patient. The MAE results for the synthetic CT image relative to the planning CT image are summarized in Table 1.

Table 1: Synthetic CT evaluation

[0045] A source of error in the Synthetic CT may be misalignment of structures between the

CT and MR images, and the areas where misalignment may occur most frequently may include the body/air and bone/tissue interfaces. Overall, the average MAE over the entire volume for all synthetic CT images was ~75 HU with a standard deviation of 12 HU. Calculated HU value differences between synCTs and SIM-CTs were 2.0 ±8.1 HU and 11.9 ± 46.7 HU for soft tissue structures and femoral bones, respectively.

[0046] As an example of treatment planning, the D99, mean dose, and max dose to the prostate for both the original treatment plan and the synthetic CT-based treatment plan may be evaluated for each exemplary patient. The results are presented in Table 2.

CT (Gv) SvnCT (Gv) Difference (Gv) % Diff

I ncorrecled External

Prostate D99 69.96 ± 13.12 69.35 ± 14.26 0.62 ±2.04 2.31 ±3.27

IV slalc Max Doxo .35 .S2 1.14 o."S

Prostate Mean Dose 70.91 ± 13.26 71.16 ± 13.48 ^■0.24 ± O.i 1.22 ±0.58 Body Max Dose ^"3.45 : l.i.Sh ^"4.1Λ Ι4.Γ U.- l .(15 I .⁽in

Corrected External

Prostate D99 70.10 ± 11.59 70.22 ± 11.70 -0.12 ±0.64 0.75 ±0.35

I'awial.- Max ^"3.4s : 12.41 y : 1 .5^" -⁽».! I : d.5^" o.i .34

Prostate Mean Dose 71.36 ± 11.87 71.54 ± 11.98 -0.15 ±0.47 0.54 ±0.33

> Max l>o>e 74 12.X3 75.19 Ι ΛΜ -n.2y . u.45 o. ii .35

Table 2: Plan dose comparison for prostate

[0047] The largest error may be seen in the D99 of Patient 2 where the prostate contoured on the Tl weighted image was shifted closer to the rectum than the contour based on the planning CT. This shift caused a portion of the MR prostate to fall in a region with a large dose gradient. This change in prostate position was seen in that the DSC for Patient 2 was 0.58, which was the lowest of the patients and was notably lower than the average DSC of 0.72. Overall, even with Patient 2 included, the D99 to the patient differed on average by only 0.25% from the original plan while both the max and mean dose to the prostate and the max dose to the body remained within 0.2% of the original plan doses.

[0048] Figures 4A and 4B illustrate an exemplary comparison of the dose distributions for the planning-CT based treatment plan and the synthetic CT-based plan of an exemplary patient (labeled as Patient 11). Figure 4A includes dose distributions for a plan based on the simulation CT, while Figure 4B includes dose distributions based on the synthetic CT. Starting from the top and going clockwise, the images shown are the transverse, sagittal and coronal slices. As expected, the largest dosimetric differences between the two plans were caused by differences in the shapes of the external body contours between the planning and synthetic CT images. Using a common external body contour resulted in a reduction in both the average and standard deviation of the prostate dose differences between the two plans by greater than 50%.

[0049] Because of the large differences in bladder and rectal filling between the patients' CT and MR imaging sessions, only the prostate and bony structures are included in the representative dose-volume histogram (DVH) displayed in Figure 5, which shows good agreement in the prostate DVH between the original and synthetic CT-based plans. More specifically, Figure 5 depicts DVHs for the synthetic CT-based (solid lines) and planning CT-based (dotted lines) scans, for the prostate, left femora, and right femora of a representative patient.

[0050] Figure 6 illustrates an exemplary system 100 for the generation of synthetic CT images. The system 100 may include one or more image sources 102 from which CT and MR images may be acquired, a computing device 104 configured to generate the synthetic CT images using image data acquired from the image source 102, and a display device 106 configured to display the synthetic CT images for use by an operator. The system 100 may take many different forms and include multiple and/or alternate components and facilities. While an exemplary system 100 is shown in Figure 6, the exemplary components as illustrated is not intended to be limiting, and additional or alternative components and/or implementations may be used.

[0051] The image source 102 may include one or more devices from which images may be acquired for use by the system 100. The image source 102 may include one or more devices for generation of images, such as a Philips Brilliance Big Bore for acquisition of CT images and/or a Philips Healthcare Panorama open-platform MRI scanner for acquisition of MR images. Additionally or alternately, the image source 102 may include a database server or other data store from which previously captured CT and MR images may be acquired for use by the system 100.

[0052] The computing device 104 may be a personal computer, a portable computer, computing mainframe, or other available computing device with sufficient processing power. The computing device 104 may be configured to execute programs on one or more processors of the computing device 104, where the programs are stored on one or more memory devices included in or accessible by the computing device 104. The computing device 104 may further include network or other interface hardware configured to allow the computing device 104 to communicate with the image source 102 directly or indirectly such as over a communication network. The computing device 104 may also be connected to input hardware configured to receive input to be provided to the computing device 104. Exemplary input hardware may include buttons or other user controls for capturing input from a user of the computing device 104. The computing device 104 may also be connected to output hardware such as one or more display devices 106 to provide visual output, one or more speakers to provide audio output, and one or more haptic devices to provide haptic feedback to users of the device. The input hardware and output hardware may be used by the computing device 104 to provide a user interface between the computing device 104 and operators of the computing device 104.

[0053] Synthetic CT generation software may be one example of an application program stored on a memory of the computing device (e.g., as software, firmware, etc.). The computing device may include various modules that provide or support the various application functions or services. For example, as illustrated the computing device may include an image registration module 108, a region assignment module 110, a weight optimization module 112, a synthetic CT generation module 116, and a treatment planning module 118. Although one example of the modularization of the synthetic CT generation software is illustrated and described, it should be understood that the operations thereof may be provided by fewer, greater, or differently named modules. Specifics of the operation of these modules are discussed below with reference to the workflow Figures 7 and 8.

[0054] With reference to the exemplary validation workflow illustrated in Figure 7, the image registration module 108 may be configured to acquire MR and CT images, rigidly register the images, and resample the images to have the same resolution and/or dimensions. As one possibility the image registration module 108 may perform registration by employing a mutual information algorithm such as the one provided in the MATLAB® computer software distributed by Mathworks of Natick, MA. The image registration module 108 may further be configured to utilize an image processing library to normalize and rebin the images to have uniform dimensions and spatial resolution. For instance, the image registration module 108 may interpolate the images onto a common image grid to ensure uniform dimensions and resolution, and may normalize features (e.g., peaks) in the MR intensity value histograms to preselected reference values to minimize intensity variation occurring between scans.

[0055] The region assignment module 110 may be configured, in one aspect, to mask the received normalized and rebinned images to remove non-body voxels. For example, the region

21 assignment module 110 may use a technique such as the technique given in Johansson and discussed in detail above. As another aspect, the region assignment module 110 may be configured to use CT and MR image criteria to assign image voxels to regions. Region assignment may be performed by the region assignment module 110 using various techniques, such as one or more of those discussed above¹¹'¹²'¹⁸'¹⁹. As one possibility, a region for each voxel may be determined as being one of a set of region types, such as the five region types of: air, bone, soft tissue, muscle, and water. With respect to bone identification, in an example a segmentation routine may utilize an ultra-short echo time (UTE) sequence, alone or in combination with the other available images, to identify bone tissue segments.

[0056] The weight optimization module 112 may be configured to receive the region- identified images from the region assignment module 110, and to generate optimized weights 114 for use in generation of synthetic CT images. As one aspect, the weight optimization module 112 may utilize the images to optimize weights 114 for each of the regions as discussed above with respect to Eq. 2. An example procedure for weight optimization and optimizing a number of patients used during weight optimization is discussed above with respect to Figure 2A. For example, the weight optimization module 112 may optimize voxel intensity value weights 114 for each of the region types of voxels, to minimize a difference between (i) actual voxel intensity values of the CT image and (ii) synthetic CT voxel values determined according to weighting the voxels of the MR images according to voxel intensity value weights of corresponding region types. The resultant weights 114 may be stored in a database or other data store for use in creation of synthetic CT images.

[0057] The synthetic CT generation module 116 may be configured to use the optimized weights 114 and MR images to generate a synthetic CT image. For example, the synthetic CT generation module 116 may utilize the voxel intensity value weights 114 for each of the region types of voxels and the voxels of the MR images to compute synthetic CT voxel values to include in a synthetic CT image. As another aspect, the synthetic CT generation module 116 may rigidly register the resulting synthetic CT image to the original CT image using a mean square algorithm, and may utilize a mean absolute error (MAE) metric to evaluate the synthetic CT images (e.g., according to Eq. 3 as discussed above). In some cases, if the generated synthetic CT image is within an acceptable threshold error amount from the planning CT, then the optimized weights 114 may be accepted by the system 100 and stored for later use in generation of synthetic CT images.

[0058] Additionally or alternately, the treatment planning module 118 may be configured to generate a first treatment plan using the planning CT and a second treatment plan using the synthetic CT. As one possibility, the planning may be performed using Eclipse®. The treatment planning module 118 may be further configured to compare differences between the plan created from planning CT and the plan created from the synthetic CT. For example, if the plans differ within an acceptable threshold (e.g., 1%, 5%), then the optimized weights 114 may be accepted as producing synthetic CTs useful for treatment applications.

[0059] With reference to the exemplary clinical workflow illustrated in Figure 8, the image registration module 108 may be configured to acquire MR images for a patient, without requiring CT images of the patient. The image registration module 108 may rigidly register the images, and, if the images differ in resolution and/or dimensions, resample and/or interpolate the images to have the same resolution. The region assignment module 110 may be configured to mask and assign regions to the images, as discussed above. The synthetic CT generation module 116 may be configured to retrieve previously-computed optimized weights 114 (such as determined above with respect to the workflow of Figure 7) and to use the weights 114 to generate a synthetic CT image from the acquired MR images. For example, the synthetic CT generation module 116 may utilize the previously-generated voxel intensity value weights 114 for each of the region types of voxels, and the voxels of the MR images, to compute synthetic CT voxel values to include in a synthetic CT image.

[0060] The treatment planning module 118 may be configured to generate a treatment plan using the synthetic CT image (e.g., using Eclipse® as discussed above). Accordingly, dose calculations may be performed using MR images of a patient, without requiring the acquisition of CT scans in the dose calculation workflow. Moreover, the treatment planning module 118 may be further configured to generate a digitally reconstructed radiograph utilizing the synthetic CT image, further providing for patient positioning, also without requiring the acquisition of CT scans.

[0061] Accordingly, the above-described system 100 and computer-implemented methods utilize an approach for generating a synthetic CT image from a set of MR images. As discussed in detail above, the disclosed approach provides good agreement overall between generated synthetic CT images and their corresponding planning CT images, both qualitatively and in terms of MAE. The breakdown of HU difference by organ type illustrates the agreement between the CT images of within 1% in soft tissue and within 5% in bony structures. The largest errors may be identified at the body/air and bone/tissue interfaces. To minimize errors in the weight 114 optimization, the approach may further include achieving geometric agreement between the MR and CT images, for example by minimizing a time gap between MR and CT imaging that may result in noticeable differences in internal patient geometry. Because of differences in bladder and rectal filling, air pockets in the CT and MR images may be misaligned, which may cause MAE to be relatively larger in voxels assigned to the air region.

[0062] Moreover, beam angles, leaf positions, and monitor units from the CT simulation- based treatment plans may be kept constant for synthetic CT-based treatment plans in order to eliminate variability that may result from using different optimizations. The prostates contoured using the synthetic CT images in the exemplary set of patients were generally smaller than the original CT images, and only about half showed a fairly good match between the two contours in terms of DSC. It may be difficult to differentiate how much of the error may be due to differences in synthetic CT values calculated using the exemplary approach and how much may be caused by other factors such as different effective radiation path lengths to the prostate. Nevertheless, using the synthetic CT images, the average D99, max dose, and mean dose to the prostate as well as maximum dose to the body remained substantially within 0.5% of their original planned values.

[0063] Thus, a system and computer-implemented method may utilize the disclosed approach to generate synthetic CT images from a set of standard MR images. Treatment plans generated using these synthetic CT images with the parameters that had been used in the CT simulation-based plans provided dosimetric results for the prostate that differ by an acceptable margin (e.g., in the exemplary set of patients less than 1%) from those produced in the original treatment plans. The synthetic CT images generated using this approach may be used for dose calculation in MR-only treatment planning, without requiring the use of specialized MR sequences such as UTE.

[0064] It should be noted that the workflows of Figures described herein are exemplary only, and that the functions or steps of the method could be undertaken other than in the order described and/or simultaneously as may be desired, permitted, and/or possible.

[0065] While exemplary embodiments are described above, it is not intended that these embodiments describe all possible forms of the invention. Rather, the words used in the specification are words of description rather than limitation, and it is understood that various changes may be made without departing from the spirit and scope of the invention. Additionally, the features of various implementing embodiments may be combined to form further embodiments of the invention.

REFERENCES

1. Beavis AW, Gibbs P, Dealey RA, et al. Radiotherapy treatment planning of brain tumours using MRI alone. Br J Radiol 1998; 71 : 544-548.

2. Stanescu T, Hans HS, Pervez N, et al. A study on the magnetic resonance imaging (MRI)- based radiation treatment planning of intracranial lesions. Phys Med Biol 2008; 53: 3579- 3593.

3. Buhl SK, Duun-Christensen AD, Kristensen BH, et al. Clinical evaluation of 3D/3D MRI- CBCT automatching on brain tumors for online patient setup verification - a step towards MRI-based treatment planning. Acta Oncol 2010; 49: 1085-1091.

4. Greer PB, Dowling JA, Lambert JA, et al. A magnetic resonance imaging-based workflow for planning radiation therapy for prostate cancer. Med JAust 2011; 194: S24-S27. Kapanen M, Collan J, Beule A, et al. Commissioning of MRI-only based treatment planning procedure for external beam radiotherapy of prostate. Magn Reson Med 2013; 70: 127-135. Graves EE, Pirzkall A, Nelson SJ, et al. Registration of magnetic resonance spectroscopic imaging to computed tomography for radiotherapy treatment planning. Med Phys 2008; 28: 2489-2496. Brunt JN. Computed tomography-magnetic resonance image registration in radiotherapy treatment planning. Clin Oncol 2010; 22: 688-697. Tanaka H, Hayashi S, Ohtakara K, et al. Usefulness of CT-MRI fusion in radiotherapy planning for localized prostate cancer. J Radiat Res 2011; 52: 782-788. Devic S. MRI simulation for radiotherapy treatment planning. Med Phys 2012; 39: 6701- 6711. Karlsson M, Karlsson MG, Nyholm T, et al. Dedicated magnetic resonance imaging in the radiotherapy clinic. Int J Radiat Oncol, Biol, Phys 2009; 74; 644-651. Jonsson JH, Karlsson MG, Karlsson M, et al. Treatment planning using MRI data: an analysis of the dose calculation accuracy for different treatment regions. Radiat Oncol 2010; 5: 62. Chen L, Price RA Jr., Wang L, et al. MRI-based treatment planning for radiotherapy: Dosimetric verification for prostate IMRT. Int J Radiat Oncol, Biol, Phys 2004: 60: 636-647. Wang D, Strugnell W, Cowin G, et al. Geometric distortion in clinical MRI systems: Part II: correction using a 3D phantom. Magn Reson Imaging 2004; 22: 1223-1232. Langlois S, Desvignes M, Constans JM, et al. MRI geometric distortion: a simple approach to correcting the effects of nonlinear gradient fields. J Magn Reson Imaging 1999; 9: 821- 831. Karger CP, Hoss A, Bendl R, et al. Accuracy of device-specific 2D and 3D image distortion correction algorithms for magnetic resonance imaging of the head provided by a manufacturer. Phys Med Biol 2006; 51 : N253-261. Doran SJ, Charles-Edwards L, Reinsberg SA, et al. A complete distortion correction for MR images: I. Gradient ward correction. Phys Med Biol 2005; 50: 1343-1361. Karotki A, Mah K, Meijer G, et al. Comparison of bulk electron density and voxel-based electron density treatment planning. JAppl Clin Med Phys 2011; 12: 97-104. Lee YK, Bollet M, Charles-Edwards G, et al. Radiotherapy treatment planning of prostate cancer using magnetic resonance imaging alone. Radiother Oncol 2003; 66: 203-216. Rank CM, Tremmel C, Hunemohr N, et al. MRI-based treatment plan simulation and adaptation for ion radiotherapy using a classification-based approach. Radiat Oncol 2013; 8: 51. Dowling JA, Lambert J, Parker J, et al. An atlas-based electron density mapping method for magnetic resonance imaging (MRI)-alone treatment planning and adaptive MRI-based prostate radiation therapy. Int J Radiat Oncol, Biol, Phys 2012; 83; e5-el 1. Johansson A, Karlsson M, Nyholm T. CT substitute derived from MRI sequences with ultrashort echo time. Med Phys 2011; 38: 2708-2714. Jonsson JH, Johansson A, Soderstrom K, et al. Treatment planning of intracranial targets on MRI derived substitute CT data. Radiother Oncol 2013; 108: 118-122.

Claims

WHAT IS CLAIMED IS:

1. A system comprising:

a computing device configured to

acquire, of a treatment area of a patient, a planning computed tomography image and at least one magnetic resonance image;

assign voxels of the at least one magnetic resonance image to tissue region types from a set of region types; and

optimize voxel intensity value weights for each of the region types to minimize a difference between (i) actual voxel intensity values of the planning computed tomography image and (ii) synthetic computed tomography voxel values determined by weighting the voxels of the region types of the at least one magnetic resonance image according to the voxel intensity value weights of the corresponding region types.

2. The system of claim 1, wherein the computing device is further configured to:

generate a synthetic computed tomography image according to the at least one magnetic resonance image and the voxel intensity value weights; and

validate the synthetic computed tomography image against the planning computed tomography image.

3. The system of claim 1, wherein the set of region types includes tissue region types of: air, bone, soft tissue, muscle, and water.

4. The system of claim 1, wherein the computing device is further configured to optimize the voxel intensity value weights according to:

Ei is a synthetic computed tomography image intensity value for voxel i, Mk_,i is an intensity of voxel i of magnetic resonance image k,

r(i) is the region type assigned to voxel i, and

Wk_,i is a weight assigned to voxel i from magnetic resonance image k according to the region type r(i) assigned to voxel i.

5. The system of claim 1 , wherein the planning computed tomography and magnetic resonance images include previously captured images retrieved from an image data store.

6. The system of claim 1, wherein the computing device is further configured to rigidly register the planning computed tomography with the at least one magnetic resonance image.

7. The system of claim 1, wherein the computing device is further configured to at least one of (i) interpolate the images onto a common image grid to ensure uniform dimensions and resolution, and (ii) normalize features in the at least one magnetic resonance image to preselected reference values to minimize intensity variation occurring between scans.

8. The system of claim 1, wherein the at least one magnetic resonance image includes images captured using at least two of: a Tl -weighted sequence, a T2-weighted sequence, and a balanced turbo field echo sequence.

9. The system of claim 1, wherein the at least one magnetic resonance image includes images captured according to a standard of care for the treatment area of the patient being imaged.

10. A system, comprising :

a computing device configured to

assign voxels of magnetic resonance images of a treatment area of a patient to tissue region types from a set of region types; and

generate a synthetic computed tomography image of the treatment area of the patient according to the magnetic resonance images and optimized voxel intensity value weights for each of the tissue region types, the optimized voxel intensity value weights computed according to a minimization of a difference between (i) actual voxel intensity values of training computed tomography images and (ii) magnetic resonance voxel intensity values of training magnetic resonance images weighted according to the optimized voxel intensity value weights.

11. The system of claim 10, wherein the computing device is further configured to generate a treatment plan based on the synthetic computed tomography image.

12. The system of claim 10, wherein the computing device is further configured to generate a digitally reconstructed radiograph based on the synthetic computed tomography image.

13. A computer-implemented method comprising :

assigning voxels of magnetic resonance images of a treatment area of a patient to tissue region types from a set of region types; and

generating a synthetic computed tomography image of the treatment area of the patient according to the magnetic resonance images and optimized voxel intensity value weights for each of the tissue region types, the optimized voxel intensity value weights computed according to a minimization of a difference between (i) actual voxel intensity values of training computed tomography images and (ii) magnetic resonance voxel intensity values of training magnetic resonance images weighted according to the optimized voxel intensity value weights.

14. The method of claim 13, further comprising generating a treatment plan for the patient based on the synthetic computed tomography image, without requiring a planning computed tomography image.

15. The method of claim 14, further comprising validating the synthetic computed tomography image against the planning computed tomography image.

16. The method of claim 14, wherein the planning computed tomography and magnetic resonance images include previously captured images retrieved from an image data store.

17. The method of claim 13, further comprising generating a digitally reconstructed radiograph based on the synthetic computed tomography image.

18. The method of claim 13, wherein the set of region types includes tissue region types of: air, bone, soft tissue, muscle, and water.

19. The method of claim 13, further comprising optimizing the voxel intensity value weights according to:

Ei is a synthetic computed tomography image intensity value for voxel i, M_kj is an intensity of voxel i of magnetic resonance image k,

r(i) is the region type assigned to voxel i, and

Wk,i is a weight assigned to voxel i from magnetic resonance image k according to the region type r(i) assigned to voxel i.

20. The method of claim 13, further comprising at least one of (i) interpolating the images onto a common image grid to ensure uniform dimensions and resolution, and (ii) normalizing features in the at least one magnetic resonance image to preselected reference values to minimize intensity variation occurring between scans.