Abstract
One of the challenges in PET/MRI is the derivation of an attenuation map to correct the PET image for attenuation. Different methods have been suggested for deriving the attenuation map from an MR image. Because the low signal intensity of cortical bone on images acquired with conventional MRI sequences makes it difficult to detect this tissue type, these methods rely on some sort of anatomic precondition to predict the attenuation map, raising the question of whether these methods will be usable in the clinic when patients may exhibit anatomic abnormalities. Methods: We propose the use of the transverse relaxation rate, derived from images acquired with an ultrashort echo time sequence to classify the voxels into 1 of 3 tissue classes (bone, soft tissue, or air), without making any assumptions on patient anatomy. Each voxel is assigned a linear attenuation coefficient corresponding to its tissue class. A reference CT scan is used to determine the voxel-by-voxel accuracy of the proposed method. The overall accuracy of the MRI-based attenuation correction is evaluated using a method that takes into account the nonlocal effects of attenuation correction. Results: As a proof of concept, the head of a pig was used as a phantom for imaging. The new method yielded a correct tissue classification in 90% of the voxels. Five human brain PET/CT and MRI datasets were also processed, yielding slightly worse voxel-by-voxel performance, compared to a CT-derived attenuation map. The PET datasets were reconstructed using the segmented MRI attenuation map derived with the new method, and the resulting images were compared with segmented CT-based attenuation correction. An average error of around 5% was found in the brain. Conclusion: The feasibility of using the transverse relaxation rate map derived from ultrashort echo time MR images for the estimation of the attenuation map was shown on phantom and clinical brain data. The results indicate that the new method, compared with CT-based attenuation correction, yields clinically acceptable errors. The proposed method does not make any assumptions about patient anatomy and could therefore also be used in cases in which anatomic abnormalities are present.
A new generation of medical imaging scanners, which will combine the high sensitivity of functional imaging by PET with the wide range of applications of MRI scanners, is currently being developed (1–4). Combined PET/MRI offers several advantages over PET/CT, of which the following relate to the particular strength of the MRI method. First, it is possible to acquire anatomic images with better soft-tissue contrast than CT. Second, the radiation dose to the patient is significantly smaller. Third, other applications such as spectroscopy or molecular imaging are also possible with MRI. A PET/MRI scanner will also be a useful research tool, for example, enabling researchers to validate molecular MRI protocols against PET, which can be considered the current gold standard in molecular imaging.
One of the difficulties in the design of a multimodality PET/MRI scanner is the derivation of an attenuation map to correct the PET image for attenuation. The necessity of good attenuation correction has been demonstrated by Huang et al. (5) and, more explicitly, in the case of quantitative PET. PET scanners use a transmission source integrated in the PET gantry to acquire an attenuation map. PET/CT scanners acquire a CT image and then rescale Hounsfield units to 511-keV linear attenuation coefficients (6). In the case of a fully integrated PET/MRI scanner, the small space available inside the bore of the magnet and the presence of a high magnetic field render both options unfeasible.
The attenuation map reflects the tissue-density distribution across the imaging volume. The high density of some biologic tissues—for instance, cortical bone—is mainly caused by atoms with a high atomic number (compared with hydrogen, carbon, oxygen, and other abundant atoms in tissue) such as calcium. This density is the direct source of contrast in a CT image, making the use of CT for attenuation correction straightforward. The MRI signal is governed by the proton density and relaxation mechanisms and has no direct correlation with the tissue density as do CT images. The challenge in MRI-based attenuation correction is finding a way to determine tissue density from a set of MR images, for example, by separating different tissue classes and assigning the corresponding linear attenuation coefficients.
Some methods have been suggested for deriving attenuation maps from MR images. These methods can be divided into 2 main classes: template-based and segmentation-based. The template-based methods start from an MRI template and an attenuation map template (7). The MRI template is registered to a patient MR image using a nonrigid procedure. The same nonrigid transformation is applied to the attenuation map template, resulting in a deformed template attenuation map, which theoretically predicts the attenuation map that would be acquired with a PET transmission scan of the patient.
With segmentation-based methods, the attenuation map is derived directly from the MR image intensity (8,9). Most methods use a 2-step approach: first the image is segmented into tissue classes with a known linear attenuation coefficient (mainly bone, soft tissue, and air), after which the voxels belonging to these tissue classes are assigned the corresponding linear attenuation coefficient. However, the low signal intensity of cortical bone in images acquired with conventional MRI sequences makes it difficult to distinguish bone from air, although there is a strong difference between the attenuation coefficients of bone and air. Therefore, most segmentation-based methods also use anatomic background knowledge to determine whether a voxel contains bone or air.
Recently, a new method has been proposed that uses a combination of pattern recognition and atlas registration to predict the attenuation map from the MR image and a set of aligned MRI and CT datasets (10). First, atlas MRI datasets are coregistered to the MR image of the patient for whom the attenuation map should be derived. The resulting spatial transformations are then applied to the corresponding CT atlas datasets. Subsequently, a pseudo-CT dataset is generated for each voxel as a weighted contribution of intensities in the coregistered CT atlas images. The pseudo-CT is then used for attenuation correction.
A common property of all described methods is that the prediction of the attenuation map is somehow based on anatomic reference data. Because there can be significant variability in patient anatomy, the reliability of these methods could be an issue in routine clinical practice. One example is the frontal sinus in the skull (with significant interpatient variations in size and volume). A method that does not rely on anatomic reference data would be applicable to all patients without anatomic restrictions.
An approach in which the different tissue types can be discriminated solely on the basis of MR image contrast would be preferable over the described methods because it would render the use of anatomic reference data obsolete. The low signal intensity of cortical bone in conventional MRI, and hence the low contrast between air and cortical bone, is caused by the lower water content of this tissue than soft tissue and the fast transverse relaxation rate (R2) (short transverse relaxation time constant [T2]). Because the relaxation of protons in cortical bone occurs too quickly, the MRI signal has disappeared before it is sampled.
Ultrashort echo time (UTE) sequences were specifically designed to visualize tissues with a short T2, such as tendons or ligaments (11). However, cortical bone has an even shorter relaxation time than these tissues. Therefore, it is necessary to investigate the feasibility of visualizing cortical bone with these sequences. We first report the results of an investigation into the proton density and relaxation time of cortical bone. We then propose an image-processing method that enhances contrast between the relevant tissue classes and yields the attenuation map. As a proof of principle, the segmentation of a UTE MR image of a phantom with realistic tissue composition is compared with a segmented CT image of the same phantom. In a final step, the clinical applicability of the proposed method is evaluated on brain PET/CT and MRI datasets of 5 patients.
MATERIALS AND METHODS
MRI Properties of Cortical Bone
The feasibility of visualizing cortical bone with MRI depends on the T2 and proton density of this tissue type. To investigate these properties, relaxometry measurements were performed on samples of cortical bone. Five samples were taken from the shaft of a bovine femur bone acquired from a local slaughterhouse. The bovine femur bone was chosen because it has a thick cortical shell allowing large sample sizes to be obtained.
The samples were analyzed with a 0.5-T relaxometer (Minispec mq20; Bruker Optics). The proton density was quantified as a value relative to the proton density of water. The signal amplitude 170 μs after radiofrequency excitation of different volumes of water was measured. Subsequently, the signal amplitude of the bone samples was measured with the same receiver gain settings. The volume of the samples was derived from a micro-CT scan.
The T2 of each sample was measured with a multiple spin-echo sequence with Carr-Purcell-Meiboom-Gill phase alternation developed in-house (12). With this sequence, 50 consecutive echoes are acquired after excitation with a 90 radiofrequency pulse, yielding 50 data points on the transverse relaxation curve. By fitting an exponential decay function exp(−t/T2) to the acquired data points, the sample averaged T2 is determined. The interecho time spacing was 170 μs.
Acquisition of Phantom UTE MR and CT Images
MR images were acquired on an Achieva 3.0-T system (Philips). The UTE sequence uses radial sampling of k-space, resulting in a spheric reconstructed field of view (13). A hard radiofrequency block pulse is used for excitation, and a low flip angle (10°) makes it possible to keep the length of this excitation pulse much shorter than 100 μs. As the goal of a UTE sequence is to start acquisition as quickly as possible after excitation, the free induction decay (FID) is sampled after excitation rather than first refocusing the proton spins and sampling the resulting gradient or spin echo. For simplicity's sake, we will call the FID the first echo in what follows. The first-echo time TE1 is defined as the interval between the end of the radiofrequency excitation pulse and the start of FID sampling. Depending on the receive coil used, this interval can result in echo times of 70−150 μs—echo times that are much faster than those of conventional MRI sequences (which usually have echo times of 1 ms or higher). After sampling over the sampling time, a gradient is used to refocus the spins and sample a second echo at TE2. The same sampling time is used for the second echo.
A realistic phantom was needed for a proof of principle, resulting in MR and CT images with signal intensities comparable to actual human data. The phantom also had to contain all 3 tissue classes considered for segmentation (bone, soft tissue, and air). The head of a pig was chosen because its tissue composition is comparable to that of a human head.
A combination of Philips Flex coils was used to achieve a field of view large enough to cover the entire phantom. The first echo time was chosen to be as short as possible (0.072 ms for the coils, which was also shorter than any other receive coil) to minimize T2 relaxation before the start of sampling. The sampling time was 1.0 ms to achieve an optimal signal-to-noise ratio for tissues with a short T2 (the T2 of cortical bone is approximately 1.5 ms). The second-echo time was also chosen to optimize R2 fitting (TE2 = 1.747 ms) (14). In addition, the shortest possible repetition time (4 ms) was chosen to reduce total image acquisition time. The images were reconstructed to a 176 × 176 × 176 voxel matrix, with a 1.875-mm isotropic voxel dimension. This resulted in a field of view covering a sphere with a 165-mm radius. Constant-level-appearance coil sensitivity correction was used to manage the highly nonuniform sensitivity of the Flex coils. The total imaging time was approximately 3.5 min.
A CT scan of the pig head phantom was acquired on a Gemini TF PET/CT scanner (Philips) using the following parameters: 120 kVp, 200 mAs, and 1.5-mm slice thickness. The in-plane resolution was 1.17 mm. The total reconstructed matrix was 512 × 512 × 264 voxels.
Acquisition of Clinical PET/CT and MRI Datasets
The applicability of the method on clinical data was tested on 5 clinical PET/CT and MRI datasets. Five epilepsy patients for whom a brain PET/CT and MRI examination was planned also received a UTE MRI scan. The PET/CT was acquired on a Philips Gemini TF PET/CT. A low-dose CT suitable for attenuation correction was acquired with standard settings (50 mAs, 120 kVp). The slice thickness was 5 mm. The in-plane resolution was 1.17 mm. The reconstructed matrix was 512 × 512 × 40 voxels. 18F-FDG PET images were acquired with the standard brain PET protocol. The acquisition time was 7 min. The images were reconstructed to a 128 × 128 × 90 matrix, with 2-mm isotropic voxel dimension. After the PET/CT examination, the patient was transferred to the MRI department, where a standard brain MRI protocol was performed and a brain UTE scan was obtained.
A Philips 8-channel SENSE head coil was used for the acquisition because it has a better uniformity and is more practical in use than the Flex coils. Other echo times were used in this experiment, because the minimum echo time achievable with the head coil was 0.140 ms, which was used as the first-echo time. The second-echo time was 1.8 ms. The same sampling and repetition time were used as before. This yielded a total acquisition time of 6 min. The images were reconstructed to a 224 × 224 × 224 matrix, with 1.3-mm isotropic voxel dimension. The field of view was a sphere with a 145.6-mm radius. The acquisition parameters of the phantom and human datasets are summarized in Table 1.
Image Processing
The image-processing protocol consists of the same 3 steps for both the phantom and the patient data: first the MR, CT, and PET images are coregistered. Then the R2 map of the MR image is calculated, and a mask is applied to correct for voxels containing air. In a final step, the MR and CT images are segmented. The workflow of the image processing is summarized in Figure 1. All images were processed with methods implemented in the Insight Toolkit (15).
The R2 is the inverse of the spin–spin relaxation time T2. High R2 values are expected in cortical bone, and low R2 values are expected in soft tissue. R2 can be estimated from 2 MR images acquired at different echo times but with other parameters kept the same. If the intensity of the MRI signal in a certain voxel at time 0 is I0, the intensity at the 2 echo times will be:Eq. 1Eq. 2where I1 and I2 depict the signal intensity at TE1 and TE2, respectively. The R2 map shows the calculated R2 value of each voxel in the image. This value can be derived easily from Equations 1 and 2:Eq. 3
R2 is suitable for distinguishing cortical bone and soft tissue, because of the significant difference in relaxation rate between both tissue types. The distinction between air and tissue is more difficult because voxels containing air can have a large calculated R2 caused by artifacts (e.g., ringing) or noise in the first-echo image and the better quality of the second-echo image. This distinction is improved by smoothing the images before further processing, but some imaging artifacts cannot be corrected with that method, possibly causing certain voxels containing air to have a nonzero intensity in the first-echo image but a zero or close-to-zero intensity in the second-echo image. In Equation 3, these properties lead to a high R2. It is, therefore, necessary to apply a mask to the data, which aims at setting all voxels containing air to 0 in the R2 map and is done by creating a binary mask from the first-echo image.
The mask is derived in 2 steps. First, a region-growing approach is used to determine the outer contour of the patient (or phantom). Eight seeds are used, placed at the outer corners of the image. Starting from these seeds, all neighboring voxels that have an intensity below a certain value (outer threshold intensity) are set to 0. Next, all voxels that have an intensity below a different threshold (inner threshold intensity) but are not necessarily connected to the region that was grown in the first step are also set to 0. All remaining voxels are set to 1. The mask is derived in 2 steps because the high intensity of the border of the patient, mainly caused by fat in the skin, allows for a higher threshold outside than inside the patient.
A voxel-by-voxel multiplication of the binary mask with the R2 map is used to calculate the corrected R2 map. Because this image now has zero intensity in voxels containing air, medium intensity in voxels containing soft tissue, and high intensity in voxels containing cortical bone, the image is comparable to a CT image.
The final step in the derivation of the attenuation map is the segmentation of the corrected R2 map. This segmentation is done with a simple mapping M of the R2 value to linear attenuation coefficients:
The threshold value of 0.5 was derived from the results of the sample measurements. The attenuation coefficients of bone and soft tissue were obtained by calculating the average attenuation coefficient of bone and soft tissue in the CT-derived attenuation maps of 20 patients in our hospital. The resulting image can be used as the attenuation map.
To quantify the accuracy of the segmentation, the assigned tissue class from the MRI segmentation was compared with the actual (segmented CT) tissue class for each voxel. The CT dataset was segmented using conventional value ranges for bone, soft tissue, and air Hounsfield units, and the voxels were assigned the linear attenuation coefficients corresponding to the segmented tissue class.
The comparison of assigned tissue classes does not give an intuitive representation of the effect of the segmentation errors on final PET image quality, because the voxel-by-voxel evaluation takes into account only local differences between both attenuation maps. In PET reconstruction, the value in every voxel is affected by the attenuation values of all voxels on all lines of response that intersect with this voxel. To determine the effect the segmentation errors would have on the resulting PET images, we have simulated the acquisition and reconstruction process for each clinical dataset.
First, the PET data were forward-projected along all lines of response in the Philips Gemini scanner. The projected data were attenuated by the total attenuation factor for that line of response in the CT-based attenuation map. The reconstruction was then performed twice: once using the CT-based and once using the MRI-based attenuation map. The maximum likelihood expectation maximization algorithm was used, and the reconstruction was stopped after 20 iterations. The average percentage difference between both reconstructed images in a region of interest containing the complete brain was then calculated. The percentage difference was defined for each voxel as the ratio of the difference between both PET images to the value in the CT-based reconstruction.
MRI Properties of Cortical Bone
Figure 2 shows measured signal amplitude versus volume of the cortical bone samples and 5 samples of distilled water. Linear regression was performed on both sets, yielding approximate expressions for the signal intensity I as a function of the volume V:Eq. 4Eq. 5
A unit volume of cortical bone will thus provide approximately 29% () of the signal intensity of a unit volume of water. The signal intensity in voxels of cortical bone will be significantly lower than in voxels of water but still easily detectable.
The T2 values measured in the samples were in the range of 1.37−1.70 ms, with an average over 5 samples of 1.51 ms. This average corresponds to values for R2 between 0.59 and 0.73 ms−1. A suitable threshold value for discriminating soft tissue from cortical bone on the R2 map can now be chosen, for example, 0.5 ms−1.
Phantom CT and MR Images
The MR images acquired with the UTE imaging sequence at both echo times are reconstructed separately. Transverse slices of the first and second-echo image of the pig head MRI dataset are shown in Supplemental Figure 1 (supplemental materials are available online only at http://jnm.snmjournals.org). The images have been smoothed first. In areas of cortical bone (the jawbone is indicated), the signal intensity is medium to high in the first-echo image and much lower in the second-echo image. Soft tissue has a medium to high signal content in both images, whereas air has a low signal intensity in both images. These observations agree with the expectations. The quality of the first-echo image is not good, showing many artifacts even after smoothing. This poor quality will affect the R2 map and derivation of the air mask.
The uncorrected R2 map is shown in Figure 3A. The high R2 value of cortical bone is clear, but cavities containing air also show voxels with a high R2. The corrected R2 map (obtained by multiplying the R2 map with the air mask) (Fig. 3B) shows distinctive intensity values in bone, soft tissue, and air regions.
The same transverse slice of the segmented UTE MR image is shown in Figure 3C, together with the segmented CT image (Fig. 3D). To provide a quantitative impression of the quality of segmentation, a voxel-by-voxel comparison was done between the segmented UTE MR and CT images. A region of interest was chosen inside the head to avoid deterioration of the results by the oversegmentation of skin. A large majority of the voxels is assigned to the correct tissue class: 89% for bone, 91% for soft tissue and 82% for air. The discrimination between tissue (either soft or bone) and air is correct in 97% of the voxels.
Patient PET/CT and MR Images
The images from the 5 clinical datasets were processed in the same manner as the phantom images. Figure 4 shows transverse slices through the eye sockets of 1 patient. This region contains complex bone structures. The uncorrected R2 map in Figure 4A shows significant artifacts. The corrected version, however, shows a much better image (Fig. 4B). The segmented MR (Fig. 4C) and CT (Fig. 4D) images are also shown.
A volume of interest containing the head from the top to the base of the skull was used for evaluation. The voxel-by-voxel comparison of the segmented MR and CT images yielded slightly worse results than the phantom images, with an overall accuracy of between 85% and 95% for all patients. The discrimination between tissue and air was at least 95% in all cases. The PET images reconstructed with CT-based attenuation correction and MRI-based attenuation correction showed average differences of around 5% in the complete brain. Figure 5 shows a transverse (Fig. 5A) and sagittal (Fig. 5B) slice of the percentage difference map of the same patient as shown in Figure 4. The largest deviation in the complete brain of this patient was an underestimation of 34% in the border of the right eye socket. Maximum deviations in other datasets were also between 20% and 40%. The same transverse slice of both reconstructed PET images is shown in Figure 6.
DISCUSSION
The proton density and T2 measurements that were described in this article suggest that the visualization of cortical bone with MRI is possible if the signal is acquired shortly after excitation. This visualization is possible with UTE sequences. We measured T2 values of cortical bone, which are significantly longer than the values cited in other work (usually around 0.5 ms). These differences could be because we used samples of bovine cortical bone, although this explanation seems unlikely. Our measurements are also supported by the R2 values measured in the clinical datasets. One could argue that cortical bone should then also be visible with conventional MRI sequences, which is clearly not the case. This can be explained by the low proton density in addition to the T2 (which is still short). The combination of these 2 properties in a tissue requires a UTE acquisition scheme to yield a usable MRI signal.
The results of our phantom and clinical study demonstrate that UTE-derived R2 maps, corrected with an air mask, are promising for MRI-based attenuation correction. The different image-processing steps generate contrast between air, soft tissue, and bone. The soft tissue or bone discrimination is better using the R2 map rather than using a subtraction of both echo images. The R2 map is quite sensitive to noise in voxels containing air. The distinction between tissue and air, therefore, relies completely on the first-echo image acquired by the UTE sequence.
In the phantom images, it is clear that major bone structures such as the jawbones are correctly segmented. The air cavity in the middle is also detected. Some segmentation errors are clear, mainly an overdetection of bone in the central part of the phantom and in the skin. The latter is possibly in part related to the high collagen content and thickness of porcine skin. Collagen is also present in cortical bone and is in part responsible for the short T2 of cortical bone. Clinical images demonstrated that although the human skin also contains a lot of collagen, it is thinner so it is expected not to cause the same amount of error.
The comparison of PET images reconstructed with CT-based and with MRI-based attenuation correction suggests that the segmentation errors do not cause large errors in the resulting PET image. The percentage difference image shows that most errors in the reconstructed image can be found at the borders of the skull, which is obviously where segmentation is most difficult. One should therefore be careful when looking at the reconstructed PET images in the vicinity of bone. Overestimation or underestimation of bone content could, respectively, cause artificial hot spots or tumors to become invisible in extreme cases. None of the evaluated clinical datasets exhibited effects that large. This is also visible in the reconstructed PET images: the 34% difference in the vicinity of the orbits did not yield artifactual appearances on the MRI attenuation-corrected PET images.
Although the discrimination is adequate in the described cases, it is clear that some improvement should still be made before the method can be used in routine clinical practice. A first consideration is the low quality of the first-echo image. This low quality is in part caused by the acquisition scheme: sampling the FID yields only half the k-space data of a normal gradient echo and thus leads to incomplete sampling, which can be counteracted only by increasing the sampling density. However, this would inevitably mean increasing the imaging time. It has also recently come to our attention that by carefully tuning certain acquisition parameters some of the artifacts can be removed. Unfortunately, the clinical data had already been acquired by the time we received this information.
The 3-dimensional radial acquisition has some advantages (e.g., less sensitive to motion artifacts) but also leads to a spheric field of view. Although this spheric field of view is not a problem in brain imaging, where the sphere will not lead to truncation of the attenuation map if properly positioned, it could lead to truncation in whole-body imaging. The use of 2-dimensional radial or spiral acquisition could provide a partial solution as this would mean that a cylindric field of view could be acquired, thus making avoiding truncation much more straightforward.
A final remark can be made about the lengthy scan times of UTE sequences. In the described study, UTE imaging time was 6 min. This is quite long considering that this would mean that for each bed position 6 min of total MRI study time would be needed for acquiring the attenuation map. However, this imaging time could be significantly reduced by acquiring at a lower resolution. Although lower-resolution acquisitions may be difficult in the brain where many small bone structures and air cavities are present, it is conceivable that a much lower resolution could be used in whole-body imaging. We have not yet investigated this possibility.
Because the presented technique truly visualizes cortical bone, it is possible to use this method to estimate an attenuation map even when the considered patient has nonstandard anatomic features. This is an important advantage over most other methods, which do not explicitly visualize cortical bone but derive its location in an indirect manner, mainly by making some assumption about the general anatomy of the patient. Although the other methods will work fine in patients with an anatomy that can be fitted to a standardized anatomic template, clinical examination is often performed on patients with deviating anatomy. Therefore, it is important to make as few assumptions as possible about the anatomy of the patient when reconstructing medical images.
CONCLUSION
We have demonstrated the feasibility of using the R2 map derived from images acquired with a UTE sequence to estimate the attenuation map from MR images. The performance was evaluated on a realistic phantom CT and MRI dataset and on a patient brain PET/CT and MRI dataset. Between 85% and 95% of the voxels were assigned to the correct tissue class. By reconstructing PET images with both MRI- and CT-based attenuation correction, we showed that the effect of the segmentation errors on the resulting PET images is 5% on average. Maximum deviations varied between 20% and 40% and were mainly found close to bone structures. These deviations did not cause artifactual lesions. Our method relies only on actual patient data and creates contrast between air, soft tissue, and cortical bone based on these data. The applicability in whole-body examinations still has to be investigated.
Acknowledgments
This work was supported by the EU FP7 project HYPERimage (grant agreement 201651). This work was also supported by BELSPO, through the project NIMI.
Footnotes
-
COPYRIGHT © 2010 by the Society of Nuclear Medicine, Inc.
References
- Received for publication May 9, 2009.
- Accepted for publication February 12, 2010.