Realistic Image Synthesis Assignment 3 - Rendering Fluids
July 31, 2015
Hikaru Ikuta∗
Department of Information Physics and Computing, The University of Tokyo(a) t = 0(b) t = 4(c) t = 10(d) t = 22

Figure 1: Fluid simulation results.

Abstract
For my final project, I have implemented volume rendering forheterogeneous participating media to render scenes with realisticsmokes. To generate physically accurate density distributions ofsmokes, I have implemented fluid simulations based on the workby J.Stam [2003].
Keywords: heterogeneous participating media, fluid simulation,
path tracing

1Introduction

For my final project, I wanted to explore on implementing
computer-generated models to create sophisticated scenes. To thisend, I focused on implementing volume rendering for heterogeneous participating media, to render realistic smokes.

1.1

Summary of accomplishments

The major features I have implemented are: (i) Volume renderingfor heterogeneous participating media, and (ii) Fluid simulation.
These features will be explained in detail in this paper.
I have also implemented the following minor features: (i) Path tracing for Lambertian and specular surfaces, (ii) Importance sampler
for image based lighting, and (iii) Texture mapping for Lambertiansurfaces.

For this project, we needed a simulator for three dimensional fluids. Therefore, one main challenge of this project was to constructa three dimensional fluid simulator, based on the two dimensionaldescription. A physical understanding of the simulation method
was required to construct a working three dimensional fluid simulator based on their work. Since the [Stam 2003] focused on theimplementation of the algorithm, this was not as straightforward asmentioned in the paper. In this paper, we will study the physicaland mathematical aspects required to extend the method to a three
dimensional one.
This paper is organized as follows. First, related works of fluid simulation and volume rendering are introduced. Then, the methods
and simulation settings are described in detail. Finally, the resultsof the implementation and the discussion are given.

2

Related Work

2.1

Fluid Simulation

The fluid simulation used in this work uses a model where the fluid
space is divided into meshes, and the fluid is modeled by vector
fields and density fields. Another typical method of fluid simulation
is Smoothed Particle Hydrodynamics (SPH), as in works as [Mlleret al. 2003]. In SPH, fluids are modeled as particles, and smoothingkernels are used to obtain a continuous representation of the fluid.

Until the previous assignment, I have been working with Kento Masui, who has implemented the following features: (i) SIMD raytriangle intersection, and (ii) Multithread rendering.

One feature of SPH is that it is easy to create a more sophisticatedand accurate model, such as surface tensions in [Mller et al. 2003].
On the other hand, the vector field-based approach by [Stam 2003]
uses a more simple model and is more straightforward to implement, which is a reason why we used it for our purpose.

1.2

2.2

Challenges

The density distribution of the smokes were generated by a fluid
simulator I have implemented based on the work by J.Stam [2003].
The method introduced in [Stam 2003] simulates two dimensional
Newtonian fluids with incompressible flow, as it will be explainedlater.
∗ e-mail:hikaru ikuta@ipc.i.u-tokyo.ac.jp

Volume Rendering

The volume rendering algorithm is based on ray marching. Raymarching is a general method for evaluating values along a ray,such as line integrals, intersecting implicit surfaces such as fractals[McGuire 2014], etc. The method used in this project was inspiredby [Delalandre et al. 2010]. Our method uses ray marching to evaluate the line integral appearing in the transmittance within participating media. Other methods of volume rendering include photonmapping [Jensen and Christensen 1998], virtual ray lights [Nováket al. 2012]. In these methods, information of emitted lights are
first stored, and the radiance in the scene is then estimated using
the stored information. Among these methods, ray marching is themost simple and straightforward method, which is why we havechosen these methods for our project.

3Methods

3.1

Fluid Simulation

This equation can be interpreted as a system of linear equationsρ(t+1) (i, j, k)ρ(t+1) (i + 1, j, k)

 (t+1)ρ(t+1) (i − 1, j, k)
ρ
(i,j+1,k)1 ρ(t+1) (i, j − 1, k)
ρ(t+1) (i, j, k + 1)
ρ(t+1) (i, j, k − 1)∆t  ∆x2+6∆x2 ∆t1111

1= ρ(t) (i, j, k).

The fluid simulator is based on the Navier-Stokes equations, given
as follows:∂u= −(u · ∇)u + ν∇2 u + f∂t∂ρ= −(u · ∇)ρ + κ∇2 ρ + S,∂t(5)By collecting this equation for all i, j, and k, we obtain the linear
equation of the form

(1)Aρ(t+1) = ρ(t) ,(2)where ρ(t+1) and ρ(t) are vectors constructed by ”straighteningout” the i, j, k components of ρ(t+1) (i, j, k) and ρ(t) (i, j, k), respectively, into column vectors. Since A and ρ(t) are known, calculating the time evolution of diffusion reduces to solving this linear
system of equations for ρ(t+1) .

where u is the velocity, ρ is the density, ν is the viscosity, κ isthe diffusion coefficient, f is the external force, and S is the mass
source. In this implementation, we consider ν and κ as constants,
which describe the properties of the fluid. Fluids where ν can beconsidered as a constant is called as Newtonian fluids. Here, u andρ fully describes the state of the fluid. The terms f and S can beused as inputs to control the fluid.

(6)

The sparse structure of A can be used for efficient calculation of thesolution. In this method, an iterative method known as the GaussSeidel Method is used, as it will be explained later.

The simulation basically consists of 2 steps: evolution of the density, and evolution of the velocity. We consider a (N + 2) × (N +2) × (N + 2)-divided mesh, and calculate the time evolution of thedensity and velocity field within.

The third term of Eq. 2, S, can be implemented trivially by simply
incrementing ρ by S∆t in each time step. We also will omit thedescription of the first term, since it is throughly explained in theoriginal paper.

3.1.1 Evolution of Density

3.1.2 Evolution of Velocity

The time evolution of fluid density is described by Eq. 2. We will
first observe the second term, which is the diffusion term. Naively
discretizing Eq. 2 with the forward difference operator yields thedifference equation

As inicated in the original paper, due to the similarity of Eqs. 1 and2, evolution of velocity can be done using the exact same formulations as the evolution of density. The difference in the velocity stepis that some transformations are applied to satisfy mass conservation, for more realistic results. This appears as a preprocessing stepfor u in the algorithm.

∆tκ(∆x2(t)(t)ρ (i + 1, j, k) + ρ (i − 1, j, k)+ρ(t+1) (i, j, k) = ρ(t) (i, j, k) +(t)ρ(i, j + 1, k) + ρ

(t)

(t)

(t)

ρ(i, j, k + 1) + ρ−6ρ(t) (i, j, k) ,This step can be explained as follows. We have the assumption thatthe flow of the fluid is incompressible, i.e.

(i, j − 1, k)+∇ · u = 0.

(i, j, k − 1)(3)

where ∆t and ∆x are the time steps and mesh widths, respectively.
Eq. 3 provides an iteration for evolving ρ over time, where the density of the next step is an explicit function of the previous steps.
However, as it is mentioned in the paper, this algorithm is an unstable algorithm, i.e. ρ can diverge depending on the simulation
settings, such as values of N and κ.
The method avoids this by using a backward difference operator forthe discretization, yielding the difference equation∆tκ(∆x2
ρ(t+1) (i + 1, j, k) + ρ(t+1) (i − 1, j, k)+ρ(t) (i, j, k) = ρ(t+1) (i, j, k) −ρ(t+1) (i, j + 1, k) + ρ(t+1) (i, j − 1, k)+(t+1)ρ(t+1)(i, j, k + 1) + ρ−6ρ
(i, j, k) .
(t+1)

(i, j, k − 1)
(4)

(7)

It is known that combined with the Navier-Stokes equations Eqs.
1-2, this is equivalent to the equation∂ρ+ u · ∇ρ = 0,∂t(8)which is also known as the continuity equation, an equation that
describes mass conservation in a differential form. Therefore, wehave that incompressibility is equivalent to mass conservation.
However, the computed values of the velocity field does not necessarily satisfy Eq. 7, due to computational errors, discretization
errors, etc. The original paper mentions that the violation of massconservation hinders the quality of the simulation to a visually recognizable level. Therefore, we must somehow force the simulation
results to satisfy mass conservation.
In the algorithm, this is accomplished by decomposing the computed velocity field to a mass-conserving part and non-massconserving part, and subtracting the non-mass-conserving partaway from the computed velocity field. In the original paper, it
is mentioned that a result called Hodge decomposition is used.

Let u′ be the computed velocity field, where ∇ · u′ := g 6= 0. It is
known that every three dimensional vector field can be decomposedusing a vector potential and a scalar potential. Therefore, there always exists some vector potential χ and scalar potential φ, suchthat u can be decomposed asu′ = ∇ × χ + ∇φ.

(9)Algorithm 1 Gauss-Seidel Method
1: while Convergence criteria are not satisfied do2:
for i ← 1, 2, · · · n do(t+1)

xi
4:
end for
5: end while3:← bi −P(t)
(t+1)
− n
j=i+1 aij xj
j=1 aij xj

Pi−1

By taking the divergence of Eq. 9, we have
3.1.4 Simulation Settings∇ · u′ = ∇ · ∇ × χ + ∇2 φ,(10)

thereforeg = ∇2 φ.

(11)

Eq. 11 is known as the Poisson Equation. Suppose that we havefound solved Eq. 11 and found φ. We can then construct a newvelocity fieldu′′ := u′ − ∇φ,(12)

The initial density distribution is shown in Fig. 1 (a). Density wasdistributed randomly inside a given rectangular region on the bottom of the fluid space. The initial velocity distribution was givenso that all the voxels in the region had an almost uniform upwardvelocity with some fluctuations. The inputs f and S were set to 0 at
all times.

3.2

Volume Rendering

3.2.1 Model Construction

which must satisfy∇ · u′′ = ∇ · (u′ − ∇φ) = ∇ · (∇ × χ) = 0.

(13)

Therefore, the new velocity field u′′ satisfies Eq. 7, ultimately
meaning that it satisfies mass conservation. Therefore, we can use
this transformation from u′ to u′′ to force the computed velocity
field to satisfy the mass conservation law.

The density distrbution ρ was used for constructing the smoke
model. The smoke is characterized by the distribution of its extinction coefficient σt and scattering coefficient σs . The relationshipbetween ρ and σt was given so that ρ serves as a scaling factor ofσt :σt = ρ/800.
(15)Given σt , we gave σs as

In the algorithm, this transformation, i.e. solving Eq. 11 and using Eq. 12 is done before computing the diffusion and advection
terms for the velocity field. The Poisson Equation Eq. 11 is solved
by discretization, which yields a linear system of equations. TheGauss-Seidel method is used for solving the obtained linear systemof equations, where the sparse structure is used in this case as well.

σs = σt /2.

(16)

By definition, the relationship σt > σs must hold. This relashonship states that the albedo of the medium is 0.5 for all wavelengths,
giving the medium a gray color when illuminated with white light.
Note that our medium does not emit any radiance.

3.1.3 Gauss-Seidel Method

3.2.2 Distance sampling

The Gauss-Seidel Method is an iterative method for solving linear
systems of equations. This method aims to solve the linear system
Ax = b,(14)where A ∈ Rn×n , x ∈ Rn , b ∈ Rn , n ∈ N. For large coefficientmatrices A, it is generally computationally expensive to calculate
the inverse of A and multiplying it to x, which has a complexityof O(n3 ) with simple methods such as Gauss-Jordan elimination.
The Gauss-Seidel Method avoids this computational cost, by exploiting the structure of the matrix, such as sparsity. This method
only requires O(n2 m) time, where m is the number of iterations.
Let aij , xi , bi be the components of A, x, and b, respectively. Let(t)the tth iteration of xi be xi . Then, starting from i = 1, The GaussSeidel Method is then described as Alg. 1. A key characteristic forthis algorithm is that each component can depend on componentsof the same time step. Specifically, the first component xt+1
only1
depends on the previous steps, xki , but the second component xt+12depends on xt+1
1 , which is on the same time step. This is expressedby the summation depending on the index i. Thanks to this feature, very few memory space is required for the algorithm, since nocopies of x are required for the iteration.
The convergence criteria can be chosen arbitrarily. In this project,
we simply used a fixed number of iterations, specifically 20 iterations, for the convergence criterion.

For evaluating the transmittance of a given line segment within themedia, our method uses ray marching to evaluate the line integral
in the equation of the transmittance.
To take account of scattering, distance sampling is required to determine which point the ray was scattered inside the media. Dinstancesampling was done as shown in Alg. 2. The main idea of this sampling method is to focus on the ratio of radiance that arrives at agiven point, and regard that ratio as the probability of the light being able to pass that point. To sample distance, we first sample arandom number prob ∈ [0, 1] from a uniform distribution. Then,
we calculate the transmittance T along the ray segment that intersects with the medium. Our distance sampling method then findsthe point where T = prob. If T < prob at all points, we regard
that the ray has passed through the medium.

4Results and Discussion

4.1

Fluid Simulation

The density distribution obtained by the fluid simulation is shownin Figs. 1 a-b. The renderer for the simulation was implemented
with OpenGL. Each voxel is assigned with grayscale diffuse colors
and alpha values proportional to the density of the smoke. Thevoxels were then rendered with alpha blending, with a red-coloredlighting.

/aii

Algorithm 2 Distance Sampling
1: Intersect medium bounding box and ray2: segment ← intersection line segment of medium boundingbox and ray
3: prob ← random([0, 1])
4: S ← 0
5: for p ← points along segment do
6:
S ← S + σt (p)ds
7:
if prob < exp(−S) then8:
return p as the scattering point
9:
end if
10: end for
11: (Ray has passed through medium)

4.2

Volume Rendering and Path Tracing

The result of the entire scene is shown in Fig. 3. Notice that thefloor color is shaded upon the teapot, indicating that global illumination is accomplished by path tracing. The smoke is also renderedsuccessfully.

Figure 3: The final image, demonstrating the results of the IBL
importance sampler.

Acknowledgements
I would like to thank Kento Masui, for working together on theprevious assignments, and sharing the excitement of working on
this project. I enjoyed learning C++ techniques, and using them foraccomplishments in this project.

References
D ELALANDRE , C., G AUTRON , P., M ARVIE , J.-E., AND
F RANÇOIS , G. 2010. Single scattering in heterogenous participating media. In ACM SIGGRAPH 2010 Talks, ACM, New
York, NY, USA, SIGGRAPH ’10, 14:1–14:1.
J ENSEN , H. W., AND C HRISTENSEN , P. H. 1998. Efficient simulation of light transport in scenes with participating media usingphoton maps. In Proceedings of the 25th Annual Conference
on Computer Graphics and Interactive Techniques, ACM, New
York, NY, USA, SIGGRAPH ’98, 311–320.
Figure 2: The final image, demonstrating the results of volume rendering.

5Conclusion

The aim for this final assignment was to implement and use
computer-generated models for sophisticated scenes. To this end,I have successfully implemented volume rendering for heterogeneous participating media, based on ray marching and path tracing.
I have also successfully implemented a 3 dimensional fluid simulator, based on the work by J.Stam [2003]. Using the simulator, I
have generated a density distribution of smoke, and finally rendereda scene containing smokes illuminated with image-based lighting
and point lights.
For future works, photon mapping for enhancing volume rendering
can be considered.

L ANDIS , H., 2002. Global illumination in production. ACM SIGGRAPH 2002 Course #16 Notes, July.
M C G UIRE , M., 2014. Numerical methods for ray tracing implicitlydefined surfaces. Williams College 2014 Course CS371 Notes,September.
M LLER , M., C HARYPAR , D., AND G ROSS , M., 2003. Particlebased fluid simulation for interactive applications.
N OV ÁK , J., N OWROUZEZAHRAI , D., DACHSBACHER , C., AND
JAROSZ , W. 2012. Virtual ray lights for rendering scenes withparticipating media. ACM Trans. Graph. 31, 4 (July), 60:1–60:11.
S TAM , J. 2003. Real-time fluid dynamics for games. Proceedingsof the Game Developer Conference.
W IKIPEDIA. Gauss-seidel method.

W IKIPEDIA. Incompressible flow.