Pedro Henrique Penna

Márcio Castro

Patricia Plentz

Henrique C Freitas

François Broquedis

Jean-François Méhaut

BinLPT: A Novel Workload-Aware Loop Scheduler for Irregular Parallel Loops

Workload-aware loop schedulers were introduced to deliver better performance than classical strategies, but they present limitations on workload estimation, chunk scheduling and integrability with applications. Targeting these challenges, in this work we propose a novel workload-aware loop scheduler that is called BinLPT and it is based on three features. First, it relies on some user-supplied estimation of the workload of the target parallel loop. Second, BinLPT uses a greedy bin packing heuristic to adaptively partition the iteration space in several chunks. The maximum number of chunks to be produced is a parameter that may be fine-tuned. Third, it schedules chunks of iterations using a hybrid scheme based on the LPT rule and on-demand scheduling. We integrated BinLPT in OpenMP, and we evaluated its performance in a large-scale NUMA machine using a synthetic kernel and 3D N-Body Simulations. Our results revealed that BinLPT improves performance over OpenMP's strategies by up to 45.13% and 37.15% in the synthetic and application kernels, respectively.

Introduction

Evenly distributing the workload among the working threads of an irregular application is a NP-Hard minimization problem known as the Multiprocessor Scheduling Problem. This is a challenge to both academic and industry communities, and it is a recurring subject of research in High Performance Computing (HPC). In shared memory parallel applications, this problem emerges when scheduling iterations of parallel loops [START_REF] Polychronopoulos | Guided Self-Scheduling: A Practical Scheduling Scheme for Parallel Supercomputers[END_REF]. In this scenario, the problem is referred as the Loop Scheduling Problem (LSP), and it can be reduced to the assignment of independent loop iterations such that their load is evenly distributed and the scheduling overhead is minimized.

Several loop scheduling strategies have been proposed to address the previous problem [START_REF] Hurson | Parallelization of DOALL and DOACROSS Loops -A Survey[END_REF], and they mainly rely on two techniques. In the first, called on-demand scheduling, iterations are scheduled to threads on-the-fly at runtime, so that both load imbalance and runtime variations may be dynamically handled. In the second technique, called chunk-size tuning, iterations are scheduled in optimally sized batches (i.e. chunks) so that (i) scheduling overheads are mitigated, (ii) load imbalance is further amortized and (iii) iteration affinity is exploited. When coupled together, on-demand scheduling and chunk-size tuning may indeed deliver reasonable performance to a wide range of scenarios. Nevertheless, these techniques do not consider any knowledge about the underlying workload of the target parallel loop, and thus scheduling strategies built upon them naturally turn out to be suboptimal [START_REF] Penna | Design Methodology for Workload-Aware Loop Scheduling Strategies Based on Genetic Algorithm and Simulation[END_REF]]. To address this weakness, workload-aware strategies were introduced [START_REF] Banicescu | Performance of scheduling scientific applications with adaptive weighted factoring[END_REF][START_REF] Kejariwal | History-Aware Self-Scheduling[END_REF][START_REF] Wang | Knowledge-Based Adaptive Self-Scheduling[END_REF]]. These strategies rely on some knowledge about the workload to adaptively fine-tune chunk sizes to further amortize load imbalance, and thus deliver superior performance. Although these strategies present better performance gains than workload-unaware strategies (or blind strategies), existing workload-aware loop scheduling strategies still face some drawbacks that should be addressed.

First, these strategies rely on profiling and statistical regression techniques, and thus are inherently designed to well-behaved workloads. To tackle irregular loops, whose workload varies drastically, some alternatives for estimating the workload on-the-fly are necessary. Usually, the loop scheduling strategy itself and the workload-estimation technique should be loosely coupled. This way, HPC engineers may plug into their solutions the workload-estimator that best fits their needs. However, existing knowledge-based strategies do not provide this flexibility. Moreover, existing workload-aware loop scheduling strategies fail to apply their knowledge about the underlying workload of the target irregular loop when scheduling chunks of iterations. They only rely on-demand scheduling technique, which leads to scalability problems [START_REF] Bull | Feedback Guided Dynamic Loop Scheduling: Algorithms and Experiments[END_REF]]. Finally, no workload-aware strategy is integrated in a publicly available parallel Application Programming Interface (API) or library, since such integration is usually not trivial [START_REF] Banicescu | On the Scalability of Dynamic Scheduling Scientific Applications with Adaptive Weighted Factoring[END_REF]].

In this context, the main goal of this work is to propose a novel workload-aware loop scheduling strategy for irregular parallel loops in which iterations are independent from one another. This new strategy overcomes the aforementioned weaknesses in workload prediction and chunk-scheduling and was integrated in a widely-used API for parallel programming. In summary, this work delivers the following main contributions to the state-of-the-art:

• A novel workload-aware loop scheduling strategy entitled BinLPT. To enable superior performance and flexibility, our strategy is based on three features: (i) a user-supplied estimation about the workload of the target irregular loop; (ii) a greedy bin packing heuristic to adaptively partition the iteration space into several chunks; and (iii) a hybrid scheme based on the Longest Processing Time (LPT) rule and on-demand scheduling [START_REF] Graham | Bounds on Multiprocessing Timing Anomalies[END_REF]]. • An integration of our novel workload-aware loop scheduling strategy into the OpenMP runtime system of GCC.

Our implementation is open-source and publicly available for download, thereby enabling any parallel application that relies on this programming abstraction to seamlessly use our strategy.

To evaluate BinLPT, we carried out a performance analysis using a synthetic kernel and a 3D N-Body Simulations application kernel. We ran experiments on a large-scale Non-Uniform Memory Access (NUMA) machine and we contrasted the performance of BinLPT against the default strategies that are shipped with OpenMP. The remainder of this work is organized as follows. In Section 2, we introduce the LSP and classical loop scheduling strategies. In Section 3, we detail the workload-aware loop scheduling strategy proposed in this work. In Section 4, we present our evaluation methodology. In Section 5, we discuss the experimental results. In Section 6, we contrast related work with ours. In Section 7, we draw the conclusions of our work and discuss some future works.

Background

In this section we introduce the background on which this work relies. First, we present the LSP, and then we discuss the classical loop schedulers.

The Loop Scheduling Problem

The LSP is an instance of the NP-Hard minimization problem for multiprocessor scheduling and it down to partitioning x into n non-overlapping chunks and assigning disjoint sets of these n chunks to the p threads; such that the following are minimized: (i) the number of n chunks that are used to partition the loop iteration space x; and the maximum load imbalance between any pair of threads in p. This formulation shows the relation of the four core variables of the LSP: (i) the loop iteration space x; (ii) the load of iterations w k ; (iii) the number of chunks n in which x will be partitioned; and (iv) the number of threads p that will process the n chunks in parallel.

Additional variables may also play an important role in a real-world context. For instance, the scheduling overhead is an important concern for strategies that assign chunks of iterations to idle threads on-the-fly. If contention in synchronization structures is costly, the irregular parallel loop may face scalability issues [START_REF] Fang | Dynamic Processor Self-Scheduling for General Parallel Nested Loops[END_REF]]. The performance of memory-intensive irregular loops can also be severely impacted by data locality, so memory affinity can be exploited by loop schedulers when through the use of large chunks [START_REF] Markatos | Using Processor Affinity in Loop Scheduling on Shared-Memory Multiprocessors[END_REF].

Loop Scheduling Strategies

Loop scheduling strategies boil down to one of the following two approaches: static, in which loop iterations are assigned to the threads of the parallel application at compiletime; and runtime, in which scheduling decisions are made at runtime [START_REF] Kejariwal | History-Aware Self-Scheduling[END_REF]]. Static scheduling strategies introduce no runtime overhead, but (i) they are only possible on parallel loops which can have their bounds somehow determined at compile time, and (ii) they are suitable only for parallel loops which feature a compile-time predictable workload. In contrast, runtime strategies are employed to address parallel loops that either do not meet the aforementioned compile-time requirement, or perform computation on a workload that is known only at runtime. In this work, we address the problem of scheduling irregular parallel loops in which the workload is known only at runtime, and in the following paragraphs we discuss the classical scheduling strategies on this scenario.

Pure Dynamic Scheduling (PDS) assigns iterations to threads in unit-sized chunks and on-demand. Whenever a thread becomes idle, an iteration is assigned to it. This strategy achieves good load balancing but at the price of a possibly high runtime overhead.

Chunk Self-Scheduling (CSS) works like PDS, but instead of assigning iterations one by one, it assigns iterations in equally-sized chunks [START_REF] Fang | Dynamic Processor Self-Scheduling for General Parallel Nested Loops[END_REF]]. Small chunk sizes deliver good load balancing, but they likely introduce prohibitive runtime overheads. In contrast, large chunk sizes avoid this problem, but may increase load imbalance. When the chunk size is fine-tuned, near-optimum load balancing is achieved [START_REF] Balasubramaniam | Towards the Scalability of Dynamic Loop Scheduling Techniques via Discrete Event Simulation[END_REF], and when it equals to one, this scheduling strategy degenerates to PDS.

Guided Self-Scheduling (GSS) also assigns chunks of iterations to threads on demand, but it dynamically changes their size at runtime (the size of the next chunk is given by the number of remaining iterations divided by the number of threads) [START_REF] Polychronopoulos | Guided Self-Scheduling: A Practical Scheduling Scheme for Parallel Supercomputers[END_REF]. The idea of having a decreasing chunk size is to offer a compromise between achieving good load balancing while reducing runtime overhead.

Factoring Self-Scheduling (FSS) works similarly to GSS, but it differs in the way that chunk sizes are determined [START_REF] Hummel | Factoring: a method for scheduling parallel loops[END_REF]]. To address the scenarios in which GSS does not perform so well, FSS computes the next chunk size by dividing a subset of the remaining loop iterations (usually half) evenly among the threads. FSS introduces no significant runtime overhead compared to GSS, and it may deliver better performance.

Among the aforementioned loop scheduling strategies, OpenMP offers builtin support for FSS and CSS. The OpenMP community pragmatically refers to them as Guided and Dynamic, respectively. Therefore, we will refer to these strategies using the latter notation, unless otherwise stated.

The BinLPT Loop Scheduler

In this section, we present our novel workload-aware loop scheduling strategy. First, we discuss the internals of BinLPT, and then we detail our strategy algorithmically. We implemented BinLPT in libGOMP and we made the enhanced version of this runtime system publicly available 1 under the GPL v3 License. Therefore, any parallel application that is built on top of OpenMP may seamlessly use our scheduling strategy.

Strategy Internals

BinLPT operates in two phases to deliver load balance to irregular parallel loops, namely chunk partitioning and chunk scheduling. The heuristics used in each phase are indeed the key features that enable the superior performance of BinLPT.

In the chunk partitioning phase, BinLPT splits the iteration space into chunks so as to amortize load imbalance while minimizing the number of chunks that are produced. In this way, runtime scheduling overheads can be reduced and iteration affinity may be exploited efficiently. Indeed, this sub-problem could be optimally solved in pseudo-polynomial time using a dynamic programming algorithm for the Linear Partition Problem. Nevertheless, since loop ranges may grow asymptotically, the overhead incurred by this algorithm makes its use prohibitive. Therefore, BinLPT relies on a workloadaware adaptive technique that takes as input a user-supplied threshold k and works as follows. First it computes the average load ω avg for a chunk based on the workload information and k. Next, it uses a greedy bin packing heuristic that bundles into a single chunk the maximum number of iterations whose overall load does not exceed ω avg .

In the chunk scheduling phase, the goal is to come up with a chunk/thread assignment that minimizes load imbalance. Therefore, BinLPT relies on a hybrid scheduling scheme that works as follows. Initially, chunks are statically scheduled to threads using the LPT rule: which assigns the heaviest chunks to the least overloaded threads, and then iteratively assigns lighter chunks to more heavily loaded threads. Next, threads are unblocked and start computing. Then, whenever a thread finishes computing all its chunks, it steals a chunk from other thread that still has work left. This simple scheme optimally handles load imbalance created by both predictable and unpredictable phenomena. Static scheduling based on LPT ensures a 4/3-approximation scheduling solution to load imbalance incurred by the workload. On the other hand, on-demand scheduling ensures that unpredictable phenomena, such as communication latencies, external load interference, and poor workload estimation, are optimally tackled in a 2-approximative fashion [START_REF] Graham | Bounds on Multiprocessing Timing Anomalies[END_REF].

Algorithm 1 BinLPT loop scheduling strategy.

1: function BINLPT(A, k, n) 2: C ← COMPUTE-CHUNKS(A, k) 3:
SORT(C, descending order) 4: for i from 0 to n do 5:

T i ← 0 6: for i from 0 to | C | do 7:
T j ← min T 8:

P Tj ← P Tj ∪ { c i } 9: T j ← T j + ω( c i ) 10: return P 11: function COMPUTE-CHUNKS(A, k) 12: j ← 0 13: C ← empty multiset 14: c 0 ← empty sequence 15: ω avg ← ij ∈A w j k 16: for i from 0 to | A | do 17: c j ← ( c j , A i ) 18: if ω( c j ) > ω avg then 19: C ← C ∪ { c j } 20: j ← j + 1 21: return C

Strategy Design

The BinLPT loop scheduling strategy is outlined in Algorithm 1. In the pictured notation, means that iteration A i is added to c j . It takes as input three parameters: an array that gives a load estimation of each iteration in the target parallel loop (A), the maximum number of chunks to generate (k) and the number of working threads (n). Then it returns a multiset (P ) that states the thread/chunk assignment (i.e. P j is a set containing all chunks assigned to thread j). The algorithm starts by computing chunks according to the greedy bin packing heuristic detailed in the previous section (line 13). Then, it sorts the produced chunks according to their loads (line 14). Next, chunks are statically scheduled following the LPT rule (lines 15 to 20). Later, during the execution, whenever a thread becomes idle, it steals chunks from other threads.

Evaluation methodology

To evaluate BinLPT, we employed a synthetic kernel and an application kernel. The former was used to assess the upper bound performance of our strategy, whereas the latter was employed to enable such analysis in a realistic scenario.

The synthetic kernel performs embarrassingly parallel computations on private variables as proposed in [START_REF] Bull | Feedback Guided Dynamic Loop Scheduling: Algorithms and Experiments[END_REF]], and thus it benchmarks the load balancing performance of a scheduler. The kernel takes as input four parameters: the iteration space size, the input workload w, the scheduling strategy, and the number of working threads. The application kernel performs N-Body Simulations and it was chosen for two main reasons. First, it has great importance to the scientific community as a whole, since it finds application on different domains, such as Computation Fluid Dynamics and Molecular Dynamics [START_REF] Springel | Simulations of the formation, evolution and clustering of galaxies and quasars[END_REF]]. Second, it is a typical irregular kernel that is frequently studied within the context of loop scheduling [START_REF] Banicescu | On the Scalability of Dynamic Scheduling Scientific Applications with Adaptive Weighted Factoring[END_REF][START_REF] Wang | Knowledge-Based Adaptive Self-Scheduling[END_REF]]. The N-Body Simulations kernel that we considered (code-named LavaMD) was extracted from the Rodinia Benchmarks Suite [START_REF] Che | Rodinia: A benchmark suite for heterogeneous computing[END_REF].

Table 1 summarizes the parameters we used in each set of experiments. The workloads are frequently studied by related works, and they were generated using the tool proposed by [START_REF] Penna | Design Methodology for Workload-Aware Loop Scheduling Strategies Based on Genetic Algorithm and Simulation[END_REF]]. Problem sizes were chosen so as to reflect the full processing capacity of the experimental platform. Baseline strategies were selected to be consistent with related works. Chunk sizes were selected based on earlier experiments that revealed them to be the optimal values.

For the synthetic kernel, we adopted a full factorial experimental design, thereby resulting in 75 different scenarios. For these experiments, we set the number of threads to 192 to reflect the full computational power of the experimental platform. For the application kernel, on the other hand, we adopted a fractioned experimental design, where we varied the number of threads from 24 to 192, with a constant step of 24. We considered 27 different scenarios for each experimental configuration. We carried out five replications of each configuration to account for the inherent variance of the measures. For each replicate, the actual order in which individual runs were executed was randomly determined. This approach ensures that experimental results and errors are independent and identically distributed (i.i.d) random variables. In our experiments, the maximum relative standard deviation error (σ/µ) observed was below to 1.0%.

We considered the following performance metrics 2 [Luke et al. 1998, Cariño and[START_REF] Cariño | Dynamic load balancing with adaptive factoring methods in scientific applications[END_REF]]: (a) Parallel Time, which is the overall execution time of the parallel loop; (b) Cost, which is the aggregate time spent to execute the parallel loop, and thus quantifies the waste of processor time; (c) Performance, which is the ratio of the total amount of work to the parallel time; (d) Coefficient of Variance (C.o.V), which is the ratio between the standard deviation and the mean execution time of the threads; (e) Slowdown, which is the ratio between the execution time of the slowest thread to the fastest one.

All experiments were performed on a SGI Altix UV 2000 machine, which has 24 cache coherent NUMA nodes interconnected through SGI's NUMAlink6 (bidirectional). Each node has an Intel Xeon E5 Sandy-Bridge processor and 32 GB of DDR3 memory. Overall, this platform features 192 physical cores and 768 GB of memory. In our experiments, hyper-threading was disabled and we used a first-touch memory allocation strategy coupled with a compact thread affinity policy to mitigate runtime NUMA effects.

Experimental Results

First, we show results for the synthetic kernel, and then we analyze results for the application kernel. All experimental results are publicly available for download 3 .

Synthetic Kernel Results

The actual number of chunks produced by BinLPT varies according to the workload itself and it be fine-tuned by its k parameter. For instance, BinLPT,288 produces at most k = 2 We will refer hereafter to these metrics using capitalization along with their symbol. 3 Experimental results available at: https://doi.org/10.6084/m9.figshare.4742272 288 chunks, pragmatically. On the other hand, for the Guided and Dynamic strategies, the number of chunks that are generated depends on the iteration space | x|. For the former strategy, the number of chunks grows proportionally to O(log | x|), whereas for the latter strategy it grows with O(| x|). In both latter strategies, the granularity of the chunk-sizes may be fine-tuned according to a parameter b. For instance, Dynamic,1 will cause Dynamic to use unit-sized chunks (b = 1). On the other hand, Guided,2 instructs Guided to generate chunks which are not smaller than 2 (b = 2).

Figure 1(a) presents the number of chunks generated by each strategy for an Exponential-generated workload, when varying the size of the iteration space (i.e. | x|) at a constant ratio (2×). We observed similar behaviors for the other workloads (i.e. Gamma-and Gaussian-based), and thus we omitted them due to space limitations. Overall, the results show that the number of chunks produced by BinLPT are far fewer than the ones produced by both Guided and Dynamic, regardless of the values assigned to parameters k and b. The number of chunks generated by Guided grows linearly and for Dynamic grows exponentially. Since the scheduling overhead of runtime on-demand scheduling techniques depends on the number of chunks [START_REF] Cariño | Dynamic load balancing with adaptive factoring methods in scientific applications[END_REF], it turns out that BinLPT is the most scalable scheduling strategy for the chosen k values. Indeed, if we used larger values for k, the number of chunks produced by BinLPT could be bigger. However, in the performance analysis that follows, we reveal that the values used for k are sufficient for BinLPT to deliver superior performance.

Figure 1(b) presents Performance (ρ) results for a Exponential-generated workload, when varying the size of the iteration space (i.e. | x|) at the same constant ratio of 2×. Overall, BinLPT achieved the best results. The highest Performance (ρ) observed for BinLPT,288 was for the scenario with 768 iterations, where it delivered 34.75% superior Performance (ρ) than the best configuration for the other two strategies (Dynamic,1). The highest Performance (ρ) observed for BinLPT,384 was also in the same scenario, where it delivered 37.65% superior Performance (ρ).

The worst Performance (ρ) observed for BinLPT was for the scenario with 6144 loop iterations and k = 288, where we noted 6.09% of Performance (ρ) degradation. Although BinLPT,288 did perform slightly worse, the weight of chunks generated by BinLPT was not fine-grained enough to amortize load imbalance in this scenario. This behavior is confirmed by Figure 1(a), which shows that BinLPT generated much less chunks than the other strategies for this scenario (BinLPT,288 generated 266 chunks, whereas Guided,1, Dynamic,1 and Dynamic,2 generated 4.33×, 23.09× and 11.54× more chunks, respectively). Thus, a fair comparison between BinLPT and the other strategies accounts for the equivalence between the maximum number of chunks k produced by our strategy and the number of chunks generated by them. Recall that BinLPT,288 splits the iteration space in at most 288 variable-size chunks, whereas Dynamic,1 produces 6144 unit-sized chunks (b = 1). Therefore, for instance, a fair comparison would be in the scenario with 768 iterations, between BinLPT,384 (at most 384 chunks) against Dynamic,2 (768/2 = 384 chunks). In this case, BinLPT would deliver 45.13% better Performance (ρ). We observed similar results for the Gamma-and Gaussian generated workloads (they were omitted due to space limitations), in which BinLPT delivered 29.94% and 32.81% better Performance (ρ) over dynamic,2, respectively.

Application Kernel Results

Figure 2(a) presents Parallel Time (τ ) results for the N-Body Simulations application when using 192 threads and varying both the chunk size and the workload. In this plot, the chunk size for the BinLPT means that the parameter k of our strategy was chosen so as it would lead a fair-comparison with Guided and Dynamic in each scenario. Overall, the results unveiled that BinLPT delivers better Parallel Time (τ ) regardless the scenario.

Figure 2(b) presents Cost (γ) results when varying the number of threads, for a scenario with the following configurations (worst-case scenario for BinLPT): Guided,3 Dynamic,3 and BinLPT,384. When analyzing the results, we observed that BinLPT delivers near constant scalability in a large-scale NUMA machine, for all the three workloads. On the other hand, the Cost (γ) scales up quasi-linear for Guided and linear for Dynamic. Nevertheless, it is worthy to note that Cost (γ) results suggest that BinLPT performs better on the Exponential and Gaussian results. Indeed, we confirmed this finding with Performance (ρ) results showed in Figure 2(c). In the case of Exponential-and Gaussian-generated workloads, BinLPT delivers up to 37.15% (168 threads) and 34.45% (192 threads) better Performance (ρ), when considering Dynamic as baseline. In, contrast, for the Gamma workload, we noted 30% of improvement.

So far, we considered metrics that rely on the performance of the slowest thread executing the parallel loop. Thus, we now analyze the other two metrics that consider the performance of all threads: C.o.V (λ) and Slowdown (χ). When we analyzed C.o.V (λ), we found out that BinLPT presented lower values than Guided and Dynamic on the Exponential-generated workload; and it showed up worst C.o.V (λ) results for the other two distributions. Indeed, these results at glance suggest that BinLPT does not deliver load balancing. Nevertheless, since our strategy did present better overall Cost (γ) and Performance (ρ) results than Guided and Dynamic in all the three workloads regardless the number of threads, we drawn the following conclusion. Even though there is a great difference between the execution times of threads in BinLPT, the actual absolute difference is smaller than it is for Guided and Dynamic. We confirmed this finding by plotting histograms of execution times of threads, when running with 96 to 192 threads. However, we had to omit these plots due to space limitations. This conclusion is also pictured in Figure 3(b), which presents Slowdown (χ) results. This plot evidences that the difference between the slowest and fastest threads in BinLPT is actually smaller than in Guided and Dynamic for all thread configurations. Considering Dynamic as baseline, and when using 192 threads, gains are 45.54%, 26.32% and 12.62% for the Exponential-Gamma-and Gaussian-generated workloads, respectively.

Related Work

Targeting time-step applications with irregular parallel loops, Banicescu [START_REF] Banicescu | On the Scalability of Dynamic Scheduling Scientific Applications with Adaptive Weighted Factoring[END_REF] proposed Adaptive Weighted Factoring (AWF). In this strategy, the chunk size is dynamically adapted after each step in the application. The newly computed chunk size decreases across the iteration space, likewise in FSS. However, the performance of the threads during the last time-step and their accumulative performance during all the previous ones is additionally considered in this adjustment. To evaluate the performance of AWF, two inhouse applications were studied: (i) Laplace's Equation Solver on an unstructured grid using Jacobi's method; and (ii) N-Body Simulations. Pure Static Scheduling (PSS) and FSS strategies were considered as baselines. Experiments were carried out on a syntheticallyloaded homogeneous cluster, and the results unveiled that AWF may achieve up to 46% better performance than the baseline strategies. Due to the notable performance of AWF, extensions have been proposed to enable its use on non-iterative applications as well [START_REF] Cariño | Dynamic load balancing with adaptive factoring methods in scientific applications[END_REF]. Nevertheless, the enhanced version of this strategy presented a performance that is comparable to the one achieved by FSS.

To address a broader class of applications, [START_REF] Kejariwal | History-Aware Self-Scheduling[END_REF] proposed History-Aware Self-Scheduling (HSS). Unlike AWF, HSS relies on statistical information collected offline via profiling to carry out a smarter scheduling. Based on this extra knowledge, at every scheduling round, HSS computes chunk sizes in a decreasing fashion like FSS, but also considering the load of previously executed iterations and their corresponding actual loads. To assess the performance of HSS, irregular parallel loops extracted from the Standard Performance Evaluation Corporation (SPEC) Benchmarks were studied, and the FSS and AWF strategies were considered as baselines. Experiments were carried out on an in-house simulator, and the results unveiled that HSS may outperform baseline strategies up to 18%.

Based on a similar offline profiling-guided approach to HSS, Wang et al. [START_REF] Wang | Knowledge-Based Adaptive Self-Scheduling[END_REF]] introduced Knowledge-Based Adaptive Self-Scheduling (KASS). This strategy works on two phases: a static partitioning phase, and a dynamic scheduling phase. In the first phase, a knowledge-based approach is used to partition iterations of the parallel loop into local work queues of threads, which makes the total workload to be equally distributed to the threads, approximately. In the second phase, iterations on local work queues are scheduling with decreasing sizes likewise in FSS. Each thread gets a chunk from its local queue to execute, and when it finishes the execution of all the chunks in its local queue, it steals chunks from other threads. To evaluate the performance of KASS, two scenarios were studied: (i) parallel loops extracted from the SPEC Benchmarks; and (ii) and three in-house application kernels, namely Over-Relaxation, Jacobi Iteration and Transitive Closure. The GSS, FSS, Trapezoid Self-Scheduling (TSS) and Affinity Self-Scheduling (AFS) strategies were considered as baselines. Experiments were carried out on a Symmetric Multiprocessing (SMP) machine, and the results unveiled that for the parallel loops, KASS is up to 16.9% faster than the baseline strategies. On the other hand, for application kernels, KASS achieved up to 21% better performance.

The workload-aware loop scheduling strategy that we proposed in this work differs from the related strategies discussed above in several points, thereby delivering contributions to the state-of-the-art. First, unlike all these strategies, BinLPT does not rely on a particular workload-estimation technique. Indeed, the HPC engineer is free to couple BinLPT with the one that yields to the best workload estimation for the application. Consequently, the applicability of our strategy is not restricted to time-step applications such as AWF nor to applications that present well-behaved workloads like HSS and KASS, which rely on offline profiling and online regression techniques. Second, even though the aforementioned strategies do use their estimations on the workload to partition the iteration space in several chunks, they lack in using this knowledge to actually schedule chunks of iterations. Alternatively, BinLPT uses a hybrid scheduling scheme based on the LPT rule and on the on-demand scheduling technique. The former handles workload imbalance in a 4/3-approximative fashion, while the latter deals with unpredictable phenomena in a 2-approximation optimally [START_REF] Graham | Bounds on Multiprocessing Timing Anomalies[END_REF]]. Finally, existing strategies lack on integrability with applications. For instance, their source-code is not available for download, and their reported algorithmic description is not detailed enough to enable an in-house implementation of them. Our solution, on the other hand, is open-source and is built into GCC's OpenMP runtime.

Conclusions

In this work we proposed a novel workload-aware loop scheduling strategy called BinLPT. To enable superior performance and flexibility, our strategy is based on three features. First, it relies on some user-supplied estimation of the workload of the target irregular loop. Such estimation may be derived either from the problem structure or through online/offline profiling, thus enabling maximum flexibility. Second, BinLPT uses a greedy bin packing heuristic to adaptively partition the iteration space in several chunks. The number of chunks to be produced is a parameter of our strategy that may be fine-tuned. Third, it schedules chunks of iterations using a hybrid scheme based on the LPT rule and on-demand scheduling.

We integrated BinLPT into OpenMP, and we evaluated its performance in a largescale NUMA machine using a synthetic kernel and a 3D N-Body Simulations application kernel. We considered several workloads and we contrasted the performance of our strategy against other strategies available in OpenMP (Guided and Dynamic). Our results unveiled that BinLPT achieves up to 45.13% and 37.15% better performance in the synthetic and application kernels, respectively. As future work, we intend to enhance BinLPT so that it also accounts for data locality when scheduling chunks; and to extend our strategy to emerging manycore platforms which feature a distributed shared memory.

Figure 1 .

 1 Figure 1. Synthetic kernel benchmarking results for Exponential workload.

  Time (τ ) vs. chunk size when running with 192 threads.

Figure 2 .

 2 Figure 2. Experimental results for the N-Body Simulations application kernel.

Figure 3 .

 3 Figure 3. Experimental results for the N-Body Simulations application kernel.

Table 1 . Parameters for experiments.

 1

	Parameters	Synthetic Kernel	Application Kernel
	Workload PDF Exponential, Gamma and Gaussian Exponential, Gamma and Gaussian
	Loop Size	{384, 768, 1536, 3072, 6144}	11 × 11 × 11
	Chunk Size	Guided {1}, Dynamic {1,2}, BinLPT {288, 384}	Guided {1,2,3}, Dynamic {1,2,3}, BinLPT {384, 576, 768}

www.github.com/lapesd/libgomp

Acknowledgment

This work was supported by FAPEMIG, FAPESC, CAPES, CNPq, CNRS and INRIA.


Dobrina Boltcheva
email: dobrina.botlcheva@univ-lorraine.fr

Bruno Lévy

Surface reconstruction by computing restricted Voronoi cells in parallel

Keywords: surface reconstruction, point cloud, restricted Voronoi diagram

We present a method for reconstructing a 3D surface triangulation from an input point set. The main component of the method is an algorithm that computes the restricted Voronoi diagram. In our specific case, it corresponds to the intersection between the 3D Voronoi diagram of the input points and a set of disks centered at the points and orthogonal to the estimated normal directions. The method does not require coherent normal orientations (just directions). Our algorithm is based on a property of the restricted Voronoi cells that leads to an embarrassingly parallel implementation. We experimented our algorithm with scanned point sets with up to 100 million vertices that were processed within few minutes on a standard computer. The complete implementation is provided.

Introduction

This article deals with surface reconstruction, a problem that may be formulated as follows: given a set of points sampled from a surface, recover the original surface from which those points were sampled. This general problem is motivated by numerous applications in reverseengineering, prototyping, visualization, or computer vision since a growing variety of scanning devices nowadays provides measurements of objects in the form of point sets.

In practice, the common surface reconstruction pipeline usually involves several steps: After the acquisition with a laser scanner, for instance, the scans from the different views are first aligned into a common global coordinate system. Then they are filtered to reduce the input noise and remove some outliers. Then the surface normals at the points are estimated (and possibly oriented) and finally, a surface reconstruction algorithm builds a triangular mesh approximating the underlying surface.

Contribution: We propose a reconstruction method that is intended to be a part of the global reconstruction pipeline. It is designed to take as input an already filtered point set and connect the input points with triangles by computing their restricted Voronoi diagram.

The algorithm is similar in spirit to Localized Cocone [START_REF] Dey | Localized cocone surface reconstruction[END_REF] and Tangential Delaunay Complex [START_REF] Boissonnat | Manifold reconstruction using tangential Delaunay complexes[END_REF]. As in these local algorithms, it does not compute any 3D data structure like a grid or a triangulation of the ambient space. The main difference is that our algorithm computes the intersection between the 3D Voronoi diagram of the points and a set of disks centered at the points and orthogonal to the estimated normal directions. We perform an iterative disk clipping in parallel for every point. As a consequence, our algorithm is both embarrassing parallel and memory efficient. It can handle large point sets, with up to 100 million points, in few minutes on an off-the-shelf PC.

Moreover, our algorithm dose not require consistent normal orientations, like in the Poisson methods [START_REF] Kazhdan | Screened poisson surface reconstruction[END_REF], nor a scale value for every input sample, like in the Floating Scale Surface method [START_REF] Fuhrmann | Floating scale surface reconstruction[END_REF]. It solely relies on estimated normal directions.

In addition, we also propose a manifold extraction algorithm and a set of post-processing heuristics to recover a clean mesh, specifically a 2-pseudo-manifold with or without boundaries [START_REF] Munkres | Elements of algebraic topology[END_REF].

Main limitations: Although our reconstruction algorithm is designed for already filtered point sets (where the points are assumed to be located exactly on the underlying surface), it is also useful in practice for real datasets (produced as unions of range maps where the points create some sort of "solid/thick surface"). For such datasets, the filtering becomes a mandatory pre-processing step which is a difficult issue and is beyond our current scope. Here, we use one or two iterations of a standard smoothing approach (projection onto the average plane of neighbors) to smooth the raw datasets which is often sufficient as confirmed by our experiments, see Figures 13,[START_REF] Manson | Streaming surface reconstruction using wavelets[END_REF]19 and 20. Note however that this smoothing is very moderate and does not wash out the small details, as shown on Fig. 16.

Related work

There are two classes of approaches which differ as to whether they approximate or interpolate the input points.

Approximation methods. They aim at fitting a smooth surface to the input point set. These methods are wellequipped to handle various imperfections in the data such as noise, occlusion and registration errors. A recent stateof-the-art and a benchmark can be found in [START_REF] Berger | A benchmark for surface reconstruction[END_REF][START_REF] Berger | State of the Art in Surface Reconstruction from Point Clouds[END_REF]. Many of these algorithms [START_REF] Hoppe | Surface reconstruction from unorganized points[END_REF][START_REF] Boissonnat | Smooth surface reconstruction via natural neighbour interpolation of distance functions[END_REF] compute a distance field by estimating the tangent plane at every point and computing closest distances using these tangent planes. The method of VRIP [START_REF] Curless | A volumetric method for building complex models from range images[END_REF] takes advantage of the range scans acquired through laser triangulation to construct a volumetric signed distance field which merges the scans in a least-squares sense. The other approaches range from computing an indicator function [START_REF] Kazhdan | Poisson surface reconstruction[END_REF][START_REF] Kazhdan | Screened poisson surface reconstruction[END_REF][START_REF] Manson | Streaming surface reconstruction using wavelets[END_REF][START_REF] Alliez | Voronoibased variational reconstruction of unoriented point sets[END_REF][START_REF] Fuhrmann | Floating scale surface reconstruction[END_REF], to locally fitting functions and moving least-squares methods [START_REF] Alexa | Computing and rendering point set surfaces[END_REF][START_REF] Ohtake | Sparse low-degree implicit surfaces with applications to high quality rendering, feature extraction, and smoothing[END_REF][START_REF] Nagai | Smoothing of partition of unity implicit surfaces for noise robust surface reconstruction[END_REF]. In particular, the Poisson surface reconstruction method [START_REF] Kazhdan | Poisson surface reconstruction[END_REF] solves for an approximate indicator function of the inferred solid, whose gradient best matches the input normals. The resulting scalar function is represented in an adaptive octree and is iso-contoured using a variant of the marching cubes algorithm [START_REF] Lorensen | Marching cubes: A high resolution 3d surface construction algorithm[END_REF]. Later, a Screened Poisson algorithm [START_REF] Kazhdan | Screened poisson surface reconstruction[END_REF] has been introduced which resolves the over-smoothing problem and offers the ability to reconstruct meshes with boundaries. Recently, a Floating Scale Surface reconstruction method [START_REF] Fuhrmann | Floating scale surface reconstruction[END_REF] has been proposed for multi-video-stereo applications that operates on large, redundant and noisy point sets but relies on a scale value for every input point which is estimated from the depth maps. Most of these methods require point sets with consistently oriented normals to estimated the signed-distance function. However, normal estimation in the presence of imperfect data remains an open problem [START_REF] Mitra | Estimating surface normals in noisy point cloud data[END_REF][START_REF] Dey | Normal estimation for point clouds: A comparison study for a voronoi based method[END_REF].

Interpolation methods. They elaborate upon Voronoi-Delaunay concepts and produce an interpolating surface in the form of a triangulation which uses a subset of the input points as vertices [START_REF] Edelsbrunner | Triangulating topological spaces[END_REF]. Often these algorithms extract this triangulation from the Delaunay triangulation of the input point set. The theory of -sampling [START_REF] Amenta | Surface reconstruction by voronoi filtering[END_REF] provides a solid background to study the properties of these algorithms. Extensive research effort has built on this theory, producing several Restricted Delaunay based methods such as Cocone [START_REF] Amenta | A simple algorithm for homeomorphic surface reconstruction[END_REF], Super Cocone [START_REF] Dey | Localized cocone surface reconstruction[END_REF] and Power Crust [START_REF] Amenta | The power crust[END_REF] algorithms. Many other extensions have been compiled in a survey [START_REF] Cazals | Delaunay triangulation based surface reconstruction: Ideas and algorithms[END_REF] and a monograph [START_REF] Dey | Curve and Surface Reconstruction: Algorithms with Mathematical Analysis[END_REF]. Among the Delaunay-based methods, there are the greedy approaches ( [START_REF] Cohen-Steiner | A greedy delaunay-based surface reconstruction algorithm[END_REF], [START_REF] Kuo | A delaunay-based region-growing approach to surface reconstruction from unorganized points[END_REF]), the Natural Neighbors ( [START_REF] Boissonnat | Smooth surface reconstruction via natural neighbour interpolation of distance functions[END_REF]), the "convection" algorithm [START_REF] Allegre | A flexible framework for surface reconstruction from large point sets[END_REF], the sampling theorems introduced in [START_REF] Boissonnat | Provably good sampling and meshing of surfaces[END_REF], and the WRAP algorithm together with its analysis in [START_REF] Edelsbrunner | Surface Reconstruction by Wrapping Finite Sets in Space[END_REF][START_REF] Ramos | Geometric and topological guarantees for the wrap reconstruction algorithm[END_REF].

To avoid the high computational cost associated with computing the 3D Delaunay triangulation of the point set, several local computation methods were proposed. For example, The Ball Pivoting algorithms (BPA) [START_REF] Bernardini | The ball-pivoting algorithm for surface reconstruction[END_REF][START_REF] Angelo | A new mesh-growing algorithm for fast surface reconstruction[END_REF] are local, optimized variants of the Alpha Shapes [START_REF] Edelsbrunner | Three-dimensional alpha shapes[END_REF]. These algorithms can be parallelized to handle large datasets [START_REF] Digne | An Analysis and Implementation of a Parallel Ball Pivoting Algorithm[END_REF]. The method introduced by Digne et al. [START_REF] Digne | Scale space meshing of raw data point sets[END_REF] improves the noise-resilience of the BPA by applying it to a scale-space version of the input point set.

Some local variants of the Cocone algorithm were also developed to mitigate the scalability problem. These algorithms are local in the sense that no 3D data structure like a grid or a triangulation of the ambient space is computed and the final triangulation is obtained by gluing local triangulations around each sample point. Funke and Ramos [START_REF] Funke | Smooth-surface reconstruction in nearlinear time[END_REF] demonstrated that the Cocone algorithm can be modified so that no global computation of the 3D Delaunay triangulation of the entire point set is required. They showed that with a locally uniform -sampling the Cocone of each point is a local object and it is completely determined by nearby samples. Dumitriu et al. [START_REF] Dumitriu | On the locality of extracting a 2-manifold in r3[END_REF] used this framework to design a surface reconstruction algorithm with correctness guarantees. Although the algorithm is quite complicated, they have produced an implementation [START_REF] Dumitriu | How much geometry it takes to reconstruct a 2-manifold in r3[END_REF]. Amenta and Kil [START_REF] Kil | Gpu-assisted surface reconstruction on locally-uniform samples[END_REF] employed similar techniques and designed a parallel algorithm operating on the GPU.

Alternatively, Dey et al. [START_REF] Dey | Localized cocone surface reconstruction[END_REF] developed an octree-based version of the Cocone algorithm such that the 3D Delaunay triangulation is only computed on small clusters of the initial point set. By adopting the assumption that the local sampling density is bounded by a constant [START_REF] Attali | Complexity of the delaunay triangulation of points on surfaces: the smooth case[END_REF], they showed that the theoretical guarantees of the original Cocone algorithm are maintained. However, when the input point set does not meet the sampling requirements, it is necessary to use some heuristics in the manifold extraction step.

Another promising idea to deal with the scalability problem is to restrict the computations to planes [START_REF] Gopi | Surface reconstruction based on lower dimensional localized delaunay triangulation[END_REF][START_REF] Gopi | A fast and efficient projection-based approach for surface reconstruction[END_REF] or to the tangent space [START_REF] Boissonnat | Manifold reconstruction using tangential Delaunay complexes[END_REF]. The tangential complex is obtained by gluing local (Delaunay) triangulations around each sample point and the neighborhoods are "reconciled" with a "Star-Splaying" approach [START_REF] Shewchuk | Star splaying: An algorithm for repairing delaunay triangulations and convex hulls[END_REF]. Although, this has been the first theoretical algorithm able to reconstruct a smooth closed manifold in a time depending only linearly on the dimension of the ambient space, to the best of our knowledge, there is no existing implementation.

Computing the restricted Voronoi diagram in parallel

The general work-flow of our algorithm is summarized in Figure 1. The input is a set of 3D points, with normals or not. Note that no consistent orientation of the normals is required but only directions. The output of the algorithm is a triangular mesh connecting the input points which approximates the surface. Specifically, the output triangulation is a compact 2-pseudo-manifold (see Section 3.5), it is orientable (no Moebius strip) and may contain several connected components and boundaries.

The next sections describe each step in detail, as well as the parameters they depend on. The first three steps, kd-tree construction §3.1, point set smoothing §3.2 and normals estimation §3.3 use classical methods. Our main contribution is the candidate triangles step ( §3.4), that computes the intersection between the Voronoi cells of the input points with a set of disks centered on the input points. It shares some similarities with the theoretical method of Funke and Ramos [START_REF] Funke | Smooth-surface reconstruction in nearlinear time[END_REF], the Localized Cocone algorithm [START_REF] Dey | Localized cocone surface reconstruction[END_REF], the Tangential Delaunay Complex [START_REF] Boissonnat | Manifold reconstruction using tangential Delaunay complexes[END_REF] and the Star Splaying [START_REF] Shewchuk | Star splaying: An algorithm for repairing delaunay triangulations and convex hulls[END_REF]. Instead of computing the Delaunay triangulation then restricting it, we directly compute the restricted Voronoi cells (intersections between the Voronoi cells and the disks), and deduce the restricted Delaunay triangles from the combinatorial information of the restricted Voronoi vertices (see Section 3.4, Figure 2, Figure 4 below). This results in an algorithm which is simple and easily parallelizable.

We also introduce a new manifold mesh extraction method ( §3.5) which avoids the trimming step used in previous work (the most critical because it can end up with an empty mesh) [START_REF] Amenta | A simple algorithm for homeomorphic surface reconstruction[END_REF]. Starting from a subset of candidate triangles which forms a clean and orientable 2-manifold with boundaries, we add the remaining triangles, one by one, in an arbitrary order, only if they do not break the topological properties of the initial mesh. The algorithm terminates with the largest possible non-empty manifold mesh.

Kd-tree construction

The basic geometric operation in our algorithm is finding the nearest neighbors of a point. To optimize this operation, we organize the input points in a kd-tree [START_REF] Samet | Foundations of Multidimensional and Metric Data Structures[END_REF].

To do this, one may use the Approximate Nearest Neighbors (ANN) library [START_REF] Mount | ANN: A library for approximate nearest neighbor searching[END_REF] which, with set to 0, computes the exact nearest neighbors for every point. But the memory consumption can be significantly reduced by replacing the kd-tree used in ANN with a balanced binary tree, with links that are completely implicit. This makes it possible to process large point sets (up to 50 million points) on a desktop PC with 16 Gb RAM. Our implementation is given in the supplemental material.

Smoothing

Since our algorithm interpolates the input points, it requires a clean and accurately registered point set. Thus, if the input set is noisy such as raw scanning data, it has to be smoothed before the reconstruction. Note however that any point cloud smoothing technique can be used at this step.

Our current implementation uses the projection onto the best approximating plane of the k-nearest neighbors, [START_REF] Hoppe | Surface reconstruction from unorganized points[END_REF]. Specifically, for a point p and its k nearest neighbors {p i } k i=1 , we find the fitting plane n T x = c for p by minimizing the error term e(n, c) = k i=1 (n T p i -c) 2 under the constraint n T n = 1. Then, the point p is projected onto this plane. We compute the best fitting plane for each point in parallel. Once all the points are projected, we update the kd-tree.

Parameters. This step uses two parameters: the number of nearest neighbors nb neighbors and the number of smoothing iterations smooth iter. We used 30 neighbors for all datasets, and 1 or 2 smoothing iterations for the raw scanned data.

Normal estimation

If the input points come without normals, our algorithm estimates them locally with the best fitting plane as in the smoothing step. Note however that the algorithm does not need to orient the normals, it only needs a normal direction at every point which is used to create a disk in the subsequent step of the algorithm.

Candidate Triangles

Given a point set {p i } n i=1 with estimated normals {n i } n i=1 and a radius r, our algorithm computes the intersection between the Voronoi cells of the points V (p i ) and the disks D r i of radius r centered on the points and orthogonal to the normals.

During this step, we fist build a disk D r i at every point p i in the input point set P orthogonal to n i and with usergiven radius r. In our implementation, we use a polygon (with 10 vertices) that approximates the disk. In contrast with previous algorithms, we explicitly build the restriction of the Voronoi cell

V (p i ) to the disk D r i at p i , V (p i ) ∩ D r
i , as shown on Fig. 2 and Fig. 4. Thus we collect a set of restricted Voronoi vertices. The set T of candidate triangles is then defined by the triangles which are dual to the restricted Voronoi vertices (more on this below).

Note that the disk radius has to be chosen in function of the sampling rate and it entails a trade-off between accuracy and computation speed. A very small radius leads to uncovered holes in the final mesh (2(A) and (B)), while a very large one increases the computation time (since the algorithm looks for neighbors in a larger space).

The basic operation in this step of the algorithm consists of clipping the disk D r i with the Voronoi cell V (p i ) of every point p i . Here we use ideas that are similar to those presented in the Localized Cocone algorithm [START_REF] Dey | Localized cocone surface reconstruction[END_REF] which allows us to build locally the candidates triangles by using only the k-nearest neighbors. More precisely, we use the clipping algorithm and the local characterization of the restricted Voronoi cells introduced in [START_REF] Levy | Variational Anisotropic Surface Meshing with Voronoi Parallel Linear Enumeration[END_REF] which we recall in the following in order to make the explanation self-contained.

Computing Voronoi cells through iterative clipping

Let us recall that the Voronoi cell V (p i ) corresponds to the intersection of the halfspaces : V (p i ) = j =i Π + (i, j), where Π + (i, j) denotes the halfspace bounded by the bisector of (p i , p j ) that contains p i . Note that the bisectors between p i and all the other points are involved in the definition above, whereas only a small subset of them corresponds to actual Voronoi edges (see Figure 3). We can classify the bisectors Π + (i, j) into two sets: contributing and non-contributing. Clipping a Voronoi cell by a non-contributing bisector does not change the result. Therefore, the intersection of any superset of the contributing bisectors corresponds to the Voronoi cell. Let p j1 , p j2 , . . . p jn-1 denote the vertices sorted by increasing distance from p i . Let V k denote the intersection of the k first halfspaces and let R k denote its p i -centered radius:

V k (p i ) = k l=1 Π + (i, j l ) R k = max{d(p i , x)|x ∈ V k (p i )}.
Data: the index i of the point x i , the disk D r centered on it, the number of neighbors n and the set of points P Result: Restricted Voronoi Facet at p i :

RV

(p i ) = V (p i ) ∩ D r i RV ← D r i R k ← max{d(p, p i )|p ∈ RV } k ← 1 while d(x i , x j k ) < 2R k f and k < n do RV ← RV ∩ Π + (i, j k ) R k ← max{d(p, p i )|p ∈ RV } k ← k + 1 end Algorithm 1:
Computes a restricted Voronoi cell as the intersection between a tangential disk and the Voronoi cell of a point.

The configuration in Figure 3 suggests that for all j such that d(p i , p j ) > 2R k , the bisector Π + (i, j) is noncontributing, i.e. V k (p i ) ⊂ Π + (i, j). This observation can be formally proven as follows:

Consider p ∈ V k (p i ) and p j such that d(p i , p j ) > 2R k . By definition of R k , d(p, p i ) < R k . We have d(p i , p) + d(p, p j ) ≥ (p i , p j ) (trian- gular inequality) and d(p i , p j ) > 2R k , therefore d(p, p j ) > R k > d(p, p i ) and p ∈ Π + (i, j).
As a direct consequence:

d(p i , p j k+1 ) > 2R k ⇒ V k = V (p i ).
We call radius of security the first value of R k that satisfies this condition. Note however that this observation does not have a practical value in the case of the unrestricted Voronoi diagram, since some cells are unbounded and have infinite R k (for an infinite cell, all the bisectors are considered, leading to prohibitive computation time). However, the observation can be clearly used to compute the restricted Voronoi cells RV (p i ) = V (p i ) ∩ D r i , ∀p i ∈ P , since they are finite (thay are contained in D r i ). We use Algorithm 1 to build the restricted Voronoi cell for every point, which is the intersection between the disk D r i , centered on it and orthogonal to the estimated normal and the Voronoi cell of the point. Note that every restricted Voronoi vertex q is the intersection between the disk D r i and the bisectors Π i,j1 and Π i,j2 . We deduce from the indices j 1 and j 2 the Delaunay triangle i, j 1 , j 2 associated with each restricted Voronoi vertex, see Fig. 4. We denote by T the so-defined set of candidate triangles.

Parameters. This step uses one parameter: the disk radius (in % of the bounding box diagonal). Note that when the radius r is larger than the width of the Voronoi cells, the disks D r touch all the sides of their Voronoi cells (like on Figure 2 (C)), and increasing the radius no longer changes the result. Thus, in our experimental results, we set the radius to a large value, r=5% of the bounding box diagonal. For very large datasets (> 10M vertices), since Voronoi cells are smaller, we decrease the size (r=0.5% bbox diagonal) to make computations faster because the algorithm looks for the neighbors within a smaller space.

Structure of the set of candidate triangles

For applications that just need to visualize the surface, providing the set of candidate triangles may be sufficient. However, most of the applications in Geometry Processing require clean manifold meshes with coherent orientations for further processing. This is why we present hereafter a manifold extraction method which uses a set of heuristics to build a valid mesh from the soup of triangles T .

We first need to take a closer look on the structure of T . The set T of candidate triangles is composed of the restricted Delaunay triangles T 3 plus some other triangles, as explained below.

The set T 3 of restricted Delaunay triangles is defined as:

T3 = {(i, j, k) |(D r i ∩ V (pj) ∩ V (p k ) = ∅)and (D r j ∩ V (pi) ∩ V (p k ) = ∅)and (D r k ∩ V (pi) ∩ V (pj) = ∅)}. (1)
In other words, they correspond to triples of restricted Voronoi cells (RV (p i ), RV (p j ), RV (p k ) that are mutually in contact (see Figures 2 and4). These triangles can be easily found, by generating all the index triples that correspond to the restricted Voronoi vertices, sorting them (using, for instance, the lexicographic order) and keeping the triples that appear three times in the sorted sequence. The so-computed T 3 triangles are displayed in Figure 5-A. Note that in the sorted sequence, there are also triples that appear once or twice. The corresponding set of triangles is referred to as T 1,2 and is defined by:

T12 = {(i, j, k) | D r i ∩ V (pj) ∩ V (p k ) = ∅} -T3. (2)
This means that a triangle (i, j, k) is in T 12 whenever the restricted Voronoi cell of p i "sees" p j and p k but not both conversely. Some configurations are shown in Figure 6 that correspond to nearly co-spherical points. Since the restricted Voronoi vertices are computed independantelly for every sample point, the numerical computations lead to vertices whose cooridantes differ slightly du to the finite computer precision. Therfore, there are restricted Voronoi vertices with different cardianlities -the majority are computed for exactly three neighoring samples, but there are also restricted vertices obtained only from 1 or 2 sample points.

The so-defined T 12 triangles are displayed in red in Figure 5-B. As expected, they are not as "reliable" as the T 3 triangles, in the sense that they are likely to generate non-manifold configurations (see for instance the four interconnected red triangles that form a "sliver" in the top right part of the neck in Figure 5-B). However, one can observe that they contain interesting information, and could be used to fill the holes in the T 3 triangles (Figure 5-A). Thus, the idea is to start from the T 3 triangles, then carefully insert the triangles from the T 12 set, one triangle at a time, making sure that the updated mesh remains manifold. If it is not the case, the triangle is rejected. The associated algorithm is detailed in the next section.

Manifold mesh extraction

We first build a clean orientable 2-manifold mesh with a subset of the T 3 triangles. Then we fill the holes iteratively, with the remaining triangles, if they do not jeopardize the topological properties of the output.

The mesh structure that we build incrementally is a manifold mesh where every edge can only be adjacent to two facets, edges with only one incident triangle are allowed, as well as triangles adjacent only by one vertex (this is a non-manifold vertex whose neighborhood is not a topological 2-sphere). The output mesh (also known as a 2-pseudo-manifold [START_REF] Munkres | Elements of algebraic topology[END_REF]) is orientable and may have boundaries and several connected components.

We initialize the output mesh with all T 3 triangles. Then we remove all the triangles incident to non-manifold edges and non-manifold vertices by excess. Recall that a nonmanifold edge is incident to more than 2 triangles, while a 'by-excess' non-manifold vertex has a closed loop of triangles (clean umbrella) in its neighborhood plus additional triangles. At this point, we also orient the initial mesh and check that there are no Moebius strips within it. If the algorithm finds a Moebius strip, it removes the last triangle that closes the loop. Note that the initial mesh can however exhibit triangles incident to each other only by a vertex.

We then fill the holes of the initial mesh with the remaining triangles from T 12 . The T 12 triangles are tested against the following criteria, from the simplest to the most time consuming. Only the triangles that pass all the tests are inserted in the final surface.

• Geometric test: ensures that the normals of the neighboring triangles agree, meaning that they do not make a sharp angle (less than π/3).

• Combinatorial test I: ensures that the triangle is properly connected to the current mesh. Every new triangle should be either incident to two edges of existing triangles, or to one edge of an existing triangles and one isolated point.

• Combinatorial test II: checks for non-manifold edges and tests whether the three candidate edges are manifold.

• Combinatorial test III: checks that inserting the new triangle do not generate 'by-excess' non-manifold vertices.

• Combinatorial test IV: checks if in the neighborhood of the triangle the same connected component appears with two opposite orientations, then connecting the triangle would create a Moebius strip. The triangle is rejected if it is incident to the same connected component with two different orientations.

The output of the algorithm is the largest possible triangular mesh that connects the input points. It is orientable 2-manifold (without non-manifold edges) but it can have non-manifold vertices, such as a pinched torus. It possibly has several connected components and boundaries.

Parameters. This step uses one parameter which is the maximum deviation angle between the triangle normals max N angle. We always use the default value set to π/3.

Post-processing

At this point, after the manifold extraction step, the resulting mesh may still contain a number of isolated holes, caused mainly by insufficiently sampled sharp features (see e.g. the top of the ears of the horse in Figure 11 C). To eliminate these, we add an extra post-processing step, that fills some of the holes.

The algorithm takes an user-defined maximum hole size parameter, in terms of number of edges. Then it looks for simply connected holes and fills them using a classical loop-split algorithm [START_REF] Liepa | Filling holes in meshes[END_REF]. This algorithm splits the hole recursively until it gets holes composed of exactly 3 edges and fills them with the corresponding triangle.

In order to ensure that the holes are simple loops, we first search and remove the configurations where the holes contain "bridges"(see Fig. 7). We call bridge triangles those that are visited more than once when we walk around the border edges of a hole. Note that this definition dose not capture bridges wider than two triangles but we did not encounter any, in practice.

During this last post-processing step, we can also remove some remaining small components, in function of their size (in number of facets).

Parameters. This step uses the following two parameters which give the size of the holes to be filled and the connected components to be removed. In our experiments, we have always used:

• max hole edges (=500) : Fill holes with a smaller number of edges

• min comp facets (=10) : Remove small components (in facet number)

Experimental results and Discussion

This section details our empirical results and gives comparisons with four state-of-the-art approaches implemented by their authors (ScaleSpace reconstruction [START_REF] Digne | An Implementation and Parallelization of the Scale-Space Meshing Algorithm[END_REF], Screened Poisson [START_REF] Kazhdan | Screened poisson surface reconstruction[END_REF], SuperCocone [START_REF] Dey | Localized cocone surface reconstruction[END_REF] and Wavelet reconstruction [START_REF] Manson | Streaming surface reconstruction using wavelets[END_REF]). We evaluated the methods both in an artificial scenario with point sets sampled from a known reference surface ( §4.1) and with a database of models with sizes ranging from a few thousand to 100 million points ( §4.2).

Comparison -Hausdorff distance to reference surface

We first evaluate experimentally whether our method reconstructs the same surface as previous work, under different sampling conditions. We created a point set by sampling a given reference surface and measured the Hausdorff distance between the reference surface and the output of the reconstruction methods. The samplings are shown in Figure 9. We tested our method both without and with post-processing (holes filling). The results are reported in Table 1. Our method gives an output with a deviation from the reference surface similar to the state-of-the-art methods.

Real-scale data

We experimented our algorithm with datasets from public repositories [START_REF]Aim shapes repository[END_REF]53,[START_REF]Farman institute 3d point sets[END_REF][START_REF]Statue model repository[END_REF][START_REF]Smart multimedia[END_REF][START_REF]Large-scale point cloud classification benchmark[END_REF] and the data used in the original articles of the four state-of-the-art methods we compared with. We conducted different experiments, to understand the behavior of our method and the previous work with point sets of different sizes (from a few thousands to 100 million points) and with point sets with different structures (non-uniform sampling, small details, noise, overlapping scans).

The table shown on Figure 8 lists some of the models and reports the timings for the algorithms. We used implementations written by the authors of the cited methods, and the same computer for all the tests [Processor: Intel Core i7-4900MQ CPU @ 2.80GHz, RAM:16Gb]. We measured the timings by using the time function in Linux. Thus, they indicate all computation steps (loading and saving included) for all methods. In particular, they comprise the pre-and post-processing steps of our method. Size-time Table 1: For each algorithm, the table reports the two one-sided distances d(M 1, M 2) and d(M 2, M 1) between M 1 (the reference mesh) and M 2 (the mesh built by the algorithm). Distances are expressed as percentages of the diagonal of the mesh bounding box. curves for all methods and datasets are shown in Figure 10.

To be able to analyze and compare qualitatively the behavior of the algorithm, we deactivated the post-processing step (hole filling). In the following (unless specified otherwise), we compare the output of our manifold extraction step with the outputs of the four methods cited before.

Non-uniform sampling

Figure 11 shows a horse sampled with 100K points with a varying sampling density. ScreenedPoisson reconstructs a perfect surface in 19 sec, while Wavelets fails to reconstruct this type of data. ScaleSpace has difficulties finding the correct scale in this dataset, and takes more than 3000 seconds. The result has holes when the triangles are larger than the estimated scale. SuperCocone divides the input points into several subsets using an octree. It sometimes fails to recover triangles that cross octree cells boundaries. Our method misses some triangles resulting into small holes that can be filled in with our post-processing step.

Small details

Figure 12 compares the results on a clean point set having many small bumpy features (the feathers of the eagle), [START_REF]Statue model repository[END_REF]. Poisson faithfully reconstructs the features. Wavelets performs similar to Poisson. ScaleSpace reconstructs a continuous manifold mesh that also follows the bumpy features. SuperCocone and our algorithm give similar results, with more missing triangles than ScaleSpace. However, the missing triangles in our reconstruction are mostly small holes that can be efficiently removed by our post-processing step.

Noise

Real datasets such as raw scans produced as unions of range maps usually exhibit alignment errors. In this data, the points are not located exactly on the underlying surface and form some sort of "solid surface". As it has been stressed previously, our method requires a clean point set since it connects the input points with triangles. The method typically fails when registration errors create two layers of points in overlapping regions of different scans (see Fig. 15). However, managing the over-sampling and the noise of the dataset is an active research subject which is beyond our scope. To accommodate the raw scans in our database we have used 1 or 2 iterations of our simple pre-processing smoothing step which have worked in most cases, as shown by our experiments. Note that the smoothing time is included in the timings given in Figure 8.

In Figure 13, we show the behavior of the algorithms for a noisy point set (original registered dataset for the Stanford bunny [53]). For this data, Poisson and Wavelets both cancel the noise and recover the fine details. ScaleSpace manages to extract a manifold mesh that exactly follows the details of the data, even when it is very bumpy. Our algorithm outputs a result that is very similar to that of SuperCocone. The same dataset with one iteration of pre-processing smoothing is shown in Figure 14. All algorithms had the same smoothed input and gave similar results. Notice the aligned missing triangles in SuperCocone at boundaries between adjacent octree cells.

We stress, however, that our method needs a very moderate smoothing which does not wash out the small details. With the Tanagra model (see Fig. 16), we reconstruct the stamp with the same resolution as Screened Poisson and Scale Space, after 2 iterations of smoothing.

Impact of pre and post-processing steps

Figures 17 and18 show raw LiDAR scans with respectively 14M and 12M input points (from [START_REF]Smart multimedia[END_REF]) and illustrate the impact of the pre-processing smoothing and the postprocessing hole filling steps of our method.

The mesh in Fig. 17, built without pre-processing (in 200s), is quite noisy since it reconstructs the data noise with 123 connected components and 7250 border edges. The mesh generated with 1 iteration of pre-processing (in 123s) is smoother and has only 15 connected components and 1585 border edges. Figures 18 demonstrates the same behavior. Our method manages to reconstruct the raw scans but the resulting meshes have many holes that are left to be filled in by the post-processing step. On this data the smoothing step takes 10s and the post-processing step -30s. The raw candidate triangles generation is done in 30s while the manifold extraction takes between 20s, if the points are smoothed, and 100s otherwise. The total reconstruction times range from 90s to 190s.

Large datasets

For the large datasets, existing Delaunay-based methods (ScaleSpace and SuperCocone) start encountering memory limitations. ScaleSpace crashes on all examples. SuperCocone gives an incorrect result on one of them ('facade-valde-grace') since many triangles are missing and crashes on the others ('chateau-rives' and 'rosetta'). The results are displayed in Figures 19 and Figure 20. A rendering of 'rosetta' is available in the supplemental material.

Results of the 'facade-val-de-grace' are shown in Figure 19. For this dataset, Poisson has problems with some incoherent normals, causing bubbles to appear. Wavelets have similar behavior. It may not be difficult to fix, but the need for having consistent normal orientations is an issue for using approximation based methods. SuperCocone missed many triangles, probably because the implementation starts to encounter memory problems.

We also show, in Figure 20, the 'chateau-rives' example with 33M points. For this dataset, our algorithm outputs a valid result with 65M facets in 361s. Among the other methods, only Wavelets build an approximating coarse mesh with 3M facets in 2050 seconds, with octree depth=10 (See Fig. 20). [START_REF]Large-scale point cloud classification benchmark[END_REF], with respectively 14M and 10M of points. The market place mesh in Fig. 22 has 20M of triangles that have been generated in 200s, with 1 iteration of smoothing (done in 15s). The saint Gallen mesh in Fig. 21 has 28M of triangles generated in 280s, with 1 iteration of smoothing (done in 20s). These meshes are far from being perfect but can be used as initial reconstructions for some downstream geometry processing applications.

Finally, we tested our method with 100 million points (sampled randomly on the Stanford bunny) which does not fit in 16Gb of RAM. Using a 32Gb computer, our method generates 600M candidate triangles in 90s, and then takes 140s to output the raw manifold mesh. With post-processing, the complete algorithm takes 533s and the final mesh has 200M facets. Among the other methods that we tested, only Poisson and Wavelets can also process this dataset. With octree depth=10, Poisson takes 1531s to output a mesh with 3M vertices and 6M facets and Wavelets takes 387 seconds for a mesh with 10M vertices and 5M facets. With octree depth=12, Poisson takes 5809s to output a mesh with 48M vertices and 96M facets and Wavelets takes 762 seconds to build a mesh with 162M Figure 17: A raw LiDAR scan (14M points) of a sandcastle (note the shovel). Top, our mesh reconstructed with no pre-processing (built in 200s) Bottom, our result reconstructed after 1 iteration of smoothing (built in 123s).

vertices and 81M facets. With depth=13, both do not fit in 32 Gb.

Conclusion and Discussion

In this paper, we have presented a surface reconstruction algorithm from a set of points in 3D based on a simple observation about the restricted Voronoi cells. Our algorithm has a trivial implementation, is embarrassingly parallel and economic in terms of memory resources. Therefore it scales up very well to datasets with tens of millions vertices while keeping computation times smaller than 5 minutes. Its simplicity makes it possible to embed it into devices such as smart-phones and cameras with a CPU (an Android ARMV7 implementation is provided in the supplemental material). Compared with the state-of-the-art methods, it is not as good as Poisson and Wavelets for filling the holes, as expected since these two methods build smooth approximation meshes. Concerning the interpolation methods, our method has a behavior similar to the Cocone and ScaleSpace even if the latter may have a better set of candidate triangles, when it can estimate the correct scale (see e.g. the eagle dataset).

The algorithm outputs exactly the Delaunay triangulation of the input points restricted to the union of the D r i disks (Section 3.4). In the experiments, we observed that the behavior of our algorithm is very similar to the Cocone algorithm (Table 1). because the disks that we use do not deviate much from the co-cones. From a theoretical point of view, it would be interesting to characterize the configurations where the outputs of both algorithms match exactly. This probably involves some variants of -sampling conditions.

Regarding timings, our algorithm is the fastest, by two order of magnitudes, in some cases. However, one needs to take into account that the SuperCocone's implementation that we obtained from its authors is sequential, and could probably be easily parallelized. We think that our algorithm will still be faster since it does not need to compute any 3D Delaunay triangulation (it directly computes the Restricted 2D Delaunay triangulation).

Unlike Poisson reconstruction and its variants, our algorithm does not need coherently oriented normals, it solely uses normal directions. This can be interesting when such information is not available or when the input point set has a shape that is too complicated for greedy orientation algorithms, such as in Fig. 19. Our algorithm requires a moderate smoothing when the data presents alignment errors, but it does not require a per-sample resolution es-Figure 19: A raw scanned point set with 30M vertices. On this dataset, some normals are not correct causing some bubbles in the Poisson reconstruction (built in 856s). SuperCocone gives a result in 5011 seconds, but has many missing triangles (maybe it starts reaching memory limits). Our algorithm (with 2 iterations of smoothing) takes 340s timate as in [START_REF] Fuhrmann | Floating scale surface reconstruction[END_REF].

Observing the relative timings of the different steps of our algorithm, manifold extraction and post-processing clearly dominate the core of the algorithm. However, they are needed by most geometry processing applications that rely on a clean manifold mesh. For future work, there is room for improvement in these steps. Our implementation of hole filling is not well optimized, not parallel, and we think it could probably be made two or three times faster.

The complete implementation in C++ is available in our open source GEOGRAM library [START_REF] Geogram | [END_REF]. A raw point set with 33M points. On the top, the Wavelet algorithm builds, with depth=10, a coarse approximation mesh with 6M vertices. On the bottom, our algorithm builds a detailed mesh with 65M vertices. The total computation takes 360 seconds -49s for 2 iterations of smoothing and normal estimation, 78s reconstruction, 61s surface extraction, 114s post-processing and 58s loading and saving.

Figure 1 :

 1 Figure 1: Overview of the algorithm.

Figure 2 :

 2 Figure 2: We compute the Voronoi diagram of the input point set restricted to a union of disks centered on the points and orthogonal to the estimated normals. From A to C: the radius of the disks is increased from 0.2% to 4% of the bounding box diagonal. All the restricted Voronoi cells (colored polygons) are computed in parallel.

Figure 3 :

 3 Figure 3: The Voronoi cell V (p i ) is the intersection of the halfspaces defined by the bisectors of the segments [p i , p j ]. Once a point (here p j 8 ) is further away twice than the current Voronoi cell at p i , it no longer contributes to the cell.

Figure 4 :

 4 Figure 4: The Delaunay triangles (p i , p j 1 , p j 2 ) are deduced from the indices of the bisectors that generate the Voronoi vertices q.

Figure 5 :

 5 Figure 5: A: the set of triangles seen from three restricted Voronoi cells (T 3 ), with some holes (yellow). B: the set of triangles seen from one or two Voronoi cells (T 1,2 , in red) fills some holes, but have non-manifold configurations (four red triangles in the right). C: result after carefully inserting the T 1,2 triangles one by one. D: suppressing the remaining holes.

Figure 6 :

 6 Figure 6: Configurations of 4 co-circular sample points (black dots) with their partial restricted Voronoi facets (sketched with black segments). The dual T 12 triangles are shown with different colors corresponding to their cardianlities.

Figure 7 :

 7 Figure 7: (a) Simple bridge. (b) Double bridge.

Figure 8 :

 8 Figure 8: Comparison of processing times for the methods. Timings are in seconds and model sizes in million vertices. We used implementations written by the authors of the cited methods, and the same computer for all the tests. For Poisson and Wavelets we used an octree with depth 10. ScaleSpace and SuperCocone adapt the depth of their octrees (from 6 to 10). KO means that the result has too many holes or is mostly unrelated with the data or the computation time explodes as compared to data of similar size.

Figure 9 :

 9 Figure 9: The used sampling patterns. A: homogeneous; B: homogeneous regularized; C: curvature adaptive

Figure 10 :

 10 Figure 10: Performance of the tested algorithms. Up, timings obtained on the Small data. Bottom, timings obtained on the Large data.

Figure 11 :

 11 Figure 11: Comparison of the algorithms with a small dataset (100K vertices) with high variations of sampling density.

Figure 12 :

 12 Figure 12: Comparison of the algorithms with a clean dataset with small-scale bumpy features.

Figure 13 :

 13 Figure 13: A noisy dataset (original scanner data of the Stanford Bunny). Poisson both filters out the noise and recovers the fine features. ScaleSpace successfully extracts a closed manifold mesh. SuperCocone and our algorithm recover a smaller number of triangles.

Figure 14 :

 14 Figure 14: The same dataset as in Figure 13 with one iteration of smoothing applied. The methods give similar results.

Figure 15 :

 15 Figure 15: A: overlapping scans; B: our reconstruction without smoothing; C: our reconstruction with one iteration of smoothing (data from [4])

Figure 16 :

 16 Figure 16: Reconstruction of the Tanagra sculpture (detail). Left: Screened Poisson reconstruction (octree depth=12). Right: our result (with 2 iterations of smoothing) where the fine-scale details are preserved.

Figure 18 :

 18 Figure 18: A raw LiDAR scan (12M points) from [49]. Ours results in function of pre-and post-processing steps.

Figure 20 :

 20 Figure 20:

Acknowlagement. B. Levy's work was partially funded by INRIA Exploragram grant and ANR MAGA and ROOT grants.