MRI/CT - CiteSeerX

6 downloads 0 Views 1MB Size Report
Mar 14, 2008 - at the Robert Steiner MR Unit, Hammersmith Hospital, London. Pixel-by-pixel ..... Hoffman EJ, Cutler PD, Digby WM, Mazziotta JC. 3-D phantom ...
Journal of Nuclear Medicine, published on March 14, 2008 as doi:10.2967/jnumed.107.041871

PET Image Denoising Using a Synergistic Multiresolution Analysis of Structural (MRI/CT) and Functional Datasets Federico E. Turkheimer1,2, Nicolas Boussion3, Alexander N. Anderson2, Nicola Pavese2, Paola Piccini1,2, and Dimitris Visvikis3 1Department

of Clinical Neuroscience, Division of Neuroscience and Mental Health, Imperial College, London, United Kingdom; Clinical Sciences Centre, Hammersmith Hospital, London, United Kingdom; and 3INSERM, U650, Laboratoire du Traitement de l’Information Medicale (LaTIM), CHU Morvan, Brest, France

2MRC

PET allows the imaging of functional properties of the living tissue, whereas other modalities (CT, MRI) provide structural information at significantly higher resolution and better image quality. Constraints for injected radioactivity, technologic limitations of current instrumentation, and inherent spatial uncertainties on the decaying process affect the quality of PET images. In this article we illustrate how structural information of matched anatomic images can be used in a multiresolution model to enhance the signal-to-noise ratio of PET images. The model states a flexible relation between function and structure in the brain and replaces high-resolution information of PET images with appropriately scaled MRI or CT local detail. The method can be naturally extended to other functional imaging modalities (SPECT, functional MRI). Methods: The methodology is based on the multiresolution property of the wavelet transform (WT). First, the coregistered structural image (MRI/CT) is downgraded to the resolution of the PET volume through appropriate filtering. Second, a redundant version of the WT is applied to both volumes. Third, a linear model is applied to the set of local coefficients of both image volumes and resulting parameters are recorded. The overall set of linear coefficients is then modeled as a mixture of multivariate gaussian distributions and fitted through a k-means algorithm. Finally, the local wavelet coefficients of the PET image are substituted by the corresponding values of the MRI/CT set calibrated according to the resulting clustering. The methodology was validated on digital simulated images and clinical data to evaluate its quantitative potential for individual as well as group analysis. Results: Application to real and simulated datasets shows very effective noise reduction (15% SD) while resolution is preserved. Conclusion: The methodology is robust to errors in the coregistration parameters, practical to implement, and computationally fast. Key Words: PET; sensitivity; multimodality; PET/CT; PET/MRI; fusion; synergistic; denoising; wavelets J Nucl Med 2008; 49:657–666 DOI: 10.2967/jnumed.107.041871

Received Mar. 20, 2007; revision accepted Dec. 20, 2007. For correspondence contact: Federico E. Turkheimer, PhD, Department of Clinical Neuroscience, Division of Neuroscience and Mental Health, Cyclotron Bldg., Room 236, Hammersmith Hospital, DuCane Rd., London W12 0NN, UK. E-mail: [email protected] COPYRIGHT ª 2008 by the Society of Nuclear Medicine, Inc.

T

he application of PET to image radiotracer distribution in humans is limited largely by scanner sensitivity (1). Although in the recent past counting statistics have been steadily increasing with the introduction of new technologies such as 3-dimensional PET (2) and new scintillator materials such as lutetium oxyorthosilicate (3), increased sensitivity has been outstripped by improvements in the spatial resolution of the detectors. In the near future, further improvements in count statistics may be brought by improvements in scanner technology (1). However, poor signal-to-noise ratio (SNR) is likely to distinguish PET from CT and MRI in the years to come. Computational approaches may play an important role in PET noise-reduction as long as the increase in the SNR is not associated with a resolution loss and their implementation is practical and fast. In particular, PET denoising approaches based on the wavelet transform (WT) (4–6) have consistently been reported to increase accuracy and precision of PET images in a wide variety of contexts (7). The WT is a recently introduced mathematic tool for the treatment of signals with nonstationary behavior (e.g., a hammer blow, a plane flyover noise, etc.) (8). The counterpart of the WT is the Fourier transform that achieves optimal encoding of periodic signals. In PET, the WT is applied to the spatial distribution of the radiotracer after reconstruction, and it converts the original pixel values in wavelet coefficients that represent signal intensity at different locations at different resolutions. The resulting multiscale representation has interesting properties (5), among which the one relevant to the present work is the ability to concentrate signal in a few coefficients at different scales while the noise remains homogeneously distributed. This separation between signal and noise may be exploited to increase image SNR if one is able to remove the noise in the wavelet space before transforming back the data into images. Alpert et al. (7) have shown through accurate simulations of dynamic PET sequences that in the case of neuroreceptor studies one could ideally achieve a 2-fold increase in the SNR of the 4-dimensional

PET DENOISING

WITH

MRI/CT WAVELETS • Turkheimer et al.

jnm041871-tp n 3/12/08

657

dataset and 1.5 SNR improvements in the derived parametric map. However, ideal performance implies the knowledge of the distribution of the ‘‘true’’ wavelet coefficients as opposed to the distribution of those due to noise. Recently, we have shown that, through a multiscale description of the physiology of the organ imaged by PET, a hierarchical filter can be derived that is able to produce SNR improvements in parametric images close to the predicted optimal mentioned above (9). In this work we further explore the possibility of using information on the object imaged by the PET procedure to enhance image quality without resolution loss. In particular, we make use of a coregistered structural/anatomic image and postulate a stochastic model between the wavelet distribution of the structural image (MRI or CT) and the wavelet distribution of the functional image (PET). This work builds on previously published results on the use of coregistered MRI or CT datasets for the partialvolume correction of PET images (10). By substituting the high-resolution wavelet coefficients of the PET volume with those of the MRI volume, Boussion et al. were able to increase the resolution of the PET volume, whereas the appropriate scaling between PET and MRI/CT coefficients at the finest resolution was obtained by assuming a global linear relationship between PET and MRI/CT wavelets at lower resolutions (10). In the procedure introduced here, we adopt the inverse approach where, after degradation of the matched MRI/CT volume to the resolution of the PET, we substitute the PET wavelets with the corresponding MRI/CT wavelets after appropriate scaling. Differently from before, however, we relax the assumption of a simple scalar relation between structure and function and postulate a stochastic model to relate structural and functional information. The performance of the resulting procedure was evaluated through simulations using a digital phantom and in the clinical setting, which obviously entails the presence in real terms of factors such as coregistration errors and unknown relations between structure and function. Clinically, the SNR improvement brought by the new algorithm was assessed in whole-body PET/CT and as improvement in detection of microglial activity in a group of Huntington’s disease (HD) patients imaged with PET and [11C]-(R)PK11195.

MATERIALS AND METHODS The computational procedure consists of 3 steps. First, the coregistered MRI/CT volume is degraded to the resolution of the PET scanner. In all instances we applied to the MRI/CT volume an isotropic gaussian kernel so that final image resolution matched the one of the PET scanner (simulated or real). Second, the WT is applied to both PET and MRI/CT volumes. Third, the wavelet coefficients of the PET volume are replaced by those of the MRI/CT volume after appropriate scaling. The scaling between MRI/CT and PET coefficients is obtained by application

658

THE JOURNAL

OF

of an appropriate stochastic model. Detailed implementation is as follows:

WT Regression Model The application of the WT to tomographic images has been described in detail previously (5) and will be summarized here. The dyadic wavelet transform (DWT) consists of the iterative application of 2 conjugate filters, the high-pass filter (H) and the low-pass filter (L) and subsequent decimation. At each iteration, the output of filter H is decimated and stored as the wavelet coefficients for that resolution, whereas the decimated output of filter L, the residuals, is passed to the following stage for subsequent filtering and decimation. The inversion of the procedure returns the original data-vector. The 2-dimensional DWT is obtained by separate application of the DWT to rows and columns. This operation generates 4 quadrants. The HH quadrant (diagonal details) collects the decimated output of the H filter applied to both rows and columns of the data matrix. The HL quadrant (horizontal details) contains the coefficients resulting from the application of the H filter to the rows and the L filter to the columns. The LH quadrant (vertical details) is obtained when the L filter is applied to the matrix rows and the H filter to the matrix columns. Finally, the LL quadrant contains the residuals obtained from the L filter operated on both rows and columns. The coefficients in the quadrants HH, HL, and LH represent the 2-dimensional DWT for the particular resolution. Formally, let I(s) be the image matrix that, for ease of comprehension, we assume to be square N·N (note that the application of the DWT requires dimension N to be a power of 2) and is indexed by s:s 5 (x,y). We define the DWT operator W and W21 its inverse. Application of the DWT to the image produces: IðwÞ 5 WðIðsÞÞ 5 Ui ½HHi ðwÞ HLi ðwÞ LHi ðwÞ:

i 5 1; 2; . . . ; D Eq. 1

The index w:w 5 (x,y) spans the dimension N/2i·N/2i, i is an index of the resolution levels, U is the union operator, and D is the top wavelet resolution level such that 2D 5 N. The 3 coefficients for quadrants HHi(w), HLi(w), and LHi(w) correspondent to pixel s can be recovered from the correspondence: w 5 bs=2i c;

Eq. 2

where bAc is the floor operator that rounds the elements of A to the nearest integer less than or equal to A. In the context of tomographic images however, the traditional DWT is usually discarded in favor of redundant WT approaches (5,10). In this work we used the cycle-spinning DWT (WTCS) that applies iteratively the WT to a series of shifts, column wise and row wise, of the original image (5). This procedure generates a WT representation that preserves the original signal content of the image to the expense of a bigger storage matrix. A detailed description of the WTCS is contained in (5); here the procedure is only sketched. The first-resolution level is obtained by application of the DWT to the original image and to the 3 images obtained, respectively, by a single pixel shift applied to the rows, then to the columns, then to both. The second resolution level is generated by the DWT applied to the 4 residual quadrants LL obtained from the previous step and to

NUCLEAR MEDICINE • Vol. 49 • No. 4 • April 2008

jnm041871-tp n 3/12/08

their shifts (3 each for a total of 16 transforms), etc. If we define the WTCS operator as Wcs and Wcs21 is its inverse then: IðwÞ 5 Wcs ðIðsÞÞ 5 Uij ½HHij ðwÞ HLij ðwÞ LHij ðwÞ: i 5 1; 2; . . . ; D j 5 1; 2; . . . ; 4i :

Eq. 3

Consider now a functional image F(s) (PET, SPECT) and the correspondent structural image S(s) (CT, MRI). Application of the WTCS generates: FðwÞ 5 Wcs ðFðsÞÞ;

Eq. 4

Equation 6b is still valid and linear least squares can still be applied as long as an appropriate variance transformation (e.g., logarithm) is applied to the functional data. At this stage we make the second assumption that the relationship between function and structure varies throughout the image but that there may be homogeneities in this relation that can be exploited to improve the SNR of the functional image. Under this assumption, the most general statistical model P(ai,bi) for the set of coefficients ai(w) and bi(w) is a mixture of multivariate normal distributions, that is: Pðai ; bi Þ 5 Sk Nk ðmk ; s2k Þ=K:

and SðwÞ 5 Wcs ðSðsÞÞ:

Eq. 5

The procedure developed here models the relation between wavelet coefficients of a functional volume F(w) and those of structural volume S(w) one resolution at a time. Therefore, what follows must be intended to be applied to F(w) and S(w) for a fixed index i. The procedure should be applied only at the finest resolution levels (the first 3 in this work), as noise is not an issue at lower resolutions. For a resolution i, we postulate for every pixel s a linear relationship between the correspondent wavelets in the functional image and those in the structural image—in other words:

F2;2N-2k

LHij ðwÞS 1bi ðwÞ1eðwÞ: i

w 5 bs=2 c:

Eq. 6a

The linear relation between the wavelet of the functional image and the one of the structural image is the first order approximation used by Boussion et al. (10); in their work this relation was established globally on the image whereas here it is applied locally. In Equation 6a the linear vectors [. . .]F and [. . .]S collect the wavelet coefficients of the different subbands of resolution ‘‘i’’ for image F and S corresponding to the pixel s, while ai(w) and bi(w) are the coefficients of the linear model and e(w) is a homogeneous zeromean gaussian process. Representation in Equation 6a assumes that, given the higher noise of the functional image, vector [. . .]S can be considered as noiseless. The model in Equation 6a can be resolved by application of linear least squares to:

Eq. 7

In Equation 7, each bivariate normal probability function Nk() is parameterized by its mean vector mk 5 [ak;b k ] and relative variance vector sk2. The model in Equation 7 can be estimated through a clustering K-means algorithm, where the number of clusters K can be calculated by progressive application of the algorithm to increasing values of k and comparative testing of the residuals of model k toward the more complex model k11 through the pseudo F test for nested models (12). The pseudo F can be constructed as:

½HHij ðwÞ HLij ðwÞ LHij ðwÞF 5 ai ðwÞ  ½HHij ðwÞ HLij ðwÞ j 5 1; 2; . . . ; 4I

k 5 1; 2; . . . K:

ðRSSk 2 RSSk 2 1 Þ RSSk 5 : N2k k11  21 N2k21 k

Eq. 8

In Equation 8, RSSi is the residual sum of squares of model i and N is the total number of (a,b) pairs. The threshold for the F test should be selected suitably small; in this instance, the P value 0.01 was adopted. Once the model in Equation 7 is estimated, the wavelet coefficients at Equation 6 can be estimated by incorporating the corresponding estimated parameters mk so that: ½HHij ðwÞ HLij ðwÞ LHij ðwÞFE 5 ak  ½HHij ðwÞ HLij ðwÞ LHij ðwÞS k: 1b Eq. 9 Let FE(w) be the functional volume obtained by application of the regression model for several wavelet resolution levels. The denoised function image can be obtained by application of the inverse DWTCS as:

½HHij ðwÞ HLij ðwÞ LHij ðwÞF 5 ai ðwÞ  ½HHij ðwÞ HLij ðwÞ

FE ðsÞ 5 Wcs21 ðFE ðwÞÞ:

S

Eq. 10

LHij ðwÞ 1bi ðwÞ: j 5 1; 2; . . . ; 4I w 5 bs=2i c:

Eq. 6b

In Equation 6b ai(w) and bi(w) are sample realizations of ai(w) and bi(w). Because the model in Equation 6b is local, the assumption on the gaussian nature and homogeneity of e(w) can accommodate both stationary and nonstationary noise conditions that depend on the particular reconstruction algorithm (11). With iteratively reconstructed images, with a large number of iterations, the nature of the noise may become non-gaussian and maintain non-gaussian properties also in the wavelet domain. If the latter is the case, the model in

Validation I: Synthetic Data Given the specified model in Equation 7, the performance of the algorithm will vary according to the image patterns and to the structural/functional relation. To assess a possible worse-case scenario, we based the simulation study on a realistic but complex pattern, the brain, with variable and mismatched structural/functional relations. We considered a digital 128 · 128 · 10 brain Hoffman phantom (13) that was used as the structural image (Fig. 1A). The functional ½Fig: 1 image was obtained by adding to the phantom different patterns of varying intensities and sizes (Fig. 1B). To achieve realistic condi-

PET DENOISING

WITH

MRI/CT WAVELETS • Turkheimer et al.

jnm041871-tp n 3/12/08

659

FIGURE 1. Digital structural (A) and functional (B) volumes used for synthetic data simulations. Structural volume is a Hoffman brain phantom (128 · 128 · 10). Functional image (ground truth) (C) was obtained by adding local changes in intensities of various sizes not necessarily matching the underlying structural patterns (arrows) and then smoothed with a gaussian filter with FWHM 5 5 mm to simulate a typical neurologic image.

tions, we created local differences in intensities and mismatches between the structural and functional image (note brighter areas in cerebellum as indicated by arrows in Fig. 1B). The image was then smoothed with a gaussian filter of 5-mm full width at half maximum (FWHM). Gaussian noise of varying SD and with the same smoothness (SD) (5%, 7.5%, 10%, 15%, 20%) was then added, and 100 realizations for each noise level were created and processed with the wavelet-based denoising algorithm. A second set of simulations was also devised with the same settings but with a proportional noise model, where SD was not uniform across the image but proportional to the underlying signal. The standard SURE wavelet filter, purposely devised for PET (5) and considered as a state-of-the-art denoising filter, was used for comparison purposes. Performance measures were bias (difference between true and estimated pixel value) and SD of the estimates (at the pixel level) averaged either across the whole image or the regions of mismatch between the simulated functional and structural images. Validation II: Phantom Data from GATE Simulator Although gaussian stationary and nonstationary noise conditions may account for a fraction of PET studies (filtered backprojection [FBP] reconstruction and iterative reconstruction with precorrections), we explored further the performance of the denoising

660

THE JOURNAL

OF

algorithm for emission images using accurate simulation of a cylindric phantom. Images consisted of a simplified numeric version of the IEC phantom 61675-1. The latter consists of a 20-cm diameter by 20-cm long cylinder, containing 6 spheres of different diameters (37, 28, 22, 17, 13, and 10 mm). The numeric version was produced as a set of 64 contiguous planes of 64 · 64 square pixels of 4 · 4 mm in size, and an acquisition with the Philips Allegro PET scanner was subsequently simulated using GATE, a Monte Carlo– based simulator (14). A total of 60 million coincidences were simulated considering a sphere/cylinder activity concentration ratio of 5:1. Images were reconstructed using the OPLEM algorithm (11 iterations) (15). Validation III: Clinical Whole-Body [18F]FDG PET/CT Images The proposed denoising technique was first tested on clinical whole-body PET/CT studies to assess the portability of the algorithm considering images of different characteristics. Whole-body PET/CT images of 4 patients acquired using the Discovery LS (GE Healthcare) were analyzed. Patients were scanned at 55–60 min after injection of an average 365 MBq of [18F]FDG. The whole-body PET acquisition protocol comprised 2-dimensional emission scans of 3 min per bed position, whereas the CT scans (tube acquisition

NUCLEAR MEDICINE • Vol. 49 • No. 4 • April 2008

jnm041871-tp n 3/12/08

parameters of 140 keV and 80 mA) were acquired under normal breathing conditions. PET images were reconstructed with CTbased attenuation correction (16), using the ordered-subsets expectation maximization (OSEM) algorithm (128 · 128 · 205 matrices; voxel size of 4.35 · 4.35 · 4.25 mm3). The number of iterations and subsets used for the reconstruction have been previously optimized for an improved SNR (17). Images were analyzed quantitatively by calculating intensity inside regions of interest (ROIs) manually drawn over either identified lesions or large homogeneous regions (lungs, liver, and other soft tissues). ROIs were thus separated into 2 different groups. In the first one containing large ROIs, noise (defined as the SD of the voxel intensity) was measured before and after denoising. In the second set containing ROIs surrounding small lesions, the SNR was computed as the ratio between mean value and the SD of the voxel intensity inside the specific ROIs and SNR was calculated before and after denoising. Only uniform ROIs around lesions were considered, where uniformity was defined as SD , 10% relative to the mean intensity. Validation IV: [11C]-(R)-PK11195 Group Study In this second part of the validation, the technique was evaluated in a clinical setting that not only contained realistic patterns, varying anatomic/functional correlates, and coregistration errors but also was characterized by a controlled context where patient diagnosis was 100% accurate and the disease pattern was known. For this study we selected a set of 22 [11C]-(R)-PK11195 scans, 12 of which were performed on HD patients with overt symptoms and 10 were age-matched controls (18). [11C]-(R)-PK11195 is a radiotracer selective for the peripheral benzodiazepine receptor that is the most consistently upregulated protein in reactive microglia (19). Microglia, the resident immune cells of the brain, react to any insult in the central nervous system, from stroke to neurodegeneration to neoplasia; this makes [11C]-(R)-PK11195 a widely used general indicator of brain disease (20). In the case of HD, microglia have a pattern of activity that has been well characterized and has a defined anatomic localization in the globus pallidus as well as caudate and putamen (21). Further activity has also been detected in the cortex but less colocalized within defined anatomic structures (21). This bivalent relation between functional activation and anatomic localization represents an ideal test-bed for the methodology introduced in this work. Diagnosis of the disease cohort was 100% accurate as all HD patients had genetically proven disease with an expanded CAG repeat in the IT15 gene on chromosome 4. All [11C]-(R)-PK11195 studies were performed on an ECAT EXACT HR11 (Siemens 966; Siemens Medical Solutions, Inc.) camera (23.4-cm axial field of view, 95 transaxial planes, spatial resolution of 4.8-mm FWHM [transaxial] and 5.6-mm FWHM [axial]). A transmission scan was acquired before the emission scan using a single rotating photon point source of 150 MBq of 137Cs for subsequent attenuation correction and scatter correction. Images were then reconstructed with the reprojection algorithm with the ramp and Colsher filters set to Nyquist frequency (22). The ramp filter offers the highest resolution at the expense of higher noise levels and is therefore the best case for the application of the methodology. [11C]-(R)-PK11195 was manufactured by Hammersmith Imanet Ltd. A mean dose of 294 MBq was injected as a bolus and 3-dimensional sinograms of emission data were then acquired over 60 min as 18 time frames. Volumetric T1-weighted MR images were obtained on a 1.0-T Picker HPQ scanner (Picker International) at the Robert Steiner MR Unit, Hammersmith Hospital, London.

Pixel-by-pixel maps of binding potential (BP) were generated using the simplified reference tissue model (SRTM) (23). Extraction of the reference region was performed using a supervised clustering algorithm that selects gray matter pixels devoid of specific [11C]-(R)-PK11195 binding (24). Parametric images were then coregistered to their respective T1-weighted MR image using as reference the relative sum of the activity images obtained as the sum of all the frames weighted by the estimated number of true events in each frame. At this stage the BP images were denoised using the proposed synergistic fusion technique. Both original and denoised images were subsequently spatially normalized into standard Montreal Neurologic Institute (MNI) space. Statistical parametric mapping SPM2 (Functional Imaging Laboratory, Wellcome Department of Imaging Neuroscience, University College London, London) was used for the image transformations. Normalized parametric images in HD patients and healthy subjects were then contrasted to create mean increase maps of microglial activation in HD using SPM2. Pixel-by-pixel application of the 1-tailed Student t test generated maps of z scores; in the generation of these maps all thresholds (intensity threshold, F statistic threshold, and Student t test threshold) were removed. The average of the z scores in ROIs comprising caudate, putamen, and globus pallidus was used to score the sensitivity of the algorithm under evaluation toward the original parametric maps. The same contrast for both denoised and raw parametric maps was also visualized as z-score maps to allow qualitative evaluation of the results. RESULTS

The synergistic methodology described in the previous sections was implemented using Matlab (The Mathworks Inc.) on a Sun Ultra10 workstation. Computationally, the technique required ,10-min processing time to denoise a 128 · 128 · 100 emission volume. Synthetic Data

Figure 2 allows the qualitative appreciation of the appli- ½Fig: 2 cation of the procedure on a single noisy realization obtained with 10% constant noise. Note that the synergistic filter image is less noisy but still sharp, suggesting no resolution loss, which, instead, is apparent in the SURE filtered image. Results of the simulation studies are summarized in Table 1 in ½Table 1 terms of bias and noise reduction. The 2 rows indicate the average results for the whole phantom (‘‘Global’’) and for those specific areas of the phantom where mismatch was introduced between the structural and functional simulated images (‘‘Mismatch’’). The additional line illustrates the performance of the SURE filter globally (the SURE does not use structural information). There was a remarkable consistency of the results among noise levels and noise models indicating a consistent reduction of the variability that ranged from ;12% on the entire pattern with a minimum of 9.8% on those cerebellar areas where we introduced a mismatch between the structural and functional patterns. Bias on average was quite low reaching a maximum of 1% throughout the 5 simulation studies. This was primarily due to signal loss at the edges with a maximum overall value of 4% indicating small resolution loss; the local regression model could not be effective in the high-resolution frequencies

PET DENOISING

WITH

MRI/CT WAVELETS • Turkheimer et al.

jnm041871-tp n 3/12/08

661

FIGURE 2. Results of simulation study (in this case, with stationary 10% noise level). One realization (10 planes) and output of synergistic and of SURE wavelet filters are shown.

(edges) but still preserved the lower frequency information. The SURE performance clearly illustrates the properties of standard wavelet filters: the significant noise reduction (up to 23% with the constant noise model) was counterbalanced by an average bias of 8% that was due to resolution loss, which around the edges caused signal decreases of .35% (data not shown). Results with the proportional noise model were almost identical in terms of bias and noise reduction and are not reported.

GATE Simulated Phantom

Data and results for the phantom simulation are illustrated in Figure 3, which shows also the line profiles through the 6 ½Fig: 3 spheres for the emission image before and after denoising. As it is apparent from the profiles, application of the structural denoising caused no appreciable loss of resolution. Noise reduction was measured as the percentage reduction of the SD in a square region, 20 · 20 pixels, placed in the center of each of the 64 planes of the cylinder. The resulting reduction

TABLE 1 Constant Noise Model Mean absolute bias of synergistic denoising approach for simulation studies expressed as fraction of signal Noise levels Characteristic 5% 8% 10% 15% Global 0.003 0.002 0.000 0.003 Mismatch 0.010 0.009 0.008 0.007 SURE* 0.076 0.070 0.066 0.055 Mean noise change of synergistic denoising approach compared with raw images for simulation studies

20% 0.005 0.006 0.046

Noise levels Characteristic Global Mismatch SURE*

5% 20.123 20.097 20.227

8% 20.124 20.098 20.231

10% 20.124 20.098 20.230

*Results for the wavelet SURE filter.

662

THE JOURNAL

OF

NUCLEAR MEDICINE • Vol. 49 • No. 4 • April 2008

jnm041871-tp n 3/12/08

15% 20.124 20.099 20.230

20% 20.124 20.099 20.230

FIGURE 3. Data and results for simulated phantom (1 plane only). (A–C) Original digital image (A) and GATE simulated image before (B) and after (C) denoising. (D) Composite of profiles (black lines on images) through hot spheres for raw (dotted) and filtered (continuous) image.

in SD due to denoising was 31% (maximum, 39.9%; minimum, 19.5%).

Clinical Whole-Body [18F]FDG PET/CT Images

Preliminary analysis of the relation between functional and structural wavelets was performed by visual inspection of the coefficients when plotted, ones against the other ones, and confirmed a very good linear distribution. Variance inhomogeneity was also verified using logarithmically transformed functional data but was discarded as unsuitable for this application.

Figure 4 presents 2 examples of images before and after ½Fig: 4 denoising by synergistic PET/CT fusion. Despite a slight smoothing effect across the whole image, the qualitative appearance of images is not degraded by the operation. Furthermore, the denoising is having the same effect across tissues with high and low uptakes, contributing this way to keeping the global homogeneity of processed images. With regard to the quantitative analysis, absolute intensity values were on average reduced by ,5% in the different ROIs considered, demonstrating a small loss of resolution as already seen in the simulated datasets. The noise reduction in the set of large ROIs inside homogeneous tissues, such as lungs or liver, is illustrated in Figure 5A. Fifteen ROIs were ½Fig: 5 drawn for each patient leading to a total of 60 different regions of mean size 5,783 6 976 mm2. Denoising was significant, as noise reduction reached 21% 6 8.2% with minimal and maximal values of 6.2% and 36%, respectively. SNRs before and after denoising in ROIs related to small lesions are given in Figure 5B (mean ROI size, 199 6 90 mm2, 12 lesions in total). The corresponding SNR increase in terms of percentage is also presented, with a mean value of 45.7% 6 22.5% (minimum, 17.1%; maximum, 83.4%). [11C]-(R)-PK11195 Group Study

Quantitatively, the noise reduction introduced by the methodology was measured in the group study as the average of the z scores on the striatum—that is, the area with known pathologic insult and related intense microglial activity. The z-score average was 1.160 when the contrast was generated using raw BP images and 1.355 when the de-noised images were used. Considering that the signal intensity remains practically constant after processing, this result corresponds to a reduction in noise of 15.5%, which is right above the lower bound predicted by the simulations. Figure 6 illustrates the overall result of the comparison ½Fig: 6 between the 2 groups of scans, HD subjects versus controls, in terms of z-score maps for the raw BP maps (Fig. 6A) and those that were denoised (Fig. 6B). Localization of activation pattern is very similar, if not identical, for both raw and

FIGURE 4. Application of waveletsynergistic methodology to [18F]FDG whole-body PET images in oncology. Anatomic data are provided by CT images acquired on a dedicated PET/CT scanner. Coronal slices of 2 patients are presented (from top to bottom): PET before denoising (left), CT (middle), and PET after denoising (right).

PET DENOISING

WITH

MRI/CT WAVELETS • Turkheimer et al.

jnm041871-tp n 3/12/08

663

FIGURE 5. Results of application of the methodology to 4 [18F]FDG whole-body images. (A) Percentage of noise decrease in large and homogeneous areas. Fifteen ROIs per patient (60 ROIs in total) were drawn in lungs and liver. (B) SNR in 12 lesions before and after wavelet denoising. Percentage of SNR increase is presented as well (right-hand side, y-axis) to help in appreciating the level of improvement achieved.

denoised maps, but the latter display an evident increase in sensitivity as the signal is more intense and better delineated without obvious evidence of resolution loss. DISCUSSION

The synergistic use of imaging modalities is meant to overcome the intrinsic limitations and to enhance the specific advantages of the different approaches to allow both anatomic and functional correlation for visualization, increased precision, and accuracy for quantification. In the past, fusion between functional and structural images has been made practical—particularly for brain applications—by the availability of software packages that are able to coregister with good precision 2 modalities acquired independently (25). It has recently become more convenient for different imaging applications with the introduction of PET/CT scanners (26) and the upcoming development of PET/MRI tomographs (27). So far, the quantitative application of image fusion has focused on the use of structural information for partialvolume correction (28). The methodology presented here, which stems from a previous approach to partial-volume correction (10), is intended to use fused structural information to denoise emission tomography images. The technique is of general use and we have shown its applicability to both PET/CT and PET/MRI coregistered volumes. As a totally automatic and fast postprocessing step it also represents a

664

THE JOURNAL

OF

very practical approach. The resulting performance was consistent across several simulated and clinical datasets and across modalities. The study on the synthetic data and the digital phantom predicted a noise reduction between 12.5% and 31%, on average, depending on the complexity of the underlying pattern, without significant loss of resolution. In real terms, this amounts to a 25%–45% reduction in sample size for a group study or to a similar reduction in the injected dose. In the clinical brain study involving the use of [11C]-(R)-PK11195 and PET on control subjects and HD patients, the resulting noise reduction amounted to ;15%. In the PET/CT application of thorax FDG studies, noise reduction was 20%, on average, with a significant increase in SNR (45%). However, on average, absolute intensity changes were similar in magnitude to those seen in the simulated datasets and limited to ,5% irrespective of the activity level in the ROIs considered. The variability in noise suppression was dependent on the varying number of clusters detected by the procedure (data not shown) and, therefore, by the homogeneity in the relation between function and structure in the image (the greater the homogeneity, the smaller the number of the clusters, the higher the noise reduction). The ability of this filter to reduce noise without resolution loss stems from the use of additional information in the form of the structural image. This feature sets this filter apart from other postreconstruction filters as well as from other wavelet filters previously proposed for emission images (5). With such filters as the SURE, which was used as a reference in the

NUCLEAR MEDICINE • Vol. 49 • No. 4 • April 2008

jnm041871-tp n 3/12/08

RGB

FIGURE 6. Results in form of z-score maps for group comparison between an HD group and a control group in study of neuroinflammation measured with PET [11C]-(R)-PK11195. No smoothing was applied. (A) Result when raw parametric maps were used. (B) Outcome when denoised maps were used. Intense signal in basal ganglia and diffuse microglial activation in cortical gray matter are pathologic hallmarks of the disease.

simulation study of this work, no extra information besides the emission image itself is available and the achieved variance reduction is counterbalanced with a resolution loss that is generally not acceptable for clinical use. The synergistic filter is also more generally applicable than the one in (9) which applies to dynamic acquisitions for irreversible tracers only. CONCLUSION

ACKNOWLEDGMENTS

This study was funded in part by the EC-FP6-project DiMI, LSHB-CT-2005-512146A, and the Brittany Region Research Program ‘‘QUANTEP.’’ This project was funded by a grant from Hammersmith Imanet (GE-Healthcare). Hammersmith Imanet Ltd. manufactured 11C-(R)-PK11195 and maintains all PET scanning equipment. REFERENCES

Wavelet models are effective computational tools for the synergistic fusion of images of different modalities. We have demonstrated that, in the context of PET/CT and MRI/PET, these models can deliver a substantial increase in sensitivity of emission images with minimal resolution loss.

1. Cherry SR. Of mice and men (and positrons): Advances in PET imaging technology. J Nucl Med. 2006;47:1735–1745. 2. Cherry SR, Dahlbom M, Hoffman EJ. 3D PET using a conventional multislice tomograph without septa. J Comput Assist Tomogr. 1991;15:655–668. 3. Melcher CL, Schweitzer JS. Cerium-doped lutetium oxyorthosilicate: a fast, efficient new scintillator. IEEE Trans Nucl Sci. 1992;39:502–505.

PET DENOISING

WITH

MRI/CT WAVELETS • Turkheimer et al.

jnm041871-tp n 3/12/08

665

4. Turkheimer FE, Banati RB, Visvikis D, Aston JA, Gunn RN, Cunningham VJ. Modeling dynamic PET-SPECT studies in the wavelet domain. J Cereb Blood Flow Metab. 2000;20:879–893. 5. Turkheimer FE, Brett M, Visvikis D, Cunningham VJ. Multiresolution analysis of emission tomography images in the wavelet domain. J Cereb Blood Flow Metab. 1999;19:1189–1208. 6. Turkheimer FE, Brett M, Aston JA, et al. Statistical modeling of positron emission tomography images in wavelet space. J Cereb Blood Flow Metab. 2000;20:1610–1618. 7. Alpert NM, Reilhac A, Chio TC, Selesnick I. Optimization of dynamic measurement of receptor kinetics by wavelet denoising. Neuroimage. 2006;30:444–451. 8. Mallat SG. A theory of multiresolution signal decomposition: the wavelet representation. IEEE Trans Pattern Anal Mach Intell. 1989;11:673–693. 9. Turkheimer FE, Aston JA, Asselin MC, Hinz R. Multi-resolution Bayesian regression in PET dynamic studies using wavelets. Neuroimage. 2006;32:111–121. 10. Boussion N, Hatt M, Lamare F, et al. A multiresolution image based approach for correction of partial volume effects in emission tomography. Phys Med Biol. 2006;51:1857–1876. 11. Wilson DW, Tsui BMW. Noise properties of filtered-backprojection and ML-EM reconstructed emission tomography images. IEEE Trans Nucl Sci. 1993;40: 1198–1203. 12. Sparks DN. Euclidean cluster analysis. J R Stat Soc Ser C Appl Stat. 1973; 22:126–130. 13. Hoffman EJ, Cutler PD, Digby WM, Mazziotta JC. 3-D phantom to simulate cerebral blood flow and metabolic images for PET. IEEE Trans Nucl Sci. 1990; 37:616–620. 14. Lamare F, Turzo A, Bizais Y, Le Rest CC, Visvikis D. Validation of a Monte Carlo simulation of the Philips Allegro/GEMINI PET systems using GATE. Phys Med Biol. 2006;51:943–962. 15. Reader AJ, Ally S, Bakatselos F, et al. One-pass list-mode EM algorithm for high-resolution 3-D PET image reconstruction into large arrays. IEEE Trans Nucl Sci. 2002;49:693–699.

666

THE JOURNAL

OF

16. Visvikis D, Costa DC, Croasdale I, et al. CT-based attenuation correction in the calculation of semi-quantitative indices of [18F]FDG uptake in PET. Eur J Nucl Med Mol Imaging. 2003;30:344–353. 17. Visvikis D, Cheze-LeRest C, Costa DC, Bomanji J, Gacinovic S, Ell PJ. Influence of OSEM and segmented attenuation correction in the calculation of standardised uptake values for [18F]FDG PET. Eur J Nucl Med. 2001;28:1326– 1335. 18. Pavese N, Gerhard A, Tai YF, et al. Microglial activation correlates with severity in Huntington disease: a clinical and PET study. Neurology. 2006;66:1638– 1643. 19. Gebicke-Haerter PJ. Microarrays and expression profiling in microglia research and in inflammatory brain disorders. J Neurosci Res. 2005;81:327–341. 20. Banati RB. Visualising microglial activation in vivo. Glia. 2002;40:206–217. 21. Sapp E, Kegel KB, Aronin N, et al. Early and progressive accumulation of reactive microglia in the Huntington disease brain. J Neuropathol Exp Neurol. 2001;60:161–172. 22. Spinks TJ, Jones T, Bloomfield PM, et al. Physical characteristics of the ECAT EXACT3D positron tomograph. Phys Med Biol. 2000;45:2601–2618. 23. Gunn RN, Lammertsma AA, Hume SP, Cunningham VJ. Parametric imaging of ligand-receptor binding in PET using a simplified reference region model. Neuroimage. 1997;6:279–287. 24. Turkheimer FE, Edison P, Pavese N, et al. Reference and target region modeling of [11C]-(R)-PK11195 brain studies. J Nucl Med. 2007;48:158–167. 25. Woods RP, Mazziotta JC, Cherry SR. MRI-PET registration with automated algorithm. J Comput Assist Tomogr. 1993;17:536–546. 26. Beyer T, Townsend DW, Brun T, et al. A combined PET/CT scanner for clinical oncology. J Nucl Med. 2000;41:1369–1379. 27. Shao Y, Cherry SR, Farahani K, et al. Simultaneous PET and MR imaging. Phys Med Biol. 1997;42:1965–1970. 28. Aston JA, Cunningham VJ, Asselin MC, Hammers A, Evans AC, Gunn RN. Positron emission tomography partial volume correction: estimation and algorithms. J Cereb Blood Flow Metab. 2002;22:1019–1034.

NUCLEAR MEDICINE • Vol. 49 • No. 4 • April 2008

jnm041871-tp n 3/12/08