Emmanuel Agullo
email: emmanuel.agullo@inria.fr

Cédric Augonnet
email: cedric.augonnet@inria.fr

Jack Dongarra

Mathieu Faverge
email: faverge@eecs.utk.edu

Julien Langou
email: julien.langou@ucdenver.edu

Hatem Ltaief
email: ltaief@eecs.utk.edu

Stanimire Tomov
email: tomov@eecs.utk.edu

Dongarra

LU Factorization for Accelerator-based Systems

Keywords: High Performance Computing, Dense Linear Algebra, LU Factorization, Hybrid Architecture, Multiple GPU Accelerators, Multicore, Tile Algorithm, Numerical Accuracy

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

I. INTRODUCTION

The LU factorization (or decomposition) of a matrix A consists of writing that matrix as a matrix product A = LU where L is lower triangular and U is upper triangular. It is a central kernel in linear algebra because it is commonly used in many important operations such as solving a nonsymmetric linear system, inverting a matrix, computing a determinant or an approximation of a condition number. In most cases, the decomposition is the computationally dominant step. For instance, solving a linear system (Ax = b) of n equations can be performed as an LU factorization (A = LU ), requiring θ(n 3 ) operations, and two triangular substitutions (Ly = b and U x = y), requiring θ(n 2 ) operations. The performance of the decomposition is thus critical and needs to be optimized especially for the emergent platforms. This is why a large effort was recently dedicated to accelerate it on multicore platforms [START_REF] Buttari | A class of parallel tiled linear algebra algorithms for multicore architectures[END_REF], [START_REF] Chan | Supermatrix Out-of-Order Scheduling of Matrix Operations for SMP and Multi-Core Architectures[END_REF]. In this paper, we propose to accelerate the LU factorization on a multicore node enhanced with multiple GPU accelerators. We follow a methodology previously employed in the context of the Cholesky factorization [START_REF] Agullo | GPU Computing Gems[END_REF] and QR factorization [START_REF] Agullo | QR Factorization on a Multicore Node Enhanced with Multiple GPU Accelerators[END_REF] that we apply to the tile LU decomposition algorithm [START_REF] Buttari | A class of parallel tiled linear algebra algorithms for multicore architectures[END_REF]. We bring four contributions. First, we present the design of new CPU/GPU hybrid kernels for performing an LU factorization on a GPU associated to a CPU core. Second, we study the impact on performance of the looking variants as well as the storage layout (row or column major) in presence of pivoting. Third, we adapt and tune our kernels for the recent NVIDIA Tesla S1070 and S2050 (Fermi) Computing Processors and we present experimental results on two machines composed of such multiple GPUs. And last but not least, we present the impact of the number of tiles on the numerical accuracy.

The reminder of the paper is organized as follows. In Section II, we present the LU factorization, some related work and our experimental environment. The LU algorithm we rely on, so-called tile LU factorization, splits the computational routines into tasks of fine granularity that can be run concurrently on different GPUs. We describe the individual tasks, the kernels, and their optimized hybrid CPU/GPU implementation in Section III. In sections IV, V and VI, we show how we schedule those operations through two different runtime environments and discuss the subsequent performance. We briefly discuss the numerical accuracy of tile LU factorization in Section VII to motivate the use of an efficient double precision factorization. We conclude and present our future work in Section VIII.

II. BACKGROUND

A. Standard LU Algorithm

LAPACK and SCALAPACK are the current de facto standard libraries for performing advanced dense linear algebra operations. The LU factorization implemented in those packages is designed as a high-level algorithm relying on basic block routines from the Basic Linear Algebra Subprograms (BLAS) and the Basic Linear Algebra Communication Subprograms (BLACS) for LAPACK and SCALAPACK, respectively. The matrix is conceptually split in blocks of columns, so called panels. At each column of the panel being factorized, a partial pivoting scheme is used. It consists of finding the element of highest magnitude (socalled pivot) within the column, swapping the corresponding rows and then dividing the other column elements by that value. This process is repeated for each column within the panel. The factorization then proceeds by an update procedure of the trailing submatrix at each step of the factorization until all the subsequent panels are processed. This panel-update sequence, called blocking, enables the use of the level-3 BLAS matrix multiplication (gemm) leading to a better data reuse, which is critical on cache-based architectures. However, a synchronization point is required in-between, similar to a fork-join paradigm. Furthermore, the parallelism only occurs within each phase (panel or update) and is expressed at the level of BLAS and BLACS, which ultimately leads to a limited performance.

B. Tile LU Algorithm

To alleviate this bottleneck, so-called tile algorithms permit to break the panel factorization and the update of the trailing submatrix into tasks of fine granularity. The update phase can now be initiated while the corresponding panel is still being factorized. Initially proposed for updating factorizations [START_REF] Golub | Matrix Computations[END_REF], the idea of tile algorithms consists of splitting the panels in square submatrices (so-called tiles). Using updating techniques to obtain tasks of fine granularity was first exploited in the context of out-of-memory (often called out-of-core) factorizations [START_REF] Yip | Fortran subroutines for out-of-core solutions of large complex linear systems[END_REF]. Two implementations of the tile LU factorization for multicore architectures were more recently proposed [START_REF] Buttari | A class of parallel tiled linear algebra algorithms for multicore architectures[END_REF], [START_REF] Chan | Supermatrix Out-of-Order Scheduling of Matrix Operations for SMP and Multi-Core Architectures[END_REF]. Tile LU algorithm annihilates matrix elements by tiles instead of rectangular panels as in LAPACK. In this paper, we use the high-level algorithm presented in [START_REF] Buttari | A class of parallel tiled linear algebra algorithms for multicore architectures[END_REF], which relies on four kernels and whose first panel-update sequence is unrolled in Figure 1.

C. Related work

Various projects provide GPU kernels for linear algebra. The CULA library implements many BLAS and LAPACK kernels [START_REF] Humphrey | CULA: hybrid GPU accelerated linear algebra routines[END_REF], but those are limited to a single GPU and problems that fit into the memory. Kurzak et al. accelerated the QR decomposition of the PLASMA library with a single GPU by designing non standard BLAS kernels. The PLA-PACK [START_REF] Fogué | Retargeting plapack to clusters with hardware accelerators flame working note #42[END_REF] and the DPLASMA [START_REF] Bosilca | A unified HPC environment for hybrid manycore/GPU distributed systems[END_REF] linear algebra libraries both target accelerator-based clusters. Both of them rely on their own runtime systems that provide intra-node dynamic scheduling capabilities.

Different runtime systems were designed to support accelerator-based platforms. The StarSs [START_REF] Ayguadé | An Extension of the StarSs Programming Model for Platforms with Multiple GPUs[END_REF] language is an annotation-based language executed either on CPUs, GPUs or Cell processors, respectively, using SMPSs, GPUSs or CellSs. Diamos and Yalamanchili also propose to use features that are similar to those available in StarPU in the Harmony runtime system [START_REF] Diamos | Harmony: an execution model and runtime for heterogeneous many core systems[END_REF]. Sequoia [START_REF] Fatahalian | Sequoia: Programming the memory hierarchy[END_REF] statically maps hierarchical applications on top of clusters of hybrid machines, and the Charm++ runtime system was extended to support GPUs [START_REF] Jetley | Scaling Hierarchical N -body Simulations on GPU Clusters[END_REF].

D. Experimental environment

The experiments have been conducted on two different hybrid platforms in single and double precision arithmetics. The first machine is a quad-socket quad-core host machine based on an AMD Opteron(tm) Processor operating at 2.4 GHz. The cache size per core is 512 KB and the size of the main memory is 32 GB. A NVIDIA Tesla S1070 graphical card is connected to the host via PCI Express 16x adapter cards. It is composed of four GPUs C1060 with two PCI Express connectors driving two GPUs each. Each GPU has 4 GB GDDR-3 of memory and 30 processors (240 cores), operating at 1.44 GHz. The computational time ratio between single and double precision arithmetics on this Tesla S1070 is 6. We will refer to this machine as Tesla-C1060. The second machine is a dual-socket hexa-core host machine based on INTEL Nehalem Processors operating at 2.67 GHz. Each core has 256 KB of L2 cache and each socket has 12 MB of L3 cache. The size of the main memory is 48 GB. Three NVIDIA Fermi C2050 cars are connected to the host with a 16x PCI bus. Each GPU has 3 GB of GDDR-5 of memory and 14 processors (448 cores), operating at 1.15 GHz. ECC was disabled on these cards to have as much memory as possible. The computation in double precision arithmetic has been fixed in this last GPU generation and the ratio compared to single precision arithmetic is now two. We will refer to this machine as Fermi-C2050.

III. KERNEL IMPLEMENTATION

A. Design

The tile LU algorithm of interest is based on four computational kernels, described as follows:

GETRF This kernel performs the standard LU factorization of a tile as in LAPACK and generates a permutation vector according to the pivot selection;

GESSM This kernel applies to a tile from the left the permutations followed by a triangular solve using the permutation vector and the lower triangular tile computed by the GETRF kernel;

TSTRF This kernel performs the block LU factorization on a pair of tiles, one on top of the other, where the upper tile (also the diagonal tile) is already upper triangular. A permutation vector is also produced as well as a unit, lower triangular tile as a result of swapping rows in order to prevent overwriting of existing data;

SSSSM This kernel applies the permutations and the transformations computed from the TSTRF kernel to two stacked tiles from the left. A triangular solve is applied to update the top tile using the unit, lower triangular tile from TSTRF. A matrix-matrix multiplication is then applied to update the bottom tile using the full tile computed by TSTRF.

Critical for the performance of the overall algorithm, SSSSM is the most compute intensive kernel and asymptotically dominates all three other kernels in terms of flops. We integrated MAGMA BLAS matrix-matrix multiplication (gemm) and triangular matrix solver (trsm), specifically tuned for the tile LU algorithm. Moreover, we developed a hybrid TSTRF kernel which takes advantage of look-ahead technics, similar to the one-sided hybrid LU factorizations in MAGMA for single GPU. The other two kernels, GETRF and GESSM, were adapted and derived from the MAGMA library. Those kernels have been implemented in single and double precision arithmetics.

B. Tuning

Two parameters must be selected to tune the tile-LU algorithm: the size of the tiles, and the internal blocking used by the different kernels. Empirical tuning was performed on the most compute intensive kernel (SSSSM). While it is crucial to optimize the most time consuming kernel, large internal blocking results in a fast SSSSM kernel but slows down the overall performance due to the overhead brought by the extra flops [START_REF] Agullo | Comparative study of one-sided factorizations with multiple software packages on multi-core hardware[END_REF]. Several sets of parameters were thus tested to find the combination that not only maximizes the performance of the kernel, but also that of the overall algorithm.

Once tuned for a specific hardware, an efficient runtime environment to ensure a proper scheduling across the host (multicore system) and the device (GPU accelerators) becomes paramount.

IV. STATIC SCHEDULING

A. Design

Originally implemented to schedule the Cholesky and QR factorizations on the Cell processor [START_REF] Kurzak | QR factorization for the CELL processor[END_REF], the hand-coded static scheduler has been extended to handle the dispatch of the hybrid kernels (from Section III) across all CPU-GPU pairs available on the system in one dimensional cyclic fashion. Figure 2 shows how the column-wise partitioning is achieved. This static runtime imposes a linear scheduling order on all the kernels during the factorization. This order enforces the execution of a predetermined subset of kernels on a particular CPU-GPU pair.

Two distinct global progress table are required to first ensure numerical correctness and second, to get high performance. The first progress table keeps track of dependencies among the different tasks at different steps of the factorization. A dependency check is performed before executing each kernel by examining the local copy of the progress table. The hosting CPUs stall with busy waiting on volatile variables until the corresponding dependencies are satisfied, which simultaneously triggers the release of the kernel to the designated CPU-GPU pair.

The second progress table is required to avoid unnecessary copies from/to CPU to/from device GPU and therefore, optimize the performance of the overall algorithm. In fact, at any specific time of the computation, a same data tile can reside both on the host's and/or device's main memory (MM). The issue is to determine the consistency of the data. This problem is well-know in cache coherency protocols. The second progress table helps to keep track of:

• which data tile is considered as final output to ensure data consistency with the other processing units (PUs).

If the tile is located on the device's MM, it is copied back to the host's MM for a true shared access with the other PUs available on the system.

• which data tile is still transient and has not yet reached its final state. This tile cannot be shared among other PUs (false-sharing). This progress table is incremented at each step of the factorization. Each CPU thread hosting a GPU probes the progress table about the state of the data tile it is about to access. If the data tile state is final, no transfer is necessary and the computation proceeds. If its state is transient, the CPU thread will initiate a transfer to its GPU in order to acquire the up-to-date data tile prior to execute its scheduled task.

This turns out to considerably decrease the number of communications involved between a CPU-GPU pair, which is critical given that the PCI bus is two orders of magnitude less efficient than the computational power of the accelerators.

V. ALGORITHMIC OPTIMIZATIONS

We introduce four incremental optimizations integrated into the application to improve the overall performance. First, we manually adjust the order of the loops, and then we optimize the compute kernels so that they fit to the specificities of GPUs.

A. Looking Variants

We study in this section the effect of three looking variants, which have an influence on the overall scheduling.

The right-looking variant (RL) consists in updating first the trailing submatrix, located on the right side of the current panel, before the execution of the next panel starts. This looking variant generates many concurrent tasks which can potentially run in parallel.

The left-looking variant (LL), also called the lazy variant, applies all subsequent updates generated from the left side to the current panel before proceeding with the next panel. Since the updates are restricted to the panel only, data reuse is maximized while at the same time parallelism can get limited.

The methodology of this static scheduling allows for pipelined execution of factorizations steps, which usually provides similar benefits to dynamic scheduling to some extent. For example, the execution of the inefficient Level 2 BLAS operations can be overlapped by the efficient Level 3 BLAS operations. This phenomenon has been successfully demonstrated for the scheduling of one-sided factorizations in the context of homogeneous multicore [START_REF] Kurzak | Scheduling dense linear algebra operations on multicore processors[END_REF]. However, when tackling hybrid system components, i.e., multicore associated with GPU accelerators, a tremendous gap in terms of sequential performance may exist between the hybrid CPU-GPU kernels and the GPU kernels as seen in Section III-B. Therefore, one of the main disadvantages of the static scheduling is a potential suboptimal scheduling, i.e., stalling in situations where work is available.

Figures 3(a) and 3(b) clearly describe this drawback of the static scheduling on four CPU-GPU pairs. The dark purple colors represent the panel tasks (GETRF and TSTRF) and the light green colors are the update kernels (GESSM and SSSSM). The panel tasks are hybrid and the GPU needs the CPU to perform the Level 2 BLAS operations while the update kernels are highly efficient Level 3 BLAS operations performed on the GPU only. Figure 3(a) shows the RL variant with lots of stalls represented by white empty spaces. The panel tasks indeed become a bottleneck and the updates tasks cannot proceed until the completeness of the panel tasks. Figure 3(b) presents the LL variant. The scheduling contains less gaps but still suffers from the lack of parallelism, especially in the beginning. And this inefficiency is even more exacerbated by the slow panel hybrid kernels.

A new looking variant has then been implemented to alleviate this bottleneck combining the previous LL version with a breadth-first search task execution (BF-LL). Each CPU-GPU pair applies all subsequent transformations, once for all, (update and panel tasks) on a particular tile on the current panel before proceeding with the tile below it. The obtained trace as seen in Figure 3(c) is very compact and dense as shown in the bottom trace. On Figure 4, the BF-LL variant therefore clearly outperforms the other variants for all problem sizes. Noteworthy to mention that some natural load imbalance starts to appear toward the end of the traces. This explains the irregularities observed on the different curves of Figure 4. tially i.e., they swap rows one after another. This operation can be extremely slow on GPUs. For that reason, we have implemented a single GPU kernel that applies all permutations within a block at once. Applying this technique to the BF-LL variant improves the performance by a few percents on Figure 5.

C. Triangular Solve

The second optimization consists in splitting the call to triangular solve (trsm) on the GPU by a triangular matrix inversion (trtri) on the CPU followed by (1) a triangular matrix-matrix multiplication (trmm) and (2) a matrix-matrix multiplication (gemm) both on the GPU to further improve the performance of the most compute intensive GPU kernel. The same modification has been proposed by Volkov and Demmel [17]. This affects the overall stability analysis of the algorithm but, as long as a reasonable pivoting strategy is used (so that the norm of L is not too large), we can guarantee that the scheme stays normwise stable. The benefit in terms of performance is however very significant: on Figure 5, the asymptotic speed is now 821 Gflop/s compared to 638 Gflop/s before.

D. Storage Layout

One critical part of the standard LU factorization is the pivoting scheme. The LU algorithm as implemented in LAPACK uses partial pivoting which generates swapping between rows. The tile LU factorization uses a pairwise pivoting scheme which implies pivoting between a pair of tiles. This generates more pivoting than the standard LU factorization and also presents some stability issues which are explained in Section VII. The storage layout optimization enables to improve the memory accesses on the GPU when pivoting by transposing the matrix accordingly. The swapping of rows, which are now contiguous in memory, can use coalescent memory accesses to substantially improve the overall performance of the application. Implementation and experimentation were therefore carried out with tiles in two storage layouts -column and row major correspondingly. On Figure 5, we obtain about 882 Gflop/s by converting the matrix into a Row-Major layout. Previous work for example, has shown that the computational cost for pivoting can be prohibitively high for performance. Figure 6 shows the strong scalability of the SGETRF algorithm with respect to the number of GPU-CPU pairs. We obtain a 1.9 speedup on 2 GPUs and a 3.6 speedup on 4 GPUs. While the successive optimizations ensure that enough parallelism is created, the scalability is limited by the contention on the bus. As shown on the trace on Figure 3(c), this statically scheduled algorithm also suffers from some load imbalance which affects its scalability.

E. Scalability

VI. DYNAMIC SCHEDULING

The dynamic scheduler described in this section proposes a more productive and portable approach to schedule the tile LU factorization across a hybrid system.

A. Design

The application to be scheduled can be split into multiple tasks. We use the StarPU [START_REF] Augonnet | StarPU: a Runtime System for Scheduling Tasks over Accelerator-Based Multicore Machines[END_REF] runtime system to schedule these tasks dynamically simultaneously on CPU cores and on CUDA devices. StarPU is a library that schedules tasks on accelerator-based platforms. It automates data management and provides a flexible framework to design scheduling policies. The two main steps required to port an application on StarPU are first to register all data to StarPU, and then to split the computation into tasks that access registered data. StarPU provides performance portability by making it possible to concentrate on purely algorithmic concerns and on kernel optimization rather than dealing with low-level and non portable problems such as data management.

1) Data management: Once a piece of data has been registered, StarPU ensures that it is available for the different tasks that need to access it, either from a CPU or from a GPU. StarPU keeps track of the location of the different data replicates and make sure that they remain coherent: when a task modifies a piece of data, all other replicates are automatically invalidated.

2) Task scheduling: The task structure proposed by StarPU contains a multi-versioned kernel, named codelet, and the list of data which should be accessed by the task, as well as the access mode (e.g., read or read/write). The codelet structure encapsulates the different implementations of the task on the different types of processing units: when StarPU assigns the task to a certain type of processing unit (e.g., CUDA or OpenCL), the appropriate implementation is called. StarPU will not assign a task to a processing unit for which the codelet is not implemented. It is therefore not required to give an implementation of the codelet for each type of processing unit.

Task dependencies are either explicit or implicitly derived from data dependencies. In the tile-LU algorithm, StarPU automatically detects dependencies so that we simply unfold the DAG describing the tile-LU algorithm in a sequential fashion.

StarPU provides various scheduling strategies to address a wide scope of applications. For the tile-LU algorithm, we used the HEFT-TMDP-PR policy. This strategy relies on performance predictions based on previous executions, which make it possible to accurately distribute the load between the heterogeneous processing units. This strategy also minimizes data transfer overhead by initiating asynchronous data transfers as soon as a task is attributed to a processing unit prior to its actual execution. The total amount of data transfers is also reduced by predicting data transfer time, and therefore avoiding to schedule tasks on processing units that do not have a local copy of the input data. More details on the HEFT-TMDP-PR policy are available in a previous study [START_REF] Augonnet | StarPU: a Runtime System for Scheduling Tasks over Accelerator-Based Multicore Machines[END_REF].

B. Performance results

Figure 7 present performance results obtained with either the static or the dynamic schedulers. The three optimizations of sections V-D, V-B and V-C are applied again. Dynamically scheduled tasks are submitted following the Right-Looking ordering. We therefore show the speed of both the Right-Looking (RL) and the Breadth-First Looking (BF-LL) statically scheduled algorithms. It is worth noting that data transfer time is not accounted in the static algorithms. Contrary to the static algorithm that was manually tuned to minimize the amount of data transfers and to balance the load, the dynamic scheduler has to take decisions even before all tasks have been submitted. For small problems, static algorithms thus perform better than the dynamic one. For larger problems, the dynamic scheduler is able to distribute the load better than the RL algorithm even though tasks are submitted in the same order in both cases.

The static algorithms are unable to handle problems larger than the size of the GPUs, while StarPU transparently handles larger problems. In the future, we plan to unfold the dynamically scheduled DAG following the BF-LL order in order to guide the dynamic scheduler. This would make it possible to benefits from the expertise introduced within the static algorithm while preserving code maintainability.

Performance portability is also illustrated on the curve obtained when using all CPUs in conjunction with the GPUs on Figure 7.

VII. NUMERICAL ACCURACY

We have mentioned above that the use of double precision may be critical when relying on the Tile LU factorization. We now illustrate this claim. The purpose of this section is not to accomplish a full study of the numerical stability and accuracy of the Tile LU algorithm; this would deserve an entire manuscript. Our goal is to justify that Tile LU algorithm can be used in practice on the matrices we usually consider and to highlight its limits.

Pivoting is necessary during the LU factorization process in order to ensure the stability of the algorithm. The standard way to perform pivoting is partial pivoting: for each column, one pivot is used to eliminate all the others elements of this column. The pivot is chosen to be an element in the column with the largest magnitude. This scheme is, for example, implemented in LAPACK and ScaLAPACK and it is practically stable; in the sense that pathological matrices exist for which it is not stable; but, these matrices have not been observed in practical applications. During the Tile LU algorithms, the partial pivoting scheme is altered in order to enable a pipeline in the update of the columns. Indeed, for a single column, several pivots can be used during the elimination process. The pivoting scheme for Tile LU is a block variant of pairwise pivoting and we call it block pairwise pivoting. Stability of pairwise pivoting has been studied in [START_REF] Trefethen | Average-case stability of Gaussian elimination[END_REF]. The authors experimentally showed that pairwise pivoting is less stable than partial pivoting. This motivates this experimental section on the stability of block pairwise pivoting. A preliminary study has been proposed in [START_REF] Quintana-Ortí | Updating an LU factorization with pivoting[END_REF] but this only addresses the case of two tiles. In practical applications, our algorithms will have tens to hundreds of tiles.

Our main result can be stated as follows. We experimentally observe that, when compared to partial pivoting, the stability of block pairwise pivoting is independent of the size of the tiles and solely depends on the number of tiles in a column. As the number of tiles increases, the accuracy of the factorization deteriorates. We observe that, on random matrices tile LU performed on 1, 000-by-1, 000 tile matrix loses at most 2-digit of accuracy with respect to partial pivoting.

Here, we consider ten matrices of order ranging from 1024 to 10240 that were randomly generated with a uniform distribution. For each matrix, Figure 8 shows the infinity norm of the residual ( Axb ∞) obtained for Tile LU relatively to the residual obtained with LAPACK. LAPACK LU algorithm implements a standard partial pivoting scheme consisting of choosing the element of highest magnitude as pivot for the current column. When using one single tile, Tile LU exactly performs a partial pivoting scheme. Therefore, Tile LU and LAPACK lead to the exact same residual in that case (t = 1, leftmost point of Figure 8). When the number of tiles per row (the matrix is composed of t×t tiles if there are t tiles per row) increases, Figure 8 shows that the numerical accuracy decreases. One decimal digit of accuracy is lost (ratio equal to 10) when 10 to 30 tiles are used, depending on the matrices considered. On the right-most part of the graph (the scale is logarithmic), we study the loss of accuracy when considering up to an improbable number of 1, 000 tiles that would not fit in the memory of our target machines if we were considering a tile size optimized for performance. Only two matrices out of ten loose a second decimal digit of accuracy (ratio equal to 100). One of our main observation is that the stability is influenced by the number of tiles used in the algorithm, and is not influenced by the size of the tiles. If this cannot be considered as a statistical result, it still shows that the degradation of Tile LU on numerical accuracy is not necessarily dramatic. Indeed, an iterative refinement procedure often recovers several decimal digit per iteration. Therefore, a Tile LU factorization in double precision is likely to lead to a very acceptable numerical accuracy. On the contrary, in single precision, the numerical accuracy of the standard LU factorization with partial factorization is often barely high enough to yield to the solution. Therefore, the computation of the solution of a linear system may fail because of one or two lost decimal digits. As a partial conclusion, for the matrix sizes and tile sizes we consider, we have shown that the obtained accuracy is critical in single precision but is likely to be acceptable in double precision, emphasizing the importance of efficient numerical solvers in double precision. On Figure 9, the dynamically scheduled DGETRF kernel reaches 500 Gflop/s on our Fermi-based system, which is half of the speed obtained by the single precision algorithm on the same platform. The different results that we have presented in single precision are also directly applicable to the double precision implementation.

VIII. CONCLUSION

We have presented the design and implementation of a new hybrid algorithm for performing the tile LU factorization on a multicore node enhanced with multiple GPUs. This work completes our study of one-sided factorizations on such platforms. We have shown that we could reach very high performance thanks to new CPU/GPU hybrid kernels. These kernels have been specifically tuned to fully exploit the potential of the NVIDIA Tesla S1070 and the new generation Fermi-based S2050 GPU accelerators. Tile LU achieves 1 Tflop/s in single precision arithmetic and 500 Gflop/s in double precision on the Fermi-based system. The numerical accuracy of the LU factorization is very sensitive to the underlying algorithm. From a high-level point view, we rely on the Tile LU factorization algorithm. We have exhibited that the numerical accuracy decreases when the number of tiles increases. For the matrix sizes and tile sizes we consider, we have shown that this degradation is critical in single precision but acceptable in double precision, highlighting the importance of using efficient double precision hardware and software.

This work aims at unifying PLASMA (initially designed for homogeneous multicore architectures) and MAGMA (initially designed for single-GPU architectures) and the resulting software will be incorporated in those libraries. The purpose is to constitute a library equivalent in functionality to LAPACK. We plan to pursue this work with two-sided factorizations (Hessenberg reduction, tridiagonalization and bidiagonalization). We are also investigating memory bound algorithms such as the solution step that follows a onesided factorization. For those algorithms, minimizing data movement is critical. This is still a complex and open problem on heterogeneous machines.

Figure 1 .

 1 Figure 1. Tile LU Factorization: first panel factorization and corresponding updates.

Figure 2 .

 2 Figure 2. Task partitioning across six processing units.

  (a) Right Looking (RL). (b) Left Looking (LL). (c) Breadth-First (BF-LL).

Figure 3 .

 3 Figure 3. Trace of the different looking variants.

Figure 4 .Figure 5 .

 45 Figure 4. Looking variants for static SGETRF on Tesla-C1060

Figure 6 .

 6 Figure 6. Scalability of static SGETRF on Tesla-C1060

Figure 7 .

 7 Figure 7. Static vs. Dynamic scheduling on the 3 GPUs of Fermi-C2050

Figure 8 .

 8 Figure 8. Numerical accuracy of the residual ( Axb ∞, logarithmic y-axis) when the number of tiles per row (t, logarithmic x-axis) varies, for random matrices A of various orders N in double precision. The residual is normalized to the residual obtained with LAPACK.

Figure 9 .

 9 Figure 9. Scalability of dynamic DGETRF on Fermi-C2050


Pierre-Yves Laont

Adrien Bousseau

George Drettakis

Rich Intrinsic Image Separation for Multi-View Outdoor Scenes

Keywords: Intrinsic images, image-guided propagation, multi-view stereo, mean-shift algorithm

Intrinsic images aim at separating an image into its reectance and illumination components to facilitate further analysis or manipulation. This separation is severely ill-posed and the most successful methods rely on user indications or precise geometry to resolve the ambiguities inherent to this problem. In this paper we propose a method to estimate intrinsic images from multiple views of an outdoor scene without the need for precise geometry or involved user intervention. We use multiview stereo to automatically reconstruct a 3D point cloud of the scene.

Although this point cloud is sparse and incomplete, we show that it provides the necessary information to compute plausible sky and indirect illumination at each 3D point. We then introduce an optimization method to estimate sun visibility over the point cloud. This algorithm compensates for the lack of accurate geometry and allows the extraction of precise shadows in the nal image.

We nally propagate the information computed over the sparse point cloud to every pixel in the photograph using image-guided propagation. Our propagation not only separates reectance from illumination, but also decomposes the illumination into a sun, sky and indirect layer. This rich decomposition allows novel image manipulations as demonstrated by our results.

1 Introduction

Editing material properties and lighting conditions is a common image manipulation task that requires signicant expertise to achieve consistent results.

The main diculty in such manipulations resides in the fact that a pixel color aggregates the eect of both material and lighting, so that standard color manipulations are likely to aect both components.

A typical usage scenario would be to take an interesting photograph, and manipulate its content after the capture, e.g., to change the color of an object, modify the lighting of the scene, or add in virtual objects. In this paper we present a method to achieve this just by taking a few extra photographs, which represents a minimal additional cost. We focus on outdoor scenes and separate photographs into a material layer (also called reectance) and several illumination layers that describe the contributions of sun, sky and indirect lighting (Fig. 1e-h). This intrinsic images decomposition [START_REF] Barrow | Recovering intrinsic scene characteristics from images[END_REF] allows easy editing of each component separately, and subsequent compositing with consistent lighting. We illustrate applications of our decomposition in Fig. 1. We rst alter the oor material with a grati while maintaining consistent shadows (b); we then add a virtual object in the scene with consistent lighting (c); and we nally change the lighting color and blur the shadows to simulate sunset (d). Please also see the video which illustrates the image manipulation process.

Estimating intrinsic images from photographs is an ill-posed problem that has been traditionally addressed with restrictive assumptions on reectance and illumination variations [START_REF] Land | Lightness and retinex theory[END_REF], signicant user intervention [START_REF] Bousseau | User-assisted intrinsic images[END_REF], or multiple images under varying lighting [START_REF] Weiss | Deriving intrinsic images from image sequences[END_REF] captured over several hours for outdoor scenes. In addition, most existing methods only generate a single illumination layer, grouping the eects of sun, sky and indirect lighting. Inverse global illumination methods [START_REF] Yu | Inverse global illumination: recovering reectance models of real scenes from photographs[END_REF][START_REF] Debevec | Estimating surface reectance properties of a complex scene under captured natural illumination[END_REF] also estimate the material properties of a scene but require very accurate geometric models, either built by hand or acquired with complex laser scanners.

In our approach we address the above shortcomings by providing an illu-

RR n°7851

mination decomposition with minimal user intervention, resulting in intrinsic image layers for sun, sky and indirect illumination. We achieve this by combining sparse geometric reconstruction [START_REF] Snavely | Photo tourism: exploring photo collections in 3D[END_REF][START_REF] Furukawa | Accurate, dense, and robust multiview stereopsis[END_REF] with image-guided propagation [START_REF] Levin | Colorization using optimization[END_REF],

thus leveraging their respective strengths. We exploit the automatically reconstructed 3D information to compute lighting information for a subset of pixels, and use image-guided propagation to decompose the photographs into intrinsic images.

Our algorithm takes as input a small number of photographs of the scene captured at a single time of the day, along with an environment map. From this lightweight capture we use recent computer vision algorithms to reconstruct a sparse 3D point cloud of the scene. The point cloud only provides an imprecise and incomplete representation of the scene. However, we show that this is sufcient to compute plausible sky and indirect illumination at each reconstructed 3D point, without complex inverse global illumination computations. The coarse geometry is however unreliable for sun illumination that typically contains highfrequency features such as cast shadows. We introduce a new parameterization of reectance with respect to sun visibility that we integrate in an optimization algorithm to robustly identify the 3D points that are in shadow. We developed an optimization inspired by mean shift [START_REF] Comaniciu | Mean shift: a robust approach toward feature space analysis[END_REF] where we use asymmetric regions of inuence and constrain the evolution of the estimates.

Image-guided propagation algorithms are typically used to propagate user scribbles [START_REF] Levin | Colorization using optimization[END_REF][START_REF] Bousseau | User-assisted intrinsic images[END_REF]; we show how to use these algorithms to propagate the illumination information computed at 3D points over the image pixels. Our approach generates intrinsic images of similar quality as scribble-based approaches, with only a small amount of user intervention for capture and calibration. In addition, our ability to separate sun, sky and indirect illumination opens the door for advanced image manipulations, as demonstrated in Fig. 1b-d.

To summarize, this paper makes the following contributions:

We show how to compute sky, indirect, and sun (ignoring cast shadows) illumination at automatically reconstructed 3D points, using incomplete and imprecise geometry and a small set of input images.

We introduce an algorithm to reliably identify points in shadow based on a new parameterization of the reectance with respect to sun visibility.

Our algorithm compensates for the lack of accurately reconstructed and complete 3D information.

We show how to propagate reectance, sun, sky and indirect illumination to all pixels in an image, without user intervention or involved inverse global illumination computation. We achieve this by using the illumination values computed at 3D points as constraints for image propagation algorithms.

After a presentation of previous work, the denition of our image formation model and a description of capture, the structure of our paper follows the three contributions described above.

Inria 2 Related Work

Intrinsic images. Several methods have been proposed to estimate reectance and illumination from a single image. This decomposition is severely ill-posed and can only be solved with additional knowledge or assumptions about the scene content. The Retinex algorithm [START_REF] Land | Lightness and retinex theory[END_REF][START_REF] Berthold | Robot Vision. McGraw-Hill Higher Education[END_REF] assumes smooth illumination and piece-wise constant reectance, while Li and Yeo [START_REF] Li | Intrinsic image decomposition using a local and global sparse representation of reectance[END_REF] assume that neighboring pixels with similar chromaticity have the same reectance and that the image is composed of a small set of reectances. Tappen et al. [START_REF] Tappen | Recovering intrinsic images from a single image[END_REF] train a classier to discriminate image derivatives due to reectance and illumination, and Shen et al. [START_REF] Shen | Intrinsic image decomposition with non-local texture cues[END_REF] introduce texture constraints to ensure that pixels with similar texture share the same reectance. These various approaches produce encouraging decompositions on isolated objects, as evaluated by the ground-truth dataset of Grosse et al. [START_REF] Grosse | Ground-truth dataset and baseline evaluations for intrinsic image algorithms[END_REF]. However, the automatic decomposition of complex outdoor

images remains an open challenge, in part because most existing methods assume monochrome lighting while outdoor scenes are often lit by a mixture of colored sun and sky light.

Weiss [START_REF] Weiss | Deriving intrinsic images from image sequences[END_REF] demonstrates how multiple images of a scene under dierent illuminations can be factored into a reectance image and a set of illumination images. Sunkavalli et al. [START_REF] Sunkavalli | Factored time-lapse video[END_REF] decompose similar image sequences of outdoor scenes into a shadow mask and images illuminated only by skylight or sunlight.

By also capturing an environment map at multiple times of the day, Matusik et al. [START_REF] Matusik | Progressivelyrened reectance functions from natural illumination[END_REF] estimate the reectance eld of an outdoor scene that can then be used for relighting. These methods assume a xed viewpoint and varying illumination (i.e., timelapse sequences), while our method relies on images captured under dierent viewpoints and xed illumination. The main advantage of our capture approach is to reduce the acquisition time to a few minutes while timelapses typically require at least several hours of capture to cover as many lighting directions as possible. Most related to our capture strategy is the system of Melendez et al. [START_REF] Melendez | Relightable Buildings from Images[END_REF] that relights buildings reconstructed from multiple photographs. However their method necessitates the additional acquisition of ash / no-ash image pairs to capture a material exemplar for every material of a building. Users then need to associate a material exemplar to each texture region of the reconstructed building.

Bousseau et al. [START_REF] Bousseau | User-assisted intrinsic images[END_REF] and Shen et al. [START_REF] Shen | Intrinsic images using optimization[END_REF] rely on user scribbles to disambiguate reectance from illumination while Okabe et al. [START_REF] Okabe | Single-view relighting with normal map painting[END_REF] propagate a sparse set of user-specied normals over the image. Given an intrinsic image decomposition, Carroll et al. [START_REF] Carroll | Illumination decomposition for material recoloring with consistent interreections[END_REF] propose a user-assisted decomposition to isolate the indirect contribution of each colored material in the scene. Our method shares similarities with these user-assisted approaches but we propagate illumination values automatically computed from a sparse point cloud instead of user indications.

Inverse rendering. Inverse rendering methods [START_REF] Yu | Recovering photometric properties of architectural scenes from photographs[END_REF][START_REF] Yu | Inverse global illumination: recovering reectance models of real scenes from photographs[END_REF][START_REF] Loscos | Interactive virtual relighting of real scenes[END_REF][START_REF] Debevec | Estimating surface reectance properties of a complex scene under captured natural illumination[END_REF] recover the reectance and illumination of a scene by inverting the rendering equation. These methods require an accurate 3D model of the scene that is either modeled manually or acquired with expensive laser scanners and then solve an often costly global illumination system. In contrast our method is robust to incomplete geometry and handles sparse point clouds automatically reconstructed from a few photographs of the scene. In addition, we are able to estimate indirect illumination at every pixel without the need to solve inverse global illumination.

Nonetheless, inverse rendering methods allow applications that are beyond the RR n°7851 scope of this paper, such as free viewpoint navigation and dynamic relighting.

Image-based propagation. Image guided interpolation methods have been introduced by Levin and Lischinski to propagate colors [START_REF] Levin | Colorization using optimization[END_REF] and tonal adjustments [START_REF] Lischinski | Interactive local adjustment of tonal values[END_REF]. In this paper we use the propagation algorithm of Bousseau et al. [START_REF] Bousseau | User-assisted intrinsic images[END_REF] that was originally designed to propagate user indications for intrinsic image decompositions. This algorithm is inspired by the matting Laplacian of Levin and colleagues [START_REF] Levin | A closed-form solution to natural image matting[END_REF] that has been used to decompose an image into a foreground and background layer, and to recover white balanced images under mixed lighting [START_REF] Hsu | Light mixture estimation for spatially varying white balance[END_REF].

Shadow removal. Our work is also related to shadow removal methods [START_REF] Finlayson | Removing shadows from images[END_REF][START_REF] Finlayson | Intrinsic images by entropy minimization[END_REF][START_REF] Mohan | Editing soft shadows in a digital photograph[END_REF][START_REF] Wu | Natural shadow matting[END_REF][START_REF] Shor | The shadow meets the mask: Pyramidbased shadow removal[END_REF]] that aim at identifying and removing cast shadows in an image, either automatically or with user assistance. While our method is also able to identify cast shadows, our main goal is to extract a reectance image, as well as smooth illumination variations. We also separate the contribution of sun, sky and indirect illumination, which enables novel image manipulations.

Image Formation Model

We assume Lambertian surfaces and model the image values at each pixel as the product between the incident illumination and the object reectance R.

Formally, the radiance I towards the camera at each non-emissive, visible point corresponding to a pixel is given by the rendering equation [START_REF] Kajiya | The rendering equation[END_REF]:

I = R * Ω cos θ ω L(ω)dω (1)
where we integrate over the entire hemisphere Ω at the visible point, L(ω) is the incoming radiance in direction ω, θ ω is the angle between the normal at the visible point and direction ω. Capital bold letters represent RGB color values and * denotes per-channel multiplication.

For our purposes, we will separate out the incoming radiance into three components: the radiance due to the sun, that due to the sky and that due to indirect lighting. To simplify notation, we dene two subsets of the hemisphere:

Ω sky , i.e., the subset of directions in which the visible point sees the sky, and Ω ind the subset of directions in which another object is visible, and thus contributes to indirect lighting. We however explicitly represent the sun visibility v sun , rst because precise computation of v sun is necessary to capture sharp shadows, and second because estimating v sun robustly is one of our main contributions.

We can now re-write Equation 1:

I = R * v sun max(0, cos θ sun )L sun + Ω sky cos θ ω L sky (ω)dω + Ω ind cos θ ω L ind (ω) dω
where L sun , L sky and L ind are radiance from the sun, the sky and indirect lighting respectively, θ sun is the angle between the normal at the visible point and the sun modeled as a directional light source, and θ ω is the angle between the normal and the direction of integration ω over the corresponding hemisphere. The scalar v sun ∈ [0, 1] models the visibility of the sun (0 for completely hidden, 1 for fully visible).

Inria

We next dene the following simplied quantities at each pixel:

S sun = v sun max(0, cos θ sun )L sun = v sun Ŝsun (2)
S sky =

Ω sky cos θ ω L sky (ω) dω (3)
S ind =

Ω ind cos θ ω L ind (ω) dω. ( 4
)
where Ŝsun corresponds to the sun illumination when cast shadows are ignored.

This allows us to dene a simplied image formation model, which we use from now on:

I = R * (S sun + S sky + S ind ) (5) = R * S total ( 6
)
where R is the object RGB reectance. S sun , S sky and S ind are the RGB incident illumination (or shading) from the sun, sky and indirect lighting respectively.

Our rst goal is to extract the reectance R and illumination S total from this image formation model. We demonstrate how to make this problem tractable by leveraging the sparse geometric information generated by multiview stereo algorithms to compute Ŝsun , S sky and S ind at each reconstructed 3D point (Section 5). We then introduce a new algorithm to estimate the sun visibility v sun precisely despite the approximate available geometry (Section 6), which will allow us to obtain all illumination components using image propagation (Section 7).

4

Capture and Reconstruction

Our method relies on a relatively lightweight capture setup composed of a digital camera (preferably on a tripod), a photographer's gray card and a simple reective sphere 1 to capture the environment map (see inset). No other special capture or measurement hardware is required.

Photography

We rst capture a few ordinary low-dynamic range photographs (LDR) which we use to perform approximate geometric reconstruction of the scene. This set of photographs should have a good coverage of the scene from dierent viewpoints and sucient overlap between neighboring viewpoints to facilitate multiview stereo reconstruction. The number of photographs required to obtain an acceptable reconstruction depends on the complexity of the scene and the presence of image features. We captured between 10 and 31 LDR photographs for each of the scenes presented in our results.

We then capture two high-dynamic range (HDR) images of the front and side of the reective sphere, placed in the scene, to obtain an angular environment 1 The reective sphere does not need to be as accurate as in the apparatus described in [START_REF] Debevec | Estimating surface reectance properties of a complex scene under captured natural illumination[END_REF]; in practice we used an inexpensive pétanque ball. map of the scene. We do this using the standard HDR assembly technique of Debevec et al. [START_REF] Debevec | Recovering high dynamic range radiance maps from photographs[END_REF].

We nally capture linear HDR images of the viewpoints for which we want to estimate the intrinsic images. Recall that we capture all images in one session, at a single time of day.

Geometry and Illuminant Calibration

local 3D point cloud using the PCA method of Hoppe et al. [START_REF] Hoppe | Surface reconstruction from unorganized points[END_REF]. We discard 3D points in too sparse or too degenerate local neighborhoods. We consider a neighborhood as degenerate when the rst singular value of the PCA is twice greater than the second singular value.

We recover a geometric proxy from the oriented point cloud using the scale space meshing method of Digne et al. [START_REF] Digne | Scale space meshing of raw data point sets[END_REF]. This automatic approach produces detailed accurate meshes in regions where the point cloud is dense, while leaving holes in areas where the point cloud is irregularly sampled. We use the automatic hole lling tool available in MeshLab [START_REF] Cignoni | Meshlab: an open-source 3d mesh processing system[END_REF] to further improve the reconstruction in those areas. We also experimented with the Poisson reconstruction of Kazhdan et al. [START_REF] Michael Kazhdan | Poisson surface reconstruction[END_REF]. Although this algorithm generates a closed surface with no holes, we found that it tends to produce bumpy surfaces due to the irregular sampling of the point cloud. We nevertheless used the Poisson reconstruction on the Rocks scene for which the point cloud is dense and uniform.

Some manual interaction is also required to calibrate the sun and sky illumination. We rst identify the sun position, orient the environment map and label sky pixels. We also take two photographs of a gray card, in the sun and in shadow, which are used to estimate the color transfer function of the reective sphere and the radiance L sun of the sun. The details of this process are described in the Appendix.

The output of our capture is: a moderately dense point cloud reconstruction of the scene which we will refer to as the PMVS points (Fig. 2b), a very approximate geometric proxy, which often contains signicant geometric errors (Fig. 2c), the direction and radiance of the sun and a correctly aligned and scaled HDR environment map containing the sky and distant indirect radiance (Fig. 4c).

Geometry-Based Computation

We describe here how to compute sun, sky and indirect illumination values for each PMVS point. These points have been generated using multiview stereo and also have normals (Section 4.2).

We rst compute sun illumination Ŝsun ignoring cast shadows, i.e., unoc- cluded sunlight. We already know the required quantities for this computation, i.e., the normal at the point and the sun radiance and direction. Treating sun visibility requires much higher precision than the one provided by the proxy, and is treated separately in Section 6.

We then compute sky and indirect illumination at each PMVS point. Fig. 3 illustrates this computation that we detail bellow. In a nutshell, we use the HDR environment map to compute both sky illumination and distant indirect illumination, while we compute the near-eld indirect illumination from the proxy geometry. Note however that we do not need to know the reectance of the proxy for this step, we instead use the captured photographs to recover the necessary outgoing radiance over the reconstructed geometry.

Assigning radiance to the proxy geometry. We rst assign radiance to each PMVS point by looking up its pixel values in each HDR image where it appears (Fig. 3). We assign the average value as the outgoing radiance of the point, which is assumed constant in all directions. We then assign the radiance Our approach shares similarities with the techniques described by Yu and Malik [START_REF] Yu | Recovering photometric properties of architectural scenes from photographs[END_REF], but while their method was designed for accurate geometry constructed manually, we handle sparse incomplete geometry reconstructed automatically.

Finally we compute an initial guess of the sun visibility v init sun at each PMVS point by tracing a ray against the proxy geometry in the direction of the sun.

Note however that this visibility test is very sensitive to errors in the reconstructed proxy, as illustrated in Fig. 2c. This fact underlines the importance of accurately estimating v sun and we next show how to rene this initial estimate.

Estimating Sun Visibility at PMVS Points

One key contribution of our approach is a novel algorithm for identifying visibil- where I is the RGB pixel value in the image we wish to process, S sky , S ind and Ŝsun are the illumination values of the corresponding PMVS point computed in Section 5, and the division is per-channel. With this parameterization, varying v sun in [0, 1] generates a curve of candidate reectances in RGB space. Our goal in this section is to nd for each PMVS point, the position on its candidate curve corresponding to its reectance R (or correspondingly visibility v sun ).

The intuition of our approach is that multiple PMVS points sharing the same reectance will generate intersecting curves in color space; their (shared) reectance will be the color where the candidate curves intersect. This is illustrated on Fig. 5a-b. By nding these intersections we can deduce the value of v sun for the PMVS point corresponding to each curve. However, imprecision in the capture process and in the geometry-based computation prevent the curves from perfectly intersecting in color space. In addition, multiple intersections can occur along a curve, which gives multiple candidates for the visibility term. We address these issues with a robust iterative procedure inspired by the mean shift algorithm.

Overview. Mean shift [START_REF] Fukunaga | The estimation of the gradient of a density function, with applications in pattern recognition[END_REF][START_REF] Comaniciu | Mean shift: a robust approach toward feature space analysis[END_REF] is a non-parametric mode-seeking algorithm that aims to locate the maxima of a density function, given a set of data points. First, a kernel (or window) is placed at each data point; it represents the region Inria of inuence of this point. In an iterative process, each kernel is then moved in a direction that increases the local density, computed as the weighted average of nearby data points. The process stops when all kernels have reached a stationary point (or mode).

In our approach, we dene an asymmetric region of inuence for each candidate curve. We use mean shift iterations to maximize an energy that measures the overlap among pairs of curves, and iteratively update the estimated reectances while constraining them to lie on their candidate curves. After convergence, for each curve we obtain the reectance (and corresponding visibility) that tends to maximize the number of PMVS points sharing a similar reectance, and correspondingly tends to minimize the number of reectances in the scene. This algorithm assumes that the scene is composed of a sparse set of reectances shared by multiple points, which is a common assumption recently used in image segmentation [START_REF] Omer | Color lines: Image specic color representation[END_REF] and white balance algorithms [START_REF] Hsu | Light mixture estimation for spatially varying white balance[END_REF].

Region of inuence. Equation 9 denes the candidate reectances of a

PMVS point as a rational curve that is parameterized non-uniformly by v sun . In order to obtain a uniform parametrization we rst approximate each curve c as a piecewise linear curve in CIE L*a*b* space, and parameterize it by arc length from the shadowed end of the curve: t in [0, 1] so that R(v sun = 0) = R(t = 0).

We use this uniform parametrization to compute distances along the curve. We chose to work in CIE L*a*b* space because it denes a perceptually uniform distance metric.

We then dene the inuence of curve c on a point x of color space as:

A xc = h ⊥ d 2 ⊥ (x, c) σ 2 ⊥ h d 2 (x, c) σ 2 (10)
where

d 2 ⊥ (x, c) = x -proj 3D (x, c) 2
is the squared distance in color space between x and its projection on curve c (i.e., the distance perpendicular to the curve).

d 2 (x, c) = (t c -proj t (x, c)) 2
is the squared dierence between the position of the current reectance estimate along the curve t c and the arc length proj t (x, c) corresponding to the projection of x on the curve (i.e., the distance along the curve).

h ⊥ and h are Gaussian kernel proles with the form h(x) = e -x controlled by the standard deviations σ 2 ⊥ and σ 2 (larger σ values correspond to a wider region of inuence).

The rst term in Equation 10 compensates for curves that do not exactly intersect, by dening a region of inuence around the curve with a Gaussian fallo orthogonal to the curve. The second factor makes our algorithm robust to false intersections, i.e., intersections of curves that in fact do not share the same reectance. The Gaussian kernel h reduces the inuence of the false intersections far from t c as t c converges toward the most likely reectance.

The regions of inuence of two curves are illustrated in Fig. 5c.

RR n°7851

Laont & Bousseau & Drettakis Energy Formulation. We then dene an energy that measures the overlap between the regions of inuence of pairs of curves:

E = V   c∈C c ′ =c A xc A xc ′   dx (11)
where we integrate over the entire color space V . We evenly discretize the 3D color space into a set S of samples, and rewrite this energy as a discrete sum:

E total = s∈S E sample (s). ( 12
)
where the energy of a sample E sample (s) accumulates the contribution of each pair of curves intersecting nearby:

E sample (s) = c∈C c ′ =c A sc A sc ′
With this formulation, two curves will contribute to the energy of a sample only if they (almost) intersect near this sample.

Derivation of the energy gradient. We seek the positions of the reectance estimates along the curves t c for c ∈ C that maximize E total :

argmax tc E total (13)
i.e. that are located at the zeros of the gradient function.

The derivative of E total with respect to t c is given by: ∂E total

∂t c = 4 σ 2 s∈S c ′ =c (proj t (s, c) -t c ) A sc A sc ′ (14)
Setting Equation 14to 0 for all c ∈ C gives:

∂E total ∂v c = 0 ⇔   s∈S c =c ′ A sc A sc ′     t c - s∈S proj t (s, c) c ′ =c A sc A sc ′ s∈S c ′ =c A sc A sc ′   = 0
which is analogous to the form obtained in the mean-shift algorithm [START_REF] Comaniciu | Mean shift: a robust approach toward feature space analysis[END_REF].

Iterative process. We dene our iterative procedure recursively by computing at iteration i the weights A i sc (i.e., the regions of inuence) using the estimates t i c for all curves c ∈ C, then updating the estimates t i+1 c for each curve c using a weighted average of the projection of nearby samples on c:

t i+1 c = s∈S proj t (s, c) c ′ =c A i sc A i sc ′ s∈S c ′ =c A i sc A i sc ′ (15)
Inria Finally, we replace all the curves belonging to each cluster by one representative curve: the curve closest to the median of this cluster's feature vectors.

The initial guess of the sun visibility v init sun is assigned to the median visibility of each cluster. This process greatly simplies and cleans up the set of curves, and lowers the computational cost of our algorithm. As shown in Fig. 6d, the curves corresponding to groups of PMVS points with similar reectances (nearly) intersect after clustering.

After clustering, we run the optimization to maximize E total . Upon convergence, for each cluster we obtain the position of the estimated reectance along its representative curve t nal . This position value is assigned to all curves belonging to this cluster, from which we can deduce the nal estimated reectance of each PMVS point corresponding to these curves.

At the end of this process, we obtain a list of PMVS points for which the position along the curve t, the sun visibility v sun and the reectance R have been estimated.

Implementation. In practice we use a truncated kernel prole h ⊥ (Equation 10) so that h ⊥ (x) = 0 when x > λ (we use λ = 3). This means that each curve will only inuence a limited number of samples; for each curve c we can precompute the indices of these samples, as well as their orthogonal distance s -proj 3D (s, c) and position of their projection along the curve proj t (s, c).

For all our results, we evenly discretized the CIE L*a*b* color space into 60 × 36 × 36 samples, with L* in [START_REF] Carroll | Illumination decomposition for material recoloring with consistent interreections[END_REF]95], a* in [-25, 25], b* in [-25, 25]. We used a xed bandwidth for the curve clustering using 6D mean-shift clustering, as well as for the 3D mean-shift clustering for orientation separation.

The clustering is the most costly part of the algorithm and takes from 25 seconds to a few minutes with our Matlab implementation, depending on the number of PMVS points. Once the clustering has been performed, the iterative optimization takes around 10 seconds, which allowed us to test many parameters for σ 2 ⊥ and σ 2 .

We found that the algorithm produces good results for a wide range of parameters, and that the best values were scene-dependent. In our experiments we found that large values of σ 7

Estimating Illumination at Each Pixel

In previous steps, we have used multiview stereo methods to generate a sparse set of 3D points on which we compute the illumination values Ŝsun , S sky and S ind (Section 5) along with the visibility of the sun v sun (Section 6). We next show how to leverage image-guided propagation methods to assign reectance and illumination values to all pixels in the input photographs. We rst show how to propagate the total illumination S total and then describe a method to subsequently separate the contribution of each lighting component (i.e., sun, sky and indirect).

Image Guided Propagation

We use the intrinsic images algorithm of Bousseau et al. [START_REF] Bousseau | User-assisted intrinsic images[END_REF] that was designed to propagate user indications for separating reectance and illumination in a single image. This algorithm makes the intrinsic image decomposition tractable by assuming that the reectance values in a pixel neighborhood lie in a plane in color space. This planar reectance assumption translates to a set of linear equations so that the illumination image is expressed as the minimizer of a least-square energy

argmin Stotal ST total M Stotal (16)
where the vector Stotal stacks the pixels of the estimated illumination image and the matrix M encodes the planar reectance assumption (see the paper by Bousseau et al. [START_REF] Bousseau | User-assisted intrinsic images[END_REF] for the complete derivation). For colored illumination the optimization is solved for each color channel separately.

In their original paper, Bousseau et al. constrain the least-square system with user indications. Users can specify the value of Stotal over a few pixels or indicate that several pixels share the same illumination or reectance. In our approach we use instead the illumination and visibility estimated at PMVS points to constrain the optimization. We express these constraints as an additional quadratic penalty

argmin Stotal ST total M Stotal + w p∈P ( Sp total -S p total ) 2 (17)

RR n°7851

where P is the set of pixels covered by PMVS points and S p total = v p sun Ŝp sun + S p sky + S p ind their illumination values computed in Sections 5 and 6. The weight w controls the importance of the constraints, in practice we use w = 1 to give equal importance to the constraints and the propagation model. Fig. 7 shows the results of the image guided propagation. This result shows the power of our approach. We exploit the information provided by the sparse and imprecise geometry to provide constraints automatically to the algorithm of Bousseau et al. [START_REF] Bousseau | User-assisted intrinsic images[END_REF], thus eliminating the need for user scribbles. A visualization of the constraints for all scenes can be found in the supplemental material.

Light Source Separation

Given the estimated illumination image Stotal , we wish to separate the contribu- tion of each illumination component Ssun , Ssky and Sind . Inspired by previous work on white balance under mixed lighting [START_REF] Hsu | Light mixture estimation for spatially varying white balance[END_REF], we show how to express our light source separation as two successive matting problems. In a rst step we decompose the illumination into a sun component and a diuse component Sdi that includes the contribution of both sky and indirect lighting. In a second step we decompose the diuse component into its two terms.

We rst express each illumination term S as the product between a scalar intensity s = S and a chromaticity C = S/ S :

Stotal = s sun C sun + s sky C sky + s ind C ind = s sun C sun + s di C di . (18)
Denoting α = s sun /(s sun + s di ) we express the illumination image values at each pixel as a mixture between two values weighted by α:

Stotal = α(s sun + s di )C sun +(1 -α)(s sun + s di )C di (19) = αF + (1 -α)B (20)
We can now recover Ssun = αF and Sdi = (1 -α)B by solving a standard matting problem. We compute α at each PMVS point from the illumination values estimated in Sections 5 and 6. We then propagate α over the image Stotal using the matting Laplacian algorithm of Levin et al. [START_REF] Levin | A closed-form solution to natural image matting[END_REF]. Finally, given α at every pixel and the known S sun and S di at each PMVS point, we solve for the sun and diuse illumination images with the following least-square optimization: The rst term of this functional ensures that the decomposition explains the input illumination Stotal and follows the estimated ratio α. The second term adds a smoothness regularization on each component while the third term constrains the solution to agree with the illumination values computed at PMVS points. We used λ = 0.1 and w = 0.01 for all our results.

argmin F,B i∈I Si total -α i F i + (1 -α i )B i 2 + λ (F i x ) 2 + (F i y ) 2 + (B i x ) 2 + (B i y ) 2 + w p∈P α p F p -S p sun 2 + (1 -α p )B p -S p
As a second step we apply the same matting approach to further separate the diuse illumination Sdi as the sum of the sky illumination Ssky and the indirect illumination Sind . We show in Fig. 9 and supplemental materials the results of this decomposition. The overall decomposition takes around 90 seconds to compute for a 3.2 megapixel image, where 30 seconds are necessary to compute Stotal and 60 seconds to perform the two subsequent separations.

Inpainting This overall process gives a satisfactory decomposition in the scenes we have tested. However, a small residual border can remain around the hard shadow boundary (see Fig. 8 (left)), a common artifact of shadow removal [START_REF] Finlayson | Intrinsic images by entropy minimization[END_REF][START_REF] Wu | Natural shadow matting[END_REF]. We identify these shadow boundaries as pixels located on sun illumination discontinuities but not on normal discontinuities. We rst propagate normals from the PMVS points over the image using the method of Okabe et al. [START_REF] Okabe | Single-view relighting with normal map painting[END_REF]. We then run an edge detector over the sun illumination image and label edge pixels that do not correspond to edges in the normal image. We remove the labeled pixels and their immediate neighbors from the reectance, sky and indirect illumination images and use inpainting to ll in the holes. This post-process takes 60 seconds on average. Fig. 8 illustrates the result of this process.

8

Results and Discussion

In Fig. 9 we show results on four dierent scenes. For all results, we show reectance, sun, sky and indirect illumination layers. The number of photographs used for each scene is shown in Table 2. The estimated illumination at PMVS points and additional views for each scene are shown in the supplemental material.

The rst scene shows the base of a statue on a square; we took HDR images RR n°7851 The second scene is challenging, since the reconstruction process is unable to capture details of the railings and does not reconstruct the vegetation casting the main shadows (see also Fig. 2d). In this view, despite the intricate geometric Next we compare our approach to three state-of-the-art methods in Fig. 10.

All these methods take a single image as input. The user-assisted approach of Bousseau et al. [START_REF] Bousseau | User-assisted intrinsic images[END_REF] produces result of a similar quality to ours, but requires a signicant number of user indications (between 25 and 105 scribbles). The automatic method of Shen et al. [START_REF] Shen | Intrinsic image decomposition with non-local texture cues[END_REF] is able to extract most of the illumination variations but colored shadows remain in the reectance image. These colored residuals are due to the variation in color between sun and sky illumination, which violates the gray illumination assumption of this method. Residual shadows are also present in the reectance estimated with the automatic method of Shen et al. [START_REF] Shen | Intrinsic images using optimization[END_REF], as well as reectance residuals in the illumination image (tiled oor in the statue scene). Although this method can support user scribbles, the authors reported that scribbles did not improve the result signicantly in these examples. All the results of these comparisons have been kindly provided by the respective authors.

We also provide in Fig. 11 a comparison to the result obtained simply by using the proxy to compute sun, sky and indirect illumination using PBRT, and then inverting Equation 5 to obtain reectance. As we can see, for a given image, there are many holes due to the incomplete proxy. In addition, the sun illumination is completely erroneous, due to the lack of reconstruction of the surrounding objects. In contrast, our method correctly captures these sun shadows, and removes them in the reectance layer as well.

We illustrate applications of our decomposition in Fig. 1. We rst alter the reectance of the ground to insert a grati while maintaining consistent shadows (b). We then add a virtual object in the scene with consistent lighting and shadows (c). We used PBRT to render the dinosaur surrounded by the captured environment map, and its shadow cast on a horizontal plane. Finally we simulate a sunset by changing the color and intensity of each illumination component separately; in addition, our decomposition allows us to blur shadows without aecting the reectance of the scene (d). All these manipulations can be performed easily in image-editing software with layer support; the supplementary video shows how we created the images in Fig. 1b-d with Adobe Photoshop.

Discussion A drawback of our current method is the need for the reective sphere to capture the environment map. We could investigate tting a sky model to sky pixels visible in the other input photographs, similar to [START_REF] Yu | Recovering photometric properties of architectural scenes from photographs[END_REF], and extracting information about non-sky surroundings in a similar manner. We RR n°7851 will investigate ways to simplify the calibration process of the sphere so that it can be done once and subsequently used for any scene.

Since we were interested in separating sun lighting from other sources, we have not shown overcast scenes. There is no fundamental reason that our method would not work with overcast scenes; the decomposition would simply rely on the S sky and S ind values and ignore the sun illumination component.

Our method is robust to holes in the proxy, since in all visibility calculations rays that do not hit the proxy will hit the environment map and get a plausible color. However, our method can fail if the initial guess for sun visibility is completely wrong. In the supplemental material we show a case where a spurious object appeared in a single view, and was thus not reconstructed at all, resulting in this type of failure. Similarly, if objects have very dark reectance the PMVS reconstruction procedure does not provide a sucient number of points, resulting in errors. Our optimization for sun visibility exploits the redundant information provided by points of the same reectance that are in sun light and shadow. Nevertheless our method can also handle reectances that are only in light or shadow. Most often in such case the proxy initialization will result

in the correct answer. The only case which could potentially cause errors is when other reectance curves incorrectly intersect with or inuence the curve of this material. Finally, the candidate curves in our sun visibility method may not intersect if the illuminant calibration is too imprecise. We are investigating ways to improve calibration as well as a way to automate the process.

Conclusion

We have presented a method to estimate rich intrinsic images for outdoor scenes.

In addition to reectance, our algorithm generates a separate image for the sun, sky and indirect illumination. Our method relies on a lightweight capture (10-31 photographs in the scenes shown here) to estimate a coarse geometric RR n°7851

representation of the scene. This geometric information allows us to estimate illumination terms over a sparse 3D sampling of the scene. We then introduce an optimization algorithm to rene the inaccurate estimations of sun visibility. While incomplete, we demonstrate that this sparse information provides the necessary constraints for an image-guided propagation algorithm that recovers the reectance and illumination components at each pixel of the input photographs.

Our intrinsic image decomposition allows users of image manipulation tools to perform consistent editing of material and lighting in photographs. An interesting direction of future work is to adapt the method of [START_REF] Carroll | Illumination decomposition for material recoloring with consistent interreections[END_REF] to provide an alternative way of computing the illumination components. It will be interesting to investigate the tradeos between complexity, speed and quality obtained from such an adaptation compared to our approach.

Our current method works independently for each view. Enforcing consistent intrinsic image properties between views is an interesting future research direction. With such consistency, our method will open the way for dynamically relightable environments and free-viewpoint navigation for Image Based Rendering systems [START_REF] Buehler | Unstructured lumigraph rendering[END_REF], in addition to the applications demonstrated in this paper.

An important step for this goal is the generation of plausible shadow motion in the sun illumination layer in order to simulate moving light sources.

  Décomposition en images intrinsèques riches à partir de plusieurs vues d'une scène extérieure Résumé : Nous présentons une méthode capable de décomposer les photographies d'une scène en trois composantes intrinsèques la réectance, l'illumination due au soleil, l'illumination due au ciel, et l'illumination indirecte. L'extraction d'images intrinsèques à partir de photographies est un problème dicile, généralement résolu en utilisant des méthodes de propagation guidée par l'image nécessitant de multiples indications utilisateur. Des méthodes récentes en vision par ordinateur permettent l'acquisition facile mais approximative d'informations géométriques d'une scène à l'aide de plusieurs photographies selon des points de vue diérents. Nous développons un nouvel algorithme qui nous permet d'exploiter cette information bruitée et peu able pour automatiser et améliorer les algorithmes d'estimation d'images intrinsèque par propagation. En particulier, nous développons une nouvelle approche par optimisation an d'estimer les ombres portées dans l'image, en peaunant une estimation initiale obtenue à partir des informations géométriques reconstruites. Dans une dernière étape nous adaptons les algorithmes de propagation guidée par l'image, en remplaçant les indications utilisateurs manuelles par les données d'ombre et de réectance déduite du nuage de points 3D par notre algorithme. Notre méthode permet l'extraction automatique des images intrinsèques à partir de multiples points de vue, permettant ainsi de nombreux types de manipulations d'images. Mots-clés : Images intrinsèques, propagation guidée par l'image, stereo multivues, algorithm mean-shift

Figure 1 :

 1 Figure 1: Starting from multiple views of the scene (a), our method decomposes photographs into four intrinsic components the reectance (e), the illumination due to sun (f ), the illumination due to sky (g) and the indirect illumination (h). Each intrinsic component can then be manipulated independently for advanced image editing applications (b-d).

  RR n°7851 (a) Input photograph (b) Reconstructed PMVS points (c) Proxy shadows (d) Shadow conguration

Figure 2 :

 2 Figure 2: Multiview geometry provides us with a sparse oriented point cloud (b, red pixels). Compared to the original image (a), the initial guess of the sun shadows computed using the proxy reconstructed from the point cloud is very inaccurate due to geometric errors, and incomplete (c). In particular, a bush casting a large shadow on the oor (d) is not reconstructed.

  of the closest PMVS point to each vertex of the proxy, and interpolate these RR n°7851 (a) Input photograph (b) Colored proxy (c) Environment map

Figure 4 :

 4 Figure 4: We estimate sky and indirect illumination at PMVS points using (b) a colored proxy, obtained by projecting and interpolating the radiance of PMVS points on the reconstructed geometry, and (c) an environment map, shown in a latitude-longitude parameterization where the red curve separates pixels contributing to sky illumination and those contributing to distant indirect illumination.

  ity v sun with respect to the sun for each PMVS point. From our image formation model in Equations 2 and 5 we express the reectance at each PMVS point as RR n°7851 (a) Selected PMVS points (b) Candidate (c) Regions of inuence (d) Regions of inuence reectance curves at initialization after one iteration

Figure 5 :

 5 Figure 5: Multiple PMVS points sharing the same reectance will generate intersecting curves in color space. (a) We selected two PMVS points with the same reectance but dierent illuminations (red and blue squares). (b) The corresponding candidate reectance curves (nearly) intersect at one end, which corresponds to the reectance of both points. Diamond markers on the curves correspond to the initial guess for visibility (randomly set in this example). (c) Each curve aects a region of the surrounding color space, according to Equation 10 . Regions of inuence at initialization are illustrated as isosurfaces with varying opacity. Regions closer to the curve and the current visibility estimate are more aected. (d) After one iteration, the visibility estimates have moved towards the intersecting end of the curves, increasing the overlap between the regions of inuence (i.e., the energy E total in Equation12).

Figure 7 :

 7 Figure 7: Separation of Reectance and Total Illumination. The optimization procedure described in Equation 17 enforces illumination constraints at PMVS points (b) and propagates them to all pixels, in order to separate an input photograph (a) into total illumination (c) and reectance (d).

di 2 where

 2 I is the image domain, P is the set of pixels covered by PMVS points and F x , F y , B x and B y are the x and y derivatives of F and B computed with nite dierences. Inria (a) Before inpainting (b) After inpainting

Figure 8 :

 8 Figure 8: The reectance, sky illumination (here, close up) and indirect illumination estimated by our algorithm can contain residual variations along hard shadow boundaries (a). We use inpainting to remove these artifacts (b).

Figure 9 :

 9 Figure 9: Results of our decomposition. We adjusted the brightness of images for illustration purpose (the scaling factors used can be found in the supplementary materials). For each scene, the sun illumination is usually much more intense than sky illumination on average, and the sky illumination is more intense than indirect illumination on average.

  Inria conguration, shadows are successfully removed from the reectance layer, and the three other layers show good results. In the third scene, a staircase is shown. There are some residual artifacts at the shadow boundaries because the vegetation moved in the breeze during HDR acquisition. In the fourth scene, a rock wall is shown. Notice how the indirect layer well represents lighting in the cracks between rocks where neither sun nor skylight is present. The third and fourth scene also show that our method succeeds in creating plausible reectance, sun, sky and indirect layers.

11 :

 11 Figure 11: Comparison between the decomposition estimated directly from the geometric proxy (rst row), and our results (second row). We obtain the proxy reectance by dividing the input image by the sum of the illumination components. Holes and inaccuracies in the proxy translate to artifacts and residual shadows in the reectance.

2

  ⊥ can compensate for calibration errors that prevent reectance curves to perfectly intersect, while large values of σ 2 can compensate for the erroneous initialization of v init sun provided by the approximate proxy. However, σ values should remain small enough to prevent curves of dierent materials to inuence one another. Table1summarizes the parameter values used for the examples in this paper.

	Statue	Rocks	Ramp	Stairs
	(σ 2 , σ 2 ⊥ ) (0.1, 1) (0.01, 1) (0.3, 1) (0.01, 20)
	(0.1, 5) (0.01, 5) (0.3, 5)

Table 1 :

 1 Sets of region of inuence parameters used on the four scenes. When specied, the two sets of parameters produced equally good results.

	Inria

Table 2 :

 2 Number of photographs captured for the geometric reconstruction, for each scene in this paper.

We apply Structure from motion using Bundler[START_REF] Snavely | Photo tourism: exploring photo collections in 3D[END_REF] and the patch-based multiview stereo (PMVS) algorithm[START_REF] Furukawa | Accurate, dense, and robust multiview stereopsis[END_REF] on the set of LDR+HDR photographs using the respective authors' publicly available implementations 2 . The result of this process is an oriented point cloud (3D positions and normals), and calibrated cameras (extrinsic and intrinsic parameters). The process also returns whether each point is visible from each camera. We rectify the images for radial distortion, and we consider only visible PMVS points in each image.We found that the PMVS normals were too noisy for illumination computation. We thus estimate normals at each 3D point by tting a plane on the[START_REF] Bousseau | User-assisted intrinsic images[END_REF] Available at http://phototour.cs.washington.edu/bundler/ and http://grail.cs.washington.edu/software/pmvs/, as of 2011/10/10 Inria

Laont & Bousseau & DrettakisWe nally align the environment map and sun with the reconstructed scene.

Appendix: Illuminant Calibration

In this appendix, we describe the details of the illuminant calibration step for the sky and the sun.

First, because our model separates sun light from sky light, we need to remove sun pixels from the environment map. We use inpainting to ll-in the saturated sun pixels from their neighbors and dene the sun position as the barycenter of the inpainted pixels. Since our model also separates sky light from indirect light, we use a standard color selection tool to label sky pixels that will contribute to the sky illumination, while other pixels (building, trees) will contribute to indirect lighting.

Second, the environment map only captures a scaled version of the scene lighting since the sphere is not perfectly specular. We need to compensate for this scaling factor.

In our system the environment map is used to compute both the sky illumination S sky and part of the indirect illumination S ind , the other part being computed from the geometric proxy (see Section 5 for more details). We estimate the color transfer function of the reective sphere K (represented as a RGB vector) by taking a photograph of a neutral gray card with known reectance placed in sun shadow. We intentionally place the card at a position where we expect its geometry to be well reconstructed. From the image formation model we have

where S env . denotes the illumination terms computed from the environment map and S proxy . the ones computed from the geometric proxy. We use this equation to solve for the unknown K.

We similarly recover the sun radiance L sun by taking a second picture of the gray card placed in sunlight (v sun = 1). From this picture we have

where L sun is the only unknown.

RR n°7851

To do so we manually mark a vertical edge of the reconstructed geometry and rotate the environment map and sun until the cast shadow of the virtual edge is aligned with that in the photograph.

We could envisage other approaches for this calibration. Cabral et al. [START_REF] Cabral | Relighting photographs of tree canopies[END_REF] use a compass to compute the environment map orientation with respect to the view direction in each photograph. In Yu and Malik [START_REF] Yu | Inverse global illumination: recovering reectance models of real scenes from photographs[END_REF], a set of photographs of the horizon are taken, and a sky model is tted. This involves taking several photographs, as opposed to the only two we take of the sphere, and also requires the use of two neutral lens lters to estimate sun radiance. In addition, tting parameters of sky models can require specic tuning for relighting applications, as was the case for morning and evening sky in the work of Cabral et al. [START_REF] Cabral | Relighting photographs of tree canopies[END_REF].

Inria