Editing Radiosity (computer graphics) (section)

== Mathematical formulation ==
The basic radiosity method has its basis in the theory of [[heat|thermal radiation]], since radiosity relies on computing the amount of light energy transferred among surfaces.  In order to simplify computations, the method assumes that all scattering is [[Lambert's cosine law|perfectly diffuse]].  Surfaces are typically discretized into quadrilateral or triangular [[finite element method|elements]] over which a piecewise polynomial function is defined.

After this breakdown, the amount of light energy transfer can be computed by using the known reflectivity of the reflecting patch, combined with the [[view factor]] of the two patches. This [[Dimensionless number|dimensionless quantity]] is computed from the geometric orientation of two patches, and can be thought of as the fraction of the total possible emitting area of the first patch which is covered by the second.

More correctly, radiosity ''B'' is the energy per unit area leaving the patch surface per discrete time interval and is the combination of emitted and reflected energy:

:<math>B(x)\, dA = E(x) \, dA +\rho(x) \, dA \int_{S}B(x') \frac{1}{\pi r^2} \cos\theta_x\cos\theta_{x'} \cdot \mathrm{Vis}(x,x') \,\mathrm dA'</math>

where:
* ''B(x)''<sub>i</sub> d''A''<sub>i</sub> is the total energy leaving a small area d''A''<sub>i</sub> around a point ''x''.
* ''E(x)''<sub>i</sub> d''A''<sub>i</sub> is the emitted energy.
* ''&rho;(x)'' is the reflectivity of the point, giving reflected energy per unit area by multiplying by the incident energy per unit area (the total energy which arrives from other patches).
* ''S'' denotes that the integration variable ''x' '' runs over all the surfaces in the scene
* ''r'' is the distance between ''x'' and ''x' ''
* θ<sub>x</sub> and θ<sub>x'</sub> are the angles between the line joining ''x'' and ''x' '' and vectors normal to the surface at ''x'' and ''x' '' respectively.
* Vis(''x'',''x' '') is a visibility function, defined to be 1 if the two points ''x'' and ''x' '' are visible from each other, and 0 if they are not.

[[File:Nusselt analog.svg|thumb|right|The geometrical form factor (or "projected solid angle") ''F''<sub>ij</sub>.<br /><br /> ''F''<sub>ij</sub> can be obtained by projecting the element ''A''<sub>j</sub> onto the surface of a unit hemisphere, and then projecting that in turn onto a unit circle around the point of interest in the plane of ''A''<sub>i</sub>. The form factor is then equal to the proportion of the unit circle covered by this projection.<br /><br />Form factors obey the reciprocity relation ''A''<sub>i</sub>''F''<sub>ij</sub> = ''A''<sub>j</sub>''F''<sub>ji</sub>]]
If the surfaces are approximated by a finite number of planar patches, each of which is taken to have a constant radiosity ''B<sub>i</sub>'' and reflectivity ''&rho;<sub>i</sub>'', the above equation gives the discrete radiosity equation,

:<math>B_i = E_i + \rho_i \sum_{j=1}^n F_{ij} B_j</math>

where ''F<sub>ij</sub>'' is the geometrical [[view factor]] for the radiation leaving ''j'' and hitting patch ''i''.

This equation can then be applied to each patch. The equation is monochromatic, so color radiosity rendering requires calculation for each of the required colors.

=== Solution methods ===
The equation can formally be solved as matrix equation, to give the vector solution:
:<math>B = (I - \rho F)^{-1} E \;</math>
This gives the full "infinite bounce" solution for B directly.  However the number of calculations to compute the matrix solution scales according to ''n''<sup>3</sup>, where ''n'' is the number of patches.  This becomes prohibitive for realistically large values of ''n''.

Instead, the equation can more readily be solved iteratively, by repeatedly applying the single-bounce update formula above.  Formally, this is a solution of the matrix equation by [[Jacobi iteration]].  Because the reflectivities ρ<sub>i</sub> are less than 1, this scheme converges quickly, typically requiring only a handful of iterations to produce a reasonable solution.  Other standard iterative methods for matrix equation solutions can also be used, for example the [[Gauss–Seidel method]], where updated values for each patch are used in the calculation as soon as they are computed, rather than all being updated synchronously at the end of each sweep.  The solution can also be tweaked to iterate over each of the sending elements in turn in its main outermost loop for each update, rather than each of the receiving patches.  This is known as the ''shooting'' variant of the algorithm, as opposed to the ''gathering'' variant.  Using the view factor reciprocity, ''A''<sub>i</sub> ''F''<sub>ij</sub> = ''A''<sub>j</sub> ''F''<sub>ji</sub>, the update equation can also be re-written in terms of the view factor ''F''<sub>ji</sub> seen by each ''sending'' patch ''A''<sub>j</sub>:
:<math>A_i B_i = A_i E_i + \rho_i \sum_{j=1}^n A_j B_j F_{ji}</math>
This is sometimes known as the "power" formulation, since it is now the total transmitted power of each element that is being updated, rather than its radiosity.

The [[view factor]] ''F''<sub>ij</sub> itself can be calculated in a number of ways. Early methods used a ''[[Hemicube (computer graphics)|hemicube]]'' (an imaginary cube centered upon the first surface to which the second surface was projected, devised by [[Michael F. Cohen]] and [[Donald P. Greenberg]] in 1985). The surface of the hemicube was divided into pixel-like squares, for each of which a view factor can be readily calculated analytically.  The full form factor could then be approximated by adding up the contribution from each of the pixel-like squares. The projection onto the hemicube, which could be adapted from standard methods for determining the visibility of polygons, also solved the problem of intervening patches partially obscuring those behind.

However all this was quite [[computation]]ally expensive, because ideally [[form factor (radiative transfer)|form factor]]s must be derived for every possible pair of patches, leading to a [[Quadratic function|quadratic]] increase in computation as the number of patches increased. This can be reduced somewhat by using a [[Binary space partitioning|binary space partitioning tree]] to reduce the amount of time spent determining which patches are completely hidden from others in complex scenes; but even so, the time spent to determine the form factor still typically scales as ''n'' log ''n''.  New methods include adaptive integration.<ref>G Walton, ''Calculation of Obstructed View Factors by Adaptive Integration'', [http://www.bfrl.nist.gov/IAQanalysis/docs/NISTIR-6925.pdf NIST Report NISTIR-6925], see also http://view3d.sourceforge.net/</ref>

=== Sampling approaches ===
The form factors ''F''<sub>ij</sub> themselves are not in fact explicitly needed in either of the update equations; neither to estimate the total intensity Σ<sub>j</sub> ''F''<sub>ij</sub> ''B''<sub>j</sub> gathered from the whole view, nor to estimate how the power ''A''<sub>j</sub> ''B''<sub>j</sub> being radiated is distributed.  Instead, these updates can be estimated by sampling methods, without ever having to calculate form factors explicitly.  Since the mid 1990s such sampling approaches have been the methods most predominantly used for practical radiosity calculations.

The gathered intensity can be estimated by generating a set of samples in the unit circle, lifting these onto the hemisphere, and then seeing what was the radiosity of the element that a ray incoming in that direction would have originated on.  The estimate for the total gathered intensity is then just the average of the radiosities discovered by each ray.  Similarly, in the power formulation, power can be distributed by generating a set of rays from the radiating element in the same way, and spreading the power to be distributed equally between each element a ray hits.

This is essentially the same distribution that a [[path-tracing]] program would sample in tracing back one diffuse reflection step; or that a bidirectional ray-tracing program would sample to achieve one forward diffuse reflection step when light source mapping forwards.  The sampling approach therefore to some extent represents a convergence between the two techniques, the key difference remaining that the radiosity technique aims to build up a sufficiently accurate map of the radiance of all the surfaces in the scene, rather than just a representation of the current view.