diff --git a/README.md b/README.md
index 110697c..080bae0 100644
--- a/README.md
+++ b/README.md
@@ -3,11 +3,97 @@ 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)
+* Wanru Zhao
+ * [LinkedIn](www.linkedin.com/in/wanru-zhao).
+* Tested on: Windows 10, Intel(R) Core(TM) i7-8750H CPU@2.2GHz, GTX 1070 with Max-Q Design(Personal Laptop)
-### (TODO: Your README)
+### Final
+
+
+
+
-*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.
+### Features
+- Shading kernel with BSDF features: diffuse and perfectly specular
+- Path termination with stream compaction
+- Material sort
+- First cache
+- Refraction with fresnel effects
+- Physically based depth of field
+- Stochastic sampled antialiasing
+- Direct lighting
+- Motion blur
+
+### Analysis
+#### Stream compaction
+
+In each iteration, apply stream compaction on pathsegments according their remaining bounces after each bounce, remove the pathsegments with 0 remaining bounce, and adjust the number of launched threads. When there is no pathsegments alive left, bouncing can be terminated earlier.
+
+
+
+
+#### Cache first intersection
+
+Cache first intersection at the first bounce and first iteration. Cached data is used for initializing following iterations as results for first bounce. This skips the same intersection calculation. (Since this assuming the intersectionsof first bounce are always the same, this cannot be used with antialiasing.)
+
+
+
+From above image, we can see for depth from 6 to 16, performance with cache first intersection is slightly better than naive way, for it reduces the calculation of first bounce for iteration 2 to iteration 5000. However, when bounce count increases from 16 to 20, performance of caching is worse than naive way. I think this is because when bounce times is large enough and the number of iterations is the same, the benefit from reducing first bounce calculation can not exceed the increasing cost of branching.
+
+#### Material sort
+
+Sort the rays/path segments so that rays/paths interacting with the same material are contiguous in memory before shading.
+
+
+
+The performance of material sort is much worse than naive way, since it does the sorting each bounce each iteration, which consumes a lot of time. This might be optimized by better sorting algorithm and the following optimization of shading with different materials.
+
+#### Different materials
+
+Diffuse, perfect specular, refraction (5000 iterations, 8 bounces)
+
+
+
+#### Antialiasing
+
+Perform antialiasing by jittering ray direction calculated by camera and averaging results of multiple iterations.
+
+(5000 iterations, 8 bounces)
+
+
+
+
+
+
+
+
+
+
+#### Depth of field
+
+Jitter the ray origin and direction from camera according to lens effects.
+
+DOF
+:--:
+
+(5000 iterations, 8 bounces, len radius 0.5, focal distance 7.0)
+
+
+#### Direct lighting
+
+Take the final ray directly to a random point on an emissive object.
+
+Sample final ray direction | Without direct lighting | Direct lighting
+:--:|:--:|:--:
+ |  | 
+(5000 iterations, 8 bounces) |(5000 iterations, 8 bounces)| (5000 iterations, 8 bounces)
+
+#### Motion blur
+
+Average samples at different times in the animation.
+
+Motion Blur
+:--:
+
+(5000 iterations, 8 bounces)
diff --git a/img/final/Running-Time-Comparison-WO-First-Cache.png b/img/final/Running-Time-Comparison-WO-First-Cache.png
new file mode 100644
index 0000000..ef19107
Binary files /dev/null and b/img/final/Running-Time-Comparison-WO-First-Cache.png differ
diff --git a/img/final/anti.PNG b/img/final/anti.PNG
new file mode 100644
index 0000000..369463a
Binary files /dev/null and b/img/final/anti.PNG differ
diff --git a/img/final/cornell.2018-09-29_21-38-52z.5000samp.png b/img/final/cornell.2018-09-29_21-38-52z.5000samp.png
new file mode 100644
index 0000000..2b4fa18
Binary files /dev/null and b/img/final/cornell.2018-09-29_21-38-52z.5000samp.png differ
diff --git a/img/final/cornell.2018-10-02_01-06-49z.10samp.png b/img/final/cornell.2018-10-02_01-06-49z.10samp.png
new file mode 100644
index 0000000..1f08207
Binary files /dev/null and b/img/final/cornell.2018-10-02_01-06-49z.10samp.png differ
diff --git a/img/final/cornell.2018-10-02_02-18-51z.5000samp-streamcompact.png b/img/final/cornell.2018-10-02_02-18-51z.5000samp-streamcompact.png
new file mode 100644
index 0000000..f584bf0
Binary files /dev/null and b/img/final/cornell.2018-10-02_02-18-51z.5000samp-streamcompact.png differ
diff --git a/img/final/cornell.2018-10-02_03-08-17z.5000samp-anti.png b/img/final/cornell.2018-10-02_03-08-17z.5000samp-anti.png
new file mode 100644
index 0000000..becc833
Binary files /dev/null and b/img/final/cornell.2018-10-02_03-08-17z.5000samp-anti.png differ
diff --git a/img/final/cornell.2018-10-02_03-13-11z.5000samp.png b/img/final/cornell.2018-10-02_03-13-11z.5000samp.png
new file mode 100644
index 0000000..cdf4ac3
Binary files /dev/null and b/img/final/cornell.2018-10-02_03-13-11z.5000samp.png differ
diff --git a/img/final/cornell.2018-10-02_03-23-43z.5000samp.png b/img/final/cornell.2018-10-02_03-23-43z.5000samp.png
new file mode 100644
index 0000000..0a78ea6
Binary files /dev/null and b/img/final/cornell.2018-10-02_03-23-43z.5000samp.png differ
diff --git a/img/final/cornell.2018-10-02_03-44-40z.5000samp.png b/img/final/cornell.2018-10-02_03-44-40z.5000samp.png
new file mode 100644
index 0000000..12e4f4f
Binary files /dev/null and b/img/final/cornell.2018-10-02_03-44-40z.5000samp.png differ
diff --git a/img/final/cornell.2018-10-02_07-15-51z.10000samp.png b/img/final/cornell.2018-10-02_07-15-51z.10000samp.png
new file mode 100644
index 0000000..857cb7a
Binary files /dev/null and b/img/final/cornell.2018-10-02_07-15-51z.10000samp.png differ
diff --git a/img/final/cornell.2018-10-02_07-43-16z.10000samp.png b/img/final/cornell.2018-10-02_07-43-16z.10000samp.png
new file mode 100644
index 0000000..9b810a3
Binary files /dev/null and b/img/final/cornell.2018-10-02_07-43-16z.10000samp.png differ
diff --git a/img/final/materialsort.png b/img/final/materialsort.png
new file mode 100644
index 0000000..a414d19
Binary files /dev/null and b/img/final/materialsort.png differ
diff --git a/img/final/origin.PNG b/img/final/origin.PNG
new file mode 100644
index 0000000..b99944d
Binary files /dev/null and b/img/final/origin.PNG differ
diff --git a/img/final/raybounce.png b/img/final/raybounce.png
new file mode 100644
index 0000000..8b1d174
Binary files /dev/null and b/img/final/raybounce.png differ
diff --git a/scenes/cornell.txt b/scenes/cornell.txt
index 83ff820..0857069 100644
--- a/scenes/cornell.txt
+++ b/scenes/cornell.txt
@@ -67,6 +67,7 @@ material 0
TRANS 0 10 0
ROTAT 0 0 0
SCALE 3 .3 3
+MOTION 0
// Floor
OBJECT 1
@@ -75,6 +76,7 @@ material 1
TRANS 0 0 0
ROTAT 0 0 0
SCALE 10 .01 10
+MOTION 0
// Ceiling
OBJECT 2
@@ -83,6 +85,7 @@ material 1
TRANS 0 10 0
ROTAT 0 0 90
SCALE .01 10 10
+MOTION 0
// Back wall
OBJECT 3
@@ -91,6 +94,7 @@ material 1
TRANS 0 5 -5
ROTAT 0 90 0
SCALE .01 10 10
+MOTION 0
// Left wall
OBJECT 4
@@ -99,6 +103,7 @@ material 2
TRANS -5 5 0
ROTAT 0 0 0
SCALE .01 10 10
+MOTION 0
// Right wall
OBJECT 5
@@ -107,11 +112,13 @@ material 3
TRANS 5 5 0
ROTAT 0 0 0
SCALE .01 10 10
+MOTION 0
// Sphere
OBJECT 6
sphere
material 4
-TRANS -1 4 -1
+TRANS 0 3 0
ROTAT 0 0 0
SCALE 3 3 3
+MOTION 1
\ No newline at end of file
diff --git a/scenes/cornellbox.txt b/scenes/cornellbox.txt
new file mode 100644
index 0000000..6580ab5
--- /dev/null
+++ b/scenes/cornellbox.txt
@@ -0,0 +1,190 @@
+// Emissive material (light)
+MATERIAL 0
+RGB 0.96 0.69 0.25
+SPECEX 0
+SPECRGB 0 0 0
+REFL 0
+REFR 0
+REFRIOR 0
+EMITTANCE 20
+
+// Diffuse white
+MATERIAL 1
+RGB .98 .98 .98
+SPECEX 0
+SPECRGB 0 0 0
+REFL 0
+REFR 0
+REFRIOR 0
+EMITTANCE 0
+
+// Diffuse red
+MATERIAL 2
+RGB .85 .35 .35
+SPECEX 0
+SPECRGB .85 .35 .35
+REFL 0.6
+REFR 0
+REFRIOR 0
+EMITTANCE 0
+
+// Diffuse green
+MATERIAL 3
+RGB .35 .85 .35
+SPECEX 0
+SPECRGB 0 0 0
+REFL 0
+REFR 0
+REFRIOR 0
+EMITTANCE 0
+
+// Specular white
+MATERIAL 4
+RGB .98 .98 .98
+SPECEX 0
+SPECRGB 1 1 1
+REFL 0.2
+REFR 0.8
+REFRIOR 1.5
+EMITTANCE 0
+
+// Specular Diffuse white
+MATERIAL 5
+RGB .98 .98 .98
+SPECEX 0
+SPECRGB .98 .98 .98
+REFL 0
+REFR 0.4
+REFRIOR 2
+EMITTANCE 0
+
+// Emissive material (light)
+MATERIAL 6
+RGB 0.36 0.67 0.88
+SPECEX 0
+SPECRGB 0 0 0
+REFL 0
+REFR 0
+REFRIOR 0
+EMITTANCE 20
+
+// Specular Diffuse white
+MATERIAL 7
+RGB .98 .98 .98
+SPECEX 0
+SPECRGB .98 .98 .98
+REFL 1
+REFR 0
+REFRIOR 0
+EMITTANCE 0
+
+// Camera
+CAMERA
+RES 800 800
+FOVY 45
+ITERATIONS 10000
+DEPTH 8
+FILE cornell
+EYE 0.0 5 10.5
+LOOKAT 0 5 0
+UP 0 1 0
+
+
+// Ceiling light
+OBJECT 0
+cube
+material 0
+TRANS -5 16 0
+ROTAT 0 0 0
+SCALE 3 .3 3
+MOTION 0
+
+// Floor
+OBJECT 1
+cube
+material 1
+TRANS 0 0 0
+ROTAT 0 0 0
+SCALE 16 .01 16
+MOTION 0
+
+// Ceiling
+OBJECT 2
+cube
+material 1
+TRANS 0 16 0
+ROTAT 0 0 90
+SCALE .01 16 16
+MOTION 0
+
+// Back wall
+OBJECT 3
+cube
+material 1
+TRANS 0 8 -8
+ROTAT 0 90 0
+SCALE .01 16 16
+MOTION 0
+
+// Left wall
+OBJECT 4
+cube
+material 2
+TRANS -8 8 0
+ROTAT 0 0 0
+SCALE .01 16 16
+MOTION 0
+
+// Right wall
+OBJECT 5
+cube
+material 3
+TRANS 8 8 0
+ROTAT 0 0 0
+SCALE .01 16 16
+MOTION 0
+
+// Box
+OBJECT 6
+sphere
+material 5
+TRANS -2 1.5 -2
+ROTAT 0 0 0
+SCALE 3 3 3
+MOTION 0
+
+// Box
+OBJECT 7
+sphere
+material 7
+TRANS 2 1.5 0
+ROTAT 0 0 0
+SCALE 3 3 3
+MOTION 0
+
+// Ceiling light
+OBJECT 8
+cube
+material 6
+TRANS 5 16 0
+ROTAT 0 0 0
+SCALE 3 .3 3
+MOTION 0
+
+// Sphere
+OBJECT 9
+sphere
+material 4
+TRANS -1 1.5 2
+ROTAT 0 0 0
+SCALE 3 3 3
+MOTION 0
+
+// Back wall
+OBJECT 10
+cube
+material 1
+TRANS 0 8 8
+ROTAT 0 90 0
+SCALE .01 16 16
+MOTION 0
\ No newline at end of file
diff --git a/scenes/material.txt b/scenes/material.txt
new file mode 100644
index 0000000..0f29423
--- /dev/null
+++ b/scenes/material.txt
@@ -0,0 +1,152 @@
+// Emissive material (light)
+MATERIAL 0
+RGB 1 1 1
+SPECEX 0
+SPECRGB 0 0 0
+REFL 0
+REFR 0
+REFRIOR 0
+EMITTANCE 5
+
+// Diffuse white
+MATERIAL 1
+RGB .98 .98 .98
+SPECEX 0
+SPECRGB 0 0 0
+REFL 0
+REFR 0
+REFRIOR 0
+EMITTANCE 0
+
+// Diffuse red
+MATERIAL 2
+RGB .85 .35 .35
+SPECEX 0
+SPECRGB 0 0 0
+REFL 0
+REFR 0
+REFRIOR 0
+EMITTANCE 0
+
+// Diffuse green
+MATERIAL 3
+RGB .35 .85 .35
+SPECEX 0
+SPECRGB 0 0 0
+REFL 0
+REFR 0
+REFRIOR 0
+EMITTANCE 0
+
+// Specular white
+MATERIAL 4
+RGB .98 .98 .98
+SPECEX 0
+SPECRGB .98 .98 .98
+REFL 1
+REFR 0
+REFRIOR 0
+EMITTANCE 0
+
+// Refract white
+MATERIAL 5
+RGB .98 .98 .98
+SPECEX 0
+SPECRGB .98 .98 .98
+REFL 0
+REFR 1
+REFRIOR 1.5
+EMITTANCE 0
+
+// Camera
+CAMERA
+RES 800 800
+FOVY 45
+ITERATIONS 5000
+DEPTH 8
+FILE cornell
+EYE 0.0 5 10.5
+LOOKAT 0 5 0
+UP 0 1 0
+
+
+// Ceiling light
+OBJECT 0
+cube
+material 0
+TRANS 0 10 0
+ROTAT 0 0 0
+SCALE 3 .3 3
+MOTION 0
+
+// Floor
+OBJECT 1
+cube
+material 1
+TRANS 0 0 0
+ROTAT 0 0 0
+SCALE 10 .01 10
+MOTION 0
+
+// Ceiling
+OBJECT 2
+cube
+material 1
+TRANS 0 10 0
+ROTAT 0 0 90
+SCALE .01 10 10
+MOTION 0
+
+// Back wall
+OBJECT 3
+cube
+material 1
+TRANS 0 5 -5
+ROTAT 0 90 0
+SCALE .01 10 10
+MOTION 0
+
+// Left wall
+OBJECT 4
+cube
+material 2
+TRANS -5 5 0
+ROTAT 0 0 0
+SCALE .01 10 10
+MOTION 0
+
+// Right wall
+OBJECT 5
+cube
+material 3
+TRANS 5 5 0
+ROTAT 0 0 0
+SCALE .01 10 10
+MOTION 0
+
+// Sphere
+OBJECT 6
+sphere
+material 1
+TRANS -3 3 0
+ROTAT 0 0 0
+SCALE 3 3 3
+MOTION 1
+
+// Sphere
+OBJECT 7
+sphere
+material 4
+TRANS 0 3 2
+ROTAT 0 0 0
+SCALE 3 3 3
+MOTION 1
+
+// Sphere
+OBJECT 8
+sphere
+material 5
+TRANS 3 3 0
+ROTAT 0 0 0
+SCALE 3 3 3
+MOTION 1
\ No newline at end of file
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index a1cb3fb..5172d5f 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -15,6 +15,7 @@ set(SOURCE_FILES
"preview.cpp"
"utilities.cpp"
"utilities.h"
+ "globals.h"
)
cuda_add_library(src
diff --git a/src/globals.h b/src/globals.h
new file mode 100644
index 0000000..aace433
--- /dev/null
+++ b/src/globals.h
@@ -0,0 +1,69 @@
+#pragma once
+
+#include
+#include
+#include
+#include
+#include
+#include
+
+#include "glm/glm.hpp"
+#include "glm/gtx/norm.hpp"
+
+__device__ const float Pi = 3.14159265358979323846;
+__device__ const float TwoPi = 6.28318530717958647692;
+__device__ const float InvPi = 0.31830988618379067154;
+__device__ const float Inv2Pi = 0.15915494309189533577;
+__device__ const float Inv4Pi = 0.07957747154594766788;
+__device__ const float PiOver2 = 1.57079632679489661923;
+__device__ const float PiOver4 = 0.78539816339744830961;
+__device__ const float Sqrt2 = 1.41421356237309504880;
+__device__ const float RayEpsilon = 0.001f;
+__device__ const float FloatEpsilon = 0.000002f;
+
+__host__ __device__ inline float AbsDot(glm::vec3 a, glm::vec3 b) {
+ return glm::abs(glm::dot(a, b));
+}
+
+__host__ __device__ inline bool IsBlack(glm::vec3 a) {
+ return a.x == 0.0f && a.y == 0.0f && a.z == 0.0f;
+}
+
+__host__ __device__ inline bool SameHemisphere(glm::vec3 w, glm::vec3 wp, glm::vec3 n) {
+ return glm::dot(w, n) > 0 && glm::dot(wp, n) > 0;
+}
+
+__host__ __device__ inline float AbsCosTheta(glm::vec3 w, glm::vec3 n) {
+ return glm::abs(glm::dot(w, n));
+}
+
+__host__ __device__ inline bool fequal(float a, float b) {
+ return std::abs(a - b) < FloatEpsilon;
+}
+
+__host__ __device__ glm::vec2 ConcentricSampleDisk(glm::vec2 u) {
+ glm::vec2 uOffset = 2.0f * u - glm::vec2(1.0f, 1.0f);
+ if (uOffset.x == 0 && uOffset.y == 0) {
+ return glm::vec2(0.0f, 0.0f);
+ }
+ float theta, r;
+ if (std::abs(uOffset.x) > std::abs(uOffset.y)) {
+ r = uOffset.x;
+ theta = PiOver4 * (uOffset.y / uOffset.x);
+ }
+ else {
+ r = uOffset.y;
+ theta = PiOver2 - PiOver4 * (uOffset.x / uOffset.y);
+ }
+ return r * glm::vec2(std::cos(theta), std::sin(theta));
+}
+
+__host__ __device__ glm::vec3 squareToSphereUniform(glm::vec2 sample) {
+ float r = glm::sqrt(sample.x);
+ float phi = sample.y * 2 * Pi;
+ float u = r * glm::cos(phi);
+ float v = r * glm::sin(phi);
+ float x = u * glm::sqrt(1 - glm::pow(r, 2.0f)) * 2;
+ float y = v * glm::sqrt(1 - glm::pow(r, 2.0f)) * 2;
+ return glm::vec3(x, y, 1 - 2 * glm::pow(r, 2.0f));
+}
\ No newline at end of file
diff --git a/src/interactions.h b/src/interactions.h
index 5ce3628..e8ae6ae 100644
--- a/src/interactions.h
+++ b/src/interactions.h
@@ -1,6 +1,9 @@
#pragma once
#include "intersections.h"
+#include "globals.h"
+
+#define fresnel 1
// CHECKITOUT
/**
@@ -41,6 +44,75 @@ glm::vec3 calculateRandomDirectionInHemisphere(
+ sin(around) * over * perpendicularDirection2;
}
+__host__ __device__ void spawnRay(PathSegment &path, glm::vec3 intersection, glm::vec3 normal, glm::vec3 wi)
+{
+ path.ray.origin = intersection + RayEpsilon * normal;
+ path.ray.direction = wi;
+}
+
+#ifdef BSDF
+//********************************* lambertian bsdf ************************************
+__host__ __device__ glm::vec3 diffuseF(const Material &mat)
+{
+ return mat.color * InvPi;
+}
+
+__host__ __device__ float diffusePDF(const glm::vec3 &wi, const glm::vec3 &wo)
+{
+ return SameHemisphere(wi, wo) ? AbsCosTheta(wi) * InvPi : 0.0f;
+}
+
+__host__ __device__ glm::vec3 diffuseSampleF(
+ glm::vec3 wo,
+ glm::vec3 *wi,
+ glm::vec3 normal,
+ float *pdf,
+ const Material &m,
+ thrust::default_random_engine &rng)
+{
+ *wi = calculateRandomDirectionInHemisphere(normal, rng);
+ if (wo.z < 0.0f) wi->z *= -1;
+ *pdf = diffusePDF(*wi, wo);
+ return diffuseF(m);
+}
+
+__host__ __device__ void diffuse(
+ PathSegment &path,
+ glm::vec3 intersection,
+ glm::vec3 normal,
+ const Material &m,
+ thrust::default_random_engine &rng)
+{
+
+ glm::vec3 wi = calculateRandomDirectionInHemisphere(normal, rng);
+ wi = glm::normalize(wi);
+ path.color *= m.color;
+
+ spawnRay(path, intersection, normal, wi);
+}
+#endif
+
+//__host__ __device__ glm::vec3 fresnelDielectric(glm::vec3 normal, glm::vec3 direction, float ei, float et) {
+// bool isEnter = glm::dot(normal, direction) > 0;
+// if (!isEnter) {
+// std::swap(ei, et);
+// normal = -normal;
+// }
+// float eta = ei / et;
+// float cosThetaI = glm::clamp(glm::dot(glm::normalize(normal), glm::normalize(direction)), -1.0f, 1.0f);
+// float sinThetaI = std::sqrt(std::max(0.0f, 1.0f - cosThetaI * cosThetaI));
+// float sinThetaT = eta * sinThetaI;
+//
+// if (sinThetaT >= 1.0f) return glm::vec3(1.0f);
+// float cosThetaT = std::sqrt(std::max(0.0f, 1.0f - sinThetaT * sinThetaT));
+// float rpar = (et * cosThetaI - ei * cosThetaT) / (et * cosThetaI + ei * cosThetaT);
+// float rper = (ei * cosThetaI - et * cosThetaT) / (ei * cosThetaI + et * cosThetaT);
+//
+// return glm::vec3((rpar * rpar + rper * rper) / 2.0f);
+//}
+
+
+
/**
* Scatter a ray with some probabilities according to the material properties.
* For example, a diffuse surface scatters in a cosine-weighted hemisphere.
@@ -68,7 +140,8 @@ glm::vec3 calculateRandomDirectionInHemisphere(
*/
__host__ __device__
void scatterRay(
- PathSegment & pathSegment,
+ int idx,
+ PathSegment & path,
glm::vec3 intersect,
glm::vec3 normal,
const Material &m,
@@ -76,4 +149,53 @@ void scatterRay(
// TODO: implement this.
// A basic implementation of pure-diffuse shading will just call the
// calculateRandomDirectionInHemisphere defined above.
+
+ thrust::uniform_real_distribution u0(0, 1);
+ float prob = u0(rng);
+
+ glm::vec3 wi;
+
+ if (prob < m.hasReflective) {
+
+ wi = glm::reflect(path.ray.direction, normal);
+ wi = glm::normalize(wi);
+ path.color *= m.specular.color;
+
+ }
+ else if(prob < m.hasReflective + m.hasRefractive)
+ {
+
+ bool isEnter = glm::dot(normal, -path.ray.direction) > 0;
+ float etaI = 1.0f, etaO = m.indexOfRefraction;
+ glm::vec3 n = normal;
+ if (!isEnter) {
+ std::swap(etaI, etaO);
+ n = -n;
+ }
+ float eta = etaI / etaO;
+
+ float R0 = (1.0f - eta) / (1.0f + eta);
+ R0 *= R0;
+ float mc = 1 - glm::dot(glm::normalize(normal), glm::normalize(-path.ray.direction));
+ float R = R0 + (1.0f - R0) * glm::pow(mc, 5);
+ if (u0(rng) < R) {
+ wi = glm::reflect(path.ray.direction, normal);
+ }
+ else {
+ wi = glm::refract(path.ray.direction, normal, eta);
+ }
+ wi = glm::normalize(wi);
+ path.color *= m.specular.color;
+ normal = -n;
+
+ }
+ else
+ {
+ wi = calculateRandomDirectionInHemisphere(normal, rng);
+ wi = glm::normalize(wi);
+ path.color *= m.color;
+ }
+
+ spawnRay(path, intersect, normal, wi);
+
}
diff --git a/src/main.cpp b/src/main.cpp
index fe8e85e..b73126f 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -26,6 +26,8 @@ int iteration;
int width;
int height;
+float runningTime = 0.0f;
+
//-------------------------------
//-------------MAIN--------------
//-------------------------------
@@ -70,6 +72,7 @@ int main(int argc, char** argv) {
init();
// GLFW main loop
+ runningTime = 0.0f;
mainLoop();
return 0;
@@ -135,11 +138,13 @@ void runCuda() {
// execute the kernel
int frame = 0;
pathtrace(pbo_dptr, frame, iteration);
+ runningTime += timer().getGpuElapsedTimeForPreviousOperation();
// unmap buffer object
cudaGLUnmapBufferObject(pbo);
} else {
saveImage();
+ std::cout << "running time: " << runningTime << std::endl;
pathtraceFree();
cudaDeviceReset();
exit(EXIT_SUCCESS);
diff --git a/src/pathtrace.cu b/src/pathtrace.cu
index c1ec122..4bb0cf1 100644
--- a/src/pathtrace.cu
+++ b/src/pathtrace.cu
@@ -9,12 +9,24 @@
#include "scene.h"
#include "glm/glm.hpp"
#include "glm/gtx/norm.hpp"
+#include "glm/gtc/matrix_inverse.hpp"
#include "utilities.h"
#include "pathtrace.h"
#include "intersections.h"
#include "interactions.h"
#define ERRORCHECK 1
+#define FAKE 0
+#define CACHE 0
+#define SORT 0
+#define ANTIALIASING 1
+#define DOF 1
+#define DIRECTLIGHT 0
+#define MOTION 0
+#define TIMER 1
+
+#define lenRadius 0.3f
+#define focalDistance 3.0f
#define FILENAME (strrchr(__FILE__, '/') ? strrchr(__FILE__, '/') + 1 : __FILE__)
#define checkCUDAError(msg) checkCUDAErrorFn(msg, FILENAME, __LINE__)
@@ -38,6 +50,13 @@ void checkCUDAErrorFn(const char *msg, const char *file, int line) {
#endif
}
+PerformanceTimer& timer()
+{
+ static PerformanceTimer timer;
+ return timer;
+}
+
+
__host__ __device__
thrust::default_random_engine makeSeededRandomEngine(int iter, int index, int depth) {
int h = utilhash((1 << 31) | (depth << 22) | iter) ^ utilhash(index);
@@ -75,6 +94,9 @@ static PathSegment * dev_paths = NULL;
static ShadeableIntersection * dev_intersections = NULL;
// TODO: static variables for device memory, any extra info you need, etc
// ...
+static ShadeableIntersection * dev_firstIntersects = NULL;
+static Light * dev_lights = NULL;
+
void pathtraceInit(Scene *scene) {
hst_scene = scene;
@@ -96,6 +118,11 @@ void pathtraceInit(Scene *scene) {
cudaMemset(dev_intersections, 0, pixelcount * sizeof(ShadeableIntersection));
// TODO: initialize any extra device memeory you need
+ cudaMalloc(&dev_firstIntersects, pixelcount * sizeof(ShadeableIntersection));
+ cudaMemset(dev_firstIntersects, 0, pixelcount * sizeof(ShadeableIntersection));
+
+ cudaMalloc(&dev_lights, scene->lights.size() * sizeof(Light));
+ cudaMemcpy(dev_lights, scene->lights.data(), scene->lights.size() * sizeof(Light), cudaMemcpyHostToDevice);
checkCUDAError("pathtraceInit");
}
@@ -106,7 +133,10 @@ void pathtraceFree() {
cudaFree(dev_geoms);
cudaFree(dev_materials);
cudaFree(dev_intersections);
+
// TODO: clean up any extra device memory you created
+ cudaFree(dev_firstIntersects);
+ cudaFree(dev_lights);
checkCUDAError("pathtraceFree");
}
@@ -129,13 +159,34 @@ __global__ void generateRayFromCamera(Camera cam, int iter, int traceDepth, Path
PathSegment & segment = pathSegments[index];
segment.ray.origin = cam.position;
- segment.color = glm::vec3(1.0f, 1.0f, 1.0f);
+ segment.color = glm::vec3(1.0f, 1.0f, 1.0f);
// TODO: implement antialiasing by jittering the ray
+ thrust::default_random_engine rng = makeSeededRandomEngine(iter, index, 0);
+ thrust::uniform_real_distribution u01(0, 1);
+#if ANTIALIASING
+ float xJitter = u01(rng);
+ float yJitter = u01(rng);
+ segment.ray.direction = glm::normalize(cam.view
+ - cam.right * cam.pixelLength.x * ((float)x + xJitter - (float)cam.resolution.x * 0.5f)
+ - cam.up * cam.pixelLength.y * ((float)y + yJitter - (float)cam.resolution.y * 0.5f)
+ );
+#else
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)
- );
+ );
+#endif
+#if DOF
+ if (lenRadius > 0) {
+ glm::vec2 pLens = (float)lenRadius * ConcentricSampleDisk(glm::vec2(u01(rng), u01(rng)));
+ float ft = glm::abs(focalDistance / segment.ray.direction.z);
+ glm::vec3 pFocus = segment.ray.origin + ft * segment.ray.direction;
+ segment.ray.origin += cam.right * pLens.x + cam.up * pLens.y;
+ segment.ray.direction = glm::normalize(pFocus - segment.ray.origin);
+ }
+#endif
+
segment.pixelIndex = index;
segment.remainingBounces = traceDepth;
@@ -208,6 +259,7 @@ __global__ void computeIntersections(
intersections[path_index].t = t_min;
intersections[path_index].materialId = geoms[hit_geom_index].materialid;
intersections[path_index].surfaceNormal = normal;
+ intersections[path_index].position = intersect_point;
}
}
}
@@ -222,12 +274,11 @@ __global__ void computeIntersections(
// Your shaders should handle that - this can allow techniques such as
// bump mapping.
__global__ void shadeFakeMaterial (
- int iter
- , int num_paths
+ int iter
+ , int num_paths
, ShadeableIntersection * shadeableIntersections
, PathSegment * pathSegments
- , Material * materials
- )
+ , Material * materials)
{
int idx = blockIdx.x * blockDim.x + threadIdx.x;
if (idx < num_paths)
@@ -264,6 +315,156 @@ __global__ void shadeFakeMaterial (
}
}
}
+//************************* Direct Lighting *******************************
+
+__device__ glm::vec3 sampleOnLight(
+ glm::vec2 xi,
+ const Geom& geom)
+{
+ glm::vec4 localP;
+ if (geom.type == CUBE) {
+ localP = glm::vec4(xi.x - 0.5f, -0.5f, xi.y - 0.5f, 1.0f);
+
+ }
+ else if (geom.type == SPHERE) {
+ localP = glm::vec4(squareToSphereUniform(xi), 1.0f);
+ }
+
+ glm::vec3 worldP = glm::vec3(geom.transform * localP);
+ return worldP;
+}
+
+
+__device__ glm::vec3 directLight(
+ int num_lights,
+ int num_geoms,
+ glm::vec3 intersect,
+ glm::vec3 normal,
+ Light* dev_lights,
+ Geom* dev_geoms,
+ Material* materials,
+ thrust::default_random_engine &rng)
+{
+ // random select one light
+ thrust::uniform_real_distribution u01(0, 1);
+ int lightIdx = glm::min((int)(u01(rng) * num_lights), num_lights - 1);
+ Light &light = dev_lights[lightIdx];
+
+ // sample on light
+ glm::vec3 sampleP = sampleOnLight(glm::vec2(u01(rng), u01(rng)), dev_geoms[light.geomId]);
+ glm::vec3 wi = glm::normalize(sampleP - intersect);
+ float distance = glm::distance(sampleP, intersect);
+ Ray ray;
+ ray.origin = intersect + wi * 0.0001f;
+ ray.direction = wi;
+
+ // test if blocked
+ for (int i = 0; i < num_geoms; i++) {
+ Geom &geom = dev_geoms[i];
+ float t = 0;
+ bool outside;
+ glm::vec3 normalt;
+ if (geom.type == CUBE)
+ {
+ t = boxIntersectionTest(geom, ray, intersect, normalt, outside);
+ }
+ else if (geom.type == SPHERE)
+ {
+ t = sphereIntersectionTest(geom, ray, intersect, normalt, outside);
+ }
+ if (t > 0.0f && materials[geom.materialid].emittance < 0.0001f && t < distance) {
+ return glm::vec3(0.0f);
+ }
+
+ }
+
+ return materials[light.materialId].color * materials[light.materialId].emittance * AbsDot(wi, normal);
+}
+//*************************************************************************
+
+//************************* Shading Kernal ********************************
+
+#pragma region BSDFShading
+
+__global__ void kernShading(
+ int iter,
+ int depth,
+ int num_paths,
+ ShadeableIntersection* shadeableIntersections,
+ PathSegment* pathsegments,
+ Material* materials
+#if DIRECTLIGHT
+ ,
+ int num_lights,
+ int num_geoms,
+ Light* lights,
+ Geom* geoms
+#endif
+)
+{
+ int idx = blockDim.x * blockIdx.x + threadIdx.x;
+ if (idx < num_paths) {
+
+ ShadeableIntersection intersection = shadeableIntersections[idx];
+
+ if (intersection.t > 0.0f) {
+ thrust::default_random_engine rng = makeSeededRandomEngine(iter, idx, depth);
+ //thrust::uniform_real_distribution u01(0, 1);
+
+ Material material = materials[intersection.materialId];
+ glm::vec3 materialColor = material.color;
+
+ if (material.emittance > 0.0f) {
+ pathsegments[idx].color *= (materialColor * material.emittance);
+ pathsegments[idx].remainingBounces = 0;
+ return;
+ }
+ else
+ {
+ scatterRay(idx, pathsegments[idx], intersection.position, intersection.surfaceNormal, material, rng);
+ pathsegments[idx].remainingBounces--;
+ if (pathsegments[idx].remainingBounces == 0)
+ {
+#if DIRECTLIGHT
+ pathsegments[idx].color *= directLight(num_lights, num_geoms, intersection.position, intersection.surfaceNormal, lights, geoms, materials, rng);
+#else
+ pathsegments[idx].color = glm::vec3(0.0f);
+#endif
+ }
+ }
+ }
+ else {
+ pathsegments[idx].color = glm::vec3(0.0f);
+ pathsegments[idx].remainingBounces = 0;
+ }
+ }
+}
+
+#pragma endregion BSDFShading
+
+//*************************************************************************
+
+//************************* Stream Compact ********************************
+struct isRemainBouncing
+{
+ __host__ __device__ bool operator()(const PathSegment& path)
+ {
+ return path.remainingBounces == 0;
+ }
+};
+//*************************************************************************
+
+//***************************** Sort **************************************
+struct materialSort
+{
+ __host__ __device__ bool operator()(const ShadeableIntersection &a, const ShadeableIntersection &b)
+ {
+ return a.materialId < b.materialId;
+ }
+};
+//*************************************************************************
+
+
// Add the current iteration's output to the overall image
__global__ void finalGather(int nPaths, glm::vec3 * image, PathSegment * iterationPaths)
@@ -273,7 +474,9 @@ __global__ void finalGather(int nPaths, glm::vec3 * image, PathSegment * iterati
if (index < nPaths)
{
PathSegment iterationPath = iterationPaths[index];
- image[iterationPath.pixelIndex] += iterationPath.color;
+ if (iterationPath.remainingBounces == 0) {
+ image[iterationPath.pixelIndex] += iterationPath.color;
+ }
}
}
@@ -325,60 +528,150 @@ void pathtrace(uchar4 *pbo, int frame, int iter) {
// for you.
// TODO: perform one iteration of path tracing
+#if TIMER
+ timer().startGpuTimer();
+#endif
- generateRayFromCamera <<>>(cam, iter, traceDepth, dev_paths);
+#if MOTION
+ float moveSpeed = 1.0f / (float)hst_scene->state.iterations;
+ for (int i = 0; i < hst_scene->geoms.size(); i++)
+ {
+ if (hst_scene->geoms[i].motion)
+ {
+ hst_scene->geoms[i].translation += glm::vec3(0.0f, 3.0f, 0.0f) * moveSpeed;
+ hst_scene->geoms[i].transform = utilityCore::buildTransformationMatrix(hst_scene->geoms[i].translation, hst_scene->geoms[i].rotation, hst_scene->geoms[i].scale);
+ hst_scene->geoms[i].inverseTransform = glm::inverse(hst_scene->geoms[i].transform);
+ hst_scene->geoms[i].invTranspose = glm::inverseTranspose(hst_scene->geoms[i].transform);
+ }
+ }
+ cudaMemcpy(dev_geoms, &(hst_scene->geoms)[0], hst_scene->geoms.size() * sizeof(Geom), cudaMemcpyHostToDevice);
+ checkCUDAError("motion blur error");
+#endif
+
+ generateRayFromCamera << > >(cam, iter, traceDepth, dev_paths);
checkCUDAError("generate camera ray");
+
int depth = 0;
PathSegment* dev_path_end = dev_paths + pixelcount;
int num_paths = dev_path_end - dev_paths;
+ dim3 numBlocksPixels = (pixelcount + blockSize1d - 1) / blockSize1d;
+
// --- PathSegment Tracing Stage ---
// Shoot ray into scene, bounce between objects, push shading chunks
- bool iterationComplete = false;
+ bool iterationComplete = false;
while (!iterationComplete) {
- // clean shading chunks
- cudaMemset(dev_intersections, 0, pixelcount * sizeof(ShadeableIntersection));
-
- // tracing
- dim3 numblocksPathSegmentTracing = (num_paths + blockSize1d - 1) / blockSize1d;
- computeIntersections <<>> (
- depth
- , num_paths
- , dev_paths
- , dev_geoms
- , hst_scene->geoms.size()
- , dev_intersections
- );
- checkCUDAError("trace one bounce");
- cudaDeviceSynchronize();
- depth++;
-
-
- // TODO:
- // --- Shading Stage ---
- // Shade path segments based on intersections and generate new rays by
- // evaluating the BSDF.
- // Start off with just a big kernel that handles all the different
- // materials you have in the scenefile.
- // TODO: compare between directly shading the path segments and shading
- // path segments that have been reshuffled to be contiguous in memory.
-
- shadeFakeMaterial<<>> (
- iter,
- num_paths,
- dev_intersections,
- dev_paths,
- dev_materials
- );
- iterationComplete = true; // TODO: should be based off stream compaction results.
+ // clean shading chunks
+ cudaMemset(dev_intersections, 0, pixelcount * sizeof(ShadeableIntersection));
+
+ // tracing
+ dim3 numblocksPathSegmentTracing = (num_paths + blockSize1d - 1) / blockSize1d;
+#if CACHE
+ if (iter == 1 && depth == 0) {
+ computeIntersections << > > (
+ depth
+ , num_paths
+ , dev_paths
+ , dev_geoms
+ , hst_scene->geoms.size()
+ , dev_intersections
+ );
+ checkCUDAError("trace one bounce");
+ cudaDeviceSynchronize();
+ cudaMemcpy(dev_firstIntersects, dev_intersections, pixelcount * sizeof(ShadeableIntersection), cudaMemcpyDeviceToDevice);
+ }
+ else if(iter != 1 && depth == 0){
+ cudaMemcpy(dev_intersections, dev_firstIntersects, pixelcount * sizeof(ShadeableIntersection), cudaMemcpyDeviceToDevice);
+ }
+ else {
+ computeIntersections << > > (
+ depth
+ , num_paths
+ , dev_paths
+ , dev_geoms
+ , hst_scene->geoms.size()
+ , dev_intersections
+ );
+ checkCUDAError("trace one bounce");
+ cudaDeviceSynchronize();
+ }
+#else
+ computeIntersections << > > (
+ depth
+ , num_paths
+ , dev_paths
+ , dev_geoms
+ , hst_scene->geoms.size()
+ , dev_intersections
+ );
+ checkCUDAError("trace one bounce");
+ cudaDeviceSynchronize();
+#endif
+
+ depth++;
+
+#if SORT
+ thrust::sort_by_key(thrust::device, dev_intersections, dev_intersections + num_paths, dev_paths, materialSort());
+#endif
+
+ // TODO:
+ // --- Shading Stage ---
+ // Shade path segments based on intersections and generate new rays by
+ // evaluating the BSDF.
+ // Start off with just a big kernel that handles all the different
+ // materials you have in the scenefile.
+ // TODO: compare between directly shading the path segments and shading
+ // path segments that have been reshuffled to be contiguous in memory.
+
+#if FAKE
+ shadeFakeMaterial<<>> (
+ iter,
+ num_paths,
+ dev_intersections,
+ dev_paths,
+ dev_materials);
+ iterationComplete = true; // TODO: should be based off stream compaction results.
+#else
+ kernShading << > > (
+ iter,
+ depth,
+ num_paths,
+ dev_intersections,
+ dev_paths,
+ dev_materials
+#if DIRECTLIGHT
+ ,
+ hst_scene->lights.size(),
+ hst_scene->geoms.size(),
+ dev_lights,
+ dev_geoms
+#endif
+ );
+ checkCUDAError("kern Shading error");
+
+ finalGather << > >(num_paths, dev_image, dev_paths);
+
+ if (depth == traceDepth) {
+ iterationComplete = true;
+ }
+
+ dev_path_end = thrust::remove_if(thrust::device, dev_paths, dev_paths + num_paths, isRemainBouncing());
+ num_paths = dev_path_end - dev_paths;
+ if (num_paths == 0) {
+ iterationComplete = true;
+ }
+
+#endif
}
- // Assemble this iteration and apply it to the image
- dim3 numBlocksPixels = (pixelcount + blockSize1d - 1) / blockSize1d;
- finalGather<<>>(num_paths, dev_image, dev_paths);
+#if TIMER
+ timer().endGpuTimer();
+#endif
+
+
///////////////////////////////////////////////////////////////////////////
@@ -390,4 +683,5 @@ void pathtrace(uchar4 *pbo, int frame, int iter) {
pixelcount * sizeof(glm::vec3), cudaMemcpyDeviceToHost);
checkCUDAError("pathtrace");
+
}
diff --git a/src/pathtrace.h b/src/pathtrace.h
index 1241227..f3734ab 100644
--- a/src/pathtrace.h
+++ b/src/pathtrace.h
@@ -3,6 +3,94 @@
#include
#include "scene.h"
+class PerformanceTimer
+{
+public:
+ PerformanceTimer()
+ {
+ cudaEventCreate(&event_start);
+ cudaEventCreate(&event_end);
+ }
+
+ ~PerformanceTimer()
+ {
+ cudaEventDestroy(event_start);
+ cudaEventDestroy(event_end);
+ }
+
+ void startCpuTimer()
+ {
+ if (cpu_timer_started) { throw std::runtime_error("CPU timer already started"); }
+ cpu_timer_started = true;
+
+ time_start_cpu = std::chrono::high_resolution_clock::now();
+ }
+
+ void endCpuTimer()
+ {
+ time_end_cpu = std::chrono::high_resolution_clock::now();
+
+ if (!cpu_timer_started) { throw std::runtime_error("CPU timer not started"); }
+
+ std::chrono::duration duro = time_end_cpu - time_start_cpu;
+ prev_elapsed_time_cpu_milliseconds =
+ static_cast(duro.count());
+
+ cpu_timer_started = false;
+ }
+
+ void startGpuTimer()
+ {
+ if (gpu_timer_started) { throw std::runtime_error("GPU timer already started"); }
+ gpu_timer_started = true;
+
+ cudaEventRecord(event_start);
+ }
+
+ void endGpuTimer()
+ {
+ cudaEventRecord(event_end);
+ cudaEventSynchronize(event_end);
+
+ if (!gpu_timer_started) { throw std::runtime_error("GPU timer not started"); }
+
+ cudaEventElapsedTime(&prev_elapsed_time_gpu_milliseconds, event_start, event_end);
+ gpu_timer_started = false;
+ }
+
+ float getCpuElapsedTimeForPreviousOperation() //noexcept //(damn I need VS 2015
+ {
+ return prev_elapsed_time_cpu_milliseconds;
+ }
+
+ float getGpuElapsedTimeForPreviousOperation() //noexcept
+ {
+ return prev_elapsed_time_gpu_milliseconds;
+ }
+
+ // remove copy and move functions
+ PerformanceTimer(const PerformanceTimer&) = delete;
+ PerformanceTimer(PerformanceTimer&&) = delete;
+ PerformanceTimer& operator=(const PerformanceTimer&) = delete;
+ PerformanceTimer& operator=(PerformanceTimer&&) = delete;
+
+private:
+ cudaEvent_t event_start = nullptr;
+ cudaEvent_t event_end = nullptr;
+
+ using time_point_t = std::chrono::high_resolution_clock::time_point;
+ time_point_t time_start_cpu;
+ time_point_t time_end_cpu;
+
+ bool cpu_timer_started = false;
+ bool gpu_timer_started = false;
+
+ float prev_elapsed_time_cpu_milliseconds = 0.f;
+ float prev_elapsed_time_gpu_milliseconds = 0.f;
+};
+
void pathtraceInit(Scene *scene);
void pathtraceFree();
void pathtrace(uchar4 *pbo, int frame, int iteration);
+
+extern PerformanceTimer &timer();
diff --git a/src/scene.cpp b/src/scene.cpp
index cbae043..76d4f16 100644
--- a/src/scene.cpp
+++ b/src/scene.cpp
@@ -30,6 +30,7 @@ Scene::Scene(string filename) {
}
}
}
+ loadLight();
}
int Scene::loadGeom(string objectid) {
@@ -74,7 +75,10 @@ int Scene::loadGeom(string objectid) {
newGeom.rotation = glm::vec3(atof(tokens[1].c_str()), atof(tokens[2].c_str()), atof(tokens[3].c_str()));
} else if (strcmp(tokens[0].c_str(), "SCALE") == 0) {
newGeom.scale = glm::vec3(atof(tokens[1].c_str()), atof(tokens[2].c_str()), atof(tokens[3].c_str()));
- }
+ }
+ else if (strcmp(tokens[0].c_str(), "MOTION") == 0) {
+ newGeom.motion = atoi(tokens[1].c_str()) == 1;
+ }
utilityCore::safeGetline(fp_in, line);
}
@@ -186,3 +190,16 @@ int Scene::loadMaterial(string materialid) {
return 1;
}
}
+
+
+void Scene::loadLight()
+{
+ for (int i = 0; i < geoms.size(); i++) {
+ if (materials[geoms[i].materialid].emittance > 0) {
+ Light light;
+ light.geomId = i;
+ light.materialId = geoms[i].materialid;
+ lights.push_back(light);
+ }
+ }
+}
diff --git a/src/scene.h b/src/scene.h
index f29a917..67512d7 100644
--- a/src/scene.h
+++ b/src/scene.h
@@ -16,11 +16,13 @@ class Scene {
int loadMaterial(string materialid);
int loadGeom(string objectid);
int loadCamera();
+ void loadLight();
public:
Scene(string filename);
~Scene();
std::vector geoms;
std::vector materials;
+ std::vector lights;
RenderState state;
};
diff --git a/src/sceneStructs.h b/src/sceneStructs.h
index b38b820..1048ccf 100644
--- a/src/sceneStructs.h
+++ b/src/sceneStructs.h
@@ -26,6 +26,13 @@ struct Geom {
glm::mat4 transform;
glm::mat4 inverseTransform;
glm::mat4 invTranspose;
+
+ bool motion;
+};
+
+struct Light {
+ int materialId;
+ int geomId;
};
struct Material {
@@ -71,6 +78,7 @@ struct PathSegment {
// 2) BSDF evaluation: generate a new ray
struct ShadeableIntersection {
float t;
+ glm::vec3 position;
glm::vec3 surfaceNormal;
int materialId;
};
diff --git a/src/utilities.h b/src/utilities.h
index abb4f27..461ccbb 100644
--- a/src/utilities.h
+++ b/src/utilities.h
@@ -9,6 +9,14 @@
#include
#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+
#define PI 3.1415926535897932384626422832795028841971f
#define TWO_PI 6.2831853071795864769252867665590057683943f
#define SQRT_OF_ONE_THIRD 0.5773502691896257645091487805019574556476f