Annette Scheel
email: scheel@mpi-sb.mpg.de

Marc Stamminger
email: stamminger@mpi-sb.mpg.de

Jörg P Ütz

Hans-Peter Seidel
email: hpseidel¢@mpi-sb.mpg.de

Enhancements to Directional Coherence Maps

Keywords: Directional Coherence Maps, Rendering, Interpolation

Directional Coherence Maps as proposed by Guo in '98 are a very efficient acceleration technique for ray tracing based global illumination renderers. It vastly reduces the number of pixels which have to be computed exactly by identifying regions which are suitable for interpolation. By using oriented finite elements for interpolation, the sampling density can be kept low for large regions of the target image. In this paper we describe extensions of the method. An improved object test is presented which prevents that small objects are missed. Additionally, it is shown how to handle textures efficiently, which was not possible with the original approach.

Introduction

Creating a high-quality image from a global illumination solution can take hours for complex scenes if every pixel is evaluated with high accuracy. Therefore, several methods have been developed which decide which pixels can be approximated by interpolation.

An early, widely used strategy for adaptive refinement for ray tracing was developed by [START_REF] Painter | Antialiased ray tracing by adaptive progressive refinement[END_REF]. They progressively add new image samples and store them in a k-D tree. The decision where to place new samples is based on the variance of the samples (to detect features) and on the area which is covered by them (to cover under-populated areas). They use a piecewise constant reconstruction from the cells of the k-D tree plus filtering. A large number of methods followed that tried to improve the strategies for finding discontinuities both in image and in object space and to exploit this knowledge for faster image synthesis. This includes for example methods for finding shadows [START_REF] Chin | Near real-time shadow generation using bsp trees[END_REF][START_REF] Teller | Computing the antipenumbra of an area light source[END_REF][START_REF] Drettakis | A fast shadow algorithm for area light sources using backprojection[END_REF][START_REF] Stewart | Fast computation of shadow boundaries using spatial coherence and backprojections[END_REF][START_REF] Hart | Direct illumination with lazy visibility evaluation[END_REF], discontinuity meshing [START_REF] Heckbert | Discontinuity meshing for radiosity[END_REF][START_REF] Lischinski | Combining hierarchical radiosity and discontinuity meshing[END_REF], subdividing image space according to object and shadow boundaries [START_REF] Pighin | Progressive previewing of ray-traced images using image-plane discontinuity meshing[END_REF], and adaptive interpolation of indirect light [START_REF] Ward | Irradiance gradients[END_REF].

£ currently at iMagis-GRAVIR/INRIA, Sophia-Antipolis, France

The Directional Coherence Maps of Guo [START_REF] Guo | Progressive radiance evaluation using directional coherence maps[END_REF] provide an interpolation method even for finer details of the solution. This method tries to find the main orientation of shading details (a shadow or a glossy highlight) and object boundaries and interpolates along those discontinuities by using so-called oriented finite elements. In contrast to discontinuity meshing the performance of the method is not constrained by the number and type of lightsources and complex geometry calculations are avoided. The DCM reconstruction can be used in conjunction with every global illumination algorithm which is able to compute exact results for single pixels.

The time to compute a high-quality image is usually very long; therefore it is desirable to have already meaningful early images. With progressing computation the solution should improve and finally converge to the correct solution [START_REF] Kajiya | The rendering equation[END_REF]. The DCM incorporate such a progressive refinement as well by using a hierarchical block refinement technique.

In the following, first the DCM algorithm will be described in more detail (Section 2) and then our improvements for this method will be given (Section 3).

The goal of the Directional Coherence Maps by [START_REF] Guo | Progressive radiance evaluation using directional coherence maps[END_REF] is to compute as few image points as possible exactly and to determine the other points by using some interpolation scheme.

Of course anti-aliasing is an issue for all ray-tracing based rendering systems. In order to keep things simple, we will assume in the following that anti-aliased images are obtained by creating images of multiple resolution and scaling them down to the desired image size. In the following a pixel is an image point of this enlarged high resolution image, several such pixels are averaged to obtained a single pixel value of the final image.

Overview

The image space is first subdivided into so called elementary blocks which are small squares, typically of size 8x8 pixels. At the beginning of the DCM reconstruction each elementary block is classified either as a smooth block or an edge block (see Figure 1 a) and b) and 2.2 for a description of the classification procedure). In contrast, edge blocks contain one or more discontinuities and simple bilinear interpolation would not be sufficient in this case. For those blocks all pixels on the block's border are evaluated exactly. With this information the direction of the "strongest" discontinuity is found. This is done by determining the least discrepancy direction, i.e. the direction for which the differences between border pixels that are opposite with respect to this direction are minimal (see Figure 2). Typically, eight different directions are tested. Parallel to this discontinuity so called oriented finite elements are established, which are longish quadratic polygons (see Figure 2). The interior of those elements is obtained by linear interpolation between the border values. In order not to miss any features, the four children of a smooth block undergo the same classification into smooth or edge blocks as elementary blocks (see 2.2). For edge blocks another test (see 2.3) is done before they are subdivided into four quads in order to determine if they contain only one single discontinuity which is only slightly curved (so-called simple edge blocks) or some more complex features (complex edge blocks) (see Figure 1 b) and c)). For children of complex edge blocks the oriented finite elements have to be re-calculated, whereas for simple edge blocks only a lazy boundary evaluation is done.

Here, only the center pixel (of the parent) is computed exactly and this value is compared to the result of the finite element approximation. If this does not differ significantly the old orientation of the finite elements is kept, the interpolation is just corrected with the new values. Otherwise, a new full boundary evaluation has to be done.

Classification Smooth-Edge Block

Two tests are done to determine if a block is a smooth or an edge block. The first test is a contrast test: From the exact radiance values at the four vertices the maximum and minimum luminances are determined to compute the block contrast following Mitchells def-

inition of contrast: ¤ ¦¥ ¨ § © ¤ ¤ ¦¥ ¨ § ¤
. If this contrast is above some threshold, the block is classified as an edge block. For the other test, all visible lines from object boundaries (as well as their reflection in planar mirrors) are computed by creating an object-ID index image with OpenGL. If a block is crossed by a visible line, it is classified as an edge block, too.

Classification Simple-Complex Edge Block

A block has to pass three tests before it can be confirmed to be a simple edge block -otherwise it will be classified as a complex edge block. The first test (the so-called zero order test) tries to find the number of edges crossing this block by counting the significant changes in color along the border. If there are only two changes, it can be supposed that there is only one discontinuity (see Figure 3).

The next test, the first order test, attempts to find the gradient of the discontinuity for each of the two points, where the discontinuity intersects the block border. For this purpose, the pixels around the crossing point inside the border have to be evaluated. If the gradient is nearly the same at both crossing points, the discontinuity (probably) is a straight line (see Figure 3).

The last test uses the object IDs which were already used for the smooth-edge block classification. No more than two objects are allowed on the block border for a simple edge block. In our implementation of the algorithm of Guo we achieved similar results as were described in the original paper (see Figures 4 and7). Depending on the scene's character only 6-10% of the pixels had to be evaluated exactly without a noticeable loss of quality. Most images were produced with oversampling (4 or 16 times) to obtain anti-aliased results. For large images (i.e. screen size) it might be possible that the data used for the DCM does not fit into main memory. We therefore implemented a possibility to subdivide the image into smaller parts and to do the algorithm for all sub-images independently.

The remainder of this section describes enhancements we introduced in order to handle thin objects on the finest block refinement level (see 3.1) and to improve the object tests which are used in several occasions during the algorithm (see 3.2). Furthermore, because the original algorithm was not able to handle textures efficiently, we developed an enhancement for textures which will be described in 3.3.

Handling Thin Objects

The DCM algorithm is designed to create progressively finer and better images, theoretically starting with a single elementary block containing the entire image and finally reaching pixel level. In practice it is a good starting point to use elementary blocks of size 8x8 (corresponding to Guo and our experience, too). In our implementation we then refine only once to obtain 4x4 blocks. For the next finer block size of 2x2 it would not make sense to do the complex finite element orientation procedure.

Especially very thin and longish objects like the shaft of the desk lamp of Figure 4 or the pricks of the cactus in Figure 7 which only have a thickness of 1 or 2 pixels are difficult to represent by the finite element interpolation. Blocks containing those objects will be classified as (complex) edge blocks due to the object ID tests. The finite element orientation procedure tries then to find the least discrepancy direction. Because the thin line only covers very few border pixels it will hardly contribute to overall differences. Especially in the presence of other shading details (like shadow) having a different orientation than the line the finite elements will not be oriented into the direction of the thin line. Therefore, the thin line will be washed out by the interpolation with this scheme (cp. Figure 5, left image).

We therefore use another test after the finite elements of a 4x4 block have been established. Before a value of an interior pixel is determined by interpolation it is tested if the object IDs of the border pixels correspond to the object ID of this pixel (the object IDs are already known from the classification tests 2.2).

If not, this single pixel is computed exactly. on the bottom shows which pixels are computed additionally (those not having a block structure). The total increase in exactly computed samples was only about 1% for this scene.

The images of

Improved Object Test

The original algorithm uses object IDs for the smooth-edge block as well for the simple-complex edge block classification. However, often more complex objects are composed from simple drawing primitives like triangles all consisting of the same material. In this case, it is desirable to interpolate across those primitives in order not to stress the borders between them. Therefore, we additionally compare the shaders of the objects. Only if object ID and shader differ, blocks will be classified as edge or complex edge blocks, respectively. Because shaders can be of arbitrary complexity, we decided to compare only shader IDs.

Enhancement for Textures

Texturing is a very efficient and widely used means for increasing the visual richness of synthetic renderings. By mapping images onto geometric objects, even simple objects can offer high visual detail in the final rendering, almost without increase of overall rendering time.

Such high visual detail from textures, however, re-tards the method of Guo. The method considers the texture as lighting detail, which entails fine sampling of the entire object. Compare Figure 6, where the floor is a single polygon with a parquet texture. It exhibits fine detail both due to the parquet texture and due to lighting, in this case mainly the chair's shadow.

However, texture detail is usually produced very lately, namely at the last light bounce before the light reaches the eye (unless you look onto textured surfaces via a mirror for example). It is not the incident light at the object that changes rapidly, but just its reflection.

This distinction opens the way for improving Guo's method in the presence of textured objects: if we could interpolate the incident light at the object with Guo's method and perform the reflection computation including the textures afterwards, in the above example the method can well interpolate large regions of the floor and can concentrate on the critical shadow boundaries and the glossy reflection. The right column of Figure 6 shows the result of our approach based on this idea. Although the final image has rich detail allover due to the textured floor, only the shadow boundaries and glossy reflection boundaries were sampled finely. In contrast, the left column of Figure 6 shows that the textured floor is sampled densely when the original approach is used.

This approach was chosen by [START_REF] Pighin | Progressive previewing of ray-traced images using image-plane discontinuity meshing[END_REF], but interpolating incident light is only possible if the last reflection is diffuse. To be able to handle more complex BRDFs we only delay the texturing and thus the most rapidly changing part of reflection until the final interpolation. More precisely, we split the reflected light of a pixel into two parts: the first part ! is the one that has been modulated by the object's texture, the second one #" is the constant offset unaffected by the texture. Note that e.g. for plastic a texture only modulates the diffuse reflection of an object ($ ), whereas the glossy reflection does not vary ($ "

).

Instead of !

, we consider the value before texturing, i.e. the value &% ' )( 10 2 43 65 87 @9

. Altogether, if we have an image sample on an object at texture coordinates 2 A3 B5 )7 C9

, and the texture function on the object is 0 2 43 65 87 @9 , the resulting color of the pixel is % ED " % 0 2 A3 B5 )7 C9 D " 5

The gain of this separation is that we now have two colors and "

that are about as smooth as if the scene was not textured at all. If we apply Guo's approach to this color tuple 2 4 5 8 F" G9 instead of to , less samples will become necessary.

On the downside, the interpolation becomes more complicated. First, we have to consider two colors, i.e. six color channels, for interpolation, contrast computation etc. Second, we have to evaluate the texture function 0 2 43 65 87 @9 during interpolation, which means that we have to compute the texture coordinates 2 A3 B5 )7 C9 for an interpolated point and have to evaluate the texture to obtain the modulation factor.

Six color channels instead of three is just a small memory problem. More difficult is the computation of texture coordinates. But since we know from the tests using object IDs for each interpolation point which object it sees, we can quickly intersect the corresponding eye ray with the object to obtain the texture coordinates.

Results and Conclusion

The image reconstruction method of Guo turned out to be a very efficient means for the acceleration of ray-tracing based renderers. Image discontinuities are detected reliably, and with the oriented finite elements simple discontinuities are approximated very efficiently. In particular if supersampling is used for anti-aliasing and the final images are obtained by averaging the reconstruction result, the quality loss is low.

In the following table we compared the method of Guo with a simple baseline renderer for the scene in Figure 4. The computation time was split into the time needed for pure illumination computation and for the DCM. The table shows that the time overhead introduced by the DCM is rather small but the overall time gain is high, although we used only a very simple illumination computation algorithm (300 rays are shot to the lightsources, for glossy objects 20 reflection rays are shot additionally). For a more sophisticated illumination algorithm (like for example a final gather of a finite element solution) even a bigger time gain can be expected.

Baseline vs. DCM (for Figure 4) Method Pixel Time DCM Time Illum Baseline 100 % -3672 s DCM 5.9 % 43 s 245 s

With our improved testing scheme, most complex cases where the original method fails are detected. By switching back to a full sampling for these few critical cases many small artifacts can be avoided by only a slight decrease of rendering performance. As a result, the quality of the final images, although obtained from usually less than 10% of the samples, is close to that of a fully sampled image (see Figure 7).

The table below contains a comparison of Guo's original method (refined to 4x4) and the extension for thin objects for the scene in Figure 7. The higher quality of the image was achieved by the expense of only approximately 1% more samples. The time for the DCM computation is moderately increased because of the additional test during the interpolation. The scene is only directly lit by point light sources, therefore the DCM time exceeds the time for the illumination computation. Again, with a more complex lighting algorithm the time for illumination computation will dominate over the time needed for the DCM.

Thin Objects (see Figure 7 6 and7).

In the following table we compared the improved method with Guo's original method for the scene in Figure 6. It shows that the number of exactly computed samples and hence also the computation time could be strongly reduced, while the increase in the time needed for the DCM computation is minimal. Again, with a more complex lighting algorithm an even higher time gain can be expected. On the downside, the method suffers from noisy illumination computations. The noise cannot be distinguished from 'real' detail and results in too fine sampling of the noisy regions. We are not aware of a simple solution to this problem.

Figure 1 :

 1 Figure 1: Block types

Figure 2 :

 2 Figure 2: Finding the least discrepancy direction and drawing the oriented finite elements

Figure 3 :

 3 Figure 3: Zero and first order tests

  Figures 4 and 7 were computed with this extension. The right image of Figure 5 contains a detail of Figure 7. Especially the pricks of the cactus are reproduced better than in the left image which was created by the original model. The index image

Figure 4 :

 4 Figure 4: Size 512x512, 4 times oversampling, 6% of pixels are computed exactly. The bottom images show two details of the index image (desk lamp and chair).

Figure 5 :

 5 Figure 5: Detail of the cactus. Left: original method (refined to 4x4 blocks), Right: new method, Bottom: Enlarged part of the index for the new method.

Figure 6 :

 6 Figure 6: Left: without texture extension, Right: with texture extension

  Furthermore, we could solve the problem of very dense sampling enforced by textured surfaces. By delaying the texturing to the final interpolation step, textured surfaces are handled without an increase in sampling density (see Figures

			)
	Method	Pixel Time DCM Time Illum
	Guo	7.28 %	70.82 s	52.38 s
	New	8.17 %	94.16 s	55.42 s

Acknowledgments

This research was funded by the ESPRIT Open LTR Project #35772 SIMULGEN "Simulation of Light for General Environments".


Annette Scheel

Marc Stamminger¡

Hans-Peter Seidel

Thrifty Final Gather for Radiosity

Finite Element methods are well suited to the computation of the light distribution in mostly diffuse scenes, but the resulting mesh is often far from optimal to accurately represent illumination. Shadow boundaries are hard to capture in the mesh, and the illumination may contain artifacts due to light transports at different mesh hierarchy levels. To render a high quality image a costly final gather reconstruction step is usually done, which re-evaluates the illumination integral for each pixel. In this paper an algorithm is presented which significantly speeds up the final gather by exploiting spatial and directional coherence information taken from the radiosity solution. Senders are classified, so that their contribution to a pixel is either interpolated from the radiosity solution or recomputed with an appropriate number of new samples. By interpolating this sampling pattern over the radiosity mesh, continuous solutions are obtained.

Introduction

In the past 15 years, much research has been concentrated on improvements of the radiosity method. The computation of radiosity solutions of more than one million patches is now possible on contemporary standard PCs. However, these solutions most often do not meet high-quality demands of many commercial applications. The problem is inherent to the method: because lighting detail is generated in object space, very fine tessellation is necessary to capture fine lighting detail. Thus quality is not only a time, but also a memory issue. Furthermore, the widely used linear Gouraud interpolation is prone to Mach banding, so that a human observer perceives the tessellation structure easily. In addition, long thin triangles can lead to frayed shadow boundaries, and finally the different levels at which light is transported to a patch may lead to artifacts. To some extent, these issues can be addressed by adapting subdivision to lighting discontinuities [START_REF] Lischinski | Combining hierarchical radiosity and discontinuity meshing[END_REF][START_REF] Volevich | Using the visible difference predictor to improve performance of progressive global illumination computation[END_REF] or by using higher-order interpolation (e.g. [START_REF] Gortler | Wavelet radiosity[END_REF]).

Ray-based Monte-Carlo methods usually compute the illumination at independent sample positions in image space. This point sampling allows exact lighting computations, but also makes the exploitation of coherence more difficult. As a result, illumination from big light sources or indirect light make the lighting computation expensive. For stochastic sampling as a rule 500 or 1000 samples per pixel are needed in such cases, and in some cases noise is still present.

The idea of final gathering is to combine both approaches. In a view-independent preprocess a not necessarily perfect radiosity solution is computed. Then in a second view-dependent step ray tracing is performed that recomputes parts of the illumination considered critical. One common technique is to only recompute the direct light in this pass and to add the indirect light from the finite element solution (e.g. [START_REF] Shirley | Physically Based Lighting Calculations for Computer Graphics[END_REF][START_REF] Dumont | A progressive algorithm for three point transport[END_REF]). Alternatively, the radiosity solution can be interpreted as a light source definition, and the last light bounce towards the eye is recomputed with ray tracing (e.g. [START_REF] Per | Clustering for glossy global illumination[END_REF]).

Previous Work

Compared to the huge number of publications on radiosity methods, the number of publications on a high quality rendering postprocess is surprisingly low. Often, a high quality rendering step is only a final chapter for a new radiosity algorithm, e.g. [START_REF] Lischinski | Combining hierarchical radiosity and discontinuity meshing[END_REF][START_REF] Smits | Efficient Hierarchical Radiosity in Complex Environments[END_REF][START_REF] Per | Clustering for glossy global illumination[END_REF][START_REF] Mostefaoui | Rendering inhomogeneous surfaces with radiosity[END_REF], but only very few papers address the issue in detail.

In [START_REF] Arjan | Adaptive sampling of area light sources in ray tracing including diffuse interreflection[END_REF], Kok et al. describe a method to decrease the number of shadow samples for area light sources exploiting information from a progressive radiosity pass. Their algorithm subdivides senders adaptively in order to detect completely visible or occluded parts that can be processed quickly. With a clever shadow pattern scheme, the number of samples is reduced further. To some extent, this radiosity postprocessing step is done implicitly in the radiosity computation itself by the idea of Hierarchical Radiosity [START_REF] Hanrahan | A rapid hierarchical radiosity algorithm[END_REF].

Several authors proposed the use of (hierarchical) radiosity for indirect light only, and to regenerate the usually more detailed direct illumination by ray tracing. More information from the radiosity solution is exploited in the local pass by Linschinski et al. [START_REF] Lischinski | Combining hierarchical radiosity and discontinuity meshing[END_REF]: Form factors are recomputed for each pixel and visibility for direct illumination only. Their approach is extended by Stürzlinger [START_REF] Stürzlinger | Optimized local pass using importance sampling[END_REF], where the number of samples is adapted to the relative contribution of the sources and visibility is computed in a stochastic approach. Christensen et al. [START_REF] Per | Clustering for glossy global illumination[END_REF] used a final gather step for reconstructing their radiance solutions. They gather all links with a fixed number of visibility samples, but also use the cluster hierarchy. It is pointed out that the final gather requires significant computation time and that despite their conservative resampling scheme artifacts still remain visible. Smits [START_REF] Smits | Efficient Hierarchical Radiosity in Complex Environments[END_REF] uses a fixed number of visibility tests and analytic form factor evaluations for each point. For links carrying a relative error below a certain threshold, the estimate of the link itself is taken. To avoid bias which might occur due to the threshold, Russian Roulette is used with links whose relative error is below the threshold.

The final gather step which is very briefly described for the algorithm presented in [START_REF] Mostefaoui | Rendering inhomogeneous surfaces with radiosity[END_REF] extends these ideas. For each receiver, critical senders are determined, whose contribution is to be computed exactly by final gather.

Bekaert et al. [START_REF] Bekaert | Final radiosity gather step using a monte-carlo technique with optimal importance sampling[END_REF] use the radiosity solution for defining an importance function for a Monte-Carlo ray tracing step. A similar idea later followed in [START_REF] Perez | Acceleration of monte carlo path tracing for general environments[END_REF], where a radiosity solution is used to guide a Monte-Carlo path tracer.

Finally, the photon map approach [START_REF] Jensen | Efficient simulation of light transport in scenes with participating media using photon maps[END_REF] is relevant, because it is essentially a final gathering approach, however it is based on a global photon tracing pass as preprocess. From the photon hits an approximate lighting solution can be computed quickly. The high quality images, however, are then obtained from a well optimized final gather step.

Idea

The idea of this paper is to better exploit the information from the global radiosity step for a more 'thrifty' final gather. Consider the example in Fig. 1, middle, showing a coarse radiosity solution from a candle, illuminated by a big light source on the left and a smaller one on the right. Due to the different sizes of the light sources, the shadows are very different. The left shadow has a clear outline that blurs with the distance to the object; the right shadow is blurry and completely washed out after a short distance.

The radiosity solution is obtained in a few seconds and thus of poor quality, but it is clear that for a final gather step valuable information can be obtained from it. In addition to the radiosity mesh which already captures the most important changes in illumination, the link structure contains per patch hints about shadow boundaries and the light sources responsible. Furthermore, during the radiosity computation more approximate information has been computed with each link, such as the variation of the form factor or a shadow gradient.

In previous work, this information was exploited by recomputing partially visible links only, or by using a fixed number of samples per link, thus implicitly performing importance sampling guided by the link structure resulting from radiosity. We experienced severe problems with these simple approaches: At patch boundaries, the set of sampled senders and their hierarchy level change abruptly. Even for very precise solutions, these discontinuities, although small in size, remain visible. For indirect light a large fraction of links is partially occluded, so resampling all of them is very expensive. When not resampling indirect light, scenes with dominant indirect light essentially show the radiosity solution with all its artifacts. Taking a fixed number of samples results in bad importance sampling. It is true that hierarchical radiosity creates a mesh of links with similar errors, so using a fixed number of samples per link assigns roughly similar error to each sample. However, the error measured in HR for refinement decision is not necessarily the error we are interested in for final gather. Using a minimal area parameter also breaks the uniform error assumption, in particular for critical links with high variation. Links arrive at different levels of the hierarchy, making comparisons difficult. Furthermore, the error of a link must be seen in the context of all other links arriving at the same point.

Our new method tries to address these deficiencies. For a given receiver we extract a list of potential senders and estimate whether their contribution needs to be resampled or can be interpolated from the radiosity solution. If resampling is necessary, an appropriate sample number is determined. Both decisions can take into account the magnitude of the illumination, its relation to other contributing patches, visibility gradient estimates, perceptual issues etc.

By gathering this sampling information at the vertices and interpolating it over the patches again, we obtain a continuous solution without discontinuities at the patch boundaries. Nevertheless, the solution will not necessarily have a continuous derivative, so Mach banding can still appear. However, the goal of the method is to resample the solution so accurately, that this effect should also disappear.

Our final gather accuracy is (almost) independent of the radiosity accuracy. If the radiosity solution is already sufficient, the final gather will simply interpolate the radiosity values without further sampling. On the other hand a coarse solution will lead to a more expensive final gather step. However, the method can fail in detecting lighting detail that has not been captured at all by the input radiosity solution.

Note that the goal of our method is to enhance the final gather step as such. This does not exclude the use of ray-tracing acceleration methods in image space, such as adaptive progressive refinement [START_REF] Painter | Antialiased ray tracing by adaptive progressive refinement[END_REF]. In fact, we use directional coherence maps [START_REF] Guo | Progressive radiance evaluation using directional coherence maps[END_REF] (see also results) to exploit image space coherence.

Algorithm

Our method is based on the results obtained by a hierarchical radiosity algorithm with clustering. In the final radiosity solution, every object is represented by a triangle mesh with an approximate radiosity value per triangle. The link mesh contains the approximate value of the form factor, its gradient, and visibility. The HR oracle is based on an estimated radiosity change over a receiver, multiplied by the receiver's area. The goal of the final gather step is to accurately recompute the last light bounce towards the camera. For now, we will handle diffuse reflection only.

Overview of the Algorithm

First, for each vertex under consideration, a list of contributing senders is created by gathering all links arriving at all levels and at all surrounding patches (Sect. 2.2 and Fig. 2, left). Then we determine for each sender, whether its contribution can be interpolated, based on the estimated magnitude of the change in illumination it will introduce (Sect. 2.3). Furthermore, we distinguish between variations caused by the form factor and the visibility term. For a low variation the corresponding term will be interpolated, otherwise an appropriate sample number for final gathering is computed and stored with the sender (Sect. 2.4 and Fig. 2, middle).

In order to perform final gathering at a visible object point, the surrounding vertices are considered. The contributions of all senders that are marked as 'smooth' in all vertices are interpolated from the vertices to the point with its barycentric coordinates. For all other senders, samples are spawned, where the number of samples is again interpolated from the vertices (Fig. 2, right). Sect. 2.5 describes our interpolation process that leads to a continuous final gather result.

In the following, we describe the single phases of the method in more detail, discuss alternatives and justify our choices.

Collection of Links for Vertices

In the radiosity solution links might arrive at different levels in the hierarchy at the receiver. Therefore, prior to the main algorithm, the links are 'pushed down' to leaf level on the receiver. During this push we recompute the form factor for the new links, but reuse the much more expensive visibility term. This recomputation avoids artifacts from the constant push step that typically appear when light is gathered at different hierarchy levels.

The algorithm begins by collecting the links from all patches surrounding a vertex, such that the vertex obtains a list of all senders. The information of the links is either averaged (for example for the irradiance) or the maximum is taken (for example for the visibility and form factor gradients). In this step links with visibility zero are neglected, which can be inaccurate if visibility is sampled by casting rays. However, a 'false reject' of a sender partially visible to the vertex requires that for all surrounding patches the sender was wrongly classified as completely occluded.

Classification of Senders

Usually, illumination exhibits strong coherence; hard-to-capture variations in form factor or visibility are only produced by a small portion of senders. These do not necessarily have to be light sources but can also be strongly illuminated patches creating strong indirect illumination. In either case, our algorithm tries to identify those "trouble makers" and resample them.

From the existing radiosity solution and from values computed for vertices information can be extracted to drive an oracle to detect senders with difficult contributions. Our oracle estimates the magnitude of the change in radiosity around a vertex due to form factor and visibility variations 1 . The threshold, which decides if this magnitude is small enough for interpolation, should respect human perception, because differences in color are not perceived equally well for all levels of illumination. In our implementation we use the CIELAB space [3]2 , because it is designed to be perceptually uniform and a color difference measure is defined which outputs JNDs (Just Noticeable Differences). Given a user specified JND, we take the radiosity values from the radiosity solution as a basis to determine which differences in illumination are acceptable.

Our first attempt was to decide according to the gradient of the illumination at the vertex. Although high gradients are a valuable hint for non-smooth illumination, the absolute difference turned out to be more significant. Additionally, it is much harder to find an intuitive threshold for gradients.

Variations in radiosity can be caused either by the form factor or the visibility term, in the diffuse case. For example, a light source could cause strong variations of the form factor, but if no occluders are present, only the form factor has to be recomputed, while visibility can be interpolated. Therefore, we estimate form factor and visibility variation separately and decide individually if form factor or visibility can be interpolated. However, the two functions depend on each other: if a patch is illuminated only weakly, shadows will be less obvious. On the other hand, a strong variation of illumination may be canceled out by a shadow. Therefore, if the variation of one term is measured the other term is kept constant.

We have developed two different oracles. The first, which compares illumination values of neighboring vertices, is more time consuming but turned out to be more reliable than the faster second, which re-uses link information from the radiosity solution.

Estimating Variation from Vertex Differences. The radiosity at a vertex can be computed precisely using the sample distribution which will be described in Section 2.4. The radiosity of each sender is compared with the corresponding values of the neighboring vertices ¡ £¢ to find out how fast the illumination changes within the neigh- borhood of vertex .

We determine the variation of the visibility term ¢¡ ¢ ¤£ £ ¦¥ ¨ § by computing the maximum difference between the radiosity at the center vertex due to sender © and the radiosity at all neighbor vertices ¡ ¢ due to the same sender. As explained before, the form factor is kept constant for this comparison.

¡ ¢ £ £ ¥ ¨ § is then ¡ ¢ ¤£ £ ¥ § £ £ ¥ ¨ § ! #" %$ ¢ ¥ ¡ '& (© ¨£ ¥ ¨ § 0) ¡ '& (© ¨£ ¥ ¡ ¢ 1 § 1 § 32 54 ( 1
)
where £ is the sender's radiosity, 6£ the form factor between and © , ¡ !& 7© %£ ¥ ¨ § is the percentage of © which is visible from 8 , and 2 is the reflectivity of the receiving patch. The variation of the form factor 9 £ @¥ ¨ § is computed analogously:

9 £ @¥ ¨ § 6 £ A #" %$ ¢ ¥ B£ ¥ ¨ § 0) C £ ¥ ¡ ¢ 1 § 1 § ¡ '& (© ¨£ ¥ § D2 FE (2)
For both, G¡ ¢ ¤£ £ H¥ ¨ § and ¢9 £ I¥ ¨ § , a perceptual error threshold can be used. The measure can only determine if there are strong changes expected around the vertex or not, but not if the shape of the function is suited for interpolation. Consequently the measure is rather strict: it will select those contributions which can safely be interpolated, but does not detect all contributions which would also be suitable for interpolation. Of course, it is possible that features inside a triangle are missed: we can only capture what is possible with the given radiosity mesh.

Estimating Variation from Links. The links of the radiosity solution give us a cheaper estimate for the change in illumination across a receiving patch. Again, the maximum of the variation of all patches surrounding a vertex is used as estimation of the variation at a vertex.

For the hierarchical refinement an estimation of the variation of the form factor is needed anyway. For example, we use Bounded Radiosity [START_REF] Stamminger | Bounded radiosity -illumination on general surfaces and clusters[END_REF], which gives a conservative lower and upper bound for the form factor. Thus, the change in radiosity due to the form factor 9 £ @¥ 8 P § over patch 8 due to sender © is:

9 £ I¥ 8 P § B Q £ ¥ upper £ ¥ 8 R § 0) C lower £ ¥ 8 P § S § ¡ !& (© £ ¥ 8 P § 32 5E (3)
The vertex variation can then be computed from the maximum variation of the surrounding patches, e.g. 9 £ ¥ ¨ § B Q #" %$ 'T ¥ 9 £ ¥ 8 R § 1 § .

If partial visibility between a sender and a receiver is detected in a link, the maximum possible change in radiosity goes from zero up to the unoccluded illumination. This gives a conservative bound on the variation due to visibility.

Additionally, we estimate the size U of the penumbra by taking into account geo- metric considerations (see next paragraph). With this penumbra size the bound can be tightened further: if U is larger than the patch the full cut-off will not take place in- side the patch. If the visibility produces a function with gradient V£ W £ ¥ 8 R § D2 YX %U then the change inside the patch with size © ¨`is only £ £ ¥ 8 P § 32 a© bbX cU . In summary the visibility variation G¡ ¢ £ £ ¥ 8 P § is determined by:

¡ ¢ ¤£ £ ¥ 8 R § d £ £ ¥ 8 P § 32 a© b`X %U 04 eU gf h© b £ £ ¥ 8 P § 32 54 else (4)
As before, the vertex visibility variation ¢¡ ¢ ¤£ £ H¥ pi § is then the maximum variation of the surrounding patches. Using the maximum of all surrounding patches decreases the risk of missing an occluder. Estimating the Size of the Penumbra. The size of the penumbra which might be produced by a sender is estimated using simplified geometric configurations. This was done before in a similar way by [START_REF] Soler | Fast calculation of soft shadow textures using convolution[END_REF].

The size U of the penumbra can be approximated by (see Fig. 3)

: U G ¡ £ ¢ ¤£ ¦¥ ¥ ¨ § £ § ¢ ¤£ ¦¥ ¥ ¨ § ` § © © © ) © © © © ) © © 4 ( 5
)
where £ is the maximum diameter of the sender, ¢ ¤£ ¥ ¥ ¨ § £ b § and ¢ ¤£ ¥ ¥ § ` § are the angles between the direc- tion between the sender's and the receiver's centers and the corresponding surface normals.

During the visibility computation of the radiosity algorithm we store the occluder with hit point closest to the receiver. For this closest hit point we compute the relation between distance of the sender to the blocker © © © I) © © and of the blocker to the receiver © © ) © ©

.

Addition of Sender Variations. So far, only single sender contributions have been investigated. But the contributions of different senders may influence each other and in the worst case several small variations may sum up and thus cause a significant variation.

To capture this worst case we select the candidates for interpolation in the following way: The variations of all senders of a vertex are sorted in decreasing order. Then, from the total sum of all variations, we subtract the single variations until the sum drops below the error threshold. The remaining contributions can then safely be interpolated, because even if they should all sum up, they will be below the error threshold. Since the highest variations are subtracted first, we ensure that the most important ones are not missed. On the other hand, this method may be too strict, because in many cases the contributions do not sum up and might even cancel out.

Number of Samples at a Vertex

If the classification procedure determines that it is better to re-compute a contribution exactly, the next step is to find the number of samples which should be used for each sender. This is not exactly true for classification by vertex comparison: the number of samples is determined first to compute the illumination at a vertex exactly in this case.

To determine the irradiance ¥ pi § at a point i the integral over all sending surfaces has to be evaluated: ¥ pi § 6

Senders £ ! £ ¥ ¨" § 3 ¥ Di 4 " § ¡ '& (© ¥ pi 4 " § # " 4 (6)
where ¥ ¨" § is the radiosity of the sender, ¥ pi 4 " § the point-to-point form factor, and ¡ '& (© ¥ pi 4 " § is true if point " is visible from point i .

We want to evaluate this integral for each vertex within a user given precision. As we take the radiosity solution, this means the final gather recomputes the final reflection. The form factor and the visibility term are re-evaluated by randomly choosing $ £ sample points on each sender.

In our algorithm, $ #£ is chosen in accordance to the magnitude of the irradiance at i due to this sender. From the radiosity solution we already have an estimate for this value-not for each vertex but for each patch. The relative contribution of a sender © to a receiving patch 8 compared to the total irradiance at this receiver is £ ¥ 8 R § 7X ¡ £¢ £ ¥ 8 R § , where the irradiance £ ¥ 8 P § can be taken from the links, if the form factor is stored. We use the unoccluded irradiance, because the visibility term would reduce the number of samples in the presence of occluders, although this situation needs a high number of samples to capture the shadow boundary.

Given a total sample number $ to approximate integral [START_REF] Gortler | Wavelet radiosity[END_REF] the number of samples $ £ ¥ 8 P § for sender © and patch 8 is given by $ £ ¥ 8 P § B ¡$ £ ¥ 8 P § (X ¡ £¢ £ ¥ 8 P § . To finally de- termine the number of samples $ £ for a vertex i the average of the number of samples proposed by each neighboring patch sharing the vertex is taken. Only if all senders are classified as "to be re-sampled" the full sample number $ will be used. In general, our algorithm manages with only 30% of $ (see Sect. 3).

One could argue that it would be more efficient to use a total number of samples $ which is not fixed but proportional to the aggregate irradiance, this generally means fewer samples for darker areas. In our case, we base our decision on interpolating contributions on the values computed for the vertices (see Sect. 2.3); therefore having a fixed precision for the vertices is more advantageous. Furthermore, it may be inaccurate to lower the precision for dark areas, because differences in illumination are perceived more easily there.

In many cases the sending patch will be a triangle or another primitive for which an analytic patch-to-point-form factor is known. Although this formula is also expensive to compute-it contains three " ¥¤ ¢ ¢ £ ¦¥ evaluations for a triangle for example-it might be faster than sampling if a high number of samples is used. However, our measurements (see Sect. 3) show that the time spent for form factor computation is low compared to the time spent for visibility computation.

Rendering a Triangle

After the preprocess each vertex has a list of senders with either precomputed illumination for interpolation or a set of sample positions for exact computation.

From the number of samples for the vertices the sample numbers for points inside a triangle can be determined. This is again done per sender: For a point i the number of samples for a sender © is the sum of the three vertex sample numbers for © weighted by the barycentric coordinates of point i .

Note, that we use real numbers as sample numbers to achieve smooth transitions.

Having ¦ ¨ § © samples, actually ¦ § $ samples will be evaluated. Each sample is weighted by one, except for the last which gets weight ¦ ¡ ) ! . Although this results in overhead, it ensures that the sample number is continuous inside the triangle enabling a smooth representation of illumination. Furthermore, this methods works well only if the sample positions are fixed for a sender, i.e., if the & -th sample is always located at position " ¢ on the sender.

It is possible that a sender is classified differently by two neighboring vertices. In this case, we compute the result according to the first as well as according to the second classification and interpolate the two results, weighted by the barycentric coordinates, to obtain a continuous representation.

Results

Test Series

In order to test the behavior of our method, we examined the influence of different parameters on execution time and visual quality of the result. Our test scene is shown in Fig. 4(a). Three objects on a table are lit by two area light sources, a very large one and a smaller one. Both create large variations in irradiance with respect to form factor and visibility. We used the large light source as a representative for indirect light, but by making it a light source we can generate a reference solution by ray tracing.

The large penumbrae result in big problems for standard ray tracers and on final gather approaches. Obtaining an image without visible noise requires about 1000 light source samples per pixel; our reference solution was computed with Monte-Carlo ray tracing with 2000 samples in 90 minutes. The second column shows a less accurate solution with only 56 exact samples per pixel for visibility and 23 for the form factor on average. Artifacts become visible, but the time also decreased to 289s with an overhead of 20%. Sample numbers were determined using link information instead of the more accurate criterion (see Sect. 2.3).

For the right column, the same parameters as for the first image were used, but a coarser radiosity solution (see Fig. 4(c)) was taken as input. The resulting image has almost the same quality as the first one, but rendering time went up to 985s. The number of exact samples for visibility and form factor went up to 265 and 234, respectively. The overhead went down to 7%.

Performance for More Complex Scenes

To test the performance and behavior of the final gather for more complex environments we used three scenes with different lighting characteristics. In the first scene "office" (taken from the Radiance homepage) indirect illumination is present but not very obvious (color bleeding on the wall from the table, see color page, first row). In contrast, the 578s 289s 985s second scene "work and rest" (w.a.r.) is mostly lit indirectly including shadows from indirect light (see color page, second row). Note that, although we selected 400 (importance weighted) samples and the algorithm recomputed the highest possible sample number in these regions, still small artifacts in the indirect shadows of the plants are visible. The last scene, "public library" has raised complexity (see color page, last row). It contains 83 mostly big and lengthy area light sources, casting all kinds of sharp and blurry shadows on the floor, which causes significant problems on any Monte-Carlo and final gather approach. The HR solutions contained 7,800, 53,252, and 97,986 triangles, for office, w.a.r, and publib, respectively. Table 1 shows the sampling statistics. The average proportion of recomputed samples for form factor and visibility varies from 25% for the office to 36% for the public library. The (v) and (l) after the scene denote whether the vertex comparison or the link information was used to find senders suitable for interpolation, respectively. Table 2 lists the timings for the three scenes, split up into the single computation steps. The time to classify for each sender if its contribution should be interpolated or re-sampled, depends strongly on the chosen method. As expected, reusing the link information is much faster than comparing the vertices. This is mainly because for the vertex comparison the corresponding senders for all neighbor vertices have to be found. Both methods turned out to be reliable in finding shadows and other lighting details, even for quite coarse meshes. In some cases the link information was slightly less reliable. The reliability of the link information depends also on the number of test rays used for HR. Furthermore, the bounds computed from the links on the form factor and visibility are not very tight. Therefore generally more samples are used (see Table 1), which might cancel out the time savings during the classification process.

The table column 'vertex' contains all other preprocessing steps: collection of links, determination of sample numbers, and computation of vertex irradiance. The timings here show that for highly refined radiosity solutions the vertex preprocess time slows down the algorithm.

Finally, the time for rendering the interior points was split into time for visibility and form factor computation and the remaining time which is spent for collecting the senders and interpolation. In general, by far most time is spent for visibility computation (75-80%), only a very small amount for form factors ( The images shown on the color plate were created in combination with the Directional Coherence Maps (DCM) [START_REF] Guo | Progressive radiance evaluation using directional coherence maps[END_REF]. In our case the DCMs were used to obtain antialiased images (36 fold oversampling) with only a moderate increase in computation time. The quality of the lighting reconstruction was not affected. The timings in Table 2 are without DCM, with DCM the whole computation took 508s, 1489s, and 4435s for office, w.a.r, and publib, respectively.

Conclusion and Future Work

We presented an improved, adaptive final gather step. For all pairs of receivers and senders an oracle decides whether the sender's contribution can be interpolated or needs resampling. This classification is done separately for the form factor and visibility term. The required information can be taken from the radiosity solution or it can be obtained more precisely from a full final gather at the vertices. If resampling is necessary, the number of samples is chosen adaptively. By interpolating the sample pattern, continuous solutions are obtained.

The accuracy of the final gather step is chosen independently of the radiosity solution. More sophisticated error measures can be used than in standard radiosity, for example perception based criteria, that require a more 'global view' than that of a radiosity refiner. Integration of perception is planned to be enhanced in the future, for example reducing accuracy where textures can mask out artifacts similar to [START_REF] Ramasubramanian | A perceptually based physical error metric for realistic image synthesis[END_REF] or for patches creating frequencies in a range where sensitivity for contrast goes down.

The performance of the final gather scales with the accuracy of its input: for accurate radiosity solutions it is usually faster than for a coarse radiosity solution. This is true as long as single patches of the radiosity solution are larger than a few image pixels. If this is not the case, the overhead of the sample classification becomes dominant. Here, a faster classification procedure which compares total illumination and not the illumination of each sender would be advantageous. For objects containing triangles in sub-pixel space it is probably necessary to switch from the polygonal mesh to a differ-ent representation (e.g. use face clusters [START_REF] Willmott | Face cluster radiosity[END_REF] or a 3D grid bounding complex objects [START_REF] Dischler | Radiosity including complex surfaces and geometric textures using solid irradiance and virtual surfaces[END_REF].) For very complex scenes the push step for links also might become a bottleneck both in terms of time and memory. Therefore, we would like to investigate if the final gather can work directly on the hierarchy of objects and links as created by the HR and not only at leaf level as was done in [START_REF] Mostefaoui | Rendering inhomogeneous surfaces with radiosity[END_REF].

Fig. 1 .

 1 Fig. 1. Left: a candle stick casting a shadow from a large and a small area light source. Center: radiosity solution. Right: visibility sample numbers considered appropriate by our method.

Fig. 2 .

 2 Fig. 2. Overview. Left: Collecting links from neighboring patches. Middle: Determination of sample numbers for each vertex. Right: Computation of illumination for a point inside a triangle.

Fig. 3 .

 3 Fig. 3. Size of the penumbra.

3

 3

Fig. 4 .

 4 Fig. 4. Test scene and radiosity solutions used for final gather Fig. 5 shows three final gather results with sample number 800. The first one was obtained from the radiosity solution in Fig. 4(b) with small final gather thresholds in 578s. The result (top row) is essentially indistinguishable from the reference solution. Only the difference image (bottom row, scaled by 10) reveals errors in particular in the penumbrae. On average, visibility was resampled by 118 samples per pixel, illumination by 65 (out of 800 possible). 81% of the time were spent on the exact sampling, the remaining 19% are used for sample classification etc. (algorithm overhead).The second column shows a less accurate solution with only 56 exact samples per pixel for visibility and 23 for the form factor on average. Artifacts become visible, but the time also decreased to 289s with an overhead of 20%. Sample numbers were determined using link information instead of the more accurate criterion (see Sect. 2.3).For the right column, the same parameters as for the first image were used, but a coarser radiosity solution (see Fig.4(c)) was taken as input. The resulting image has almost the same quality as the first one, but rendering time went up to 985s. The number of exact samples for visibility and form factor went up to 265 and 234, respectively. The overhead went down to 7%.

Fig. 5 .

 5 Fig. 5. Final gather results. Left: fine radiosity solution with accurate final gather. Middle: fine radiosity solution with less accurate final gather. Right: coarse radiosity solution with accurate final gather. Lower row: difference to reference solution.

Fig. 6 .

 6 Fig. 6. Results. Upper small image: number of samples (red: visib., green: formfactor). Lower small image: radiosity solution. Images sizes (top to bottom): 424x363, 512x384, 427x363.

Table 1 .

 1 Sample numbers for test scenes. Total: max. number of importance weighted samples per pixel, avg. vis: recomputed visibility samples, avg. ff: recomputed form factor samples.

	scene	total avg. vis avg. ff
	office(v) 400	101	148
	office(l)	400	130	230
	w.a.r.(l)	400	126	122
	publib(l) 1000	361	230

Table 2 .

 2 Timings in seconds for test scenes (see color plates)

	1-10%) and the overhead is

The term variation will be used in the following not in its mathematical meaning, but to describe magnitude of change in radiosity.

[START_REF] Per | Clustering for glossy global illumination[END_REF] Prior to the conversion to CIELAB a tone mapping should be applied to the color, which is not done in this implementation yet.

All tests in this section were perfomed on a Linux Pentium PC with 733MHz, those in the next sections on a 1.1GHz Pentium

Acknowledgments

The first author is supported by the ESPRIT Open LTR Project 35772 "SIMULGEN", Simulation of Light for General Environments. We would also like to thank George Drettakis for fruitful discussions and proofreading and the reviewers for their valuable comments. iMAGIS is a joint project of CNRS/INRIA/UJF/INPG.