🎯 Radiotherapy Treatment Planning System Dose Calculation

From CT Numbers to Electron Density

Every modern photon dose calculation begins with a patient-specific 3D density map. A planning CT scan assigns each voxel a Hounsfield Unit (HU) describing its X-ray attenuation relative to water, and a calibration curve converts that number into the relative electron density the dose engine actually needs to model radiation transport.

• 1–3 mm: CT slice thickness (typical simulation protocol)
• -1000 to +3000: Hounsfield range (air to dense cortical bone)
• CIRS / Gammex: Calibration phantom (tissue-equivalent plug set)
• 2–5 million: Typical voxel count (per patient CT dataset)

Hounsfield units and the CT calibration curve

The Hounsfield scale is defined so that water reads 0 HU and air reads -1000 HU, with denser materials (cortical bone) reaching +1000 to +3000 HU depending on scanner kVp. HU by itself is not directly usable by a dose engine — it depends on beam energy and scanner model, not purely on physical density.\n\nAt commissioning, physicists scan a phantom containing tissue-equivalent plugs of known physical and electron density (lung, adipose, muscle, trabecular and cortical bone) to build a scanner-specific, piecewise-linear HU-to-relative-electron-density curve. Every planning CT voxel is then mapped through this curve, giving the dose algorithm the electron density it needs to compute attenuation and scatter correctly.\n\nA more rigorous variant, the stoichiometric calibration method, additionally derives elemental composition per HU bin so that Compton and photoelectric interaction probabilities are represented more faithfully across the full range of tissues encountered clinically.

A single miscalibrated CT-to-density curve can shift computed dose by several percent across an entire treatment course, which is why AAPM TG-66 requires an independent HU-to-density calibration check at commissioning and at least annually thereafter.

Structure delineation and density segmentation

Once the density grid exists, the target volumes and organs at risk are contoured on the same CT dataset — the same voxel grid that the dose engine will later trace rays through. Contours that straddle sharp density boundaries (a lung tumour abutting aerated lung, a spine target beside cortical bone) are where dose calculation accuracy matters most.\n\nMost engines internally classify tissue into broad density classes — air/lung (density below roughly 0.3), soft tissue (around 1.0), and bone (above roughly 1.3) — because interaction physics and the width of the scatter kernel differ qualitatively between these regimes, not just quantitatively.\n\nOverriding CT density inside high-Z artifacts (dental fillings, hip prostheses, fiducial markers) is standard practice, since metal streak artifacts corrupt the HU-to-density mapping locally and can otherwise introduce large, spurious dose errors right where accuracy is needed most.

Voxel grid resolution trade-offs

The image (CT) grid and the dose calculation grid are usually distinct: the CT is acquired at sub-millimetre in-plane resolution, but the dose is calculated on a coarser grid to keep memory and computation time manageable. A finer grid resolves steep dose gradients more accurately but increases voxel count — and hence computation time — roughly with the cube of the resolution improvement.\n\nClinically, 2–3 mm dose grids are standard for conventional fractionation, while stereotactic body radiotherapy (SBRT) plans — with their small fields and very steep dose falloffs — commonly drop to 1–1.5 mm to avoid smearing out the high-gradient penumbra that defines a conformal SBRT dose cloud.

Ray Tracing and the TERMA Calculation

Before any dose can be computed, the planning system needs to know how much energy the primary photon beam releases at every point in the patient. Ray tracing along divergent beam lines through the electron density grid, attenuated according to the Beer-Lambert law, yields the TERMA — Total Energy Released per unit MAss — at each voxel.

• 6–15 MV: Linac photon energy (typical clinical beams)
• I = I₀e^-μx: Attenuation law (exponential fluence decay)
• Σρᵢ·dᵢ: Radiological path length (density-weighted depth)
• Batho, ETAR: Correction-era algorithms (1980s heterogeneity corrections)

Effective depth and the radiological path length

Ray tracing does not measure how many centimetres of geometric distance a ray has travelled — it measures radiological (water-equivalent) depth: the density-weighted line integral Σρᵢ·dᵢ along the ray. A beam crossing 3 cm of aerated lung (relative density ≈ 0.25) has effectively only traversed the equivalent of about 0.75 cm of water, while the same geometric distance through cortical bone (density ≈ 1.7–1.9) is radiologically much deeper.\n\nThis single reinterpretation of “depth” is what allows a dose engine trained on homogeneous-water beam data to extend its primary-beam attenuation model to a real, heterogeneous patient.

From correction-based to model-based algorithms

Early treatment planning systems (1980s–1990s) used correction-based algorithms: compute the dose in an equivalent homogeneous water phantom, then apply a 1D scaling factor — Batho power-law or Equivalent Tissue-Air Ratio (ETAR) — based only on the radiological depth along the central ray. These were cheap to compute but fundamentally blind to lateral scatter: they could not represent how a density change beside the ray, not just along it, redistributes dose.\n\nModel-based algorithms — convolution/superposition and Monte Carlo — replaced this by explicitly transporting energy in three dimensions, which is essential wherever density varies sharply in the plane perpendicular to the beam, such as at a lung/tumour or lung/chest-wall interface.

Correction-based algorithms can be accurate to within 2–3% in a homogeneous water-equivalent patient, yet drift to 10% or more locally near a lung interface, because they only ever correct depth dose along the central ray and never model lateral disequilibrium.

TERMA is not yet dose

TERMA marks where photons interact and release energy — but that is not where the energy is ultimately absorbed. The secondary (Compton) electrons produced at each interaction site carry that energy outward through the surrounding tissue before finally depositing it, travelling anywhere from a fraction of a millimetre in dense bone to a centimetre or more in low-density lung at megavoltage energies.\n\nConverting TERMA into absorbed dose — accounting for exactly this electron transport — is the job of the convolution/superposition step covered next.

Convolution / Superposition — From TERMA to Absorbed Dose

Secondary electrons and scattered photons released at a primary interaction site do not deposit their energy on the spot — they travel outward through the surrounding tissue before losing all of their kinetic energy. Convolution/superposition algorithms capture this by spreading a pre-computed energy deposition kernel around every interaction point and summing the results.

• Monte Carlo: Kernel source (pre-generated point kernels)
• Mackie, 1985: Kernel formulation (polyenergetic point-spread)
• several cm: Kernel lateral extent (electron + photon scatter range)
• O'Connor theorem: Density scaling rule (radiological kernel scaling)

The point dose deposition kernel

A dose deposition kernel describes, in three dimensions, what fraction of the energy released at a single interaction point is absorbed at every surrounding location. These kernels are pre-computed offline with Monte Carlo simulation of a mono-energetic photon interacting in a homogeneous water medium, then summed (weighted by the beam's energy spectrum) into a polyenergetic kernel for clinical use.\n\nThe kernel has a characteristic forward-peaked, teardrop shape — most energy is carried a short distance forward along the original photon direction by Compton electrons, with a longer, weaker tail of multiply-scattered photon dose extending outward in every direction.\n\nThe absorbed dose at any point is then the convolution integral D(r) = ∫ TERMA(r′)·K(r − r′) dr′ — summing every interaction site's kernel contribution that reaches that point. This superposition of many spread-out kernels is exactly what stage 3's animated pulses along the ray paths represent.

Collapsed-cone acceleration

A literal 3D convolution over every voxel pair is computationally expensive — effectively an operation that scales with the square of the voxel count. The collapsed-cone convolution (CCC) method, introduced by Ahnesjö in 1989, discretizes the kernel into a modest number of fixed solid-angle directions ("cones") radiating from each interaction point, and transports energy along those discrete rays instead of to every individual voxel.\n\nThis collapses a full solid-angle integral into perhaps a few dozen directions per point, cutting computation time by roughly an order of magnitude while retaining most of the three-dimensional scatter accuracy — the reason collapsed-cone convolution/superposition became the clinical workhorse algorithm from the 1990s onward.

O'Connor's density scaling theorem states that a kernel computed in water can be rescaled for a different density medium simply by shrinking or stretching its physical dimensions inversely with local density — the same principle animated here as narrower pulses in bone and wider pulses in lung.

Why heterogeneities challenge the kernel model

Density scaling works well when density changes gradually, but it is only an approximation across a sharp interface. Immediately beyond a soft-tissue-to-lung boundary, the medium has not yet accumulated enough depth for secondary electrons to reach charged-particle (electronic) equilibrium — the kernel's implicit assumption of local equilibrium briefly breaks down, and computed dose can be measurably wrong in a thin region right at the boundary.\n\nThis effect is strongest for small fields and higher photon energies, which is exactly the combination used in lung SBRT — making interface accuracy a first-order clinical concern rather than a theoretical curiosity.

Summing Beams into the Final 3D Dose Distribution

A modern treatment plan rarely uses a single beam. IMRT and VMAT plans deliver dose from many gantry angles, or from a continuously rotating arc, and each contributes its own convolved dose distribution. The planning system sums every beam's contribution voxel-by-voxel across the full 3D grid to build the total dose the patient will actually receive.

• 5–9 fields: Typical beam count (step-and-shoot IMRT)
• 1–2 arcs: VMAT arcs (continuous gantry rotation)
• 2–3 mm: Dose grid resolution (common clinical setting)
• 10% increments: Isodose line spacing (of prescription dose)

Beam summation and plan normalization

Each field or arc is calculated as its own 3D dose cube, segment by segment for step-and-shoot IMRT or control point by control point for VMAT, then combined with the monitor-unit weighting the optimizer assigned to each. The summed cube is finally normalized so that a chosen prescription metric — commonly that 95% of the target volume receives 100% of the prescribed dose — is satisfied exactly, which is what allows two very differently shaped plans to be compared on a common percentage scale.

Isodose lines and the dose-volume histogram

Isodose lines connect points of equal dose, conventionally drawn every 10% of the prescription dose (90%, 80%, 70%…), giving planners an immediate visual read of conformality and hot or cold spots. The same voxel dose grid is also binned per anatomical structure into a dose-volume histogram (DVH), which compresses the entire 3D distribution into a single 2D plot used for formal plan approval against clinical dose constraints.\n\nBoth tools depend entirely on the accuracy of the underlying voxel-by-voxel accumulation performed in this stage — any local error upstream propagates directly into the isodose shape and the DVH curve a physician signs off on.

Modern IMRT and SBRT plans are defined by steep dose gradients — sometimes falling from 100% to 20% of prescription within a few millimetres — so accurate 3D accumulation at fine grid resolution is what actually delivers that conformality in the treated patient, not just on the screen.

Grid resolution and calculation cost

Because voxel count scales with the cube of linear resolution, halving the dose grid spacing roughly multiplies computation time and memory by a factor of eight. Clinical practice therefore settles on the coarsest grid that still resolves clinically relevant gradients — 2–3 mm for the majority of conventionally fractionated plans, tightened to 1–1.5 mm for small, high-gradient SBRT targets where coarser grids would visibly round off the true dose falloff.

Algorithm Choice — Speed versus Accuracy

No single dose algorithm is free of trade-offs. Pencil Beam models are fast but can misjudge lateral scatter near density interfaces; collapsed-cone convolution/superposition balances speed and 3D accuracy for routine use; Monte Carlo simulates individual particle histories and is treated as the accuracy benchmark against which every other algorithm is validated.

• 2–3%: AAPM TG-53 tolerance (homogeneous-medium target)
• 3%/3mm, 2%/2mm: Clinical gamma criterion (dose/distance agreement)
• up to 10–15%: Pencil Beam lung-interface error (local dose, small SBRT fields)
• <60 s: GPU Monte Carlo run time (modern commercial engines)

Where simple algorithms break down

Pencil Beam convolution models the lateral spread of dose with a fixed kernel shape calibrated in water, then only corrects the depth-dose curve for the radiological path length — it never widens or narrows that lateral kernel to reflect the local density it is actually passing through. At a lung/soft-tissue or lung/bone interface this can over- or under-estimate local dose by several percent, occasionally exceeding 10–15% for the small fields typical of lung SBRT, where loss of lateral electronic equilibrium is most pronounced.\n\nThe AAPM TG-101 report on SBRT explicitly warns against relying on correction-based or simple pencil-beam algorithms for lung targets, recommending convolution/superposition or Monte Carlo instead.

AAPM Task Group 101 (SBRT) specifically flags lung tumours as a scenario where algorithm choice alone can change reported target coverage by a clinically meaningful margin — the same plan can look adequate under Pencil Beam and inadequate under Monte Carlo.

Monte Carlo as the gold standard

Monte Carlo dose calculation simulates individual photon and electron histories directly through the actual voxelized density grid, using the real physical interaction cross-sections at every step rather than a pre-scaled kernel. As the number of simulated histories grows, the result converges to the true answer, with statistical uncertainty falling proportionally to 1/√N.\n\nHistorically this made Monte Carlo impractical for daily clinical use — a fully converged calculation could take hours on a single CPU core. GPU-accelerated implementations, now standard in several commercial planning systems, have cut that to well under a minute for a full clinical plan, making Monte Carlo a practical first-line algorithm rather than only an offline verification tool.

Commissioning, tolerances, and independent verification

AAPM Task Group 53 (1998) established the framework clinics still use to commission a dose algorithm: benchmark it against measurement in homogeneous phantoms (target agreement roughly 2–3%) and against heterogeneous phantoms or independent Monte Carlo, where wider tolerances are accepted precisely because the physics is harder.\n\nBeyond initial commissioning, most departments run an independent second-check calculation (an independent point-dose or monitor-unit verification using a different algorithm) for every patient plan, and participate in periodic external audits — such as those run by IROC — that irradiate a phantom under the clinic's own plan and compare measured to computed dose, catching systematic algorithm or calibration errors before they ever reach a patient.