Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC2861831

Formats

Article sections

- Abstract
- 1. Introduction
- 2. The planogram transform
- 3. Linogram filtered backprojection algorithm from multiple views
- 4. Exact reconstruction and characterization of error of the PFDR
- 5. PFDR algorithm in detector coordinates
- 6. Numerical experiments
- 7. Discussion and conclusion
- References

Authors

Related links

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/045008PMCID: PMC2861831

NIHMSID: NIHMS191512

Email: ude.notgnihsaw.u@klpmahc

See other articles in PMC that cite the published article.

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.

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 2*L* = 195.3 mm by 2*H* = 144.9 mm which is composed of a 94 by 70 array of detector pixels with a distance of *T _{d}* = 2.1 mm between the centers of two adjacent detector pixels. The distance between two opposing detector heads is 2

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).

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.

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.

In the following discussion, we will refer to the support of a function, $h:{\mathbb{R}}^{n}\to \mathbb{R}$, as the set

$$\mathrm{supp}\left(h\right)\equiv \{x\in {\mathbb{R}}^{n}\phantom{\rule{thinmathspace}{0ex}}:h(x)\ne 0\}.$$

The inner product will be denoted by ·, ·. The complement of a set **A** will be denoted by **A**^{c}. The indicator function of a set **A** will be denoted by

$${1}_{\mathbf{A}}\left(x\right)\equiv \{\begin{array}{cc}1,\hfill & x\in \mathbf{A},\hfill \\ 0,\hfill & x\notin \mathbf{A}.\hfill \end{array}\phantom{\}}$$

Binary subscripts mark variables for which we have taken the Fourier transform. For example, if $h\in {L}^{1}\left({\mathbb{R}}^{2}\right)$, then the Fourier transform of *h* in the second variable is

$${\widehat{h}}_{01}(x,Y)\equiv \mathcal{F}\left\{h\right(x,\cdot \left)\right\}\left(Y\right)\equiv {\int}_{\mathbb{R}}h(x,y){\mathrm{e}}^{-2\pi \mathrm{i}yY}\mathrm{d}y.$$

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).

Let $f\in {L}^{1}\left({\mathbb{R}}^{n}\right)$ be real-valued and A be a nonsingular matrix and B = (A^{T})^{−1}. Then

$$\mathcal{F}\left\{f\right(Ax\left)\right\}\left(X\right)={\mid \mathrm{det}\phantom{\rule{thinmathspace}{0ex}}A\mid}^{-1}\mathcal{F}\left\{f\right(x\left)\right\}\left(BX\right).$$

Moreover, if A is a rotation, then

$$\mathcal{F}\left\{f\right(Ax\left)\right\}\left(X\right)=\mathcal{F}\left\{f\right(x\left)\right\}\left(AX\right).$$

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

$$\mathrm{supp}\left(f\right)\subseteq {\mathbf{S}}_{f}\equiv \left\{\right(x,y,z):{x}^{2}+{y}^{2}<{a}^{2},\phantom{\rule{thinmathspace}{0ex}}\mid z\mid <c\},$$

where *a* < min*(R, L)* and *c* ≤ *H*. 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 *c* ≤ *H*.

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

$${u}_{0}=\frac{1}{2}({s}_{A}-{s}_{B}),\phantom{\rule{1em}{0ex}}{v}_{0}=\frac{1}{2R}({s}_{A}+{s}_{B})$$

(1)

$${u}_{1}=\frac{1}{2}({t}_{A}+{t}_{B}),\phantom{\rule{1em}{0ex}}{v}_{1}=\frac{1}{2R}({t}_{A}-{t}_{B}).$$

(2)

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

$$g({u}_{0},{v}_{0},{u}_{1},{v}_{1})\equiv Pf({u}_{0},{v}_{0},{u}_{1},{v}_{1})\equiv {\int}_{\mathbb{R}}f({u}_{0}-{v}_{0}y,y,{u}_{1}-{v}_{1}y)\mathrm{d}y,$$

(3)

which is the 3D extension of the linogram transform:

$$g(u,v,z)\equiv {\int}_{\mathbb{R}}f(u-vy,y,z)\mathrm{d}y$$

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

$$\Vert \frac{\mathrm{d}}{\mathrm{d}y}\left[\begin{array}{c}{u}_{0}-{v}_{0}y\hfill \\ y\hfill \\ {u}_{1}-{v}_{1}y\hfill \end{array}\right]\Vert =\Vert \left[\begin{array}{c}-{v}_{0}\hfill \\ 1\hfill \\ -{v}_{1}\hfill \end{array}\right]\Vert =\sqrt{1+{v}_{0}^{2}+{v}_{1}^{2}}.$$

(4)

Given supp(*f*) **S**_{f}, the support of *g* is contained in the butterfly-shaped set given by

$$\mathrm{supp}\left(g\right)\subseteq {\mathbf{S}}_{g}\equiv \left\{\right({u}_{0},{v}_{0},{u}_{1},{v}_{1})\in {\mathbb{R}}^{4}:\mid {u}_{0}\mid \le a\sqrt{1+{v}_{0}^{2}},\phantom{\rule{thinmathspace}{0ex}}\mid {u}_{1}\mid \le c+a\mid {v}_{1}\mid \}.$$

Due to the finite size of the detectors, the planogram transform cannot be measured for all $({u}_{0},{v}_{0},{u}_{1},{v}_{1})\in {\mathbb{R}}^{4}$. The measured domain is the diamond cross diamond shape given by

$$\begin{array}{cc}\hfill & \mathbf{M}\equiv {\mathbf{M}}_{L}\times {\mathbf{M}}_{H}\hfill \\ \hfill & {\mathbf{M}}_{L}\equiv \left\{\right({u}_{0},{v}_{0})\in {\mathbb{R}}^{2}:\mid {u}_{0}\mid \le L-R\mid {v}_{0}\mid \}\hfill \\ \hfill & {\mathbf{M}}_{H}\equiv \left\{\right({u}_{1},{v}_{1})\in {\mathbb{R}}^{2}:\mid {u}_{1}\mid \le H-R\mid {v}_{1}\mid \}.\hfill \end{array}$$

We define the measured planogram transform by

$${g}^{m}({u}_{0},{v}_{0},{u}_{1},{v}_{1})=\{\begin{array}{cc}g({u}_{0},{v}_{0},{u}_{1},{v}_{1}),\hfill & ({u}_{0},{v}_{0},{u}_{1},{v}_{1})\in \mathbf{M},\hfill \\ 0,\hfill & ({u}_{0},{v}_{0},{u}_{1},{v}_{1})\notin \mathbf{M}.\hfill \end{array}\phantom{\}}$$

After observing the intersection of **S**_{g} and **M**, we see that there is no truncation in the *u*_{0} and *u*_{1} directions for

$$\mid {v}_{0}\mid \le {v}_{m0}\equiv \frac{RL-a\sqrt{{R}^{2}+{L}^{2}-{a}^{2}}}{{R}^{2}-{a}^{2}},\phantom{\rule{1em}{0ex}}\mid {v}_{1}\mid \le {v}_{m1}\equiv \frac{H-c}{R+a}.$$

(5)

In other words for *v*_{0} ≤ *v*_{m0} and *v*_{1} ≤ *v*_{m1}, the 2D projections given by *g ^{m}*(·,

$${\mathbf{M}}_{m0}\equiv \left\{\right({u}_{0},{v}_{0},{u}_{1},{v}_{1})\in \mathbf{M}:\mid {v}_{0}\mid \le {v}_{m0}\}$$

and the restricted planogram by

$${g}^{m0}({u}_{0},{v}_{0},{u}_{1},{v}_{1})=\{\begin{array}{cc}g({u}_{0},{v}_{0},{u}_{1},{v}_{1}),\hfill & ({u}_{0},{v}_{0},{u}_{1},{v}_{1})\in {\mathbf{M}}_{m0},\hfill \\ 0,\hfill & ({u}_{0},{v}_{0},{u}_{1},{v}_{1})\notin {\mathbf{M}}_{m0}.\hfill \end{array}\phantom{\}}$$

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}_{\psi}({u}_{0},{v}_{0},{u}_{1},{v}_{1})\equiv {\int}_{\mathbb{R}}f({u}_{0}\phantom{\rule{thinmathspace}{0ex}}\mathrm{cos}\phantom{\rule{thinmathspace}{0ex}}\psi -({v}_{0}\phantom{\rule{thinmathspace}{0ex}}\mathrm{cos}\phantom{\rule{thinmathspace}{0ex}}\psi +\mathrm{sin}\phantom{\rule{thinmathspace}{0ex}}\psi )t,{u}_{0}\phantom{\rule{thinmathspace}{0ex}}\mathrm{sin}\phantom{\rule{thinmathspace}{0ex}}\psi -({v}_{0}\phantom{\rule{thinmathspace}{0ex}}\mathrm{sin}\phantom{\rule{thinmathspace}{0ex}}\psi -\mathrm{cos}\phantom{\rule{thinmathspace}{0ex}}\psi )t,{u}_{1}-{v}_{1}t)\mathrm{d}t.$$

A more convenient way to express *g _{ψ}* is

$$\begin{array}{cc}\hfill & {g}_{\psi}({u}_{0},{v}_{0},{u}_{1},{v}_{1})={\int}_{\mathbb{R}}{f}_{-\psi}({u}_{0}-{v}_{0}y,y,{u}_{1}-{v}_{1}y)\mathrm{d}y,\hfill \\ \hfill & {f}_{-\psi}(x,y,z)=f(x\phantom{\rule{thinmathspace}{0ex}}\mathrm{cos}\phantom{\rule{thinmathspace}{0ex}}\psi -y\phantom{\rule{thinmathspace}{0ex}}\mathrm{sin}\phantom{\rule{thinmathspace}{0ex}}\psi ,x\phantom{\rule{thinmathspace}{0ex}}\mathrm{sin}\phantom{\rule{thinmathspace}{0ex}}\psi +y\phantom{\rule{thinmathspace}{0ex}}\mathrm{cos}\phantom{\rule{thinmathspace}{0ex}}\psi ,z).\hfill \end{array}$$

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

$$\begin{array}{cc}\hfill g& ({u}_{0},{v}_{0},{u}_{1},{v}_{1})=\frac{1}{\mathrm{cos}\phantom{\rule{thinmathspace}{0ex}}\psi +{v}_{0}\phantom{\rule{thinmathspace}{0ex}}\mathrm{sin}\phantom{\rule{thinmathspace}{0ex}}\psi}\hfill \\ \hfill & \times {g}_{\psi}\left(\frac{{u}_{0}}{\mathrm{cos}\phantom{\rule{thinmathspace}{0ex}}\psi +{v}_{0}\phantom{\rule{thinmathspace}{0ex}}\mathrm{sin}\phantom{\rule{thinmathspace}{0ex}}\psi},\frac{{v}_{0}\phantom{\rule{thinmathspace}{0ex}}\mathrm{cos}\phantom{\rule{thinmathspace}{0ex}}\psi -\mathrm{sin}\phantom{\rule{thinmathspace}{0ex}}\psi}{\mathrm{cos}\phantom{\rule{thinmathspace}{0ex}}\psi +{v}_{0}\phantom{\rule{thinmathspace}{0ex}}\mathrm{sin}\phantom{\rule{thinmathspace}{0ex}}\psi},{u}_{1}-\frac{{u}_{0}{v}_{1}\phantom{\rule{thinmathspace}{0ex}}\mathrm{sin}\phantom{\rule{thinmathspace}{0ex}}\psi}{\mathrm{cos}\phantom{\rule{thinmathspace}{0ex}}\psi +{v}_{0}\phantom{\rule{thinmathspace}{0ex}}\mathrm{sin}\phantom{\rule{thinmathspace}{0ex}}\psi},\frac{{v}_{1}}{\mathrm{cos}\phantom{\rule{thinmathspace}{0ex}}\psi +{v}_{0}\phantom{\rule{thinmathspace}{0ex}}\mathrm{sin}\phantom{\rule{thinmathspace}{0ex}}\psi}\right).\hfill \end{array}$$

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 *v*_{m0} 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}_{\psi}={N}_{\psi}\left(a\right)\equiv \mathrm{min}\{n\in \mathbb{N}:2{\phi}_{m0}n\ge \pi \}$$

detector positions, where _{m0} tan^{−1}(*v*_{m0}); or equivalently

$${N}_{\psi}=\mathrm{min}\left\{n\in \mathbb{N}:{v}_{m0}\ge \mathrm{tan}\left(\frac{\pi}{2n}\right)\right\}.$$

If we write the equation for *v*_{m0} in terms of *a*, we have

$$a=\frac{L-R{v}_{m0}}{\sqrt{1+{v}_{m0}^{2}}}$$

which is a decreasing function for 0 ≤ $0\le {v}_{mo}\le {\scriptstyle \frac{L}{R}}$, so

$$\begin{array}{cc}\hfill {N}_{\psi}& =\mathrm{min}\left\{n\in \mathbb{N}:a\le \frac{L-R\phantom{\rule{thinmathspace}{0ex}}\mathrm{tan}\left(\frac{\pi}{2n}\right)}{\sqrt{1+{\mathrm{tan}}^{2}\left(\frac{\pi}{2n}\right)}}\right\}\hfill \\ \hfill & =\mathrm{min}\left\{n\in \mathbb{N}:a\le L\phantom{\rule{thinmathspace}{0ex}}\mathrm{cos}\left(\frac{\pi}{2n}\right)-R\phantom{\rule{thinmathspace}{0ex}}\mathrm{sin}\left(\frac{\pi}{2n}\right)\right\}.\hfill \end{array}$$

(6)

For example, using *a* = 60 mm, we will need three gantry positions (*N _{ψ}* = 6 detector positions), spaced ${\scriptstyle \frac{180\xb0}{6}}=30\xb0$ apart to fulfill Orlov’s condition. Thus, in this case we will have

$$\{{g}_{\psi}:\psi =0\xb0,30\xb0,60\xb0,90\xb0,120\xb0,150\xb0\}.$$

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

$$\begin{array}{cc}\hfill {g}_{\psi}(u,v,z)& \equiv {\int}_{\mathbb{R}}f(u\phantom{\rule{thinmathspace}{0ex}}\mathrm{cos}\phantom{\rule{thinmathspace}{0ex}}\psi -(v\phantom{\rule{thinmathspace}{0ex}}\mathrm{cos}\phantom{\rule{thinmathspace}{0ex}}\psi +\mathrm{sin}\phantom{\rule{thinmathspace}{0ex}}\psi )t,u\phantom{\rule{thinmathspace}{0ex}}\mathrm{sin}\phantom{\rule{thinmathspace}{0ex}}\psi -(v\phantom{\rule{thinmathspace}{0ex}}\mathrm{sin}\phantom{\rule{thinmathspace}{0ex}}\psi -\mathrm{cos}\phantom{\rule{thinmathspace}{0ex}}\psi )t,z)\mathrm{d}t\hfill \\ \hfill & ={\int}_{\mathbb{R}}{f}_{-\psi}(u-vy,y,z)\mathrm{d}y.\hfill \end{array}$$

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

Let $f\in {L}^{1}\left({\mathbb{R}}^{3}\right)$ and g_{ψ} be defined as above. Then we have that

$${\widehat{g}}_{\psi ,100}(U,v,z)={\widehat{f}}_{110}\left(U\right(\mathrm{cos}\phantom{\rule{thinmathspace}{0ex}}\psi -v\phantom{\rule{thinmathspace}{0ex}}\mathrm{sin}\phantom{\rule{thinmathspace}{0ex}}\psi ),U(\mathrm{sin}\phantom{\rule{thinmathspace}{0ex}}\psi +v\phantom{\rule{thinmathspace}{0ex}}\mathrm{cos}\phantom{\rule{thinmathspace}{0ex}}\psi ),z).$$

For *ψ* = 0, we have ${\widehat{g}}_{100}(U,v,z)={\widehat{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.

Let $f\in {L}^{1}\left({\mathbb{R}}^{3}\right)$ and g_{ψ} be given as above. Then

$$f(x,y,z)={\int}_{\mathbb{R}}{\int}_{\mathbb{R}}\mid U\mid {\widehat{g}}_{\psi ,100}(U,v,z){\mathrm{e}}^{2\pi \mathrm{i}[x\phantom{\rule{thinmathspace}{0ex}}\mathrm{cos}\phantom{\rule{thinmathspace}{0ex}}\psi +y\phantom{\rule{thinmathspace}{0ex}}\mathrm{sin}\phantom{\rule{thinmathspace}{0ex}}\psi -v(x\phantom{\rule{thinmathspace}{0ex}}\mathrm{sin}\phantom{\rule{thinmathspace}{0ex}}\psi -y\phantom{\rule{thinmathspace}{0ex}}\mathrm{cos}\phantom{\rule{thinmathspace}{0ex}}\psi \left)\right]U}\mathrm{d}U\phantom{\rule{thinmathspace}{0ex}}\mathrm{d}v.$$

We have

$$\begin{array}{cc}\hfill f(x,y,z)& ={\int}_{\mathbb{R}}{\int}_{\mathbb{R}}{\widehat{f}}_{110}(X,Y,z)\phantom{\rule{thinmathspace}{0ex}}{\mathrm{e}}^{2\pi \mathrm{i}[xX+yY]}\mathrm{d}X\phantom{\rule{thinmathspace}{0ex}}\mathrm{d}Y\hfill \\ \hfill & ={\int}_{0}^{\infty}{\int}_{0}^{2\pi}r{\widehat{f}}_{110}(r\phantom{\rule{thinmathspace}{0ex}}\mathrm{cos}\phantom{\rule{thinmathspace}{0ex}}\phi ,r\phantom{\rule{thinmathspace}{0ex}}\mathrm{sin}\phantom{\rule{thinmathspace}{0ex}}\phi ,z)\phantom{\rule{thinmathspace}{0ex}}{\mathrm{e}}^{2\pi \mathrm{i}r[x\phantom{\rule{thinmathspace}{0ex}}\mathrm{cos}\phantom{\rule{thinmathspace}{0ex}}\phi +y\phantom{\rule{thinmathspace}{0ex}}\mathrm{sin}\phantom{\rule{thinmathspace}{0ex}}\phi ]}\mathrm{d}r\phantom{\rule{thinmathspace}{0ex}}\mathrm{d}\phi \hfill \\ \hfill & ={\int}_{\mathbb{R}}{\int}_{0}^{\pi}\mid r\mid {\widehat{f}}_{110}(r\phantom{\rule{thinmathspace}{0ex}}\mathrm{cos}\phantom{\rule{thinmathspace}{0ex}}\phi ,r\phantom{\rule{thinmathspace}{0ex}}\mathrm{sin}\phantom{\rule{thinmathspace}{0ex}}\phi ,z)\phantom{\rule{thinmathspace}{0ex}}{\mathrm{e}}^{2\pi \mathrm{i}r[x\phantom{\rule{thinmathspace}{0ex}}\mathrm{cos}\phantom{\rule{thinmathspace}{0ex}}\phi +y\phantom{\rule{thinmathspace}{0ex}}\mathrm{sin}\phantom{\rule{thinmathspace}{0ex}}\phi ]}\mathrm{d}r\phantom{\rule{thinmathspace}{0ex}}\mathrm{d}\phi .\hfill \end{array}$$

Now let

$$r\left[\begin{array}{c}\hfill \mathrm{cos}\phantom{\rule{thinmathspace}{0ex}}\phi \hfill \\ \hfill \mathrm{sin}\phantom{\rule{thinmathspace}{0ex}}\phi \hfill \end{array}\right]=U\left[\begin{array}{c}\hfill \mathrm{cos}\phantom{\rule{thinmathspace}{0ex}}\psi -v\mathrm{sin}\phantom{\rule{thinmathspace}{0ex}}\psi \hfill \\ \hfill \mathrm{sin}\phantom{\rule{thinmathspace}{0ex}}\psi +v\mathrm{sin}\phantom{\rule{thinmathspace}{0ex}}\psi \hfill \end{array}\right].$$

(7)

Then *r* d*r* d = *U* d*U* d*v* and by taking the dot product of (7) with $\left[{\scriptstyle \begin{array}{c}\hfill x\hfill \\ \hfill y\hfill \end{array}}\right]$, we have

$$r[x\phantom{\rule{thinmathspace}{0ex}}\mathrm{cos}\phantom{\rule{thinmathspace}{0ex}}\phi +y\phantom{\rule{thinmathspace}{0ex}}\mathrm{sin}\phantom{\rule{thinmathspace}{0ex}}\phi ]=U[x\phantom{\rule{thinmathspace}{0ex}}\mathrm{cos}\phantom{\rule{thinmathspace}{0ex}}\psi +y\phantom{\rule{thinmathspace}{0ex}}\mathrm{sin}\phantom{\rule{thinmathspace}{0ex}}\psi -v(x\phantom{\rule{thinmathspace}{0ex}}\mathrm{sin}\phantom{\rule{thinmathspace}{0ex}}\psi -y\phantom{\rule{thinmathspace}{0ex}}\mathrm{cos}\phantom{\rule{thinmathspace}{0ex}}\psi \left)\right].$$

Therefore,

$$\begin{array}{cc}\hfill f(x,y,z)=& {\int}_{\mathbb{R}}{\int}_{\mathbb{R}}{\widehat{f}}_{110}\left(U\right(\mathrm{cos}\phantom{\rule{thinmathspace}{0ex}}\psi -v\phantom{\rule{thinmathspace}{0ex}}\mathrm{sin}\phantom{\rule{thinmathspace}{0ex}}\psi ),U(\mathrm{sin}\phantom{\rule{thinmathspace}{0ex}}\psi +v\phantom{\rule{thinmathspace}{0ex}}\mathrm{cos}\phantom{\rule{thinmathspace}{0ex}}\psi ),z)\hfill \\ \hfill & \times {\mathrm{e}}^{2\pi \mathrm{i}[x\phantom{\rule{thinmathspace}{0ex}}\mathrm{cos}\phantom{\rule{thinmathspace}{0ex}}\psi +y\phantom{\rule{thinmathspace}{0ex}}\mathrm{sin}\phantom{\rule{thinmathspace}{0ex}}\psi -v(x\phantom{\rule{thinmathspace}{0ex}}\mathrm{sin}\phantom{\rule{thinmathspace}{0ex}}\psi -y\phantom{\rule{thinmathspace}{0ex}}\mathrm{cos}\phantom{\rule{thinmathspace}{0ex}}\psi \left)\right]U}\mid U\mid \mathrm{d}U\phantom{\rule{thinmathspace}{0ex}}\mathrm{d}v\hfill \\ \hfill =& {\int}_{\mathbb{R}}{\int}_{\mathbb{R}}\mid U\mid {\widehat{g}}_{\psi ,100}(U,v,z)\phantom{\rule{thinmathspace}{0ex}}{\mathrm{e}}^{2\pi \mathrm{i}[x\phantom{\rule{thinmathspace}{0ex}}\mathrm{cos}\phantom{\rule{thinmathspace}{0ex}}\psi +y\phantom{\rule{thinmathspace}{0ex}}\mathrm{sin}\phantom{\rule{thinmathspace}{0ex}}\psi -v(x\phantom{\rule{thinmathspace}{0ex}}\mathrm{sin}\phantom{\rule{thinmathspace}{0ex}}\psi -y\phantom{\rule{thinmathspace}{0ex}}\mathrm{cos}\phantom{\rule{thinmathspace}{0ex}}\psi \left)\right)]U}\mathrm{d}U\phantom{\rule{thinmathspace}{0ex}}\mathrm{d}v.\hfill \end{array}$$

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.

Let $f\in {L}^{1}\left({\mathbb{R}}^{3}\right)$ and supp(*f*) **S**_{f}. Also let N_{ψ} defined as in (6) and ${\psi}_{k}={\scriptstyle \frac{\pi}{{N}_{\psi}}}k$ for k = 0, 1, …, N_{ψ} − 1. Then

$$\begin{array}{cc}\hfill f(x,y,z)=& \sum _{k=0}^{{N}_{\psi}-1}{\int}_{-{v}_{m0}}^{{v}_{m0}}\frac{1}{\mu \left(v\right)}{\int}_{\mathbb{R}}\mid U\mid {\widehat{g}}_{{\psi}_{k},100}^{m0}(U,v,z)\hfill \\ \hfill & \times {\mathrm{e}}^{2\pi i[x\phantom{\rule{thinmathspace}{0ex}}\mathrm{cos}\phantom{\rule{thinmathspace}{0ex}}{\psi}_{k}+y\phantom{\rule{thinmathspace}{0ex}}\mathrm{sin}\phantom{\rule{thinmathspace}{0ex}}{\psi}_{k}-v(x\phantom{\rule{thinmathspace}{0ex}}\mathrm{sin}\phantom{\rule{thinmathspace}{0ex}}{\psi}_{k}-y\phantom{\rule{thinmathspace}{0ex}}\mathrm{cos}\phantom{\rule{thinmathspace}{0ex}}{\psi}_{k}\left)\right]U}\mathrm{d}U\phantom{\rule{thinmathspace}{0ex}}\mathrm{d}v,\hfill \end{array}$$

where the sum is over all the gantry rotations and

$$\mu \left(v\right)=\sum _{k=0}^{{N}_{\psi}-1}{1}_{\left\{v:\mid \frac{v\phantom{\rule{thinmathspace}{0ex}}\mathrm{cos}\phantom{\rule{thinmathspace}{0ex}}{\psi}_{k}-\mathrm{sin}\phantom{\rule{thinmathspace}{0ex}}{\psi}_{k}}{\mathrm{cos}\phantom{\rule{thinmathspace}{0ex}}{\psi}_{k}+v\phantom{\rule{thinmathspace}{0ex}}\mathrm{sin}\phantom{\rule{thinmathspace}{0ex}}{\psi}_{k}}\mid \le {v}_{m0}\right\}}\left(v\right)$$

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

The proof follows easily from theorem 3.2.

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:{L}^{1}\left({\mathbb{R}}^{4}\right)\to {L}^{1}\left({\mathbb{R}}^{3}\right)$ by

$$\begin{array}{cc}\hfill & {\widehat{\left(Kh\right)}}_{110}({U}_{0},\phantom{\rule{thinmathspace}{0ex}}{V}_{0},\phantom{\rule{thinmathspace}{0ex}}z;{v}_{1\mathrm{max}})\equiv {\int}_{-{v}_{1\mathrm{max}}}^{{v}_{1\mathrm{max}}}{\widehat{h}}_{1100}\left({U}_{0},{V}_{0},z-\frac{{V}_{0}}{{U}_{0}}{v}_{1},{v}_{1}\right)\mathrm{d}{v}_{1}\hfill \\ \hfill & {\widehat{\left(Kh\right)}}_{110}({U}_{0},\phantom{\rule{thinmathspace}{0ex}}{V}_{0},\phantom{\rule{thinmathspace}{0ex}}z)\equiv \underset{{v}_{1\mathrm{max}}\to \infty}{\mathrm{lim}}{\widehat{\left(Kh\right)}}_{110}({U}_{0},\phantom{\rule{thinmathspace}{0ex}}{V}_{0},\phantom{\rule{thinmathspace}{0ex}}z;{v}_{1\mathrm{max}}).\hfill \end{array}$$

For $f\in {L}^{1}\left({\mathbb{R}}^{3}\right)$ and supp(*f*) **S**_{f}, the PFDR algorithm rebins the planogram data into a stack of linograms according to

$$\begin{array}{cc}\hfill & {\widehat{g}}_{110}^{\mathrm{reb}}({U}_{0},\phantom{\rule{thinmathspace}{0ex}}{V}_{0},\phantom{\rule{thinmathspace}{0ex}}z;{v}_{1\mathrm{max}})\equiv \frac{1}{N\left(-\frac{{V}_{0}}{{U}_{0}},z;{v}_{1\mathrm{max}}\right)}{\widehat{\left(K{g}^{m}\right)}}_{110}({U}_{0},\phantom{\rule{thinmathspace}{0ex}}{V}_{0},\phantom{\rule{thinmathspace}{0ex}}z;{v}_{1\mathrm{max}})\hfill \\ \hfill & N\left(-\frac{{V}_{0}}{{U}_{0}},z;{v}_{1\mathrm{max}}\right)\equiv {\int}_{-{v}_{1\mathrm{max}}}^{{v}_{1\mathrm{max}}}{1}_{\mathbf{M}\cup {\mathbf{S}}_{g}^{c}}\left(0,0,z-\frac{{V}_{0}}{{U}_{0}}{v}_{1},{v}_{1}\right)\phantom{\rule{thinmathspace}{0ex}}\mathrm{d}{v}_{1}.\hfill \end{array}$$

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 *v*_{1max} controls the maximum co-polar acceptance angle to be rebinned. The term *N*(−*V*_{0}/*U*_{0},*z*; *v*_{1max}) 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 *v*_{1max} ≤ *v*_{m1}, then *N*(−*V*_{0}/*U*_{0},*z*; *v*_{1max}) = 2*v*_{1max} and in this case PFDR algorithm is given by

$$\begin{array}{cc}\hfill {\widehat{g}}_{110}^{\mathrm{reb}}({U}_{0},\phantom{\rule{thinmathspace}{0ex}}{V}_{0},\phantom{\rule{thinmathspace}{0ex}}z;{v}_{1\mathrm{max}})& \equiv \frac{1}{2{v}_{1\mathrm{max}}}{\widehat{\left(K{g}^{m}\right)}}_{110}({U}_{0},\phantom{\rule{thinmathspace}{0ex}}{V}_{0},\phantom{\rule{thinmathspace}{0ex}}z;{v}_{1\mathrm{max}})\hfill \\ \hfill & =\frac{1}{2{v}_{1\mathrm{max}}}{\int}_{-{v}_{1\mathrm{max}}}^{{v}_{1\mathrm{max}}}{\widehat{g}}_{1100}^{m}\left({U}_{0},{V}_{0},z-\frac{{V}_{0}}{{U}_{0}}{v}_{1},{v}_{1}\right)\phantom{\rule{thinmathspace}{0ex}}\mathrm{d}{v}_{1}.\hfill \end{array}$$

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

$${P}^{\ast}h(x,y,z)\equiv {\int}_{\mathbb{R}}{\int}_{\mathbb{R}}h(x+{v}_{0}y,{v}_{0},z+{v}_{1}y,{v}_{1})\mathrm{d}{v}_{0}\phantom{\rule{thinmathspace}{0ex}}\mathrm{d}{v}_{1}.$$

Let $h\in {L}^{1}\left({\mathbb{R}}^{4}\right)$. Then

$${\widehat{\left({P}^{\ast}h\right)}}_{100}({U}_{0},y,z)={\widehat{\left(Kh\right)}}_{110}({U}_{0},-y{U}_{0},z).$$

(8)

We have

$$\begin{array}{cc}\hfill {\widehat{\left(Kh\right)}}_{110}({U}_{0},-y{U}_{0},z)\equiv & {\int}_{\mathbb{R}}{\widehat{h}}_{1100}({U}_{0},-y{U}_{0},z+{v}_{1}y,{v}_{1})\mathrm{d}{v}_{1}\hfill \\ \hfill =& {\int}_{\mathbb{R}}{\int}_{\mathbb{R}}{\int}_{\mathbb{R}}h({u}_{0},{v}_{0},z+y{v}_{1},{v}_{1}){\mathrm{e}}^{2\pi \mathrm{i}[{v}_{0}y-{u}_{0}]{U}_{0}}\mathrm{d}{u}_{0}\phantom{\rule{thinmathspace}{0ex}}{\mathrm{d}}_{0}\phantom{\rule{thinmathspace}{0ex}}\mathrm{d}{v}_{1}\hfill \\ \hfill =& {\int}_{\mathbb{R}}{\int}_{\mathbb{R}}{\int}_{\mathbb{R}}h(x+{v}_{0}y,{v}_{0},z+y{v}_{1},{v}_{1}){\mathrm{e}}^{-2\pi \mathrm{i}x{U}_{0}}\mathrm{d}x\phantom{\rule{thinmathspace}{0ex}}\mathrm{d}{v}_{0}\phantom{\rule{thinmathspace}{0ex}}\mathrm{d}{v}_{1}\hfill \\ \hfill =& {\widehat{(P\ast h)}}_{100}({U}_{0},y,z).\hfill \end{array}$$

Thus, the K operator computes the backprojection over *v*_{1} and the substitution *V*_{0} = −*yU*_{0} followed by an inverse Fourier transform in *U*_{0} variable computes the backprojection over *v*_{0}. The rebinned measured data with *v*_{1max} = *H/R* are related to the backprojection operation of the measured data by the relation

$$\begin{array}{cc}\hfill {\widehat{g}}^{reb}({U}_{0},-y{U}_{0},z;H\u2215R)=& \frac{1}{N(y,z;H\u2215R)}{\widehat{\left(K{g}^{m}\right)}}_{110}({U}_{0},-y{U}_{0},z;H\u2215R)\hfill \\ \hfill =& \frac{1}{N(y,z;H\u2215R)}{\widehat{(P\ast {g}^{m})}}_{110}({U}_{0},y,z).\hfill \end{array}$$

As a numerical algorithm, the *K* operator applied to measured data is actually ${\widehat{\left(K{g}^{m}\right)}}_{110}({U}_{0},{V}_{0},z)$ and not ${\widehat{\left(K{g}^{m}\right)}}_{110}({U}_{0},-y{U}_{0},z)$. The distinction is subtle, but in the discrete case the calculation of $\widehat{\left(K{g}^{m}\right)}({U}_{0},-y{U}_{0},z)$ requires 1D interpolations in Fourier space. This can be implemented accurately and efficiently with the CHIRPZ algorithm (Magnusson 1993). Calculating the backprojection over *v*_{0} 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.

Let $f\in {L}^{1}\left({\mathbb{R}}^{3}\right)$, supp(f) **S**_{f}, and let g^{m0} be the restricted planogram as defined above. Then

$$\begin{array}{cc}\hfill {\widehat{g}}_{\psi ,1100}^{m0}({U}_{0},{V}_{0},{u}_{1},{v}_{1})=& \frac{1}{\mid {U}_{0}\mid}{\int}_{\mathbb{R}}{\int}_{\mathbb{R}}{\mathrm{e}}^{2\pi \mathrm{i}[-({V}_{0}\u2215{U}_{0})Y+({u}_{1}+{v}_{1}({V}_{0}\u2215{U}_{0})\left)Z\right]}\hfill \\ \hfill & \times {\widehat{f}}_{-\psi ,111}({U}_{0},Y,Z){1}_{\left\{\right({U}_{0},Y,Z):\mid Y-{v}_{1}Z\mid \le {v}_{m0}\mid {U}_{0}\mid \}}({U}_{0},Y,Z)\mathrm{d}Y\phantom{\rule{thinmathspace}{0ex}}\mathrm{d}Z\hfill \\ \hfill & \times {1}_{\mathit{M}\cup {\mathit{S}}_{g}^{c}}(0,0,{u}_{1},{v}_{1})\hfill \end{array}$$

(9)

$$\begin{array}{cc}\hfill =& {\phantom{\mid}{\widehat{f}}_{-\psi ,100}({U}_{0},y,{u}_{1}-{v}_{1}y)\ast \frac{\mathrm{sin}\left(2\pi {v}_{m0}{U}_{0}y\right)}{\pi {U}_{0}y}\mid}_{y=-{V}_{0}\u2215{U}_{0}}\hfill \\ \hfill & \times {1}_{\mathit{M}\cup {\mathit{S}}_{g}^{c}}(0,0,{u}_{1},{v}_{1}).\hfill \end{array}$$

(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.

Let $f\in {L}^{1}\left({\mathbb{R}}^{3}\right)$ and supp(f) **S**_{f} with $a=L\phantom{\rule{thinmathspace}{0ex}}\mathrm{cos}\left({\scriptstyle \frac{\pi}{2{N}_{\psi}}}\right)-R\phantom{\rule{thinmathspace}{0ex}}\mathrm{sin}\left({\scriptstyle \frac{\pi}{2{N}_{\psi}}}\right)$. Let

$$\left[\begin{array}{c}\hfill {x}_{{\psi}_{k}}\hfill \\ \hfill {y}_{{\psi}_{k}}\hfill \end{array}\right]\equiv \left[\begin{array}{cc}\mathrm{cos}\phantom{\rule{thinmathspace}{0ex}}{\psi}_{k}\hfill & \mathrm{sin}\phantom{\rule{thinmathspace}{0ex}}{\psi}_{k}\hfill \\ -\mathrm{sin}\phantom{\rule{thinmathspace}{0ex}}{\psi}_{k}\hfill & \mathrm{cos}\phantom{\rule{thinmathspace}{0ex}}{\psi}_{k}\hfill \end{array}\right]\left[\begin{array}{c}\hfill x\hfill \\ \hfill y\hfill \end{array}\right],$$

where ${\psi}_{k}={\scriptstyle \frac{\pi}{{N}_{\psi}}}k$ for k = 0, 1, …, N_{ψ} − 1. Then $\left\{{g}_{{\psi}_{k}}^{m0}:k=0,1,\dots ,{N}_{\psi}-1\right\}$ contains a partition of the set of projections over all azimuthal angles. Define

$$\begin{array}{cc}\hfill {f}^{\mathrm{PFDR}}(x,y,z)\equiv & \sum _{k=0}^{{N}_{\psi}-1}{\int}_{-{v}_{m0}}^{{v}_{m0}}{\int}_{\mathbb{R}}\frac{\mid {U}_{0}\mid}{N\left({y}_{{\psi}_{k}},z;{v}_{1\mathrm{max}}\right)}{\widehat{\left(K{g}_{{\psi}_{k}}^{m0}\right)}}_{100}({U}_{0},{v}_{0},z;{v}_{1\mathrm{max}})\hfill \\ \hfill & \times {\mathrm{e}}^{2\pi \mathrm{i}[{x}_{{\psi}_{k}}+{v}_{0}{y}_{{\psi}_{k}}]{U}_{0}}\mathrm{d}{U}_{0}\mathrm{d}{v}_{0},\hfill \end{array}$$

where g^{m0} is defined as in section 1. Then

$$\begin{array}{cc}\hfill {f}^{\mathrm{PFDR}}(x,y,z)=& {\int}_{{\mathbb{R}}^{3}}\mathrm{d}{U}_{0}\mathrm{d}Y\mathrm{d}Z{\widehat{f}}_{111}({U}_{0},Y,Z){\mathrm{e}}^{2\pi \mathrm{i}[x{U}_{0}+yY+zZ]}\hfill \\ \hfill & \times [\sum _{k=0}^{{N}_{\psi}-1}\frac{1}{N({y}_{{\psi}_{k}},z;{v}_{1\mathrm{max}})}{\int}_{-{v}_{1\mathrm{max}}}^{{v}_{1\mathrm{max}}}{1}_{\mathit{M}\cup {\mathit{S}}_{g}^{c}}(0,0,z+{y}_{{\psi}_{k}}{v}_{1},{v}_{1})\phantom{]}\hfill \\ \hfill & \times \phantom{[}{1}_{\left\{\right({U}_{0},Y,Z):\mid Y-{v}_{1}Z\mid \le {v}_{m0}\mid {U}_{0}\mid \}}({U}_{0}\phantom{\rule{thinmathspace}{0ex}}\mathrm{cos}\phantom{\rule{thinmathspace}{0ex}}{\psi}_{k}+Y\phantom{\rule{thinmathspace}{0ex}}\mathrm{sin}\phantom{\rule{thinmathspace}{0ex}}{\psi}_{k},Y\phantom{\rule{thinmathspace}{0ex}}\mathrm{cos}\phantom{\rule{thinmathspace}{0ex}}{\psi}_{k}-{U}_{0}\phantom{\rule{thinmathspace}{0ex}}\mathrm{sin}\phantom{\rule{thinmathspace}{0ex}}{\psi}_{k},Z)\mathrm{d}{v}_{1}].\hfill \end{array}$$

Let

$${A}_{\psi}=\left[\begin{array}{ccc}\mathrm{cos}\phantom{\rule{thinmathspace}{0ex}}\psi \hfill & -\mathrm{sin}\phantom{\rule{thinmathspace}{0ex}}\psi \hfill & \hfill 0\hfill \\ \mathrm{sin}\phantom{\rule{thinmathspace}{0ex}}\psi \hfill & \mathrm{cos}\phantom{\rule{thinmathspace}{0ex}}\psi \hfill & \hfill 0\hfill \\ 0\hfill & 0\hfill & \hfill 1\hfill \end{array}\right].$$

Then with the help of theorem 4.2 we have

$$\begin{array}{cc}\hfill {\int}_{-{v}_{m0}}^{{v}_{m0}}& {\int}_{\mathbb{R}}\frac{\mid {U}_{0}\mid}{N({y}_{\psi},z;{v}_{1\mathrm{max}})}{\widehat{\left(K{g}_{\psi}^{m0}\right)}}_{100}({U}_{0},{v}_{0},z;{v}_{1\mathrm{max}}){\mathrm{e}}^{2\pi \mathrm{i}[{x}_{\psi}+{v}_{0}{y}_{\psi}]{U}_{0}}\mathrm{d}{U}_{0}\mathrm{d}{v}_{0}\hfill \\ \hfill =& {\int}_{\mathbb{R}}\frac{\mid {U}_{0}\mid}{N({y}_{\psi},z;{v}_{1\mathrm{max}})}{\widehat{\left(K{g}_{\psi}^{m0}\right)}}_{110}({U}_{0},-{y}_{\psi}{U}_{0},z;{v}_{1\mathrm{max}}){\mathrm{e}}^{2\pi \mathrm{i}{x}_{\psi}{U}_{0}}\mathrm{d}{U}_{0}\hfill \\ \hfill =& {\int}_{\mathbb{R}}\frac{\mid {U}_{0}\mid}{N({y}_{\psi},z;{v}_{1\mathrm{max}})}{\int}_{-{v}_{1\mathrm{max}}}^{{v}_{1\mathrm{max}}}{\widehat{g}}_{\psi ,1100}^{m0}({U}_{0},-{y}_{\psi}{U}_{0},z+{y}_{\psi}{v}_{1},{v}_{1}){\mathrm{e}}^{2\pi \mathrm{i}{x}_{\psi}{U}_{0}}\mathrm{d}{v}_{1}\mathrm{d}{U}_{0}\hfill \\ \hfill =& \frac{1}{N({y}_{\psi},z;{v}_{1\mathrm{max}})}{\int}_{-{v}_{1\mathrm{max}}}^{{v}_{1\mathrm{max}}}\mathrm{d}{v}_{1}{1}_{\mathbf{M}\cup {\mathbf{S}}_{g}^{c}}(0,0,z+{y}_{{\psi}_{k}}{v}_{1},{v}_{1})\hfill \\ \hfill & \times {\int}_{{\mathbb{R}}^{3}}{\widehat{f}}_{-\psi ,111}({U}_{0},Y,Z){1}_{\left\{\right({U}_{0},Y,Z):\mid Y-{v}_{1}Z\mid \le {v}_{m0}\mid {U}_{0}\mid \}}({U}_{0},Y,Z){\mathrm{e}}^{2\pi \mathrm{i}[{x}_{\psi}{U}_{0}+{y}_{\psi}Y+zZ]}\mathrm{d}{U}_{0}\phantom{\rule{thinmathspace}{0ex}}\mathrm{d}Y\phantom{\rule{thinmathspace}{0ex}}\mathrm{d}Z\hfill \\ \hfill =& \frac{1}{N({y}_{\psi},z;{v}_{1\mathrm{max}})}{\int}_{-{v}_{1\mathrm{max}}}^{{v}_{1\mathrm{max}}}\mathrm{d}{v}_{1}{1}_{\mathbf{M}\cup {\mathbf{S}}_{g}^{c}}(0,0,z+{y}_{{\psi}_{k}}{v}_{1},{v}_{1})\hfill \\ \hfill & \times {\int}_{{\mathbb{R}}^{3}}{\widehat{f}}_{-\psi ,111}({U}_{0},Y,Z){1}_{\left\{\right({U}_{0},Y,Z):\mid Y-{v}_{1}Z\mid \le {v}_{m0}\mid {U}_{0}\mid \}}({U}_{0},Y,Z)\hfill \\ \hfill & \times {\mathrm{e}}^{2\pi \mathrm{i}\langle {A}_{\psi}^{T}{(x,y,z)}^{T},{({U}_{0},Y,Z)}^{T}\rangle}\mathrm{d}{U}_{0}\phantom{\rule{thinmathspace}{0ex}}\mathrm{d}Y\phantom{\rule{thinmathspace}{0ex}}\mathrm{d}Z\hfill \\ \hfill =& \frac{1}{N({y}_{\psi},z;{v}_{1\mathrm{max}})}{\int}_{-{v}_{1\mathrm{max}}}^{{v}_{1\mathrm{max}}}\mathrm{d}{v}_{1}{1}_{\mathbf{M}\cup {\mathbf{S}}_{g}^{c}}(0,0,z+{y}_{{\psi}_{k}}{v}_{1},{v}_{1})\hfill \\ \hfill & \times {\int}_{{\mathbb{R}}^{3}}{\widehat{f}}_{111}\left(A\psi {({U}_{0},Y,Z)}^{T}\right){1}_{\left\{\right({U}_{0},Y,Z):\mid Y-{v}_{1}Z\mid \le {v}_{m0}\mid {U}_{0}\mid \}}({U}_{0},Y,Z)\hfill \\ \hfill & \times {\mathrm{e}}^{2\pi \mathrm{i}\langle {(x,y,z)}^{T},{A}_{\psi}^{T}{({U}_{0},Y,Z)}^{T}\rangle}\mathrm{d}{U}_{0}\phantom{\rule{thinmathspace}{0ex}}\mathrm{d}Y\phantom{\rule{thinmathspace}{0ex}}\mathrm{d}Z\hfill \\ \hfill =& \frac{1}{N({y}_{\psi},z;{v}_{1\mathrm{max}})}{\int}_{-{v}_{1\mathrm{max}}}^{{v}_{1\mathrm{max}}}\mathrm{d}{v}_{1}{1}_{\mathbf{M}\cup {\mathbf{S}}_{g}^{c}}(0,0,z+{y}_{{\psi}_{k}}{v}_{1},{v}_{1})\hfill \\ \hfill & \times {\int}_{{\mathbb{R}}^{3}}{\widehat{f}}_{111}({U}_{0},Y,Z){1}_{\left\{\right({U}_{0},Y,Z):\mid Y-{v}_{1}Z\mid \le {v}_{m0}\mid {U}_{0}\mid \}}\left({A}_{\psi}^{T}{({U}_{0},Y,Z)}^{T}\right){\mathrm{e}}^{2\pi \mathrm{i}[x{U}_{0}+yY+zZ]}\mathrm{d}{U}_{0}\phantom{\rule{thinmathspace}{0ex}}\mathrm{d}Y\phantom{\rule{thinmathspace}{0ex}}\mathrm{d}Z\hfill \end{array}$$

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

$$\begin{array}{cc}\hfill & \sum _{k=0}^{{N}_{\psi}-1}\frac{1}{N({y}_{{\psi}_{k}},z;{v}_{1\mathrm{max}})}{\int}_{-{v}_{1\mathrm{max}}}^{{v}_{1\mathrm{max}}}{1}_{\left\{\right({U}_{0},Y,Z):\mid Y-{v}_{1}Z\mid \le {v}_{m0}\mid {U}_{0}\mid \}}\left({A}_{\psi}^{T}{({U}_{0},Y,Z)}^{T}\right)\hfill \\ \hfill & \times {1}_{\mathbf{M}\cup {\mathbf{S}}_{g}^{c}}(0,0,z+{y}_{{\psi}_{k}}{v}_{1},{v}_{1})\mathrm{d}{v}_{1}.\hfill \end{array}$$

(11)

When we restrict the integration of *v*_{1} to *v*_{1} ≤ *v*_{1max} ≤ *v*_{m1}, we have $N(y,z;{v}_{1\mathrm{max}})=\frac{1}{2{v}_{1\mathrm{max}}}$ and this filter is independent of *y _{ψk}* =

$$\frac{1}{2{v}_{1\mathrm{max}}}\sum _{k=0}^{{N}_{\psi}-1}{\widehat{h}}^{\mathrm{PFDR}}\left({A}_{{\psi}_{k}}^{T}{({U}_{0},Y,Z)}^{T};{v}_{1\mathrm{max}}\right),$$

(12)

where

$${\widehat{h}}^{\mathrm{PFDR}}({U}_{0},Y,Z;{v}_{1\mathrm{max}})\equiv {\int}_{-{v}_{1\mathrm{max}}}^{{v}_{1\mathrm{max}}}{1}_{\left\{\right({U}_{0},Y,Z):\mid Y-{v}_{1}Z\mid \le {v}_{m0}\mid {U}_{0}\mid \}}({U}_{0},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 *v*_{1max}.

When *v*_{1max} → 0, *f*^{PFDR} = *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.

Let $f\in {L}^{1}\left({\mathbb{R}}^{3}\right)$, supp(f) **S**_{f}, and let g^{m0} be the restricted planogram as defined above. Then for v_{1} ≤ v_{m1} we have

$${\widehat{g}}_{\psi ,1010}^{m0}({U}_{0},{v}_{0},Z,{v}_{1})={\widehat{f}}_{-\psi ,111}({U}_{0},{v}_{0}{U}_{0}+{v}_{1}Z,Z){1}_{\{{v}_{0}:\mid {v}_{0}\mid \le {v}_{m0}\}}\left({v}_{0}\right).$$

For a proof see Brasse *et al* (2004).

The next theorem establishes a connection between the direct data ${\widehat{g}}_{1010}^{m0}({U}_{0},{v}_{0},Z,0)$ and the rebinned date ${\widehat{\left(K{g}^{m0}\right)}}_{101}({U}_{0},{v}_{0},Z;{v}_{1\mathrm{max}})$.

Let $f\in {L}^{1}\left({\mathbb{R}}^{3}\right)$ with supp(f) **S**_{f} and v_{1max} ≤ v_{m1} be given. Define

$${\widehat{h}}^{\mathrm{PFDR}}({U}_{0},Y,Z;{v}_{1\mathrm{max}})\equiv {\int}_{-{v}_{1\mathrm{max}}}^{{v}_{1\mathrm{max}}}{1}_{\left\{\right({U}_{0},Y,Z):\mid Y-{v}_{1}Z\mid \le {v}_{m0}\mid {U}_{0}\mid \}}({U}_{0},Y,Z)\mathrm{d}{v}_{1}$$

just as we have above. Then for v_{0} ≤ v_{m0} we have that

$${\widehat{\left(K{g}_{\psi}^{m0}\right)}}_{101}({U}_{0},{v}_{0},Z;{v}_{1\mathrm{max}})={\widehat{f}}_{-\psi ,111}({U}_{0},{v}_{0}{U}_{0},Z){\widehat{h}}^{\mathrm{PFDR}}({U}_{0},{v}_{0}{U}_{0},Z;{v}_{1\mathrm{max}})$$

(13)

$$={\widehat{g}}_{\psi ,1010}^{m0}({U}_{0},{v}_{0},Z,0){\widehat{h}}^{\mathrm{PFDR}}({U}_{0},{v}_{0}{U}_{0},Z;{v}_{1\mathrm{max}}).$$

(14)

Since *v*_{1max} ≤ *v*_{m1}, ${1}_{\mathbf{M}\cup {\mathbf{S}}_{g}^{c}}(0,0,{u}_{1},{v}_{1})=1$ for *v*_{1} ≤ *v*_{1max}. Now we use theorem 4.2 with the substitutions *V*_{0} = −*yU*_{0} and *u*_{1} = *z* + *v*_{1}*y* to establish the relation

$$\begin{array}{cc}\hfill {\widehat{g}}_{\psi ,1100}^{m0}({U}_{0},-y{U}_{0},z+{v}_{1}y,{v}_{1})=& \frac{1}{\mid {U}_{0}\mid}{\int}_{{\mathbb{R}}^{2}}{\widehat{f}}_{-\psi ,111}({U}_{0},Y,Z)\hfill \\ \hfill & \times {1}_{\left\{\right({U}_{0},Y,Z):\mid Y-{v}_{1}Z\mid \le {v}_{m0}\mid {U}_{0}\mid \}}{\mathrm{e}}^{2\pi \mathrm{i}(yY+zZ)}\mathrm{d}Y\phantom{\rule{thinmathspace}{0ex}}\mathrm{d}Z\hfill \end{array}$$

and thus

$$\begin{array}{c}\hfill {\widehat{\left(K{g}_{\psi}^{m0}\right)}}_{110}({U}_{0},-y{U}_{0},z;{v}_{1\mathrm{max}})={\int}_{-{v}_{1\mathrm{max}}}^{{v}_{1\mathrm{max}}}{\widehat{g}}_{\psi ,1100}^{m0}({U}_{0},-y{U}_{0},z+y{v}_{1},{v}_{1})\mathrm{d}{v}_{1}\hfill \\ \hfill =\frac{1}{\mid {U}_{0}\mid}{\int}_{{\mathbb{R}}^{2}}{\widehat{f}}_{-\psi}({U}_{0},Y,Z){\widehat{h}}^{\mathrm{PDFR}}({U}_{0},Y,Z;{v}_{1\mathrm{max}}){\mathrm{e}}^{2\pi \mathrm{i}(yY+zZ)}\mathrm{d}Y\phantom{\rule{thinmathspace}{0ex}}\mathrm{d}Z.\hfill \end{array}$$

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

$$\begin{array}{cc}\hfill & \frac{1}{\mid {U}_{0}\mid}{\widehat{f}}_{-\psi ,111}({U}_{0}.Y,Z){\widehat{h}}^{\mathrm{PFDR}}({U}_{0},YmZ;{v}_{1\mathrm{max}})\hfill \\ \hfill & ={\int}_{\mathbb{R}}{\widehat{\left(K{g}_{\psi}^{m0}\right)}}_{111}({U}_{0},-y{U}_{0},Z;{v}_{1\mathrm{max}}){\mathrm{e}}^{-2\pi \mathrm{i}yY}\mathrm{d}y\hfill \\ \hfill & =\frac{1}{\mid {U}_{0}\mid}{\int}_{\mathbb{R}}{\widehat{\left(K{g}_{\psi}^{m0}\right)}}_{111}({U}_{0},{V}_{0},Z;{v}_{1\mathrm{max}}){\mathrm{e}}^{-2\pi \mathrm{i}{V}_{0}(Y\u2215{U}_{0})}\mathrm{d}{V}_{0}\hfill \\ \hfill & =\frac{1}{\mid {U}_{0}\mid}{\widehat{\left(K{g}_{\psi}^{m0}\right)}}_{101}({U}_{0},Y\u2215{U}_{0},Z;{v}_{1\mathrm{max}}).\hfill \end{array}$$

Making the substitution *v*_{0} = *Y/U*_{0} we have

$$\begin{array}{cc}\hfill {\widehat{\left(K{g}_{\psi}^{m0}\right)}}_{101}({U}_{0},{v}_{0},Z;{v}_{1\mathrm{max}})=& {\widehat{f}}_{-\psi ,111}({U}_{0},{v}_{0}{U}_{0},Z){\widehat{h}}^{\mathrm{PFDR}}({U}_{0},{v}_{0}{U}_{0},Z;{v}_{1\mathrm{max}}),\hfill \\ \hfill =& {\widehat{g}}_{\psi ,1010}^{m0}({U}_{0},{v}_{0},Z,0){\widehat{h}}^{\mathrm{PFDR}}({U}_{0},{v}_{0}{U}_{0},Z;{v}_{1\mathrm{max}}).\hfill \end{array}$$

The filter ${\widehat{h}}^{\mathrm{PFDR}}(X,Y,Z;{v}_{1\mathrm{max}})$ is explicitly derived in the next proposition.

Define

$${\widehat{h}}^{\mathrm{PFDR}}(X,Y,Z;{v}_{1\mathrm{max}})\equiv {\int}_{-{v}_{1\mathrm{max}}}^{{v}_{1\mathrm{max}}}{1}_{\left\{\right(X,Y,Z):\mid Y-{v}_{1}Z\mid \le {v}_{m0}\mid X\mid \}}(X,Y,Z)\phantom{\rule{thinmathspace}{0ex}}\mathrm{d}{v}_{1}.$$

Then

$${\widehat{h}}^{\mathrm{PFDR}}(X,Y,Z;{v}_{1\mathrm{max}})={1}_{\left\{\right(X,Y,Z):\mid Y\mid -{v}_{m0}\mid X\mid \le {v}_{1\mathrm{max}}\mid z\mid \}}(X,Y,Z)\times \{\begin{array}{cc}2{v}_{1\mathrm{max}},\hfill & Z=0,\hfill \\ 2{v}_{1\mathrm{max}},\hfill & \mid Y\mid \le {v}_{m0}\mid X\mid -{v}_{1\mathrm{max}}\mid Z\mid ,\hfill \\ \frac{{v}_{m0}\mid X\mid -\mid Y\mid +{v}_{1\mathrm{max}}\mid Z\mid}{\mid Z\mid},\hfill & -\mid Y\mid \le {v}_{m0}\mid X\mid -{v}_{1\mathrm{max}}\mid Z\mid \le \mid Y\mid ,\hfill \\ \frac{2{v}_{m0}\mid X\mid}{\mid Z\mid},\hfill & {v}_{m0}\mid X\mid -{v}_{1\mathrm{max}}\mid Z\mid \le -\mid Y\mid .\hfill \end{array}\phantom{\}}$$

Furthermore, when v_{0} ≤ v_{m0}

$$\begin{array}{cc}\hfill \frac{\mid {U}_{0}\mid}{{\widehat{h}}^{\mathrm{PFDR}}({U}_{0},{v}_{0}{U}_{0},Z;{v}_{1\mathrm{max}})}& =\{\begin{array}{cc}\frac{\mid {U}_{0}\mid}{2{v}_{1\mathrm{max}}},\hfill & {v}_{1\mathrm{max}}\mid Z\mid \le ({v}_{m0}-\mid {v}_{0}\mid )\mid {U}_{0}\mid ,\hfill \\ \frac{\mid {U}_{0}Z\mid}{{v}_{1\mathrm{max}}\mid Z\mid +({v}_{m0}-\mid {v}_{0}\mid )\mid {U}_{0}\mid},\hfill & \mid {v}_{1\mathrm{max}}\mid Z\mid -{v}_{m0}\mid {U}_{0}\mid \mid \le \mid {v}_{0}{U}_{0}\mid ,\hfill \\ \frac{1}{2{v}_{m0}}\mid Z\mid ,\hfill & ({v}_{m0}+\mid {v}_{0}\mid )\mid {U}_{0}\mid \le {v}_{1\mathrm{max}}\mid Z\mid ,\hfill \end{array}\phantom{\}}\hfill \\ \hfill & =\mathrm{max}\left(\frac{\mid {U}_{0}\mid}{2{v}_{1\mathrm{max}}},\frac{\mid {U}_{0}Z\mid}{{v}_{1\mathrm{max}}\mid Z\mid +({v}_{m0}-\mid {v}_{0}\mid )\mid {U}_{0}\mid},\frac{1}{2{v}_{m0}}\mid Z\mid \right)\hfill \end{array}$$

is a continuous function on $\left\{\right({U}_{0},{v}_{0},Z)\in {\mathbb{R}}^{3}:\mid {v}_{0}\mid \le {v}_{m0}\}$.

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.

Under the hypotheses of theorem 4.5 we have that

$$f(x,y,z)=\sum _{k=0}^{{N}_{\psi}-1}{\int}_{-{v}_{m0}}^{{v}_{m0}}\frac{1}{\mu \left({v}_{0}\right)}{\int}_{{\mathbb{R}}^{2}}\frac{\mid {U}_{0}\mid}{{\widehat{h}}^{\mathrm{PFDR}}({U}_{0},{v}_{0}{U}_{0},Z;{v}_{1\mathrm{max}})}{\widehat{\left(K{g}_{{\psi}_{k}}^{m0}\right)}}_{101}({U}_{0},{v}_{0},Z;{v}_{1\mathrm{max}})\times \phantom{\rule{thickmathspace}{0ex}}{\mathrm{e}}^{2\pi \mathrm{i}\left[\right({x}_{{\psi}_{k}}+{v}_{0}{y}_{{\psi}_{k}}){U}_{0}+zZ]}\mathrm{d}{U}_{0}\phantom{\rule{thinmathspace}{0ex}}\mathrm{d}Z\phantom{\rule{thinmathspace}{0ex}}\mathrm{d}{v}_{0},$$

where

$$\mu \left(v\right)=\sum _{k=0}^{{N}_{\psi}-1}{1}_{\left\{v:\mid \frac{v\phantom{\rule{thinmathspace}{0ex}}\mathrm{cos}\phantom{\rule{thinmathspace}{0ex}}{\psi}_{k}-\mathrm{sin}\phantom{\rule{thinmathspace}{0ex}}{\psi}_{k}}{\mathrm{cos}\phantom{\rule{thinmathspace}{0ex}}{\psi}_{k}+v\phantom{\rule{thinmathspace}{0ex}}\mathrm{sin}\phantom{\rule{thinmathspace}{0ex}}{\psi}_{k}}\mid \le {v}_{m0}\right\}}\left(v\right).$$

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, *U*_{0}, is replaced by

$$\frac{\mid {U}_{0}\mid}{{\widehat{h}}^{\mathrm{PFDR}}({U}_{0},{v}_{0}{U}_{0},Z;{v}_{1\mathrm{max}})}.$$

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 *v*_{m1} = 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 *v*_{m1}. 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 *v*_{m0}.

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 *v*_{1}), filtering the stack of rebinned linograms, and then backprojecting over *v*_{0}. The planogram filtered backprojection algorithm filters the data and then performs 2D backprojection over both *v*_{0} and *v*_{1}. The filter used by the planogram filtered backprojection algorithm is given by

$$\frac{1}{\sum _{k=0}^{{N}_{\psi}-1}\widehat{p}\left({A}_{-{\psi}_{k}}{({U}_{0},{v}_{0}{U}_{0}+{v}_{1}{U}_{1},{U}_{1})}^{T}\right)},$$

where

$$\frac{1}{\widehat{p}({U}_{0},{v}_{0}{U}_{0}+{v}_{1}Z,Z)}=\{\begin{array}{cc}\frac{\mid {U}_{0}\mid}{2{v}_{1\mathrm{max}}},\hfill & \mid {U}_{0}\mid {v}_{m0}-\mid Z\mid {v}_{1\mathrm{max}}\ge \mid {v}_{0}{U}_{0}+{v}_{1}Z\mid ,\hfill \\ \frac{\mid {U}_{0}Z\mid}{{v}_{1\mathrm{max}}\mid Z\mid +{v}_{m0}\mid {U}_{0}\mid -\mid {v}_{0}{U}_{0}+{v}_{1}Z\mid},\hfill & \mid {v}_{1\mathrm{max}}\mid Z\mid -{v}_{m0}\mid {U}_{0}\mid \mid \le \mid {v}_{0}{U}_{0}+{v}_{1}Z\mid ,\hfill \\ \frac{\mid Z\mid}{2{v}_{m0}},\hfill & \mid {U}_{0}\mid {v}_{m0}-\mid Z\mid {v}_{1\mathrm{max}}\le -\mid {v}_{0}{U}_{0}+{v}_{1}Z\mid \hfill \end{array}\phantom{\}}$$

Note that for $\left\{\right({U}_{0},{v}_{0},Z)\in {\mathbb{R}}^{3}:\mid {v}_{0}\mid \le {v}_{m0}\}$

$${\phantom{\mid}\frac{\mid {U}_{0}\mid}{{\widehat{h}}^{\mathrm{PFDR}}({U}_{0},{v}_{0}{U}_{0},Z;{v}_{1\mathrm{max}})}=\frac{1}{\widehat{p}({U}_{0},{v}_{0}{U}_{0}+{v}_{1}Z,Z)}\mid}_{{v}_{1}=0}.$$

The PFDRX algorithm is computationally more efficient than the planogram filtered backprojection algorithm by implementing a more efficient method of backprojecting over *v*_{1} 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.

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({s}_{A},{s}_{B},{t}_{A},{t}_{B})\equiv {\int}_{\mathbb{R}}f\left(\frac{1}{2}({s}_{A}-{s}_{B})-\frac{1}{2R}({s}_{A}+{s}_{B})y,y,\frac{1}{2}({t}_{A}+{t}_{B})-\frac{1}{2R}({t}_{A}-{t}_{B})y\right)\mathrm{d}y.$$

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

$$\begin{array}{cc}\hfill & {q}^{m}({s}_{A},{s}_{B},{t}_{A},{t}_{B})\equiv q({s}_{A},{s}_{B},{t}_{A},{t}_{B}){1}_{\left\{\right({s}_{A},{s}_{B},{t}_{A},{t}_{B})\in {[-L,L]}^{2}\times {[-H,H]}^{2}\}}({s}_{A},{s}_{B},{t}_{A},{t}_{B}),\hfill \\ \hfill & {q}^{m0}({s}_{A},{s}_{B},{t}_{A},{t}_{B})\equiv {q}^{m}({s}_{A},{s}_{B},{t}_{A},{t}_{B}){1}_{\left\{\right({s}_{A},{s}_{B})\in {\mathbb{R}}^{2}:\frac{1}{2R}({s}_{A}+{s}_{B})\le {v}_{m0}\}}({s}_{A},{s}_{B}),\hfill \end{array}$$

respectively.

Now let

$$\begin{array}{cc}\hfill & {A}_{0}\equiv \left[\begin{array}{cc}1\hfill & \hfill R\hfill \\ -1\hfill & \hfill R\hfill \end{array}\right],\phantom{\rule{1em}{0ex}}{B}_{0}\equiv {\left({A}_{0}^{T}\right)}^{-1}=\left[\begin{array}{cc}\frac{1}{2}\hfill & \hfill \frac{1}{2R}\hfill \\ -\frac{1}{2}\hfill & \hfill \frac{1}{2R}\hfill \end{array}\right]\hfill \\ \hfill & {A}_{1}\equiv \left[\begin{array}{cc}1\hfill & \hfill R\hfill \\ 1\hfill & \hfill -R\hfill \end{array}\right].\hfill \end{array}$$

Then we have that

$${A}_{0}\left[\begin{array}{c}\hfill {u}_{0}\hfill \\ \hfill {v}_{0}\hfill \end{array}\right]=\left[\begin{array}{c}\hfill {s}_{A}\hfill \\ \hfill {s}_{B}\hfill \end{array}\right]\phantom{\rule{1em}{0ex}}\text{and}\phantom{\rule{1em}{0ex}}{A}_{1}\left[\begin{array}{c}\hfill {u}_{1}\hfill \\ \hfill {v}_{1}\hfill \end{array}\right]=\left[\begin{array}{c}\hfill {t}_{A}\hfill \\ \hfill {t}_{B}\hfill \end{array}\right]$$

for the detector coordinates (*s _{A},s_{B},t_{A},t_{B}*) and the planogram coordinates (

$$q({A}_{0}{({u}_{0},{v}_{0})}^{T},{A}_{1}{({u}_{1},{v}_{1})}^{T})=g({u}_{0},{v}_{0},{u}_{1},{v}_{1}).$$

Then by theorem 1.1 we have that

$$\begin{array}{cc}\hfill {\widehat{g}}_{1100}({U}_{0},{V}_{0},{u}_{1},{v}_{1})& =\mathcal{F}\left\{g\right(\cdot ,\cdot ,{u}_{1},{v}_{1}\left)\right\}({U}_{0},{V}_{0})\hfill \\ \hfill & =\mathcal{F}\left\{q\right({A}_{0}{(\cdot ,\cdot )}^{T},{A}_{1}{({u}_{1},{v}_{1})}^{T}\left)\right\}({U}_{0},{V}_{0})\hfill \\ \hfill & ={\mid \mathrm{det}\phantom{\rule{thinmathspace}{0ex}}{A}_{0}\mid}^{-1}\mathcal{F}\left\{q\right(\cdot ,\cdot ,{A}_{1}{({u}_{1},{v}_{1})}^{T}\left)\right\}\left({B}_{0}{({U}_{0},{V}_{0})}^{T}\right)\hfill \\ \hfill & =\frac{1}{2R}{\widehat{q}}_{1100}({B}_{0}{({U}_{0},{V}_{0})}^{T},{A}_{1}{({u}_{1},{v}_{1})}^{T}).\hfill \end{array}$$

(15)

The basis for the PFDR algorithm lies in the relation

$${\widehat{g}}_{1100}({U}_{0},{V}_{0},{u}_{1},0)={\widehat{g}}_{1100}\left({U}_{0},{V}_{0},{u}_{1}-{v}_{1}\frac{{V}_{0}}{{U}_{0}},{v}_{1}\right),$$

which holds for ${v}_{1}\in \mathbb{R}$. By the above relationship and (15), we have that

$${\widehat{q}}_{1100}\left({B}_{0}\left[\begin{array}{c}\hfill {U}_{0}\hfill \\ \hfill {V}_{0}\hfill \end{array}\right],{A}_{1}\left[\begin{array}{c}\hfill {u}_{1}\hfill \\ \hfill 0\hfill \end{array}\right]\right)={\widehat{q}}_{1100}\left({B}_{0}\left[\begin{array}{c}\hfill {U}_{0}\hfill \\ \hfill {V}_{0}\hfill \end{array}\right],{A}_{1}\left[\begin{array}{c}{u}_{1}-{v}_{1}\frac{{V}_{0}}{{U}_{0}}\hfill \\ {v}_{1}\hfill \end{array}\right]\right).$$

(16)

Now let (*X _{A},X_{B}*) be such that

$${B}_{0}\left[\begin{array}{c}\hfill {U}_{0}\hfill \\ \hfill {V}_{0}\hfill \end{array}\right]=\left[\begin{array}{c}\hfill {X}_{A}\hfill \\ \hfill {X}_{B}\hfill \end{array}\right]$$

and thus

$$\left[\begin{array}{c}\hfill {U}_{0}\hfill \\ \hfill {V}_{0}\hfill \end{array}\right]={B}_{0}^{-1}\left[\begin{array}{c}\hfill {X}_{A}\hfill \\ \hfill {X}_{B}\hfill \end{array}\right]={A}_{0}^{T}\left[\begin{array}{c}\hfill {X}_{A}\hfill \\ \hfill {X}_{B}\hfill \end{array}\right]=\left[\begin{array}{cc}1\hfill & -1\hfill \\ R\hfill & R\hfill \end{array}\right]\left[\begin{array}{c}\hfill {X}_{A}\hfill \\ \hfill {X}_{B}\hfill \end{array}\right]=\left[\begin{array}{c}{X}_{A}-{X}_{B}\hfill \\ R({X}_{A}+{X}_{B})\hfill \end{array}\right].$$

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

$${\widehat{q}}_{1100}({X}_{A},{X}_{B},{u}_{1},{u}_{1})={\widehat{q}}_{1100}\left({X}_{A},{X}_{B},{u}_{1}-{X}_{B}\left(\frac{2R{v}_{1}}{{X}_{A}-{X}_{B}}\right),{u}_{1}-{X}_{A}\left(\frac{2R{v}_{1}}{{X}_{A}-{X}_{B}}\right)\right).$$

(17)

Thus, the PFDR algorithm in detector coordinates is given by

$${q}_{110}^{\mathit{reb}}({X}_{A},{X}_{B},z)\equiv \frac{1}{N\left(-R\frac{{X}_{A}+{X}_{B}}{{X}_{A}-{X}_{B}},z;{v}_{1\mathrm{max}}\right)}\times {\int}_{-{v}_{1\mathrm{max}}}^{{v}_{1\mathrm{max}}}{\widehat{q}}_{1100}^{m}\left({X}_{A},{X}_{B},z-{X}_{B}\left(\frac{2R{v}_{1}}{{X}_{A}-{X}_{B}}\right),z-{X}_{A}\left(\frac{2R{v}_{1}}{{X}_{A}-{X}_{B}}\right)\right)\mathrm{d}{v}_{1}.$$

(18)

Suppose that we wish to perform the PFDR algorithm on the measured data set and we have *N _{s}* detector pixels in the transaxial direction and

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.

By letting *z* = *u*_{1} and $t={\scriptstyle \frac{2R{v}_{1}}{{X}_{A}-{X}_{B}}}$ in equation (17) we have

$${\widehat{q}}_{1100}({X}_{A},{X}_{B},z,z)={\widehat{q}}_{1100}({X}_{A},{X}_{B},z-{X}_{B}t,z-{X}_{A}t).$$

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

$${\widehat{q}}_{1100}^{m}({X}_{A},{X}_{B},z,z)={\widehat{q}}_{1100}^{m}({X}_{A},{X}_{B},z-{X}_{B}t,z-{X}_{A}t)-{\int}_{0}^{t}\frac{\mathrm{d}}{\mathrm{d}{t}^{\prime}}{\widehat{q}}_{1100}^{m}({X}_{A},{X}_{B},z-{X}_{B}{t}^{\prime},z-{X}_{A}{t}^{\prime})\mathrm{d}{t}^{\prime},$$

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 *v*_{0} ≤ *v*_{m0} and *v*_{m1} ≤ *v*_{m1}, 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 ${\widehat{h}}^{\mathrm{PFDR}}({U}_{0},{v}_{0}{U}_{0},Z){\mid}_{{U}_{0}=0}$, we may not filter the rebinned data ${\widehat{\left(K{g}^{m0}\right)}}_{101}({U}_{0},{v}_{0},Z;{v}_{m1})$. In theorem 4.7, exact reconstruction is performed by filtering the rebinned data by

$$\frac{\mid {U}_{0}\mid}{{\widehat{h}}^{\mathrm{PFDR}}({U}_{0},{v}_{0}{U}_{0},Z;{v}_{1\mathrm{max}})}.$$

The *U*_{0} term prevents this filter from blowing up near *U*_{0} = 0 and thus provides a stable filter.

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.

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, *g*^{m0} for *v*_{1} ≤ *v*_{m1}. In the PFDR+FBP algorithm one can rebin the measured data set, *g ^{m}*, or just the restricted data set,

The restricted data set for *v*_{m0} = 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

$${\mathrm{e}}^{-2\pi \mathrm{i}y{v}_{1}Z}{\widehat{f}}_{-\psi ,101}(X,y,Z).$$

Then for each *v*_{0} we used the CHIRPZ (Magnusson 1993) algorithm to calculate a Fourier transform in the *y* dimension to output frequency samples spaced by *T*_{v0}*X* units, where *T*_{v0} is the sampling distance in the *v*_{0} variable. Thus, we have calculated

$${\widehat{f}}_{-\psi ,101}(X,{v}_{0}X+{v}_{1}Z,Z)={\widehat{g}}_{\psi}(X,{v}_{0},Z,{v}_{1}).$$

To prevent aliasing, one should multiply ${\widehat{g}}_{\psi}(X,{v}_{0},Z,{v}_{1})$ by

$${1}_{\left\{Y\in \mathbb{R}:\mid Y\mid \le \frac{1}{2Ty}\right\}}({v}_{0}X+{v}_{1}Z),$$

where ${\scriptstyle \frac{1}{2{T}_{y}}}$ is the Nyquist frequency of our image and *T _{y}* is the sampling distance in the

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

$$\begin{array}{cc}\hfill & \frac{\mid {U}_{0}\mid}{{\widehat{h}}^{\mathrm{PFDR}}({U}_{0},{v}_{0}{U}_{0},Z;{v}_{1\mathrm{max}})}(1-{\left(2{T}_{u0}{U}_{0}\right)}^{8})(1-{\left(2{T}_{u1}Z\right)}^{8}),\hfill \\ \hfill & \frac{1}{\sum _{k=0}^{{N}_{\psi}-1}\widehat{p}\left({A}_{-\psi}{({U}_{0},{v}_{0}{U}_{0}+{v}_{1}{U}_{1},{U}_{1})}^{T}\right)}(1-{\left(2{T}_{u0}{U}_{0}\right)}^{8})(1-{\left(2{T}_{u1}Z\right)}^{8}),\hfill \end{array}$$

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

$$\mid {U}_{0}\mid (1-{\left(2{T}_{u0}{U}_{0}\right)}^{8}),$$

where *T*_{u0} and *T*_{u1} are the sampling rates for *u*_{0} and *u*_{1}, respectively. The terms (1−(2*T*_{u0}*U*_{0})^{8}) and (1 − (2*T*_{u1}*Z*)^{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 ${g}_{110}^{\mathrm{reb}}({U}_{0},{V}_{0},z;{v}_{1\mathrm{max}})$ for all *z* ≤ *c* + *av*_{1max}. This maximum value for *z* was chosen by empirical observations.

From theorem 4.4, we have

$${\int}_{{\mathbb{R}}^{2}}{g}^{m0}({u}_{0},{v}_{0},{u}_{1},{v}_{1})\phantom{\rule{thinmathspace}{0ex}}\mathrm{d}{u}_{0}\phantom{\rule{thinmathspace}{0ex}}\mathrm{d}{u}_{1}={\int}_{{\mathbb{R}}^{3}}f(x,y,z)\phantom{\rule{thinmathspace}{0ex}}\mathrm{d}x\phantom{\rule{thinmathspace}{0ex}}\mathrm{d}y\phantom{\rule{thinmathspace}{0ex}}\mathrm{d}z$$

for *v*_{0} ≤ *v*_{m0} and *v*_{1} ≤ *v*_{m1}. 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

$$\frac{{(1+{v}_{0}^{2}+{v}_{1}^{2})}^{3\u22152}}{{T}_{v0}{T}_{v1}}.$$

Note that

$${(1+{v}_{0}^{2}+{v}_{1}^{2})}^{-3\u22152}\phantom{\rule{thinmathspace}{0ex}}\mathrm{d}{v}_{0}\phantom{\rule{thinmathspace}{0ex}}\mathrm{d}{v}_{1}=\mathrm{cos}\phantom{\rule{thinmathspace}{0ex}}\theta \phantom{\rule{thinmathspace}{0ex}}\mathrm{d}\theta \phantom{\rule{thinmathspace}{0ex}}\mathrm{d}\phi ,$$

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

$${(1+{v}_{0}^{2}+{v}_{1}^{2})}^{-1\u22152}{(1+{v}_{0}^{2}+{v}_{1}^{2})}^{3\u22152}=1+{v}_{0}^{2}+{v}_{1}^{2}.$$

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 _{ψ}* :

Image reconstruction computation times (s) with *v*_{1max} =tan(15°) and voxel side length of *d*/2, where *d* = 2.1 mm is the detector sampling distance.

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.

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 *v*_{m0} = 0.2691 ≈ tan(15°). To fulfill the Orlov conditions, we need *N _{ψ}* = 6 detector positions (see equation (6)) at rotation angles

Figure 7 displays noiseless reconstructions with *v*_{1max} = tan(15°), and figure 8 displays plots of profiles through these images. Figure 9 displays noiseless reconstructions with ${v}_{1\mathrm{max}}={\scriptstyle \frac{H}{R}}\approx \mathrm{tan}(29\xb0)$ and figure 10 displays plots of profiles through these images. Note that the maximum oblique LOR measured has a value of ${v}_{1}={\scriptstyle \frac{H}{R}}$. 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.

Reconstructions using the PFDR+FBP, PFDRX and Brasse algorithms with *v*_{1max} = tan(15°) from simulated noiseless data. Axial slices are on the top row, and coronal slices are on the bottom row.

Reconstructions using the PFDR+FBP, PFDRX and Brasse algorithms with *v*_{1max} = tan(29°) from simulated noiseless data. Axial slices are on the top row, and coronal slices are on the bottom row.

Reconstructions using the PFDR+FBP, PFDRX and Brasse algorithms with *v*_{1max} = 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 *v*_{1max} = tan(15°) and *v*_{1max} = tan(29°), respectively.

Root mean squared error (RMSE) for large hot (left) and cold (right) disks of images reconstructed from simulated noiseless data with *v*_{1max} = tan(15°).

Root mean squared error (RMSE) for large hot (left) and cold (right) disks of images reconstructed from simulated noiseless data with *v*_{1max} = tan(29°).

The relative *L*^{2} error of images reconstructed from the noisy data is shown in table 3. The relative *L*^{2} error is defined by ${\scriptstyle \frac{{\Vert \stackrel{~}{f}-f\Vert}_{2}}{{\Vert f\Vert}_{2}}}$, where *f* is the true image and $\stackrel{~}{f}$ is the reconstructed image.

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.

Signal-to-noise and contrast recovery coefficient measurements of images reconstructed from simulated noisy data with *v*_{1max} = 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.

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 *v*_{1max} = tan(15°). These images are shown in figure 16.

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 *v*_{0} 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.

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).

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.

- 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.

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |