prof. Ioan Mircea Risteiu

assoc.prof Ileana

Mihail Pascariu

Digital Alba Iulia. The System Integrate for Regional E-government SIREgo Integrating Information Technology and Social Science Research for Effective

Keywords: regional development, information technology, social science, e-government, GIS, digital entities

The Digital Government Society (DGS), or Intelligent Society is a global multi-disciplinary state particular "to society" contributors and practitioners engaged in and committed to democratic digital government. Digital (or electronic) government fosters the use of information and technology to support and improve public policies and government operations, engage citizens, and provide comprehensive and timely government services. DGS equips its members with a professional support network focused on both education and effective practices that nurture technical, social, and organizational transformation in the public sector. The DGS welcomes members from all sectors, endorses diverse, multi-, and interdisciplinary research undertakings relevant to both theory and practice, and strongly encourages practitioner-researcher exchanges at local, regional, national, and international levels. The purpose of this project is to delimit the Digital Entities, as "images of reality", made of Data, (the bits -zeros and ones-put on a storage system, Information (the attributes used to assign semantic meaning to the data), until Knowledge (the structural relationships described by a data model or semantic relationships between attributes). Every digital entity requires information and knowledge to correctly interpret and display

I. INTRODUCTION 1.1 Local/regional government vs communication technology

As literature explains, e-Government explains that the government's best practice use of information and communication technology (ICT) in order to exchange information and services with citizens, businesses, and other arms of government. e-Government includes legislature, judiciary, or administration, culture and recreation, education and promotion in order to improve internal efficiency, the delivery of public services, or processes of democratic governance. The primary delivery models are Government-to-Citizen or Government-to-Customer (G2C), Government-to-Business (G2B) and Government-to-Government (G2G). The important benefits of e-government include improved efficiency, convenience, and better accessibility of public services. [START_REF] Baller | Structural covariates of U.S. county homicide rates: Incorporating spatial effects[END_REF] has proposed a structure of the three modules. Each one has a data distribution as is shaped in the figure 1.

Available information technology

While e-government is often thought of as "online government" or "Internet-based government"many non-Internet based "electronic government" technologies can be used in this context, including telephone, fax, PDA, SMS text messaging, MMS, and 3G, GPRS, WiFi, WiMAX and Bluetooth. Other technologies can include CCTV, tracking systems, RFID, biometric identification, road traffic management and regulatory enforcement, identity cards, smart cards and other NFC applications; polling station technology (where non-online evoting is being considered), TV and radio-based delivery of government services, email, online community facilities, newsgroups and electronic mailing lists, online chat, and instant messaging technologies. There are also some technologyspecific sub-categories of e-government, such as mgovernment (mobile government), u-government (ubiquitous government), and g-government (GIS/GPS applications for e-government). There are many considerations and potential implications of implementing and designing e-government, including disintermediation of the government and its citizens, impacts on economic, social, and political factors, and disturbances to the status quo in these areas. In countries such as the United Kingdom, there is interest in using electronic government to re-engage citizens with the political process. In particular, this has taken the form of experiments with electronic voting, aiming to increase voter turnout by making voting easy. The UK Electoral Commission has undertaken several pilots, though concern has been expressed about the potential for fraud with some electronic voting methods [START_REF] Bailey | Interactive Spatial Data Analysis[END_REF].

Data quality ensuring

Data quality can be assessed through data accuracy (or error), precision, uncertainty, compatibility, consistency, completeness, accessibility, and timeliness as recorded in the lineage data. Spatial error refers to the difference between the true value and the recorded value of non-spatial and nontemporal data in a database. Attribute error is more complicated than other types of spatial errors. It is related to scale of measurements. The upper figure is extremely consistent from this point of view.

On the other hand, at one scale of measurement, the difference may be regarded as error while not at another scale. For example, an elevation of 497 m recorded in the database with its true value being 492 m will be considered erroneous at the ratio and interval scales but accurate in a general category such as an elevation class between 450 and 500 m which is at the nominal scale. However, sometimes the true value is not known, error can not be evaluated. Under such circumstances, uncertainty is used. Statistically, we use the average from multiple measurements to estimate the true value and the standard deviation of the multiple values as an indicator of the level of uncertainty.

Avoiding inaccurate information

In order to avoid situations like we particularly mentioned a method which is used in remote sensing image classification, where a relatively complete procedure exists for classification error analysis could be used here. Related to this, we -determine the sampling method for ground truth data collection; Methods include systematic sampling, random sampling, stratified sampling and systematic unaligned sampling, etc. -determine the sample size; -determine the attribute of sample location; this is usually done by field survey or the use of more accurate data sources such as aerial photo interpretation; -compare sample data with classification data and establish the contingency matrix; -calculate various types of errors from the contingency matrix.

This can be applied to the determination of any type of attribute error at the nominal measurement scale, and based on this, we are building a new data processing structure for the designed GIS.

II. BUILDING A DIGITAL GOVERNMENT SOCIETY

The development of specialized software for spatial data analysis has seen rapid growth since the lack of such tools. The main objective of the software is to provide the user with a natural path through an empirical spatial data analysis exercise, starting with simple mapping and geovisualization, moving on to exploration, spatial autocorrelation analysis, and ending up with spatial regression. Let's start with a basic data processing example. By using a single database variable, as elevation positioning, we have detected a logical inconsistency. The standard error difference between area mean and population mean caused bias in estimating population mean.

The first issue of this project is to delimit the Digital Entities, as "images of reality", made of Data, (the bits -zeros and ones-put on a storage system, Information (the attributes used to assign semantic meaning to the data), until knowledge (the structural relationships described by a data model or semantic relationships between attributes). Every digital entity requires information and knowledge to correctly interpret and display.

The second issue is to build an HCIbased application with common facilities:

- Depending on the size and complexity of the Web application, the user may be able to handle anywhere from 10 to thousands of concurrent users.

compatibility
The number of concurrent connections to your Web server(s) will ultimately have a direct impact on your site's performance. Therefore, designer performance objectives must include two dimensions:

-the speed of a single user's transaction; -the amount of performance degradation related to the increasing number of concurrent users on your Web servers.

On the other hand, load management refers to the method by which simultaneous user requests are distributed and balanced among multiple servers (Web, ColdFusion, DBMS, file, and search servers). Effectively balancing load across your servers ensures that they do not become overloaded and eventually unavailable. There are several different methods that you can use to achieve load management: hardware-based solutions; softwarebased solutions, including round-robin Internet DNS or third-party clustering packages; Hardware and software combinations; each option has its own distinct merits.

II.2 The context for GIS

In many respects, GeoDa is a reinvention of the original SpaceStat package (Anselin 1992), which by now has become quite dated, with only a rudimentary user interface, an antiquated architecture and performance constraints for medium and large data sets. The software was redesigned and rewritten from scratch, around the central concept of dynamically linked graphics. This means that different "views" of the data are represented as graphs, maps or tables with selected observations in one highlighted in all. In that respect, GeoDa is similar to a number of other modern spatial data analysis software tools, although it is quite distinct in its combination of user friendliness with an extensive range of incorporated methods. A few illustrative comparisons will help clarify its position in the current spatial analysis software landscape. Linking in these implementations is constrained by the architecture of the GIS, which limits the linking process to a single map (in GeoDa, there is no limit on the number of linked maps). In this respect, GeoDa is similar to other freestanding modern implementations of ESDA, such as the cartographic data visualizer, or cdv (Dykes 1997), GeoVISTA Studio (Takatsuka and Gahegan 2002) and STARS (Rey and Janikas 2004). These all include functionality for dynamic linking, and to a lesser extent, brushing. They are built in open source programming environments, such as Tkl/Tk (cdv), Java (GeoVISTA Studio) or Python (STARS) and thus easily extensible and customizable. In contrast, GeoDa is (still) a closed box, but of these packages it provides the most extensive and flexible form of dynamic linking and brushing for both graphs and maps. For all this platforms the main subject is focused on common spatial autocorrelation and statistics. With the original implementation for all users, the data is assumed to be local, on a heterogeneous network, or obtainable via an external program. With the inclusion of DEI's capabilities into applications, it is possible to access data via gateways located virtually anywhere on the Internet. The data can be accessed piecewise allowing large datasets to be supported.

II.3 The issue today's developers

The design of present projects consists of an interactive environment that combines maps with statistical graphs, using the technology of dynamically linked windows. It is geared to the analysis of discrete geospatial data, i.e., objects characterized by their location in space either as points (point coordinates) or polygons (

III. NEW IMPLEMENTATION IN DATA REPRESENTATION

In [START_REF] Luc Anselin | An Introduction to Spatial Data Analysis, Spatial Analysis Laboratory[END_REF] there is shown the situation when the polygon data requires more information to represent the spatial aspect of the DataSet; the non-spatial aspect of the data being represented is the same. Instead of a MeshCell type of PointCell, now the type is Polygon. The Field no longer has 'vertices' dependency, but now has 'cell' dependency. Each 'row' of non-spatial data (one element from the data side of each Field) is dependent on a cell, not just one vertex. The different types of Connections that IGOR provides manage how the vertices make up the Polygon, allowing for more efficient representation of the vertices in memory. In the closed representation (figure 2) is a simple diagram that helps to see how the usual structures represent the polygon data [START_REF] Luc Anselin | An Introduction to Spatial Data Analysis, Spatial Analysis Laboratory[END_REF]. This structure helps the designers to understand the data processing procedure. Each data entity is queried through different standards and criteria.

IV. PROPOSED PROCESSING METHOD

Spatial modeling and mapping of local/regional data is a relatively represented, and is placing higher demands on data quality than traditional, point-based analyses. Previous experience with spatial modeling projects has demonstrated the problems experienced with insufficiently qualitycontrolled (QC) data sets. A single outlier has the capability of causing spatial abnormalities, sometimes quite large and obvious. The outlier may not appear unusual in a time series plot for a single station, but spatial analyses can often show the inconsistency of the outlier in relation to values from nearby stations.

While a spatial analysis can be an effective way of identifying outliers, the detection and correction of such errors after the spatial data set has been produced can be time-consuming and inconsistent.

Each spatial grid must be visually examined for possible outliers, and special calculations must be made to help show outliers where they may be otherwise hidden, such as areas of complex terrain. When possible errors are found, they must be individually verified, and the erroneous grid reproduced without the offending data point.

Temporal relationships between pairs of data coordinates have been calculated in anomaly space, that is, the variations remaining after the means had been subtracted off. While this is a generally a more stable way to express continuous and dynamic data, it requires that a stable, long-term mean value be available.

For data supplying, because of heterogeneous primary data sources, we have used:

-import/ export facilities; -sessions for executing executable script for o Standard queries o Real time designing, running and killing temporary databases

IV.1 Integration of the real time error detection into Internet-based application

This part of the project is focused in combining real time error detection with temporary databases by running server side scripting programs, and client side fast processing enabling. Today, networkenabled consumer devices such as mobile phones or PDAs allow new user groups to take advantage of the Internet. Besides the limitations posed by these new devicesincluding restricted screen real estate, small bandwidth, numeric keyboard, etc.there are also new opportunities for application developers. For instance, information services may be designed to consider the users' locations, time of the day, etc.; portability and usability of mobile devices allow for instant and ubiquitous access to information; and, content and access structures can be delivered considering the users' interests and preferences.

-Structured documents and style sheets: This approach implies that content is stored independently of the target language by using for instance the Extensible Markup Language (XML, [START_REF] Luc Anselin | An Introduction to Spatial Data Analysis, Spatial Analysis Laboratory[END_REF]) as a declarative specification language. XSL (Extensible Stylesheet Language, [START_REF] Bailey | Interactive Spatial Data Analysis[END_REF]) style sheets can then be used to define the replacement of XML element types with tags specific to the target language. The transformation is done dynamically either on the server side upon page requests or by clients capable of doing the transformation themselves.

This The modeling techniques upon which they are based are targeted at modeling functional behavior and the principle of separation of concerns is applied by organizing an application as a set of cooperating functional units (where, depending on the particular programming paradigm, a functional unit may be a procedure, a module, an object, a class, etc) see the figure 4.

IV.2 Viewing data on different computing systems

First approach for building a behavioral talk between user-side and server-side is focused on dynamic creation of content.

State-of-the-art applications which support personalized presentation of content demand for dynamic content generation. This is supported on the one hand by parameterizable application manager and on the other hand by the possibility to create the content tree at runtime by the application (figure 5). This procedure is based on checking code following user setup. The applied method is a form of automatic code generation where both the starting and the end code are written in the same language. By putting the aspect oriented programming algorithm into the client-server application for GIS management, the logical diagram looks like in figure 6. The generated beans of information became apparent that not only primitive widgets needed to be supported but also more sophisticated ones (figure 7). With these so-called "group widgets" it is possible to represent coherent user interactions.

Figure 1

 1 Figure 1 Implemented structures

Figure 2

 2 Figure 2 Data representation for errors detection

Figure 3 Figure 4 Figure 5

 345 Figure 3 Static and instance table

Figure 6 Figure 7

 67 Figure 6 Diagram structure for error checking

  All we know about the importance of having Web servers that scale well, but what exactly scalability is difficult to define. Simply, scalability is a Web server's ability to maintain a site's availability, reliability, and performance as the amount of simultaneous Web traffic, or load, hitting the Web server increases. The major issues that affect Web site scalability include: "Performance" and "Load management". Web application architects and developers must design and code an application with performance in mind. Once the application is built, various administrators/users can tune performance by setting specific flags and options on the database, the operating system, and often the application itself to achieve peak performance.

			and applications has to be used for
			access facilities
			-fast data sharing, is demanded for a
			fast access to data and complex
			sharing
			II.1 Scalability in www
	file format with neutral
	developing environment, with open
	compatibilities	with	actual
	applications: ESRI, AutoDesk, Word,
	Excel, PDF, PS, PRT;
	-open property database for with
	maintaining data property ;
	-multiple file management, where the
	integrated management of the files

Following the construction and tuning efforts, quality assurance testers should test and measure an application's performance prior to deployment to establish acceptable quality benchmarks. If all of these efforts are performed well, consequently you are able to better diagnose whether the Web site is operating within established operating parameters when reviewing the statistics generated by Web server monitoring and logging programs.

  The use of dynamic linking and brushing as a central organizing technique for data visualization has a strong tradition in exploratory data analysis (EDA), going back to the notion of linked scatterplot brushing (Stuetzle 1987), and various methods for dynamic graphics outlined inCleveland and McGill (1988). In geographical analysis, the concept of "geographic brushing" was introduced by Monmonier (1989) and made operational in the Spider/Regard toolboxes of Haslett, Unwin and associates (Haslett et al. 1990, Unwin 1994). Several modern toolkits for exploratory spatial data analysis (ESDA) also incorporate dynamic linking, and, to a lesser extent, brushing. Some of these rely on interaction with a GIS for the map component, such as the linked frameworks combining XGobi or XploRe with ArcView(Cook et al. 1996(Cook et al. , 1997, Symanzik et al. 2000), the SAGE toolbox, which uses ArcInfo (Wise et al. 2001), and the DynESDA extension for ArcView (Anselin 2000), GeoDa's immediate predecessor.

  polygon boundary coordinates). The current version adheres to ESRI's shape file as the standard for storing spatial information. It contains functionality to read and write such files, as well as to convert ASCII text input files for point coordinates or boundary file coordinates to the shape file format. It uses ESRI's MapObjects LT2 technology for spatial data access, mapping and querying.

		The analytical
	functionality is implemented in a modular fashion,
	as a collection of C++ classes with associated
	methods. In broad terms, the functionality can be
	classified into six categories:
	 spatial data manipulation and utilities: data
	input, output, and conversion
	 data	transformation:	variable
	transformations and creation of new
	variables
	 mapping: choropleth maps, cartogram and
	map animation
	 EDA: statistical graphics  spatial autocorrelation: global and local
	spatial autocorrelation statistics, with
	inference and visualization
	 spatial regression: diagnostics and
	maximum likelihood estimation of linear
	spatial regression models

  Based of our experience in remote database controlling, we have applied procedures for instance table controlling, like: designing, running and killing temporary tables or databases. Such procedure is represented in next figure and is framed in aspect oriented programming technique.

				ng code: Using this approach the
				presentation is calculated at the client
				side. Pages transmitted consist of
				content and scripting code. At
				runtime, the scripting code is executed
				and	compiles	an	appropriate
				presentation of the content by
				generating the required target format
				(e.g., WML or XHTML). This
				approach is not feasible in situations
				where for instance mobile devices set
				restrictions in code size or computing
				power.
				-Pre-generation: It is common practice
				on the WWW to pre-generate static
				Web pages so that they can be
				accessed at prime time [11]. While the
				WAP architecture in general does not
				prevent this approach it does not meet
				two	major	challenges:	firstly,
				browsers for new mobile devices are
				even more heterogeneous than today's
				Web browsers; secondly, most E-
				commerce applications are expected
				to deliver dynamic content with a high
				update rate thus leading to a humble
				lifespan for a single static page.
				First step in this implementation is focused on
				designing static and dynamic associated tables.
				Based of our experience in remote database
				controlling, we have applied procedures for
				instance table controlling, like: designing, running
				and killing temporary tables or databases. Such
				procedure is represented in next figure (figure 3)
				and is framed in aspect oriented programming
	idea	is	a	technique.
	simplification of using Composite
	Capability/Preference	Profiles
	(CC/PP, [13]) as proposed by the
	World Wide Web Consortium (W3C).
	-S
	c
	r
	i
	p
	t
	i

 -------+-----------+-----+----- -------+-----------+-----+----- -------+-----------+-----+----- -------+-----------+-----+----- -------+-----------+-----+----- -------+-----------+-----+----- -------+-----------+-----+-----
V. CONCLUSIONS

The main task of the formulated project was focused in creating the basement for maximum mobility to remotely monitor and control, by allocating rights and programming facilities to the both mobile and desktop computing systems.

The harder part of the project was related to the synchronizing and interfacing the proposed setup with the existent processing implementation. By using aspect oriented programming that is derivate from object and service-based programming, we have facilitated running applications:

-at the engineering level model integration into force simulated environments with component repository building blocks; -for developing tools and interfaces to rapidly configure models into distributed federation from mobile devices; -for defining meta-model with data structures needed to maintain semantic consistency between models; -for exposing interfaces in legacy simulation systems to fully leverage their capabilities; -remote controlling decision making system.


L Alonso
email: alonso@loria.fr

F Cuny
email: cuny@loria.fr

S Petitjean
email: petitjea@loria.fr

J.-C Paul
email: paul@loria.fr

S Lazard
email: lazard@loria.fr

E Wies
email: wies@loria.fr

Loria -Inria Lorraine

The Virtual Mesh: A Geometric Abstraction for Efficiently Computing Radiosity

Keywords: I.3 [Computing Methodologies]: Computer Graphics; I.3.7 [Computer Graphics]: Three-Dimensional Graphics and Realism Algorithms, Performance Radiosity, wavelets, illumination, curves & surfaces, virtual mesh

In this paper, we introduce a general-purpose method for computing radiosity on scenes made of parametric surfaces with arbitrary trimming curves. By contrast with past approaches that require a tessellation of the input surfaces (be it made up of triangles or patches with simple trimming curves) or some form of geometric approximation, our method takes fully advantage of the rich and compact mathematical representation of objects. At its core lies the virtual mesh, an abstraction of the input geometry that allows complex shapes to be illuminated as if they were simple primitives. The virtual mesh is a collection of normalized square domains to which the input surfaces are mapped while preserving their energy properties. Radiosity values are then computed on these supports before being lifted back to the original surfaces.

To demonstrate the power of our method, we describe a high-order wavelet radiosity implementation that uses the virtual mesh. Examples of objects and environments, designed for interactive applications or virtual reality, are presented. They prove that, by exactly integrating curved surfaces in the resolution process, the virtual mesh allows complex scenes to be rendered more quickly, more accurately and much more naturally than with previously known methods.

INTRODUCTION

Since radiosity techniques have interesting properties that make them useful for new applications (simulation of light to give a filmic "look" to computer animations, precomputed environments for virtual reality), a direct rendering of scenes made of complex shapes becomes highly desirable. Scenes created with modeling software typically include parametric surfaces such as NURBS, cylinders, spheres, toroidal patches, arbitrary planar primitives and their CSG combinations. Moreover, the limit surfaces of many popular subdivision schemes (such as Loop or Catmull-Clark) are known to be parametric surfaces -see, e.g., [START_REF] Stam | Exact evaluation of Catmull-Clark subdivision surfaces at arbitrary parameter values[END_REF]]. This paper presents a novel and powerful method for the direct radiosity-based rendering of such surfaces.

One assumption of the original radiosity algorithm is that every object of a scene is made of a collection of planar patches [START_REF] Goral | Modelling the interaction of light between diffuse surfaces[END_REF]. Unfortunately, few of the Fig. 1. This geometry room has been illuminated with a high-order wavelet radiosity solution using the virtual mesh. The curved objects have been modeled with the SGDL modeling kernel (http://www.sgdl.com).

later improvements have entirely freed themselves of this hypothesis. For instance, in hierarchical radiosity [START_REF] Hanrahan | A rapid hierarchical radiosity algorithm[END_REF]], lower-level elements are assumed to be planar polygons. Galerkin radiosity [START_REF] Heckbert | Finite element methods for global illumination[END_REF][START_REF] Zatz | Galerkin radiosity: a higher-order solution method for global illumination[END_REF]] is an exception, but its hierarchical versions have been essentially limited to planar primitives (including wavelet radiosity [START_REF] Gortler | Wavelet radiosity[END_REF][START_REF] Schröder | Wavelet algorithms for illumination computations[END_REF]). Thus, with all known global illumination methods, incorporating curved patches is always achieved at the cost of some geometric approximation or important restriction on the shapes allowed: either the input surfaces are tessellated in a number of planar elements until a certain flatness is achieved (which can result in an enormous increase in the number of initial surfaces) or their geometry is approximated for the evaluation of energy transfers (as in the original progressive radiosity solutions -see [START_REF] Arquès | A new radiosity approach for regular objects: application to ruled surfaces[END_REF] and the references therein). In both cases, the integration of trimmed curved patches is not a natural extension of the polygonal case.

This paper presents the virtual mesh, a novel way of computing radiosity on curved surfaces with arbitrary trimming curves that combines the power of adaptive subdivision and functional analysis while avoiding the tessellation of input surfaces. We effectively bridge the gap between general parametric patches and a hierarchical Galerkin formulation with efficient and compact function bases. Furthermore, since our method does not resort to approximations of curved patches, it eliminates almost all geometric sources of error in the algorithm and thus avoids the resolution of a perturbated balance equation [START_REF] Arvo | A framework for the analysis of error in global illumination algorithms[END_REF]]. The virtual mesh is thus a powerful solution to the problem of incorporating curved surfaces in radiosity algorithms. With it, scenes made of complex shapes can be illuminated both quickly and very accurately, as Figure 1 illustrates.

The essence of radiosity in a functional setting implies that calculations must be performed over normalized domains (i.e., unit squares or unit right triangles). This is especially important for constructing the basis of functions, defining optimal quadrature rules and maintaining coherence between all the levels in the hierarchy. Accordingly, our virtual mesh method uses an abstract representation, made of normalized virtual supports, of the input scene to perform radiosity calculations which would otherwise be intractable. This really is a three-step process (see Figure 2):

1. map the input surfaces (curved patches and complex polygons1 ) to normalized planar domains, 2. extend the radiosity function to the virtual mesh, the collection of these supports, and compute radiosity using the virtual mesh as the new input geometry, 3. read off the intensity values for the original surfaces from the values computed on the virtual mesh.

A key challenge of step 1 is to find appropriate mapping functions, i.e., functions that preserve energy properties of the graphics primitives. Given an input surface, such mappings, called energy-preserving, leave the ratio of the areas of two surface elements invariant 2 . The use of energy-preserving mappings allows hierarchical Galerkin methods with high-order basis functions to be used without having to recompute many new attributes at each level of the hierarchy. This invariance property has a simple mathematical expression when the input surfaces can be globally parameterized, so we focus on trimmed parametric surfaces in the rest of the paper. This is only a mild restriction since parametric objects are those that can be most easily displayed in computer graphics.

The virtual mesh can be plugged into any advanced variation of radiosity algorithms, including discontinuity meshing, clustering and hierarchical subdivision, and in a variety of finite element methods. To demonstrate its efficiency, we have developed and implemented a high-order wavelet radiosity algorithm that uses our new tool. Since the cost of tessellation increases with the order of the wavelet basis, the virtual mesh is especially applicable to higher-order wavelets. Our implementation can illuminate scenes made of two major types of primitives: complex polygons (and in particular ordinary polygons, the most common primitives used by present-day modelers) and low-degree trimmed algebraic patches, a primitive set known to encompass 95 percent of conventional, unsculptured parts in industrial environments [START_REF] Requicha | Solid modeling: a historical summary and contemporary assessment[END_REF]. It can easily be extended to higher-order surfaces as soon as the appropriate energy-preserving mappings can be numerically evaluated.

The paper is organized as follows. Section 2 describes related work on radiosity with curved primitives. The next two sections give the high-level ideas behind the virtual mesh in radiosity-based global illumination, independently of any resolution method: Section 3 shows how to build the virtual mesh and Section 4 how to use it in a radiosity algorithm. Section 5 then describes how we have incorporated the virtual mesh in a hierarchical radiosity implementation which uses high-order wavelet functions. Section 6 demonstrates the efficiency of the proposed technique on several test cases, before concluding.

RELATED WORK

In radiosity algorithms [START_REF] Cohen | Radiosity and Realistic Image Synthesis[END_REF][START_REF] Sillion | A clustering algorithm for radiosity in complex environments[END_REF], one is led to work over very simple supports. This is made necessary because representing functions over general domains is not straightforward. The solution presented in the literature has been to use the tensor product of two 1D bases [START_REF] Gortler | Wavelet radiosity[END_REF][START_REF] Zatz | Galerkin radiosity: a higher-order solution method for global illumination[END_REF]]. Often, orthonormal polynomials like the Legendre or Jacobi polynomials are used as 1D bases. Unfortunately, this construction may only be carried out if the support is a simple polygon, i.e., a triangle or a parallelogram. For more general supports, there does not seem to be a simple way of finding 2D basis functions.

In the course of the resolution, one is also led to compute integrals on geometric supports as well as the mean value of the residual. Usually, such computations can only be carried out with numerical methods. Such methods (known as quadrature rules) involve sampling the kernel of the radiosity equation at various points and approximating the integrals as weighted sums of kernel samples. For general supports, there is no automatic way of finding good quadrature points that ensure reconstruction up to a desired accuracy. Also, to avoid recalculation of quadrature points and weights for different types of simple polygons, it is desirable to map the input geometry to normalized domains.

For all these reasons, radiosity-based rendering has largely been limited until now to scenes made of simple polygons. Exceptions to this general rule rely, to illuminate curved surfaces, on some kind of tessellation or geometric approximation [Schäfer • 5 1997;[START_REF] Stamminger | Bounded radiosity -Illumination on general surfaces and clusters[END_REF] or introduce errors in the resolution by working in the parametric domain without preserving areas [START_REF] Christensen | Global illumination of glossy environments using wavelets and importance[END_REF].

To incorporate curved objects into a hierarchical framework, Schäfer [START_REF] Schäfer | Hierarchical radiosity on curved surfaces[END_REF]] propose a three-dimensional extension of adaptive meshing, called objectbased meshing, which avoids the initial subdivision of curved parts. First, each curved object is wrapped in a very crude bounding box. Then, local improvements of the polygonal approximation can be computed at any stage of the algorithm by generating new vertices on the faces of the bounding box and projecting them on the surface of the original object. Since patches at a given level of the hierarchy are no longer the exact unions of child patches and have a different orientation, push-pull operations have to be modified. Also, the influence of geometric approximations on the result may be difficult to evaluate. [START_REF] Stamminger | Bounded radiosity -Illumination on general surfaces and clusters[END_REF]] take a different approach with their bounded radiosity method. Their algorithm takes into account the full extent of the interacting objects without resorting to point sampling. The refinement criterion, however, is only a function of upper and lower bounds on form factors, computed using an estimate of the range of normals and boxes bounding the geometry. This may produce unnecessary subdivisions, inducing a lack of efficiency of the hierarchical operations. Note also that this method (as well as the one of Schäfer) assumes the lower-level elements of complex surfaces to be polygons, so it does not work on the exact geometry of the input surfaces.

Christensen et al. [START_REF] Christensen | Global illumination of glossy environments using wavelets and importance[END_REF]] present an importance-driven algorithm for the global illumination of glossy environments using piecewise-constant wavelets. A four-dimensional wavelet representation is used to represent the spatially-and angularly-varying radiance distributions across surfaces in the environment. Support for curved surfaces, a feature of the implementation, is limited to surfaces which can be parameterized on the unit square. Handling of parametric surfaces with general trimming curves is thus unclear. Also, only convex patches are allowed, since interactions of an object with itself are not taken into account. In addition, the spatial representation of wavelets on the unit square uses the natural parameters of the patch and so may be heavily biased. Indeed, a large part of the parameter domain may well represent a small part of the patch area, depending on the parameterization of the patch. Since the authors use a uniform sampling of this parameter domain to approximate the transport coefficients, there is no simple way to control the density of sample points needed to sufficiently cover the curved patch and to bound the approximation error. Similar comments apply to the refinement oracle and to the uniform subdivision of the parameter domain, which may corrupt the radiosity computations. In other words, the physical meaning of the quantities computed is lost.

The classical hierarchical radiosity algorithms have a complexity of O(k 2 + n), where k is the number of initial patches and n the number of final patches. Since the O(k 2 ) term has a strong influence on the resolution process, several algorithms have sought to eliminate it by grouping input surfaces into clusters [Smits et al. 1994;[START_REF] Sillion | A unified hierarchical algorithm for global illumination with scattering volumes and object clusters[END_REF][START_REF] Christensen | Clustering for glossy global illumination[END_REF][START_REF] Willmott | Face cluster radiosity[END_REF] or simplifying the geometry [START_REF] Rushmeier | Geometric simplification for indirect illumination calculations[END_REF]. Though they could be used to illuminate curved patches, clustering algorithms have only been used in polygonal scenes up to now. And as indicated by a recent study [START_REF] Hasenfratz | A practical analysis of clustering strategies for hierarchical radiosity[END_REF], applying clustering algorithms to tessellated curved surfaces without preserving topology would be largely inefficient.

To sum up, the best known approaches for globally illuminating curved surfaces suffer from the following shortcomings:

• the traditional meshing approach applied to curved surfaces may lead to severe artifacts and continuity problems. For instance, it can generate a large number of triangles or slivers which are difficult to handle in visibility computations and can cause Z-buffer artifacts when the radiosity solution is visualized; • all "smart" (i.e., other than straight tessellation) hierarchical radiosity methods for illuminating curved patches introduce errors of a geometrical or physical kind which may deeply affect the accuracy of the radiosity solution; • since the hierarchical subdivision of curved parts is not the natural extension of the subdivision of flat regions, either geometric values (area, range of normals) have to be recomputed and stored for each patch of the hierarchy [START_REF] Schäfer | Hierarchical radiosity on curved surfaces[END_REF]], which slows down the simulation, or refinement may be very inefficient (depending on the parameterizations used) because the subdivision of the parameter domain does not correspond to a quadtree of the surface [START_REF] Christensen | Global illumination of glossy environments using wavelets and importance[END_REF]]; • all past methods are limited in the range of shapes that they are capable of illuminating.

Our contributions

By contrast, the virtual mesh method has many advantages over these previous attempts at illuminating general parametric surfaces:

• genericity -It is a natural and strong solution to the rendering of curved objects in a radiosity setting and it provides a currently missing link between hierarchical Galerkin radiosity and parametric surfaces (both convex and concave) with arbitrary trimming curves. • accuracy -By keeping a virtual but exact image of the input geometry throughout the resolution process, the virtual mesh avoids most of the geometric errors that plague radiosity algorithms when applied to complex shapes. • flexibility -Since the virtual mesh is a "normalized" representation of the input geometry with a set of bijective mappings that allow to go back and forth between real and virtual supports, our method can be combined with any advanced variation of the basic radiosity method.

To illustrate this flexibility, we present an implementation of the virtual mesh in a hierarchical Galerkin radiosity algorithm. Our radiosity program takes advantage of the very good compression ratios provided by high-order wavelets for storing radiosity. Recent work [START_REF] Cuny | A novel approach makes higher order wavelets really efficient for radiosity[END_REF]] has shown that high-order wavelets are more efficient and more powerful than piecewise-constant wavelets for the illumination of real-world scenes (using a shooting algorithm and getting rid of links [START_REF] Stamminger | Getting rid of links in hierarchical radiosity[END_REF][START_REF] Winkler | Expérimentation d'algorithmes de calcul de radiosité à base d'ondelettes[END_REF]]), disproving a common belief.

The virtual mesh can also be combined with a clustering strategy, to allow for a more efficient manipulation of complex shapes inside clusters and accelerate the computations for the local energy transfers between close surfaces. Motivated by applications where physical accuracy is critical (as in lighting design, where

• 7

an exact measure of energy prevails), we have decided not to experiment with clustering in this paper. Indeed, error control is very difficult with clustering and its practical use may result in unexpected behaviors, as recently emphasized in [START_REF] Hasenfratz | A practical analysis of clustering strategies for hierarchical radiosity[END_REF]].

• efficiency -Since it gives radiosity algorithms the ability to work with compact representations of real-world scenes, the virtual mesh induces important benefits in terms of computation time, memory consumption and visual quality, compared to the illumination of tessellated objects. In turn, these improvements allow for high-accuracy wavelet radiosity solutions to be computed on large scenes with a reasonable time and memory consumption.

BUILDING THE VIRTUAL MESH

As we have illustrated, the main reason why, in the past, radiosity-based rendering has been performed essentially on scenes made of simple polygons is because of the repeated need to do computations on simple normalized domains. Rather than prematurely approximate complex surfaces via triangulation, we use a high-level representation of the scene, the virtual mesh. We show that this representation allows for general planar primitives and trimmed parametric surfaces to be handled as naturally and efficiently as traditional primitives. This section and the next are dedicated to the principles of the proposed technique, independently of any particular resolution algorithm of the radiosity equation. We start here by giving a thorough presentation of the construction of the virtual mesh for parametric patches. Section 4 will then show how this representation can be used to efficiently and accurately compute radiosity on trimmed curved patches and complex polygons.

The main idea of the virtual mesh is to compute radiosity on a set of unit virtual supports representing the (possibly curved) input surfaces. For this to have a physical meaning, the mappings involved must be energy-preserving. (As we shall see, this invariance property is crucial when using non-constant basis functions in a hierarchical Galerkin setting -Section 5 -and allows for important speedups -Section 6). In mathematical terms, since radiosity is a measure of power by unit area, a mapping is energy-preserving if it leaves the ratio of the areas of two surface elements invariant, i.e., if its Jacobian3 is constant. Now, a direct mapping of a curved or complex polygonal primitive S onto a unit square is not always possible, so we separate the problem at hand into three different tasks (see Figure 3):

(1) find a constant Jacobian mapping Ψ -1 that maps a curved patch S to the plane (Section 3.1),

(2) enclose the image S of the patch with a parallelogram S (Section 3.2),

(3) map the parallelogram S onto a unit square S (Section 3.3).

Obviously, if S is already planar, the first step is skipped.

S S' S S' S Λ Ψ -1 -1
(2)

(3)

(1) Fig. 3. Mapping a trimmed curved patch to a unit square domain.

Mapping curved patches to the plane

Finding a mapping with a constant Jacobian which "flattens" a curved surface is a key issue in our method. It allows calculations to be performed on a simpler geometry while preserving energy properties of the input patch. This section examines this issue.

3.1.1 Intuitively. Let S be a smooth curved patch parameterized by

x(θ, ϕ) = (x(θ, ϕ), y(θ, ϕ), z(θ, ϕ)), (θ, ϕ) ∈ I, I ⊆ R 2 .
θ and ϕ are called the natural parameters of S.

In what follows, subscripts indicate partial derivatives. The isoparameter curves

C θ : x(θ = cst, ϕ), C ϕ : x(θ, ϕ = cst),
where cst means constant, define a net of curves on S. At each parameter pair (θ 0 , ϕ 0 ) ∈ I, the vector x θ (θ 0 , ϕ 0 ) is tangent to C θ0 at x 0 = x(θ 0 , ϕ 0 ) and the speed of the position vector along C θ0 at x 0 is x θ (θ 0 , ϕ 0 ) . Similar observations hold for C ϕ0 . The quantity x θ (θ 0 , ϕ 0 ) × x ϕ (θ 0 , ϕ 0 ) dθ dϕ measures the area of a differential element around point x 0 . The image of this differential element by the "natural coordinates" mapping

χ -1 : (x, y, z) -→ (θ, ϕ) is a rectangle of area dθ dϕ. Thus χ is energy-preserving if x θ (θ 0 , ϕ 0 ) × x ϕ (θ 0 , ϕ 0 ) = cst, ∀(θ 0 , ϕ 0 ) ∈ I.
There is no reason for this to be the case. Thus, to find an energy-preserving mapping, we need to reparameterize S by the parameters (u, v) in such a way that

x u (θ 0 , ϕ 0 ) × x v (θ 0 , ϕ 0 ) = cst, ∀(θ 0 , ϕ 0 ) ∈ I.
(1)

It is always possible to construct an energy-preserving mapping Ψ : (u, v) → (x, y, z) by taking for instance u = θ and adjusting the speed of the position vector along C ϕ in such a way that the constraint of Equation ( 1) is satisfied at each point of S.

3.1.2 Formally. What the above shows is that there exists at least one energypreserving mapping that sends S to a suitable parameter domain. Now, there are

• 9 actually many qualitatively different energy-preserving mappings for a given curved surface S and some may be easier to deal with computationally than others.

The general case is to consider θ = f (u, v) and ϕ = g(u, v). Actually, we have observed that it is usually not necessary to have both u and v depend on θ and ϕ. So assume for now that θ = f (u), ϕ = g(u, v). Let us determine the constraint that f and g must satisfy for Ψ to be energy-preserving. Since

x u (θ, ϕ) = f (u)x θ (θ, ϕ) + g u (u, v)x ϕ (θ, ϕ), x v (θ, ϕ) = g v (u, v)x ϕ (θ, ϕ),
the Jacobian of Ψ is [START_REF] Kaplan | Advanced Calculus[END_REF]]:

J Ψ (u, v) = x u (θ, ϕ) × x v (θ, ϕ) = |f (u)g v (u, v)| x θ (θ, ϕ) × x ϕ (θ, ϕ) .
Now, Ψ and Ψ -1 are energy-preserving if the Jacobian J Ψ is constant, say equal to 1/c, c = 0. In virtue of the theorem on the derivatives of inverse functions, we get (up to sign)

∂v ∂ϕ (θ, ϕ) = 1 g v (u, v) = c|f (u)| x θ (θ, ϕ) × x ϕ (θ, ϕ) .
Integrating with respect to ϕ and replacing θ and ϕ by f (u) and g(u, v) respectively leads to

v = c|f (u)| g(u,v) 0 x θ (f (u), t) × x ϕ (f (u), t) dt. (2)
Any choice of f and g satisfying the constraint of Equation ( 2) gives an energypreserving mapping Ψ -1 : (x, y, z) → (u, v) and its inverse.

For the different operations we need to carry out (computing the virtual mesh from the input geometry, forwarding visibility queries to the virtual mesh, reconstructing the radiosity function on curved patches), we want to be able to go back and forth between (x, y, z), (θ, ϕ) and (u, v). Thus, we need to know the inverse of all the different mappings involved. This is not necessarily as difficult as it may seem. Indeed, all these mappings are bijective (except possibly on a set of measure zero) so that, in case the inversion cannot be achieved exactly, a numerical algorithm will do a fine job with arbitrary precision. Also, the parameterizations we shall consider in this paper all involve circular or hyperbolic functions, for which the inverse mapping (x, y, z) → (θ, ϕ) can be computed exactly.

3.1.3

Example: a ring torus. Consider for instance the doughnut-like shape of Figure 4.a. Mathematically, it corresponds to a patch on the ring torus given by the parameterization

   x(θ, ϕ) = cos θ (cos ϕ + 2), y(θ, ϕ) = sin θ (cos ϕ + 2), z(θ, ϕ) = sin ϕ, θ ∈ 0, π 2 , ϕ ∈ [0, π).
Equation ( 2) rewrites as: Clearly, choosing f to be a linear function of u will get rid of the |f (u)| term. Let us for instance choose f (u) = u/2 and c = 1. Equation ( 3) boils down to a constraint on g for Ψ : (u, v) → (x, y, z) to be energy-preserving:

v = c|f (u)| g(u, v) + 1 2 sin g(u, v) . (3)
v = g(u, v) + 1 2 sin (g(u, v)).
Here, it turns out that g does not depend on u, so we may write g(u, v) = g(v).

With the above equation, g can be determined numerically, from which we can deduce the energy-preserving mapping Ψ:

Ψ : (u, v) -→ cos u 2 (cos g(v) + 2), sin u 2 (cos g(v) + 2), sin g(v) . Since v = g(v) + 1 2 sin g(v) = ϕ + 1 2 sin ϕ, the inverse mapping Ψ -1 is: Ψ -1 : (x, y, z) -→ 2θ, ϕ + 1 2 sin ϕ .
The result of applying Ψ -1 to our ring torus patch gives the planar surface displayed in Figure 4.b. This should be compared with Figure 4.c which shows the result of applying the mapping (x, y, z) → (2θ, ϕ), the Jacobian of which is non-constant, which sends a 3D point to its natural coordinates (up to a constant factor).

Appendix A gives additional examples of energy-preserving mappings for a few curved surfaces of interest.

Enclosing complex polygonal primitives

What the previous section has essentially proved is that a curved input patch S can be replaced with a planar primitive S . Now, this primitive may be somewhat arbitrary (as in Figure 3): it may have curved boundaries, holes. . . In general, it is not possible to map a curved patch directly onto a unit square in such a way that the mapping is bijective, except in important but special cases (an example of which is the torus patch of Figure 4). An intermediate step is to bound the planar primitive S by a parallelogram S . We do it in such a way as to minimize the voids (i.e., the area of S \ S ). This way, we also minimize the regions where the radiosity function has to be extended (see Section 4.1) and thus the amount of unnecessary calculations.

How does one get a good bounding box? If S has a polygonal outer boundary with n vertices, we start by computing its convex hull in O(n). Then we use the • 11 minimal enclosing parallelogram algorithm of Schwarz et al. [START_REF] Schwarz | Minimal enclosing parallelogram with application[END_REF] which has a complexity linear in the number of vertices. The resulting enclosing parallelogram S indeed has the property that it minimizes the voids. If S has a curved boundary, we simply discretize this boundary, compute a minimal enclosing parallelogram of the discrete curve with the algorithm of Schwarz et al. and use this parallelogram as the bounding box of S . Note that the overhead of computing an enclosing parallelogram for each input patch -something that has to be done only once at the beginning of the algorithm -is negligible compared to the overall computation time (i.e., less than one thousandth of the total).

Mapping quadrilaterals to unit squares

The last step of the chain is to map the parallelogram S found previously onto a unit square S. This is an easy task, which can be done using an affine mapping Λ -1 whose inverse can be written as follows:

(u , v ) = Λ(ũ, ṽ) = A + ũ--→ AB + ṽ--→ AD,
where the points A, B, C, D are the vertices of the quadrilateral. As can readily be observed, the Jacobian of this mapping is constant: it is equal to the area of S .

USING THE VIRTUAL MESH

Now that we have shown how to build the virtual mesh, our aim in this section is to discuss its use in a radiosity algorithm. The main issue here is that, generally speaking, not all points of the unit virtual support S correspond to points of the original S. Thus, if we want to compute radiosity on S, we need to extend the definition of the radiosity function on S to S (or rather to S, the 3D equivalent of S). Section 4.1 presents our simple extension of the radiosity function to the virtual mesh. Then, Section 4.2 shows how to compute radiosity on the virtual mesh, knowing that it represents a geometry that is larger than the original geometry. Once radiosity values have been computed on the virtual mesh up to the desired convergence rate, they can be lifted to the original mesh using the inverse mappings S → S → S and subsequently passed to the rendering engine as in traditional radiosity algorithms.

Defining radiosity on S

In general, there are "more" points on S than on S. We examine here how radiosity can be defined on those extra points.

Let S be the pre-image of S by (Ψ • Λ) -1 (see Figure 3). Call a point of S real if it has a pre-image on S, virtual otherwise (i.e., if belongs to S \ S). The extension of the radiosity function to virtual supports should follow two basic rules:

• The extended function is the same as the original one on real points of S.

• The extended function should be as simple and continuous as possible on virtual points of S and on the boundaries between real and virtual points.

The first rule is essential since we want to deduce radiosity values on S from radiosity values on S. The second rule is no less important. Indeed, in hierarchical radiosity (Section 5), the radiosity value at any point of S (e.g., at a virtual point) can influence the radiosity value at any other point of S (e.g., at a real point). Thus, the radiosity function must be computed with the same precision at all points of S.

Since discontinuities of the derivatives of the radiosity function are only due to discontinuities of the visibility function [START_REF] Heckbert | Discontinuity meshing for radiosity[END_REF]], we do not want our extension to introduce new visibility discontinuities. To this aim, we introduce an extension ṽ of the visibility function v as follows. Visibility between a point p of the environment and a point q of S is defined as the visibility between p and q, where q is the point of S that is closest to q: ṽ(p, q) = v(p, q), qq = min q ∈S qq .

The radiosity function over S is then defined in a classical way, but using ṽ as visibility function (see Figure 5):

B(x) = E(x) + ρ(x)
S B(y)G(x, y)ṽ(x, y) dA.

Even though this may not be the best extension of the radiosity function, it is very simple and has proved to be very reliable in practice.

Computing radiosity on S

Now that we have properly defined the radiosity function at all points of the virtual mesh, we are almost ready to compute the extended radiosity function and we must do so in a meaningful way. The general idea is to work with our high-level representation whenever possible and go back to the original geometry only when needed (i.e., for visibility queries). The main issue here is that since we work with a virtual geometry that represents more points than the original geometry, we must be careful not to propagate the wrong quantity of energy.

Let us analyze the three different cases, i.e., when S is a receiver, a blocker or an emitter.

S is a receiver. When S acts as a receiver, we use S exactly as if this virtual support was the receiver. Real and virtual points are handled in exactly the same way and with the same accuracy. This way, the radiosity function is consistent over the whole virtual support and can be used in a hierarchical algorithm without the virtual points disturbing the real points of S. S is a blocker. When S acts as a blocker, we must use its real geometry, and not S, in visibility calculations. We proceed as follows.

Let Ṡ be the untrimmed equivalent of S (i.e., if S is defined by x(θ, ϕ), (θ, ϕ) ∈ I, the surface Ṡ is defined by x(θ, ϕ), (θ, ϕ) ∈ R 2 ). First, we compute the intersection of the visibility ray with Ṡ:

• If S is curved, computing the intersection is not too hard since the surface is given in parametric form. If there is an intersection, we project it using (Ψ • Λ) -1 . If the projection falls outside of the virtual support S, we are done since this means that the 3D intersection point is not on S. • If S is planar, then Ṡ is a plane so computing the intersection point (there is almost always one) is easy. Project this point using Λ -1 . As before, if the projection falls outside of the virtual support S, we are done.

To determine if the projected point is a real or a virtual point of S (i.e., if S is a blocker for the ray considered or not), we use a trapezoidal map algorithm [Boissonnat and Yvinec 1998] which we have slightly improved to take curved arcs into account. The rough idea is to project the trimming curves of S (the projected curves come in two flavors: line segments and curved arcs) and to divide the map area (here: the virtual support S) into a number of "slabs" by drawing vertical lines through the ends of each line segment and through the inflections and extrema in the u and v directions of each curved arc. It is then possible with a fast binary search to find which slab the projected point is in and consequently whether S is a blocker for the current visibility ray (see the example of Figure 6).

S is an emitter. When S acts as an emitter, we need to be careful because both real and virtual points of its virtual support have received energy. In other words, since only the real points have a physical meaning, S has accumulated more energy than it should have. To avoid repropagating this extra energy, the idea is to continue working with the extended radiosity function on S but to use as weights of quadrature points only a fraction of what they are over real regions, as in Figure 7. (Adjusting quadrature weights to the virtual support is similar in spirit to the shadow masks used by Zatz to smooth the shadow discontinuities out of the radiosity distribution seen by the Galerkin method [START_REF] Zatz | Galerkin radiosity: a higher-order solution method for global illumination[END_REF]].)

To compute the interaction between two basis functions (one on the emitter and one on the receiver) using a quadrature rule, we are lead to compute the interaction between a function on S and a point on the receiver. Assume for now that the whole of S shoots energy. We need to evaluate integrals of the form emitter f (p) dA, where dA is an element of area on S, which are approximated as weighted sums

p i ∈ emitter w i f (p i ),
where w i is the weight associated to quadrature point p i . Since we work with constant Jacobian mappings, a possibility is to place a regular n × n grid on S, take the quadrature points at the center of each cell and compute

S f (p) dA ≈ n i=1 n j=1 w ij f (p i,j ),
where p i,j = 2i-1 2n , 2j-1 2n and w i,j = 1 n 2 is the area of the cell in which p i,j falls. Now, a simple way to adjust this quadrature rule to the real part of the virtual support and to avoid energy overpropagation is to take as quadrature points pi,j = p i,j and to use as quadrature weight of pi,j the area of the "real" part of the cell in which it falls, i.e. the part that corresponds to real points of S (as in Figure 7.b). Clearly, as the size n of the grid is increased, the total area of the mixed cells (i.e., cells having both real and virtual points) will decrease as a function of 1 n 2 , so that the right amount of energy will be emitted.

VIRTUAL MESH IN HIGH-ORDER WAVELET RADIOSITY

Since the virtual mesh uncouples the input geometry from the representation of radiosity, the general method presented in the previous section can be efficiently combined with any advanced radiosity algorithm. In particular, it fits nicely in a hierarchical Galerkin radiosity algorithm. We report here an implementation of a high-order wavelet radiosity algorithm with the virtual mesh.

Algorithm overview

Our hierarchical radiosity algorithm paired with the virtual mesh is outlined in Algorithm 1. It is based on the well-known shooting method with Southwell relaxation. As can readily be seen, this algorithm has a lot in common with a classical hierarchical radiosity algorithm, apart from some steps which we shall detail in the following sections.

The interested reader is referred to [START_REF] Cuny | A novel approach makes higher order wavelets really efficient for radiosity[END_REF] for the details of our implementation. Suffice it to say for the purpose of the present paper that the refinement oracle is similar to the one of [START_REF] Bekaert | Error control for radiosity[END_REF]

(estimate

The Virtual Mesh

• 15

Algorithm 1 Overview of HR with the virtual mesh.

1: for all i ∈ {1, . . . , n} do 2:

Si ←-VirtualSupport(S i )

3:

Bi ←-Ẽi , ∆ Bi ←-Ẽi 4: end for 5: while k S k ,λ ∆ Bk > ε do 6:

Choose i such that Si has max unshot energy ∆ Bi ∆ Bi ←-0 13: end while 14: for all i ∈ {1, . . . , n} do 15:

UpdateHierarchy( Si )

16:

Draw (S i ) 17: end for errors on the radiosity function using control points on the receiver) and that we use shooting rather than the more conventional gathering approach because it is what allows to get rid of links (as explained in [START_REF] Stamminger | Getting rid of links in hierarchical radiosity[END_REF]] and [START_REF] Winkler | Expérimentation d'algorithmes de calcul de radiosité à base d'ondelettes[END_REF]]) and also because it is what makes higher-order wavelets good and efficient for radiosity.

We use the most general wavelet basis functions known in the radiosity field, i.e. high-order wavelets. We have observed that for linkless radiosity, and with the refinement oracle described in [START_REF] Cuny | A novel approach makes higher order wavelets really efficient for radiosity[END_REF]], a previously noted annoying behavior of high-order tree-like wavelets (namely, they can produce ringing discontinuities when the refinement oracle is not high enough [START_REF] Willmott | An empirical comparison of progressive and wavelet radiosity[END_REF]) is rarely a problem.

When working over curved surfaces, which may exhibit large variations of curvatures and normals, it is desirable to increase the order of the wavelet functions to better capture the details of the radiosity function and reduce the number of final patches. For instance, we have determined that when the input scene is made of quadric surface patches, the best visual results are obtained using linear or quadratic wavelets. And as shown in [START_REF] Cuny | A novel approach makes higher order wavelets really efficient for radiosity[END_REF]], using linear (M 24 ) or quadratic (M 3 ) wavelets in a radiosity algorithm is actually faster and less memory-consuming than using Haar (M 1 ) wavelets, for the specific case of a shooting algorithm that does not store links.

As said in Section 4, we use the original geometry of the surface only when needed, i.e., for visibility queries. This means that the hierarchy of subdivisions of the surface is constructed over the virtual mesh and not over the original surface. In Section 5.2, we show that the constant Jacobian property is crucial for the push-pull coefficients computed on the virtual mesh to be meaningful for the original surface and provides a very natural answer to the extension of wavelet radiosity to curved surfaces and complex polygons.

Two points will then remain to be discussed, which correspond to steps underlined in red in Algorithm 1. Section 5.3 shows how to handle self-shooting surfaces (line 8 of the algorithm, a surface can interact with itself) and Section 5.4 details the management of the hierarchy (line 9).

Push-pull coefficients

In this section, we prove that the use of energy-preserving mappings greatly simplifies the computation of push-pull coefficients over general supports and allows for a key property of classical push-pull coefficients (namely, their independence on the level of the hierarchy) to be transfered to general parametric surfaces.

To study this property, recall that 1D wavelet functions are generated from a single function φ by dilations and translations [START_REF] Schröder | Wavelet algorithms for illumination computations[END_REF]]:

φ m i (x) = √ 2 m φ(2 m x -i).
Let S m i be the support of φ m i . Let f be the projection of a univariate function f defined on [0, 1] in a tree wavelet basis. Then we have:

f m (x) = i α m i φ m i (x),
where m is the level of representation and the α m i are weights associated to φ m i . The hierarchical properties of wavelets allow information to be transfered from one level of the hierarchy to the next or the previous by computing a set of push-pull coefficients h m i,j and g m i,j :

α m+1 j = i g m i,j α m i , α m-1 i = j h m i,j α m j . (4)
The classical case. Suppose given f m and α m i . Since level m + 1 is finer than level m, f m can also be represented at level m + 1. In other words, there exist coefficients α m+1 i such that

α m i φ m i (x) = α m+1 i φ m+1 i (x).
Multiplying both sides by φ m+1 j (x) and taking the inner product on the space of functions yields:

α m i φ m i , φ m+1 j S m+1 j = α m+1 i φ m+1 i , φ m+1 j S m+1 j . ( 5
)
The classical wavelet functions form an orthonormal basis. Thus:

φ m+1 i , φ m+1 j S m+1 j = δ ij ,
where δ ij is the Kronecker delta. We end up with

α m+1 j = α m i φ m i , φ m+1 j S m+1 j .
The Virtual Mesh

• 17

With this and another similar calculation, we conclude that:

g m i,j = φ m i , φ m+1 j S m+1 j , h m i,j = φ m i , φ m+1 j S m i .
Since classical wavelet functions are defined on normalized domains, one can prove that the push-pull coefficients do not depend on the level of representation m. For instance, for h m i,j :

h m i,j = S m i φ m i (x)φ m+1 j (x) dx = 2 m √ 2 S m i φ(2 m x -i)φ(2 m+1 x -j) dx, = √ 2 [0,1] φ(y)φ(2y + 2i -j) dy, (y = 2 m x -i), = φ, φ 1 2i-j [0,1] = h 0 0,2i-j .
Similar properties can be proved for 2D wavelet functions. Suppose that these functions are defined as tensor products of 1D wavelet functions (other definitions are possible, but lead to the same conclusions):

φ m i,j (x, y) = φ m i (x)φ m j (y).
A scalar function f defined over [0, 1] 2 (e.g., the radiosity function) can be projected in this basis of functions:

f (x, y) = i,j
α m i,j φ m i,j (x, y).

As before, a set of push-pull coefficients can be defined which allow for information to be transfered between levels of the hierarchy:

h m (i,j),(k,l) = φ m i,j , φ m+1 k,l S m i,j , g m (i,j),(k,l) = φ m i,j , φ m+1 k,l S m+1 k,l
.

And as before, these coefficients do not depend on m. For instance:

h m (i,j),(k,l) = S m i,j φ m i,j (x, y)φ m+1 k,l (x, y) dx dy, = S m i φ m i (x)φ m+1 k (x) dx S m j φ m j (x)φ m+1 l (x) dy , = h m i,k h m j,l = h 0 0,2i-k h 0 0,2j-l = h 0 (0,0),(2i-k,2j-l) .
The general case. The above properties are very desirable to avoid heavy and costly calculations each time a hierarchy is created. We now show that to have similar properties in the general case, we need to work with constant Jacobian mappings.

When the support is not normalized (say, in 1D, a curved arc), the basis functions φm i are not orthonormal and they are different from the classical wavelet functions. In this case, combining Equations ( 4) and ( 5) gives a definition of the push-pull coefficients in matrix form:

H m = (A m ) -1 B m , G m = (C m ) -1 D m ,
where

H m = ĥm i,j , G = ĝm i,j , A m = φm i , φm j Ŝm i , C m = φm i , φm j Ŝm j , B m = φm i , φm+1 j Ŝm i , D m = φm i , φm+1 j Ŝm+1 j .
We have a similar definition in the 2D case which involves coefficients of the form:

φm i,j , φm k,l Ŝm i,j , φm i,j , φm k,l Ŝm k,l , φm i,j , φm+1 k,l Ŝm i,j , φm i,j , φm+1 k,l Ŝm+1 k,l
.

Let S be a parametric surface. Recall the notations of Section 3:

S Ψ -1 ---→ S , S Λ -1 ---→ S,
where Ψ -1 : (x, y, z) → (u, v) sends S to the parametric domain (u, v), S is a parallelogram enclosing the image S of S by Ψ -1 and Λ -1 : (u , v ) → (ũ, ṽ) sends S to [0, 1] 2 . Let Ŝm i,j = (Ψ • Λ)(S m i,j ). We aim to see how push-pull coefficients on S, or more precisely on S = (Ψ • Λ)( S), relate to those of S. We start by computing the following coefficient:

φm i,j , φm+1 k,l Ŝm i,j = Ŝm i,j φm i,j (θ, ϕ) φm+1 k,l (θ, ϕ) dA, = Šm i,j φm i,j (u, v) φm+1 k,l (u, v)J Ψ (u, v) du dv, = S m i,j
φ m i,j (ũ, ṽ)φ m+1 k,l (ũ, ṽ)J Ψ (u, v)J Λ| S m i,j dũ dṽ,

where Λ| S m i,j is the restriction of Λ to S m i,j . As seen in Section 3.3, J Λ is a constant function of ũ, ṽ so it can be taken out of the integral (J Λ| S m i,j depends on m though, since it is related to the area of S m i,j ). Now, if Ψ is energy-preserving, J Ψ is constant and can be taken out of the integral. In this case,

φm i,j , φm+1 k,l Ŝm i,j = J Ψ J Λ| S m i,j S m i,j φ m i,j (u, v)φ m+1 k,l (u, v) du dv = J Ψ J Λ| S m i,j h m (i,j),(k,l) .
Similarly,

φm i,j , φm k,l Ŝm i,j = J Ψ J Λ| S m i,j δ i,k δ j,l .
This means that the formula defining the push-pull coefficients reduces in the energypreserving case to a single term, ĥm (i,j),(k,l) = h m (i,j),(k,l) , which we know to be independent of m.

Wrapping up. What we just proved is an important result. First, it gives yet another hint that the virtual mesh with constant Jacobian mappings is the natural way of extending traditional radiosity techniques to complex primitives. Second, it means that computing push-pull coefficients on a parametric surface is very simple when using energy-preserving projections, while they involve a double sum of coefficients in the general case. Third, and most importantly, the push-pull coefficients ĥm (i,j),(k,l) do not depend on m. They can thus be precomputed and stored at the top level of the hierarchy, which simplifies the process of redistributing energy within the hierarchy. We have determined that this speeds up the process by a factor of 3 for linear wavelets and by a factor of 7 for quadratic wavelets. The explanation is simple: the number of basis functions increases at each level of the hierarchy as a function of the order n of the wavelet basis M n and the number of push-pull coefficients is the square of the number of basis functions. Working with constant Jacobian mappings thus has a deep impact on performance, since many new hierarchies are created all along the radiosity simulation.

It could be argued that the extra cost ascribable to the use of mappings with non-constant Jacobian must be small for Haar (M 1 ) wavelets. This is true, since it essentially amounts to recomputing the surface area of a patch. But as we shall see in the results section, the trade-off (in terms of computation time, number of final patches and visual accuracy) is clearly in favor of using M 2 or M 3 wavelets for illuminating objects with varying curvatures and normals. For such bases, using energy-preserving mappings is very important.

Self-shooting surfaces

Since with the virtual mesh method the input geometry is kept intact, a curved input patch may very well shoot energy onto itself. More precisely, every curved primitive that has at least a concave "island" will self-shoot. For every such surface, we create a decoy, i.e., a new surface with the same geometry but zero energy levels.

Handling self-shooting surfaces is then a two-step process. First, the surface shoots energy on its decoy as if it was a different surface. Then the energy levels, residuals and mesh subdivisions of the decoy are copied back to the original surface.

Managing the hierarchy

Suppose the image of a curved patch by a constant Jacobian mapping is the support of a wavelet approximation of the radiosity function. To benefit from the simple transfer coefficient computation of its virtual support, we proceed as if the virtual support was the real support of the function. All calculations are done on a unit square and so the hierarchical meshing is very natural (see Figure 9). Indeed, a simple quadtree algorithm will produce four identical squares (representing four regions of identical area on S) at each level which are assumed to be the new virtual supports of the parts of the original shape that fall over them. During subdivision, two special kinds of meshes may appear: empty meshes, which will simply not be considered for emission in future steps, and plain meshes, which can be managed with a classical hierarchical subdivision.

Subdividing the virtual support of a shape rather than the shape itself avoids

• • •
empty mesh plain mesh boundary mesh the use of complicated meshing schemes that would produce polygons which would be hard to manage. And the use of the virtual mesh does not simply reduce the number of initial primitives, it also reduces the number of final patches, as will appear clearly in the next section.

RESULTS

The wavelet radiosity algorithm with the virtual mesh presented above has been implemented. We now present results obtained with our implementation for various test scenes. Our implementation of the virtual mesh can currently handle complex polygons with arbitrary trimming curves, low-degree curved patches and developable surfaces.

We have used several test cases to illustrate different advantages and properties of the virtual mesh5 :

• the opera, to compare the illumination of scenes with complex polygons to premeshed scenes; • the teapot, to compare the illumination of scenes with curved patches to tessellated scenes; • the chess set, to compare the illumination of scenes with curved patches with or without constant Jacobian mappings; • the Reality Center TM and the Soda Hall, to show that our implementation can also handle real-size scenes with a reasonable time and memory consumption.

All test cases show that the virtual mesh does make a considerable difference performance-wise and memory-wise. We have decided not to use any form of clustering for our tests. Our motivations for this were twofold. First, although clustering offers an elegant theoretical solution for reducing the complexity of the radiosity algorithm, its practical use may behave unexpectedly and error control is difficult. Clustering is good to have a quick rough simulation of a scene but not to have simulations that are everywhere very accurate, of the kinds needed in lighting design. Second, we wanted to determine what performance gains could be ascribed to the use of the virtual mesh alone and not to other factors. In addition, since the virtual mesh can also be used in conjunction with clustering, we have good reasons to believe that the benefits of using the virtual mesh observed without clustering would be similar with clustering. The results were collected on a PC Pentium III 633 MHz with 500 MB of free memory and 500 MB of disk cache. We took care to use similar parameter settings for each test case (in particular, same minimum subdivision size and subdivision oracle threshold). Iterations of the radiosity method are stopped when a user-defined convergence rate is reached, defined as the ratio

E i -E r E i ,
where E i is the initial energy and E r is the energy remaining to propagate after the current iteration. Finally, the memory results we give are for the difference between the total memory use after simulation minus the resident set size of the process.

Complex vs. simple polygons

In architectural modeling, many polygons having non-trivial shapes appear. Up to now, these primitives had to be tessellated to be handled by radiosity algorithms. With the virtual mesh, each of these polygons is illuminated as a unique primitive. Our first test case illustrates this property and its impact on visual quality and computation time.

Our first scene, the opera, is a faithful modeling of an opera house of the classical style designed by Emmanuel Héré (18 th century). The original model has 17,272 complex polygons. When only simple polygons are allowed, the number of input primitives rises to 32,429 (see Figure. 10.a). Note that the triangulation algorithm used is based on the Open GL tessellator which may generate slivers. Clearly, better polygon triangulators exist (like the constrained Delaunay algorithm). But they generally produce a much larger number of triangles than Open GL. Since the performance of any radiosity algorithm depends heavily on the number of input surfaces, we have decided to stick with the Open GL tessellator so as to have meaningful comparisons. Qualitatively speaking, Figure 10 shows that the results are visually more pleasing when the virtual mesh is used to handle complex polygons. Clearly, separating the representation of the radiosity function from the geometry and using simple and computationally convenient domains to build the hierarchy yields a much better approximation of radiosity. When the tessellated model is used for the simulation, false geometric discontinuities due to the triangulation may appear, as is the case in Figure 10.c along the shadows at the bottom left of the right window frame. Also, problems due to shapes not well suited to visibility computations can be seen above the window: slivers are missed by the visibility algorithm.

Quantitatively speaking, Table 1 shows that the virtual mesh, in addition to improve the visual quality and accuracy of the simulation, also has a strong impact on performance. Indeed, for the same convergence rate of 99 %, the overall computation time is divided by 2. And as shown by Figure 11, the closer to the "ideal" solution one wants to be, the better it gets. In addition, the number of surfaces that have shot energy during the simulation is lower by more than 25 %: this is an important point which saves visibility computation time.

Note that the gain in terms of memory usage (45 %) is less important than what it is for running time. This is easily explained by the fact that subdivisions are not simply due to the shape of the input primitives but also to the complexity of the radiosity function to represent. Here, many subdivisions are due to the complex patterns of shades and shadows around the many windows of the model, thus generating many final meshes and increasing memory usage.

• 23 a.

b.

M 1 M 2 M 3 c. d.

Curved vs. tessellated primitives

With the teapot test scene (Figure 12), we illustrate the fact that it is faster and visually more realistic to illuminate scenes made of parametric surface patches directly using the virtual mesh rather than via triangulation. The teapot is made of 12 quadric surface patches (ellipsoids, cylinders, cones and hyperboloids of one sheet) and 7 complex polygons. Note that some of the patches, notably around the handle, have complex trimming curves. The teapot is illuminated by a white point light source which reflects on a planar red plate, so that there is a lot of indirect lighting. Comparison is made with a tessellation of the teapot. We have chosen a reasonable mesh size (1,392 triangles and parallelograms -see Figure 12.a) so that the teapot looks good enough when viewed from some distance: it is a good compromise between a very coarse mesh (a few hundred primitives) and a mesh fine enough (of the order of 10,000 triangles) to grab the full curvedness of the teapot. Each vertex of the tessellation is equipped with the real normal of the teapot at that point, so that normals can be interpolated during the simulation.

The results of our simulations, all made for the same convergence rate of 99 %, are given in Table 2. For the same wavelet basis (M 3 ), the illumination of the non-tessellated version of the teapot is more than twice faster than that of the tessellated one. And as illustrated by Figs. 12.c and .d, the visual result is more pleasing when the input surfaces are quadric patches. To obtain a similar visual result, the tessellated teapot would need a very large number of input patches (and thus of final patches, inducing a greater memory usage).

Table 2 shows that the number of final patches is not very different between the two simulations. This is easily explained by the fact that most subdivisions occur along the sharp shadow boundary that appears on the lid of the teapot and that many subdivisions are needed to "grab" the curve, both for the tessellated and the quadrics teapot. However, memory consumption is significantly lower for the illumination with curved patches and the virtual mesh: memory usage is linked to the number of hierarchies created during the simulation and is thus influenced also by the number of input primitives.

Note that the artifacts that one can observe on the tessellated teapot simulation (Fig. 12.c) are due to an acceleration of the algorithm: before considering a patch for subdivision, the oracle tests whether its orientation prevents it from being seen by the emitter by looking at the interpolated normal at a number of Gauss quadrature points. The one artifact that appears near the shadow boundary on the lid is thus the result of a small patch being marked as entirely not visible from the light source, while in fact its top left corner is visible from the source.

The teapot test case also gives another confirmation of the findings of Cuny et al. [START_REF] Cuny | A novel approach makes higher order wavelets really efficient for radiosity[END_REF]]: the benefits of higher-order wavelets outweigh their costs.

Constant vs. non-constant Jacobian mappings

The chess set test scene is used to illustrate the impact of the constant Jacobian property on the performance of the wavelet radiosity/virtual mesh algorithm. The scene is made of a chessboard and its chess pieces illuminated by a desk lamp. The chess set alone is made of 555 surface patches, 264 of which are curved (some of them having non-trivial trimming curves, as the knight of Fig. 13). Figure 13 shows a radiosity illumination of the scene with a final pass of ray tracing (3 bounces) that has a small impact on the overall performance (8 seconds for a 256 × 256 image).

We have compared the performance of two versions of the virtual mesh: using constant Jacobian mappings and using natural coordinates mapping of the input primitives. The main difference lies in the necessity, in the latter case, to compute the push-pull coefficients at each level of the hierarchy, as explained in Section 5.2. Our results are summarized in Table 3. They show that the constant Jacobian requirement is very important. We have determined that the algorithm does not Table 3. Results for the chess set test scene using linear wavelets, for a 90 % convergence rate.

even converge below 5 × 5 Gauss quadrature points per support for coefficient estimation and at least 10×10 points are needed to have a satisfactory visual result. Note that it could have been worse for other types of mappings: if a large region in parameter space represents only a small region in equal-area parameter space, then approximation errors can only be compensated by significantly increasing the number of quadrature points. But this has a deep impact on performance, as Table 3 indicates: when a 30×30 "grid" is used, which is the one we take for the computation of the top-level coefficients in the constant Jacobian version, the running time is almost four times slower for the natural coordinates simulation and memory usage is significantly larger.

Environments with real-world complexity

In this final section, we show that the virtual mesh can also efficiently and robustly handle scenes of high geometric complexity. The Reality Center TM test scene is a faithful modeling of the SGI immersive visualization facility located in Cortaillod, Switzerland. It is made of 2,561 surfaces, more than half of which are curved (mainly spherical, conical, cylindrical and toroïdal patches), illuminated by 41 point light sources. Figure 14 shows that this relatively low number of primitives is enough to represent a complex environment.

And we can zoom in on any part of the simulation while keeping a good visual quality. To reach the same visual result with simple polygons, the scene would need to be tessellated in more than 100,000 triangles, substantially hampering the performance. Table 4 summarizes the numerical results of this simulation on our PC Pentium III.

For larger scenes, the PC may not have enough memory for very accurate simulations. To demonstrate that the virtual mesh has no problem with such scenes, we have developed a parallel version of our wavelet radiosity algorithm on a distributed shared-memory machine, the SGI Origin 2000 with 4 195 MHz R10000 MIPS processors and 4 Gb of main memory (the details of the algorithm can be found in [START_REF] Paul | Partitioning and scheduling large radiosity computations in parallel[END_REF]). Our test scene is the 4 th floor6 of the famous Soda Hall (the computer science building on the campus of the University of California at Berkeley) with its furniture -see Figure 15 for an overview of the illuminated floor. While the original model is made of 272,450 simple polygons, allowing simple curved patches and complex polygons reduces this number to 70,645. Table 4 shows that this falloff does make a difference performance-and memory-wise when the virtual mesh is used, as already observed. The fact that the memory consumption is quite high is a reflection of the fact that there is a lot of data replication in the parallel algorithm.

Note that the images of Figure 1 were made from a simulation of the Soda Hall scene augmented with a special room (room 420A) furnished with several sorts of curved objects (teapot, chess set, benches, jars, desk lamps. . . ).

CONCLUSION

In this paper, we have presented a general method for separating the structure of the mesh from the complexity of the geometry in radiosity algorithms. At its core is the virtual mesh, an abstraction of the input geometry that is both simple and computationally convenient. The virtual mesh completely eliminates the effect of often highly tessellated but otherwise simple surface shapes on the algorithm. The method works for both trimmed parametric surfaces and complex polygons. We have shown how to interface the virtual mesh with a wavelet radiosity algorithm. Our implementation provides a hitherto missing link between hierarchical Galerkin radiosity and general trimmed curved surfaces. The virtual mesh has been shown to yield very good performances in terms of running time, memory usage and visual quality.

Our method is a new step towards making radiosity robust enough to be employed on useful scenes. Still, several problems remain to be addressed:

• The virtual mesh is not limited to Galerkin radiosity. Indeed, it is a very general method and can be combined with any improved radiosity algorithm. In the future, it will probably see applications to clustering algorithms and to handle glossy surfaces. We are currently working on coupling the virtual mesh with discontinuity meshing, to benefit from the best of the two methods and go one step further in improving the accuracy and realism of the simulations.

• Some triangle meshes (like a laser scan of the Venus de Milo) do not represent simple surface shapes: the underlying surfaces are curved but without a compact mathematical representation. For such surfaces, our method cannot be applied in its current form. However, we are currently investigating ways of globally parameterizing sets of triangles, as discussed in [START_REF] Lévy | Non-distorted texture mapping for sheared triangulated meshes[END_REF]]. The two main issues here are how to compute energy-preserving mappings for a given set of triangles and how to split a triangle mesh into useful sets of triangles.

• Up to recently, the future of parametric surfaces in computer graphics was open to question. Smoothed versions of polygonal surfaces seemed to take hold. But recently, parametric surfaces received a renewed interest since it was shown that subdivision surfaces are piecewise-parametric (Stam proved this for the Loop and Catmull-Clark subdivision schemes [START_REF] Stam | Exact evaluation of Catmull-Clark subdivision surfaces at arbitrary parameter values[END_REF]. We now need to evaluate how practical our method is for such general surfaces as NURBS and subdivision surfaces. One short-term work is to develop equal-area parameterizations for is thus energy-preserving for the spheroid. We obtain a mapping for the sphere (a = b) as a limit case when e → 0:

(x, y, z) -→ (θ, sin ϕ).

A.2 Flat surfaces

A flat surface is a surface for which the Gaussian curvature vanishes everywhere. Generalized cylinders, generalized cones and tangent developables are popular examples of flat surfaces. A generalized cylinder is defined by the following parametric equation:

x(θ, ϕ) = pϕ + y(θ),

where p is an arbitrary point and y(θ) is the parametric equation of the base curve C of the cylinder. The energy-preserving constraint of Equation ( 2) writes down as

v = c|f (u)| y (θ) × p g(u, v).
If C is planar, such that its equation is y(θ) = (y 1 (θ), y 2 (θ), 0) T , and p = (0, 0, 1) T (right generalized cylinder), then this constraint becomes: v = c|f (u)| y 2 1 (θ) + y 2 2 (θ)g(u, v). For instance, if C is a cardioid, parameterized by y(θ) = (a cos θ(1 + cos θ), a sin θ(1 + cos θ), 0) T , then we can take J Ψ = 2 √ 2a and the mapping (x, y, z) -→ (u, v) = ( √ 1 -cos θ, ϕ)

The Virtual Mesh

• 31 is energy-preserving.

A tangent developable is parameterized by

x(θ, ϕ) = y(θ) + ϕy (θ), with y(θ) is the equation of a base curve C as above. Equation ( 2) is here

v = c 2 |f (u)| y (θ) × y (θ) g(u, v) 2 .
If C is again a cardioid, then Ψ -1 : (x, y, z) -→ (θ+sin θ, ϕ 2 ) is an energy-preserving mapping such that J Ψ = 3a 2 2 .

Fig. 2 .

 2 Fig. 2. Overview of the illumination of a trimmed torus patch by four point light sources. a. The input scene. b. The outside of the torus patch mapped to a virtual 2D support. c. The illumination of the 2D support. d. Global illumination of the scene, with radiosity values on the torus deduced from the illumination of its virtual 2D support.

Fig. 4 .

 4 Fig. 4. Mapping a torus patch with a curve drawn on it to the domain [0, π] 2 . a. A ring torus patch. b. A mapping of this patch with constant Jacobian. c. A mapping with non-constant Jacobian.

Fig. 5 .

 5 Fig. 5. Smoothly extending the radiosity function.

Fig. 6 .

 6 Fig. 6. Trapezoidal map of an arrangement of line segments and curved arcs.

Fig. 7 .

 7 Fig. 7. Example of adjusting quadrature weights to virtual supports when S is considered an emitter. a. Classical 2D quadrature points and weights. b. Modified weights on virtual supports.

  HR-Transfer (S i -→ Sj , ∆ Bi ) 10: Bj ←-Bj + Ẽj , ∆ Bj ←-∆ Bj + Ẽj 11:

• 19 Fig. 8 .

 198 Fig. 8. A self-shooting surface shoots on a decoy.

Fig. 9 .

 9 Fig. 9. Hierarchical subdivision of the virtual support of the patch delimited by the curve drawn on the torus of Figure 4.a.

Fig. 10 .

 10 Fig. 10. The opera scene. a. A close-up on the initial mesh of the opera: with simple or complex primitives. b. c. Radiosity solution with (b.) and without (c.) the virtual mesh.

Fig. 12 .

 12 Fig. 12. The teapot test scene. a. The tessellated model. b. The teapot modeled with quadrics. c. Illumination of a. with different wavelet bases. d. Illumination of b. with different wavelet bases.

  Fig. 13. Illumination of the chess set with constant Jacobian mappings.

• 27 Fig. 14 .

 2714 Fig. 14. Simulation of the Reality Center TM test scene.

Fig. 15 .

 15 Fig. 15. Aerial view of the illumination of a floor of the Soda Hall.

  So one can take JΨ = a 2 √ a 2 + b 2 , f (u) = u and g(u, v) = √ v as solution. The mapping (x, y, z) -→ (u, v) = (θ, ϕ 2 )is thus energy-preserving. Sphere and spheroid. A prolate spheroid is defined by sin θ sin ϕ, y = a cos θ sin ϕ, z = b cos ϕ, with b > a. Equation (2) writes down as v = c 2 ab |f (u)| sin g 1 -e 2 sin 2 g + 1 e arcsin (e sin g) , where e = 1 -a 2b 2 > 0 is the ellipticity of the spheroid. To determine an energy-preserving mapping for this surface, one can take J Ψ = ab 2 and u = θ. The mapping (x, y, z) -→ θ, sin ϕ 1 -e 2 sin 2 ϕ + 1 e arcsin (e sin ϕ)

Table 1 .

 1 Results for the opera test scene using Haar wavelets, for a 99 % convergence rate. Fig. 11. Evolution of the convergence rate as function of time for the opera test scene.

	geometry	# input prim.	# final meshes	computation time (s)	memory usage (Mb)	# shooting surfaces
	simple polygons	32,429	819,109	952		203.5	9,358
	complex polygons	17,272	550,159	482		140.3	6,776
	gain	87 %	49 %		98 %		45 %	38 %
	conv.	100
	rate			complex polygons
		99			triangles
		95
		90
		80
		50
		20
		0	100	200	300	400	time

Table 2 .

 2 Results for the teapot scene with various wavelet bases, with or without tessellation of the input geometry, for the same global error and the same convergence rate (99 %).

	geometry	wavelet basis	# input primitives	# final meshes	running time (s)	memory usage (Mb)	# shooting surfaces
	simple polygons	M 1		52,346	795	28.1	1,300
	with interpolation	M 2	1,392	25,475	558	23.6	1,378
	of normals	M 3		18,217	605	22.4	1,410
	quadrics and	M 1		39,895	583	9.2
	constant Jacobian	M 2	19	20,920	341	7.8	21
	mappings	M 3		16,537	262	8.1

Table 4 .

 4 Results for the Reality Center TM (linear wavelets) and Soda Hall (quadratic wavelets) test scenes, for a 95 % convergence rate.

	scene	geometry	# input surfaces	# final meshes	running time (s)	memory usage (Mb)
	Reality Center TM on PC Pentium III	low-degree patches	2,561	119,486	1,510	157.1
		simple polygons	272,450	1,124,585	44,585	3,425
	Soda Hall on SGI	complex polygons	199,550	927,428	30,232	2,612
	Origin 2000	quadrics and complex polygons	70,645	827,62	13,411	1,210

In the rest of this paper, "simple polygons" is used to denote triangles and/or parallelograms and "complex polygons" to denote arbitrary generalized polygons, i.e. possibly non-convex, multiply connected polygons with holes and curved trimming boundaries.

Energy-preserving mappings are also known as equal-area projections in cartography.

See[Kaplan 1984, pp. 98-99 and 238-245] for precise mathematical definitions and properties.

The notation Mn refers to the basis of multi-wavelets of order n, i.e., the wavelets whose smoothing functions are tensor products of the first n Legendre polynomials[START_REF] Gortler | Wavelet radiosity[END_REF]].

Additional examples and images available at http://www.vsp-tech.com.

Available at http://www.cs.berkeley.edu/~kofler.

ACKNOWLEDGMENTS

The authors wish to thank the anonymous reviewers for useful comments and suggestions.

• 29 cubic Bézier patches.

APPENDIX A. EQUAL-AREA MAPPINGS FOR CURVED PATCHES

Following the ideas developed in Sections 3 and 4, the virtual mesh can be applied to a variety of complex shapes. The only thing that needs to be done to add a new shape to the dictionary of primitives that our wavelet radiosity algorithm can illuminate is to determine an energy-preserving mapping for this shape. We illustrate here several simple but important cases: natural quadrics and flat surfaces.

A.1 Natural quadrics

Because of their fairly important modeling power and modest complexity, quadrics (i.e., algebraic surfaces of degree 2) are important in geometric modeling and visualization. Some quadrics are more important than others. Patches of natural quadrics (planes, cones, spheres and cylinders) make up, according to [START_REF] Nourse | Natural quadrics in mechanical design[END_REF], 85 % of all mechanical pieces.

This section shows how to compute energy-preserving mappings for natural quadrics. Constant Jacobian mappings for the remaining types of quadrics can be found in [START_REF] Cuny | Radiosité à base d'ondelettes sur des surfaces paramétriques[END_REF]].

Cylinder. Consider a right circular cylinder given in canonical form by:

Using the same notations as in Section 3.1, where θ = f (u) and ϕ = g(u, v),

Thus we can take J Ψ = 1/c = 1, g(u, v) = ϕ = v and f (u) = θ = u. In other words, the "natural coordinates" mapping (x, y, z) -→ (θ, ϕ) is in this case energypreserving.

Cone. A circular cone can be described by the equation

The constant Jacobian constraint of the mapping (x, y, z) → (u, v) boils down to