= Physically Realistic Volume Rendering =

== Group Members ==
==== James Hegarty ====
==== Qiqi Wang ====

== 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.

|| attachment:Yosemite1.jpg || attachment:Yosemite2.jpg || attachment:Yosemite3.jpg ||

== 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.

attachment:Mesh1.png

attachment:Mesh2.png

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.

attachment:flow1.png

attachment:flow2.png

== Ray tracing through unstructured mesh (Qiqi Wang) ==

In PBRT, we added a file attachment: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:
{{{struct Face {
        int id;
        vector<Point> vertices;
        Cell *c0, *c1;
        Vector faceNormal;
        void CalculateNormal();
        bool Intersect( Point o, Vector d, float& t ) const;
        Point Center() const;
};

struct Cell {
        int id;
        vector<Face*> faces;
        vector<float> scalars;
        bool Intersect( Point o, Vector d, float& t0, float& t1,
                        Face*&f0, Face*&f1 ) const;
};
}}}

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:
{{{bool Face::Intersect( Point o, Vector d, float& t ) const {
    const float EPS = 0.0001;
    // calculate intersection point and t
    float total = Dot( d, faceNormal );
    float sect1 = Dot( vertices[0] - o, faceNormal );
    if (total == 0) return false;
    t = sect1 / total;
    // calculate total angle around itersection point
    float angle = 0;
    for (int i=0; (size_t)i<vertices.size(); i++) {
         int ip = (i+1) % vertices.size(); // next vertex
         Vector arm1 = vertices[i ] - o;
         Vector arm2 = vertices[ip] - o;
         // project to the foot of the ray
         Vector unitd = Normalize(d);
         arm1 -= unitd * Dot(arm1,unitd);
         arm2 -= unitd * Dot(arm2,unitd);
         // right on top
         if (arm1.Length() == 0 || arm2.Length() == 0) return true;
         // compute the angle of this side
         arm1 = Normalize(arm1);
         arm2 = Normalize(arm2);
         float sine = Dot( Cross(arm1, arm2), unitd );
         float cosine = Dot(arm1, arm2);
         float side_angle = atan2( sine, cosine );
         // add to the total angle
         angle += side_angle;
    }
    // intersect within the face if the total angle is +- 2*PI
    if (2*M_PI-EPS < fabs(angle) && fabs(angle) < 2*M_PI+EPS)
        return true;
    // otherwise the total angle must be 0
    else if ( fabs(angle) > EPS ) {
        Error("\nTotal angle equal to %f in intersection\n", angle);
        cout << "o = " << o << " d = " << d << "\n";
        for (int i=0; (size_t)i<vertices.size(); i++) cout << vertices[i] << endl;
        return false;
    }
    return false;
}

bool Cell::Intersect( Point o, Vector d, float& t0, float& t1, Face*&f0, Face*&f1 ) const {
    t0 = 1E10; t1 = -1E10;
    bool iIntersect = false;
    // find intersection with all faces.
    for (int i=0; (size_t)i<faces.size(); i++) {
        float t;
        if (faces[i]->Intersect( o, d, t )) {
            iIntersect = true;
            if (t <= t0) { t0 = t; f0 = faces[i]; }
            if (t >  t1) { t1 = t; f1 = faces[i]; }
        }
    }
    if (iIntersect && t0 == t1) {
        if (f0 != f1) {
            // decide which face is more on the side of the ray direction
            if (rand() > rand()) {
                Face* ff = f0;
                f0 = f1;
                f1 = ff;
            }
        }
        else {
            Error("\nt0 == t1 and f0 == f1 in Cell::Intersect, t0=%f, t1=%f\n", t0, t1);
            cout << "CellId = " << id << " o " << o << " d " << d << endl;
            cout << "f0 = " << f0->id << " f1 = " << f1->id << endl;
            return false;
        }
    }
    return iIntersect;
}}}

|| attachment:verify1.png || attachment:verify2.png || attachment:verify3.png ||

== Modeling of light scattering ==

||attachment:rainbow1.jpg||attachment:rainbow2.png||attachment:rainbow3.jpg||

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