DEFORMATION-AWARE NEURAL NETWORKS
Year: 2019

Abstract:
We propose a novel approach for deformation-aware neural networks that learn the weighting and synthesis of dense volumetric deformation fields. Our method specifically targets the space-time representation of physical surfaces from liquid simulations. Liquids exhibit highly complex, non-linear behavior under changing simulation conditions such as different initial conditions. Our algorithm captures these complex phenomena in two stages: a first neural network computes a weighting function for a set of pre-computed deformations, while a second network directly generates a deformation field for refining the surface. Key for successful training runs in this setting is a suitable loss function that encodes the effect of the deformations, and a robust calculation of the corresponding gradients. To demonstrate the effectiveness of our approach, we showcase our method with several complex examples of flowing liquids with topology changes. Our representation makes it possible to rapidly generate the desired implicit surfaces. We have implemented a mobile application to demonstrate that real-time interactions with complex liquid effects are possible with our approach.

heading: 1 INTRODUCTION
text: Learning physical functions is an area of growing interest within the research community, with applications ranging from physical priors for computer vision problems Kyriazis & Argyros (2013), over robotic control Schenck & Fox (2017), to fast approximations for numerical solvers Tompson et al. (2017). While the underlying model equations for many physics problems are known, finding solutions is often prohibitively expensive for phenomena on human scales. At the same time, the availability of model equations allows for the creation of reliable ground truth data for training, if enough computational resources can be allocated.
Water, and liquids in general, are ubiquitous in our world. At the same time, they represent an especially tough class of physics problems, as the constantly changing boundary conditions at the liquid-gas interface result in a complex space of surface motions and configurations. In this work we present a novel approach to capture parametrized spaces of liquid behavior that is based on space-time deformations. We represent a single 3D input surface over time as a four-dimensional signed-distance function (SDF), which we deform in both space and time with learned deformations to recover the desired physical behavior. To calculate and represent these deformations efficiently, we take a two-stage approach: First, we span the sides of the original parameter region with precomputed deformations, and infer a suitable weighting function. In a second step, we synthesize a dense deformation field for refinement. As both the parameter weighting problem and the deformation synthesis are highly non-linear problems, we demonstrate that neural networks are a particularly suitable solver to robustly find solutions.
We will demonstrate that it is possible to incorporate the non-linear effects of weighted deformations into the loss functions of neural networks. In particular, we put emphasis on incorporating the influence of deformation alignment into the loss gradients. This alignment step is necessary to ensure the correct application of multiple consecutive deformations fields. The second stage of our
algorithm is a generative model for deformation fields, for which we rely on a known parametrization of the inputs. Thus, in contrast to other generative models which learn to represent unknown parametrization of data sets Radford et al. (2016), our models are trained with a known range and dimensionality to parameter range, which serves as input.
Once trained, the models can be evaluated very efficiently to synthesize new implicit surface configurations. To demonstrate its performance, we have implemented a proof-of-concept version for mobile devices, and a demo app is available for Android devices in the Google Play store. Our approach generates liquid animations several orders of magnitude faster than a traditional simulator, and achieves effective speed up factors of more than 2000, as we will outline in Sec. 5. The central contributions of our work are:
• A novel deformation-aware neural network approach to very efficiently represent large collections of space-time surfaces with complex behavior.
• We show how to compute suitable loss gradient approximations for the sub-problems of parameter and deformation inference.
• In addition we showcase the high performance of our approach with a mobile device implementation that generates liquid simulations interactively.

heading: 2 RELATED WORK
text: Capturing physical behavior with learning has a long history in the field of learning. Early examples targeted minimization problems to determine physical trajectories or contact forces Bhat et al. (2002); Brubaker et al. (2009), or plausible physics for interacting objects Kyriazis & Argyros (2013; 2014). Since initial experiments with physically-based animation and neural networks Grzeszczuk et al. (1998), a variety of new deep learning based works have been proposed to learn physical models from videos Battaglia et al. (2016); Chang et al. (2016); Watters et al. (2017). Others have targeted this goal in specific settings such as robotic interactions Finn et al. (2016), sliding and colliding objects Wu et al. (2015; 2016), billiard Fragkiadaki et al. (2015), or trajectories in height fields Ehrhardt et al. (2017). The prediction of forces to infer image-space motions has likewise been targeted Mottaghi et al. (2016a;b), while other researchers demonstrated that the stability of stacked objects can be determined from single images Lerer et al. (2016); Li et al. (2016). In addition, the unsupervised inference of physical functions poses interesting problems Stewart & Ermon (2017). While these methods have achieved impressive results, an important difference to our method is that we omit the projection to images. I.e., we directly work with three-dimensional data sets over time.
In the context of robotic control, physics play a particularly important role, and learning object properties from poking motions Agrawal et al. (2016), or interactions with liquids Schenck & Fox (2017) were targeted in previous work. Learning physical principles was also demonstrated for automated experiments with reinforcement learning Denil et al. (2016). Recently, first works have also addressed replacing numerical solvers with trained models for more generic PDEs Farimani et al. (2017); Long et al. (2017). In our work we target a more narrow case: that of surfaces deforming based on physical principles. However, thanks to this narrow scope and our specialized loss functions we can generate very complex physical effects in 3D plus time.
Our method can be seen as a generative approach representing samples from a chaotic process. In this context, the learned latent-space representation of regular generative models Masci et al. (2011); Rezende et al. (2014); Goodfellow et al. (2014); Radford et al. (2016); Isola et al. (2017) is replaced by the chosen parametrization of a numerical solver. Our model shares the goal to learn flow physics based on examples with other methods Ladicky et al. (2015); Tompson et al. (2017); Chu & Thuerey (2017), but in contrast to these we focus on 4D volumes of physics data, instead of localized windows at single instances of time. Alternative methods have been proposed to work with 4D simulation data Raveendran et al. (2014); Thuerey (2017), however, without being able to capture larger spaces of physical behavior. Due to its focus on deformations, our work also shares similarities with methods for optical flow inference and image correspondences Bailer et al. (2016); Dosovitskiy et al. (2015); Ranjan & Black (2016); Ilg et al. (2016). A difference to these approaches is that we learn deformations and their weighting in an unsupervised manner, without explicit ground truth data. Thus, our method shares similarities with spatial transformer networks (STNs) Jaderberg et al. (2015), and unsupervised approaches for optical flow Meister et al. (2017).
However, instead of aligning two data sets, our method aims for representing larger, parametrized spaces of deformations.

heading: 3 LEARNING DEFORMATIONS
text: We first explain our formulation of the simulation parameter space, which replaces the latent space of other generative models such as autoencoders Masci et al. (2011), GANs Radford et al. (2016), or auto-regressive models Van Oord et al. (2016). Given a Navier-Stokes boundary value problem with a liquid-gas interface, we treat the interface over time as a single space-time surface. We work with a set of these space-time surfaces, defined over a chosen region of the N-dimensional simulation parameters α. We assume the parameters to be normalized, i.e. α ∈ [0, 1]N . In practice, α could contain any set of parameters of the simulation, e.g. initial positions, or even physical parameters such as viscosity. We choose implicit functions φ(α) ∈ R4 → R to represent specific instances, such that Γ(α) = { x ∈ R4;φ(α,x) = 0 } is the space of surfaces parametrized by α that our generative model should capture. In the following, φ and ψ will denote four-dimensional signed distance functions. We will typically abbreviate φ(α) with φα to indicate that φα represents a set of constant reference surfaces. While we will later on discretize all functions and operators on regular Cartesian grids, we will first show the continuous formulation in the following.
Deforming Implicit Surfaces Representing the whole set φα is very challenging. Due to bifurcations and discretization artifacts φα represents a process that exhibits noisy and chaotic behavior, an example for which can be found in Fig. 1. Despite the complex changes in the data, our goal is to find a manageable and reduced representation. Thus, we generate an approximation of φα by deforming a single initial surface in space and time. We apply space-time deformations u : R4 → R4 to the initial surface ψ0(x). As a single deformation is limited in terms of its representative power to produce different shapes, we make use of N sequential, pre-computed deformations, thus ui with i ∈ [1 · · ·N ], each of which is scaled by a scalar weight parameter βi, whose values are normally between 0 and 1. This gives us the freedom to choose how much of ui to apply. The initial surface deformed by, e.g., u1 is given by ψ0(x− β1u1).
The sequence of deformed surfaces by sequentially applying all pre-computed ui is given by ψi(x,β) = ψi−1(x − βiui). It is crucial to align such sequences of Eulerian deformations with each other. Here we employ the alignment from previous work Thuerey (2017), which we briefly summarize in the following, as it influences the gradient calculation below. Each deformation i relies on a certain spatial configuration for the input surface from deformation i− 1. Thus, when applying ui−1 with a weight βi−1 < 1, we have to align ui correspondingly. Given the combined deformation vsum(x,α) = ∑N i=1 βiu ∗ i (x), with intermediate deformation fields u ∗ i−1(x) = ui−1(x − u∗i (x)),
we compute an inversely weighted offset field as vinv(x,α) = − ∑N i=1(1 − βi)u∗i (x). This offset field is used to align the accumulated deformations to compute the final deformation as vfin(x + vinv(x,β),β) = vsum(x,β). vfin now represents all weighted deformations βiui merged into a single vector field. Intuitively, this process moves all deformation vectors to the right location for the initial surface ψ0, such that they can be weighted and accumulated.
To achieve the goal of representing the full set of target implicit surfaces φα our goal is to compute two functions: the first one aims for an optimal weighting for each of the deformations in the sequence, i.e. β, while the second function computes a final refinement deformation w after the weighted sequence has been applied. We will employ two neural networks to approximate the two functions, which we will denote as fp, and fd below. Both functions depend only on the simulation parameters space α, i.e., fp(α) = β, and fd(α) = w.
Splitting the problem into fp and fd is important, as each of the pre-computed deformations weighted by fp only covers a single trajectory in the space of deformed surfaces. In the follow-
ing, we employ an optical flow solve from previous work to pre-compute deformations between the inputs Thuerey (2017), which deform the input for the extrema of each original parameter dimension αi. E.g., for a two-dimensional parameter space this yields two deformations along the sides of the unit cube. This first step robustly covers rough, large scale deformations. As a second step, we employ fd, which is realized as a generative CNN, to infer a deformation for refining the solution. Below we will explain the resulting equations for training. The full equations for applying the deformations, and a full derivation of our loss functions can be found in the supplemental materials. To shorten the notation, we introduce the helper function D(xi,α), which yields a deformed set of coordinates in R4 depending on α that incorporates the deformation sequence weighted by fp(α), and a refinement deformation from fd(α).
We express the overall goal in terms of minimizing the L2 distance between the deformed and the target implicit surfaces for all possible values in the parameter space α, using β and w as degrees of freedom:
argmin β,w L,L =
∫ ‖ψ0(D(xi,α))− φα‖22 dα . (1)
Our work addresses the problem of how to compute weighting of the deformations and on synthesizing the refinement field. The main difficulty lies in the non-linearity of the deformations, which is why we propose a novel method to robustly approximate both functions with NNs: fp will be represented by the parameter network to compute β, and we make use of a deformation network that to generate w. We employ relatively simple neural networks for both functions. Key for training them is encoding the effect of deformations in the loss functions to make the training process aware of their influence. Hence, we will focus on describing the loss functions for both networks and the corresponding discretized gradients in the following.
Learning Deformation Weights For training the NNs we propose the following objective function, which measures the similarity of a known reference surface φα and the corresponding, approximated result ψ0(x,α) for a parameter value α. We introduce the numerical equivalent of the continuous L2 loss from Eq. (1) as
L = 1
2 ∑ i (ψ0(D(xi,α))− φα(xi))2 ∆xi , (2)
which approximates the spatial integral via the sum over all sample points i with corresponding evaluation positions xi ∈ R4, where ∆xi = ‖xi − xi−1‖ is constant in our case. This corresponds to a regular grid structure, where the loss is accumulated per cell. The central challenge here is to compute reliable gradients for D, which encapsulates a series of highly non-linear deformation steps. We first focus on inferring β, with w = 0.
The gradient of Eq. (2) with respect to one component of the deformation parameters βj is then given by
d dβj L = ∑ i ( − ∂ ∂βj vfin(xi + vinv(xi,β),β) · ∇ψ0(xi − vfin(xi,β)) ) (ψ0(xi,β)− φα(xi)) ,
(3) where the first term on the sum over i in parentheses represents the gradient of the deformed initial surface. Here we compute the derivative of the full deformation as ∂∂βj vfin(xi + vinv(xi,β),β) = u∗i (xi). The offset by vinv on the left hand side indicates that we perform a forward-advection step for this term. Details of this step, and for the full derivation of the gradient are given in Appendix B.1.
After parameter network + Deformation network Reference surface
w
↵
A trained NN with this loss functions yields an instance of fp, with which we can infer adjusted deformation weights fp(α) = β.
Learning to Generate Deformations Based on β, we apply the deformation sequence ui. The goal of our second network, the deformation network fd, is to compute the refinement deformation w. In contrast to the pre-computed ui, fd(α) = w now directly depends on α, and can thus capture the interior of the parameter space. Given the initial surface ψ0 deformed by the set of βiui, which we will denote as ψ̃ below, the refinement deformation is applied with a final deformation step as ψ(x) = ψ̃ (x−w(x, α)). In order to compute the gradient of the deformation loss, we introduce the indicator function χj(x) for a single deformation vector wj of w. We found it useful to use a fine discretization for the implicit surfaces, such as ψ, and lower resolutions for w. Hence, each discrete entry wj can act on multiple cells of ψ, which we enumerate with the help of χj . Now the derivative of Eq. (1) for a fixed β with respect to a single deformation vector wj of w is given by
d dwj L = − ∑ i χj(xi)∇ψ̃(xi −w(xi,α)) ( ψ̃(xi,α)− φα(xi) ) . (4)
The full derivation of this gradient is given in Appendix B.2. Our approach for deformation learning can be regarded as an extension of STNs Jaderberg et al. (2015) for dense, weighted fields, and semi-Lagrangian advection methods. The parameter network corresponds to an STN which learns to combine and weight known deformation fields. The deformation network, on the other hand, resembles the thin plate spline STNs, where a network generates an offset for each cell center, which is then used to sample a deformed image or grid. Note that in our case, this sampling process corresponds to the semi-Lagrangian advection of a fluid simulation.
Training Details For fp we use a simple structure with two fully connected layers, while fd likewise contains two fully connected layers, followed by two or more four-dimensonal de-convolution layers. All layers use ReLU activation functions. Details can be found in App. B, Fig. 12.
In practice, we also found that a small amount of weight decay andL2 regularization of the generated deformations can help to ensure smoothness. Thus, the loss function of the deformation network, with regularization parameters γ1 and γ2 is
Lt = L+ γ1||θ||2 + γ2||w||2 , (5) where θ denotes the network weights. In addition, regular SDFs can lead to overly large loss values far away from the surface due to linearly increasing distances. Thus, we apply the tanh() function to the SDF values, in order to put more emphasis on the surface region.
Special care is required for boundary conditions, i.e, the sides of our domain. Assuming constant values outside of the discretized domain, i.e. ∂ψ(x)/∂n = 0 for all x at the domain sides leads to vanishing gradients ∇ψ(x) = 0 in App. B, Eq. (4). We found this causes artificial minima and maxima in the loss function impeding the training process. Hence, we extrapolate the SDF values with ∂ψ(x)/∂n = ±1 in order to retrieve non zero gradients at the domain sides. To train both networks we use stochastic gradient descent with an ADAM optimizer and a learning rate of 10−3. Training is performed separately for both networks, with typically 1000 steps for fd,
1
Initial surface Reference surface+ Parameter network
0
↵
and another ca. 9000 steps for fd. Full parameters can be found in App. B, Table 2. As training data we generate sets of implicit surfaces from liquid simulations with the FLIP method Bridson (2015). For our 2D inputs, we use single time steps, while our 4D data concatenates 3D surfaces over time to assemble a space-time surface. Working in conjunction, our two networks capture highly complex behavior of the fluid space-time surface φα over the whole parameter domain. We will evaluate the influence of the networks in isolation and combination in the following.

heading: 4 EVALUATION
text: In order to evaluate our method, we first use a two-dimensional parameter space with two dimensional implicit surfaces from a liquid simulation. An overview of the space of 2156 training samples of size 1002 can be found in the supplemental materials. For our training and synthesis runs, we typically downsample the SDF data, and use a correspondingly smaller resolution for the output of the deformation network, see Appendix C.2, Table 2. The effect of our trained networks in terms of loss reductions is shown on the left side of Fig. 4 under Liquid 2D. As baseline we show the loss for the undeformed surface w.r.t. the test data samples. For this 2D data set, employing the trained parameter network reduces the loss to 59.4% of the initial value. Fig. 3 shows the surface of an exemplary result. Although our result does not exactly match the target due to the constraints of the pre-computed deformations, the learned deformation weights lead to a clear improvement in terms of approximating the target.
The inferred deformation of our deformation network further reduces the surface loss to 26.6% of its initial value, as shown in Fig. 4. This is equivalent to a reduction to 44.8% compared the result after applying the weighted deformations. Note that the graphs in this figure correspond to an ablation study: starting with the loss for the undeformed surface, over using only the parameter network, to deformations for a flat surface, to our full algorithm. An example surface for these two-dimensional cases can be found in Fig. 5. This figure compares the surface after applying weighted and inferred deformations, i.e. our full method (right), with a surface deformed by only by deformations weighted by the parameter network (left). The NN deformation manages to reconstruct the two arm in the center of the surface, which the pre-computed deformations fail to capture. It is also apparent
drop_res070000_025000_066666: frames 2,10,25
that despite the improvement, this surface does not reach the tip of the arms. This is caused by regularization over the varying set of target surfaces, leading to an averaged solution for the deformation. Additional examples for this two dimensional setup can be found in the supplemental video.
4D Surface Data Next we consider complete space-time data sets in four dimensions. with a three dimensional parameter space α. The three parameter dimensions are x- and y-coordinates of the initial drop position, as well as its size. We use a total of 1764 reference SDFs with an initial resolution of 1004, which are down-sampled to a resolution of 404. To illustrate the capabilities of the deformation network, we start with a completely flat initial surface as ψ0, and train the deformation network to recover the targets. As no pre-computed deformations are used for this case, we do not train a parameter network. The flat initial surface represents an especially tough case, as the network can not rely on any small scale details in the reference to match with the features of the targets. Despite this difficulty, the surface loss is reduced to 30% of the initial loss purely based on the deformation network. A set of visual examples can be seen in Fig. 6. Due to the reduced resolution of the inferred deformation w.r.t. the SDF surface, not all small scale features of the targets are matched. However, the NN manages to reconstruct impact location and size very well across the full parameter space. Additional 2D and 4D results can be found in the supplemental materials.
Once we introduce a regular initial surface for ψ0, in this case the zero parameter configuration with the drop in one corner of the domain with the largest size, our networks perform even better than for the 2D data discussed above. The weighted deformations lead to a loss reduction to 55.6% of the initial value, and the learned deformation reduce the loss further to 22.2% of the baseline loss (Fig. 4). An example is shown in Fig. 7. In contrast to the flat surface test, the network deformation can now shift and warp parts of ψ0, such as the rim of the splash of Fig. 7 to match the targets.

heading: 5 ADDITIONAL RESULTS WITH IMPLICIT SURFACES IN 4D
text: Our method yields highly reduced representations which can be used to very efficiently synthesize simulation results. To demonstrate the representational capabilities and the performance of our method, we have integrated the evaluation pipeline for our trained networks into an Android application. As our method yields an implicit surface as 4D array, visualizing the resulting animation is very efficient. We render slices of 3D data as implicit surfaces, and the availability of a full
3D representations makes it possible to add curvature-based shading and secondary particle effects on the fly. In this context, please also consider the supplemental materials, which contain these sequences in motion. They are available at: https://ge.in.tum.de/publications/ 2017-prantl-defonn/.
Performance One of the setup available in this app is the liquid drop setup with 3D parameter space described above. With this setup a user can release drops of liquid at varying positions and of varying size. An example reference and generated result can be found in Fig. 8. For this liquid drop setup, evaluating the network takes 69ms on average, and assembling the final deformation field another 21.5ms. We use double buffering, and hence both tasks are evaluated in a background thread. Rendering the implicit surface takes an average of 21ms, with an average frame rate of 50 fps. The original simulation for the drop setup of Fig. 8 took 530 seconds on average with a parallel implementation to generate a single 4D surface data point.Assuming a best-case slowdown of only 4x for the mobile device, it would require more than 32 minutes to run the original simulation there. Our app generates and renders a full liquid animation in less than one second in total. Thus, our algorithm generates the result roughly 2000 times faster than the regular simulation. Our approach also represents the space of more than 1700 input simulations, i.e., more than 17GB, with less than 30MB of storage.
Stairs A next setup, shown in Fig. 9, captures a continuous flow around a set of obstacles. Liquid is generated in one corner of the simulation domain, and then flows in a U-shaped path around a wall, down several steps. In the interactive visualization, green arrows highlight in- and outflow regions. The three dimensional parametrization of this setup captures a range of positions for the wall and two steps leading to very different flow behaviors for the liquid. In this case the data set consists of 1331 SDFs, and our app uses an output resolution of 504. The corresponding loss measurements can be found in the right graphs of Fig. 4. As with the two previously discussed data sets, our approach leads to very significant reductions of the surface loss across the full parameter space, with a final residual loss of 31.3% after applying the learned deformation. Due to larger size of the implicit surfaces and the inferred deformation field, the performance reduces to a frame rate of 30 fps on average, which, however, still allows for responsive user interactions.
Discussion Our approach in its current form has several limitations that are worth mentioning. E.g., we assume that the space of target surfaces have a certain degree of similarity, such that a single surface can be selected as initial surface ψ0. In addition, our method currently does not make use of the fact that the inputs are generated by a physical process. E.g., it would be highly interesting for future work to incorporate additional constraints such as conservation laws, as currently our results can deviate from an exact conservation of physical properties. E.g., due to its approximating nature our method can lead to parts of the volume disappearing and appearing over time. Additionally, the L2 based loss can lead to rather smooth results, here approaches such as GANs could potentially improve the results.

heading: 6 CONCLUSIONS
text: We have presented a novel method to generate space-time surfaces with deformation-aware neural networks. In particular, we have demonstrated the successful inference of weighting sequences of aligned deformations, and the generation of dense deformation fields across a range of varied inputs. Our method exhibits significant improvements in terms surface reconstruction accuracy across the full parameter range. In this way, our networks can capture spaces of complex surface behavior, and
allow for real-time interactions with physics effects that are orders of magnitudes slower to compute with traditional solvers.
Beyond liquid surfaces, our deformation networks could also find application for other types of surface data, such as those from object collections or potentially also moving characters. Likewise, it could be interesting to extend our method in order to infer deformations for input sets without an existing parametrization.

heading: A DEFORMATION ALIGNMENT
text: As before, φα denotes the reference signed distance functions (SDFs) of our input parameter space, while ψ denotes instances of a single input surface, typically deformed by our algorithm. We will denote the single initial surface without any deformations applied with ψ0. Here, we typically use the zero point of our parameter space, i.e., ψ0 = φα0 , with α0 = 0. Hence, we aim for deforming ψ0 such that it matches all instances of φα as closely as possible.
For the pre-computed, end-point deformations, it is our goal to only use a single deformation for each dimension of the parameter space α. Thus u1 will correspond to α1 and be weighted by β1, and we can apply β1u1 to compute a deformation for an intermediate point along this dimension. Given the sequence of pre-computed deformations {u1,u2, . . . ,uN} and a point in parameter space {β1, . . . , βN} a straight-forward approach is to apply each deformation sequentially
ψ1(x,β) = ψ0(x− β1u1) ψ2(x,β) = ψ1(x− β2u2)
... ψN (x,β) = ψN−1(x− βNuN ). (6)
However, there are two disadvantages to this approach. The main problem is that the deformations ui are only meaningful if applied with βi = 1.
Thus, if a previous deformation wasn’t applied fully with a weight of 1, each subsequent deformation will lead to an accumulation of deformation errors. The second disadvantage of this simple method
is the fact that many advection steps have to be applied in order to arrive at the final deformed SDF. This also affects performance as each advection step introduces additional computations, and scattered memory accesses. This is best illustrated with the an example, shown in Fig. 10. In row 1) the first two figures in red show two pre-computed end-point deformations (red arrows). The first one (u1) moves a drop to the right, while u2 changes its size once it has reached its final position. Images with deformations show the source in dark blue, and the deformed surface in light blue. In this example, the two deformations should be combined such that the horizontal position and the drop size can be independently controlled by changing β1 and β2. E.g., on the top right, a correct solution for β1 = 0.5, β2 = 1 is shown. Row 2) of Fig. 10 shows how these deformations are applied in previous work: The second deformation acts on wrong parts of the surface, as the drop has not reached its left-most position for β1 = 0.5. The undesirable result in this case is a partially deformed drop, shown again in the middle row on the right.
We present an alternative approach which aligns the deformation fields to the final position of the deformation sequence. Then, all aligned deformation fields can simply be accumulated by addition, and applied to the input in a single step. To do this, we introduce the intermediate deformation fields:
u∗N (x) = uN (x),
u∗N−1(x) = uN−1(x− u∗N (x)), u∗N−2(x) = uN−2(x− u∗N (x)− u∗N−1(x)),
... u∗1(x) = u1(x− u∗N (x)− u∗N−1(x) . . .− u∗2(x)). (7)
Each u∗i is moved by all subsequent deformations u ∗ j , j ∈ [i + 1 · · ·N ], such that it acts on the correct target position under the assumption that βi = 1 (we will address the case of βi 6= 1 below). The Eulerian representation we are using means that advection steps look backward to gather data, which in our context means that we start with the last deformation uN to align previous deformations. Using the aligned deformation fields u∗ we can include β and assemble the weighted intermediate fields
vsum(x,β) = N∑ i=1 βiu ∗ i (x) (8)
and an inversely weighted correction field
vinv(x,β) = − N∑ i=1 (1− βi)u∗i (x). (9)
The first deformation field vsum represents the weighted sum of all aligned deformations, weighted with the correct amount of deformation specified by the deformation weights βi. The second deformation vinv intuitively represents the offset of the deformation field vsum from its destination caused by the β weights. Therefore, we correct the position of vsum by this offset with the help of an additional forward-advection step calculated as:
vfin(x + vinv(x,β),β) = vsum(x,β), (10)
This gives us the final deformation field vfin(x,β). It is important to note that the deformation vsum for a position x is assembled at a location x′ that is not known a-priori. It has to be transported to x with the help of vinv, as illustrated in Fig. 11. This correction is not a regular advection step, as the deformation is being ’pushed’ from x + vinv(x,β) to x. In order to solve this advection equation we use an inverse semi-Lagrangian step, inspired by algorithms such as the one by Lentine et al. Lentine et al. (2011), pushing values forward with linear interpolation. As multiple values can end up in a single location, we normalize their contribution. Afterwards, we perform several iterations of a ”fill-in” step to make sure all cells in the target deformation grid receive a contribution (we simply extend and average deformation values from all initialized cells into uninitialized regions).
The deformed SDF is then calculated with a regular advection step applying the final, aligned deformation with
ψ(x,β) = ψ0(x− vfin(x,β)). (11)
Based on our correction step from Eq. (10) this method now respects the case of partially applied deformations. As the deformations u∗(x) are already aligned and don’t depend on β, we can precompute them. To retrieve the final result it is now sufficient to sum up all deformations in vsum and vinv, then apply one forward-advection step to compute vfin, and finally deform the input SDF by applying semi-Lagrangian advection. While our method is identical with alignment from previous work Thuerey (2017) for βi = 1, it is important for practical deformations with weights βi 6= 1. Our method is illustrated in row 3) of Fig. 10. In the bottom left we see the deformation field vsum from previous work. It is also the starting point of our improved alignment, but never applied directly. Our method corrects vsum by transforming it into vfin, bottom center, which acts on the correct spatial locations. In this example, it means that the expanding velocities from u2 are shifted left to correctly expand the drop based on its initial position. Our method successfully computes the intended result, as shown in the bottom right image.
This algorithm for aligning deformations will be our starting point for learning the weights β. After applying the weighted deformations, we adjust the resulting surface we an additional deformation field generated by a trained model. In the following, we will derive gradients for learning the weighting as well as the refinement deformation.

heading: B LEARNING DEFORMATIONS
text: As outlined in the main document, we aim for minimizing the L2 loss between the final deformed surface and the set of reference surfaces φα, i.e.:
L = 1
2 ∑ i (ψ0(D(xi,α))− φα(xi))2 , (12)
where D(xi,α) denotes the joint application of all weighted and generated deformation fields.
B.1 LEARNING DEFORMATION WEIGHTING
We will first focus on the parameter network to infer the weighting of the pre-computed deformation fields based on the input parameters α. Thus, the NN has the goal to compute β(α) ∈ RN = (β1(α), . . . , βN (α)) in order to minimize Eq. (12). The application of the deformations weighted by β includes our alignment step from Sec. A, and hence the neural networks needs to be aware of its influence. To train the parameter network, we need to specify gradients of Eq. (12) with respect to the network weights θi. With the chain rule we obtain ddθlij L = dβ dθlij dL dβ . Since the derivative of the network output βi with respect to a specific network weight θlij is easily calculated with backpropagation Bishop (2006), it is sufficient for us to specify the second term. The gradient of Eq. (12) with respect to the deformation parameter βi is given by
d dβi L = ∑ j d dβi ψ(xj ,β) [ψ(xj ,β)− φα(xj)] , (13)
where we have inserted Eq. (11). While the second term in the sum is easily computed, we need to calculate the first term by differentiating Eq. (11) with respect to βi, which yields
d dβi ψ(x,β) = − d dβi vfin(x,β) · ∇ψ0(x− vfin(x,β)). (14)
As the gradient of ψ0 is straight forward to compute, ddβivfin(x,β) is crucial in order to compute a reliable derivative. It is important to note that even for the case of small corrections vinv(x,β), Eq. (10) cannot be handled as another backward-advection step such as vfin(x,β) = vsum(x − vinv(x,β),β). While it might be tempting to assume that differentiating this advection equation will produce reasonable outcomes, it can lead to noticeable errors in the gradient. These in turn quickly lead to diverging results in the learning process, due to the non-linearity of the problem.
The correct way of deriving the change in vfin(x,β) is by taking the total derivative of vsum(x,β) = vfin(x + vinv(x,β),β) with respect to βi:
d dβi vsum(x,β)
= ∂
∂βi vfin(x + vinv(x,β),β) + JV (x + vinv(x,β),β)
∂
∂βi vinv(x,β), (15)
where, JV (x + vinv(x,β),β) denotes the Jacobian of vfin with respect to x, evaluated at x + vinv(x,β). Rearranging Eq. (15) and inserting vsum and vinv yields
∂
∂βi vfin(x + vinv(x,β),β) (16)
= d
dβi vsum(x,β)− JV (x + vinv(x,β),β)
∂
∂βi vinv(x,β)
= d
dβi N∑ i=1 βiu ∗ i (x) + JV (x + vinv(x,β),β) ∂ ∂βi N∑ i=1 (1− βi)u∗i (x)
= [1− JV (x + vinv(x,β),β)]u∗i (x). (17)
We note that the Jacobian in the equation above has small entries due to the smooth nature of the deformations vfin. Thus, compared to the unit matrix it is small in magnitude. Note that this relationship is not yet visible in Eq. (15). We have verified in experiments that JV does not improve the gradient significantly, and we thus set this Jacobian to zero, arriving at
∂
∂βi vfin(x + vinv(x,β),β) ≈ u∗i (x), (18)
where the u∗ are the deformation fields aligned for the target configuration from Eq. (7). We use Eq. (18) to estimate the change in the final deformation fields for changes of the i-th deformation parameter. We see that this equation has the same structure as Eq. (10). On the left-hand side, we have ∂∂βivfin, evaluated at x + vinv(x,β), whereas u ∗ i on the right-hand side is evaluated at x. To calculate ddβivfin(x,β) then, we can use the same forward-advection algorithm, which is applied to the correction in Eq. (10). With this, we have all the necessary components to assemble the gradient from Eq. (13) for training the parameter network with back-propagation.
B.2 LEARNING TO GENERATE DEFORMATIONS
Our efforts so far have been centered around producing a good approximation of φα, with a set of given end-point deformations {u0, . . . ,un}. The performance of this method is therefore inherently constrained by the amount of variation we can produce with the deformation inputs. To allow for more variation, we propose to generate an additional space-time deformation field w(α), that changes with the simulation parameters α. Once again, we model this function with a neural network, effectively giving the network more expressive capabilities to directly influence the final deformed surface.
For this network we choose a structure with a set of four-dimensional deconvolution layers that generate a dense space-time deformation field. We apply the trained deformation with an additional
advection step after applying the deformations weighted by the parameter network:
ψ̃(x) = ψ0 (x− vfin(x,β(α))) , (19) ψ(x) = ψ̃ (x−w(x,α)) . (20)
Thus, the deformation network only has to learn to refine the surface ψ̃ after applying the weighted deformations, in order to accommodate the nonlinear behavior of φα.
As input, we supply the deformation network with the simulation parameters α = (αi, . . . , αN ) as a vector. The output of the network are four-component vectors, with the resolutionRx×Ry×Rz×Rt. Note that in general the SDF resolution and the deformation resolution do not need to be identical. Given a fixed SDF resolution, we can use a smaller resolution for the deformation, which reduces the number of weights and computations required for training. Thus in practice, each four-dimensional vector of the deformation acts on a region of the SDF, for which we assume the deformation to be constant. Therefore, we write the deformation field as
w(x,α) = ∑ j ξj(x) wj(α), (21)
where ξj(x) is the indicator function of the j-th region on which the four-dimensional deformation vector wj(α) acts. This vector is the j-th output of the deformation network.
For training, we need to calculate the gradient of the loss-function Eq. (12) with respect to the network weights. Just like in the previous section, it is sufficient to specify the gradient with respect to the network outputs wi(α). Deriving Eq. (12) yields
d dwi L
= ∑ j d dwi ψ(x) (ψ(x)− φα(x))
= ∑ j d dwi ψ̃ (x−w(x,α)) (ψ(x)− φα(x))
=− ∑ j Xi(xj)∇ψ̃(xj −w(xj ,α)) (ψ(xj ,α)− φα(xj)) . (22)
Thus, we can calculate the derivative by summation over the region that is affected by the network output wi. The gradient term is first calculated by evaluating a finite difference stencil on ψ̃(xj) and then advecting it with the corresponding deformation vector w(xj ,α). The other terms in Eq. (22) are readily available. Alg. 1 summarizes our algorithm for training the deformation network. In particular, it is important to deform the input SDF gradients with the inferred deformation field, in order to calculate the loss gradients in the correct spatial location for backpropagation.
ALGORITHM 1: Training the deformation network Data: training samples from φα Result: trained deformation network weights Θ for each training sample {α̃, φ̃} do
evaluate neural network to compute β(α̃) load reference SDF φ̃, initial SDF ψ0 calculate vfin(xi, β(α̃)) ψ̃ = advect ψ0 with vfin calculate ∇ψ̃ evaluate neural network to compute wi(α̃) assemble w(xi) from wi(α̃,Θ) according to Eq. (21) advect ψ̃ with w advect ∇ψ̃ with w for each wi do
calculate the gradient ddwiL according to Eq. (22) end backpropagate ddwiLt from Eq. (5) to adjust Θ
end

heading: C ADDITIONAL EVALUATION
text: C.1 2D DATA SET
In the following, we explain additional details of the evaluation examples. For the two dimensional data set, we use the SDFs extracted from 2D simulations of a drop falling into a basin. As simulation parameters we choose α1 to be the size of the drop, and α2 to be its initial x-position, as shown in Fig. 13. From this simulation we extract a single frame at t = 30, which gives us a two-dimensional parameter-space α = (α1, α2), where each instance of α has a corresponding two-dimensional SDF. In order to train the networks described in section 3, we sample the parameter domain with a regular 44 × 49 grid, which gives us 2156 training samples, of which we used 100 as a validation set.
Fig. 14 shows the validation loss and the training loss over the iterations both for parameter learning and for deformation learning. We observe that in both cases the learning process reduces the loss, and finally converges to a stable solution. This value is lower in the case of deformation training, which can be easily explained with the increased expressive capabilities of the deformation network.
We verified that the solution converged by continuing training for another 36000 steps, during which the change of the solution was negligible.
As mentioned above, it might seem attractive to use a simpler approximation for the forward advection in Eq. (13), i.e., using a simpler, regular advection step. However, due to the strong non-linearity of our setting, this prevents the network from converging, as shown in Fig. 15.
The effect of our deformation network approach is illustrated in Fig. 5. This figure compares our full method (on the right) with several other algorithms. A different, but popular approach for non-linear dimensionality reduction, which can be considered as an alternative to our method, is to construct a reduced basis with PCA. Using the mean surface with four eigenvectors yields a similar reduction to our method in terms of memory footprint. We additionally re-project the different reference surfaces into the reduced basis to improve the reconstruction quality of the PCA version. However, despite this the result is a very smooth surface that fails to capture any details of the behavior of the parameter space, as can be seen in the left column of Fig. 5.
The next column of this figure (in pink) shows the surfaces obtained with the learned deformation weights with our parameter network (Fig. 12 top), but without an additional deformation network. As this case is based on end-point deformations, it cannot adapt to larger changes of surface structure in the middle of the domain. In contrast, using our full pipeline with the deformation network yields surfaces that adapt to the varying behavior in the interior of the parameter space, as shown on the right side of Fig. 5. However, it is also apparent that the deformations generated by our approach do not capture every detail of the references. The solution we retrieve is regularized by the varying reference surfaces in small neighborhoods of α, and the networks learns an averaged behavior from the inputs.
debug note: frames 7, 19, 31
C.2 4D DATA SETS
Below we give additional details for our results for the 4D data sets and experiments presented in the main document.
Liquid Drop As our first 4D test case, we chose a drop of liquid falling into a basin. As our simulation parameters we chose the x- and y-coordinates of the initial drop position, as well as the size of the drop. We typically assume that the z-axis points upwards. To generate the training data, we sample the parameter space on a regular grid, and run simulations, each with a spatial resolution of 1003 to generate a total of 1764 reference SDFs. Here, ψ0 contains a 4D SDF of a large drop falling into the upper right corner of the basin. In Fig. 17 we show additional examples how the introduction of the deformation network helps to represent the target surface across the full parameter range.
The advantages of our approach also become apparent when comparing our method with a direct interpolation of SDF data-sets, i.e., without any deformation. Our algorithms requires a single fullresolution SDF, three half resolution deformations, and the neural network weights (ca. 53.5k). While a single 404 SDF requires ca. 2.5m scalar values, all deformations and network weights require ca. 2m scalars in total. Thus our representation encodes the full behavior with less storage than two full SDFs. To illustrate this point, we show the result of a direct SDF interpolation in Fig. 18. Here we sample the parameter space with 8 SDFs in total (at all corners of the 3D parameter space). Hence, this version requires more than 4x the storage our approach requires. Despite the additional memory, the direct interpolations of SDFs lead to very obvious, and undesirable artifacts. The results shown on the right side of Fig. 18 neither represent the initial drop in (a), nor the resulting splash in (b). Rather, the SDF interpolation leads to strong ghosting artifacts, and an overall loss of detail. Instead of the single drop and splash that our method produces, it leads to four smoothed, and repeated copies. Both the PCA example above, and this direct SDF interpolation illustrate the usefulness of representing the target surface in terms of a learned deformation.
For the falling drop setup, our video also contains an additional example with a larger number of 14 pre-computed deformations. This illustrates the gains in quality that can be achieved via a larger number of deformation fields. However, in this case the parameter and deformation network only lead to negligible changes in the solution due to the closely matching surface from the pre-computed deformations.
Stairs Our second test setup illustrates a different parameter space that captures a variety of obstacle boundary conditions parametrized with α. Our first two simulation parameters are the heights of
two stair-like steps, while the third parameter is controlling the position of a middle divider obstacle, as illustrated in Fig. 19. The liquid flows in a U-shaped manner around the divider, down the steps. For this setup, we use a higher overall resolution for both space-time SDFs, as well as for the output of the deformation network. Performance details can be found in Table 1.
Fig. 20 depicts still frames captured from our mobile application for this setup. With this setup the user can adjust stair heights and wall width dynamically, while deformations are computed in the background. While this setup has more temporal coherence in its motion than the drop setup, the changing obstacle boundary conditions lead to strongly differing streams over the steps of the obstacle geometry. E.g., changing the position of the divider changes the flow from a narrow, fast stream to a slow, broad front.
Parameter 1 - raise corner Parameter 2 - lower platform Parameter 3 - wall width
Figure 19: The geometric setup of the three deformations of our stairs setup from 20 are illustrated in this figure.
Figure 20: These screens illustrate our stairs setup running in our mobile application. From left to right, the middle divider is pulled back, leading to an increased flow over the step in the back. In the right-most image, the left corner starts to move up, leading to a new stream of liquid pouring down into the outflow region in the right corner of the simulation domain.

Review summaries:
This is an application paper on dense volumetric synthesis of liquids and smoke. Given densely registered 4D implicit surfaces (volumes over time) for a structured scene, a neural-network based model is used to interpolate simulations for novel scene conditions (e.g. position and size of dropped water ball). The interpolation model composes two components -- given these conditions, it first regresses weights combining a set of precomputed deformation fields, and then a second model regresses dense volumetric deformation corrections -- these are helpful as some events are not easily modeled with a set of basis deformations. 

This paper introduces a deep learning approach for physical simulation. The approach combines two networks for synthesizing 4D data that represents 3D physical simulations. Here the first network outputs an initial guess, and the second network adds details. The first network utilizes a set of precomputed deformations, while the weights can be set to generate different output shapes. The precomputed deformations are applied in a recurrent manner. The second network is a variant of STN.