Infrared Image Synthesis. C. Garnier et al. [GCFM*99] described an IR sensor model
developed from study of physical effects involved in IR image acquisition process. Their approach consists of a combination and an extension of current camera models used in visible and infrared image synthesis, and they also merges ray tracing and post-processing techniques.
Zhangye Wang et al. [WWBP01] established an infrared model for ground targets, such as tank. They considered the effect of inner heat source, friction and various environmental factors and drew infrared images of target at various states by the Computer Graphics techniques. Zhangye Wang et al. [WPLJ01] and Zhaoyi Jiang et al. [JWJP04] proposed new IR image synthesis model accounting for meteorological, environmental, material and artificial factors. The energy equilibrium equation is built based on the principle of heat transfer and infrared physics and a finite difference method is adopted to solve the equations.
Zhaoyi Jiang et al. [JWP03] proposed a method to constitute the dynamic infrared scene by combining scene geometric modeling with multi-spectral image. Multi-spectral IR textures are correspondent to materials of object surface. The actual attenuated IR radiation intensity is expressed by Phong illumination model between 10℃and 50℃.
Figure 2.1 The synthesis procedure of infrared targets and the background [WWBP01].
10
Molecular Dynamic for Conduction. Molecular dynamic (MD) methods are now generally
accepted approaches for simulating molecular scale models of matter. The essence of MD simulation methods is simply stated: numerically solve the N-body problem of classical mechanics. H. Xue and C. Shu [XS99] presented investigation to deals with the equilibration of heat conduction simulation in a very thin film using MD. David M. Harrild et al. [HPH00]
described the novel application of a Finite Volume method (FVM) derived from computational fluid dynamics to the field of computational cardiac electrophysiology. They applied this method to simulate conduction in an arbitrarily shaped or complex region.
Jiaocheng Ma et al. [MXJ08] proposed real-time mathematical 2D heat transfer and solidification model. This model was presented for billet continuous casting of low carbon steel and solved by FVM. The GPU-based implementation is faster than that of CPU-based MD simulation. Juekuan Yang et al. [YWC07] presented an implementation of MD simulation on modern graphics processing units (GPU). This GPU algorithm was used to calculate the thermal conductivities of solid argon and to reduce the total computational time of conduction simulation at high performance. These methods mentioned above are used to analyze and solve heat flux propagation on a single solid object but they are difficult to represent interaction with other objects.
Figure 2.2 The mesh used to incorporate a complex conduction environment and the result of stimulation [HPH00].
11
Fluid dynamic. The solution to a fluid dynamics problem typically involves calculating
various properties of the fluid, such as velocity, pressure, density, and temperature, as functions of space and time. Wei Li et al. [LWK03] presented a physically-based flow simulation which supports complex boundary conditions running on the general-purpose graphics hardware. For computing the flow field, they accelerated the computation of the Lattice Boltzmann Method (LBM) on GPU, by grouping particle packets into 2D textures and mapping the Boltzmann equations completely to the rasterization and frame buffer operations.
To handle complex, moving and deformable boundaries, they proposed a generic voxelization algorithm of the boundaries using depth peeling, and extended it to a dynamic boundary generation method that converts any geometric boundary to LBM boundary nodes on-the-fly.
Zhe Fan et al. [FKZQ*09] presented a simulation and visualization system for the thermal fluid dynamics inside a pressurized water reactor of a nuclear power plant when cold water is injected into the reactor vessel. They employed a hybrid thermal lattice Boltzmann method, which has the advantages of ease of parallelization and ease of handling complex simulation boundaries.
Figure 2.3 Particles advected in a 2D flow field based on the D2Q9 LBM model [LWK03].
12
3. Overview
In this chapter, we will briefly describe our method and each chapter afterward.
We propose our voxelization algorithm in chapter 4. The energy of molecules on object surface may not be the same in dynamic state of heat transfer. The objects in simulation should be divided into voxels as tiny as well. We generate voxelized textures from objects by depth peeling in axis-aligned orthogonal views. These textures which stored energy and geometrical information are able to compute variation of heat in parallel through GPU. Since the method is only valid to surfaces, we fill the voxels in multiple layered textures with large weight of heat capacity for simulating volume of object. By comparing with full volume resolution, our methods compress voxel number into 10~25%.
Chapter 5 includes the details of heat transfer simulating. In order to rapidly generate visually reasonable scene in simulation of radiation, we ignore the interactive reflection between objects with low temperature difference. We implement shadow map and Phong shading model while low intensity heat source. In addition, screen space ambient occlusion (SSAO) [Mitt07, RGS09] is used to enhance detail in high temperature difference.
To simulate conduction, we separate the render pass into pruning, contact detection and spreading phases. First, we use the linear-time CULLIDE [GRLM03] algorithm to prune collision-free objects in image space which are unnecessary for computation of heat exchange.
Once a collision is detected, the contact surfaces make heat pass through from hot to cool side.
We apply Layered Depth Images (LDI) extendedmethod [FBAF08] which processes contacts between objects bounded by triangular surfaces. The last step is spreading. Each object
13
spreads heat through itself until the distribution of heat is equilibrium.
Computation with computational fluid dynamics (CFD) usually involves intensive computation in general convection processing. Due to our real-time and visually acceptable requirement, we represent the shape of object by combination of cylinders. Heat transfer between fluid and such simple shape can get approximate distribution [BJO98, KCY06]. By appling this approximation we can speed up the calculation between surface and fluid from O(n3) to O(n2logn) which n defined boundary of the scene.
Our experiment and result are presented in chapter 6. According to performance and visual result, we provide optimal parameters for this simulation system. At last, we summarize our approach and discussion in chapter 7.
14
4. Voxelization of 3D Model
Since the temperatures between various parts of an object may not be the same, we should divide the object into regular volumes, so called voxels, with its own thermal energy for simulating heat transfer in dynamic state. A voxel is a volume element, comprising a value on the regular grid in three dimensional space. In our system, the variations of temperatures on object surfaces are displayed by voxels in continuous domain.
4.1 GPU-Based Voxelization
Figure 4.1 From left to right: mesh of penguin, large voxel size and small voxel size voxelized model by slice method. The smaller voxel size results more precise shape of mesh and more storage data used.
An intuitive voxelization approach is the slicing method [FC00]. This method sets distance of near and far clips planes as a single voxel size and renders only the geometry falling into the slab between the two clip planes. The clip planes are shifted to generate subsequent slices until entire volume covered with these slices. Our simulation method mostly interacts around surfaces of objects which are sparse in a volume slice in most case. In other words, only a small percentage of voxels are intersected by the boundary surfaces. There is no need to voxelize the empty space that corresponds to non-boundary voxels. Furthermore, only self
15
heat spreading involves interior volume and we present how we reduce the amount of voxel and keep visually-similar approximation in chapter 4.3.
Figure 4.2 The flow chart of GPU-based voxelization algorithm.
The GPU-based voxelization algorithm we used is the method proposed by Wei Li et al.
[LFWK03]. They avoid a slicing method by the idea of depth peeling [Ever01] used for order-independent transparency. The depth layers in the scene are stripped away with successive rendering passes. The procedure can be divided into three steps:
1. At first, the scene is rendered normally and the layer of nearest fragments as voxels is obtained. During the second rendering pass, each fragment compares depth value with the depth texture obtained from the depth buffer of the previous pass.
True
Mesh Data
Depth peeling
Occlusion Query
Voxelized Texture
End state
False
16
2. The fragment is discarded if either depth test is not pass or depth value greater than that depth texture.
3. The process continues until no fragment is farther away than the corresponding pixel in the depth texture. This condition is best determined by using a hardware occlusion query, which returns the number of pixels written to the frame buffer.
We apply the peeling process three times. Each time, the image plane is orthogonal to one of the major axes. The viewport is translated so that the layer images do not missing voxels but are tiled as tightly as possible.
Figure 4. 3 Voxelized textures of “Happy Buddha”. There are 18 layers in x axis, 20 in y and 8 in z; the left shows depth value in the first three textures of z axis.
We apply depth peeling and save layer images with depth value in textures named voxelized texture. There are three series of voxel data from each orthogonal view in these textures as shown in Figure 4.3. The 3D position of voxel can be restored as texture coordinate and depth in voxelized texture. As a result of peeling from three orthogonal views, some of the voxels may be rendered more than once. The replication does not affect the accuracy but we can reduce it to save more space and computing times.
Y Axis
Z Axis
X Axis
17
4.2 Optimization of Voxelized Texture
We get voxelized textures by three-view depth peeling algorithm but these voxelized texture data can further be optimized for acceleration of simulation system.
The improvement mainly occurs at overlapping areas in voxelized textures if the areas can be rendered from two or three orthogonal views. We take a sphere for example. Figure 4.4 shows the first voxelized textures from x and z axis. These textures which are rendered as hemisphere of the sphere overlap at the same quarter with another one. Though all voxels are unique in single orthogonal view, the same area may be rendered repeatedly by each peeling process in the worst case. That means we may simulate on the same voxel twice or more but it is not necessary. Therefore, we should store each voxel only once in all voxelized textures.
Figure 4.4 The first layer images from orthogonal view on X and Z axis overlap on cross area.
To keep the structure of voxel integrity, stored areas in each axis-aligned texture must be independent of other orthogonal views. Then, we consider the attribute of surface normal to check each fragment on voxelized texture whether the absolute normal weight of view axis is greater than other two axes. For instance, if n is normal vector of surface and absolute normal
X Axis Z Axis
Overlapping Area
18
weight of x-axis ∥nx∥ is greater than ∥ny∥ and ∥nz∥, the corresponding fragment is preserved as a voxel in x-axis peeling process.
Figure 4.5 Voxelized texture of sphere rendered from x axis: (a) depth value; (b) normal of surface; (c) repainting red on fragment if ∥nx∥ is greatest, green for y axis and blue for z;
(d) preserved area of (c).
After removing replication of voxelized textures, more empty space may be generated due to image space quantization at view change. As shown in Figure 4.6 (a) and (b), some of these textures contain only few voxels or completely empty. We pack texels to fill empty space before its texture as forward as possible and remove textures without any voxel data. Through packing textures process, the number of layers is decreased which affects the number of simulation computing times. We compare the textures of “Happy Buddha” with/without data packing. The number of layers decreases to a half.
The last step of optimization is to build connection of voxels. We do not know whether these voxels connect others around after packing. Even though we can find the connection by restored 3D position of voxels, it has to be precomputed to reduce traversal neighbor texels every time.
19
(a)
(b)
(c)
Figure 4 6 Voxelized texture of “Happy Buddha” rendered from z axis with surface normal, from left to right: (a) original, (b) stored and (c) packed texture. After packing, the number of layers in x is 8, 11 in y and 6 in z axis.
We summarize the procedure of optimizing voxelized textures. First we remove the replicate voxels by comparing absolute normal weight of orthogonal axis. The voxels with maximum absolute normal weight in viewing axis should be preserved. We pack the preserved voxels as forward as possible in voxelized textures for fill the empty spaces. After packing, we reduce the layer if layer of texture is completely empty. At last step we reconstruct the connection of voxels and save connecting information in corresponding textures. Figure 4.7 represents the sample of optimizing procedure.
20
Figure 4.7 The example of optimization of voxelized texture.
4.3 Multi-Level Voxelized Texture
Our voxelization algorithm generates only a single shell of target model. However, not all of objects fit in this case; it may be solid or more complicated structure inside. We will deal with solid and uniform material for conduction simulating. For this reason, we build voxelized textures inside the shell.
Figure 4.8 Voxelized texture of sphere slice. From left to right is single shell, multi-level and irregular multi-level. The irregular multi-level varies voxel size with different levels.
Because of overlap occurring on some voxels, the conduction weight between voxels should be adjusted.
We can build inner textures easily from original object and shell textures. If we construct all Level 1
21
voxels inside of object with the same size of shell voxels, the number of textures explosively increases and it cannot take our advantage of acceleration anymore. The voxels inside should be considered with larger capacity than those outside, and then we only need to build fewer textures for interior voxels.
We use a half resolution to record interior voxel data. Such as mipmap texture, each level means one more layer inside and texture size smaller than level before. Figure 4.8 shows the voxelized texture of sphere slice by single shell, multi-level and irregular multi-level. The voxel size is presented as conduction weight between different levels of textures in simulation.
Besides, we must be careful about increasing voxel size, and total capacity must be kept in original volume. We discuss multi-level and irregular multi-level structure in chapter 5 and 6 in detail.
We construct multi-level structure on full volume voxel data but our voxelized texture. Since adding one level means to add the same size of our voxelized texture, it is difficult to divide a large amount of voxels inside into three axis-aligned textures. However, we can build the irregular multi-level structure in voxelized texture. There are few irregular voxels inside because of the voxel size increasing exponentially each level.
22
5. Simulation of Heat Transfer
5.1 Radiation from Heat Source
According to Stefan–Boltzmann law, the power emitted per unit area of the surface of a black body is directly proportional to the fourth power of its absolute temperature. That is
𝑞
"= 𝜎𝜀𝑇
4(1)
Where q” is the total power radiated per unit area, T is the temperature in the Kelvin scale, ε is the thermal emissivity and σ = 5.67×10−8W·m−2·K−4 is the Stefan–Boltzmann constant. It means that anything emits thermal radiation if it is higher than absolute zero in temperature.
Besides, thermal radiation may be absorbed, reflected or transmitted. It is difficult to apply radiosity algorithm with each feature in real-time. For the reason, we approximate parts of processes of total radiation and simulation with the following assumptions.
At first, the heat source, such as the sun, emits thermal energy to the scene. In general case, the environment temperature is about 300K and the sun is about 5800K. During the daytime, the sun emit energy about 1000W·m-2 to Earth. In contract, the energy emits from other object can usually be omitted because of low temperature difference. We assume the max temperature difference in the scene to be 100K, if there is no additional heat source, and the max radiation is about 31.15W·m-2(distance between objects: 10m). We project this energy to cubes with different materials whose volume is 1m3 and assume cubes absorbing energy fully as black body. The variations of temperature are shown as Table 5.1. Actually, these materials cannot completely absorb the energy and the variations of temperature may be close to zero.
23
As this result, we only calculate radiation from heat source to the scene but from other objects.
Table 5.1 The variations of temperature by fully absorbing 31.15Wm-2 on 1m3 cube.
Material Density kg/m3 Heat Capacity(25℃)
J/(kg·K)
Gasoline 0.73×103 2220 1.92×10-5
Since the reflection term is based on intensity of heat source but variation of temperature, it cannot simplify as emission. At assumption of the radiation under a single wavelength, objects between Ai and Aj, the intensity of radiation from patch dAi to dAj is:
dq
i→j= I
icosθ
idA
idA
jcosθ
jR
ij2(2)
24
Ii is the intensity of Ai, θi andθj is included angles between 𝑑𝐴 and normal of patches 𝑖𝑑𝐴𝑗 and Rij is length between dAi and dAj. If the distance is too far from heat source to object surface or heat source is very small as a point, such as the sun, the intensity of radiation can be formed as the diffuse term in Phong reflection model.
I
d= k
d(L∙N)i
d(3)
kd is diffuse reflection constant, L is the direction vector from the point on the surface toward each light source, N is the normal at this point on the surface and id is intensity of light source.
The ambient term in Phong model is to simulate accumulated inter-reflection in the scene but it is unidirectional and too rough for complicated scene. The reflection of radiation is similar with global illumination which contains direct and indirect light computing. We can get the first reflection from heat source easily but it’s costly to compute indirect reflection repeatedly until thermal equilibrium. For simulating dynamic scene in real-time, we use Phong model and shadow map to render direct light and ambient occlusion map to approximate indirect light if the scene is complicated.
We assume that all materials in the scene are opacity to simplify the problem. For these assumptions, we only deal with emission and reflection of radiation from heat source as the scene by using Phong model and shadow map, optional indirect reflection by ambient occlusion if the scene is complicated. We render the emission of heat source to voxelized textures represented each object in the scene. The reflection is rendered in screen independently and blended with final results.
25
5.2 Heat Propagation
The area irradiated by heat source gains energy and the heat is propagated in the whole solid object until energy balance. In each rendering pass, we spread the energy of voxels to connected neighbor for each voxelized texture.
𝑞 = 𝑘 ∙ 𝐴
𝐿 ∆𝑇 (4)
In conduction equation (4), k is the thermal conductivity of material, ΔT is variation of temperature, A is contact area and L is the distance between two elements. For a single object, we can get A, L and ΔT from voxelized texture and set the same k for all voxels.
In conduction equation (4), k is the thermal conductivity of material, ΔT is variation of temperature, A is contact area and L is the distance between two elements. For a single object, we can get A, L and ΔT from voxelized texture and set the same k for all voxels.