Fluid Simulation with an L0 Based Optical Flow Deformation

: Fluid simulation can be automatically interpolated by using data-driven ﬂuid simulations based on a space-time deformation. In this paper, we propose a novel data-driven ﬂuid simulation scheme with the L 0 based optical ﬂow deformation method by matching two ﬂuid surfaces rather than the L 2 regularization. The L 0 gradient smooth regularization can result in prominent structure of the ﬂuid in a sparsity-control manner, thus the misalignment of the deformation can be suppressed. We adopt the objective function using an alternating minimization with a half-quadratic splitting for solving the L 0 based optical ﬂow deformation model. Experiment results demonstrate that our proposed method can generate more realistic ﬂuid surface with the optimal space-time deformation under the L 0 gradient smooth constraint than the L 2 one, and outperform the state-of-the-art methods in terms of both objective and subjective quality.


Introduction
Fluid simulations are increasingly important in computer graphics and have made an impact in real-time games. The simulation algorithms typically include physics-based methods and data-driven approaches. The physics-based fluid solvers [1][2][3][4] for working with the large amounts of simulation data are challenging, as the simulation data is typically just passed on to a rendering stage and any required change means re-starting a new simulation from scratch [5]. A category of data-driven fluid simulations [6] interpolates between two existing fluid animations based on a space-time deformation. Once the matches of two 4D shapes are computed, a user can freely choose any version in between the two extremes efficiently, without starting a new simulation.
Matching two fluid surfaces means finding the corresponding points between the source surface and the target surface, and this procedure needs to deform the source onto the target. To choose the suitable overall match, the deformation should be regularized with some constraints. A novel optical flow approach [5] is to register deforming space-time source surfaces and target surfaces for the robust registration and blend the complex volumetric phenomena by minimizing energy function with the L 2 norm smoothness regularizer and the Tikhonov regularizer, which is inherently based on closest distances in four-dimensional Euclidean space. The L 2 norm based smoothness regularization in [5] can be enforced on the gradient of the deformations to suppress the differences of the neighbors, however, this method can not always require an ideal deformation for the large differences of the source and target surfaces. An example of this is demonstrated in Figure 1, where the droplets of water with different shapes are released at different positions. As the errors of the deformation run into the registration procedure, the inaccurate matching caused by the sub-optimal deformation will create ambiguous corresponding points between the surface of the pool and and that of the droplet, which brings out the artifacts in the interpolated results.  [5]. While final output retrieves large scale motions, the artifacts between the droplet and surface are noticeable.
We propose an L 0 based optical flow deformation optimization framework to reduce the undesirable artifacts by employing sparse smoothness constraint, and then to match the appropriate points of each shape. Because some large differences between source and target shapes have non-zero deformation gradient values while other regions keep zero, which indicates the sparse property, we present a new measure via L 0 norm as smoothness regularization constraint. The sparse gradient counting measure can confine the number of deformation gradient changes, which allows large gradient magnitude by nature. In addition, to retrieve the motion optimally, this sparse property shows that different degrees of smoothness deformation priors should be introduced to different regions. As a consequence, our proposed L 0 based method can retrieve an optimal deformation to yield better motion than the L 2 based method [5] and avoid suffering from the error accumulation.
A preliminary version of this work appears in [7]. In this paper, we have made efforts to improve our paper in terms of the approach, theoretical analysis, and performance evaluation. We add more scenarios for evaluation and discuss the convergence, execution time, limitation of our proposed method. We demonstrate that our proposed L 0 based optical flow deformation method can robustly and automatically match the surfaces of the fluid simulation data produced by semi-Lagrangian advection method than the L 2 ones [5,6]. Our proposed method with the sparsity measure can outperform both the user guidance registration and the third data-point support in terms of accuracy with small details and moderate correspondences. Our proposed method can generate more realistic fluid surface with the optimal space-time deformation under the L 0 gradient smooth constraint than the L 2 one, and outperform the state-of-the-art methods in terms of both objective and subjective quality.
The rest of this paper is organized as follows. Section 2 provides a brief survey of related work. Section 3 introduces the background of the fluid flow and the optical flow solve. Section 4 presents the proposed fluid simulation method with the L 0 based optical flow deformation. In Section 5, experiments are conducted for our proposed method in comparison with the state-of-the-art methods. The conclusion is given in Section 6.

Related Work
In this section, we briefly introduce the related work about the fluid simulation method and the motion reconstruction method with sparse constraints.
Fluid simulation has a long history in computer graphics research. One of the influential works in this area is published in Foster and Metaxas [4], and the stable advection approach introduced by [3] has made a significant progress of fluid simulation for visual effects. The FLIP method [8] which combines the particle and the grid data, and the purely particle based SPH approach [9] have been widely used in movie productions. Numerous extensions and improvements have been proposed, e.g., to introduce small-scale turbulent details [10], or to address the problem emerged in reality from fluid-air and fluid-solid interactions with the surface tension force and the adhesion force [11].
Several works have used guide behaviors to make liquid surfaces follow a desired shape or motion [12,13]. Unlike this modification approaches on existing fluid simulation, the work-flow that concentrates on align pre-computed space-time surfaces and interpolates between them for intermediate motions is mainly based on non-rigid ICP algorithms [6], and optical flow solve [5].
The non-rigid ICP algorithm has been applied to reconstruct the temporally coherent dynamic objects from data captured by the multi-view acquisition system [14] or the monocular depth camera [15]. This registration technique has also been applied to the tracking of changing topology shapes such as liquid surfaces [16]. These method will run into the problem of local minima that is associated with non-rigid ICP. Raveendran et al. [6] roughly prevents the problem with user-guided correspondences.
Regular optical flow aims at recovering the motion of dynamic objects, especially for consecutive frames. Up to now, the Horn-Schunck algorithm [17] is extended to estimate high quality flow field. The connection between the optical flow and fluid flow are established for typical flow visualizations [18]. The optical flow techniques have also been used to acquire the fluid velocities based on physical priors [19]. Optical flow estimation commonly minimizes an energy function to obtain an optimal result, where the smoothness regularization is incorporated as a convex regularizer. The optical flow deformation approach [5] is essentially similar to the non-rigid iterative closest point algorithm that incorporates user guidance [6] in term of the registration process. Likewise, this optical flow deformation approach robustly recovers a match among the inputs based on local smoothness constraints and obtains a deformation field to align the surfaces. The L 2 based smoothness regularization [5] may cause the error propagation and the accumulated errors. In this paper, we utilize a specific sparse regularization as a constraint to improve the results.
In fact, sparse representation [20] has successfully applied in field of image processing [21], computer vision [22], and graphics [15]. Especially, the sparsity priors have been utilized in the optical flow estimation and motion reconstruction. For example, a new sparsity prior of the flow field in certain domains is introduced to estimate dense flow fields, and a stronger additional sparsity constraint on the flow gradients is incorporated into the model to cope the measurement noises [22]. While [23] decomposes the motion field into sparse and non-sparse components to restrain the over-sparsity. In the motion tracking [15], the L 0 based motion regularizer is incorporated into the non-rigid deformation to propose the L 0 -L 2 non-rigid motion tracking method to efficiency stop the error propagation.
As we know, the points on the same surface part mostly share similar motion transformations, thus, the gradient field of the deformation exhibits the sparse property. The straightforward solution to exploit the sparsity is integrated the L 0 norm regularization to the energy function. However, the optimization problem with the L 0 norm minimization is non-convex problem, which can not be easily solved. In image smoothing problems [21], an L 0 gradient minimization method is proposed to solve the similar issue, where an alternating optimization strategy with half-quadratic splitting is utilized. In this paper, we propose an L 0 based optical flow deformation model and present an effective approach by adopting a similar alternating optimization strategy with half quadratic splitting. In addition, we conduct experiments on the fluid droplets and dam break scenarios to evaluate our proposed method.

Background
For brevity, we will give a brief overview of fluid flow, and then introduce the previous method of optical flow solve [5]. Throughout this paper we will use symbol (e.g., u) to denote continuous value, and bold symbol (e.g., u) to denote the discrete counterpart. The symbol with subscript (e.g., φ r ) is usually used to represent the derivative of the corresponding variable. Thus, a subscript will not denote the component of a vector except for the special statement.
Fluid Flow. The fluid flow is governed by incompressible Navier-Stokes equations [24], and a set of partial different equations that are supposed to hold throughout the fluid ∂v ∂t where v is used in fluid for the velocity of the fluid, ∂v ∂t is the commonly used material derivative w.r.t. time t for advection, f is the external body forces(to make the fluid behave in some desired way), the letter ρ stands for the fluid density, P is the force due to pressure that the fluid exerts on anything, and γ is the "kinematic viscosity".
Following the description of [19], the observable concentration of a marker fluid is denoted with φ and the advection of a scalar field φ based on a velocity field v can be expressed as where s is a source term, including internal and external sources. Assume that the marker fluid and the ambient fluid have approximately the same mass density, the flow will perform a pure advection when the negligible diffusivity (γ = 0), incompressibility (∇ · v = 0), and s = 0, the standard form of advection is given by Equation (2) is equivalent to the brightness-constancy constraint of optical flow as introduced by Horn and Schunck algorithm [17].
Optical Flow solve. The work of [5] uses the L 2 norm based optical flow solve to perform a registration step. In this setting, two animations from different simulation setups are taken as the inputs: a source and a target. Each animation is a sequence of surfaces over time, where the surfaces can be concatenated as time slices of a 4D volume. The registration mainly transforms source to further fit it into target through deformations u.
A flow simulations can be represented by 4D signed-distance values φ, and the parameter r is an external control knob to modify the simulation from source to target. The regular optical flow usually recovers the motion of dynamic object from two images taken at different times. Thuerey [5] instead uses optical flow to recover the motion among surfaces of different simulations. The algorithm employs the brightness-constancy assumption to recover the flow with motion discontinuities. Please refer to [5] for further details.
Applying the optical flow alone is underdetermined because of the non-linearity, so that various additional constraints or priors are imposed to regularize the problem. The Tikhonov regularizer term restricts the deformation to be as small as possible. E smooth defines a L 2 regularizer which smoothes the consistent motions along spatial as well as temporal dimensions. Given the optical flow constraint and some common regularizers, the estimation of deformation field u is achieved by minimizing the following energy min where β S and β T are the weights directly controlling the importance of each energy term. The energy to be minimized in a least-squares sense is given by where the data derivative φ r corresponds to the time derivation φ t for optical flow applications. The deformation vector u of (4) is devoted to move a points to its corresponding position from the source to the target, and is calculated over four dimensions (m = 4). The spatial gradient ∇φ is computed with the target surface. The minimization of E total is performed by an hierarchical scheme. The hierarchical procedure is used to recursively re-sample and deform the simulation data on coarser resolutions. The first deformation u 1 calculated by optical flow deformation model transforms source inputs closer to target inputs. Then another solve is performed for the remaining difference, yielding the second deformation u 2 . The deformation is continually refined with residual iterations until the refinement does not yield gains in quality anymore. We exactly follow [5] to align multiple consecutive Eulerian deformations so that they can be merged into a single, smooth linear deformation from the source to the target.

Method
The 4D registrations and interpolation pipeline were processed automatically as illustrated in Figure 2. The optical flow solve computed a dense field of four-component deformation vectors u, which mapped points of input A onto input B. The residual optical flow was performed until the refinement did not yield gains in quality. Then the deformation sequences from each iteration were aligned into a single deformation field. It was observed that the source surface should ideally deform to end up on the target surface. Based on this goal, the projection step [16] was adopted to retrieve the fine details at the surface that prevented by regularization. It is notable that the projection only gave good results when a suitable overall match was obtained by optical flow step. After deformations were calculated, the interpolation step [5] was further utilized to produce intermediate surface animations. Given two input sequences ψ 1 , ψ 2 with deformations u 1→2 and u 2→1 , the linear interpolation could significantly generate the final quality.

L 0 Based Optical Flow Deformation
The goal of optical flow in this work is to compute deformations and to recover the discontinuous fluid motions. In the following, a sparsity constraint on the deformation gradients is incorporated into the optical flow estimation model to reduce the undesirable errors.
The L 2 based optical flow solve [5] for the certain inputs produced by semi-Lagrangian advection algorithm cannot robustly and accurately match surfaces in the nearest-neighbor fashion. Especially in the inputs with little similarities because of the surface flickering, the deformations will try to match pieces in proximity but not necessarily the whole target shape. We propose an L 0 based deformation regularizer over the existing optical flow solve to implicitly define the appropriate correspondences and to prevent motion changes on other regions.
Specifically, the deformation is denoted by u to transform the source surface closer into the target surface. We formulate the L 0 based deformation by minimizing the following energy function where a new smoothness regularizer E smooth is introduced to the general optical flow solve to replace the L 2 norm smoothness term. E smooth measures the deformation differences between all pairs of neighboring points, and is given by where S(u) is the L 0 norm for the deformation gradient ∇u. The gradient ∇u p = {∂ x u p , ∂ y u p , ∂ z u p , ∂ t u p } for each point p of the deformation field is calculated based on the two 4D SDFs between neighboring points in space and time. The gradient magnitude |∂u p | is defined as the sum of gradient magnitudes over four dimensions. Our smoothness constraint with L 0 norm regularization is used to control the number of non-zero gradients globally and present the sparsity of deformation gradients.
We estimated the optimized deformation by minimizing the energy function The E smooth term of L 0 norm regularizer in (5) is a non-convex problem. Due to the computationally intractability of the L 0 norm, the half-quadratic splitting strategy utilized in image smoothing achieves a high qualitative effect [21]. Inspired by the strategy described in [21], we split the objective function (8) into two sub-problems by introducing auxiliary variables into the energy function. We argue that this method found appropriate matches and retrieved the surface motions robustly.

Solver
We introduced an auxiliary variable v p , corresponding to ∇u p , and reformulated the energy function as and β A is an automatically adapting parameter to control the similarity between the gradients ∇u p and their corresponding variables v p . To solve this problem, we alternatively fix variables u p to solve v p and fix variables v p to solve u p in each iteration.

Minimizing v
The terms not involving v p in (9) are removed to estimate v p As ∇u p is pre-fixed, (10) has a close form solution

Minimizing u
As v p is pre-fixed, the u p estimation can be represented in a quadratic function (12) formulates a L 2 based minimization. The solution is similar to the method of L 2 based optical flow [5]. Setting the gradient of E(u) to be zero gives a closed solution for (12), which can be expressed as a linear equation The ψ r was computed with ψ 1 − ψ 2 , where the source and target inputs are denoted by ψ 1 and ψ 2 respectively. The spatial gradient ∇ψ is denoted with the target surface. L denotes the discrete Laplacian.
The optical flow based on L 0 regularizer is summarized in Algorithm 1. Initially, the algorithm was performed with two SDF inputs ψ and a zero deformation field u = 0. The function deform(ψ, u, α) applied the deformation u to input ψ with weight α ∈ [0...1]. A single, smooth linear deformation from source to target was required rather than deformation sequences produced by the iteration runtime. We followed [5] to align multiple deformations with the alignVelocity(u) function.
The reliable error metric [5] was adopted to quantitatively evaluate whether the current optical flow solve led to an improvement. More emphasis was put on the surface region rather than the volume itself. Given two inputs ψ 1 and ψ 2 (which we define them in a discrete setting), the error metric becomes where V is the volume of a cell. The indicator function h detecting misalignment surfaces is defined as where s 1 and s 2 represent two input SDFs ψ 1 and ψ 2 .

Algorithm 1 Optical flow based on L 0 regularizer
Input: inputs ψ 1 and ψ 2 , u, smoothing weight β S , parameters In our proposed method, we also utilized the commonly used hierarchical procedure to recursively re-sample and deform the data on coarser resolutions. The deformation was continually refined with residual iterations until the refinement did not yield gains in quality anymore, e.g., the ratio of the error of two consecutive iterations was greater than a given threshold for each resolution.

Experimental Results
We conducted the experiments on the fluid simulation of both the flying droplets and dam break scenarios. The standard fluid simulation, surface matching process and the interpolation step were run on the extensible fluid simulation framework, called mantaflow [25]. The fluid surface sequences of both the source and target inputs were generated by running the standard fluid simulation. When the deformations were calculated, new fluid animations could be produced in the interpolation weight range from 0 to 1. We utilized our proposed L 0 based fluid simulation method and the state-of-the-art method [5] to calculate the deformation and further generate the interpolation results as the outputs, respectively. The interpolation step yielded a smooth transition between source inputs and target inputs, and generated new outputs along a range of input positions. In the following, we compare our L 0 based refinement results with those of the state-of-the-art [5] to demonstrate the accuracy of our proposed method.

Convergence Analysis
We took the flying droplets as an example to study the convergence property of our proposed L 0 based optical flow solve method in Figure 3. The error metric error in [5] was computed between the target SDFs and the deformed one at different iterations. We observe that the errors descended with the increase of iterations so that our proposed algorithm was convergent. In order to speed up the convergence, we monotonically increased β A by multiplying the parameter κ = 2. In addition, Figure 3 shows the influence of the error curves with different parameters κ. We let the parameter κ varying from 1.2 to 2. From the convergence curve, we can observe that the best performance was generally achieved at κ = 2. Thus, in our experiments with a resolution of 64 4 , we set β A to be 1E − 3 in the first iteration and increase it by multiplying parameter κ = 2 after each iteration when error was reduced.

Deformation Results
We evaluate the quality of the deformation through the visual quality of the interpolated results, where the uni-directional nearest-neighbor interpolation is used, such as the deformation direction from ψ 1 to ψ 2 .
We consider two flying droplets with different initial shapes and different positions. They drop in the pool at different times to create different water ripples. We try to deform the source one to the target one by our proposed fluid simulation method and state-of-the-art method [5] respectively, with the interpolation weight as 1. Figure 4 shows the interpolated results of the L 2 based and the L 0 based deformation methods respectively. In order to easily observe, the sectional view images at a fixed frame are shown along with the interpolated results. We can observe that our results are closer to the target input ψ 2 and with less artifacts than those of the L 2 based method. In addition, the interpolated results of these two different methods with the deformation from the input ψ 2 to the input ψ 1 are shown in Figure 5. The L 2 based deformation method may create some ambiguous correspondences between the input ψ 1 to the input ψ 2 , and the undesired ripples appear at the surface due to the suboptimal matching. Our proposed method can alleviate the artificial motions brought by the L 2 based method effectively.
Three examples of interpolated results at different frames for flying droplets are demonstrated in Figure 6, where the droplets were falling into the pool and causing the varying splashes. Figure 7 demonstrates the interpolation results of the droplets falling into the pool and the varying splashes. We observe that our interpolated results could preserve the large-scale structures and recover the small-scale details. Our L 0 based method could reduce the artificial motions brought by the L 2 based method, which demonstrates that our L 0 based method could calculate moderate correspondences and a good deformation field and could further reduce the undesired errors of complex surfaces of the liquid.
We take the dam break as an example to demonstrate the accuracy of our optimization method. In this set of animations, a block of water was released at one side of a pool, and this formed a wave that swept the pool. As shown in Figure 8, we produced fluid animations for two different widths for this block of water. The interpolated results were generated with the interpolation weight as 1 to compare the generated surfaces with the target surfaces visually. We compared the results of our proposed L 0 based method with those of [5]. As shown in the right columns of Figure 8, the L 2 based method caused the undesired ripples at the surface, and the L 0 optimization method captured the accurate motions of the waves. Therefore, the surface of our results based on the L 0 regularization was closer to the target one than that based on the L 2 optimization. Figure 9 shows the whole range of the interpolation results of dam break deformed from ψ 1 to ψ 2 with different interpolation weights. The inputs had inconsistent behaviors in this static frame. Compared with the L 2 based regularizer which restricted the deformation locally across the surface region, the L 0 based regularizer effectively recovered the small-scale details. The animation simulation of all cases by using the L 2 based method [5] and our L 0 methods can be found in the supplemental video. The accompanying video shows a direct comparison to illustrate the difference between the results of [5] and our method.

Evaluation
The error metric mentioned in Section 4.2 could quantitatively and significantly detect the differences between the deformed surface and the target surface. We utilized the error metric (14) to calculate the errors between the target surfaces and the generated surfaces which were deformed with the L 2 based and the L 0 based deformation methods respectively. The error values can be found in Table 1, which includes two aforementioned scenarios (droplet, dam break). We show the results with the deformation direction from ψ 1 to ψ 2 and the inverse one. We can observe that our proposed L 0 based fluid simulation method could obtain less errors than the L 2 based method [5], and generate the desired correspondences. Our experiments used data-sets with a resolution of 64 × 64 × 64 × 64. We utilized the commonly used hierarchical procedure to recursively re-sample and deform the data on coarser resolutions, including the resolutions of 16 × 16 × 16 × 16, 32 × 32 × 32 × 32, and 64 × 64 × 64 × 64, denoted by 16 4 , 32 4 , and 64 4 , respectively. The changes of resolution were very useful to reduce the inaccuracies of the optical flow solve. Table 2 lists the error between the deformed surface and the target surface, the factor (the ratio of the error of this iteration to that of last iteration), and the time to calculate the deformation during each iteration. In our experiment, we stopped the iterations when the factors for the resolution of 16 4 , 32 4 , and 64 4 were greater than 70%, 95%, and 95%, respectively. Particularly for the resolution of 16 4 , the iterations stopping criterion should exceed three times. Note that the deformation was not refined if the gains in quality could not be yielded anymore. Thus, we did not use the deformation obtained by the last iteration. From the Table 2, we can see the execution time of the deformation by our L 0 based optical flow method with the half quadratic splitting approach was less than that of the L 2 based optical flow method.

Limitation
We take an example to illustrate in which case the interpolation method yielded undesirable results. Figure 10 shows the interpolation results of the water hitting different obstacles. The obstacle in the pool led to different liquid motions for the splash near the obstacle. The surface of both the L 0 and the L 2 based interpolation results demonstrate that part of the water was pushed into the obstacle. Even though the L 0 based result outperformed the L 2 based result, both them were undesirable results. Because the large differences of the source and the target surface had difficulty finding the corresponding points between the source and target surface, the inaccuracy deformation calculation led to incorrect matching and further resulted in undesirable results and artifacts. Through our experiments, we observed that either the L 0 or the L 2 based interpolation method failed when the fluid had relatively complex structures and large differences between the source and the target. Therefore, there is still room to improve the fluid simulation approach.

Conclusions
In this paper, we proposed a L 0 -norm based deformation for fluid simulation to obtain a good alignment of two inputs. The main contribution of our work is a L 0 -norm optimization approach to align fluid simulation with changing topologies. The L 0 based optical flow takes advantage of the sparsity property of the fluid motion of neighboring frames to constrain the deformation. Moreover, we propose a novel effective approach to solve the L 0 -based optimization problem by making use of L 0 gradient minimization, which can globally control the non-zero gradients to approximate the deformation in a sparsity-control manner. Experimental results demonstrate that our method outperforms the state-of-the-art L 2 based optical flow [5] and accurately recovers motions between the space-time surfaces of different simulations.
There are two different applications of the interpolation method. The time step of the input data influences the continuity of the animations visually. We can make use of our proposed fluid simulation scheme to produce a highly continuous behavior for arbitrary neighboring sequences with a similar topology. Alternatively, we can utilize different scenarios to interpolate new fluid animations frame-by-frame. Our proposed approach is a data-driven optical flow solve, however, the underlying physics is ignored. In the future, we will incorporate with other physical optimization of the surface motions, such as the constraints for mass conservation.