James Gregson

Ivo Ihrke

Nils Thuerey

Wolfgang Heidrich

From Capture to Simulation -Connecting Forward and Inverse Problems in Fluids

Keywords: I.3.3 [Computer Graphics]: Picture/Image Generation-Digitizing and Scanning fluid-capture, velocity estimation, optimization

Velocity reconstruction & advection

: Low-resolution captures obtained by tomographic scanning (left) are used as inputs to our method which estimates physically plausible dense velocity fields. Such velocity fields fully determine the fluid state and can be applied in a variety of applications including fluid super-resolution (right) allowing capture to be integrated into pipelines for visual effects simulation.

Introduction

Fluid simulation and fluid capture are active research areas in computer graphics. Due to difficulties in acquiring full flow descriptions, the use of captured data in simulations has been limited [START_REF] Wang | Physically Guided Liquid Surface Modeling from Videos[END_REF][START_REF] Li | Water Surface Modeling from a Single Viewpoint Video[END_REF].

We aim at decreasing this gap by showing a connection between proximal methods and fluid tracking as well as simulation. Proximal methods are a class of recent, successful numerical methods from the image processing/computer vision community that split difficult multi-objective optimization schemes into simpler subproblems. These sub-problems are expressed and accessed only by application-specific solvers known as proximal operators. This approach encourages modularity and the development of highly optimized solvers for each sub-problem.

Our key finding is that the pressure solve employed in many fluid simulators can be interpreted as a proximal operator. This enables its incorporation as a physical constraint into proximal methods. We demonstrate the flexibility of this approach by applying it to the problem of fluid tracking where it allows for an elegant inclusion of prior knowledge into the velocity estimation procedure. In particular, we couple the Navier-Stokes equations into the estimation procedure, which enables significant improvements in tracking accuracy as compared to state-of-the-art techniques.

Having an accurate velocity field of the flow phenomena enables flexible new applications, such as applying fluid super-resolution techniques to the captured flows, as well as re-simulations in modified domains. In this way we demonstrate a meaningful integration of captured data into fluid simulation pipelines. The specific contributions of this paper are

• establishing the relationship between fluid pressure solves and proximal operators,

• a physically constrained, multi-scale tracking method for fluid densities to reconstruct temporally coherent velocity fields,

• re-simulation of the captured fluids for adding detail or changing the domain, and

• initial results for solving other mixed forward and inverse problems using the same proximal operator framework.

We believe many more such applications can be developed based on the same methodology. As such, we are confident that recent convex optimization methods will prove to be a very powerful tool in the field of computational fluid dynamics.

Related Work

Fluid Capture aims at the measurement of the full fluid properties of a flow existing in the real world. The flow is typically characterized by its physical properties such as advected densities and fluid velocities. Capture work is made difficult by the fact that optically observable properties are further influenced by other factors than those governing the flow. These phenomena include e.g. emission, scattering, absorption, fluorescence, and refraction.

For this reason, capture technologies in computer graphics and computer vision have first concentrated on phenomenological descriptions. [START_REF] Hasinoff | Photo-consistent Reconstruction of Semitransparent Scenes by Density-sheet Decomposition[END_REF], as well as [START_REF] Ihrke | Image-Based Tomographic Reconstruction of Flames[END_REF] reconstruct emission volumes of flames. The captured information relates to the chemical composition, density, and burn rate of the combustible but not to physical flow properties.

Smoke scanning [START_REF] Hawkins | Acquisition of Time-Varying Participating Media[END_REF][START_REF] Gu | Compressive Structured Light for Recovering Inhomogeneous Participating Media[END_REF]] recovers scattering properties or smoke densities, but again, no information relevant to flow. [START_REF] Atcheson | Time-resolved 3D Capture of Non-stationary Gas Flows[END_REF] recover dynamic refractive index volumes, a quantity related to fluid density, while [START_REF] Gregson | Stochastic Tomography and its Applications in 3D Imaging of Mixing Fluids[END_REF] investigate fluid mixing by observing dye concentrations. The transport of both observed quantities is directly governed by the underlying fluid velocities, which are therefore suitable data for velocity tracking. However, all references above only recover per-frame information without explicit temporal correlation. Fluid velocities are not computed.

Methods described in the physics literature often involve sophisticated hardware setups such as scanning with the help of laser and electron beams, or complex optical arrangements. Examples include Particle Image Velocimetry (PIV), (Planar) Laser Induced Fluorescent Imaging (PLIF), Electron Beam Fluorescence (EBF), or Schlieren-based systems. In some cases, such as PIV, the flow has to be physically seeded with particles to enable its visualization. An extensive overview is given by [START_REF] Tropea | Handbook of Experimental Fluid Mechanics[END_REF]. For Schlieren systems, both classical and computational, the reader is referred to [START_REF] Settles | Schlieren and Shadowgraph Techniques[END_REF][START_REF] Panigrahi | Schlieren and Shadowgraph Methods in Heat and Mass Transfer[END_REF]. A relatively recent addition to the fluid capture toolbox is TomoPIV [START_REF] Elsinga | Tomographic Particle Image Velocimetry[END_REF]], a method to reconstruct instantaneous 3D PIV data. None of the capture methods mentioned above, except for EBF, is able to capture fluid velocity directly.

Velocity Estimation for 3D flows has only received sparse attention due to its inherent difficulties. In graphics, 3D velocity reconstruction have been used to guide fluid simulations of water surfaces. [START_REF] Wang | Physically Guided Liquid Surface Modeling from Videos[END_REF] capture free-flowing water columns, whereas [START_REF] Li | Water Surface Modeling from a Single Viewpoint Video[END_REF] aim at generating dynamic height field water surfaces from video.

The fluid mechanics literature has more extensive coverage, even though full 3D tracking is rare and temporal coherence has not been addressed. Early work aimed at compensating PIV plane scanning artifacts by performing an operation similar to rolling shutter compensation on acquired volume data sets [START_REF] Van Vliet | Time-Resolved, 3D, Laser-Induced Fluorescence Measurements of Fine-Structure Passive Scalar Mixing in a Tubular Reactor[END_REF]]. However, only a coarse mean flow was computed by crosscorrelation.

The equivalence between optical flow and fluid flow was realized early, but its formal conditions have not been established until recently [START_REF] Liu | Fluid Flow and Optical Flow[END_REF]. The authors put particular emphasis on the image space projection of 3D transport mechanisms and associated boundary effects. The full 3D transport equation was inverted by [START_REF] Su | Scalar Imaging Velocimetry Measurements of the Velocity Gradient Tensor Field in Turbulent Flows. I. Assessment of Errors[END_REF] in the context of scanned PLIF data.

The acquisition of fluid velocities using optical flow techniques rests on physical priors, such as divergence-free constraints for incompressible flows [START_REF] Yuan | Discrete Orthogonal Decomposition and Variational Fluid Flow Estimation[END_REF]Vlasenko and Schnörr 2009]. These methods typically require complex discretization steps due to their higher order regularizers. Moreover, the discretized operators must be compatible, i.e. they must preserve the continuous case vector identities in the discrete setting. As an alternative to these complex numerical techniques, basis function expansions that only cover the appropriate sub-spaces have been used, e.g. divergencefree wavelets [START_REF] Kadri Harouna | Divergence-free Wavelets and High Order Regularization[END_REF]]. This method allows for an elegant incorporation of the constraint and a natural multi-scale representation, but inter-connecting with other parts of a fluid simulation pipeline requires conversion. All examples discussed in this paragraph have been developed for (quasi-)2D flows.

In contrast to the fluid mechanics literature, we show that pressure projection can be naturally incorporated into a velocity estimation scheme as a convex regularizer. This allows for a simpler handling of physical constraints in the estimation, while maintaining compatibility and modularity of the implementation. We tightly couple the velocity estimation and fluid simulation steps to predict good initial guesses for the flow velocities in subsequent time steps, effectively coupling a fluid simulation with global history into the estimation procedure. Most similar to this approach, [START_REF] Barbu | Joint Estimation of Volume and Velocity in TomoPIV[END_REF] have recently coupled two frame estimates in a discrete TomoPIV setting, however, without incorporating physical constraints as we do.

Fluid simulation has the goal of calculating the motion of a flow based on a set of boundary conditions and a physical model. Such simulations have a long history in research. One of the influential papers in the area, e.g., introduced the staggered grid discretization [START_REF] Harlow | Numerical Calculation of Time-dependent Viscous Incompressible Flow of Fluid with Free Surface[END_REF], which is the basis for many works in computer graphics. While Foster and Metaxas [1996] were the first to use fluid simulations in graphics, the introduction of the stable fluids approach [START_REF] Stam | Stable Fluids[END_REF]] has lead to the wide-spread use of fluid simulation for visual effects. These simulations have since become a crucial tool in digital movie productions, especially the FLIP method [START_REF] Zhu | Animating Sand as a Fluid[END_REF]] and the purely particle based SPH approaches [Müller et al. 2003] are widely used. Recent works make it possible to simulate advanced rigid-body interactions with particle-based flows [START_REF] Akinci | Versatile Rigid-fluid Coupling for Incompressible SPH[END_REF], propose adaptive discretizations for gains in efficiency [START_REF] Ando | Highly Adaptive Liquid Simulations on Tetrahedral Meshes[END_REF], and show how to handle complex incompressible materials [START_REF] Stomakhin | A Material Point Method for Snow Simulation[END_REF]]. The single phase fluid solver we use is in line with those commonly used in graphics. A good overview is, e.g., given in the book by R. [START_REF] Bridson | Fluid Simulation for Computer Graphics[END_REF]. Our particular choice of algorithms, such as the use of BFECC [START_REF] Kim | Flow-Fixer: Using BFECC for Fluid Simulation[END_REF]] is described in more detail below.

In a flow modification context, previous work has demonstrated the generation of flow detail based on physical models [START_REF] Kim | Wavelet Turbulence for Fluid Simulation[END_REF][START_REF] Narain | Fast Animation of Turbulence Using Energy Transport and Procedural Synthesis[END_REF]], which we show to be applicable to captured data using our estimation framework.

Our method can be considered as a simulation that is guided by observed data. Previous work has employed non-linear optimizations with the adjoint method to guide fluids [START_REF] Mcnamara | Fluid Control Using the Adjoint Method[END_REF], while other works have used forces-based approaches targeting liquid surfaces [START_REF] Shi | Taming Liquids for Rapidly Changing Targets[END_REF][START_REF] Thuerey | Detail-Preserving Fluid Control[END_REF]]. The guiding of high resolution simulations based on a coarse input flow has, e.g., been investigated by [START_REF] Nielsen | Guiding of Smoke Animations through Variational Coupling of Simulations at Different Resolutions[END_REF], and we will explore the connection between guiding and our ADMM framework in Sec-tion 6.

Background and Notation

Incompressible Fluid Flow can be modeled by the incompressible Navier-Stokes equations. These equations describe the timeevolution of the fluid velocity u = (u, v, w) as a function of itself, pressure P , density ρ, viscous forces induced by the viscous stress tensor τ as well as external body forces f , subject to the condition that the flow be mass-conserving (divergence-free):

∇ • u = 0 (1) ut + u • ∇ u + ∇P/ρ = (∇ • τ )/ρ + f (2)
Throughout this paper we will use symbols with arrows (e.g. u) to denote continuous vector fields, and bold symbols (e.g. u) to denote the corresponding discretized field, with all components rearranged into a single column vector.

In a general flow, the advection of a scalar field φ based on a velocity field u is given by

φt + ∇ • ( uφ) -α∇ 2 φ = s,
where φ refers to an observable concentration of a marker fluid such as a dye or smoke, α is the diffusivity and s is a source term. The marker fluid and the ambient fluid are assumed to have approximately the same mass density. In the case of s = 0, negligible diffusivity (α = 0), and incompressibility (∇ • u = 0), the equation reduces to

φt + u • ∇φ = 0, (3)
which corresponds to a pure advection. We note that Eq. 3 is equivalent to the brightness constancy constraint of optical flow as introduced by [START_REF] Horn | Determining Optical Flow[END_REF].

Eqs. 2 and 3, in conjunction with the constraint in Eq. 1, describe the time-evolution of a non-diffusive incompressible fluid flow carrying a marker density. A difficulty in fluid simulation is that there is no time-evolution equation for fluid pressure P . When evolving the fluid velocities, the gradient of the scalar pressure field can only contribute to divergent portions of the updated velocity field. Modern incompressible fluid simulations advance the velocity forward in time, omitting the ∇P term. This is followed by a pressure-projection step in which a pressure field is computed that exactly counteracts divergent components of velocity introduced during time stepping.

This pressure projection step is the central component for the new proximal operator that we introduce below.

Convex Optimization is an area that has seen rapid progress in recent years, and has found many applications in areas such as imaging and inverse methods. A common form of convex optimization problem is given as

min u F (u) + G(u), (4)
where F is often an ill-posed linear least squares data term F (u) = b -Au 2 2 , while G is a convex regularizer or prior. To solve this kind of problem, we can turn the unconstrained optimization problem of Equation 4 into a constrained one by introducing a slack variable v:

min u,v F (u) + G(v) subject to u = v.
(5)

This splitting approach is advantageous because it separates the two component terms F and G in Eq. 4, and allows us to use an array of recent optimization methods to solve the problem. Many recent non-linear solvers are proximal methods [START_REF] Parikh | Proximal Algorithms[END_REF], which are based on the so-called proximal operators for F and G:

prox λF (z) := arg min u F (u) + 1 2λ z -u 2 2 (6) and prox λG (z) := arg min v G(v) + 1 2λ z -v 2 2 . (7)
These two operators are individual minimizers for F and G respectively, which remain close to some reference point z. In this way the proximal operator can be interpreted as a trust-region method with weight 1 λ and has several benefits, among them • the proximal operators for F and G can be independently derived and implemented, making it possible to easily combine high performance implementations of different F and G.

• prox λF is strongly convex even when F is weakly convex (e.g. under-determined for least-squares problems)

• generally improves conditioning of quadratic F by adding a small multiple ( 1 2λ ) of identity. Given the proximal operators for a given problem, we can choose between many different algorithms to solve Equation 5. One of the most popular and easy to implement choices is the Alternating Direction Method of Multipliers (ADMM, [START_REF] Parikh | Proximal Algorithms[END_REF]), which can be expressed with proximal operators as follows:

Algorithm 1 Alternating Direction Method of Multipliers 1: procedure ADMM(prox λF , prox λG ) 2:
while k < maxIters do 3: // Update primary variable u 4:

u k+1 ← prox λF (v k -q k ) 5: 6: // Update slack variable v 7: v k+1 ← prox λG (u k+1 + q k ) 8: 9:
// Update Lagrange multipliers 10:

q k+1 ← q k + u k+1 -v k+1 11:
end while 12: return u k 13: end procedure Here k is an iteration index, v is a vector of dual (slack) variables, and q is a vector of Lagrange multipliers that is initialized to zero. For convex F and G, and reasonable choices of λ, ADMM causes the value of u k to converge to a minimizer of the original optimization problem of Equation 4 [START_REF] Parikh | Proximal Algorithms[END_REF].

Pressure Projection as Proximal Operator

We now show that pressure projection as used in many incompressible flow solvers is equivalent to a proximal operator for the indicator function of divergence-free velocity fields CDIV . Intuitively, this operator corresponds to a projection ΠC DIV (•) onto the set of divergence-free vector fields. This follows from the Helmholtz decomposition, which states that every vector field u * , subject to appropriate boundary conditions, can be uniquely decomposed into a curl-free component, ∇P , defined by the gradient of a scalar function (pressure in our case) and a divergence-free component, ∇× a, defined by the curl of a vector potential a.

The vector identities ∇ × (∇P ) = 0 and ∇ • (∇ × a) = 0 further show that these fields are orthogonal. Incompressible flow solvers commonly exploit this fact by advancing the incompressible Navier-Stokes equations forward in time using all terms from Equation 2 except for the pressure gradient ∇P/ρ to obtain an intermediate velocity field u * . Since the pressure gradient term can only contribute to the divergent components of the velocity field, u * can be made divergence-free by computing a pressure field that exactly counteracts the divergent component introduced by time integration, i.e.:

P = ∇ 2 -1 (-∇ • u * ) (8) u = u * -∇P/ρ, (9)
where ∇ 2 -1 indicates Poisson integration for the pressure. Solving Equations 8-9 in sequence results in a velocity field u that is divergence-free. From the orthogonality of the divergence-free and curl-free component and from the uniqueness of the decomposition it is apparent that the update is the minimum norm change to the velocity field required to make it divergence-free, and therefore, a projection onto the subspace of divergence-free vector fields CDIV .

This implies that the pressure projection step commonly employed in fluid simulators can alternatively be stated as a simple minimization problem. In terms of the proximal operator framework, we may write:

prox λDIV (u * ) = arg min u∈C DIV u -u * 2 2 = ΠC DIV (u * ). ( 10
)
Equation 10 shows that the pressure projection step is a convex operation and thus may be used in convex optimization algorithms. Moreover, the projection onto the sub-space of divergence-free vector fields is a proximal operator, enabling its application as a regularizer in proximal optimization algorithms. We note that Equation 10 may practically be implemented either through a pressure solve as in Equations 8-9 or as a projection operation as in Eq. 10.

The implementation using the pressure solve step benefits from the availability of optimized solvers in the fluid simulation community. Moreover, since a scalar update function is computed, only one third of the variables are needed in the solve as compared to implementing Eq. 10 directly.

Fluid Tracking

As mentioned in Sect. 3, the optical flow brightness constancy formulation and the transport equation for inviscid, incompressible flow are equivalent. In order to perform fluid tracking, we can therefore use the proximal operator derived in the previous section with the 3D equivalent of an optical flow data term (Eq. 3) to formulate incompressible fluid tracking as a minimization problem

min u EP C (u) + αESM (u) + βEKE(u) E F T (u)=F (u) + γEDIV (u) G(u)
.

Here, EP C is the optical-flow photoconsistency (or fluid transport), term, ESM is a smoothness term, EKE is a kinetic energy penalty, and EDIV is a divergence-free term as in the previous section. We have already indicated how this minimization problem maps to the one from Equation 4. In order to apply ADMM to this problem, we have to define the individual terms, and derive the proximal operators for F and G. As we show in the following, these are a velocity estimation problem EF T (u) = F (u) and a pressure projection

EDIV (u) = G(u), respectively.
Technically, we propose several changes to established optical flow algorithms. These modifications significantly enhance the velocity reconstructions compared to previous work, and we evaluate the resulting gain in quality in Section 6.1. The main differences between fluid tracking and standard optical flow can be summarized as:

• the suppression of spurious velocities, Sect. 5.1,

• a divergence-free multi-scale scheme incorporating advection of fluid velocities, Sect. 5.1.1 and 5.2, and

• a sequential reconstruction pipeline that ensures temporal coherence of the flow by coupling the Navier-Stokes equations into the estimation procedure, Sect. 5.3.

Velocity Estimation

First, we derive the velocity estimation subproblem. We use a standard ℓ2 objective function with a second-order smoothness term and a kinetic energy constraint. The data term is chosen according to Eq. 3:

EP C ( u) = 1 2 Ω (φt -u • ∇φ) 2 dΩ, ( 12
)
where φt is the time-derivative of the volume density and ∇φ is its spatial gradient.

We use a voxel basis to represent the continuous space fields, and discretize the image gradient ∇φ with the first-order finite difference matrix [∇Φ] to obtain

EP C (u) = 1 2 (Φt -[∇Φ]u) 2 (13) 1 2 u T [∇Φ] T [∇Φ]u -Φ T t [∇Φ]u, (14)
where the second line omits terms that are constant with respect to u, and therefore not relevant to the minimization problem.

In flow tracking applications we typically encounter regions with very little variation in density (analogous to untextured image regions in optical flow applications). In order to extrapolate velocities to these regions we introduce a smoothness prior as a penalty. This can be expressed in continuous space as a minimization over the spatial gradients of the components of the flow field, i.e.

ESM ( u) = j Ω ∇uj 2 2 dΩ, (15)
which after discretization becomes

ESM (u) = 1 2 j (Dju) 2 = 1 2 j u T Lju, (16)
where D j is a discrete derivative matrix for the j th velocity component, and Lj = D T j Dj is the corresponding discrete Laplacian.

Finally we introduce a Tikhonov regularization term that is physically motivated as a kinetic energy penalty to favor minimally energetic flows:

EKE(u) = 1 2 u T u (17)
The overall objective function F (u) for the velocity estimation subproblem is then the weighted sum of Equations 13, 16 and 17 as indicated in Equation 11. Since all terms of the fluid tracking problem EF T (u) = F (u) are quadratic functions, the corresponding proximal operator turns out to be a simple linear least-squares problem:

prox λF (v) = (I + λAF T ) -1 (v -λbF T ) (18)
AF T = [∇Φ] T [∇Φ] + α j Lj + βI bF T = [∇Φ] T Φt
Note that (I + λAF T ) is a positive definite matrix since it is encoding a proximal operator. Thus, Conjugate Gradient is a feasible solver for this linear system. We use it together with the ILDLTT preconditioner (incomplete Cholesky factorization), caching the preconditioner and warm-starting the solver every time the proximal operator is called by ADMM.

The velocity estimation method derived in this section can be interpreted as a variant on the well known [START_REF] Horn | Determining Optical Flow[END_REF] optical flow algorithm. Instead of using a photo-consistency argument, we have exploited the link to the transport equation (Eq. 3) and incorporated physical constraints on the divergence of the velocity field and the kinetic energy. This last prior can be very important for fluids: without it, Horn-Schunck optical flow will find the smoothest possible solution. Although this is desirable for optical flow which generally has large simple motions, for fluids it is preferable to look for low-energy solutions which helps to reduce spurious velocities and allows lower smoothness penalties to be used. For a visual comparison of the effect, see Fig. 2.

Scale Splitting of Velocity

In preparation for introducing a multi-scale algorithm for the fluid tracking problem, we revisit the derivation of the velocity estimation step, and rewrite it in terms of a low frequency estimate uin and a high frequency perturbation ∆ u: u = uin + ∆ u. In the multi-scale algorithm below, we obtain uin by upsampling a low-resolution solution, and optimize for EF T (∆u) at each scale. Following the same derivation as above, we can derive a new proximal operator, which is identical to the one from Equation 18, except for the right hand side vector bF T , which is now: bF

T = [∇Φ] T Φt + α j L T j u in + βu in

Pressure Projection

To apply the proximal operator from Equation 10 in the context of our splitting of the fluid velocity into an initial estimate and a perturbation, u = uin + ∆ u, we modify the proximal operator slightly to account for the fact that the velocity estimate uin may not be divergence-free:

Algorithm 2 Pressure projection proximal operator, prox λG 1: procedure prox λG (∆u * ) 2:

// u in is coarse scale velocity estimate 3:

u * = u in + ∆u * 4: u = prox λDIV (u * )
⊲ Project u via Equations 8 and 9

5:

return u -u in 6: end procedure We base our implementation on the staggered discretization with co-located velocities from [START_REF] Ando | Highly Adaptive Liquid Simulations on Tetrahedral Meshes[END_REF]]. We use it to discretize anisotropic, rectangular grids, but our algorithm is oblivious to the different types of pressure solvers, as long as pressure projection and optical flow energy are discretized consistently.

Extension to Multiscale

We now describe a multi-scale extension that allows the method to recover large motions. Each scale of the method computes an update ∆ u to the velocity field computed by coarser scales that improves the fit to the input data, returning the sum as output. The multi scale algorithm makes use of three auxiliary functions:

• SmoothAndCubicDownsample(φ, η, σ) blurs φ by a Gaussian of width σ and cubicly downsamples by a factor of η

• CubicUpsample(φ, η) cubicly upsamples by a factor of η

• Advect(φ, u, ∆t) advects φ through the velocity field u by ∆t units of time, using BFECC path tracing with 3rd order spatial interpolation and min-/max-limiting Algorithm 3 Multiscale fluid tracking, 0 < η < 1 is downsampling factor, σ is blur for downsampling, lmax is the number of levels 1: procedure FLUIDTRACKINGMS(φ1, φ2, ∆t, η, σ, l)

2: if l < lmax then 3: φ C 1 ← SmoothAndCubicDownsample(φ1, η, σ) 4: φ C 2 ← SmoothAndCubicDownsample(φ2, η, σ) 5: u C ← FluidTrackingMS(φ C 1 , φ C 2 , ∆t, η, σ, l + 1) 6: u ← 1 η * CubicUpsample(u C , 1 η ) 7: end if 8: φ * 2 ← Advect(φ2, u, -∆t) 9:
// Solve for ∆u using ADMM 10:

// using u in = u in Algorithm 2.

11:

∆u ← ADM M (prox λF , prox λG , φ1, φ * 2 )

12:

return u + ∆u 13: end procedure

The multi-scale algorithm conceptually works by warping the second density backwards in time along the estimated velocity field u computed at the next coarser scale, so that the Taylor expansion used to discretize the fluid tracking time step is valid. Note that ADMM in line 11 of Algorithm 3 only computes an update to the velocities u. Note that u is only used as a constant offset for prox λG to enforce a divergence free total velocity field in Algorithm 2 and is otherwise unused by ADMM.

Our algorithm largely follows the procedure used by [START_REF] Meinhardt-Llopis | Horn-Schunck Optical Flow with a Multi-Scale Strategy[END_REF], except we use ADMM to additionally enforce the divergence-free constraint and use higher-order path tracing for the advection step to improve the results. However, a crucial topic for accurately tracking fluids is temporal coherence of the flow, which is closely related to conservation of momentum. We outline our modifications to embed multi-scale tracking into the fluid simulation context in the following section.

Temporal Coherence

In order to improve tracking quality we first compute an estimate of the next velocity, and then use the fluid tracking algorithm of the previous section to compute an update of these estimated velocities. Based on a velocity from the previous time step, we estimate the next velocity by warping it with itself, and projecting it to make it divergence free (essentially a forward fluid simulation step with operator splitting). This gives us an estimate for the current time step, for which we only have the captured densities φ t . As the constraints we are solving for are convex, using the velocity estimate as initial guess has little influence on the solution. Instead, we warp φ t by -∆t with the velocity estimate to align it with φ t-1 as much as possible, and use our fluid tracker to compute the residual change of the velocities. The final velocity at time t is thus given by the tracking solution plus the velocity estimate. This algorithm is summarized in Algorithm 4. This tight coupling of fluid simulation estimates and tracking has several implications: most importantly, it significantly simplifies the problem the fluid tracking step has to solve. Instead of relying on it to solve the complex, non-linear Navier-Stokes dynamics from Eq. 2, we give it a much simpler problem, namely to only estimate changes in initial conditions, and similar effects that are not covered by the simulation step. This coupling also serves as an implicit temporal coherence prior since the method searches for low-energy updates to the velocities predicted by forward simulation. Even though there is no direct expression of temporal smoothness, velocity fields that are temporally smooth imply lower-energy updates and so are favored over those that 'pop'. This complements the spatial prior, which extrapolates the motions computed in high-texture regions to areas without detail.

Unfortunately the coupling of tracking and simulation also implies that tracking is now a sequential process, in contrast to algorithms that only take two density images as input. In fact, the complete simulation history is contributing to the current time estimate.

Algorithm 4 Fluid tracking, compute next velocity field given previous state φ t-1 , u t-1 and next density φ t 1: procedure COMPUTEVELOCITY(φ t-1 , φ t , u t-1 , ∆t)

2: ũt ← Advect(u t-1 , u t-1 , ∆t) 3:
// project to be div. free with u in = 0 in Alg.

Results and Discussion

In this section we provide quantitative evaluations of the fluid tracking performance, and show several applications of using, resimulating, and modifying the resulting flow fields. All results are based on real captured data of mixing fluids from [START_REF] Gregson | Stochastic Tomography and its Applications in 3D Imaging of Mixing Fluids[END_REF]. These datasets are based on tomographic reconstructions of dye concentrations in mixing fluids. These are very challenging datasets for a number of reasons. First, they were acquired with a very small number of cameras, which limits the amount of interior texture detail that our flow tracking could rely on. Second, as outlined by [START_REF] Gregson | Stochastic Tomography and its Applications in 3D Imaging of Mixing Fluids[END_REF], their method suffers from self-shadowing and other photometric calibration problems, which means that the brightness constancy assumption holds only approximately and on short time scales. Finally, each frame is reconstructed individually, so that the data suffers from small amount of temporal flicker. However, despite these challenges, the data is representative of the best captures of fluid densities described in the literature so far, which motivates our decision to use these datasets in our experiments. Note that single-phase phenomena in air or water can be regarded as equivalent if they are correctly matched in terms of their parametrization and Reynolds number. As such we will treat the flows captured in water as rising smoke in some of the examples below.

Tracking Validation

We compare our proposed method with two common approaches for estimating the motion between two still frames. The first is a multi-scale extension of the well known Horn-Schunck ("H&S") technique that is designed to recover large motions. This choice is motivated by its continued competitiveness when implemented in a modern way [START_REF] Sun | A Quantitative Analysis of Current Practices in Optical Flow Estimation and the Principles behind Them[END_REF]]. The second is an extension to the multi-scale Horn-Schunck method "H&S div. penalty" that incorporates a divergence penalty. This is to mimic the performance of divergence-free methods such as [START_REF] Kadri Harouna | Divergence-free Wavelets and High Order Regularization[END_REF].

Our proposed method is referred to as "ADMM full" and we include a version without kinetic energy constraint "ADMM no kin.". The reconstruction parameters are provided in Table 1.

Our validation experiment consists in generating a simulated ground truth sequence of a 3D incompressible flow (resolution 50 × 100 × 50) transporting a density through a simple domain.

The sequence consists of 118 frames. The central slice 25 through the simulated density volume is shown in Fig. 3 (upper left).

To evaluate the different methods, we investigate two error measures that are adapted from the optical flow community. The velocity norm error measures the norm of the vector difference between the estimated and the ground truth vectors. The angular error measures the directional accuracy of the motion. We weight the angular error with the norm of the velocities to help distinguish large-scale trends from finer features such as a slight shifting of eddies that produce high angular errors even though the overall flow behavior is consistent.

Figure 3 shows a visual comparison of the angular error, velocity norm error, and a color-coded visualization of the velocity vectors. All plots show the same volume slice 25 and have a common scale. We see that our proposed method using tightly coupled simulation and velocity estimation significantly outperforms the other two methods. The error is significantly lower, and our method is the only one that successfully captures the overall structure of the velocity field. This is even more apparent in the supplemental video footage, where the Horn-Schunck and divergence-penalized methods show severe temporal artifacts.

The temporal behavior of the different algorithms is also shown in Fig. 4. The Horn-Schunck derived techniques are susceptible to incoherent estimates whereas our proposed method shows a smooth error evolution that is due to cumulative tracking error (left). The large angular errors at the beginning of the sequence stem from a random velocity initialization of the ground truth simulation (right): Initially, there is no density in the simulation domain which could be used by the algorithm do identify these "hidden" velocities. Around frame 60, the simulation has reached the full domain and the optical flow obtains its highest accuracy. Whereas the directional data recovered by the Horn-Schunck derived methods are Capture H&S H&S div. ADMM The H&S solvers compute a solution to the linear system, corresponding to Eq. 11 using CG, whereas ADMM uses the algorithms described in the paper. The multi-scale parameters are common for all methods: 3 levels, σ =0.5, η =0.65.

almost random, our proposed algorithm computes reasonable estimates.

To test the effect of the kinetic energy regularizer, Fig. 4 also includes a version "ADMM no kinetic" of our proposed algorithm with deactivated kinetic energy term. The numerical performance is worse than the full method for most of the sequence. It is interesting to note that once the data has reached the full domain, the kinetic energy term has a slightly negative effect: when data is available everywhere over-regularization occurs. An adaptive determination of the corresponding parameter may increase the accuracy of our algorithm in future implementations.

We also tested the performance of the method on flows with solid obstacles and with varying degrees of data corruption. Input data was generated for a domain with 100x200x100 cells, consisting of a spherical dye inflow with randomized initial velocities and a spherical obstacle that triggers turbulent mixing (shown in Figure 6). We enforce the solid boundary conditions via the variational approach of [START_REF] Batty | A fast variational framework for accurate solid-fluid coupling[END_REF] in the pressure solve and as a penalty in the velocity estimation subproblem.

We evaluated the fluid tracking performance at the full resolution and on a half-resolution grid with varying levels of blur added to the input. This was done to approximate the quality of data obtained by optical tomography capture setups (e.g. as used by [START_REF] Gregson | Stochastic Tomography and its Applications in 3D Imaging of Mixing Fluids[END_REF]) which typically yield fairly low spatial resolution with substantial filtering. The velocity norm errors and angular errors are plotted in Figure 6.

The results in Figure 6 show that reasonable velocity predictions may be obtained even from coarse input volumes with moderate levels of blur. High levels of blur cause the reconstructions to diverge more quickly over time in both metrics.

Validation with Resimulation

Having developed a method to estimate the velocity fields of incompressible flows and assessed its accuracy on simulated data, we are now able to restart fluid simulations using initial conditions computed from captured data. In this way we are trying to answer a question that is interesting for scientific applications of our method: based on our velocity reconstructions, how close can we get to the real phenomenon with a simulation? We begin with some discussion on what can be expected in such a situation.

The results in Section 6.1 demonstrate that tightly coupling simulation and fluid tracking significantly improves the accuracy of the computed velocity fields. Passive advection of marker fields through the velocity fields computed by our method confirms that results reproduce the medium and coarse scales of the velocity fields well, Fig. 5. We see that the majority of significant flow features are reproduced in the advected marker fields, despite being transported over several seconds through turbulent mixing. In comparison, an advection with the velocities computed by the Horn-Schunck derived techniques are not able to reproduce the behavior of the capture. Temporal tracking artifacts lead to flicker in the resimulations, which is most easily seen in the supplemental video. It is well known that turbulent flows exhibit considerable energy transfer between spatial scales in the velocity field, including scales that are significantly finer than those that can be resolved by flow tracking. We consequently expect that the resimulated behavior matches better for flows without effects such as buoyancy that cause high frequency behavior. In the former case, numerical and physical viscosity rapidly dampens out high frequencies that cannot be estimated accurately, an effect that is advantageous for longer resimulation similarity. This expectation is confirmed in the resimulated "bloom" capture, Figure 7 (a), which involves dye being poured into still water. Although high-frequencies are quickly lost, the medium and large scale behavior persists for over a second of turbulent mixing since there are no density gradients driving the flow. The energy transfer in this case is primarily coarse-to-fine, and our resimulation gives a motion close to the real flow.

The opposite effect can be observed for the "smoke" capture, shown in Figure 7 (b). Here the buoyancy term results in kinetic energy being introduced in high spatial frequencies of the flow, resulting in energy being transferred from fine-to-coarse scales. The captured data is under-resolved for this type of flow and diverges more quickly. As such it highlights a limitation that our method shares with many other methods for simulating fluids: when our solver and the chosen discretization are unable to represent the physics correctly, we obtain results that deviate from the desired behavior, although they can sometimes be improved by optimizing simulation parameters.

Stylistic Capture Modifications

A benefit of computing velocity fields for capture data is that now the full flow parameters are known and can be used for artistic modifications. This measure makes captured flows amenable for use in fluid simulations. In the following, we demonstrate exemplary changes that can be applied.

Resolution Enhancement can be achieved by two means: The first is a resimulation of the captured velocity fields in conjunction with passive advection that carries initial frame densities forward in time. The availability of full flow information enables a resampling in both space and time. For an illustration of this approach, see Figs. 8 and9. The apparent resolution increase can be appreciated by comparing the columns 'Passive Advection' and 'Tomography'. The advected densities are not guaranteed to match the capture partially because we model the density inflow conditions from a single frame of the sequence, which can lead to changes of the overall shape (e.g. on the left side of Figure 9).

The second approach for apparent resolution increase is the application of a synthetic turbulence model. As an example, we apply wavelet turbulence [START_REF] Kim | Wavelet Turbulence for Fluid Simulation[END_REF] in the column labeled 'Superresolved' (calculated with mantaflow [START_REF] Pfaff | MantaFlow[END_REF]).

Since the original characteristics of the flow are not changed, the

Frame 40

Frame 65

Frame 97 Domain Change can be affected by running the resimulation with advection (discussed above) up to a desired time-frame, after which the simulation is switched to a dynamic setting that takes into account the new boundary conditions. This way, the computed velocity fields can be used to transfer captured flows into modified domains. Figure 10 shows the 'smoke' dataset resimulated after adding a large sphere to the domain. The resimulation is robust despite a significant modification of the domain. Guided Simulation can be achieved by exploiting a second connection to optimization techniques from image processing: Given an initial simulation with velocity fields u, and a target state u ref that may e.g. result from flow tracking, we seek an updated state u + ∆u of the original simulation that reproduces the coarse scale behavior of the target flow while maintaining some characteristics of the original simulation.

We may formulate the coarse scale requirement as a low frequency blur affected by a matrix B on the simulation data. With these prerequisites, we can express a new optimization that seeks the necessary update ∆u as:

arg min ∆u 1 2 (B∆u -B(u ref -u)) 2 2 (19) +αESM (∆u) + βEKE(∆u) + γEDIV (u + ∆u),
i.e. the update ∆u, when blurred, should mimic the differences in the original and the target flows, as before subject to physical constraints. The data term now resembles a deblurring problem, but, as outlined above, has a different interpretation in fluid problems. The regularizers, again, ensure physicality of the flow. The proximal operators only change in the data term, which can be derived in a similar manner to Eq. 18.

The fluid guiding scheme is then performed by interleaving a simulation and the guiding step Eq. 19 until the full sequence is computed. An example of this application is shown in Fig. 11. Note that for this example we are effectively running the ADMM framework twice, first to reconstruct the guiding velocities, and then a second time to guide the dynamic simulation.

Limitations

The present work has several limitations that bear consideration. The most significant one is that the method depends on the flows having time-varying detail that can be used to infer the underlying velocities. Consequently the performance of the method suffers when such texture is unavailable. This occurs most frequently in laminar flows which develop a static or near-static distribution of the marker fluid. This limits the applicability of the approach in such areas as laminar boundary layers.

Less severe limitations include the fact that through the use of simulation as a prediction step frames must be reconstructed in sequence, preventing reconstructing multiple frames in parallel. Finally the weights of the smoothness and kinetic energy penalties must be tuned for a given flow and capture setup to achieve a good balance between reconstructed detail and noise levels.

Conclusions & Future Work

In this paper we have shown that the pressure solve that is integral to many fluid simulators can be cast as a proximal operator. This formulation allows us to implement highly optimized versions of this operator, and re-use it in another context using a proximal optimization framework such as ADMM.

There are many applications for this approach. In this paper we have focused primarily on tracking fluid velocities based on densities captured by prior art fluid imaging methods. We have shown that these tracked fluid velocities can serve as the basis for resimulation, for example to add details not captured in the original data (fluid super-resolution), or to change the domain of the fluid flow. We have also shown some initial results of a guided simulation application of the same framework, which is akin to a deblurring problem in image processing. These applications show how the proximal operator framework can be used to connect forward and inverse problems in fluids, and therefore provides an avenue for using captured fluid data in graphics applications.

We believe that the applications we have explored so far are only touching upon the possibilities offered by this approach. In general, casting the fluid problems into this convex optimization framework using proximal operators allows us to separately implement and optimize different components of both forward and inverse problems, and to combine them in a modular and efficient way. We believe that it is possible to implement additional local constraints on the flow field as proximal operators, for example to allow for flexible boundary conditions such as the wall separation boundaries proposed by [START_REF] Batty | A fast variational framework for accurate solid-fluid coupling[END_REF]. This approach could be used to implement a valve, or to allow fluid flows to peel of the surface of a solid in a fluid simulation. We believe that many more such applications are possible.

Figure 2 :

 2 Figure 2: Comparison of 3D reconstructed velocity fields to ground truth data (left) for the full ADMM method that includes kinetic and divergence-free penalties (middle left), the same method without kinetic energy penalty (β = 0, middle right), and multi-scale Horn-Schunck flow (β = 0, γ = 0, no ADMM, far right). The figure shows three different time steps (20, 25 and 30). The full method clearly reproduces the features of the ground truth velocity pattern, whereas the method without kinetic energy penalty invents new structures (lower right corner). Standard Horn-Schunck flow produces smooth flow fields that only on average correspond to the ground truth. For visualization, the 3D velocity vectors are converted to spherical coordinates which are then mapped to HSV color space. The figure shows a central slice of the 3D domain.

  2

Figure 3 :

 3 Figure 3: A visual comparison between Horn-Schunck (H&S), divergence-penalized H&S, and our proposed method. From top to bottom: magnitude of vector difference w.r.t. ground truth, angular error, absolute HSV color-encoded direction plots. The figure shows the central slice of our simulation volume.

Figure 4 :

 4 Figure 4: Ground truth comparison for a 3D simulation for different methods. We show velocity norm error (left) and norm-weighted angular error (right) plots. The Horn-Schunck-based two-frame methods suffer from temporal artifacts. The direction estimates of the flow vectors are unsatisfactory. In contrast, our proposed method shows very good temporal coherence and significantly improved directional estimates.

Figure 5 :

 5 Figure 5: Passive advection of densities in a reconstructed timevarying velocity field. From left to right: captured densities, result based on Horn-Schunck reconstruction, H&S with divergence penalty, and our proposed algorithm.

Figure 6 :

 6 Figure 6: Fluid tracking for a flow interacting with a solid obstacle. Top row: dye enters from a spherical inflow and impinges upon a spherical obstacle. A slice through the domain shows dye concentration (fringe color), ground-truth velocity vectors (black) and reconstructed velocity vectors (white). Bottom row: velocity norm and angular errors for native resolution ("original") and half-resolution reconstructions with varying levels of blur (σ blur = [0.0, 0.5, 1.0] voxels) applied to the input density fields.

Figure 7 :

 7 Figure 7: Resimulation of the bloom and smoke datasets at 10 frame intervals. Despite fine-scaled turbulent mixing, the simulator preserves many features of the flow for over a second. In contrast, the smoke capture with its buoyancy induced instabilities causes more rapid divergence of the simulated flow.

Figure 8 :

 8 Figure 8: Proposed fluid-tracking applied to the 'smoke' capture consisting of dye and alcohol rising under buoyancy in still water. Left column: tomographically reconstructed concentration fields, center column: passive advection using fluid tracking velocities, right column: additional synthetic turbulence. The reconstruction was performed on a 50 × 150 × 50 grid, and took 2.35 minutes per frame on average. global behavior of the capture remains intact.

Figure 9 :

 9 Figure 9: Proposed fluid-tracking applied to the 'bloom' capture consisting of dye poured into still water. Left column: tomographically reconstructed concentration fields, center column: passive advection using fluid tracking velocities, right column: additional synthetic turbulence. The flow velocities were reconstructed on a 100 × 150 × 100, taking 12.5 minutes on average per frame.

Figure 10 :

 10 Figure 10: The ability to estimate fluid velocities allows us to restart simulations in modified domains as compared to the capture setting. The example uses the 'smoke' data set.

Figure 11 :

 11 Figure11: As an application of the guided simulation approach, we add high-frequency curl noise to a resimulation of the smoke dataset in order to affect its visual appeal. Since no measures are taken to avoid affecting the global flow behavior, the overall lowfrequency characteristics of the captured flow is lost (left). Incorporating the guided simulation approach helps to recover the lowfrequency motions (right).

Table 1 :

 1 Optimization parameters for ground truth comparisons.

'n.a.' stands for 'not applicable', 'h.c.' for 'hard constraint'.

Acknowledgements

This was in part supported by NSERC and the GRAND NCE, Ivo Ihrke has been supported by the German Research Foundation (DFG) through Emmy-Noether Grant IH 114/1-1.


Serap Kırbız

Alexey Ozerov

Antoine Liutkus

Laurent Girin

PERCEPTUAL CODING-BASED INFORMED SOURCE SEPARATION

Keywords: Informed source separation, source coding, perceptual models

Informed Source Separation (ISS) techniques enable manipulation of the source signals that compose an audio mixture, based on a coder-decoder configuration. Provided the source signals are known at the encoder, a low-bitrate sideinformation is sent to the decoder and permits to achieve efficient source separation. Recent research has focused on a Coding-based ISS framework, which has an advantage to encode the desired audio objects, while exploiting their mixture in an information-theoretic framework. Here, we show how the perceptual quality of the separated sources can be improved by inserting perceptual source coding techniques in this framework, achieving a continuum of optimal bitrateperceptual distortion trade-offs.

INTRODUCTION

Active listening of music and audio consists in interacting with the different audio objects that compose a mixture signal. This goes from generalized karaoke to respatialization and audio scene content modification. Since only the mixture signal is usually available at the user's side, signal processing techniques must be used to recover the individual audio objects from the mix. However, blind or semi-blind methods [START_REF] Comon | Handbook of Blind Source Separation: Independent Component Analysis and Blind Deconvolution[END_REF] are still not efficient enough on complex mixtures to be embedded in large audience products, and solutions based on a coder-decoder configuration have been recently proposed. Provided the source signals are known at the encoder, a low-bitrate side-information is sent to the decoder -in addition to the mix signal itself -and used to achieve efficient source manipulation. In this line, the MPEG Spatial Audio Object Coding (SAOC) technology has been derived from multichannel spatial coding and makes use of spatial cues [START_REF] Engdegård | MPEG spatial audio object coding -the ISO/MPEG standard for efficient coding of interactive audio scenes[END_REF]. In parallel, the coder-decoder configuration has This work is partly funded by the French National Research Agency (ANR) as a part of the DReaM project (ANR-09-CORD-006-03) and partly supported by the European Research Council (Grant 278025), the Emergence(s) program from the City of Paris, and LABEX WIFI (Laboratory of Excellence within the French Program "Investments for the Future") under references ANR-10-LABX-24 and ANR-10-IDEX-0001-02 PSL*.

been introduced in the source separation framework, leading to informed source separation (ISS) systems [START_REF] Liutkus | Informed source separation through spectrogram coding and data embedding[END_REF]. In a general manner, the side-information consists of models of mixing parameters and source power spectral densities (PSD) which are used to build demixing filters at the decoder [START_REF] Liutkus | Informed source separation through spectrogram coding and data embedding[END_REF].

Although efficient, these methods have their performance bounded by the best estimates that can be provided by the separation method and source model. Hence, a large bitrate spent on the separation parameters does not improve consequently the quality of the estimates. Recently, the ISS problem was reformulated in [START_REF] Ozerov | Coding-based informed source separation: Nonnegative tensor factorization approach[END_REF] so as to consistently benefit from additional bitrate as in source coding [START_REF] Gray | Source coding theory[END_REF][START_REF] Zhao | On entropy-constrained vector quantization using Gaussian mixture models[END_REF]. Coding-based ISS (CISS) was proposed, based on a probabilistic model as in source coding. However, instead of using a prior distribution which does not make use of the mixture, CISS encodes the sources relying on their posterior distribution given the mixture, whose smaller entropy leads to reduced bitrates. CISS thus provides a flexible scheme that "unifies" source separation and multichannel source coding, and allows for fine and extended tuning of the trade-off between side-information bitrate and separation quality. Now, CISS can be further improved. In particular, no perceptual considerations have been introduced so far in this framework, though it is well-known that audio coding techniques benefit a lot from the perceptual approach [START_REF] Van De Par | A perceptual model for sinusoidal audio coding based on spectral integration[END_REF]. In the present paper, we propose and experiment a first series of strategies to introduce perceptual models and criteria in the CISS scheme. We report the improvements resulting from those new developments.

This paper is structured as follows. First, we rapidly present the CISS framework in Section 2. Section 3 presents the integration of perceptual models into the CISS framework. Finally, the proposed method is evaluated in Section 4.

CODING-BASED ISS

All signals are processed in a Time-Frequency (TF) representation, for instance the Modified Discrete Cosine Transform (MDCT), since it is critically sampled and thus a good candidate for source coding. Let

s ωt = [s 1ωt • • • s Jωt ]
⊤ be the J × 1 column vector gathering all J sources for TF bin (ω, t),

where ω ∈ [1, N ω ] and t ∈ [1, N t ] respectively denote the frequency and frame index, and • ⊤ denotes transposition. In the present study, the mixture is assumed to be a linear combination of the sources, with 1 × J mixing vector A:

x ωt = As ωt , ∀ (ω, t) .

(

This single-channel case is considered here for simplicity of presentation, but the model can be extended to handle many complex multichannel mixing scenarios [START_REF] Liutkus | Low bitrate informed source separation of realistic mixtures[END_REF].

The sources are assumed to be independent and to follow a Local Gaussian Model (LGM) [START_REF] Liutkus | Gaussian processes for underdetermined source separation[END_REF]. Under this assumption, all TF bins of s j are independent and normally distributed, thus yielding the following prior distribution:

s ωt ∼ N (0, C ss,ωt ) , ( 2
)
where N denotes the Gaussian distribution and

C ss,ωt = diag{v 1ωt , • • • , v Jωt } is the prior source covariance matrix.
The variance v jωt of source j at TF bin (ω, t) can be understood as its PSD. The LGM is hence parametrized by the J × N ω × N t sources PSD tensor V = {V jωt } j,ω,t . Assuming that a model V of V is transmitted to the decoder along with the mixture x and the mixing coefficients, the (modelled) variance of the mixture at TF bin (ω, t) is

C xx,ωt = AC ss,ωt A ⊤ = ∑ J j=1 A 2 j Vjωt . ( 3
)
It is straightforward to show that the posterior distribution of s ωt given x ωt , V and A writes:

s ωt | x ωt , V, A ∼ N (µ ωt , K ωt ) , (4)
where:

µ ωt = C ss,ωt A ⊤ C -1 xx,ωt x ωt , K ωt = C ss,ωt -C ss,ωt A ⊤ C -1
xx,ωt AC ss,ωt .

(

A first trend of research on ISS [START_REF] Liutkus | Informed source separation through spectrogram coding and data embedding[END_REF] has focused on using V and A at the decoder to build Wiener demixing filters, which exactly corresponds to estimate the sources s ωt by their posterior mean µ ωt in [START_REF] Gray | Source coding theory[END_REF]. This leads to a computationally efficient ISS method with good but bounded separation quality. To reach arbitrary distortion on separated sources at the cost of higher bitrate, CISS has been proposed in [START_REF] Ozerov | Coding-based informed source separation: Nonnegative tensor factorization approach[END_REF] and has the important advantage of being optimal from an informationtheoretical point of view. The main idea of CISS is to apply traditional (lossy) source-coding [START_REF] Gray | Source coding theory[END_REF] on s ωt , with the important difference that the considered distribution is the posterior distribution (4) instead of the classical prior distribution, for instance [START_REF] Engdegård | MPEG spatial audio object coding -the ISO/MPEG standard for efficient coding of interactive audio scenes[END_REF]. This difference leads to an important decrease of the source coding bitrate, due to the mutual information between the sources s and the mixture x which is available at the decoder. Given this fundamental difference, the practical implementation follows the usual fundamentals of source coding (see [START_REF] Ozerov | Coding-based informed source separation: Nonnegative tensor factorization approach[END_REF] for a detailed description). In short, the Algorithm 1 Coding-based ISS for squared error. For all TF bins (ω, t) (which are omitted here for conciseness): 1. Compute µ and K as in (5).

Compute the eigenvalue decomposition

K = Udiag {[λ 1 , • • • , λ D ]} U T ,
where U T is the KLT.

Compute

z = U T (s -µ) .
4. Quantize each dimension of z using a uniform quantizer of constant step-size ∆ s to yield quantized z. Using an arithmetic coder as an entropy coder [START_REF] Zhao | On entropy-constrained vector quantization using Gaussian mixture models[END_REF], the effective codeword length (in bits) is given by:

- ∑ J j=1 log 2 ∫ zj +∆s/2 zj -∆s/2 N (z j | 0, λ j ) dz j ,
where N (• | 0, λ j ) is the probability density function of the zero-mean normal distribution with variance λ j .

5. Quantized source vector is given by s = Uz + µ.

sources s ωt are encoded using model-based (for instance (4)) constrained entropy quantization based on scalar quantization in the mean-removed Karhunen-Loeve Transform (KLT) domain [START_REF] Zhao | On entropy-constrained vector quantization using Gaussian mixture models[END_REF]. The user chooses a step-size ∆ s to control the ratedistortion trade-off and Algorithm 1 is applied, while using an arithmetic coder for entropy coding. Finally, the problem of estimating V and transmitting it to the decoder has to be addressed (the cost of transmitting A is assumed to be negligible). Due to its important number of entries, factorization and compression techniques are required to model and encode V.

Previous studies [START_REF] Liutkus | Informed source separation through spectrogram coding and data embedding[END_REF][START_REF] Ozerov | Coding-based informed source separation: Nonnegative tensor factorization approach[END_REF][START_REF] Liutkus | Low bitrate informed source separation of realistic mixtures[END_REF] have concentrated on Nonnegative Tensor Factorization (NTF) [START_REF] Nikunen | Multichannel audio upmixing based on non-negative tensor factorization representation[END_REF], which drastically reduces the number of parameters by assuming that V is the superposition of a relatively small number K of rank-1 tensors:

Vjωt = ∑ K k=1 W ωk H tk Q jk , (6)
where W, H and Q are N ω ×K, N t ×K and J ×K nonnegative matrices, respectively. Learning these model parameters is achieved through the maximum likelihood (ML) strategy that can be shown equivalent to the following optimization criterion [START_REF] Févotte | Nonnegative matrix factorization with the Itakura-Saito divergence. With application to music analysis[END_REF]:

{W ⋆ , H ⋆ , Q ⋆ } = argmin {W,H,Q} ∑ j,ω,t d IS ( V jωt | Vjωt ) , (7)
where d IS is the Itakura-Saito divergence (see, e.g., [START_REF] Févotte | Nonnegative matrix factorization with the Itakura-Saito divergence. With application to music analysis[END_REF] for a definition). This optimization problem can be tackled with classical NTF multiplicative updating [START_REF] Liutkus | Gaussian processes for underdetermined source separation[END_REF], i.e., Algorithm 2 (with P jωt = 1). Once {W, H, Q} are estimated, they are quantized to be transmitted to the decoder. It has been demonstrated in [START_REF] Ozerov | Coding-based informed source separation: Nonnegative tensor factorization approach[END_REF] that nearly optimal performance can be achieved by uniform quantization of log W, log H and log Q. The respective quantization steps are provided in [START_REF] Ozerov | Coding-based informed source separation: Nonnegative tensor factorization approach[END_REF]. The resulting indices are entropy encoded, e.g. by a Gaussian Mixture Model-based arithmetic coding as in [START_REF] Ozerov | Coding-based informed source separation: Nonnegative tensor factorization approach[END_REF].

PERCEPTUAL CISS

It is well known from audio coding that optimizing the Mean Squared Error (MSE) is far to be the best strategy to achieve the best perceptual quality possible at a given bitrate. Instead, audio coding is more efficient when the distortion measure to be minimized between compressed and original data includes some perceptual weighting [START_REF] Ozerov | Flexible quantization of audio and speech based on the autoregressive model[END_REF]. Here, we show how this principle can be exploited in the CISS framework to yield an optimal CISS strategy that minimizes both the bitrate and the perceptual distortion. Actually, we propose to integrate perceptual models into the CISS scheme either at the NTF model estimation step (Algorithm 2), or at the posterior source coding step (Algorithm 1), or at both steps. The corresponding perceptually motivated CISS-NTF schemes are referred to as CISS-PNTF, PCISS-NTF, and PCISS-PNTF, respectively.

Perceptual Models and associated Weighting Tensor

Implementing a psychoacoustic model into the CISS framework amounts in defining a perceptual weighting tensor P = {P jωt } j,ω,t (P jωt ≥ 0) that gives the perceptual importance of each TF bin (ω, t) of each source j. Depending on the perceptual model considered, the resulting coefficients P jωt are computed differently, but all further processing may then be understood as optimizing a perceptually weighted quadratic distortion, which is classical in the audio coding literature. More precisely, while parameter learning presented in section 3.2 below is amenable to a weighted log-likelihood optimization, perceptual source coding is achieved by identifying P with a sensitivity matrix, as presented in [START_REF] Van De Par | A perceptual model for sinusoidal audio coding based on spectral integration[END_REF][START_REF] Petkov | The sensitivity matrix for a spectral auditory model[END_REF]. We first rapidly present the two auditory perceptual models that we tested in this study. The first model, called here the Par model was proposed by van de Par et al. in [START_REF] Van De Par | A perceptual model for sinusoidal audio coding based on spectral integration[END_REF] and the second model is defined in ITU-R BS.1387 [START_REF]Method for Objective Measurements of Perceived Audio Quality[END_REF].

The Par model [START_REF] Van De Par | A perceptual model for sinusoidal audio coding based on spectral integration[END_REF] is a psychoacoustic masking model based on the peripheral bandpass filtering properties of the human auditory system (HAS). It is based on frequency masking only and does not take into account temporal masking. It evaluates the distortion-to-masker ratio within each auditory filter which can be used for the derivation of a masking threshold and provides a measure for distortion detectability.

In the presence of a source power spectrogram |s jωt |2 , the weighting tensor is obtained as [START_REF] Petkov | The sensitivity matrix for a spectral auditory model[END_REF]:

P 2 jωt = C s L N ω ∑ i |h ω | 2 |g i ω | 2 1 Nω ∑ ω ′ |h ω ′ | 2 |g i ω ′ | 2 |s jω ′ t | 2 + C a , (8)
where h ω is the transfer function of the outer-middle ear filter, g i ω is that of the i-th gamma-tone filter, L is the effective duration of the corresponding block, and C s and C a are some model calibration constants (see [START_REF] Van De Par | A perceptual model for sinusoidal audio coding based on spectral integration[END_REF][START_REF] Petkov | The sensitivity matrix for a spectral auditory model[END_REF] for details).

The ITU-R model recommended in ITU-R BS.1387 is a spectro-temporal model and relies on both spectral and temporal masking [START_REF]Method for Objective Measurements of Perceived Audio Quality[END_REF]. The quantitative significance of an auditory object within a mixture can be measured by its perceived loudness through the HAS. The loudness of one frame is modelled by calculating the excitation using perceptually motivated frequency scale (Bark scale) and critical bandwidth, compressing the excitation, and integrating over frequency [START_REF]Method for Objective Measurements of Perceived Audio Quality[END_REF].

The loudness coefficients Pjbt are calculated in each critical band b of ITU-R BS. 1387 model as [START_REF]Method for Objective Measurements of Perceived Audio Quality[END_REF]

Pjbt = c ( E τ b τ b E 0 ) 0.23   ( 1 -τ b + τ b Ẽs jbt E τ b ) 0.23 -1   , ( 9
)
where τ b is the threshold index, E t b is the excitation threshold, E 0 is a constant and Ẽs jbt are excitation patterns [START_REF]Method for Objective Measurements of Perceived Audio Quality[END_REF]. The weighting coefficients calculated in the critical bands are extended in the linear frequency scale of the STFT such that P jωt = Pjbt , for the frequency components ω lying in the b-th critical band [START_REF] Kırbız | Perceptually enhanced blind single-channel music source separation by non-negative matrix factorization[END_REF]. It can be seen that, loudness is affected by parameters other than sound pressure, including frequency and duration.

Perceptually weighted model estimation (PNTF)

In this study, instead of the ML approach [START_REF] Van De Par | A perceptual model for sinusoidal audio coding based on spectral integration[END_REF], we propose to estimate the model parameters through a perceptually weighted NTF, which has recently been considered both for multichannel audio coding, and upmixing [START_REF] Nikunen | Multichannel audio upmixing based on non-negative tensor factorization representation[END_REF]. Perceptuallyweighted NTF amounts to estimate {W, H, Q} by minimizing the following cost function instead of (7):

{W ⋆ , H ⋆ , Q ⋆ } = argmin {W,H,Q} ∑ jωt PjωtdIS ( Vjωt | Vjωt ) . (10)
The rationale behind the use of ( 10) is that all TF bins are not of equal perceptual importance. Hence, focusing on parameters that provide good approximation for perceptually important TF bins may be a good way to optimize performance. The optimization procedure for estimating the parameters through multiplicative updates is given in Algorithm 2.

Perceptual posterior source coding (PCISS)

Let P = { Pjωt } j,ω,t be the weighting tensor coefficients computed as in Section 3.1, its quantized version or its estimation 1 , using Vjωt as the sources PSD. To apply this weighting Algorithm 2 Weighted-NTF algorithm

• Inputs: V jωt = |s jωt |
2 , P jωt and K ∈ N

• Initialize W, H and Q randomly. Iterate several times:

Q jk ← Q jk ( ∑ ωt P jωt W ωk H tk V jωt V -2 jωt ∑ ωt P jωt W ωk H tk V -1 jωt ) , ( 11
)
W ωk ← W ωk ( ∑ jt P jωt H tk Q jk V jωt V -2 jωt ∑ jt P jωt H tk Q jk V -1 jωt ) , ( 12
)
H tk ← H tk ( ∑ jω P jωt W ωk Q jk V jωt V -2 jωt ∑ jω P jωt W ωk Q jk V -1 jωt ) . ( 13
)
matrices for quantization, i.e., to optimize a weighted distortion instead of the MSE, the entries of distribution ( 4) should be modified as follows

s ′ ωt = diag{ Pωt }s ωt , ( 14
) µ ′ ωt = diag{ Pωt }µ ωt , ( 15
) K ′ ωt = diag{ Pωt }K ωt diag{ Pωt } ⊤ , ( 16
)
where Pωt = [ P1ωt , . . . , PJωt ] ⊤ . After this modification the source vectors should be quantized exactly as described in Algorithm 1. The quantized source vectors will be weighted. Thus, the unweighted coefficients should be obtained by sωt = diag{ Pωt } -1 s′ ωt .

EXPERIMENTS

In this section, we evaluate the performance of the proposed perceptually motivated CISS-NTF method for two different auditory models. For this purpose, a set of 14 excerpts sampled at 44.1kHz from the QUASI database 2 is used. Each excerpt is a linear instantaneous mixture of 5 to 10 sources with a duration of approximately 30s long. The performance is evaluated using the Normalized Signal to Distortion Ratio (NSDR, in dB) and the Normalized Perceptual Similarity Measure (NPSM, between 0 and 1) which measure the improvement of the SDR [START_REF] Vincent | Performance measurement in blind audio source separation[END_REF] of BSSEval and PSM of PEMO-Q [START_REF] Huber | PEMO-Q -a new method for objective audio quality assessment using a model of auditory perception[END_REF], respectively, over a "do nothing separation", where the mixture itself is considered as a source estimate. NSDR and NPSM are calculated as NSDR j ( ŝj , sj , x) = SDR j ( ŝj , sj ) -SDR j (x, sj ), [START_REF] Huber | PEMO-Q -a new method for objective audio quality assessment using a model of auditory perception[END_REF] NPSM j ( ŝj , sj , x) = PSM j ( ŝj , sj ) -PSM j (x, sj ) , (18) where sj , ŝj and x denote the original source signal, estimated source signal and the mixture signal in time domain, respectively. The simulations produced many (rate, NSDR) and (rate, NPSM) pairs. Then, for each small range of rates, the pairs corresponding to the highest NSDR and NPSM are selected. The resulting (rate, NSDR) and (rate, NPSM) planes were then smoothed using the locally weighted scatterplot smoothing (LOESS) method to produce the rate/performance curves. The results can be found in Fig. 1 and audio excerpts may be downloaded on the webpage dedicated to this paper 3 .

Although we evaluated the performance using both Par and ITU-R models for extracting the weighting tensor, we reported the results corresponding to the higher NPSM. As such we report the results with the ITU-R model for PNTF and with the Par model for PCISS. If we use the loudness matrix of ITU-R model for weighting the distortion, it forces the louder components to have less distortion, while the distortion in the regions with low energy is increased. Thus, we observe an improvement in terms of NSDR around 0.5 dB for CISS-PNTF and PCISS-PNTF compared to CISS-NTF and PCISS-NTF, respectively. Therefore, applying a perceptual weight on the NTF scheme leads to an improvement in terms of NSDR. At the same time, it keeps the NPSM scores unchanged, which is a bit deceiving. On the other hand, we can see that PCISS-NTF and PCISS-PNTF are outperformed by CISS-NTF and CISS-PNTF in terms of NSDR, especially for bitrates above 6.5 kbps, up to approximately 0.5 dB. This is trivial and predictable since NSDR is optimized with the MSE minimization obtained by the CISS scheme. The important point here is that the NPSM scores are improved consistently by approximately 0.01-0.02 by using PCISS with the Par model.

CONCLUSION

In this study, we have shown how perceptual models can be included in coding-based informed source separation (CISS), so as to yield an information-theoretical optimal way to transmit audio objects when their mixtures are known at the decoder. This work extends the state of the art about CISS that previously only considered squared-error as a distortion measure.

Including perceptual weighting in informed source separation was achieved from two different perspectives. First, we showed how the estimation of the Nonnegative Tensor Factorization source model can benefit from perceptual weighting at no cost in terms of bitrate. Second, we have demonstrated that recent results on source coding using perceptual distortions could be extended to the case of posterior sourcecoding, i.e., when both the coder and the decoder share some side-information.

This approach is the first Audio Object Coding approach that permits to include perceptual constraints on the audio sources to be transmitted through their downmix. Indeed, even if other approaches from the state of the art do rely on residual perceptual coding [START_REF] Engdegård | MPEG spatial audio object coding -the ISO/MPEG standard for efficient coding of interactive audio scenes[END_REF], the perceptual models they consider are not relative to the audio objects themselves, but rather to the errors performed during estimation, which is different and suboptimal. We have demonstrated through our evaluation that the system we propose permits to recover isolated sources that are of better perceptual quality with a bitrate of only a few kilobits per seconds and per source.

Fig. 1 .

 1 Fig. 1. Performance of the proposed perceptual CISS-NTF scheme.

Note that in contrast to weighed model estimation presented in section 3.2, where the weighting tensor needs only to be computed at the encoder side, for the perceptual posterior source coding the weighting tensor needs to be known on both the encoder and the decoder sides. Thus, the weighting tensor should be either transmitted or estimated somehow. In this work we estimate it by replacing source power spectrograms |s jωt

| 2 (e.g., in[START_REF] Liutkus | Low bitrate informed source separation of realistic mixtures[END_REF]) by |µ jωt | 2 , with µ ωt from (5).

http://www.tsi.telecom-paristech.fr/aao/en/ 2012/03/12/quasi/

http://www.gipsa-lab.grenoble-inp.fr/~laurent. girin/demo/PCISS-EUSIPCO2014.zip