804 IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 30, NO. 3, MARCH 2011 Simultaneous Reconstruction of Activity and Attenuation for PET/MR André Salomon*, Andreas Goedicke, Bernd Schweizer, Til Aach, Senior Member, IEEE, and Volkmar Schulz Abstract—Medical investigations targeting a quantitative analysis of the position emission tomography (PET) images require the incorporation of additional knowledge about the photon attenuation distribution in the patient. Today, energy range adapted attenuation maps derived from computer tomography (CT) scans are used to effectively compensate for image quality degrading effects, such as attenuation and scatter. Replacing CT by magnetic resonance (MR) is considered as the next evolutionary step in the field of hybrid imaging systems. However, unlike CT, MR does not measure the photon attenuation and thus does not provide an easy access to this valuable information. Hence, many research groups currently investigate different technologies for MR-based attenuation correction (MR-AC). Typically, these approaches are based on techniques such as special acquisition sequences (alone or in combination with subsequent image processing), anatomical atlas registration, or pattern recognition techniques using a data base of MR and corresponding CT images. We propose a generic iterative reconstruction approach to simultaneously estimate the local tracer concentration and the attenuation distribution using the segmented MR image as anatomical reference. Instead of applying predefined attenuation values to specific anatomical regions or tissue types, the gamma attenuation at 511 keV is determined from the PET emission data. In particular, our approach uses a maximum-likelihood estimation for the activity and a gradientascent based algorithm for the attenuation distribution. The adverse effects of scattered and accidental gamma coincidences on the quantitative accuracy of PET, as well as artifacts caused by the inherent crosstalk between activity and attenuation estimation are efficiently reduced using enhanced decay event localization provided by time-of-flight PET, accurate correction for accidental coincidences, and a reduced number of unknown attenuation coefficients. First results achieved with measured whole body PET data and reference segmentation from CT showed an absolute mean difference of 0.005 cm 1 ( 20%) in the lungs, 0.0009 cm 1 ( 2%) in case of fat, and 0.0015 cm 1 ( 2%) for muscles and blood. The proposed method indicates a robust and reliable alternative to other MR-AC approaches targeting patient specific quantitative analysis in time-of-flight PET/MR. Manuscript received July 27, 2010; revised October 28, 2010; accepted November 15, 2010. Date of publication November 29, 2010; date of current version March 02, 2011. Asterisk indicates corresponding author. *A. Salomon is with Philips Research Europe, Molecular Imaging Systems Department, 52066 Aachen, Germany, and also with RWTH Aachen University, Institute of Imaging and Computer Vision, 52072 Aachen, Germany (e-mail: andre.salomon@philips.com). A. Goedicke and B. Schweizer are with Philips Research Europe, Molecular Imaging Systems Department, 52066 Aachen, Germany (e-mail: andreas. goedicke@philips.com). T. Aach is with RWTH Aachen University, Institute of Imaging and Computer Vision, 52072 Aachen, Germany (e-mail: til.aach@lfb.rwth-aachen.de). V. Schulz is with Philips Research Europe, Molecular Imaging Systems Department, 52066 Aachen, Germany, and also with RWTH/UKA Aachen University, Institute of Experimental Molecular Imaging, 52072 Aachen, Germany (e-mail: vschulz@ukaachen.de). Color versions of one or more of the figures in this paper are available online at http://ieeexplore.ieee.org. Digital Object Identifier 10.1109/TMI.2010.2095464 Index Terms—Attenuation correction, consistency conditions, positron emission tomography, reconstruction. I. INTRODUCTION TTENUATION correction in position emission tomography (PET) has become standard in many clinical applications and enables, e.g., a quantitative analysis of local tracer concentration by using standardized tracer uptake values (SUV) as a measure [1], [2]. Typically, a dedicated computed tomography (CT)-scanner is used (e.g., in a hybrid PET-CT scanner setup) to provide the anatomical information required to derive attenuation maps for quantitative PET reconstruction. There is also a trend to further improve the quality of diagnosis by replacing CT by another imaging modality offering information beyond just anatomical imaging: magnet resonance (MR) Tomography. However, besides the expected advantages of this new hybrid imaging modality, MR has the downside of not providing an easy access to the gamma photon attenuation. As an example, MR imaging of air-filled and cortical bone structures normally result in similar image intensity values although the attenuation for 511 keV photons is substantially different. This is because the MR signal strengths depend mainly on the induced change of nuclei spin energy levels and not on the density of the weakly bound electrons, which contribute to the attenuation mostly via Compton- and Rayleigh-scattering. There are recent attempts to overcome this limitation by including anatomical atlas information [3]–[5]. However, due to the fact that some tissue densities vary considerably among the population (e.g., lungs), the robustness of such approaches is probably not accurate enough for broad clinical use. Hofmann et al. [3] proposed a pattern recognition algorithm operating on a data base of MR and corresponding CT images. Although they achieved very promising results for human brain imaging, the accuracy for human whole body data has still to be investigated. An alternative way to overcome this substantial limitation for whole body PET/MR is to incorporate additional knowledge about the photon attenuation solely derived from the PET emission data, e.g., using so called consistency conditions [6] or simultaneous statistical reconstruction [7], [8]. The basic idea of this concept is that under ideal conditions (i.e., neglecting noise, etc.), each such measured PET sinogram can be explained by a unique attenuation-activity distribution pair [9]. If the activity distribution is a priori known, this approach leads to accurate attenuation map estimates as shown by Nuyts et al. [7], but in all other cases, where both the activity and the attenuation distribution have to be determined, inherent crosstalk between both has A 0278-0062/$26.00 © 2010 IEEE SALOMON et al.: SIMULTANEOUS RECONSTRUCTION OF ACTIVITY AND ATTENUATION FOR PET/MR shown to limit the overall success [7], [10]. This is caused by some general limitations of the linearization which is used for solving this “ill-posed” (see [11]) problem as shown by Bronnikov [12] and Krol et al. [13] regarding SPECT imaging. These limitations, i.e., circular symmetry, zero-line integrals and minimum spread of observation angles, resolution limit, and source influence (also referred to as crosstalk) in principle hold for PET imaging, too, but with ideal spatial and timing resolution of the scanner, the problem of simultaneous reconstruction can be solved by the Potato Peeler approach as used by Sidky et al. [14]. Using this, the limitations caused by circular symmetry and source influence vanish with a corresponding improvement of temporal and spatial resolution. Subsequently, as shown in [15], the incorporated time-of-flight (ToF) information enables the reduction of the source influence even at limited temporal resolutions of currently 580 ps FWHM. In the following, we propose a simultaneous reconstruction algorithm where attenuation coefficients of previously defined connected anatomical regions are adapted instead of performing modifications on a voxel-by-voxel basis. By that, the number of unknown attenuation coefficients is significantly reduced for the subsequent optimization process, leading to an increased robustness of the algorithm as well as additional reduction of crosstalk artifacts (mostly caused by source influence). Furthermore, this approach does not incorporate additional anatomical a priori information such as organ models and thus, it is not affected by the problem of individual tissue attenuation variability and unique patient anatomy. II. RECONSTRUCTION METHOD USING DATA CONSISTENCY AND ANATOMICAL INFORMATION A. Outline of the Reconstruction A simplified flow chart of the proposed method is shown in Fig. 1. Two data sets are fed into the iterative optimization scheme: the PET (list-mode) projection data and anatomical patient-specific information. The anatomical information is here provided by the segmentation into several regions as captured in the MR data. Thus, during a pre-processing step, anatomical structures are segmented and initial attenuation coefficients ( -values) from, e.g., water are assigned. As it is not even required at this point to accurately distinguish between the different tissue classes (e.g., lung, cortical bone, or soft tissue), the initial value for each region is typically set to the attenuation of cm ). Furthermore, it is also not water at 511 keV ( required to have one organ represented by one segment. The iterative part of the method starts with an attenuation corrected reconstruction of the activity distribution. In order to eliminate most of the cross-dependencies between activity and attenuation during the subsequent update of the attenuation, we here suggest performing several iterations of ML-EM incorporating the ToF information. In the second step, a gradient-ascent approach, realized by a ML algorithm, finally derives -map correction values for the subsequent update of the regions. Afterwards, the process starts again using the newly generated attenuation map until a predefined stopping criterion is fulfilled 805 Fig. 1. Iterative reconstruction scheme. (e.g., until a lower limit of the mean squared deviations between estimated and measured sinogram has been reached). B. Activity and Attenuation Estimation In the following, we denote the estimated activity in iteration of image voxel as and the estimated attenuation co. The attenuation ( -)coefficient in efficient of region as , respectively. Note that a sinogram voxel is referred to as is equivalent to a sorted list of coincident events containing spatial and additional ToF information given by the time delay between the detection of two coincident photons. Consequently, each index of a sinogram entry, denoted as , indexes both the spatial position of the corresponding detector crystals and the ToF information. In the following, the measured number of coincidences in sinogram entry is denoted as . The activity distribution can be obtained by using an update scheme as proposed by Wang et al. [16] which can be expressed in standard ML-EM-like notation as (1) with the estimated nonattenuated activity and the estimated attenuation in sinogram entry [7], where is the interaction length of the line-of-response corresponding to sinogram entry with voxel . The entries of the system matrix denote the sensitivity of voxel with respect to sinogram entry and depend on the ToF information but not the attenuation. The estimated number of accidental coincidences , depending on both the estimated activity and the attenuation, is considered later in detail. As shown in the following section, it can be determined, e.g., using the results from a Monte-Carlo (MC) simulation and an estimated distribution of random coincidences. For the latter, 806 IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 30, NO. 3, MARCH 2011 there are numerous methods dealing with the estimation and correction of random coincidences, e.g., using a simple model or by measuring delayed channel random estimates [17]–[19]. In the following section, we propose a method based on the emission data only, i.e., without using delayed channel techniques. The -coefficient for each region is updated according to mean correction values (2) where represents a correction factor derived from the the number of voxels forming region . emission data and Since the attenuation for each true coincidence is independent of its point of origin, the calculation of the correction factors for the attenuation distribution is similar to the reconstruction algorithms in CT [20] and does not incorporate any ToF characteristics . Thus, a (relaxed) gradient-ascent approach realized by the back projection according to (3) as proposed by Nuyts et al. in [7] can be used, where parameter is a relaxation coefficient, and is the diameter of the detector array. The iterative optimization loop is performed, until, e.g., the consistency error , defined as (4) does not significantly decrease anymore. Deriving an estimation for , which implicitly depends on the activity and attenuation is a crucial point for the success of the proposed method. As this has a significant influence on the overall accuracy, a deis provided tailed description of our approach to assess in the following section. C. Estimation of Accidental Coincidences %– % In clinical PET imaging, a significant amount of detected coincidences either incorporates scattered photons or originates from multiple decays events (random coincidences). An incorrect estimation of these undesired coincidences (accidentals) has an adverse influence on both the convergence behavior and the achievable accuracy and thus on the reconstructed activity and attenuation. It is therefore essential that the probability distribution of the scattered coincidences and the random coincidences , whose entries are denoted as and , as well as their total amount and have to be accurately estimated for their correction in each simultaneous reconstruction step. Accordingly, the total amount of accidentals in each sinogram entry can be formulated as (5) With respect to inhomogeneous detector crystal efficiencies and spatial variations in the dose and attenuation distribution, might significantly vary in sinogram entry , as investigated in [19] and [21]. An example for this is the detector shielding which reduces the amount of random coincidences near the upper and lower edges of the axial PET field-of-view (FoV) by about 25% in clinical studies. The relation between the spatial distributions of random coincidences and random single photons can be described as (6) where is the estimated number of random single photons in the detector crystal of the first incident of sinogram entry and of the second incident, respectively. Note that the distribution of the ToF information in (6) is assumed to be uniform inside and zero outside the corresponding coincidence window, since the detector hardware uses a similar scheme to limit the maximum delay time between two coincident photon events. is almost indeWith respect to scatter, we assume that pendent of the total decay rate in the patient. With this assumpcan be estimated by, e.g., a sequential Monte-Carlo tion, (MC) simulation of single decays and accelerated equivalent techniques as proposed in [22]–[27]. In contrast to this, the estimate for the total number of scattered coincidences , namely the actual fraction between scattered and true coincidences (scatter-fraction), is also dependent on the actual decay rate of each study [28]. In our implementation, we considered out-of-axial FOV scatter by extending the simulated FOV (continuously if not estimated in neighboring scan FOVs). Also the aforementioned lead shielding and resulting inhomogeneity of the random singles distribution is taken into account. The quality of the MC simulation results strongly depends on the current attenuation and activity distribution. While a first order estimation for the attenuation map can be easily derived from a homogeneously water-filled body shape, a similar modelling approach for the initial activity distribution typically will not be appropriate. As the activity distribution is normally nonhomogeneously distributed in the body, a scatter estimate derived from such a model would be too inaccurate for further use. Consequently, since no reliable start values for and can be inferred, is set to zero. For initialization of the randoms estimation, we further assume an uniform distribution for and set to % of the total number of measured coincidences. Starting with these values for the estimated accidentals, a first reconstruction of the patient-related activity distribution is performed which then can be used for the subsequent MC-based estimation of . Afterwards, based on the assumption that each true coincidence has its origin inside the patient, a reference of only the accidental coincidences can be obtained in a dedicated sinogram subset containing only those sinogram entries whose corresponding lines- or- tubes-of-response do not cross the patient body. In consequence, coincidences with SALOMON et al.: SIMULTANEOUS RECONSTRUCTION OF ACTIVITY AND ATTENUATION FOR PET/MR are caused by scatter and multiple decay processes and thus provide the basis of a further separation in random and scattered coincidences. Both have an inherently different spatial behaviour: the density distribution of scattered coincidences in typically decreases fast with increasing distance between the corresponding lines-of-response and the scattering object, while the density of the concurrent distribution from random coincidences slightly increases towards the borders of the axial PET FoV [19], [21]. Subsequently, the random coincidences are almost isolated from both the scattered and true coincidences in a which is defined by subset (7) By setting to a value satisfying (8) the following equation holds for : (9) Under this assumption, even if the current scatter estimate is erroneous, a fairly accurate estimation for the total number of random coincidences can be calculated by (10) Otherwise, if instead of would be used in the sums of (10), errors in the estimated scatter would have a more severe influence on which then may lead to a cyclic convergence behavior. Subsequently, the spatial distribution of random detected singles can be updated using (11) denotes all sinogram entries in subset which where correspond to detector crystal . The variation of activity and attenuation and thus its corresponding random single photon rates along the patient’s z-axis (e.g., between head and torso) is significantly higher than its local variation inside each detector ring. Furthermore, for a typical patient study, the number of detected coincidences does not allow for an accurate (i.e., noise-free) estimation of the underlying probability distribution. After estimating the number and distribution of random coincidences with (6) and (10), one can assess the total number of scattered coincidences by (12) 807 Fig. 2. IEC NEMA phantom with fillable spheres and an empty cylindrical insert. In order to reduce the influence of deviations between the esti, mated and the actual number of random coincidences on (12) uses only that part of the sinogram which is significantly influenced by scattered coincidences. Summarizing, not only the activity and attenuation, but also the scatter and random distribution are iteratively estimated in our approach. The scatter distribution in sinogram subset differs significantly from the distribution of random coincidences which allows to isolate a distribution of random coincidences by only using the aforementioned sinogram subset . Unlike the attenuation and activity distribution, the distributions of the scatter and random coincidences mainly consist of low spatial frequencies as found by Accorsi et al. in [23] and [24] which also allows for a robust simultaneous estimation of all four distributions. Next, the influence of the estimation of the accidental coincidences on the resulting attenuation coefficients is further examined and the convergence behaviour of the entire method is evaluated. III. PERFORMANCE ON MEASURED PHANTOM DATA For an evaluation of the applied estimation of the scattered and random coincidences, measured PET/CT data was used. Reference distributions, e.g., for scatter and random events, were derived from these PET data and compared to the simulation results. A. Evaluation of the Activity Reconstruction The estimation of activity, attenuation, and accidental (random scatter) coincidences as well as the attenuation dependent sensitivity were investigated. Accordingly, test data from a NEMA IEC PET phantom have been acquired using a Philips Gemini TF PET/CT scanner and FDG as radiation isotope (see Figs. 2 and 3). The phantom was equipped with four hot spheres and two cold spheres, all of different diameters, located around an (nonactive) air-filled cylinder. The activity ratio between the hot spheres and the background was chosen as 5:1. During a total scan time of 90 s, about 20.4e6 events were detected. Using a multilabel region growing approach (see e.g., [29]) the CT-map of the phantom has been segmented in 461 regions of various sizes. All attenuation coefficients except the surrounding air were initially set to cm ( -coefficient for water at 511 keV). The influence of accidental coincidences on the accuracy of the method was investigated by comparing the correction factors after a single iteration assuming the reference (CT-based) attenuation map. Here, without any compensation for scattered- and random coincidences, the correction factors showed 808 IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 30, NO. 3, MARCH 2011 Fig. 3. Transverse (upper) and coronal (lower) view of the CT map (left image) and reference activity (right image) of the measured IEC phantom. Fig. 5. Left column: Back-projected measured coincidences of sinogram subset and CT-based attenuation map (first row: transverse view, second and fourth row: sagittal view). Third and fifth row of left column: Corresponding curves of the estimated (dashed line) and measured (solid line) coincidences. Right column: Difference between the estimated and measured back-projected coincidences of sinogram subset . I I K Fig. 4. Correction values in transverse (upper) and coronal (lower) view without compensation for scattered and random coincidences (left images) and with Monte-Carlo scatter estimation and compensation for random coincidences (right images). high negative values outside as well as inside the phantom [see Fig. 4 (left)]. As a consequence, the algorithm tried to compensate the missing correction of accidental counts during further iteration steps by reducing the -coefficients to values of cm instead of 0.095 cm . The reason for that is the high level of accidental coincidences in PET-imaging which, in this case, ranged around 25% for scattered and 17% for random coincidences, meaning that only 58% of the detected events were true coincidences in this particular case. After including (Monte-Carlo based) scatter compensation for Compton- and Rayleigh-scattering and correction for random coincidences using (1), only minor correction factors were found [see Fig. 4 (right)] and the derived -coefficients stayed close to the expected value of water at 511 keV ( cm ) with maximum absolute deviations of 0.0006 cm . Referring to the accuracy of the estimation of accidental coincidences, Fig. 5 shows the back-projection of all measured (mostly accidental) coincidences in sinogram subset (i.e., ) overlayed by the CT-based attenuation map for orientation, as well as the difference between the estimated and measured back-projected accidental coincidences. Additionally, two numbered sections through the back-projection are presented. For each of these sections, Fig. 5 also shows line plots of the back-projected estimated accidental coincidences and the measured reference of subset . Generally, the results in can be extrapolated to the entire sinogram subset space , because of the low spatial frequencies in the distribution of accidental coincidences [23]. Thus, the insignificant deviations between the measurements and estimation as presented in Fig. 5 demonstrate the reliability of the entire estimation for this particular case. B. Multiple Iterations Evaluation For a first evaluation of the entire method sequentially measured PET/MR data was used containing about 6.88e7 detected coincidences (see MR intensity image in Fig. 6, first column). As the Plexiglas (PMMA) is not visible in the T1 weighted MR image, a dilation filtering of the region segmentation was performed to extend the corresponding attenuation map by the phantom wall thickness ( mm). Furthermore, since the patient table of the PET/MR significantly differs from the PET/CT system, a dedicated transmission scan was performed and the resulting table attenuation was added to each attenuation map estimate. All regions of the segmented attenuation map (further denoted as region label map, see Fig. 6, second row) were initially assigned to the -coefficient for water at 511 keV. The distortions observed in the MR image and the region label map near the edges of the axial FoV are due to the truncated RF exitation pulse (Fig. 6, right column). The computed correction factors and attenuation maps at four different iteration steps indicate a SALOMON et al.: SIMULTANEOUS RECONSTRUCTION OF ACTIVITY AND ATTENUATION FOR PET/MR Fig. 6. First row: T1 weighted MR image of the IEC phantom. Uniform background could not be produced as observed in [30]. Second row: corresponding region label map. Fig. 7. Intermediate estimates for the attenuation maps (upper row) and the resulting local correction factors (lower row) in iteration n f ; ; ; g. = 0 1 2 17 stable and consistent convergence (see Fig. 7), and remaining attenuation map deviations of less than 0.0021 cm in the for the large water-filled region and less than 0.003 cm air-filled cylinder have been found. However, some larger deviations occurred in the regions corresponding to the PMMA since those voxels are not surrounded by active material and thus are affected by general limitations of the method, i.e., “minimum spread of observation angles” as described by Krol et al. [13]. IV. PERFORMANCE ON CLINICAL PATIENT DATA For the following evaluation study on clinical patient data, we assume that the acquisition and subsequent segmentation of MR images in future combined PET/MR systems lead to a similar structural accuracy as the segmentation of CT maps used in this work. Thus, in order to investigate the general feasibility and accuracy of the proposed method to extract information about the 809 attenuation from the emission (PET) data, it is sufficient for the segmentation algorithm to take the structural information from a CT instead of an MR device. Accordingly, only the anatomical segmented information, but not the attenuation coefficients as measured by the CT were used for the iterative optimization process. PET data at five bed positions covering the whole chest region coincidences per scan position), have been acquired ( resulting in a reconstructed FoV of 587 587 518 mm . In a first step, a reference attenuation map was generated after registering the CT and PET volume. The bottom row of Fig. 8 shows the measured CT-based attenuation map of this study. The anatomical reference information derived via an automatic segmentation of the sequentially acquired CT map (using the aforementioned multi-label region growing approach [29]) provided an almost perfect registration between both modalities and consisted of about 3000 closed regions. Except of the patient table which position and attenuation distribution is a priori given, the attenuation coefficient for each region was initially set equivalent to water at 511 keV cm , see 1st row in Fig. 8). Analogous to ( the previously discussed reconstruction of measured phantom data, about 11 iterations of the entire method were performed until convergence was reached. Already after the first iteration, differences between lungs, bones and soft tissue were observed which might allow for a coarse classification of each region in soft tissue, lung/air, and bone. With more iterations (second row of Fig. 8 shows the result after the eleventh iteration), the estimation further converged towards the reference CT-map, ending up in a good agreement between both distributions. Corresponding to these different iteration steps, Fig. 8 shows the reconstructed activity distributions and Fig. 9 displays the intermediate correction values for the attenuation within the presegmented regions. Note that near the axial boundaries of the PET FoV the detector sensitivity is limited and thus not enough signal is available to derive reliable information for the proposed gradient-ascent approach. Thus, for the evaluation of (2), some upper and lower slices of the correction fields have been cropped. Fig. 10 depicts the patient- and scanner-specific sensitivity distribution, which can be described as [denominator in (1)]. Note that the accuracy of the reference attenuation map and thus its segmented regions is affected by organ motion effects introduced by the different scan times for both modalities. While (fast) CT image acquisition typically takes below 4 s per bed position, a PET scan requires about 90 s per bed position. Thus, the PET image always represents a motion averaged view whereas the fast CT scan typically represents a certain physiological state. Since we did not apply any motion correction during the processing, we neglected this effect accepting local inaccuracies due to respiratory and heart motion like the deviations at the lower edge of both lungs (marked with black arrows in Fig. 9). The overall convergence of the proposed algorithm is illustrated by the curve in Fig. 11, which shows the root summed squared error for iterations . The curve smoothly approaches a lower limit value, mainly resulting 810 IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 30, NO. 3, MARCH 2011 Fig. 8. Attenuation map and reconstructed activity distributions. Upper row: Initial estimation, second row: estimation after 11 iterations, third row: CT-based reference attenuation and corresponding reconstructed activity, and last row: difference between 11th iteration and CT-based reference. White arrows point to region averaged  coefficients. (from left to right: transverse, coronal, and sagittal transection). Fig. 10. Sensitivity distribution, represented as gray-scale from black to white, of five bed positions ranging between 0.00019 and 0.10162 based on detector geometry and CT-based attenuation of 511 keV (from left to right: transverse, coronal, and sagittal transection). Fig. 9. Correction field after first iteration (first row), after 11 iterations (second row), and after using CT-based attenuation map (third row). Blue values indicate an attenuation reduction, red values an increase of corresponding attenuation. from Poisson-noise in the list-data and inaccuracies of the reconstruction, such as an imperfect modelling of accidental coincidences and segmented regions. Furthermore, a clustering of the attenuation values into intervals of different tissue classes according to Fig. 12 was ap- Fig. 11. Curve plots the root summed squared error (RSSE) between the measured and the estimated sinogram (left plot) and the the percentaged number of correctly classified voxels inside the patient depending on the number of iterations performed. plied to illustrate, how many attenuating voxels (i.e., voxels inside the patient body) were correctly classified. As indicated by the curve shown in Fig. 11, an early optimum was already approached after only three iterations for this particular case. Finally, a quantitative analysis of the remaining (region averaged) error between measured (CT-based) and estimated attenuation SALOMON et al.: SIMULTANEOUS RECONSTRUCTION OF ACTIVITY AND ATTENUATION FOR PET/MR 811 Fig. 12. Reduced value space for -coefficients. Fig. 13. Region size dependent absolute error distribution between the estimated attenuation ( ) and the low-dose CT based -map ( ). map (Fig. 8) has shown an absolute mean difference of 0.005 cm % in the lungs, 0.0009 cm % in case of fatty tissue and 0.0015 cm % for muscles and blood pool. The averaged deviations between the estimated activity using the reconstructed and the reference attenuation are as follows: Bony tissue : %, blood and muscles : %, fatty tissue : %, and lung tissue : %. As shown in Fig. 13 the accuracy depends from the region sizes and indicate that in this particular case, regions greater than 1000 voxels ( ml) provide absolute errors lower than 0.02 cm . Similar good results (not explicitly shown here) have been found for studies with obese patients. Investigations with simulated patient data assuming different values for the temporal resolution between 300 ps and 10 000 ps FWHM indicated that already a resolution of about 1500 ps leads to a good converml. More accurate temporal gence of regions larger than resolutions ( ps) led to both increased convergence speed and more accurate attenuation coefficients. Finally, first evaluation studies in combination with measured MR data from a sequential Philips PET/MR scanner were performed. Fig. 14 shows slices of one of five investigated patient studies illustrating a couple of problems occurring using anatomical information from MR instead of CT. One is the reduced quantitative region homogeneity observed for MR images, affecting threshold based concepts aiming to separate regions of different tissue classes, such as the ribs and the lungs, as well as the backbone and the surrounding tissue. In four out of five investigated cases, the latter pair could not be accurately separated by the applied region growing approach. For the good results achieved in the present case, extensive fine-tuning of threshold parameters was required to avoid, e.g., junctions between the lungs and the backbone since both result in similar MR signal responses. Moreover, misalignment between PET Fig. 14. From left to right: T1-weighted MR image, corresponding segmented label map, initial attenuation estimate, and estimated attenuation after 12 iterations. Dashed ellipses denote cropped MR - images by the limited MR FOV corresponding label map extents are assigned to red pixels near the FOV edges in the second column. and MR was observed in some of the study cases, e.g., where the MR signal coil was removed prior to the PET scan. Another problem with MR data is caused by the limited MR FOV which especially in case of adipose patient studies and near the patient arms led to cropping artifacts as marked by the dashed ellipses in Fig. 14. In order to reduce these artifacts, an iterative ML-EM reconstruction was performed without attenuation correction (see resulting activity distribution in Fig. 15, first column). Then, using the derived activity distribution, a threshold based binary segmentation (threshold was set to the first minimum in the histogram of the activity image) was performed for each slice. Subsequently, morphological filtering was applied (see Fig. 15, second column) in order to close small defects and to fit the noncropped body contour. Here, a dilation with a spherical filter kernel (support radius: two voxels) and an erosion, also with spherical support (radius: six voxels), was used. Subsequently, the missing parts of the region label map outside the MR FOV were also segmented and finally added to the segmentation (see second column in Fig. 14—red regions near the edges of the MR FOV correspond to extended regions). After some iterations of the simultaneous reconstruction, reasonable attenuation coefficients were observed in the regions of main interest such as the back bone, the soft-tissue, and the lungs. However, some details could not be reconstructed correctly such as the attenuation of the ribs due to the respiratory motion during the PET scan. Also the attenuation estimate in the cranial bone contained inaccurate values caused by slight misalignments between the PET and MR scan of the head. In our implementation, the computation time needed for each iteration ranged between 30 and 60 min per bed position on a dual quad-core processor. However, using variance reduction techniques or computer hardware related acceleration, e.g., based on graphics processing units (GPUs), the algorithmic 812 IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 30, NO. 3, MARCH 2011 REFERENCES Fig. 15. From left to right: Reconstructed activity without attenuation correction, corresponding segmented body volume, initial attenuation corrected activity estimate, and reconstructed activity after 12 iterations. Dashed lines denote the borders of the extended MR FOV—corresponding label map extents are assigned to red pixels near the FOV edges in the second column. performance can be increased by large factors ranging between 10 and 100 as shown in [31]. V. CONCLUSION The proposed algorithm provides a robust and generic method to generate attenuation maps derived from ToF PET emission and properly segmented nonfunctional (MR, CT) data. Using a system modelling including scatter- and random coincidence compensation, a good convergence behavior was demonstrated. In our tests, remaining minor deviations between the estimated and measured (CT-based) attenuation distributions had no relevant effects on the finally reconstructed activity distributions. Especially for tissue classification, the observed fast algorithmic convergence encourages the combination of this method with other approaches in order to increase the reliability of an, e.g., emulated CT-based attenuation map [4]. A significant advantage of our method over others is that partial volume errors are partially compensated by the attenuation values estimated instead of being strengthened as in tabular methods which tend to generate high local deviations in the attenuation coefficients near tissue boundaries. This leads to a substantially higher degree of freedom for MR-sequence development and thus supports the constant attempt to further reduce the acquisition time. Summarizing, the proposed approach provides a promising and feasible method for attenuation map generation in combined PET/MR imaging. However, as also shown in this document, the MR data have to provide a certain (high) quality level in order to enable an accurate segmentation of different tissue regions. Accordingly, future investigations must be focused on novel MR sequences and segmentation methods. [1] C. C. Meltzer, P. E. Kinahan, P. J. Greer, T. E. Nichols, C. Comtat, M. N. Cantwell, M. P. Lin, and J. C. Price, “Comparative evaluation of MR-based partial-volume correction schemes for PET,” J. Nucl. Med. Mol. Imag., vol. 40, pp. 2053–2065, 1999. [2] M. A. Lodge, J. D. Lucas, P. K. Marsden, B. F. Cronin, M. J. O’Doherty, and M. A. Smith, “A PET study of 18 FDG uptake in soft tissue masses,” Eur. J. Nucl. Med. Mol. Imag., vol. 26, pp. 22–23, 1999. [3] M. Hofmann, F. Steinke, V. Scheel, G. C. Charpiat, J. M. Brady, B. Schölkopf, and B. J. Pichler, “MR-based PET attenuation correction—Method and validation,” IEEE NSS/MIC, vol. 1–2, 2007. [4] M. Hofmann, F. Steinke, V. Scheel, G. Charpiat, J. Farquhar, P. Aschoff, M. Brady, B. Schölkopf, and B. J. Pichler, “MRI-based attenuation correction for PET/MRI: A novel approach combining pattern recognition and atlas registration,” J. Nucl. Med., vol. 48, pp. 1875–188, 2008. [5] T. Beyer, M. Weigert, H. H. Quick, U. Pietrzyk, F. Vogt, C. Palm, G. Antoch, S. P. Müller, and A. Bockisch, “MR-based attenuation correction for torso-PET/MR imaging: Pitfalls in mapping MR to CT data,” Eur. J. Nucl. Med. Mol. Imag., vol. 35, pp. 1142–1146, 2008. [6] A. Welch, “Attenuation correction in PET using consistency information,” IEEE Trans. Nucl. Sci., vol. 45, no. 6, pp. 3134–3141, Dec. 1998. [7] J. Nuyts, P. Dupont, S. Stroobants, R. Bennick, L. Mortelmans, and P. Suetens, “Simultaneous maximum a posteriori reconstruction of attenuation and activity distributions from emission sinograms,” IEEE Trans. Med. Imag., vol. 18, no. 5, pp. 393–403, May 1999. [8] G. Glatting, M. Wuchenauer, and S. N. Reske, “Simultaneous iterative reconstruction for emission and attenuation images in positron emission tomography,” Med. Phys., vol. 27, pp. 2065–2071, 2000. [9] F. Natterer, “Computerized tomography with unknown source s,” J. Appl. Math, vol. 43, pp. 1201–1202, 1983. [10] V. Dicken, “A new approach towards simultaneous activity and attenuation reconstruction in emission tomography,” Inverse Prob., vol. 15, pp. 931–960, 1999. [11] J. Hadamard, “Sur les problèmes aux dérivées partielles et leur signification physique,” Princeton Univ. Bull., pp. 49–52, 1902. [12] A. Bronnikov, “Approximate reconstruction of attenuation map in SPECT imaging,” IEEE Trans. Nucl. Sci., vol. 42, no. 5, pp. 1483–1488, Oct. 1995. [13] A. Krol, J. Bowsher, S. Manglos, D. Feiglin, M. Tornai, and F. Thomas, “An EM algorithm for estimating SPECT emission and transmission parameters from emission data only,” IEEE Trans. Med. Imag., vol. 20, no. 3, pp. 218–232, Mar. 2001. [14] E. Sidky and X. Pan, “Mathematical formulation of the potato peeler perspective,” IEEE NSS / MIC, vol. 2, pp. 722–725, 2002. [15] A. Salomon, V. Schulz, B. Schweizer, A. Goedicke, and T. Aach, “Simultaneous reconstruction of activity and attenuation in multi-modal ToF-PET,” in Fully 3D Conf. Proc. Book, 2009, pp. 339–342. [16] W. Wang, Z. Hu, E. Gualtieri, M. Parma, E. Walsh, D. Sebok, Y. Hsieh, C. Tung, X. Song, J. Griesmer, J. Kolthammer, L. Popescu, M. Werner, J. Karp, and D. Gagnon, “Systematic and distributed time-of-flight list mode PET reconstruction,” IEEE NSS / MIC, vol. 3, pp. 1715–1722, 2006. [17] D. Brasse, P. E. Kinahan, C. Lartizien, C. Comtat, M. Casey, and C. Michel, “Correction methods for random coincidences in fully 3D whole-body PET: Impact on data and image quality,” J. Nucl. Med., vol. 46, pp. 859–867, 2005. [18] C. W. Stearns, D. L. McDaniel, S. G. Kohlmyer, P. R. Arul, B. P. Geiser, and V. Shanmugam, “Random coincidence estimation from single event rates on the discovery ST PET/CT scanner,” IEEE NSS / MIC, vol. 5, pp. 3067–3069, 2003. [19] M. E. Casey and E. J. Hoffman, “Quantitation in positron emission computed tomography: 7. A technique to reduce noise in accidental coincidence measurements and coincidence efficiency calibration,” J. Comput. Assist. Tomogr., vol. 10, pp. 845–850, 1986. [20] J. Nuyts, B. De Man, P. Dupont, M. Defrise, P. Suetens, and L. Mortelmans, “Iterative reconstruction for helical CT: A simulation study,” Phys. Med. Biol., vol. 43, pp. 729–737, 1998. [21] R. D. Badawi, M. P. Miller, D. L. Bailey, and P. K. Marsden, “Randoms variance reduction in 3D PET,” Phys. Med. Biol., vol. 44, pp. 941–954, 1999. [22] C. C. Watson, M. E. Casey, C. Michel, and B. Bendriem, “Advances in scatter correction for 3D PET/CT,” IEEE NSS / MIC, vol. 5, pp. 3008–3012, 2004. SALOMON et al.: SIMULTANEOUS RECONSTRUCTION OF ACTIVITY AND ATTENUATION FOR PET/MR [23] R. Accorsi, L. Adam, M. Werner, and J. Karp, “Implementation of a single scatter simulation algorithm for 3D PET: Application to emission and transmission scanning,” IEEE NSS / MIC, vol. 2, pp. 816–820, 2002. [24] R. Accorsi, L. Adam, M. Werner, and J. Karp, “Optimization of a fully 3D single scatter simulation algorithm for 3D PET,” Phys. Med. Biol., vol. 49, pp. 2577–2598, 2004. [25] C. Holdsworth, C. Levin, T. Farquhar, M. Dahlbom, and E. Hoffman, “Investigation of accelerated Monte Carlo techniques for PET simulation and 3-D PET scatter correction,” IEEE NSS / MIC, vol. 3, pp. 1500–1504, 1999. [26] C. Holdsworth, C. Levin, M. Janecek, M. Dahlbom, and E. Hoffman, “Performance analysis of an improved 3D PET Monte Carlo simulation and scatter correction,” IEEE NSS / MIC, vol. 2, pp. 13/53–13/57, 2000. [27] C. Holdsworth, R. Badawi, P. Santos, A. Van den Abbeele, E. Hoffman, and G. El Fakhri, “Evaluation of a Monte Carlo scatter correction in clinical 3D PET,” IEEE NSS / MIC, vol. 4, pp. 2540–2544, 2003. 813 [28] S. Surti, A. Kuhn, M. E. Werner, A. E. Perkins, J. Kolthammer, and J. S. Karp, “Performance of Philips Gemini TF PET/CT scanner with special consideration for its time-of-flight imaging capabilities,” J. Nucl. Med., vol. 48, pp. 471–448, 2007. [29] T. Aach and H. Dawid, “Region oriented 3D-segmentation of NMRdatasets: A statistical model-based approach,” in Proc. Vis. Commun. Image Proc., SPIE, 1990, vol. 1360, pp. 690–701. [30] B. Zhang, D. Pal, Z. Hu, N. Ojha, T. Guo, G. Muswick, C. Tung, and J. Kaste, “Attenuation correction for MR table and coils for a sequential PET/MR system,” IEEE NSS / MIC, pp. 3303–3306, 2009. [31] J. L. Herraiz, S. España, S. García, R. Cabido, A. S. Montemayor, M. Desco, J. J. Vaquero, and J. M. Udias, “GPU acceleration of a fully 3D iterative reconstruction software for PET using CUDA,” IEEE NSS / MIC, 2009.