Ice Melting Simulation
4.3 GPU-based Ray Tracing With Meshes and Metaballs
To improve the realism, we combine meshes and metaballs in the same scene, and im-plement special ray tracing for these two parts at the same time. We use OpenCL [Khr09]
to operate GPU-based parallel computing ability. The kernel function utilized on GPUs allows multiple threads computing simultaneously. Here we set a thread as a ray, and it is in charge of calculating a color of a pixel on the screen.
OpenCL doesn't support calling the same function in a kernel recursively, so we pro-pose an iterative ray tracing method and use a stack to save newly generated reflection ray.
We use a stack to save the reflection ray separated from the incident ray. The intensity of ray is reduced, and the ray terminates if its intensity is smaller than the threshold. After
terminating the ray, another ray will be popped from the stack and the intersection will be calculated for it.
The color I of pixel including three parts: color of incident ray Ic, color of reflection ray Is, and color of refraction ray It. The relation of these three parts as shown in below:
I = α∗ Ic+ (1− α) ∗ It+ Is, α∈ [0, 1], (4.7) where α is the transparency of the intersected object. When α is closer to 0, the object is more transparent.
Figure 4.9: Sketch of ray tracing.
Figure 4.9 is an illustration of ray tracing. The blue area is a transparent object, and the orange bar is the solid object. The different thicknesses of the red arrow represents different intensity of the ray. The computation of the ray is terminated when the intensity of the ray is smaller than the default intensity threshold.
Furthermore, the flow chart of ray tracing algorithm is shown in Figure 4.10. At first, we compute the direction D of incident ray from the view target. We try to find the intersection between the ray and objects, where objects may be meshes or metaballs.
In order to find the intersection, we use the ray-metaball intersection to find the nearest intersection point Hmt, and we test intersection between ray and BVH construction of
Figure 4.10: Flow chart of ray tracing.
Hmesh and find out the one which is closer to the view target. We have to compute the ray-isosurface intersection test if Hmtis closer. It is worth noting that when we compute the ray-isosurface intersection using the method illustrated in Section 4.2, the algorithm needs to be interrupted once the Hmesh is closer than the absolute position of metaball interval. This mechanism avoids wrong results from occurring as shown in Figure 4.11.
We compute the normal vector n and color Icfor the intersection point, and then we use the following equation to obtain reflection ray Rs:
Rs = D− 2 ∗ (D · n) ∗ n. (4.8)
The ray Rs will be pushed into stack. If the intersected object is transparent, we calculate the refraction angle ϕ using the Snell's refraction law:
C sin θ = Ctsin ϕ, (4.9)
where C is the refractive index of incident medium, Ct is the the refractive index of re-fractive medium, and θ is the angle between n and D. We assume the incident medium is air and the refractive medium is water (or ice). We can use the angle ϕ to obtain refraction ray Rt. Thus, we replace the incident ray D to refraction ray Rtand enter next iteration.
Figure 4.11: The special case occurs when the intersection Hmtis closer to the view target, but the intersection of isosurface is farer than Hmesh.
The ray would be terminated if the intensity of the ray is smaller than the default inten-sity threshold. In this situation, we pop one ray from the stack and enter the next iteration until there is no element in the stack. Finally, we accumulate all the colors computed by the branch rays into the final color of the pixel.
Chapter 5
Experiments and Discussions
In this section, we show the experimental results from our ice melting method and ray tracing method. We implemented the ice melting simulation and ray tracing using GPUs with OpenCL compute capability 1.1. The experiments were performed on an Intel(R) Core(TM) i7-2600 CPU @ 3.40GHz 3.40GHz with 16.0GB RAM, and NVIDIA GeForce GTX 670. Our development environment was Microsoft Windows 7 Home Premium, and Microsoft Visual Studio 2008. The Nvidia GPU Computing SDK version was 4.20.
We used four benchmarks in our experiments, and we will introduce them in next subsection. These four benchmarks are separately used in different experiments. Below we simply illustrate our experiments:
1. We compared the execution time of ice melting for the same model with different resolutions of voxelization. We want to know the impact of the execution time for different resolutions of voxelization on our system.
2. We compared time used for marching cubes construction and ray tracing between the methods proposed in [IUD10] and ours. The realism are also compared. Iwasaki et al. [IUD10] employed marching cubes to construct both ice model and water droplet, and we separated these two parts and used metaballs to represent the water particles. We are interested in the performance and result of these two methods.
3. We continued to compare the realism between the ice melting methods proposed in
[IUD10] and ours. The method proposed in [IUD10] transforms ice particle to water particle when the particle absorbs enough energy. Then, water droplet appears in place. Our method transfers virtual water volume from high to low via ice particles, and generates water particle when the amount of virtual water volume is too large.
4. We compared the rendering calculation timing between the methods with and with-out using BVH construction for metaballs. The paper [KSN08] is good for rendering opaque objects using depth peeling. However, the method suffers bad performance when the metaballs are transparent. We used LBVH to build a hierarchy construc-tion, which speeds up the execution time of ray-metaball intersection test.
5. We compared the rendering results between different heat sources acting on the ice.
We have two types of heat source: Rectangle and Sphere. Our ice melting simulation method employs two types of heat source. The different shape of heat sources will radiate energy to different directions. We are interested in the rendering results for different heat sources.
5.1 Benchmarks
In our experiments, we perform four benchmarks: Cube, Sphere, Bunny and Dragon.
Figure 5.1, 5.2, 5.3, and 5.4 show a series of snapshots for these four benchmarks.
Table 5.1 shows the model complexities, and the number of the initial ice particles that belongs to the model.
Figure 5.1: A series of snapshots for Cube.
Figure 5.2: A series of snapshots for Sphere.
Figure 5.3: A series of snapshots for Bunny.
Figure 5.4: A series of snapshots for Dragon.
Benchmarks # Tri. # Vert. # Ice Particles # Grids of MC
Cube 8 12 13248 32*32*32
Sphere 960 482 7399 32*32*32
Bunny 69668 34837 22795 64*64*64
Dragon 200000 50000 53351 64*64*64
Table 5.1: Model complexities and number of ice particles. The columns of table from left to right are number of triangles for input benchmark, number of vertices for input benchmark, number of ice particles, and number of grids of marching cubes (MC).
We converted a closed polygonal model to a voxel model using model voxelization.
However, model voxelization fails if the model is unclosed. Nevertheless, we can use the method proposed in [NT03] to voxelize unclosed models.
The experimental results of ice melting simulation for these four benchmarks are shown in Table 5.2, including average execution time, number of generated water par-ticles, and number of triangles constructed by marching cubes. Each benchmark was ex-ecuted for 1000 frames.
# Water # MC Heat Phase Virtual Water Motion of MC
Particles Tri. Transfer Transition Volume Transfer Particles Construction Sum
Cube 3348.5 9190.6 12.24 1.01 0.57 23.16 29.54 66.52
Sphere 3795.6 6770.7 6.43 0.98 0.38 15.99 21.94 45.72
Bunny 6211.6 25777.0 33.25 1.40 1.15 54.62 124.54 214.96
Dragon 17390.1 35161 144.58 3.87 4.14 263.81 268.86 685.26
Table 5.2: The experimental results of ice melting simulation. The first two columns are average number of water particles and number of triangles constructed by marching cubes.
The right part is the average execution time for each step. (msec)
Water Particle Virtual Water Creation of Gradually Reduce Decomposing Volume Computing Water Particles Virtual Water Volume Sum
Cube 82.0 311.59 34.96 139.89 568.44
Sphere 97.72 179.65 22.12 76.86 376.35
Bunny 131.60 636.22 68.49 311.23 1147.54
Dragon 598.91 2531.07 251.34 762.68 4144.0
Table 5.3: Breakdown of virtual water volume transfer. (usec)
The step of the virtual water volume transfer includes four subparts, water particle
Compute Pressure Compute Force Compute Position Sum
Cube 0.36 0.03 22.77 23.16
Sphere 0.37 0.03 15.59 15.99
Bunny 2.12 0.04 52.46 54.66
Dragon 11.46 0.06 252.29 263.81
Table 5.4: Breakdown of motion of particles. (msec)
reducing virtual water volume, as shown in Table 5.3. These four parts were implemented on CPUs. The breakdowns of motion of particles are shown in Table 5.4. In our ice melting simulation, all the steps were implemented on GPUs except for the virtual water volume transfer and marching cubes construction.
Figure 5.5: Execution time of each step of ice melting simulation relates to the number of water particles. The benchmark used here is Cube.
In Table 5.2, we can see that there are two main time-consuming parts in our ice melt-ing method. Those are motion of particles computation and marchmelt-ing cubes construction.
The execution time used for computing the motion of water particles is proportional to the number of water particles. In the other words, the more water particles we have, the more time it consume, as shown in Figure 5.5. In Figure 5.5, the water particles gradually increase with time step, and we observe that only the execution time of motion of particles increases when the number of water particles increases.
Our method constructs the marching cubes on CPU. If we split the finer grids for marching cubes, we need more time on constructing marching cubes and getting finer polygonal model. We can utilize the GPU-based method proposed in [SEL11] and spend less time on marching cubes construction.
Performance of rendering procedure is shown in Table 5.5. Ray tracing method is used and the resolution is 720∗ 720. Our ray tracer uses reflection ray, refraction ray, and shadow ray. It tests for ray-triangle intersection and ray-metaball intersection. There are three parts, BVH construction for triangles, BVH construction for metaballs, and ray tracing. We compared the Table 5.5 and the left part of Table 5.2. Ray tracing consumes more time if the number of water particles and triangles are larger.
BVH Construction BVH Construction
Table 5.5: Execution time of rendering. (msec)
5.2 Experiment 1
In experiment 1, we compared the ice melting execution time for the same model with different level of voxelization. We split the scene into the following number of voxels:
16∗ 16 ∗ 16, 32 ∗ 32 ∗ 32, and 64 ∗ 64 ∗ 64. We voxelized the polygonal model into voxel model in each level. In this experiment, we used the benchmark Sphere, and the initial ice particles in each level of voxelization were 1089, 7399, and 54057. We calculated the execution time of ice melting simulation (excluding the marching cubes construction).
The experimental result is shown in Figure 5.6.
In the Figure 5.6, we can clearly see that the level of voxelization is the biggest factor affecting performance. The experiment with 64∗ 64 ∗ 64 voxels needs more execution time than others, even though the GPU parallel computing has been used.
Figure 5.6: Execution time of ice melting simulation with different level of voxelization.
The benchmark used here is Sphere.
5.3 Experiment 2
Our ice melting method is based on Iwasaki et al. [IUDN10]. They support two ways to render melting objects. One of the rendering methods uses marching cubes to construct both ice model and water particles. However, it loses realism when the particles are smaller than the grid of the marching cubes, or too close to the ice model. We separated the ice model and the water particles and use metaballs to represent the water particles.
We implemented the method proposed in [IUDN10]. Table 5.6 compared the timing of marching cubes construction and timing of ray tracing between the method proposed in [IUDN10] and ours. The benchmark we used for this comparison is Bunny.
Without Metaball With Metaballs
Table 5.6: Model complexity and timing comparisons (in msec) of marching cubes con-struction and rendering with and without metaballs. The model we use is Bunny.
The Table 5.6 shows that our method (with metaballs) spends less time on constructing marching cubes, but our ray tracing method spend a lot of time on ray-metaball intersection test. The performance of our method is worse than the method without using metaballs.
However, our method maintains integral water droplet, while their method cannot. We use Figure 5.7 and Figure 5.8 to explain.
Figure 5.7: This is a snapshot of the Bunny ear. In the same frame, the left hand side is the method without using metaballs, and the right hand side is our method. The red and yellow circles are different parts between two pictures.
Figure 5.8: This is also a snapshot of the Bunny ear. The red and yellow circles are different parts between two pictures.
We observe the result figure (Figure 5.7 and Figure 5.8) and find some differences between our method and the method without using metaball. Some water particles in their method disappears because the field value of the field scalar is smaller than the threshold
but in their method, the surface of water particle is angular (yellow circles). Their method can be improved by using finer grids of marching cubes, but it will severely affect the performance.
5.4 Experiment 3
In the experiment 3, we compare the result with the method proposed in [IUDN10].
Their method converts an ice particle into a water particle when the particle absorbed enough energy. The water particle replaces the original ice particle in the same place. Our method transfers the water volume from high to low via the ice particles, and generates water particles when the amount of water volume is too large.
In the real world, we observe that the ice on the floor melts and shrinks smoothly, as shown in Figure 5.9. We wouldn't see a water droplet appears directly on the surface of the ice cube. On the contrary, water will emerge on the floor around the ice cube. We compared this phenomenon between their method and ours, as shown in Figure 5.10.
Figure 5.9: The snapshots of real melting ice cube. www.youtube.com/ watch?
v=WgjksZoznuA
In the Figure 5.10, the top three pictures are a series of snapshots of melting ice cube using their method, and below is our results. We can see that their method directly creates water particles on the surface of ice cube. Our method transfers water volume via the ice particles and creates water particles on the floor.
5.5 Experiment 4
The method proposed by Kanamori et al. [KSN08] is good for rendering opaque object using depth peeling. This method can find the order of intersection between meta-balls and a ray in a short time. However, it suffers bad performance when the metameta-balls are transparent. In order to test the ray-metaball intersection efficiently, we use LBVH to build a hierarchy construction and speed up the execution time. The method of LBVH construction for metaballs is similar to the method that we introduced in Section 4.1. We
Figure 5.10: The top three pictures are the ice melting results proposed in [IUDN10], and the bottom three pictures are our results.
calculated the axis-aligned bounding box (AABB) for each metaball and we use the spher-ical center of metaball as its representative point. Then, we can construct Morton code for all metaballs. We compared the ray tracing performance between the method using and not using BVH construction for metaballs.
Table 5.7: Execution time of our ray tracing method with CPU-based BVH construction, with GPU-based BVH construction, and without BVH construction for metaballs. The model we use is Dragon. (msec)
In the Table 5.7, we observe that when the BVH construction for the metaballs is not used, we need more time on ray tracing procedure. We created BVH construction on CPU to improve the performance for ray tracing, but CPU-based BVH construction will still spend more time. Finally, we decided to use the LBVH construction to construct the
BVH for metaballs on GPUs. Compared with others, the timing of ray tracing for our method is 4.8x faster than the method without BVH for metaballs, and the timing of BVH construction is about 29.3x faster than the CPU-based BVH construction.
5.6 Experiment 5
We employed two types of geometry as heat sources in the ice melting simulation:
rectangle and sphere. The different types of heat source will radiate energy in different ways. The rectangle will uniformly emit the thermal photons along the normal direction of the rectangle. The photons in the sphere uniformly arrangement on the surface of the sphere and emit along the direction which is perpendicular to the surface.
We are interested in the rendering results for different heat sources. We show the results in Figure 5.11.
Figure 5.11: The top three pictures are the rendering results for the heat source Sphere, and the below three pictures are the rendering results for the heat source Rectengle.