diff --git a/README.md b/README.md
index f044c821..8163c78b 100644
--- a/README.md
+++ b/README.md
@@ -1,13 +1,159 @@
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)
+* Xiaoyue Ma
+ * [LinkedIn](https://www.linkedin.com/in/xiaoyue-ma-6b268b193/)
+* Tested on: Windows 10, i7-12700H @ 2.30 GHz 16GB, GTX3060 8GB
-### (TODO: Your README)
+# Overview
+This project implements a CUDA-powered pathtracing denoiser that utilizes geometry buffers (G-buffers) to steer a smoothing filter. Drawing inspiration from the study "[Edge-Avoiding A-Trous Wavelet Transform for Fast Global Illumination Filtering](https://jo.dreggn.org/home/2010_atrous.pdf)", this denoiser enhances the appearance of pathtraced images by producing smoother results with fewer samples per pixel.
-*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.
+
+ 
+
+By iteratively applying sparse blurs of increasing sizes, one can approximate the effects of a Gaussian filter. Instead of using a large filter, this method utilizes a smaller one but spaces out the samples that pass through it. Interestingly, to achieve a greater blur size, it doesn't require enlarging the filter; instead, it demands more iterations. This approach optimizes performance while still achieving the desired blurring effect.
+
+
+ 
+
+
+# Features
+## Core features
+##### **A-Trous**
+
+A-Trous Filtering enhances our denoiser's efficiency by applying iterative sparse blurs, achieving large filter results with a smaller filter size.
+
+|No Filter | Filter Size = 16 | Filter Size = 64 |
+|:-----: | :-----: |:-----: |
+| |  |  |
+
+##### **Edge-Avoiding A-Trous**
+
+A-Trous Filtering can blur vital details, but using G-buffer information, we adjust blurring on sharp edges to preserve key image elements, optimizing denoising.
+
+|No Filter | A-Trous (64) | A-Trous with Edge-Avoiding (64) |
+|:-----: | :-----: |:-----: |
+| |  |  |
+
+##### **G-Buffer**
+
+Using normal, position, and time intersection data as weights, we minimize edge blurring during application. This data is viewable via the "Show GBuffer" option in the GUI.
+
+| Time to Intersect |Normal | Position |
+|:-----: | :-----: |:-----: |
+| |  |  |
+
+## Extra Feature
+
+#### Gaussian Filtering
+
+
+The Gaussian Filter calculates a pixel's new color by averaging its neighbors, giving more weight to those nearby. In my tests, it generated a more blurred image when edge-avoiding was off and a marginally noisy one when on.
+
+|No Filter | A-Trous (16) | Gaussian (16)|
+|:-----: | :-----: |:-----: |
+| |  |  |
+
+
+| A-Trous Edge-Avoiding(64) | Gaussian Edge-Avoiding(64)|
+|:-----: | :-----:|
+ |  |
+
+
+
+# Performance Analysis
+
+### **How much time denoising adds to the renders**
+
+The denoiser activates after the path tracer's image rendering, with its runtime being influenced by image resolution and filter size, not scene complexity. My tests show that for an 800x800 image with an 80x80 filter, the denoising time remains consistent regardless of iteration count.
+
+
+ 
+
+
+
+### **How denoising influences the number of iterations needed to get an "acceptably smooth" result**
+
+Perceptions of "smoothness" differ among individuals and can be influenced by various image factors. In the **'cornell_ceiling_light'** test, a smooth appearance was achieved at 300 iterations without denoising. With denoising, only 150 iterations were needed, marking a 100% reduction. While the benefit of denoising can depend on the scene, it notably reduces required iterations overall.
+
+No Denoising, 300 iterations | Denoised, 150 iterations
+:----------:|:-----------:
+ | 
+
+### **How denoising at different resolutions impacts runtime**
+
+The runtime for denoising rises with resolution, but the increase isn't linear. For instance, even though a significant resolution jump (from 200x200 to 800x800) is made, the runtime only multiplies by seven. As the resolution increases, the denoising process requires more time, attributed to the higher pixel count and added A-trous filter iterations. However, this growth in runtime doesn't scale proportionally with the resolution.
+
+
+ 
+
+
+
+### **How varying filter sizes affect performance**
+
+Predictably, an increase in filter size results in an extended denoising runtime because more neighboring pixels get sampled to determine each pixel's new color. In a chart derived from the test at an 800x800 resolution, while there's a relationship between additional time and filter size, it's not strictly linear.
+
+
+ 
+
+
+### **How visual results vary with filter size -- does the visual quality scale uniformly with filter size?**
+
+For images from the test scene at 10 iteration, denoising sees a marked improvement from 32 to 64 and a discernible one from 16 to 32. However, further increments offer limited visual benefits.
+
+Filter Size 4 | Filter Size 16 | Filter Size 32 | Filter Size 64 |
+:-----:|:-----:|:-----:|:-----:|
+ |  |  |  |
+
+### **How effective/ineffective is this method with different material types**
+
+This technique excels with diffuse materials, but often results in a softer appearance as the "roughness" diminishes from color smoothing. Its efficacy is diminished with reflections, causing them to blur noticeably and reducing the material's shine.
+
+Diffuse | Specular | Imperfect Specular
+:----------:|:-----------:|:-----------:
+ |  | 
+
+
+### **How do results compare across different scenes**
+**For example, between `cornell.txt` and `cornell_ceiling_light.txt`. Does one scene produce better denoised results? Why or why not?**
+
+The efficacy of denoising is scene-dependent. In the cornell_ceiling_light setup, it shows promising results, whereas in the standard cornell scene, it's less effective. Uniformly lit scenes with limited color variations favor denoisers, as consistent illumination aids in rapid light emission computations. In contrast, darker scenes with infrequent bright spots introduce more noise, posing challenges to the denoising process.
+
+Cornell Scene | Light Cornell Scene
+:----------:|:-----------:
+ | 
+
+
+### **A-Trous Filtering vs. Gaussian Filtering**
+
+In performance comparisons, A-Trous Filtering notably surpasses Gaussian Filtering. A-Trous runtime grows almost linearly with increasing filter size or resolution, whereas Gaussian's runtime rises exponentially.
+
+
+ 
+
+
+### Confuse
+
+In this paper, the right implementation code is like:
+```
+glm::vec3 ntmp = gBuffer[idxtmp].nor;
+t = nval - ntmp;
+dist2 = glm::max((glm::dot(t, t) / (stepWidth * stepWidth)), 0.0f);
+float n_w = glm::min(glm::exp(-dist2 / n_phi), 1.0f);
+```
+
+But I found that just applied
+```
+glm::vec3 ntmp = gBuffer[idxtmp].nor;
+t = nval - ntmp;
+dist2 = glm::dot(t, t);
+float n_w = glm::min(glm::exp(-dist2 / n_phi), 1.0f);
+```
+will make the denoise effect better.
+
+Original Version | Simple Code Version
+:----------:|:-----------:
+ | 
\ No newline at end of file
diff --git a/img/10AEdge16.png b/img/10AEdge16.png
new file mode 100644
index 00000000..52114708
Binary files /dev/null and b/img/10AEdge16.png differ
diff --git a/img/10AEdge16extra.png b/img/10AEdge16extra.png
new file mode 100644
index 00000000..3b497c31
Binary files /dev/null and b/img/10AEdge16extra.png differ
diff --git a/img/10AEdge32.png b/img/10AEdge32.png
new file mode 100644
index 00000000..2f01a882
Binary files /dev/null and b/img/10AEdge32.png differ
diff --git a/img/10AEdge4.png b/img/10AEdge4.png
new file mode 100644
index 00000000..877d949c
Binary files /dev/null and b/img/10AEdge4.png differ
diff --git a/img/10AEdge64.png b/img/10AEdge64.png
new file mode 100644
index 00000000..f1057a66
Binary files /dev/null and b/img/10AEdge64.png differ
diff --git a/img/10AnoEdge16.png b/img/10AnoEdge16.png
new file mode 100644
index 00000000..0ff6e975
Binary files /dev/null and b/img/10AnoEdge16.png differ
diff --git a/img/10AnoEdge32.png b/img/10AnoEdge32.png
new file mode 100644
index 00000000..247faa5a
Binary files /dev/null and b/img/10AnoEdge32.png differ
diff --git a/img/10AnoEdge64.png b/img/10AnoEdge64.png
new file mode 100644
index 00000000..7737d127
Binary files /dev/null and b/img/10AnoEdge64.png differ
diff --git a/img/10GEdge64.png b/img/10GEdge64.png
new file mode 100644
index 00000000..a5565fae
Binary files /dev/null and b/img/10GEdge64.png differ
diff --git a/img/10GnoEdge16.png b/img/10GnoEdge16.png
new file mode 100644
index 00000000..52332b05
Binary files /dev/null and b/img/10GnoEdge16.png differ
diff --git a/img/10cornell.png b/img/10cornell.png
new file mode 100644
index 00000000..1d7cac42
Binary files /dev/null and b/img/10cornell.png differ
diff --git a/img/10cornell1.png b/img/10cornell1.png
new file mode 100644
index 00000000..2ce3aa78
Binary files /dev/null and b/img/10cornell1.png differ
diff --git a/img/10ori.png b/img/10ori.png
new file mode 100644
index 00000000..83715091
Binary files /dev/null and b/img/10ori.png differ
diff --git a/img/150de.png b/img/150de.png
new file mode 100644
index 00000000..0aa0e6f0
Binary files /dev/null and b/img/150de.png differ
diff --git a/img/300ori.png b/img/300ori.png
new file mode 100644
index 00000000..dd7ca898
Binary files /dev/null and b/img/300ori.png differ
diff --git a/img/alg.png b/img/alg.png
new file mode 100644
index 00000000..42e02ece
Binary files /dev/null and b/img/alg.png differ
diff --git a/img/alg2.png b/img/alg2.png
new file mode 100644
index 00000000..064cd1bb
Binary files /dev/null and b/img/alg2.png differ
diff --git a/img/chart1.png b/img/chart1.png
new file mode 100644
index 00000000..a9433c6f
Binary files /dev/null and b/img/chart1.png differ
diff --git a/img/chart2.png b/img/chart2.png
new file mode 100644
index 00000000..ec444e08
Binary files /dev/null and b/img/chart2.png differ
diff --git a/img/chart3.png b/img/chart3.png
new file mode 100644
index 00000000..06917dd3
Binary files /dev/null and b/img/chart3.png differ
diff --git a/img/chart4.png b/img/chart4.png
new file mode 100644
index 00000000..0f580096
Binary files /dev/null and b/img/chart4.png differ
diff --git a/img/diffuse.png b/img/diffuse.png
new file mode 100644
index 00000000..bd13158a
Binary files /dev/null and b/img/diffuse.png differ
diff --git a/img/imperfact.png b/img/imperfact.png
new file mode 100644
index 00000000..3b9d8daf
Binary files /dev/null and b/img/imperfact.png differ
diff --git a/img/nor.png b/img/nor.png
new file mode 100644
index 00000000..639a60ab
Binary files /dev/null and b/img/nor.png differ
diff --git a/img/pos.png b/img/pos.png
new file mode 100644
index 00000000..4a1d887c
Binary files /dev/null and b/img/pos.png differ
diff --git a/img/time.png b/img/time.png
new file mode 100644
index 00000000..4ebac440
Binary files /dev/null and b/img/time.png differ
diff --git a/scenes/cornell.txt b/scenes/cornell.txt
index 83ff8202..8703405d 100644
--- a/scenes/cornell.txt
+++ b/scenes/cornell.txt
@@ -48,11 +48,29 @@ REFR 0
REFRIOR 0
EMITTANCE 0
+MATERIAL 5
+RGB .527 .804 .976
+SPECEX 0
+SPECRGB .98 .98 .98
+REFL 0
+REFR 0
+REFRIOR 0.71
+EMITTANCE 0
+
+MATERIAL 6
+RGB .527 .804 .976
+SPECEX 0
+SPECRGB .98 .98 .98
+REFL 1
+REFR 1
+REFRIOR 1.33
+EMITTANCE 0
+
// Camera
CAMERA
RES 800 800
FOVY 45
-ITERATIONS 5000
+ITERATIONS 10
DEPTH 8
FILE cornell
EYE 0.0 5 10.5
diff --git a/scenes/cornell_ceiling_light.txt b/scenes/cornell_ceiling_light.txt
index 15af5f19..5b11471e 100644
--- a/scenes/cornell_ceiling_light.txt
+++ b/scenes/cornell_ceiling_light.txt
@@ -48,6 +48,24 @@ REFR 0
REFRIOR 0
EMITTANCE 0
+MATERIAL 5
+RGB .527 .804 .976
+SPECEX 0
+SPECRGB .98 .98 .98
+REFL 0
+REFR 0
+REFRIOR 0.71
+EMITTANCE 0
+
+MATERIAL 6
+RGB .527 .804 .976
+SPECEX 0
+SPECRGB .98 .98 .98
+REFL 1
+REFR 1
+REFRIOR 1.33
+EMITTANCE 0
+
// Camera
CAMERA
RES 800 800
diff --git a/src/interactions.h b/src/interactions.h
index 144a9f5b..74edb819 100644
--- a/src/interactions.h
+++ b/src/interactions.h
@@ -1,14 +1,80 @@
#pragma once
#include "intersections.h"
+#include
+#include
+#include
+#define TWO_PI 6.2831853071795864769252866
+#define IMPERFECT_SPECULAR 1
+#define hasSubsurface 0
+
+class BSSRDF {
+ typedef glm::vec3 Vec3;
+
+public:
+ float scattering_coefficient;//0.001 - 5.0
+ float absorption_coefficient;//0.1 - 10.0
+ float refraction_index;//1.0 - 3.0
+
+ __host__ __device__ BSSRDF(float scattering, float absorption, float refraction)
+ : scattering_coefficient(scattering), absorption_coefficient(absorption), refraction_index(refraction) {}
+
+ __host__ __device__ float fresnel(float incident_angle) const {
+ float r0 = (1 - refraction_index) / (1 + refraction_index);
+ r0 = r0 * r0;
+ return r0 + (1 - r0) * pow(1 - incident_angle, 5);
+ }
+
+ __host__ __device__ float phase_function(float cos_scatter_angle) const {
+ float g = 0.5;
+ float denominator = 1.0 + g * g - 2.0 * g * cos_scatter_angle;
+ return (1.0 - g * g) / (denominator * sqrt(denominator));
+ }
+
+ __host__ __device__ float single_scattering(const Vec3& scattering_direction, const Vec3& normal) const {
+ float theta = glm::dot(scattering_direction, normal);
+ float adjusted_scattering = scattering_coefficient * (0.5 + 0.5 * theta);
+ return adjusted_scattering;
+ }
+
+ __host__ __device__ float multiple_scattering(const Vec3& view_direction, const Vec3& new_point, const Vec3& original_point, const Vec3& light_direction) const {
+ float distance = glm::length(new_point - original_point);
+ float decay = exp(-absorption_coefficient * distance);
+
+ float cos_scatter_angle = glm::dot(-view_direction, light_direction);
+ float phase_factor = phase_function(cos_scatter_angle);
+
+ float result = decay * phase_factor;
+ return glm::clamp(result, 0.0f, 1.0f);
+ }
+
+ __host__ __device__
+ float evaluate(const Vec3& light_direction, const Vec3& view_direction, const Vec3& point, const Vec3& original_point) const {
+ float incident_angle = glm::dot(light_direction, view_direction);
+ float fresnel_term = fresnel(incident_angle);
+ float single_scatter = single_scattering(-view_direction, light_direction);
+ float multiple_scatter = multiple_scattering(-view_direction, point, original_point, light_direction);
+
+ // Ensure energy conservation
+ float totalScatter = single_scatter + multiple_scatter;
+ if (totalScatter > 1.0f) {
+ single_scatter /= totalScatter;
+ multiple_scatter /= totalScatter;
+ }
+
+ return fresnel_term * (single_scatter + multiple_scatter);
+ }
+
+};
+// CHECKITOUT
/**
* Computes a cosine-weighted random direction in a hemisphere.
* Used for diffuse lighting.
*/
__host__ __device__
glm::vec3 calculateRandomDirectionInHemisphere(
- glm::vec3 normal, thrust::default_random_engine &rng) {
+ glm::vec3 normal, thrust::default_random_engine& rng) {
thrust::uniform_real_distribution u01(0, 1);
float up = sqrt(u01(rng)); // cos(theta)
@@ -23,9 +89,11 @@ glm::vec3 calculateRandomDirectionInHemisphere(
glm::vec3 directionNotNormal;
if (abs(normal.x) < SQRT_OF_ONE_THIRD) {
directionNotNormal = glm::vec3(1, 0, 0);
- } else if (abs(normal.y) < SQRT_OF_ONE_THIRD) {
+ }
+ else if (abs(normal.y) < SQRT_OF_ONE_THIRD) {
directionNotNormal = glm::vec3(0, 1, 0);
- } else {
+ }
+ else {
directionNotNormal = glm::vec3(0, 0, 1);
}
@@ -40,23 +108,199 @@ glm::vec3 calculateRandomDirectionInHemisphere(
+ sin(around) * over * perpendicularDirection2;
}
+__host__ __device__
+float calculateFresnelReflectance(const glm::vec3& rayDir, glm::vec3& normal, float indexOfRefraction) {
+ float cosThetaI = glm::dot(rayDir, normal);
+ float R1 = indexOfRefraction;
+ float R2 = 1.0f;
+
+ if (cosThetaI < 0.0f) {
+ cosThetaI = -cosThetaI;
+ std::swap(R1, R2);
+ normal = -normal;
+ }
+
+ //Schlick
+ float R0 = pow((R1 - R2) / (R1 + R2), 2);
+ return R0 + (1 - R0) * pow(1 - cosThetaI, 5);
+}
+
+
+__host__ __device__
+float sampleHG(float g, thrust::default_random_engine& rng) {
+ thrust::uniform_real_distribution u01(0, 1);
+
+ float random_val = u01(rng);
+ if (std::abs(g) < 1e-3)
+ return 2 * random_val - 1;
+
+ float term = (1 - g * g) / (1 - g + 2 * g * random_val);
+ return (1 + g * g - term * term) / (2 * g);
+}
/**
- * 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.ray.direction = newDirection;
- pathSegment.ray.origin = intersect + (newDirection * 0.0001f);
+ PathSegment& pathSegment,
+ glm::vec3 intersect,
+ glm::vec3 normal,
+ const Material& m,
+ thrust::default_random_engine& rng,
+ Geom* lights, int lightNum) {
+
+ if (pathSegment.remainingBounces <= 0) return;
+
+ thrust::uniform_real_distribution u01(0, 1);
+ thrust::uniform_real_distribution u02(0, 1);
+ float random = u01(rng);
+ float random2 = u02(rng);
+ glm::vec3 newDir;
+
+ if (m.hasReflective && m.hasRefractive) {
+ float cosThetaI = glm::dot(pathSegment.ray.direction, normal);
+ if (cosThetaI < 0.0f) {
+ cosThetaI = -cosThetaI;
+ normal = -normal;
+ }
+
+ float R = calculateFresnelReflectance(pathSegment.ray.direction, normal, m.indexOfRefraction);
+
+ if (random < R) {
+ newDir = glm::reflect(pathSegment.ray.direction, normal);
+ }
+ else {
+ float eta = cosThetaI > 0 ? 1.0f / m.indexOfRefraction : m.indexOfRefraction;
+ newDir = glm::refract(pathSegment.ray.direction, normal, eta);
+ }
+
+ // Russian Roulette with color-based termination
+ float continueProbability = glm::max(glm::max(pathSegment.color.r, pathSegment.color.g), pathSegment.color.b);
+ if (u01(rng) > continueProbability) {
+ pathSegment.color = glm::vec3(0.0f);
+ pathSegment.remainingBounces--;
+ return;
+ }
+ // Adjust color due to the probability to counteract the average reduction
+ pathSegment.color /= continueProbability;
+#if IMPERFECT_SPECULAR
+ // Imperfect specular reflection/refraction
+ newDir += calculateRandomDirectionInHemisphere(normal, rng) * 0.3f;
+ newDir = glm::normalize(newDir);
+#endif
+ pathSegment.color *= m.specular.color * m.color;
+ }
+ else if (m.hasReflective) {
+#if IMPERFECT_SPECULAR
+ float fresnelReflectance = calculateFresnelReflectance(-pathSegment.ray.direction, normal, m.indexOfRefraction);
+ if (random < fresnelReflectance)
+ {
+ newDir = glm::reflect(pathSegment.ray.direction, normal);
+ pathSegment.color *= m.specular.color;
+ }
+ else
+ {
+ newDir = calculateRandomDirectionInHemisphere(normal, rng);
+ }
+ pathSegment.color *= m.color;
+#else
+ newDir = glm::reflect(pathSegment.ray.direction, normal);
+ pathSegment.color *= m.specular.color * m.color;
+#endif
+ }
+ else if (m.hasRefractive) {
+ float cosTheta = glm::dot(-pathSegment.ray.direction, normal);
+ bool entering = cosTheta > 0;
+ float eta = entering ? 1.0f / m.indexOfRefraction : m.indexOfRefraction;
+ newDir = glm::refract(pathSegment.ray.direction, normal, eta);
+ if (glm::length(newDir) < 0.01f) { // Check for total internal reflection
+ newDir = glm::reflect(pathSegment.ray.direction, normal);
+ }
+ pathSegment.color *= m.specular.color * m.color;
+ }
+ else if (hasSubsurface) {
+ BSSRDF subsurfaceModel(1.2f, 1.0f, 1.3f);
+ thrust::uniform_real_distribution u01(0, 1);
+ float random2 = u01(rng);
+ bool outside = glm::dot(pathSegment.ray.direction, normal) > 0;
+ //outside = true;
+ if (!outside) { // Inside the object
+ float fresnelReflectance = calculateFresnelReflectance(pathSegment.ray.direction, normal, subsurfaceModel.refraction_index);
+ if (random2 < fresnelReflectance) {
+ glm::vec3 inverse_normal = -normal;
+ glm::vec3 insideDirection = calculateRandomDirectionInHemisphere(inverse_normal, rng);
+ glm::vec3 new_color = pathSegment.color * subsurfaceModel.single_scattering(insideDirection, inverse_normal) * m.color;
+ pathSegment.color = new_color;
+ pathSegment.ray.direction = insideDirection;
+ pathSegment.ray.origin = intersect + insideDirection * 0.001f;
+
+ }
+ else {
+ // Multiple Scattering
+ glm::vec3 scattering_direction = calculateRandomDirectionInHemisphere(normal, rng);
+
+ thrust::uniform_real_distribution u02(0, lightNum - 1);
+ Geom light = lights[u02(rng)];
+ glm::vec3 lightPos = glm::vec3(light.translation * glm::vec3(u01(rng), u01(rng), u01(rng)));
+
+ glm::vec3 new_origin = intersect + scattering_direction * 0.02f;
+ glm::vec3 lightDir = glm::normalize(lightPos - intersect);
+
+ glm::vec3 view_direction = glm::normalize(pathSegment.ray.origin - intersect);
+ glm::vec3 new_color = pathSegment.color * subsurfaceModel.multiple_scattering(view_direction, new_origin, intersect, lightDir) * m.color;
+ //glm::vec3 new_color = pathSegment.color * exp(-(intersect - pathSegment.ray.origin) * m.specular.exponent) * m.color;
+ pathSegment.color = new_color;
+ pathSegment.ray.origin = new_origin;
+ pathSegment.ray.direction = scattering_direction;
+ }
+
+ }
+ else { // Outside the object
+ glm::vec3 incomingDir = -pathSegment.ray.direction;
+ float fresnelEffect = subsurfaceModel.fresnel(glm::dot(incomingDir, normal));
+ glm::vec3 refractedDir = glm::refract(incomingDir, normal, 1.0f / subsurfaceModel.refraction_index);
+
+ glm::vec3 new_color = pathSegment.color * fresnelEffect * subsurfaceModel.single_scattering(refractedDir, intersect) * m.color;
+ pathSegment.color = new_color;
+
+ glm::vec3 new_origin = intersect + glm::normalize(pathSegment.ray.direction) * 0.001f + normal;
+ pathSegment.ray.origin = new_origin + refractedDir * 0.001f;
+ pathSegment.ray.direction = refractedDir;
+ }
+ pathSegment.remainingBounces--;
+ return;
+ }
+ else {
+ newDir = calculateRandomDirectionInHemisphere(normal, rng);
+ pathSegment.color *= m.color;
+ }
+
+ glm::vec3 newOrigin = intersect + 0.001f * newDir;
+
+ pathSegment.ray.origin = newOrigin;
+ pathSegment.ray.direction = newDir;
+ pathSegment.remainingBounces--;
+ return;
}
diff --git a/src/main.cpp b/src/main.cpp
index 4092ae4a..f542b20f 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -164,10 +164,14 @@ void runCuda() {
int frame = 0;
pathtrace(frame, iteration);
}
-
- if (ui_showGbuffer) {
+ if (ui_denoise) {
+ denoise(ui_colorWeight, ui_normalWeight, ui_positionWeight, ui_filterSize);
+ showDenoisedImage(pbo_dptr, iteration);
+ }
+ else if (ui_showGbuffer) {
showGBuffer(pbo_dptr);
- } else {
+ }
+ else {
showImage(pbo_dptr, iteration);
}
diff --git a/src/pathtrace.cu b/src/pathtrace.cu
index 23e5f909..1913bfba 100644
--- a/src/pathtrace.cu
+++ b/src/pathtrace.cu
@@ -16,6 +16,11 @@
#define ERRORCHECK 1
+#define GAUSS_ON_ATROUS_OFF 0
+#define AVOID_EDGE 1
+#define SHOW_POS_BUFFER 0
+#define SHOW_NOR_BUFFER 0
+
#define FILENAME (strrchr(__FILE__, '/') ? strrchr(__FILE__, '/') + 1 : __FILE__)
#define checkCUDAError(msg) checkCUDAErrorFn(msg, FILENAME, __LINE__)
void checkCUDAErrorFn(const char *msg, const char *file, int line) {
@@ -73,7 +78,7 @@ __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;
+ float timeToIntersect = gBuffer[index].t * 255.0;
pbo[index].w = 0;
pbo[index].x = timeToIntersect;
@@ -82,6 +87,35 @@ __global__ void gbufferToPBO(uchar4* pbo, glm::ivec2 resolution, GBufferPixel* g
}
}
+__global__ void gbufferToPBONormal(uchar4* pbo, glm::ivec2 resolution, GBufferPixel* gBuffer) {
+ int x = (blockIdx.x * blockDim.x) + threadIdx.x;
+ int y = (blockIdx.y * blockDim.y) + threadIdx.y;
+
+ if (x < resolution.x && y < resolution.y) {
+ int index = x + (y * resolution.x);
+ glm::vec3 nor = gBuffer[index].nor * 255.f;
+ pbo[index].w = 0;
+ pbo[index].x = abs(nor.x);
+ pbo[index].y = abs(nor.y);
+ pbo[index].z = abs(nor.z);
+ }
+}
+
+__global__ void gbufferToPBOPosition(uchar4* pbo, glm::ivec2 resolution, GBufferPixel* gBuffer) {
+ int x = (blockIdx.x * blockDim.x) + threadIdx.x;
+ int y = (blockIdx.y * blockDim.y) + threadIdx.y;
+
+ if (x < resolution.x && y < resolution.y) {
+ int index = x + (y * resolution.x);
+ glm::vec3 pos = gBuffer[index].pos * 25.5f;
+ pbo[index].w = 0;
+ pbo[index].x = abs(pos.x);
+ pbo[index].y = abs(pos.y);
+ pbo[index].z = abs(pos.z);
+
+ }
+}
+
static Scene * hst_scene = NULL;
static glm::vec3 * dev_image = NULL;
static Geom * dev_geoms = NULL;
@@ -90,7 +124,9 @@ 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
-// ...
+static glm::vec3* dev_denoised_image = NULL;
+static glm::vec3* dev_denoised_buffer = NULL;
+static Geom* lights = NULL;
void pathtraceInit(Scene *scene) {
hst_scene = scene;
@@ -114,6 +150,13 @@ void pathtraceInit(Scene *scene) {
cudaMalloc(&dev_gBuffer, pixelcount * sizeof(GBufferPixel));
// TODO: initialize any extra device memeory you need
+ cudaMalloc(&dev_denoised_image, pixelcount * sizeof(glm::vec3));
+ cudaMemset(dev_denoised_image, 0, pixelcount * sizeof(glm::vec3));
+ cudaMalloc(&dev_denoised_buffer, pixelcount * sizeof(glm::vec3));
+ cudaMemset(dev_denoised_buffer, 0, pixelcount * sizeof(glm::vec3));
+
+ cudaMalloc(&lights, scene->lightNum * sizeof(Geom));
+ cudaMemcpy(lights, scene->lights.data(), scene->lights.size() * sizeof(Geom), cudaMemcpyHostToDevice);
checkCUDAError("pathtraceInit");
}
@@ -126,7 +169,9 @@ void pathtraceFree() {
cudaFree(dev_intersections);
cudaFree(dev_gBuffer);
// TODO: clean up any extra device memory you created
-
+ cudaFree(dev_denoised_image);
+ cudaFree(dev_denoised_buffer);
+ cudaFree(lights);
checkCUDAError("pathtraceFree");
}
@@ -177,7 +222,7 @@ __global__ void computeIntersections(
float t;
glm::vec3 intersect_point;
- glm::vec3 normal;
+ glm::vec3 nor;
float t_min = FLT_MAX;
int hit_geom_index = -1;
bool outside = true;
@@ -207,7 +252,7 @@ __global__ void computeIntersections(
t_min = t;
hit_geom_index = i;
intersect_point = tmp_intersect;
- normal = tmp_normal;
+ nor = tmp_normal;
}
}
@@ -220,7 +265,7 @@ __global__ void computeIntersections(
//The ray hits something
intersections[path_index].t = t_min;
intersections[path_index].materialId = geoms[hit_geom_index].materialid;
- intersections[path_index].surfaceNormal = normal;
+ intersections[path_index].surfaceNormal = nor;
}
}
}
@@ -231,6 +276,8 @@ __global__ void shadeSimpleMaterials (
, ShadeableIntersection * shadeableIntersections
, PathSegment * pathSegments
, Material * materials
+ , Geom* lights
+ , int lightNum
)
{
int idx = blockIdx.x * blockDim.x + threadIdx.x;
@@ -258,7 +305,7 @@ __global__ void shadeSimpleMaterials (
else {
segment.color *= materialColor;
glm::vec3 intersectPos = intersection.t * segment.ray.direction + segment.ray.origin;
- scatterRay(segment, intersectPos, intersection.surfaceNormal, material, rng);
+ scatterRay(segment, intersectPos, intersection.surfaceNormal, material, rng, lights, lightNum);
}
// If there was no intersection, color the ray black.
// Lots of renderers use 4 channel color, RGBA, where A = alpha, often
@@ -282,6 +329,8 @@ __global__ void generateGBuffer (
if (idx < num_paths)
{
gBuffer[idx].t = shadeableIntersections[idx].t;
+ gBuffer[idx].nor = shadeableIntersections[idx].surfaceNormal;
+ gBuffer[idx].pos = pathSegments[idx].ray.origin + pathSegments[idx].ray.direction * gBuffer[idx].t;
}
}
@@ -389,7 +438,9 @@ void pathtrace(int frame, int iter) {
num_paths,
dev_intersections,
dev_paths,
- dev_materials
+ dev_materials,
+ lights,
+ hst_scene->lightNum
);
iterationComplete = depth == traceDepth;
}
@@ -409,6 +460,156 @@ void pathtrace(int frame, int iter) {
checkCUDAError("pathtrace");
}
+__global__ void ATrousDenoise(float c_phi, float n_phi, float p_phi, glm::ivec2 resolution, GBufferPixel* gBuffer, int stepWidth,
+ glm::vec3* origin_image, glm::vec3* denoised_image)
+{
+ int x = blockDim.x * blockIdx.x + threadIdx.x;
+ int y = blockDim.y * blockIdx.y + threadIdx.y;
+ if (x >= resolution.x || y >= resolution.y) return;
+
+ int idx = x + y * resolution.x;
+
+ glm::vec3 cval = origin_image[idx];
+ glm::vec3 nval = gBuffer[idx].nor;
+ glm::vec3 pval = gBuffer[idx].pos;
+
+ float cum_w = 0.0f;
+ glm::vec3 sum = glm::vec3(0.f);
+
+ static float kernel[5][5] = {
+ 0.00296902, 0.0133062, 0.0219382, 0.0133062, 0.00296902,
+ 0.0133062, 0.0596343, 0.0983203, 0.0596343, 0.0133062,
+ 0.0219382, 0.0983203, 0.162103, 0.0983203, 0.0219382,
+ 0.0133062, 0.0596343, 0.0983203, 0.0596343, 0.0133062,
+ 0.00296902, 0.0133062, 0.0219382, 0.0133062, 0.00296902,
+ };
+ for (int i = -2; i < 2; i++) {
+ for (int j = -2; j < 2; j++) {
+ glm::ivec2 uv;
+ uv.x = x + i * stepWidth;
+ uv.y = y + j * stepWidth;
+ uv = glm::clamp(uv, glm::ivec2(0, 0), glm::ivec2(resolution.x - 1, resolution.y - 1));
+
+ int idxtmp = uv.x + uv.y * resolution.x;
+ glm::vec3 ctmp = origin_image[idxtmp];
+ float weight = 1.0f;
+
+#if AVOID_EDGE
+
+ glm::vec3 t = cval - ctmp;
+ float dist2 = glm::dot(t, t);
+ float c_w = glm::min(glm::exp(-dist2 / c_phi), 1.0f);
+
+ glm::vec3 ntmp = gBuffer[idxtmp].nor;
+ t = nval - ntmp;
+ dist2 = glm::max((glm::dot(t, t) / (stepWidth * stepWidth)), 0.0f);
+ dist2 = glm::dot(t, t);
+ float n_w = glm::min(glm::exp(-dist2 / n_phi), 1.0f);
+
+ glm::vec3 ptmp = gBuffer[idxtmp].pos;
+ t = pval - ptmp;
+ dist2 = glm::dot(t, t);
+ float p_w = glm::min(glm::exp(-dist2 / p_phi), 1.0f);
+
+ weight = c_w * n_w * p_w;
+
+#endif
+ float curKernel = kernel[i+2][j+2];
+ sum += ctmp * weight * curKernel;
+ cum_w += weight * curKernel;
+ }
+ }
+ denoised_image[idx] = sum / cum_w;
+}
+
+__global__ void GaussDenoise(glm::ivec2 resolution, GBufferPixel* gBuffer,
+ float c_phi, float n_phi, float p_phi, int filtersize, glm::vec3* origin_image, glm::vec3* denoised_image)
+{
+ int x = blockDim.x * blockIdx.x + threadIdx.x;
+ int y = blockDim.y * blockIdx.y + threadIdx.y;
+
+ if (x >= resolution.x || y >= resolution.y) return;
+
+ int idx = x + y * resolution.x;
+
+ glm::vec3 cval = origin_image[idx];
+ glm::vec3 nval = gBuffer[idx].nor;
+ glm::vec3 pval = gBuffer[idx].pos;
+
+ glm::vec3 sum = glm::vec3(0.f);
+ float cum_w = 0.0f;
+
+ int bound = (filtersize - 1) / 2;
+ int s = filtersize / 6;
+
+ for (int i = -bound; i <= bound; ++i)
+ {
+ for (int j = -bound; j <= bound; ++j)
+ {
+ glm::ivec2 uv;
+ uv.x = x + i;
+ uv.y = y + j;
+
+ float weight = 1.0f;
+
+ if (uv.x >= 0 && uv.x < resolution.x && uv.y >= 0 && uv.y <= resolution.y)
+ {
+ int idxtmp = uv.x + uv.y * resolution.x;
+ glm::vec3 ctmp = origin_image[idxtmp];
+
+#if AVOID_EDGE
+
+ glm::vec3 t = cval - ctmp;
+ float dist2 = glm::dot(t, t);
+ float c_w = glm::min(glm::exp(-dist2 / c_phi), 1.0f);
+
+ glm::vec3 ntmp = gBuffer[idxtmp].nor;
+ t = nval - ntmp;
+ dist2 = glm::max(glm::dot(t, t), 0.0f);
+ float n_w = glm::min(glm::exp(-dist2 / n_phi), 1.0f);
+
+ glm::vec3 ptmp = gBuffer[idxtmp].pos;
+ t = pval - ptmp;
+ dist2 = glm::dot(t, t);
+ float p_w = glm::min(glm::exp(-dist2 / p_phi), 1.0f);
+
+ weight = c_w * n_w * p_w;
+
+#endif
+ float curKernel = (1.f / (2.f * PI * s * s)) * exp(-(i * i + j * j) / (2.f * s * s));
+ sum += ctmp * weight * curKernel;
+ cum_w += weight * curKernel;
+ }
+ }
+ }
+ denoised_image[idx] = sum / cum_w;
+}
+
+void denoise(float c_phi, float n_phi, float p_phi, float filterSize)
+{
+ Camera& cam = hst_scene->state.camera;
+ glm::ivec2 resolution = cam.resolution;
+
+ const dim3 blockSize2d(8, 8);
+ const dim3 blocksPerGrid2d(
+ (cam.resolution.x + blockSize2d.x - 1) / blockSize2d.x,
+ (cam.resolution.y + blockSize2d.y - 1) / blockSize2d.y);
+
+ cudaMemcpy(dev_denoised_buffer, dev_image, resolution.x * resolution.y * sizeof(glm::vec3), cudaMemcpyDeviceToDevice);
+
+#if GAUSS_ON_ATROUS_OFF
+ GaussDenoise << > > (resolution, dev_gBuffer, c_phi, n_phi, p_phi, filterSize, dev_denoised_buffer, dev_denoised_image);
+#else
+ for (int stepWidth = 1; 4 * stepWidth < filterSize; stepWidth *= 2)
+ {
+ ATrousDenoise << > > (c_phi, n_phi, p_phi, resolution, dev_gBuffer, stepWidth, dev_denoised_buffer, dev_denoised_image);
+ cudaDeviceSynchronize();
+ std::swap(dev_denoised_buffer, dev_denoised_image);
+ }
+#endif
+ cudaMemcpy(hst_scene->state.image.data(), dev_denoised_image, resolution.x * resolution.y * sizeof(glm::vec3), cudaMemcpyDeviceToHost);
+}
+
// CHECKITOUT: this kernel "post-processes" the gbuffer/gbuffers into something that you can visualize for debugging.
void showGBuffer(uchar4* pbo) {
const Camera &cam = hst_scene->state.camera;
@@ -418,7 +619,15 @@ void showGBuffer(uchar4* pbo) {
(cam.resolution.y + blockSize2d.y - 1) / blockSize2d.y);
// CHECKITOUT: process the gbuffer results and send them to OpenGL buffer for visualization
- gbufferToPBO<<>>(pbo, cam.resolution, dev_gBuffer);
+#if SHOW_POS_BUFFER
+ gbufferToPBO_Position << > > (pbo, cam.resolution, dev_gBuffer);
+
+#elif SHOW_NOR_BUFFER
+ gbufferToPBO_Normal << > > (pbo, cam.resolution, dev_gBuffer);
+
+#else
+ gbufferToPBO << > > (pbo, cam.resolution, dev_gBuffer);
+#endif
}
void showImage(uchar4* pbo, int iter) {
@@ -431,3 +640,14 @@ const Camera &cam = hst_scene->state.camera;
// Send results to OpenGL buffer for rendering
sendImageToPBO<<>>(pbo, cam.resolution, iter, dev_image);
}
+
+void showDenoisedImage(uchar4* pbo, int iter) {
+ 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);
+
+ // Send results to OpenGL buffer for rendering
+ sendImageToPBO << > > (pbo, cam.resolution, iter, dev_denoised_image);
+}
diff --git a/src/pathtrace.h b/src/pathtrace.h
index 9e12f440..d826815f 100644
--- a/src/pathtrace.h
+++ b/src/pathtrace.h
@@ -8,3 +8,5 @@ void pathtraceFree();
void pathtrace(int frame, int iteration);
void showGBuffer(uchar4 *pbo);
void showImage(uchar4 *pbo, int iter);
+void showDenoisedImage(uchar4* pbo, int iter);
+void denoise(float c_phi, float n_phi, float p_phi, float stepWidth);
\ No newline at end of file
diff --git a/src/preview.cpp b/src/preview.cpp
index 3ca27180..02d03c82 100644
--- a/src/preview.cpp
+++ b/src/preview.cpp
@@ -196,8 +196,8 @@ void drawGui(int windowWidth, int windowHeight) {
ImGui::NewFrame();
// Dear imgui define
- ImVec2 minSize(300.f, 220.f);
- ImVec2 maxSize((float)windowWidth * 0.5, (float)windowHeight * 0.3);
+ ImVec2 minSize(350.f, 280.f);
+ ImVec2 maxSize((float)windowWidth * 0.6, (float)windowHeight * 0.8);
ImGui::SetNextWindowSizeConstraints(minSize, maxSize);
ImGui::SetNextWindowPos(ui_hide ? ImVec2(-1000.f, -1000.f) : ImVec2(0.0f, 0.0f));
@@ -225,6 +225,9 @@ void drawGui(int windowWidth, int windowHeight) {
ImGui::Separator();
+ ImGui::Text("Application average %.3f ms/frame (%.1f FPS)", 1000.0f / ImGui::GetIO().Framerate, ImGui::GetIO().Framerate);
+ ImGui::Separator();
+
if (ImGui::Button("Save image and exit")) {
ui_saveAndExit = true;
}
diff --git a/src/scene.cpp b/src/scene.cpp
index cbae043c..79a1d069 100644
--- a/src/scene.cpp
+++ b/src/scene.cpp
@@ -84,6 +84,12 @@ int Scene::loadGeom(string objectid) {
newGeom.inverseTransform = glm::inverse(newGeom.transform);
newGeom.invTranspose = glm::inverseTranspose(newGeom.transform);
+ if (newGeom.materialid == 0)
+ {
+ this->lights.push_back(newGeom);
+ this->lightNum++;
+ }
+
geoms.push_back(newGeom);
return 1;
}
diff --git a/src/scene.h b/src/scene.h
index f29a9171..e63ecdb0 100644
--- a/src/scene.h
+++ b/src/scene.h
@@ -23,4 +23,6 @@ class Scene {
std::vector geoms;
std::vector materials;
RenderState state;
+ std::vector lights;
+ int lightNum = 0;
};
diff --git a/src/sceneStructs.h b/src/sceneStructs.h
index da7e558a..e2e3404e 100644
--- a/src/sceneStructs.h
+++ b/src/sceneStructs.h
@@ -79,4 +79,6 @@ struct ShadeableIntersection {
// What information might be helpful for guiding a denoising filter?
struct GBufferPixel {
float t;
+ glm::vec3 nor;
+ glm::vec3 pos;
};