Physically Realistic Volume Rendering
Group Members
James Hegarty
Qiqi Wang
Final Image
Added / Modified Code
Overview
In this project we are planning on rendering a scene based on a physically accurate light scattering model for volumetric materials. Such materials can include fire, smoke, vapor, and clouds. Special effects such as rainbows can be also rendered with this model. In order to accurately model the distribution of volume materials in complex geometry, a physically based fluid calculation must be performed, resulting in an unstructured mesh. A ray tracing algorithm for unstructured meshes needs to be developed, and light reflection and scattering models from these volume materials needs to be built.
The following pictures are from www.ImageMontage.com , I believe they are taken at the Nevada fall in Yosemite. The waterfall and rainbow is the type of things we want to render.
|
|
|
Modeling of volume material distribution (Qiqi Wang)
In order to model the distribution of volume materials accurately, we will use the commercial flow solver Fluent. In the waterfall scene, we use commercial Commputational Fluid Dynamics software Fluent 6.2 to calculate the flow of air around the waterfall, as well as the density of water drops in the air. The size of the mesh is 1000 ft by 1000 ft by 6000 ft, the waterfall alone is 100 ft tall. The following pictures demonstrates the mesh generated for the Computataional Fluid Dynamics calculation. The first picture shows the surface mesh, the second picture shows the volume mesh on the symmetric section plane.
This is a typical unstructured mesh used for Computational Fluid Dynamics calculation around complex geometry. First of all, tt has huge mesh size variation. The cells inside and near the waterfall is small, and the cells away from it is ten thousand times larger in volume. Also, it has three different types of cells. Inside the waterfall, octohedral cells are used; away from the waterfall, hexahedral cells are used; pyramid cells are used as a layer between octohedral cells and hexahedral cells.
The velocity field and the density of water drops are computed by Fluent 6.2, a most up-to-date commercial Computational Fluid Dynamics software. The small volume in the middle of the step is used to model the water fall, and the rest of the volume is used to model the air. The flow field is calculated using incompressible Navier-Stokes equation with k-e turbulence model. We model the driving force by the flow of water inside the waterfall as a constant downward source term in the momentum equation. This source term drives the flow field in the entire volume. We model the floating water drops in air as transporting and reacting chemical species, whose density is also solved by Fluent 6.2. Inside the waterfall volume, the water drop density is generated in a constant rate. In the rest of the volume, the water drops evaporates into vapor and convects away from the waterfall. The following two images shows the flow field on the symmetric cut plane and the isosurface of 20% of the maximum water drops concentration. Our rainbow image is rendered based on this calculated density field of the water drops.
Ray tracing through unstructured mesh (Qiqi Wang)
In PBRT, we added a file volumemesh.cpp in the volumes folder to deal with unstructured mesh ray tracing. The bulk of this code is four classes (VolumeMesh, Mesh, Cell, Face) and mesh file I/O routines. In our data structure of storing the mesh, each face is an object of the Face class, which stores position of all its vertices, its face normal and two pointers to its adjacent cells; each cell is an object of the Cell class, which stores pointers to all its faces and the water drops density in the cell. The definition of the classes are as follows:
1 struct Face {
2 int id;
3 vector<Point> vertices;
4 Cell *c0, *c1;
5 Vector faceNormal;
6 void CalculateNormal();
7 bool Intersect( Point o, Vector d, float& t ) const;
8 Point Center() const;
9 };
10
11 struct Cell {
12 int id;
13 vector<Face*> faces;
14 vector<float> scalars;
15 bool Intersect( Point o, Vector d, float& t0, float& t1,
16 Face*&f0, Face*&f1 ) const;
17 };
Both Face and Cell class have an Intersect method, which determines whether a given ray intersect with the Face / Cell, and calculate the intersection point(s). To calculate the intersection point of a ray with a face, we just need to compute the intersection of the ray with the face plane, which can be done in 2 dot products and 1 division. Determination of whether the ray intersects with the face is much nore difficult. Our algorithm is to calculate the total angle of the face edges with respect to the ray. If the total angle is plus or minus 2*pi, then the ray intersects with the face; if the total angle is 0, then the ray does not intersect with the face. The total angle cannot be any other value. For robustness reason that will be explained later, we project the face onto a plane that is normal to the ray before calculating the total angle. It takes roughly 4 dot products, 1 cross product, 4 normalization and 1 arctan2 operations for each edge of the face to determine intersection. To calculate the intersection points of a ray with a cell, we test intersection of the ray with all its faces. Only 0 or 2 faces should have an intersection with the ray assuming the cell is convex. The code for the two Intersect methods are:
1 bool Face::Intersect( Point o, Vector d, float& t ) const {
2 const float EPS = 0.0001;
3 // calculate intersection point and t
4 float total = Dot( d, faceNormal );
5 float sect1 = Dot( vertices[0] - o, faceNormal );
6 if (total == 0) return false;
7 t = sect1 / total;
8 // calculate total angle around itersection point
9 float angle = 0;
10 for (int i=0; (size_t)i<vertices.size(); i++) {
11 int ip = (i+1) % vertices.size(); // next vertex
12 Vector arm1 = vertices[i ] - o;
13 Vector arm2 = vertices[ip] - o;
14 // project to the foot of the ray
15 Vector unitd = Normalize(d);
16 arm1 -= unitd * Dot(arm1,unitd);
17 arm2 -= unitd * Dot(arm2,unitd);
18 // right on top
19 if (arm1.Length() == 0 || arm2.Length() == 0) return true;
20 // compute the angle of this side
21 arm1 = Normalize(arm1);
22 arm2 = Normalize(arm2);
23 float sine = Dot( Cross(arm1, arm2), unitd );
24 float cosine = Dot(arm1, arm2);
25 float side_angle = atan2( sine, cosine );
26 // add to the total angle
27 angle += side_angle;
28 }
29 // intersect within the face if the total angle is +- 2*PI
30 if (2*M_PI-EPS < fabs(angle) && fabs(angle) < 2*M_PI+EPS)
31 return true;
32 // otherwise the total angle must be 0
33 else if ( fabs(angle) > EPS ) {
34 Error("\nTotal angle equal to %f in intersection\n", angle);
35 cout << "o = " << o << " d = " << d << "\n";
36 for (int i=0; (size_t)i<vertices.size(); i++) cout << vertices[i] << endl;
37 return false;
38 }
39 return false;
40 }
41
42 bool Cell::Intersect( Point o, Vector d, float& t0, float& t1, Face*&f0, Face*&f1 ) const {
43 t0 = 1E10; t1 = -1E10;
44 bool iIntersect = false;
45 // find intersection with all faces.
46 for (int i=0; (size_t)i<faces.size(); i++) {
47 float t;
48 if (faces[i]->Intersect( o, d, t )) {
49 iIntersect = true;
50 if (t <= t0) { t0 = t; f0 = faces[i]; }
51 if (t > t1) { t1 = t; f1 = faces[i]; }
52 }
53 }
54 if (iIntersect && t0 == t1) {
55 if (f0 != f1) {
56 // decide which face is more on the side of the ray direction
57 if (rand() > rand()) {
58 Face* ff = f0;
59 f0 = f1;
60 f1 = ff;
61 }
62 }
63 else {
64 Error("\nt0 == t1 and f0 == f1 in Cell::Intersect, t0=%f, t1=%f\n", t0, t1);
65 cout << "CellId = " << id << " o " << o << " d " << d << endl;
66 cout << "f0 = " << f0->id << " f1 = " << f1->id << endl;
67 return false;
68 }
69 }
70 return iIntersect;
71 }
The Mesh class manages all the cells and faces. It also manages a 'Current Point' and 'Current Cell' that can be efficiently moved by tracing a ray through the mesh volume. The ray tracing algorithm through the volume mesh is implemented with the MoveTo method of this class.
1 class Mesh {
2 public:
3 Mesh( string filename, int numScalars );
4 Point CurrentPoint() { return currentPt; }
5 const Cell* CurrentCell() { return &cells[currentCell]; }
6 bool MoveTo( Point target, vector<const Cell*>& trajectory, vector<float> t0, vector<float> t1 );
7 BBox bbox() const { return BBox(p0, p1); }
8 private:
9 Point currentPt;
10 int currentCell;
11 int numFaces, numCells;
12 Face* faces;
13 Cell* cells;
14 Point p0, p1;
15 void CalculateBBox();
16 };
17
18 bool Mesh::MoveTo( Point target, vector<const Cell*>& trajectory, vector<float> vt0, vector<float> vt1 )
19 {
20 const Cell *c = CurrentCell();
21 Point o = currentPt;
22 Vector d = target - currentPt;
23 float t0 = 0, t1 = 0;
24 int iter = 0, maxiter = 10000;
25 while (c != NULL) {
26 if (iter++ > maxiter) break;
27 Face *f0 = NULL, *f1 = NULL;
28 bool isect = c->Intersect(o,d,t0,t1,f0,f1);
29 if (!isect) {
30 Error("Failed to find intersection for cell %d\n", c-cells);
31 cout << "o = " << o << ", d = " << d << "\n";
32 return false;
33 }
34 assert( f0!=NULL && f1!=NULL );
35 // save trajectory
36 trajectory.push_back(c);
37 vt0.push_back(t0);
38 vt1.push_back(t1);
39 if (t1 >= 1) break;
40 // move to next cell
41 if (f1->c0 == c) c = f1->c1;
42 else if (f1->c1 == c) c = f1->c0;
43 else
44 Error("Mesh structure error detected in Mesh::MoveTo, Face %d with Cell %d\n", f1-faces, c-cells);
45 }
46 if (iter > maxiter) {
47 Error("Mesh::MoveTo, maximum iterations %d reached.\n", maxiter);
48 return false;
49 }
50 // modify current point
51 if (c != NULL) {
52 currentCell = c - cells;
53 currentPt = target;
54 return true;
55 }
56 // return value indicates whether we are able to move to the point
57 else return false;
58 }
The VolumeMesh class is a subclass of DensityRegion. It manages a Mesh class. It overrides the Density function and returns the density of any given point by use Mesh::MoveTo to move the 'Current Point' to the target point and return the density of the 'Current Cell'.
1 float VolumeMesh::Density(const Point &Pobj) const {
2 if (!extent.Inside(Pobj)) return 0;
3 vector<const Cell*> trajectory;
4 vector<float> t0, t1;
5 if (mesh->MoveTo( Pobj, trajectory, t0, t1 ))
6 return trajectory.back()->scalars[0];
7 // return 1;
8 else return 0;
9 }
Finaly, the mesh file I/O routines are used in the constructor of Mesh class. It counts the number of cells and faces in the mesh, allocates memory space and reads the file content into our data structure.
Robustness Issue
In rendering each image, millions of rays sent into the volume mesh, and each ray intersect with thousands of faces and cells. Because of the huge number of intersections, extremely rare situations can occur. Rays that goes right through an edge creates two numerical problems: (1) Rays that goes right through an edge of a cell might be found to intersect with 1 or 3 faces of that cell. (2) The two intersect points of a cell can be two close to determine which is farther. Using double precision numbers instead of single precision numbers reduces the number of such numerical problems, but cannot solve them completely. In order to deal with that, we need to twick the algorithm. First of all, when calculate the total angle around the ray, we project the face into a plane that is orthogonal to the ray. This makes the calculated angle of each edge with respect to a given ray independent of the face. Also, when the two intersection points are close, we use a randomized algorithm to determine which face is farther away in order to prevent dead lock.
Verification of ray tracing algorithm
The following three images are rendered to verify that we have a correct unstructured mesh ray tracing algorithm. In the first image, we set the light emmision factor to be constant within the entire volume. The rendered result matches the shape of the mesh in the previous section. In the second and third image, we set the light emmision factor to be proportional to the water drops density calculated by Fluent 6.2. We see the shape matches the isosurface of waterdrops density shown in the previous section.
|
|
|
Modeling of light scattering
|
|
|
We will use a physically-based model to simulate the scattering of light in volumes. Our model will cover scattering with opaque particles such as smoke, in addition to transparent particles such as water. In the case of transparent particles, light of different wavelengths has different refraction indices inside the particle, and internal reflection occurs, leading light that enters the water droplet from behind the observer to reflect back at the observer spread out horizontally. This effect is commonly known as a rainbow. Our model will not involve ray tracing of individual particles, because this would require a prohibitive amount of calculation time. Instead, we will make simplifying assumptions based on the light and camera angle and properties of the particle being rendered.
While PBRT does already include some models for volume rendering, they do not account for transparent particles and angle dependent effects, like we want to render. We also suspect that there may be the possibility of improving PBRT's efficiency in rendering our volumetric meshes.
Sources
Rainbows, Halos, and Glories by Robert Greenler