Computed Tomography (CT)

Computed tomography (CT) builds a cross-sectional image of an object from X-ray measurements taken along many straight lines through it. A conventional radiograph records, at each point of the film or detector, the total attenuation along one ray, so structures at different depths are superimposed. CT instead measures the attenuation along a very large number of rays crossing a thin slice from many directions, about a million per picture in Hounsfield’s 1979 description, and a computer solves for the attenuation coefficient at every point of the slice. Hounsfield, who built the first clinical scanner, summarized the gain in his 1979 Nobel lecture: depth is no longer superimposed, soft tissues that a radiograph cannot separate become distinguishable, and the reconstructed values can be measured quantitatively relative to water (Hounsfield, 1979).

For computer vision, CT is unusual among sensors in two ways. The image is not recorded but computed, by solving an inverse problem whose physics and sampling decide which artifacts appear. And its voxel values are calibrated: a CT number in Hounsfield units refers to the attenuation of water and air, so intensities mean the same thing across scans in a way camera pixel values do not.

How It Works

Attenuation and the line integral

An X-ray beam passing through matter loses photons to photoelectric absorption and to Compton scattering, the two dominant mechanisms at diagnostic energies of about 20–150 keV (Kak & Slaney, 1988, §4.1). For a narrow beam of monoenergetic photons crossing a slice with linear attenuation coefficient $\mu(x, y)$, the number of photons $N_d$ reaching the detector relates to the number $N_\text{in}$ entering by the Beer–Lambert law, which Kak and Slaney write as

$$N_d = N_\text{in} \exp\left(-\int_\text{ray} \mu(x, y)\, ds\right) \quad\Longleftrightarrow\quad \int_\text{ray} \mu(x, y)\, ds = \ln \frac{N_\text{in}}{N_d},$$

with $ds$ an element of length along the ray (Kak & Slaney, 1988, ch. 4, Eq. 7). Taking the logarithm of the measured intensity, relative to an unattenuated reference, turns each detector reading into a line integral of $\mu$. Hounsfield described the same step: absorption is the logarithm of the ratio of X-ray intensity to detector reading (Hounsfield, 1979).

The linear model holds only under assumptions that real scanners break. Kak and Slaney list them: monoenergetic photons, detectors insensitive to scatter, a very narrow beam with no partial-volume effects, and a very small aperture. The X-ray tube emits a spectrum, and $\mu$ depends on energy, so the measurement becomes

$$N_d = \int S_\text{in}(E) \exp\left(-\int_\text{ray} \mu(x, y, E)\, ds\right) dE,$$

where $S_\text{in}(E)$ is the incident photon spectrum (Kak & Slaney, 1988, ch. 4, Eq. 9). Most tissues attenuate low-energy photons more strongly, so the beam’s mean energy rises as it travels, which is called beam hardening. Its consequences appear under Artifacts below.

Projections, sinograms and the Radon transform

Kak and Slaney parametrize a line by its angle $\theta$ and signed distance $t$ from the origin, $x\cos\theta + y\sin\theta = t$, and define the projection at angle $\theta$ as the set of line integrals

$$P_\theta(t) = \int_{-\infty}^{\infty}\int_{-\infty}^{\infty} f(x, y)\, \delta(x\cos\theta + y\sin\theta - t)\, dx\, dy,$$

which, as a function of both $\theta$ and $t$, is the Radon transform of $f$ (Kak & Slaney, 1988, ch. 3, Eq. 3). Stacking the projections as rows or columns, one per angle, gives an image called a sinogram, because a single point in the object traces a sinusoid $t = x_0\cos\theta + y_0\sin\theta$ across it. The transform is named after Johann Radon, whose 1917 paper solved the problem of recovering a function from its integrals along lines; Cormack learned of Radon’s work only in the early 1970s, after solving the problem independently (Cormack, 1979).

A set of parallel rays at one angle is a parallel projection; rays diverging from a single source point to a row of detectors form a fan-beam projection (Kak & Slaney, 1988, §3.1). The mathematics is simplest for parallel projections, and fan-beam data can be re-sorted into parallel form or reconstructed with fan-beam versions of the same algorithms.

The central slice theorem

The key to reconstruction is what Kak and Slaney call the Fourier slice theorem, also known as the central slice or projection-slice theorem: the one-dimensional Fourier transform of a parallel projection at angle $\theta$ equals the two-dimensional Fourier transform of the object along a line through the origin at the same angle (Kak & Slaney, 1988, §3.2). With the Fourier transforms of the object, $F(u, v)$, and of the projection, $S_\theta(w)$, both defined with the kernel $e^{-j2\pi(\cdot)}$,

$$S_\theta(w) = F(w\cos\theta,\; w\sin\theta).$$

It follows that projections at enough angles fill the 2D frequency plane along radial lines, and an inverse 2D Fourier transform recovers the object.

Filtered back-projection

Rewriting the inverse 2D Fourier transform in polar coordinates gives filtered back-projection, which Kak and Slaney describe as used in almost all applications of straight-ray tomography in 1988 (Kak & Slaney, 1988, ch. 3, Eqs. 24–35). Each projection is filtered with the ramp $|w|$, the Jacobian of the polar change of variables, and the filtered projection is then back-projected: its value at $t$ is added to every pixel on the line $x\cos\theta + y\sin\theta = t$, smeared back across the image. Summing over angles completes the reconstruction. Without the filter, back-projection over-weights the low frequencies, where all the radial lines meet, and returns a blurred image (Kak & Slaney, 1988, §3.3.1).

The ramp grows with frequency, and so does its amplification of noise. Multiplying it by a smooth window lowers the noise at the cost of resolution, the trade-off between resolution and noise that Hounsfield (1979) described. Fewer views leave gaps between the radial lines at high frequencies, and the reconstruction shows thin streaks: the filtered back-projection of a point from a few views is a star-shaped pattern (Kak & Slaney, 1988, §5.1.1). This is aliasing in the angular direction; Kak and Slaney’s rule of thumb asks for roughly as many projections as rays per projection. Medical image reconstruction derives the formula from the slice theorem, compares filters under photon noise and gives the view-count rule with its conditions.

Iterative and learned reconstruction

Filtered back-projection has no model of the photon statistics and needs complete, evenly spaced views. Iterative methods instead discretize the slice and solve the large, sparse system of ray equations by repeated projection, as in the algebraic reconstruction technique (ART) and its relatives SIRT and SART (Kak & Slaney, 1988, ch. 7); model-based iterative reconstruction adds a model of the data statistics and prior information about the image (Adler & Öktem, 2018). Learned methods either post-process an FBP image with a network, as FBPConvNet does for sparse-view data (Jin et al., 2017), or build the forward projector into the network, as learned primal–dual reconstruction does (Adler & Öktem, 2018). Medical image reconstruction describes these methods, low-dose CT denoising and the evaluation of learned reconstructions, including hallucinated structure.

Scanner Geometries

Hounsfield (1979) described the three scanner types in use at the time, all of which produce roughly the same set of readings as parallel sets or fans at many angles:

• Translate–rotate. The source and a small set of detectors traverse the patient, collecting parallel rays, then rotate a step and repeat. His example had 30 detectors and took 18 s per slice.
• Rotate-only fan beam. A fan of X-rays falls on an arc of about 300–500 detectors, and the tube and detectors rotate together around the patient, taking a slice in about 3 s. Kak and Slaney (1988, §4.1.5) call this the rotate–rotate or third-generation geometry.
• Stationary detector ring. About 700–1000 detectors form a fixed ring and only the tube rotates, also in about 3 s. Kak and Slaney call this the rotate–fixed or fourth-generation geometry.

The geometry determines which faults show up as which artifacts. In a third-generation scanner, a miscalibrated detector corrupts the same ray in every view, and the correlated errors back-project into a ring; in a fourth-generation scanner, a bad detector costs one projection, which barely shows, but source instabilities produce rings instead (Kak & Slaney, 1988, §4.1.5). Fourth-generation detectors must accept photons from many directions, so they cannot be collimated against scatter as well.

Two later changes made CT volumetric. In helical (spiral) CT, the table moves continuously while a slip-ring gantry rotates, so the source traces a helix relative to the patient; the American Association of Physicists in Medicine (AAPM) dates its introduction to the early 1990s (AAPM, 2008). As of 2023 the FDA describes a single rotation as typically taking a second or less, with the detectors arranged in rows (FDA, page dated 2023). Multiple-detector-row CT (MDCT) adds detector rows along the table axis $z$: the AAPM report dates dual-slice systems to 1994, four-channel systems to late 1998, 16-channel systems to the end of 2001 and 64-channel systems to mid-2004, and states that all MDCT scanners of its time used the third-generation geometry with rotation times as short as 0.33 s (AAPM, 2008). The table advance is described by the pitch, the table travel per rotation divided by the total nominal beam collimation $N \times T$ for $N$ channels of width $T$. Pitch 1 means contiguous beams, below 1 overlap and above 1 gaps; the report records a period when some manufacturers divided by a single channel width instead, giving values $N$ times larger, until the IEC definition was re-established (AAPM, 2008).

With many detector rows, the rays are no longer confined to parallel planes but form a cone. Cone-beam reconstruction treats each row as a fan tilted out of the plane of rotation. Kak and Slaney (1988, §3.6.2) present a cone-beam filtered back-projection algorithm for a circular orbit based on the analysis of Feldkamp, Davis and Kress (1984). In the plane of rotation it is identical to the fan-beam algorithm, and away from that plane their reconstructions degrade. Only points illuminated from all directions, a sphere around the centre whose radius depends on the cone angle, can be reconstructed properly.

Data It Produces

A reconstructed CT slice is a 2D array of calibrated values, and a scan is a stack of slices forming a 3D volume. Kak and Slaney (1988) mention reconstruction matrices of 256 × 256, 320 × 320 or 512 × 512 for the scanners of their time. Each voxel has a physical size: in-plane pixel spacing set by the field of view and matrix, and slice thickness and spacing set by the acquisition and reconstruction. Slices are often thicker than the in-plane spacing, and the resulting anisotropic voxels matter for any 3D filter, network or measurement. With MDCT, detector channels can be combined after acquisition: the AAPM report’s example presents one 4 × 1.25 mm acquisition as 1.25, 2.5 or 5 mm slices (AAPM, 2008).

CT images are stored in DICOM. The stored integers are mapped to output units by a linear rescale, output $= m \cdot \text{SV} + b$, with the Rescale Slope $m$ and Rescale Intercept $b$ given in the header. For original CT images that are not localizers or multi-energy acquisitions, the standard requires the output units to be Hounsfield units (NEMA, 2026, C.8.2.1). Vision code that reads raw pixel arrays without applying the rescale works with arbitrary offsets rather than HU.

Hounsfield Units

A scanner reconstructs a linear attenuation coefficient, but what it reports is a CT number in Hounsfield units (HU), defined relative to water (Kak & Slaney, 1988, ch. 4, Eq. 12):

$$H = \frac{\mu - \mu_\text{water}}{\mu_\text{water}} \times 1000 .$$

Water is 0 HU by definition, and air, whose attenuation is taken as $\mu = 0$, is −1000 HU. Each HU is 0.1% of water’s attenuation. Kak and Slaney give the typical range of reported values as −1000 to 3000. Hounsfield’s lecture describes the same scale: air at −1000 at the bottom, water chosen to be zero because its absorption is close to that of tissue, and about 1000 levels either side of water, with bone at the top (Hounsfield, 1979).

HU are calibrated but not absolute physical constants. Since $\mu$ depends on photon energy and the beam is polychromatic, a scanner effectively measures attenuation at an effective energy, the monochromatic energy at which a material would show the measured coefficient. Kak and Slaney cite McCullough et al.’s finding that the original EMI head scanner had an effective energy of 72 keV at a tube voltage of 120 kV, and stress that the concept holds only when all rays have the same exit spectrum (Kak & Slaney, 1988, §4.1.2). The same tissue can therefore have different CT numbers on different scanners or tube voltages. A real scanner rarely gives exactly 0 for water and −1000 for air; if the response is linear, measured values $H_m$ can be corrected with the measured water and air values (Kak & Slaney, 1988, ch. 4, Eq. 13):

$$H = \frac{H_m - H_{m,\text{water}}}{H_{m,\text{water}} - H_{m,\text{air}}} \times 1000 .$$

Energy dependence can also be exploited. Kak and Slaney describe how scans at two tube spectra separate the photoelectric and Compton contributions to attenuation, after Alvarez and Macovski and after Brooks, which gives some chemical identification of materials and removes beam-hardening artifacts in principle (Kak & Slaney, 1988, §4.1.2–4.1.3). This is the basis of dual-energy CT.

Windowing

A CT image spans thousands of HU, while a display shows a few hundred grey levels. Hounsfield (1979) described the solution in his lecture: the range of displayed tones is restricted to a small part of the scale, a “window” that can be raised for tissue such as the heart or lowered for detail in the lungs, and narrowing it increases the visible difference between organs such as the liver and its neighbours. A window is specified by its centre (level) $c$ and width $w$. DICOM defines the default linear mapping from input value $x$ to output $y$ in $[y_\text{min}, y_\text{max}]$ (NEMA, 2026, C.11.2.1.2.1):

$$y = \begin{cases} y_\text{min} & x \le c - 0.5 - \frac{w-1}{2}, \\[4pt] y_\text{max} & x > c - 0.5 + \frac{w-1}{2}, \\[4pt] \left(\dfrac{x - (c - 0.5)}{w - 1} + 0.5\right)(y_\text{max} - y_\text{min}) + y_\text{min} & \text{otherwise.} \end{cases}$$

Windowing is a display mapping, not a change to the data: within the window it multiplies differences, contrast and noise alike, by the same factor, and outside it clips. In vision pipelines the same operation can appear as preprocessing, when HU are clipped to a range and rescaled before entering a network.

Practical Example

The examples simulate a parallel-beam scanner with scikit-image 0.25.2, NumPy 2.3.5 and SciPy 1.16.3, using the Shepp–Logan phantom bundled with scikit-image (skimage.data.shepp_logan_phantom, loaded from a file inside the package). The function’s docstring cites Shepp and Logan’s 1974 paper and does not say which variant the image is. The original phantom, as Kak and Slaney (1988, Table 3.1) tabulate it, has ellipse values of 2.0, −0.98, −0.02 and 0.01, which give a skull of 2.0 and nearly uniform interior values of 1.00–1.04. The bundled image’s grey levels are 0, 25, 51, 76, 102 and 255 out of 255. These match, to 8-bit rounding, the amplitudes 1.0, −0.8, −0.2 and 0.1 of the higher-contrast version that Toft (1996, Table B.1) uses and describes as chosen for better contrast, and that Adler and Öktem (2018) call the modified Shepp–Logan phantom.

The phantom’s values have no physical unit, and they are not Hounsfield units. To simulate photon counting, the code adopts an illustrative scale of its own: pixel values are read as attenuation coefficients in cm⁻¹ over a 25.6 cm field, so that the value 0.2 plays the role of water. The first script prints the phantom’s grey levels, varies the number of views on noiseless data, then adds a small dense insert, a disc 0.77 cm across with $\mu = 10$ cm⁻¹ (a value of our choosing, not a model of a particular metal), and compares the purely linear data with photon counts, Poisson counts from $10^5$ incident photons per ray followed by the logarithm, for three noise seeds. Reconstruction uses iradon with its default ramp filter.

import numpy as np
from skimage.data import shepp_logan_phantom
from skimage.transform import radon, iradon

mu = shepp_logan_phantom()                 # 400 x 400, values in [0, 1]
print("grey levels x 255:", np.unique(np.round(mu * 255)).astype(int))

# Illustrative physical scale (our assumption, not part of the phantom):
# read pixel values as attenuation in 1/cm and give the 400-pixel field 25.6 cm.
dx = 25.6 / 400                            # cm per pixel; 0.2 plays the role of water
n = mu.shape[0]
yy, xx = np.mgrid[:n, :n]
flat = (slice(276, 316), slice(238, 278))  # 40 x 40 patch of value 0.2

def fbp(sino, theta):
    return iradon(sino, theta=theta, filter_name="ramp") / dx

# 1. Fewer views, noiseless data: streaks from angular undersampling.
print("views   SD in flat patch (noiseless)")
for k in [720, 360, 90, 30]:
    theta = np.linspace(0.0, 180.0, k, endpoint=False)
    rec = fbp(radon(mu, theta=theta) * dx, theta)
    print("%5d   %.4f" % (k, rec[flat].std()))

# 2. A dense insert: a 0.77 cm disc with mu = 10 /cm, in 360 views.
theta = np.linspace(0.0, 180.0, 360, endpoint=False)
dense = mu.copy()
dense[(yy - 200) ** 2 + (xx - 200) ** 2 <= 6**2] = 10.0
N0 = 1e5
for label, obj in [("no insert", mu), ("dense insert", dense)]:
    p = radon(obj, theta=theta) * dx
    sds = []
    for seed in [0, 1, 2]:
        counts = np.random.default_rng(seed).poisson(N0 * np.exp(-p))
        rec = fbp(-np.log(np.maximum(counts, 1) / N0), theta)
        sds.append(rec[flat].std())
    print("%-12s  linear data: SD %.4f   photon counts: SD %.4f-%.4f "
          "(seeds 0-2), rays with < 10 photons: %d"
          % (label, fbp(p, theta)[flat].std(), min(sds), max(sds),
             (N0 * np.exp(-p) < 10).sum()))

grey levels x 255: [  0  25  51  76 102 255]
views   SD in flat patch (noiseless)
  720   0.0043
  360   0.0083
   90   0.0335
   30   0.0689
no insert     linear data: SD 0.0083   photon counts: SD 0.0149-0.0154 (seeds 0-2), rays with < 10 photons: 0
dense insert  linear data: SD 0.0087   photon counts: SD 0.0345-0.0378 (seeds 0-2), rays with < 10 photons: 2857

With 400 detector bins, Kak and Slaney’s rule of about $\pi/2$ views per detector bin asks for about 630 views. The variation in a region that should be uniform roughly doubles from 720 to 360 views and keeps growing as views are removed, with no noise in the data at all: these are the streaks of angular undersampling. With linear, noiseless data the dense insert barely changes the uniform patch, since filtered back-projection is linear and the insert’s contribution is reconstructed in place. With photon counts, the rays through the insert receive almost no photons, their logarithms are dominated by Poisson noise and by the floor of one count that the code imposes, and back-projection spreads those errors along the rays as streaks across the slice. Over the three seeds the patch’s variation is more than twice that of the scan without the insert. Photon starvation is only one of the mechanisms behind streaks from dense objects; this simulation leaves out beam hardening and scatter.

The second script adds a faint disc to the phantom, 0.004 above its surroundings, reconstructs it with the Hann filter (the ramp multiplied by a Hann window), and converts the result to illustrative Hounsfield units by treating the value 0.2 as water and 0 as air. Under that mapping the disc is +20 HU. It then displays the image through two windows with the DICOM formula.

import numpy as np
from skimage.data import shepp_logan_phantom
from skimage.transform import radon, iradon

dx, mu_water, N0 = 25.6 / 400, 0.2, 1e5   # illustrative scale, as above
mu = shepp_logan_phantom()
n = mu.shape[0]
yy, xx = np.mgrid[:n, :n]
r2 = (yy - 285) ** 2 + (xx - 125) ** 2
lesion, ring = r2 <= 12**2, (r2 >= 18**2) & (r2 <= 28**2)
mu = mu + 0.004 * lesion                  # faint disc: +2% of "water"

theta = np.linspace(0.0, 180.0, 360, endpoint=False)
p = radon(mu, theta=theta) * dx
counts = np.random.default_rng(0).poisson(N0 * np.exp(-p))
rec = iradon(-np.log(np.maximum(counts, 1) / N0), theta=theta,
             filter_name="hann") / dx

# Illustrative HU: phantom value 0.2 is treated as water, 0 as air.
hu = 1000.0 * (rec - mu_water) / mu_water
print("HU: background %.1f, lesion %.1f, background SD %.1f"
      % (hu[ring].mean(), hu[lesion].mean(), hu[ring].std()))

def window(x, c, w, ymin=0.0, ymax=255.0):
    # DICOM PS3.3 C.11.2.1.2.1, default LINEAR function
    y = ((x - (c - 0.5)) / (w - 1) + 0.5) * (ymax - ymin) + ymin
    y = np.where(x <= c - 0.5 - (w - 1) / 2, ymin, y)
    return np.where(x > c - 0.5 + (w - 1) / 2, ymax, y)

head = mu > 0
print("window (C, W)    contrast  noise SD  contrast/noise  head pixels clipped")
for c, w in [(1500, 6000), (0, 100)]:
    g = window(hu, c, w)
    con, sd = g[lesion].mean() - g[ring].mean(), g[ring].std()
    clipped = np.mean((g[head] == 0) | (g[head] == 255))
    print("(%4d, %4d)  %10.2f %9.2f %15.2f %19.1f%%"
          % (c, w, con, sd, con / sd, 100 * clipped))

HU: background -0.4, lesion 18.3, background SD 23.6
window (C, W)    contrast  noise SD  contrast/noise  head pixels clipped
(1500, 6000)        0.80      1.00            0.79                 0.0%
(   0,  100)       46.62     59.39            0.78                30.9%

The reconstruction recovers the disc’s offset to within a couple of HU, though its pixel-level noise is larger than the offset; the disc stands out on average over its few hundred pixels. The wide window shows the whole head without clipping, but compresses the disc to less than one grey level. The narrow window spreads the same 20 HU over dozens of grey levels, and it spreads the noise by the same factor, so the ratio of contrast to noise is unchanged; the price is that nearly a third of the head, everything outside −50 to +50 HU, is clipped to black or white.

Artifacts

CT artifacts are systematic errors produced when the data violate the assumptions of the reconstruction. They matter to vision models because they are structured, often oriented along lines through the slice, and can resemble anatomy.

• Beam hardening. Because the beam hardens along its path, rays through a lot of bone emerge with a higher mean energy and appear less attenuating than the linear model expects. Kak and Slaney (1988, §4.1.3) show two effects in simulation: CT numbers raised near the skull, which in the early years was mistaken for grey matter of the cortex, and wide dark streaks between dense bones. Corrections linearize the measurements using calibration on a water-like absorber, re-project a thresholded bone image to correct the rays through bone, or use dual-energy data.
• Scatter. Scattered photons add a roughly constant intensity to the detectors. Where the primary beam is strong this hardly matters, but where it is weak the relative error is large, and because the error depends on the view it reconstructs as streaks (Kak & Slaney, 1988, §4.1.4).
• Undersampling. Too few views produce thin streaks radiating from high-contrast structures, as in the example; undersampling within each projection, where sharp edges exceed the detector sampling rate, produces streaks too (Kak & Slaney, 1988, §5.1).
• Rings. A miscalibrated detector in a rotate–rotate scanner produces concentric rings (Kak & Slaney, 1988, §4.1.5).
• Motion. Data for one slice take time to collect, and an object that moves meanwhile is inconsistent between views. Kak and Slaney describe the blurring of the heart and gated CT, which keeps only the projections taken at the same phase of the cardiac cycle recorded by an electrocardiogram (§4.1.5). Faster rotation and wider detectors in later scanners shorten the acquisition (AAPM, 2008).
• Metal. Cormack (1979) noted that a piece of metal in the body distorts a CT slice containing it beyond recognition, and that sharp attenuation changes at bone–soft-tissue interfaces produce overshoots. It follows from the mechanisms above that metal, far denser than bone, exaggerates all of them, adding photon starvation on the rays through it; the example shows that part alone.
• Partial volume. A detector element and a voxel have finite size, so a voxel straddling two tissues reports a mixture. The ideal line-integral model assumes a beam narrow enough to avoid this (Kak & Slaney, 1988, §4.1.1).

Dose

CT uses ionizing radiation, and dose constrains how it is designed and used. Image noise comes from the limited number of photons reaching the detectors; Hounsfield (1979) noted that improvements in picture grain would be proportional to dose for a given resolution, and that industrial uses have no such constraint. The trade-offs between dose, noise and resolution are the reason for the low-dose and sparse-view research described in medical image reconstruction: Jin et al. (2017) point out that reducing the number of views by a factor of twenty reduces the dose by the same factor.

Scanners report standardized dose indices, which the AAPM report defines (AAPM, 2008):

• CTDI$_\text{vol}$ (volume CT dose index, mGy) is the weighted CTDI measured in a standard cylindrical phantom divided by the pitch. It represents the average dose within the scanned volume for that phantom, not the dose to a particular patient, and does not depend on scan length.
• DLP (dose–length product, mGy·cm) is CTDI$_\text{vol}$ times the scan length, which reflects the total energy delivered by the scan.
• Effective dose (mSv) weights organ doses by their radiation sensitivity, to express a non-uniform exposure as an equivalent whole-body one. The report stresses that it describes a standard model body, not any individual, and should be quoted with no more than one or two significant figures.

The U.S. Food and Drug Administration states that effective doses from diagnostic CT procedures are typically estimated at 1 to 10 mSv, and that estimates for a procedure can vary by a factor of 10 or more with the type of procedure, patient size, and the CT system and its operating technique (FDA, page dated 2017). The FDA recommends that imaging teams use the lowest dose that yields image quality adequate for diagnosis and intervention (FDA, page dated 2023).

Processing the Data

Several properties of CT data shape the vision methods built on it.

• Use HU, not raw values. After the DICOM rescale, intensity-based steps such as clipping and thresholding transfer between scans better than they do between photographs, within the limits set by effective energy and scanner calibration.
• Respect the physical spacing. Resample or account for anisotropic voxels before 3D processing and when measuring sizes and distances. Surface meshes of bone or organs are commonly extracted from CT volumes with marching cubes, which was developed for medical slice data.
• Treat noise as structured. CT noise concentrates at high frequencies (Hounsfield, 1979), its level follows the reconstruction filter, and it follows from the photon statistics that it varies with position, since rays through more material carry fewer photons. Assumptions behind generic image noise models and image denoising methods should be checked against this.
• Expect scanner and protocol shift. The reconstruction filter, slice thickness, tube voltage (through the effective energy), dose and reconstruction method each change the noise, texture or CT numbers, and with them the inputs a model sees.
• Segment in 3D, in physical units. Delineating organs and lesions in CT volumes is the subject of medical image segmentation; the calibrated intensities and known spacing are what distinguish it from segmenting photographs.
• Registration with other modalities. CT intensities are not a function of MR intensities, which is why multimodal registration uses statistical similarity measures such as mutual information.

Strengths

• Cross-sectional images without superimposition of structures at different depths, with soft-tissue contrast that projection radiography lacks (Hounsfield, 1979).
• Calibrated values in Hounsfield units.
• Fast volumetric acquisition with helical and multi-row scanners; Hounsfield’s lecture already described CT’s use in radiotherapy planning.
• A well-understood forward model, which makes physics-based simulation, as in the example, a practical source of training and test data.

Limitations

• Ionizing radiation, so dose limits the number of photons and with it noise and resolution.
• Soft tissues differ little in attenuation, and a narrow window makes the differences visible only as far as noise allows. Intravenous contrast agents are used to improve visualization, with their own risk of reactions (FDA, page dated 2023).
• Artifacts from beam hardening, scatter, metal, motion and undersampling.
• CT numbers that depend on the X-ray spectrum and the scanner, so “calibrated” does not mean identical across sites.

History

The mathematics came first and unnoticed. Radon solved the reconstruction of a function from its line integrals in 1917; the problem was solved again in statistics by Cramér and Wold in 1936, in radio astronomy by Bracewell and in electron microscopy by DeRosier and Klug in 1968, as Cormack recounted in his Nobel lecture (Cormack, 1979). Cormack, a physicist in Cape Town, came to the problem in 1956 while supervising the use of radioisotopes at Groote Schuur Hospital, where he saw that radiotherapy planning assumed a homogeneous body. He worked out a solution over the following years, tested it in 1963 on a phantom of aluminium and Lucite scanned with a cobalt-60 gamma-ray beam, and published in 1963 and 1964 to virtually no response.

Hounsfield, an engineer at EMI’s Central Research Laboratories, had the idea for what became the EMI scanner in 1967, while exploring automatic pattern recognition (NobelPrize.org, Hounsfield biographical). His laboratory rig moved a gamma-ray source across the object on a lathe bed, rotating the object 1° after each sweep, and took nine days per picture; an X-ray tube cut this to nine hours, and the rig imaged preserved and fresh brains and pig sections. A clinical prototype was installed at Atkinson Morley’s Hospital in London, where the first patient was scanned in 1972 and the image showed a cyst in the brain (Hounsfield, 1979). Body scanners followed, and by 1979 Hounsfield put the number of hospitals operating CT at about 1000. The 1979 Nobel Prize in Physiology or Medicine was awarded jointly to Allan M. Cormack and Godfrey N. Hounsfield “for the development of computer assisted tomography” (NobelPrize.org, 1979).

By 1988 Kak and Slaney could describe almost all production scanners as fan-beam machines reconstructed by filtered back-projection; helical and multi-row scanning followed, as dated under Scanner Geometries.

Comparison with Other Sensors

• Projection radiography uses the same X-ray physics but records one projection, so depth is superimposed; CT trades a higher dose for cross-sections and calibrated values (the FDA’s table gives typical effective doses of 0.02 mSv for a chest radiograph and 7 mSv for chest CT; FDA, page dated 2017).
• Magnetic resonance imaging, which Kak and Slaney also treat as a tomographic modality (§4.4), forms images from magnetic resonance signals rather than X-ray transmission, so its contrast reflects different tissue properties.
• Positron emission tomography and single-photon emission CT reconstruct the distribution of a radioactive tracer from the photons it emits, with the same family of algorithms after correcting for attenuation (Kak & Slaney, 1988, §4.2).
• Industrial CT applies the same principles to objects such as rocket motors and oil-well cores (Kak & Slaney, 1988, §4.1.6), free of the dose limits of medical CT.
• Cameras and LiDAR measure light reflected from surfaces, so they see only the outside of objects; CT measures transmission and sees through them.