• No results found

Efficient implementation of the Particle Level Set method

N/A
N/A
Protected

Academic year: 2021

Share "Efficient implementation of the Particle Level Set method"

Copied!
68
0
0

Loading.... (view fulltext now)

Full text

(1)LiU-ITN-TEK-A--10/050--SE. Efficient implementation of the Particle Level Set method John Johansson 2010-09-02. Department of Science and Technology Linköping University SE-601 74 Norrköping, Sweden. Institutionen för teknik och naturvetenskap Linköpings Universitet 601 74 Norrköping.

(2) LiU-ITN-TEK-A--10/050--SE. Efficient implementation of the Particle Level Set method Examensarbete utfört i medieteknik vid Tekniska Högskolan vid Linköpings universitet. John Johansson Handledare Doug Roble Examinator Jonas Unger Norrköping 2010-09-02.

(3) Upphovsrätt Detta dokument hålls tillgängligt på Internet – eller dess framtida ersättare – under en längre tid från publiceringsdatum under förutsättning att inga extraordinära omständigheter uppstår. Tillgång till dokumentet innebär tillstånd för var och en att läsa, ladda ner, skriva ut enstaka kopior för enskilt bruk och att använda det oförändrat för ickekommersiell forskning och för undervisning. Överföring av upphovsrätten vid en senare tidpunkt kan inte upphäva detta tillstånd. All annan användning av dokumentet kräver upphovsmannens medgivande. För att garantera äktheten, säkerheten och tillgängligheten finns det lösningar av teknisk och administrativ art. Upphovsmannens ideella rätt innefattar rätt att bli nämnd som upphovsman i den omfattning som god sed kräver vid användning av dokumentet på ovan beskrivna sätt samt skydd mot att dokumentet ändras eller presenteras i sådan form eller i sådant sammanhang som är kränkande för upphovsmannens litterära eller konstnärliga anseende eller egenart. För ytterligare information om Linköping University Electronic Press se förlagets hemsida http://www.ep.liu.se/ Copyright The publishers will keep this document online on the Internet - or its possible replacement - for a considerable time from the date of publication barring exceptional circumstances. The online availability of the document implies a permanent permission for anyone to read, to download, to print out single copies for your own use and to use it unchanged for any non-commercial research and educational purpose. Subsequent transfers of copyright cannot revoke this permission. All other uses of the document are conditional on the consent of the copyright owner. The publisher has taken technical and administrative measures to assure authenticity, security and accessibility. According to intellectual property law the author has the right to be mentioned when his/her work is accessed as described above and to be protected against infringement. For additional information about the Linköping University Electronic Press and its procedures for publication and for assurance of document integrity, please refer to its WWW home page: http://www.ep.liu.se/. © John Johansson.

(4) Efficient implementation of the Particle Level set method. John Johansson1 Department of Science and Technology Link¨oping University. A thesis submitted for the degree of Master of Science September 10, 2010. 1. johjo729@student.liu.se.

(5) placeholder.

(6) Abstract. The Particle Level set method is a successful extension to Level set methods to improve the volume preservation in fluid simulations. This thesis will analyze how sparse volume data structures can be used to store both the signed distance function and the particles in order to improve access speed and memory efficiency. This Particle Level set implementation will be evaluated against Digital Domains current Particle Level set implementation. Different degrees of quantization will be used to implement particle representations with varying accuracy. These particles will be tested and both visual results and error measurments will be presented. The sparse volume data structures DB-Grid and Field3D will be evaluated in terms of speed and memory efficiency..

(7) placeholder.

(8) Contents 1 Introduction 1.1 Surface tracking . . . . . . . 1.2 The Particle Level set method 1.3 Digital Domain . . . . . . . . 1.4 Sparse volume data structures 1.5 Motivation . . . . . . . . . . 1.6 Aim . . . . . . . . . . . . . . 1.7 Structure of the report . . . .. . . . . . . .. . . . . . . .. . . . . . . .. . . . . . . .. . . . . . . .. . . . . . . .. . . . . . . .. . . . . . . .. . . . . . . .. . . . . . . .. . . . . . . .. . . . . . . .. . . . . . . .. . . . . . . .. . . . . . . .. . . . . . . .. . . . . . . .. . . . . . . .. . . . . . . .. . . . . . . .. . . . . . . .. . . . . . . .. . . . . . . .. . . . . . . .. . . . . . . .. 1 1 2 2 2 2 3 3. 2 Background 2.1 Implicit functions . . . . . . . . . 2.1.1 Implicit surfaces . . . . . 2.1.2 Geometric properties . . . 2.1.3 Signed distance functions 2.2 Level set methods . . . . . . . . 2.2.1 Discretization . . . . . . . 2.2.2 Advection . . . . . . . . . 2.2.3 Spatial discretization . . . 2.2.4 Temporal discretization . 2.2.5 Renormalization . . . . . 2.2.6 Optimizations . . . . . . . 2.3 Particle Level set method . . . . 2.3.1 Overview . . . . . . . . . 2.3.2 Massless particles . . . . . 2.3.3 Initialization . . . . . . . 2.3.4 Advection . . . . . . . . . 2.3.5 Error correction . . . . . 2.3.6 Radius adjustment . . . . 2.3.7 Reseeding . . . . . . . . . 2.4 Sparse volume data structures . . 2.4.1 DB-Grid . . . . . . . . . . 2.4.2 Field3D . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . .. 5 5 5 6 6 7 7 10 10 12 14 15 16 16 16 17 17 18 19 19 19 20 22. . . . . . . .. v.

(9) CONTENTS. CONTENTS. 3 Implementation 3.1 Current implementation . . . . . . . . . . 3.2 Storing the particles . . . . . . . . . . . . 3.2.1 Class Interface . . . . . . . . . . . 3.2.2 Particle pool . . . . . . . . . . . . 3.2.3 Multiple arrays . . . . . . . . . . . 3.2.4 STL containers . . . . . . . . . . . 3.2.5 Traversing particles . . . . . . . . 3.3 Particle representation . . . . . . . . . . . 3.3.1 Normal implementation . . . . . . 3.3.2 Semi-Compact implementation . . 3.3.3 Compact implementation . . . . . 3.4 Dynamic seeding condition . . . . . . . . 3.5 Putting the Particle Level set together . . 3.5.1 Grids . . . . . . . . . . . . . . . . 3.5.2 Level set advection . . . . . . . . . 3.5.3 Level set renomalization . . . . . . 3.5.4 Particle advection . . . . . . . . . 3.5.5 Radius adjustment and Reseeding. . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . .. 23 23 23 24 25 25 26 27 28 28 28 29 30 30 31 31 31 31 31. 4 Results 4.1 Evaluating DB-Grid and Field3D . . . . . . . . 4.1.1 Random access . . . . . . . . . . . . . . 4.1.2 Sequential access . . . . . . . . . . . . . 4.2 Particle Level set implementation . . . . . . . . 4.2.1 Deformation test . . . . . . . . . . . . . 4.2.2 Smart seeding . . . . . . . . . . . . . . . 4.2.3 Cache behaviour and threading . . . . . 4.3 Particle implementations . . . . . . . . . . . . . 4.3.1 Deformation Test . . . . . . . . . . . . . 4.3.2 Error measurement . . . . . . . . . . . . 4.4 Comparing Particle Level set implementations .. . . . . . . . . . . .. . . . . . . . . . . .. . . . . . . . . . . .. . . . . . . . . . . .. . . . . . . . . . . .. . . . . . . . . . . .. . . . . . . . . . . .. . . . . . . . . . . .. . . . . . . . . . . .. . . . . . . . . . . .. . . . . . . . . . . .. . . . . . . . . . . .. . . . . . . . . . . .. . . . . . . . . . . .. . . . . . . . . . . .. . . . . . . . . . . .. 33 33 33 33 35 35 36 36 38 38 38 40. 5 Discussion. . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . .. 41. 6 Conclusion 43 6.1 Future work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 References. 46. A Hamilton-Jacobi ENO. 51. B Hamilton-Jacobi WENO. 53. C Additional results. 55. vi.

(10) List of Figures 2.1 2.2 2.3 2.4 2.5 2.6. The level set of a 2D implicit sphere. . . . . . . . . . . . . . . . . Discrete representations in two dimensions. . . . . . . . . . . . . Trilinear interpolation . . . . . . . . . . . . . . . . . . . . . . . . DB-Grid example in 2D. . . . . . . . . . . . . . . . . . . . . . . . Scheme of the block array . . . . . . . . . . . . . . . . . . . . . . The Stanford Dragon represented with level sets using DB-Grid.. . . . . . .. . . . . . .. . . . . . .. . . . . . .. 3.1 3.2. Particle pool structure. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 Multiple arrays structure. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26. 4.1 4.2 4.3 4.4 4.5 4.6. Deformation test performed on Level set. . . . . Deformation test performed on Particle Level set. Deformation test performed on Particle Level set, Normal particle implementation . . . . . . . . . . Semi-compact particle implementation . . . . . . Compact particle implementation . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . smart seeding is used . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . .. . . . . . .. . 6 . 8 . 9 . 21 . 21 . 22. . . . . . .. 35 36 37 38 39 39. C.1 Higher resolution Particle Level set deformation test. . . . . . . . . . . . . . 55. vii.

(11)

(12) List of Tables 4.1 4.2 4.3 4.4 4.5 4.6 4.7 4.8. DB-Grid and Field3D random access test. . . . . . . . . . . . . . . . . . . . DB-Grid and Field3D sequential access test. . . . . . . . . . . . . . . . . . . Statistics from the Deformation test with and without smart seeding. . . . . Particle Level set advection with varying block size and number of threads. Statistics from the Deformation test with different particle implementations. Comparing error of the different particle implementations. . . . . . . . . . . Comparing Particle Level set implementations. . . . . . . . . . . . . . . . . Time breakdown of the new Particle Level set implementation. . . . . . . .. ix. 34 34 36 37 38 39 40 40.

(13)

(14) Chapter 1. Introduction Since computer generated special effects, or Visual effects as it is called, started to get popular in the early 90’s, fluids such as water, fire and smoke has been one of the most common effects. The reason is that it is easier and cheaper to simulate the fluid than to shot the real thing. It is also easier to direct the fluid and make it behave like the director wants. In order to create fluid effects that look natural and realistic, the movement and apperance is based on equations that describe the motion of fluid flow. These equations are known as the Navier-Stokes equations. When simulating liquid using these equations, the surface has to be tracked using a special method.. 1.1. Surface tracking. One of the first successful methods of tracking propagating surfaces was the Level set methods, introduced by Osher and Sethian in 1988 [26]. This method tracks the surface using a signed distance function, that is represented using a Eulerian grid. The method is widely used in other fields such as computational physics, image processing, computer vision, medical image analysis and 3D reconstruction [25]. The method has however problems conserving volume when undergoing deformation. Another Eulerian method of tracking the liquid surface is the Volume of Fluid (VOF) [11], where each cell in a grid stores how much percent water it contains. The problem with this method is that geometric properties can not be computed in the same straigtforward way as with Level set methods. A different approach of tracking a surface is the Lagrangian representation, where particles are used to track the surface. This was first introduced in [10], where they used marker particles to track the surface. The method tracks the interface with a much higher accuracy and preserves volume better than grid-based methods, but it is harder to calculate geometric properties such as gradient and curvature. Extracting a smooth surface mesh that is used in rendering is easier on a grid than from a large number of particles. Recent research has introduced two new methods that uses particles to track the fluid, Smoothed Particle Hydrodynamics (SPH) [20] and Fluid-Implicit-Particle (FLIP) [36]. In these methods the particles are used to track the whole volume of the fluid instead of just the surface.. 1.

(15) 1.2 The Particle Level set method. 1.2. 1. INTRODUCTION. The Particle Level set method. The Particle Level set (PLS) method [5] is an interface capturing method that uses the advantages of both Eulerian and Lagrangian methods. The main representation is a signed distance function evolved using Level set methods and it uses auxiliary marker particles to accurately track the surface and correct the signed distance function. The method has been very popular and succefully used in both research and in Visual Effects. In research it was initially used in fluid simulation to track water surfaces [7], but has since then also been used to track surfaces when simulating multiple interacting fluids [18], simulating interaction between rigid objects, deformable objects and fluids [30] and simulating water with splashes and spray [19]. To create spray, the authors used escaped particles to trigger a separate particle simulation based on the SPH method. Examples of when the PLS was used in Visual Effects includes simulation of a melting character in Terminator 3: Rise of the Machines [34], when simulating water pouring into a sinking ship in Poseidon [9] and when simulating water flooding New York City in Day After Tomorrow [13].. 1.3. Digital Domain. Digital Domain is a Visual Effects studio located in Venice, California. They have been involved in 80+ feature films and have earned three Academy awards for Best Visual Effects and three Scientific and Technical Achievement awards. Their fluid simulation engine fsim has earned them one of the Scientific and Technical Achievement awards. It is a grid-based solver that can simulate both liquids and fire. When simulating liquid it uses the Particle Level set method to track surfaces.. 1.4. Sparse volume data structures. In the last couple of years, scientists at Digital Domain have developed an efficient data structure for sparse volume data, the DB-Grid [22]. A similar data structure, Field3D [35], has been developed at Sony Pictures Imageworks and made public as open source. Both these data structures are very capable of representing the data used for the level set method and the Particle Level set method. These new data structures are based on the simple concept of dividing the grid domain into blocks and only allocating memory in blocks where data is explicit set. According to [21], DB-Grid should theoretically be faster than any previous sparse level set data structure.. 1.5. Motivation. The current Particle Level set implementation at Digital Domain is considered one of the slowest and most memory consuming parts of the fluid simulation. The company wants a Particle Level set implementation that is based on one of the new sparse volume data structures. This would make it more efficient both in terms of speed and memory consumption, which would make it possible to perform simulations at higher resolutions with faster execution speed.. 2.

(16) 1. INTRODUCTION. 1.6. 1.6 Aim. Aim. The first step is to evaluate the two new sparse data structures, DB-Grid and Field3D. The evaluation should test their performances in random access and sequential access. The goal is to determine which one is best in terms of speed and memory efficiency. The data structure that proves to be superior should then be used as the foundation for a new C++ implementation of the Particle Level set method. The aim of the implementation is to achieve a working prototype where both the signed distance function and the particles are stored using the sparse volume data structure. Regarding the storage of particles, a general purpose particle storage structure is requested to be developed. To reduce the memory footprint of the particles, different levels of particle quantization should be implemented and tested. The final implementation needs to be tested against the regular Level set solution and Digital Domains current Particle Level set implementation.. 1.7. Structure of the report. This report is structured as following: • Chapter 2 describes the theory that is needed in order to implement a Particle Level set solution and describe the two sparse volume data structures. • Chapter 3 will give specific details on our implementation. • Chapter 4 presents results from the tests that was carried out with the implementation. • Chapter 5 will discuss the implementation and results from the tests. • Chapter 6 concludes the report and lists what can be done to improve the implementation.. 3.

(17)

(18) Chapter 2. Background This chapter will describe the theoretical base that is needed for the implementation in the next chapter. It will describe the math behind implicit functions and signed distance functions, how Level set methods are used to evolve these functions over time and how the Level set methods are extended with particles to improve accuracy. The last section will describe concept and details of new sparse volume data structures.. 2.1. Implicit functions. In Level set methods, a signed distance function is used to represent the surface we want to track [26]. A signed distance function is a three dimensional implicit function, where the surface is defined as 0. Mathematical functions can be classified as either implicit or explicit. An implicit function is defined as f (x, y, z) = k, f : <3 → < (2.1) where the function f depends on x, y, z and equals a scalar value k. An explicit function on the other hand, is a function that depends on x, y and explicitly gives the value of z. z = f (x, y). (2.2). As an example, the implicit function of a sphere with radius r can be expressed as: f (x, y, z) = x2 + y 2 + z 2 − r2. 2.1.1. (2.3). Implicit surfaces. A implicit surface is a contour in the domain Ω given by a three dimensional implicit function. The surface can be extracted by only selecting the set of points where φ (x, y, x) = c. The constant c is called the iso value and the set of the values that is selected is called a iso-surface or level set of the implicit function φ (x, y, x). The zero level set or interface is denoted as Γ. The constant c can be varied to select different subsets of the function values, but the most common value to use as the iso value c is 0. This makes it easy to classify a value in the domain Ω because if the φ value is > 0 it is located outside the. 5.

(19) 2.1 Implicit functions. 2. BACKGROUND. surface and if it is < 0 it is located inside the surface. The implicit function with c = 0 can be defined as φ (~x) > 0 if ~x ∈ Ω+ φ (~x) = 0 if ~x ∈ Γ (2.4) φ (~x) < 0 if ~x ∈ Ω−. <0 >0. =0. Figure 2.1: The level set of a 2D implicit sphere.. 2.1.2. Geometric properties. Implicit functions has very convinient tools for analytically calculating geometric properties such as gradient, normal and mean curvature. The gradient of the implicit function is defined as   ∂φ ∂φ ∂φ , , , (2.5) Oφ = ∂x ∂y ∂z which means the partial derivative in all dimensions. The normal of the implicit function is expressed as the normalized gradient ~ = Oφ N |Oφ|. (2.6). The mean curvature of the implicit function is defined as the divergence of the normal   Oφ κ=O· (2.7) |Oφ| This can be expanded to the following expression of first and second partial derivatives of φ: κ = (φ2x φyy − 2φx φy φxy + φ2y φxx + φ2x φzz − 2φx φz φxz (2.8) +φ2z φxx + φ2y φzz − 2φy φz φyz + φ2z φyy )/ |Oφ|3. 2.1.3. Signed distance functions. A general distance function is defined as d (~x) = min (|~x − ~xI |). 6. for all ~xI ∈ Γ,. (2.9).

(20) 2. BACKGROUND. 2.2 Level set methods. which means that d (~x) gives the Euclidean distance to the closest point ~xI on the interface Γ. A signed distance function is a subset of implicit functions. It can be expressed as |φ (~x)| = d (~x). It has all the properties of a implicit function plus a couple of new that are unique for signed distance functions. One property is the Eikonal equation, kOφk = 1.. (2.10). It assures that the distances are uniformly distributed and that the derivative is continuous, which is needed in order to have a gradient perpendicular to the interface. Another property is the closest point tranform, a function that takes a point ~x in space and returns the closest point on the interface. ~ ~xC = ~x − φ (~x) N. 2.2. (2.11). Level set methods. Most geometry used in Computer Graphics and Computational Physics is represented as surfaces. The most common surface representation is the polygonal mesh, an explicit representation where points on the surface are stored and connected to form an enclosed surface. Another type of surface representation is the implicit representation, where an implicit function is sampled on points in space. Level set methods enables implicit surfaces, defined as signed distance functions, to be time dependent and dynamic. Level set methods was introduced by Osher and Sethian in 1988 [26], where they used it to track propagating interfaces in time dependent physics problems. To solve the timedependent equations that expresses a moving implicit surface, level set methods use a Hamilton-Jacobi approach. It has since then been improved with higher order spatial and temporal integration schemes and memory efficient data structures. For a thorough review of Level set methods, consult the books by Osher and Fedkiw [25] and Sethian [32]. The motion of the zero level set can be defined by a speed function F (x, ~n, φ, ...) = ~n ·. dx Oφ dx = · dt kOφk dt. (2.12). d~ x dt. is projected on the normal. where ~x is a point on the interface and the partial derivative ~n.. 2.2.1. Discretization. When working with spatially defined functions in a discrete environment such as a computer, the values has to be sampled. There are two general representations that can be used, the Lagrangian representation and the Eulerian representation.. 7.

(21) 2.2 Level set methods. 2. BACKGROUND. Lagrangian representation In the Lagrangian representation, a number of points are tracked in space that each store properties. The advantages of this representation are that the whole computational domain does not need to be stored and it can store the position of each point with very high precision. Drawbacks with this representation are that it can be hard to compute geometrical properties and if the points are connected, it can be hard to manage topology changes. Particle simulations and polygonal meshes uses the Lagrangian representation.. y. y. x. x. (a) Lagrangian. (b) Eulerian. Figure 2.2: Discrete representations in two dimensions.. Eulerian representation The Eulerian representation stores information at points with uniform distance between eachother, forming a regular grid. The distance between each point decides how accurate the sampling will be. The distance is denoted ∆x, ∆y, ∆z and to simplify the implementation the distance can be the same for all dimensions ∆x = ∆y = ∆z, forming a Cartesian grid. Knowing the distance between all the points makes it easy and fast to calculate geometric properties such as gradient and curvature. This is one main reason why signed distance functions are best represented using the Eulerian representation. There are some problems with this representation that has to be considered. It suffers from high numerical dissipation and requires large amounts of memory to store. Interpolation When a signed distance function are sampled on a discrete Eulerian grid, function values can only be accessed on grid points. In order to obtain a function value that is located within a grid cell, interpolation has to be used. The most basic interpolation method is called Nearest neighbour and as the name suggests the value of the nearest grid point is used. This method does not give good results though, in order to achieve decent results more sophisticated methods are needed. A common and straightforward method is called trilinear interpolation, where three linear interpolations are performed, one for each dimension. There are also more advanced interpolation methods, such as cubic. 8.

(22) 2. BACKGROUND. 2.2 Level set methods. interpolation and hermite interpolation, but the accuracy these methods provide are not needed for this application. Trilinear interpolation Trilinear interpolation is as described earlier three linear interpolations, one for each dimension. It will interpolate the value at position x, y, z using the eight surrounding grid values. The grid cell i, j, k which is used as base for the interpolation is determined by. Ci,j+1,k+1. Ci+1,j+1,k+1. Ci,j,k+1. Ci+1,j,k+1. x,y,z. Ci,j+1,k. Ci+1,j+1,k. Ci,j,k. Ci+1,j,k. Figure 2.3: Trilinear interpolation flooring the position x, y, z, and these two coordinates are used to calculate the interpolation weights i = bxc α = x − i j = byc β = y − j (2.13) k = bzc γ = z − k The first step is to perform linear interpolation in the x dimension, reducing 8 values to 4. X00 = (1 − α) Ci,j,k + αCi+1,j,k X10 = (1 − α) Ci,j+1,k + αCi+1,j+1,k (2.14) X01 = (1 − α) Ci,j,k+1 + αCi+1,j,k+1 X11 = (1 − α) Ci,j+1,k+1 + αCi+1,j+1,k+1 The second step is the y dimension, where 4 values are linearly interpolated to 2. Y0 = (1 − β) X00 + βX10 Y1 = (1 − β) X01 + βX11. (2.15). The third and final dimension where 2 values are weighted together to one value. Z = (1 − γ) Y0 + γY1. 9. (2.16).

(23) 2.2 Level set methods. 2.2.2. 2. BACKGROUND. Advection. The interface of a signed distance function can be evolved over time using a external ~ . The simplest way to describe the motion of an evolving surface is the velocity field V Lagrangian formulation, where points on the surface are moved using the corresponding velocity vector in the velocity field. This formulation is expressed as the ordinary differential equation d~x ~ (~x) =V (2.17) dt for all points ~x with φ (~x) = 0. This means all points ~x on the interface will be moved ~ (~x). The problem with this formulation is that if the interface is according to the vector V using a Lagrangian represention, topological changes are hard to handle and geometrical properties are non-trivial to approximate. Therefore an Eulerian formulation is better suited for advection of a signed distance function. The Eulerian formulation is a partial differential equation (PDE) that uses φ to both represent and evolve the interface. ∂φ ~ =0 + Oφ · V ∂t. (2.18). ∂φ + kOφk F (x, ~n, φ, ...) = 0 ∂t. (2.19). where O is the gradient operator. The gradient is expressed as ~ =u Oφ · V ˆ. ∂φ ∂φ ∂φ + vˆ +w ˆ ∂x ∂y ∂z. (2.20). Equation 2.18 and 2.19 are commonly known as the Level set equations. For the Lagrangian ~ only needs to be defined at the interface but for the formulation 2.17 the velocity field V Eulerian formulation it has to be defined for the whole computational domain Ω. In order to solve the Eulerian PDE 2.18, the spatial and temporal partial derivatives of φ must be numerically approximated. Solutions to this will be described in the to following sections.. 2.2.3. Spatial discretization. The gradient is one of the most important properties of a implicit function. Analytically, it is calculated by solving a partial differental equation for each dimension. But in the discrete case, the signed distance function is sampled on a Cartesian grid, where a numerical scheme is needed to approximate the derivative. A collection of numerical schemes that is well suited for this is the finite difference (FD) schemes. There are different FD schemes that gives different accuracy on the result, all with different computational requirements and implementation complexity. The schemes presented here will only be shown for the x dimension, the scheme is applied in the same way for y and z. Upwind differencing The simplest scheme is the first order accurate forward and backward upwind difference. These schemes are derived by using the first two terms of the Taylor expansion. The. 10.

(24) 2. BACKGROUND. 2.2 Level set methods. forward scheme uses the value at the current grid point i and at the next grid point i + 1 to approximate the derivative. φi+1 − φi ∂φ ≈ Dx+ φ = ∂x ∆x. (2.21). ∆x denotes the grid spacing in the current dimension. The backward scheme is almost the same as the forward scheme, but uses the value at the current grid point i and at the previous grid point i − 1 instead. ∂φ φi − φi−1 ≈ Dx− φ = ∂x ∆x. (2.22). Central difference A slightly more accurate scheme is the second order accurate central difference. It uses values at grid points on both sides of the current grid point i. This scheme is actually derived as a combination of a forward and a backward upwind step. φi+1 − φi−1 ∂φ ≈ Dx φ = ∂x 2∆x. (2.23). The problem with second order accurate central differencing is that the value at the central grid point is not evaluated. If this point should contain a value very different from its neighbors, the approximated derivative would be wrong. There is also a second order accurate central differencing scheme for approximation of the second order partial derivative. φi+1 − 2φi + φi−1 ∂2φ ≈ Dx+ Dx− φ = (2.24) 2 ∂x ∆x2 Higher order schemes The numerical solutions of upwind differencing can be consideraby improved by replacing the first order upwind schemes with higher order schemes. Two schemes that has proven to give good numerical solutions for level set methods are the Hamilton-Jacobi methods ENO (essentially nonoscillatory) [27] and WENO (weighted essentially nonoscillatory) [14]. Hamilton-Jacobi ENO is a third order accurate upwind scheme that approximates a Newton polynomial and its derivative. Hamilton-Jacobi WENO is a third to fifth order accurate upwind scheme that approximates three ENO polynomials and weights them together. The schemes are described in detail in Appendix A and B. Numerical stability Osher and Fedkiw [25] mentiones thatThe Level set equations (2.18 and 2.19) has two important properties, it is a hyperbolic PDE and it is a Hamilton-Jacobi equation. Being a hyperbolic PDE means the function does not depend on derivatives of order greater than one. It also gives a special behaivor where they propagate information along certain directions. Because of this, special care has to be applied when approximating spatial derivatives. To approximate the partial derivative φnx , upwind differencing has to be. 11.

(25) 2.2 Level set methods. 2. BACKGROUND. performed in the reverse direction of the velocity field flow. This can be divided into the following cases  + n  φx , V < 0 n φnx = φ− (2.25) x, V > 0   0, Vn =0. 2.2.4. Temporal discretization. In order to evolve the zero level set over time, the temporal derivative has to be approximated and used to advance the interface. There are a number of different schemes for integration over time and all schemes can be classified as one of two types, implicit or explicit. Implicit methods are unconditionally stable but are generally hard to implement. Explicit schemes on the other hand are often straightforward to implement but has stability issues. Forward Euler The simplest time integration scheme is forward Euler. It is an explicit scheme and is only first order accurate. The reason it is popular is because it is computationally fast and easy to implement. The general formulation is derived from the first two terms of an Taylor expansion. yn+1 = yn + hf (tn , yn ) (2.26) Forward Euler uses the current value yn and adds one timestep h times the function f (tn , yn )to reach the next value yn+1 . The Forward Euler formulation can be applied to the Level set equation, which gives the following PDE for integrating the zero level set in time. φn+1 − φn + Oφn · ~un = 0 ∆t. (2.27). TVD Runge-Kutta Total Variation Diminishing (TVD) is a property of certain temporal integration schemes. The total variation of a discrete quantity u can be calculated as T V (u) =. X. |uj+1 − uj |. (2.28). j. and a method is considered to TVD if the following condition is fullfilled. T V (un+1 ) ≤ T V (un ). (2.29). Osher and Shu [33] introduced Runge-Kutta time integrators that had the TVD property. They used the fact that convex combinations are TVD for non-negative coefficients and Runge-Kutta integrators are convex combinations of Forward Euler steps. The first order TVD RK is identical to forward Euler, presented in section 2.2.4.. 12.

(26) 2. BACKGROUND. 2.2 Level set methods. TVD RK2 The second order TVD Runge-Kutta is a convex combination of two Forward Euler steps. The method is also known as the midpoint rule or modified Euler and is an explicit scheme. Begin by taking one Euler step to tn + ∆t: φn+1 − φn + Oφn · ~un = 0 ∆t. (2.30). which is used for the second Euler step to tn + 2∆t: φn+2 − φn+1 + Oφn+1 · ~un+1 = 0 ∆t. (2.31). Finally an averaging step that takes the convex combination of the initial data and the result of 2.31, gives the result for time tn + ∆t: 1 1 φn+1 = φn + φn+2 2 2. (2.32). TVD RK3 The third order TVD Runge-Kutta is proposed to be used for level set time integration by Enright et al. [5]. The scheme can be broken down into three Forward Euler steps and two convex combinations. Begin with a Forward Euler step to tn + ∆t: φn+1 − φn + Oφn · ~un = 0 ∆t. (2.33). which is used for a second Euler step to tn + 2∆t: φn+2 − φn+1 + Oφn+1 · ~un+1 = 0 ∆t. (2.34). 1 3 1 φn+ 2 = φn + φn+2 4 4. (2.35). Average to tn + 12 ∆t using:. Take a third Euler step to tn + 23 ∆t 3. 1. 1 1 φn+ 2 − φn+ 2 + Oφn+ 2 · ~un+ 2 = 0 ∆t. (2.36). followed by a second averaging step: 3 1 2 φn+1 = φn + φn+ 2 3 3. This produces a third order accurate solution.. 13. (2.37).

(27) 2.2 Level set methods. 2. BACKGROUND. CFL condition The CFL (Courant-Friedreichs-Lewy) condition [4] is a condition to make sure that time integration of the Level set equation 2.18 is numerically stable. By dividing the cell width with the greatest velocity field component, the calculated timestep ∆t will make sure the zero level set will only be moved at most one cell width. ∆t <. ∆x max |~u (~xi,j,k )|. (2.38). This contition can be enforced by introducing a CFL number α. The multidimensional CFL condition can then be written as   ~ (~xi,j,k )| |~u (~xi,j,k )| |~v (~xi,j,k )| |w =α (2.39) ∆t max ∆x ∆y ∆z where α is a value ∈ (0, 1). α is according to Osher and Fedkiw [25], commonly selected as a value in the range [0.5, 0.9].. 2.2.5. Renormalization. According to [28], in order to maintain a signed distance function it needs to be renormalized every time the interface has been evolved over time. This is done by solving the initial value formulation of the eikonal equation:  ∂φ = sgn φ0 (1 − kOφk) ∂t. (2.40). Osher and Fedkiw [25] points out that the sign can be replaced by a smoothing function S (φ) that improves the numerical solution. φ S (φ) = q φ2 + kOφk2 (∆x)2. (2.41). A way of calculating the norm of the gradient vector kOφk was presented in [31], where they used Godunov’s method for solving PDEs. The square of the partial derivative is given by    − , 0)2 , min (φ+ , 0)2 , if F > 0  max max (φ  x x    2 + , 0)2 , if F < 0 φ2x = max min (φ− (2.42) , 0) , max (φ x x    0 F =0 and the total norm is expressed as kOφk =. q φ2x + φ2y + φ2z. (2.43). Solving equation 2.40 for the whole signed distance function can be very computation heavy. Therefore, Sethian introduced the Fast Marching Method to efficiently solve the equation. This method is described in detail in the book Level Set Methods and Fast Marching Methods [32].. 14.

(28) 2. BACKGROUND. 2.2.6. 2.2 Level set methods. Optimizations. Since the level set methods was introduced, recent research has not only improved the discretization schemes, there has also been work done trying to reduce the amount of data needed to be computed and efficient ways of storing the significant parts of the signed distance function. Narrow band method The narrow band technique is an optimization to the level set methods where computations are only performed within a distance δ from the interface. The band between −δ and δ is called the narrow band. The condition can be expressed as abs (φi,j,k ) < δ. (2.44). The method was later improved by Peng et. al. [28] where they divide the narrow band into several bands. β defined the computational band and γ defined a extended band that was used as a smooth transition to reduce oscilliations at the boundaries. They used the following cutoff function to achieve the smooth transition.. c (φ) =.    1. if |φ| ≤ β.  0. otherwise. (|φ|−γ)2 (2|φ|−3β+γ) (γ−β)3 . if β < |φ| ≤ γ. (2.45). Efficient data structures One of the biggest problem when working with signed distance functions is the large memory footprint. Sampling a signed distance function on a 5123 grid will require at least 512 MB of memory. For each operation 134 million elements has to be processed. To overcome this a couple of data structures has been developed, that takes advantage of the nature of the signed distance function. One approach of storing signed distance functions more efficiently was to use an Octree [17]. This make it possible to store data with a varying detail level, depending on the relevance of the data. Data near the interface could have higher resolution than data far away from it. One drawback with this approach is that it is complicated to implement and it takes O (log (n)) time to random access a value. A different approach than the Octree was the DT-Grid (Dynamic tubular grid) by Nielsen and Museth [23]. It is a hierarchical compression technique that recursively compresses each dimension using a similar technique to what is used to compress sparse matrices. This technique only stores the signed distance function within the narrow band, reducing both memory consumption and the amount of elements that has to be processed. Houston et. al. [12] combined the DT-Grid concept with the compression technique Run Length Encoding to form the Hierarchical Run Length Encoding method.. 15.

(29) 2.3 Particle Level set method. 2.3. 2. BACKGROUND. Particle Level set method. The Particle Level set method is an extension to Level set methods. It uses marker particles as a complementary to the Eulerian grid-based signed distance function. The particles are used to correct the signed distance function when its solution is weak. The Lagrangian approach tracks the interface more accurate and it makes it possible to capture details that are smaller than what is possible with only a grid-based solution. This higher accuracy makes it possible to preserve the volume of the level set better when is gets deformed. The concept of using marker particles was first introduced by Foster and Fedkiw [8]. They used particles on the inside of the interface to correct the signed distance function when its solution was weak. The method was improved by Enright et. al. [5] where they used particles on both sides of the interface. The method was further improved by Enright et. al. in [6], where they introduced Semi-Lagrangian advection and fast marching techniques to improve efficiency.. 2.3.1. Overview. The Particle Level set method can be divided into several well-defined steps: • Initialize particles around the zero level set. • For each timestep ∆t : 1. Advect the zero level set one timestep. 2. Advect all particles one timestep. 3. Correct the signed distance function using escaped particles. 4. Renormalize the signed distance function. 5. Correct again with escaped particles. 6. Adjust the radius of each particle. 7. Add and remove particles where needed.. 2.3.2. Massless particles. To track the interface, particles are used as a complementary to the signed distance function. Because particles are only used for tracking, there is no need to store any mass. Each particle is represented with a position and a radius. The radius vary between the particles and is used to store the distance to the zero level set. This means that the particles represent a sampled Lagrangian version of the signed distance function. The radius is limited according to rmin = 0.1 max (∆x, ∆y, ∆z) (2.46) rmax = 0.5 max (∆x, ∆y, ∆z) Particles are classified as positive or negative, where all positive particles represent the outside of the interface and the negative represent the inside. It is also necessary to keep track of the particles that are considered as escaped.. 16.

(30) 2. BACKGROUND. 2.3.3. 2.3 Particle Level set method. Initialization. Initialization is done by inserting a specified number of particles in the cells near the zero level set. This is done for all cells that has a corner within 3∆x of the zero level set. Enright et. al. [5] suggests inserting 32 positive and 32 negative particles into each cell. Each particle is randomly placed within the cell and an attraction step is performed to attract it to a band in its corresponding side of the interface. The band is defined as bmin ≤ |φ| ≤ bmax , where bmin = rmin and bmax = 3∆x. The attraction step requires a target, in this case φgoal . The value is randomly selected within the range [bmin , bmax ]. It steps the position ~xp closer to the band according to ~ (~x) . ~xnew = ~xp + λ (φgoal − φ (~x)) N. (2.47). λ = 1 when the attraction starts and for each iteration λ is halved. This procedure is performed until φ (~xp ) is within the band. When a particle is correctly placed in its corresponding band, its radius is set to the distance to the zero level set. The distance is given by interpolating the φ value at the position of the particle. To make sure the radius does not get to large or small, it is fitted to the correct range using   if sp φ(~x) > rmax rmax rp = sp φ(~x) if rmin ≤ sp φ(~x ) ≤ rmax (2.48)   rmin if sp φ(~x) < rmin . where sp is the sign of the particle.. 2.3.4. Advection. The zero level set and the particles are advected independently. Particles are advected using the Lagrangian formulation described in section 2.18. In [5] they suggest using the third order TVD Runge-Kutta integration scheme (section 2.2.4) for both level set and particle advection. In [6] it is proposed that the semi-Lagrangian method for level set advection and second order TVD Runge-Kutta integration for particles will give equally good results with faster computational time. Semi-Lagrangian advection The semi-Lagrangian advection method is a first order accurate advection scheme. The scheme handles both spatial and temporal integration. Spatial integration is done with linear interpolation according to n φn+1 i,j,k = αβγφi+1,j+1,k+1 + (1 − α) βγφni,j+1,k+1 + α (1 − β) γφni+1,j,k+1 + αβ (1 − γ) φni+1,j+1,k + (1 − α) (1 − β) γφni,j,k+1 + (1 − α) β (1 − γ) φni,j+1,k + α (1 − β) (1 − γ) φni+1,j,k + (1 − α) (1 − β) (1 − γ) φni,j,k. 17. (2.49).

(31) 2.3 Particle Level set method. 2. BACKGROUND. The interpolation weights are calculated by tracing a Lagrangian particle backwards in time in the velocity field flow as expressed in   (i − r) ∆x − ui,j,k ∆t ∆t , α= r = i − ui,j,k (2.50) ∆x ∆x   (j − s) ∆y − ui,j,k ∆t ∆t s = j − ui,j,k , β= (2.51) ∆y ∆y   (k − t) ∆z − ui,j,k ∆t ∆t , γ= t = k − ui,j,k (2.52) ∆z ∆z Enright et. al. [6] states that the method is unconditionally stable because of the linear interpolation and the Lax-Richtmyer theorem gives that it will converge to a correct solution. Being unconditinally stable means the timestep is not bound by the CFL condition.. 2.3.5. Error correction. After both the zero level set and the particles have been advected the particles are used to correct the level set interface where its solution is weak. Because the particles has been advected independent of each other, the particles carry a more accurate solution of the advection. Therefore only the particles that does not match the zero level set are used in the correction. This is done by first determining all particles that is considered escaped. Escaped means the particle has appeared on the wrong side of the interface by more than its radius. For positive particles, this means the interpolated φ value at its center is ≤ 0 and for negative particles it means φ value > 0. For each escaped particle a local distance function φp is calculated at each of the eight cell corners surrounding the particle. φp (~x) = sp (rp − |~x − x~p |). (2.53). where sp is the sign of the particle, ~x is the position of the cell corner, rp is particle radius and ~xp is particle position. If the value φp deviates from φ there is a possiblility of an error in the level set. All the positive escaped particles are used to correct the φ > 0 region of the signed distance function and all the negative escaped particles are used for the φ ≤ 0 region. The calculated value φp is compared to the value stored at the cell corner and the maximum is used. For positive particles this can be expressed as  (2.54) φ+ = max φp , φ+ ∀p∈E +. and for negative particles it can be expressed as φ− = min. ∀p∈E −. φp , φ−. . (2.55). It is not certain that φ+ and φ− will agree. They are merged to a final value with ( φ+ if |φ+ | ≤ |φ− | φ= (2.56) φ− if |φ+ | > |φ− |. 18.

(32) 2. BACKGROUND. 2.3.6. 2.4 Sparse volume data structures. Radius adjustment. When the signed distance function has been corrected, all particles must be adjusted to fit the new location in respect to the interface. This is done by setting the radius of each particle to its current distance to the zero level set, φ (~xp ). It is done in the same way as when the particles were initialized, by clamping the distance value using equation 2.48. There are some exceptions, when a particle is closer to the interface than rmin . In that case, the radius is unchanged. Particles that remain escaped after the correction has its radius set to ±rmin , depending on what side of the interface they are located.. 2.3.7. Reseeding. When the zero level set evolves over time, the shape and size of the surface will change. In order to maintain the desired number of particles in each cell, reseeding has to be performed. Each cell that contains particles has to be examined to determine if there is any need to add or delete particles. If there are to few particles, a positive and a negative particle is added until the right amount is achieved. The particles are added in the same way as in the initialization step (described in section 2.3.3). If there are to many particles in a cell, all particles in that cell that are not escaped are sorted based on their distance to the interface. Then the particle with the greatest distance are deleted until the cell contains the right amount of particles. It can be costly to examine each cell for every timestep of advection. The simplest solution to avoid this is to only perform reseeding after a fixed number of timesteps. Another more sophisticated solution is to measure how much the surface area has changed over time. If the difference is greater than a predefined threshold, reseeding is performed. The surface area of the zero level set is given by Z δ (φ) |Oφ| d~x (2.57) where δ (φ) is a numerically smeared out delta function. With first order accuracy it can be expressed as    0   φ < − 1 1 + 2 cos πφ − ≤ φ ≤  (2.58) δ (φ) = 2    0 <φ where  = 1.5∆x.. 2.4. Sparse volume data structures. In section 2.2.6 the need for efficient data structures to store signed distance functions is described. A couple of data structures dedicated to signed distance functions are briefly presented. In recent years research at two Visual Effects studios has resulted in two new general purpose volume data structures, DB-Grid and Field3D. Both these data structures are based on a block concept which provides sparse storage and should theoretically be more efficient in terms of speed and memory consumption compared to previous data structures.. 19.

(33) 2.4 Sparse volume data structures. 2.4.1. 2. BACKGROUND. DB-Grid. DB-Grid (Dynamic Blocked Grid) is a highly efficient volume data structure, developed by Ken Museth at Digital Domain. The data structure is briefly described in [21] and [22]. It is specifically designed to handle sparse volume data. Therefore it is perfectly suited for representing signed distance functions, where for most operations the only useful information needed is located at and within a close distance to the interface. Only storing information on and around the interface instead of the whole domain considerably reduces the amount of data that is needed to be processed and stored. Fast access to data is also important when working with signed distance fields. DB-Grid provides both random access in constant time and sequential access in linear time, which is highly desirable. Design concept The data structure divides the grid domain into equally sized blocks and all blocks have the same size in all dimensions. The block size is a number power of 2, usally 8 or 16. Blocks are represented as its own block object and stores an array that constains its elements. The element array can be allocated as a two dimensional array, where the second dimension is used to access different layers of data. Blocks are accessed through the block array, an array that stores pointers to all blocks. The blocks are classified as either • Dirty block • Outside block • Inside block where the data structure contains multiple Dirty block, one Outside block and one optional Inside block. The Dirty blocks are the blocks that stores explicitly set data. The Outside block is the key to make the data structure sparse. The elements of this block are all set to a predefined fill value. When the data structure is initialized, all the pointers in the block array will point at the Outside block. This means that everywhere in the grid domain the fill value will be returned, but there will only be one block of data allocated. When data is set explicitly, block objects classified as Dirty blocks will be allocated. A nice benefit of allocating data this way is that the data will be cache coherent. As long as a block is containing explicitly set data, the block will remain in memory. As soon as the data is removed, the block object will be deallocated and the block pointer will point at the Outside block again. Bit operators are used to quickly convert Cartesian grid coordinates to a block coordinate and a element coordinate, by transforming three dimensional coordinates into one dimensional linear offsets. Both the block coordinate and the element coordinate is a linear offset into each corresponding array. Only two dereferences and the coordinate conversion is needed for random access of a data value.. 20.

(34) 2. BACKGROUND. 2.4 Sparse volume data structures. 0. 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11. 12. 13. 14. 15. Figure 2.4: Shows a two dimensional example of a implicit sphere represented using the DB-Grid concept. Green blocks represent Dirty blocks, blocks outside the interface are Outside blocks and the block inside is the Inside block.. Block array 0. 1. 2. 3. Block 0. Outside block. 4. 5. Block 1. 6. 7. Block 2. 8. 9. Block 3. 10. 11. Block 4. 12. 13. Block 5. 14. 15. Block 6. Block 7. Inside block. Figure 2.5: This scheme shows the block array and how it points to the different blocks. Fast iterators The most common way of accessing data when working with signed distance functions is sequential access. To provide this in a convinient and efficient way, iterators are used. These iterators have a sparse behavior, which means they do not step throught all elements in the grid, only where explicit set data are stored. To keep track of explicit data, bit masks are used. One bit mask for the block array and one for the element array in each block object. This enables sequential access in linear time. When travering the grid using a sparse iterator it also provides a convinient method that will return the grid coordinate at the current position. The grid coordinate can easily be transformed into world coordinates using a transformation matrix provided by the DB-Grid. Level set customization This concludes the general purpose volume data structure. When used for signed distance functions a small extension is added. To complement the Outside block, the Inside block is used. It represents the inside of the zero level set and contains the same fill value as the outside block, but with a negative sign. In order to correctly assign inside and outside, a sign correction algorithm has to be run.. 21.

(35) 2.4 Sparse volume data structures. 2. BACKGROUND. Earlier in this section, different layers of data was described. When working with level sets this is a very useful feature. In can be used to store the level set solution at different time t when performing time integration.. Figure 2.6: The Stanford Dragon represented as a signed distance function using DB-Grid. Green blocks show the Dirty blocks and the yellow box represents the whole grid domain.. 2.4.2. Field3D. Field3D was developed at Sony Pictures Imageworks and it was made open source in 2009 [35]. It is a complete volume data API developed in C++ that can handle both dense and sparse volume data. It also supports the file format HDF5, a format developed by Nasa for storage of volume data. In this thesis we only focus on the SparseField of Field3D. This data structure uses the same concept as DB-Grid, by diving the grid domain into equally sized blocks. Each block is represented as a block object and each block object contains an element array and a single value. When the data structure is initialized, the domain is divided into blocks and a block object will be allocated for each. However, the element array will not be allocated. Instead only the single value will be set that represents the whole subdomain of the block. It is not until the user explicitly sets data into the block the element array will be allocated. This gives the data structure the ability to use two levels of detail within the grid domain. It provides a sparse behavior which is very well suited for representing signed distance functions where the interesting data is on and close to the interface. The Field3D API provides iterators which makes it easy to traverse all the data and transformation to transform grid coordinates to world coordinates etc.. 22.

(36) Chapter 3. Implementation This chapter will describe how a Particle Level set solution was implemented. It is based on the theory described in chapter 2. The implementation will use the advantages of the sparse volume data structures described in 2.4, to both store the signed distance function and the particles. There will also be suggested a couple of particle implementations, with varying accuracy and memory requirements.. 3.1. Current implementation. The current implementation of the Particle Level set method is an integral and important part of Digital Domains Award winning fluid simulation package fsim. Implemented in the mid 2000 it does not feature any efficient data structures. It is considered one of the slowest parts of the fluid simulation along with the Projection solver. The signed distance function is stored in a dense three dimensional grid, where the grid is implemented as a linear array. It is extended with the narrow band solution where a mask is used to keep track of the interface and its surroundings. Particles are stored in two arrays, one array with positive particles and one with negative particles. Each particle object contains • Position in world coordinates, 3 x float • Radius, 1 x float • Escaped, 1 x boolean The particles are mapped to its corresponding grid cell using a special mapper class. This is computational heavy because every time the particles are moved, it has to remap all the moved particles.. 3.2. Storing the particles. This section will describe the implementation of the particle storage. An efficient way to store and locate particles is a very important part of the Particle Level set implementation. There was a number of guidelines from Digital Domain regarding the implementation of the particle container.. 23.

(37) 3.2 Storing the particles. 3. IMPLEMENTATION. • Based on a sparse volume data structure • General purpose • Cache coherent • Minimize memory allocation • Support parallellism • Query particles from cell • Query cell from particle • Particle quantization Two different approaches of storing particles will be described, one uses an array for each element in the sparse volume and the other uses a single array and all elements points to its area in the array. We came up with the concept of the multiple arrays approach, the other approach was inspired by the solution of Nielsen and Museth [24]. Because of time limitations only the approach with multiple arrays was implemented and tested.. 3.2.1. Class Interface. The particle container was implemented as a general purpose framework that can be used for any problems that deals with particles. Therefore it is necessary to have a logical and simple interface. The following methods are suggested • Add particle(cell coordinate, particle) • Move particle(cell coordinate from, cell coordinate to) • Remove particle(cell coordinate, particle id) • Clear particles(cell coordinate) • Get Particle Count(cell coordinate) • Get Total Particle Count • Clear field One thing that is common for both approaches is that position of each particle is stored in local cell coordinates. The cell coordinate space is defined as [0, 1) for each dimension. This means that the particle itself does not know where in the grid it is located. When a particle gets advected it can be moved at max one cell according to the CFL condition 2.2.4. This means that the particle can have a position in the range [−1, 2) before it gets moved to its new cell.. 24.

(38) 3. IMPLEMENTATION. 3.2.2. 3.2 Storing the particles. Particle pool. The particle pool approach is similar to the particle bin method used in [24]. All particles are stored in a big array and each element in the volume data structure stores a linear offset and the number of particles the element contains. A predefined value determines the maximum number of particles an element can have. The particle array allocates as many position as there are allocated elements in the grid times the maximum number of particles per element. When accessing particles for a particular element, particle iterators will be returned that only accesses the particles that belongs to that element. The start iterator will point at the particle that is located where the linear offset points in the particle array and the end iterator will point at the linear offset added by the number of particles in that cell. This implementation can execute in parallell by diving the particle array into intervals and let each thread work on one interval. Pros • Allocated memory is coherent. • Particles can be traversed without travering the grid. Cons • Complex implementation. • The particle array needs to be resized each time particles are added to a new element or all particles are moved from an element. • Storing an offset and size for each element requires extra memory.. Grid elements. Particle array. Offset: 0. Offset: 5. Offset: 9. Size: 5. Size: 4. Size: 6. 0. 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. Figure 3.1: Particle pool structure.. 25. 11. 12. 13. 14.

(39) 3.2 Storing the particles. 3. IMPLEMENTATION. Grid elements. Particle arrays. 0. 0. 0. 1. 1. 1. 2. 2. 2. 3. 3. 3 4. 4. 5. Figure 3.2: Multiple arrays structure.. 3.2.3. Multiple arrays. This approach uses one array per element to store the particles that are geometrically located within its grid cell. To avoid allocating arrays for each element of a block, a pointer is stored. If particles are added to a particular element, an array is allocated and the particles are inserted. Each element of a block stores a boost shared pointer to an array and the array stores the particle objects. boost is a utility library for C++ and the shared pointer is a pointer implementation that provides reference counting which eases memory management. If an object created using the shared pointer is not referenced anywhere in the implementation, it will automatically be deleted. This implementation can execute in parallell by letting each thread work on different arrays or different blocks of arrays. Pros • Easy to add particles to new an element or remove all particles from an element. • Straightforward to implement. • The arrays can be adjusted to fit the cache memory in a threaded implementation. Cons • Arrays could be allocated at different locations in the memory. • Small memory overhead for each array that is created. • In order to traverse particles, the grid needs to be traversed.. 26.

(40) 3. IMPLEMENTATION. 3.2.4. 3.2 Storing the particles. STL containers. In both particle storage approaches, an array is needed to store the particle objects. The C++ Standard Library (STL) provides a couple of container classes, all with different advantages. The classes are all templated, which means any type of object can be stored. All the classes provides methods for safe access, add object to front, add object to back and get size. They also provide iterators for simple traversal of the data. Vector The Vector [3] is the STL implementation of a standard linear array. When a Vector object is created it internally allocates memory for a predefined number of the templated object. If more objects are added than the Vector has pre-allocated for, it has to internally create a new array with enough space. The stored objects then has to be copied to the new array and the old array is deleted. This makes the Vector useful in a case where a large number of objects are to be stored and the number of objects does not change over time. The allocated data is memory coherent. List The List [2] is a doubly linked list. Each element is stored as a node and these nodes are linked forward and backward, forming a chain. It starts out with one node and for each new element that is added, a new node is created and linked to chain. Elements can not be accessed directly by their index number, instead the list has to iterate to the desired index. The List is useful when the amount of objects to store is constantly changing and objects needs to be inserted anywhere in the list. Because new nodes are allocated anywhere in memory, the allocated memory will not be coherent. It can move and remove nodes fast by only changing the pointers that link the nodes together. Deque Deque is acronym for Double-Ended Queue [1]. The Deque internally allocates data by grouping elements into blocks and each time more elements needs to be allocated, a new block of elements is added. The blocks are linked by pointers. This is similar to a List, but in this case each node is an array of subelements. Nothing of this is exposed to the programmer, giving it the impression of a regular array. The Deque is useful when storing a small number of objects and you want the ability to quickly allocate space for a couple of more objects. It will provide memory coherency because each block of elements is allocated together in memory.. 27.

(41) 3.3 Particle representation. 3.2.5. 3. IMPLEMENTATION. Traversing particles. Most of the operations performed on particles in the Particle Level set method applies to all particles. Therefore it is convinient to have a simple way of traversing all particles. Both volume data structures provides sparse iterators that will allow us to iterate through only grid elements that contains particles. They also has an easy method of finding out the grid coordinate (i, j, k) of the current element. The grid coordinate can be transformed to world space coordinates (I, J, K) using a grid-to-world transform provided by the volume data structure. The world space grid coordinate can the be used to compute the world space position of the particle according to ~x = I~ + ~xp ∆x. (3.1). where ~xp is the particle position in local space.. 3.3. Particle representation. One of the bottlenecks when implementing the Particle Level set method is the amount of particles needed. The surface area can grow very fast when the zero level set is evolved over time, which means the particle amount increases considerably. One thing we have explored in this thesis is to implement the particle representation with different accuracy and see what concequences it has on the resulting signed distance function. As a minimum, a particle must be represented with a position, a radius and it must keep track if it has escaped and if it is a positive or negative particle. In all three implementations the radius is used to keep track of if it is a positive or negative particle. This is easily done by letting the radius be positive for positive particles and negative for negative particles. It is important that each particle implementation has the same interface. This makes it easy to switch between the different implementations without having to change any other code. Here is a simple and straightforward interface that was used for all the particle implementations • setPosition( float, float, float ) • getPosition() • setRadius( float ) • getRadius() • setEscaped( boolean ) • hasEscaped(). 28.

(42) 3. IMPLEMENTATION. 3.3.1. 3.3 Particle representation. Normal implementation. This is a very simple solution that uses floats to store position and radius. It gives high precision on each property, but it also leads to a large memory footprint. • Position: 3 floats (3 x 32 bits) • Radius: 1 float (1 x 32 bits) • Escaped: 1 boolean (1 x 8 bits) Total memory usage per particle is 17 bytes (136 bits). The implementation is very straightforward and no special considerations has to be done in respect to position and radius range.. 3.3.2. Semi-Compact implementation. The semi-compact solution replaces floats with unsigned short integers which gives half the precision of a float. Also the escaped property is not stored as a boolean. Instead one bit of the radius variable is used to keep track of this property. • Position: 3 unsigned short integers (3 x 16 bits) • Radius: 1 shared unsigned short int (15 bits) • Escaped: 1 shared unsigned short int (1 bit) The total memory usage per particle is 8 bytes (64 bits), less than half of the normal implementation. In this implementation, values has to be stored as integers. This means the float values that is passed in to the particle object has to be quantizied and scaled to fit the integer value. The quantization is done similar to what is suggested in [24]. An unsigned short integer is stored using 16 bits precision, which gives 216 = 65536 levels of quantization. This amount of levels means a coordinate x can take a quantizised value in the range [0, 65536]. As mentioned in section 3.2 all particle coordinates are limited to the range [−1, 2). The coordinate x is mapped from the local space to quatized space xq using this equation   65535 (x + 1) xq = (3.2) 3 where it first translates the value to the range [0, 3), multiplies by the number of quantization levels and divides by the length of the range. It can be mapped back to local space using 3xq x= −1 (3.3) 65535 The radius are mapped in a similar way, but are limited to the the range [−0.5, 0.5]. This range is determined by rmax (equation 2.46), assuming that ∆x = 1 in local space. The sign of the radius determines if the particle is positive or negative. The mapping equation is rq = d65534 (r + 0.5)e (3.4). 29.

(43) 3.3 Particle representation. 3. IMPLEMENTATION. where the radius is translated to [0, 1] and multiplied by the number of quantization levels - 1. This is because the first bit is used for the escaped property. The quantizised radius can be mapped back to local space using rq & 65534 r= − 0.5 (3.5) 65534 where & is the bitwise AND operator. The AND operation will make the equation consider all bits but the one to the far right. The last part of this particle implementation is to be able to set the first bit of the radius variable to 1 or 0 depending on if the particle is escaped or not. To set the bit to 1, meaning the particle has escaped, this equation is used rq = rq |1. (3.6). It uses the OR operator to only modify the first bit. To set the escaped property to 0: rq = rq & 65534. 3.3.3. (3.7). Compact implementation. The compact implementation is very similar to the semi-compact implementation in terms of mapping coordinate and radius and setting the escaped property. The difference is that it uses unsigned char instead of unsigned short integer to store the data. unsigned chars have only 8 bits precision which makes this particle implementation half the size of semicompact and less than a forth of the normal implementation. Each variable will only have 28 = 256 quantization levels. • Position: 3 unsigned char (3 x 8 bits) • Radius: shared unsigned char (7 bits) • Escaped: shared unsigned char (1 bit) The total memory usage per particle is 4 bytes (32 bits). Just as with the semi-compact implementation, the coordinate x has to be translated and scaled in order to be transformed to quantized space.   255 (x + 1) (3.8) xq = 3 The quantized coordinate is mapped back to local space 3xq x= −1 (3.9) 255 The same applies for the radius rq = d254 (r + 0.5)e. (3.10). and when mapping back from quantized space. rq & 254 r= − 0.5 (3.11) 254 where & is the bitwise AND operator. The escaped property bit is set to 1 using the OR operator: rq = rq |1 (3.12) and to set it back to 0: rq = rq & 254 (3.13). 30.

(44) 3. IMPLEMENTATION. 3.4. 3.4 Dynamic seeding condition. Dynamic seeding condition. In order to make the seeding of particles a bit more efficient, we came up with a dynamic seeding scheme where the number of particles to be seeded decreased the furter away from the interface the cell was located. We wanted to see if we could reduce the number of particles without worsen the solution of the signed distance function correction. The following equation was used to determine the number of particles in the current cell nc =. nmax ∆x |φmin |. (3.14). nmax is the maximum number of particles in a cell, ∆x is the grid spacing and φmin is the smallest φ value of the eight that belongs to the current cell. To ensure no division by zero, this was only used where φmin ≥ ∆x.. 3.5. Putting the Particle Level set together. This section describes of the Particle Level set implementation was put together. Some parts are threaded and for the threading the open source package threadpool was used.. 3.5.1. Grids. The signed distance function is stored in a DB-Grid where each allocated element stores a float. A secondary DB-Grid is used to store particles. The particle storage method used is the Particle Array method (section 3.2.3) based on Deques. Why DB-Grid was used in favor of Field3D is discussed in chapter 5. To simplify the implementation, the same grid spacing is used in all dimensions which gives ∆x = ∆y = ∆z.. 3.5.2. Level set advection. The level set advection uses a threaded implementation where each thread is working on a block. Before the execution, a queue of blocks is created for each thread. Spatial derivatives are calculated using WENO (section 2.2.3) and the level set is time integrated using TVD-RK3 (section 2.2.4). Here, the layers in DB-Grid are very useful because in TVD-RK3 three timesteps are calculated. Each can be stored in a different layer and then merged together.. 3.5.3. Level set renomalization. The renormalization is done by solving equation 2.40, where the Godunov scheme (section 2.2.5) is used to calculate the norm of the gradient vector. The time integration is approximated with Forward Euler (section 2.2.4). The timestep used is ∆t = 0.5∆x [25].. 31.

(45) 3.5 Putting the Particle Level set together. 3.5.4. 3. IMPLEMENTATION. Particle advection. Particle advection is one of the parts of the Particle Level set implementation where threading is implemented. Particles are advected using TVD-RK3, where each thread is advecting one block. However it is only the position of the particles that is advected, the particles are not moved to the correct cell. Instead all particles that has a position (~x < 0 or ~x ≥ 1) is moved to the correct cell as a subsequent step. This cannot be done in parallell because of locking, when multiple treads are trying to move to or from the same Deque.. 3.5.5. Radius adjustment and Reseeding. The radius adjustment was also threaded and each thread handled one block. The reseeding is a compution heavy method and according to [5] should not be perform for each timestep. Instead it is performed when the surface area has changed with 5% or 10%. To keep track of the surface area a variable stores the value at initialization. For each timestep the current area is compared to the initial area and if the area has changed too much, reseeding is performed and the initial area variable is reset. To estimate the surface area of the zero level set, we used a crude approximation by counting the number of dirty elements in the grid and dividing this value by the width of the narrowband.. 32.

References

Related documents

Det går att exportera alla respondenterna i en enkät till en lista i Excel, innehållande personuppgifter, kategorier och inloggningsuppgifter för den aktuella enkäten..

ANI (Asian News International), Our policy is to provide safety to all residents in Afghanistan including minorities: Top Taliban leader to Sikhs, 23 December 2021, Factiva,

Information on cases of kidnappings of adult civilians for the purpose of forced recruitment by Taliban groups, in the urban area of Islamabad, between 2012

23 In 2019, UNAMA reported on a general decrease of casualties but noted increase of targeted attacks on certain groups in Afghanistan including on Shia

Between 1 January 2020 and 6 May 2021, ACLED collected 490 violent incidents in the Northwest region of Cameroon, 213 of which were coded as battles, 20 as explosions/remote

79 Reuters, Over 2 million people displaced by conflict in Ethiopia's Tigray region - local official, 12 December 2020, url; BBC News, Tigray crisis: Ethiopia region at risk of

A 25 March 2021 UN report stated that, since 26 February 2021, ‘Palestinian militants’ had ‘launched six rockets and one incendiary balloon from Gaza into Israel’, while the

6 France24, Afghanistan’s media enters the unknown under Taliban rule, 24 August 2021, url; TOLOnews, Afghan Media Activity Faces Sharp Decline: Report, 3 October