Georges Nader
email: georges.nader@inria.fr

Instant Transport Maps on 2D Grids

Keywords: CCS Concepts:, Computing methodologies → Image processing; Optimal Transport, Nonlinear solver, Re-meshing, Sampling

Fig. 1. Our fast mass-transport solver enables many applications such as adaptive sampling, surface remeshing, heightfield morphing and caustic design with interactive performance. From left to right: a painting of Van Gogh (A Wheatfield with Cypresses), Max-Planck 3D model courtesy of Max-Planck Institut für Informatik, and volcano heightmaps courtesy of University of Otago.

In this paper, we introduce a novel and extremely fast algorithm to compute continuous transport maps between 2D probability densities discretized on uniform grids. The core of our method is a novel iterative solver computing the L 2 optimal transport map from a grid to the uniform density in the 2D Euclidean plane. A transport map between arbitrary densities is then recovered through numerical inversion and composition. In this case, the resulting map is only approximately optimal, but it is continuous and density preserving. Our solver is derivative-free, and it converges in a few cheap iterations. We demonstrate interactive performance in various applications such as adaptive sampling, feature sensitive remeshing, and caustic design.

INTRODUCTION

The notion of optimal transport (OT) was first suggested in 1781 by Monge as an engineering task to move a pile of sand into a hole with minimal effort. This problem has been formalized as the problem of finding a transport map T between two probability measures µ and ν that minimizes a certain cost c [START_REF] Kantorovich | On the transfer of masses (in Russian)[END_REF]]. Since then, optimal transport has appeared in various form in numerous application domains [START_REF] Villani | Optimal Transport: Old and New[END_REF]]. In the context of computer graphics, typical examples include color transfer [START_REF] Morovic | Accurate 3D Image Colour Histogram Transformation[END_REF], mesh generation [START_REF] Delzanno | An optimal robust equidistribution method for two-dimensional grid adaptation based on MongeśKantorovich optimization[END_REF]], BRDF design [START_REF] Bonneel | Displacement interpolation using Lagrangian mass transport[END_REF], sampling [START_REF] Fernando De Goes | Blue Noise Through Optimal Transport[END_REF][START_REF] Qin | Wasserstein Blue Noise Sampling[END_REF], shape morphing [START_REF] Lévy | A Numerical Algorithm for L2 Semi-Discrete Optimal Transport in 3D[END_REF]], inter surface mapping [START_REF] Mandad | Variance-Minimizing Transport Plans for Inter-surface Mapping[END_REF], area preserving maps [START_REF] Zhao | Area-Preservation Mapping using Optimal Mass Transport[END_REF]], caustic design [START_REF] Schwartzburg | Highcontrast Computational Caustic Design[END_REF], centroidal convex decomposition [START_REF] Shi-Qing Xin | Centroidal Power Diagrams with Capacity Constraints: Computation, Applications, and Extension[END_REF], etc.

All these applications thus rely on the availability of efficient solvers of the underlying OT problem. Its choice and design heavily depend on the domain (e.g., Euclidean space versus an arbitrary manifold), the cost function, and the specificities of the application. The L 2 norm is by far the most used and well understood cost function. This choice provides welcome properties such as smoothness and unicity of the solution.

Previous Work on OT Solvers

In this context, several very different solving strategies have been developed, and we refer to recent surveys [START_REF] Lévy | Notions of optimal transport theory and how to implement them on a computer[END_REF][START_REF] Peyré | Computational Optimal Transport[END_REF][START_REF] Solomon | Optimal Transport on Discrete Domains[END_REF]] for a more detailed overview.

The original formulation, as proposed by [START_REF] Kantorovich | On the transfer of masses (in Russian)[END_REF], consists in first relaxing the problem of computing a bijective map T : Ω → Ω by rather computing a transport plan π : Ω × Ω → R. Then, by discretizing both domains as Dirac masses, the problem boils down to a standard linear programming (LP) task. Since the number of unknowns is the square of the number of Dirac masses, such an approach has for a long time been limited to very small problems with up to a few thousands of Dirac points. More recently, through the use of continuous kernels as entropic regularizers and sparse representation of the transport plan, costly LP solver can be advantageously replaced by faster iterative scaling procedures [START_REF] Cuturi | Sinkhorn Distances: Lightspeed Computation of Optimal Transport[END_REF][START_REF] Solomon | Convolutional Wasserstein Distances: Efficient Optimal Transportation on Geometric Domains[END_REF]]. Such approaches do not produce a bijective map but a more general transport plan, which might be a strength for some applications, but also a limitation for all applications seeking for a bijective map.

Another popular class of solvers computes the transport map between continuous and point-wise measures through a convex optimization of the weights of a power-diagram [START_REF] Aurenhammer | Minkowski-Type Theorems and Least-Squares Clustering[END_REF]]. The convergence of such semi-discrete solvers can be greatly improved through multiresolution [START_REF] Lévy | A Numerical Algorithm for L2 Semi-Discrete Optimal Transport in 3D[END_REF][START_REF] Mérigot | A Multiscale Approach to Optimal Transport[END_REF] or Newton iterations [START_REF] Fernando De Goes | Blue Noise Through Optimal Transport[END_REF]. Similar to the previous discrete setting, one difficulty is that when applied to a pair of continuous source-target, one has to be discretized as Dirac masses and a continuous transport map has to be approximated using some kind of interpolation [START_REF] Schwartzburg | Highcontrast Computational Caustic Design[END_REF].

A third line of work proposes to directly compute the transport map between continuous densities. In this continuous setting, the OT problem boils down to solving a Monge-Ampère equation which is a nonlinear elliptic PDE. The non-linearity can be relaxed using computational fluid dynamics at the expense of an additional virtual time [START_REF] Benamou | The Monge-Kantorovitch mass transfer and its computational fluid mechanics formulation[END_REF][START_REF] Papadakis | Optimal Transport with Proximal Splitting[END_REF]]. More recently, a finite-difference scheme with accurate boundary condition and Newton iterations has been proposed by [START_REF] Benamou | Numerical Solution of the Optimal Transportation Problem Using the Monge-Ampère Equation[END_REF].

Whereas all the solvers of each category relies on very different discretizations, formulations, and numerical tools, it is surprising to see that the fastest solvers of each category exhibit about the same performance in terms of speed: processing typical images requires order of minutes. Achieving such a level of performance made several novel applications of OT possible. However, without significantly faster OT solvers, those novel tools can hardly be integrated within end user applications where the parameters of the given tool often need to be tuned in an interactive manner to explore the space of solutions. For the same reason, such tools can hardly be used as building blocks of higher level processing that would involve the repeated computation of OT maps.

Contribution

In order to achieve the 2 orders of magnitude speedup required to enable interactive applications, we first observe that for some graphics applications such as sampling, area-preserving mapping and re-meshing, one of the densities can be assumed to be uniform, thus opening the door for severe simplification and optimization opportunities. Moreover, in the more general case of a pair of nonuniform densities, we argue than many applications do not require a truly optimal solution regarding the transport cost as long as: 1) the densities are preserved and 2) that the transport map remains as smooth as the optimal L 2 norm solution. These two observations lead to the first central idea of this paper, which consists in focusing on computing as quickly as possible optimal L 2 transport maps between a non-uniform source and a uniform target. A more general density preserving and smooth map between arbitrary densities can then be obtained through numerical inversion and composition.

Under this hypothesis, we propose a very fast numerical method for computing the L 2 optimal transport map between a 2D Euclidean grid (i.e., an image) toward a uniform density. Our approach lies in the third category of solvers working in the continuous setting that has the main advantage of directly computing a transport map suitable for composition and fast solution space exploration for interpolation problems (e.g., color transfer, shape morphing, etc.). To achieve very high performance we make several technical contributions.

Firstly, we show that the resulting nonlinear Monge-Ampère equation can be advantageously decomposed to leverage a fixed Laplacian-like linear operator whose contribution dominates the equation to be solved. This decoupling will allow us to design an iterative solver based on a single prefactorizaton and very cheap iterations.

Secondly, we show how to adequately discretize our PDE and differential operators so that:

(1) The input density grid is properly interpreted as a piecewise constant function without implicit smoothing.

(2) The discrete equations are consistent with the initial energy preserving constraints.

(3) The discrete problem boils down to a problem of simultaneous nonlinear equations instead of a more classical minimization of a nonlinear least-square energy. The first two properties will lead us to differential operators that depart from standard finite difference/element discretizations on regular grids. The third property might look insignificant at first, but as we will show it plays a fundamental role in the efficiency of our approach.

Last but not least, we propose a novel acceleration heuristic for solving nonlinear simultaneous equations that takes inspiration from nonlinear conjugate-gradient methods for nonlinear energy minimization.

Overall, each iteration of our solver requires only a few evaluations of residual vectors (about 3 or 4) only, and does not even require the computation of the Jacobian matrix nor gradient vectors. As a result each iteration is extremely fast, and as we will show our novel solver outperforms standard Newton iterations while avoiding the numerical difficulties of Newton iteration in non-convex area or in case of badly conditioned Jacobians. This also makes our approach very simple both conceptually and in terms of implementation.

As illustrated in Figure 1, we demonstrate our solver on various applications involving either one uniform density (remeshing, sampling) or a pair of non-uniform densities (shape morphing, goal based caustics). For all these applications, the solution has been computed in about 1s on a single core CPU.

In this paper we will focus on 2D uniform grids, leaving the extension of our work to 3D and arbitrary domains for future work.

BACKGROUND 2.1 Solving Simultaneous Nonlinear Equations

As sketched in the introduction, our discretization will eventually lead us to the problem of solving a set of nonlinear simultaneous equations of the form: r(x) = 0 , (1) with as many unknowns x as equations r. The literature on this class of problems is much more scarce than on the optimization of nonlinear energies which receive of a lot of attention even in the computer graphics community [START_REF] Shahar | Accelerated Quadratic Proxy for Geometric Optimization[END_REF][START_REF] Rabinovich | Scalable Locally Injective Mappings[END_REF].

Of course, such a problem can easily be cast as an optimization problem by considering the L 2 norm of the residual vector r:

e (x) = 1 2 ∥r(x)∥ 2 , (2)
which could then be blindly minimized using a classical solver such as Newton, quasi-Newton, or nonlinear conjugate-gradient methods [START_REF] Nocedal | Numerical Optimization[END_REF]. This would, however, be a lose idea because we would loose information in the process, and the underlying squaring would make the problem much more anisotropic and thus difficult to solve [START_REF] Broyden | A Class of Methods for Solving Nonlinear Simultaneous Equations[END_REF]].

Gauss-Newton solvers can exploit the specific structure of e, but the most classical approach to solve (1) is through Newton root finding algorithm, which is directly derived from a first-order approximation of r [START_REF] Householder | Principles of Numerical Analysis[END_REF]]. In practice it is often coupled with a line-search algorithm leading to the following iterative algorithm:

d k = -J r (x k ) -1 r k (3) x k +1 = x k + α k d k (4)
where k denotes the iteration number, d k is the search direction, J r is the Jacobian of r, r k is a shortcut for r(x k ), and α k ∈ [0, 1] is found to prevent divergence in case of a poor local approximation.

The line-search algorithm finding α k requires an error function as e in equation ( 2). It is classically implemented using some variants of the Wolfe conditions [START_REF] Nocedal | Numerical Optimization[END_REF]. This Newton method suffers from two severe difficulties. The first one is the requirement of solving a new, non-symmetric, linear system of equations from scratch at each iteration, which is extremely costly especially with a large number of unknowns. This goes without saying that computing the Jacobian J r itself might be difficult and costly. The second one is that if the initial guess is not close enough to the solution, then the search direction is not necessarily a descent direction and the iterations will either stall or diverge. Even though the conditions for convergence are well known for decades [START_REF] Householder | Principles of Numerical Analysis[END_REF]], modifying the Jacobian to ensure convergence remains a tricky and costly task.

The most well known alternative to this Newton method are the two quasi-Newton rules proposed by Broyden [START_REF] Broyden | A Class of Methods for Solving Nonlinear Simultaneous Equations[END_REF]]. This approach follows the same algorithm than the above Newton iterations but using an approximation M of the inverse of the Jacobian that is updated at each iteration such that 1) it satisfies the secant equation between two successive solutions M k (r k -r k -1 ) = α k -1 d k -1 , and 2) the new approximate Jacobian (or its inverse) is the closest to the previous one. In practice, this boils down to rankone updates, and just like the famous L-BFGS method [START_REF] Nocedal | Numerical Optimization[END_REF], this strategy can easily be implemented with a bounded number of recurrence vectors. This approach makes only use of derivative information from the past. On the one hand this a strength because no additional information from the problem at hand has to be computed, but on the other hand all the differential information at hand might be obsolete when working on high-dimensional and highly anisotropic problems.

In Sections 4.1 and 4.2 we will devise a novel approach to this class of problem, which can be understood as either a conjugategradient like approach for simultaneous linear equation, or a novel quasi-Newton update with only one recurrence term. Our approach exhibits much higher convergence rate with only one additional evaluation.

Optimal Transport & the Monge-Ampère equation

When working with continuous domains for both the source and target densities, the OT problem can be described as follows. Let u and v be two probability densities supported on U and V respectively, where U and V are compact subsets of R d . Let M be the set of maps transporting the source u into the target v:

M = {T : U → V ; u = v (T )det(J T )} , ( 5
)
where J T is the Jacobian matrix of T . In this paper, we seek for the map T that minimize the L 2 cost of transportation:

c (T ) = 1 2 U ∥x -T (x)∥ 2 u (x) dx . ( 6
)
As shown by [START_REF] Brenier | Polar factorization and monotone rearrangement of vector-valued functions[END_REF], this problem is well posed and it admits a unique minimizer which is the gradient of a convex scalar potential:

φ : U → R , T (x) = ∇φ(x) . ( 7
)
Plugging equation ( 7) within the density preserving equation ( 5) yields the elliptic Monge-Ampère PDE:

det(H φ (x)) = u (x) v (∇φ(x)) , ∀x ∈ U , (8)
where H φ is a Hessian matrix of φ. Solving directly this equation in the general case is a challenging task, not only because it is highly nonlinear but also because we have to consider additional constraints such as the convexity of φ and the so called OT boundary condition ∇φ(U ) = V which is also highly nonlinear. To the best of our knowledge, the most advanced numerical solver for this equation has been proposed by [START_REF] Benamou | Numerical Solution of the Optimal Transportation Problem Using the Monge-Ampère Equation[END_REF].

The boundary condition is approximated through a signed distance field leading to an Hamilton-Jacobi equation. Both PDE are then discretized through an adequate finite-difference scheme enforcing for convexity of the solution and Newton iterations. Despite its high accuracy, handling typical size images requires order of minutes of computations.

OUR SIMPLIFIED OPTIMAL TRANSPORT MODEL

In this section, we show how to derive from the general Monge-Ampère optimal transport formalism a simplified problem by introducing adequate restrictions. The non-linearity of equation 8 comes from three terms:

(1) the determinant, (2) the quotient between u and v, (3) the composition of v with the unknown map T .

Based on this observation, our first and major simplification consists in imposing that the target density v is constant. Without lack of generality, we can further assume that the area of the domain Y of v is unit, in which case v = 1, and the PDE to solve reduces to:

det(H φ (x)) = u (x) , ∀x ∈ U , ( 9
) u i , j i j C i , j T (C i, j ) T ( x i+ 1 2 , j+ 1 2 ) h
Input density u Warped grid Fig. 2. In the discrete setting, the continous L 2 OT problem consists in displacing the vertices of the regular grid so that the area of each cell corresponds to the target density.

which is considerably simpler to solve, provided that the remaining nonlinear determinant can be efficiently handled as we will show in Section 4. At a first glance this seems to be a severe restriction, but as we will demonstrate in the result section, many applications naturally deal with a uniform density as either the source or target. For instance, to adapt a uniform distribution to a non-uniform density u, we compute the transport map T u→1 using our simplified problem and then apply its inverse T -1 u→1 numerically (Section 6). For applications requiring a non-uniform source u and target v, we can compute a density preserving map T u→v by 1) computing the maps T u→1 and T v→1 using our simplified problem, and 2) composing them numerically:

T u→v = T -1 v→1 • T u→1 . ( 10
)
Of course, this only gives an approximation of the optimal solution as the composed map T u→v is generally not a minimizer of the L 2 cost function ( 6), but it is straightforward to verify that it does provide a valid transport map (i.e. T u→v ∈ M), which we show is sufficient for many performance-demanding applications.

Our second simplification consists in imposing that both the source and target share the same support (i.e., U = V ), and U must be a rectangle so that it admits a simple discretization through a regular grid. This restriction is naturally achieved when working on images, otherwise it needs to be enforced by padding the bounding rectangle with zeros. This choice considerably simplifies the handling of the non standard OT boundary condition since it reduces to linear Neumann boundary conditions [START_REF] Delzanno | An optimal robust equidistribution method for two-dimensional grid adaptation based on MongeśKantorovich optimization[END_REF][START_REF] Haker | Optimal Mass Transport for Registration and Warping[END_REF].

More precisely, we assume that the input density u is given as a piecewise constant function defined over a 2D uniform grid consisting in n cells of size h × h with values u = (u 1 , . . . , u n ) T . With such a discretization, it is natural to assume that T is piecewise bilinear over each cell such that the image of u by T is a warped regular grid with quadrilateral cells (Figure 2). Integrating our density preserving equation det(J T ) = u over each cell leads to a system of n equations:

area(T (C i, j )) -h 2 u i, j = 0 , ∀ cell C i, j , (11)
which simply states that the area of the warped cells must be proportional to the source density u.

DESIGN OF OUR EFFICIENT SOLVER

In this section we devise an efficient solving strategy for our simplified PDE eq. ( 9). In order to efficiently handle the remaining nonlinear term (i.e., the determinant), we propose to write the transport map as a displacement by the gradient of a different scalar potential ψ : 12) Such a decomposition is rather common in the literature where -ψ is called the Kantorovich potential. It is related to the initial convex potential φ by φ(x) = 1/2∥x∥ 2 + ψ (x), and their Hessian are related by: H φ = I + H ψ . When working in R 2 , plugging this identity in our simplified PDE (9) leads, after some trivial refactorizations, to the following equation:

T (x) = x + ∇ψ (x) . (
∆ψ + det(H ψ ) + 1 = u . ( 13
)
This rewriting is one of the central ingredient of our fast solver. Indeed, the initial nonlinear determinant has now been decomposed as the sum of a linear Laplacian operator and another nonlinear determinant. Since the latter measures some kind of distortions of the displacement and that the displacement is expected to be locally smooth in most area, this expression is expected to be greatly dominated by the Laplacian term. This assumption is verified in Figure 3 on the Julia set example. For this rather complicated example, the integral of the magnitude of the Laplacian is about ×8 times larger than the one of the determinant. Such an insight immediately suggests a global/local alternate solving approach. To this end, let us introduce a dual variable ζ for the determinants that is initialized to zero. Let the index k denotes the iteration number. Then the general idea would be to iterate through the following two steps:

(1) Solve for ψ k assuming ζ k -1 is constant using the following Poisson equation:

∆ψ k = u -ζ k -1 -1 . ( 14
)
(2) Update the determinants ζ k := det(ψ k ). This strategy is appealing because the Laplacian operator can be prefactored only once so that the repeated linear problems ( 14) can be efficiently carried out during the iterations. However, since the problem to solve is not convex, this naive approach is unlikely to converge even after re-writing the steps as proper proximal operators. Moreover, the proximal operators associated to such steps involve significant computation overhead.

Before presenting our actual solver, let us precise the discretization of equation ( 13). As we will detail in Section 5, in order to yield a consistent discretization with as many unknowns as equations, we discretize the unknown function ψ on the dual grid of the input image with nodal values defined at the cell centers. Then, following the discretization made in equation ( 11) with one integrated equation per cell, we end up with a set of equations of the following form:

r(ψ ) = Lψ + h 2 (q(ψ ) + 1 -u) = 0 , ( 15
)
where r is the vector of residual, L is a sparse symmetric matrix corresponding to the linear Laplacian term, ψ is the vector of the unknowns, and q is the vector holding the determinants of the Hessian of ψ per cell. The actual form of L and q will be detailed in Section 5, whereas solving as efficiently as possible ( 15) is the subject of the next two subsections.

Our Fast Linearized Solver

Since we took care at building a proper system of simultaneous equations, our problem (15) could directly be solved using the Newton method and its variants described in Section 2.1. As we will show in the result section, in our case those methods are rather slow and even fails to converge for several inputs. We overcome these issues by taking advantage of the insights we developed at the beginning of this section, while casting the naive alternate strategy into a more robust line-search based algorithm similar to Newton iterations described in Section 2.1. To this end, the idea is to use the solution of equation ( 15) as a target (instead of directly picking it as the next solution ψ k +1 ), and then to seek for a good enough solution ψ k +1 between the current solution ψ k and this target. After some algebraic manipulations, this search direction d is more directly obtained as the solution of:

d k = -L -1 h 2 (q k + 1 -u) -ψ k = -L -1 (Lψ k + h 2 (q k + 1 -u)) = -L -1 r k . ( 16
)
The next solution ψ k +1 is obtained through a line-search using equation ( 4) and the same error function e defined in equation ( 2). Line-search algorithms are classically implemented using some variants of the Wolfe conditions [START_REF] Nocedal | Numerical Optimization[END_REF]]. However, those conditions involve the computation of the gradient of the error function e only to check whether the given approximate solution is acceptable or not. In our context, this is a waste of computation effort because the gradient of e is not needed for other purposes. We thus implemented an accurate line-search algorithm combining quadratic interpolations of the error function within the current search range, and a fallback to a golden-section search in case the interpolating quadratic is concave or its minimum leads to an invalid guess. We stop bracketing when the error is reduced and that the distance between two successive estimates of α is below some threshold t α (e.g., t α = 10 -2 ).

Algorithm 1: Our iterative solver with conjugate directions.

Initialization: build and factorize the matrix -L ψ 0 := 0, r 1 := 0 do 1: Solve

d k = -L -1 r k ▷ sec. (4.1) 2: Improve search direction, ▷ sec. (4.2) d k = d k + βd k -1 ▷ eq. (21) 3: Minimize ∥r(ψ k -1 + α k d k )∥ 2 using line-search, return new solution ψ k +1 = ψ k + α k d k and residual r k +1 . while ∥r k +1 ∥ 2 < t e /h 2
Our solver iterations are summarized in Algorithm 1 where we can for now ignore the step 2 (i.e., β = 0) that will be described later. We stop our iterations when the norm of the residual is below a given user defined tolerance t e that is scaled by the cell area h 2 to be agnostic from the grid resolution.

Let us emphasize that even though equation ( 16) resembles to Newton iterations with some arbitrary approximation of the Jacobian, it is actually stronger than that because it gives us a direction of projection towards the hyperplane of solutions defined by the Poisson equation ( 14).

Acceleration Through Conjugate Directions

The current version of our solver is particularly fast during the first iterations where L is a good approximation of the true Jacobian. However, for difficult inputs for which the determinants of ψ start to have a significant influence, the convergence is expected to slow down. In order to maintain a higher speed of convergence, we propose to take inspiration from nonlinear conjugate-gradient (NLCG) methods by updating the search direction with respect to the previous one as follows:

d k = d k + β k d k -1 , ( 17
)
with d k = -L -1 r k the initial search direction guess, and β a scalar coefficient appropriately chosen to yield a better search direction (Figure 4). However, NLCG methods have been developed to minimize nonlinear energies and none of the numerous heuristics that have been proposed so far makes sense in our context [START_REF] Hager | A survey of nonlinear conjugate gradient methods[END_REF]. We thus designed a novel heuristic dedicated to the problem of solving simultaneous equations. The general idea consists in enforcing two successive directions to be conjugate with

x k -1 x k x k +1 dk d k d k -1
Fig. 4. Illustration of one iteration of our solver accelerated with conjugate directions.

x k respect to the metric induced by r. In other words, we would like that the following relation hold:

x k x k -1 d k d k -1 5 e -3 d k
d T k J T k J k d k -1 = 0,
where we used the shortcut J k = J r (ψ k ) for brevity. This is easily accomplished by using the following formula for β:

β = - d T k-1 J T k J k d k d T k -1 J T k J k d k -1 . ( 18
)
At a first glance this formula is not very appealing as it involves the full Jacobian matrix at the current iteration. However, let us observe that a product of the form J k d is nothing else but the scaled derivatives of r in the direction d. We can thus efficiently compute those quantities using finite difference approximations:

J k d k ≈ r(ψ k + ε d k ) -r k ε = r k +ε -r k ε . ( 19
)
The derivatives along the previous search direction can directly be estimated for free from the previous residual as follows:

J k d k-1 ≈ r k -r k-1 α k -1 . ( 20
)
This leads to the following final formula for β:

β = - α k -1 (r k -r k -1 ) T (r k +ε -r k ) ε ∥r k -r k -1 ∥ 2 , ( 21
)
which, overall, requires only one additional evaluation of the residual that is cheap to compute. According to our experiments, we found that taking ε = α k -1 to estimate J k d k -1 as in eq. ( 20) led to slightly faster convergence than computing a more accurate estimate, likely because of the implicit smoothing of the error function it provides. A similar observation has been made by [START_REF] Broyden | A Class of Methods for Solving Nonlinear Simultaneous Equations[END_REF] with his quasi-Newton method. On the other hand, when estimating J k d k , since our error function e is extremely anisotropic, ε must not be chosen too large, and in practice we found that taking ε ∈ [10 -6 , 10 -2 ] has little impact on the convergence rate. We thus fixed it to 10 -3 in our implementation. Figure 4 depicts a typical step of this algorithm, which according to our experiments is almost always successful in finding the local minima in the given 2D subspace (see also supplemental material). In contrast, as depicted in Figure 5, the true Jacobian often yields very poor search directions from which only very small progress can be achieved.

∇ ψ ( x i+ 1 2 , j+ 1 2 ) ψ i , j
x 0

x 1 Fig. 6. Illustration of our dual grid discretization. Nodal values of ψ are defined at cell centers x i, j , and gradients are computed at cell's corners

x i + 1 2 , j+ 1 2
from the 4 neighbor nodal values.

DISCRETIZATION

This section details how we derived our discrete equation ( 15) from its continuous form ( 13), how to treat boundaries, and some related implementation details.

Let us start by showing how to properly construct the pseudo-Laplacian matrix L and vector of Hessian determinants q so that eq. ( 15) is rigorously equivalent to our initial area-based density preserving eq. ( 11). In other words, we must make sure that the following identity holds:

       . . . area(T (C i, j )) . . .        = Lψ + h 2 (q(ψ ) + 1) . ( 22
)
With such a goal in mind, it is pretty clear that using a standard firstorder finite difference/element scheme with nodes defined either on cell corners or on the dual grid to discretize our PDE ( 13) is doomed to fail. The main challenge to ensure a consistent discretization is to make sure that the discrete gradient and Laplacian differential operators are consistent to each other, that is, the Laplacian discretization must be based on the discrete gradient operator used to compute the displacements of the cell corners from ψ . To achieve highest performance as possible, those operators must further be based on as compact as possible stencils.

As already outlined, we achieve these goals using a finite-difference scheme based on the discretization of the potential ψ on the dual grid, that is, with one nodal value per cell center, as depicted in Figure 6.

Gradient & Laplacian

To compute the warped cell T (C i, j ) and its area we need to compute the image of each corner from the gradient of ψ :

T x i+ 1 2 , j+ 1 2 = x i+ 1 2 , j+ 1 2 + ∇ψ x i+ 1 2 , j+ 1 2 , ( 23
)
that we compute through averaged second-order central differences:

∇ψ x i+ 1 2 , j+ 1 2 :=       ψ i +1, j -1 +ψ i +1, j -ψ i, j -ψ i, j -1 2h ψ i, j +ψ i +1, j -ψ i, j-1 -ψ i +1, j -1 2h       . ( 24
)
We then define our consistent discrete Laplacian operator by composing the previous discrete gradient with an analogously defined divergence operator leading to the following 5-point stencil mask at cell i, j:

h 2 ∆ψ (x i,j ) := ψ i-1,j-1 +ψ i+1,j-1 +ψ i+1,j+1 +ψ i-1,j+1 -4ψ i,j 2 , ( 25
)
which is different than the usual 5-point Laplacian stencil. Those operators are summarized below as stencil masks:

∂• ∂x 0 ∂ • ∂x 1 h 2 ∆• 1 2h -1 1 -1 1 1 2h 1 1 -1 -1 1 2        1 0 1 0 -4 0 1 0 1       
In practice, it is easy to verify that the exact same formula for the Laplacian can be derived by extracting the linear part of the determinant of J T at the cell center, which is nothing else than the relative signed area of the displaced cell: det(J T ) = area(T (C i, j ))/h 2 , thus guaranteeing a consistent discretization.

Hessian determinants

As for the Laplacian, the vector q(ψ ) of determinants of the Hessian of ψ is equivalently obtained either as the quadratic part of det(J T ) or as the signed area of the quadrilateral defined by the gradient vectors of ψ at the four respective cell corners. For each cell of index c, it as the form:

q c = det(J T (x c )) = ψ T c Qψ c , ( 26
)
where Q is a 9 × 9 matrix and ψ c denotes the 9 stacked values of ψ for the cell c and its 8 neighbors.

Boundary Conditions

Since we are working on a rectangular domain, our boundary condition T (δU ) = δU boils down to linear Neumann conditions:

∇ψ (x) T N (x) = 0 , ∀x ∈ ∂U , ( 27
)
where N denotes the normal of the boundary ∂U of the domain U .

In our discrete settings, this is easily accomplished by considering ghost nodes (Figure 6) along the boundaries with enforced mirrored values, e.g., ψ -1, j = ψ 0, j . With such a boundary condition, our Poisson problem (14) needs to satisfy the compatibility condition U u = 1 to be well posed, which is indeed the case.

Implementation details

Our solver requires only the gradient and Laplacian operators from which we can assemble the Laplacian matrix L and efficiently compute the residual vector r at each iteration through the area of the cells displaced by the the gradient of ψ (eq. 11).

The determinants of the Hessians given in equation ( 26) are required only to assemble the full Jacobian J r for Newton or Gauss-Newton iterations.

It is worth noting that the matrix L is symmetric, but negative semidefinite with the constant vector 1 = [1, 1, . . . , 1] as singular vector. In order to enable a fast Cholesky factorization, we fix one arbitrary node to 0 prior to factorizing -L, and project the result in the subspace orthogonal to 1.

MAP INVERSION AND COMPOSITION

Numerical Inversion. Algorithm 1 is designed to efficiently compute the optimal map T u→1 from a non-uniform distribution u towards a uniform one. This map is represented as a set of displaced bi-linear quads. To compute the inverse map T 1→u (y), we need to search the point x such that T u→1 (x) = y. To this end, we first search for the displaced quad containing y using a Bounding Volume Hierarchy (BVH), and then analytically compute its bi-linear barycentric coordinates to recover x. Figure 7 shows the forward map T u→1 computed using Algorithm 1, and the inverse map T 1→u obtained by numerically inverting a regular grid. As expected, after inverting a uniform grid, the resulting density is proportional to the density of the input.

Numerical Composition. This procedure can also be used to compute a transport map T u→v (x) between a pair of non-uniform densities u and v by composing simpler maps as explained in eq. (10). To this end, we first pre-compute the two forward maps T u→1 and T v→1 , and then perform a direct bi-linear interpolation to compute y = T u→1 (x) followed by the above inversion procedure to compute T -1 v→1 (y).

Numerical Interpolation. Interpolating between two densities u and v according to a parameter α is then easily accomplished by:

T α,u→v (x) = (1 -α )x + αT u→v (x) . ( 28
)
This map permits to compute intermediate densities when starting from the density u, but in some cases, it is desirable to compute intermediate densities from a uniform one, which is accomplished as follows:

T α,1→uv = T α,u→v • T -1 u→1 = (1 -α )T -1 u→1 + αT -1 v→1 . ( 29
)
Figure 1 includes a direct use of it to interpolate between two heightfields interpreted as densities. This equation generalizes to n-way barycenters as in Figure 12.

RESULTS

Raw Performance

Our tests were conducted on an Intel i7-7820X CPU at 3.6GHz with 16GB of memory. Our prototype uses Cholmod with nested dissection as fill-in reordering to pre-factorize the Laplacian matrix and solve the linear systems, and Eigen for other matrix and vector operations. In our current implementation, only the pre-factorization step is multi-threaded, while the rest runs on a single core. Table 1 reports the computation cost of each part of our algorithm for different target difficulties and grid resolutions. The prefactorization step of the matrix L depends only on the grid resolution. When working on multiple problems having the same sizes, this step has to be done only once and can thus be neglected. In some applications, the Cholesky factors could even be cached on disk to completely save this step across multiple sessions.

The optimization itself depends on two factors: the cost of one iteration, and the number of required iterations that will be analyzed later. Each iteration mostly involves two sparse triangular solves to compute the search direction candidate d k , a few residual evaluations for the line-search, and one more to estimate β. Since triangular solves are about ×5 times more expensive than computing one residual vector, it makes sense to spend a little more time in the line-search to find a pretty accurate 1D minima, hence reducing the overall number of iterations. In practice, we found that setting a tolerance t α = 10 -2 for the line-search to be a good tradeoff. The cost of one iteration ranges from 10ms to 700ms for grids ranging from about 65k to 4M cells.

The second to last column reports the time to apply the inverse map to a uniform grid having the same resolution as the input grid, including the construction of the BVH. It represents less that 10% of the overall computation time. Moreover, the inverse map can be applied to various grids and samplings without having to recompute the forward map nor to assemble the BVH.

From this table, we see that our prototype implementation is already able to maintain interactive performance for working grids up to about 512 2 cells.

Convergence rate

Beside raw performance, the speed of our algorithm is also conditioned by the number of iterations, which itself depends on the stopping criterion and some notion of difficulty of the target density.

As explain in Algorithm 1, we stop our iteration when a certain level t e of relative accuracy has been reached. Whereas its actual choice likely depends on the application, we found that taking t e = 10 -5 to be a good, rather conservative, starting point as the differences of the maps become hardly visible afterwards (Figure 8). We used this setting for all the results of this paper.

As discussed in Section 4, the rate of convergence of our solver is expected to be correlated with the ratio between the magnitudes of the Laplacians and Hessian determinants. This behavior is depicted Figure 9 on the four disks binary example where we changed the value of the dark area to indirectly adjust this ratio. We can also observe that the working grid resolution has no impact on the convergence rate. Though this experiment might suggest that the convergence rate is somehow correlated with the contrast of the input density, this is not always the case. For instance, if the input density is constant along one axis, then our solver will find the exact solution at the first iteration regardless of the actual values. Moreover, as can be seen on the previous plots, the number of required iterations quickly stalls as the contrast increases.

Table 2 shows convergence rate for increasing grid resolution against a density with known analytic solution. On this example, our solver achieves the same accuracy as the more general and sophisticated solver of [START_REF] Benamou | Numerical Solution of the Optimal Transportation Problem Using the Monge-Ampère Equation[END_REF] while being orders of magnitude faster.

Table 2. Convergence rate and timings (in seconds) for increasing grid resolutions on the problem given in section 6.1 of [START_REF] Benamou | Numerical Solution of the Optimal Transportation Problem Using the Monge-Ampère Equation[END_REF] with L 2 errors computed against the known analysis solution. This table also reports the performance of Benamou et al. solver and the semi-discrete solver of Geogram with the maximal and average Euclidean-distance between their Laguerre centroids and our displaced quad centers (section 7.5). 1).

Comparison with classical nonlinear solvers

Figure 10 compares our original solver (with and without our conjugatedirection acceleration technique) against the damped Newton algorithm described in Section 2.1 combined with our line-search procedure, L-BFGS with 3 recurrence terms, and the Dogleg solver implemented in Ceres.

As can be seen, our solvers completely outperform both L-BFGS and Dogleg in terms of number of iterations, and since our iterations are much cheaper, the comparisons are even more impressive in terms of computation times. For the Julia set example, as expected, the Newton method performs remarkably well in terms of number of iterations, but since each Newton iteration requires the assembly of the Jacobian matrix and to solve a new linear and non-symmetric problem from scratch, each iteration is considerably more expensive and our solvers outperform Newton iterations. Moreover, for more difficult densities as the inverse Julia one, the true Jacobian often leads to poor search directions, as depicted in Figure 5, making the solver nearly stall until it finally find a good direction. Finally, Figure 10 also shows that for difficult densities, our conjugate-direction acceleration technique yields huge improvements, about ×3 for t e = 10 -5 on this example, and more for higher accuracy. Ours, 3.5s [START_REF] Solomon | Convolutional Wasserstein Distances: Efficient Optimal Transportation on Geometric Domains[END_REF]], 1200s

Fig. 12. Bilinear barycenters from the four image corners of the left panel as inputs. Left: using our solver and a generalization of eq. ( 29). Right: convolutional Wasserstein barycenters [START_REF] Solomon | Convolutional Wasserstein Distances: Efficient Optimal Transportation on Geometric Domains[END_REF]] (using the Matlab code provided by the authors). In both cases, transport maps/plans have computed on 299 × 299 images after inversion.

Composition & Interpolation

Figure 11 compares the transport map between pairs of non-uniform densities obtained using our solver and numerical composition, to the L 2 optimal maps. The later has been obtained by modifying our solver to directly solve the full Monge-Ampère equation. If v is non null on the domain, this is easily accomplished by replacing the right-hand-side of equation ( 11) by the integral of u/(v •T ) over each cell when computing the residual vector. This is of course much more expensive than computing an approximately optimal map through composition because 1) the solver has no knowledge on the variations of v, and 2) each iteration now involves costly integral computations. Our composition approximation yields very similar Wasserstein distances but tends to exhibits some undulations.

The same observation can be made in Figure 12 comparing our fast interpolation method to convolutional Wasserstein barycenters [START_REF] Solomon | Convolutional Wasserstein Distances: Efficient Optimal Transportation on Geometric Domains[END_REF]]. This effect is especially visible when the łSž splits in half. On the other hand, entropic regularization yields over blurred results.

Comparison with semi-discrete solvers

Figure 8(c) compares the forward map obtained using our solver with the one obtained by the semi-discrete Newton solver implemented in Geogram [START_REF] Lévy | Geogram[END_REF]]. To match our settings, we generate one weighted Dirac mass per grid center, and define the continuous target as a single quad with a uniform density. This configuration by-passes the very costly integrals over the Laguerre cells which are required in the case of a non-uniform target density v. It can be seen that our displaced quads overlap very well the Laguerre cells of the proved semi-discrete solution. This observation is confirmed in Table 2: both the maximal and average Euclidean-distances between the Laguerre centroids and our displaced quad centers are much smaller than the discretization step.

Such a result was expected since, as shown by [START_REF] Gu | Variational principles for Minkowski type problems, discrete optimal transport, and discrete Monge-Ampere equations[END_REF], semi-discrete solvers can be directly derived from the Monge-Ampère equation ( 8) through a discretization of the input density u as a set of Dirac masses. Furthermore, it is easy to show that the scalar potential ψ optimized by our solver is related to the power-diagram weights w optimized by the semi-discrete solver by: ψ = -2w [START_REF] Lévy | A Numerical Algorithm for L2 Semi-Discrete Optimal Transport in 3D[END_REF]]. Since the discretizations are different, such an equality is only expected to be reached for an infinitely fine grid.

Semi-discrete solvers optimize a convex energy (e.g., eq. ( 2) in [START_REF] Mérigot | A Multiscale Approach to Optimal Transport[END_REF]]) whose gradient is analogue to our residual vector r comming from equation 11, but with the major difference that the target cells are polygonal Laguerre cells instead of simple quads. This difference justifies why the semi-discrete energy is not applicable in our case and yields to Jacobians of the residual exhibiting very different structures and properties. In the semi-discrete case, the Jacobian is symmetric positive semi-definite and similar to a weighted Laplacian [START_REF] Fernando De Goes | Blue Noise Through Optimal Transport[END_REF][START_REF] Gu | Variational principles for Minkowski type problems, discrete optimal transport, and discrete Monge-Ampere equations[END_REF]]. It has a standard Laplacian structure that might change during the optimization. A staggered triangulation has been suggested by de Goes et al. [2014] but edge flips, and thus structure changes, might still be required to reach the true optimal transport solution. In our case of a regular grid, it is nonetheless interesting to observe that edge flips can only occur within each pair of diagonals of each quad, and the edges of the Laguerre cells can only takes 4 directions (0, π /4, π /2, -π /4). In contrast, our quadrilateral cell edges can take arbitrary orientations. In our case, the Jacobian is not symmetric but it has a fixed structure. We split this Jacobian as a fixed Laplacian-like matrix plus the Jacobian of a nonlinear term (which is never computed by our solver). Therefore, our Laplacian matrix should not be confused with the weighted Laplacian that appears in the semi-discrete settings. Even at the first iteration, for which our Laplacian is exactly our Jacobian and the Laguerre cells degenerate to the same uniform grid as us, the semi-discrete Laplacian boils down to a standard Laplacian mask whereas ours is aligned on the diagonal (section 5.1).

Finally, in terms of performance, Table 2 reveals that for comparable settings, our solver is orders of magnitude faster. We further measured that in the Geogram implementation, half of the time is spent in updating the Laguerre diagram, whereas the other half is spent in solving the underlying sparse linear problems. Of course, this high performance gain comes at a loss of generality, in particular semi-discrete solvers can work on arbitrary target domains with arbitrarily distributed Diracs as inputs.

APPLICATIONS

This section aims to evaluate our approach in the context of various computer-graphics applications that make a direct use of our OT solver.

2D Adaptive Sampling

The generation of random samples matching a given density is an essential ingredient in many applications such as Monté-Carlo integration or stippling. This is still an ongoing research area where the main challenge is to find the best tradeoff between speed and sampling quality. The algorithm we implemented is straightforward. Given a target density u, we pre-compute the forward map T u→1 using our solver, build a BVH on the displaced quads, and then apply the inverse map T -1 u→1 to a stream of samples uniformly distributed within the [0, 1] 2 domain. The cost of generating one adapted sample is thus the sum of generating the random sample itself plus the cost of one BVH search that has a O (log(n)) complexity, with n the number of grid cells. For the actual generation, we implemented a simple tiling strategy of a pre-computed uniform pattern of 1024 samples, but for practical Monté-Carlo integration some recent work can stream high-quality uniform samples with very high efficiency [START_REF] Perrier | Sequences with Low-Discrepancy Blue-Noise 2-D Projections[END_REF]]. This procedure is analogue to the classical inverse transform sampling method based on conditional CDFs [START_REF] Pharr | Physically Based Rendering: From Theory to Implementation[END_REF]].

This procedure is illustrated in Figure 13 in the context of stippling, in which case the input density corresponds to the inverse of the input image. Compared to discontinuous conditional CDFs, our L 2 -optimal map is continuous and smooth (for a continuous density) and thus, it much better preserves the intrinsic quality of the input uniform distribution, thus yielding to much lower reconstruction error. Compared to a sophisticated blue noise generator such as the one of [START_REF] Fernando De Goes | Blue Noise Through Optimal Transport[END_REF], our result exhibits some structures in areas of strong anisotropic stretching, but our approach is orders of magnitude faster, it can handle a stream of samples, and by construction, it does converge to the target density. In applications for which the generated samples can be globally processed, we found that very few Lloyd relaxation steps (typically 3 as in fig. 13(b)) are enough to break those structures and reach a very high quality with a very small overhead. On this particular example, the Lloyd steps improve the reconstruction in sparse areas (highlights), but tends to blur the sharpest details.. As detailed in Table 3, our approach can generate 1M samples in half a second on a single CPU core, though we believe that our current BVH implementation has plenty of room for speed improvements.

Inverted mesh Uniform mesh & inversion 3D reconstruction

Fig. 14. Remeshing pipeline. From an input density defined on the parametric domain, we first compute the map T u→1 , then invert a uniform mesh to obtain the desired vertex density. The topology of the uniform mesh is preserved.

Surface remeshing with constrained connectivity

Feature adaptive surface remeshing is a well studied problem, but none of the methods proposed so far permits to constrain the connectivity of the output mesh. This is however required by some applications such as the generation of meshes suitable for hardware tessellation [START_REF] Lambert | Multi-Resolution Meshes for Feature-Aware Hardware Tessellation[END_REF], or to make the mesh connectivity implicit for streaming or compression purposes. To this end, our general approach consists in working in the 2D parametric space [START_REF] Botsch | Polygon Mesh Processing[END_REF], and take advantage that our transport maps are bijective (i.e., foldover free) to deform a given uniform mesh according to some density metric.

Our proposed remeshing pipeline is presented in Figure 14. Starting from a 3D mesh with a rectangular 2D parametrization, we generate a regular grid on which we define a target density u. We then compute the map T u→1 using our solver, and apply its inverse to the vertices of a given uniform mesh. Finally, we recover the 3D positions of the inverted vertices through linear interpolation.

The density u should reflect surface curvatures while taking into account parametric distortions. To this end, we follow the normal lifting technique [START_REF] Nivoliers | Anisotropic and feature sensitive triangular remeshing using normal lifting[END_REF]] that consists in embedding the surface mesh in a 6D space composed of the 3D vertex coordinates p i and scaled normal coordinates γ n i . The parameter γ controls the relative importance given to the normals. The target distribution is then obtained by computing the area of each cell of the working grid in this 6D space. The 3D position and normal coordinates of the cell corners are obtained from the corresponding surface point in parametric space trough linear interpolation. Some results are given in Figure 15 for different values of γ , and a user edited target. Thanks to the high performance of our solver, it is possible to adjust the value γ and paint the target density with instant feedback (see the accompanying video). Finally, once the transport map has been computed we can generate several level of details instantly by simply applying the inverse map to meshes with different resolutions.

Freeform optic design

The problem of caustic design consists in finding a height-field surface producing through refraction of a given light source a desired image on a given receiver (Figure 1). The state-of-the-art approach to this problem [START_REF] Schwartzburg | Highcontrast Computational Caustic Design[END_REF] consists in first computing an OT map between the incoming light source and target densities. This gives an initial solution that is then iteratively improved until the map yields a normal field that is truly integrable. They reported that the initial step requires from 15mn to 4h for 1.5M samples. In this context, computing a strictly L 2 optimal map is again a waste of effort as the resulting normal field wont be integrable anyway. Our solver permits to cut down this computation time to the order of 1s, thus enabling interactive feedback during prototyping (see accompanying video). Results using our OT solver to establish the mapping followed by a direct Poisson integration of the resulting normal field are given in Figure 16.

DISCUSSIONS AND FUTURE WORK

This paper shows that computing optimal transport maps from a continuous probability density to the uniform density is a useful problem for various applications, including some involving pairs of non-uniform densities. By focusing on this simplified OT problem, we devised a dedicated solver based on 1) an original finite-difference discretization of the resulting Monge-Ampère equation, and 2) a new solver for nonlinear simultaneous equations that, for our problem at hand, outperforms classical Newton, quasi-Newton, and Gauss-Newton solvers by orders of magnitude. Even though we designed this new nonlinear solver in our OT context, it would be interesting to investigate how it would perform on other problems involving nonlinear simultaneous equations or even nonlinear least-square problems for which our conjugate-direction acceleration technique can directly be applied.

As with any numerical solvers, our approach comes with its sources of errors. The discretization of the input density on a uniform grid is exact for an input image but requires a lossy conversion for other kind of inputs (e.g., an analytic functions or a piecewise linear triangular mesh). It would thus be very interesting to generalize our approach to be able to work on adaptive grids or even directly on triangular meshes via the use of higher-order FEM discretizations.

Our solver assumes that the image of a quad cell remains a quad which is only an approximation classically admitted when solving OT through the Monge-Ampère equation [START_REF] Benamou | Numerical Solution of the Optimal Transportation Problem Using the Monge-Ampère Equation[END_REF]. This source of error can be reduced by subdividing the grid, or by investigating the use of higher-order elements.

Regarding performance, our experiments revealed that our specialized solver is already orders of magnitude faster that concurrent, but more general, alternatives when applied on similar settings. Our prototype implementation is currently mono-threaded whereas the computation of the residual vectors can trivially be parallelized on either CPU or GPU, and the triangular solves could also benefits from GPU acceleration [START_REF] Li | On Parallel Solution of Sparse Triangular Linear Systems in CUDA[END_REF]]. For difficult densities, it would also be interesting to experiment with a coarse-to-fine multi-resolution strategy, which, if successful, could drastically improve the overall performance of our solver.

In terms of generalization, it would be highly desirable to extend our approach to work on non-uniform grids while retaining the same level of accuracy and speed. This is especially the case for our remeshing application for which the input density is initially defined on a non uniform triangular mesh.

Another interesting question is whether our approach could be extended to 3D. Compared to the 2D setting, a first difference is that the decomposition we made to leverage a Laplacian will led to additional quadratic terms, but those wont be a problem as long as the Laplacian dominate. The most notable difference is that direct sparse solvers, as required to make our approach extremely fast, do not behave as well on 3D grids because of a much larger fill-in. It is thus unclear whether our approach will still outperform Newton iterations by orders of magnitude.

Fig. 3 .

 3 Fig.3. Magnitude of the Laplacian and Hessian determinant of the scalar potential ψ as obtained using our L 2 optimal transport solver from the Julia set image to a constant density.

Fig. 5 .

 5 Fig. 5. Newton step for the second iteration on the inverse Julia problem. Iso-curves of the error e and the local quadratic approximation are shown in gray and brown respectively.

Fig. 7 .

 7 Fig. 7. Visualization of the forward map computed by our solver, and the inverse map obtained after numerical inversion of a uniform grid. Original image courtesy of Fredo Durand.

Fig. 8 .Fig. 9 .

 89 Fig. 8. Closeups of maps obtained for the Julia example with superimpositions of our converged map in orange (24 iterations, e = 4 • 10 -13 ) on top of, from left-to-right: our maps obtained after 1 and 4 iterations, and the Laguerre cells obtained by a reference semi-discrete solvers (together with the Laguerre centroids and the centroids of our quads).

Fig. 10 .

 10 Fig.10. Evolution of the L 2 norm residual with respect to the number of iterations and execution time for the Julia (top row) and inverse Julia (bottom row) targets (respectively, rows 3 and 4 in Table1).

Fig. 11 .

 11 Fig. 11. Comparison of the T u→v maps computed by solving the full Monge Ampère equation (eq. 8) versus the one obtained via composition (eq. 10). Iso-curves in the input domains and their images under the T u→v maps, and the respective transport costs c (T ) are shown for two pairs of distributions. Top row: Gaussian example of Benamou et al. [2014]. Bottom row: a pair of anisotropic BRDF lobes.

Fig. 13 .

 13 Fig. 13. (a) 50k samples generated by numerically inverting a set of tiled uniform blue noise samples with respect to the forward T u→1 map. (b) 3 Lloyd steps applied to the numerical inversion output. (c) Sampling through numerical inversion of conditional CDFs of the same uniform pattern as in (a). (d) 50k blue noise samples generated with [de Goes et al. 2012]. Top-right insets show 256×256 images reconstructed through density estimation from which we have computed RMSE from the original image. Original photograph courtesy of Anton Hooijdonk.

  Fig. 15. Left: Feature sensitive remeshing obtained by setting the normal lifting parameter γ to 1.5. Middle: setting γ to 0 will yield in an isotropic remeshing. Right: remeshing using a user defined input density. 3D model courtesy of cvbtruong from Turbosquid.

Fig. 16 .

 16 Fig. 16. Surface lense generation for the target v and a uniform (left) and non-uniform (right) incoming light sources. Top row shows close ups of the forward map T v →1 and composite map T v →u . In both cases, simulated caustics (bottom row) closely ressemble the initial target v.

Table 1 .

 1 Detailed performance of our OT solver for various input densities.

	Grid Resolution (1/h)	Factorization of -L (s)	Linear solve (step 1) (ms)	Optimization Conj. Dir. (step 2) (ms) Line-search (step 3) (ms) # iterations	Total opt. time (s)	T -1 u→1 map (s)	Total Time (s)
	256	0.16	5.9	0.6	1.1 26 0.20 0.043 0.43
	512	0.64 25.7 2.9	5.2 25 0.84	0.18	1.66
	1024 2.66	115 13.5 27.1 25 3.90	0.78	7.34
	2048 11.2	518 68.9 105 25 17.3	3.60	32.1
	256	0.16	6.1 0.58 0.97 15 0.11	0.04	0.31
	450	0.50	20	1.8	6	4 0.11	0.12	0.73
	450	0.50	20	2.2	4	25 0.64	0.14	1.38
	512	0.65 25.9 2.6	5.7	5 0.17	0.16	0.98

Table 3 .

 3 Processing time for generating density-adapted samples.

	50k	250k 500k	1M

ACM Transactions on Graphics, Vol. 37, No. 6, Article 249. Publication date: November 2018.

iteration, e =

2.3e -3 4 iterations, e = 8e -6Laguerre cells

ACKNOWLEDGMENTS

This work has been supported by the ANR RichShape project (ANR-14-CE24-0004). We would like to thank the authors of [Solomon et al. 2015] and[de Goes et al. 2012] for making their code publicly available, Bruno Lévy for pointing us to his public implementation of a semi-discrete solver in Geogram, and the anonymous reviewers for their insightful comments. The 3D models and images are provided cour Additionally, we would like to acknowledge the support from the SUTD Digital Manufacturing and Design (DManD) Centre, supported by the Singapore National Research Foundation.


George Pechlivanoglou
email: george@pechlivanoglou.com

C N Nayeri

C O Paschereit2

S Y Mpo

O Tatin

Backflow flaps as stall control elements on wind turbines

Keywords:

Wind turbine blades suffer from static and dynamic stall effects during their normal operation. At the same time strong secondary flows (cross-flows) exist at the inner and the outermost regions of the blade. Most of the modern wind turbine blades are equipped with flow control devices such as vortex generators (VGs) in an effort to avoid stall and stabilize their aerodynamic performance. The current invention and research work proposes a solution to the unsteady stall issues and at the same time offers a performance boosting potential. The solution comprises a passively actuated backflow flap positioned in span-wise or chord-wise orientation depending on the local flow characteristics. Wind tunnel tests in combination with extensive field preliminary work as a part of the VG development efforts of the authors form the ground work for a new and innovative aerodynamic solution for wind turbines.

INTRODUCTION

Wind turbine blades traditionally were characterized by low drag "clean" configurations aiming at peak efficiency. However in recent years it has become apparent that the exposed rotorblade surface suffers from heavy rain erosion [START_REF] Hugh Keegan | On erosion issues associated with the leading edge of wind turbine blades[END_REF] as well as surface roughness and irregularities [START_REF] Pechlivanoglou | The effect of distributed roughness on the power performance of wind turbines[END_REF]. Many blade manufacturers therefore tend to shift their focus from the ultimate efficiency towards performance reliability. Furthermore there is a rapid increase in rotor diameters in an effort to increase the capacity factors of wind turbines and reduce the Levelized Cost of Energy while trying to maintain constant mechanical loads. This leads to slender blade plan-form designs which utilize more cambered airfoils in order to produce the same or even more lift from a given blade area.

In order to avoid stall and its adverse effects in AEP performance as well as with respect to aeroelastic loads many manufacturers try various designs of vortex generators [START_REF] Castro | The use of a new fatigue tool (albdes) to analyse the effects of vortex generators on wind turbines[END_REF]. These energize the boundary layer and avoid stall. The authors have done considerable research and development efforts in this field and more information on VGs and their function and placement can be found at the publications of Pechlvanoglou et al. [START_REF] Mueller-Vahl | Vortex generators for wind turbine blades: A combined wind tunnel and wind turbine parametric study[END_REF], [START_REF] Pechlivanoglou | Passive and active flow control solutions for wind turbine blades[END_REF]. The main disadvantage of vortex generators is that they generate drag. This is especially present in the situations where the VG effect is not needed. There they act as efficiency inhibitors rather than performance boosters. The aim of the current research work is to propose and analyze an alternative to VGs for the blade regions where VGs could be replaced.

CONCEPT AND INITIAL RESEARCH

The basic idea for the current research was first expressed by Prof. Bechert et al. from HFI TU Berlin in the 90's [START_REF] Dw Bechert | Aeroflexible oberflächenklappen als rückstrombremsen nach dem vorbild der deckfedern des vogelflügels[END_REF]. After examining the behavior of birds when gliding at high angle of attack (e.g. when landing), it was apparent that many species were able to actually land without flapping their wings. When trying to explore the way those birds avoid stall, it became apparent that they were equipped with a stall control solution on their wings (see Fig. 1). A line of feathers located close to the trailing edges of the bird wings is able to lift during the initial development of stall and thus restrict the stall severity and enhance the generation of lift. This mechanism allows certain bird species to achieve controllable landings with very low speeds without the need of energy consuming wing flapping actions.

Further research on this topic was performed by Meyer [START_REF] Kurt | Experimentelle untersuchungen von rückstromklappen auf tragflügeln zur beeinflussung von strömungsablösungen[END_REF] and Schatz [START_REF] Schatz | Numerische Simulation der Beeinflussung instationärer Strömungsablösung durch frei bewegliche Rückstromklappen auf Tragflügeln[END_REF] at the institute of fluid dynamics and technical acoustics (ISTA/HFI) of TU Berlin. The work of Meyer focused around the HQ17 airfoil and glider aircraft applications. Extensive wind tunnel tests demonstrated the feasibility of this concept and it was shown that a slotted version of the design, similar to bird feathers was the configuration that achieved the highest lift increase with a subsequently minimal drag penalty (see Fig. 2 and3).

The proposal of Meyer at al. was the utilization of such devices on aircraft as a means of lift enhancement during landing. The solution was tested on a glider and performed very well. However to the best of the authors' knowledge it was never implemented on a serial basis for any other aircraft.

Stall control in the field of wind energy

Wind turbine blades have to operate in severe atmospheric conditions under extreme turbulence, gusts and wind sheer. Rotorblade manufacturing is also heavily driven by the extreme cost-reduction pressure of the competitive electricity market. Automated high precision production is therefore out of the scope for wind turbine blade manufacturers. Therefore wind turbine rotorblades have to be able to operate effectively without stall under extreme conditions while being simple and very tolerant to manufacturing deviations. A common method in the direction of "reliable aerodynamics" for wind turbines is the use of vortex generators.

The company SMART BLADE GmbH TM in collaboration with the authors has developed and installed retrofit VGs (Fig. 4) on more than 2000 wind turbines as a means of performance boosting. The performance improvement is achieved by the elimination of the flow separation from the inner 1/3 of the blade which usually operates under static stall. Depending on the blade design details as well as the local site characteristics, the retrofit VGs achieve a performance enhancement of 1-3% over the Annual Energy Production (AEP) of the turbine [START_REF] Mueller-Vahl | Vortex generators for wind turbine blades: A combined wind tunnel and wind turbine parametric study[END_REF], [START_REF] Ekaterinaris | Know-blade task-2 report aerodynamic accessories[END_REF] and [START_REF] Pechlivanoglou | Vortex generators for wind turbine blades: Wind tunnel tests, field simulations and structural analysis[END_REF].

A drawback however of the VG solution is the fact that it consists of rigid units that are able to enhance performance only in certain operating regimes. Outside their operational envelope, vortex generators simply add drag. The net benefit from the use of VGs on blade regions that experience constant stall is positive, however other stall control solutions (e.g. backflow flaps [START_REF] Pechlivanoglou | Rotorblatt[END_REF]) are explored in an effort to further improve the performance of outboard blade sections, which do not experience severe permanent stall but are currently equipped with VGs (see Fig. 5).

WIND TUNNEL TESTING

Wind Tunnel Setup Wind tunnel measurements were carried out in the wind tunnel facility of ISTA/HFI, TU Berlin. The closed loop wind tunnel is powered by a 500 kW speed regulated motor. The test section has an overall length of 10m, a cross-sectional area of 2 × 1.4 m 2 and a contraction ratio of 6.25 : 1 (see Fig. 6). The turbulence level of the free stream is below 0.5% in the test section. The airfoil was a NACA 63 (3) 618 with a chord length of 600 mm. The airfoil is mounted to an external 6-component balance, which is decoupled from the wind tunnel. The test section is equipped with splitter walls to minimize the influence of boundary layer effects while the wing is mounted horizontally extending over the entire span (a=1.554m) between the splitter walls. The gaps (≈ 0.5 mm) between test wing and splitter walls are covered with thin adhesive tape to reduce the pressure losses at the tips of the airfoil. The sealing tape is not in contact with the splitter walls and thus the force measurements are not impaired. All tests are carried out at a wind speed of 29 m s . The resulting Reynolds number based on the chord length was Re ≈ 1.1 • 10 6 .

Generally all the wind tunnel measurements performed at the facilities of HFI TU Berlin undergo a solid blockage and wake blockage correction according to Barlow et al. [START_REF] Pope | Low Speed Wind Tunnel Testing[END_REF] and Meyer [START_REF] Dw Bechert | Aeroflexible oberflächenklappen als rückstrombremsen nach dem vorbild der deckfedern des vogelflügels[END_REF]. However in the case of the backflow flap tests on the NACA 63 (3) 618 test wing some of the solidity requirements of the blockage correction were causing accuracy problems. Considering that the correction does not give reliable results for solidity levels higher than 7.5% [START_REF] Pope | Low Speed Wind Tunnel Testing[END_REF], which is the case for our test wing at the AoA range where backflow flaps are operational, it was decided to avoid corrections all-together. For an actual blade design process, corrections would be necessary, but for our comparative study, a total absence of polar corrections is considered non-problematic.

Vortex Generators Extensive multi-parameter testing of vortex generators has been performed by various researchers (e.g. [START_REF] Wa Timmer | Summary of the delft university wind turbine dedicated airfoils[END_REF], [START_REF] Mueller-Vahl | Vortex generators for wind turbine blades: A combined wind tunnel and wind turbine parametric study[END_REF] and [START_REF] Fuglsang | Development of the riso wind turbine airfoils[END_REF]), however the specific shape and positioning of the SMART BLADE vortex generators was tested on the NACA 63 (3) -618 in order to examine the combination of lift enhancement at high AoA in combination with the re-

) b y t h e a u t h o r s o n t h e N A C A 6 3 (3 ) 6 1 8 t e s t w i n g a t t h e T U B e r l i n w i n d t u n n e l . I n p r i n c i p l e t h e b a c k fl o w fl a p s w e r e p o s i t i o n e d a t t h e t r a i l i n g e d g e o f t h e w i n g . W h e n m u l t i p l e b a c k fl o w fl a p c o n fi g u r a t i o n s w e r e t e s t e d t h o s e w e r e p l a c e d u ps t r e a m o f t h e t r a i l i n g e d g e a n d u p t o t h e p o s i t i o n o f m a x i m u m t h i c k n e s s o f t h e t e s t w i n g . F r o m a l l t h e d i ffe r e n t b a c k fl o w fl a p m a t e r i a l s , s i z e s , p os i t i o n a n d o t h e r c o m b i n a t i o n s t h e c u r r e n t p a p e r p r e s e n t s t h e p e r f o r m a n c e r e s u l t s o f t h e f o l l o w i n g c o n fi g u r a t i o n s :

1 . B a s e l i n e a i r f o i l e x h i b i t i n g f r e e t r a n s i t i o n (i . e . n o t t r i p p e d ). What was especially interesting in the current case was the performance with double rows of backflow flaps, where the lift was further enhanced and the flow separation was delayed even more (see Fig. 11). The outcome of that case was an almost neutral behavior with respect to drag at low AoA in combination with a drag reduction at higher AoA due to the flow separation delay and the subsequent reduction of the pressure drag component. On the negative side, the pre-stall lift seems to be slightly adversely affected by the presence of the backflow flap. The reason for that seems to be the slight modification of the Kutta condition of the airfoil at the pre-stall region. More specifically, the existence of the compliant backflow flap material at the trailing edge region of the suction sice of the airfoil leads to a less-than-optimal trailinge edge shape. The overall behavior therefore of this solution was similar with that of the VGs, but without the associated drag penalty.

. S i n g l e b a c k fl o w fl a p w i t h a w i d t h

T h e r e s u l t s o f t h e w i n d t u n n e l m e a s u r e m e n t s w e r e i n g e ne r a l a g r e e m e n t w i t h t h e r e s u l t s o f M e y e r [ 7 ] a n d s h o w e d t h e e x p e c t e d b e h a v i o r o f s t a l l d e l a y a n d e n h a n c e d l i f t a t h i g h A o A . A t t h i s p o i n t i t i s w o r t h n o t i n g t h a t t h e w i n d t u n n e l t e s t p r e s e n t e d i n t h i s s t u d y a s w e l l a s t h o s e o f M e y e r [ 7 ] w e r e p e r f o r m e d a t t h e s a m e w i n d t u n n e l f a c i l i t y w i t h t h e s a m e i n s t r u m e n t a t i o n a n d g e n e r a l b o u n d a r y c o n d i t i o n s . T h e d i

The addition of a Gurney flap (see Fig. 9) at the pressure side enhanced even more the lift through the entire range of inflow angles with a very small overall drag penalty, thus proving that the existence of a Gurney flap with its associated local circulation modification does not have an adverse effect on the backflow flap physical mechanism. This makes this combined solution quite attractive for wind turbine applications.

SMARTVIZ TM FIELD ANALYSIS

Description of the method The final part of the current research was the investigation of the aerodynamic performance for a real life wind turbine blade in the field. This was possible with the use of the proprietary experimental flow visualization technique, called SMARTviz jointly developed by SMART BLADE GmbH and HFI/ISTA TU Berlin. This method is based on the well known use of flow tufts to visualize the near-wall flow. However the use of modern high resolution digital cameras and image processing capabilities allows the automatic position recognition of each tuft and the calculation of the statistics of its behavior in time. In addition to that, parameters such as the instantaneous wind speed and wind direction at the turbine are recorded and synchronized with the images thus forming a comprehensive quantitative aerodynamic performance database for the blade. More details on the SMARTviz method can be found in the relevant publications of Vey et al. [START_REF] Vey | Extracting quantitative data from tuft flow visualizations on utility scale wind turbines[END_REF], [START_REF] Vey | Experimental and numerical investigations of a small research wind turbine[END_REF] and [START_REF] Vey | Utility scale wind turbine yaw from a flow visualization view[END_REF].

The wind turbine selected for this experimental campaign was a 1.5MW pitch regulated, variable speed utility scale turbine. The turbine was located at the state of Lower Saxony in Germany and was equipped with 37m rotorblades which were not equipped with any form of flow control add ons (see Fig. 12). The condition of the turbine and its rotorblades was good without signs of leading edge erosion. More than 6000 flow tufts were installed along the entire span of the blade and the turbine was then put back in operation. Several thousands of images were captured until enough images were collected for every wind-bin of the operating range of the turbine. At this point it is important to note that only limited images were captured for wind speeds above the nominal wind speed. The reason for that is the fact that above the rated wind speed the pitch system rotates the blades towards reduced inflow angles thus eliminating stall from the outboard blade regions (i.e. this region is irrelevant for backflow flaps and/or VGs).

Field Results

The experimental field results are analyzed by means of various statistical modes. However one of the prevailing methods is that of the standard deviation of the flow tuft position [START_REF] Haans | Stall in yawed flow conditions: A correlation of blade element momentum predictions with experiments[END_REF] (i.e. tuft activity). On the color scale (see Fig. 13), red regions denote stall (i.e. high tuft activity) and blue regions denote attached flow (i.e. low tuft activity). The right side graph is a form of synthetic flow streamlines on the blade surface. It is the result of the LIC (Linear In- tegral Convolution) [START_REF] Cabral | Imaging vector fields using line integral convolution[END_REF] method and it is very useful for the identification of the cross-flow zones on the blade. After the analysis of the flow tuft behavior, the date showed some distinctive patterns. The inner 1/3 of the blade span exhibited the typical static stall behavior for the entire range of wind speeds below the nominal. This is the region that is typically treated with vortex generators in order to re-energize the flow and reduce stall.

The outboard part of the blade exhibited a small region of static stall at the trailing edge. Due to the low frequency of the image capture (2 images per rotation) it was not possible to examine the presence of dynamic stall, especially in the event of gust or temporary yaw misalignment. However the onset of static stall at the outboard part of the blade anyway denotes that this region of the blade could be benefited by the installation of backflow flaps. In essence the stall delay effect of backflow flaps would act beneficially, preventing the loss of lift (and therefore the loss of AEP) due to T.E. stall. Furthermore the stabilization of lift would reduce the loading variations throughout the rotation of the blade and would therefore be beneficial with respect to the fatigue loading of the blade. The effects of such devices would therefore be overall beneficial, probably justifying the added device cost and installation & maintenace cost.

DISCUSSION & CONCLUSIONS

The current paper describes the efforts of the authors to identify experimentally the real-life aerodynamic shortcomings of rotorblades and to develop methods to rectify the blade performance. Wind tunnel tests have been performed on a relatively common wind turbine airfoil (NACA 63 (3) 618) with vortex generators as well as with bio-inspired backflow flaps. The results of the wind tunnel tests are very encouraging for the use of backflow flaps on wind turbines either as retrofit solutions or integral design elements. It was shown that the performance of backflow flaps is similar to that of VGs with respect to lift enhancement, but without the inevitable drag penalty of VGs for out-of-design operation points.

The aerodynamic field analysis also demonstrated the potential use of such solutions as retrofit elements for the reduction of trailing edge stall of existing rotorblades. Future work in this project will include field experiments by means of SMARTviz of the turbine equipped with VGs and back-flow flaps. The realization of the backflow flaps is in process together with a large materials corporation and will have the form of custom-built multi-layered tape. Special fiber reinforcements are integrated in the tape material in order to enhance the fatigue life of the tape and thus reduce maintenance and replacement costs.

Furthermore, additional efforts for full characterization of this solution are necessary in order to be able to simulate the turbine performance by means of aeroelastic tools. Finally the last part of the planned future experiments include the use of backflow flaps in the dynamic pitching setup of the laminar wind tunnel of TU Berlin in order to study the dynamic response of backflow flaps.

Figure 1 .

 1 Figure 1. Backflow feathers lifted by the trailing edge separation. [Source: Stewart Baird]

Figure 2 .Figure 3 .

 23 Figure 2. Segmented backflow flaps mounted near the trailing edge of the HQ-17 airfoil. [Source: Meyer[7]]

Figure 4 .

 4 Figure 4. 3D model of the SMART BLADE vortex generator unit. [Source: SMART BLADE GmbH]

Figure 5 .

 5 Figure 5. View of an outboard section of a SIEMENS blade with trailing edge serrations and VGs. [Source: SIEMENS]

Figure 6 .Figure 7 .Figure 8 .

 678 Figure 6. Schematic representation of the airfoil test section of the ISTA/HFI TU BERLIN large wind tunnel

  o f 1 0 % c m o u n t e d a t 9 0 % c v i a fl e x i b l e h i n g e a n d w i t h a d e fl e c t i o n l i m i t o f 6 0 • . T h e d e fl e c t i o n w a s l i m i t e d b y m e a n s o f Dynema T M s t r i n g s c o n n e c t i n g t h e b a c k fl o w fl a p a n d t h e t r a i l i n g e d g e o f t h e w i n g w i t h a d i s t a n c e o f 1 /1 0 o f t h e w i n g ' s s p a n .

Figure 9 .Figure 10 .

 910 Figure 9. S c h e m a t i c v i e w o f G u r n e y fl a p s [ 1 6 ]

  ffe re n c e s i n t h e p e r f o r m a n c e o f t h e H Q -1 7 a i r f o i l o f M e y e r a n d the NACA airfoil presented in this study has most to do with the inherent airfoil characteristics rather than the different implementation of the backflow flap solution.

  ow Flaps B ackfl . F. + G urney F. P orous B ackfl . F.

Figure 11 .

 11 Figure 11. Experimental results of various backflow flap configurations mounted on the NACA 63 (3) 618 test wing, measured at TU Berlin.

Figure 12 .

 12 Figure 12. Picture of the turbine with the flow tufts installed on one of its blades.

Figure 13 .

 13 Figure 13. Experimental SMARTviz field results showing on the left side the standard deviation[20] (i.e. activity) of tufts. Red regions denote stall (i.e. high tuft fluctuations) and blue regions denote attached flow. The right side graph is the result of LIC (Linear Integral Convolution)[21] method showing the streamlines and respective cross-flow zones.

ACKNOWLEDGMENTS

The authors would like to acknowledge the contribution of Dr.-Ing. S. Vey of HFI TU Berlin and of Dipl.-Ing. G. Weinzierl of SMART BLADE for the entire SMARTviz field project.

Interactive Techniques, SIGGRAPH '93, pages 263-270, New York, NY, USA, 1993. ACM.