Hybrid resampling method for fan beam spect

A method of resampling of a sinogram from a fan beam SPECT includes the steps of angular resampling of the sinogram using frequency modulation, and thence transfers resampling of the sinogram using linear interpolation. This permits the use of parallel geometry algorithms in quantitative work without the expected degradation of resolution or the introduction of artifacts.

TECHNICAL FIELD 
The present invention relates generally to reconstruction of an object from 
fan beam data obtained from sinograms, and more particularly to an 
improved method of resampling for fan beam SPECT for quantitative work. 
BACKGROUND OF THE INVENTION 
Fan beam reconstruction methods for single photon computerized tomography 
(SPECT) of the prior art may be divided into two general categories: (1) a 
direct approach, and (2) a resampling approach. In the direct approach, 
the fan beam data are the input to a fan beam back projector, which can be 
derived from the integral for the parallel beam back projector by a change 
of variables. One disadvantage to the direct approach is that the ramp 
filter convolution must be replaced by a direct evaluation of the 
convolution. This process, which is implemented efficiently by Fourier 
transforms (FT) in standard parallel beam geometry, is almost always less 
efficient because of the need for direct evaluation of the convolution. 
The resampling approach seeks to apply a resampling or rebinning technique 
which results in a parallel beam data set. The advantage of the resampling 
approach is in the fact that most tomographic algorithms have been 
developed for parallel geometry. The disadvantage of the resampling 
procedure, which is usually by linear interpolation, is that the procedure 
may degrade resolution and introduce artifacts. 
SUMMARY OF THE INVENTION 
It is therefore a general object of the present invention to provide 
quantitative fan beam SPECT using a hybrid resampling algorithm which 
permits use of parallel geometry algorithms without the expected 
degradation of resolution or introduction of artifacts. 
These and other objects will be apparent to those skilled in the art. 
The method of resampling of a sinogram from a fan beam SPECT of the present 
invention includes the steps of angular resampling of the sinogram using 
frequency modulation, and thence transfers resampling of the sinogram 
using linear interpolation.

DESCRIPTION OF THE PREFERRED EMBODIMENT 
The assumption is usually made that the projection image of a SPECT fan 
beam collimator closely approximates the X-CT fan beam line integral. 
However, the inventor herein has found that this is not strictly true. 
There is an additional factor due to geometric effects that causes the 
sensitivity to vary across the collimator face. Thus, while the resampling 
by itself corrects for the distance dependent component of sensitivity, it 
does not correct for the transverse component. As discussed in more detail 
hereinbelow, this transverse component can be corrected by a post 
processing of the resampled projections. However, on balance, the direct 
methods do not correct for the transverse variation in sensitivity. 
The varying fan beam collimator resolution has been characterized utilizing 
collimator modulation transfer function (MTF) analysis. With line source 
measurements, it has been found that both the perpendicular FWHM 
(full-width-at-half-max) and parallel FWHM were invariant over the field 
of view. The parallel, or axial, FWHM varied linearly with distance at 
approximately the same rate as a parallel hole collimator of the same 
description. The perpendicular (or transverse) FWHM had an additional 
increase that was a monotonically increasing function of the source to 
collimator distance, and was directly related to the diverging beam 
geometry. As discussed below, the depth dependent resolution of the fan 
beam collimator transforms via resampling to a parallel beam depth 
dependent model. For this reason, a frequency filter based upon the 
frequency distance relation (FDR) is possible. 
For a fan beam collimator of thickness s.sub.2, focal length F with 
aperture function g.sub.1 and detector plane function g.sub.2, and a point 
source of strength f.sub.0 at a distance s.sub.1 from the aperture plane, 
and a transverse fan beam distance u.sub.0 from the optical axis of the 
collimator, the total counts per second tc obtained is: 
##EQU1## 
The tilde indicates the 2-D Fourier transform of the aperture function. The 
first fraction expresses the effect of collimator thickness on 
sensitivity. The second fraction is the term expressing the distance 
dependence, and the third fraction expresses the transverse dependence. 
The distance dependence is rectified by the resampling. The transverse 
component can be removed either at preprocessing or post processing. (See 
FIGS. 1 and 2). 
Let g.sub.1 and g.sub.2 be circular apertures of radius R.sub.1 and R.sub.2 
respectively. The 2-D Fourier transform of the circular aperture is the 
bess-sinc function and along the central ray, the result is: 
##EQU2## 
For purposes of measuring sensitivity only, the space sampled by a fan beam 
bore can be replaced with a strip integral of uniform width. As a fan beam 
line integral is traced back to the focus, the spatial density of the 
strip integrals increases, so the sensitivity increases. From geometrical 
reasoning along, this overlap is proportional to (F+s.sub.2)/(F-s.sub.1). 
At present, both resampling and direct methods are based on line 
integrals. No method of either type has been proposed that corrects for 
the transverse component in sensitivity. 
Both resampling and direct methods are based on line integrals. This means 
that the data is independent of geometry. The sampling of the sinogram for 
fan beam scanning is a two-dimensional coordinate transformation. The fan 
beam data is obtained from a symmetric fan beam collimator in a circular 
orbit. The data consists of N.times.N planar projections, taken at M 
equiangular stops in 360 degrees. 
The data is a fan beam data set f(u,.theta.) where u is the transverse 
coordinate, and .theta. is the angular coordinate. For definiteness, let 
us say that the acquisition is collected in a counterclockwise rotation, 
head view. The sinogram has been sampled at equal intervals in u and 
.theta., and is given by: 
EQU .function.(u.sub.j,.theta..sub.k), u.sub.j =(j-N/2+1/2).DELTA.u, j=0, . . . 
, N-1, and .theta..sub.k =k.DELTA..theta., k=0, . . . , M-1.(Eq. 3) 
The quantity .DELTA.u is the pixel size, and 
##EQU3## 
is the angle between projections. 
The fan beam resampling problem for noncircular orbits is intrinsically two 
dimensional. In other words, x=x(u.sub.j,.theta..sub.k) and 
.phi.=.phi.(u.sub.j,.theta..sub.k). On the other hand, if the orbit is 
circular, the sampling uncouples, and the coordinate systems for the flat 
fan beam collimator are related by: 
##EQU4## 
In FIG. 2, the Radon (parallel) space coordinates are (x,.phi.). The 
parameter R is the distance from the center of rotation to the focal point 
of the collimator. We will consider only the case for circular orbits, 
with R constant. 
For circular orbits, the two dimensional (2-D) resampling can be separated 
in two one dimensional (1-D) resamplings that can be performed in either 
order. In the process of the present invention, an algorithm is used that 
first performs the angular resampling using frequency modulation, followed 
by transverse resampling using linear interpolation. 
For angular resampling, first take the Fourier transform along the columns 
of the sinogram, with respect to .theta., and determine the phase factors 
exp (imtan.sup.-1 (u.sub.j /F)), where m is the angular frequency. The 
resampling algorithm is as follows: 
1. Determine .DELTA..phi.=.DELTA..theta.. 
2. Take the discrete Fourier transform (DFT) with respect to .theta.. We 
obtain: 
##EQU5## 
3. Use frequency modulation based on Equation 2 to obtain a translation in 
sinogram space: 
##EQU6## 
If the rotation is clockwise, it is easy to see that: 
##EQU7## 
4. Inverse transform to obtain the resampling: 
##EQU8## 
After multiplying by this phase factor and inverse transforming, the result 
is an intermediate sinogram that is a function of u and .phi.. This 
intermediate sinogram is defined as: 
EQU .function..sub.l 
(u.sub.j,.phi..sub.k)=.function.(u.sub.j,.theta.(.phi..sub.k)).(Eq. 8) 
If M is an integral power of 2, use the Fast Fourier Transform (FFT) for 
steps 2 and 4. If M is not an integral power of 2, use the arbitrary 
factors algorithm that is based upon the FFT and the Z transform, or a 
prime factor FFT. Zero fill should not be used because the sinogram is a 
periodic function of .theta.. 
Also note that the transfer function of angular resampling, (Equation 6), 
results in an MTF that is an ideal low pass filter. This helps explain the 
improved resolution of this method over other types of resampling. 
The transverse resampling is accomplished using linear interpolation. The 
interpolation algorithm is: 
1. Determine .DELTA.x from .DELTA.x=R/N. (Eq. 9) 
2. For x.sub.j =(j-N/2+1/2).DELTA.x, j=0, . . . , N-1, find the interval 
u.sub.J,u.sub.J+1 ! such that u.sub.J .ltoreq.x.sub.j &lt;u.sub.J+1. 
3. The interpolated result is given by: 
##EQU9## 
The resampled sinogram is 
EQU .function..sub.2 (x.sub.j,.phi..sub.k)=r.sub.j .function..sub.1 
(u.sub.J,.phi..sub.k)+(1-r.sub.j).function..sub.1 
(u.sub.J+1,.phi..sub.k)=r.sub.j 
.function.(u.sub.J,.theta.(x.sub.j,.phi..sub.k))+(1-r.sub.j).function.(u.s 
ub.J+1,.theta.(x.sub.j,.phi..sub.k)) (Eq. 11) 
A necessary step for the utilization of noise dependent frequency space 
filters such as the 2-D or 3-D Wiener filters is the accurate 
determination of the variance of the DFT of noisy data. The raw data is a 
radionuclide sinogram image, so each image element is an independent 
Poisson random variable. First, facts are presented that will be needed 
for the determination. The results are described as 1-D continuous 
functions, for simplicity. These results correctly generalize to 2-D and 
3-D discrete functions. 
1. Both the DFT and the continuous Fourier transform (CFT) are utilized. 
The DFT and its inverse are defined by: 
##EQU10## 
The CFT and its inverse are defined by: 
##EQU11## 
The set X is the domain of the variable x, that for our purposes is a 1-D 
DFT, a 1-D CFT or 2-D versions of these. An "overloaded" notation will 
also be used to describe a one dimensional DFT, a one dimensional CFR or 
two dimensional versions of these. Thus, .function.(x.sub.i,.omega..sub.k) 
is a 1-D DFT with respect to the second variable, .function.(.omega.,v) is 
a 2-D CFT, etc. 
The only difficulty with this notation is that the function arguments must 
be stated explicitly. Recall that, because of periodicity, the CFT with 
respect to the angular variable of the sinogram is identical with the DFT, 
with .omega.=.omega..sub.k =k, and k an integer. 
2. The expectation or mean of an ensemble of random variables is defined as 
.function..sub.l =(.function..sub.l). The operator(.) is the expectation 
operator based on the underlying probability density function of the 
stochastic process. Similarly, for the Fourier transform, we define 
.function.(.omega.)=(.function.(.omega.)), and 
.function.(.omega.)=.function.(.omega.). (The mean of the FT is equal to 
the FT of the mean). 
3. The quantity z* denotes the complex conjugate of z. 
4. The covariance of the FT of a r.v. is defined by: 
##EQU12## 
The variance of the FT is then obtained from the covariance by setting 
.omega..sub.1 =.omega..sub.2 : 
EQU Var(g(.omega.)).ident.def Cov(g(.omega.)g(.omega.)). 
5. If g(x) is an independent Poisson random variable at each position x, 
then: 
EQU Cov(g(x.sub.1)g(x.sub.2))=g(x.sub.1).delta.(x.sub.1 -x.sub.2) with 
Cov(g(.omega..sub.1)g(.omega..sub.2))=g(.omega..sub.1 -.omega..sub.2), and 
Var(g(.omega.))=g(0). (Eq. 13) 
6. Let 
##EQU13## 
the convolution of g with a linear operator A. Then: 
EQU Cov(h(.omega..sub.1)h(.omega..sub.2))=A(.omega..sub.1)A(.omega..sub.2)g(.om 
ega..sub.1 -.omega.)=A(.omega..sub.1)A(.omega..sub.2)g(.omega..sub.1 
-.omega.). (Eq. 14) 
If A is the shift operator h(x)=Ag=g(x-x.sub.0), then 
A(.omega.)=exp(-i.omega.x.sub.0), and 
EQU Cov(h(.omega..sub.1)h(.omega..sub.2))=exp(-i(.omega..sub.1 
-.omega..sub.2)x.sub.0)h(.omega..sub.1 -.omega..sub.2), and 
Var(h(.omega.))=h(0). (Eq. 15) 
Comparing this with Equation 13, we see that FTI using the angular 
resampling algorithm of Equations 5-8 leaves the variance of FT unchanged. 
Thus, we expect to see no change, on average, to the variance of the FT 
after FTI filtering. The Monte Carlo simulations verify that conclusion. 
To obtain the result for the combined resampling, we note from Equation 15 
that frequency modulation results in a shifted Poisson random variable. In 
other words: 
EQU Cov(h(.omega..sub.1)h(.omega..sub.2))=exp(-i(.omega..sub.1 
-.omega..sub.2)x.sub.0)h(.omega..sub.1 -.omega..sub.2) implies that 
EQU (10b) Cov(h(x.sub.1 -x.sub.0)h(x.sub.2 -x.sub.0))=h(x.sub.1 -x.sub.0) 
.delta.(x.sub.1 -x.sub.0 -x.sub.2). (Eq. 16) 
Thus, we can consider the output of the FTI angular resampling to be a 
shifted ensemble of independent Poisson random variables. Hence, we can 
apply Equation 16 (below), with interpolation weights given in Equation 10 
to determine the combined effect of angular FTI and transverse linear 
interpolation. 
It is also clear that the Modulation Transfer Function (MTF) of frequency 
modulation is equal to one, so frequency modulation is an ideal low pass 
filter. 
Another result needed is the variance of a function resulting from 
multiplication by a real scalar, a. From Equation 12, the result is: 
EQU Var(ag(.omega.))=a.sup.2 g(0). (Eq. 17) 
If the action is a vector multiply of the image g by a vector a, with 
h(x.sub.j)=a.sub.j g(x.sub.j), the variance of the FT is: 
##EQU14## 
The Parseval summation formula for the 2-D DFT is: 
##EQU15## 
The simulated fan beam projections were generated for a "Liver-Spleen" 
phantom. The pixel size was chosen to be 5.85 mm, so that truncated 
projections could be avoided. The resampled pixel size is then 2.70 mm, 
which is comparable to parallel ray simulations. The processing steps 
were: 
1. Resample the fan beam projection data to parallel beam projection data, 
using the two step procedure outlined above. 
2. Correct for the transverse component of sensitivity using Equation 22 
below. 
3. Generate contours based upon an algorithm that back projects edge 
information to form the contour of a convex body outline. 
4. Premultiply. This operation is the first step in the intrinsic 
attenuation correction. This algorithm accurately accounts for the 
"shadow" of an off-center convex attenuating body. This procedure makes it 
possible to accurately apply exponential back projection to complete the 
reconstruction. 
5. Take the 2-D Fourier transform the sinogram, then remove the "far-field" 
component of the sinogram. This filter is described as follows: 
##EQU16## 
sinogram frequency. 
6. Complete the reconstruction by exponential back projection, by either a 
direct image space method or Fourier-Bessel reconstruction using the CHT 
algorithm. 
The 2-D Wiener prefilter is the convolution of a gaussian, representing 
collimator resolution, with (modified Bessel) K.sub.o function for scatter 
compensation. The coordinates (x,y) are the Cartesian coordinates of the 
planar images. The parameters g.sub.x and g.sub.y are the 
full-width-at-half-max (FWHM) of the collimator resolution with respect to 
x and y, respectively. The parameter .beta. is related to the FWHM of the 
scatter function. The parameter SP is the scatter-to-primary ration. The 
symbol ** represents 2-D convolution. The filter kennel is: 
##EQU17## 
The evaluation of the K.sub.o function in frequency space avoids having to 
evaluate the removable singularity of the K.sub.o function at zero. This 
filter has been implemented for nonisotropic pixel sizes that is required 
by the resampling. The FT variance is utilized as the estimate of the 
noise power spectrum term in a Wiener filter previously described and 
applied to a 3-D resolution recovery filter based on the FDR for parallel 
geometry. 
A critical question is determining the probability distribution function of 
the FT noise variance. The formulation for a nonuniform distribution of 
independent Poisson random variables is needed, a nontrivial number of 
which will have a mean at or near zero. This was approximated with a 
nonuniform distribution of gaussian variables. It is shown that this model 
accurately predicts the experimentally derived NPS of the data. The 
important fact for this application is that the noise power spectrum is a 
constant. 
In FIG. 3, the total counts are shown from a resampled point source scan. A 
point source 1 cm in diameter was filled with about 400 .mu.Ci of Tc-99m 
and placed about 5.5 cm from the center of rotation. The radius of 
rotation of the camera was set at 19.45 cm. The matrix size was 128, with 
pixel size 3.56 mm. There were 120 projections in 360 degrees, in 
step-and-shoot mode, with 20 sec. per stop on a Picker 3000 system with 
UHRFAN beam collimators. The data set was then decay corrected, and 
processed by the 2-step resampling method described above, except that no 
transverse sensitivity correction was performed. The resampled data was 
then analyzed by totaling the counts in each projection image. FIG. 3 
shows that the number of counts in each projection is uniform, but there 
is still a small variation. In FIG. 4, the resampled counts from the study 
in FIG. 3 is plotted against the residual transverse sensitivity given in 
Eq. 1: 
##EQU18## 
where ro=5.5 cm, and R=34.5 cm. The normalization constant K is chosen so 
that the total counts over all views is the same. There is good visual 
agreement between the predicted and measured total counts. Statistical 
tests also support the null hypothesis. A chi-square test for the first 60 
projections (59 degrees of freedom) yielded a value of 4100, indicating 
significant differences. On the other hand, a Kolmogorov-Smirnov test 
indicated no significant difference at the 20% rejection level, with a 
t-value of 0.0048, and a t-critical of 0.098. A t-test (two-tailed, 
unequal variances, 119 degrees of freedom) yielded a similar result: There 
was no significant difference, with a t-value of 1.times.10.sup.-12 and a 
t-critical of 1.29 at the 20% rejection level. 
The post processing correction is given by: 
##EQU19## 
where u(x) is given by Equation 4. 
Although others have determined a 3-D FDR filter for diverging collimators, 
their method was designed to work with a direct method of reconstruction, 
as a prefiltering step. In the method of this invention, I will apply the 
3-D FDR to the fan beam data after resampling. Thus, the diverging beam 
resolution as a function of distance needs to be transformed to the 
resampled parallel beam geometry. The object to be reconstructed is: 
EQU O(x,z,.phi.)=o(s,t,z:.phi.). 
The unblurred projection f.sub.ll in parallel geometry is given by: 
##EQU20## 
The unblurred fan beam projection f.sub.fan is: 
##EQU21## 
where 
##EQU22## 
The blurred 3-D projection sinogram f.sub.1.sup.b (s,z,.phi.) for parallel 
geometry is: 
##EQU23## 
where c(s, t, z) is the blur function, and typically has a FWHM that 
varies linearly with the source-to-detector distance. With the 
substitution of Equation 17 into Equation 18, we obtain: 
##EQU24## 
But the factor (R-t)/(F+s.sub.2)=(F-s.sub.1)/(F+s.sub.2) is just the 
additional depth dependent factor in the fan beam transfer function. The 
2-D Fourier transform of Equation 28 is: 
##EQU25## 
The collimator resolution in parallel geometry transforms into fan beam 
geometry. We obtain the same Equation 29 if we use another form of 
Equation 28: 
##EQU26## 
the 2-D Fourier transform is unchanged, so Equation 28 and Equation 30 are 
equivalent. 
If b has a Fourier transform: 
##EQU27## 
This is the MTF of the fan beam collimator transfer function with the back 
plane distance B=O. Therefore the blur function to use with the resampled 
projection data is the parallel model given by Equation 27. 
Equation 28 also yields a proof of the distance dependent sensitivity of a 
diverging beam collimator. If o is a delta function centered at 
(s.sub.0,t.sub.0,z.sub.0) so that: 
EQU o(.omega..sub.1,t,.omega..sub.2 ;.phi.)=exp(-i2.pi.(.omega..sub.1 s.sub.0 
+.omega..sub.2 z.sub.0)).delta.(t-t.sub.0). (Eq. 33) 
then the sensitivity at a distance t.sub.o in the rotated frame is 
proportional to: 
##EQU28## 
The first fraction is due to the change of variables between l and t, and 
is almost constant for realistic collimators. The second fraction 
(F+s.sub.2)(R-t.sub.0) is the distance dependent sensitivity term given in 
Equations 21 and 22. 
To apply the F DP to the blurred projection data, first resample the data 
to obtain f.sub.2 (x,z,.phi.) and take the 3-D-FT of f.sub.2 : 
##EQU29## 
where t=-m/2.pi.v yields a stationary phase approximation to Equation 27. 
Then 
##EQU30## 
is used in a Wiener filter to obtain an estimate of .function..sub.1 
(m,.omega..sub.1,.omega..sub.2) from .function..sub.b 
(m,.omega..sub.1,.omega..sub.2). 
For the 3-D FDR, the blur function b is a gaussian, with: 
##EQU31## 
The Wiener filter is: 
EQU .function..sub.1 
(m,.omega..sub.1,.omega..sub.2)=W(m,.omega..sub.1,.omega..sub.2).function. 
.sub.b (m,.omega..sub.1,.omega..sub.2). (Eq. 36) 
with: 
##EQU32## 
and where the estimate of the projection power spectrum 
.vertline..function..sub.e (m,.omega..sub.1,.omega..sub.2).vertline..sup.2 
is obtained from: 
##EQU33## 
The function metz is a Metz filter with kernel 
##EQU34## 
and exponent X. The parameter X is set to 12. The quantity 
.vertline..pi..vertline..sup.2 is the estimate of the FT variance of the 
modified projection data. 
The same study that was used for the sensitivity study was also used for a 
study of resolution. The FWHM of the resampled point source images was 
measured as a function of distance from the detector, and compared to a 
parallel collimator model based on the same length (34.9 mm), width (1.4 
mm) sepal thickness (0.15 mm) and shape (regular hexagonal hole). The 
results are shown in FIG. 5. I utilized the collimator MTF function for a 
regular hexagon. I used the quadrature rule to combine the various 
contributions to the calculated FWHM. The FWHM of the 1 cm point source 
projection and the intrinsic resolution (3.5 mm) were combined with the 
FWHM of a point source to yield an estimate of the measured FWHM. A mean 
interaction depth of 8.9 mm was incorporated in the calculated FWHM. The 
utility function used to measure FWHM demonstrated a lot of variability, 
but the two curves are fairly close. Also, the resolution study with the 
cold rod phantom and the 3-D FDR resolution filter, discussed below, 
demonstrate that the result of Equation 32 is valid. 
For overall quantitative accuracy with uniform attenuation, the 
Liver-Spleen phantom for fan beam geometry, that consists of attenuating 
ellipses with varying amounts of activity within each ellipse, is shown. 
The concentrations were in the ratios 0:2:5:10. The quantitative accuracy 
using the resampling method has been preserved. A simulated fan beam 
acquisition with the following parameters was used: Focal length: 53.5 cm, 
focal point radium of rotation 34.82, matrix size 128, pixel size 0.585 
cm, and 256 projections in 360 degrees. The large pixel size was necessary 
to include the liver-spleen phantom in the field-of-view. Thus, the 
phantom dimensions are not realistic from a physical point of view. The 
large size of the phantom, however, means that it is more difficult to 
obtain an accurate attenuation correction. The large numbers of 
projections were necessary because the effective transverse pixel length 
was about 0.3 cm, and the liver-spleen phantom was 30 cm wide from right 
to left. The reconstruction of the liver-spleen phantom using the 
resampling algorithms is illustrated in Table 2. 
TABLE 2 
______________________________________ 
Region Simulated Activity 
Measured Activity 
______________________________________ 
Air 0 0.069 .+-. 0.005 
Background 2 2.01 .+-. 0.04 
Spleen 5 5.03 .+-. 0.045 
Liver 10 9.95 .+-. 0.08 
Lesion 0 -0.049 .+-. 0.095 
______________________________________ 
To determine how resampling algorithm affected resolution, a point source 
in air was scanned on the Picker 3000 equipped with UHRFAN collimators. 
The point source was placed 5.5 cm from the center of rotation, and the 
cameras were placed 19.3 cm from the center of rotation. About 300 .mu.Ci 
of Tc-99m was placed in a spherical point source 1 cm i.d., and the point 
source was placed on the headrest. The matrix size was 128, with pixel 
size 3.56 mm. One-hundred twenty projections were acquired in step and 
shoot mode. The projections were decay corrected and then reconstructed 
using the direct method, and again using the resampling algorithm with 
back projection. No low pass filtering was used. 
The FWHM and FWTM of the reconstructed point source images were then 
measured. 
Very little difference in the reconstructions was found, except that the 
3-D FDR method resulted in a noticeable resolution recovery. In Table 3, 
the FWMH and FWTM measurements from the two methods are shown. There is 
also very little difference in the measurements. The 3-D FDR 
reconstruction has minimal artifacts, and the best resolution of the three 
methods was tested. 
TABLE 3 
______________________________________ 
FWHM,mm FWTM,mm 
______________________________________ 
Direct method 11.03 20.42 
Resampling method 
11.28 20.28 
3-D FDR Filter 9.70 16.46 
______________________________________ 
To test the overall practicality and effectiveness of the combined 
protocol, a scan was made of a standard SPECT phantom with hot spheres, 
cold spheres, and cold rod inserts in a uniform background of about 1 
.mu.Ci/ml of Tc-99m. The sphere sizes were 31.8 mm, 25.4 mm, 19.1 mm, 15.9 
mm, 12.7 mm and 9.5 mm in diameter. The rod sizes were 12.7 mm, 11.1 mm, 
9.5 mm, 7.9 mm, 6.4 mm and 4.6 mm. The acquisition parameters were 120 
projections in 360 degrees, continuous mode, 128 matrix with a 3.56 mm 
pixel size, and 15s per stop. There were about 120,000 counts in each 
projection. 
The protocol consisted of the following steps. First, the raw data was 
decay corrected, scaled the output to block exponent form, and the 
estimate of the variance of the DFT was created. Next, the data was 
resampled using hybrid resampling method, and the FT variance was updated 
at each step. Then a 2-D stationary Wiener filter or the 3-D FDR Wiener 
filter was applied to the data set. A "noise multiplier" of 0.5 was used 
with both data sets. No attempt was made to find an optimal noise 
multiplier. Then the filtered projections were back projected using the 
CHT-Bellini method to obtain the final result. 
The conventional methodology consisted of direct fan beam back projection, 
followed by a 3-D Wiener filter for resolution recovery and scatter 
correction. The default value of 0.5 was used for the noise multiplier. 
The last step of the conventional methodology was attenuation compensation 
using the post attenuation correction (one iteration). 
For the hot/cold spheres study, the results of the conventional method, the 
3-D FDR/resampling method, and the 2-D prefilter/resampling method were 
about the same, with the 3-D FDR yielding perhaps the best qualitative 
result. The 2-D prefilter showed somewhat more noise amplification, while 
the conventional 3-D post filter overfiltered the cold spheres, and was 
not able to resolve the smallest hot sphere, as were the other two 
methods. 
The resampled fan beam reconstructions of the cold rod study show good 
resolutions, and are relatively free of artifacts. On the other hand, the 
conventional reconstruction was artifactual, and could not reliably detect 
even the large 12.7 mm rods. This protocol cannot be used to detect small 
cold focal lesions. The data set was very ill posed because of the low 
count level, a minimal number of projections and the continuous scanning. 
The 3-D FDR reconstruction shows the best resolution, with the 9.5 mm rods 
clearly resolved, and good uniformity, but with a slight ringing artifact 
at the center. This may be due to an inaccurate FDR filter at the low 
frequencies, where the stationary phase approximation is not very 
accurate. It is also characteristic of over correction for resolution, and 
is a form of aliasing. Aliasing is not well understood in resampling, 
because the frequency map is transformed by the resampling as well. This 
type of artifact is also characteristic of an inaccurate center of 
rotation correction, but because the 2-D prefiltered reconstruction did 
not show a similar result, this is not likely to be the source of error. 
The 3-D FDR reconstruction also showed a slight "scalloping" artifact at 
the periphery of the reconstruction. This was due to angular aliasing 
because of angular understampling. If an effective resolution of about 6 
mm is assumed for the 3-D FDR and the diameter of the FoV of 25 cm, the 
120 projections were below the number recommended to avoid angular 
aliasing 22, 23!. The 3-D FDR filter, because it results in additional 
resolution, will need to have M larger than with more conventional 
reconstruction methods. 
In summary, there has been shown that there is a depth dependent factor and 
transverse factor in the determination of fan beam sensitivity. 
Straightforward geometric arguments determine the nature of the 
depth-dependent fan beam sensitivity, and how it can be rectified with 
resampling to parallel geometry. The transverse component closely follows 
the predicted values, and is removed by a post processing procedure. The 
sensitivity point source study validated this result. The variation, 
however, is not great; it is less than 5% for a FoV 11 cm in diameter. 
The depth dependent resolution also depends on a similar geometric factor, 
and there is no transverse variation in resolution. It has been shown that 
the depth dependent fan beam resolution can be transformed to an 
equivalent parallel beam model with reasonable accuracy. 
Also, resampling allowed the use of existing and effective methodology for 
compensation of uniform attenuation: The contouring algorithm, the 
premultiply algorithm, and the CHT back projection algorithm for parallel 
geometry and uniform attenuation were utilized without modification. The 
3-D FDR algorithm and the 2-D prefilter had to be modified for anisotropic 
pixel sizes. 
An interesting observation is that frequency modulation is an ideal low 
pass filter with maximum bandwidth, whereas linear interpolation 
necessarily cannot be. 
The reconstructed point source study demonstrated that resolution is not 
much affected by the resampling algorithm. This is unexpected, because 
linear interpolation is still needed for this method. Perhaps the greatest 
errors in resampling methods based on 2-D interpolation came from the 
angular interpolation. The 3-D FDR filter gave the best improvement in 
resolution. 
The 3-D FDR method developed here is effective and is fairly efficient, the 
greatest numerical cost being the computation of 3-0D FFTs for large 
clinical data sets. In its present configuration, the small amount of 
processing time required will not preclude it from clinical application.