\documentclass[a4paper,11pt]{article}

% Definitions for the document
\def \Course {Interactive Computer Graphics }
\def \LectureNo {14}
% font definitions used in vector algebra
\newcommand{\Vector}[1]{\mathbf{#1}}
\newcommand{\UVector}[1]{\textbf{\em #1}}
\newcommand{\Matrix}[1]{\mathcal{#1}}
% standard layout for Duncan's lecture notes and tutorial material
\usepackage[margin=2cm, nohead]{geometry}
\usepackage{times}
\usepackage{wrapfig}
\usepackage{subfigure}
\usepackage{graphicx}
\usepackage{fancyhdr}
\usepackage{amsmath}
\usepackage{amsfonts}
\pagestyle{fancy}
\fancyhead{}
\fancyfoot{}
\renewcommand{\footrulewidth}{0.4pt}
\renewcommand{\headrulewidth}{0pt}
\lfoot{\footnotesize{\Course Lecture \LectureNo}}
\rfoot{\footnotesize{\thepage}}
% end of standard layout

\begin{document}

\section*{Lecture \LectureNo: Computational Issues in Radiosity}

%\begin{wrapfigure}{r}{90mm}
%  \vspace{-10mm}
%  \begin{center}
%     \includegraphics[width=90mm]{Figures/GraphicsFigure01_01.jpg} 
%  \end{center}
%  \vspace{-5mm}
%  \label{fig:GraphicsFigure01_01}
%  \caption{The first figure}
%\end{wrapfigure}
The radiosity method has stages that require intensive use of cpu time, and involve approximate methods that can lead to unwanted visual artefacts. The problems that are associated with each stage of the method are sumarised as follows:
\begin{enumerate}
\item Divide the graphics world into discrete patches \\
\indent Meshing strategies, meshing errors
\item Compute form factors by the hemicube method \\
\indent Alias errors
\item Solve the matrix equation for the radiosity of each patch. \\
\indent Computational strategies
\item Average the radiosity values at the corners of each patch \\
\indent Interpolation approximations
\item Compute a texture map of each point or render directly \\
\indent At least this stage is relatively easy
\end{enumerate}

\subsection*{Alias Errors}
Computation of the form factors will involve alias errors. These are equivalent to those errors that occur in texture mapping, due to discrete sampling of a continuous environment. However, these errors are perhaps the least of the problems we will encounter. Errors in the form factors are a secondary effect, as each patch radiosity will be determined by a large number of other patches. Thus each form factor will make only a small contribution to the result.  The alias errors can be reduced, if necessary, by increasing the sampling of the hemi-cube, but this in turn puts up the computational demand. 

\subsection*{Form Factor Computation}
Form factors have a reciprocal relationship:
\begin{equation}
\nonumber
F_{ij} = \frac{ cos\phi_i \; cos\phi_j |A_j|}{\pi r^2} 
\end{equation}
\begin{equation}
\nonumber
F_{ji} = \frac{ cos\phi_i \; cos\phi_j |A_i|}{\pi r^2} 
\end{equation}
\begin{equation}
\nonumber
F_{ji} = \frac{ F_{ij} |A_i|}{|A_j|} 
\end{equation}
Thus, providing we can store the form factors, only half the patches need be computed. However, the number of form factors will be very large, and in early solutions for radiosity it was not possible to pre-compute and save them. To give an idea of the size of the problem, we note that for 60,000 patches, there are 3,600,000,000 form factors. We only need to store half of these (reciprocity), but we will need four bytes for each, hence 7Gbytes are needed. This memory requirement is going up with the square of the number of patches. As many of them are zero we can save space by using an indexing scheme. An index could be created with one bit per form factor. If a bit is zero this implies that the form factor zero and not stored. All non zero form factors can be stored in a one dimensional array. Using such a scheme reduces the storage requirement, possibly to one quarter. However, it is easy to conceive of scenes where it is not possible to pre-compute and store the form factors.

\subsection*{Matrix Inversion}
Inverting the matrix can be done by the Gauss Siedel method. This is not really a matrix method as such, but rather an iterative scheme based on the form factor equations. Each patch equation has the form:
\begin{equation}
\nonumber
 B_i = E_i +R_i \sum B_j F_{ij}
\end{equation}
We use the iteration
\begin{equation}
\nonumber
 B_i^k = E_i +R_i \sum B_j^{k-1} F_{ij}
\end{equation}
The initial values, $B_i^0$, may be set to zero. At the first iteration the emitted light energy is distributed to those patches that are illuminated, in the next cycle, those patches illuminate others and so on. The Gauss Siedel inversion is stable and converges fast since the $E_i$ terms are constant and correct at every iteration, and all $B_i$ values are positive.

\subsection*{Progressive refinement}
The nature of the Gauss Siedel method allows a partial solution to be rendered as the computation proceeds. Without altering the method we could render the image after each iteration, allowing the designer to stop the process and make corrections quickly. This may be particularly important if the scene is so large that we need to re-calculate the form factors every time we need them.

If the image is progressively refined in this way, it will start dark and gradually illuminate as the iterations proceed. For this reason it is sometimes thought advantageous to add in a constant ambient term which is reduced as each iterative epoch is completed. 

\subsection*{Shooting and Gathering}

\begin{wrapfigure}{r}{30mm}
 \vspace{-5mm}
 \begin{center}
    \includegraphics[width=30mm]{Figures/GraphicsRadiosityGatheringPatch.jpg} 
 \end{center}
 \vspace{-5mm}
\end{wrapfigure}
The Gauss Siedel inversion can be modified to make it faster by making use of the fact that it is essentially distributing energy around the scene.  By choosing a good order of evaluation we can improve the quality of a partial solution and reduce the number iterations required for convergence.

The evaluation of one $B$ value using one row of the matrix in the normal Gauss Siedel way is like a process of gathering. The value of $B$ is estimated from the $B$ values of all the other patches:
 \begin{equation}
\nonumber
 B_i^k = E_i +R_i \sum B_j^{k-1} F_{ij}
\end{equation}

\begin{wrapfigure}{r}{30mm}
 \vspace{-10mm}
 \begin{center}
    \includegraphics[width=30mm]{Figures/GraphicsRadiosityShootingPatch.jpg} 
 \end{center}
 \vspace{-5mm}
\end{wrapfigure}
Suppose, in one iteration, $B_i$ changes by $\Delta B_i$. Every other patch in the scene will change (providing its form factor is not zero). The magnitude of this change can be written:
\begin{equation}
\nonumber
 \Delta B_j = R_jF_{ij} \Delta B_i
\end{equation}        
Which can be simply changed into an iterative updating scheme:
\begin{equation}
\nonumber
  B_j^k = B_j^{k-1} + R_jF_{ij} \Delta B_i^{k-1}
\end{equation} 
This is the process of shooting, and is equivalent to evaluating the matrix column wise. The use of shooting allows us to choose an evaluation order that ensures fastest convergence. The patches with the largest change $\Delta B$ (called the unshot radiosiy) are distributed first. In an iteration, the unshot radiosity of a chosen patch is reduced to zero, and that of other patches is incrementally increased. The process starts by shooting the emitting patches, since at the first iteration $\Delta B_i = Ei$ . 

\subsection*{Interpolation Strategies}
Visual artefacts do occur with interpolation strategies, but may not be significant for small patches. Figure~\ref{fig:Interpolation} illustrates the different interpolation strategies and the results they will produce. As noted previously the computed values alone will produce undesirable visual artefacts. However, the use of a cubic interpolation scheme, though more accurate than a purely linear scheme, will probably not produce any significant improvement.
\begin{figure}
\begin{minipage}[b]{0.48\linewidth}
\centering
\includegraphics[height=43mm]{Figures/GraphicsRadiosityInterpolation.jpg}
\caption{Interpolation Strategies}
\label{fig:Interpolation}
\end{minipage}
\hspace{5mm}
\begin{minipage}[b]{0.48\linewidth}
\centering
\includegraphics[height=40mm]{Figures/GraphicsRadiosityMeshing.jpg}
\caption{$D^0$ Meshing Problems}
\label{fig:Meshing}
\end{minipage}
\end{figure}
\subsection*{Meshing}
Meshing is the process of dividing the scene into patches, and causes the worst problems in terms of visual quality. Meshing artefacts are scene dependent. The most obvious are called $D^0$ artefacts, and are caused by caused by discontinuities in the radiosity function. These could be due to shadows or to surface albedo (ratio of incident to reflected light) or texture. Discontinuities of this kind cause problems when accompanied by bad patching. Figure~\ref{fig:Meshing} illustrates the problem. Here a single polygon is divided into regular square patches. A shadow boundary crosses, resulting in a serious visual alias in the computed radiosity. Interpolation will reduce the visual effect of the alias, but will smooth out the hard edge of the shadow.

\subsection*{Discontinuity Meshing (a priori)}
The idea behind discontinuity meshing is to compute discontinuities in advance, and align the patches with them. Additionally, the discontinuities are made known to the interpolation procedure so that hard boundaries are not smoothed out. Many of the discontinuities are scene dependent. These in particular will include object boundaries, which we would naturally expect to be patch boundaries as well, and albedo discontinuities, possibly specified in texture maps. Shadows can also be determined, and this can be done to a process similar to ray tracing.  Books contain a lot of algorithmic detail on this subject, but we will not cover it here.

\subsection*{Adaptive Meshing (a posteriori)}
The idea of adaptive meshing is to recompute the mesh as the radiosity calculation proceeds. Places where there are strong discontinuities in the radiosity can be found by comparing values at adjacent patches. There are essentially two approaches:
\begin{enumerate}
\item put more patches (elements) into areas with high discontinuity or
\item move the mesh boundary to coincide with the greatest change.
\end{enumerate}

\subsection*{Subdivision of Patches (h refinement)}

Using the h-refinement method we compute the radiosity at the vertices of the coarse grid of patches, as shown in Figure~\ref{fig:HRefinement}. If the discontinuities exceed a threshold we subdivide the patch into elements, and recompute those elements. They can be checked for discontinuity and subdivided further. When a patch is divided into elements the radiosity of each element is computed using the original radiosity solution for all other patches. The assumptions behind this are that: 
\begin{enumerate}
\item the radiosity of a patch is equal to the sum of the radiosity of its elements, and, 
\item the distribution of radiosities among elements of a patch do not affect the global solution significantly.
\end{enumerate}
These assumptions are a reasonable, and they save re-computing the complete radiosity solution after each subdivision.
\begin{figure}[ht]
\begin{minipage}[b]{0.45\linewidth}
\centering
\includegraphics[height=50mm]{Figures/GraphicsRadiosityHRefinement.jpg}
\caption{Adaptive Meshing using h-refinement}
\label{fig:HRefinement}
\end{minipage}
\hspace{0.5cm}
\begin{minipage}[b]{0.45\linewidth}
\centering
\includegraphics[height=50mm]{Figures/GraphicsRadiosityRRefinement.jpg}
\caption{Adaptive Meshing using r-refinement}
\label{fig:RRefinement}
\end{minipage}
\end{figure}
\subsection*{Patch Refinement (r refinement)}
Here we compute the radiosity at the vertices of the coarse grid, and move the patch boundaries closer together if they have high radiosity changes. Unlike the h-refinement solution it is necessary to recompute the entire radiosity solution each refinement. However the method should make more efficient use of patches by shaping them correctly. Hence a smaller number of patches could be used. In Figure~\ref{fig:RRefinement} we see that the uniform patches tend to get larger, and the ones around the discontinuities tend to get smaller, hence some other forms of re-patching may be required as the algorithm proceeds.

\subsection*{Adding Specularities}
We noted that specularities (being viewpoint dependent) cannot be calculated by the standard radiosity method. However, they do form an important aspect of visual realism, and it is desirable to add them where possible. This can be done after the radiosity solution has been obtained, by means of ray tracing. The complete ray tracing solution is not required, just the specular component in the viewpoint direction. Thus only one primary ray per pixel and one secondary ray in the reflected direction are required to achieve this. Further ray tracing can be used to enhance radiosity solutions with effects such as reflections and translucency.


\end{document}