LSO background radiation as a transmission source using time of flight

The scanner used in these experiments is a Siemens Biograph mCT with 4 rings of detectors. The coincidence signal from Lu-176 decay is not seen in the standard mCT for a few reasons. The coincidence window on the mCT is ~4 ns which equates to a spatial acceptance diameter of ~60 cm. This removes most of the longer lines of response in the scanner. The line of responses that are less than the 30 cm in length are further discriminated by the lower level discriminator (LLD), which is set to 435 keV. In order to detect the decay from Lu-176, some modifications are made to the standard mCT. The coincidence window is increased to ~6.6 ns which results in a chord length of 99 cm and larger than the physical diameter of the scanner of 86 cm. The LLD is lowered to ~160 keV. The constant fraction discriminators (CFD) thresholds are also lowered to a value of 160 keV. In the processing firmware for the detectors, multiple energy windows were added to discriminate between the original emission 511 keV photons from a positron annihilation event and the two gammas (307 and 202 keV) from Lu-176. These additional energy windows were centered on 307 and 202 keV. The 307 keV window had an upper threshold of 350 keV and a lower bound of 250 keV. The 202 keV window had an upper bound of 250 keV and a lower bound of 160 keV. The events within these energy windows are tagged in listmode data and are used in the rebinner for energy discrimination.

The signal from Lu-176 as transmission data is measured by recording the beta in the originating detector. This beta will ionize its energy very locally in the LSO material and is accepted if it has enough energy to trigger the CFD. If one of the 307 or 202 keV gammas (the 88 keV gammas are neglected) escapes the originating detector, it has to traverse the field of view and be absorbed by an opposing detector to be recorded as a coincidental event. The PET scanner records the event's positions to create a line of response and records a time difference for the two events. By knowing the spatial positions of the two pixels that recorded the line of response, a look up table can be created to relate the chord length of the measured line of response (LOR) to the time of flight of the traversing gamma. Using this relation, a transmission coincidence time window can be created for each LOR (figure 2).

Zoom In Zoom Out Reset image size

The events are further processed by qualifying the energy of the gamma and using the information from the time of flight to get directionality of the line of response. Since the beta is emitted and captured locally in the LSO, the beta event happens first in time. Therefore the detector element with the smaller time stamp of the two events is where the beta emission occurred. The beta only has to trigger the CFD whereas the gamma must deposit enough energy to fall into one of the added energy window for the Lu-176 events. An additional qualification was added to the 202 keV sinogram to minimize emission backscatter contamination into the transmission data. If emission data is assumed to be 511 keV, the 180º backscattered energy deposited in the first hit detector would be ~340 keV and the Compton scattered gamma would carry an energy of ~170 keV. The time of flight of such an event is consistent with the Lu-176 events. A condition is added where if the second event (generally a gamma from a Lu-176 event) is in the 202 keV window, the first signal (beta from Lu-176) must not fall into the 307 keV window. If all conditions are satisfied, the event is added to one of the 2 sinograms for 307 and 202 keV events. In order to create the attenuation map from transmission data, a blank scan is performed. The blank scans are acquired for long time periods with no objects in the FOV to obtain fairly noiseless blank scans. They are rebinned in the same manner as described above and also separated into 2 separate sinograms depending on the detected gamma's energy. The typical rates found in the blank scans are ~400 events per minute per chord for the 307 keV energy window and ~250 events per minute per chord for the 202 keV window.

The attenuation maps are reconstructed with an ML-TR iterative algorithm with quadratic regularization that models the transmission data statistics (Fessler 1995, Nuyts et al 1998). The reconstruction of attenuation maps (rather than attenuation factor estimation) is necessary for a few reasons. First, LSO background data are low count nature and inconsistent due to presence of ignored scatter components. Therefore additional constraint of attenuation by regularized image reconstruction partially removes these inconsistencies. The attenuation map in PET is necessary to produce a scatter component, use in activity image reconstruction. Finally, LSO background attenuation map coefficients need to be scaled to 511 keV for activity reconstruction.

The prompt data with a mean value $\bar{p}$, spatial projection index i, can be modeled by combining LSO background source attenuated projections with random count mean values:

$$\begin{eqnarray}&&{{\bar{p}}_{i}}\left(\mu \right)={{a}_{i}}{{B}_{i}}+{{\bar{r}}_{i}},{{a}_{i}}={{e}^{-\underset{j}{\sum} {{c}_{ij}}{{\mu}_{j}}}}\end{eqnarray} \tag{ 1 }$$

Here, µ is the attenuation coefficient distribution, identified at the image voxel with index j. B represents blank projection data. c_ij is the system matrix of non-TOF line-integral projectors. Mean random counts are $\bar{r}$and a is the attenuation factor. Scatter contribution from LSO background transmission source as well from emission activity sources was ignored.

The attenuation map and emission activity images are estimated by maximizing the objective function:

$$\begin{eqnarray}&&L\left(\mu \right)-\beta U\left(\mu \right)=\underset{i}{\sum} \left({{p}_{i}}\ln {{\bar{p}}_{i}}\left(\mu \right)-{{\bar{p}}_{i}}\left(\mu \right)\right)-\beta \underset{j}{\sum} \underset{k\in {{N}_{j}}}{\sum} {{\left({{\mu}_{k}}-{{\mu}_{j}}\right)}^{2}}\end{eqnarray} \tag{ 2 }$$

The first term L represents the Poisson Likelihood function to match modeled $\bar{p}$ and measured p prompt data. The second term U represents a simple quadratic penalty (also called Gaussian prior) on attenuation map values. This simplest prior imposes smoothing on the attenuation map image in order to reduce noise influence. We used six closest neighbors in the voxel j neighborhood N_j in a prior design. A regularization parameter β was chosen using a trial and error approach for a particular data set. The update equation represents a scaled gradient ascent algorithm, where a gradient of (1) can be recognized in the nominator:

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \mu _{j}^{\left(0\right)}=0 \\ \mu _{j}^{\left(n+1\right)}=\mu _{j}^{\left(n\right)}+\alpha \frac{\underset{i}{\sum} {{C}_{ij}}{{B}_{i}}a_{i}^{\left(n\right)}\left(1-\frac{{{p}_{i}}}{{{B}_{i}}a_{i}^{\left(n\right)}+{{{\bar{r}}}_{i}}}\right)-\beta \frac{\partial U\left({{\mu}^{\left(n\right)}}\right)}{\partial {{\mu}_{j}}}}{D\underset{i}{\sum} {{C}_{ij}}\frac{{{\left({{B}_{i}}a_{i}^{\left(n\right)}\right)}^{2}}}{{{B}_{i}}a_{i}^{\left(n\right)}+{{{\bar{r}}}_{i}}}+\beta \frac{{{\partial}^{2}}U\left({{\mu}^{\left(n\right)}}\right)}{\partial \mu _{j}^{2}}}, a_{i}^{\left(n\right)}={{e}^{-\underset{j}{\sum} {{C}_{ij}}\mu _{j}^{\left(n\right)}}} \\ \end{array}\end{eqnarray} \tag{ 3 }$$

Here, α is the relaxation parameter used to accelerate convergence (we used a 1.5 value), D is the transaxial diameter of the reconstructed image support, and n is the iteration number. The summation over the voxel j index is a forward projector operation used in the computation of attenuation factor a_i.

When using the attenuation maps from the Lu-176 transmission data for corrections to 511 keV emission data, scaling the attenuation map to 511 keV energies is performed. The raw attenuation data is acquired separately for the 202 keV and 307 keV events. The reconstruction of the attenuation maps are also performed separately for the 2 transmission data sets. This allows for the scaling of the attenuation values to 511 keV with separate ratios of the total attenuation coefficients of water at the values of 511 keV and 202 or 307 keV. For this work, when emission and transmission data is collected simultaneously, only the 307 keV transmission data is used.

Figure 3 demonstrates that the scaling of the reconstructed attenuation map from the transmission data from Lu-176 is similar to that derived from a CT scan. There is a loss of fine structure detail using the transmission method. How this translates into PET emission data corrections will be investigated.

Zoom In Zoom Out Reset image size

After the initial estimation of the attenuation map from LSO background transmission imaging, attenuation is refined in simultaneous attenuation and activity reconstruction. In this work we used the MLACF algorithm presented in (Panin et al 2012) which extends the algorithm (Defrise et al 2014) for the case on non-zero background events. The scatter component is supposed to be known (since it cannot be updated without the attenuation map reconstruction) and the scaling parameter is defined by the initial estimate of the activity image. In the presence of a background the maximum likelihood estimate of the attenuation at fixed activity has no closed form expression. Therefore the algorithm approximates the attenuation by the maximizer of the non-TOF likelihood of the TOF summed data y_i. This approximate ML estimate of the attenuation factor is equal to the ratio between the background corrected measured data, and the current modeled projection data (equal to the forward projection of the current activity image). Inserting this ratio into the original TOF Poisson likelihood eliminates the attenuation unknown so that the cost function becomes a function of the activity only. Optimization of that function is done using a monotone iterative algorithm with complexity similar to the commonly-used MLEM.

The attenuation of lines-of-response, passing outside the support of the activity image cannot be recovered in MLACF. For such LORs, MLACF uses known attenuation from transmission scan.

TOF prompt data y (511 keV window) with spatial projection (LOR) index i and TOF bin index t can be modeled by combining the modeled projection $\bar{p}$ from the emission object f, corrected for scanner efficiency by a normalization array N and for attenuation a. The background events have a known mean $\bar{b}$, equal to the sum of the estimated efficiency corrected scatter and of the estimated randoms. The Poisson Likelihood objective function has the following form:

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} L\left(\mathbf{f},\mathbf{a}\right)=\underset{i,t}{\sum} \left({{y}_{it}}\ln \left({{a}_{i}}{{{\bar{p}}}_{it}}\left(f\right)+{{{\bar{b}}}_{it}}\right)-{{a}_{i}}{{{\bar{p}}}_{it}}\left(f\right)-{{{\bar{b}}}_{it}}\right), \\ {{{\bar{p}}}_{it}}\left(f\right)=\underset{j}{\sum} \frac{{{c}_{it,j}}}{{{N}_{i}}}{{f}_{j}} \\ \end{array}\end{eqnarray} \tag{ 4 }$$

where c_it,j is the TOF system matrix. In the following, quantities without index t denote TOF summed quantities, e.g.

$$\begin{eqnarray}&&{{c}_{i,j}}=\underset{t}{\sum} {{c}_{it,j}}\text{},\text{}{{y}_{i}}=\underset{t}{\sum} {{y}_{it}}\text{},\text{}{{\bar{p}}_{i}}=\underset{t}{\sum} {{\bar{p}}_{it}}\text{},\text{}{{\bar{b}}_{i}}=\underset{t}{\sum} {{\bar{b}}_{it}}\end{eqnarray} \tag{ 5 }$$

The attenuation factors that maximize the non-TOF likelihood of the TOF summed data y_i at given activity are given by:

$$\begin{eqnarray}&&{{\rm{a}}_{\rm{i}}}{\rm{ = }}\left\{ {\begin{array}{*{20}{l}} {{{\rm{a}}_{{{\rm{i}}^{\rm^{\prime}}}}}{\rm{,known\ attenuation}}}\\ {\frac{{{{\rm{N}}_{\rm{i}}}\left( {{{\rm{y}}_{\rm{i}}}{\rm{ - }}{{{\rm{\bar b}}}_{\rm{i}}}} \right)}}{{{{{\rm{\bar p}}}_{\rm{i}}}}}{\rm{ = }}\frac{{{{\rm{N}}_{\rm{i}}}{{{\rm{\tilde y}}}_{\rm{i}}}}}{{{{{\rm{\bar p}}}_{\rm{i}}}}}{\rm{ > 0}}} \end{array}} \right.\end{eqnarray} \tag{ 6 }$$

Here prime spatial projection indices denote LORs that have known attenuation. Substituting (6) into (4) leads to an objective function (the reduced log- likelihood as similarly defined in Defrise et al 2014) that can be maximized by the following iterative algorithm:

$$\begin{eqnarray}&&f_{j}^{\left(n+1\right)}=\frac{f_{j}^{\left(n\right)}}{\underset{{{i}^{\prime}}}{\sum} \frac{{{c}_{{{i}^{\prime}},j}}{{a}_{{{i}^{\prime}}}}}{{{N}_{{{i}^{\prime}}}}}+\underset{i}{\sum} \frac{{{c}_{i,j}}}{{{N}_{i}}}\frac{{{y}_{i}}}{\bar{p}_{i}^{\left(n\right)}}}\left\{\underset{{{i}^{\prime}},t}{\sum} \frac{{{c}_{{{i}^{\prime}}t,j}}}{{{N}_{{{i}^{\prime}}}}}\frac{{{y}_{{{i}^{\prime}}t}}}{\bar{p}_{{{i}^{\prime}}t}^{\left(n\right)}+{{a}_{{{i}^{\prime}}}}^{-1}{{{\bar{b}}}_{{{i}^{\prime}}t}}}+\underset{i,t}{\sum} \frac{{{y}_{it}}}{{{N}_{i}}}\frac{\left({{c}_{it,j}}{{{\tilde{y}}}_{i}}+{{c}_{i,j}}{{{\bar{b}}}_{it}}\right)}{{{{\tilde{y}}}_{i}}\bar{p}_{it}^{\left(n\right)}+{{{\bar{b}}}_{it}}\bar{p}_{i}^{\left(n\right)}}\right\},\\ &&\bar{p}_{it}^{\left(n\right)}\equiv {{p}_{it}}\left({{f}^{\left(n\right)}}\right)\end{eqnarray} \tag{ 7 }$$

Here n is iteration number. We initialize the MLACF iteration (7) with the activity image reconstructed with one OSEM iteration using an initial estimate of the attenuation from LSO background transmission

A phantom study was performed to demonstrate the ability to separate the emission and transmission data. A uniform phantom (20 cm diameter × 26 cm in height) with ~1 mCi of activity was placed on the bed next to a cold CT calibration phantom. Figure 4 shows the sinograms acquired simultaneously. The emission sinogram (bottom left) shows a 10 min acquisition and clearly shows the emission events with the shadowing of the cold CT phantom. The transmission sinogram (right figures) shows the collected data also for 10 min and separated by the gamma energies. In both transmission sinograms, both phantoms in the field of view are visible as well as the bed that the phantoms are placed on. All sinograms in the figures are displayed with the same polarity of greyscale.

Zoom In Zoom Out Reset image size

From figure 4, one can observe a noisier 202 keV sinogram when compared to the 307 keV sinogram. A study was performed to see if the transmission data was contaminated by the emission data and if the contamination is a function of emission source activity. 3 uniform phantoms (20 cm diameter × 26 cm in height) with varying activity from no activity to 0.5 mCi and 2.2 mCi were placed in the geometric center of the PET field of view. These phantoms were measured for 1 h and rebinned into the 2 transmission trues sinograms (prompts – delays). A profile across all the sinogram elements with summing over 100 angles and summing of all axial planes are shown in figure 5 for 202 and 307 keV transmission sinograms.

Zoom In Zoom Out Reset image size

There is some contamination in the 307 keV transmission data around the edge of the object which can be seen as a small mismatch between the profiles. The 202 keV transmission data has more visible contamination from emission data. The tails of the sinogram increases in events further away from the center of the field of view. This may be a result of when the backscattering angle approaches 90° and the energies are close to equal to each other and both events fall into the 202 keV window. The change of the profiles in the region outside the objects boundaries in the cold phantom when compared to the region outside the object in the 307 keV profiles could be a result of the 307 keV photons scattering in the object into the 202 keV energy window.

A study of human subjects was performed in order to observe how well the method performs for the complex structure of the human body. The scan duration times were set to 10 min to simulate fairly realistic scan times and minimize movement from the volunteers. Each volunteer was placed on the bed and inserted into the PET field of view with no activity in or around the PET scanner. A corresponding blank was acquired for 36 h and used for reconstructions of the attenuation maps. The reconstruction algorithm parameters used were 10 iterations with 24 subsets with some regularization for all human volunteer studies. No corrections were performed to the data that corrects for object originating physical effects to the transmission gamma such as scatter or attenuation.

Figure 6 shows two volunteer scans of the head with a carbon fiber head holder in the field of view. From the figure the sinus are visible and the head holder and outline of the head are well defined. Some high density regions are also visible such as parts of the skull and teeth.

Zoom In Zoom Out Reset image size

Figure 7 are reconstructed images of the torso regions of males with weights of ~70 kg. The study of the torso region was performed with both arms up and arms down in a relaxed position. Arms are usually up in a clinical scan unless it is not practical. The arms down case can suffer from truncation of the CT. Using transmission data from Lu-176, the field of view for the attenuation maps are matched to the PET field of view. From both studies, body outline is resolved and internal details such as lungs and the heart are visible. The two studies also differ as the arms down case had the volunteer laying directly on the carbon fiber bed and the arms up case had the volunteer laying on a foam mat between the body and the bed. The bed is easily seen in both images but only having the whole bed visible when the volunteer is lying on a foam mat that separates the body from the bed. The patient's arms up while lying on the foam mat is the typical clinical procedure for imaging in the torso region.

Zoom In Zoom Out Reset image size

Figure 8 shows a larger volunteer of weight of ~ 180 kg that would experience truncation within the CT FOV. This study was performed with arms down and 10 min. The circle illustrates the CT's 50 cm FOV and demonstrates even if this study was performed with arms extended overhead, truncation would still occur to this volunteer. It is observed that the body contour is still resolved and some internal structures are visible such as the lungs and heart, but not as clear as the smaller volunteers' case.

Zoom In Zoom Out Reset image size

Work was performed to show the PET reconstructed images if the attenuation maps from the simultaneous transmission scan were used for the corrections to the PET emission data. The first case was a striatal head phantom with a fillable water cavity for addition of activity to the phantom. The phantom was filled with 4 mCi of F-18 and placed in a carbon fiber head holder. The phantom was scanned using a standard head-neck protocol for duration of 10 min. A corresponding CT was performed before the PET scan was acquired. The PET scan was the same for all three cases where the emission data was rebinned for 10 min into time of flight sinograms. Transmission based attenuation maps were reconstructed with a blank scan with duration of 36 h and the same reconstruction protocol used in the cold human studies above. Figure 9 shows three cases of interest, CT corrected, 1 h Lu-176 transmission and 10 min of Lu-176 transmission data. The 1 h collect of Lu-176 data was performed to show the quality of the method with higher statistics and compare it to a more reasonable time for patient studies. The attenuation correction and scatter corrections were performed using the associated attenuation map and the PET emission reconstruction performed was OP-OSEM with time of flight using 2 iterations and 24 subsets.

Zoom In Zoom Out Reset image size

The PET emission data shows little difference between all three cases. Uniformity in all three cases also shows little differences demonstrating that the attenuation maps of this simulated brain scan are good enough to perform the corrections to the emission data even with 10 min of data.

An image quality phantom was scanned to extend the study to a torso sized object. The phantom was filled with Ge-68 in an epoxy matrix and had an activity approximately 2 mCi at the time of the measurement. The phantom has 6 spheres with 4 hot (4x activity concentration from background) and 2 cold spheres with a cold Styrofoam filled cylinder in the center of the phantom. The phantom was placed in the centered to the bore and set on top of the bed in a foam holder for this particular phantom. A CT was performed before the phantom was moved into the PET field of view. The listmode acquisition was performed for 30 min. The emission data was rebinned for 10 min acquisition time and all 30 min of transmission data for the 307 keV photons were rebinned for transmission data. The longer transmission scan was performed to account for the exclusion of 202 keV transmission data due to observed contamination of the 202 keV transmission data from the emission data.

Figure 10 shows the derived attenuation maps from the CT scan and from 30 min of transmission data from Lu-176. From the attenuation maps, it is observed that the cold cylinder in the Lu-176 transmission image is not well resolved. The corresponding PET emission images show some artifacts that come from having residual values in the cylinder that should be empty. This problem is not seen in the cold spheres because the spheres are filled with epoxy. The cross talk between the emission data and the attenuation map puts activity in the region where there should be air and no activity.

Zoom In Zoom Out Reset image size

The attenuation map created using the Lu-176 decays have some inaccuracies but are still useful for larger objects. These attenuation maps generally define the boundaries of the object being scanned fairly well. Using the Lu-176 attenuation maps, the scatter correction can be performed on the object and the resulting scatter correction sinogram is inputted to the MLACF algorithm. The Lu-176 attenuation map can also be used as a starting image for the attenuation estimate of MLACF. The image quality phantom's data from the previous section was reconstructed using the MLACF algorithm with 5 iterations and 24. The resulting emission and attenuation maps are shown in figure 11.

Zoom In Zoom Out Reset image size

The attenuation map shown in figure 11 is an estimate from MLACF and shows that the internal structures are well defined with respect to the attenuation map derived from just Lu-176 data (figure 10). There is some disadvantage to the MLACF attenuation map in that there is no estimation for the line of responses that have no emission data. The bed and the shell of the phantom are not recovered but do not seem to vary much from the starting image. The center hole is recovered and the emission reconstruction now has no emission contamination in the center cold region of the phantom as seen in figure 10. Although the images are similar, the comparison is challenging as the two emission images are reconstructed using different algorithms with different objective functions and convergence rates. An observation is that the uniform regions do appear uniform with no visible artifacts. The sphere recovery is similar between the two cases and the iterations were selected to try to achieve similar noise structure between the two cases.