Florent Duguet

George Drettakis

Robust Epsilon Visibility

Keywords: CR Categories: I.3.7 [Computer Graphics]: Three-Dimensional Graphics and Realism-Color, shading, shadowing, and texture I.3.7 [Computer Graphics]: Three-Dimensional Graphics and Realism-Visible line/surface algorithms Illumination, Shadow Algorithms, 3D Visibility, Robust Visibility Predicates, Epsilon Visibility

Figure 1: (Far left) A Dandelion model of 62,000 polygons. Using robust ε visibility, the output shadow mesh has 270,000 polygons for ε 10 4 . From left to right: shadows on the leaves and ground (top view) for ε 10 4 , and ε 10 2 , resulting in 220,000 polygons, with no visible difference. Rightmost two images: zoom of the ground for ε 10 4 and ε 10 2 ; some missing features are now visible.

Introduction

Computing visibility is central for many aspects of computer graphics. From the most basic visible surface determination, to shadows and global illumination calculations, a major part of computational effort is spent on visibility. Although discrete approximations or sampling such as the z-buffer or ray-casting are often used in practice, analytic methods in which a continuous and precise solution is computed, can be very useful. In particular, for shadow calculations [START_REF] Lischinski | Discontinuity meshing for accurate radiosity[END_REF][START_REF] Stewart | Fast computation of shadow boundaries using spatial coherence and backprojections[END_REF][START_REF] Drettakis | A fast shadow algorithm for area light sources using back projection[END_REF], or global illumination [START_REF] Durand | Fast and accurate hierarchical radiosity using global visibility[END_REF], the importance and utility of analytic visibility methods has been demonstrated. For virtual environments or games, graphics engines can now handle large polygon counts, even for low-end platforms. In certain cases, for games engines or virtual reality, texture memory may be a scarce resource, precluding its use for shadow representation. As a result, subdividing the input geometry into (partially) shadowed and lit sub-polygons using one of these methods could well be the best solution for displaying high-quality shadows for these applications. Note that such subdivision is a view-independent solution with shadows, with no additional shadow processing per frame. We will refer to this subdivision of the input scene as a shadow mesh in what follows. The shadow mesh can be created either for hard shadows (such as from a directional or point source) or for soft shadows from an area source.

In practice however, analytic methods have not been used because they suffer from robustness problems, algorithmic complexity and/or memory restrictions. These problems render them unusable for the type of scene used in virtual reality, games or other interactive applications. Such scenes are geometrically complex, and typically contain a large number of degeneracies: objects touch, edges are often aligned or coplanar etc. In addition, since floating point numbers are used both for the modeling phase and visibility computation, below a certain threshold, small features of the model or of the resulting shadows can lead geometric algorithms to fail. Object connectivity is not always given in such models and cannot always be reconstructed, and intersecting polygons are often present in the models. All of these properties can lead to robustness problems for geometric algorithms.

Motivation

To create shadow meshes, in particular for soft shadows, discontinuity meshing approaches [START_REF] Heckbert | Discontinuity meshing for radiosity[END_REF][START_REF] Drettakis | A fast shadow algorithm for area light sources using back projection[END_REF][START_REF] Lischinski | Discontinuity meshing for accurate radiosity[END_REF] intersect shadow boundary surfaces, or swaths, with scene geometry. This approach is inevitably unrobust for large and degenerate scenes, since numerical problems quickly lead to loss of connectivity in the shadow mesh, or result in geometrical errors due to small features.

We can observe however that shadow boundaries are delimited by extremal stabbing lines [START_REF] Teller | Computing the antipenumbra of an area light source[END_REF]. Two examples are shown in red in Fig. 2 (a). One way to generate shadow boundaries is to compute extremal stabbing lines, and then join them to form the actual swaths which constitute shadow boundaries. A simpler problem is thus solved, since in essence we only perform line-or ray-casting. This results in a stabler approach, as was shown in the Visibility Skeleton [START_REF] Durand | The Visibility Skeleton: A Powerful and Efficient Multi-Purpose Global Visibility Tool[END_REF]. For non-degenerate scenes, the number of possible configurations of connectivity between extremal stabbing lines and swaths is finite and small. As a result, Durand et al., [START_REF] Durand | The Visibility Skeleton: A Powerful and Efficient Multi-Purpose Global Visibility Tool[END_REF], developed a catalog of swaths and their neighboring lines to establish connectivity, based on a small set of non-degenerate configurations of edges, vertices and faces. The result of this construction, among other applications, is the shadow mesh. Unfortunately, this approach will fail for complex scenes, such as those encountered in games or VR applications. Consider Fig. 2(b), where vertices of several leaves of a tree are collinear. Building shadow adjacencies based on a finite set of configurations is no longer possible, since a catalog capable of treating all cases would be infinitely big. Despite treatment of certain cases for the test scenes used in [START_REF] Durand | Fast and accurate hierarchical radiosity using global visibility[END_REF], a complete solution is lacking, as mentioned by the authors ( [START_REF] Durand | Fast and accurate hierarchical radiosity using global visibility[END_REF], pages 164 -166). This fact is also clearly demonstrated by the small size and specific (i.e., custombuilt) nature of the test scenes used in all previous analytic approaches [START_REF] Drettakis | A fast shadow algorithm for area light sources using back projection[END_REF][START_REF] Stewart | Fast computation of shadow boundaries using spatial coherence and backprojections[END_REF][START_REF] Lischinski | Discontinuity meshing for accurate radiosity[END_REF][START_REF] Durand | The Visibility Skeleton: A Powerful and Efficient Multi-Purpose Global Visibility Tool[END_REF][START_REF] Durand | Fast and accurate hierarchical radiosity using global visibility[END_REF], in which the largest test scenes used were under 2,000 polygons. If edges or vertices are almost aligned or coplanar (see Fig. 2(c)) problems will appear, since resulting small shadow mesh features become sources of numerical instability, such as small, sliver triangles, edge side-test problems etc. Finally, many models are inaccurate, and may contain intersecting polygons or lack connectivity information. To use the catalog approach, expensive and numerically unstable preprocessing would have to be performed to re-facetize the input scene based on these intersections, and then to reconstruct all connectivity information around edges.

Contributions

To robustly compute shadow boundaries for real-world scenes, we develop an approach based on generalized extremal stabbing lines. Candidate extremal stabbing lines are proposed based on native generators for non-degenerate configurations [START_REF] Durand | The Visibility Skeleton: A Powerful and Efficient Multi-Purpose Global Visibility Tool[END_REF], i.e., VV , V EE and E4, with V a vertex and E an edge. We attach additional generators to the candidate line, using ε methods to robustly determine whether a single feature (edge or vertex) is a generator or a blocker. The set of attached generators ensures that degenerate stabbing lines are correctly constructed, and allows us to robustly establish shadow boundary connectivity. We show that our ε method is consistent, by avoiding undesired propagation of ε contact and by ensuring consistent shadow-boundary connectivity.

The use of ε methods requires special attention when encountering faces (polygons), smaller than the ε threshold. We introduce robust blocker predicates both for the case when face connectivity is available, and for models lacking this information. We also treat models with intersecting polygons as input.

The resulting structure is sufficient for the generation of a shadow mesh. An interesting advantage of our methods is that when the value of ε used in the robust computation is increased, small shadow features will disappear, reducing the size of the shadow mesh.

These methods have allowed us to implement robust shadow mesh computation for soft and directional shadow boundaries for complex, real-world scenes, such as those shown in Fig. 1, 14 -19.

Previous Work

Very early work in computer graphics attempted to compute shadows using polygon clipping methods [START_REF] Weiler | Hidden surface removal using polygon area sorting[END_REF]. These were the precursors of analytic visibility methods which have inspired our work.

It is beyond the scope of this paper to survey all work on visibility. An excellent survey of the whole spectrum of visibility methods in computer graphics can be found in F. Durand's Ph.D. thesis [START_REF] Durand | 3D Visibility: analytical study and applications[END_REF]. For directional and point light sources, shadow volumes [START_REF] Crow | Shadow algorithms for computer graphics[END_REF] allow the computation of shadow boundaries. They suffer however from robustness problems and computational expense, both in the preprocess and for rendering since they are very polygon-fill intensive [17]. Various other methods have been developed, for example the methods of Campbell and Fussell [START_REF] Campbell | Adaptive mesh generation for global diffuse illumination[END_REF] or Chin and Feiner [START_REF] Chin | Fast object-precision shadow generation for areal light sources using BSP trees[END_REF] based on BSP-trees, which also run into robustness or memory problems due to the geometric operations and data structures involved.

Discontinuity meshing approaches were introduced by Heckbert [13] and Lischinski et al. [START_REF] Lischinski | Discontinuity meshing for accurate radiosity[END_REF]. These methods compute shadow discontinuity surfaces (line swaths), by intersecting them with the scene, with the resulting robustness problems previously discussed in Sect. 1.1. The ideas developed in the original aspect graph vision literature [START_REF] Koenderink | The internal representation of solid shape with respect to vision[END_REF][START_REF] Gigus | Computing the aspect graph for the line drawings of polyhedral objects[END_REF], and introduced to graphics by Teller [START_REF] Teller | Computing the antipenumbra of an area light source[END_REF] have led to algorithms computing the complete discontinuity mesh [START_REF] Drettakis | A fast shadow algorithm for area light sources using back projection[END_REF][START_REF] Stewart | Fast computation of shadow boundaries using spatial coherence and backprojections[END_REF]. These include the backprojection data structure representing the visible part of the source for rapid computation of penumbral lighting. These methods suffer from the same robustness problems.

We have already mentioned the Visibility Skeleton [START_REF] Durand | The Visibility Skeleton: A Powerful and Efficient Multi-Purpose Global Visibility Tool[END_REF], which computes extremal stabbing lines by casting them into the scene, and uses a finite catalog to establish swath connectivity. This approach was used to compute shadows, but also in a global illumination algorithm which adapts the radiosity mesh using perceptual criteria to insert discontinuities [START_REF] Durand | Fast and accurate hierarchical radiosity using global visibility[END_REF].

Many approximate or discrete methods have been developed for shadow computation (e.g., shadow maps [START_REF] Williams | Casting curved shadows on curved surfaces[END_REF]) and ray-tracing is evidently one of the most widespread point-sampling approaches used in practice. For soft shadows, convolution methods [START_REF] Soler | Fast calculation of soft shadow textures using convolution[END_REF] give convincing results, which may contain some error however, depending on the configuration. All these methods use a discrete buffer, and thus resolution and resulting aliasing are the predominant problems. More recent work [START_REF] Agrawala | Efficient image-based methods for rendering soft shadows[END_REF] presents an image-based approach comprising an efficient interactive solution with limited control of sampling and an expensive, view-dependent ray-tracing solution.

We adopt an ε geometry approach [START_REF] Salesin | Epsilon geometry: Building robust algorithms from imprecise computations[END_REF], albeit without perturbing the input data, using interval arithmetic to robustly treat degenera-cies for our visibility computations. Interval approaches have been extensively used in computer graphics, for example for implicit functions and CSG [START_REF] Duff | Interval arithmetic and recursive subdivision for implicit functions and constructive solid geometry[END_REF] or by Snyder [START_REF] Snyder | Interval analysis for computer graphics[END_REF] for a variety of problems such as ray-tracing, interference detection etc. To our knowledge, ε methods have only been used to a limited extent for 3D visibility or shadow problems, notably by Bala et al. [START_REF] Bala | Radiance interpolants for accelerated boundederror ray tracing[END_REF] for ray-tracing.

Generalized Visibility Events

As mentioned in the introduction, we will treat scenes containing degeneracies, such as those used in real-world applications. We first define our new framework, which is based on a general definition of extremal stabbing lines (ESL) and line swaths. This definition overcomes the limitations of a catalog of visibility events [START_REF] Durand | The Visibility Skeleton: A Powerful and Efficient Multi-Purpose Global Visibility Tool[END_REF], since it will not depend on a finite enumeration of specific configurations. We first define the basic entities we use.

A generator is an edge or a vertex of the scene 1 . For edges we distinguish silhouette edges with respect to a given line direction, in the traditional sense of the term, i.e., that the edge is attached to a single face, or that one of the two faces connected to the edge is back-facing with respect to the line. A silhouette vertex is defined depending on its local blocking properties, i.e., whether it blocks an extremal stabbing line or not, depending on its neighboring faces. This is defined in detail in Section 4.1 (see also Fig. 6).

E 2 V 2 V 1 block down block up l 2 l 1 E 1 Figure 3: An extremal stabbing line l 1 , with generators V 1 V 2 E 1 E 2 ,
and up and down blockers block up and block down , and a swath between ESL l 1 and ESL l 2 , sharing E 2 and V 2 .

A blocker is a face, or non-silhouette edge or vertex touched by a line.

A generalized extremal stabbing line (ESL) is a data structure, containing the corresponding line equation (also called supporting line), the set of associated generators, and its up and down blocker. A degenerate ESL, l 1 , its generators and blockers are shown in Fig. 3.

The limits of a critical line swath, or simply swath, are its up and down blockers, and the start and end ESL's (see Fig. 3).

To establish shadow boundaries, we need to connect stabbing lines with swaths. Since we can no longer use a finite catalog, we will define when it is possible to connect pairs of ESL's. The intuition is that there needs to be a sufficient number of shared generators to define a planar (e.g., edge-vertex) or quadric (edgeedge-edge) swath joining the two ESL's. The process we use is procedural, and detailed in Section 4. [START_REF] Campbell | Adaptive mesh generation for global diffuse illumination[END_REF].

In what follows we consider visibility events between an "emitter" and a "receiver" pair. This is a logical choice for shadow computations, and allows us to define distances and sizes in a meaningful manner for our ε algorithms. This can be seen as a lazy evaluation of parts of the visibility skeleton [START_REF] Durand | The Visibility Skeleton: A Powerful and Efficient Multi-Purpose Global Visibility Tool[END_REF]. In theory, the Visibility Skeleton could be computed using our algorithm, if we computed and stored all emitter-receiver pairs. This would however be impractical for any scene of reasonable size, due to the (at least) quadratic complexity in memory. Our ε approach is based on distance, defined across the shaft formed by the emitter and receiver. Thus merging operations vary depending on the chosen receiveremitter pair; querying the entire structure would require special manipulation. However, for all the applications using the visibility skeleton to date [START_REF] Durand | The Visibility Skeleton: A Powerful and Efficient Multi-Purpose Global Visibility Tool[END_REF][START_REF] Durand | Fast and accurate hierarchical radiosity using global visibility[END_REF], (shadows, discontinuity mesh, global illumination), emitter-receiver pair computations suffice, so we do not consider this restriction significant.

Enumeration and Validation of ESL's and Swaths

The first step in our approach is the enumeration and validation of extremal stabbing lines (ESL). We construct ESL's, associating all appropriate generators as we proceed, as well as start and end blockers, if they exist. The second step is the identification and validation of swaths connecting ESL pairs. ESL's are "cast" into the scene to identify potential generators. We will use ε methods to consistently treat "almost" coplanar or "almost" collinear events, during this process. Before the enumeration and validation, we present the specifics of our ε methods.

ε methods

While casting lines, we need to identify whether a scene element, face (polygon), edge or vertex interacts with the line. If there is an interaction, it can be either a blocking interaction, i.e., we consider that the line is blocked, or can be a generator interaction, i.e., the line grazes the element, for example an edge, and thus we consider that the edge is a generator of the line. Consistent attribution of generators is a key to correct connectivity of shadow boundaries, since ESL's are connected by looping over generators.

We perform all calculations to determine whether an interaction exists, and its type, in floating point arithmetic, up to a precision of a pre-determined ε threshold.

2e 2e The ε parameter is given in the same units as the model of the scene. Thus if the scene is defined for example in meters, ε will also be defined in meters. Thus, the algorithm reacts as expected to scaling operations, i.e., if the entire scene is scaled by a coefficient of say 2, then the result will be the same in the topological sense, as if ε was divided by 2.

E (a) v 1 v 2 (b) (c)
We use an interval arithmetic approach with floating point values for upper and lower bounds. We define ε contact in a natural manner: two elements are in ε contact, if and only if their distance is below ε. We say that an edge or vertex of the scene is in interaction with a line if and only if its distance to the line is below ε. For a vertex p and line we project the vertex onto the line to find the closest point q, and for a line and an edge, p and q are the pair of closest points of the two lines. An equivalent definition is that a vertex or an edge interacts with the line if it hits the associated fat element, i.e., sphere or shaft. This ε contact approach is central to our algorithm.

As a result of our ε approach, scene elements are transformed. A line is transformed into a shaft, which we call a fat line or fat ray (Fig. 4(a)). This shaft is the cylinder defined by the spheres of any two points centered on the line.

A vertex is a sphere of radius ε, which we will call a fat vertex, Fig. 4(b), and an edge is a cylindrical "shaft" (similar in spirit to that of [START_REF] Haines | Shaft culling for efficient ray-traced radiosity[END_REF]) between the two spheres of its extremal vertices, which we will call a fat edge, Fig. 4(c).

Since faces are only considered as blockers in the context of our approach, they remain as is.

There are five kinds of interactions between a line and a face, shown in Fig. 5. The first two cases, no-stab and full-stab, are handled trivially (Fig. 5(a) and (b)), since the first case results in no interaction, and the second results in the line being blocked. The third and fourth (Fig. 5(c) and (d)) are followed by a blocker/generator test described below, in order to check if the element is a blocker or a generator for the line. The last case is more complex and is the essence of our ε approach (Fig. 5(e)).

In Section 5, we introduce techniques to robustly determine whether this configuration blocks the line.

The blocker/generator tests for edges and vertices proceed as follows. Since the line may not pass exactly through the edge, for an edge-stab we form a plane Π with the edge (points A and B in Fig. 6(a)) and normal to a direction perpendicular to the line with its origin on the edge. We can always define this plane since we do not have a vertex-stab, and thus the edge is not on the line. We then compute the position of the faces attached to the edge with respect to this plane. If faces are present on both sides, the edge is a blocker, otherwise, the edge is a generator (Fig. 6(b)). This is equivalent to determining whether we have a silhouette edge with respect to the line.

For the vertex stab, let Π be a plane orthogonal to the line, for instance a plane passing through the vertex. We compute the orthogonal projections of the faces surrounding the vertex onto Π. Each face around the vertex can be seen as a slice of a pie. We then merge the slices in contact. If the slices merge into a whole pie, then the vertex is a blocker, otherwise, it is a generator (see Fig. 6(c)).

Enumeration and Validation of ESL's

We have now defined the framework for generalized visibility events and our ε visibility approach, which are the basic elements required to define our new, robust, algorithm. We start with the enumeration of non-degenerate or generic ESL's, that is VV , V EE and E4, as defined in previous work [START_REF] Teller | Computing the antipenumbra of an area light source[END_REF][START_REF] Durand | The Visibility Skeleton: A Powerful and Efficient Multi-Purpose Global Visibility Tool[END_REF]. For VV and V EE the line equation is defined as in [START_REF] Durand | The Visibility Skeleton: A Powerful and Efficient Multi-Purpose Global Visibility Tool[END_REF]. For E4 ESL's we adopt the method previously proposed by Teller [START_REF] Teller | Visibility Computations in Densely Occluded Polyhedral[END_REF]. An alternative would be to use the bisection approach described in [START_REF] Durand | The Visibility Skeleton: A Powerful and Efficient Multi-Purpose Global Visibility Tool[END_REF].

The choice of which generators are considered as native depends on the order in which we evaluate them: We demonstrate below that a consistent result is obtained with our method independently of the order of evaluation.

The algorithm proceeds with an ESL casting process. First, generic ESL's are formed by enumeration of native generators. For each such generic line, we propose a candidate ESL and cast it through the scene. The ESL-casting process has a dual function: we find and attach all additional generators to form the appropriate set attached to a degenerate ESL, and we test visibility of the stabbing line to decide whether the candidate ESL is valid. Finally, redundant ESL's are also eliminated.

ESL Casting

A candidate ESL is defined by a line and a set of native generators (VV , V EE, E4). The ESL casting process is very similar to a traditional ray-casting algorithm, but is based on "fat rays" as described previously. The casting process identifies the elements, either face, edge or vertex, which interact with the ESL, as described in the previous section. The fat ray emanates from the emitter, and we visit elements in order of increasing distance from the origin of the fat ray. This is achieved cheaply by the traversal routine of the acceleration structure described below; a sort is performed on elements in the interior of a cell of the acceleration structure.

For each such element, we test if it is a blocker, using the ε stabbing process. If it is, we flag the element as the down blocker of the extremal stabbing line, and the casting process is stopped. If not, we add the element to the set S of stabbed elements associated with the ESL, which is initially empty. Note that the native generators are not contained in S at the outset of the casting process.

Once the ray is stopped, we compare the set of stabbed elements S, and the set of native generators N. If N S, then all the native generators are touched by the candidate ESL, which is then validated. The set of stabbed elements also constitutes the set of generators associated with the ESL. This process is illustrated in Fig. 7. To accelerate swath validation, references to the ESL's are stored on the scene edges and vertices. Even though the combinatorial criterion is satisfied with respect to l 1 , the topological criterion is not, and thus l 3 is maintained as a separate ESL.

em itt er em itt er f v 1 v 2 S={v 1 , v 2 , e 1 , e 2 } e 1 e 2 v 1 v 2 e 1 e 2 f S={v 1 } v 1 v 2 S={} em itt er (a) (b) (c)
v 2 v 4 v 1 v 3 v 1 v 2 v 4 v 1 v 2 v 4 v 3 v 3 >e 2e l 1 l 2 l 1 l 1 l 3
We use standard acceleration techniques for fat-ray casting. Faces are simply added to the acceleration structure cells if they are either inside the cell or at a distance less than ε. We use a recursive grid for faces and subdivide based on a 3 Ô n criterion, where n is the number of polygons at each level. Most other acceleration techniques for ray tracing could also be used for our ESL casting.

Elimination of Redundant ESL's

As mentioned above, we consider the order in which we traverse generators to be irrelevant. For this to give a consistent result, we need a method which will eliminate redundant ESL's when they are "sufficiently close" to an existing ESL, with respect to the ε chosen. This elimination step must be done carefully, to avoid propagation of ε across generators, which would result in the "merge" of lines which in reality are geometrically far apart.

We thus define a consistent elimination method for ESL's, based on two criteria: a combinatorial and a topological criterion.

An ESL l c satisfies the combinatorial criterion with respect to another ESL l e , if and only if the set of its native generators is completely contained in the set of generators of the other ESL l e .

An ESL l c satisfies the topological criterion with respect to another ESL l e , if and only if the distance between the two supporting lines is less than ε.

The distance between two lines used for this test is defined with respect to the emitter-receiver shaft. In particular it is defined as the maximum distance (over the shaft), between any point of one line with respect to the other line. If an ESL satisfies both criteria with respect to another, existing ESL, then it is redundant and is eliminated. The entire process is illustrated in Fig. 8.

Note that the test can be performed before the actual ESL casting. We thus compare the candidate line to existing ESL's preventing unnecessary computation implied by the actual casting algorithm. We only search a subset of existing ESL's, since we first need to satisfy the first criterion.

Enumeration and Validation of Line Swaths

Once we have generated all the ESL's, we create the swaths which join them, if required by the application. Since we have foregone the catalog, we use a procedural construction. Swath creation involves two steps, enumeration, to determine whether a swath can exist between two ESL's, and validation, which verifies visibility with respect to the rest of the scene.

Swaths are found by looping over all edges and vertices in the scene. We first loop over all pairs of edges EE to treat the special case of coplanar edges. We also collect any other edges which are coplanar to this pair; We discuss details of this case later. At the same time we find any potential EEE swaths. We then loop over all V and E pairs.

For each potential swath, EEE or EV , we examine the extremal stabbing lines shared by each generator, E or V . If there is a shared pair of ESL's, we create a candidate swath, which is then tested for visibility.

We need to validate the candidate swath, that is treat visibility with respect to the other elements of the scene. In the catalog approach [START_REF] Durand | The Visibility Skeleton: A Powerful and Efficient Multi-Purpose Global Visibility Tool[END_REF], this step is not necessary, since it is included in the local connectivity information. An example of a swath which must be discarded is shown in Fig. 9.

Since all ESL's have been computed, there cannot be a visibility discontinuity along the swath, otherwise it would be represented by an ESL. Thus, to validate the swath, it suffices to sample visibility by a line in the middle of the line set, that is the line through the point in the middle of the edge portion concerned by the swath, which we call the midline.

The midline is then cast in the same manner as an ESL into the scene, using the acceleration structures. The swath is valid if the midline is not blocked before reaching the receiver.

A special case occurs for coplanar edges, as was pointed out in [START_REF] Heckbert | Discontinuity meshing for radiosity[END_REF]. This special case corresponds to a line set of dimension 2, as opposed to other swaths (EEE and EV ) which are of dimension 1, since they only have one degree of freedom (along the E in an EV for example). We first collect all the ESL's contained in the plane, and then find the extremal lines, in the sense of the emitterreceiver pair, shown in red in Fig. 10. We then sort the ESL's across the generating emitter edge, and perform the visibility test for each

Swath Connectivity Consistency

To demonstrate consistency, we adopt the approach of the visibility skeleton [START_REF] Durand | The Visibility Skeleton: A Powerful and Efficient Multi-Purpose Global Visibility Tool[END_REF], in which ESL's are the nodes of a graph, connected by the swaths which are arcs.

We want to demonstrate that if we merge a certain number of nodes together, due to our ε approach, the graph remains consistent. By consistent, we mean that all arcs which could potentially exist in the new configuration are preserved.

Consider the example shown in Fig. 11 (left). If we increase ε sufficiently, l 1 and l 2 will merge, Fig. 11 (right).

Once the merge has taken place, the graph becomes as shown in Fig. 11, upper right. In practice, all generators of l 2 are attached to l 1 , assuming that l 1 was created first. We refer to the new "merged" node as l m . We need to show that all arcs which were connected to all the nodes before merging are taken into account. Arcs which connected the two merged nodes no longer exist. We need to show that arcs originating at another node and arriving at one of the nodes merged into l m , will now be connected to l m . This is simple to show, since the generator set of the merged node l m is the union of the generator sets of all constituent nodes (l 1 and l 2 ). Consider any given external node l e (e.g., l 3 ) which was linked to one of the constituent nodes before the merge. The merged node l m thus contains the necessary generators and we are able to create an arc between l e and l m , when looping over generators for swath enumeration.

When creating swaths between merged nodes, an edge is necessarily shared between them. This edge is used for the midline computation, thus guaranteeing a consistent result, since it will reflect the original structure of the shadow boundary. (Right) The ε threshold is higher, and thus l 1 and l 2 will merge. A merged node l m now takes the place of l 1 and l 2 in the graph.

Robust Blocker Predicates for ε Visibility

The use of fat lines, or an ε approach, implies that all geometric computation is performed up to a certain precision. Once the threshold is chosen some elements may be of size smaller than the given ε. For instance, consider the example shown in Figure 12(a) where line l interacts with many small faces. Since the faces are smaller than the ε chosen, they should not be considered individually.

e e l l In particular, we need a robust way of determining whether the line is blocked by multiple faces.

For this, we develop two solutions: the multi-face structure, which assumes that the scenes contain faces with connectivity information; and the blocker-fan, which treats multiple interactions without connectivity information, and handles surface contact.

If the original model does not contain connectivity information we reconstruct it as much as possible, by finding vertices at the same position and determining faces sharing an edge. The blockerfan is used for all remaining cases, due either to numerical imprecision or touching faces.

The Multi-face Structure

If an ESL has multiple interactions with a face (e.g., Fig. 5(e) or Fig. 12(a)), we apply the multi-face approach. The multi-face is a data structure containing the faces interacting with the fat line, and the contour of these faces used for the blocker test. The idea of the multi-face is related to the application of face-clustering to Upper row: The line is not blocked, both due to depth and to slice coverage. Middle row: The line is blocked, since it satisfies both criteria. Lower row: Special case of surface contact. The valid region is defined by the surfaces in contact. visibility presented by Leblanc and Poulin [START_REF] Leblanc | Guaranteed occlusion and visibility in cluster hierarchical radiosity[END_REF], but is developed and used in a very different context.

For each surrounding face which is touched by the line, in the ε sense (i.e., distance less than ε), we add the face to the multi-face. All connected faces touched by the line are added if the edge is not a silhouette with respect to the line. Once the set is complete, we identify the boundary edges of this set.

To determine if the multi-face blocks a given fat line, we check if the projection of the boundary edges onto a plane orthogonal to the line is a complete loop around the line. If it is, then the multiface is a blocker; Otherwise, it is not. In addition, the line must be entirely contained in the loop in the ε sense (see Fig. 12). This test is performed by traversing the neighboring faces to the line, adjacent to the face(s) interacting with the line.

The Blocker-fan structure

For unconnected polygons we use a different approach. We collect a "fan" of independent polygons encountered along the line of sight of the fat ray, as it is propagated through the scene. The elements encountered are either edges or vertices, and they form slices with respect to the shaft, as shown in red in Figure (13, middle column). We also consider an ε depth for each slice, in the direction of the line. The depth is taken around the point of intersection of the edge or the vertex with the line.

As with the multi-face, we project the faces onto a perpendicular plane. For each face, we take the edges outside the fat line, and create angular slices (see Fig. 13). If the slices collected cover a whole section of the shaft, and they overlap in depth, then we have a blocker. Note that, as with fat-ray casting (see Sec. 4.2.1), we visit elements in order of distance on the line.

The blocker fan also handles surface contact. In particular for touching objects and a line passing "between" them, we perform the same slice operation (see Fig. 13, lower row), and verify that the slices overlap in depth. Note that the blocker fan works together with fat ray casting, so that the routines must be called by the casting process.

Implementation and Results

We have implemented a system based on the algorithms presented above. In addition to the methods described in the previous sections, we compute intersections between polygons. Instead of refacetizing the input polygons along these discontinuities, we consider them as special generators and blockers, rather than true polygonal elements. They are thus attached to the original input geometry.

We have also implemented an octree acceleration structure for the hourglass formed by two edges when looping over generators as described in [START_REF] Durand | The Visibility Skeleton: A Powerful and Efficient Multi-Purpose Global Visibility Tool[END_REF].

All the results and tests we present have been run on scenes containing unmodified objects we have either found on the internet or which are used in real applications. All timings are on a Pentium III 1 Ghz PC under Linux, with a GNU C++ implementation. The example scenes we use are a Dandelion (Dand) (Fig. 1), a Cactus (Fig. 15), "Vizzy the skeleton" (Fig. 16), the Agave plant (Fig. 19 right), the "Big" scene (Fig. 14), which is the largest scene tested, and a very degenerate aligned cube scene (Fig. 17). For soft shadows we have also used the Tree of "Big", containing 33,000 input polygons (Fig. 19, middle).

Shadow Computation

To compute shadow boundaries for directional light sources we only need to enumerate ESL's in the direction of the light source, and we only have planar swaths. We use all algorithms developed; that is the general definition of ESL, ε methods, multiface, blocker fan and intersecting polygons. Results for directional light sources are shown in Fig. 141516.

We have also computed a solution for an extremely degenerate model, which is a cube of cubes, in which everything is aligned, shown in Figure 17. This computation took 20 minutes.

To compute soft shadows, we compute the ESL's and the swaths for each such emitter-receiver pair, as described previously. Note that we only compute visibility events with one generator on the source; in the full discontinuity mesh [START_REF] Drettakis | A fast shadow algorithm for area light sources using back projection[END_REF][START_REF] Stewart | Fast computation of shadow boundaries using spatial coherence and backprojections[END_REF] other events also exist, but we choose to ignore them. We made this choice for computational efficiency and because the impact of these events on shadows is usually minimal. The algorithm could easily be adapted to compute them all, at the cost of an expensive additional search for ESL's within the shaft.

To create images with shadows, we insert all discontinuity lines into a Constrained Delaunay Triangulation (CDT), and then sample the light source. Nonetheless, standard constrained Delaunay triangulation packages are numerically unrobust. For soft shadows, the CDT packages we tried failed for very complex models. We show however the discontinuity edges for an area light source computed by our algorithm, without actually computing a CDT and the resulting shadows. We computed the discontinuity edges on all objects for the Vizzy and Tree models. Due to memory limitations of our current implementation, we restricted the computation to the floor only for Agave. These results are shown in Figure 19.

Figure 18: Soft shadow on the rose and on a 69,000 polygon bunny.

For the rose (1,700 polygons) and the bunny (69,000 polygons) models, the CDT was able to produce a soft shadow mesh, Fig. 18.

Statistics

As we can see from Table 1, we can treat models of significant complexity (from 11,000 to 121,000 input polygons), and the resulting shadow meshes multiply the polygon count by 2-5 times for a directional source; this a reasonable overhead for high-quality shadows. The Cactus is an exception where we have a factor of 13 increase. For soft shadows, the shadow mesh for the rose is 49,500 polygons and required 2h, while for the bunny it is about 95,000 polygons and required 3h; the respective times for computing the discontinuity edges only were 8m and 1h46m. Our soft shadow computation could be significantly optimized by using a hierarchical data structure on the scene graph and using more appropriate acceleration structures. For the discontinuity edges only, the times are 1h43m for Vizzy, 17h30m for Tree and 2h40m for Agave (floor only). Other than the unoptimized implementation, the reason for the high computational cost for Tree is also the large number of E4 events (91,000) compared for example to Vizzy (446). It is also interesting to note the number of swaths, 464,000 for Tree compared to 34,000 for Vizzy. It is interesting to note the much higher number of ESL's for the Cactus which has multiple interactions between spikes, compared to the Vizzy model which is quite sparse and does not contain multiple interactions. The percentage of degenerate ESL's varies, but is always more than 10%, and in some cases the number of degenerate ESL's is at the same level as the non-degenerate lines, even for small ε (see Table 2) which underlines the importance of our approach. Similarly, the number of multi-face operations on the various models shows that the systematic treatment we propose is indispensable. The running times do not always decrease with the increase in ε, due to the potential additional complexity of certain steps of our algorithm such as the blocker-fan computation.

Discussion and Conclusion

We have presented a novel, systematic approach to robust shadow boundary computation on real-world models. Our method first redefines visibility events in a generalized framework, thus taking into account typical degenerate configurations, notably collinear and coplanar scene features. We then develop an algorithm to robustly compute generalized extremal stabbing lines using ε methods, and a procedural approach to connect shadow boundaries. We demonstrated that our approach is consistent both in terms of the usage of ε and in preservation of shadow boundary connectivity. Finally, we present robust blocker predicates for ε computations, both for connected models, and models which do not have connectivity information. We also handle intersecting polygons. Our results show that we can robustly compute both directional and soft shadow boundaries on complex models.

Our method greatly improves the usability of analytic shadow algorithms for real scenes. We consider the idea of consistent feature merging to be of particular interest, since it may lead to an interesting development of hierarchical visibility methods. Currently, the multi-face approach requires connectivity information. If connectivity information is missing, interesting problems arise, since, for large ε values we may wish to merge small features which are relatively far apart, such as the bones of the Vizzy model Fig. 16. We find this line of research interesting, and inevitably dependent on application requirements in terms of output.

Evidently, our approach is limited by the inherent computational complexity of visibility events. It is known that the worst-case computational complexity of the full discontinuity mesh is O´n 4 µ for the E4 ESL's, with n the number of objects/edges, multiplied by mesh construction cost which can be O´n 2 µ on each object. In cases where the number of E4 events is large, such as the Tree example, this effect becomes significant, although the worst case bounds are still very pessimistic. In our approach ESL enumeration and elimination are brute-force, but cheap, while the more expensive ESL validation step is output dependent. Developing optimal algorithms is an interesting research direction, and should follow the spirit of sweep approaches (e.g., [START_REF] Stewart | Fast computation of shadow boundaries using spatial coherence and backprojections[END_REF]).

As with any robust algorithm, it is interesting to discuss possible cases of failure. The use of ε predicates means that we will never block an ESL which should not be blocked. For the same reason however, we may enumerate and validate an ESL which should have been blocked. This is not problematic since no swath will be connected to this ESL, and therefore the output will remain consistent. In its current form, we do not guarantee topological consistency of shadows for very large ε values; shadows simply disappear after a certain point. More sophisticated application-dependent approaches could be developed to guide the way simplified swaths are created.

Our approach could be directly applied to the global illumination method of [START_REF] Durand | Fast and accurate hierarchical radiosity using global visibility[END_REF], since all operations are performed in that method on a link basis, which is exactly equivalent to the emitter-receiver pair we use here.

Figure 2 :

 2 Figure 2: (a) Shadow boundaries are delimited by extremal stabbing lines (in red). (b) In complex models, stabbing lines may contain multiple collinear vertices. (c) If vertices are almost collinear, we choose to consider them as being on the same line.

Figure 4 :

 4 Figure 4: (a) A fat line or fat ray, (b) A fat vertex (c) A fat edge.

Figure 5 :

 5 Figure 5: (a) No-stab, the line does not interact with the triangle (b) Full-stab, the line is blocked by the triangle (c) Edge-stab, the line touches the edge but nothing else (d) Vertex-stab the line touches the vertex but nothing else (e) Multi-stab: the fat line touches multiple elements.

Figure 6 :

 6 Figure6: (a) Determining whether an edge is a blocker for an edgestab (b) Top, the edge is blocked, bottom the edge is a generator (c) Top the vertex is a generator, bottom, the vertex is a blocker.

Figure 7 :

 7 Figure 7: The ESL casting process. (a) Before casting, the line is formed by the native generators v 1 and v 2 , S is empty and the line originates at the emitter. (b) The ESL encounters the generator v 1 , and in this configuration, the ESL is invalidated by the face f . (c) Here, the ESL encounters generators v 1 v 2 e 1 e 2 . This is a valid ESL, with down blocker f , since S contains the native generators.

Figure 8 :

 8 Figure 8: Elimination of redundant ESL's. (a) The native generators of l 1 are v 1 and v 2 , and we add v 3 and v 4 . (b) Native generators of l 2are v 2 and v 3 , and we add v 1 . l 1 and l 2 satisfy the combinatorial criterion and the topological criterion, and thus l 2 is eliminated. (c) Native generators of l 3 are v 3 and v 4 , and we add v 2 . Even though the combinatorial criterion is satisfied with respect to l 1 , the topological criterion is not, and thus l 3 is maintained as a separate ESL.

Figure 9 :

 9 Figure 9: Swath validation. (a) The swath ev is valid since the midline does not touch any other object. (b) The swath e 2 v is invalid, since the midline is blocked by the face f , and the receiver has not been encountered.

Figure 10 :

 10 Figure 10: Coplanar edge treatment and shadow mesh segment insertion.

Figure 11 :

 11 Figure11:(Left) Original configuration with a small ε. The corresponding graph is shown in the upper left. (Right) The ε threshold is higher, and thus l 1 and l 2 will merge. A merged node l m now takes the place of l 1 and l 2 in the graph.

Figure 12 :

 12 Figure 12: (a) A fat line covers several faces, and the fat line is blocked since it is contained inside the multi-face contour, in black. (b) The multi-face does not block the fat line, since it is not contained in the multi-face contour.

Figure 13 :

 13 Figure 13: Blocker-fan: (Left) Perspective view of the unconnected faces encountered by the line. (Middle) The slices, on a plane perpendicular to the line. The red lines show the angular slices. (Right) View of the slice along the line direction to verify depth overlap. Upper row: The line is not blocked, both due to depth and to slice coverage. Middle row: The line is blocked, since it satisfies both criteria. Lower row: Special case of surface contact. The valid region is defined by the surfaces in contact.

Figure 14 :

 14 Figure 14: "Big" scene, overview and details of shadows on flower, leaves and ground.

Figure 15 :

 15 Figure 15: Two views of the cactus model.

Figure 16 :

 16 Figure 16: An overview and a zoom of "Vizzy the skeleton" with a directional light.

Figure 17 :

 17 Figure 17: The cube of cubes degenerate model and resulting shadows (directional light source parallel to the cubes).

Figure 19 :

 19 Figure 19: Discontinuity edges for an area light source for three complex models (Vizzy, Tree and Agave -floor only).

Table 1 :

 1 Statistics for test scenes with directional shadows, polys: number of input polygons, time: time to compute shadow mesh,

	Scene/ε	Vizzy/1e 4	Cactus/1e 4	Dand/1e 4	Agave/1e 4	Big1e 4
	polys	11,000	34,000	61,000	114,000	121,000
	time	0m35s	1h10m	2h40m	0h04m	0h10m
	smesh	31,600	428,000	287,000	153,000	329,500
	ESL	4,300	45,000	75,000	20,000	19,000
	dESL	535	19,600	7,700	5,000	4,300
	mface	254	17,000	1,200	7,700	3,100

smesh: size of shadow mesh, ESL: number of generic ESL's, dESL: number of degenerate ESL's, mface: number of ESL's implying multi-faces.

Table 2 :

 2 Additional statistics for test scenes for varying ε.

In contrast to[START_REF] Durand | The Visibility Skeleton: A Powerful and Efficient Multi-Purpose Global Visibility Tool[END_REF], we do not consider faces as generators. In the context of a finite catalog, such a consideration is appropriate, since it allows a more concise presentation of possible configurations.

Acknowledgments

Thanks to P. Poulin for reading an early version and M. Stamminger for his many comments and for making the video. This work was partially supported by the INRIA Action de Recherche Coopérative (ARC Vis3D) on 3D Visibility. The initial ideas leading to this work originated in F. Durand's Ph.D. thesis. We would like to thank him for the many fruitful discussions and exchanges.


Eric Paquette

The Simulation of Paint Cracking and Peeling

Keywords: Deteriorations, surface imperfections, paint, cracks, multi-layer surfaces, natural phenomena

Weathering over long periods of time results in cracking and peeling of layers such as paint. To include these effects in computer graphics images it is necessary to simulate crack propagation, loss of adhesion, and the curling effect of paint peeling. We present a new approach which computes such a simulation on surfaces. Our simulation is inspired by the underlying physical properties. We use paint strength and tensile stress to determine where cracks appear on the surface. Cracks are then propagated through a 2D grid overlaid on the original surface, and we consider elasticity to compute the reduction of paint stress around the cracks. Simulation of the adhesion between the paint and the underlying material finally determines how the paint layer curls as it peels from the surface. The result of this simulation is rendered by generating explicit geometry to represent the peeling curls. We provide user control of the surface properties influencing the propagation of cracks. Results of our simulation and rendering method show that our approach produces convincing images of cracks and peels.

Introduction

Cracks and peels in paint are a common everyday phenomenon. For example, after a few years of exposure to the elements, cracks appear on painted wooden doors, often resulting in peeling of the paint layer away from the underlying wood. These visual effects can be very significant (see for example the left images of Figures 8,9, and 10), and we believe that it is important to include them in computer generated images. As with all aging and weathering effects, simulating and rendering such phenomena results in images which have a much higher degree of realism. In typical synthetic environments, used for example in films or video games, many painted objects exist as well as objects which are made from layered surfaces, and thus it is important that they include the simulation of cracking and peeling.

Previous approaches in computer graphics have either only simulated cracks and their propagation without treating multi-layer phenomena and peeling [START_REF] Gobron | Crack pattern simulation based on 3D surface cellular automata[END_REF][START_REF] Hirota | Generation of crack patterns with a physical model[END_REF], or presented tools to control the location of detached paint areas without treating the creation and propagation of cracks [START_REF] Wong | A geometry dependent texture generation framework for simulating surface imperfections[END_REF]. Consequently, no complete method currently exists, which allows a user to simply request the generation of cracks and peels on a surface.

Our approach is inspired by the underlying physical phenomena, but uses a simplified model. It allows a fast and simple simulation process, provides good control to obtain the desired effects, and gives convincing results.

Our current implementation operates on a planar surface, and assumes that the outer paint layer is infinitesimally thin. We generate cracks on the surface, based on the simulation of tensile stress in the paint layer. Cracks are created and propagated on the surface using a 2D grid to access local surface properties. These properties can be represented in different ways, either procedurally, as a constant, or a texture. For example, if we wish to reproduce a specific type of crack pattern, textures can be used to define surface properties. The cracks are then guided by texture value and gradient. Simple intersection of cracks is handled, and fine details of cracks are maintained separately, thus increasing simulation speed. The simplicity of our crack propagation algorithm allows us to compute a simulation efficiently, giving the user fine control and the ability to easily experiment to obtain the desired result.

In addition to crack propagation, we generate and render the curling effect of peeling on our two-layer representation based on the simulation of adhesion between these layers. To our knowledge this is the first treatment of adhesion and curling for peeling in computer graphics. This effect adds important and convincing visual detail.

Curling is computed at an appropriate level of detail of the crack path as paint loses its adhesion with the underlying surface. As a result, peeling can occur in different directions as simulation proceeds, which corresponds to observed physical reality. It is important to note that we do not modify the underlying geometry to add this curl-ing effect, but simply superimpose additional geometry on the underlying object.

The entire approach, crack propagation and peeling, is relatively simple to use with standard rendering packages.

Previous Work

Simulating "imperfections" or "weathering" is a relatively young area in computer graphics. To our knowledge, the first such approach was that of Becket and Badler in 1990 [START_REF] Becket | Imperfection for realistic image synthesis[END_REF] which described a general system treating scratches, stains, splotches, and rust. Dorsey and colleagues have developed a series of methods treating a number of different phenomena, such as patinas, weathering due to flow of liquids, and weathering of stone (a survey is presented in [START_REF] Dorsey | Digital materials and virtual weathering[END_REF]). These models are typically based on involved physical simulations, with relatively expensive computation times and the need to determine physical parameters. Such models are indispensable for correct and accurate simulation of weathering phenomena, and have produced images of striking realism and beauty. Nonetheless, their usage can be cumbersome for some applications, and more simplistic, albeit less accurate, approaches can also be useful. Approximate models can often be more appropriate, for the computer graphics industry [START_REF] Warniers | Dirty pictures[END_REF] for example.

In this paper we treat generation and propagation of cracks on a two-layer model. Dorsey and Hanrahan [START_REF] Dorsey | Modeling and rendering of metallic patinas[END_REF] developed a model based on multiple layers to simulate metallic patinas, in the context of texture generation tools. Sophisticated shading models were used, resulting in great realism. Merillou et al. [START_REF] Merillou | Surface scratches: Measuring, modeling and rendering[END_REF] presented a two-layer model to simulate scratches on a surface. Their approach does not add geometry to represent the scratches, but simply defines an appropriate anisotropic reflectance function to give the correct visual impression of scratches.

Another aspect of our work deals with the formation and propagation of cracks which can then peel. Terzopoulos and Fleischer [START_REF] Terzopoulos | Modeling inelastic deformation: Viscoelasticity, plasticity, fracture[END_REF] developed an early physically-based model of fracture, as part of the treatment of inelastic behavior in general. O'Brien and Hodgins [START_REF] O'brien | Graphical modeling and animation of brittle fracture[END_REF] developed a model for fracture on an adaptively subdivided tetrahedral mesh. A more complete physically-based simulation is used, and cracks are propagated in 3D. This model concentrates on the dynamics and the motion of fracture, yielding impressive visual results, but at a non-negligible computational cost.

Hirota et al. [START_REF] Hirota | Generation of crack patterns with a physical model[END_REF] simulated cracks on a surface caused by surface contraction. Their approach is based on physics to a certain extent, and is modeled as a network of springs. Cracks appear as springs break, and thus their paths are constrained by the spring network. The nodes are squares or hexagons, resulting in different types of crack shapes. The visual quality of the cracks is directly affected by aliasing across the network cells. Neff and Fiume [START_REF] Neff | A visual model for blast waves and fracture[END_REF] simulated the generation of cracks due to a blast wave. A seed crack is used and cracks are propagated, releasing energy as they proceed. This enables the creation of new crack branches. A simple random variation scheme on the branching direction generates visually convincing results.

For our crack propagation method, we have been inspired by the ideas presented in Gobron and Chiba [START_REF] Gobron | Crack pattern simulation based on 3D surface cellular automata[END_REF]. Their method however deals only with cracks and does not simulate adhesion or peeling. As with our approach, cracks are propagated across a 2D grid, using a relaxation method. Certain aspects however of their simulation are more simplistic, as their approach does not treat peeling (adhesion, curling, etc.). Since between 100,000 to a million cells are required for the simulation of their cracks, computation time is quite expensive, ranging from 20 minutes to several hours.

Finally, in the general system presented by Wong et al. [START_REF] Wong | A geometry dependent texture generation framework for simulating surface imperfections[END_REF], peeling is treated along with rust and dust. This approach simply modifies the surface color based on a threshold function for the choice of the appropriate layer color, and thus the geometric effect of cracks, loss of adhesion, and surface curling are not simulated. No simulation of crack creation and propagation is performed, and the user must manually position "peeling light sources".

Physical Basis

Paint and other thin layers are applied on many surfaces for both decorative and protective purposes. The condition of the paint layer over years of exposure depends on many factors including the paint itself (e.g., latex or oil), the surface on which it is applied (e.g., material such as wood or metal, roughness, penetration [START_REF] Long | Treatise on Coatings[END_REF]), as well as environmental factors (e.g., UV light, oxygen, pollutants). All these factors affect paint differently during its lifespan [START_REF] Fischer | Organic Protective Coatings[END_REF]. Four properties are of major importance to paint cracking and peeling: stress, strength, elasticity, and adhesion to the base surface.

One factor is inherent to every kind of paint: as paint dries, it shrinks, producing stress in the paint layer [START_REF]Adhesion and Adhesives Fundamentals and Practice[END_REF]. In this context, thinner paint layers are preferable to thicker ones since they induce less stress. A dirty or wet base surface and insufficient stirring of the paint also results in weaker paint layers [START_REF] Cyril | Paint Technology Manuals[END_REF]. The method of application (e.g., brush, spray, roller) also affects the layer quality [START_REF] Cyril | Paint Technology Manuals[END_REF]. Even after paint has dried, its properties are not fixed as it will slowly become weaker because of interaction with the elements of its environment. The principal cause of damage is moisture [START_REF] Sward | Paint Testing Manual -physical and chemical examination of paints, varnishes, lacquers, and colors[END_REF], but many other envi-ronmental factors [START_REF]Aspects of Adhesion -proceedings[END_REF][START_REF] Cyril | Paint Technology Manuals[END_REF] such as water, UV light, temperature variations (contraction or dilation of the base surface), pollutants, abrasion, and impacts all weaken the paint layer.

When paint is weakened, the physical properties of the layer change. Paint is to some extent elastic, that is it can be stretched thus redistributing surface stress. As years pass, the elasticity of paint is likely to decrease. Long term exposure to pollutants for example can greatly decrease elasticity. The paint also has a specific strength, which determines the greatest stress it can handle before it starts to break. As is the case for elasticity, exposure to harsh elements or UV light can reduce the strength of the paint, making it prone to develop cracks and peels. When a crack forms, it is readily visible since it enables the paint on either side to shrink back to a state with less stress. This release of tensile stress occurs perpendicularly to the crack [START_REF] Gobron | Crack pattern simulation based on 3D surface cellular automata[END_REF]. Adhesion is the most important property of paint. We define it as the ability of paint to attach itself to the underlying surface. Moisture, base surface contraction, or rust can all greatly reduce the adhesion of paint. When the paint shrinks as cracks release tensile stress, a shearing force is introduced at the interface between the paint and the base surface, pulling the paint away from its initial position. As paint loses adhesion with the base surface, it can either peel by slowly curling away [START_REF] Fischer | Organic Protective Coatings[END_REF][START_REF] Cyril | Paint Technology Manuals[END_REF], or flake [START_REF] Cyril | Paint Technology Manuals[END_REF] and eventually fall off the base surface.

As we can see, peeling is a complex phenomenon that spans over years. It is not fully understood, since most of the (physical) research focuses on preventing its occurrence [START_REF] Fischer | Organic Protective Coatings[END_REF]. Moreover, most of its parameters are difficult to measure. In this paper, we develop a model inspired by, rather than strictly simulating, the physics of the phenomenon. Our goal is ease of control, allowing the user to obtain effects that are visually appealing. In our simulation process, the paint strength and stress determine where cracks appear on the surface. Reduction of stress around the cracks is computed by considering the elasticity. Finally, curling of the paint layer as it peels from the surface is computed with the adhesion. In the next two sections, we detail the structures and algorithms of our empirical peeling model.

Simulation

Overview

Our simulation model is simple and summarized in Figure 1. A surface is defined by two layers that we will refer to as base layer and paint layer. 1 When tensile stress (Figure 1-i) is too high compared to the paint strength, 1 While the top layer can be of a different material than paint (e.g., wallpaper), it is simpler to refer to it as paint. a crack can appear (Figure 1-ii) and expand on the paint layer according to the local surface properties. By allowing the paint to shrink (Figure 1-iii) to a lower energy state, cracks reduce the tensile stress (Figure 1-iv) in the paint layer, producing an increase in shearing stress (Figure 1-v), and can cause adhesion to the base layer to be lost (Figure 1-vi). Thereafter, the paint is free to peel away from the base layer (Figure 1-vii).

Our simulation process is presented in Figure 2. Several cracks evolve on a surface and can intersect each other. To produce an image, the cracks and their parameters are passed to the rendering module (described in Section 5) where their shape is computed and rendered. This can be done at the end of the simulation or at time step intervals.

The surface properties are represented by a 2D grid structure. Each grid cell stores its local paint properties (adhesion, stress, etc.) and bilinear interpolation is used to compute properties for any location on the surface. Every surface property can have a value that varies with respect to direction on the surface. We represent such directional variations by storing samples for a predefined number of directions. 2 The resolution of the grid is responsible for many features (propagation speed, crack details, efficiency, etc.) that will be explained in the next sections.

Cracks 1.

New Cracks and Propagation

We assume that if nothing else held it in place (e.g., neighbor cells, adhesion to base layer), the paint layer in a cell would shrink to a stable lower energy state. The tensile stress σ is the difference in energy (force) between the current state and the stable state of the paint layer.

The tensile break strength S b is the maximal force the paint layer can handle before it breaks. Thus, if the stress σ exceeds the strength S b , the paint will break. A new crack appears within the grid cell with the highest stress to strength ratio σ/S b (Figure 3(a)).

The creation and propagation of cracks are strongly influenced by the stress and strength properties of paint. These properties can be dependent on the direction on the surface. For instance, wood fibers in the base layer can induce directional variations in stress and strength. When creating a crack, the orientation of the initial segment is perpendicular to the orientation of the maximal stress to strength ratio σ/S b . Since we represent the directionality with samples for fixed directions, we use random distributions based on the sample values to compute the orientation of the maximal value of a property.

A propagation step moves each crack end point by a linear segment of length smaller or equal to the grid size, thus guaranteeing that we do not miss any grid cell contribution along this segment (we typically use a uniform distribution between 0.5 and 1.0 times the grid cell size). This propagation defines the path followed by a crack (Figure 3(b)). When propagating cracks, the stress σ is compared to another paint property, the tensile crack strength S c . A crack stops when the local stress σ is below the strength S c , or when it intersects another one (see Section 4.5). As when computing the orientation of new cracks, each of the two end points of a crack has a tendency to propagate in the direction perpendicular to the local highest stress to strength ratio σ/S c . We also mix the computed orientation with the previous crack direction to prevent changes in direction which are too abrupt.

The top arrow in Figure 2 presents the main simulation loop, showing how the cracks and paint are adjusted in each step. We propagate the cracks until extinction of movement, before creating a new crack (bottom arrow in Figure 2). This simulates the fact that it often requires less energy to propagate a crack than to create a new one [START_REF] Smith | Fast and controllable simulation of the shattering of brittle objects[END_REF]. It also prevents the creation of many separate cracks within regions of high stress, and results in more continuous and realistic cracks. This strategy also greatly reduces the number of intersections that must be computed. In the same spirit, our crack strength S c is lower than the break strength S b to ensure that cracks will not stop after only a few propagation steps (we typically use a S b that is equal to 1.5 times S c ).

Relaxation of Tensile Stress

As we saw in Section 3, cracks reduce tensile stress in the paint layer. We simulate this relaxation process by reducing the tensile stress of the cells close to each crack segment.

The elastic relaxation distance property D r determines how far from the segment we should release stress. A tensile deformation δ of the paint layer, corresponds to the release in tensile stress. This is shown in Figure 4 The paint tensile stress released by a crack should be roughly perpendicular to the crack. The direction of relaxation for a point on the paint layer is interpolated from the two bisector vectors for each segment (see Figure 4(a)). Compared to Gobron and Chiba [START_REF] Gobron | Crack pattern simulation based on 3D surface cellular automata[END_REF] where only the direction perpendicular to the segment is considered, our approach avoids unwanted overlap (when consecutive segments form a concave curve) or cells that can receive no relaxation (in the convex case).

The set of cells close to the crack end points will be used to compute the paint tensile stress on the next propagation step. Since they will serve in the crack propagation computation, relaxing their stress would cause the crack to release stress ahead of itself. This would be incorrect with respect to the real physical phenomenon, and could cause the crack to stop for no apparent reason. These cells are thus excluded from the relaxation computation.

As described so far, the crack propagation has resulted in a set of piecewise linear segments indicating the central path of the crack. However a crack is an opening in the paint layer. As the paint layer is relaxed by the cracks, it shrinks locally by the deformation lengths δ. This causes the crack to widen, and the base layer to become visible. The deformation length δ of each relaxed cell is added to the crack width property of the crack segment(s) responsible for this relaxation. For each crack segment, an independent width is computed for each side of the crack. This results in two sets of corresponding piecewise linear segments with independent crack width values defining the crack sides. The crack widths are displayed in Figure 4(a) as segments on both sides of the crack path; linear interpolation is used to ensure continuity from one segment to the next.

Adhesion Modification

We next explain how we determine when curling of the peels can occur. As the paint layer shrinks by relaxation, it causes a lateral deformation of the layer, which induces a shearing stress τ counteracted by the paint adhesion. The adhesion strength S a represents the maximal force (due to shrinking) that the paint layer can tolerate before adhesion to the base layer is broken. If the shearing stress τ increases above the adhesion strength S a in any direction, the paint layer in this cell is considered as "detached" from its base layer for all directions. This computation only accounts for deformation local to the cell, without considering deformation of adjacent cells. It also is a binary operator, so even if adhesion to the base layer would be lost in the vicinity of the crack, when considering the whole cell, adhesion to the base layer can remain.

It should be noted that when the relaxation of tensile stress (Section 4.3) encounters cells detached from the base layer, it sums the deformation lengths δ, but extends the relaxation to a larger region by considering the detached cells as having a null size. The linear decrease of the relaxation deformation is thus moved further away from the crack, due to the loss of adhesion.

When adhesion is lost, the paint can curl away from the base layer. To compute this peeling, each crack segment keeps, in an adhesion width, the distance to the closest cell (with respect to the region defined by its bisectors) that is still attached to the base layer (Figure 4(b)).

Intersection

As a crack propagates through the grid, every cell within a distance of half a cell is marked with its crack ID. We do not allow two cracks to traverse the same cell. Independent cracks are therefore always separated by a minimal distance related to the grid resolution. Self-intersection of a crack cannot be detected by the crack ID but can be handled by separately testing for intersection of the crack with itself.

When a crack enters a previously marked cell, the new segment is guaranteed not to cross over a previous crack segment. At this point the intersecting crack has not yet reached the intersected (previous) crack (see Figure 5) but we force the new crack to intersect the previous crack, as if the weakened paint layer suddenly broke. We thus need to connect the two cracks. The intersection point is the "closest" point on one segment of the intersected crack. First we eliminate the segments of the intersected crack for which the intersecting segment end point is not contained in the space confined within the intersected segment and bisector directions. For the remaining segments, we identify their respective intersection point and keep the segment with the closest intersection point.

We then add two segments to the intersecting crack to connect it to the intersected crack. Elastic relaxation of tensile stress is computed with respect to the added segments to determine the crack width of the intersecting crack where the intersection occurred. After relaxation, the intersected crack is split at the intersection point with respect to the intersecting crack width. The newly formed (opened) crack smoothly joins each side of the intersecting segment.

Segment Merging and Crack Detail

As the paint loses adhesion to the base layer, elastic relaxation reaches further away from the crack and the peels also curl further. Similarly to elastic relaxation, the peels curl away from the base layer in the direction perpendicular to the crack path (in our case, the directions of the interpolated bisectors of Figure 4(a)). As the relaxation or peeling reach further away from the segment, they should consider a larger portion of the crack. The computations thus have to consider the crack at a coarser level of detail as the adhesion width increases.

As explained earlier, the crack path traverses the paint layer with changes in direction that are not too abrupt. When little adhesion is lost, the adhesion width W a vectors (Figure 4(b)) are short and the peel will curl very locally. When the loss of adhesion extends much further, the adhesion width W a vectors are much longer and they can even intersect each other, yielding unexpected results as illustrated in Figure 6(a). It will also be difficult to create geometry for curling such peels. We solve these problems by merging together adjacent segments. This is similar to considering the crack path at a coarser level of detail, disregarding high frequency variations of the crack path (see Figure 6(b)). Merge operations are controlled by a set of metrics that determine if segments should be merged or not. These metrics consider different factors such as the segment lengths with respect to their adhesion width, and the difference in direction between consecutive segments.

Of course, we still want to keep the high frequencies of the crack path to compute the geometry of the peels, even though we do not use them to compute the relaxation or the peels curling directions. We do this by keeping detail records that store the difference (in crack and adhesion widths, as well as position) between the merged segments and the coarser level segment (see Figure 6(b)). This information is used when generating the geometry and for rendering the final peels.

Interaction

We implemented this technique in an interactive system in which the user can specify the different properties and simulation parameters. Running on an Octane2 with a MIPS R12k 400 MHz CPU, the propagation is interactive (many new segments per second) for up to about 30 cracks. This enables the user to observe the evolution and change the parameters early in the simulation process.

The specification of paint layer parameters can be done with properties constant over the entire surface, isotropic properties, or directional properties. We found that textures to control the location of high and low strength are particularly useful (in Figures 9 and10, the second image is the texture used to define the tensile strength property). To guide the crack to higher stress regions and to retrieve directional information from the texture, gradient operators are computed.

Rendering

Once the paths and widths of the cracks have been determined by the simulation, we pass this information to the rendering phase. It generates the peeling geometry on top of the original geometry, locally along the crack paths, and computes the final image with proper shading. The final rendering phase could be done with any highquality renderer; we use Alias |Wavefront's Maya™.

After the crack is formed, new portions of the paint layer (sides and underneath the crack) get exposed to environmental attacks (e.g., oxygen, UV light). The differential stress between the top and the bottom within paint layer will cause this layer to curl. Instead of simulating this process in detail, we simply rely on a controllable approximation shown in Figure 7.

We compute the outward curling of the paint layer using a spiral-shaped geometry. The end of the curl forms an angle θ with the base surface (Figure 7(a)). This angle is user controlled and vanishes to zero at the end points of a crack.

To generate and render the effect of peeling, we create a mesh of micro-polygons with elements smaller than a user specified size (Figure 7(b)). A crack segment is subdivided in as many sub-segments as needed to respect this maximal size. It is also subdivided at each location where detail information is available (Section 4.6). Along the crack, the geometry is built generating one quadrilateral strip per sub-segment. To ensure smooth shading of the peels (Figure 7(c)), the surface normals of the micropolygons are averaged and stored at every mesh vertex location.

This technique is simple and efficient, and can thus be generated on the fly in most standard rendering environments such as Maya™, or could be generated in a preprocess to create the peels geometry into a file added to the final scene. The approach has the advantage that the additional geometry is "glued" to the base surface, and thus does not require remeshing or modifying the objects. Some care must however be taken to add a small offset with respect to the original surface to avoid "surface acne" artifacts.

In the rendering process we also need to compute the shading of the base geometry. It must show the base layer color inside a crack and the paint color outside. This corresponds to evaluating if the point to be rendered is inside or outside a crack. We use a 2D grid structure to accelerate this test. We implemented this inside a plugin shader with user input for paint and base layer shaders and output of the appropriate color for the point being shaded. This is very flexible, allowing the use of any shader, 2D/3D texture, or other utilities to send the desired color to our shader. Shadows and other illumination effects are handled by the renderer.

Results

Figures 8,9, and 10 present real photographs (left) and synthetic results (rightmost two) of our system for a garage door, wood shutters, and an indoor wall. When the peels are viewed from a closer viewpoint (right), the shading and shadowing clearly show how both the curl geometry and the base layer appearing inside the crack are important to the visual quality of the results. In Figures 9 and 10, the second image from the left shows the tensile break strength S b texture used, with darker representing lower strength and lighter representing higher strength.

Table 2 summarizes the statistics for these three scenes. Our simulation system is time and memory efficient, requiring about 35 MB of memory for the presented simulations. The computation of the elastic relaxation is the most time consuming part (about 80% of the total simulation time). The synthetic images were rendered with high quality anti-aliasing at a resolution of 640×480 pixels.

We have presented a simple, fast, and easy to use approach for crack generation and peeling on layers such as paint. If desired, the user may control the generation of cracks with textures to describe the different surface properties. The crack propagation is performed on the surface, with a 2D grid to query and update the various surface properties. As cracks propagate, they remove constraints in the paint layer, allowing it to slightly shrink. We simulate both this shrinking and the related loss of adhesion between the paint and the base layer. The result of the crack propagation simulation is passed to the rendering phase, which generates geometry to represent the curling process of peeling and uses the appropriate shader when inside or outside of the cracks. Our approach produces visually appealing results, which are close to observed phenomena.

In future work, we will treat pieces of peels which break off, as often happens in reality. To add further realism, we should be able to treat more than two layers, for example when paint is added over wallpaper, etc. Another important direction is the treatment of curved surfaces, which complicates both the crack propagation algorithm, in particular the determination of the directions of propagation, and the generation and rendering of curl geometry.

Finally, our overall goal is a complete system to treat a wide variety of deterioration effects, in an integrated and easy to control manner. This will require careful determination of parameter choices, which will undoubtedly be a combination of automatic and user-assisted solutions.

Figure 1 :Figure 2 :

 12 Figure 1: Peeling model.

Figure 3 :

 3 Figure3: Snapshots from our system for creation and propagation of cracks. (a) The initial stress to strength ratio σ/S b scaled such that a ratio of 1.0 corresponds to the size of a cell and presented as lines for the four sampled orientations. A new crack is created perpendicular to the maximal ratio. (b) The propagation of a crack (note that σ/S c > σ/S b and that no relaxation was computed on the stress σ used to compute the presented ratios).

Figure 4 :

 4 Figure 4: Snapshots from our system for relaxation and adhesion. (a) The cumulated deformation lengths δ caused by the relaxation. The bisectors on both sides of the crack represent the crack width and determine the region affected by each segment. (b) The shearing stress to adhesion strength ratio τ /S a . Grey cells have lost adhesion to the base layer, and arrows on each side of the crack represent how far adhesion to the base layer is lost with respect to the crack.

Figure 5 :

 5 Figure5: Intersection between two cracks, the intersected (previous) crack and the intersecting (new) crack. The identified end point corresponds to the intersecting crack end point before the intersection is detected. Cells with a light grey background are those marked by the intersected crack ID, while the cells with a dark grey are those marked by the intersecting crack.

Figure 6 :

 6 Figure 6: Merge of segments at a coarser level of detail. (a) Before merging, the peeling regions of some segments overlap. (b) After merging and recomputation of the adhesion widths, peeling regions respect the linear behavior of the crack path. The detail records for the crack path position (black dot on grey circle) and the adhesion width (white dot on grey circle) are also shown.

Figure 7 :

 7 Figure 7: Geometry of the peels. (a) Cross-section of the curling geometry with three sections and angles θ = 90 o , φ = 90 o /(1 + 2 + 3) = 15 o . (b) Generated mesh. (c) Rendered result.

Table 1 :

 1 propagate on the surface and are represented Paint and crack properties. The equations are simple so their effect is intuitive to understand. All the user specified properties can be constant or vary over the surface. Constant terms should be added in the equations of δ and τ so that the units match, but were not included for clarity.

	Paint Properties	Units	Source	Use
	σ Tensile stress	force / length user specified	creation, propagation, relaxation
	S b Tensile break strength	force / length user specified	creation
	S c Tensile crack strength	force / length user specified, S c < S b	propagation
	δ Tensile deformation	length	δ = σ.cellSize.(1 -distance Dr	) relaxation
	D r Elastic relaxation distance length	user specified	relaxation
	τ Shearing stress	force / length τ = δ/cellSize	adhesion
	S a Adhesion strength	force / length user specified	adhesion
	Crack Properties
	W c Crack width	length	δ i	rendering
	W a Adhesion width	length	distance	merging, rendering

by connected segments forming piecewise linear curves. They are independent of the 2D grid representation which only serves to query the local surface properties. Typically, crack segments are smaller than the grid cell size. The surface and crack properties used throughout Sections 4 and 5 are presented in Table

Table 2 :

 2 Statistics for the scenes of Figures 8, 9, and 10. † Simulations were computed on an Octane2 with a single MIPS R12k 400 MHz CPU. ‡ Rendering was done using the multiprocessing capabilities of Maya™ on an Onyx 3400 with 16 MIPS R12k 400 MHz CPUs.

		Door	Shutters	Wall
	Simulation time †	75 min	20 min	3 min
	Number of cracks	200	200	30
	Number of segments 2900	1500	700
	After merge	1200	850	250
	Grid resolution	50 k	36 k	30 k
	Rendering time ‡	3 min	3 min	3 min
	Number of polygons 400 k	150 k	130 k

Since these properties are symmetrical, they need to be stored only within 180 o . We typically use four directions (0 o , 45 o , 90 o , 135 o ).

Acknowledgements

Our work was supported by grants and scholarships from FCAR, NSERC, and MRI-MEQ. We thank Alias |Wavefront for the Maya™ system, the reviewers, and everyone from LIGUM, REVES, and iMAGIS. We would like to thank A. Reche for his help with the photography.

Real