In contrast to previous cases, the application of harmonic functions does not require the subdivision of surfaces into planar polygons, but deals with the original geometry. This property makes it especially useful when the view-dependent rendering phase uses ray-tracing.
Suppose surface A is defined parametrically by a
position vector function, 
, where parameters u and v are in
the range of [0,1].
Let a representative of the basis functions be:
( 
substitutes 
for notational simplicity). Note that the basis functions have two indices,
hence the sums should also be replaced by double summation in
equation 1.79. Examining the basis
functions carefully, we can see that the goal is the calculation of the
Fourier series of the radiosity distribution.
In contrast to the finite element method, the basis functions are now non-zero almost everywhere in the domain, so they can approximate the radiosity distribution in a wider range. For that reason, approaches applying this kind of basis function are called global element methods.
In the radiosity method the most time consuming step is the evaluation of the integrals appearing as coefficients of the linear equation system (equation 1.79). By the application of cosine functions, however, the computational time can be reduced significantly, because of the orthogonal properties of the trigonometric functions, and also by taking advantage of effective algorithms, such as Fast Fourier Transform (FFT).
In order to illustrate the idea, the calculation of
for each k,l is discussed.
Since 
, it can be regarded as a function defined
over the square 
. Using the equalities of surface integrals, and
introducing the notation
for surface element magnification, we get:
Let us mirror the function
onto coordinate system axes u and v, and repeat the resulting
function having its domain in 
infinitely in both directions with period 2. Due to mirroring and periodic
repetition, the final function 
will be even and periodic with
period 2 in both
directions. According to the theory of the Fourier series, the function can be
approximated by the following sum:
All the Fourier coefficients 
can be calculated by a
single, two-dimensional FFT. (A D-dimensional FFT of N samples
can be computed by taking
number of one-dimensional FFTs [Nus82] [PFTV88].)
Since 
if 
,
this Fourier series and the definition of the basis functions
can be applied to equation 1.98, resulting in:
Consequently, the integral can be calculated in closed form, having replaced the original function by Fourier series. Similar methods can be used to evaluate the other integrals. In order to compute
J(u,v) must be Fast Fourier Transformed.
To calculate
the Fourier transform of
is needed. Unfortunately the latter requires a 4D FFT which involves many operations. Nevertheless, this transform can be realized by two two-dimensional FFTs if g(p,p') can be assumed to be nearly independent of either p or p', or it can be approximated by a product form of p and p' independent functions.
Finally, it should be mentioned that other global function bases can also be useful. For example, Chebyshev polynomials are effective in approximation, and similar techniques to FFT can be developed for their computation.
