Example-based Synthesis of Turbulence by Flow Field Style Transfer
Syuhei Sato∗
UEI Research

Yoshinori Dobashi
Hokkaido University / UEI Research

Tomoyuki Nishita
UEI Research / Hiroshima Shudo University
local patch

Abstract

global patchXu(p)The computational cost of physically-based fluid simulation forgenerating realistic animations is expensive. Therefore, it takes along time for an animator to create the desired animation. To reduce such costs, several methods have been developed that employ
an approach in which turbulence is synthesized as a post process.
Since, with this approach, global motion can be obtained using lowresolution fluid simulation, realistic animations with the desired behavior can be created at low cost. In the field of image editing,style transfer methods, which can efficiently achieve the desired
stylization by transferring features of an example image to a userspecified image, have been proposed. We combine these concepts
and present a novel style transfer method for turbulence by reusing
the existing fluid animations. Our system allows the user to transfer
turbulent motion from one flow field to other flow fields. We do
this by extending example-based image synthesis methods to theflow field: a patch-based synthesis and an optimization-based texture synthesis. Our method is designed such that the resulting flow
fields satisfy the incompressibility condition. Our method can intuitively and easily create high-resolution fluid animations with thedesired turbulent motion.

p

xu(p)
p

B

b

B
b

Ut

ut

Figure 1: Definition of local (the orange squares) and global (theblue square) patches. For visual clarity, we show only one global
patch, but actually the global patch is defined for each position p.

of our method, with which the desired turbulence can be intuitively
and easily transferred onto a user designed low-resolution simulation.

Keywords: fluid simulation, reusing existing fluid animations, incompressibility, texture synthesis, patch-based synthesis

2

Concepts: •Computing methodologies → Animation; Physicalsimulation;

Stam developed a semi-Lagrangian method [Stam 1999], which isan unconditionally stable solver for the Navier-Stokes equations.
Since this method was introduced, many methods have been proposed for simulating various phenomena. Although realistic animations can be produced with these methods, the user must repeat
the simulation and search for good parameter settings until the desired animation is obtained.

1

Introduction

In many entertainment applications, such as movies and computer
games, physically-based fluid simulation has generally been used.
However, one of the problems of using this is the expensive computational cost, which makes it tedious to adjust the parameters tocreate the desired fluid animation. The animator must repeat thefluid simulation multiple times to find an appropriate set of parameters that can produce a satisfactory result. To avoid such costly
simulation, post-processing approaches for synthesizing turbulence
have been proposed: e.g., adding plausible turbulence [Kim et al.
2008], guide-based simulation [Nielsen and Christensen 2010].Using these methods, the time for creating the desired animation canbe significantly reduced.
In this paper, we also propose a method for synthesizing turbulence but take a different approach to those used in the methods
previously presented. Our method provides a tool for transferring
turbulent motion between two input flow fields. The turbulencein a high-resolution flow field can be transferred onto a user designed low-resolution flow field. This transfer is done by extending an example-based image synthesis approach to the flow field.
Firstly, we apply a patch-based synthesis to transfer turbulent motion while preserving the global distribution of the turbulence. Secondly, optimization-based texture synthesis is applied to the boundary of each patch to iron out the discontinuities and satisfy the incompressibility condition. We call our approach turbulence styletransfer. In this paper, we apply this method to smoke simulation
defined on a uniform grid. We show several examples of the use∗ e-mail:syuhei sato@dwango.co.jp

Related Work

Several post-processing approaches for synthesizing turbulence
have been proposed. Kim et al. [2008] synthesized detailed turbulent motion onto low-resolution simulation using wavelet noise.
Neilsen et al. [2010] controlled a high-resolution fluid simulationaccording to a user designed low-resolution simulation by solving aminimization problem. Chu et al. [2017] proposed a method basedon a convolutional neural network. By contrast, we reuse existingflow fields with the desired turbulent motion, and then the turbulence is transferred onto any other flow field. Therefore, the desired
turbulence can be intuitively and easily synthesized just by specifying such flow field.
In the literature on texture synthesis, many optimization-based
methods have been proposed. One popular method was proposed
by Kwatra et al. [2005]. This method repeats a nearest neighborsearch and solving linear system, and then a texture sufficiently
similar to an example texture can be synthesized. Kwatra et al.
[2007] also proposed a method for synthesizing the texture on amesh representing a water surface. Narain et al. [2007] extended
the method of Kwatra et al. [2007] such that details can be added
to a water surface from example images. When multiple images,
such as those representing water, foam, and so on, are given to thissystem, an appropriate image is automatically selected according
to the local features of the water mesh. In our method, we focus
on smoke defined on a grid. Therefore, those methods that focus
on water meshes are not suitable for our purpose. Jamriska et al.
[2015] proposed an appearance transfer method for 2D fluid ani-global searchsearch regionxuplocal search
Xu
zu

b

search region
qh

pB
B

ql
b

Zu
target fieldsource fieldtarget fieldsource field

Figure 2: Two-level patch search. A similar region is globally searched using the global patches (the blue squares), and then locally searchedusing the local patches (the orange squares). For the global and local search, a search region is shown by the red dotted squares.

mations based on the above optimization-based texture synthesis.
This method distinguishes between the features of the boundary
and the inner region of the fluid using an alpha mask. Similar tothis method, we create a mask that represents the presence of fluid
from a density field and use it for patch detection.
Ma et al. [2009] synthesized the small-scale detail in a motion
field from an example image using Kwatra’s method [Kwatra et al.
2005]. This motion field is then added to an input low-resolution
motion field; however, the detailed motion added is uniform overthe entire field. In fluids, the size and distribution of turbulence
varies. Therefore, this method [Ma et al. 2009] is not suitable forour purpose. We take into account features of the flow field bycalculating the similarity using the velocity and density fields tosynthesize plausible turbulence.

3

Our Method

In this section, we describe our example-based turbulence synthesis. For the sake of clarity, we describe the method for a 2D flow
field. Two input flow fields are required and we assume that the
number of frames for these flow fields is the same, though themethod is not restricted to this assumption. In the following description, we use lower and upper case symbols to represent variables related to the low and high resolution grids, respectively.

3.1

Input Flow Fields

The inputs to our method are a low-resolution flow field with thedesired global motion (target field) and a high-resolution flow fieldwith the desired turbulent motion (source field). The target andsource fields consist of velocity fields (ut (i) and Us (i)) and density masks (dt (i) and Ds (i)), where i = 1, 2, · · · , I represents theframe number, and I is the number of frames. dt (i) and Ds (i)are set to 1 for the grid points where the smoke density is larger
than Dth , and set to 0 otherwise. Dth is a threshold for the density
which is chosen so that the smoke region is extracted. The grid sizeand grid interval for the target and source fields are n × m, N × M
and ∆x, ∆X, respectively (n < N , m < M and ∆x > ∆X).
From these inputs, we evaluate the similarity of the flow and theshape of the smoke.
The source field is firstly downsampled so that the grid size becomes the same as that of the target field (n × m), and the lowfrequency components of the velocity and density fields of thesource field are obtained, denoted by us (i) and ds (i) respectively.
Next, ut (i), dt (i) and us (i) are linearly upsampled to the resolution of the source field (N × M ), denoted as Ut (i), Dt (i) andUsl (i), respectively. The high-frequency components of the source
field Ush (i) are then computed by Us (i) − Usl (i).

Next, we subdivide the target field into small square regions, which
we call patches. Ut (and Dt ) are subdivided into patches, shown
as orange squares in Fig. 1, with grid size B × B and grid interval∆X. We denote the patches at p in the velocity and density fields
by Xu (p) and Xd (p), respectively. We call these local patches.
In addition, for each center p, a square region whose size and gridinterval are b × b and ∆x (b > B × ∆X/∆x) is defined on ut anddt (the blue square in Fig. 1). These patches are denoted as xu (p)and xd (p) for the velocity and density fields, respectively, and wecall these global patches. The global patches are used for the global
search over the entire simulation space while the local patches areused for a local search inside a global patch for a similar region (seeSection 3.2). In the case of 3D flow fields, we define the patches
using voxels with a grid size of B 3 or b3 .

3.2

Patch-based Turbulence Synthesis

Our method adds high-frequency turbulent motion to the target fieldby searching for the best patch in the source field that is most similar to each of the patches in the target field. To efficiently detect
the patch, we develop a two-level search algorithm using two resolutions (low and high) as shown in Fig. 2. First, for each local
patch, we check the existence of the smoke using Dt (i). When thesmoke exists in a patch, the following two-level search process is
executed. That is, a similar region is globally detected using the
global patches on a low-resolution grid, and the local search is thendone inside the region of the detected global patch. However, wefound that repeating this two-level search process at each frame independently produced temporal artifacts in the resulting flow field.
We address this problem as follows. The global search is done only
once at the frame when the smoke firstly appears in the patch and,at the subsequent frames, the global patch is simply translated sothat its center coincides with the center of the optimal local patch
computed in the previous frame. After searching for and finding
the best local patch, we use this to synthesize the high-frequency
components Uh (i).
In the global search process, we find the best global patch in thesource patch that is most similar to each of the global patches inthe target field (see left in Fig. 2). We search the low-resolution (ordownsampled) source field for the position ql such that the following energy function is minimized.
El (p, ql ) = ||xu (p) − zu (ql )||2 + α||xd (p) − zd (ql )||2 , (1)where α is a user specified coefficient for adjusting the influenceof the second term, the density difference, on the above equation.
zu (ql ) and zd (ql ) are b × b size patches whose center is locatedat ql in the velocity and density fields, us and ds , respectively, inthe downsampled source field. By the above search process, we canprocessing region Ω

using the velocity gradients. In order to satisfy the incompressibility of the flow, the divergence of the velocity is also introduced tothe energy function. More concretely, we minimize the following
energy function Eb at the boundary region Ω, in order to remove
the discontinuities between the local patches.

baundary patchr

C

C X(r)

Eb (r, s) =X||∇X(r) − ∇Z(s)||2 + γ(∇ · X(r))2 ,

(3)

r,s

Uh

Figure 3: Definition of a processing region Ω (the green shaded
regions) and a boundary patch (the green square).

find the best region in the source field corresponding to each of the
global patches in the target region.
Next, in the local search process, we search the region in the source
field for the local patch that is most similar to it in the corresponding global patch in the target field (see right in Fig. 2). To take into
account temporal coherence, we introduce the high-frequency components of the velocity field Ûh (i) obtained by advecting the highfrequency components synthesized in the (i−1)-th frame Uh (i−1)
according to the target field Ut (i − 1), and its patch is denoted as
Xû (p). Then, we search the search region in the source field (thered dotted square in Fig. 2) for the center position qh of the local patch whose high-frequency components are similar to those ofXû (p) by minimizing the following energy function, Eh :where γ is a relative weighting coefficient.X(r) and Z(s) are C ×
C size patches with centers located at r and s on Uh and Ush ,respectively. These patches are called boundary patch (the green
square in Fig. 3). To minimize Eq.(3), we adopt an ExpectationMaximization like approach, the same as Kwatra et al. [2005].
In the E-step, X(r) is calculated by solving the linear system ofEq. (3), while keeping Z fixed at the value obtained in the M-step.
This is done by setting the derivation of Eq. (3) with respect to Uh
to zero. To obtain a numerical solution to the linear system, we
use the conjugate gradient (CG) method. This linear system is very
similar to that used in the previous methods [Sato et al. 2016].Next,
in the M-step, we search Z(s) on Ush keeping X fixed at the value
estimated in the E-step, such that Eq. (3) is minimized. Since X
is fixed, the second term of Eq. (3) becomes constant in this Mstep. Therefore, the energy Eb is calculated using only the firstterm of Eq. (3). We alternately update X and Z by repeating theabove E-step and M-step. When the energy transition Eb is sufficiently small, we terminate this smoothing process, and obtain theboundary-smoothed high-frequency components Uh . Finally, theresultant velocity field U is obtained as Uh + Ut . To visualize theflow, the smoke density is advected along with U.

Eh (p, qh ) = ||Xu (p) − Zu (qh )||2

4

+ α||Xd (p) − Zd (qh )||22

+ β||Xû (p) − Zû (qh )|| ,(2)

where β is a user specified coefficient for adjusting the influence ofthe third term on Eq. 2. Zu (qh ), Zd (qh ) and Zû (qh ) are B × B
size patches for Usl , Ds and Ush , respectively, with the center at
qh . Using the above search, the most similar patch is obtained fromthe source field for each local patch of the target field, and then we
synthesize the high-frequency components of the velocity field Uhusing these local patches. However, since we have ignored the spatial continuity between local patches, the velocities are discontinuous at the boundary. Wherefore, we eliminate these discontinuities
by a smoothing process described in the next section.

3.3

Smoothing Velocities between Patches

In order to make the velocities smooth at the boundaries between
local patches, we apply an optimization-based texture synthesis
approach [Kwatra et al. 2005] to Uh . However, since Kwatra’s
method does not take into account the incompressibility of the fluid,
we cannot apply this method to the flow fields in a straightforwardmanner. We modify the energy function in the original method[Kwatra et al. 2005] to satisfy the incompressibility. The detailsof this process and our modification are described in the following.
First, a processing region, Ω, with a predefined width is createdaround the boundaries between each local patch on Uh (the green
shaded regions in Fig. 3), where the width is specified by the user.
In this process, only the high-frequency components (Uh and Ush )are considered. The high-frequency components contain the smallscale turbulent motion and vorticity, which are captured in the spatial derivative of the velocity field, but smooth flow (e.g., laminarflow) is not contained. Therefore, we formulate our energy function

ResultsFigs. 4, 5 show results created using our method. We used a desktop
PC with an Intel Core i7-5820K CPU, 32GB memory to computeall the examples. α, β and γ are set to 0.004, 0.5 and 10, respectively. The videos corresponding to these examples can be found in
the supplementary material.
Fig. 4 shows results synthesized using various source fields withdifferent turbulent motion, in 2D flow field. Grid size of target andsource fields is 24 × 32 and 192 × 256. Fig. 4 (a) shows the target
field linearly upsampled to 192 × 256 grid size. Figs. 4 (b1), (c1)and (d1) are source fields. Figs. 4 (b2), (c2) and (d2) are our results
created by transferring the turbulent motion of (b1), (c1) and (d1)to (a). The computation time for our method and 192 × 256 size
fluid simulation is 0.7 and 1.0 seconds per frame on average.
Fig. 5 shows various results in 3D flow field. Grid size of target andsource fields is 16 × 16 × 24 and 128 × 128 × 192. Fig. 5 (a) shows
the target field linearly upsampled to the size of the source field.
Figs. 5 (b1) and (c1) are source fields. Figs. 5 (b2) and (c2) are ourresults created by transferring the turbulent motion of (b1) and (c1)to (a). The computation time for our method and 128 × 128 × 192
size fluid simulation is 270 and 50 seconds per frame on average.
In this 3D case, our method is 5 times slower than the simulation.
The most expensive process in our method is the search process.
However, since the search process for each patch is independent,
we expect a significant acceleration would be achieved by using
the GPU to execute the search process for each patch in parallel.
Therefore, fast computing can expect by implementing the method
on GPU. As shown in our 2D and 3D results, the turbulent motion
is similar to each source field, and the global motion is preserved tobe same with the target field.

(b1) source 1(a) target(b2) result 1(c1) source 2(c2) result 2(d1) source 3(b1) source 1(c1) source 2(b2) result 1(c) result 2(d2) result 3

Figure 4: Experimental results in 2D flow field.

5Conclusions

We have proposed a method for synthesizing turbulent motion byreusing existing flow fields while preserving the incompressibilityof the fluid. Our method used patches to synthesize turbulent motion while preserving the global distribution of the turbulence. In
addition, we applied an optimization-based texture synthesis approach to smooth the velocity at the boundary of each patch. Furthermore, the energy function for the texture synthesis is modified
such that the incompressibility is satisfied.
One of the limitations of our method is that it is difficult to apply it
to flow fields with largely different directions and speeds, because
we evaluate the velocities directly as described in Sec. 3.2. To address this problem, we plan to evaluate flow fields after rotating andnormalizing them. We also plan to extend our method to other typesof fluid, such as fire, cloud, and water.

Acknowledgements(a) target

Figure 5: Experimental results in 3D flow field.

K WATRA , V., A DALSTEINSSON , D., K IM , T., K WATRA , N.,
C ARLSON , M., AND L IN , M. C. 2007. Texturing fluids. IEEE
Transactions on Visualization and Computer Graphics 13, 5,
939–952.
M A , C., W EI , L., G UO , B., AND Z HOU , K. 2009. Motion fieldtexture synthesis. ACM Transactions on Graphics 28, 5, Article
110.
NARAIN , R., K WATRA , V., L EE , H.-P., K IM , T., C ARLSON , M.,
AND L IN , M. C. 2007. Feature-guided dynamic texture synthesis on continuous flows. In Proceedings of the 18th Eurographics
conference on Rendering Techniques, 361–370.

This work was supported by JSPS KAKENHI Grant Number
JP15H05924.

N IELSEN , M. B., AND C HRISTENSEN , B. B. 2010. Improved
variational guiding of smoke animations. Computer Graphics
Forum 29, 2, 705–712.

References

S ATO , S., D OBASHI , Y., AND N ISHITA , T. 2016. A combining
method of fluid animations by interpolating flow fields. In Proceedings of SIGGRAPH Asia 2016 Technical Briefs, Article 4.

C HU , M., AND T HUEREY, N. 2017. Data-driven synthesis ofsmoke flows with cnn-based feature descriptors. ACM Transactions on Graphics 36, 4, Article 14.
JAMRISKA , O., F ISER , J., A SENTE , P., L U , J., S HECHTMAN , E.,
AND S YKORA , D. 2015. Lazyfluids: appearance transfer for
fluid animations. ACM Transactions on Graphics 34, 4, Article
92.
K IM , T., T HUREY, N., JAMES , D., AND G ROSS , M. 2008.
Wavelet turbulence for fluid simulation. ACM Transactions on
Graphics 27, 3, Article 3.
K WATRA , V., E SSA , I., B OBICK , A., AND K WATRA , N. 2005.
Texture optimization for example-based synthesis. ACM Transactions on Graphics 24, 3, 795–802.

S TAM , J. 1999. Stable fluids. In Proceedings of ACM SIGGRAPH
1999, Annual Conference Series, 121–128.