Chapter 6: Cone-Beam CT and Feldkamp (FDK) Reconstruction

The fan beam of the previous chapter images one slice with a one-dimensional detector. Extending the detector to a two-dimensional flat panel makes the X-ray beam spread as a cone instead of a fan, and a single rotation of the source captures projection data covering the entire object at once. This is cone-beam CT, widely used in dental CT, angiography systems (C-arms), and industrial CT. The geometry is the $z$-extension of the fan beam: the source travels the circle of radius $D_{\mathrm{SO}}$ in the $z=0$ plane, parameterized by the source angle $\beta$, and the detector becomes two-dimensional with an in-plane coordinate $u$ (unit vector $\hat e_u$) and a rotation-axis coordinate $v$ (unit vector $\hat e_v$). As before, computations use a virtual detector scaled down to pass through the isocenter (origin).

$$S(\beta) = D_{\mathrm{SO}}(\cos\beta,\ \sin\beta,\ 0), \qquad P(u,v) = u\,\hat e_u + v\,\hat e_v, \quad \hat e_u = (-\sin\beta, \cos\beta, 0),\ \hat e_v = (0,0,1)$$

The circular orbit is nevertheless mechanically simple, and most practical scanners use it. Its de-facto standard reconstruction is the Feldkamp–Davis–Kress (FDK, 1984) algorithm, the natural 3D extension of fan-beam FBP. It has three steps. (1) Weight each detector pixel by the cosine factor $D_{\mathrm{SO}}/\sqrt{D_{\mathrm{SO}}^2+u^2+v^2}$, which compensates for the path obliquity of tilted rays. (2) Convolve each detector row (the 1D data along $u$ at fixed $v$) with the ramp filter. (3) Backproject in 3D with the distance weight $(D_{\mathrm{SO}}/L)^2$. Here $L$ is the depth from the source along the central-ray direction, and $(u^*, v^*)$ is where the ray through voxel $(x,y,z)$ hits the virtual detector; the backprojection bilinearly interpolates the value there and sums over all $\beta$.

$$\tilde p(\beta, u, v) = \frac{D_{\mathrm{SO}}}{\sqrt{D_{\mathrm{SO}}^2 + u^2 + v^2}}\; p(\beta, u, v), \qquad q(\beta, \cdot, v) = \tilde p(\beta, \cdot, v) * h$$ $$f(x,y,z) = \frac{\pi}{N_\beta} \sum_{\beta} \left(\frac{D_{\mathrm{SO}}}{L}\right)^{\!2} q\bigl(\beta,\, u^*,\, v^*\bigr)$$ $$L = D_{\mathrm{SO}} - x\cos\beta - y\sin\beta, \qquad u^* = \frac{D_{\mathrm{SO}}\,(-x\sin\beta + y\cos\beta)}{L}, \qquad v^* = \frac{D_{\mathrm{SO}}\, z}{L}$$

In the central plane $z=0$ we get $v^* = 0$, and the cosine weight, the filter, and the backprojection all reduce exactly to the fan-beam FBP formulas. The central slice of FDK therefore agrees exactly with 2D fan-beam FBP. Away from the plane, the rays tilt more and more, and the data missing under Tuy's condition takes its toll. The top and bottom edges of structures that vary quickly along $z$ become blurred and their values sink (intensity drops). That is the typical cone artifact, and the simulations below let you verify it directly.

Simulation: cone-beam geometry in 3D

The cone-beam acquisition geometry is shown in 3D. Rotate the viewpoint by left-dragging and zoom with the mouse wheel. Move the source along its orbit with the $\beta$ slider or the play button, and change the source distance with the $D_{\mathrm{SO}}$ slider to see the cone half-angle respond. What is displayed is the physical detector on the far side of the source (source–detector distance $= 2\cdot D_{\mathrm{SO}}$, magnification 2); the reconstruction math uses the virtual detector scaled down to the isocenter (half-width 1.4). Both represent the same ray bundle.

Simulation: running the FDK reconstruction

A 3D Shepp-Logan phantom is rasterized, forward-projected with a cone beam, and reconstructed with FDK, all inside your browser. Choose the resolution $N$ (the volume has $N^3$ voxels), the number of projections, and the reconstruction filter, then press Run. The computation uses WebGPU when available and falls back to the CPU (Web Worker) otherwise. The result appears in the three-plane viewer below.

Simulation: orthogonal three-plane viewer

The reconstructed volume is displayed in three orthogonal planes: axial (fixed $z$), coronal (fixed $y$), and sagittal (fixed $x$). In each column the top image is the FDK reconstruction and the bottom one is the true phantom; the slider moves the slice position. The slice RMSE compares against the truth inside the FOV cylinder ($r \le 0.9$, $|z| \le 0.9$).

The quickest way to see cone artifacts is to move the axial $z$ slider away from the center ($z = 0$). Near the center the RMSE is smallest and the slice matches the truth almost perfectly. Within the extent of the object, the RMSE grows monotonically with $|z|$ as the tops and bottoms of the ellipsoids blur and sink in value (beyond the top of the phantom, $|z| \approx 0.8$, the slice is nearly empty and the RMSE drops again). The coronal and sagittal views show it at a glance: the upper and lower edges of the structures (edges along $z$) are blurred and darkened, while in-plane edges (along $x$ and $y$) stay sharp even in the outer slices. That contrast is worth a close look.

The limit of a circular orbit

Cone-beam CT collects projections of a whole volume in a single rotation thanks to its 2D detector, but the circular orbit violates Tuy's condition everywhere except the central plane, so the data is inherently incomplete. The FDK algorithm is the approximate reconstruction that extends fan-beam FBP to 3D in three steps: cosine weighting, row-wise ramp filtering, and distance-weighted backprojection. Its central slice agrees exactly with 2D fan-beam FBP, and quality degrades with $|z|$ through cone artifacts. Everything so far has been analytic, one-pass reconstruction. The next chapter changes perspective and solves reconstruction as an optimization problem.

Voxels and partial-volume effects

Each value in a 3D image is an average over a voxel with finite width and height. If one voxel contains more than one material, as at a bone–soft-tissue boundary, its CT number falls between the component values. This is the partial-volume effect. Thinner sections and smaller in-plane pixels separate finer structures, but fewer photons contribute to each voxel and noise rises. In-plane pixel size depends on the field of view and matrix; longitudinal resolution also depends on detector collimation, helical pitch, reconstruction method, and section thickness.

References

• Feldkamp LA, Davis LC, Kress JW. Practical cone-beam algorithm. Journal of the Optical Society of America A 1, 612–619 (1984).
• Tuy HK. An inversion formula for cone-beam reconstruction. SIAM Journal on Applied Mathematics 43, 546–552 (1983).