Project code https://github.com/Remyuu/GAMES101-Homework
This article is divided into two parts:Path tracing code implementationandMicro-texture model.
We are HW.5 Constructed the Whitted-Style Ray Tracing algorithm ray tracing project, HW.6 The BVH acceleration structure is used to speed up the intersection process. This time, we built Path Tracing's ray tracing and used multi-threading to accelerate rendering. Finally, the micro-surface model is used to provide the project with a more rough material.
In addition, it should be noted that the content of this article on microsurface models is mainly derived from Ref.5 , mainly talks about the basic theory and code implementation of the Cook-Torrance model.
This article basically explains the entire content of the framework. If there are any errors in the content, please point them out. This project is about rendering a CornellBox scene, and the final effect is roughly as shown below:
Parameter 1: {SSP:64, res:{784, 784}, parallel: false, RussianRoulette = 0.8}, render time:{4101 seconds},
Parameter 2: {SSP:64, res:{784, 784}, parallelism: true, RussianRoulette = 0.8, cookTorrance, PDF = GGX}, render time: {3415 seconds}
Seven Frameworks of OperationDownload (Please forgive the slow download of self-built small water pipe)
Project Flow – main.cpp
As usual, we start the analysis from the main function.
The process of this project is very simple:Set up the scene and render it.Let’s take a closer look.
First initialize the Scene object and set the scene resolution.
// Change the definition here to change resolution Scene scene(784, 784);Create four materials - red, green, white and light. These materials use the DIFFUSE type and set different diffuse coefficients (Kd).
Material* red = new Material(DIFFUSE, Vector3f(0.0f)); red->Kd = Vector3f(0.63f, 0.065f, 0.05f); ... Material* light = new Material(DIFFUSE, (8.0f * Vector3f (0.747f+0.058f, 0.747f+0.258f, 0.747f) + 15.6f * Vector3f(0.740f+0.287f,0.740f+0.160f,0.740f) + 18.4f *Vector3f(0.737f+0.642f,0.737f+0.159f,0.737f))); light->Kd = Vector3f (0.65f);Create the Cornell scene objects and then add them to the scene.
MeshTriangle floor("../models/cornellbox/floor.obj", white); ... MeshTriangle light_("../models/cornellbox/light.obj", light); scene.Add(&floor); .. . scene.Add(&light_);Build a BVH to accelerate collision detection between rays and objects in the scene.
scene.buildBVH();Finally, a renderer object r is created to render the scene and the rendering time is recorded.
Renderer r; auto start = std::chrono::system_clock::now(); r.Render(scene); auto stop = std::chrono::system_clock::now();The above is the general process of the project.
Object abstract base class – Object
The Object class defines all the basic behaviors that an object needs in the ray tracing algorithm. It uses pure virtual functions to indicate that this is an interface that needs to be inherited by specific object classes (MeshTriangle, Sphere, and Triangle) and implement these methods.
For detailed instructions, see the following code comments:
Object() {} virtual ~Object() {} virtual bool intersect(const Ray& ray) = 0;// Used to determine whether a ray intersects with the object virtual bool intersect(const Ray& ray, float &, uint32_t &) const = 0;// Also used to detect whether a ray intersects with an object, but this function also returns a parameterized representation of the intersection and the intersection index virtual Intersection getIntersection(Ray _ray) = 0;// Returns the intersection information between the ray and the object virtual void getSurfaceProperties(const Vector3f &, const Vector3f &, const uint32_t &, const Vector2f &, Vector3f &, Vector2f &) const = 0;// This function is used to obtain the properties of the object surface, such as the surface normal, texture coordinates, etc. virtual Vector3f evalDiffuseColor(const Vector2f &) const =0;// Evaluates the diffuse color of an object at a specific texture coordinate virtual Bounds3 getBounds()=0;// Returns the bounding box of the object virtual float getArea()=0;// Returns the surface area of the object. The calculation method for each shape can be different virtual void Sample(Intersection &pos, float &pdf)=0;// Sample a point from the surface of the object for light source sampling. The `pos` parameter is the information of the sampling point, and `pdf` is the probability density function value of the point. virtual bool hasEmit()=0;// Determine whether the object emits light, that is, whether it is a light source.Based on this class, we also created three specific object classes: MeshTriangle, Sphere and Triangle. These three classes are subclasses of the object class Object and are used to represent different geometric shapes in three-dimensional space.
Since we need to render the picture by ray tracing, we need to implement an important operation called intersect, which is used to detect whether a ray intersects with an object.
In addition, each class has a data member m of the Material type, which represents the material of the object. The material defines the object's color, texture, emitted light and other properties.
In Triangle.hpp, rayTriangleIntersect uses the Möller-Trumbore algorithm to determine whether a ray intersects a triangle in 3D space, and if so, it can calculate the exact location of the intersection. For a detailed code explanation, please see my other articlearticle , the main steps are written in the code comments.
bool rayTriangleIntersect(const Vector3f& v0, const Vector3f& v1, const Vector3f& v2, const Vector3f& orig, const Vector3f& dir, float& tnear, float& u, float& v){ // First, calculate the vectors of the two sides of the triangle (edge1 and edge2), and then calculate a new vector pvec based on the vector product (outer product) of the ray direction dir and the edge edge2. Vector3f edge1 = v1 - v0; Vector3f edge2 = v2 - v0; Vector3f pvec = crossProduct(dir, edge2); float det = dotProduct(edge1, pvec); if (det == 0 || det < 0) return false; // Then, by calculating the dot product (inner product) of pvec and edge edge1, a determinant value is obtained. If this value is 0 or negative, it means that the ray is parallel to the triangle or the ray is in the opposite direction of the triangle, and false should be returned at this time. Vector3f tvec = orig - v0; u = dotProduct(tvec, pvec); if (u < 0 || u > det) return false; // Then, calculate the cross product of tvec and edge1 to get qvec, and calculate its dot product with dir to get v. If v is less than 0 or u+v is greater than det, return false. Vector3f qvec = crossProduct(tvec, edge1); v = dotProduct(dir, qvec); if (v < 0 || u + v > det) return false; float invDet = 1 / det; // Finally, if all tests pass, the ray intersects the triangle. Calculate the depth tnear of the intersection point, as well as the barycentric coordinates (u, v) inside the triangle. tnear = dotProduct(edge2, qvec) * invDet; u *= invDet; v *= invDet; return true; }Due to space constraints, we will only focus on these three categories.
Triangle
Next, a Triangle object represents a triangle in three-dimensional space.
Constructor
Triangle(Vector3f _v0, Vector3f _v1, Vector3f _v2, Material* _m = nullptr) : v0(_v0), v1(_v1), v2(_v2), m(_m) { e1 = v1 - v0; e2 = v2 - v0; normal = normalize(crossProduct(e1, e2)); area = crossProduct(e1, e2).norm()*0.5f; }Each triangle has three vertices (v0, v1, and v2), two edge vectors (e1 and e2), a normal vector (normal), an area (area), and a material pointer (m). In the triangle constructor, the edge vector, normal vector, and area are calculated based on the three vertices input. The area is calculated by half the modulus of the cross product of e1 and e2.
Triangle related operations
There are three functions (Sample, getArea and hasEmit) that are directly overridden.
... void Sample(Intersection &pos, float &pdf){ float x = std::sqrt(get_random_float()), y = get_random_float(); pos.coords = v0 * (1.0f - x) + v1 * (x * (1.0f - y)) + v2 * (x * y); pos.normal = this->normal; pdf = 1.0f / area; } float getArea(){ return area; } bool hasEmit(){ return m->hasEmission(); }- The Sample function randomly samples a point on the surface of a triangle and returns:
- Information about this point (including position and normal vector)
- Probability density function value (pdf) of the sampling point
- The getArea function returns the area of the triangle.
- The hasEmit function checks whether the triangle's material emits light.
MeshTriangle
The MeshTriangle class is also a subclass of the Object class. It represents a 3D model or mesh composed of many triangles. In other words, a MeshTriangle object may contain many Triangle objects. For example, you can use 12 triangles (2 triangles per face, a total of 6 faces) to represent a cube model.
The MeshTriangle class also includes some additional functions, such as calculating the model's AABB bounding box, BVHAccel objects, etc.
Constructor
The following pseudo code succinctly describes the construction process of MeshTriangle:
MeshTriangle(string filename, Material* mt) { // 1. Load the model file loader.LoadFile(filename); // 2. Create a Triangle object for each face and store for (each face in model) { Triangle tri = create triangles (face vertices, mt); triangles.push_back(tri); } // 3. Calculate the bounding box and total area of the model bounding_box = calculate bounding box (all vertices of the model); area = calculate total area (all triangles); // 4. Create a BVH for fast intersection testing bvh = create BVH (all triangles); }The constructor accepts a file name filename and a material mt, and then uses objl::Loader to load the 3D model. After loading the model, it traverses all the triangles of the model and creates corresponding Triangle objects. At the same time, the bounding box of the entire model and the total area of all triangles are calculated and stored.
In the framework, we use the objl::Loader class to read the .obj file. Call loader.LoadFile(filename) to complete the loading. Then access loader.LoadedMeshes to obtain the loaded 3D model data.
objl::Loader loader; loader.LoadFile(filename); ... auto mesh = loader.LoadedMeshes[0];When loading, we noticed an assertion, which checks whether a model has only one mesh. If there are multiple meshes or no mesh, an assertion error is triggered.
assert(loader.LoadedMeshes.size() == 1);Initializes the minimum and maximum values of the model's vertices, which are used to calculate the axis-aligned bounding box of the 3D model.
Vector3f min_vert = Vector3f{std::numeric_limits ::infinity(), std::numeric_limits ::infinity(), std::numeric_limits ::infinity()}; Vector3f max_vert = Vector3f{-std::numeric_limits ::infinity(), -std::numeric_limits ::infinity(), -std::numeric_limits ::infinity()};Next we need to understand the data structure of the mesh in the objl library, that is, what is stored in the objl::Loader loader.
- MeshName: stores the name of the mesh
- Vertices: Stores all vertex data in the model, including position, normal and texture coordinates.
- Indices: Stores all face (usually triangle) data in the model. Each face consists of a set of indices pointing to vertices in Vertices.
- MeshMaterial: Stores all material data in the model, including diffuse color, specular color, texture and other properties.
The vertex information of each triangle is stored continuously in Vertices, so we construct a Triangle with three vertices as a group, and then set the AABB.
for (int i = 0; i < mesh.Vertices.size(); i += 3) { std::array face_vertices; for (int j = 0; j < 3; j++) { auto vert = Vector3f(mesh.Vertices[i + j].Position.X, mesh.Vertices[i + j].Position.Y, mesh.Vertices [i + j].Position.Z); face_vertices[j] = vert; min_vert = Vector3f(std::min(min_vert.x, vert.x), std::min(min_vert.y, vert.y), std::min(min_vert.z, vert.z)); max_vert = Vector3f(std::max(max_vert.x, vert.x ), std::max(max_vert.y, vert.y), std::max(max_vert.z, vert.z)); } triangles.emplace_back(face_vertices[0], face_vertices[1], face_vertices[2], mt); } bounding_box = Bounds3(min_vert, max_vert);Finally, the areas of all triangles are calculated and the BVH acceleration structure is constructed. In the loop, the code first stores the pointers of all triangles into ptrs, then calculates the sum of the areas of all triangles. Then all triangle pointers are passed into the BVH constructor.
std::vectorMesh triangle related operations
1. Calculation of patch attributes
First isCalculation of patch attributes — getSurfaceProperties, where the following properties need to be calculated:
- The normal vector of a triangleN: This is easy to do. Just find the two sides of the triangle and do a cross product.
- Texture Coordinatesst: The texture coordinates of the triangle vertices are interpolated, and uv is the barycentric coordinate of the intersection point inside the triangle. The following is a detailed description.
void getSurfaceProperties(const Vector3f& P, const Vector3f& I, const uint32_t& index, const Vector2f& uv, Vector3f& N, Vector2f& st) const{ const Vector3f& v0 = vertices[vertexIndex[index * 3]]; const Vector3f& v1 = vertices[vertexIndex[index * 3 + 1]]; const Vector3f& v2 = vertices[vertexIndex[index * 3 + 2]]; Vector3f e0 = normalize(v1 - v0); Vector3f e1 = normalize(v2 - v1); N = normalize(crossProduct(e0, e1)); const Vector2f& st0 = stCoordinates[vertexIndex[index * 3]]; const Vector2f& st1 = stCoordinates[vertexIndex[index * 3 + 1]]; const Vector2f& st2 = stCoordinates[vertexIndex[index * 3 + 2]]; st = st0 * (1 - uv.x - uv.y) + st1 * uv.x + st2 * uv.y; }We have also seen this function in Triangle, but it is different in MeshTriangle.
As their names imply, Triangle represents an independent triangle, while MeshTriangle represents a group of interconnected triangles, that is, a triangular mesh. The getSurfaceProperties of Triangle directly uses the vertex and texture coordinate information stored in the Triangle class. Since our MeshTriangle has multiple triangles, we know the triangle to be obtained through the parameter index.
For the calculation of texture coordinates, first we know that UV coordinates are used in the process of mapping 2D textures onto 3D models. Use the uv coordinates to weight the st coordinates of the triangle vertices to obtain the st coordinates of the intersection. This is called interpolation. These coordinates define the corresponding position of each vertex of the 3D model on the 2D texture.
In this function,
- st0, st1, and st2 are the texture coordinates corresponding to the triangle vertices. These coordinates specify the position of the vertices in the texture map.
- 1 – uv.x – uv.y, uv.x and uv.y correspond to the weights of the three vertices of the triangle. In other words, if you are at a vertex of the triangle, the weight of that vertex is 1, and the weights of other vertices are 0.
- In addition, stCoordinates is defined in advance by the artists, so programmers don’t need to worry about it.
According to the Möller Trumbore algorithm, it is determined that st0 corresponds to (1 – uv.x – uv.y), st1 corresponds to uv.x, and st2 corresponds to uv.y.
To summarize what this function does: multiply the texture coordinate st of each vertex by its corresponding weight, and then add them together. This is an interpolation method that can be used to find the texture coordinates of any point within the triangle.
2. The diffuse color of the uv coordinate on a specific material
The evalDiffuseColor function is used to calculate the diffuse color of a given 2D texture coordinate on a specific material.
When light hits an object's surface, it will reflect according to certain rules, which are affected by the material of the object's surface. Diffuse color is a way to describe this reflection effect, which represents the object's ability to reflect light.
The effect diagram when the pattern item is set to 0, default and 1 respectively:
Vector3f evalDiffuseColor(const Vector2f& st) const { float scale = 5; float pattern = (fmodf(st.x * scale, 1) > 0.5) ^ (fmodf(st.y * scale, 1) > 0.5); return lerp( Vector3f(0.815, 0.235, 0.031), Vector3f(0.937, 0.937, 0.231), pattern); }Collision point information class – Intersection
Collision information structure
This class is used to store information about the intersection of light and objects. I have written the function of each item in the comments below for your reference.
struct Intersection { Intersection(){ happened=false; coords=Vector3f(); normal=Vector3f(); distance= std::numeric_limits ::max(); obj =nullptr; m=nullptr; } // happened indicates whether an intersection actually occurred. // If the ray did not hit any object, then happened will be false. bool happened; // coords Represents the coordinates of the intersection. //If happened is true, then coords will contain the exact position where the light intersects the object. Vector3f coords; // coords represents the texture coordinates of the intersection. //It is used to get the texture of the object surface at the intersection position Information. Vector3f tcoords; // normal represents the normal vector at the intersection. // The normal vector is a vector perpendicular to the surface of an object. It is used to determine the orientation of the object and plays a key role in lighting calculations. Vector3f normal; // emit Indicates the emission value of the light source at the intersection. //If the object at the intersection is a light source, this vector is non-zero. Vector3f emit; //Indicates the distance from the origin of the light to the intersection. double distance; //Points to the object that the light hits . Object* obj; // Points to the material of the object at the intersection, including the object's color, smoothness, reflectivity, etc. Material* m; };Get collision information
This function is actually in Triangle, Sphere and BVHAccel, but the function is inseparable from the Intersection structure. At the same time, for the sake of typesetting, it is simply placed in this chapter.
Intersection in Triangle class
First, paste the source code directly:
inline Intersection Triangle::getIntersection(Ray ray) { Intersection inter; if (dotProduct(ray.direction, normal) > 0) return inter; double u, v, t_tmp = 0; Vector3f pvec = crossProduct(ray.direction, e2) ; double det = dotProduct(e1, pvec); if (fabs(det) < EPSILON) return inter; double det_inv = 1. / det; Vector3f tvec = ray.origin - v0; u = dotProduct(tvec, pvec) * det_inv; if (u < 0 || u > 1) return inter; Vector3f qvec = crossProduct(tvec, e1) ; v = dotProduct(ray.direction, qvec) * det_inv; if (v < 0 || u + v > 1) return inter; t_tmp = dotProduct(e2, qvec) * det_inv; if (t_tmp < 0) { return inter; } inter.distance = t_tmp; inter.coords = ray(t_tmp); inter.happened = true; inter.m = m; inter.normal = normal; inter .obj = this; return inter; }Renderer Class – Renderer
The operation process of the Renderer is very simple: it loops to generate an image for each pixel on the screen.
Here are some simple instructions:
- spp (samples per pixel): The number of samples per pixel, which indicates how many rays the ray tracing algorithm will cast in one pixel.
- framebuffer: A one-dimensional array of size width*height that stores the color of each pixel in the scene.
- scene.castRay(Ray(eye_pos, dir), 0) function: Calls the ray tracing algorithm. This function will shoot a ray and calculate the color of the object that this ray hits in the scene. This color value is then added to the color of the corresponding pixel in the framebuffer.
So the key point of rendering the Renderer class is this line castRay. And castRay is in the Scene class, which is also the focus of this article.
Some notes on intersect
We see a lot of intersect functions in the project, which makes people dizzy. For example, when we first see intersect(), it directly returns true in the triangle class. But in fact, we will wonder, shouldn’t it contain judgment logic? Instead of directly returning 0 or 1. There are also Intersect() in BVHAccel, IntersectP() in Bounds3, and the relationship between getIntersection of specific objects, etc.
In this regard, the process is briefly explained.
According to the process, first in the castRay() function, we call BVHAccel's Intersect(), and then BVHAccel's Intersect() will find the AABB of the leaf node in the BVH data structure. Then it will call Bounds3's IntersectP to determine whether the bounding box of the node intersects with the light.
- If they do not intersect: return the default constructor of the Intersection class (empty collision data structure);
- If the current node is a leaf node: directly call the getIntersection method of the object to calculate the intersection of the light and the object.
Next, the getIntersection method will return the Intersection data structure corresponding to their object type. The following are the getIntersection methods of Sphere, MeshTriangle and Triangle respectively.
Intersection getIntersection(Ray ray){ Intersection result; ... return result; } Intersection getIntersection(Ray ray){ Intersection intersec; if (bvh) intersec = bvh->Intersect(ray); return intersec; }The getIntersection method of Triangle is marked with override, which means that the function declaration is inside the class definition, while the function definition is outside the class definition. Therefore, the part that is actually executed is inline:
inline Intersection Triangle::getIntersection(Ray ray){ ... }Scene
Before talking about castRay, let's take a quick look at the general shape of the Scene.
This class contains all the objects and lights in the scene, as well as some rendering parameters, such as the scene's width, height, field of view, background color, etc. Next, we will explain the main parts of this class:
Member variables:
- width and height are the pixel width and height of the scene, fov is the field of view of the camera, and backgroundColor is the background color of the scene.
- objects and lights store the objects and light sources in the scene respectively. objects is a vector of pointers to Object type, and lights is a vector containing smart pointers\ type.
- bvh is a pointer to BVHAccel, which is used to store the scene's bounding volume hierarchy (BVH) structure to speed up the intersection calculation between rays and objects.
- maxDepth and RussianRoulette are used to control the details of the path tracing algorithm (we will come back to this in a moment).
Member functions:
- The Add function is used to add objects or light sources to the scene.
- The HandleAreaLight function is used to handle ray tracing of area light sources.
- The reflect and refract functions are used to calculate the direction of reflection and refraction of light, respectively. The refract function implements Snell's law and is used to calculate the direction of refraction of light when it enters different media. Two situations need to be handled specially: light is emitted from the inside of the object, and light is emitted from the outside of the object.
- The fresnel function is used to calculate the Fresnel equation and obtain the ratio of reflected light to transmitted light. That is, the ratio of reflection to transmission when light passes through the interface of two different media.
- The functions of intersect, buildBVH, castRay, sampleLight and trace are described in detail later.
So far, this Scene class encapsulates all the data and operations of a ray tracing rendering scene. Next, we will talk about some details of this class in detail.
In Scene.cpp, we first have the buildBVH() function, which creates a bounding volume hierarchy (BVH).
The intersect(const Ray &ray) const function is then used to find the intersection of the ray defined by ray with any object in the scene. It does this for efficiency by using the previously constructed BVH.
sampleLight(Intersection &pos, float &pdf) const The function is used to randomly sample the luminous objects in the scene. This function first calculates the sum of the areas of all luminous objects, and then randomly selects a position in this area as a sampling point.
void Scene::sampleLight(Intersection &pos, float &pdf) const{ float emit_area_sum = 0; for (uint32_t k = 0; k < objects.size(); ++k) { if (objects[k]->hasEmit() ){ emit_area_sum += objects[k]->getArea(); } } float p = get_random_float() * emit_area_sum; emit_area_sum = 0; for (uint32_t k = 0; k < objects.size(); ++k) { if (objects[k]->hasEmit()){ emit_area_sum += objects[k]->getArea( ); if (p <= emit_area_sum){ objects[k]->Sample(pos, pdf); break; } } } }enter: Intersection refers to pos, and floating point numbers refer to pdf.
These two parameters are output parameters, that is, this method will change their values. The object of type Intersection is used to store the intersection information between the light and the object, and pdf represents the probability density of selecting this intersection.
Output: This method has no return value, its result is returned by changing the values of pos and pdf.
Specific process:
- First, this method calculates the sum of the areas of all emitting objects, emit_area_sum, through the first loop. For each object, it first calls the hasEmit method to check whether the object can emit light. If it can, it calls the getArea method to get the area of the object and adds it to emit_area_sum.
- Next, it generates a random number p between 0 and emit_area_sum. This random number p is used to randomly select an emitter object.
- Then, it goes through the second loop to select the emitter. In this loop, it first checks whether each object can emit light. If it can, it adds the area of the object and then checks whether p is less than or equal to the current emit_area_sum. If so, then this object is the selected object. Then, it calls the Sample method of this object, randomly selects a point on this object, and updates the values of pos and pdf. Finally, it jumps out of the loop and ends this method.
Path Tracing Implementation – castRay()
Finally, we come to castRay. We are going to use the Path Tracing algorithm to implement the ray tracing function.
// Implementation of Path Tracing Vector3f Scene::castRay(const Ray &ray, int depth) const { // TO DO Implement Path Tracing Algorithm here }Next, let's introduce the key processes: light source sampling and indirect lighting calculation.
In other words, we can roughly divide the PT algorithm into two parts: direct illumination and indirect illumination.
First some basic definitions:
- p: The collision point is the intersection of the light and the object in the scene.
- wo: outgoing direction, which is the direction reflected from the intersection point to the camera.
- wi: incoming direction, which is the direction in which the light is reflected from the light source to the intersection point.
The following is pseudo code that details the above process:
shade(p, wo) // First calculate the direct light source sampleLight(inter, pdf_light) Get x, ws, NN, emit from inter Shoot a ray from p to x If the ray is not blocked in the middle L_dir = emit * eval( wo, ws, N) * dot(ws, N) * dot(ws, NN) / |xp|^2 / pdf_light // Then calculate the indirect light source L_indir = 0.0 Test Russian Roulette with probability RussianRoulette wi = sample(wo, N) Trace a ray r(p, wi) If ray r hit a non-emitting object at q L_indir = shade(q, wi) * eval(wo, wi , N) * dot(wi, N) / pdf(wo, wi, N) / RussianRoulette Return L_dir + L_indirMake sure you have a clear understanding of how Path Tracing is implemented, then start writing code.
Vector3f Scene::castRay(const Ray &ray, int depth) const { const float EPLISON = 0.0001f; Intersection p_inter = intersect(ray); if (!p_inter.happened) return Vector3f();// Default construction zero vector - blackif (p_inter.m->hasEmission())// Is it self-luminousreturn p_inter.m->getEmission();// Directly return the emission color// Get the intersection between ray and object plane Intersection x_inter; float pdf_light = 0.0f; // Sample light source at intersection point sampleLight(x_inter, pdf_light); // Get x, ws, N, NN, emit from inter Vector3f p = p_inter.coords;// Object intersection coordinates Vector3f x = x_inter.coords;// Light source intersection coordinates Vector3f ws_dir = (x - p).normalized();// float ws_distance = (x - p).norm(); // Distance from object to light source Vector3f N = p_inter.normal.normalized(); // Normal of object intersection Vector3f NN = x_inter.normal.normalized(); // Normal of light source intersection Vector3f emit = x_inter.emit; // Color vector of light source intersection // Shoot a ray from p to x Vector3f l_dir(0.0f), l_indir(0.0f); // See the following description Ray ws_ray(p, ws_dir); // Make a Ray Intersection from p to the light point ws_ray_inter = intersect(ws_ray); // Then intersect // If the ray is not blocked in the middle // That is, check whether the straight line path from object p to light source x is blocked by other objects if(ws_ray_inter.distance - ws_distance > -EPLISON) { // See below for details l_dir = emit * p_inter.m->eval(ray.direction, ws_ray.direction, N) * dotProduct(ws_ray.direction, N) * dotProduct(-ws_ray.direction, NN) / (ws_distance * ws_distance) / pdf_light; } // Test Russian Roulette with probability RussianRoulette if(get_random_float() <= RussianRoulette) { // See below for details Vector3f wi_dir = p_inter.m->sample(ray.direction, N).normalized(); Ray wi_ray(p_inter.coords, wi_dir); // If ray r hit a non-emitting object at q Intersection wi_inter = intersect(wi_ray); // If a collision is detected and the collision point is not emitting, start calculating indirect lighting if (wi_inter.happened && (!wi_inter.m->hasEmission())) { // See the description below for details l_indir = castRay(wi_ray, depth + 1) * p_inter.m->eval(ray.direction, wi_ray.direction, N) * dotProduct(wi_ray.direction, N) / p_inter.m->pdf(ray.direction, wi_ray.direction, N) / RussianRoulette; } } return l_dir + l_indir; }A few points need to be explained:
- ray or wo_ray: This is the ray we are dealing with, called the outgoing ray, which means the ray that propagates from a certain point to a certain direction. In the castRay() function, ray is the ray passed as a parameter.
- intersect: will call the overloaded method located in bvh, as follows:
c++ Intersection BVHAccel::Intersect(const Ray& ray) const { Intersection isect; if (!root) return isect; isect = BVHAccel::getIntersection(root, ray); return isect; }
- p_inter: This isLight and SurfaceSpecifically, p_inter is an object of type Intersection.
- x_inter: This isLight and light sourcesThe intersection information. Similar to p_inter, x_inter is also an object of Intersection type.
- l_dir (direct illumination): This is the lighting contribution from the light source directly hitting the surface. In path tracing, this is typically obtained by sampling the light source from the surface point and calculating the direct illumination.
- l_indir (indirect illumination): This is the indirect lighting contribution reflected by the environment onto the surface. In path tracing, this is typically calculated by sampling the BRDF of the surface and recursively tracing reflected rays.
- ws_ray: This is the ray from point p to point x. When calculating direct lighting, a new ray needs to be sent towards the light source to check if the object can directly see the light source (i.e. there are no other objects in the way). This ray is ws_ray.
Here we will focus on the BRDF formula for direct lighting. The formula here takes into account the geometry term and the PDF of the light source sampling, so the formula is the same as mine.Previous Articles The formula given in is different.
The mathematical formula for direct lighting is as follows: Radiance under direct lightingRadiance emitted by the light sourceCosine of the angle between the incident light ray and the surface normalCosine of the angle between the outgoing light ray and the normal vector of the light sourceAttenuation functionProbability density function for selecting the direction of the light sourceLdir=Li⋅fr(wi,wo,N)⋅cos(θ)⋅cos(θ′)r2⋅p(wi)where,Ldir: radiance under direct lightingLi: radiance emitted by the light sourcefr(wi,wo,N):BRDFcos(θ): cosine of the angle between the incident light ray wi and the surface normal Ncos(θ′): cosine of the angle between the outgoing light ray −wsdir and the normal vector NN of the light sourcer2: Attenuation functionp(wi): Probability density function for selecting the direction of the light source (PDF) Converted into code is:
l_dir = emit * p_inter.m->eval(ray.direction, ws_ray.direction, N) * dotProduct(ws_ray.direction, N) * dotProduct(-ws_ray.direction, NN) / (ws_distance * ws_distance) / pdf_light;The eval function is as follows, which describes how the light intensity changes with the change of the direction of the incident light and the outgoing light. The input of this function is the incident direction of the light ray.direction, the reflection direction ws_ray.direction, and the surface normal N, and the output is a value that measures the intensity of the reflected light. In the code, diffuse BRDF = ρ/π.
... float cosalpha = dotProduct(N, wo); if (cosalpha > 0.0f) { Vector3f diffuse = Kd / M_PI; return diffuse; } else return Vector3f(0.0f); ...After calculating direct illumination, we start to calculate indirect illumination using Russian roulette. The process is as follows:
In Scene.hpp we define the value of RussianRoulette, and we will calculate indirect lighting only when the random number obtained is less than RussianRoulette. In other words, it is possible that the indirect lighting of a pixel is not calculated at all.
If we need to calculate indirect lighting now, we generate a new ray wi_ray for it. The direction of this ray is sampled based on the surface material of the current intersection and the original ray direction. This part implements material-based importance sampling. In the code, we doUniform SamplingThe specific calculation method is shown in the following code:
... // uniform sample on the hemisphere float x_1 = get_random_float(), x_2 = get_random_float(); float z = std::fabs(1.0f - 2.0f * x_1); // r - distance from the origin to the hemisphere; phi - angle in polar coordinates float r = std::sqrt(1.0f - z * z), phi = 2 * M_PI * x_2; Vector3f localRay(r*std::cos(phi), r*std::sin(phi), z); return toWorld(localRay, N); ...Finally, let's talk about the calculation of indirect lighting l_indir. The calculation formula for indirect lighting is: Lo=shade(q,−wi)⋅fr⋅cosine⋅1pdf(wi) But here
l_indir = castRay(wi_ray, depth + 1) * p_inter.m->eval(ray.direction, wi_ray.direction, N) * dotProduct(wi_ray.direction, N) / p_inter.m->pdf(ray.direction, wi_ray .direction, N) / RussianRoulette;castRay(wi_ray, depth + 1): This is a recursive call, which represents emitting a new ray from the current intersection point and obtaining the lighting contribution generated by the reflection of the ray on all objects.
- Why divide by RussianRoulette at the end?
Because when using Russian Roulette, we randomly generate a number, and if this number is greater than a threshold (here is the RussianRoulette variable), we terminate the ray tracing. When we terminate the tracing of a certain ray path, we are actually abandoning all possible subsequent reflections of this ray path, and these abandoned reflections may contribute to the final ray information.
Therefore, to compensate for the effect of terminating the trace, we amplify the intensity of the remaining rays by a factor of 1 / RussianRoulette. This ensures that the expected values of the intensity of all retained ray paths are equal to their actual intensities, thus ensuring the unbiased nature of the ray tracing algorithm. I will discuss this in myPrevious article It is also mentioned in.
At this point, we have basically explained the entire project.Previous Articles There are involved.
Finally, we achieved the following rendering:
Next, I will talk about how to use multi-threaded computing to speed up rendering.
Multithreaded acceleration
A quick introduction to C++ multithreading
For readers who are not familiar with C++ multithreading, here is an example to help you get started quickly.
#include #include // This is the function we will run in both threads void printMessage(std::string message) { std::cout << message << std::endl; } int main() { // Create and run two threads std::thread thread1(printMessage, "Hello from thread 1"); std::thread thread2(printMessage, "Hello from thread 2"); // Wait for both threads to finish thread1.join(); thread2.join(); return 0; }Before our main function ends, we need to wait for all threads to end, so we must call the join() function. If the thread does not end, but the main program ends early, this may cause unexpected results.
Deploy multithreading
A relatively simple method is to use OpenMP, which will not be described here.
The other is the most direct method, manual block + mutex, the code is as follows:
// change the spp value to change sample ammount int spp = 4; int thread_num = 16; int thread_height = scene.height / thread_num; std::vector threads(thread_num); std::cout << "SPP: " << spp << "\n"; std::mutex mtx; float process=0; float Reciprocal_Scene_height=1.f/ (float)scene.height; auto castRay = [&](int thread_index) { int height = thread_height * (thread_index + 1); for (uint32_t j = height - thread_height; j < height; j++) { for (uint32_t i = 0; i < scene.width; ++i) { // generate primary ray direction float x = (2 * (i + 0.5) / ( float)scene.width - 1) * imageAspectRatio * scale; float y = (1 - 2 * (j + 0.5) / (float)scene.height) * scale; Vector3f dir = normalize(Vector3f(-x, y, 1)); for (int k = 0; k < spp; k++) framebuffer[j*scene.width+i] += scene.castRay(Ray(eye_pos, dir), 0) / spp; } mtx.lock( ); process = process + Reciprocal_Scene_height; UpdateProgress(process); mtx.unlock(); } }; for (int k = 0; k < thread_num; k++){ threads[k] = std::thread(castRay,k); std::cout << "Thread[" << k << "] Started:" << threads[k] .get_id() << "\n"; } for (int k = 0; k < thread_num; k++){ threads[k].join(); } UpdateProgress(1.f);But surprisingly, on my macOS, the speed of single thread and multi-thread (8t) is the same. My configuration is as follows:
- Model Name: MacBook Pro – 14-inch, 2021
- Chip: Apple M1 Pro
- Total Number of Cores: 8 (6 performance and 2 efficiency)
- Memory: 16G
- System Version: macOS 13.2.1 (22D68)
16.threads vs 1.thread
Optimize progress display
Let each thread track its own progress and find that each thread is indeed calculating the corresponding area synchronously.
auto castRay = [&](int thread_index) { // Create a local progress variable for this thread float thread_progress = 0.f; int height = thread_height * (thread_index + 1); for (uint32_t j = height - thread_height; j < height; j++) { for (uint32_t i = 0; i < scene.width; ++i) { // ... Rest of your code ... } // Update this thread's progress thread_progress += Reciprocal_Scene_height; std::cout << "Thread[" << thread_index << "] Progress: " << thread_progress << "\n"; } };
After testing, using 2 threads on my device is the best speed, which is about 70% faster. This statement is also confirmed in my task manager (the program sets the number of threads to 2):
It is speculated that the scheduling strategy of the M1pro chip is that when encountering a large number of parallel calculations, it will prioritize starting only two large cores to collaborate, rather than starting all cores like Windows (generally speaking).
Microfacet Models
- The microsurface theory
As early as 1967, physicists Torrance and Sparrow proposed the theory of microfacet models in their paper "Theory for Off-Specular Reflection From Roughened Surfaces".
The application in computer graphics is attributed to the work of Robert L. Cook and Kenneth E. Torrance. In their 1982 paper "A Reflectance Model for Computer Graphics", they introduced the microsurface model to the field of computer graphics, providing a physically based method for simulating various materials in the real world.
- What is a microsurface model?
The microsurface model is a theoretical model used in computer graphics to simulate the reflection and refraction of light from rough surfaces. The basic assumption of this model is that a complex microsurface can be replaced by a simplified macrosurface, and the scattering function (BSDF) of the macrosurface can match the overall directional scattering behavior of the microsurface. In other words, the microsurface model attempts to simulate how light scatters on a rough surface through statistical methods, rather than trying to accurately simulate the details of each microsurface.
The microsurface model provides a concrete, physically based form for the BRDF, which is a core component of the lighting model.
- Why the Microsurface Model?
To quote the paper "A Reflectance Model for Computer Graphics":
why images rendered with previous models often look plastic
In other words, the model proposedOne of themThe purpose is to solve the "plastic feeling".
definition
In the previous shading model, we have the formula for the total reflectance r: r=ka+∑i=0nlc⋅(rd+rs)With,rd=kd⋅(n→⋅l→)rs=ks⋅(h→⋅n→)p $r_d, r_s$ represents the diffuse reflectance and specular reflectance respectively.
Cook-Torrance provides a variation of this equation: r=ka+∑i=0nlc⋅(n→⋅l→)⋅(d⋅rd+s⋅rs) This equation introduces two new concepts:
- The diffuse reflection term $r_d$ and the specular reflection term $r_s$ are controlled by the variables d and s respectively.
- The dot product of the normal vector and the light vector is separated from the diffuse reflectance, making it part of the sum. This allows the diffuse reflectance to be a constant, and this definition can be modified later if necessary.
In addition, the specular reflectivity $r_s$ contains a divisor of $(n⋅l)$, which cancels out the $(n⋅l)$ in the sum, which may be very confusing to the reader. Let me explain why.
- About $s$ and $d$
Since $s$ and $d$ are used to control the balance between diffuse and specular terms, there is the following relationship between the two: s+d=1 Therefore, we generally ignore $d$ and simplify the Cook-Torrance formula: r=ka+∑i=0nlc⋅(n→⋅l→)⋅((1−s)⋅rd+s⋅rs) It is not difficult to understand the relationship between $s$ and $d$: According to the law of conservation of energy, the total light energy reflected by a material (including diffuse and specular reflection) cannot exceed the light energy it absorbs. Therefore, if a material reflects a lot of light energy, then it must reflect less light energy diffusely, and vice versa. For example, a very smooth surface (such as a mirror or metal) will have high specular reflection but almost no diffuse reflection. A very rough surface (such as a brick wall or cloth) will have high diffuse reflection but almost no specular reflection.
- About $r_s$
rs=D*G*F4*(n→⋅l→)*(n→⋅v→)
Note that the literature and other places sometimes use $π$ in the denominator of the formula instead of $4$. This is an error that was made in the Cook-Torrance paper and has been repeated in many places. The correct derivation of the formula should use $4$. However, this is a small difference in a constant factor so it is not particularly important if not done exactly.
In the Cook-Torrance specular reflectance (r_s) formula, D, G, and F are three functions that can be selected in different forms. They represent the distribution function (Distribution), the geometric attenuation function (Geometry), and the Fresnel reflection function (Fresnel).
Microsurface
The microfacet model considers the surface of an object as composed of countless tiny face elements, each of which can have its own normal. Such a model can better simulate the details of the surface of an object, including roughness and smoothness.
The surface normal and the normal of each microfacet may not be the same. For very smooth surfaces, such as a perfect mirror, all microfacets face the same direction, which is the surface normal $n$. However, for matte or rough surfaces, such as matte paint or stone, the directions faced by microfacets are random. Or, on rough surfaces, these microfacets may occlude other facets and cast shadows on the facets.
Cook-Torrance Model
The Cook-Torrance model attempts to explain the above three phenomena:
- The distribution of microfacets (via the Normal Distribution Function D): Different microfacets face different directions, depending on the roughness of the surface. This distribution describes the proportion of microfacets that reflect light for an observer.
- Mutual occlusion and shadow effects between microfacets (through geometric attenuation function G): This function describes the mutual occlusion and shadow effects between microfacets. On a rough surface, some microfacets may be occluded by other facets or cast shadows on other facets.
- The interaction between light and microfacets (through the Fresnel reflection function F): This determines how much light will be reflected, absorbed, or transmitted after it hits a microfacet.
Although the contribution of G and F functions to the final rendering result may be small, they are still important components of the model. In the physically based BRDF model, the distribution function D is a very critical part because it determines how light is reflected from various parts of the object's surface, thus affecting the visual effect of the object. Commonly used distribution functions include Beckmann distribution, GGX distribution, etc. Choosing different distribution functions can simulate the visual effects of different types of materials.
Fresnel term
This part is the Cook-Torrance model used to explain the Fresnel effect.
In graphics, the Schlick approximation is generally used: F=F0+(1−F0)∗(1−(v→⋅h→))5 where
- $F_0$ is the reflectivity when the light is incident vertically, which is calculated from the refractive index n of the object
- v and h are the view direction and half-distance vector (the average direction of the incident and reflected rays) respectively.
This formula assumes that the reflectivity changes linearly with the incident angle, but the reflectivity increases faster in the edge area (where the incident angle is close to 90 degrees), so a fifth power term is added to better simulate this phenomenon.
Finally, this formula can be implemented in C++ as follows:
double Fresnel_Schlick(double n, const Vector3d &v, const Vector3d &h) { double F0 = ((n - 1) * (n - 1)) / ((n + 1) * (n + 1)); double cosTheta = std ::max(dot(v, h), 0.0); return F0 + (1.0 - F0) * std::pow(1.0 - cosTheta, 5); }Geometric attenuation term
In the Cook-Torrance model, geometric attenuation is usually defined as the minimum of two attenuation factors, corresponding to the two effects. These two attenuation factors can be calculated by some simple formulas, called the Torrance-Sparrow geometric shielding model. G=min(1,2(h→⋅n→)(n→⋅v→)(v→⋅h→),2(h→⋅n→)(n→⋅l→)(v→⋅h→)) where,
- V is the viewing direction
- H is the half vector (average of the light direction and the view direction)
- N is the surface normal
- n_dot_v and n_dot_l are the dot products of the normal with the view direction and the light direction respectively.
This model is more suitable when the surface roughness is low (that is, the surface is closer to complete mirror reflection), and can be expressed in C++ as follows:
double v_dot_h = dotProduct(V, H); double n_dot_h = dotProduct(N, H); double G1 = 2 * n_dot_h * n_dot_v / v_dot_h; double G2 = 2 * n_dot_h * n_dot_l / v_dot_h; double G = clamp(0, 1, std::min(G1, G2)); // Note that the clamp method is overridden in the framework, so there is no error here~D Function
Finally, let's talk about the most critical and most contributing item, the microscopic surface item of specular reflection.
The specular part of the Blinn-Phong model is one possible choice of D function that will make the specular part have Blinn-Phong-like properties, but it is not the only choice. In fact, there are many possible NDFs, each with different surface roughness and reflectivity properties.
The most commonly used NDF is probably the GGX (Generalized Trowbridge-Reitz) model, which works well for a wide range of smooth to rough surfaces. Another common choice is the Beckmann model, which generally performs better for surfaces with medium to high roughness.
Here we use the isotropic Beckmann model: D = e−tan2αm2πm2cos4α where,
- $m$ represents the parameter of surface roughness
- $\alpha$ is the angle between the half-vector $H$ and the surface normal $N$, calculated using $acos(\text{n_dot_h})$.
double m = (type == MaterialType::MICROFACET_DIFFUSE) ? 0.6 : 0.2; double alpha = acos(n_dot_h); double D = exp(-pow(tan(alpha)/m, 2)) / (M_PI*m*m *pow(cos(alpha), 4));As the surface becomes very rough, the orientation of these microfacets becomes more and more random, causing the reflected light to become more scattered.
Cook and Torrance recommended the Beckmann distribution. But now more people useGGX Distribution, I will talk about it in depth later.
Code implementation
The Cook-Torrance formula we need to implement is: r=ka+∑i=0nlc⋅(n→⋅l→)⋅((1−s)⋅rd+s⋅rs) We have just given the specific codes for the three terms of DGF. We can directly write the $r_s$ term in the Cook-Torrance formula:
Vector3f Material::cookTorrance(const Vector3f &wi, const Vector3f &wo, const Vector3f &N) { auto V = wo; auto L = wi; auto H = normalize(V + L); auto type = m_type; double n_dot_v = dotProduct(N , V); double n_dot_l = dotProduct(N, L); if(!(n_dot_v > 0 && n_dot_l > 0)) return 0; double n_air = 1, n_diff = 1.2, n_glos = 1.2; double n2 = (type == MaterialType::MICROFACET_DIFFUSE) ? n_diff : n_glos; double r0 = (n_air-n2)/(n_air+ n2); r0*=r0; double F = r0+(1-r0)*pow(1 - n_dot_v, 5); double v_dot_h = dotProduct(V, H); double n_dot_h = dotProduct(N, H); double G1 = 2 * n_dot_h * n_dot_v / v_dot_h; double G2 = 2 * n_dot_h * n_dot_l / v_dot_h; double G = clamp(0, 1, std::min(G1, G2)); double m = (type == MaterialType::MICROFACET_DIFFUSE) ? 0.6 : 0.2; double alpha = acos(n_dot_h); double D = exp(-pow(tan(alpha)/m, 2)) / (M_PI*m*m*pow (cos(alpha), 4)); auto ans = F * G * D / (n_dot_l * n_dot_v * 4); return ans; }$r_d $ directly uses $\frac{1}{\pi}$. This assumption is based on the properties of Lambertian reflection. Lambertian reflection is an ideal, completely diffuse reflection surface.
Then $k_a$ represents the reflection coefficient of ambient light. Ambient light is used to simulate the global illumination in the scene. Here, ambient light is directly ignored, that is, $k_a=0$.
case MICROFACET_DIFFUSE: { float cosalpha = dotProduct(N, wo); if (cosalpha > 0.0f) { auto ans = Ks * cookTorrance(wi, wo, N) + Kd * eval_diffuse(wi, wo, N); return ans; } else return Vector3f(0.0f); break; }Additionally, the two coefficients for the material are set as follows:
Material *white_m = new Material(MICROFACET_DIFFUSE, Vector3f(0.0f)); white_m->Kd = Vector3f(0.725f, 0.71f, 0.68f); white_m->Ks = Vector3f(1,1,1) - white_m-> Kd;Here is the final rendered result. The ball specifically uses the microsurface model MICROFACET_DIFFUSE, SSP=64, and resolution 784, 784:
Reference
- Fundamentals of Computer Graphics 4th
- GAMES101 Lingqi Yan
- "A Reflectance Model for Computer Graphics" - 357290.357293
- https://zhuanlan.zhihu.com/p/152226698
- https://graphicscompendium.com/gamedev/15-pbr
- "Theory for Off-Specualr Reflection From Roughened Surfaces"
- https://hanecci.hatenadiary.org/entry/20130511/p
EN 
CN
Leave a Reply