Finish Refraction Ray & Anti-Aliasing (CIS565-Fall-2021#2)

* remove or move "using namespace std" into source

* finished part of core features

- Ideal diffuse surface model for shading
- Ray continuation with stream compaction
- Thrust sort for memory coalescence

* add `static` qualifier to device arrays

* increase hidpi scene resolution

* add `depth` to RNG generator in `shadeMaterial`

* add `const` qualifiers

* move extract material IDs to `computeIntersections`

* remove redundant `iterationComplete` variable

* finish core shading features

- perfect specular reflection
- ideal diffuse shading
- transition between half-reflection & half-diffuse

* increase lighting effects

* add Cornell variant -- multiple balls

* fix floating point issue in `getPointOnRay`

- avoid new ray colliding with the surface by increasing the floating point offset

* decrease light intensity

* move constants to static_config.h

* fix incorrect light rendering & reflections

* fix specular material color value in multiball scene

* add caching first depth ray trace intersections

* add sample img for core feature demo

* WIP on refraction

* add bloopers

* add glass material to cornell.txt

* WIP on debugging refraction code

* multiple improvements

- always normalize ray direction before calculation
- always use true_normal for reflection & diffusion

* add specular blue to cornell.txt

* add stochastic sampled anti-aliasing

* add bloopers in anti-aliasing

* fix ray refraction bug

- bug due to some magic floating point number
- thanks to https://piazza.com/class/kqefit06npb37y?cid=129

* remove one blooper

* add reflection img; add simple README

Author: shineyruan

README.md

@@ -18,3 +18,17 @@ Finished path tracing core features:
 - radix sort by material type
 - path continuation/termination by Thrust stream compaction 
 
+Finished Advanced Features:
+- Refraction with Fresnel effects using Schlick's approximation
+- Stochastic sampled anti-aliasing
+
+### Ray Refraction for Glass-like Materials
+|            Perfect Specular Reflection             |               Glass-like Refraction                |
+| :------------------------------------------------: | :------------------------------------------------: |
+| ![](img/cornell.2021-10-04_02-10-06z.5000samp.png) | ![](img/cornell.2021-10-04_01-57-31z.5000samp.png) |
+
+### Stochastic Sampled Anti-Aliasing
+|                 1x Anti-Aliasing (Feature OFF)                  |                        4x Anti-Aliasing                         |
+| :-------------------------------------------------------------: | :-------------------------------------------------------------: |
+| ![](img/cornell.2021-10-04_01-14-01z.5000samp-antialias-1x.png) | ![](img/cornell.2021-10-04_01-07-24z.5000samp-antialias-4x.png) |
+

img/cornell.2021-10-02_21-54-13z.5000samp.png

@@ -48,6 +48,26 @@ REFR        0
 REFRIOR     0
 EMITTANCE   0
 
+// Glass white
+MATERIAL 5
+RGB         1 1 1
+SPECEX      0
+SPECRGB     1 1 1
+REFL        0
+REFR        1
+REFRIOR     1.5
+EMITTANCE   0
+
+// Specular blue
+MATERIAL 6
+RGB         .35 .35 .85
+SPECEX      0
+SPECRGB     .35 .35 .85
+REFL        1
+REFR        0
+REFRIOR     0
+EMITTANCE   0
+
 // Camera
 CAMERA
 RES         800 800

img/cornell.2021-10-03_21-15-11z.5000samp.png

@@ -42,6 +42,23 @@ __host__ __device__ glm::vec3 calculateRandomDirectionInHemisphere(
          sin(around) * over * perpendicularDirection2;
 }
 
+/**
+ * @brief Use Schlick's approximation to calculate Fresnel factor based on given
+ * incident angle & index of refraction. Assumes ray always shoots from air to
+ * material (n_1 = 1).
+ *
+ * Reference: https://en.wikipedia.org/wiki/Schlick's_approximation
+ *
+ * @param cos_theta:    cos(theta), theta is incident angle
+ * @param idx_refract:  index refraction of material
+ * @return float
+ */
+__host__ __device__ float schlicks(const float cos_theta,
+                                   const float idx_refract) {
+  float r0 = powf((1.0f - idx_refract) / (1.0f + idx_refract), 2.0f);
+  return r0 + (1.0f - r0) * powf(1.0f - cos_theta, 5.0f);
+}
+
 /**
  * Scatter a ray with some probabilities according to the material properties.
  * For example, a diffuse surface scatters in a cosine-weighted hemisphere.
@@ -78,17 +95,41 @@ __host__ __device__ void scatterRay(PathSegment &pathSegment,
     thrust::uniform_real_distribution<float> u01(0, 1);
     const float random_sample = u01(rng);
 
-    if (random_sample < m.hasReflective) {
+    // check if the current ray is in material
+    glm::vec3 in_direction = glm::normalize(pathSegment.ray.direction);
+    bool is_inMaterial     = glm::dot(in_direction, normal) > 0.0f;
+
+    if (random_sample < m.hasRefractive) {
+      float idx_refraction_ratio =
+          (is_inMaterial) ? m.indexOfRefraction : (1.0f / m.indexOfRefraction);
+
+      // get the incident angle
+      float cos_angle = glm::abs(glm::dot(in_direction, normal));
+      // calculate Fresnel factor
+      float fresnel_factor = schlicks(cos_angle, m.indexOfRefraction);
+      // As refraction is determined by both reflection & transmission,
+      // probabilistically determine whether the next ray is reflective ray or
+      // transmission ray
+      glm::vec3 out_refracted_dir =
+          glm::refract(in_direction, normal, idx_refraction_ratio);
+      const float prob = u01(rng);
       pathSegment.ray.direction =
-          glm::reflect(pathSegment.ray.direction, normal);
+          (prob < fresnel_factor || glm::length(out_refracted_dir) < EPS)
+              ? glm::reflect(in_direction, normal)
+              : out_refracted_dir;
+      pathSegment.color *= m.specular.color;
+    } else if (random_sample < m.hasReflective) {
+      pathSegment.ray.direction = glm::reflect(in_direction, normal);
       pathSegment.color *= m.specular.color;
     } else {
       pathSegment.ray.direction =
           calculateRandomDirectionInHemisphere(normal, rng);
     }
 
     --pathSegment.remainingBounces;
-    pathSegment.ray.origin = intersect + EPS * normal;
+    // 0.001f is a really WIERD magic number here. We cannot unify this offset
+    // with EPS since that would BREAK the refraction renderer.
+    pathSegment.ray.origin = intersect + 0.001f * pathSegment.ray.direction;
     pathSegment.color *= m.color;
     pathSegment.color =
         glm::clamp(pathSegment.color, glm::vec3(0.f), glm::vec3(1.0f));

img/cornell.2021-10-04_01-07-24z.5000samp-antialias-4x.png

@@ -81,10 +81,14 @@ static PathSegment *dev_paths                   = NULL;
 static ShadeableIntersection *dev_intersections = NULL;
 // TODO: static variables for device memory, any extra info you need, etc
 // ...
-static int *dev_materialIDs       = NULL;
-static int *dev_materialIDBuffers = NULL;
+static int *dev_materialIDs        = NULL;
+static int *dev_materialIDBuffers  = NULL;
+static glm::vec3 *dev_image_buffer = NULL;
+
+// first-bounce intersection caching
 #ifdef CACHE_INTERSECTIONS
 static ShadeableIntersection *dev_intersections_cache = NULL;
+static int *dev_materialIDs_cache                     = NULL;
 #endif
 
 void pathtraceInit(Scene *scene) {
@@ -95,7 +99,7 @@ void pathtraceInit(Scene *scene) {
   cudaMalloc(&dev_image, pixelcount * sizeof(glm::vec3));
   cudaMemset(dev_image, 0, pixelcount * sizeof(glm::vec3));
 
-  cudaMalloc(&dev_paths, pixelcount * sizeof(PathSegment));
+  cudaMalloc(&dev_paths, ANTIALIAS_FACTOR * pixelcount * sizeof(PathSegment));
 
   cudaMalloc(&dev_geoms, scene->geoms.size() * sizeof(Geom));
   cudaMemcpy(dev_geoms, scene->geoms.data(), scene->geoms.size() * sizeof(Geom),
@@ -106,15 +110,22 @@ void pathtraceInit(Scene *scene) {
              scene->materials.size() * sizeof(Material),
              cudaMemcpyHostToDevice);
 
-  cudaMalloc(&dev_intersections, pixelcount * sizeof(ShadeableIntersection));
-  cudaMemset(dev_intersections, 0, pixelcount * sizeof(ShadeableIntersection));
+  cudaMalloc(&dev_intersections,
+             ANTIALIAS_FACTOR * pixelcount * sizeof(ShadeableIntersection));
+  cudaMemset(dev_intersections, 0,
+             ANTIALIAS_FACTOR * pixelcount * sizeof(ShadeableIntersection));
 
   // TODO: initialize any extra device memeory you need
-  cudaMalloc((void **)&dev_materialIDs, pixelcount * sizeof(int));
-  cudaMalloc((void **)&dev_materialIDBuffers, pixelcount * sizeof(int));
+  cudaMalloc((void **)&dev_materialIDs,
+             ANTIALIAS_FACTOR * pixelcount * sizeof(int));
+  cudaMalloc((void **)&dev_materialIDBuffers,
+             ANTIALIAS_FACTOR * pixelcount * sizeof(int));
+  cudaMalloc(&dev_image_buffer, pixelcount * sizeof(glm::vec3));
 #ifdef CACHE_INTERSECTIONS
   cudaMalloc((void **)&dev_intersections_cache,
-             pixelcount * sizeof(ShadeableIntersection));
+             ANTIALIAS_FACTOR * pixelcount * sizeof(ShadeableIntersection));
+  cudaMalloc((void **)&dev_materialIDs_cache,
+             ANTIALIAS_FACTOR * pixelcount * sizeof(int));
 #endif
 
   checkCUDAError("pathtraceInit");
@@ -129,8 +140,10 @@ void pathtraceFree() {
   // TODO: clean up any extra device memory you created
   cudaFree(dev_materialIDs);
   cudaFree(dev_materialIDBuffers);
+  cudaFree(dev_image_buffer);
 #ifdef CACHE_INTERSECTIONS
   cudaFree(dev_intersections_cache);
+  cudaFree(dev_materialIDs_cache);
 #endif
 
   checkCUDAError("pathtraceFree");
@@ -150,22 +163,40 @@ __global__ void generateRayFromCamera(Camera cam, int iter, int traceDepth,
   int y = (blockIdx.y * blockDim.y) + threadIdx.y;
 
   if (x < cam.resolution.x && y < cam.resolution.y) {
-    int index            = x + (y * cam.resolution.x);
-    PathSegment &segment = pathSegments[index];
-
-    segment.ray.origin = cam.position;
-    segment.color      = glm::vec3(1.0f, 1.0f, 1.0f);
+    int index = x + (y * cam.resolution.x);
 
-    // TODO: implement antialiasing by jittering the ray
+    // primary ray per pixel
+    PathSegment &segment     = pathSegments[index];
+    segment.ray.origin       = cam.position;
+    segment.color            = glm::vec3(1.0f, 1.0f, 1.0f);
+    segment.pixelIndex       = index;
+    segment.remainingBounces = traceDepth;
     segment.ray.direction =
         glm::normalize(cam.view -
                        cam.right * cam.pixelLength.x *
                            ((float)x - (float)cam.resolution.x * 0.5f) -
                        cam.up * cam.pixelLength.y *
                            ((float)y - (float)cam.resolution.y * 0.5f));
 
-    segment.pixelIndex       = index;
-    segment.remainingBounces = traceDepth;
+    // implement antialiasing by jittering the ray
+    // sub-sampled extra rays per pixel
+    int pixelcount = cam.resolution.x * cam.resolution.y;
+    for (int i = 1; i < ANTIALIAS_FACTOR; ++i) {
+      PathSegment &extra_segment     = pathSegments[i * pixelcount + index];
+      extra_segment.ray.origin       = cam.position;
+      extra_segment.color            = glm::vec3(1.0f, 1.0f, 1.0f);
+      extra_segment.pixelIndex       = index;
+      extra_segment.remainingBounces = traceDepth;
+      thrust::default_random_engine rng =
+          makeSeededRandomEngine(iter, index, i);
+      thrust::uniform_real_distribution<float> u01(0, 1);
+      extra_segment.ray.direction = glm::normalize(
+          cam.view -
+          cam.right * cam.pixelLength.x *
+              ((float)x + u01(rng) - (float)cam.resolution.x * 0.5f) -
+          cam.up * cam.pixelLength.y *
+              ((float)y + u01(rng) - (float)cam.resolution.y * 0.5f));
+    }
   }
 }
 
@@ -263,14 +294,25 @@ __global__ void shadeMaterial(
   }
 }
 
-// Add the current iteration's output to the overall image
-__global__ void finalGather(int nPaths, glm::vec3 *image,
-                            PathSegment *iterationPaths) {
+// Add the current iteration's output to the image buffer
+__global__ void finalGather(int nPaths, glm::vec3 *img_buffer,
+                            const PathSegment *iterationPaths) {
   int index = (blockIdx.x * blockDim.x) + threadIdx.x;
-
   if (index < nPaths) {
     PathSegment iterationPath = iterationPaths[index];
-    image[iterationPath.pixelIndex] += iterationPath.color;
+    atomicAdd(&img_buffer[iterationPath.pixelIndex][0], iterationPath.color[0]);
+    atomicAdd(&img_buffer[iterationPath.pixelIndex][1], iterationPath.color[1]);
+    atomicAdd(&img_buffer[iterationPath.pixelIndex][2], iterationPath.color[2]);
+  }
+}
+
+// Average the accumulative subpixel values in image buffer & add it to final
+// image
+__global__ void addToImage(int pixelcount, glm::vec3 *image,
+                           const glm::vec3 *img_buffer) {
+  int index = (blockIdx.x * blockDim.x) + threadIdx.x;
+  if (index < pixelcount) {
+    image[index] += (img_buffer[index] / (1.0f * ANTIALIAS_FACTOR));
   }
 }
 
@@ -326,15 +368,14 @@ void pathtrace(uchar4 *pbo, int frame, int iter) {
   checkCUDAError("generate camera ray");
 
   int depth            = 0;
-  int num_active_paths = pixelcount;
+  int num_active_paths = ANTIALIAS_FACTOR * pixelcount;
 
   // --- PathSegment Tracing Stage ---
   // Shoot ray into scene, bounce between objects, push shading chunks
-
   while (num_active_paths > 0) {
     // clean shading chunks
     cudaMemset(dev_intersections, 0,
-               pixelcount * sizeof(ShadeableIntersection));
+               ANTIALIAS_FACTOR * pixelcount * sizeof(ShadeableIntersection));
 
     // --- Tracing Stage ---
     dim3 numblocksPathSegmentTracing =
@@ -345,6 +386,18 @@ void pathtrace(uchar4 *pbo, int frame, int iter) {
       cudaMemcpy(dev_intersections, dev_intersections_cache,
                  pixelcount * sizeof(ShadeableIntersection),
                  cudaMemcpyDeviceToDevice);
+      cudaMemcpy(dev_materialIDs, dev_materialIDs_cache,
+                 pixelcount * sizeof(int), cudaMemcpyDeviceToDevice);
+      if (num_active_paths - pixelcount > 0) {
+        dim3 numBlocksAntialiasTracing =
+            (num_active_paths - pixelcount + blockSize1d - 1) / blockSize1d;
+        computeIntersections<<<numBlocksAntialiasTracing, blockSize1d>>>(
+            depth, num_active_paths - pixelcount, dev_paths + pixelcount,
+            dev_geoms, hst_scene->geoms.size(), dev_intersections + pixelcount,
+            dev_materialIDs + pixelcount);
+        checkCUDAError("anti-alias extra rays trace one bounce");
+        cudaDeviceSynchronize();
+      }
     } else {
       computeIntersections<<<numblocksPathSegmentTracing, blockSize1d>>>(
           depth, num_active_paths, dev_paths, dev_geoms,
@@ -356,6 +409,8 @@ void pathtrace(uchar4 *pbo, int frame, int iter) {
         cudaMemcpy(dev_intersections_cache, dev_intersections,
                    pixelcount * sizeof(ShadeableIntersection),
                    cudaMemcpyDeviceToDevice);
+        cudaMemcpy(dev_materialIDs_cache, dev_materialIDs,
+                   pixelcount * sizeof(int), cudaMemcpyDeviceToDevice);
       }
     }
 #else
@@ -401,9 +456,14 @@ void pathtrace(uchar4 *pbo, int frame, int iter) {
   }
 
   // Assemble this iteration and apply it to the image
+  cudaMemset(dev_image_buffer, 0, pixelcount * sizeof(glm::vec3));
+  dim3 numBlocksSubPixels =
+      (ANTIALIAS_FACTOR * pixelcount + blockSize1d - 1) / blockSize1d;
+  finalGather<<<numBlocksSubPixels, blockSize1d>>>(
+      ANTIALIAS_FACTOR * pixelcount, dev_image_buffer, dev_paths);
   dim3 numBlocksPixels = (pixelcount + blockSize1d - 1) / blockSize1d;
-  finalGather<<<numBlocksPixels, blockSize1d>>>(pixelcount, dev_image,
-                                                dev_paths);
+  addToImage<<<numBlocksPixels, blockSize1d>>>(pixelcount, dev_image,
+                                               dev_image_buffer);
 
   // Send results to OpenGL buffer for rendering
   sendImageToPBO<<<blocksPerGrid2d, blockSize2d>>>(pbo, cam.resolution, iter,