Project 3: Li Zheng

README.md

@@ -3,11 +3,29 @@ CUDA Path Tracer
 
 **University of Pennsylvania, CIS 565: GPU Programming and Architecture, Project 3**
 
-* (TODO) YOUR NAME HERE
-* Tested on: (TODO) Windows 22, i7-2222 @ 2.22GHz 22GB, GTX 222 222MB (Moore 2222 Lab)
+* Li Zheng
+  * [LinkedIn](https://www.linkedin.com/in/li-zheng-1955ba169)
+* Tested on: Windows CUDA10, i5-3600 @ 3.59GHz 16GB, RTX 2060 6GB (personal computer)
 
-### (TODO: Your README)
+![render_cornell](img/part1.PNG)  
 
-*DO NOT* leave the README to the last minute! It is a crucial part of the
-project, and we will not be able to grade you without a good README.
+## Part 1 - Core Features
+In this part, I implete
+- A shading kernel with BSDF evaluation for ideal Diffuse surfaces and perfectly specular-reflective surfaces.
+- Path continuation/termination using Stream Compaction from Project 2.
+- Toggleable meaterial sort to make rays/pathSegments/intersections contiguous in memory.
+- Toggle first bounce cache for re-use across all subsequent iterations
 
+### Sort Material
+Test on iterations of 5000 and depth of 8 with first bounce cache. The rendering with material sort takes 402.3s compared with 403.9s without material sort. The material sort is more efficient, because different materials take different time to complete. Compared with reflective material, diffusive material takes more time to complete BSDF evaluations. It’s better to make path segments continuous in memory by material type. Then the threads of the same warp are more likely to finish at the same time, avoiding the wait.
+
+### First Bounce Cache
+![cache_first_bounce](img/cache_first_bounce.PNG)  
+The diagram demonstrates the time to complete 5000 iterations of path tracing with different max depths. With the max depth increases, the time consumption of first-bounce-cache method is less than the case without the cache. It is predictable since the cache saves the time to calculate the first bounce of each iteration. 
+
+## Part 2 - Advance Features
+Implement refraction and depth-of-field
+![refraction and depth-of-field](img/part2.PNG)  
+
+## Part 3 - Octree 
+Tried to implement Octree, but there are some bugs haven't been fix.

models/12190_Heart_v1_L3.mtl

@@ -0,0 +1,14 @@
+# 3ds Max Wavefront OBJ Exporter v0.97b - (c)2007 guruware
+# File Created: 05.12.2011 16:42:58
+
+newmtl Heart
+	Ns 55.0000
+	Ni 1.5000
+	d 1.0000
+	Tr 0.0000
+	Tf 1.0000 1.0000 1.0000 
+	illum 2
+	Ka 0.6980 0.0510 0.0510
+	Kd 0.6980 0.0510 0.0510
+	Ks 0.4140 0.4140 0.4140
+	Ke 0.0000 0.0000 0.0000

scenes/cornell.txt

@@ -53,7 +53,7 @@ CAMERA
 RES         800 800
 FOVY        45
 ITERATIONS  5000
-DEPTH       8
+DEPTH       4
 FILE        cornell
 EYE         0.0 5 10.5
 LOOKAT      0 5 0
@@ -115,3 +115,12 @@ material 4
 TRANS       -1 4 -1
 ROTAT       0 0 0
 SCALE       3 3 3
+
+// Arrow
+OBJECT 7
+mesh
+../models/12190_Heart_v1_L3.obj
+material 2
+TRANS      -1 1 -5
+ROTAT       0 0 0
+SCALE       .4 .4 .4 

src/interactions.h

@@ -1,6 +1,8 @@
 #pragma once
 
 #include "intersections.h"
+#include "glm/glm.hpp"
+#include "glm/gtx/norm.hpp"
 
 // CHECKITOUT
 /**
@@ -76,4 +78,47 @@ void scatterRay(
     // TODO: implement this.
     // A basic implementation of pure-diffuse shading will just call the
     // calculateRandomDirectionInHemisphere defined above.
+
+	if (m.emittance > 0.0f) {
+		pathSegment.color *= (m.color * m.emittance);
+	}
+	if (pathSegment.remainingBounces <= 0) {
+		return;
+	}
+	// pure-diffusive
+	/*float reflectiveRatio;
+	if (m.hasReflective) {
+		reflectiveRatio = 0.8;
+	}
+	else {
+		reflectiveRatio = 0.2;
+	}
+
+	thrust::uniform_real_distribution<float> u01(0, 1);
+	float randomNum = u01(rng);
+	*/
+	if (!m.hasReflective) {
+		glm::vec3 diffuseDirection = calculateRandomDirectionInHemisphere(normal, rng);
+		pathSegment.ray.direction = diffuseDirection;
+		//float lightTerm = glm::abs(glm::dot(normal, diffuseDirection));
+		pathSegment.color *= m.color;
+	}
+	if (m.hasReflective) {
+		glm::vec3 reflectDirection = glm::reflect(pathSegment.ray.direction, normal); 
+		pathSegment.ray.direction = reflectDirection;
+		float lightTerm = pow(glm::abs(glm::dot(normal, reflectDirection)), m.specular.exponent);
+		pathSegment.color *= (m.color * lightTerm);
+	}
+	if (m.hasRefractive) {
+		float n = m.indexOfRefraction;
+		glm::vec3 refractDirection = glm::refract(pathSegment.ray.direction, normal, n);
+		pathSegment.ray.direction = refractDirection;
+		float F0 = (n - 1) * (n - 1) / (n + 1) / (n + 1);
+		float lightTerm = glm::abs(glm::dot(normal, refractDirection));
+		float F = F0 + (1 - F0) * pow((1 - lightTerm), 5);
+		pathSegment.color *= (m.color * F);
+	}
+
+	pathSegment.ray.origin = intersect + 0.01f * normal;
+	pathSegment.remainingBounces--;
 }

src/intersections.h

@@ -142,3 +142,166 @@ __host__ __device__ float sphereIntersectionTest(Geom sphere, Ray r,
 
     return glm::length(r.origin - intersectionPoint);
 }
+
+__host__ __device__ float meshIntersectionTest(Geom mesh, Triangle *triangles, Ray r,
+	glm::vec3 &intersectionPoint, glm::vec3 &normal, bool &outside) {
+
+	glm::vec3 ro = multiplyMV(mesh.inverseTransform, glm::vec4(r.origin, 1.0f));
+	glm::vec3 rd = glm::normalize(multiplyMV(mesh.inverseTransform, glm::vec4(r.direction, 0.0f)));
+
+	for (int i = mesh.tIndexStart; i <= mesh.tIndexEnd; i++)
+	{
+		Triangle triangle = triangles[i];
+		glm::vec3 p0 = triangle.vertices[0];
+		glm::vec3 p1 = triangle.vertices[1];
+		glm::vec3 p2 = triangle.vertices[2];
+		glm::vec3 n = glm::cross(p1 - p0, p2 - p0); // normal to check intersection inside a triangle, could opposite to triangle.normal 
+		float d = -glm::dot(triangle.normal, p1);
+
+		if (glm::abs(glm::dot(triangle.normal, rd)) < EPSILON) {
+			continue;
+		}
+		float t = (-d - glm::dot(triangle.normal, ro)) / (glm::dot(triangle.normal, rd));
+		if (t < 0) {
+			continue;
+		}
+		glm::vec3 p = ro + t * rd;
+
+		if (glm::dot(glm::cross(p1 - p0, p - p0), n) < 0) {
+			continue;
+		}
+		if (glm::dot(glm::cross(p2 - p1, p - p1), n) < 0) {
+			continue;
+		}
+		if (glm::dot(glm::cross(p0 - p2, p - p2), n) < 0) {
+			continue;
+		}
+		if (glm::dot(rd, triangle.normal) <= 0) {
+			outside = true;
+		}
+		else {
+			outside = false;
+		}
+
+		glm::vec3 objspaceIntersection = p;
+		intersectionPoint = multiplyMV(mesh.transform, glm::vec4(objspaceIntersection, 1.f));
+		normal = glm::normalize(multiplyMV(mesh.invTranspose, glm::vec4(triangle.normal, 0.f)));
+		return glm::length(r.origin - intersectionPoint);
+	}
+	return -1;
+
+}
+
+// Octree
+__host__ __device__ void octreeIntersection(Geom mesh, Triangle *triangles,
+	int *sortTringles, OctreeNode_cuda *octreeVector, int nodeIdx,
+	Ray &r, glm::vec3 &intersectionPoint, glm::vec3 &normal,
+	bool &outside, glm::vec3 &ro, glm::vec3 &rd, float &t) {
+
+	if (nodeIdx == -1) {
+		return;
+	}
+	OctreeNode_cuda root = octreeVector[nodeIdx];
+	float tmin = -1e38f;
+	float tmax = 1e38f;
+	float amin, amax;
+	for (int xyz = 0; xyz < 3; ++xyz) {
+		if (xyz == 0) {
+			amin = root.xmin;
+			amax = root.xmax;
+		}
+		else if (xyz == 1) {
+			amin = root.ymin;
+			amax = root.ymax;
+		}
+		else {
+			amin = root.zmin;
+			amax = root.zmax;
+		}
+		float rdxyz = r.direction[xyz];
+		{
+			float t1 = (amin - r.origin[xyz]) / rdxyz;
+			float t2 = (amax - r.origin[xyz]) / rdxyz;
+			float ta = glm::min(t1, t2);
+			float tb = glm::max(t1, t2);
+			glm::vec3 n;
+			n[xyz] = t2 < t1 ? +1 : -1;
+			if (ta > 0 && ta > tmin) {
+				tmin = ta;
+			}
+			if (tb < tmax) {
+				tmax = tb;
+			}
+		}
+	}
+
+	if (tmax >= tmin && tmax > 0) {
+		if (root.triangleStart != -1 && root.triangleEnd != -1) {
+			for (int i = root.triangleStart; i <= root.triangleEnd; i++) {
+
+				int triangleIdx = sortTringles[i];
+				Triangle triangle = triangles[triangleIdx];
+				glm::vec3 p0 = triangle.vertices[0];
+				glm::vec3 p1 = triangle.vertices[1];
+				glm::vec3 p2 = triangle.vertices[2];
+				glm::vec3 n = glm::cross(p1 - p0, p2 - p0); // normal to check intersection inside a triangle, could opposite to triangle.normal 
+				float d = -glm::dot(triangle.normal, p1);
+
+				if (glm::abs(glm::dot(triangle.normal, rd)) < EPSILON) {
+					continue;
+				}
+				float t = (-d - glm::dot(triangle.normal, ro)) / (glm::dot(triangle.normal, rd));
+				if (t < 0) {
+					continue;
+				}
+				glm::vec3 p = ro + t * rd;
+
+				if (glm::dot(glm::cross(p1 - p0, p - p0), n) < 0) {
+					continue;
+				}
+				if (glm::dot(glm::cross(p2 - p1, p - p1), n) < 0) {
+					continue;
+				}
+				if (glm::dot(glm::cross(p0 - p2, p - p2), n) < 0) {
+					continue;
+				}
+				if (glm::dot(rd, triangle.normal) <= 0) {
+					outside = true;
+				}
+				else {
+					outside = false;
+				}
+
+				glm::vec3 objspaceIntersection = p;
+				intersectionPoint = multiplyMV(mesh.transform, glm::vec4(objspaceIntersection, 1.f));
+				normal = glm::normalize(multiplyMV(mesh.invTranspose, glm::vec4(triangle.normal, 0.f)));
+				t = glm::length(r.origin - intersectionPoint);
+				return;
+			}
+		}
+		else {
+			octreeIntersection(mesh, triangles, sortTringles, octreeVector, root.tlf, r, intersectionPoint, normal, outside, ro, rd, t);
+			octreeIntersection(mesh, triangles, sortTringles, octreeVector, root.tlb, r, intersectionPoint, normal, outside, ro, rd, t);
+			octreeIntersection(mesh, triangles, sortTringles, octreeVector, root.trf, r, intersectionPoint, normal, outside, ro, rd, t);
+			octreeIntersection(mesh, triangles, sortTringles, octreeVector, root.trb, r, intersectionPoint, normal, outside, ro, rd, t);
+			octreeIntersection(mesh, triangles, sortTringles, octreeVector, root.blf, r, intersectionPoint, normal, outside, ro, rd, t);
+			octreeIntersection(mesh, triangles, sortTringles, octreeVector, root.blb, r, intersectionPoint, normal, outside, ro, rd, t);
+			octreeIntersection(mesh, triangles, sortTringles, octreeVector, root.brf, r, intersectionPoint, normal, outside, ro, rd, t);
+			octreeIntersection(mesh, triangles, sortTringles, octreeVector, root.brb, r, intersectionPoint, normal, outside, ro, rd, t);
+		}
+	}
+	return;
+}
+
+__host__ __device__ float meshIntersectionTest_octree(
+	Geom mesh, Triangle *triangles, int *sortTringles,
+	OctreeNode_cuda *octreeVector, Ray r,
+	glm::vec3 &intersectionPoint, glm::vec3 &normal, bool &outside) {
+
+	glm::vec3 ro = multiplyMV(mesh.inverseTransform, glm::vec4(r.origin, 1.0f));
+	glm::vec3 rd = glm::normalize(multiplyMV(mesh.inverseTransform, glm::vec4(r.direction, 0.0f)));
+
+	float t = -1;
+	octreeIntersection(mesh, triangles, sortTringles, octreeVector, 0, r, intersectionPoint, normal, outside, ro, rd, t);
+	return t;
+}

src/main.cpp

@@ -42,7 +42,7 @@ int main(int argc, char** argv) {
 
     // Load scene file
     scene = new Scene(sceneFile);
-
+	
     // Set up camera stuff from loaded path tracer settings
     iteration = 0;
     renderState = &scene->state;
@@ -71,7 +71,6 @@ int main(int argc, char** argv) {
 
     // GLFW main loop
     mainLoop();
-
     return 0;
 }