Project 4: (Charles) Zixin Zhang

README.md

@@ -3,11 +3,93 @@ CUDA Denoiser For CUDA Path Tracer
 
 **University of Pennsylvania, CIS 565: GPU Programming and Architecture, Project 4**
 
-* (TODO) YOUR NAME HERE
-* Tested on: (TODO) Windows 22, i7-2222 @ 2.22GHz 22GB, GTX 222 222MB (Moore 2222 Lab)
+* (Charles) Zixin Zhang
+* Tested on: Windows 11, 11th Core i7, 3060 Laptop GPU
 
-### (TODO: Your README)
+![pathTracer](images/pathTracer.gif)
 
-*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.
+
+
+# Results
+
+Figure 1: Zoomed-in "So many balls" Scene: 
+
+Denoised: 
+
+![](Denoised.png)
+
+Noised: 
+
+![](Noised.png)
+
+---
+
+Figure 2: "So many balls" Scene:
+
+Denoised: 
+
+![](Denoised2.png)
+
+Noised: 
+
+![](Noised2.png)
+
+
+---
+
+Figure 3: A Simple Cornell Box:  
+
+![](comp.png)
+
+
+This denoiser is achieved by implementing the paper "Edge-Avoiding A-Trous Wavelet Transform for fast Global Illumination Filtering," by Dammertz, Sewtz, Hanika, and Lensch. 
+
+# Performance Analysis
+
+ In my implementation, denoising is performed during the last iteration. We perform 10 iterations on the ray traced images and the 11th iteration is the denoising step. 
+
+## Denoising Performance
+
+|             | How much time does denoising take in ms? |
+| ----------- | ----------- |
+| Figure 1      | 52       |
+| Figure 2   | 43        |
+| Figure 3   | 44        |
+
+Since we apply denoising once, this technique is very efficient at eliminating noise in the scene. 
+
+
+## Comparison with Pure Ray Traced Images
+
+|             | How many iterations do ray traced images need to achieve a smooth result? |
+| ----------- | ----------- |
+| Figure 1      | ~1000     | 
+| Figure 2   | ~5000        | 
+| Figure 3   | ~5000        | 
+
+As shown in the table, if we took more samples per pixel to eliminate the noise, it would take more iterations and time to achieve a similar result. 
+
+## Resolution
+
+Denoising took a significant amount of extra time when the resolution went from 1080p to 4k. It is expected since our run time depends linearly on the pixel count in the scene. 
+
+|             | 720p | 1080p | 4K |
+| ----------- | ----------- | ----------- | ----------- |
+| Figure 1      | 10ms  |             ~52ms         |  120ms     |
+| Figure 2   |   ~16ms   |        ~43ms        |~215ms|
+| Figure 3   |  ~17ms    |        ~44ms        |~215ms|
+
+
+# Changes 
+
+In order to programmatically generate a lot of balls (121 balls in figure 1 and 2), scene files are only used to specify the camera specifications. All geometric shapes and materials are specified at run time in the actual code. 
+
+Fuzziness is added to the material properties such as we can have [fuzzy reflection](https://raytracing.github.io/books/RayTracingInOneWeekend.html#metal/fuzzyreflection). 
+
+Figure 1 and 2 is genreated by `scenes/manyBalls.txt`. 
+
+# References
+
+- [_Ray Tracing in One Weekend Series by Peter Shirley_](https://raytracing.github.io/books/RayTracingInOneWeekend.html)
+- [_Edge-Avoiding A-Trous Wavelet Transform for fast Global Illumination Filtering by Dammertz, Sewtz, Hanika, and Lensch_](https://jo.dreggn.org/home/2010_atrous.pdf)
 

scenes/cornell_ceiling_light.txt

@@ -50,9 +50,9 @@ EMITTANCE   0
 
 // Camera
 CAMERA
-RES         800 800
+RES         1280  720
 FOVY        45
-ITERATIONS  10
+ITERATIONS  1000
 DEPTH       8
 FILE        cornell
 EYE         0.0 5 10.5

scenes/manyBalls.txt

@@ -0,0 +1,10 @@
+// Camera
+CAMERA
+RES         1280 720
+FOVY        45
+ITERATIONS  10
+DEPTH       8
+FILE        sphere
+EYE         13 -2 3
+LOOKAT      0 0 0
+UP          0 1 0

src/interactions.h

@@ -40,23 +40,122 @@ glm::vec3 calculateRandomDirectionInHemisphere(
         + sin(around) * over * perpendicularDirection2;
 }
 
+__host__ __device__
+glm::vec3 calculateRandomDirectionInSphere(thrust::default_random_engine& rng)
+{
+    thrust::uniform_real_distribution<float> u01(-1.f, 1.0f);
+    while (true)
+    {
+        glm::vec3 p = glm::vec3(u01(rng), u01(rng), u01(rng));
+        float length = glm::length(p);
+        if (length >= 1.0f) { continue; }
+        return p;
+    }
+}
+
+__host__ __device__
+float reflectance(float cosine, float refIdx)
+{
+    // Use Schlick's approximation for reflectance.
+    float r0 = (1.0f - refIdx) / (1.0f + refIdx);
+    r0 *= r0;
+    return r0 + (1.0f - r0) * glm::pow(1.0f - cosine, 5.0f);
+}
+
 /**
- * Simple ray scattering with diffuse and perfect specular support.
+ * Scatter a ray with some probabilities according to the material properties.
+ * For example, a diffuse surface scatters in a cosine-weighted hemisphere.
+ * A perfect specular surface scatters in the reflected ray direction.
+ * In order to apply multiple effects to one surface, probabilistically choose
+ * between them.
+ *
+ * The visual effect you want is to straight-up add the diffuse and specular
+ * components. You can do this in a few ways. This logic also applies to
+ * combining other types of materias (such as refractive).
+ *
+ * - Always take an even (50/50) split between a each effect (a diffuse bounce
+ *   and a specular bounce), but divide the resulting color of either branch
+ *   by its probability (0.5), to counteract the chance (0.5) of the branch
+ *   being taken.
+ *   - This way is inefficient, but serves as a good starting point - it
+ *     converges slowly, especially for pure-diffuse or pure-specular.
+ * - Pick the split based on the intensity of each material color, and divide
+ *   branch result by that branch's probability (whatever probability you use).
+ *
+ * This method applies its changes to the Ray parameter `ray` in place.
+ * It also modifies the color `color` of the ray in place.
+ *
+ * You may need to change the parameter list for your purposes!
  */
+
 __host__ __device__
 void scatterRay(
-		PathSegment & pathSegment,
-        glm::vec3 intersect,
-        glm::vec3 normal,
-        const Material &m,
-        thrust::default_random_engine &rng) {
-    glm::vec3 newDirection;
-    if (m.hasReflective) {
-        newDirection = glm::reflect(pathSegment.ray.direction, normal);
-    } else {
-        newDirection = calculateRandomDirectionInHemisphere(normal, rng);
+    PathSegment& pathSegment,
+    glm::vec3 intersect,
+    glm::vec3 normal,
+    bool outside,
+    const Material& m,
+    thrust::default_random_engine& rng,
+    const ShadeableIntersection& sInter) {
+
+    thrust::uniform_real_distribution<float> u01(0, 1);
+
+    if (m.hasReflective)
+    {
+        pathSegment.ray.origin = intersect + EPSILON * normal;
+        glm::vec3 reflectedDir = glm::reflect(glm::normalize(pathSegment.ray.direction),
+            normal);
+        // add fuzziness
+        reflectedDir += m.fuzziness * calculateRandomDirectionInSphere(rng);
+
+        if (glm::dot(reflectedDir, normal) > 0.0f)
+        {
+            pathSegment.color *= m.specular.color;
+            pathSegment.ray.direction = glm::normalize(reflectedDir);
+            pathSegment.remainingBounces--;
+        }
+        // for big sphere or grazing rays, we may scatter below the 
+        // surface. In this case, terminate this segment. 
+        else {
+            // NOTE: this line is necessary to prevent the white boundary
+            // if we terminate the ray path because reflected ray goes below
+            // the surface, this path's contribution should be set to black (0.f)
+            pathSegment.color *= glm::vec3(0.f);
+            pathSegment.remainingBounces = 0;
+        }
+    }
+    else if (m.hasRefractive)
+    {
+        float refractionRatio = outside ? (1.f / m.indexOfRefraction) : m.indexOfRefraction;
+        // refractive rays always shoots towards negative normal direction: that's why we use subtraction
+        // since intersect falls slightly short to the object it's hitting, 
+        // we need a bigger EPSILON so that reflective rays are shoot 
+        // from a point that is not occluded by the surface
+        glm::vec3 unitRayDir = glm::normalize(pathSegment.ray.direction);
+        float cosTheta = fmin(glm::dot(-unitRayDir, normal), 1.0f);
+        float sinTheta = sqrt(1.0f - cosTheta * cosTheta);
+        bool cannotReflect = refractionRatio * sinTheta > 1.0f;
+        glm::vec3 newRayDir;
+        if (cannotReflect || reflectance(cosTheta, refractionRatio) > u01(rng))
+        {
+            pathSegment.ray.origin = intersect + EPSILON * normal;
+            newRayDir = glm::reflect(unitRayDir, normal);
+        }
+        else {
+            pathSegment.ray.origin = intersect - EPSILON * 100.0f * normal;
+            newRayDir = glm::refract(unitRayDir, normal, refractionRatio);
+        }
+        pathSegment.color *= m.color;
+        pathSegment.ray.direction = glm::normalize(newRayDir);
+        pathSegment.remainingBounces--;
+    }
+    else {
+        pathSegment.ray.origin = intersect + EPSILON * normal;
+        glm::vec3 diffuseDir = glm::normalize(calculateRandomDirectionInHemisphere(normal, rng));
+        pathSegment.color *= m.color;
+        pathSegment.ray.direction = diffuseDir;
+        // diffuse always scatter
+        pathSegment.remainingBounces--;
     }
 
-    pathSegment.ray.direction = newDirection;
-    pathSegment.ray.origin = intersect + (newDirection * 0.0001f);
 }

src/intersections.h

@@ -19,11 +19,13 @@ __host__ __device__ inline unsigned int utilhash(unsigned int a) {
     return a;
 }
 
+// CHECKITOUT
 /**
  * Compute a point at parameter value `t` on ray `r`.
  * Falls slightly short so that it doesn't intersect the object it's hitting.
  */
 __host__ __device__ glm::vec3 getPointOnRay(Ray r, float t) {
+    // return r.origin + (t - .0001f) * glm::normalize(r.direction);
     return r.origin + (t - .0001f) * glm::normalize(r.direction);
 }
 
@@ -34,6 +36,7 @@ __host__ __device__ glm::vec3 multiplyMV(glm::mat4 m, glm::vec4 v) {
     return glm::vec3(m * v);
 }
 
+// CHECKITOUT
 /**
  * Test intersection between a ray and a transformed cube. Untransformed,
  * the cube ranges from -0.5 to 0.5 in each axis and is centered at the origin.
@@ -44,9 +47,9 @@ __host__ __device__ glm::vec3 multiplyMV(glm::mat4 m, glm::vec4 v) {
  * @return                   Ray parameter `t` value. -1 if no intersection.
  */
 __host__ __device__ float boxIntersectionTest(Geom box, Ray r,
-        glm::vec3 &intersectionPoint, glm::vec3 &normal, bool &outside) {
+    glm::vec3& intersectionPoint, glm::vec3& normal, bool& outside) {
     Ray q;
-    q.origin    =                multiplyMV(box.inverseTransform, glm::vec4(r.origin   , 1.0f));
+    q.origin = multiplyMV(box.inverseTransform, glm::vec4(r.origin, 1.0f));
     q.direction = glm::normalize(multiplyMV(box.inverseTransform, glm::vec4(r.direction, 0.0f)));
 
     float tmin = -1e38f;
@@ -87,6 +90,7 @@ __host__ __device__ float boxIntersectionTest(Geom box, Ray r,
     return -1;
 }
 
+// CHECKITOUT
 /**
  * Test intersection between a ray and a transformed sphere. Untransformed,
  * the sphere always has radius 0.5 and is centered at the origin.
@@ -97,8 +101,8 @@ __host__ __device__ float boxIntersectionTest(Geom box, Ray r,
  * @return                   Ray parameter `t` value. -1 if no intersection.
  */
 __host__ __device__ float sphereIntersectionTest(Geom sphere, Ray r,
-        glm::vec3 &intersectionPoint, glm::vec3 &normal, bool &outside) {
-    float radius = .5;
+    glm::vec3& intersectionPoint, glm::vec3& normal, bool& outside) {
+    float radius = 0.5f;
 
     glm::vec3 ro = multiplyMV(sphere.inverseTransform, glm::vec4(r.origin, 1.0f));
     glm::vec3 rd = glm::normalize(multiplyMV(sphere.inverseTransform, glm::vec4(r.direction, 0.0f)));
@@ -121,21 +125,34 @@ __host__ __device__ float sphereIntersectionTest(Geom sphere, Ray r,
     float t = 0;
     if (t1 < 0 && t2 < 0) {
         return -1;
-    } else if (t1 > 0 && t2 > 0) {
+    }
+    else if (t1 > 0 && t2 > 0) {
         t = min(t1, t2);
-        outside = true;
-    } else {
+        // outside = true;
+    }
+    else {
         t = max(t1, t2);
-        outside = false;
+        // outside = false;
     }
 
     glm::vec3 objspaceIntersection = getPointOnRay(rt, t);
 
     intersectionPoint = multiplyMV(sphere.transform, glm::vec4(objspaceIntersection, 1.f));
     normal = glm::normalize(multiplyMV(sphere.invTranspose, glm::vec4(objspaceIntersection, 0.f)));
+
+    if (sphere.scale[0] < 1.0f)
+    {
+        normal = -normal;
+    }
+
+    outside = glm::dot(r.direction, normal) < 0;
+    // make it so that normals always point against the incident ray
+    // If the ray is outside the geometry, the normal will point outward, 
+    // but if the ray is inside the geometry, the normal will point inward.
     if (!outside) {
         normal = -normal;
     }
 
+
     return glm::length(r.origin - intersectionPoint);
 }