Correction of Partial-Volume Effect for PET Striatal Imaging: Fast Implementation and Study of Robustness

Principles of GTM Method.

The principles of the method are briefly recalled here, and further considerations can be found in the original article by Rousset et al. (11). We use the notations and assumptions of that article. The objects being acquired in PET are supposed to consist of N tissue structures D_i with homogeneous activity T_i (i = 1,…, N). The true activities T_i can be estimated from a set of N ROIs ROI_j, (j = 1,…, N) of volume V_j and mean activity t_j, drawn on the PET image. According to the fifth equation of Rousset et al., one can write: Eq. 1 where x is the space coordinate in the image and RSF_i(x) (where RSF is regional spread function) gives the contribution of structure D_i to the PET measure at x. Finally, if one considers t the vector of measurements t_j, T the vector of the true homogeneous activities T_i, and the N × N GTM with terms ω_ij defined as: Eq. 2 we can rewrite Equation 1 as t = GTM·T.

The ROI_j may be any regions and they must be different from one another to provide a full-rank GTM matrix. These regions may have some voxels in common. Although not required by the GTM method, an ROI_j set for which each ROI includes mainly 1 given tissue structure should be chosen for practical numeric considerations. If one takes the set of regions ROI_j to be equal to the set D_i, the diagonal terms of the matrix can be interpreted as the recovery coefficients of each structure and the off-diagonal terms as spillover coefficients. If the number of regions ROI is less than N, the problem is underconstrained and does not allow the identification of the true concentration. If the number of ROIs is greater than N, the problem is overconstrained and the true concentration may be estimated in the least squares sense (12).

Taking for granted the hypothesis that the objects being acquired consist of a finite number of structures with a homogeneous activity, the a priori information is introduced through delineation of the tissue structures D_i corresponding to the homogeneous objects from the MR image. This task is performed in the discrete sampling of the MR image and implies the resampling of either the ROI_j function or the RSF_i function in Equation 2 to account for the intermodality registration.

Computation of RSF.

For the current acquisition conditions, we will consider that the center of mass for the object studied (brain with striatum) is at about the center of both the MRI scanner and the PET scanner. This choice reduces the generality of the work but focuses the study on usual acquisition conditions to infer operative results. Estimation of the different RSF maps appears to be crucial. Two implementations of this estimation have been proposed in the literature. The first is the original one given by Rousset et al. (11), which simulates exhaustively the interactions between the emitting objects and the scanner by using a transaxial PET simulator (19,20) (sinogram GTM approach); the second is an estimation obtained by spatial filtering in the image space with a gaussian kernel (image GTM approach). The sinogram GTM approach was demonstrated in 2 dimensions, and because most PET acquisitions are now performed in 3D mode, an extension to 3 dimensions is needed. In this 3D implementation, we accounted for most of the acquisition and reconstruction phenomena and considered the sinogram GTM implementation as the reference. Conversely, we chose to implement the image GTM approach with the aim of obtaining a fast and easily operated method.

Sinogram GTM Approach.

In the original implementation of the GTM method, the RSF_i corresponding to the tissue structure D_i was estimated with a 2D noise-free analytic simulator. For the 3D extension of this implementation, we used the projection module of a 3D analytic simulator originally written for whole-body scanning on the ECAT HR+ scanner (Siemens/CTI PET Systems, Knoxville, TN) (21) and adapted for this study (22). Compared with previous studies (11,19,20), 2 new effects are accounted for in our analytic simulator that act on the spatial resolution of the scanner operated in 3 dimensions. First, the projection module explicitly considers any 3D lines of response (LORs) defined by any crystal pair in coincidence, instead of only 2D transaxial LORs. These LORs are convolved with the intrinsic point spread function (PSF) of the 2 crystals; this corresponds to 2 successive convolutions with a 2D crystal PSF that models its response and the depth of interaction effect. Second, the blurred LORs are resampled as they are recorded on the scanner.

The simulated sinogram corresponding to each D_i structure is reconstructed to form the RSF_i with the same parameters as those used for the PET data to which the GTM algorithm will be applied. Attenuation is not simulated because it may be considered constant in the brain soft tissues and does not contribute to any acquisition effect in our application.

Image GTM Approach.

In this approach, the intrinsic PSF of the scanner and the effects of the reconstruction are compounded in a global PSF in the image space. From the above remarks, it appears that the resolution varies in the field of view (FOV) of the PET scanner and that the global PSF should be modeled with a nonstationary kernel. Yet we modeled the global PSF with a stationary 3D elliptic gaussian kernel for 3 reasons. First, many contributions, such as that of Mawlawi et al. (17), report results with this simple model; second, the striatal regions are always located near the center of the PET FOV; and third, this choice eventually enables operative clinical implementation. The RSF_i(x) functions are obtained as the convolution of support of D_i with the global image PSF. The choice of a stationary 3D gaussian kernel allowed us to use a fast and accurate recursive technique to perform the gaussian filtering. We used the ideas of Deriche (23,24), who proposed approximation of the gaussian filter with an exponentially based filter family that can be implemented through a recursive algorithm. This made the gaussian filtering fast and its computation time independent of the gaussian kernel width.

Sinogram Approach.

The shape of the intrinsic PSF was estimated with a 2-mm-diameter collimated pencil beam of γ-rays coming from an ¹⁸F source. For a ring detector geometry, we defined the intrinsic PSF h_α(r,l) for a crystal belonging to ring r_c, at azimuthal location l_c along this ring, and where α is the angle between the direction of an LOR and the crystal normal direction. The PSF h_α(r,l) is parameterized with the sum of 2 gaussian functions: Eq. 3 to better reproduce the tails (A < 1, ε_r > 1, and ε_l > 1) of the individual detector PSF. The transverse width σ₁(α) of the gaussian functions increases when the angle α increases to account for the resolution degradation toward the edge of the FOV. The PSF widths σ₁(α) and σ_r were tuned by matching simulated and measured data of a reconstructed 3-mm-diameter sphere filled with ¹⁸F and located at various radial distances from the center of the scanner FOV. We did not use a 1-mm pointlike sphere, as the analytic simulator is not meant to simulate objects much smaller than the crystal detector width (of the order of 4 mm for the ECAT HR+ scanner).

Image Approach.

The image GTM approach uses a convolution with a 3D gaussian kernel. We acquired experimental data for a 1-mm-diameter pointlike sphere filled with ¹⁸F. It was located 5 cm from the center of the scanner FOV, and the data were reconstructed using the same protocol as for the clinical PET data that are candidate to PVE correction. Considering the reconstructed sphere as the impulse response of the system (acquisition and reconstruction process), we fitted the measured profiles with a stationary gaussian kernel using different widths for the axial and transaxial directions.

Monte Carlo Simulation.

We simulated 2 reference datasets: an ¹⁸F-l-dopa–like dynamic study and an ¹¹C-raclopride–like dynamic study. The numeric head phantom of Zubal et al. (25) was used; it consists of a volume of manually labeled anatomic regions segmented from an actual MRI acquisition. As is sketched in Figure 1, we selected 7 ROIs that usually show uptake in ¹⁸F-l-dopa and ¹¹C-raclopride studies to create the reference anatomic volume: caudate (right and left), putamen (right and left), white matter, gray matter, and skin/muscle structures. The numeric phantom is a 256 × 256 × 120 voxel volume with a voxel size of 1.1 × 1.1 × 1.4 mm.

• Download figure
• Open in new tab
• Download powerpoint

FIGURE 1.

The 7 compartments considered in Monte Carlo simulations—skin/muscle, gray matter, white matter, right caudate, left caudate, right putamen, and left putamen—shown with arbitrary gray levels in axial, coronal, and sagittal views. Plot shows TACs used for each structure in ¹⁸F-l-dopa Monte Carlo simulation.

From an actual protocol using ¹⁸F-l-dopa tracer with nine 10-min frames, a physiologic TAC was derived for each structure. From these TACs and the numeric phantom, we generated 1 PET dynamic study with a Monte Carlo simulator (26) for the ECAT HR+ scanner. This simulator is, of course, distinct from the analytic simulator used in the sinogram implementation of the GTM algorithm. We reconstructed the data according to the brain protocol used in our laboratory: The filtered backprojection algorithm that we use is 3D reprojection reconstruction (27), with a Hann apodizing window and the Nyquist frequency cutoff (0.2 mm⁻¹). The simulation did not consider random events, but attenuation and scatter were modeled and were corrected in the reconstruction process (28). The reconstructed simulated image size was a 128 × 128 × 63 voxel volume with an isotropic voxel size of 2.43 × 2.43 × 2.43 mm and contained 9 frames. The simulation was performed with realistic parameters for geometry (size and relative position of the structures in the FOV) and statistics (statistics for emission corresponding to a 333-kBq ¹⁸F-l-dopa injection).

Similarly, an actual ¹¹C-raclopride study enabled us to obtain a second set of TACs. The simulated dynamic study contained 23 frames of various durations ranging from 30 s to 5 min (statistics for emission corresponding to a 370-kBq ¹¹C-raclopride injection).

The Monte Carlo simulator has already been validated, as reported by Reilhac et al. (26). Nevertheless, a Monte Carlo simulation of a 3-mm-diameter sphere was performed for 2 reasons: first, to check the concordance between the PSF from the measured data and the PSF from the simulated data; second, to assess the goodness of fit of the gaussian kernel chosen in the image implementation approach.

Brain Phantom Acquisition.

We used an anthropomorphic head phantom. This phantom (Radiology Support Devices, Long Beach, CA) consists of a main cavity containing 4 small chambers that can be filled independently. The shape of the cavities mimics the cortex volume and the putamen and caudate nuclei of the striatum. The main cavity may be embedded in a synthetic head skull. The cavities are made of plastic, whereas the skull is a bonelike material that accounts for realistic acquisition conditions (shape, complexity, and attenuation). A PET dynamic acquisition was performed with the ECAT HR+ scanner. The internal cavities were filled with a solution of ¹⁸F (414 Bq·mL⁻¹), and the internal chambers received a solution of ¹¹C (332 Bq·mL⁻¹). The protocol consisted of 20 acquisitions of 1 min each and 1 late acquisition of 10 min. The day after, a transmission scan of 20 min was performed. The acquisition conditions corresponded to a typical ¹⁸F-l-dopa brain protocol. The images were reconstructed with the brain algorithm described above. The decay correction step of the reconstruction process was disabled to synthetically create a temporal evolution of activity in the image. Samples from the ¹⁸F and ¹¹C solutions were taken and counted in a well counter for cross-calibration. The anthropomorphic phantom underwent MRI on a Signa 1.5-T scanner (General Electric Medical Systems, Milwaukee, WI). The main cavity only was filled with water to yield good contrast between brain and striatum with a T1-weighted spoiled gradient-recalled acquisition in a steady state. This set of data enabled us to validate the whole process of the GTM correction method: Indeed, we registered the PET and MRI studies with an algorithm based on mutual information (29) and used the MRI acquisition to manually delineate structures of homogeneous activity. The walls of the internal cavities, which could be seen in the MRI acquisition, were not considered as part of the striatal nuclei.

Factors of Errors for Correction.

If one takes aside the noise amplification effect inherent in any deconvolution method, some specific factors of error modify the GTM correction performance even if the hypotheses of the method are fulfilled. The a priori information used by the GTM method involves delineation of the homogeneous D_i structures either with a visual criterion (as in this work) or with a criterion of an automatic method. We chose to infer the regions ROI_j in the MR image coordinates and from the D_i structures with j = i: This can be performed through the selection of voxels of D_i-corresponding RSF_i with value above half the maximum of RSF_i or through the morphologic closing of D_i structures. This choice implies resampling of the PET data to account for the registration between the PET image and MR image. Consequently, the PET–MRI registration and the segmentation of the structures D_i are considered parts of the PVE correction in the GTM method.

In considering the organization of our implementation and the sources of errors described above, one notices 2 points. First, the registration errors between the PET image and MR image introduce errors that affect only the measures t_j; the registration error does not modify the GTM matrix for a given choice of the structures D_i and regions ROI_j. Second, from Equation 2, it appears that potential error in the delineation of the structure of homogeneous activity will bring some artifacts in the computation of RSF_i and consequently in the ω_ij terms.

Robustness with Respect to Registration Error.

The misregistration is a major source of errors for the GTM method. Many multimodality registration algorithms for cerebral studies are currently available; most estimate a rigid-body transform between the images to be registered. In recent developments in registration methods, the matching criterion is most often a similarity measure at the voxel level and especially the mutual information criterion (29,30) that extends the early ideas of Woods et al. (31). Estimation of the registration error for a given algorithm is an active domain of study in itself. The signal-to-noise ratio, resolution in the images, and contrast in the images significantly affect the accuracy of the registration. In the comparative study reported by West et al. (32), a median error is defined to compare different registration algorithms, and typical error values are given for different conditions. Following this work, a median error ranging from 1.5 to 5 mm is expected for the studies currently performed in our laboratory based on images with the characteristics described in the Phantom Simulation and Acquisition section. This error range is confirmed by a study that we performed previously (33). Yet a simple visual quality check enables detection of a shift error greater than 2.5 mm. Consequently, the registration errors we studied corresponded to shifts ranging from 0.5 to 2.5 mm at the level of the striatum.

To be realistic, we generated rigid transforms that simulated the registration error of a given amplitude. We separated our simulations into 2 categories and considered either the translation-only error or the rotation-only error. For a given norm of translation, we generated 21 unit vectors that were directed toward randomly chosen points on the surface of a hemisphere, which were nontrivial directions of the image. For a given value of rotation error, we inferred the 21 corresponding rotation-only rigid transforms. To explore the range of current residual error yielded by any registration method (0.5–2.5 mm), we made the norm of the 21 translations vary from 0.5 to 2.5 mm and the angle of the 21 rotations vary from 1° to 5°. The effects of the registration error were studied with the ¹⁸F-l-dopa Monte Carlo PET simulation.

Robustness with Respect to Segmentation Error.

Missegmentation of the structures D_i is another major source of errors. The segmentation error cannot be modeled or quantified easily. To study the effect of missegmentation, we varied the shape of the structures D_i (34). The missegmentation process was obtained by delineating different sets of D_i structures. We chose 1 structure that varied—the right caudate—and left the others unchanged. To describe the segmentation error with a 1-dimensional parameter, we generated a set of structures D_{Right Caudate} and endowed them with good set properties. We generated several sequences of right caudate shapes D_{Right Caudate}^k, k = (1,…, 70) that followed the property D_{Right Caudate}^k ⊂ D_{Right Caudate}^{k + 1}. With this condition, we may parameterize the segmentation error of a given D_{Right Caudate} structure by its error size relative to the true one. Given the right caudate nucleus available in the Zubal phantom (3,000 voxels), we varied the structure from a 1-voxel-thick eroded shape to a 1-voxel-thick dilated shape. The relative error size varied from −30% to +30%. The whole process is sketched in Figure 2. We started from the initial right caudate and obtained the internal layer (respectively, the external layer), which is between the original shape and the 1-voxel-thick eroded shape (respectively, the 1-voxel-thick dilated shape). Each of the obtained thin layers was split automatically into small patches of 35 voxels each. Then, we progressively added the patches to go from the small eroded shape to the dilated one. Each time we added a patch, we checked that the object still contained only 1 connected component. To get relevant results, we defined 30 different sequences of evolutions. To quantify the effects of a segmentation error on the GTM correction procedure and compare them with the effects of other sources of error, we used the ¹⁸F-l-dopa–simulated PET data.

• Download figure
• Open in new tab
• Download powerpoint

FIGURE 2.

Flow chart illustrating building of family of missegmented right caudate. (A) Original caudate. On left, 1-voxel eroded core is shown in light gray. On right, 1-voxel-thick internal layer is shown in gray and 1-voxel-thick external layer, from original shape, is shown in dark gray. (B) Addition of one tile after another. (C) Tessellated caudate.

Influence on GTM Matrix.

Table 2 reports the 2 GTM matrices obtained for the numeric phantom used in the simulations. The structures D_i corresponded to the tissue structures with homogeneous activity that were considered in the simulation, and the regions ROI_j were the voxels with a value greater than 50% of the maximum value of the RSF_j(x) corresponding to D_j for j = (1,…, N). The diagonal terms reveal the recovery coefficients, and one can notice that smaller values are yielded by the image implementation, compared with the sinogram implementation, for the small central structures. If one considers the spillover coefficients greater than 0.15, which correspond to pairs of regions that actually interacted, the relative difference between the 2 approaches is less than 4%. For 3D PET, the computation time is about 1,000 times longer in the sinogram approach than in the image approach. To fulfill the comparison, we applied both implementations to the simulated phantoms (dopa and raclopride) and to the brain phantom acquisition. The D_i structures were taken from the Zubal model (Fig. 1) for the Monte Carlo simulation and from the MR image for the realistic brain acquisition. The ROIs were derived as described above. Two sets of PVE-corrected values were obtained, one for each implementation (sinogram and image).

Monte Carlo Simulation.

The PVE-corrected TACs of concentration of ligand for the dopa and raclopride studies were compared with the known curves introduced into the simulation. We used the apparent recovery coefficient (ARC), defined as the ratio of the recovered activity to the known activity for each structure, as a figure of merit to compare the 2 implementations. In Figure 4, we plotted the 3 ARC values (uncorrected, GTM image, and GTM sinogram) for each frame of the simulated dynamic PET scan. The mean and SD values of the ARC for the 2 studies are reported in Table 3. The mean and SD values are consistent with those reported by Rousset et al. (11); the values vary between 93% and 104% for the mean and between 0.1% and 3.4% for the SD. The SD values indicate a slight advantage for the image implementation.

• Download figure
• Open in new tab
• Download powerpoint

FIGURE 4.

Plot of ARC with time that indexes PET dynamic simulations. Regional ARCs shown are uncorrected (solid line), sinogram GTM PVE corrected (dashed line), and image GTM PVE corrected (dotted line). First row displays curves on left caudate and second row on right putamen. First column displays curves from ¹⁸F-dopa simulation and second column from ¹¹C-raclopride simulation.

Brain Phantom Acquisition.

We fitted a monoexponential function to the corrected TAC of each structure to obtain the isotope half-life and concentration of the isotope. As criteria for validation, we used ARC; isotope half-life recovery; and initial activity concentration recovery, defined by the ratio of the recovered values to the calibrated (counted) ones. The values of half-life recovery, initial activity concentration recovery, and ARC for each implementation are given in Table 4 and can be compared with the values obtained when no correction is performed. The values for half-life recovery or initial activity concentration recovery obtained after correction vary from 89% to 100% for the sinogram approach and from 93% to 99% for the image approach. These values are in accordance with values reported in the literature (11). Small differences between the 2 implementations may be noted. The ARC mean criterion indicates better raw performance for the sinogram approach; the ARC SD is smaller for the image approach. The 2 implementations therefore do not induce the same bias on the correction: The bias is constant with contrast for the image implementation, whereas the bias varies with a small amplitude for the sinogram implementation.

The PVE recovery capability of the method is reproducible for each structure regardless of the implementation. No definitive difference may be found to distinguish between the 2 implementations for this realistic striatal phantom.

Robustness with Respect to Registration Error.

For each simulated misregistration (translation or rotation) and for each structure, we computed the mean of the ARC (mARC) across all the frames of the dopa Monte Carlo dynamic study. To characterize these mARC values for a given amplitude of misregistration, we computed the mean and SD values of the mARC over all the drawings (21 transforms) for each given amplitude. Figures 5 and 6 show these values for the right caudate in the dopa simulation. The mARC mean value remains stable, whereas the SD increases from 2%–3% (with no registration error) to 6%–8% for a 2.5-mm translation error or a 5° rotation error.

• Download figure
• Open in new tab
• Download powerpoint

FIGURE 5.

Robustness of GTM method with respect to translation error between PET and MRI: example of right caudate structure with ¹⁸F-dopa simulation study. Evolution of mean (%) and SD of mARC value is plotted against amplitude of translation residual error. One mARC value is obtained for each randomly generated translation error (21 drawings per amplitude).

• Download figure
• Open in new tab
• Download powerpoint

FIGURE 6.

Robustness of GTM method with respect to rotation error between PET and MRI: example of right caudate structure with ¹⁸F-dopa simulation study. Evolution of mean (%) and SD of mARC value is plotted against amplitude of rotation residual error. One mARC value is obtained for each randomly generated rotation error (21 drawings per amplitude).

Robustness with Respect to Segmentation Error.

As described previously, the mean and SD values of the mARC were computed over all the configurations (30 configurations) for a given size of the missegmentation object. In Figure 7, these values are plotted for the right caudate in the dopa simulation. The mARC mean values vary from +7% to −7%, which corresponds to, respectively, a configuration with an undersegmented object and a configuration with an oversegmented object. The SD value of mARC remains stable. There is no spurious variation in the behavior of the correction method; the bias follows a piecewise decrease. For the unchanged objects, the GTM algorithm keeps performing correctly. Results for the unchanged objects, not reported here, showed that the bias of TACs after correction for those objects remained constant regardless of the shape of the right caudate.

• Download figure
• Open in new tab
• Download powerpoint

FIGURE 7.

Robustness of GTM method with respect to segmentation error: example of right caudate structure with ¹⁸F-dopa simulation study. Evolution of mean (%) and SD of mARC value is plotted against amplitude of segmentation residual error (relative error size). One mARC value is obtained for each randomly generated segmentation error (30 drawings per amplitude).