PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
Inverse Probl. Author manuscript; available in PMC 2010 April 30.
Published in final edited form as:
Inverse Probl. 2010 March 25; 26(4): 45008.
doi:  10.1088/0266-5611/26/4/045008
PMCID: PMC2861831
NIHMSID: NIHMS191512

Advancements to the planogram frequency–distance rebinning algorithm

Abstract

In this paper we consider the task of image reconstruction in positron emission tomography (PET) with the planogram frequency–distance rebinning (PFDR) algorithm. The PFDR algorithm is a rebinning algorithm for PET systems with panel detectors. The algorithm is derived in the planogram coordinate system which is a native data format for PET systems with panel detectors. A rebinning algorithm averages over the redundant four-dimensional set of PET data to produce a three-dimensional set of data. Images can be reconstructed from this rebinned three-dimensional set of data. This process enables one to reconstruct PET images more quickly than reconstructing directly from the four-dimensional PET data. The PFDR algorithm is an approximate rebinning algorithm. We show that implementing the PFDR algorithm followed by the (ramp) filtered backprojection (FBP) algorithm in linogram coordinates from multiple views reconstructs a filtered version of our image. We develop an explicit formula for this filter which can be used to achieve exact reconstruction by means of a modified FBP algorithm applied to the stack of rebinned linograms and can also be used to quantify the errors introduced by the PFDR algorithm. This filter is similar to the filter in the planogram filtered backprojection algorithm derived by Brasse et al. The planogram filtered backprojection and exact reconstruction with the PFDR algorithm require complete projections which can be completed with a reprojection algorithm. The PFDR algorithm is similar to the rebinning algorithm developed by Kao et al. By expressing the PFDR algorithm in detector coordinates, we provide a comparative analysis between the two algorithms. Numerical experiments using both simulated data and measured data from a positron emission mammography/tomography (PEM/PET) system are performed. Images are reconstructed by PFDR+FBP (PFDR followed by 2D FBP reconstruction), PFDRX (PFDR followed by the modified FBP algorithm for exact reconstruction) and planogram filtered backprojection image reconstruction algorithms. We show that the PFDRX algorithm produces images that are nearly as accurate as images reconstructed with the planogram filtered backprojection algorithm and more accurate than images reconstructed with the PFDR+FBP algorithm. Both the PFDR+FBP and PFDRX algorithms provide a dramatic improvement in computation time over the planogram filtered backprojection algorithm.

1. Introduction

Positron emission tomography (PET) is a medical imaging method that utilizes positron-emissing pharmaceuticals to assess the physiology and patho-physiology of patients by mapping the distribution of the tracer in their body (Wernick and Aarsvold 2004, Bendriem and Townsend 1998, Natterer 1986, Natterer and Wübbeling 2001). This task is accomplished by detecting the anti-colinear annihilation photons with specialized systems. The virtual line connecting a given pair of detector elements is called a line of response (LOR). The number of coincidences measured by a given pair of detector elements can be modeled by the Radon transform which is a line integral of the concentration of the radiotracer along the LOR determined by the given detector element pair.

This work is motivated by a new positron emission mammography/tomography (PEM/PET) system (Raylman et al 2006, 2007, 2008) developed at West Virginia University. The system comprises two pairs of parallel planar detectors. Only the panels parallel to each other are in coincidence. Each detector plane has an effective size of 2L = 195.3 mm by 2H = 144.9 mm which is composed of a 94 by 70 array of detector pixels with a distance of Td = 2.1 mm between the centers of two adjacent detector pixels. The distance between two opposing detector heads is 2R = 264 mm. The detector gantry is capable of rotation (in step-and-shoot mode) to acquire data at all azimuthal angles.

High quality, rapid image reconstruction algorithms for this system are desired for the detection and quantitation of small lesions (less than 5 mm). If a lesion is detected, an image-guided biopsy is performed while the patient remains on the imaging table. An efficient image reconstruction algorithm is desired to minimize the time the patient must remain on the imaging table.

Image reconstruction algorithms can be developed by analytic inversion of the Radon transform which is parameterized into sets of projections. A projection (figure 1) is the 2D set of LORs that intersect the image for a fixed direction (fixed polar and azimuthal angle).

Figure 1
Sketch of a (left) nontruncated projection and a (right) truncated projection.

Sufficient conditions for image reconstruction from the PET projection data have been characterized by Orlov (1976). Suppose we plot the angular orientation of all the measured projections on a unit sphere (known as the Orlov sphere) and denote this set by Ω. The Orlov condition states that the measured data are sufficient if there is no great circle on the unit sphere that does not intersect Ω. If a proper subset of Ω satisfies the Orlov condition then the problem is over determined, i.e. there is redundancy in the measurements. Figure 2 displays two Orlov sphere examples.

Figure 2
Examples of sets of projections (shown in black) plotted on the Orlov sphere. Both Ω0 and Ω1 represent sufficient data, but Ω1 ⊃ Ω0 contains additional redundant data.

Image noise is related to the number of measured coincident photon pairs. To optimize the signal-to-noise ratio (SNR), one should use as much of the data as possible. The parameterization of LORs in three-dimensional Euclidean space (i.e. the acquired data set) is a four-dimensional set. Image reconstruction in PET is an overdetermined problem because a four-dimensional data set is used to determine a three-dimensional object (the source of redundancy in the Orlov condition). We will refer to the data whose LORs have a co-polar angle of zero as the direct data and the remaining data as the oblique data. The Orlov sphere of the direct data is shown in the left plot of figure 2. One could reconstruct the image by only using the (three-dimensional) direct data instead of the entire (four-dimensional) set of data, but this reconstruction scheme results in a suboptimal SNR.

Many fully 3D reconstruction algorithms have been developed for PET imaging. These algorithms reconstruct the image directly from the entire 4D set of measured data. Examples of such algorithms are 3D filtered backprojection (Colsher 1980) and OSEM (Hundson and Larkin 1994). Brasse et al (2004) have developed a 3D filtered backprojection algorithm in the planogram coordinate system.

A rebinning algorithm averages over the oblique (redundant) data and outputs a three-dimensional direct data set. Standard 2D reconstruction algorithms can then be used to reconstruct the image, slice-by-slice, from this rebinned direct data. Image reconstruction with rebinning algorithms is more computationally efficient than exact analytic reconstruction methods such as 3D filtered backprojection. Another advantage of rebinning algorithms over 3D filtered backprojection algorithms is their ability to use data associated with truncated projections. A truncated projection is a projection in which some of the LORs in the projection are not measured (figure 1) due to the finite size of the detectors. Filtered backprojection algorithms must either throw out data associated with a truncated projection or include an extra step to compensate for the missing data such as the implementation of a reprojection algorithm (Kinahan and Rogers 1989).

The PFDR (Champley et al 2008) and Kao (Kao et al 2004) algorithms are the only rebinning algorithms that have been developed specifically for PET scanners with parallel planar detectors. Several rebinning algorithms have been developed for cylindrical PET scanners such as Fourier rebinning (FORE) (Defrise et al 1997), exact Fourier rebinning (FOREX) (Defrise et al 1997) and Fourier rebinning based on John’s equation (FORE-J) (Defrise and Liu 1999). Single slice rebinning (SSRB) (Daube-Witherspoon and Muehllehner 1987) is a general rebinning algorithm that can be applied to any PET geometry, but is less accurate than the PFDR and FORE algorithms.

The content of this paper is concerned with further advancements to the planogram frequency–distance rebinning (PFDR) algorithm (Champley et al 2008) which is an approximate rebinning algorithm derived in the planogram coordinate system. This rebinning algorithm is based on the frequency-distance relationship (Champley et al 2008, Bernardi et al 2007). The planogram coordinate system is a native coordinate system for PET scanners with parallel planar detectors. In the original PFDR paper (Champley et al 2008) the authors developed the PFDR algorithm and provided numerical experiments using both simulated data and data from the PEM/PET system described above. The numerical experiments provided a qualitative and quantitative comparison of images reconstructed with the 2D FBP algorithm from data rebinned with the PFDR and SSRB algorithms. Computation times for the PFDR algorithm were also stated.

This paper is organized as follows: in the next section we describe the planogram coordinate system which is the cross product of two linogram coordinate systems (Edholm and Herman 1987). Section 3 briefly describes the 2D filtered backprojection algorithm in linogram coordinates for variable number of views (rotations of the detector gantry). The main results of this paper are contained in section 4, where we show that if we reconstruct our image by implementing the PFDR algorithm followed by the linogram filtered backprojection algorithm from section 3, the result can be expressed as a filtered version of the true image. An explicit formula for this filter is derived. This filter can be used to appropriately modify the ramp filter in the linogram filtered backprojection algorithm to achieve exact reconstruction or explicitly characterize the errors in images reconstructed from the PFDR algorithm followed by filtered backprojection. This method for exact reconstruction requires complete projections which must be completed with a reprojection algorithm. We also show how the filter used in the Brasse algorithm (Brasse et al 2004) is related to the filter we derive for exact reconstruction. In section 5 we develop the PFDR algorithm in detector coordinates. Expressing the PFDR algorithm in detector coordinates also enables us to provide a comparative analysis of the PFDR and Kao algorithms. Numerical experiments using both simulated and measured data from the PEM/PET scanner are outlined in section 6. We finish with some concluding remarks in section 7.

1.1. Notation

In the following discussion, we will refer to the support of a function, h:RnR, as the set

supp(h){xRn:h(x)0}.

The inner product will be denoted by left angle bracket·, ·right angle bracket. The complement of a set A will be denoted by Ac. The indicator function of a set A will be denoted by

1A(x){1,xA,0,xA.}

Binary subscripts mark variables for which we have taken the Fourier transform. For example, if hL1(R2), then the Fourier transform of h in the second variable is

h^01(x,Y)F{h(x,)}(Y)Rh(x,y)e2πiyYdy.

This notation was introduced in Brasse et al (2004) and is necessary because in the following it will be difficult to identify the variables for which we have taken the Fourier transform based on the arguments.

We will use the following theorem throughout this paper. For a proof see Folland (1999).

Theorem 1.1

Let fL1(Rn) be real-valued and A be a nonsingular matrix and B = (AT)−1. Then

F{f(Ax)}(X)=detA1F{f(x)}(BX).

Moreover, if A is a rotation, then

F{f(Ax)}(X)=F{f(x)}(AX).

2. The planogram transform

Let f (x, y, z) denote the concentration of the PET radiotracer with

supp(f)Sf{(x,y,z):x2+y2<a2,z<c},

where a < min(R, L) and cH. In practice, it is often the case that c > H, but since none of the coincident photon pairs are measured for annihilations beyond the axial extent of the scanner, we may assume without loss of generality that cH.

The signals measured by the interaction of annihilation photons in the detector module are decoded by the electronics, and the estimated interaction position is given in detector (s, t) coordinates (figure 3). We will parameterize the PET data in planogram coordinates which can be expressed as a linear transformation of detector coordinates (figure 3) by

u0=12(sAsB),v0=12R(sA+sB)
(1)

u1=12(tA+tB),v1=12R(tAtB).
(2)
Figure 3
Parameterization of the LORs in the planogram coordinate system.

The PET data can be modeled by the planogram transform (Brasse et al 2004) given by

g(u0,v0,u1,v1)Pf(u0,v0,u1,v1)Rf(u0v0y,y,u1v1y)dy,
(3)

which is the 3D extension of the linogram transform:

g(u,v,z)Rf(uvy,y,z)dy

developed by Edholm (Edholm and Herman 1987). The variables (v0,v1) describe the angular orientation of the LOR, so g(·, v0, ·, v1) is a 2D parallel beam projection. Observe that the planogram transform (3) is a parameterization of line integrals of f weighted by

ddy[u0v0yyu1v1y]=[v01v1]=1+v02+v12.
(4)

Given supp(f) [subset, dbl equals] Sf, the support of g is contained in the butterfly-shaped set given by

supp(g)Sg{(u0,v0,u1,v1)R4:u0a1+v02,u1c+av1}.

Due to the finite size of the detectors, the planogram transform cannot be measured for all (u0,v0,u1,v1)R4. The measured domain is the diamond cross diamond shape given by

MML×MHML{(u0,v0)R2:u0LRv0}MH{(u1,v1)R2:u1HRv1}.

We define the measured planogram transform by

gm(u0,v0,u1,v1)={g(u0,v0,u1,v1),(u0,v0,u1,v1)M,0,(u0,v0,u1,v1)M.}

After observing the intersection of Sg and M, we see that there is no truncation in the u0 and u1 directions for

v0vm0RLaR2+L2a2R2a2,v1vm1HcR+a.
(5)

In other words for [mid ]v0[mid ]vm0 and [mid ]v1[mid ]vm1, the 2D projections given by gm(·, v0, ·, v1) are not truncated. Now we define the restricted domain by

Mm0{(u0,v0,u1,v1)M:v0vm0}

and the restricted planogram by

gm0(u0,v0,u1,v1)={g(u0,v0,u1,v1),(u0,v0,u1,v1)Mm0,0,(u0,v0,u1,v1)Mm0.}

2.1. Gantry rotations

Complete 3D tomographic reconstruction is not possible with one detector orientation of the PEM/PET scanner since we do not have projection data over all azimuthal angles (i.e. Orlov’s condition (Orlov 1976) is not met). To acquire full tomographic data, the PEM/PET scanner is rotated around the scanner axis (z-axis). If the detector gantry is rotated azimuthally counterclockwise by ψ, the planogram data are given by

gψ(u0,v0,u1,v1)Rf(u0cosψ(v0cosψ+sinψ)t,u0sinψ(v0sinψcosψ)t,u1v1t)dt.

A more convenient way to express gψ is

gψ(u0,v0,u1,v1)=Rfψ(u0v0y,y,u1v1y)dy,fψ(x,y,z)=f(xcosψysinψ,xsinψ+ycosψ,z).

We will use the convention g = gψ=0 and f = f−ψ=0. One can show that

g(u0,v0,u1,v1)=1cosψ+v0sinψ×gψ(u0cosψ+v0sinψ,v0cosψsinψcosψ+v0sinψ,u1u0v1sinψcosψ+v0sinψ,v1cosψ+v0sinψ).

From this last relation, we see that the sampling of the data is quite different for each view; the transformation of the sampling lattice is a nonlinear mapping.

Note that due to the dependence of vm0 on a, the number of detector positions necessary (for this particular geometry and rotation direction) for tomographic reconstruction (to satisfy the Orlov condition) depends on a. Specifically, for a given a, one needs

Nψ=Nψ(a)min{nN:2φm0nπ}

detector positions, where [var phi]m0 [equivalent] tan−1(vm0); or equivalently

Nψ=min{nN:vm0tan(π2n)}.

If we write the equation for vm0 in terms of a, we have

a=LRvm01+vm02

which is a decreasing function for 0 ≤ 0vmoLR, so

Nψ=min{nN:aLRtan(π2n)1+tan2(π2n)}=min{nN:aLcos(π2n)Rsin(π2n)}.
(6)

For example, using a = 60 mm, we will need three gantry positions (Nψ = 6 detector positions), spaced 180°6=30° apart to fulfill Orlov’s condition. Thus, in this case we will have Nψ = 6 sets of data:

{gψ:ψ=0°,30°,60°,90°,120°,150°}.

3. Linogram filtered backprojection algorithm from multiple views

After application of the PFDR algorithm, we are left with a stack of data in linogram coordinates for each detector position. In this section we explore image reconstruction with the linogram filtered backprojection algorithm from multiple views, i.e. data from multiple gantry rotations. For a detailed analysis of the specific case where only two detector positions are used, see Magnusson (1993). The rotational linogram transform is given by

gψ(u,v,z)Rf(ucosψ(vcosψ+sinψ)t,usinψ(vsinψcosψ)t,z)dt=Rfψ(uvy,y,z)dy.

We now state the projection-slice theorem for the general linogram transform.

Proposition 3.1

Let fL1(R3) and gψ be defined as above. Then we have that

g^ψ,100(U,v,z)=f^110(U(cosψvsinψ),U(sinψ+vcosψ),z).

Proof

For ψ = 0, we have g^100(U,v,z)=f^110(U,vU,z), see Magnusson (1993). For ψ ≠ 0, the result follows by application of theorem 1.1 to the last relation.

The image, f, can be recovered from one detector orientation if the detector is of infinite length, as suggested by the following theorem.

Proposition 3.2

Let fL1(R3) and gψ be given as above. Then

f(x,y,z)=RRUg^ψ,100(U,v,z)e2πi[xcosψ+ysinψv(xsinψycosψ)]UdUdv.

Proof

We have

f(x,y,z)=RRf^110(X,Y,z)e2πi[xX+yY]dXdY=002πrf^110(rcosφ,rsinφ,z)e2πir[xcosφ+ysinφ]drdφ=R0πrf^110(rcosφ,rsinφ,z)e2πir[xcosφ+ysinφ]drdφ.

Now let

r[cosφsinφ]=U[cosψvsinψsinψ+vsinψ].
(7)

Then [mid ]r[mid ] dr d[var phi] = [mid ]U[mid ] dU dv and by taking the dot product of (7) with [xy], we have

r[xcosφ+ysinφ]=U[xcosψ+ysinψv(xsinψycosψ)].

Therefore,

f(x,y,z)=RRf^110(U(cosψvsinψ),U(sinψ+vcosψ),z)×e2πi[xcosψ+ysinψv(xsinψycosψ)]UUdUdv=RRUg^ψ,100(U,v,z)e2πi[xcosψ+ysinψv(xsinψycosψ))]UdUdv.

The next theorem generalizes the above result to the case of detectors of finite length and is the basis for multiview linogram filtered backprojection algorithm.

Proposition 3.3

Let fL1(R3) and supp(f) [subset, dbl equals] Sf. Also let Nψ defined as in (6) and ψk=πNψk for k = 0, 1, …, Nψ − 1. Then

f(x,y,z)=k=0Nψ1vm0vm01μ(v)RUg^ψk,100m0(U,v,z)×e2πi[xcosψk+ysinψkv(xsinψkycosψk)]UdUdv,

where the sum is over all the gantry rotations and

μ(v)=k=0Nψ11{v:vcosψksinψkcosψk+vsinψkvm0}(v)

is the number of gantry rotations that have measured LORs in this particular direction.

The proof follows easily from theorem 3.2.

4. Exact reconstruction and characterization of error of the PFDR

The PFDR algorithm (Champley et al 2008) provides a method to reconstruct an approximation of the distribution of the radiotracer, f, with any 2D reconstruction algorithm. In this section we provide additional analysis of the PFDR algorithm including a method for exact reconstruction and an error analysis of reconstructing images with the PFDR algorithm followed by the filtered backprojection algorithm outlined in the previous section.

For notational convenience, we define a linear operator K:L1(R4)L1(R3) by

(Kh)^110(U0,V0,z;v1max)v1maxv1maxh^1100(U0,V0,zV0U0v1,v1)dv1(Kh)^110(U0,V0,z)limv1max(Kh)^110(U0,V0,z;v1max).

For fL1(R3) and supp(f) [subset, dbl equals] Sf, the PFDR algorithm rebins the planogram data into a stack of linograms according to

g^110reb(U0,V0,z;v1max)1N(V0U0,z;v1max)(Kgm)^110(U0,V0,z;v1max)N(V0U0,z;v1max)v1maxv1max1MSgc(0,0,zV0U0v1,v1)dv1.

In the above we have the PFDR algorithm applied to the measured data, but the same procedure can be applied to the restricted data. The parameter v1max controls the maximum co-polar acceptance angle to be rebinned. The term N(−V0/U0,z; v1max) is used to normalize the varying number of contributions from each of the oblique projections resulting from the finite size of the detectors. Note that if v1maxvm1, then N(−V0/U0,z; v1max) = 2v1max and in this case PFDR algorithm is given by

g^110reb(U0,V0,z;v1max)12v1max(Kgm)^110(U0,V0,z;v1max)=12v1maxv1maxv1maxg^1100m(U0,V0,zV0U0v1,v1)dv1.

The next theorem establishes a connection with the K operator and the backprojection operator. The backprojection operator is the (formal) adjoint of the planogram transform which is given by

Ph(x,y,z)RRh(x+v0y,v0,z+v1y,v1)dv0dv1.

Theorem 4.1

Let hL1(R4). Then

(Ph)^100(U0,y,z)=(Kh)^110(U0,yU0,z).
(8)

Proof

We have

(Kh)^110(U0,yU0,z)Rh^1100(U0,yU0,z+v1y,v1)dv1=RRRh(u0,v0,z+yv1,v1)e2πi[v0yu0]U0du0d0dv1=RRRh(x+v0y,v0,z+yv1,v1)e2πixU0dxdv0dv1=(Ph)^100(U0,y,z).

Thus, the K operator computes the backprojection over v1 and the substitution V0 = −yU0 followed by an inverse Fourier transform in U0 variable computes the backprojection over v0. The rebinned measured data with v1max = H/R are related to the backprojection operation of the measured data by the relation

g^reb(U0,yU0,z;HR)=1N(y,z;HR)(Kgm)^110(U0,yU0,z;HR)=1N(y,z;HR)(Pgm)^110(U0,y,z).

As a numerical algorithm, the K operator applied to measured data is actually (Kgm)^110(U0,V0,z) and not (Kgm)^110(U0,yU0,z). The distinction is subtle, but in the discrete case the calculation of (Kgm)^(U0,yU0,z) requires 1D interpolations in Fourier space. This can be implemented accurately and efficiently with the CHIRPZ algorithm (Magnusson 1993). Calculating the backprojection over v0 by numerical integration gives slightly better results (data for comparison not shown).

The next theorem is the basis for the PFDR algorithm for restricted data and will be used in the subsequent analysis.

Theorem 4.2

Let fL1(R3), supp(f) [subset, dbl equals] Sf, and let gm0 be the restricted planogram as defined above. Then

g^ψ,1100m0(U0,V0,u1,v1)=1U0RRe2πi[(V0U0)Y+(u1+v1(V0U0))Z]×f^ψ,111(U0,Y,Z)1{(U0,Y,Z):Yv1Zvm0U0}(U0,Y,Z)dYdZ×1MSgc(0,0,u1,v1)
(9)

=f^ψ,100(U0,y,u1v1y)sin(2πvm0U0y)πU0yy=V0U0×1MSgc(0,0,u1,v1).
(10)

For a proof see Champley et al (2008).

A rebinning algorithm such as PFDR enables one to (approximately) reconstruct an image slice-by-slice with a 2D reconstruction algorithm. The next theorem derives the result of reconstructing our image from the rebinned data with the 2D filtered backprojection algorithm (theorem 3.3). We will refer to the discrete implementation of this theorem as the PFDR+FBP algorithm.

Theorem 4.3

Let fL1(R3) and supp(f) [subset, dbl equals] Sf with a=Lcos(π2Nψ)Rsin(π2Nψ). Let

[xψkyψk][cosψksinψksinψkcosψk][xy],

where ψk=πNψk for k = 0, 1, …, Nψ − 1. Then {gψkm0:k=0,1,,Nψ1} contains a partition of the set of projections over all azimuthal angles. Define

fPFDR(x,y,z)k=0Nψ1vm0vm0RU0N(yψk,z;v1max)(Kgψkm0)^100(U0,v0,z;v1max)×e2πi[xψk+v0yψk]U0dU0dv0,

where gm0 is defined as in section 1. Then

fPFDR(x,y,z)=R3dU0dYdZf^111(U0,Y,Z)e2πi[xU0+yY+zZ]×[k=0Nψ11N(yψk,z;v1max)v1maxv1max1MSgc(0,0,z+yψkv1,v1)]×[1{(U0,Y,Z):Yv1Zvm0U0}(U0cosψk+Ysinψk,YcosψkU0sinψk,Z)dv1].

Proof

Let

Aψ=[cosψsinψ0sinψcosψ0001].

Then with the help of theorem 4.2 we have

vm0vm0RU0N(yψ,z;v1max)(Kgψm0)^100(U0,v0,z;v1max)e2πi[xψ+v0yψ]U0dU0dv0=RU0N(yψ,z;v1max)(Kgψm0)^110(U0,yψU0,z;v1max)e2πixψU0dU0=RU0N(yψ,z;v1max)v1maxv1maxg^ψ,1100m0(U0,yψU0,z+yψv1,v1)e2πixψU0dv1dU0=1N(yψ,z;v1max)v1maxv1maxdv11MSgc(0,0,z+yψkv1,v1)×R3f^ψ,111(U0,Y,Z)1{(U0,Y,Z):Yv1Zvm0U0}(U0,Y,Z)e2πi[xψU0+yψY+zZ]dU0dYdZ=1N(yψ,z;v1max)v1maxv1maxdv11MSgc(0,0,z+yψkv1,v1)×R3f^ψ,111(U0,Y,Z)1{(U0,Y,Z):Yv1Zvm0U0}(U0,Y,Z)×e2πiAψT(x,y,z)T,(U0,Y,Z)TdU0dYdZ=1N(yψ,z;v1max)v1maxv1maxdv11MSgc(0,0,z+yψkv1,v1)×R3f^111(Aψ(U0,Y,Z)T)1{(U0,Y,Z):Yv1Zvm0U0}(U0,Y,Z)×e2πi(x,y,z)T,AψT(U0,Y,Z)TdU0dYdZ=1N(yψ,z;v1max)v1maxv1maxdv11MSgc(0,0,z+yψkv1,v1)×R3f^111(U0,Y,Z)1{(U0,Y,Z):Yv1Zvm0U0}(AψT(U0,Y,Z)T)e2πi[xU0+yY+zZ]dU0dYdZ

The last equality follows from theorem 1.1.

Thus, the PFDR+FBP algorithm reconstructs a filtered version of the true object, where the (shift-variant) 3D filter is given in frequency space by

k=0Nψ11N(yψk,z;v1max)v1maxv1max1{(U0,Y,Z):Yv1Zvm0U0}(AψT(U0,Y,Z)T)×1MSgc(0,0,z+yψkv1,v1)dv1.
(11)

When we restrict the integration of v1 to [mid ]v1[mid ]v1maxvm1, we have N(y,z;v1max)=12v1max and this filter is independent of yψk = y cos ψkx sin ψk and z (i.e. the filter is stationary) and is equal to

12v1maxk=0Nψ1h^PFDR(AψkT(U0,Y,Z)T;v1max),
(12)

where

h^PFDR(U0,Y,Z;v1max)v1maxv1max1{(U0,Y,Z):Yv1Zvm0U0}(U0,Y,Z).

This filter is illustrated in figure 4 and explicitly derived in proposition 4.6. The filter is equal to one except inside a cone aligned in the axial direction. The aperture of this cone is given by the parameter v1max.

Figure 4
Axial slices of the filter 12v1maxk=1Nψ1h^PFDR(AψkT(U0,Y,Z)T;v1max) given by (12) with v1max = 0.5.

When v1max → 0, fPFDR = f as shown in theorem 3.3.

The following theorem can be used to implement an efficient reprojection algorithm in the planogram coordinate system (Brasse et al 2004) by employing only FFT operations, similar to the 3D Fourier reprojection (3D-FRP) algorithm (Matej and Lewitt 2001) derived in sinogram coordinates. We will refer to this method of reprojection as the planogram Fourier reprojection (PFRP) algorithm. The theorem is also used in the proof of theorem 4.5.

Theorem 4.4

Let fL1(R3), supp(f) [subset, dbl equals] Sf, and let gm0 be the restricted planogram as defined above. Then for [mid ]v1[mid ] ≤ vm1 we have

g^ψ,1010m0(U0,v0,Z,v1)=f^ψ,111(U0,v0U0+v1Z,Z)1{v0:v0vm0}(v0).

For a proof see Brasse et al (2004).

The next theorem establishes a connection between the direct data g^1010m0(U0,v0,Z,0) and the rebinned date (Kgm0)^101(U0,v0,Z;v1max).

Theorem 4.5

Let fL1(R3) with supp(f) [subset, dbl equals] Sf and v1max ≤ vm1 be given. Define

h^PFDR(U0,Y,Z;v1max)v1maxv1max1{(U0,Y,Z):Yv1Zvm0U0}(U0,Y,Z)dv1

just as we have above. Then for [mid ]v0[mid ] ≤ vm0 we have that

(Kgψm0)^101(U0,v0,Z;v1max)=f^ψ,111(U0,v0U0,Z)h^PFDR(U0,v0U0,Z;v1max)
(13)

=g^ψ,1010m0(U0,v0,Z,0)h^PFDR(U0,v0U0,Z;v1max).
(14)

Proof

Since v1maxvm1, 1MSgc(0,0,u1,v1)=1 for [mid ]v1v1max. Now we use theorem 4.2 with the substitutions V0 = −yU0 and u1 = z + v1y to establish the relation

g^ψ,1100m0(U0,yU0,z+v1y,v1)=1U0R2f^ψ,111(U0,Y,Z)×1{(U0,Y,Z):Yv1Zvm0U0}e2πi(yY+zZ)dYdZ

and thus

(Kgψm0)^110(U0,yU0,z;v1max)=v1maxv1maxg^ψ,1100m0(U0,yU0,z+yv1,v1)dv1=1U0R2f^ψ(U0,Y,Z)h^PDFR(U0,Y,Z;v1max)e2πi(yY+zZ)dYdZ.

By taking the Fourier transform of the variables (y, z) of both sides of the above equality we have that

1U0f^ψ,111(U0.Y,Z)h^PFDR(U0,YmZ;v1max)=R(Kgψm0)^111(U0,yU0,Z;v1max)e2πiyYdy=1U0R(Kgψm0)^111(U0,V0,Z;v1max)e2πiV0(YU0)dV0=1U0(Kgψm0)^101(U0,YU0,Z;v1max).

Making the substitution v0 = Y/U0 we have

(Kgψm0)^101(U0,v0,Z;v1max)=f^ψ,111(U0,v0U0,Z)h^PFDR(U0,v0U0,Z;v1max),=g^ψ,1010m0(U0,v0,Z,0)h^PFDR(U0,v0U0,Z;v1max).

The filter h^PFDR(X,Y,Z;v1max) is explicitly derived in the next proposition.

Proposition 4.6

Define

h^PFDR(X,Y,Z;v1max)v1maxv1max1{(X,Y,Z):Yv1Zvm0X}(X,Y,Z)dv1.

Then

h^PFDR(X,Y,Z;v1max)=1{(X,Y,Z):Yvm0Xv1maxz}(X,Y,Z)×{2v1max,Z=0,2v1max,Yvm0Xv1maxZ,vm0XY+v1maxZZ,Yvm0Xv1maxZY,2vm0XZ,vm0Xv1maxZY.}

Furthermore, when [mid ]v0[mid ] ≤ vm0

U0h^PFDR(U0,v0U0,Z;v1max)={U02v1max,v1maxZ(vm0v0)U0,U0Zv1maxZ+(vm0v0)U0,v1maxZvm0U0v0U0,12vm0Z,(vm0+v0)U0v1maxZ,}=max(U02v1max,U0Zv1maxZ+(vm0v0)U0,12vm0Z)

is a continuous function on {(U0,v0,Z)R3:v0vm0}.

The proof of the above proposition is straight-forward and thus will be omitted for the sake of brevity. Now we are ready to state our theorem for exact reconstruction with the PFDR algorithm which is the main result of this paper. We will refer to the discrete implementation of this theorem as the PFDRX algorithm.

Theorem 4.7

Under the hypotheses of theorem 4.5 we have that

f(x,y,z)=k=0Nψ1vm0vm01μ(v0)R2U0h^PFDR(U0,v0U0,Z;v1max)(Kgψkm0)^101(U0,v0,Z;v1max)×e2πi[(xψk+v0yψk)U0+zZ]dU0dZdv0,

where

μ(v)=k=0Nψ11{v:vcosψksinψkcosψk+vsinψkvm0}(v).

Proof

The proof immediately follows from theorems 4.5 and 3.3.

Thus one can obtain exact reconstruction by the use of the PFDR algorithm followed by a modified linogram filtered backprojection, where the usual ramp filter, [mid ]U0[mid ], is replaced by

U0h^PFDR(U0,v0U0,Z;v1max).

Theorem 4.7 enables one to obtain an exact reconstruction from the set of nontruncated projections. Often the support of f extends to the ends of our axial field of view. In this case vm1 = 0. The PFDRX algorithm along with 3D FBP algorithms may use a reprojection algorithm to complete the unmeasured projections (Kinahan and Rogers 1989) to extend vm1. A reprojection algorithm takes an initial estimate of the image and re-projects the data onto virtual detectors to estimate the missing data. Brasse et al (2004) suggest using theorem 4.4 as a means of efficient reprojection. Reprojection could also be used to extend vm0.

4.1. Connection to the planogram filtered backprojection algorithm developed by Brasse et al

There is a lot of similarity between the PFDRX algorithm and the planogram filtered backprojection algorithm developed by Brasse et al (2004). The PFDRX algorithm reconstructs the exact image by first rebinning the algorithm into a stack of linograms (which is equivalent to a backprojection over v1), filtering the stack of rebinned linograms, and then backprojecting over v0. The planogram filtered backprojection algorithm filters the data and then performs 2D backprojection over both v0 and v1. The filter used by the planogram filtered backprojection algorithm is given by

1k=0Nψ1p^(Aψk(U0,v0U0+v1U1,U1)T),

where

1p^(U0,v0U0+v1Z,Z)={U02v1max,U0vm0Zv1maxv0U0+v1Z,U0Zv1maxZ+vm0U0v0U0+v1Z,v1maxZvm0U0v0U0+v1Z,Z2vm0,U0vm0Zv1maxv0U0+v1Z}

Note that for {(U0,v0,Z)R3:v0vm0}

U0h^PFDR(U0,v0U0,Z;v1max)=1p^(U0,v0U0+v1Z,Z)v1=0.

The PFDRX algorithm is computationally more efficient than the planogram filtered backprojection algorithm by implementing a more efficient method of backprojecting over v1 and filtering a three-dimensional set instead of a four-dimensional set of data. The next section outlines a method to further improve the computational complexity of the PFDR algorithm where applied to the measured data.

5. PFDR algorithm in detector coordinates

In this section we derive the PFDR algorithm in detector coordinates to improve the computational efficiency. The relations between the detector coordinates and planogram coordinates are given by equations (1) and (2) and drawn in figure 3. Define the PET data transform in detector coordinates by

q(sA,sB,tA,tB)Rf(12(sAsB)12R(sA+sB)y,y,12(tA+tB)12R(tAtB)y)dy.

Similar to the planogram transform, we define the measured and restricted data in detector coordinates by

qm(sA,sB,tA,tB)q(sA,sB,tA,tB)1{(sA,sB,tA,tB)[L,L]2×[H,H]2}(sA,sB,tA,tB),qm0(sA,sB,tA,tB)qm(sA,sB,tA,tB)1{(sA,sB)R2:12R(sA+sB)vm0}(sA,sB),

respectively.

Now let

A0[1R1R],B0(A0T)1=[1212R1212R]A1[1R1R].

Then we have that

A0[u0v0]=[sAsB]andA1[u1v1]=[tAtB]

for the detector coordinates (sA,sB,tA,tB) and the planogram coordinates (u0,v0,u1,v1). The planogram coordinate and detector coordinate transforms are related by

q(A0(u0,v0)T,A1(u1,v1)T)=g(u0,v0,u1,v1).

Then by theorem 1.1 we have that

g^1100(U0,V0,u1,v1)=F{g(,,u1,v1)}(U0,V0)=F{q(A0(,)T,A1(u1,v1)T)}(U0,V0)=detA01F{q(,,A1(u1,v1)T)}(B0(U0,V0)T)=12Rq^1100(B0(U0,V0)T,A1(u1,v1)T).
(15)

The basis for the PFDR algorithm lies in the relation

g^1100(U0,V0,u1,0)=g^1100(U0,V0,u1v1V0U0,v1),

which holds for v1R. By the above relationship and (15), we have that

q^1100(B0[U0V0],A1[u10])=q^1100(B0[U0V0],A1[u1v1V0U0v1]).
(16)

Now let (XA,XB) be such that

B0[U0V0]=[XAXB]

and thus

[U0V0]=B01[XAXB]=A0T[XAXB]=[11RR][XAXB]=[XAXBR(XA+XB)].

Then from (16) we arrive at the PFDR relation in detector coordinates

q^1100(XA,XB,u1,u1)=q^1100(XA,XB,u1XB(2Rv1XAXB),u1XA(2Rv1XAXB)).
(17)

Thus, the PFDR algorithm in detector coordinates is given by

q110reb(XA,XB,z)1N(RXA+XBXAXB,z;v1max)×v1maxv1maxq^1100m(XA,XB,zXB(2Rv1XAXB),zXA(2Rv1XAXB))dv1.
(18)

Suppose that we wish to perform the PFDR algorithm on the measured data set and we have Ns detector pixels in the transaxial direction and Nt in the axial direction. The transformation from detector coordinates (sA,sB) to linogram coordinates (u0,v0) transforms the data set from a square block to a diamond-shaped block (ML) with zeros placed within the data set in a checkerboard pattern as shown in figure 5. In other words the sampling in linogram coordinates is on an interlaced lattice. Thus, in detector coordinates we have a data set of Ns2 elements and in linogram coordinates we have a data set of size (2Ns − 1)2. For fixed (tA,tB) an FFT of the data set q(·,·,tA,t) requires Ns2log(Ns2) operations and an FFT of the data set g(·, ·,u1,v1) requires (2Ns − 1)2log((2Ns − 1)2) operations. Computation of the FFT in detector coordinates provides a computational complexity gain by a factor of about 4. Also note that in the absence of zero padding in the FFT, the PFDR algorithm must rebin (2Ns − 1)2 different frequencies in planogram coordinates and only Ns2 in detector coordinates. Therefore, we expect that implementing the PFDR algorithm applied to the measured data set in detector coordinates requires about four times fewer operations than the PFDR algorithm in planogram coordinates.

Figure 5
Detector coordinate samples (left). Linogram samples (right). The ×s mark the samples and the dots mark zeros in the data. Here we have shown 36 samples which would require a data size of 11 × 12 = 121 in linogram coordinates.

5.1. Comparison of PFDR and the Kao rebinning algorithm

By letting z = u1 and t=2Rv1XAXB in equation (17) we have

q^1100(XA,XB,z,z)=q^1100(XA,XB,zXBt,zXAt).

This relation, which we call the PFDR relation, is only an equality in the case of ideal data (Champley et al 2008). In the case of measured or restricted data this relation is merely approximate. The Kao rebinning algorithm (Kao et al 2004) is described in our notation by

q^1100m(XA,XB,z,z)=q^1100m(XA,XB,zXBt,zXAt)0tddtq^1100m(XA,XB,zXBt,zXAt)dt,

Thus, one can perform exact rebinning on the entire set of measured data with the Kao algorithm. Unfortunately in practice this algorithm (Kao et al 2004) distorts the data and increases noise levels. This increase in noise level is likely due to the use of finite difference methods to approximate derivatives in the correction term.

With nontruncated projections, i.e. with [mid ]v0[mid ]vm0 and [mid ]vm1[mid ]vm1, one can perform exact reconstruction with the PFDRX algorithm by theorem 4.7. However, we must note that this does not imply that we have developed a method for exact rebinning. Consider equation (14). Since h^PFDR(U0,v0U0,Z)U0=0, we may not filter the rebinned data (Kgm0)^101(U0,v0,Z;vm1). In theorem 4.7, exact reconstruction is performed by filtering the rebinned data by

U0h^PFDR(U0,v0U0,Z;v1max).

The [mid ]U0[mid ] term prevents this filter from blowing up near U0 = 0 and thus provides a stable filter.

6. Numerical experiments

In this section we describe the numerical experiments that were performed to test the PFDR+FBP and PFDRX algorithms. We start by discussing the implementation details of the PFRP, PFDR+FBP, PFDRX and Brasse algorithms. The next two subsections discuss the experiments performed on simulated and measured data, respectively.

6.1. Implementation details

The PFDR algorithm was implemented in detector coordinates (18). The FFTW library (Frigo and Johnson 2005) was used to compute the FFTs. The accuracy of the PFDR algorithm depends on the amount of zero padding used. For our experiments, we used FFTs of size 144 × 144 for the 94 × 94 array of data. By definition, the PFDRX algorithm can only be applied to the restricted data set, gm0 for [mid ]v1[mid ]vm1. In the PFDR+FBP algorithm one can rebin the measured data set, gm, or just the restricted data set, gm0, but in either case the FBP algorithm only uses the projections such that [mid ]v0[mid ]vm0. We found that using only the restricted data set in the rebinning step of the PFDR+FBP algorithm resulted in images in lower contrast as compared to images where the measured data set was rebinned. For this paper we only evaluated the PFDR+FBP algorithm applied to the measured data set.

The restricted data set for vm0 = tan(15°) represents about 59% of the data. This loss of data is certain to result in a reduction in the SNR, but for parallax considerations (Wernick and Aarsvold 2004) in the PEM/PET scanner this data is discarded anyway (Raylman et al 2008).

The planogram Fourier reprojection (PFRP) algorithm was implemented as follows. First we calculated an initial estimate of our image. One can use only the direct data for this initial estimate, but we found that better results could be achieved by reconstructing the initial estimate with the PFDR+FBP algorithm with a small polar acceptance angle of 5°. This method seemed to provide a nice compromise between accuracy and noise. For each data set, gψ, we then computed fψ. Then we calculated samples of f^ψ,101(X,y,Z) with the FFT algorithm. Now for each v1 and z we calculated the product

e2πiyv1Zf^ψ,101(X,y,Z).

Then for each v0 we used the CHIRPZ (Magnusson 1993) algorithm to calculate a Fourier transform in the y dimension to output frequency samples spaced by Tv0X units, where Tv0 is the sampling distance in the v0 variable. Thus, we have calculated

f^ψ,101(X,v0X+v1Z,Z)=g^ψ(X,v0,Z,v1).

To prevent aliasing, one should multiply g^ψ(X,v0,Z,v1) by

1{YR:Y12Ty}(v0X+v1Z),

where 12Ty is the Nyquist frequency of our image and Ty is the sampling distance in the y dimension. Finally we calculate an inverse 2D FFT algorithm to get gψ(u0,v0,u1,v1) which is our reprojected data.

To minimize ringing artifacts introduced by Gibb’s phenomenon, the filters were modified slightly. For the PFDRX and Brasse algorithms, we filtered the rebinned data by

U0h^PFDR(U0,v0U0,Z;v1max)(1(2Tu0U0)8)(1(2Tu1Z)8),1k=0Nψ1p^(Aψ(U0,v0U0+v1U1,U1)T)(1(2Tu0U0)8)(1(2Tu1Z)8),

respectively, and for the PFDR+FBP, we filtered the rebinned data by

U0(1(2Tu0U0)8),

where Tu0 and Tu1 are the sampling rates for u0 and u1, respectively. The terms (1−(2Tu0U0)8) and (1 − (2Tu1Z)8) are used to apodize the digital filters and were optimized empirically.

The axial support of the rebinned data extends past the axial field of view. The PFDRX algorithm filters along this dimension, so to avoid ringing artifacts caused by Gibb’s phenomenon, we calculated g110reb(U0,V0,z;v1max) for all [mid ]z[mid ]c + av1max. This maximum value for z was chosen by empirical observations.

From theorem 4.4, we have

R2gm0(u0,v0,u1,v1)du0du1=R3f(x,y,z)dxdydz

for [mid ]v0[mid ]vm0 and [mid ]v1[mid ]vm1. This relationship is known as the zeroth-order Radon consistency condition (Natterer 1986), and we use this to enforce the total number of counts in the reconstructed image to be the same as the number of counts averaged over all the 2D projections.

To account for the effects of detector sensitivity due to the solid angle of a detector pair as seen from the point at the center of the scanner, we multiplied the data by

(1+v02+v12)32Tv0Tv1.

Note that

(1+v02+v12)32dv0dv1=cosθdθdφ,

where θ is the co-polar angle and [var phi] in the azimuthal angle. Along with the planogram correction (reciprocal of (4)), the data were weighted by an overall factor of

(1+v02+v12)12(1+v02+v12)32=1+v02+v12.

The reconstruction software was written in C++ and run on a 3.06 GHz Intel Core 2 Duo processor. The images were all reconstructed from six sets of data: {gψ : ψ = 0°, 30°, 60°, 90°, 120°, 150°}. Since the PFDR+FBP, PFDRX and Brasse reconstruction algorithms process each set of data independently, parallelizing the reconstruction algorithms is a trivial matter. One can parallelize the data onto up to six processors to speed the computation. We parallelized the data onto the two processors on our machine. Image reconstruction times for each of the algorithms (with v1max = tan(15°)) we tested are shown in table 1 for images of voxel side length 2.1/2 mm (image is 141 × 115 × 115) and table 2 for images of voxel side length 2.1/4 mm (image is 281 × 233 × 233). All images shown below have voxels of side length 2.1/4 mm. When the reconstructions were serialized onto one processor, the computation times were almost exactly twice as long. Thus, we expect the reconstruction time to be reduced to a third of what was reported in tables tables11 and and22 when the reconstruction is parallelized onto six processors.

Table 1
Image reconstruction computation times (s) with v1max =tan(15°) and voxel side length of d/2, where d = 2.1 mm is the detector sampling distance.
Table 2
Image reconstruction computation times (s) with v1max =tan(15°) and voxel side length of d/4, where d = 2.1 mm is the detector sampling distance.

6.2. Simulated data

Simulated data enables controlled experiments with known ground-truth to provide quantitative analysis (accuracy and noise properties) of image reconstruction algorithms. Digital phantoms were reconstructed with the PFDR+FBP, PFDRX and Brasse algorithms. We implemented analytic ray-tracing algorithms to simulate data from the PEM/PET scanner whose specifications were described in the introduction. Each LOR was weighted by the solid angle that the detector crystal pair subtends to provide a rough estimate of the detector sensitivity. Physical effects such as detector response, scatter, attenuation, etc were not incorporated in our simulation. Where noted, Poisson noise was applied to the simulated data.

Our phantom (figure 6) consists of a warm background with hot and cold cylindrical regions (contrast between hot cylinders and background is 2:1 and cold regions have no activity). The small cylinders have a diameter and height of 2.4 mm, and the large cylinders have a diameter of 9.6 mm and a height of 3.6 mm. The phantom extends to the top edge of the axial field of view. This is similar to breast imaging where the breast extends to the top edge of the axial field of view where it meets the chest wall.

Figure 6
Phantom used for numerical experiments. Axial slice at z = 2.4 (left). Sagittal slice at x =−36.686 (center). Coronal slice at y = 0 (right). The black lines profiles in figures figures88 and and1010.

The support of our phantom is contained within a cylinder of radius 60 mm. Thus, from equation (5) we have that vm0 = 0.2691 ≈ tan(15°). To fulfill the Orlov conditions, we need Nψ = 6 detector positions (see equation (6)) at rotation angles ψ = 0°, 30°, 60°, 90°, 120°, 150°. Since the phantom extends to ends of the axial extent of the scanner, vm1 = 0.

Figure 7 displays noiseless reconstructions with v1max = tan(15°), and figure 8 displays plots of profiles through these images. Figure 9 displays noiseless reconstructions with v1max=HRtan(29°) and figure 10 displays plots of profiles through these images. Note that the maximum oblique LOR measured has a value of v1=HR. Reconstructions with noisy data are shown in figure 13. Approximately 777 million coincidences were used in the reconstruction algorithms. This count level is higher than what one should expect from the PEM/PET scanner, but was chosen to illustrate the differences between the three reconstruction algorithms.

Figure 7
Reconstructions using the PFDR+FBP, PFDRX and Brasse algorithms with v1max = tan(15°) from simulated noiseless data. Axial slices are on the top row, and coronal slices are on the bottom row.
Figure 8
Profiles of images reconstructed from the PFDR+FBP, PFDRX and Brasse algorithms with v1max = tan(15°) from simulated noiseless data. Figure 6 shows the location of these profiles.
Figure 9
Reconstructions using the PFDR+FBP, PFDRX and Brasse algorithms with v1max = tan(29°) from simulated noiseless data. Axial slices are on the top row, and coronal slices are on the bottom row.
Figure 10
Profiles of images reconstructed from the PFDR+FBP, PFDRX and Brasse algorithms with v1max = tan(29°) from simulated noiseless data. Figure 6 shows the location of these profiles.
Figure 13
Reconstructions using the PFDR+FBP, PFDRX and Brasse algorithms with v1max = tan(15°) from simulated noisy data. Axial slices are on the top row and coronal slices are on the bottom row.

The root mean squared error (RMSE) of the large hot and cold cylinders reconstructed with the noiseless data is plotted in figures figures1111 and and1212 for v1max = tan(15°) and v1max = tan(29°), respectively.

Figure 11
Root mean squared error (RMSE) for large hot (left) and cold (right) disks of images reconstructed from simulated noiseless data with v1max = tan(15°).
Figure 12
Root mean squared error (RMSE) for large hot (left) and cold (right) disks of images reconstructed from simulated noiseless data with v1max = tan(29°).

The relative L2 error of images reconstructed from the noisy data is shown in table 3. The relative L2 error is defined by f~f2f2, where f is the true image and f~ is the reconstructed image.

Table 3
Relative L2 error of images reconstructed from simulated noisy data with v1max = tan(15°).

Noise properties were evaluated in the images reconstructed from noisy data. The signal-to-noise ratio (SNR) was measured by calculating the ratio of the mean to the standard deviation in regions of interest in the warm background of the phantom. We chose regions of interest that were localized within one axial slice of the image to evaluate how the SNR varies axially. The SNR is plotted versus the axial (z) position of the region of interest in figure 14.

Figure 14
Signal-to-noise and contrast recovery coefficient measurements of images reconstructed from simulated noisy data with v1max = tan(15°). Signal-to-noise ratio (SNR) in the background versus depth (left). The contrast recovery coefficient (CRC) ...

The contrast recovery coefficients (CRC) for the 2.4 mm hot cylinders of images reconstructed from simulated noisy data are plotted in figure 14.

The root mean squared error (RMSE) of the large hot and cold cylinders reconstructed from noisy data is plotted in figure 15.

Figure 15
Root mean squared error (RMSE) for large hot (left) and cold (right) disks of images reconstructed from simulated noisy data with v1max = tan(15°).

6.3. PEM/PET data

We collected data from our PEM/PET scanner (Raylman et al 2006, 2007, 2008) to demonstrate the quality of images reconstructed with the PFDR+FBP and PFDRX algorithms from real data. The micro-Derenzo cold rod phantom was placed on its side (long axis of the cold rods aligned in the transaxial direction) and scanned with the PEM/PET scanner. A total of about 66 million counts were collected and images were reconstructed with v1max = tan(15°). These images are shown in figure 16.

Figure 16
The micro-Derenzo cold rod phantom was placed on its side in the PEM/PET scanner to provide challenging axial variation for a rebinning algorithm. This phantom has cylindrical rods of diameters 1.2, 1.6, 2.4, 3.2, 4.0 and 4.8 mm with no activity. Images ...

7. Discussion and conclusion

We have developed a theoretically exact and numerically stable method for analytic image reconstruction with the PFDR algorithm that is more computationally efficient than filtered backprojection. We refer to this as the PFDRX algorithm (theorem 4.7). After one completes the 2D projection data by a reprojection algorithm, the PFDRX algorithm can be implemented in three steps. First one applies the PFDR algorithm to the four-dimensional set of restricted data to rebin the data into a three dimensional set that consists of a stack of linogram data. Next one applies a 2D filter to the rebinned data. Finally one back-projects the rebinned data over v0 to form the image.

In this paper we have also developed more supporting theory to the PFDR rebinning algorithm. Theorem 4.3 characterizes the error in images reconstructed by first rebinning the data with the PFDR algorithm and then with the 2D filtered backprojection in linogram coordinates. We refer to this reconstruction procedure as the PFDR+FBP algorithm. In section 5 we derived the PFDR algorithm in detector coordinates. We expect that implementing the PFDR algorithm in detector coordinates requires about four times fewer operations than the PFDR algorithm in planogram coordinates when used to rebin the measured data.

This work is the first to present the implementation of the Fourier reprojection algorithm in planogram coordinates based on theorem 4.4. This method was first suggested as a means of reprojection by Brasse et al (2004). Similar efficient reprojection algorithms in sinogram coordinates can be implemented with 3D-FRP (Matej and Lewitt 2001) or FOREPROJ (Liu et al 1999).

A comparative analysis of the PFDR+FBP, PFDRX and Brasse (filtered backprojection algorithm in planogram coordinates) algorithms has been evaluated with simulated noise-free and noisy data. From figures figures77--10,10, we see that in the case of noise-free data that the Brasse algorithm seems to provide the most qualitatively accurate images, followed by the PFDRX algorithm and then the PFDR+FBP algorithm. The three algorithms rank in the opposite order when it comes to computation speed as shown in tables tables11 and and2.2. The three algorithms provide a tradeoff between accuracy and noise as illustrated in figure 17. The PFDR+FBP and PFDRX algorithms provide a dramatic reduction in computation time (tables (tables11 and and2)2) for a small reduction in accuracy (table 3 and figure 15) versus the Brasse algorithm.

Figure 17
Tradeoff between accuracy and speed of computation between the PFDR+FBP, PFDRX and Brasse algorithms.

Figures Figures77--1010 also illustrate how the accuracy of the PFDR+FBP algorithm seems to degrade with increased axial acceptance angle, while the accuracy of the PFDRX and Brasse algorithms remain essentially the same when applied to the noise-free data.

Reconstructions with simulated noisy data are shown in figure 13. One can still note the distortions present in the images reconstructed with the PFDR+FBP algorithm. Differences in noise texture between the three algorithms can also be seen at the edges of the field of view along the axial direction. The apparent reduced noise level at the edges of the field of view in the PFDRX and Brasse algorithms is likely due to the use of data estimation via reprojection. Overall, the noise levels seem to be the largest in the images reconstructed with the PFDRX algorithm. This apparent increase in noise is verified quantitatively in figure 14. The slight increase in noise is likely due to the fact that one must correct for approximations made in the rebinning step by implementing a sharper filter.

Table 3 provides a quantitative measure of the global error of the PFDR+FBP, PFDRX and Brasse algorithms applied to noisy data. These measurements indicate that the PFDRX algorithm provides more quantitative accuracy than the PFDR+FBP algorithm, but less accuracy than the Brasse algorithm. To measure the quantitative accuracy in specific regions of interest, we plotted the RMSE error in the large hot and cold cylinders in figure 15. These metrics further support the assertion of the order of the rank in accuracy of the three algorithms tested.

Significant differences cannot be seen in the measurements of the contrast recovery coefficients of the small hot cylinders among the three reconstruction algorithms. Although these plots do indicate that the Brasse algorithm provides the best contrast and the PFDR+FBP provides the worst contrast, for a fixed z, the data differ by less than 10%.

Reconstructions of PEM/PET measured data agree with our findings with simulated data. The images reconstructed from the PEM/PET data with the PFDRX algorithm displays better contrast (more accuracy) than the images reconstructed with the PFDR+FBP algorithm (figure 16).

Acknowledgments

This work is sponsored by the National Institutes of Health under grants CA74135 and CA94196. The authors would like to thank Larry Pierce for his valuable comments and David Brasse for providing us with a copy of his planogram reconstruction code.

References

  • Bendriem B, Townsend DW. The Theory and Practice of 3D PET. Springer; Berlin: 1998.
  • Bernardi ED, Mazzoli M, Zito F, Baselli G. Evaluation of frequency–distance relation validity for FORE optimization in 3-D PET. IEEE Trans. Nucl. Sci. 2007;54:1670–8.
  • Brasse D, Kinahan P, Clackdoyle R, Comtat C, Defrise M, Townsend D. Fast fully 3D image reconstruction in PET using planograms. IEEE Trans. Med. Imaging. 2004;23:413–25. [PubMed]
  • Champley KM, Defrise M, Clackdoyle R, Raylman RR, Kinahan PE. Planogram rebinning with the frequency–distance relationship. IEEE Trans. Med. Imaging. 2008;27:925–33. [PMC free article] [PubMed]
  • Colsher J. Fully three-dimensional positron emission tomography. Phys. Med. Biol. 1980;25:103–15. [PubMed]
  • Daube-Witherspoon M, Muehllehner G. Treatment of axial data in three-dimensional PET. J. Nucl. Med. 1987;28:1717–24. [PubMed]
  • Defrise M, Kinahan P, Townsend D, Michel C, Sibomana M, Newport D. Exact and approximate rebinning algorithms for 3-D PET data. IEEE Trans. Med. Imaging. 1997;16:145–58. [PubMed]
  • Defrise M, Liu X. A fast rebinning algorithm for 3-D positron emission tomography using John’s equation. Inverse Problems. 1999;15:1047–65.
  • Edholm PR, Herman GT. Linograms in image reconstruction from projections. IEEE Trans. Med. Imaging. 1987;6:301–7. [PubMed]
  • Folland GB. Real Analysis. Wiley; New York: 1999.
  • Frigo M, Johnson SG. The design and implementation of FFTW3. Proc. IEEE. 2005;93:216–31. (special issue on ‘Program Generation, Optimization and Platform Adaptation’)
  • Hundson H, Larkin R. Accelerated image reconstruction using ordered subsets of projection data. IEEE Trans. Med. Imaging. 1994;13:601–9. [PubMed]
  • Kao CM, Pan X, Chen CT. An exact fourier rebinning algorithm for 3D PET imaging using panel detectors. Phys. Med. Biol. 2004;49:2407–23. [PubMed]
  • Kinahan P, Rogers JG. Analytic 3D image reconstruction using all detected events. IEEE Trans. Nucl. Sci. 1989;36:964–8.
  • Liu X, Defrise M, Michel C, Sibomana M, Comtat C, Kinahan P, Townsend D. Exact rebinning methods for three-dimensional pet. IEEE Trans. Med. Imaging. 1999;18:657–64. [PubMed]
  • Magnusson M. PhD dissertation. Department of Electrical Engineering, Linkping University; Linkping, Sweden: 1993. Linogram and other direct methods for tomographic reconstruction.
  • Matej S, Lewitt RM. 3D-FRP: direct fourier reconstruction with fourier reprojection for fully 3D PET. IEEE Trans. Nucl. Sci. 2001;48:1378–85.
  • Natterer F. The Mathematics of Computerized Tomography. Wiley; New York: 1986.
  • Natterer F, Wübbeling F. Mathematical Methods in Image Reconstruction. Society for Industrial & Applied Mathematics; Philadelphia: 2001.
  • Orlov SS. Theory of three-dimensional image reconstruction: I. Conditions for a complete set of projections. Sov. Phys.—Crystallogr. 1976;20:429–33.
  • Raylman R, Majewski S, Kross B, Popov V, Proffitt J, Smith M, Weisenberger A, Wojcik R. Development of a dedicated positron emission tomography system for the detection and biopsy of breast cancer. Nucl. Instrum. Methods Phys. Res. A. 2006;569:291–5.
  • Raylman R, et al. Initial testing of a positron emission mammography/tomography (PEM/PET) breast imaging and biopsy device. J. Nucl. Med. 2007;48:436.
  • Raylman R, et al. The positron emission mammography/tomography breast imaging and biopsy system (PEM/PET): design, construction and phantom-based measurements. Phys. Med. Biol. 2008;53:637–53. [PubMed]
  • Wernick MN, Aarsvold JN. Emission Tomography. Academic; New York: 2004.