Radiometry and Photometry

Radiometry measures light as electromagnetic radiation, in physical units: joules, watts, square metres and steradians. Photometry measures the same light weighted by the spectral sensitivity of a standard human observer, in lumens, candelas and lux. Both describe how much light leaves a source, how it spreads through space, how much falls on a surface and how much a surface sends towards a camera. For computer vision the central quantity is radiance, the power per unit projected area per unit solid angle: it is what a surface “looks like” from a direction, it does not change along a ray in empty space, and it is what a linear camera pixel measures up to constants set by the lens and the exposure.

Image formation follows the whole chain from light source to stored pixel and uses radiance, irradiance and the BRDF without defining them in detail. This article defines the quantities and the geometry behind them: solid angle and projected solid angle, the ladder from flux to radiance, the cosine and inverse-square laws and when the inverse-square law fails, Lambertian surfaces, the reason a matte wall is equally bright from near and far, Planck’s law as a reference spectrum, and the photometric units built on the SI definition of the candela. Its example computes the irradiance of a finite disc against the point-source laws, renders a floor and a patch with a pinhole camera to show radiance invariance, and converts two spectra to lumens and lux.

Definition

Radiometry and photometry share one set of geometric quantities; photometry weights each by the luminous efficiency of the eye, $V(\lambda)$, and a constant.

Radiometric quantity	Symbol	SI unit	Photometric counterpart	Unit
Radiant flux (power)	$\Phi$	W	Luminous flux	lm
Radiant intensity	$I = d\Phi/d\omega$	W sr⁻¹	Luminous intensity	cd = lm sr⁻¹
Irradiance (incident), exitance (leaving)	$E = d\Phi/dA$, $M$	W m⁻²	Illuminance, luminous exitance	lx = lm m⁻²
Radiance	$L = d^2\Phi/(dA\cos\theta\, d\omega)$	W m⁻² sr⁻¹	Luminance	cd m⁻²

The units of radiant intensity and radiance are those of the SI Brochure (BIPM, 2019, Table 6, p. 136), which notes that the steradian is usually kept in radiometric units (Table 6, note a), where it keeps units of flux and of intensity apart (p. 133); the lumen (cd sr) and lux (lm/m²) are in its Table 4 (pp. 133–134). The pairings in the last two columns follow ISETCam’s ieLuminanceFromEnergy.m, which states that one formula turns spectral radiance into luminance (cd/m²), spectral irradiance into illuminance (lux), spectral radiant intensity into luminous intensity (cd) and spectral power into luminous flux (lm). Each radiometric quantity also has a spectral version per unit wavelength, such as spectral radiance $L_\lambda$ in W m⁻² sr⁻¹ nm⁻¹; multispectral satellite imaging works with band radiances in W m⁻² sr⁻¹ µm⁻¹.

In vision the word “intensity” usually means a pixel value. In radiometry it means only radiant intensity, power per steradian from a source small enough to treat as a point; this article uses it only in that sense.

Intuition

Think of light as a flow of energy along rays. Flux counts the total energy per second through some surface. Irradiance asks how densely that flux lands on a surface: the same lamp gives a lower irradiance on a desk when it is lifted, because the same flux spreads over a larger area. Radiant intensity asks how much flux a small source sends into a given cone of directions. Radiance asks the finest question: how much power travels along a narrow bundle of rays, per unit cross-section of the bundle and per unit spread of its directions.

Radiance is what makes a surface look bright. A matte wall photographed from 2 m and from 8 m gives the same pixel values: at 8 m each pixel sees a patch of wall sixteen times larger, but the lens collects only a sixteenth as much light from each unit of wall. The wall’s image shrinks; its brightness does not. Only when a source is smaller than a pixel does its image dim with distance, because then all its light lands in one pixel and the inverse-square law takes over. Stars and distant lamps behave like points; walls, faces and skies behave like extended sources.

Photometry answers a different question: how bright the light is to a person. Two lamps of equal radiant power can differ several-fold in lumens, because the eye is far more sensitive in the middle of the visible spectrum than at its ends, and insensitive outside it, where a 2856 K blackbody sends most of its power (the example below computes that 10.5% of it lies between 380 and 780 nm).

Formal Definition

Solid angle and projected solid angle

A solid angle is the area a cone of directions cuts out of a sphere, divided by the squared radius. One steradian subtends an area $r^2$ on a sphere of radius $r$, and the full sphere subtends $4\pi$ sr (BIPM, 2019, Table 4, note c). In spherical coordinates about a surface normal, with polar angle $\theta$ and azimuth $\phi$,

$$d\omega = \sin\theta\, d\theta\, d\phi, \qquad d\Omega = \cos\theta\, d\omega,$$

where $\Omega$ is the projected solid angle (Nicodemus et al., 1977, note to Table 1, p. 11). The hemisphere above a surface has solid angle $2\pi$ sr but projected solid angle $\pi$ sr; the factor $\pi$ that appears throughout radiometry is this projected hemisphere. A small surface of area $dA'$ at distance $r$, tilted by $\theta'$ from the line of sight, subtends $d\omega = dA' \cos\theta' / r^2$ (our statement of the standard geometry).

Radiance and irradiance

Nicodemus et al. (1977, Eq. 4, p. 5) write the flux that a beam of radiance $L_i$ arriving from direction $\theta_i$ within $d\omega_i$ delivers to a surface element as

$$d\Phi_i = L_i \cos\theta_i \, d\omega_i \, dA_i = dE_i \, dA_i .$$

Read backwards, this is the definition of radiance as flux per unit area projected perpendicular to the beam and per unit solid angle, $L = d^2\Phi/(dA\cos\theta\,d\omega)$. Integrating over all incoming directions gives the irradiance,

$$E = \int_{\text{hemisphere}} L_i(\theta, \phi) \cos\theta \, d\omega = \int L_i \, d\Omega .$$

The factor $\cos\theta$ is Lambert’s cosine law for incidence: a beam striking a surface obliquely spreads over a larger area. A uniform sky of radiance $L$ gives $E = \pi L$.

Point sources and the inverse-square law

A point source of radiant intensity $I$ illuminates a surface element at distance $r$, with angle of incidence $\theta$, through the solid angle $dA \cos\theta/r^2$, so (our derivation from the definitions above)

$$E = \frac{I \cos\theta}{r^2}.$$

The $1/r^2$ is the inverse-square law. Levoy (Image Formation, slide 53) gives the conservation-of-energy picture behind it, for light passing through a lens aperture: doubling the distance to the sensor makes the cone’s cross-section four times larger, so the light per unit area falls fourfold. For a point source at height $h$ above a horizontal floor, a floor point at horizontal distance $x$ has $r = h/\cos\theta$, and the irradiance is

$$E(x) = \frac{I h}{(h^2 + x^2)^{3/2}} = \frac{I}{h^2}\cos^3\theta ,$$

one cosine from incidence and two from the growing distance. If the source is itself a small flat Lambertian emitter facing down, its intensity is not isotropic but falls as $I(\theta) = I_0\cos\theta$, because its projected area shrinks with angle; the floor then receives $E = I_0 \cos^4\theta / h^2$. These are not the $\cos^4\alpha$ of lens vignetting, which image formation gives for the light inside a camera, although the factors have the same origin.

Extended sources

The inverse-square law holds exactly only for a point. For a uniform Lambertian disc of radius $R$ and radiance $L$, facing down at height $h$, integrating $L \cos\theta_e \cos\theta_i / r^2$ over the disc gives on the axis

$$E_{\text{disc}}(0) = \frac{\pi L R^2}{R^2 + h^2},$$

and at horizontal offset $x$, with $s = h^2 + x^2$,

$$E_{\text{disc}}(x) = \frac{\pi L}{2}\left(1 - \frac{s - R^2}{\sqrt{(s + R^2)^2 - 4x^2R^2}}\right).$$

Both are our derivation, checked numerically in the example by quadrature and by Monte Carlo integration. Expanding the ratio of the disc to a Lambertian point with the same on-axis intensity (our derivation) gives $E_{\text{disc}}/E_{\text{far}} - 1 = (R/h)^2(2\tan^2\theta - 1)\cos^4\theta + O\big((R/h)^4\big)$, so to first order the point model overestimates on the axis and underestimates beyond $\tan\theta = 1/\sqrt{2}$ (35.264°); with the next term, the sign change moves to $\tan\theta \approx 1/\sqrt{2} + (5\sqrt{2}/36)(R/h)^2$. A point with the same on-axis intensity $I_0 = \pi R^2 L$ predicts $I_0/h^2$, which exceeds the exact value by the factor $(R^2 + h^2)/h^2$. The relative error of the inverse-square law on the axis is therefore exactly $R^2/h^2$: 25% at a distance of one diameter ($h = 2R$), 4% at 2.5 diameters, 1% at 5 diameters. Rules of thumb that the inverse-square law holds beyond some multiple of the source size are this calculation in another form; no source read for this article states one.

Lambertian surfaces

A perfectly diffuse or Lambertian surface reflects with the same radiance in all directions, whatever the illumination. Nicodemus et al. (1977, Appendix C, Eqs. C1–C4, p. 43) show that this requires a constant bidirectional reflectance distribution function (BRDF), $f_{r,d} = L_r/E_i$ in sr⁻¹, and that the fraction of incident flux reflected, the diffuse reflectance $\rho_d$, is $f_{r,d}$ times the projected solid angle of the hemisphere:

$$\rho_d = \pi f_{r,d}, \qquad f_{r,d} = \frac{\rho_d}{\pi}, \qquad L_r = \frac{\rho_d E_i}{\pi}, \qquad M = \pi L_r = \rho_d E_i .$$

The exitance of a Lambertian surface is $\pi$ times its radiance, not $2\pi$ times, because the outgoing radiance is weighted by $\cos\theta$ over the hemisphere. The same monograph uses $1/\pi$ sr⁻¹ as the BRDF of a perfect (lossless) Lambertian reflector when it defines reflectance factors (p. 11). Visionbook (§5.2.1) writes its Lambertian model as outgoing ray power $a\,\ell_{\text{in}}(\mathbf{n}\cdot\mathbf{p})$ (its Eq. 5.2), in its own units of ray power and with no $1/\pi$; in §5.3, its integral of that law over incoming directions prints $\cos(\mathbf{n}\cdot\mathbf{p})$ where $\mathbf{n}\cdot\mathbf{p}$, itself the cosine, is meant, and refers to the Lambertian law as Equation 5.1, which is the general BRDF equation (the Phong paragraph does the same). When radiance is the unit, the factor $1/\pi$ is needed. The lidar range equation uses the same $\Gamma/\pi$ for a Lambertian target, and multispectral satellite imaging the same $\pi$ when it converts radiance to reflectance. The wider family, from Phong to microfacet BRDFs, is the subject of reflectance models; photometric stereo surveys the ones used for shape recovery.

Conservation of radiance

In empty space radiance does not change along a ray (our derivation). Take two small areas $dA_1$ and $dA_2$ a distance $r$ apart, tilted by $\theta_1$ and $\theta_2$ from the line joining them. The flux leaving $dA_1$ towards $dA_2$ is $L_1\, dA_1 \cos\theta_1 \cdot dA_2\cos\theta_2/r^2$, and the flux arriving at $dA_2$ from $dA_1$ is $L_2\, dA_2 \cos\theta_2 \cdot dA_1\cos\theta_1/r^2$. With no loss in between these are equal, so $L_2 = L_1$. The factor $dA_1 dA_2 \cos\theta_1\cos\theta_2/r^2$ is symmetric in the two areas; it is the geometric “throughput” of the pencil of rays.

A camera exploits this. A pixel’s footprint on a surface grows as $d^2$ with distance $d$, while the solid angle of the lens seen from the surface falls as $1/d^2$, so the flux per pixel is fixed by the surface’s radiance and by the camera alone. Image formation gives the result, the irradiance on the sensor $E = L\,(\pi/4)(d_a/f)^2\cos^4\alpha$ for an aperture of diameter $d_a$; aperture and f-number explains the $1/N^2$ in it, and exposure integrates the irradiance over time. Methods that read pixel values as measurements, such as shape from shading, need images whose values are proportional to scene radiance, which is why they ask for raw or linearized data.

Spectral radiance and Planck’s law

A blackbody at temperature $T$ emits the spectral radiance given by Planck’s law (standard physics, stated here as our statement):

$$B_\lambda(T) = \frac{2hc^2}{\lambda^5}\,\frac{1}{\exp\!\big(hc/(\lambda k T)\big) - 1},$$

with Planck’s constant $h$, the speed of light $c$ and Boltzmann’s constant $k$. Thermal and infrared cameras apply the same law in the infrared. Integrated over wavelength and multiplied by $\pi$, it gives the Stefan–Boltzmann exitance $\sigma T^4$; the example checks this numerically with the CODATA constants in scipy.constants. Conventions differ between libraries: colour-science’s planck_law (v0.4.7) takes the wavelength in metres and returns radiance per metre of wavelength (its docstring gives the unit as W sr⁻¹ m⁻², without the wavelength; its sd_blackbody multiplies by $10^{-9}$ to get per nanometre), and uses the exitance constant $c_1 = 3.741771 \times 10^{-16}$ W m² (printed W/m² in its docstring) divided by $\pi$, with $c_2 = 1.4388 \times 10^{-2}$ m K.

Photometry

The luminous efficiency function and the candela

Photometric quantities are radiometric quantities weighted by the photopic luminous efficiency function $V(\lambda)$, the standard observer’s relative sensitivity in daylight, and scaled by a constant $K_m$:

$$\Phi_v = K_m \int V(\lambda)\, \Phi_\lambda(\lambda)\, d\lambda .$$

Both colour-science (luminous_flux, by the trapezoid rule, with the “CIE 1924 Photopic Standard Observer” as its default $V$) and ISETCam (ieLuminanceFromEnergy.m, as $683 \sum E\, V \Delta\lambda$, with $\Delta\lambda$ the sample spacing, or 10 nm by default for a single wavelength) compute exactly this. The same weighting gives the tristimulus value $Y$ of the CIE 1931 observer: color spaces notes that $Y$ is defined to be proportional to luminance, and ISETCam’s ieXYZFromEnergy.m returns $Y$ in cd/m² for spectral radiance and in lux for spectral irradiance. Wyman, Sloan and Shirley (2013) give analytic fits to the 1931 colour-matching functions; their multi-lobe fit to $\bar{y}$ is a sum of two piecewise Gaussians (Eq. 4, Table 1), and the example uses it in place of the tabulated curve.

The constant comes from the SI. In the current SI the candela is defined by fixing the luminous efficacy of monochromatic radiation of frequency $540 \times 10^{12}$ Hz, $K_{cd}$, at exactly 683 lm/W (BIPM, 2019, §2.2, p. 123, and p. 131). The effect is that one candela is the luminous intensity, in a given direction, of a source emitting monochromatic radiation of that frequency with a radiant intensity of $1/683$ W/sr in that direction (p. 131). The frequency corresponds to a vacuum wavelength of $c/\nu = 555.17$ nm (our arithmetic), where the example’s fit to $\bar{y}$ reads 0.998. The brochure calls $K_{cd}$ a technical constant that links the physical power stimulating the eye to the photobiological response of a standard observer (p. 125). colour-science uses $K_m = 683$ lm/W and notes that the value has been rounded from 683.002 for practical use, citing Wyszecki and Stiles (not read here). The lumen is a candela steradian, the lux a lumen per square metre, and the candela per square metre the unit of luminance.

Efficacy and efficiency

colour-science separates two terms. Its luminous_efficiency is the dimensionless ratio $\int V\Phi_\lambda\,d\lambda / \int \Phi_\lambda\,d\lambda$; its luminous_efficacy is that ratio times 683, in lm/W. The efficacy of a source therefore depends on the range over which its radiant power is counted: a spectrum normalized only over the visible gives a far larger number than the same spectrum normalized over all wavelengths. The example normalizes over all wavelengths.

Properties

• The inverse-square law is a far-field law. It is exact for a point, and for an extended source it is wrong by $R^2/h^2$ on the axis of a disc. Off the axis a finite flat source is better approximated by a Lambertian point ($\cos^4$) than by an isotropic one ($\cos^3$).
• Radiance is invariant along rays; irradiance is not. Distance changes how much of a surface each pixel sees, not how bright it looks, until the surface is smaller than a pixel.
• Cosine factors accumulate. Incidence, emission and distance each contribute a cosine; on a floor under a point source the result is $\cos^3\theta$ or $\cos^4\theta$, and inside a camera the lens adds its own $\cos^4\alpha$.
• Everything is linear. Radiometric quantities from several sources add, which is why a linear image of a scene under two lamps is the sum of the images under each.
• Photometric quantities describe a human observer, not a camera. A sensor’s channels have their own spectral responses, which image formation notes do not match the CIE colour-matching functions; in our reading, an illuminance in lux tells a camera user how bright a scene is for people, not how many electrons a given sensor will collect.
• The factor π is projected solid angle. It turns radiance into exitance for a Lambertian surface, a uniform sky into irradiance, and a Lambertian BRDF $\rho/\pi$ into a reflectance $\rho$.

Examples

Inverse square, cosine and a finite source, seen by a camera

The script below has four parts, all exact computations on synthetic scenes with known radiometry. Units are metres, watts and steradians; the floor is the plane $z = 0$; every source and camera faces straight down; $\theta$ is the angle from the vertical at the floor point. Radiance images store the scene radiance $\rho E/\pi$ of a Lambertian floor with $\rho = 0.5$, the quantity in front of any lens.

• M1, a disc against point models. A uniform Lambertian disc of radius $R = 0.2$ m faces down at heights $h = R$ to $20R$. A point source with the same on-axis intensity $I_0 = \pi R^2 L$ gives $I_0/h^2$ on the axis; the disc’s exact on-axis value is computed by numerical quadrature. Off the axis, the disc is compared with its far-field model, a Lambertian point ($I_0\cos^4\theta/h^2$), over $0 \le x \le 2h$, and the script finds where the error changes sign.
• M2, a floor image. A 257 × 257 orthographic radiance image of the floor over $[-2h, 2h]^2$ at $h = 1$ m, under an isotropic point, a Lambertian point and the disc with $h/R = 5$.
• M3, radiance invariance. A square Lambertian patch, 0.256 m on a side and uniformly lit ($E_0 = 100$ W/m², so $L = 50/\pi$ W m⁻² sr⁻¹), is seen by a downward pinhole camera (10 mm pinhole distance, 10 µm pixels, 256 × 256, aperture radius 2.5 mm) from 2, 4 and 8 m. Each of 1024² floor samples sends $L\cos\theta_e\,dA \cdot A_{\text{ap}}\cos\theta_a/r^2$ into the aperture (a small-aperture approximation), and the power lands in the pixel to which the sample projects. Pixel values are sensor irradiances in units of $L A_{\text{ap}}/f^2$, so a pixel on the axis should read 1, and a pixel whose ray makes angle $\alpha$ with the axis should read $\cos^4\alpha$, the term image formation gives. Since $\cos^4\alpha \approx 1 - 2\tan^2\alpha$, the patch mean should fall below 1 by about $2\langle\tan^2\alpha\rangle = s^2/(3d^2)$ for a patch of side $s$ at distance $d$ (our derivation): 0.005461, 0.001365 and 0.000341 at 2, 4 and 8 m (our arithmetic).
• M4, photometry. One watt at 540 THz, by the SI definition and with Wyman et al.’s fit to $\bar{y}$; and one watt of a 2856 K blackbody spectrum, normalized over all wavelengths ($\sigma T^4/\pi$), converted to lumens with a 1 nm trapezoid and with a 10 nm bin sum (ISETCam’s convention), with the fraction of its power between 380 and 780 nm.

Seven checks verify the code: the hemisphere’s solid angle and projected solid angle by quadrature (C1); the closed-form disc irradiance against quadrature (C2) and Monte Carlo integration with $10^6$ points (C3, seed 0); the rendered images against $\cos^3$ and $\cos^4$ (C4); the pixel-count ratios as integers (C5); Planck’s law against Stefan–Boltzmann (C6); and the same Planck formula, with colour-science’s constants, against the value printed in colour-science’s docstring for 500 nm and 5500 K (C7; the formula is written from Planck’s law, not ported). The hypotheses, fixed before the pilot, were: (H1) the on-axis error of the inverse-square law equals $R^2/h^2$ to within $10^{-9}$; (H2a) for $h \ge 2R$ the largest error of the Lambertian-point model over $0 \le x \le 2h$ is on the axis; (H2b) for $h \ge 5R$ the error changes sign within 0.1° of the first-order angle $\tan\theta \approx 1/\sqrt{2} + (5\sqrt{2}/36)(R/h)^2$ given above; (H3) the floor images follow $\cos^3$ and $\cos^4$; (H4) in M3 the centre pixel is the same at all three heights to within $10^{-7}$, the pixel counts fall exactly by 4 and 16, the patch mean falls below 1 by the first-order $\cos^4$ deficit $2\langle\tan^2\alpha\rangle$ to within 10% of that deficit, and after division by $\cos^4\alpha$ every pixel reads 1 to within $10^{-6}$; (H5) the fit’s $\bar{y}$ at 555.1712 nm differs from 1 by no more than 0.0080, the square root of the largest squared error Wyman et al. report for their $\bar{y}$ fit (Table 2, 6.4e-5); (H6) the 2856 K source gives less than 68.3 lm per watt, the two integrations agree to 0.5%, and less than half its power lies between 380 and 780 nm. H5 assumes that the tabulated CIE $\bar{y}$ is 1 at that wavelength; the CIE table was not read for this article, so that premise is unverified. The parameters of M1 and M3 were chosen to show derived effects, so these parts illustrate their derivations rather than test them. The script runs in under a second with NumPy 2.3.5 and SciPy 1.16.3 (Python 3.13).

import os
os.environ["OPENBLAS_NUM_THREADS"] = os.environ["OMP_NUM_THREADS"] = "1"
import numpy as np
from scipy import integrate, optimize, constants as sc

# Units: metres, watts, steradians; wavelengths in nm where stated. Floor is z = 0, sources and
# cameras look straight down; theta is the angle from the vertical at the floor point.

def disc_closed(x, h, R, L=1.0):        # irradiance under a Lambertian disc (radiance L, radius R)
    s = h*h + x*x                        # at height h, floor point at horizontal offset x
    return 0.5*np.pi*L*(1 - (s - R*R)/np.sqrt((s + R*R)**2 - 4*x*x*R*R))

def disc_quad(x, h, R, L=1.0):           # E = L * integral of cos(e) cos(i) / r^2 over the disc
    f = lambda rho, phi: rho/(rho*rho + x*x - 2*rho*x*np.cos(phi) + h*h)**2
    return L*h*h*integrate.dblquad(f, 0, 2*np.pi, 0, R, epsabs=0, epsrel=1e-12)[0]

def disc_mc(x, h, R, n, rng, L=1.0):     # Monte Carlo over uniform points on the disc
    rho, phi = R*np.sqrt(rng.random(n)), 2*np.pi*rng.random(n)
    g = h*h/(rho*rho + x*x - 2*rho*x*np.cos(phi) + h*h)**2
    return L*np.pi*R*R*g.mean(), L*np.pi*R*R*g.std(ddof=1)/np.sqrt(n)

print("== Checks")
w = integrate.dblquad(lambda th, ph: np.sin(th), 0, 2*np.pi, 0, np.pi/2)[0]
W = integrate.dblquad(lambda th, ph: np.cos(th)*np.sin(th), 0, 2*np.pi, 0, np.pi/2)[0]
print("C1 hemisphere: solid angle %.12f (2 pi), projected %.12f (pi)" % (w, W))
R = 0.2
c2 = max(abs(disc_quad(xh*hr*R, hr*R, R)/disc_closed(xh*hr*R, hr*R, R) - 1)
         for hr in [1, 2, 3, 5, 10, 20] for xh in [0, 0.5, 1, 2])
print("C2 disc, quadrature vs closed form: max rel. diff %.1e (tol 1e-9)" % c2)
rng = np.random.default_rng(0)
zs = []
for hr, xh in [(1, 0), (1, 1), (5, 0), (5, 1)]:
    m, se = disc_mc(xh*hr*R, hr*R, R, 10**6, rng)
    zs.append((m - disc_closed(xh*hr*R, hr*R, R))/se)
print("C3 disc, Monte Carlo (1e6 points, seed 0) z-scores:", " ".join("%.2f" % z for z in zs), "(tol |z| < 4)")

print("\n== M1: a disc of radius R against point models with the same on-axis intensity I0 = pi R^2 L")
print(" h/R  on-axis err   R^2/h^2     diff      max|far/disc-1| at x/h   sign change (deg)")
for hr in [1, 2, 3, 5, 10, 20]:
    h = hr*R
    E_disc = disc_quad(0.0, h, R)                      # exact on-axis value by quadrature
    err0 = (np.pi*R*R/h**2)/E_disc - 1                 # I0/h^2 against the disc
    xs = np.linspace(0, 2*h, 2001)
    far = (np.pi*R*R)*np.cos(np.arctan(xs/h))**4/h**2  # Lambertian point: I0 cos(theta) * cos(theta)/r^2
    rel = far/disc_closed(xs, h, R) - 1
    k = np.argmax(np.abs(rel))
    g = lambda x: (np.pi*R*R)*np.cos(np.arctan(x/h))**4/h**2/disc_closed(x, h, R) - 1
    ang = np.degrees(np.arctan(optimize.brentq(g, 0.0, 2*h, xtol=1e-14)/h)) if g(0)*g(2*h) < 0 else np.nan
    print("%4d  %.8f  %.8f  %8.1e   %.6f   %.4f        %.3f" % (hr, err0, 1/hr**2, err0 - 1/hr**2, abs(rel[k]), xs[k]/h, ang))

print("\n== M2: floor radiance image (257 x 257, orthographic, x and y in [-2h, 2h]), h = 1 m, rho = 0.5")
h, rho, n = 1.0, 0.5, 257
xs = np.linspace(-2*h, 2*h, n); X, Y = np.meshgrid(xs, xs)
r2 = X**2 + Y**2 + h*h; cos = h/np.sqrt(r2)
img_iso = rho*(1.0*cos/r2)/np.pi                   # isotropic point, I = 1 W/sr: E = I cos/r^2
img_lam = rho*(1.0*cos*cos/r2)/np.pi               # Lambertian point, I(theta) = 1 W/sr * cos
img_disc = rho*disc_closed(np.sqrt(X**2 + Y**2), h, h/5, L=1/(np.pi*(h/5)**2))/np.pi  # h/R = 5
print("C4 render vs cos^3/h^2 and cos^4/h^2: max rel. diff %.1e, %.1e (tol 1e-12)" % (
      np.abs(img_iso/(rho*cos**3/h**2/np.pi) - 1).max(), np.abs(img_lam/(rho*cos**4/h**2/np.pi) - 1).max()))
c = n//2
print(" x/h   iso/iso(0)  cos^3     lam/lam(0)  cos^4     disc/disc(0)")
for j in [c, c + 32, c + 64, c + 96, c + 128]:
    print("%4.1f   %.6f    %.6f  %.6f    %.6f  %.6f" % (xs[j], img_iso[c, j]/img_iso[c, c], cos[c, j]**3,
          img_lam[c, j]/img_lam[c, c], cos[c, j]**4, img_disc[c, j]/img_disc[c, c]))

print("\n== M3: a uniformly lit Lambertian patch seen by a pinhole camera from three heights")
f, p, a = 0.010, 10e-6, 2.5e-3                     # pinhole distance, pixel pitch, aperture radius
A_ap = np.pi*a*a; E0 = 100.0; L = rho*E0/np.pi     # patch irradiance and radiance (M = pi L)
side, ds = 0.256, 0.25e-3                          # patch side and floor sampling step
g = (np.arange(int(round(side/ds))) + 0.5)*ds - side/2
GX, GY = np.meshgrid(g, g)
ii = (np.arange(256) + 0.5 - 128)*p; U, V = np.meshgrid(ii, ii)
c4 = (f*f/(U*U + V*V + f*f))**2                    # cos^4 of each pixel centre's ray
rows = []
for d in [2.0, 4.0, 8.0]:
    r2 = GX**2 + GY**2 + d*d
    dP = L*(d*d/r2)*ds*ds*A_ap/r2                  # L cos(e) dA * A_ap cos(a) / r^2 into the aperture
    iu = np.floor(f*GX/(d*p) + 128).astype(int); iv = np.floor(f*GY/(d*p) + 128).astype(int)
    img = np.zeros((256, 256)); np.add.at(img, (iv, iu), dP)
    img = img/(p*p)/(L*A_ap/f**2)                  # sensor irradiance in units of L A_ap / f^2
    on = img > 0
    rows.append((d, on.sum(), dP.sum(), img[127:129, 127:129].mean(), img[on].mean(),
                 (img[on]/c4[on]).mean(), np.abs(img[on]/c4[on] - 1).max()))
print("  d   pixels  P_patch (W)    P(2)/P(d)  (d/2)^2  centre       patch mean   mean/cos^4   max|v/cos^4-1|")
for d, npx, P, cen, mean, mc4, dev in rows:
    print("%3.0f  %6d   %.6e   %8.5f  %6.1f   %.9f  %.6f     %.9f  %.1e" % (
          d, npx, P, rows[0][2]/P, (d/2)**2, cen, mean, mc4, dev))
print("C5 pixel-count ratios %d, %d (exact 4, 16)" % (rows[0][1]//rows[1][1], rows[0][1]//rows[2][1]),
      "remainders %d, %d" % (rows[0][1] % rows[1][1], rows[0][1] % rows[2][1]))

print("\n== M4: photometry with Wyman, Sloan and Shirley's multi-lobe fit to ybar (1931), K = 683 lm/W")
def ybar(lam):                                     # lam in nm; Wyman et al. (2013), Eq. 4, Table 1
    t1 = (lam - 568.8)*np.where(lam < 568.8, 0.0213, 0.0247)
    t2 = (lam - 530.9)*np.where(lam < 530.9, 0.0613, 0.0322)
    return 0.821*np.exp(-0.5*t1*t1) + 0.286*np.exp(-0.5*t2*t2)
K, d_obs = 683.0, 2.0
lam540 = sc.c/540e12*1e9
for name, eff in [("SI definition (540 THz)", K), ("683 * ybar_fit(%.4f nm)" % lam540, K*ybar(lam540))]:
    print("1 W isotropic, %-30s %8.3f lm, %7.3f cd, %7.4f lx at %.0f m" % (
          name + ":", eff, eff/(4*np.pi), eff/(4*np.pi)/d_obs**2, d_obs))
print("ybar_fit(555.1712) - 1 = %+.5f (bound from Table 2: 0.0080)" % (ybar(lam540) - 1))

T = 2856.0
def planck(lam_m, T, c1L=2*sc.h*sc.c**2, c2=sc.h*sc.c/sc.k):   # spectral radiance, W m^-2 sr^-1 m^-1
    return c1L/lam_m**5/np.expm1(c2/(lam_m*T))
lo = integrate.quad(planck, 1e-8, 1e-6, args=(T,), epsabs=0, epsrel=1e-12, limit=200)[0]
hi = integrate.quad(planck, 1e-6, 1e-2, args=(T,), epsabs=0, epsrel=1e-12, limit=500,
                    points=[1e-5, 1e-4, 1e-3])[0]
print("C6 pi * integral of B over wavelength / (sigma T^4) - 1 = %.1e (tol 1e-9)" % (np.pi*(lo + hi)/(sc.sigma*T**4) - 1))
v = 3.741771e-16/np.pi/(500e-9)**5/np.expm1(1.4388e-2/(500e-9*5500))
print("C7 colour-science constants at 500 nm, 5500 K: %.4f (docstring 20472701909806.5...)" % v)
B_tot = sc.sigma*T**4/np.pi                        # integral of B over all wavelengths
lam1 = np.arange(360.0, 831.0, 1.0); lam10 = np.arange(360.0, 831.0, 10.0)
S1 = planck(lam1*1e-9, T)*1e-9/B_tot               # spectrum per nm, normalized to 1 W in total
S10 = planck(lam10*1e-9, T)*1e-9/B_tot
eff_trap = K*np.trapezoid(ybar(lam1)*S1, lam1)
eff_bin = K*np.sum(ybar(lam10)*S10)*10.0
lam_v = np.arange(380.0, 781.0, 1.0)
frac_vis = np.trapezoid(planck(lam_v*1e-9, T)*1e-9, lam_v)/B_tot
print("2856 K, 1 W: %.4f lm (1 nm trapezoid), %.4f lm (10 nm bin sum), rel. diff %.2e" % (
      eff_trap, eff_bin, eff_bin/eff_trap - 1))
print("2856 K: fraction of radiant power in 380-780 nm %.4f; %.4f lx at %.0f m (isotropic)" % (
      frac_vis, eff_trap/(4*np.pi)/d_obs**2, d_obs))

== Checks
C1 hemisphere: solid angle 6.283185307180 (2 pi), projected 3.141592653590 (pi)
C2 disc, quadrature vs closed form: max rel. diff 1.2e-12 (tol 1e-9)
C3 disc, Monte Carlo (1e6 points, seed 0) z-scores: -0.33 -0.33 1.11 1.30 (tol |z| < 4)

== M1: a disc of radius R against point models with the same on-axis intensity I0 = pi R^2 L
 h/R  on-axis err   R^2/h^2     diff      max|far/disc-1| at x/h   sign change (deg)
   1  1.00000000  1.00000000   4.4e-16   1.000000   0.0000        41.132
   2  0.25000000  0.25000000   2.2e-16   0.250000   0.0000        37.014
   3  0.11111111  0.11111111   5.6e-17   0.111111   0.0000        36.072
   5  0.04000000  0.04000000   3.5e-17   0.040000   0.0000        35.561
  10  0.01000000  0.01000000   8.7e-18   0.010000   0.0000        35.339
  20  0.00250000  0.00250000  -2.8e-16   0.002500   0.0000        35.283

== M2: floor radiance image (257 x 257, orthographic, x and y in [-2h, 2h]), h = 1 m, rho = 0.5
C4 render vs cos^3/h^2 and cos^4/h^2: max rel. diff 5.6e-16, 5.6e-16 (tol 1e-12)
 x/h   iso/iso(0)  cos^3     lam/lam(0)  cos^4     disc/disc(0)
 0.0   1.000000    1.000000  1.000000    1.000000  1.000000
 0.5   0.715542    0.715542  0.640000    0.640000  0.656954
 1.0   0.353553    0.353553  0.250000    0.250000  0.262547
 1.5   0.170677    0.170677  0.094675    0.094675  0.099770
 2.0   0.089443    0.089443  0.040000    0.040000  0.042069

== M3: a uniformly lit Lambertian patch seen by a pinhole camera from three heights
  d   pixels  P_patch (W)    P(2)/P(d)  (d/2)^2  centre       patch mean   mean/cos^4   max|v/cos^4-1|
  2   16384   5.092197e-06    1.00000     1.0   0.999998672  0.994570     0.999999675  3.3e-07
  4    4096   1.278255e-06    3.98371     4.0   0.999998668  0.998637     0.999999669  3.3e-07
  8    1024   3.198908e-07   15.91855    16.0   0.999998667  0.999659     0.999999667  3.3e-07
C5 pixel-count ratios 4, 16 (exact 4, 16) remainders 0, 0

== M4: photometry with Wyman, Sloan and Shirley's multi-lobe fit to ybar (1931), K = 683 lm/W
1 W isotropic, SI definition (540 THz):        683.000 lm,  54.351 cd, 13.5879 lx at 2 m
1 W isotropic, 683 * ybar_fit(555.1712 nm):    681.538 lm,  54.235 cd, 13.5588 lx at 2 m
ybar_fit(555.1712) - 1 = -0.00214 (bound from Table 2: 0.0080)
C6 pi * integral of B over wavelength / (sigma T^4) - 1 = -6.6e-12 (tol 1e-9)
C7 colour-science constants at 500 nm, 5500 K: 20472701909806.5820 (docstring 20472701909806.5...)
2856 K, 1 W: 16.4737 lm (1 nm trapezoid), 16.4744 lm (10 nm bin sum), rel. diff 4.32e-05
2856 K: fraction of radiant power in 380-780 nm 0.1053; 0.3277 lx at 2 m (isotropic)

Every check held. The solid angle of the hemisphere is $2\pi$ and its projected solid angle $\pi$ (C1), the closed-form disc irradiance agrees with quadrature to $1.2 \times 10^{-12}$ (C2), the four Monte Carlo estimates lie within 1.3 standard errors of it (C3), the rendered images match the cosine formulas to rounding (C4), the pixel counts fall by exactly 4 and 16 (C5), Planck’s law integrates to $\sigma T^4$ to $6.6 \times 10^{-12}$ (C6), and the colour-science constants reproduce its docstring value (C7).

The inverse-square law against a disc. H1 held at all six heights: the error of $I_0/h^2$ is $R^2/h^2$ to rounding, a factor of two at $h = R$, 4% at $h = 5R$ and 1% at $h = 10R$, where the disc’s diameter is a fifth of the distance. H2a held: for every height, including $h = R$, the Lambertian-point model’s largest error over the floor is on the axis, where it overestimates by $R^2/h^2$. Off the axis the error shrinks, changes sign, and the model underestimates; the first-order coefficient of the series above, $(2\tan^2\theta - 1)\cos^4\theta$, reaches 1 in magnitude on the axis and at most 1/3 elsewhere (our derivation). H2b held: for $h/R = 5$, 10 and 20 the printed sign-change angles lie within 0.003° of the first-order formula (35.563°, 35.339° and 35.283°, our arithmetic). At $h/R = 2$ and 3, where no criterion was set, they lie below it, by less than 0.1° (the formula gives 37.097° and 36.089°).

Cosine fall-off. H3 held. Along the floor the isotropic source falls to 0.354 of its centre value at $x = h$ (45°), and the Lambertian point to 0.25. The disc at $h/R = 5$ follows the Lambertian point, not the isotropic one, lying slightly above its normalized profile (by 5.2% at $x = 2h$), because its centre value is the one the far-field model overestimates. In our reading, a flat diffusing source therefore lights a floor less evenly than the isotropic $\cos^3$ suggests.

Radiance invariance. H4 held in all its parts. Moving the camera from 2 m to 8 m divides the patch’s pixel count by exactly 16 and the power the aperture collects by 15.92, yet the centre pixel reads the same value, 0.9999987, at every height, to within $5 \times 10^{-9}$. That is the conservation of radiance: each pixel sees a footprint sixteen times larger through a solid angle sixteen times smaller. The patch mean falls below 1 by close to the first-order deficits $s^2/(3d^2)$, well within the 10% the criterion allows, because the outer pixels see the patch at angles where the camera’s own $\cos^4\alpha$ applies; dividing by $\cos^4\alpha$ removes it to $3.3 \times 10^{-7}$. The same term explains why the collected power falls slightly less than $d^2$ (by 15.92 rather than 16): at 2 m more of the patch lies off the axis. The centre reads slightly below 1 for the same reason, since the four central pixels’ rays are not quite on the axis.

Lumens and lux. H5 held: the fit gives $\bar{y}(555.1712\text{ nm}) = 0.99786$, so one watt at 540 THz yields 681.5 lm instead of the SI’s 683 lm by definition, a deviation of 0.21%, well inside the fit’s reported error; a 1 W isotropic source of that frequency gives 54.35 cd and 13.59 lx at 2 m. H6 held: one watt of a 2856 K blackbody yields only 16.47 lm, 2.4% of 683, because only 10.5% of its power lies between 380 and 780 nm and, in our reading, much of that lies at the red end, where $V$ is small; the 10 nm bin sum agrees with the 1 nm trapezoid to $4.3 \times 10^{-5}$, since both the spectrum and the fit are smooth on that scale.

Process note. The plan for this example listed eight departures from the batch plan’s description, before the pilot: (1) the floor image is 257 × 257 instead of 256 × 256, so that a pixel lies on the axis and the edges at $\pm 2h$, and orthographic, with a pinhole camera only in M3; (2) the far-field match is taken as $I_0 = \pi R^2 L$, and M1 adds the off-axis comparison with the Lambertian-point model and its sign-change angle, since an isotropic point is not the disc’s far field off the axis; (3) the disc’s exact irradiance uses a closed form (our derivation), checked by quadrature and Monte Carlo, where the batch plan named only those two; (4) M2 prints the $\cos^4$ and disc profiles beside $\cos^3$; (5) M3 lights the patch uniformly rather than with the point source, computes pixel values by radiometric transfer and adds the total power, the centre value and the $\cos^4$-divided mean, so that radiance invariance is not mixed with the lens term; (6) M4 normalizes the 2856 K spectrum over all wavelengths and adds the 10 nm bin sum and the visible fraction, replacing a planning expectation that compared the result with colour-science’s tabulated value, which cannot be computed without the CIE table; (7) checks C1–C7 were added; (8) the planning research’s wording of the inverse-square hypothesis (“below 1% for h ≥ 10R”) was corrected, since the error is exactly 1% at $h = 10R$. Three planning scripts ran before the pilot. Two, symbolic, computed the series of the disc against the far-field model; in the plan’s words, “H2 was written after this script printed, and is disclosed as such”. The first version of the second one printed a coefficient as nan and stopped with an error; it was fixed by taking the series coefficients from the expanded expression, and only the fixed version’s output was used. A third script computed the predicted $\cos^4$ deficits and the H5 bound from formulas. A reachable-value script then set the check tolerances on the pilot’s set-up; in the plan’s words, it “also printed the M3 centre 2×2 values for d = 2, 4, 8 m … which is H4a’s measure. H4a is therefore written after its value was seen and is an illustration, not a test”. The post-pilot log records the editor’s instruction to disclose both, quoted: “H2 written after the SymPy series (plan01, plan02, including plan02’s failed first version) and H4a printed by reach.py (illustration)”. The pilot’s output was identical to the first run’s; no setting, measure or printed value changed after the pilot.

Common Misconceptions

• “Intensity falls off with the square of distance.” Radiant intensity, power per steradian, is a property of the source and does not change with distance. The irradiance it produces on a surface falls as $1/r^2$, and only for a point or at distances large compared with the source.
• “Distant objects look darker.” A surface’s radiance, and so its pixel value, does not depend on its distance in clear air; the example’s centre pixel is the same at 2, 4 and 8 m. Distant objects look smaller, and become dimmer only once they are smaller than a pixel.
• “A Lambertian surface emits equally in all directions.” It has the same radiance in all directions, so its radiant intensity falls as $\cos\theta$, and its exitance is $\pi L$, not $2\pi L$.
• “A Lambertian BRDF is the albedo.” It is $\rho/\pi$ sr⁻¹ (Nicodemus et al., 1977, Eq. C4). Models written in units of ray power, such as visionbook’s Eq. 5.2, leave out the $\pi$, which matters as soon as radiance is compared with irradiance.
• “Luminance is the Y channel of YCbCr.” That channel is luma, a weighted sum of gamma-encoded values; luminance is a photometric quantity in cd/m², linear in the light, as color spaces and gamma correction explain.
• “Lumens measure a lamp’s power.” They measure power weighted by human vision. In the example one watt yields 683 lm at 540 THz and 16.5 lm from a 2856 K blackbody.
• “Lux describes what a camera sees.” It describes what a standard observer would see. A camera with a different spectral response, or with sensitivity in the near infrared, collects a different signal from two scenes of equal illuminance.

Where It Is Used

Image formation is the chain from scene radiance to pixel value, and its reflection equation, lens irradiance formula and sensor model are written in the quantities defined here. Exposure integrates sensor irradiance over time into radiant exposure, aperture and f-number explains why the irradiance falls as $1/N^2$, and the image sensor article converts that energy into photons and electrons. Shape from shading and photometric stereo invert the Lambertian law to recover surface orientation, and brightness constancy assumes that a scene point keeps its value between frames, which holds for Lambertian surfaces under unchanged light because their radiance does not depend on the viewing direction (our reading). HDR imaging recovers radiance maps from bracketed exposures, and neural radiance fields and volume rendering model radiance along camera rays. Multispectral satellite imaging converts measured spectral radiance into reflectance, and the lidar range equation combines the inverse-square law with a Lambertian target. Photometric units appear where people set the scene: low-light image enhancement describes its datasets by the illuminance at the camera in lux, gamma correction specifies displays in cd/m², and color spaces builds the CIE XYZ space on the same luminous weighting through $Y$.