A new method to correct the attenuation map in simultaneous transmission/emission tomography using Gd/Ga radioisotopes.

Reconstruction of the tomographic images without attenuation correction can cause erroneously high count densities and reduced image contrast in low attenuation regions. In order to solve the problem of photon attenuation, one needs to know the attenuation coefficient for the individual patient being studied. Therefore, we made an attempt to correct the attenuation map in simultaneous transmission/emission tomography with (153)Gd/(67)Ga using maximum likelihood method using the expectation maximization (ML-EM) algorithm to correct the transmission window for both the spillover and downscatter. Spillover fraction, scatter fraction and parameters for the scatter function (A, b and c) were determined experimentally and optimized using the optimization program written in IDL based on simplex theory. All measurements were performed on a Vertex gamma camera using the anthropomorphic thorax phantom for validation of data obtained by the proposed method. It was observed that without spillover and downscatter correction, the mean counts were 19.29 in liver and 26.90 in lung, whereas after after applying the corrections, the mean counts were reduced to 3.80 and 15.10 in liver and lung, respectively, which were close to true mean counts (liver 2.15 and lung 14.89). In this proposed method, we introduced the set of F(t)(spillover) and K(t)(downscatter) to account for the variations in projection pixels (f(t) and k(t)) with the density and thickness. The F(t) and K(t) were determined using the transmission data by an iterative process. The quantitative error was reduced by 98.0% for lung and 90.0% for liver when the corrected transmission images were obtained after the subtraction of spillover and downscatter fraction.

Introduction

Single photon emission computed tomography (SPECT) is a tomographic technique that allows three-dimensional (3D) visualization of biochemical or physiological processes in the human body by external detection of photons from an administered radiopharmaceutical. SPECT reveals the function of the body rather than its structure. However, its ability to depict a true activity map depends largely on the imaging properties of the SPECT system as well as on the methods used for converting the acquired data into values reflecting the activity distribution. The physical factors which severely affect the quantitative accuracy of SPECT images include attenuation and scatter. Scattered radiation may cause deviation in measurements relative to their true value. The magnitude of the error introduced by the photon attenuation depends upon tissue thickness, tissue path length and the type of the tissue. Reconstruction of the tomographic images without attenuation correction can cause erroneously high count densities and reduced image contrast in low attenuation regions (e.g. lung), which introduce errors when images are evaluated quantitatively.-[12]

The human body is mainly constituted of low atomic number (Z) elements, i.e. soft tissue and water; therefore, the dominant interaction process in the biological system is Compton scattering. Further, these factors are related to the density and composition of the body tissues, and it is therefore crucial to have access to individual attenuation maps to account for the attenuation and scatter effects when high quantitative accuracy is desired. Fakhri et al.[3] and de Vries et al.[4] have studied the impact of scattered photons during SPECT acquisition using various 99mTc compounds and found that the scattered photons have a relatively minor impact on the reconstructed image quality. However, during 67Ga SPECT imaging, large number of scattered photons appear in lower energy acquisition window as explained by Farncombe et al.[5]

In order to solve the problem of photon attenuation, one needs to know the attenuation coefficient or attenuation map for the individual patient being studied. Techniques have been developed to correct the attenuation and scattering in SPECT for single-energy photon emitters such as 99mTc and multienergy emitters such as 201Tl where there is no spillover in the transmission window due to overlapping of the emission and transmission windows. However, these techniques are not satisfactory for generating the patient-specific attenuation maps in the case of a multi-energy emitter such as 67Ga, in which the emitted photon energy window overlaps with the transmission window for 153Gd in addition to the downscatter. Owing to the strong influence of the downscatter and the spillover in transmission map during simultaneous emission and transmission tomography, it is necessary to correct the transmission map before reconstruction of the emission images.[6] It is therefore desirable to develop a correction method to produce a corrected attenuation map that can be used for attenuation correction in the emission reconstruction. In the present study, we made an attempt to correct the attenuation map in simultaneous transmission/emission tomography with 153Gd/67Ga radioisotopes using maximum likelihood method using the expectation maximization (ML-EM) algorithm.

Materials and Methods

For the proposed study, we used the ML-EM algorithm.-[78] There are two steps in this method: (1) E-step, which computes the expected projection data for the current estimate and is equivalent to a forward projection, and (2) M-step, in which the expected value is compared to the measured projection data and then the estimated reconstructed image is updated by effectively backprojecting the ratio of the measured to estimated projections. In the algorithm, geometric response, attenuation, and scatter can be accounted for and modeled. Thus, it gives a quantitatively more accurate reconstruction.

Attenuation of emitted radiation as it passes through the subject under study is accounted using transmission source and an attenuation map is obtained for image reconstruction. Attenuation map was estimated based on measurements of attenuation in a phantom.

We developed a method to correct the transmission window for both the spillover and downscatter. We used three windows, two for emission, i.e. 67Ga (93 and 184 keV), and one for transmission, i.e. 153Gd (100 keV). Multiple energy window scatter estimation was not carried out due to constraint in the available acquisition system (Philips/ADAC). Therefore, convolution subtraction was investigated. The available system for transmission measurement consisted of two 153Gd line sources in a dual head 90° configuration.

In order to minimize the depth dependence of resolution, attenuation, spillover fraction and downscatter, the Geometric Mean (GM) of conjugate views was considered.

The observed image in the transmission window is given by

g = g + fg + g ……………………………………(1)

where gois the observed transmission image, g is the transmission image without contamination, g is the lower energy window image, ft is the spillover fraction and g is the downscatter image in transmission window.

Experimental setup to determine spillover fraction

The measurements of the spillover fraction were performed using the ADAC Vertex gamma camera. A point source having 130 MBq of 67Ga in a 5-ml syringe was placed in air at 15 cm from the face of the detector. A medium energy collimator was fitted to the gamma camera and images were acquired using a 128 × 128 matrix for 84 seconds. This duration was chosen so that it could accumulate 6 million counts in air to see the spillover fraction in the transmission window.

Spillover fraction was determined by varying the energy window setting of the lower energy window of 67Ga as shown in [Table 1]. This was done by shifting the lower energy window of 67Ga and transmission window (no 153Gd present) from their estimated photopeak energy windows, i.e. 90 (14%) and 105 (15%). The counts in the emission and transmission windows were noted. The experimental setup is shown in [Figure 1].

Table 1

Spillover fraction (ft) (%) in the transmission window at different energy window setting

Figure 1

Experimental setup to measure the spillover fraction in air

Determination of spillover fraction and downscatter in transmission window

As a first approximation, the lower energy window was used to estimate the spillover fraction in air in the transmission window. Multiplying the lower energy emission image with the spillover fraction (f) derives a spillover image and it was subtracted from the transmission image. The upper energy window was used to estimate the downscatter (g) in the transmission window by convolving the upper energy emission image using bi-exponential function.[9] rather than the monoexponential function used for photopeak scatter correction[10].

GM data were used as first approximation assuming that the scatter function is spatially invariant. Therefore, the scatter function can be written as

S = 1/N{A * exp(–br) + exp (–cr)} …………………………… (2)

where r is radial distance, A, b and c are the parameters for scatter function, and N preserves total counts by normalizing S to unity.

Then, convolved upper energy window is scaled by a scatter fraction (k) before the subtraction of estimated downscatter from the spillover corrected transmission window. The estimated corrected transmission image (g) can be given by:

g = g – {ftg + k (S × g)} …………………………………(3)

where g is the image in 93 keV window, g is the image in 184 keV window and kis a scatter fraction (in transmission window). Spillover fraction, scatter fraction and A, b and c parameters for the scatter function were determined experimentally and optimized using the optimization program written in IDL [Figure 2] based on simplex theory. The Sum of Absolute Differences (SAD) between the measured spillover and downscatter image (gmt) and estimated spillover and downscatter image (ĝmt) is calculated and can be given by:

Figure 2

Schematic diagram for determination of A, b, c, kt and ft by using optimization program

Σ|g – ĝ|

There were considerable variations in the values for k and fwith some degree of interdependency. To reduce the variation and interdependency of k and f, we fixed the A, b and c parameters (A = 49.22, b = 1.56, c = 0.122) for the scatter function and repeated the optimization procedure. The results for the repeat optimization are given in Table 2. It was found that these parameters have slight variation compared to the results shown in Table 3. However, there was still considerable variation observed for different source depths and materials. These results were evaluated graphically by plotting ktand f against μd (μd is the attenuation path length obtained by multiplying the narrow beam attenuation coefficient μ by source depth d). The results for k with fixed scatter function parameters could be described by a single buildup relation (equation (4)) which is independent of material:[11]

Table 2

Values of kt and ft when A, b and c values are fixed for different materials at various depths

Table 3

Values of A, b, c, k and f values for different materials at different thicknesses (cm)

where A1, B, and β are constants and μd is attenuation path length.

The correlation between measured k values and those predicted from the buildup equation was high (r = 0.976). The chi-squared values obtained during curve fitting reduced from 0.034 to 0.031 for the data obtained with non-fixed versus fixed scatter function parameters. It was evident from the observations that f had a linear relationship with source depth, which differed significantly for the different materials. The variation of individual values from the best fitted line was reduced when the scatter function parameters were fixed. The pattern of linear slope was similar for different materials [Figure 3]. In an attempt to further improve the fdata, another optimization was performed using fixed values for A, b, c, while kvalues were obtained from the fitted curve [Figure 4]. The results for k suggest that externally measured transmission data can be used to estimate the scatter fraction pixel by pixel. However, it is not possible to distinguish depth from material using an external transmission measurement. Therefore, a single linear equation (a’ + b’ * μd) was used to fit the f data with moderate success (r = 0.70) as shown in Figure 5. The fitted values of f and k are given in Table 4.

Figure 3

versus μd(A, b, care fixed)

Figure 4

k versus μd(A, b, c are fixed)

Figure 5

Spillover fraction (f) versus μd

Table 4

Fitted values of f and k for different μT/2 for different materials

If the point source is measured from conjugate views and the geometric mean of those views is calculated, then the spillover fraction (ft) and scatter fraction (k) are approximately equivalent to those from a source at the center of the object.[12-13] Therefore, ft and kt in the transmission window are approximately equivalent to those from a source at the center of the object irrespective of the depth. The μd can be replaced by μT/2, (f= a’ + b’ * μT/2) in the linear equation and buildup equation (1) can be rewritten as:

where T is the total thickness of the object.

In most cases, B is set to (A’ – 1), which results in kt = 0 in air, outside the attenuating object, where e = 1. These values of f and k were used when applying the convolution subtraction method in the phantom study.

In air, scattered photons are not present in the beam; therefore, spillover fraction in the transmission window due to 90 keV would be expected to remain constant. However, in a medium where high-energy photons of 67Ga interact by Compton scattering, there is loss of energy and some photons are detected in the transmission window (105 keV) as well as in the lower energy window (90 keV). Scattering depends on the density of the materials; therefore, downscatter also increases with the density and thickness of the object. The scatter due to the 93 keV photon (i.e. photopeak scatter) also increases with thickness and density of material. As a result, the ratio of counts (spillover) for the 105 keV to 90 keV windows is a complex mixture of unscattered and scattered photons. Consequently, the spillover can no longer be considered constant. The spillover fraction in the transmission window (g) is estimated by the spillover fraction f multiplied by the lower energy window (g) of 67Ga. To account for regional differences in scattering material, f is replaced by the pixel by pixel spillover fraction F.

A two-dimensional matrix of spillover fractions (F) which is the set of f for the projection pixels can be calculated using the linear equation which is derived from the fitted curve (f vs. μd) [Figure 5]. Similarly, the buildup equation establishes a relationship between kand μd as shown in Figure 4. Therefore, k also varies with density and thickness of the object. The scatter fraction k is replaced by K to account for variation in scatter fraction pixel by pixel. Therefore, the equation for the corrected transmission image can be written as:

Solution of this equation is only possible by an iterative process since F and K are functions of path length as estimated by transmission measurement.

Validation of transmission data correction

The method was validated for 67Ga/153Gd emission–transmission tomography by variable density (anthropomorphic thorax) phantom using various parameters for scatter function (S), scatter fraction (k), spillover fraction (f) and spatial window width (SWT) 13% (13% of the field of view is exposed by the transmission window at any time) for a more realistic situation. Measurements were performed on a Vertex gamma camera (Philips dual head) with the medium energy collimator. A reference scan was acquired before the phantom study. The reference scan consists of two static images that are used to calculate the μT. These images represent the blank scans (zero attenuation). The detectors were positioned in a 90° position. The imaging table was removed and the line source of 153Gd was used to acquire the reference scan in air using 105 keV (15% window width) with a 128 × 128 matrix size. A true transmission scan was acquired using 128 × 128 matrix size for three photopeak energy windows with centerlines and widths of 90 keV (14%), 105 keV (15%) and 184 keV (20%) to compare with the corrected transmission scan. Simultaneous emission and transmission scans were performed using the 67Ga/153Gd. The matrix size and window settings were kept same to acquire the true transmission map, uncorrected transmission map and emission image of thorax phantom. Two syringes (20 ml) having 50 MBq and 110 MBq of 67Ga in 15 ml solution were placed in the lung and liver of the phantom, respectively.

Results and Discussion

The aim of using thorax phantom was to assess the transmission-based correction method in a more realistic imaging situation. The results obtained are shown in Table 5, and Figures 6 and 7. The results were assessed by defining two ROIs (region of interests) over the contaminated region of lung and liver in the uncorrected transmission scan. These ROIs were applied to the uncorrected, corrected and the true transmission scans and the mean counts were recorded [Table 5]. It was observed from the mean counts that where no spillover and downscatter correction was applied in the transmission scan (uncorrected) of liver and lung, the mean counts were 19.29 and 26.90, respectively. After applying the transmission-based spillover and downscatter corrections, the mean counts were reduced to 3.80 and 15.10 in liver and lung, respectively. These mean counts were close to mean counts of liver 2.15 and lung 14.89 of the true transmission scans. Therefore, using the transmission-based correction methods for spillover and downscatter in the uncorrected transmission scan, the quantitative error was reduced by 90.0% and 98.2% in liver and lung regions, respectively. The results were also evaluated by examining the count profiles through lungs and liver of the uncorrected, corrected and true transmission scans as shown in Figures 6 and 7. Before the correction, there was significant difference between the profiles of both liver and lungs of uncorrected and true transmission scans. After performing the transmission-based downscatter and spillover corrections, the profiles through the liver and lungs of the corrected transmission scan and true transmission scan showed excellent agreement. The results for the thorax phantom demonstrated the importance of using transmission data for accurate downscatter and spillover correction in the transmission window. The average downscatter and spillover fractions varied considerably from subject to subject and were difficult to correct without transmission data. Furthermore, it is inappropriate to use a constant spillover fraction and downscatter fraction when the object has non-uniform density. The relationship between μd, f and k for this combination of radionuclides (67Ga/153Gd) can be established and used to implement transmission-based spillover and downscatter correction, which takes into account the effect of variable tissue density and the object thickness on the spillover and downscatter distribution.

Table 5

Quantitative accuracy using transmission-based downscatter and spillover corrections in the anthropomorphic thorax phan

Figure 6

Comparison between uncorrected, corrected and true transmission scans (liver) for the anthropomorphic thorax phantom

Figure 7

Comparison between uncorrected, corrected and true transmission scans (lung) for the anthropomorphic thorax phantom

Conclusion

In the proposed new method, we introduced the set of F(spillover) and K (downscatter) to account for the variations in projection pixels (f and k) with the density and thickness of the objects. The F and K were determined using the transmission data by an iterative process. The quantitative error was reduced by 98.0% for lung and 90.0% for liver in thorax phantom when the corrected transmission images were obtained after the subtraction of spillover and downscatter fractions from the contaminated transmission images. It is therefore concluded that the quantitative error in the contaminated transmission image can be considerably reduced using the proposed transmission-based scatter correction method.