T1/T2*-weighted MRI provides clinically relevant pseudo-CT density data for the pelvic bones in MRI-only based radiotherapy treatment planning

Abstract

Background and purpose. In radiotherapy (RT), target soft tissues are best defined on MR images. In several cases, CT imaging is needed only for dose calculation and generation of digitally reconstructed radiographs (DRRs). Image co-registration errors between MRI and CT can be avoided by using MRI-only based treatment planning, especially in the pelvis. Since electron density information can not be directly derived from the MRI, a method is needed to convert MRI data into CT like data. We investigated whether there is a relationship between MRI intensity and Hounsfield unit (HU) values for the pelvic bones. The aim was to generate a method to convert bone MRI intensity into HU data surrogate for RT treatment planning. Material and methods. The HU conversion model was generated for 10 randomly chosen prostate cancer patients and independent validation was performed in another 10 patients. Data consisted of 800 image voxels chosen within the pelvic bones in both T1/T2*-weighted gradient echo and CT images. Relation between MRI intensity and electron density was derived from calibrated HU-values. The proposed method was tested by constructing five “pseudo”-CT series. Results. We found that the MRI intensity is related to the HU value within a HU range from 0 to 1400 within the pelvic bones. The mean prediction error of the conversion model was 135 HU. Dose calculation based on the pseudo-CT images was accurate and the generated DRRs were of good quality. Conclusions. The proposed method enables generation of clinically relevant pseudo-CT data for the pelvic bones from one MRI series. It is simpler than previously reported approaches which require either acquisition of several MRI series or T2* maps with special imaging sequences. The method can be applied with commercial clinical image processing software. The application requires segmentation of the bones in the MR images.

Magnetic resonance imaging (MRI) is widely used in radiotherapy (RT) treatment planning due to its superior soft tissue contrast when compared with computed tomography (CT). This means that MR and CT images should be co-registered to enable dose calculation and the generation of the digitally reconstructed radiographs (DRRs). Because of the different patient setup in these images, this increases the geometrical uncertainty of the position of the target and normal tissues [1]. MRI simulators are MRI devices which have a wide or open bore, flat imaging table top and a laser system suitable for skin marking. They can be used for MRI-only based RT treatment planning since they allow using the same fixation devices and patient setup that are used in RT treatment [2,3]. The advantage is that image co-registration errors are omitted. This approach, however, may result in dose calculation errors of 1–2% in the pelvic region due to the lack of electron density information which is derived directly from the CT images.

Hounsfield unit (HU) information has been assigned for the MR images in literature. The simplest approach used a bulk density for the whole bone volume and assumed everything else as water equivalent tissue [2,4,5]. In the more elaborated techniques, density values or bone characteristics have been derived from T2* (apparent transverse relaxation time) maps [6], by determining T2 (transverse relaxation time) spectrum of bone [7] or by combining information from several MRI series [8]. These approaches have utilized several image series and/or special imaging sequences based on ultra short echo times (TE ≤ 0.1 ms). These sequences are still prototypes and their application requires extra imaging time. The essential characteristics of bones have also been evaluated simply by using anatomical T1/T2*-weighted gradient echo MRI [9].

In this study, we investigated whether there is a relation between MRI intensity and the HU-value within the pelvic bones. The aim was to find a method that could provide HU data from the MRI to generate a pseudo-CT series which could be then used in clinical RT treatment planning. T1/T2*-weighted anatomical gradient echo imaging was chosen because it provides essential characteristics of the bones and especially because T2* time correlates with bone density [9,10]. Segmentation of the bones seemed possible in these images because the intensity range in the bones was clearly different from that of the adjacent muscles. The relation between MRI intensity and the relative electron density was obtained from calibrated HU-values. The results were validated in an independent patient group. We introduced the method in a form that can be applied with commercial clinical image processing software with limited processing capability. This was possible because the conversion model was constructed and validated by choosing the voxels randomly by hand and because the pseudo CT series were generated based on discrete MRI intensity values. The pseudo-CT series were constructed for five patients to test the feasibility of the method for RT treatment planning. We investigated the achieved dose calculation accuracy and constructed sample DRRs.

Material and methods

Equipment and imaging parameters

The MRI was performed with a 1.5 T scanner GE Optima MR450w (GE Medical Systems, Milwaukee, WI, USA). The scanner has a flat table top and a wide bore (diameter 70 cm) enabling the use of similar patient immobilization devices as in external RT. A T1/T2*-weighted (spoiled steady state) three-dimensional (3D) gradient echo sequence LAVA Flex® was chosen. This is a dual echo sequence which provides out-phase and in-phase images for water and fat signals acquired at echo times (TE) of 2.1 and 4.2 ms, respectively. Fat and water images can be calculated from these images. Repetition time (TR) was 6.8 ms, flip angle 15^o, reconstruction slice thickness 2.4 mm (interpolated to 1.2 mm by zerofilling) and matrix size 428x448 (FOV 42 × 42 cm). Bandwidth (BW) was 90.9 kHz for the entire FOV. Only in-phase images were used (except for segmentation of bone) because bone marrow contains both fat and water. The images were corrected for intensity inhomogeneity using PURE® (based on calibration series) and for geometrical distortion using Grapwarp 3® both of these software provided by the vendor. The imaging area in the superior-inferior direction covered the prostate and seminal vesicles with margins ranging from 2 to 5 cm. The coverage in the anterior-posterior and left-right directions covered the whole body contour to enable dose calculation. Imaging was carried out using an eight channel phased array cardiac coil with the posterior part of the coil placed below the flat table top and the anterior part supported by two plastic semi arcs constructed in our hospital. This system setup enables the use of the same patient setup as in treatment. A drawback is a lower signal to noise ratio when compared with purely diagnostic imaging due to the greater distance between the coil and patient skin. For this kind of approach, however, the image quality has been considered sufficient [11].

The CT imaging was carried out with GE LightSpeed RT scanner (4-slice, GE Medical Systems, Waukesha, WI) at 120 kVp with slice thickness of 2.5 mm. The HU-values were calibrated to electron density values using the RMI-phantom (Gammex, Middleton, WI). The image processing was performed using MIM 5® workstation (MIM Software Inc., Cleveland, OH, USA). Treatment planning was carried out with a grid size of 2.5 mm using Eclipse® workstation (Varian Medical Systems, Palo Alto, CA, USA).

Patients and generation of data points

Twenty randomly chosen prostate cancer patients (age 59 to 82 years, mean of 70 years) prescribed to external RT of the prostate and seminal vesicles were MRI and CT imaged according to our clinical routine procedure. The patients were informed about the study and verbal agreement was obtained for the inclusion. The study did not have any effect on the treatment of these patients but one or two extra MRI series were acquired (extra imaging time < 6 minutes per series). Five patients were imaged also with slightly different values of the TE = 5.2 ms, TR = 8.3 ms and BW = 62.5 kHz to estimate the effect of imaging parameters. The axial MR and CT images were first co-registered according to bony landmarks. Ten patients were included in the generation of the model and another 10 for independent validation. We chose 40 voxels randomly within the pelvic bone area for each patient. The corresponding HU and MRI intensity values of a voxel gave a data point for the generation and validation of the model. This kind of model construction is simple and it is possible with clinically accepted image processing software with limited processing capability. The superior-inferior (S-I) coverage was from the femoral bone/trochanter major to the middle of the iliac bone. The voxels were chosen uniformly within the segmented bone volume and it was confirmed that they cover the whole observed HU range from about 0 up to 1400. The voxels were located at different kind of bone structures such as cortical bone, trabecular bone and bone marrow.

The intensity in MR images may be inhomogeneous due to positional variations in the sensitivity of the surface coil and magnetic field imperfections caused by the system or the patient. Also the skin distance of a surface coil and gains in signal acquisition may differ between the patients. Even though the inhomogeneity correction (PURE®) equalizes intensity variations within a patient, overall intensity levels may still differ between patients. To minimize the effect of this, the MRI intensity values of the voxels were normalized to the average intensity of a muscle that is located closest to the bone (see Figure 1). For this purpose, eight small ROIs were drawn in different image slices excluding contribution of fat. The ROIs were uniformly located within same S-I coverage as the voxels. The normalized data was multiplied by 150 (mean normalization factor) to achieve realistic MRI intensity level.

Figure 1. Example of a small normalization ROI (upper white circle) chosen within a muscle that is located next to the bone. The edge of the bone (lower white contour) is clearly visible enabling segmentation of the bone edge in the MR image (LAVA Flex® in-phase image).

Modeling of the relation between MR intensity and HU-value inside the bone

Small solid fragments of the bone generate susceptibility gradients within pores that are located next to them. These pores contain both water and fat. The gradients have an effect on MRI intensity through changes in T2*. The relationship is very complex and depends on the volume fraction of these solid parts, their number, thickness and separation [12,13]. We generated a simplified model to give a theoretical basis for the relation between the MRI intensity and HU-value (or electron density). Signal intensity s in gradient echo imaging is proportional to s = A(T1)·P·exp(-TE·R2*), where R2* = 1/T2*, P is proton (spin) density and A is a T1-dependent factor determined by the flip angle, TR and level of spoiling of transversal magnetization. The visible MRI signal comes from different soft tissue components (such as collagen, collagen-bound water and pore water) and is the integral over different values of T1, P and T2* within a voxel. When normalizing to the signal from a muscle (water) we obtain

where ρ(x,y,z)·dx·dy·dz is the fraction of a component having T1, P and T2* within a voxel and x, y and z are space coordinates.

The aim was to investigate the overall relation between MRI intensity and HU-values and to find out whether the Equation 1 can be applied in a simple way and with sufficient accuracy for clinical RT planning. In practice, the inverse form of this formula is preferable (i.e. HU-values should be determined from the MRI intensity values). Since mathematic forms of all the factors in Equation 1 are not known, the inversion was done by a series expansion which results in a simple polynomial model of

where A, B and C are the fitting parameters and s is MRI intensity (normalized or unnormalized). The fitting parameters are complex functions of the distributions of T1, T2*, P and tissue component fractions within a voxel. The order for the polynomial fit was chosen based on minimum value of Akaike information criterion AIC [14] calculated as

where P is the number of fitting parameters, N is the number of data points and Σ(Δy)² is the sum of the squares of the fitting errors [Dy = HU_measured-HU(s)].

Each fitting parameter value is given as the obtained value (± 1 SD). A 95% confidence limit for a parameter is ± 2 SD. The Pearson correlation coefficient (ρ) was used to quantify linear correlation between two sets of data and to approximate overall linear correlation for slightly non-linear relations.

Generation of pseudo-CTs

Five pseudo-CT series were generated to demonstrate that the obtained model could provide appropriate data for dose calculation and DRR generation. Five patients were chosen randomly from the validation group. Pelvic bones were manually contoured according to the LAVA Flex in-phase, out-phase and water-images. The pseudo-CTs were generated with clinically accepted software with limited processing capability and, therefore, a voxel-based conversion was not possible. Instead, the voxels within the segmented bone structure were divided into 16 subgroups based on their MRI intensity values using an intensity ranging tool. According to the conversion model, the MRI intensity ranges of the subgroups corresponded to the following HU ranges: below 0, 0...100, ..., 1300...1400, and above 1400. The volume that had HU values above 1400 was below 1 ml. Therefore, the upmost subgroup was ‘above 1400’. Since only very small negative HU values were observed for the patients (HU > 210), the lowest subgroup was ‘below 0’. Each pseudo-CT was generated by converting mean MRI intensities of the subgroups into HUs by using the obtained model. The rest of the body was assumed water equivalent (HU = 0). Rectum and skin contour in MRI and CT were corrected to be similar to minimize their effect on dose calculation. DRRs were generated in a treatment planning system from these pseudo-CT images.

A simple four-field ‘box’ treatment plan was calculated for 6 and 15 MV beams to demonstrate the effect of the method on dose calculation. PTV included the prostate and seminal vesicles. Dose calculation differences between a pseudo-CT and actual CT were determined for five different versions of a pseudo-CT: i) by assuming the whole body as water equivalent; ii) by using a bulk HU value for the bones which corresponded to the average obtained for the whole group of the patients; iii) by assuming an individually determined bulk HU values for the bones converted from the MRI; iv) by determining regional HU values converted from the MRI; and v) as in iv) but the HU values for the body outside the pelvic bones were corrected to correspond to that of actual CT (‘fat correction’).

Results

Relation between HU-value and relative electron density

The electron density values relative to water of 1.000, 1.039, 1.116, 1.142, 1.285, 1.473 and 1.707 of the RMI phantom resulted in measured HU-values of 10, 12, 97, 181, 516, 931 and 1401, respectively. Pearson correlation coefficient of these values was ρ = 0.996 (p < 0.0001). The AIC method suggested a linear model. The fitting parameters were A = 1.037 (± 0.013) and B = 0.00048 (± 0.00002) 1/HU. The obtained linear model was used to convert HU- values into relative electron density values.

Modeling the relation between MRI-intensity and HU-value

The correlation between MRI-intensity and HU-value was ρ = -0.862 (p < 0.0001) (total of 400 data points). The AIC method suggested a second order fit. Mean, SD and maximal error of the fit were 131, 115 and 512 HU, respectively. The fitting parameters were A = 1947 (± 71) HU, B = 211.39 (± 0.79) HU/MR_{intensity unit} and C = 0.0169 (± 0.0022) HU/MR²_{intensity unit}. The slope of a linear fit was 25.04 HU/MR_{intensity unit} demonstrating a mean steepness of the non-linear fit. The model fit is illustrated in Figure 2.

Figure 2. Relation between MRI intensity (in-phase image) and HU number obtained for 400 randomly chosen voxels within the pelvic bones. The second order polynomial fit is expressed as a solid line.

Quite similar model fit was obtained without data normalization but systematic errors for an individual patient raised from 85 HU up to 154 HU. Normalization to fat resulted in poorer fit and the correlation of data was below 0.4.

The derived absolute MRI intensity values were lower for the different set of imaging parameters. This demonstrates that the parametrization for the relation between MRI intensity and HU depends on the imaging parameters. The correlation between the MRI intensities obtained for the two parameter sets was almost perfect 0.992 (p < 0.0001).

Validation group

The correlation between the calculated and measured HUs was ρ = 0.913 (p < 0.0001) for the 400 data points. The mean, SD and maximal differences were 135, 121 and 578 HU, respectively. Prediction error of the HU did not correlate with HU-level (ρ = 0.291) within 0 ≤ HU ≤ 1400. The linear fit of the calculated and measured HU values resulted in the parameter values of A = 25 (± 18) HU and B = 0.96 (± 0.03) suggesting that there was no systematic error between the measured and calculated values. The relationship between the measured and calculated HUs is presented in Figure 3. The range of mean HU values of the pelvic bones varied between 225 and 473 with a mean ± SD of 337 ± 63 HU.

Figure 3. Relation between the measured and calculated HU-values for 400 data points obtained from the validation group. Straight solid line with slope equal to one illustrates ideal calculation.

Generated pseudo-CT images

An example of the generated DRR images is presented in Figure 4. The dose calculation results are summarized in Table I. The proposed method increased overall calculation accuracy for all of the five patients. The individual bulk HUs for the pelvic bones were calculated within ± 23 HU.

Figure 4. DRR image generated from A) actual CT image and B) pseudo-CT image with HU-values inside the bones converted from the MRI. The field size is equal in the both images. Lateral images were chosen to demonstrate clear visibility of bone edges through whole transversal section of the pelvic bones.

Table I. Maximal dose calculation differences (in percentages) between pseudo2CT and actual CT for five randomly chosen patients (equal MU counts were applied for both images).

		HU information used inside the pelvic bones
Energy	Fraction of PTV	Water only	Bulk value of 337 corresponding group average	Individual bulk value obtained from MRI	Regional HU obtained from MRI	Regional HU from MRI with body corrected^a
6 MV	V95%	1.8%	− 1.9%	− 1.1%	− 1.2%	− 0.6%
	V50%	2.0%	− 1.9%	− 0.9%	− 1.0%	− 0.4%
	V5%	1.7%	− 1.9%	− 1.2%	− 1.2%	− 0.5%
	all points within	± 3.9%	± 3.2%	± 2.2%	± 2.4%	± 1.3%
15 MV	V95%	1.9%	− 1.4%	− 0.9%	− 0.9%	− 0.4%
	V50%	1.6%	− 1.4%	− 0.7%	− 0.8%	− 0.4%
	V5%	1.5%	− 1.2%	− 0.8%	− 0.8%	− 0.4%
	all points within	± 3.7%	± 2.7%	± 2.2%	± 2.2%	± 1.2%

^aHU value for the body outside the bones was corrected for fat content (on the average to about 220) to demonstrate differences related mainly to the conversion of HU values inside the bones.

Discussion

We have demonstrated that there is relation between the intensity of T1/T2*-weighted gradient echo MR image and HU number (or relative electron density) within the pelvic bones. This relation provides patient specific bone density information both for dose calculation and the generation of verification images (DRRs) for MRI-only based RT treatment planning. Such data have earlier been derived from T2* maps [6], T2 spectrum of bone [7] or from several MRI series [8] based on special ultra short echo time (UTE) sequences. To the best of our knowledge, this is the first report of generation of bone density data from only one anatomical MRI series. The proposed method is simple and that it utilizes the same MRI series that can be used for the segmentation of bones, dose calculation and DRR generation. There is no need for multiple MRI series to replace CT. The method can be applied with commercial clinical image processing systems with limited processing capability. The provided conversion model was validated inside the pelvic bones and it should not be adopted for use outside the segmented bone structures.

The mean prediction error obtained is equal to value of 137 HU reported for a voxel based method utilizing four image series [8]. In this light, this study presents an equal alternative for HU conversion. The mean prediction and fitting errors, and dose calculations were considered acceptable for a HU range from 0 to 1400. The MRI signal from dense cortical bone (HU ≥ 1200) was dark [15] and just about 10 units above the noise level. Even though the effect of noise was most pronounced with the lowest MRI intensities, the model seemed to average the data properly. The large maximal prediction and fitting errors were most probably due to slight image co-registration errors, noise and partial volume effects. These factors may depend on the voxel size.

The results were compared with the conventional bulk density approach since it utilizes data obtained purely for the bones and it does not adjust HU based on dose calculation or on density values of body organs outside the bones (e.g. fat and muscle). Similar bulk HU values have been reported for several groups of patients (density from 1.19 to 1.22 g/cm³ corresponding to HU range from 288 to 333) [2,4,5]. Mean prediction error of the relative electron density was 0.06 (corresponding to 135 HU) which should not normally result in a dose calculation error greater than about 1% for photons. By minimizing the effect of fat we demonstrated that the uncertainty of dose calculation related to the HU conversion was acceptable. We are aiming to develop a method for regional fat correction based on the LAVA Flex® images. In addition, we are evaluating calculation accuracy achievable by using the presented method in detail for different treatment sites and techniques.

The normalization of the MRI intensity was observed important to minimize systematic errors in HU conversion. This is a drawback but no better method was found to correct for variations in the overall intensity level between patients. The normalization was also tested to one ROI but poorer results were obtained than by normalizing to the average of several ROIs. This is most likely explained by image noise and small local intensity variations left after the inhomogeneity correction (PURE®). Since PURE® is based on calibration images, also potential patient movement during the imaging may result in uncorrected intensity variations. The use of larger normalization ROIs can not solve the problem since the ROI area should be kept small in order to minimize the contribution of fat and other tissues adjacent to the muscle. As a conclusion, the normalization should be done to the average of several small ROIs (at least eight) drawn within the muscles closest to the bones on different image slices (see Figure 1).

The MRI intensity values could not be directly converted into electron densities because of the lack of a MRI compatible phantom. The uncertainty related to conversion of electron density from HU is regarded clinically acceptable and, therefore, we conclude that the presented conversion of the MRI intensity into electron density through HUs is acceptable. The relation between the MRI intensity and HU seemed to exist at least with this type of MRI imaging sequence within a moderate range of the imaging parameter values. The presented model should be tuned for the individual MRI scanner, imaging sequence and parameters and also for the CT scanner and CT imaging protocol. Changes in system performance (e.g. magnetic field shimming) due to aging may require updating the model and its validity should be established regularly.

The presented method can be considered sensitive to bone segmentation errors. The intensity conversion, however, increased the intensity of cortical bone area (and thus enhancing the bone edge) and minimized the effect of such errors. The bone edge was observed to be consistent within about 1 mm in MRI and CT. In addition, the location of the bone edge was observed to be accurate within 1 mm (about voxel size) in an on-going phantom study by our clinic. On this ground, the generated DRRs are of good quality and accurately present the structures of pelvic bones essential for treatment localization. In T2-weighted images, bone edges had intensity values very similar to those of adjacent muscles. This may limit their use for automatic segmentation of the bones supporting our choice of the MRI sequence. Our next goal is to develop an automated segmentation method for the pelvic bones in the T1/T2*-weighted MR images. The segmentation is aimed to be implemented into the same step with the segmentation of other tissues needed for the RT treatment planning.

Acknowledgements

Rita Janes, M.D., is greatly acknowledged for her valuable comments and language corrections that significantly improved the manuscript. The authors report no conflicts of interest. The authors alone are responsible for the content and writing of the paper.

Declaration of interest: The authors report no conflicts of interest. The authors alone are responsible for the content and writing of the paper.