Fran Coise Andr

Marc Le Fur

Yves Mah

Jean-Louis Pazat
email: pazat@irisa.fr

Parallelization of a Wave Propagation Application using a Data Parallel Compiler

Abstract| This paper presents the parallelization process of a Wave Propagation application using the Pandore environment. Tools are brie y described, the stress is put on the description of the parallelization by data distribution and performance results are shown.

I. Introduction

The di culty to program massively parallel architectures with distributed memory is a severe impediment to the use of these parallel machines. In the past few years, the data parallel model has been used to de ne new languages such as HPF 7], tools and compilers: the programmer is provided a familiar uniform logic address space and a sequential ow of control. He controls the distributed aspect of the computation by specifying the data distribution on the local memories of the processors. The compiler generates code according to the spmd model and the links between the code execution and the data distribution is enforced by the owner-writes rule.

To achieve good performance when following this approach, several sophisticated compilation techniques and run-time systems have been studied and integrated into environments [START_REF] Guidec | POM: a Virtual Parallel Machine Featuring Observation Mechanisms[END_REF][START_REF] Bareau | Parallel Program Performance Debugging with the Pandore II Environment[END_REF]1,[START_REF] Brezany | Automatic Parallelization for GENESIS[END_REF]12]. Among these environments, the Pandore environment allows the compilation of both HPF and Pandore programs into spmd machine independent code 2]. A series of experiments on classical kernels have already led to satisfactory results. The next step is the validation of the compiler and the run-time system on real applications; the Wave Propagation application presented here is one of them.

The paper is organized as follows: we brie y present the Pandore environment and give an overview of the C-Pandore language. The compilation schemes are then expounded and the distribution of the well-known Jacobi kernel is detailed. Finally, the di erent steps of the parallelization of the Wave Propagation application are described and the results of this experiment are discussed.

II. The Pandore Environment

The Pandore environment has been designed to facilitate the programming of data distributed applications for distributed memory computers or clusters of workstations. Figure 1 shows the components of the Pandore environment.

The source program can be written in HPF or in a dialect of C (C-Pandore) augmented with data distribution features. In the rst case, a source to source translator is used to translate a subset of HPF into C-Pandore. This translator 13] is built upon a Fortran 90/HPF Front-end and a C-code generator written at GMD-First with the Cocktail toolbox 8].

From the source program, the Pandore compiler automatically generates a machine independent spmd code according to the owner-writes rule. Optimizations are included in the compiler that is described in section IV. For performance measurements we also provide a pro ler and a post-mortem analysis tool 3].

The Pandore run-time is in charge of the management of distributed arrays 15]; it uses a generic message passing library called POM (Parallel Observable Machine 9]). This library o ers limited but ecient services that makes it possible to run the same program on a wide range of distributed memory computers such as the iPSC or the Paragon; clusters of workstations are also supported. This library is also used by other compilers and tools developed in our team.

III. The Pandore Language

The Pandore language is based on a sequential imperative language using a subset of C (excluding pointers) as a basis. We have added a small set of simple and well-de ned data distribution features in order to describe frequently used decompositions. Alignments can be easily included and we plan to handle them in a near future.

A Pandore program is a sequential program which calls distributed phases. The sequential part is in charge of all I/O operations and is executed on the host processor (if exists) or on one node of the distributed computer. Each distributed phase is spread over the processors of the target machine and is executed in parallel according to the owner-writes rule. The speci cation of a distributed phase is described similarly to the de nition of a procedure by the statement dist d-phase (distributed parameter list) fd-blockg

The distributed parameter list is used for specifying the partitioning and the mapping of the data used in the distributed phase. The array is the only data type that may be partitioned. The means to decompose an array is to split it into blocks. The number of blocks is independent of the number of processors: both block and cyclic(k) HPF distribution features 7] are handled.

The speci cation of the partitioning for a d-dimensional array is given by the construct block (t 1 ; :::; t d ) where t i indicates the size of the blocks in the i th dimension. For example: B N] N] by block (N; N=P) indicates that the array B of N N elements is decomposed into P blocks of size N N=P: the array is decomposed into P blocks of contiguous columns. The following partitioning A N] N] by block (1; N) indicates that the array A of N N elements is decomposed into blocks of size 1 N: the array is decomposed into N rows.

Then, the mapping of the blocks onto the architecture will be achieved by the compiler in a regular or cyclic way according to the mapping parameters (regular or wrapped). In Pandore, we consider only one dimensional processor arrays whose size is not speci ed in the source code but used as a parameter by the compiler. As we allow the mapping of multidimensional decompositions, it is needed to indicate the order for the mapping of blocks: (1,0) states for column rst, (0,1) states for row rst. For example oat A N] N] by block(1;N) map wrapped(1;0) maps the N rows of A cyclically onto the processors; it is equivalent to the HPF distribution CYCLIC(1,*). The distribution speci cation oat B N] N] by block(N;N=P) map wrapped(1;0) maps the blocks of N=P columns onto the processors. If there are P processors, the mapping is similar to a HPF block decomposition BLOCK(N=P , *).

The last speci cation given in the parameter list concerns the transfer mode for values between the caller and the distributed phase: allowed modes are IN, OUT and INOUT. This speci cation is similar to the Fortran90 intent attribute. Figure 2 shows an example of a distributed phase.

#define N 512 #define P 4 dist jacobi(double B N] N] by block(N,N/P) map regular(0,1) mode INOUT) double A N] N] by block(N,N/P) map regular(0,1); { int i,j; for (j=1; j<N-1; j++) for (i=1; i<N-1; i++) language, with no direct relation with distribution, to improve the ease of programming. Ordinary C functions are not allowed in the Pandore language but in addition to distributed phases, two features are o ered to the programmer: macros and closed functions. Macro declarations are similar to procedure de nitions but parameters are passed \by name" and calls to a macro are in-lined by the compiler. Closed functions are similar to C functions but cannot access global variables nor modify distributed arrays. Close functions are similar to HPF PURE functions.

A i] j] = 0.5 * B i] j] + 0.125 * (B i-1] j] + B i+1] j] + B i] j-1] + B i] j+1]); for (j=1; j<N-1; j++) for (i=1; i<N-1; i++) B i] j] = A i] j]; } SUBROUTINE JACOBI (B)
,J) + B(I+1,J) + B(I,J-1) + B(I,J+1)) END DO END DO B(1:N-2, 1:N-2) = A(1:N-2, 1:N-2) END SUBROUTINE JACOBI

IV. The Pandore Compiler

The compiler produce spmd code from the usersupplied data decomposition according to the ownerwrites rule. Two compilation schemes are embedded in the compiler. For reductions and parallel loops with one statement, a ne bounds and array references, the compiler applies an optimized scheme 14] performing loop bounds reduction and message vectorization.

For statements that cannot be optimized, the compiler relies on the well-known runtime resolution technique. Because distributed array management is a critical point to achieve good performances, an original distributed array management based on paging 15] has been developed to support both schemes.

A. Compilation of Parallel Loops

For reductions and parallel loops, the compiler generates a code that comprises two parts: a communication part and then a computation part. Only one computation code is generated so, no di erence is made between local computations (that is those involving only local data) and non-local computations (that need to fetch data owned by other processors). Avoiding this separation does not a ect the performance of the generated code since our distributed array management provides a uniform and e cient access method for local data and copies of distant data. Besides, performing this separation does not seem realistic in the general case, with regard to compilation time and code fragmentation.

Actually, the compiler generates a series of communication codes. One communication code is produced for each right hand side reference to a distributed array and is decomposed in its turn into a send part and a receive part. Loop bounds and array subscripts but also the distribution of the arrays involved in the computation are analyzed by the compiler. For a given right hand side reference to a distributed array, the associated set of data that must be exchanged between processors is characterized by a polyhedron whose scanning 11, 6] provides the spmd send code and receive code for the reference. Then, the way arrays are represented in the local memories is taken into account by the compiler so that the data to be moved from one processor to another are scanned in the appropriate direction. This permits the transfer of contiguous zones (both on the sender and the receiver) and so eliminates the need of coding/decoding and copying between message bu ers and local memories. The generation of the spmd computation code rely on the same technique: according to the analysis of the left hand side reference, the compiler constructs a polyhedron which de nes the set of iterations that must be performed on each processor.

B. Management of Distributed Arrays

The distributed array management that completes the two compilation schemes balances the memory requirements and the speed of accesses to local data. It provides a uniform representation for local data and copies of distant data. Each block of a distributed array is decomposed into pages thus an array is represented on a processor by a table of pages that contains both local pages (pages of the blocks owned by the processor) and distant pages (copies of pages owned by other processors). For a given distributed array, the direction and the size of its pages are determined by the compiler so that the global to local index transformation involves only low level operations (shifts, masks) and the size of the table of pages is minimized.

V. The Jacobi Kernel

The Jacobi Relaxation Iterative Method can be used to approximate the solution of a partial di erential equation discretized on a grid. At each time step, the current approximation is updated by computing for each grid point the weighted average of the values of the neighboring points. We focus here on the kernel of this algorithm that consists of two loop nests working on two 2-D arrays A and B. The rst loop nest one computes in array A the current approximation from the values stored in array B that represents the last approximation. The second one transfers elements of A into B.

The distributed phase corresponding to the kernel is shown in gure 2. Arrays A and B are both distributed into P groups of columns (blocks of size N N=P), one group on each of the P processors. The mapping is not signi cant here as there is only one block per processor. With such a distribution, the workload is evenly distributed among the processors and the second loop nest is executed without any communication because arrays A and B are fully aligned. Communication is needed in the rst loop nest. Indeed, computing elements of the boundaries of each block necessitates accessing elements situated on the neighboring processors. The symmetry of the accesses would permit a row-wise distribution, leading to the same cost of communication. However, as in both loop nests elements are accessed column-wise (loop j is the outer loop), locality is best exploited with a columnwise distribution 1 .

As they conform with the restrictions of the application of the polyhedron-based compilation scheme, these two loop nests are fully optimized by the Pandore compiler. Although its scope is wider, the joint use of the optimized scheme and the run-time system prove to be as e cient as more classic compilation methods such as the overlap 16]. Indeed, the following optimizations automatically apply in this case: Iteration domains are restricted, messages are fully vectorized, direct unbu ered communications are used and index conversions are reduced to the identity function (the page number is the column number and the page o set is the row number).

The performances of the produced code for various input sizes 2 are summarized in table 1. They are almost optimal for small numbers of processors. For a given array size, the number of operations needing only local data performed by a processor is inversely proportional to the number of processors. The boundaries are of xed size, so performances decrease with large numbers of processors. However, one can notice that performances remain at a good level even with small data sizes: for N = 128, the e ciency |i.e. the ratio tp ts p where tp is the time for the parallel Pandore program and ts the time for the sequential C program| reaches 67% for [START_REF] Mah | Distributed Array Management for HPF Compilers[END_REF] Several tests on classical algorithms such as matrixmatrix multiplication, Gram-Schmidt or LU factor-1 The compilation process does not perturb the loop order when restricting the computation loop domains 2 The memory size on the ipsc/2 nodes did not allow the program execution for N = 512 and P = 2 ization have already been conducted to evaluate the Pandore compiler 2]; the next step of the validation of the compiler goes through experimentations on real applications.

A. The Wave Propagation Application

We describe here the parallelization/distribution with Pandore of a wave propagation algorithm that has been developed by the French Petroleum Institute for use by seismology experts to analyze the impacts of seismic shocks. The application, whose core is about one thousand line long, simulates the wave propagation in a bounded 2D space. Waves are generated by an explosion triggered at a given point of the considered space. The program studies the temporal evolution of the waves at several points of the space where some sensors are located. It takes as input a number of simulation parameters such as the time and space steps, the frequency of the explosion source and the positions of the sensors.

The numerical algorithm follows a discretizednite element method. It corresponds to the second order time discretization (V (t = n + 1) = F(V (t n), V (t = n 1))) and to the second order discretization of the spatial partial derivatives from the continuous system of PDE describing the waves propagation in an heterogeneous domain.

The results of the algorithm are the values of the horizontal and vertical movements at each time step, for each sensor. The program is divided into two phases:

The initialization phase: it de nes the initial conditions of the explosion and the constraints associated with the nature of the propagation domain.

The computation phase: this phase consists of a main loop representing the time evolution. At each iteration step we compute the horizontal and vertical movements at time t+1 and t+2, for each point of the grid representing the 2D space.

The nal results, i.e. the movements associated to the sensors, correspond to a grid sampling. During the computation phase, for each movement, four arrays are used: U p and U m for the horizontal movements and W p and W m for the vertical ones.

The body of the loop is composed of four similar parts corresponding to the following computations:

8 > > < > > : U p (t + 1) = f(U m (t); U p (t + 1)) W p (t + 1) = g(W m (t); W p (t + 1)) U m (t + 2) = f(U p (t + 1); U m (t)) W m (t + 2) = g(W p (t + 1); W m (t))
Functions f and g comprises two series of nested loops. The rst one is a series of 2-deep loops operating on the inner part of the grid. These loops are comparable to the rst part of the Jacobi kernel presented earlier. The second series is made of several 1-deep loops operating on the upper border of the grid.

B. Distribution of the Program

We describe here the steps to transform the initial sequential program into a C-Pandore one. Only one distributed phase is needed for this application, corresponding to the initialization phase followed by the computation phase. So, using the dist construct of the C-Pandore language to encapsulate this phase, all the computation will be automatically distributed on the nodes of the target parallel computer.

Exchanges between the host (or a dedicated node) and the computing nodes are only performed at the beginning of the distributed phase (to send the simulation parameters) and at the end of the computation (to transfer the nal results giving the movements associated with the sensors). These exchanges will be automatically handled by the compiler according to the parameters speci cation of the distributed phase.

The body itself of the distributed phase has been slightly modi ed in order to exhibit parallel nested loops which conform to the conditions under which the compiler may perform loop optimization. We easily obtain that for the four main computation parts described in A.because the loops naturally appear as parallel loop nests with a ne array references and loop bounds.

The main task when writing the C-Pandore program resides in the the choice of the array decompositions. The main arrays, described in A., are 2D arrays representing the propagation space grid. They are used together with a tenth of 2D coe cient arrays of the same type. Eight other 1D arrays are used for the computation of the border of the grid.

For this algorithm, we chose a column-wise decomposition for all the 2D arrays because in the computation of the values associated with the inner part of the grid, the dependencies are similar those found in the Jacobi leading to this decomposition as one of the best choices; the 1-deep loops operating on the upper border of the grid ( rst row of the arrays) necessitates a column-wise decomposition in order to distribute the workload evenly; the 2-deep loops are column oriented so a column decomposition enforces the locality for the greatest part of the computation.

Given P, the number of processors, each (N,N) array is partitioned into blocks of size (N,N/P). So we obtain a C-Pandore program that, except for the speci cation of the distribution phase and the partitioning of the arrays, corresponds almost exactly to the sequential one.

C. Performance Results

C.1. First Results

We ran the above described version of the Wave Propagation program on the ipsc/2. The performance results are given in table 2. The overall performances are satisfactory, considering that the parallel code has been produced automatically. With a good adequation between the data size and the number of processors, an e ciency around 70% is reached. Furthermore, it can be noticed that, for a given array size, performances decrease few when adding processors. In the C-Pandore program obtained in section B. the parallel loops performing the sampling lead to communications that could be avoided. This is due only to the fact that alignment cannot be expressed yet in C-Pandore; hence, the arrays storing the movements associated with the sensors are not aligned with the arrays representing the space grid. However, a trivial manual renumbering of the sensors points is possible to make the sampling arrays aligned with the grid arrays. After this transformation, the sampling loops are executed without communication. The performances are not very much a ected by this modication as the sampling only represent 5% of the total execution time.

Avoiding Multiple Transfers

During execution some array elements are sent several times to the same processor without having been modi ed because loops are optimized independently and there is no inter-loop def-use analysis for non-local data. These unnecessary transfers may be avoided by storing the data, after the rst send, in auxiliary arrays. We experimented this {non trivial{ optimization which necessitates to declare new arrays. The performances are improved by about 5% of the total execution time. This does not appear as a signi cant gain, considering that for this application the memory cost is of great importance. In fact, the size of the memory on each node of the parallel computer severely limits the experiments that may be conducted.

VII. Conclusion

We have shown in this paper that the Pandore compiler is able to distribute e ciently a real application, without any signi cant e ort from the programmer. Indeed, the rst Pandore source of the Wave Propagation program presented in the paper is very similar to the sequential original program and produces correct performances. This tends to con rm the viability of the data-parallel approach for scienti c computing provided that general enough and non-naive compiling and run-time techniques are applied.

To handle very large applications, in particular applications that comprise multiple modules or necessitate intensive I/O, other techniques must be integrated to existing environments. For this purpose, the joint study of redistribution, procedures and separate compilation is under way in the Pandore project.

Figure 1 :

 1 Figure 1: The Pandore Environment

  INTEGER, PARAMETER :: N = 512 REAL(KIND=8), DIMENSION(0:N-1,0:N-1) :: B REAL(KIND=8), DIMENSION(0:N-1,0:N-1) :: A !HPF$ PROCESSORS PROCS(4) !HPF$ DISTRIBUTE (BLOCK, *) ONTO PROCS :: A, B INTEGER I, J DO J=1, N-2 DO I=1, N-2 A(I,J) = 0.5 * B(I,J) + 0.125 * (B(I-1

Figure 2 :

 2 Figure 2: Kernel of the Jacobi algorithm in C-Pandore and HPF Some other constructs have been added to the C

Table 1 :

 1 processors although 1=4 of the columns are exchanged in the rst loop nest. Performances results for the Jacobi Kernel

	N	128	256	512
	p time (s) e . time (s) e . time (s) e . 2 0.242 86% 0.970 87% | | 4 0.128 81% 0.495 86% 1.961 86% 8 0.068 76% 0.253 84% 0.999 85% 16 0.039 67% 0.132 80% 0.503 84% 32 0.027 48% 0.075 70% 0.261 81%
	VI. Parallelization of a Wave Propagation Application with Pandore

Table 2 :

 2 Performance results for the Wave Propaga-

	tion Application
	C.2. Further Optimizations
	Sampling Associated with the Sensors


C Bareau

Y Mah Eo

J.-L Pazat

Parallel Program Performance Debugging with the Pandore II Environment

Keywords:

In this paper, we present the overall design of Pandore II, an environment dedicated to the experimentation of distribution of sequential programs for their execution on distributed memory parallel architectures. The emphasis is then put on two performance analysis tools integrated in this environment.

INTRODUCTION

Programming distributed memory parallel computers (dmpcs for short) is a di cult task because of the management of parallel processes. The intricateness between computation and interprocess communication is often very tight, leading to very cumbersome programs. So the wish of most users is to be freed from all these \low level" details. In many cases they do not want to take into account the distributed aspects of their programs, which they also want to be machine independent.

A solution to this problem is given by the use of sequential programming languages extended with data distribution features like High Performance Fortran 7]. In this case, the compiler is in charge of generating communicating parallel processes from the sequential code and the data distribution speci cation. This solution has been validated 4], and implemented in the Pandore II environment 1]. Yet, e ciency issues necessitate to develop speci c performance analysis tools.

THE PANDORE II ENVIRONMENT

Pandore II is an environment designed for parallel execution of imperative sequential programs on dmpcs. It comprises a compiler, libraries for di erent dmpcs and execution analysis tools including a pro ler and a trace generator.

The source language

The source language for this environment is a subset of C augmented with features for data decomposition purpose { similar characteristics appear in the recently de ned High Performance Fortran language. No speci c knowledge of the target machine is required of the user: only the speci cation of data decomposition is left to his duty.

A Pandore II program is a collection of distributed phases which may be seen like functions. The speci cation of data distribution is expressed as attributes of the formal parameters. Arrays can be partitioned into rectangular blocks and their mapping on processors can be regular or cyclic. For example : dist myphase(float A N] N] by block(N,1) map wrapped(0,1) mode INOUT) speci es that the array A is to be partitioned into columns mapped cyclically onto the processors. The INOUT mode speci es that A is to be read at the beginning of the phase and written back at the end. Figure 1 shows an example of a Pandore II program.

The compiler

The Pandore II compiler generates parallel processes according to the data decomposition speci ed by the programmer. The compilation scheme is based on the locality of writes, on the host/node and the SPMD model. The host process executes the code of the main() part of the program, whereas node processes execute the SPMD code produced by the compiler from the distributed phases. A node process is in charge of its local variables: it updates their values and sends them to other processes when needed, according to the original sequential program. For example the assignment

A i] = B i + 1] + C i] is translated into: refresh(ftmp1,tmp2g, fB i+1],C i]g, owner(A i])) exec(owner(A i]), tmp1+tmp2) free(ftmp1,tmp2g)
When executing the refresh macro, owners of B i + 1] and C i] send their values to the owner of A i]. This process receives the values into bu ers tmp1 and tmp2. The exec macro is a guarded command that insures that only the owner of A i] executes the statement A i] =tmp1+tmp2. If the assigned variable is replicated on all the processors (as scalars for example), distant values are broadcasted through the network.

This basic translation scheme is not very e cient but optimization techniques based on the same model exist, for example 2] carries out a static domain analysis of loops to generate e cient code.

The runtime library

The runtime permits the execution of object code on di erent dmpcs. Its goal is to implement memory and process management, communication of data elements between processes, distributed data management and e cient index translation. It relies on machine dependent libraries and is implemented using macro de nitions.

The machine model is a fully connected network of processing elements, communicating through reliable FIFO channels. Sends are non-blocking, whereas receives are blocking. The design and the implementation of this library increase the system portability and facilitates the instrumentation of the code.

Need for performance debugging

The performances of the code generated by the Pandore II system are dependent on the appropriateness of the chosen data distribution to the algorithm but also on strategies directing the compiler and the runtime implementation. It appears that evaluating the in uence of these parameters is necessary in order to guide the system designers and to provide the user with tools helping him to distribute his data. A dynamic evaluation { as opposed to a static estimation 3, 6] { o ers the advantage of being applicable to every type of program and yields precise results.

Two techniques are used for performance measurement in the Pandore II environment: tracing and pro ling. These two techniques di er in their aim and in their implementation. Tracing permits the recording of events to which are assigned at least a type and a timestamp. With this method, the parallel activity can be recorded in order to rebuild the program behavior; we present here an extension of usual tracing for extracting information on all the potential behaviors of the program. The counterpart of this method is pro ling, whose aim is to gather enough statistics for execution analysis. A number of counters related to events are updated during program execution. In Pandore II, an enhancement of mere pro ling is used: in addition to their occurrences, the durations of events may also be cumulated 9].

TRACE GENERATION AND ANALYSIS 3.1. Principle of the Pandore II trace analyzer

The principle of performance analysis by trace generation is to study the causal structure of the program that takes into account the load balancing induced by the chosen data distribution and the compilation scheme. Indeed, from an intuitive point of view, \a distributed program is very parallel when, most of the time, most processors can actually perform an action, i.e. are not blocked waiting for a message". Now these blockings, due to the asynchronism of the communications, happen when an event is causally dependent on an event on a di erent processor; thus, a tool that makes these dependences between processors clearly visible gives an estimation of their importance in terms of e ciency. Moreover, by focusing on events rather than on relations between them, one may get an estimation of the load balance of the program.

Our tool permits to visualize whether an action creates a bottleneck or can be performed in parallel with many other actions on other processors. For this purpose, we build the lattice of all the possible behaviors of the program: it is a graph in which each vertex represents an instant of the execution and the outgoing edges are the actions that may be performed at this instant by the unblocked processors. Hence each path from the initial vertex to the nal one corresponds to one of the interleaving of the actions of the program.

For instance, from the following program, where x i] and y j] are placed on processor P1 and z i] on processor P2, Source code execution on P1: on P2:

z i]:=y j]+1; send(y j],P2); (a) receive(y j],P1); (e) x i]:=5;

x i]:=5;

(b) z i]:=y j]+1; (f) y j]:=x i]+3*z i]; receive(z i],P2); (c) send(z i],P1); (g) z i]:=z i]-3; y j]:=x i]+3*z i]; (d) z i]:=z i]-3; (h)
we get, giving a direction to each processor, the lattice of executions of gure 2(a). The main drawback of the lattice of executions is its size. Furthermore, events are not all interesting: for instance when an assignment always follows a reception, one of these two events could be abstracted without loss of precision. The user may also want to analyze only parts of its program. In that respect, we allow him to include his own observation points either in the Pandore II source program or in the generated code. In the previous example, if we choose to observe only assignments, we obtain the new lattice shown in gure 2(b).

Construction

The construction of the lattice of executions is based on results of order theory, using the fact that a distributed execution can be seen as a partial order. An algorithm 10] that computes vector timestamps coding this order has been implemented in ECHIDNA 8], a programming environment for execution of Estelle speci cations on dmpcs, networks of workstations, or by simulation on monoprocessors 1 . In order to use this tool, an Estelle-code generator has been added as new back-end to the Pandore II environment. The execution times of the Estelle code are of course very di erent from those obtained with the C code, but the notion of intrusion is here irrelevant for we focus on the very structure of the algorithm, on the causal dependences between events that ensue only from the source, not from the low-level mechanisms.

The lattice of executions is actually the lattice of the ideals of an execution seen as a partial order, therefore it can be created from only one execution. We use an algorithm described in 5]; this algorithm is online: during an execution, the observed events generate traces with timestamps coding the partial order, and the lattice is incrementally constructed from these traces.

Applications

(a) • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • (b) • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • (c) • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • Figure 3.

Matrix product lattices of executions

With our tool, the user is given an intuitive look over the performance of his distributed program. The lattice of executions highlights the bottlenecks as well as the \fully parallel" parts of the program, therefore it is possible to have an overall estimation of the performance, and thus to compare several data distributions. Furthermore, when a bottleneck appears in the lattice, the user immediately knows which events are causing it, and so which part of the program should be reconsidered.

As an example, we show in gure 3 the lattices corresponding to a matrix product on four processors with three di erent distributions. The observed events are the assignments of the program. It can be seen that in the distribution (c) there are many points with only one or two outgoing edges: parallelism is here very weak, especially at the beginning and at the end of the execution when only one processor achieves assignments, the others performing only communication actions. On the contrary, in (a), most points correspond to instants when three or four processors are \ready to work", leading to a greater amount of potential parallelism. This approach may seem unrealistic if we consider lattices with many processors; however we observe that such lattices present the same overall shape (for example the juxtaposition of the same pattern in distribution (a)), so the observations made with few processors are pertinent.

THE PANDORE II PROFILER 4.1. Instrumentation

The Pandore II pro ler allows the user to collect a number of quantitative measures on his program's execution with minimal intervention. The use of pro ling restrains the amount of storage needed; the number of counters to be updated is of the order of the number of variables declared in the source program. This pro ler has been implemented on a 32-node iPSC/2 but is easily portable. Sensors are inserted in modi ed versions of some runtime macros, thus the compiler generates a similar code whether an instrumentation is demanded or not. Lapses of time are measured with a software microsecond clock.

As most information for updating counters is available at compile time, the level of intrusion remains low (limited to a few percent execution overhead). Measurements are performed on each node and counters are brought back to the host at the end of the execution and then written down into a le that can be exploited by appropriate tools. The host code is not instrumented due to the lack of precision of time measurement on a time-shared multi-user system.

The links between the source and the evaluation results are established two di erent ways: rst the user bounds fragments of the distributed phases he wants to be evaluated by de ning some instrumentation zones, typically loop nests. Moreover, output gures are associated to objects of the source program such as arrays, scalars, conditional statements or loops.

Results

Besides execution times and the load balance, the main results produced by a pro led execution are related to communications and synchronizations. They may be classi ed in two categories: measures speci c to distributed phases (communication with host at the beginning and at the end of each phase, phase triggering) and measures concerning assignments within instrumentation zones.

These statistics give information about the e ciency of the runtime implementation especially for message passing. Moreover, with the last class of results, data distribution for a given algorithm can be evaluated. An assignment of a distributed array element in which another distributed array reference appears in the right hand side may generate a message from the owner of the right hand side to the owner of the left hand side. The aim of the measurement is to globally build a directed graph where vertices are array partitions and arcs describe the tra c between partitions. Arcs are valued by the number of messages, the transferred volume or the waiting time on reception. For example, the assignment A 3; 5] = B 4] will increase the value of the arc (B1 ! A2) if element A 3; 5] is on processor 2 and B 4] on processor 1.

In a similar way, an assignment of a replicated variable in which a distributed array reference appears in the right hand side will generate a broadcast message from the owner of the array element. For the entire execution, broadcasting is described (number of messages, transferred volume, waiting time on reception) for each pair (var,part), where the partition part represents the source and var the assigned replicated variable. For example, the assignment x = A 3; 5] will increase the value of the counters related to the pair (x; A1) if A 3; 5] is located on processor 1.

The produced results may be analyzed as such or treated by a speci c tool that can give partial and abstracted views (e.g. by selecting or grouping processors or variables). Another choice would be to rst look at the second loop nest. This would lead to the decomposition of A into groups of N=P columns (distribution b). One could also think to an intermediate solution: decomposition of A into (2 N=P; 2 N=P) blocks (distribution c).

V1 V0 V3 V2 A0 A1 A2 A3 A0 A1 A2 A3 A0 A2 A1 A3 V1 V0 V3 V2 (
Figure 4 gives the communication graphs for the three distributions. The row-wise distribution seems preferable. This is con rmed by the waiting time graphs which show strong synchronization for the column-wise and block-wise distributions.

CONCLUSION

The approach of distribution of sequential programs by data distribution is now recognized. We have presented Pandore II, a complete environment for experimenting this method. As e ciency is a key issue, there is a great need for performance evaluation. However, because of the speci city of the codes generated by systems like Pandore II, usual performance debugging tools are not well adapted. Therefore, we have designed new tools and integrated them in our environment. They are based on two complementary techniques of execution analysis (tracing and pro ling) which permit qualitative and quantitative evaluation. They are aimed to help the user to distribute his program's data and to give information to the system designers. They have already been employed to improve the compiler and the runtime; nevertheless, experimentation must be pursued in order to tune these tools as well as the compilation and runtime techniques involved in the environment.

Figure 1 .

 1 Figure 1. Example of Pandore II source program

Figure 2 .

 2 Figure 2. Lattices of executions

Figure 4 .

 4 Figure 4. Communication graphs for 3 distributions

Estelle is an ISO language for protocol speci cation.

for (i=0; i<N; i++) /* */ for (j=0; j<N; j++) /* Instrumentation zone 1 */ V[i] = f(V[i/* __________________ */ for (j=0; j<N; j++) /* */ for (i=1; i<N-1; i++) /* Instrumentation zone 2 */


Nader Salman

Mariette Yvinec

Quentin Merigot

Feature Preserving Mesh Generation from 3D Point Clouds

Keywords: Descriptors (according to ACM CCS): [Computational Geometry and Object Modeling] [I.3.5]: Curve, surface, solid, and object representations-

We address the problem of generating quality surface triangle meshes from 3D point clouds sampled on piecewise smooth surfaces. Using a feature detection process based on the covariance matrices of Voronoi cells, we first extract from the point cloud a set of sharp features. Our algorithm also runs on the input point cloud a reconstruction process, such as Poisson reconstruction, providing an implicit surface. A feature preserving variant of a Delaunay refinement process is then used to generate a mesh approximating the implicit surface and containing a faithful representation of the extracted sharp edges. Such a mesh provides an enhanced trade-off between accuracy and mesh complexity. The whole process is robust to noise and made versatile through a small set of parameters which govern the mesh sizing, approximation error and shape of the elements. We demonstrate the effectiveness of our method on a variety of models including laser scanned datasets ranging from indoor to outdoor scenes.

Introduction

Surface reconstruction consists in inferring the shape of an object from a 3D point cloud sampled on the surface. The problem of automatic surface reconstruction from point samples has received considerable amount of efforts in the past twenty years. Surface reconstruction problems occurs in different flavors according to the quality of the data which are more or less entangled with noise, and the properties of the reconstructed object whose surface can be smooth or structured with sharp features. Moreover, an implicit description of the reconstructed object is enough for some applications, whereas some other require a data structure more amenable to further processing, like a mesh.

In this paper, we focus on the case where the sampled surface is piecewise smooth and the output sought after is a high quality triangle mesh approximating the surface. Our goal is to generate from the input point cloud a mesh including a faithful representation of sharp edges, thus improving the trade-off between the accuracy and the size off the mesh. We propose for that an algorithm which experimentally proves to be fairly robust to noise.

Related work

Reconstruction and quality surface mesh generation are often tackled separately in the literature. We first review the approaches dealing with the problem of robust reconstruction of piecewise smooth surfaces, before reviewing the ap-proaches that combine reconstruction and quality surface mesh generation.

Towards piecewise smooth reconstructions.

In the quest to achieve both noise robustness and recovery of sharp features, we can identify an evolution of approaches ranging from smooth to sharp.

Smooth surfaces. As noise and sparseness are common in laser scanner datasets, we have seen in recent years a variety of approaches for robust reconstruction of smooth surfaces. A popular approach consists of searching an implicit function whose zero level set fits the input data points. The searched implicit functions can take different forms, being for instance approximated signed distance functions to the inferred surface [START_REF] Ohtake | Sparse lowdegree implicit surfaces with applications to high quality rendering, feature extraction, and smoothing[END_REF][START_REF] Guennebaud | Algebraic point set surfaces[END_REF] or approximated indicator functions of the inferred solid [START_REF] Kazhdan | Poisson Surface Reconstruction[END_REF]. Implicit functions methods often provide spectacular robustness to noise. However most of them assume that the data points come with oriented normals and are sampled from a smooth surface. The latter assumption leads to erroneous surface reconstruction at locations where the smoothness assumption is not met.

An intermediate step toward handling sharp features consists of using locally adapted anisotropic basis functions when computing the implicit function [START_REF] Dinh | Reconstructing surfaces using anisotropic basis functions[END_REF]. Adamson and Alexa [START_REF] Adamson | Anisotropic point set surfaces[END_REF] perform piecewise polynomial approximation through anisotropic moving least-squares (MLS). More recently, Oztireli et al. [OGG] used a kernel regres-Figure 1: Overview of the method at a glance. From left to right: given a point cloud sampled on a piecewise smooth surface; detect the potential sharp edges points, cluster them with respect to the underlying sharp features direction, and extract explicit polylines for each cluster; sample protecting balls on the extracted polylines; combine extracted polylines to implicit surface reconstruction through feature preserving mesh generation to output a 3D mesh where the sharp edges are explicitly tesselated and appear as a subset of the triangle edges.

sion technique to extend the MLS reconstruction to surfaces with sharp features.

While satisfactory for some applications (e.g., visualization as opposed to mesh generation), these approaches generate surfaces which are more accurate around sharp features but are still smooth. In addition, boosting the anisotropy of a function is a dangerous game when simultaneously targeting robustness to noise, as such a strategy may amplify noise instead of the true features (features and noise are ambiguous as both are high frequency). Furthermore, a purely local distinction between smooth parts and sharp edges is likely to yield fragmented edges while a global feature extraction approach can favor long edges in order to recover the true structure.

Handling sharp features. This brings out another thread of work aimed at extracting long sharp features. Pauly et al. [START_REF] Keiser R | Multi-scale feature extraction on point-sampled surfaces[END_REF] use a multi-scale approach to decide if a point is on a feature, and construct a minimum-spanning tree to recover the feature graph. Daniels et al. [START_REF] Daniels | Robust smooth feature extraction from point clouds[END_REF] extract featurecurves using a method based on a robust MLS projection operator which locally projects points onto the sharp edges, and grow a set of polylines through the projected points. Jenke et al. [START_REF] Jenke | Patch-graph reconstruction for piecewise smooth surfaces[END_REF] extract the feature lines by robustly fitting local surface patches and by computing the intersection of close patches with dissimilar normals. Both methods show satisfactory results even for defect-laden point sets. More recently, Merigot et al [START_REF] Mérigot | Robust voronoi-based curvature and feature estimation[END_REF] proposed a robust estimation of curvature tensors based on a covariance measure defined from the data points Voronoi cells. Their method provides theoretical guarantees on robustness under bad sampling and noise conditions.

One way to make a local, binary distinction between smooth parts and sharp features consists of performing a local clustering of the inferred normals to the surface [OBA * 03]. If this process reveals more than one cluster, the algorithm does not fit just one low-degree implicit function, but as many quadrics as the number of clusters. This leads to faithful local reconstruction of a sharp edge as the intersection of two implicit primitives in the case of two clusters. A corner is reconstructed as the intersection of three or more primitives.

Note in passing that the method misses tips (sharp point that are not incident to any sharp crease), cusps and darts.

In order to achieve improved robustness, Fleishman et al. [START_REF] Fleishman | Robust moving least-squares fitting with sharp features[END_REF] segment neighborhoods of points by growing regions belonging to the same part of the surface. Robustness is obtained by using a forward search technique which finds reference planes corresponding to each region. While constituting an important progress, the method requires very dense point clouds, and is rather expensive. Also, potential instabilities in the classification can create fragmented surface parts. Targeting even higher robustness, Lipman et al. [START_REF] Lipman | Datadependent MLS for faithful surface approximation[END_REF] enrich the MLS projection framework with sharp edges by defining a singularity indicator field based on the error of the MLS approximation. However, they choose to restrict their approach to a single singularity within each influence radius, which significantly limits the type of feature points that can be handled.

Coupling reconstruction and mesh generation.

One popular approach that couples reconstruction and mesh generation is the extended marching cubes introduced by Kobbelt et al. [START_REF] Kobbelt | Feature sensitive surface extraction from volume data[END_REF]. The input is assumed to be a discretized signed distance function enriched with oriented normals, the latter being used to decide whether or not a voxel contains a sharp feature. If a voxel is classified as containing a feature, some additional vertices are constructed within the voxel and placed at intersections between the planes defined by the vertices and their affiliated normals. This approach works robustly and efficiently by truly coupling feature extraction and mesh generation. Nevertheless, it relies on the marching cubes process which is well known to produce rather poor quality surface meshes.

Two approaches which are able to generate quality surfaces while preserving the features [DHOS07, JWSA08] proceed sequentially by first extracting a set of features, then using the advancing front method proposed by Schreiner et al. [START_REF] Schreiner | Direct (re) meshing for efficient surface processing[END_REF]. The latter produces quality isotropic meshes from different inputs including MLS representations, and can control the density of the mesh by pre-computing a guidance field. The algorithm uses time-consuming pre-processing of the point cloud in order to construct a guidance field, as well as a final post-processing using MLS projections.

submitted to Eurographics Symposium on Geometry Processing (2010)

A relevant alternative to the marching cubes and advancing front strategies is the Delaunay refinement paradigm [START_REF] Boissonnat | Provably good sampling and meshing of surfaces[END_REF]. In this method the surface mesh is intersection free by construction as filtered out of a 3D Delaunay triangulation. A number of additional guarantees are also provided after termination of the refinement process, such as a good shape of elements, a faithful approximation of geometry and normals, and a low complexity of the mesh. More interestingly in our context, Delaunay refinement is able to couple reconstruction and mesh generation. At the intuitive level, the refinement procedure is combined with a sensing algorithm probing an implicit surface defined from the data points. The probing is performed along Voronoi edges of sampling points which are not only longer than the short edges of the marching cubes but also data-dependent as they become more and more orthogonal to the sensed surface as the refinement process goes along. However, in [START_REF] Boissonnat | Provably good sampling and meshing of surfaces[END_REF] sharp features are not explicitly handled, hence they are not accurately represented in the output mesh. One way to circumvent this drawback is to extract the feature lines and to constrain the Delaunay refinement to preserve them. However, it is well known that the presence of constrained edges where meeting surface patches form small angles, endangers the termination of the Delaunay refinement [START_REF] Shewchuk | Mesh generation for domains with small angles[END_REF]. Recently Cheng et al. [START_REF] Cheng S.-W | Delaunay refinement for piecewise smooth complexes[END_REF][START_REF] Cheng | A Practical Delaunay Meshing Algorithm for a Large Class of Domains*[END_REF] propose to deal with this problem using the method of protecting balls set up around sharp features.

Contributions

In this paper, we argue that some modifications of the implicit smooth surface reconstruction algorithms can produce results with accuracy on par with the best current piecewise smooth reconstruction methods. Given a 3D point cloud sampled on a piecewise smooth surface, our algorithm consists of two steps. In the first step, we extract from the input 3D point cloud a set of polylines approximating the sharp edges of the object. This feature extraction is based on the robust feature estimation framework of Merigot et al. [START_REF] Mérigot | Robust voronoi-based curvature and feature estimation[END_REF]. In parallel, our algorithm runs on the data points a reconstruction process, providing an implicit surface, i.e. an implicit function whose zero level sets approximates the data points. The second step of our algorithm generates a mesh approximating the implicit surface while including a faithful representation of the extracted sharp edges. This step runs a feature preserving version of the Delaunay refinement meshing algorithm, based on the meshing algorithm of Boissonnat and Oudot [START_REF] Boissonnat | Provably good sampling and meshing of surfaces[END_REF] and using the mechanism of protecting balls introduced by Cheng et al. [START_REF] Cheng S.-W | Delaunay refinement for piecewise smooth complexes[END_REF][START_REF] Cheng | A Practical Delaunay Meshing Algorithm for a Large Class of Domains*[END_REF]. An overview of the method is illustrated in Figure 1.

The technical tools used in this algorithm are directly inspired by recent work [START_REF] Mérigot | Robust voronoi-based curvature and feature estimation[END_REF][START_REF] Cheng | A Practical Delaunay Meshing Algorithm for a Large Class of Domains*[END_REF]. However, they have not been used previously in this combination, and we argue that it is this particular synthesis of existing ideas which is the key of its success.

The main benefits of our algorithm are: Accuracy. The method provides a meshed representation of a point sampled object with an enhanced trade-off between accuracy and mesh size. Robustness. The method is robust to noise and stable under sparse sampling. Flexibility. The method offers considerable flexibility as it can be tuned with respect to the scale and the targeted size of the mesh.

Feature extraction

Feature detection

A point on a piecewise-smooth surface, is called a (sharp) feature point, if the surface has no tangent plane at this point. The locus of feature points on a piecewise-smooth surface has the structure of a graph.

We are given a point cloud P sampling a piecewise-smooth surface S. In order to allow the mesh generation algorithm to faithfully reproduce the sharp features of S, we need to extract from the input point cloud an approximation of the feature graph of S. Our first goal is to identify the data points that are close to sharp edges. To this aim, we use the method of Merigot et al. [START_REF] Mérigot | Robust voronoi-based curvature and feature estimation[END_REF] based on an analysis of the shape of Voronoi cells through their covariance matrices.

Recall that the covariance matrix of a bounded domain E of R 3 with respect to a base point p is defined by cov(E, p) = R E (xp)(x-p) T dx. The eigenvectors of this matrix capture the principal axes of E while the ratio of the eigenvalues gives information on the anisotropy of E.

The Voronoi diagram of the sampling P is the cellular decomposition of the space, where each sample p ∈ P has a cell, V (p), that is the locus of points closer to p than to any other sample in P :

V (p) = {x ∈ R d : ∀q ∈ P, x -p ≤ x -q }
By construction, the shape of Voronoi cells depends on the whole set of data, and some cells are unbounded. In order to get more local information, [START_REF] Mérigot | Robust voronoi-based curvature and feature estimation[END_REF] introduce an offset parameter R and consider for each data point p of P, the covariance matrix, M(p, R), of the intersection of the Voronoi cell V (p) with a ball, B(p, R), of radius R centered at p:

M(p, R) = cov (V (p) ∩ B(p, R), p)
The principal drawbacks of Voronoi-based methods is their high sensitivity to noise. To alleviate the effect of noise, [START_REF] Mérigot | Robust voronoi-based curvature and feature estimation[END_REF] smoothes the information contained in the covariance matrices using convolution. The method computes for each point p the convolved Voronoi covariance matrix (CVCM), M c (p, R, r), by summing the covariance matrices of all points q of P at a distance less than r to x:

M c (p, R, r) = ∑ q∈B(p,r) M(q, R)
The crucial remark is that sample points lying on smooth surfaces, called smooth points, have a pencil shaped Voronoi cell; hence their CVCM has two "small" eigenvalues and a larger submitted to Eurographics Symposium on Geometry Processing (2010)

one, with the eigenvector corresponding to the large eigenvalue directed along the surface normal. Sample points that are close to sharp edges, called edge points, have rather flat Voronoi cells; therefore there CVCM has a single "small" eigenvalue and two larger ones, and the eigenvector corresponding to the small eigenvalue is directed along the edge. Sample points close to feature graph nodes, called corner points, tend to have cone-shaped Voronoi cells thus all three eigenvalues have comparable values.

Based on those observations, the algorithm we use to detect feature points (edge and corner points) is the following. We compute the CVCM, M c (p, R, r), of each point p ∈ P, and sort the eigenvalues of this matrix by decreasing order λ 0 (p) ≥ λ 1 (p) ≥ λ 2 (p). If the ratio λ 0 (p)/λ 1 (p), is greater than a threshold t 1 , p is a smooth point. Else, if the ratio between λ 1 (p)/λ 2 (p), is greater than a threshold t 2 , p is an edge point. Else, point p is a corner point. We gather, in a set F, the pairs (p, e 2 (p)), where p is a point classified as an edge point and e 2 (p) is the eigenvector of its CVCM corresponding to the smallest eigenvalue λ 2 (p), which is supposed to be aligned with the sharp edge direction. Note that we choose not to keep in F , the pairs corresponding to corner points because of the instabilities of the data around these positions. Instead, we will recover corner points as extremities of sharp edges (Section 2.3).

Parameters:The method needs four parameters: the offset radius R, the convolution radius r and both thresholds t 1 ,t 2 . Parameter R should be chosen smaller than to the external reach of the underlying surface defined as the smallest distance from a point of the surface to the external medial axis. Parameter r has two different roles: besides smoothing the information, it determines the curvature radius under which a curved area is considered as a sharp feature. Consequently, it should be chosen depending on the noise level and on the sizing field of the targeted mesh. In all our experiments, we fixed parameter R to be one tenth of the point cloud bounding sphere radius and r to be a tenth of R, which leads to good results. Moreover, we allow the user to modify t 1 and t 2 manually, using an interactive interface.

Algorithm 1 Clustering points in F Place each point p i ∈ F in its own cluster F i with it's potential extremity points E i = ∅ for all points p i do for all p j ∈ F such that d(p i , p j ) < ρ k (p i ) do if p j and p i are NOT in the same cluster (

F j = F i ) then Compute angle α i = ∠(e 2 (p i ), p i p j ) Compute angle α j = ∠(e 2 (p j ), p i p j ) if α i ≤ θ and α j ≤ θ then Merge the clusters F i and F j Merge E i and E j if α i ≤ θ and α j > θ then E i ← p j

Feature clustering

In [START_REF] Mérigot | Robust voronoi-based curvature and feature estimation[END_REF], it is proven that the convolved covariance matrix is robust to noise, by deriving a bound on the quality of the results as a function of the Hausdorff distance between the point cloud and the sampled surface. However, the reconstruction of the sharp feature graph is left unaddressed in this work. Here, we need to approximate the sharp edges with polygonal lines. In their work [START_REF] Daniels | Robust smooth feature extraction from point clouds[END_REF] presented a technique for solving a similar problem. We however proceed differently taking advantage of the edge direction information we get from the CVCM analysis.

Our goal here is to group edge points of F into clusters, where each cluster samples a straight or bending sharp edge.

For that, we need a similarity measure that will quantify if two edge points belong to the same sharp edge.

A pair (p i , p j ) of edge points is said to be a close pair if the Euclidean distance d(p i , p j ) is less than ρ k (p i ), where

ρ k (p i ) = 1 k ∑ k l=1 |p i -p l | 2
, the sum running on the k nearest neighbors to p i in F . Two close points p i , p j ∈ F belong to the same cluster if the angles α i = ∠(e 2 (p i ), p i p j ) and α j = ∠(e 2 (p j ), p i p j ) are smaller than a given threshold θ.

This criteria is injected in an union-find algorithm (Algorithm 1) that partitions points of F into clusters F i . In addition, if only one of the angular conditions is satisfied for a close pair (p i , p j ), say α i < θ but α j ≥ θ, point p j is classified as an extremity point of cluster F i and included in a set E i of extremity points related to F i . These extremity points will be of use to elongate the constructed polylines in Section 2.3, facilitating the resolution of the polylines junction recovery problem. Moreover, our implementation prunes clusters consisting of few points (|F i | ≤ 6) and unlabels the potential extremity points associated to them.

Parameters:The identification of good clusters relies on an appropriate selection of parameters θ and k. For big values of θ the detected clusters will be thick, thus it will be hard to sort the clustered points, a step we need for recovering the features. For our purposes, θ = 15 • and k = 20 provided sufficiently thin and representative clusters to facilitate the feature recovery procedure (Figure 2), without the need to apply the thinning heuristics of [START_REF] Lee | Curve reconstruction from unorganized points[END_REF].

submitted to Eurographics Symposium on Geometry Processing (2010)

Feature recovery

In the following our aim is to recover from each cluster, provided by the previous step, a polyline that will best fit the shape of the underlying sharp edge.

For each cluster F i we build a polyline Poly(F i ) by applying iteratively a procedure succ(p) that discovers the successor of a point p (Figure 3 top). We start at a random seed point p ∈ F i . The successor s(p) of p is defined as follows: we consider the subset N(p) ⊂ F i of points included in a ball of radius rp centered at p. All points in N(p) are then projected onto the line defined by p and the edge direction e 2 (p). Points in N(p) can be divided into two subsets by inspecting the position of their projections relative to p. If p is the starting point, we choose randomly one of the subsets to be N + (p). Otherwise, let N + (p) be the subset of points whose projections are in the direction opposite to the predecessor of p. Then, s(p) is defined as the point in N + (p) with closest projection to p. The procedure succ(p) is applied iteratively until an endpoint is reached, i.e. a point for which N + (p) is empty. The procedure is then restarted from the seed point in the opposite direction.

The polyline Poly(F i ) and the edge directions attached to its vertices are then smoothed using a standard PCA based smoothing process.

To lengthen the recovered feature polylines Poly(F i ), we use the detected potential extremity points E i attached to each cluster. For each endpoint p ∈ F i , if E i is not empty, we search for neighboring points p j ∈ E i such as the Euclidean distance d(p, p j ) < rp. All such points are projected onto the line defined by p and the edge direction e 2 (p). Then, we compute the centroid c of the projected points located on the half-line defined by p and not containing the projection of N(p). Finally, point c is added to the polyline as the successor of p and attached a direction e 2 (c) = e 2 (p). Point c is considered as the new endpoint of Poly(F i ) (see Figure 3 bottom).

Our implementation maintains for each polyline Poly(F i ) a value corresponding to the polyline length that will be of use in the feature junction recovery step to prevent joining the endpoints of a single short polyline.

Parameters:The recovered polyline result is dependent on the size of the selected neighborhood. The radius parameter rp should be related to the local feature size of the polyline to be recovered. Optimally, one should study the effect of the neighborhood size and propose an optimal neighborhood-size-selection scheme. Practically, we do not have this information, hence we choose rp = ρ k (p)

where k = 5.

Junction recovery. In order to build a realistic approximation of the sharp feature graph, we need to construct appropriate junctions between the constructed polylines. The task is simplified as we have extended the polylines and thus reduced the search space for potential junctions. Note that

P oly(F i ) E i p c p F i s(p) N + (p)
Figure 3: Illustration of feature recovery step. Top: procedure succ(p) that discovers the successor, s(p), of a point p ∈ F i . Bottom: after smoothing the polyline Poly(F i ), we elongate it at each endpoint p using the detected extremity points E i .

sharp edges joined at vertices of the feature graph are classified as follows: a dart is incident to a single sharp edge; a cusp is incident to two sharp edges forming a sharp angle at the junction; a corner is incident to three or more sharp edges.

A natural solution is to merge in a single vertex, polylines endpoints that are spatially close. We could simply take as position for the merged vertex the barycenter of the merged endpoints. However, to increase the accuracy of the approximation of sharp feature junctions we proceed as follows.

For each endpoint p of Poly(F i ), we start by constructing the set J (p) of endpoints that are within a ball B(p, r junc ) of radius r junc centered at p. In case there is no endpoint in J (p) other than p itself, we tag p as a dart and continue.

Lj

L k p k p j q j q k c
Otherwise, |J (p)| ≥ 2 and we compute the junction position as follows. For every pair of endpoints p j , p k ∈ J (p), we consider the two points q j , q k which are the closest points such that q j belongs to the line L j passing through p j and directed along e 2 (p j ) and q k belongs to the line L k passing through p k and directed along e 2 (p k ). The junction point is computed as the barycenter of the set of closest points, except that each closest point that is not in the ball B(p, r junc ) is replaced by the endpoint p in the barycenter computation.

Our procedure has the advantage of being easy to implement and is computationally inexpensive since it does not introduce any plane nor surface fitting and intersection computations. An example of finding junctions from polylines including three junction cases, darts, cusps and corners is illustrated in Figure 5. Note that we do not capture tips, i.e. vertices incident to zero sharp edges but are sharp.

Parameters:The radius parameter r junc should be chosen large enough, in order to cover the empty space between polyline endpoints. On the other hand, r junc generally shouldn't be larger than the length of the shortest polyline -which would lead to both of its endpoints being merged into a single junction. In practice, we found that choosing r junc to be equal to five times the average spacing of the ending four points of a polyline leads to good results. 3 Feature preserving mesh generation Besides extracting sharp edges, our algorithm runs on the data points a reconstruction process, such as Poisson reconstruction [START_REF] Kazhdan | Poisson Surface Reconstruction[END_REF] or any variant of moving least square [START_REF] Guennebaud | Algebraic point set surfaces[END_REF], providing an implicit description of a surface approximating the data points. Such an implicit surface is represented as a function (R 3 → R) whose zero level set approximates the data points. Our goal in this section is to show how, using a feature preserving variant of a Delaunay refinement mesh generation algorithm, we obtain a surface mesh approximating the implicit surface and including a faithful representation of the sharp edges extracted in the previous section.

Delaunay refinement surface mesh generation

To generate a mesh approximating the implicit surface, we use the surface meshing algorithm of Boissonnat and Oudot [START_REF] Boissonnat | Provably good sampling and meshing of surfaces[END_REF]. This algorithm is a Delaunay refinement, based on the notion of restricted Delaunay triangulation.

Let P be a set of points in R 3 and S be a surface. The restricted Delaunay triangulation, D S (P), is the sub-complex of the Delaunay triangulation D(P) formed by the faces of D(P) whose dual Voronoi faces intersect S. It has been proven that if P is a "sufficiently dense" sample of S, D S (P) is homeomorph to the surface S and is a good approximation of S in the Hausdorff sense [ES97, AB99, BO05].

The meshing algorithm of [START_REF] Boissonnat | Provably good sampling and meshing of surfaces[END_REF] builds on this result. It maintains the restricted Delaunay triangulation D S (P) of a sampling P of the surface, and refines the sampling P until it is dense enough for the restricted Delaunay triangulation to be a good approximation of the surface. More precisely, the refinement process tracks bad facets in D S (P), i.e. facets that do not comply to the meshing criteria. Each bad facet is killed out by the insertion of a new vertex at its surface center which is the point where the Voronoi edge, dual to the facet, intersects the surface.

Such a Delaunay refinement surface meshing algorithm yields a quality surface mesh, free of self-intersections. Parameters:The parameters of the Delaunay refinement are the parameters of the criteria defining the bad facets. A facet is bad if either some facet angle is smaller than α, or some facet edge is longer than l or the distance between the center of the facet and the surface is more than d.

Feature preserving extension

To achieve an accurate representation of the extracted sharp features in the final mesh, we resort to the method of protecting balls proposed by Cheng et al. [START_REF] Cheng S.-W | Delaunay refinement for piecewise smooth complexes[END_REF] and experimented in [START_REF] Cheng | A Practical Delaunay Meshing Algorithm for a Large Class of Domains*[END_REF]. This method performs, before the Delaunay refinement, a protecting step in which each sharp feature is covered by a set of protecting balls, such that: ⊲ protecting balls are centered on sharp features and each feauture is completely covered by the union of protecting balls centered on it; ⊲ any two balls have an empty intersection, except balls with consecutive centers on a given sharp feature that intersect significatively but do not include each other's center; ⊲ any three balls have an empty intersection.

Once the protecting balls have been computed, they are regarded as weighted points and inserted altogether, as initial vertices, in a weighted Delaunay triangulation. The Delaunay refinement process is then performed using this weighted Delaunay triangulation: at each step, the refinement process computes the weighted version of a bad facet surface center. This point is affected a zero weight and inserted in the weigthed Delaunay triangulation.

Such a weight setting, together with the protecting ball properties, enforces the fact that two protecting ball centers, consecutive on a sharp feature polyline are guaranteed to remain connected by a restricted Delaunay edge. Furthermore, owing to the fact that no three balls intersect, the refinement process never inserts refinement point in the union of the protecting balls nor probes the surface in this protected region. This ensures the termination of the refinement process whatever may be the dihedral angle formed by smooth patches incident on a sharp feature.

submitted to Eurographics Symposium on Geometry Processing (2010) Protecting balls have a crucial influence on the sizing of the output mesh. On one hand protecting ball centers define the approximation of sharp features in the final mesh, and the sizing parameter l of the Delaunay refinement cannot be smaller than the spacing between protecting ball centers. On the other hand, to provide a surface mesh with facets incident to the mesh edges approximating a sharp feauture, the refinement process requires that the union of protecting balls covers the gap between the implicit surface (orange) and the extracted feature polylines (green, blue).

In practice, to compute the protecting balls, we first keep all endpoint vertices of the detected polylines Poly(F i ) and then sample points on Poly(F i ) according to some user defined uniform distance d ′ . Then the algorithm tries setting a uniform weight wp = 2/3 * d ′ on sample points. If the resulting set of balls do not comply with the above rules, it is locally fixed by adding new sample points and reducing ball radii.

Parameters:In our experiments, the mesh generation step depends on a single user-defined parameter δ. Indeed, we use as unit length the quantity avg defined as the mean of the distances from each input data point to each one of its six first nearest neighbors. Then, we fix the protecting ball parameter and Delaunay refinement parameters to:

d ′ = δ * avg, α = 25 • , l = 2 * d ′ , d = 0.6 * d ′ .

Experimental results

Our prototype is implemented in C++ using the robust primitives provided by the CGAL library [cga]. We use CGAL 3D triangulation as the core data structure to compute the Voronoi covariance matrices of the input point cloud [START_REF] Mérigot | Robust voronoi-based curvature and feature estimation[END_REF]. The CGAL library also provides us with an implementation of the Poisson reconstruction method [START_REF] Kazhdan | Poisson Surface Reconstruction[END_REF] to compute an implicit surface from the input point cloud.

Lastly CGAL provides an implementation of the Delaunay refinement surface meshing algorithm of Boissonnat and Oudot [START_REF] Boissonnat | Provably good sampling and meshing of surfaces[END_REF].

A crucial component for reaching good timings is the efficient answer to queries for the subset of a point cloud contained in a given ball. Such a query is required for computing the CVCM used to detect edge points, it is also used for clustering edge points, and finally for recovering and joining the feature polylines. We have chosen to implement these queries using the ANN library [ann]. Note that Poisson reconstruction assumes that the input points come with oriented normals. Here, we take advantage of the result from the eigenvector analysis of the CVCM. Each data point is attached a normal direction corresponding to the eigenvector of its CVCM with largest eigenvalue. The normals are oriented using the method described by Hoppe et al. in [HDD * 92] and provided by the point set processing package of CGAL.

Experiments where conducted on various datasets, we report here some of them: five "clean" synthetic models (fandisk, octa-flower, carved object, blade, block), a synthetic model with added noise (fandisk), the raw output of a Minolta TM laser scanner for an indoor model (Ramses), and the raw output of a Leica TM laser scanner for an outdoor urban scene (Church of Lans le Villard).

All results presented here have been produced using fixed values of the parameters that are mentioned in the Parameters paragraph at each of the algorithm steps, except for the three parameters t 1 ,t 2 and δ.

Accuracy versus mesh size. Figure 6, illustrates the effect of parameter δ on the output mesh resolution. Here we have chosen t 1 = 50, t 2 = 6 and protected the extracted polylines with balls obtained for respectively δ = 3.5 (a), δ = 7 for (b).

We further evaluated our method on the three synthetic datasets, shown in Figure 7. The octa-flower model is interesting as it has curvy, sharp features. The two other models carved_object, and blade show the ability of our approach to submitted to Eurographics Symposium on Geometry Processing (2010) As already mentioned, our goal is to improve the accuracy for a given mesh size, by ensuring a faithful representation of sharp edges. To illustrate our achievement with repect to this goal, we show, for each model in Figure 7, two meshes computed with the same value or parameter δ, hence the same value of Delaunay refinement algorithm parameters

(α = 25 • , l = 2 * d ′ , d = 0.6 * d ′ ).
For quantitative evaluation of the accuracy improvement, Figure 8 shows an approximation of the Hausdorff distance between generated meshes and the ground truth model for various size of the generated mesh. Note that, the feature preserving Delaunay refinement (FPE) always outperforms standard Delaunay refinement (PIF) and that the effect is all the more important for coarse meshes.

Figure 9 is another illustration of the ability of our method to enhance the trade-off between accuracy and mesh size through feature extraction and preservation. Note that the only way for standard mesh generation methods (marching cubes and Delaunay" refinement) to overcome the lack of feature lines is to generate high resolution meshes.

Complex datasets. In order to evaluate the ability of our feature extraction method to cope with noise, we perturbed the points sampled on the fandisk model. The perturbation is uniform and has an amplitude of 1% of the diagonal of the dataset bounding box. Keeping the same set of parameters (t 1 ,t 2 ) as for the unperturbed fandisk model, Figure 10 shows the detected edge points. Although the edge points are more diffused with this amount of noise, the computed feature directions remain quite stable. Moreover, despite introducing strong noise, the feature extraction step constructed almost all sharp edges.

To evaluate the practical usability of this new approach, we applied it to two laser scan datasets acquired from two different environments. The first one is an indoor laser acquisition of the Ramses model, Figure 11. The second is an outdoor laser acquisition of the Church of Lans le Villard (France), Figure 12. These datasets pose some problems compared to the synthetic models presented previously as near sharp features data are sparse and entangled with a high level of noise.

The total time required to construct the feature polylines depends on the number of detected sharp edges points. Our prototype program has not been optimized and the details shown in Table 1 are just given for indication (on a laptop with 3.5Gb memory, Intel R Core TM 2 2.7Ghz processor).

Model/ Parameters |F | Timings [sec] #points t 1 , t 2 , δ (1)
Limitations. Since feature extraction and smooth surface reconstruction are run independently, the extracted feature graph may have a certain discrepancy with respect to the smooth reconstructed surface. The feature preserving mesh generation algorithm yields a surface mesh with a faithfull representation of sharp features provided that the protecting balls cover the gap between the extracted features and the reconstructed smooth surface. This condition puts one more constraint on the parameter δ and, hence on the mesh accuracy and sizing.

The pipeline currently described does not support adaptative sampling and the generation of meshes with non uniform sizing field. This however is not a real limitation of the method but simply results from our concern to limit and simplify the parameters of the algorithm.

Conclusion

In this paper, we have presented an efficient and robust feature preserving mesh generation strategy from 3D point clouds. The approach we use builds a bridge between implicit surface reconstruction and mesh generation.

As for future work, it would be interesting to explore the possibility to automatically pre-compute the different parameters from the input point cloud. We are also planning to exploit the mesh generation algorithm flexibility to apply our method to more challenging datasets provided by multi-view stereo reconstruction. We also plan to make our code available very soon and we hope that it will be useful in the fields which require simplified and realistic geometric models, like in urban modelling.

Figure 2 :

 2 Figure 2: Result of the clustering step on the octa-flower model.

Figure 4 :

 4 Figure 4: Left: Feature recovery for the octa-flower model. Right: Closeup view of the recovered sharp edges and the recovered junction (black dot).

Figure 5 :

 5 Figure 5: Junction recovery. (a) a dart and multiple corners are recovered from the smooth feature model; (b) multiple border cusps are recovered from the half cylinder model

Figure 6 :

 6 Figure 6: Feature preserving mesh generation for the fandisk model. The threshold δ is used to adjust the size of the protective balls thus the mesh resolution; for (a) and (b), the meshed Poisson implicit function (PIF) is on the left whereas the meshes generated with our feature preserving extension using PIF are on the right. For both threshold values the extracted polylines and their endpoints are tesselated explicitly.

Figure 7 :

 7 Figure 7: Feature preserving mesh generation on synthetic data sets. From left to right: Our extracted polylines overlayed to the input point cloud; Delaunay refinement mesh of the Poisson implicit surface; output of our feature preserving mesh generation. Top: octa-flower; middle: carved_object (courtesy of P. Jenke); and bottom: blade models recover sharp edges from sparse samplings (≤ 50k points). Notice, that in the case of carved_object and blade, the original models are 3D meshes, but the vertices of those meshes are NOT lying on the feature lines. The initial triangle edges are zigzaggy along the feature lines.

Figure 8 :

 8 Figure 8: Quantitative evaluation of the accuracy versus mesh size : (a) the fandisk model, points on the curve correpond to values 10, 7, 4.5 and 3.5 of δ (b) the block model, points on the curve correpond to values 9, 6, 3 and 1.5 of δ. The Hausdorff distance plotted vertically is given in units of 1 200 of the bounding box diagonal length.

Figure 9 :

 9 Figure 9: From left to right: input point cloud and overlayed extracted feature polylines for the block model and three generated meshes from the implicit surface obtained by Poisson reconstruction: marching cubes implementation[START_REF] Kazhdan | Poisson Surface Reconstruction[END_REF], with default octree depth, (62k facets); standard Delaunay refinement[START_REF] Boissonnat | Provably good sampling and meshing of surfaces[END_REF], with δ = 0.5, (60k facets); our feature preserving Delaunay refinement, with δ = 6, (900 facets).

Figure 10 :

 10 Figure 10: Feature extraction and preserving mesh generation for the fandisk model corrupted with uniform noise of 1% of its bounding box. From top to bottom, left to right: the detected edge points (t 1 = 50,t 2 = 6); the associated feature directions; the extracted feature polylines; the output mesh with (δ = 7). Note how choosing large protective balls patches the small gaps (pointed by small arrows) between feature polylines.

Figure 11 :

 11 Figure 11: Left: input point cloud of Ramses model downsampled to about 27% of the original dataset (top), extracted feature polylines (bottom). Middle: meshed Poisson implicit function. Right: output of our feature preserving mesh generation (5k facets).

Figure 12 :

 12 Figure 12: From top to bottom, left to right: input point cloud of Church model (courtesy of INPG) downsampled to about 7.2% of the original dataset; extracted feature polylines; meshed Poisson implicit function; output of our feature preserving mesh generation.

Table 1 :

 1 Statistics of our prototype implementation.

	(2)	(3)	(4)	Fig.

submitted to Eurographics Symposium on Geometry

Processing (2010)

submitted to Eurographics Symposium on GeometryProcessing (2010)

Acknowledgements

We wish to express our thanks to the reviewers for their insightful comments. We also thank Pierre Alliez for the precious discussions, Stephane Tayeb and Laurent Rineau to provide access to their experimental meshing code, Joel Daniels and the Aim@Shape repository for the several data-sets. This work was partially supported by the ANR (Agence Nationale de la Recherche) under the "Gyroviz" project (No ANR-07-AM-013) www.gyroviz.org and by the Fondation d'entreprise EADS, contract no 3610.