Thirty unknowns instead of six

A new piece of machinery in a collection like this one is a liability until it has reproduced something already known. This one has an unusually clean way of doing that, and the way is the reason to trust everything downstream of it.

The claim #

The directional solver reduces to the radiosity solver exactly, and the reduction is what makes every other number it produces a statement about lobes rather than about new arithmetic.

• Thirty unknowns rather than six: one radiance per ordered pair of patches in a six-face box, per wavelength band.
• The transport is L(i→j) = Lᵉ(i→j) + Σₖ f_i(ω_ik → ω_ij) κ_ik L(k→i), with κ_ik = π F_ik the cosine-weighted solid angle patch k subtends at i.
• Setting f = ρ/π collapses it. Every one of the thirty radiances becomes its patch’s radiosity divided by π, and the answer matches the existing solver at the floating-point floor.
• And that is a stronger check than agreement to a tolerance. A method that agreed to three decimal places would be a different method that happened to be close.

What the unknown is #

Radiosity’s unknown is one number per patch: the total power leaving per unit area, identical in every direction. The directional solver’s unknown is one radiance per ordered pair: how much light leaves patch i in the direction of patch j, which for six faces is thirty numbers rather than six.

Each patch is treated as a point with a full bidirectional distribution on it. Light arrives at i from k along the direction ω_ik, is scattered by f_i, and leaves towards j along ω_ij. How much arrives is the radiance from k times the cosine-weighted solid angle k subtends at i, and that quantity is already computed: it is π times the form factor.

That identification is the piece of luck the whole construction rests on. F_ik is defined as the fraction of light leaving i diffusely that reaches k, and by reciprocity and the definition of a form factor, π F_ik is exactly ∫ cos θ dω over the solid angle k occupies at i. So a solver that already has form factors has, for free, the geometry a directional solver needs.

The reduction is the check that matters. There are two ways to gain confidence in a new numerical method. One is to compare it against a reference on a problem both can solve. The other is to show that it becomes the reference under a parameter setting.

The second is much stronger and it is available here. Putting f = ρ/π — a Lambertian surface, no lobe at all — makes the scattering independent of both directions, so

L(i→j) = Eᵢ/π + (ρᵢ/π) Σₖ κ_ik L(k→i)

has a right-hand side that does not depend on j. All thirty unknowns for a given i are therefore equal, and substituting L = B/π recovers the radiosity equation exactly.

That is an algebraic identity, so the numerical check is asking whether the implementation has it. It does, to 9.8 × 10⁻¹⁶ relative — which is what an identity looks like in double precision after a few thousand operations, and is not what a merely-similar method looks like.

A check at the floating-point floor rules out a class of errors that a check at a tolerance cannot. A sign error in a cosine, a transposed index, a form factor used in the wrong direction: all of those would produce agreement to a few per cent on a symmetric room and none of them survives a comparison at 10⁻¹⁶.

The reduction, approached rather than imposed #

The check above imposes the Lambertian limit exactly. The roughness sweep approaches it, and the two together are more convincing than either.

At roughness 0.8 the departure is 0.969 ΔE₀₀, at 0.6 it is 1.457, at 0.45 it is 2.051, and the sequence rises monotonically towards the glossy end. Extrapolating the other way, the sequence heads towards zero, and the exact check says it arrives.

That monotonicity is asserted rather than observed, and it could have failed. A solver with a subtle error in its lobe would very likely produce a non-monotone sequence — a bump somewhere in the middle where two errors of different sign crossed — and the assertion would report it.

Two independent checks on the same machinery, one exact and one about shape, and neither is implied by the other. The exact one says the limit is right; the monotone one says the path to it is.

That pair of curves is a second kind of check and it is worth having beside the reduction. The reduction verifies that the solver’s geometry is right. This verifies that its results are sensitive to the modelling decisions they should be sensitive to — a solver whose answer did not change when the energy accounting changed would be a solver in which the lobe was not doing anything.

The two together bracket a common failure. Machinery that reproduces a known limit and responds to nothing else is machinery whose new degrees of freedom are inert; machinery that responds to everything and reproduces no limit is machinery nobody should trust. Passing both is not proof of correctness and it eliminates the two cheapest ways of being wrong.

What a point-sampled patch commits to #

Every model has a commitment and this one’s is worth stating before any of its results are used.

A patch here is a point with a bidirectional distribution on it. The direction from patch i to patch j is the direction between their centres, and the variation of that direction across i’s own extent is ignored.

That is the same thing form factors already do — a form factor is an average over both patches — so the commitment is not new to this solver. What is new is that it now matters differently, because a Lambertian surface does not care about direction and a lobe cares about nothing else.

The consequence is a resolution limit rather than a bias. A lobe wider than the cone a patch subtends is averaged correctly; one narrower is not, and the boundary is measured at about roughness 0.15 for a cube. Below it the model returns unphysical answers rather than merely imprecise ones, which is the honest failure mode and is the one that gets noticed.

Why the lobe is averaged over a cone #

The first version of this solver evaluated f_i(ω_ik → ω_ij) at the two centre directions and it was catastrophic.

Two patch centres in a cube sit at directions where a microfacet distribution can be enormous — near a mirror direction, the distribution’s value goes as the inverse square of the roughness — and the transport then took that spike as though the whole patch were that bright. At roughness 0.05 the solver returned a chroma of thirty-four thousand and a negative lightness.

The repair is to average the lobe over the cone each patch actually subtends. Ω is read off the form factor, the cone is sampled at its centre plus three rings of twelve, and both the incoming and outgoing directions are averaged.

That is not a numerical detail; it is the difference between a model and a mistake. A patch subtends a cone and evaluating at its centre is a quadrature of one point, which is adequate for a function that varies slowly across the cone and useless for one that varies by orders of magnitude.

The lesson generalises past this solver. Any calculation that replaces an extended object by a point is performing a one-point quadrature, and whether that is acceptable depends entirely on what is being integrated — which is a property of the integrand rather than of the geometry, and is therefore not visible from the geometry alone.

What it costs to run #

The system is thirty by thirty per wavelength band and there are eighty-one bands, so a room is eighty-one Gaussian eliminations of a thirty-square matrix — a few million floating-point operations, which is nothing.

What costs is the geometry. Each of the thirty ordered pairs needs the lobe evaluated between every incoming patch and its own outgoing direction, averaged over two cones of thirty-seven directions each, which is about thirteen hundred bidirectional evaluations per pair per incoming patch. That is the bulk of the work and it is independent of wavelength, so it is computed once and reused across all eighty-one bands.

The separation of the spectral part from the geometric part is what makes the solver practical, and it is available because a Fresnel lobe is spectrally neutral — the same fact that makes a highlight the lamp rather than the paint. The body term carries the pigment’s spectrum and the lobe does not, so the geometry can be computed once in a spectrum-free way and combined with the reflectance band by band.

That is a modelling assumption as well as an optimisation: a dielectric interface’s reflectance varies weakly with wavelength through the refractive index, and treating it as constant is a simplification the solver depends on structurally.

The spread in that figure is a useful way to price the new machinery. Twenty-four additional unknowns per band, several thousand extra bidirectional evaluations, and a modelling commitment about point-sampled patches, in exchange for a quantity of about one colour difference that the previous method could not represent at all.

Whether that is worth it depends on the question. For a room’s overall cast it is a correction of a few per cent and radiosity is fine. For a question about what a surface looks like from a particular place — which is what a photograph is, and what an inspection is — it is the entire answer, and radiosity has nothing to say.

What was computed, and how #

The form factors are this collection’s own, from the closed forms for parallel and perpendicular rectangles with closure used to fix the four adjacent faces. They were verified against a numerical quadrature when they were written and the closure check is still in the suite, alongside the agreement between two independent form-factor computations.

The bidirectional distribution is the previous round’s microfacet surface with a Lambertian body, used unchanged. Reusing it rather than writing a second one is deliberate: two implementations of one object is how two implementations come to disagree, which this fleet has recorded about a different pair of measurements — most sharply about the two implementations of one reuse measure.

The reduction assertion compares every one of the thirty ordered-pair radiances against its patch’s radiosity divided by π, in every one of the eighty-one bands, and reports the worst relative difference. Comparing only the total, or only one patch, would pass on several errors that this catches.

What the extra unknowns do not buy #

It is worth being clear about what the new machinery does not add, because the list is longer than the list of what it does.

It does not add a bounce. The linear system is solved exactly, so all orders of interreflection are included, exactly as radiosity includes them. Nothing here is a truncation.

It does not add spectral resolution. The same eighty-one bands, the same reflectances, the same illuminant. A spectral solve was already what this collection does, against the three-channel shortcut a renderer would use.

And it does not add geometry. The same six faces, the same closed-form factors, the same cube. A room with furniture in it is as far out of reach as it was.

What it adds is exactly one thing: the direction light leaves a surface in. That single addition is worth between one and five colour differences and a quantity radiosity could not express, which is a good return for one axis — and it is worth noticing that the axis was available for two rounds and was not taken, because taking it meant rewriting a solver rather than parameterising one.

Where the model stops #

The patch-as-point commitment is the largest limitation and it is stated above.

The energy accounting takes the interface’s return out of the body term using the lobe’s directional-hemispherical albedo, which is the physically sensible choice and is not the only one anybody makes. The two choices disagree about the sign of the chroma result, so this is a modelling decision with consequences rather than a bookkeeping one.

And the lobe is isotropic, dielectric and single-scattering. A metallic surface has a coloured Fresnel term; a rough surface has multiple scattering between microfacets that the single-scattering model loses; and a real paint has both a specular lobe and a diffuse one with its own directional structure.

One practical note for anybody wanting to reproduce this. The reduction check is only available because the two solvers share their geometry — the same form factors, the same faces, the same closure. A directional solver built from scratch with its own geometry would produce an answer close to radiosity’s and would not reduce to it exactly, and the difference between “close” and “exact” is the whole value of the check.

So the way to build such a thing is to build it on top of the existing solver’s geometry rather than beside it, even when that constrains the formulation. Here it constrained it substantially: one direction per patch pair is a cruder discretisation than a hemispherical one, and it was chosen because it is the discretisation the form factors already carry.

The generalisation #

The habit is about how to introduce a method that supersedes one already in use.

The temptation is to demonstrate the new method on the case the old one could not handle, because that is the point of it. That demonstration establishes nothing, because there is nothing to compare against — the new answer is the only answer, and its correctness is exactly what is in question.

The move is to find the parameter setting under which the new method becomes the old one, and to check that it does, at the tightest tolerance the arithmetic supports. A method that reduces correctly has had its geometry, its indexing, its normalisation and its solver checked all at once, against an implementation somebody else already trusted.

The failure mode is to check against a tolerance that a wrong method could also meet. A check at three decimal places on a symmetric problem is nearly free, and a check at the floating-point floor is not — which is why the second one is worth arranging even when it takes a special case to arrange it.

A last observation about cost, since a fivefold increase in unknowns sounds expensive and is not. The linear algebra is negligible; the geometry dominates; and the geometry is spectrum-free, so it is amortised across every band. A room that took four seconds to solve by radiosity takes about six by directional transport in this implementation, and nearly all of the difference is the cone quadrature.

That ratio would not survive a finer discretisation. A hemispherical directional radiosity with a hundred directions per patch has ten thousand unknowns rather than thirty, and the linear algebra stops being negligible immediately — which is why the graphics literature abandoned the approach rather than refining it, and why a coarse version of it is still useful for a question about colour rather than about images.

Who found it, and when #

Directional extensions of radiosity date from Immel, Cohen and Greenberg in 1986, who discretised the hemisphere at each patch and solved for a directional radiance distribution. The formulation here is cruder — one direction per patch pair rather than a hemispherical discretisation — and is adequate because a box has six faces and a colour is an integral.

The identification of π F with the cosine-weighted solid angle is standard and is usually stated as the reason form factors sum to one in a closed enclosure. Using it to build a directional solver on top of an existing radiosity implementation appears not to be written down anywhere, probably because in graphics the natural move is to abandon patches entirely and trace paths.

Where the ladder goes next #

The machinery is checked and the first result it produces is one the Lambertian assumption could not have suggested: a lobe takes colour out of a bounce, and it does so for a reason that has nothing to do with how much light it carries.