Project 4: Helena Zhang

README.md

@@ -3,11 +3,95 @@ 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)
+* Helena Zhang
+* Tested on: Windows 11, i7-10750 @ 2.6GHz 16GB, Geforce RTX 2060 6GB
 
-### (TODO: Your README)
+### Results
+Pathtraced Cornell Box at 20 iterations, without and with Atrous Denoising
 
-*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.
+<p float="left">
+  <img src="img/20IterNaive.png" width="40%" />
+  <img src="img/20IterSmooth.png" width="40%" />
+</p>
 
+### Analysis
+
+**Performance:**
+The denoising functionality adds two additional kernels to the base pathtracing code: a gBuffer storage kernel that records a path's first bounce's position and normals, and the atrous denoising convolution kernel (function not the matrix) that will adjust the color results with its neighbors. Since the two kernels will run for every path used to render the image, the number of calls and the runtime increases as the number of paths, or resolution increases. 
+<p float="left">
+  <img src="img/threadsRuntime.png" width="40%" />
+  <img src="img/threadsRuntimeLog.png" width="40%" />
+</p>
+
+The two kernel's runtimes are compared above for 100 by 100, 200 by 200, 400 by 400, 800 by 800, and 1600 by 1600 pixels. The runtimes for both kernels increased exponentially when compared to the length of the images, but linearly when compared to the number of pixels, or threads. 
+
+Another variable used in the atrous kernel is the size of the convolution filter. The filter takes 25 pixels around the target pixel, where the farther pixel taken is +/- (2s + 1, 2s + 1) where s is gap between each adjacent pixels used to calculate the weighted average color. 
+
+<img src="img/filterSizeRuntime.png" width="80%" />
+
+The figure above compares the average runtime for each execution of the atrous kernel where each pixel used for the kernel is 1, 2, 5, and 10 pixels apart. There is a slight decrease in average runtime per kernel as the filter size increases. This is because for larger kernel sizes, the atrous kernel will find pixels that are farther away, and these neighbors may not exist for pixels near the edges of the image. Thus, for these pixels, its atrous calls will have fewer memory access, hence the shorter runtime. However, the overall runtime will still be higher for larger filter sizes, since it will have more atrous kernel calls, but each atrous call takes on average a short time. 
+
+**Visuals:**
+With the correct weights and filter size, the denoiser can smooth images with just two iteration of pathtracing:
+<p float="left">
+  <img src="img/lowestIterNaive.png" width="40%" />
+  <img src="img/lowestIterSmooth.png" width="40%" />
+</p>
+
+As filter sizes increase, more pixels with a larger range will be taken into consideration when smoothing. 
+<p float="left">
+ <img src="img/default1.png" width="51%" />
+
+  Filter size: 5*5 (No space between pixels)
+
+
+  <img src="img/default2.png" width="51%" />
+
+   Filter size: 9*9 (1 space between pixels)
+
+
+  <img src="img/default3.png" width="51%" />
+
+   Filter size: 13*13 (2 spaces between pixels)
+
+
+  <img src="img/default4.png" width="51%" />
+
+ Filter size: 17*17 (3 spaces between pixels)
+
+
+</p>
+
+The images appear more smooth as filter size increases, as the number of atrous kernel calls and kernels taken into consideration increases. However, there are diminishing returns for larger filter sizes. 
+
+<p float="left">
+  <img src="img/default5.png" width="51%" />
+
+   Filter size: 21*21 (4 spaces between pixels)
+
+
+  <img src="img/default7.png" width="51%" />
+
+   Filter size: 29*29 (6 spaces between pixels)
+
+
+  <img src="img/default10.png" width="51%" />
+
+ Filter size: 41*41 (9 spaces between pixels)
+
+</p>
+
+Finally, although the atrous denoising technique can smooth images without much additional runtime, there are still some caveats. 
+
+For certain materials and scenes, the technique doesn't perform as well. For instance:
+
+<img src="img/20IterSmooth.png" width="80%" />
+
+The render above appears mostly smooth, but there are reflected edges of the reflecting ball are blurred. 
+
+<p float="left">
+  <img src="img/20IterCornell.png" width="40%" />
+  <img src="img/25IterCornell.png" width="40%" />
+</p>
+
+In a Cornell box with less light (hence more noise), there are noisy specs (left), even though the render uses the same setting as the brighter Cornell box above it. With some adjustments to the filter size and weights, most of the specs are gone (right), but there is still noise around the light source, and the perfect reflective sphere appears blurry. 

scenes/cornell_ceiling_light.txt

@@ -50,7 +50,7 @@ EMITTANCE   0
 
 // Camera
 CAMERA
-RES         800 800
+RES         1200 1200
 FOVY        45
 ITERATIONS  10
 DEPTH       8

src/intersections.h

@@ -24,7 +24,7 @@ __host__ __device__ inline unsigned int utilhash(unsigned int a) {
  * 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 * glm::normalize(r.direction);
 }
 
 /**

src/main.cpp

@@ -23,11 +23,11 @@ int ui_iterations = 0;
 int startupIterations = 0;
 int lastLoopIterations = 0;
 bool ui_showGbuffer = false;
-bool ui_denoise = false;
-int ui_filterSize = 80;
-float ui_colorWeight = 0.45f;
-float ui_normalWeight = 0.35f;
-float ui_positionWeight = 0.2f;
+bool ui_denoise = true;
+int ui_filterSize = 5;
+float ui_colorWeight = 0.50f;
+float ui_normalWeight = 0.50f;
+float ui_positionWeight = 0.50f;
 bool ui_saveAndExit = false;
 
 static bool camchanged = true;
@@ -151,6 +151,7 @@ void runCuda() {
 
     if (iteration == 0) {
         pathtraceFree();
+        updateData(ui_denoise, ui_filterSize, ui_positionWeight, ui_normalWeight, ui_colorWeight);
         pathtraceInit(scene);
     }
 
@@ -164,11 +165,12 @@ void runCuda() {
         int frame = 0;
         pathtrace(frame, iteration);
     }
-
+    
     if (ui_showGbuffer) {
       showGBuffer(pbo_dptr);
-    } else {
-      showImage(pbo_dptr, iteration);
+    } 
+    else {
+        showImage(pbo_dptr, iteration);
     }
 
     // unmap buffer object

src/pathtrace.cu

@@ -73,24 +73,58 @@ __global__ void gbufferToPBO(uchar4* pbo, glm::ivec2 resolution, GBufferPixel* g
 
     if (x < resolution.x && y < resolution.y) {
         int index = x + (y * resolution.x);
-        float timeToIntersect = gBuffer[index].t * 256.0;
+
+        glm::vec3 positions = gBuffer[index].pos * 25.5f;
+        glm::vec3 normals = 127.0f * (gBuffer[index].nor + glm::vec3(1.0f));
 
         pbo[index].w = 0;
-        pbo[index].x = timeToIntersect;
-        pbo[index].y = timeToIntersect;
-        pbo[index].z = timeToIntersect;
+        pbo[index].x = normals.x;
+        pbo[index].y = normals.y;
+        pbo[index].z = normals.z;
     }
 }
 
+
+
 static Scene * hst_scene = NULL;
 static glm::vec3 * dev_image = NULL;
+static glm::vec3* dev_image_clean = NULL;
 static Geom * dev_geoms = NULL;
 static Material * dev_materials = NULL;
 static PathSegment * dev_paths = NULL;
 static ShadeableIntersection * dev_intersections = NULL;
 static GBufferPixel* dev_gBuffer = NULL;
 // TODO: static variables for device memory, any extra info you need, etc
-// ...
+
+__constant__ float dev_kernel[25];
+__constant__ int dev_offset[50];
+static bool denoising;
+static int fs;
+static float pw;
+static float nw;
+static float cw;
+
+const float kernel_host[25] =
+{ 0.0039, 0.0156, 0.0234, 0.0156, 0.0039,
+0.0156, 0.0625, 0.0936, 0.0625, 0.0156,
+0.0234, 0.0936, 0.1406, 0.0936, 0.0234,
+0.0156, 0.0625, 0.0936, 0.0625, 0.0156,
+0.0039, 0.0156, 0.0234, 0.0156, 0.0039 };
+const glm::ivec2 offset_host[25] =
+{ glm::ivec2(-2, -2), glm::ivec2(-1, -2), glm::ivec2(0, -2), glm::ivec2(1, -2), glm::ivec2(2, -2),
+glm::ivec2(-2, -1), glm::ivec2(-1, -1), glm::ivec2(0, -1), glm::ivec2(1, -1), glm::ivec2(2, -1),
+glm::ivec2(-2, 0), glm::ivec2(-1, 0), glm::ivec2(0, 0), glm::ivec2(1, 0), glm::ivec2(2, 0),
+glm::ivec2(-2, 1), glm::ivec2(-1, 1), glm::ivec2(0, 1), glm::ivec2(1, 1), glm::ivec2(2, 1),
+glm::ivec2(-2, 2), glm::ivec2(-1, 2), glm::ivec2(0, 2), glm::ivec2(1, 2), glm::ivec2(2, 2) };
+
+void updateData(bool denoise, int filterSize, float posWeight, float norWeight, float colWeight) {
+    denoising = denoise;
+    fs = filterSize;
+    pw = posWeight;
+    nw = norWeight;
+    cw = colWeight;
+}
+
 
 void pathtraceInit(Scene *scene) {
     hst_scene = scene;
@@ -100,6 +134,9 @@ void pathtraceInit(Scene *scene) {
     cudaMalloc(&dev_image, pixelcount * sizeof(glm::vec3));
     cudaMemset(dev_image, 0, pixelcount * sizeof(glm::vec3));
 
+    cudaMalloc(&dev_image_clean, pixelcount * sizeof(glm::vec3));
+    cudaMemset(dev_image_clean, 0, pixelcount * sizeof(glm::vec3));
+
   	cudaMalloc(&dev_paths, pixelcount * sizeof(PathSegment));
 
   	cudaMalloc(&dev_geoms, scene->geoms.size() * sizeof(Geom));
@@ -115,11 +152,16 @@ void pathtraceInit(Scene *scene) {
 
     // TODO: initialize any extra device memeory you need
 
+    cudaMemcpyToSymbol(dev_kernel, kernel_host, 25 * sizeof(float));
+    cudaMemcpyToSymbol(dev_offset, offset_host, 50 * sizeof(int));
+
+
     checkCUDAError("pathtraceInit");
 }
 
 void pathtraceFree() {
     cudaFree(dev_image);  // no-op if dev_image is null
+    cudaFree(dev_image_clean);
   	cudaFree(dev_paths);
   	cudaFree(dev_geoms);
   	cudaFree(dev_materials);
@@ -130,6 +172,69 @@ void pathtraceFree() {
     checkCUDAError("pathtraceFree");
 }
 
+
+// A-trous
+__global__ void atrous(
+    const int filterSize,
+    const float pos_weight,
+    const float nor_weight,
+    float col_weight,
+    glm::ivec2 resolution,
+    GBufferPixel* gBuffer,
+    glm::vec3* src,
+    glm::vec3* dst,
+    float iter_i) {
+
+    int x = (blockIdx.x * blockDim.x) + threadIdx.x;
+    int y = (blockIdx.y * blockDim.y) + threadIdx.y;
+
+    // pixel (x, y)
+    if (x < resolution.x && y < resolution.y) {
+        int index = x + (y * resolution.x);
+
+        glm::vec3 sum(0.0f);
+        glm::vec3 pval = gBuffer[index].pos;
+        glm::vec3 nval = gBuffer[index].nor;
+        glm::vec3 cval = src[index];
+
+        float cum_w = 0.0f;
+
+        for (size_t i = 0; i < 25; ++i) {
+            int px = x + dev_offset[2 * i] * filterSize;
+            int py = y + dev_offset[2 * i + 1] * filterSize;
+
+            if (px < 0 || px >= resolution.x || py < 0 || py >= resolution.y) {
+                continue;
+            }
+
+            int pIndex = px + (py * resolution.x);
+
+            glm::vec3 ctmp = src[pIndex];
+            glm::vec3 t = iter_i * (cval - ctmp);
+            float dist2 = glm::dot(t, t);
+            float c_w = exp(-(dist2) / pow(col_weight, (float)filterSize));
+
+            glm::vec3 ptmp = gBuffer[pIndex].pos;
+            t = pval - ptmp;
+            dist2 = glm::dot(t, t);
+            float p_w = exp(-(dist2) / pow(pos_weight, (float)filterSize));
+
+            glm::vec3 ntmp = gBuffer[pIndex].nor;
+            t = nval - ntmp;
+            dist2 = glm::max(glm::dot(t, t), 0.0f);
+            float n_w = exp(-(dist2) / pow(nor_weight, (float)filterSize));
+
+            float weight = c_w * p_w * n_w;
+            sum += ctmp * weight * dev_kernel[i];
+            cum_w += weight * dev_kernel[i];
+
+        }
+
+        sum = sum / (cum_w);
+        dst[index] = sum;
+    }
+}
+
 /**
 * Generate PathSegments with rays from the camera through the screen into the
 * scene, which is the first bounce of rays.
@@ -277,11 +382,17 @@ __global__ void generateGBuffer (
   int num_paths,
   ShadeableIntersection* shadeableIntersections,
 	PathSegment* pathSegments,
-  GBufferPixel* gBuffer) {
+  GBufferPixel* gBuffer,
+    Material* materials) {
   int idx = blockIdx.x * blockDim.x + threadIdx.x;
   if (idx < num_paths)
   {
+      ShadeableIntersection intersection = shadeableIntersections[idx];
+      PathSegment segment = pathSegments[idx];
+
     gBuffer[idx].t = shadeableIntersections[idx].t;
+    gBuffer[idx].pos = intersection.t * segment.ray.direction + segment.ray.origin;
+    gBuffer[idx].nor = intersection.surfaceNormal;
   }
 }
 
@@ -297,6 +408,7 @@ __global__ void finalGather(int nPaths, glm::vec3 * image, PathSegment * iterati
 	}
 }
 
+
 /**
  * Wrapper for the __global__ call that sets up the kernel calls and does a ton
  * of memory management
@@ -379,7 +491,7 @@ void pathtrace(int frame, int iter) {
 	cudaDeviceSynchronize();
 
   if (depth == 0) {
-    generateGBuffer<<<numblocksPathSegmentTracing, blockSize1d>>>(num_paths, dev_intersections, dev_paths, dev_gBuffer);
+    generateGBuffer<<<numblocksPathSegmentTracing, blockSize1d>>>(num_paths, dev_intersections, dev_paths, dev_gBuffer, dev_materials);
   }
 
 	depth++;
@@ -400,6 +512,19 @@ void pathtrace(int frame, int iter) {
 
     ///////////////////////////////////////////////////////////////////////////
 
+    if (denoising) {
+        const Camera& cam = hst_scene->state.camera;
+        const dim3 blockSize2d(8, 8);
+        const dim3 blocksPerGrid2d(
+            (cam.resolution.x + blockSize2d.x - 1) / blockSize2d.x,
+            (cam.resolution.y + blockSize2d.y - 1) / blockSize2d.y);
+
+        for (size_t i = 1; i <= fs; ++i) {
+            atrous << <blocksPerGrid2d, blockSize2d >> > (i, pw, nw, cw, cam.resolution, dev_gBuffer, dev_image, dev_image_clean, 1.0f / (float) iter);
+            std::swap(dev_image, dev_image_clean);
+        }
+    }
+
     // CHECKITOUT: use dev_image as reference if you want to implement saving denoised images.
     // Otherwise, screenshots are also acceptable.
     // Retrieve image from GPU

src/pathtrace.h

@@ -8,3 +8,4 @@ void pathtraceFree();
 void pathtrace(int frame, int iteration);
 void showGBuffer(uchar4 *pbo);
 void showImage(uchar4 *pbo, int iter);
+void updateData(bool denoise, int filterSize, float posWeight, float norWeight, float colWeight);