diff --git a/.gitattributes b/.gitattributes new file mode 100644 index 0000000..d41a940 --- /dev/null +++ b/.gitattributes @@ -0,0 +1 @@ +*.gif filter=lfs diff=lfs merge=lfs -text diff --git a/CMakeLists.txt b/CMakeLists.txt index d3d976c..7eb294f 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -2,6 +2,16 @@ cmake_minimum_required(VERSION 3.0) project(cis565_path_tracer) +# Crucial magic for CUDA linking +find_package(Threads REQUIRED) +find_package(CUDA 8.0 REQUIRED) + +## Because screw life +## set(CUDACC_DEFINE -D__CUDACC__) + +## Because fuck nvcc +## add_definitions(-D__CUDACC__) + set(CMAKE_MODULE_PATH "${CMAKE_SOURCE_DIR}/cmake" ${CMAKE_MODULE_PATH}) # Set up include and lib paths @@ -60,10 +70,6 @@ if (WIN32) list(APPEND CORELIBS legacy_stdio_definitions.lib) endif() -# Crucial magic for CUDA linking -find_package(Threads REQUIRED) -find_package(CUDA 8.0 REQUIRED) - set(CUDA_ATTACH_VS_BUILD_RULE_TO_CUDA_FILE ON) set(CUDA_SEPARABLE_COMPILATION ON) @@ -71,6 +77,9 @@ if(${CMAKE_SYSTEM_NAME} MATCHES "Darwin") set(CUDA_PROPAGATE_HOST_FLAGS OFF) endif() +## Because screw life +set(CUDACC_DEFINE -D__CUDACC__) + include_directories(.) #add_subdirectory(stream_compaction) # TODO: uncomment if using your stream compaction add_subdirectory(src) diff --git a/README.md b/README.md index 110697c..a56c88f 100644 --- a/README.md +++ b/README.md @@ -1,13 +1,153 @@ -CUDA Path Tracer -================ +# **University of Pennsylvania, CIS 565: GPU Programming and Architecture, Project 3:** -**University of Pennsylvania, CIS 565: GPU Programming and Architecture, Project 3** +# **CUDA Path Real Time Path Tracer** -* (TODO) YOUR NAME HERE -* Tested on: (TODO) Windows 22, i7-2222 @ 2.22GHz 22GB, GTX 222 222MB (Moore 2222 Lab) -### (TODO: Your README) -*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. + + +Tested on: Windows 10, Intel Core i7-7700HQ CPU @ 2.80 GHz, 8GB RAM, NVidia GeForce GTX 1050 + + ![Built](https://img.shields.io/appveyor/ci/gruntjs/grunt.svg) ![Issues](https://img.shields.io/github/issues-raw/badges/shields/website.svg) ![CUDA 8.0](https://img.shields.io/badge/CUDA-8.0-green.svg?style=flat) ![Platform](https://img.shields.io/badge/platform-Desktop-bcbcbc.svg) ![Developer](https://img.shields.io/badge/Developer-Youssef%20Victor-0f97ff.svg?style=flat) + + + + +- [Features](#features) + + + +- [In-Depth](#indepth) + + + +- [Blooper](#blooper) + + + + +____________________________________________________ + + + +The goal of this project was to run an algorithm that clears out all zeros from an array on the GPU using CUDA. This parallel reduction is done using the scan algorithm that computes the exclusive prefix sum. I also implemented a parallelized Radix Sort using the exclusive prefix sum algorithm developed. + + + +### Things Done + +#### Core Features + + - [x] Shading for Ideal Diffuse Surfaces + - [x] Shading for Perfectly Specular Surfaces + - [x] Early Ray Termination Using Stream Compaction + - [x] Caching of First Bounces + + #### Spicy Features + - [x] Shading for Transmissive Surfaces (maybe) + - [x] Work-efficient Stream Compaction Using Shared Memory Across Multiple Blocks + - [x] Environment Maps + - [x] Direct Lighting + - [x] Multiple Importance Sampling + - [x] Stochastic Antialiasing + - [x] Arbitrary Mesh Loading + - [ ] Kd-Tree & Stackless Kd-Tree Traversal on the GPU + - [ ] [Physically Accurate Lens Flares](https://placeholderart.wordpress.com/2015/01/19/implementation-notes-physically-based-lens-flares/) + - [ ] [Specular BRDF with Microstructure BRDF](https://people.eecs.berkeley.edu/~lingqi/publications/paper_glints2.pdf) + - [ ] Texture Mapping (maybe) + - [ ] Depth of Field/Cooler Lens Effects such as bloom + +![Cornell Box](/img/cornellOBJ.png) + +Multiple Importance Sampling, 7500 Samples, 1080x1080 px, OBJ dodecahedron mesh. + + + +### In-Depth: + +#### Transmissive Surfaces + +This uses basic `glm::refract` so it should work. For some reason it doesn't. I should change the implementation of it, but for now, maybe I can get half-credit points. To try out a scene with it simply set a material's bsdf to 2. +`TODO: Show picture` + +#### Work Efficient Scan Using Shared Memory Across Multiple Blocks + +This took me like 5 days to do, I ended up simply using thrust's partition, but well, that's programming for you. The benefit to my implementation is that it works all the way up to `2^24` as opposed to some others' implementation which fails before that. +`TODO: Show performance comparison` + +#### Environment Maps: + +![environment-map](/img/envMap2.png) + +An environment map with a purely specular sphere in the middle. + +![diffuse-environment-map](/img/envMap.png) + +An environment map with a purely diffuse sphere in the middle. + +In the sample images above, you can clearly see the effect the environment map has with specular surfaces. With the purely specular surface it sort of looks like there is some bloom effeect showing as well. I think that is really cool. + +The second image with a diffuse surface is supposed to be a 98% white sphere. But as you see, as in the real world, the environment also reflects light onto the object. The result is the beautifully shaded sphere you see. + +In the following images you will see how the environment maps affect the lighting of a cornell box scene with a specular sphere in the middle: + +![cornellBoxNoEnvMap](/img/cornellSpec.png) +Normal cornell box scene with no environment map. + +![cornellBoxNoEnvMap](/img/cornellEnvSpec.png) +Normal cornell box scene with environment map added. + +The environment map adds better lighting to the scene. + + +#### Direct Lighting + +This took a while to get completely right, it uses a light-based sample to shade the entire scene. As such, there are no reflections on objects. The result however is a much more nicely converged scene. + +`TODO: Show Pic` + +#### Multiple Importance Sampling + +When you combine the light based sampling and the bsdf-based sampling, you get what you saw in the representative image at the start: a completely converged. Here are some more pictures: + +![mis-sample1](/img/cornellTwo.2017-10-01_21-11-00z.5000samp.png) + +`TODO: Add More Images` + +#### Stochastic Anti-Aliasing: + +This Feature was fairly simple to implement yet took a lot to perfect. This feature does not work with first-bounce caching unless you jitter the ray direction after you generate the ray, which is possible and might be a `TODO`, but given the lack of time, it will probably remain an idea. Here is what happens when I try anti-aliasing with caching turned on: + +![bad-aliasing](/img/cornell.2017-10-01_18-30-30z.5000samp.png) + +With my cached first bounce, there are jagged edges surrounding everything as all the rays have been jittered, but they do not change across iterations, so the random jitterness stays and remains prominent. As such the picture looks very weird. Here's what it looks like without anti-aliasing at all: + +![no-aliasing](/img/cornellTwo.2017-10-01_09-02-39z.5000samp.png) + +The image here has very rough aliased edges. Here is what it looks like with everything fixed: + +![aliasing](/img/cornellTwo.2017-10-01_21-11-00z.5000samp.png) + +Now the image is much smoother. + +#### Arbitrary Mesh Loading + +Took me forever to get this to work, triangle intersection and intersection testing on the GPU caused a lot of problems. Eventually they were all resolved, but still. Here is the beautiful representative image once again with a stretched out cyan dodecahedron model: + +![Cornell Box](/img/cornellOBJ.png) + + +### Bloopers / Lessons Learned + +This is a blooper I got while trying to get MIS to work + +![bnw](/img/bloopers/bnw.gif) + +Another blooper I got where I was sampling the light weirdly + +![the-v](/img/bloopers/v.gif) + +My favorite Blooper, which also shows Direct Ligting Rays being sampled: + +![starry stuff](/img/bloopers/starry.gif) diff --git a/img/MIS-1000samps.PNG b/img/MIS-1000samps.PNG new file mode 100644 index 0000000..6693552 Binary files /dev/null and b/img/MIS-1000samps.PNG differ diff --git a/img/bloopers/bnw.gif b/img/bloopers/bnw.gif new file mode 100644 index 0000000..3fbe9c5 Binary files /dev/null and b/img/bloopers/bnw.gif differ diff --git a/img/bloopers/starry.gif b/img/bloopers/starry.gif new file mode 100644 index 0000000..046a9f6 Binary files /dev/null and b/img/bloopers/starry.gif differ diff --git a/img/bloopers/v.gif b/img/bloopers/v.gif new file mode 100644 index 0000000..2224621 Binary files /dev/null and b/img/bloopers/v.gif differ diff --git a/img/cornell.2017-10-01_18-30-30z.5000samp.png b/img/cornell.2017-10-01_18-30-30z.5000samp.png new file mode 100644 index 0000000..d70ce26 Binary files /dev/null and b/img/cornell.2017-10-01_18-30-30z.5000samp.png differ diff --git a/img/cornellEnvSpec.png b/img/cornellEnvSpec.png new file mode 100644 index 0000000..fd396ee Binary files /dev/null and b/img/cornellEnvSpec.png differ diff --git a/img/cornellOBJ.png b/img/cornellOBJ.png new file mode 100644 index 0000000..305d292 Binary files /dev/null and b/img/cornellOBJ.png differ diff --git a/img/cornellSpec.png b/img/cornellSpec.png new file mode 100644 index 0000000..e953b70 Binary files /dev/null and b/img/cornellSpec.png differ diff --git a/img/cornellTwo.2017-10-01_00-59-01z.5000samp.png b/img/cornellTwo.2017-10-01_00-59-01z.5000samp.png new file mode 100644 index 0000000..8b55eb4 Binary files /dev/null and b/img/cornellTwo.2017-10-01_00-59-01z.5000samp.png differ diff --git a/img/cornellTwo.2017-10-01_02-07-17z.5000samp.png b/img/cornellTwo.2017-10-01_02-07-17z.5000samp.png new file mode 100644 index 0000000..503fd73 Binary files /dev/null and b/img/cornellTwo.2017-10-01_02-07-17z.5000samp.png differ diff --git a/img/cornellTwo.2017-10-01_04-03-53z.5000samp.png b/img/cornellTwo.2017-10-01_04-03-53z.5000samp.png new file mode 100644 index 0000000..0475f8e Binary files /dev/null and b/img/cornellTwo.2017-10-01_04-03-53z.5000samp.png differ diff --git a/img/cornellTwo.2017-10-01_04-16-37z.0samp.png b/img/cornellTwo.2017-10-01_04-16-37z.0samp.png new file mode 100644 index 0000000..d2ec3f1 Binary files /dev/null and b/img/cornellTwo.2017-10-01_04-16-37z.0samp.png differ diff --git a/img/cornellTwo.2017-10-01_08-41-23z.5000samp.png b/img/cornellTwo.2017-10-01_08-41-23z.5000samp.png new file mode 100644 index 0000000..a550340 Binary files /dev/null and b/img/cornellTwo.2017-10-01_08-41-23z.5000samp.png differ diff --git a/img/cornellTwo.2017-10-01_09-02-39z.26samp.png b/img/cornellTwo.2017-10-01_09-02-39z.26samp.png new file mode 100644 index 0000000..6a137f3 Binary files /dev/null and b/img/cornellTwo.2017-10-01_09-02-39z.26samp.png differ diff --git a/img/cornellTwo.2017-10-01_09-02-39z.5000samp.png b/img/cornellTwo.2017-10-01_09-02-39z.5000samp.png new file mode 100644 index 0000000..e1ed79e Binary files /dev/null and b/img/cornellTwo.2017-10-01_09-02-39z.5000samp.png differ diff --git a/img/cornellTwo.2017-10-01_19-53-41z.5000samp.png b/img/cornellTwo.2017-10-01_19-53-41z.5000samp.png new file mode 100644 index 0000000..f2b82c3 Binary files /dev/null and b/img/cornellTwo.2017-10-01_19-53-41z.5000samp.png differ diff --git a/img/cornellTwo.2017-10-01_21-11-00z.5000samp.png b/img/cornellTwo.2017-10-01_21-11-00z.5000samp.png new file mode 100644 index 0000000..9cf3b3b Binary files /dev/null and b/img/cornellTwo.2017-10-01_21-11-00z.5000samp.png differ diff --git a/img/envMap.PNG b/img/envMap.PNG new file mode 100644 index 0000000..d5dd9cd Binary files /dev/null and b/img/envMap.PNG differ diff --git a/img/envMap2.PNG b/img/envMap2.PNG new file mode 100644 index 0000000..9990f2f Binary files /dev/null and b/img/envMap2.PNG differ diff --git a/img/environmentTest.2017-10-02_03-40-14z.5000samp.png b/img/environmentTest.2017-10-02_03-40-14z.5000samp.png new file mode 100644 index 0000000..48b08a0 Binary files /dev/null and b/img/environmentTest.2017-10-02_03-40-14z.5000samp.png differ diff --git a/img/environmentTest.2017-10-02_04-37-11z.5000samp.png b/img/environmentTest.2017-10-02_04-37-11z.5000samp.png new file mode 100644 index 0000000..4937d86 Binary files /dev/null and b/img/environmentTest.2017-10-02_04-37-11z.5000samp.png differ diff --git a/scenes/cornell.txt b/scenes/cornell.txt index 83ff820..21069bf 100644 --- a/scenes/cornell.txt +++ b/scenes/cornell.txt @@ -6,7 +6,8 @@ SPECRGB 0 0 0 REFL 0 REFR 0 REFRIOR 0 -EMITTANCE 5 +EMITTANCE 3 +BSDF -1 // Diffuse white MATERIAL 1 @@ -17,6 +18,7 @@ REFL 0 REFR 0 REFRIOR 0 EMITTANCE 0 +BSDF 0 // Diffuse red MATERIAL 2 @@ -27,6 +29,7 @@ REFL 0 REFR 0 REFRIOR 0 EMITTANCE 0 +BSDF 0 // Diffuse green MATERIAL 3 @@ -37,6 +40,7 @@ REFL 0 REFR 0 REFRIOR 0 EMITTANCE 0 +BSDF 0 // Specular white MATERIAL 4 @@ -47,6 +51,7 @@ REFL 1 REFR 0 REFRIOR 0 EMITTANCE 0 +BSDF 0 // Camera CAMERA diff --git a/scenes/cornellDraggo.txt b/scenes/cornellDraggo.txt new file mode 100644 index 0000000..8357f77 --- /dev/null +++ b/scenes/cornellDraggo.txt @@ -0,0 +1,150 @@ +// Emissive material (light) +MATERIAL 0 +RGB 0.8 1 1 +SPECEX 0 +SPECRGB 0 0 0 +REFL 0 +REFR 0 +REFRIOR 0 +EMITTANCE 10 +BSDF -1 + +// Emissive material (light) +MATERIAL 1 +RGB 1 1 0.8 +SPECEX 0 +SPECRGB 0 0 0 +REFL 0 +REFR 0 +REFRIOR 0 +EMITTANCE 10 +BSDF -1 + +// Diffuse white +MATERIAL 2 +RGB .98 .98 .98 +SPECEX 0 +SPECRGB 0 0 0 +REFL 0 +REFR 0 +REFRIOR 0 +EMITTANCE 0 +BSDF 0 + +// Diffuse red +MATERIAL 3 +RGB .85 .35 .35 +SPECEX 0 +SPECRGB 0 0 0 +REFL 0 +REFR 0 +REFRIOR 0 +EMITTANCE 0 +BSDF 0 + +// Diffuse green +MATERIAL 4 +RGB .35 .85 .35 +SPECEX 0 +SPECRGB 0 0 0 +REFL 0 +REFR 0 +REFRIOR 0 +EMITTANCE 0 +BSDF 0 + +// Specular white +MATERIAL 5 +RGB .98 .98 .98 +SPECEX 0 +SPECRGB .98 .98 .98 +REFL 1 +REFR 0 +REFRIOR 0 +EMITTANCE 0 +BSDF 1 + +// Camera +CAMERA +RES 800 800 +FOVY 45 +ITERATIONS 5000 +DEPTH 8 +FILE cornellTwo +EYE 0.0 5 10.5 +LOOKAT 0 5 0 +UP 0 1 0 + + +// Ceiling light +OBJECT 0 +plane +material 0 +TRANS 3 9.99 0 +ROTAT 90 0 0 +SCALE 3 3 1 + +// Ceiling light 2 +OBJECT 1 +plane +material 1 +TRANS -3 9.99 0 +ROTAT 90 0 0 +SCALE 3 3 1 + +// Floor +OBJECT 2 +cube +material 2 +TRANS 0 0 0 +ROTAT 0 0 0 +SCALE 10 .01 10 + +// Ceiling +OBJECT 3 +cube +material 2 +TRANS 0 10 0 +ROTAT 0 0 90 +SCALE .01 10 10 + +// Back wall +OBJECT 4 +cube +material 2 +TRANS 0 5 -5 +ROTAT 0 90 0 +SCALE .01 10 10 + +// Left wall +OBJECT 5 +cube +material 3 +TRANS -5 5 0 +ROTAT 0 0 0 +SCALE .01 10 10 + +// Right wall +OBJECT 6 +cube +material 4 +TRANS 5 5 0 +ROTAT 0 0 0 +SCALE .01 10 10 + +// Sphere +OBJECT 7 +sphere +material 2 +TRANS -1 4 -4 +ROTAT 0 0 0 +SCALE 0.1 0.1 0.1 + +// Cube Mesh +OBJECT 8 +mesh +../scenes/dragon.obj +material 4 +TRANS 0 2 -4 +ROTAT 0 0 0 +SCALE 2 2 2 \ No newline at end of file diff --git a/scenes/cornellReflect.txt b/scenes/cornellReflect.txt new file mode 100644 index 0000000..5e096cd --- /dev/null +++ b/scenes/cornellReflect.txt @@ -0,0 +1,133 @@ +// Emissive material (light) +MATERIAL 0 +RGB 1 1 1 +SPECEX 0 +SPECRGB 0 0 0 +REFL 0 +REFR 0 +REFRIOR 0 +EMITTANCE 5 +BSDF -1 + +// Diffuse white +MATERIAL 1 +RGB .98 .98 .98 +SPECEX 0 +SPECRGB 0 0 0 +REFL 0 +REFR 0 +REFRIOR 0 +EMITTANCE 0 +BSDF 0 + +// Diffuse red +MATERIAL 2 +RGB .85 .35 .35 +SPECEX 0 +SPECRGB 0 0 0 +REFL 0 +REFR 0 +REFRIOR 0 +EMITTANCE 0 +BSDF 0 + +// Diffuse green +MATERIAL 3 +RGB .35 .85 .35 +SPECEX 0 +SPECRGB 0 0 0 +REFL 0 +REFR 0 +REFRIOR 0 +EMITTANCE 0 +BSDF 0 + +// Specular white +MATERIAL 4 +RGB .98 .98 .98 +SPECEX 0 +SPECRGB .98 .98 .98 +REFL 1 +REFR 0 +REFRIOR 0 +EMITTANCE 0 +BSDF 0 + +// Specular white +MATERIAL 5 +RGB .98 .98 .98 +SPECEX 0 +SPECRGB .98 .98 .98 +REFL 1 +REFR 0 +REFRIOR 0 +EMITTANCE 0 +BSDF 1 + +// 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 + +// Floor +OBJECT 1 +cube +material 1 +TRANS 0 0 0 +ROTAT 0 0 0 +SCALE 10 .01 10 + +// Ceiling +OBJECT 2 +cube +material 1 +TRANS 0 10 0 +ROTAT 0 0 90 +SCALE .01 10 10 + +// Back wall +OBJECT 3 +cube +material 1 +TRANS 0 5 -5 +ROTAT 0 90 0 +SCALE .01 10 10 + +// Left wall +OBJECT 4 +cube +material 2 +TRANS -5 5 0 +ROTAT 0 0 0 +SCALE .01 10 10 + +// Right wall +OBJECT 5 +cube +material 3 +TRANS 5 5 0 +ROTAT 0 0 0 +SCALE .01 10 10 + +// Sphere +OBJECT 6 +sphere +material 5 +TRANS -1 4 -1 +ROTAT 0 0 0 +SCALE 3 3 3 diff --git a/scenes/cornellRefract.txt b/scenes/cornellRefract.txt new file mode 100644 index 0000000..c267ae3 --- /dev/null +++ b/scenes/cornellRefract.txt @@ -0,0 +1,144 @@ +// Emissive material (light) +MATERIAL 0 +RGB 1 1 1 +SPECEX 0 +SPECRGB 0 0 0 +REFL 0 +REFR 0 +REFRIOR 0 +EMITTANCE 5 +BSDF -1 + +// Diffuse white +MATERIAL 1 +RGB .98 .98 .98 +SPECEX 0 +SPECRGB 0 0 0 +REFL 0 +REFR 0 +REFRIOR 0 +EMITTANCE 0 +BSDF 0 + +// Diffuse red +MATERIAL 2 +RGB .85 .35 .35 +SPECEX 0 +SPECRGB 0 0 0 +REFL 0 +REFR 0 +REFRIOR 0 +EMITTANCE 0 +BSDF 0 + +// Diffuse green +MATERIAL 3 +RGB .35 .85 .35 +SPECEX 0 +SPECRGB 0 0 0 +REFL 0 +REFR 0 +REFRIOR 0 +EMITTANCE 0 +BSDF 0 + +// Specular white +MATERIAL 4 +RGB .98 .98 .98 +SPECEX 0 +SPECRGB .98 .98 .98 +REFL 1 +REFR 0 +REFRIOR 0 +EMITTANCE 0 +BSDF 0 + +// Reflective +MATERIAL 5 +RGB .98 .98 .98 +SPECEX 0 +SPECRGB .98 .98 .98 +REFL 1 +REFR 0 +REFRIOR 0 +EMITTANCE 0 +BSDF 1 + +// Refractive +MATERIAL 6 +RGB .98 .98 .98 +SPECEX 0 +SPECRGB .98 .55 .55 +REFL 1 +REFR 0 +REFRIOR 1.517 +EMITTANCE 0 +BSDF 2 + +// 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 + +// Floor +OBJECT 1 +cube +material 1 +TRANS 0 0 0 +ROTAT 0 0 0 +SCALE 10 .01 10 + +// Ceiling +OBJECT 2 +cube +material 1 +TRANS 0 10 0 +ROTAT 0 0 90 +SCALE .01 10 10 + +// Back wall +OBJECT 3 +cube +material 1 +TRANS 0 5 -5 +ROTAT 0 90 0 +SCALE .01 10 10 + +// Left wall +OBJECT 4 +cube +material 2 +TRANS -5 5 0 +ROTAT 0 0 0 +SCALE .01 10 10 + +// Right wall +OBJECT 5 +cube +material 3 +TRANS 5 5 0 +ROTAT 0 0 0 +SCALE .01 10 10 + +// Sphere +OBJECT 6 +sphere +material 6 +TRANS -1 4 -1 +ROTAT 0 0 0 +SCALE 3 3 3 diff --git a/scenes/cornellTwo.txt b/scenes/cornellTwo.txt new file mode 100644 index 0000000..762ff39 --- /dev/null +++ b/scenes/cornellTwo.txt @@ -0,0 +1,160 @@ +// Emissive material (light) +MATERIAL 0 +RGB 0.8 1 1 +SPECEX 0 +SPECRGB 0 0 0 +REFL 0 +REFR 0 +REFRIOR 0 +EMITTANCE 10 +BSDF -1 + +// Emissive material (light) +MATERIAL 1 +RGB 1 1 0.8 +SPECEX 0 +SPECRGB 0 0 0 +REFL 0 +REFR 0 +REFRIOR 0 +EMITTANCE 10 +BSDF -1 + +// Diffuse white +MATERIAL 2 +RGB .98 .98 .98 +SPECEX 0 +SPECRGB 0 0 0 +REFL 0 +REFR 0 +REFRIOR 0 +EMITTANCE 0 +BSDF 0 + +// Diffuse red +MATERIAL 3 +RGB .85 .35 .35 +SPECEX 0 +SPECRGB 0 0 0 +REFL 0 +REFR 0 +REFRIOR 0 +EMITTANCE 0 +BSDF 0 + +// Diffuse green +MATERIAL 4 +RGB .35 .85 .35 +SPECEX 0 +SPECRGB 0 0 0 +REFL 0 +REFR 0 +REFRIOR 0 +EMITTANCE 0 +BSDF 0 + +// Specular white +MATERIAL 5 +RGB .98 .98 .98 +SPECEX 0 +SPECRGB .98 .98 .98 +REFL 1 +REFR 0 +REFRIOR 0 +EMITTANCE 0 +BSDF 1 + +// Diffuse green +MATERIAL 6 +RGB .35 .65 .85 +SPECEX 0 +SPECRGB 0 0 0 +REFL 0 +REFR 0 +REFRIOR 0 +EMITTANCE 0 +BSDF 1 + +// Camera +CAMERA +RES 1080 1080 +FOVY 45 +ITERATIONS 7500 +DEPTH 8 +FILE cornellTwo +EYE 0.0 5 10.5 +LOOKAT 0 5 0 +UP 0 1 0 + + +// Ceiling light +OBJECT 0 +plane +material 0 +TRANS 3 9.99 0 +ROTAT 90 0 0 +SCALE 3 3 1 + +// Ceiling light 2 +OBJECT 1 +plane +material 1 +TRANS -3 9.99 0 +ROTAT 90 0 0 +SCALE 3 3 1 + +// Floor +OBJECT 2 +cube +material 2 +TRANS 0 0 0 +ROTAT 0 0 0 +SCALE 10 .01 10 + +// Ceiling +OBJECT 3 +cube +material 2 +TRANS 0 10 0 +ROTAT 0 0 90 +SCALE .01 10 10 + +// Back wall +OBJECT 4 +cube +material 2 +TRANS 0 5 -5 +ROTAT 0 90 0 +SCALE .01 10 10 + +// Left wall +OBJECT 5 +cube +material 3 +TRANS -5 5 0 +ROTAT 0 0 0 +SCALE .01 10 10 + +// Right wall +OBJECT 6 +cube +material 4 +TRANS 5 5 0 +ROTAT 0 0 0 +SCALE .01 10 10 + +// Sphere +OBJECT 7 +sphere +material 2 +TRANS -1 4 -4 +ROTAT 0 0 0 +SCALE 0.1 0.1 0.1 + +// Cube Mesh +OBJECT 8 +sphere +material 6 +TRANS 0 3 -2 +ROTAT 30 10 20 +SCALE 3 2 3 \ No newline at end of file diff --git a/scenes/cornellTwoCubes.txt b/scenes/cornellTwoCubes.txt new file mode 100644 index 0000000..675476a --- /dev/null +++ b/scenes/cornellTwoCubes.txt @@ -0,0 +1,141 @@ +// Emissive material (light) +MATERIAL 0 +RGB 0.8 1 1 +SPECEX 0 +SPECRGB 0 0 0 +REFL 0 +REFR 0 +REFRIOR 0 +EMITTANCE 5 +BSDF -1 + +// Emissive material (light) +MATERIAL 1 +RGB 1 1 0.8 +SPECEX 0 +SPECRGB 0 0 0 +REFL 0 +REFR 0 +REFRIOR 0 +EMITTANCE 5 +BSDF -1 + +// Diffuse white +MATERIAL 2 +RGB .98 .98 .98 +SPECEX 0 +SPECRGB 0 0 0 +REFL 0 +REFR 0 +REFRIOR 0 +EMITTANCE 0 +BSDF 0 + +// Diffuse red +MATERIAL 3 +RGB .85 .35 .35 +SPECEX 0 +SPECRGB 0 0 0 +REFL 0 +REFR 0 +REFRIOR 0 +EMITTANCE 0 +BSDF 0 + +// Diffuse green +MATERIAL 4 +RGB .35 .85 .35 +SPECEX 0 +SPECRGB 0 0 0 +REFL 0 +REFR 0 +REFRIOR 0 +EMITTANCE 0 +BSDF 0 + +// Specular white +MATERIAL 5 +RGB .98 .98 .98 +SPECEX 0 +SPECRGB .98 .98 .98 +REFL 1 +REFR 0 +REFRIOR 0 +EMITTANCE 0 +BSDF 0 + +// Camera +CAMERA +RES 800 800 +FOVY 45 +ITERATIONS 5000 +DEPTH 8 +FILE cornellTwo +EYE 0.0 5 10.5 +LOOKAT 0 5 0 +UP 0 1 0 + + +// Ceiling light +OBJECT 0 +cube +material 0 +TRANS 4 10 2.5 +ROTAT 0 0 0 +SCALE 2 .3 2 + +// Ceiling light 2 +OBJECT 1 +cube +material 1 +TRANS -3 10 -3 +ROTAT 0 0 0 +SCALE 2 .3 2 + +// Floor +OBJECT 2 +cube +material 2 +TRANS 0 0 0 +ROTAT 0 0 0 +SCALE 10 .01 10 + +// Ceiling +OBJECT 3 +cube +material 2 +TRANS 0 10 0 +ROTAT 0 0 90 +SCALE .01 10 10 + +// Back wall +OBJECT 4 +cube +material 2 +TRANS 0 5 -5 +ROTAT 0 90 0 +SCALE .01 10 10 + +// Left wall +OBJECT 5 +cube +material 3 +TRANS -5 5 0 +ROTAT 0 0 0 +SCALE .01 10 10 + +// Right wall +OBJECT 6 +cube +material 4 +TRANS 5 5 0 +ROTAT 0 0 0 +SCALE .01 10 10 + +// Sphere +OBJECT 7 +sphere +material 5 +TRANS -1 4 -1 +ROTAT 0 0 0 +SCALE 3 3 3 diff --git a/scenes/cube.mtl b/scenes/cube.mtl new file mode 100644 index 0000000..c89145d --- /dev/null +++ b/scenes/cube.mtl @@ -0,0 +1,6 @@ +newmtl initialShadingGroup +illum 4 +Kd 0.50 0.50 0.50 +Ka 0.00 0.00 0.00 +Tf 1.00 1.00 1.00 +Ni 1.00 diff --git a/scenes/cubeTest.txt b/scenes/cubeTest.txt new file mode 100644 index 0000000..716c72e --- /dev/null +++ b/scenes/cubeTest.txt @@ -0,0 +1,50 @@ +// Specular white +MATERIAL 0 +RGB .98 .98 .98 +SPECEX 0 +SPECRGB .98 .98 .98 +REFL 0 +REFR 0 +REFRIOR 0 +EMITTANCE 0 +BSDF 0 + +// Emissive Plane Light Material +MATERIAL 1 +RGB .98 .98 .98 +SPECEX 0 +SPECRGB .98 .98 .98 +REFL 0 +REFR 0 +REFRIOR 0 +EMITTANCE 5 +BSDF -1 + +// Camera +CAMERA +RES 800 800 +FOVY 45 +ITERATIONS 5000 +DEPTH 8 +FILE cubeTest +EYE 0.0 5 10.5 +LOOKAT 0 5 0 +UP 0 1 0 + +// LIGHT +OBJECT 0 +plane +material 1 +TRANS 5 5 5 +ROTAT 90 0 0 +SCALE 1 1 1 + + +// Mesh +OBJECT 1 +mesh +../scenes/cube.obj +material 0 +TRANS 2 3 0 +ROTAT 90 0 0 +SCALE 2 2 3 \ No newline at end of file diff --git a/scenes/dragon.mtl b/scenes/dragon.mtl new file mode 100644 index 0000000..9c09dab --- /dev/null +++ b/scenes/dragon.mtl @@ -0,0 +1,2 @@ +# WaveFront *.mtl file (generated by CINEMA 4D) + diff --git a/scenes/environment1.png b/scenes/environment1.png new file mode 100644 index 0000000..a272149 Binary files /dev/null and b/scenes/environment1.png differ diff --git a/scenes/environment2.jpg b/scenes/environment2.jpg new file mode 100644 index 0000000..d96f41c Binary files /dev/null and b/scenes/environment2.jpg differ diff --git a/scenes/environmentTest.txt b/scenes/environmentTest.txt new file mode 100644 index 0000000..f929be5 --- /dev/null +++ b/scenes/environmentTest.txt @@ -0,0 +1,35 @@ +// Specular white +MATERIAL 0 +RGB .98 .98 .98 +SPECEX 0 +SPECRGB .98 .98 .98 +REFL 1 +REFR 0 +REFRIOR 0 +EMITTANCE 0 +BSDF 1 + +// Camera +CAMERA +RES 800 800 +FOVY 45 +ITERATIONS 5000 +DEPTH 8 +FILE environmentTest +EYE 0.0 5 10.5 +LOOKAT 0 5 0 +UP 0 1 0 + + +// Sphere +OBJECT 0 +sphere +material 0 +TRANS 0 0 0 +ROTAT 0 0 0 +SCALE 3 3 3 + +// Environment +ENVIRONMENT +FILENAME ..\scenes\environment1.png +DIMENSIONS 1100 550 4 \ No newline at end of file diff --git a/scenes/environmentTest2.txt b/scenes/environmentTest2.txt new file mode 100644 index 0000000..9021e3f --- /dev/null +++ b/scenes/environmentTest2.txt @@ -0,0 +1,55 @@ +// Specular white +MATERIAL 0 +RGB .98 .98 .98 +SPECEX 0 +SPECRGB .98 .98 .98 +REFL 0 +REFR 0 +REFRIOR 0 +EMITTANCE 0 +BSDF 1 + +// Specular white +MATERIAL 1 +RGB .98 .98 .98 +SPECEX 0 +SPECRGB .98 .98 .98 +REFL 0 +REFR 0 +REFRIOR 0 +EMITTANCE 5 +BSDF -1 + +// Camera +CAMERA +RES 1080 1080 +FOVY 45 +ITERATIONS 1000 +DEPTH 8 +FILE environmentTest2 +EYE 0.0 5 10.5 +LOOKAT 0 5 0 +UP 0 1 0 + +// LIGHT +OBJECT 0 +plane +material 1 +TRANS 0 5 0 +ROTAT 90 0 0 +SCALE 1 1 1 + + +// Sphere +OBJECT 1 +sphere +material 0 +TRANS 0 5 0 +ROTAT 0 0 0 +SCALE 5 5 5 + + +// Environment +ENVIRONMENT +FILENAME ..\scenes\environment2.jpg +DIMENSIONS 8192 4096 3 \ No newline at end of file diff --git a/scenes/environmentTest3.txt b/scenes/environmentTest3.txt new file mode 100644 index 0000000..b436cc3 --- /dev/null +++ b/scenes/environmentTest3.txt @@ -0,0 +1,63 @@ +// Specular white +MATERIAL 0 +RGB .98 .98 .98 +SPECEX 0 +SPECRGB .98 .98 .98 +REFL 0 +REFR 0 +REFRIOR 0 +EMITTANCE 0 +BSDF 0 + +// Specular white +MATERIAL 1 +RGB .98 .98 .98 +SPECEX 0 +SPECRGB .98 .98 .98 +REFL 0 +REFR 0 +REFRIOR 0 +EMITTANCE 5 +BSDF -1 + +// Camera +CAMERA +RES 800 800 +FOVY 45 +ITERATIONS 5000 +DEPTH 8 +FILE environmentTest2 +EYE 0.0 5 10.5 +LOOKAT 0 5 0 +UP 0 1 0 + +// LIGHT +OBJECT 0 +plane +material 1 +TRANS 0 5 0 +ROTAT 90 0 0 +SCALE 1 1 1 + + +// Sphere +OBJECT 1 +sphere +material 0 +TRANS 0 5 0 +ROTAT 0 0 0 +SCALE 5 5 5 + +// Mesh +OBJECT 1 +mesh +../scenes/cube.obj +material 0 +TRANS 2 3 0 +ROTAT 90 0 0 +SCALE 2 2 3 + +// Environment +ENVIRONMENT +FILENAME ..\scenes\environment2.jpg +DIMENSIONS 8192 4096 3 \ No newline at end of file diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index a1cb3fb..e5b5ca7 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -1,5 +1,6 @@ set(SOURCE_FILES "stb.cpp" + "misHelper.h" "image.cpp" "image.h" "interactions.h" @@ -11,6 +12,8 @@ set(SOURCE_FILES "scene.cpp" "scene.h" "sceneStructs.h" + "tiny_obj_loader.cc" + "tiny_obj_loader.h" "preview.h" "preview.cpp" "utilities.cpp" @@ -19,5 +22,5 @@ set(SOURCE_FILES cuda_add_library(src ${SOURCE_FILES} - OPTIONS -arch=sm_20 + OPTIONS -arch=sm_61 ) diff --git a/src/helperKerns.h b/src/helperKerns.h new file mode 100644 index 0000000..fe5d42e --- /dev/null +++ b/src/helperKerns.h @@ -0,0 +1,16 @@ +#pragma once + +#include +#include + +//This would normally be wi.z in my CPU pathtracer, but this is not tangent space! +__host__ __device__ float CosTheta(const glm::vec3& n, const glm::vec3& wi) { + return glm::abs(glm::dot(n, wi)); +} + +__host__ __device__ bool SameHemisphere(const glm::vec3& normal, const glm::vec3& wi, const glm::vec3& wo) { + float dotWi = glm::dot(wi, normal); + float dotWo = glm::dot(wo, normal); + + return (dotWi >= 0 && dotWo >= 0) || (dotWi < 0 && dotWo < 0); +} \ No newline at end of file diff --git a/src/interactions.h b/src/interactions.h index 5ce3628..fcf99e7 100644 --- a/src/interactions.h +++ b/src/interactions.h @@ -1,6 +1,7 @@ #pragma once #include "intersections.h" +#include "misHelper.h" // CHECKITOUT /** @@ -9,36 +10,38 @@ */ __host__ __device__ glm::vec3 calculateRandomDirectionInHemisphere( - glm::vec3 normal, thrust::default_random_engine &rng) { - thrust::uniform_real_distribution u01(0, 1); - - float up = sqrt(u01(rng)); // cos(theta) - float over = sqrt(1 - up * up); // sin(theta) - float around = u01(rng) * TWO_PI; - - // Find a direction that is not the normal based off of whether or not the - // normal's components are all equal to sqrt(1/3) or whether or not at - // least one component is less than sqrt(1/3). Learned this trick from - // Peter Kutz. - - 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) { - directionNotNormal = glm::vec3(0, 1, 0); - } else { - directionNotNormal = glm::vec3(0, 0, 1); - } - - // Use not-normal direction to generate two perpendicular directions - glm::vec3 perpendicularDirection1 = - glm::normalize(glm::cross(normal, directionNotNormal)); - glm::vec3 perpendicularDirection2 = - glm::normalize(glm::cross(normal, perpendicularDirection1)); - - return up * normal - + cos(around) * over * perpendicularDirection1 - + sin(around) * over * perpendicularDirection2; + glm::vec3 normal, thrust::default_random_engine &rng) { + thrust::uniform_real_distribution u01(0, 1); + + float up = sqrt(u01(rng)); // cos(theta) + float over = sqrt(1 - up * up); // sin(theta) + float around = u01(rng) * TWO_PI; + + // Find a direction that is not the normal based off of whether or not the + // normal's components are all equal to sqrt(1/3) or whether or not at + // least one component is less than sqrt(1/3). Learned this trick from + // Peter Kutz. + + 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) { + directionNotNormal = glm::vec3(0, 1, 0); + } + else { + directionNotNormal = glm::vec3(0, 0, 1); + } + + // Use not-normal direction to generate two perpendicular directions + glm::vec3 perpendicularDirection1 = + glm::normalize(glm::cross(normal, directionNotNormal)); + glm::vec3 perpendicularDirection2 = + glm::normalize(glm::cross(normal, perpendicularDirection1)); + + return up * normal + + cos(around) * over * perpendicularDirection1 + + sin(around) * over * perpendicularDirection2; } /** @@ -47,11 +50,11 @@ glm::vec3 calculateRandomDirectionInHemisphere( * 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 @@ -68,12 +71,109 @@ glm::vec3 calculateRandomDirectionInHemisphere( */ __host__ __device__ void scatterRay( - PathSegment & pathSegment, - glm::vec3 intersect, - glm::vec3 normal, - const Material &m, - thrust::default_random_engine &rng) { - // TODO: implement this. - // A basic implementation of pure-diffuse shading will just call the - // calculateRandomDirectionInHemisphere defined above. + PathSegment & path, + ShadeableIntersection& isect, + const Material &m, + thrust::default_random_engine &rng, + float& pdf) { + + // Light + if (m.bsdf == -1) { + path.color *= (m.color * m.emittance); + path.remainingBounces = 0; + } + // Diffuse + else if (m.bsdf == 0) { + //Straight from Line 7 of my NaiveIntegrator.cpp + const glm::vec3 &n = isect.surfaceNormal; + const glm::vec3 &wo = -path.ray.direction; + + //This is "f". See Line 7 of LambertBRDF.cpp in my CPU Pathtracer + glm::vec3 accum_color = m.color * InvPi; + + //This is lamberFactor. See Line 23 of NaiveIntegrator.cpp in my CPU Pathtracer + float lambert_factor = fabs(glm::dot(n, wo)); + + //PDF Calculation + float dotWo = glm::dot(n, wo); + float cosTheta = fabs(dotWo) * InvPi; + pdf = cosTheta; + + if (pdf == 0) { + path.remainingBounces = 0; + return; + } + + glm::vec3 integral = (accum_color * lambert_factor) + / pdf; + path.color *= integral; + + //Scatter the Ray + path.ray.origin = isect.point + n*EPSILON; + path.ray.direction = calculateRandomDirectionInHemisphere(n, rng); + path.remainingBounces--; + } else if (m.bsdf == 1) { //Reflective + const glm::vec3 &n = isect.surfaceNormal; + + //Scatter the Ray + path.ray.origin = isect.point + n*EPSILON; + path.ray.direction = glm::reflect(path.ray.direction, n); + path.remainingBounces--; + pdf = 0.f; + } else if (m.bsdf == 2) { //Refractive + const glm::vec3 n = isect.surfaceNormal; + const glm::vec3 wo = -path.ray.direction; + + //Figure out which way we're going in->out or out -> in + //This is needed for incrementing the point along the normal + //and refraction + const bool entering = glm::dot(wo, n) > 0; + const float eta = entering ? 1 / m.indexOfRefraction : m.indexOfRefraction; + glm::vec3 faceforwardN = entering ? n : n; + + //Perform the Refraction + path.ray.direction = glm::refract(-wo, faceforwardN, eta); + + pdf = 0; + + //Increment based on whether or not we've changed mediums (!TIR) + const bool changedMediums = glm::dot(path.ray.direction, n) < 0.f; + const glm::vec3 increment = changedMediums ? -n*EPSILON : n*EPSILON; + + //Change the path according to calculations + path.ray.origin += increment; + path.color *= m.specular.color * 0.9f; + path.remainingBounces--; + } else { + //SHOULDN'T EVER HAPPEN + } + + path.ray.direction = glm::normalize(path.ray.direction); } + +//Picks a light at random, then gets the color at that a point on that light +__host__ __device__ glm::vec3 sample_li(const Geom& light, const Material& m, const glm::vec3& ref, thrust::default_random_engine &rng, glm::vec3 *wi, float* pdf_li) { + if (light.type == CUBE) { + //SAMPLE SHAPE + glm::vec3 shape_sample = sampleCube(light, ref, rng, pdf_li); + *wi = glm::normalize(shape_sample - ref); + + if (*pdf_li == 0 || shape_sample == ref) { + return glm::vec3(0.f); + } + + return m.color * m.emittance; + } else if (light.type == PLANE) { + //SAMPLE SHAPE + glm::vec3 shape_sample = samplePlane(light, ref, rng, pdf_li); + *wi = glm::normalize(shape_sample - ref); + + if (*pdf_li == 0 || shape_sample == ref) { + return glm::vec3(0.f); + } + + return m.color * m.emittance; + } + + return glm::vec3(0.f); +} \ No newline at end of file diff --git a/src/intersections.h b/src/intersections.h index 6f23872..c4832eb 100644 --- a/src/intersections.h +++ b/src/intersections.h @@ -45,7 +45,7 @@ __host__ __device__ glm::vec3 multiplyMV(glm::mat4 m, glm::vec4 v) { * @param outside Output param for whether the ray came from outside. * @return Ray parameter `t` value. -1 if no intersection. */ -__host__ __device__ float boxIntersectionTest(Geom box, Ray r, +__host__ __device__ float boxIntersectionTest(const Geom& box, Ray r, glm::vec3 &intersectionPoint, glm::vec3 &normal, bool &outside) { Ray q; q.origin = multiplyMV(box.inverseTransform, glm::vec4(r.origin , 1.0f)); @@ -140,5 +140,166 @@ __host__ __device__ float sphereIntersectionTest(Geom sphere, Ray r, normal = -normal; } - return glm::length(r.origin - intersectionPoint); + return glm::length(r.origin - intersectionPoint) > EPSILON ? glm::length(r.origin - intersectionPoint) : -1.f; } + + +/** +* Test intersection between a ray and a transformed sphere. Untransformed, +* the sphere always has radius 0.5 and is centered at the origin. +* +* @param intersectionPoint Output parameter for point of intersection. +* @param normal Output parameter for surface normal. +* @return Ray parameter `t` value. -1 if no intersection. +*/ +__host__ __device__ float planeIntersectionTest(Geom plane, Ray r, + glm::vec3 &intersectionPoint, glm::vec3 &normal) { + + Ray r_loc; + r_loc.origin = multiplyMV(plane.inverseTransform, glm::vec4(r.origin, 1.0f)); + r_loc.direction = multiplyMV(plane.inverseTransform, glm::vec4(r.direction, 0.0f)); + + float t = glm::dot(glm::vec3(0, 0, 1), (glm::vec3(0.5f, 0.5f, 0) - r_loc.origin)) / glm::dot(glm::vec3(0, 0, 1), r_loc.direction); + glm::vec3 p = glm::vec3(t * r_loc.direction + r_loc.origin); + + if (t > 0 && p.x >= -0.5f && p.x <= 0.5f && p.y >= -0.5f && p.y <= 0.5f) { + intersectionPoint = multiplyMV(plane.transform, glm::vec4(p,1)); + normal = glm::normalize(multiplyMV(plane.invTranspose, glm::vec4(0, 0, 1, 0))); + return t; + } + + return -1; +} + + +/*********************************************** + +MESH INTERSECTION ZONE +ENTER WITH CAUTION + +************************************************/ + +//Gets Area of a defined Triangle +__host__ __device__ float TriArea(const glm::vec3 &p1, const glm::vec3 &p2, const glm::vec3 &p3) +{ + return glm::length(glm::cross(p1 - p2, p3 - p2)) * 0.5f; +} + +//Does Barycentric Interpolation Between Triangle Points and reference point + +__host__ __device__ glm::vec3 getNormal(const Geom* tri, const glm::vec3& p) { + glm::vec3 p0 = tri->positions[0]; + glm::vec3 p1 = tri->positions[1]; + glm::vec3 p2 = tri->positions[2]; + + glm::vec3 n0 = tri->normals[0]; + glm::vec3 n1 = tri->normals[1]; + glm::vec3 n2 = tri->normals[2]; + + float A = TriArea(p0, p1, p2); + float A0 = TriArea(p1, p2, p); + float A1 = TriArea(p0, p2, p); + float A2 = TriArea(p0, p1, p); + + return glm::normalize(n0 * A0 / A + n1 * A1 / A + n2 * A2 / A); +} + +/** +* Test intersection between a ray and a transformed mesh. Untransformed, +* the mesh is cool. +* +* @param intersectionPoint Output parameter for point of intersection. +* @param normal Output parameter for surface normal. +* @return Ray parameter `t` value. -1 if no intersection. +*/ +//For some reason 561 code didn't work so I used: +//https://en.wikipedia.org/wiki/M%C3%B6ller%E2%80%93Trumbore_intersection_algorithm +__host__ __device__ float triangleIntersectionTest(const Geom& tri, Ray r_world, + glm::vec3 &intersectionPoint, glm::vec3 &normal, bool& outside) { + + Ray r; + r.origin = multiplyMV(tri.inverseTransform, glm::vec4(r_world.origin, 1.0f)); + r.direction = glm::normalize(multiplyMV(tri.inverseTransform, glm::vec4(r_world.direction, 0.0f))); + + glm::vec3 vertex0 = tri.positions[0]; + glm::vec3 vertex1 = tri.positions[1]; + glm::vec3 vertex2 = tri.positions[2]; + glm::vec3 edge1, edge2, h, s, q; + float a, f, u, v; + + edge1 = vertex1 - vertex0; + edge2 = vertex2 - vertex0; + + h = glm::cross(r.direction, edge2); + a = glm::dot(edge1, h); + + if (a > -EPSILON && a < EPSILON) + return -1; + f = 1.f / a; + s = r.origin - vertex0; + u = f * (glm::dot(s,h)); + if (u < EPSILON || u > 1.0 - EPSILON) return -1; + + q = glm::cross(s,edge1); + v = f * glm::dot(r.direction,q); + if (v < EPSILON || u + v > 1.0 - EPSILON) return -1; + // At this stage we can compute t to find out where the intersection point is on the line. + float t = f * glm::dot(edge2,q); + + //Local Normal: + glm::vec3 norm = glm::normalize(glm::cross(edge1, edge2)); + + if (t > EPSILON) // ray intersection + { + intersectionPoint = multiplyMV(tri.transform, glm::vec4(getPointOnRay(r,t),1.f)); + normal = getNormal(&tri, intersectionPoint); + outside = glm::dot(norm, r.direction) < EPSILON; + return glm::length(r_world.origin - intersectionPoint); + } + else // This means that there is a line intersection but not a ray intersection. + return -1; + +} + + +/** +* Test intersection between a ray and a transformed mesh. Untransformed, +* the mesh is cool. +* +* @param intersectionPoint Output parameter for point of intersection. +* @param normal Output parameter for surface normal. +* @return Ray parameter `t` value. -1 if no intersection. +*/ +/** +__host__ __device__ float meshIntersectionTest(Geom mesh, Ray r, + glm::vec3 &intersectionPoint, glm::vec3 &normal) { + + Ray r_loc; + r_loc.origin = multiplyMV(mesh.inverseTransform, glm::vec4(r.origin, 1.0f)); + r_loc.direction = multiplyMV(mesh.inverseTransform, glm::vec4(r.direction, 0.0f)); + + float closest_t = -1; + Triangle closestTri; + + //For every triangle: + for (int i = 0; i < mesh.tri_count; i++) { + int increment = (mesh.tri_index + i); + const Triangle tri = dev_triangles[increment]; + float tri_t = triangleIntersectionTest(tri, r_loc, intersectionPoint, normal); + if (tri_t > 0 && (tri_t < closest_t || closest_t < 0)) { + closest_t = tri_t; + closestTri = tri; + } + } + + if (closest_t > 0) + { + glm::vec3 p = glm::vec3(closest_t * r_loc.direction + r_loc.origin); + glm::vec3 n = getNormal(&closestTri, p); + return closest_t; + } + + return -1; +} +*/ + diff --git a/src/misHelper.h b/src/misHelper.h new file mode 100644 index 0000000..50a7100 --- /dev/null +++ b/src/misHelper.h @@ -0,0 +1,136 @@ +#pragma once + +#include +#include +#include "src\sceneStructs.h" +#include + + +/********************************************* +********************************************** +****** Multiple Importance ****** +****** Sampling ****** +********************************************** +**********************************************/ + +// This file has lots of important helper kernels for MIS + +__device__ __host__ float cubeArea(const Geom& light) { + return 2 * light.scale.x * light.scale.y * + 2 * light.scale.z * light.scale.y * + 2 * light.scale.x * light.scale.z; +} +__device__ __host__ float planeArea(const Geom& light) { + return light.scale.x * light.scale.y; +} + +__device__ __host__ glm::vec3 sampleCube(const Geom& light, const glm::vec3& ref, thrust::default_random_engine &rng, float* pdf) { + + //Get a sample point + thrust::uniform_real_distribution u_verts(0, 1); + glm::vec3 sample_li = glm::vec3(u_verts(rng) - 0.5f, 0, u_verts(rng) - 0.5f); + sample_li = glm::vec3(light.transform * glm::vec4(sample_li, 1)); + + glm::vec3 normal_li = glm::vec3(0, 0, -1); + + glm::vec3 wi = glm::normalize(sample_li - ref); + + //Get shape area and convert it to Solid angle + float cosT = fabs(glm::dot(-wi, normal_li)); + float solid_angle = (glm::length2(sample_li - ref) / cosT); + + *pdf = solid_angle / cubeArea(light); + + //Check if dividing by 0.f + *pdf = isnan(*pdf) ? 0.f : *pdf; + + return sample_li; +} + +__device__ __host__ glm::vec3 samplePlane(const Geom& light, const glm::vec3& ref, thrust::default_random_engine &rng, float* pdf) { + + //Get a sample point + thrust::uniform_real_distribution u_verts(0,1); + glm::vec3 sample_li_local = glm::vec3(u_verts(rng) - 0.5f, u_verts(rng) - 0.5f, 0); + glm::vec3 sample_li = multiplyMV(light.transform, glm::vec4(sample_li_local, 1)); + + glm::vec3 wi = glm::normalize(sample_li - ref); + + glm::vec3 normal_li = glm::vec3(light.invTranspose * glm::vec4(0, 0, 1, 0)); + *pdf = 1 / planeArea(light); + + //Get shape area and convert it to Solid angle + float cosT = glm::abs(glm::dot(-wi, normal_li)); + float solid_angle = (glm::length2(sample_li - ref) / cosT); + + *pdf *= solid_angle; + + return sample_li; +} + +__host__ __device__ float power_heuristic(int nf, float fpdf, int ng, float gpdf) { + float f = nf * fpdf, g = ng * gpdf; + if (fpdf == 0 && gpdf == 0) return 0.f; + + return (f*f) / + (f*f + g*g); +} + + +__host__ __device__ float pdf(int bsdf, const glm::vec3& wo, const glm::vec3& wi, const glm::vec3& n) { + if (bsdf == 0) { + //PDF Calculation + float dotWo = glm::dot(n, wo); + float cosTheta = fabs(dotWo) * InvPi; + return cosTheta; + } else if (bsdf == 1 || bsdf == 2) { + return 0.f; + } + else { + return 0.f; + } +} + +__host__ __device__ glm::vec3 f(const Material m, const glm::vec3& wo, const glm::vec3& wi) { + int bsdf = m.bsdf; + if (bsdf == 0) { + return m.color * InvPi; + } + else if (bsdf == 1 || bsdf == 2) { + return glm::vec3(0.f); + } + else { + return glm::vec3(0.3f); + } +} + +__host__ __device__ bool isBlack(const glm::vec3& vec) { + return vec[0] == 0 && vec[1] == 0 && vec[2] == 0; +} + +///__host__ __device__ float planeIntersectionTest(Geom plane, Ray r, glm::vec3 &intersectionPoint, glm::vec3 &normal) +__host__ __device__ float pdfLi(const Geom& light,const ShadeableIntersection& ref, const glm::vec3 wi) { + if (light.type == PLANE) { + //To Be Filled: + glm::vec3 isectLightPoint; + glm::vec3 normal; + //Input + Ray ray; + ray.origin = ref.point; + ray.direction = wi; + + if (planeIntersectionTest(light, ray, isectLightPoint, normal) < 0) { + return 0.f; + } + + return glm::length2(ref.point - isectLightPoint) / + (glm::abs(glm::dot(normal, -wi)) * planeArea(light)); + } else if (light.type == CUBE) { + //TODO ? + } + return 0.f; +} + +__host__ __device__ bool isSpecular(const int& bsdf) { + return bsdf == 1 || bsdf == 2; +} \ No newline at end of file diff --git a/src/pathtrace.cu b/src/pathtrace.cu index c1ec122..c119ee9 100644 --- a/src/pathtrace.cu +++ b/src/pathtrace.cu @@ -4,11 +4,15 @@ #include #include #include +#include +#include +#include #include "sceneStructs.h" #include "scene.h" #include "glm/glm.hpp" #include "glm/gtx/norm.hpp" +#include "glm/gtx/component_wise.hpp" #include "utilities.h" #include "pathtrace.h" #include "intersections.h" @@ -18,53 +22,78 @@ #define FILENAME (strrchr(__FILE__, '/') ? strrchr(__FILE__, '/') + 1 : __FILE__) #define checkCUDAError(msg) checkCUDAErrorFn(msg, FILENAME, __LINE__) + +#define CACHE_FIRST 0 + +#define MIS 1 + +#define ANTIALIAS 1 + +#define DOF 1 +#define LENS_RADIUS 3.f +#define FOCAL_DIST 12.f + void checkCUDAErrorFn(const char *msg, const char *file, int line) { #if ERRORCHECK - cudaDeviceSynchronize(); - cudaError_t err = cudaGetLastError(); - if (cudaSuccess == err) { - return; - } - - fprintf(stderr, "CUDA error"); - if (file) { - fprintf(stderr, " (%s:%d)", file, line); - } - fprintf(stderr, ": %s: %s\n", msg, cudaGetErrorString(err)); + cudaDeviceSynchronize(); + cudaError_t err = cudaGetLastError(); + if (cudaSuccess == err) { + return; + } + + fprintf(stderr, "CUDA error"); + if (file) { + fprintf(stderr, " (%s:%d)", file, line); + } + fprintf(stderr, ": %s: %s\n", msg, cudaGetErrorString(err)); # ifdef _WIN32 - getchar(); + getchar(); # endif - exit(EXIT_FAILURE); + exit(EXIT_FAILURE); #endif } + + +//Predicate functor +struct hasMoreBounces +{ + __host__ __device__ + bool operator()(const PathSegment& path) + { + return (path.remainingBounces > 0); + } +}; + + + __host__ __device__ thrust::default_random_engine makeSeededRandomEngine(int iter, int index, int depth) { - int h = utilhash((1 << 31) | (depth << 22) | iter) ^ utilhash(index); - return thrust::default_random_engine(h); + int h = utilhash((1 << 31) | (depth << 22) | iter) ^ utilhash(index); + return thrust::default_random_engine(h); } //Kernel that writes the image to the OpenGL PBO directly. __global__ void sendImageToPBO(uchar4* pbo, glm::ivec2 resolution, - int iter, glm::vec3* image) { - 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 pix = image[index]; - - glm::ivec3 color; - color.x = glm::clamp((int) (pix.x / iter * 255.0), 0, 255); - color.y = glm::clamp((int) (pix.y / iter * 255.0), 0, 255); - color.z = glm::clamp((int) (pix.z / iter * 255.0), 0, 255); - - // Each thread writes one pixel location in the texture (textel) - pbo[index].w = 0; - pbo[index].x = color.x; - pbo[index].y = color.y; - pbo[index].z = color.z; - } + int iter, glm::vec3* image) { + 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 pix = image[index]; + + glm::ivec3 color; + color.x = glm::clamp((int)(pix.x / iter * 255.0), 0, 255); + color.y = glm::clamp((int)(pix.y / iter * 255.0), 0, 255); + color.z = glm::clamp((int)(pix.z / iter * 255.0), 0, 255); + + // Each thread writes one pixel location in the texture (textel) + pbo[index].w = 0; + pbo[index].x = color.x; + pbo[index].y = color.y; + pbo[index].z = color.z; + } } static Scene * hst_scene = NULL; @@ -73,42 +102,101 @@ static Geom * dev_geoms = NULL; static Material * dev_materials = NULL; static PathSegment * dev_paths = NULL; static ShadeableIntersection * dev_intersections = NULL; +static ShadeableIntersection * dev_fst_bounce = NULL; // TODO: static variables for device memory, any extra info you need, etc // ... +static unsigned char* dev_environment = NULL; + +/********************************************* +********************************************** +****** Environment Mapping ****** +****** Functions ****** +********************************************** +**********************************************/ + +__device__ glm::vec3 getEnvMapColor(unsigned char* dev_environment, const glm::vec3& dir, const int& width, const int& height, const int& bpp) { + float x = dir.x, y = dir.y, z = dir.z; + + float u = atan2f(x, z) / (2 * PI) + 0.5f; + float v = y * 0.5f + 0.5f; + + v = 1-v; + + u *= width; + v *= height; + + int u_i = u; + int v_i = v; + + // Transform coordinates + unsigned char* r = dev_environment + bpp * (u_i + width*v_i); + + glm::vec3 color = glm::vec3(*(r + 0), *(r + 1), *(r + 2)); + + return glm::vec3(color.x, color.y, color.z) / 255.f; +} + + +void loadInEnvironment(const unsigned char* environment, const glm::ivec3 environment_dims) { + int width = environment_dims.x; + int height = environment_dims.y; + int bpp = environment_dims.z; + + // Allocate CUDA array in device memory + cudaMalloc(&dev_environment, bpp * width * height * sizeof(unsigned char)); + checkCUDAError("cudaMallocArray while mallocing texture array"); + + cudaMemcpy(dev_environment, environment, bpp * width * height * sizeof(unsigned char), cudaMemcpyHostToDevice); + checkCUDAError("cudaMempcpying while mallocing texture array"); +} + void pathtraceInit(Scene *scene) { - hst_scene = scene; - const Camera &cam = hst_scene->state.camera; - const int pixelcount = cam.resolution.x * cam.resolution.y; + hst_scene = scene; + const Camera &cam = hst_scene->state.camera; + const int pixelcount = cam.resolution.x * cam.resolution.y; - cudaMalloc(&dev_image, pixelcount * sizeof(glm::vec3)); - cudaMemset(dev_image, 0, pixelcount * sizeof(glm::vec3)); + cudaMalloc(&dev_image, pixelcount * sizeof(glm::vec3)); + cudaMemset(dev_image, 0, pixelcount * sizeof(glm::vec3)); - cudaMalloc(&dev_paths, pixelcount * sizeof(PathSegment)); + cudaMalloc(&dev_paths, pixelcount * sizeof(PathSegment)); - cudaMalloc(&dev_geoms, scene->geoms.size() * sizeof(Geom)); - cudaMemcpy(dev_geoms, scene->geoms.data(), scene->geoms.size() * sizeof(Geom), cudaMemcpyHostToDevice); + //Transfers Triangles from CPU to GPU - cudaMalloc(&dev_materials, scene->materials.size() * sizeof(Material)); - cudaMemcpy(dev_materials, scene->materials.data(), scene->materials.size() * sizeof(Material), cudaMemcpyHostToDevice); + cudaMalloc(&dev_geoms, scene->geoms.size() * sizeof(Geom)); + cudaMemcpy(dev_geoms, scene->geoms.data(), scene->geoms.size() * sizeof(Geom), cudaMemcpyHostToDevice); + checkCUDAError("copying dev_geoms"); - cudaMalloc(&dev_intersections, pixelcount * sizeof(ShadeableIntersection)); - cudaMemset(dev_intersections, 0, pixelcount * sizeof(ShadeableIntersection)); + cudaMalloc(&dev_materials, scene->materials.size() * sizeof(Material)); + cudaMemcpy(dev_materials, scene->materials.data(), scene->materials.size() * sizeof(Material), cudaMemcpyHostToDevice); + + cudaMalloc(&dev_intersections, pixelcount * sizeof(ShadeableIntersection)); + cudaMemset(dev_intersections, 0, pixelcount * sizeof(ShadeableIntersection)); - // TODO: initialize any extra device memeory you need + // TODO: initialize any extra device memeory you need + cudaMalloc(&dev_fst_bounce, pixelcount * sizeof(ShadeableIntersection)); + cudaMemset(dev_fst_bounce, 0, pixelcount * sizeof(ShadeableIntersection)); + checkCUDAError("pathtraceInit"); - checkCUDAError("pathtraceInit"); + //Load In Environment Map + if (scene->environment != NULL) { + loadInEnvironment(scene->environment, scene->environment_dims); + checkCUDAError("environmentLoading"); + } } void pathtraceFree() { - cudaFree(dev_image); // no-op if dev_image is null - cudaFree(dev_paths); - cudaFree(dev_geoms); - cudaFree(dev_materials); - cudaFree(dev_intersections); - // TODO: clean up any extra device memory you created - - checkCUDAError("pathtraceFree"); + cudaFree(dev_image); // no-op if dev_image is null + cudaFree(dev_paths); + cudaFree(dev_geoms); + cudaFree(dev_materials); + cudaFree(dev_intersections); + // TODO: clean up any extra device memory you created + cudaFree(dev_fst_bounce); + + if (dev_environment != NULL) { + cudaFree(dev_environment); + } } /** @@ -129,19 +217,99 @@ __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); +#if MIS + segment.color = glm::vec3(0.f); +#else + segment.color = glm::vec3(1.f); +#endif + segment.throughput = glm::vec3(1.f); - // TODO: implement antialiasing by jittering the ray +#if ANTIALIAS + thrust::default_random_engine rng = makeSeededRandomEngine(iter, x, y); + thrust::uniform_real_distribution u01(-0.5, 0.5); + float jitter_x = u01(rng); + float jitter_y = u01(rng); + + segment.ray.direction = glm::normalize(cam.view + - cam.right * cam.pixelLength.x * ((float)x - (float)cam.resolution.x * 0.5f + jitter_x) + - cam.up * cam.pixelLength.y * ((float)y - (float)cam.resolution.y * 0.5f + jitter_y) + ); +#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 segment.pixelIndex = index; segment.remainingBounces = traceDepth; } } +__host__ __device__ void getIntersection( + const Ray& ray + , Geom* geoms + , const int geoms_size + , ShadeableIntersection& intersection) { + + float t; + glm::vec3 intersect_point; + glm::vec3 normal; + float t_min = FLT_MAX; + int hit_geom_index = -1; + bool outside = true; + + glm::vec3 tmp_intersect; + glm::vec3 tmp_normal; + + // naive parse through global geoms + for (int i = 0; i < geoms_size; i++) + { + const Geom& geom = geoms[i]; + + //printf("Getting mesh intersection for: a %d type geom\n", geoms[i].type); + + if (geom.type == CUBE) + { + t = boxIntersectionTest(geom, ray, tmp_intersect, tmp_normal, outside); + } + else if (geom.type == SPHERE) + { + t = sphereIntersectionTest(geom, ray, tmp_intersect, tmp_normal, outside); + } + else if (geom.type == PLANE) { + t = planeIntersectionTest(geom, ray, tmp_intersect, tmp_normal); + } + else if (geom.type == TRIANGLE) { + t = triangleIntersectionTest(geom, ray, tmp_intersect, tmp_normal, outside); + } + // TODO: add more intersection tests here... triangle? metaball? CSG? + + // Compute the minimum t from the intersection tests to determine what + // scene geometry object was hit first. + if (t > 0.0f && t_min > t) + { + t_min = t; + hit_geom_index = i; + intersect_point = tmp_intersect; + normal = tmp_normal; + } + } + + if (hit_geom_index == -1) + { + intersection.t = -1.0f; + } + else + { + //The ray hits something + intersection.t = t_min; + intersection.materialId = geoms[hit_geom_index].materialid; + intersection.surfaceNormal = normal; + intersection.point = intersect_point; + } +} + // TODO: // computeIntersections handles generating ray intersections ONLY. // Generating new rays is handled in your shader(s). @@ -153,7 +321,7 @@ __global__ void computeIntersections( , Geom * geoms , int geoms_size , ShadeableIntersection * intersections - ) +) { int path_index = blockIdx.x * blockDim.x + threadIdx.x; @@ -161,108 +329,443 @@ __global__ void computeIntersections( { PathSegment pathSegment = pathSegments[path_index]; - float t; - glm::vec3 intersect_point; - glm::vec3 normal; - float t_min = FLT_MAX; - int hit_geom_index = -1; - bool outside = true; + getIntersection(pathSegment.ray, geoms, geoms_size, intersections[path_index]); + } +} - glm::vec3 tmp_intersect; - glm::vec3 tmp_normal; +// A Naiive Integrator +__global__ void shadeMaterialNaive( + int iter + , int num_paths + , ShadeableIntersection * shadeableIntersections + , PathSegment * pathSegments + , Material * materials +) +{ + int idx = blockIdx.x * blockDim.x + threadIdx.x; + if (idx < num_paths) + { + ShadeableIntersection intersection = shadeableIntersections[idx]; + PathSegment& path = pathSegments[idx]; + if (intersection.t > 0.0f) { // if the intersection exists... + // Set up the RNG + // LOOK: this is how you use thrust's RNG! Please look at + // makeSeededRandomEngine as well. + thrust::default_random_engine rng = makeSeededRandomEngine(iter, idx, path.remainingBounces); + thrust::uniform_real_distribution u01(0, 1); - // naive parse through global geoms + Material m = materials[intersection.materialId]; - for (int i = 0; i < geoms_size; i++) - { - Geom & geom = geoms[i]; + // If the material indicates that the object was a light, "light" the ray + float pdf; + scatterRay(path, intersection, m, rng, pdf); + + } + else { // If there was no intersection, color the ray black. + path.color = glm::vec3(0.0f); + path.remainingBounces = 0; + } + } +} + +// A Direct Lighting Integrator +__global__ void shadeMaterialDirect( + int iter + , int depth, int depthLimit + , int light_count + , int geoms_size + , int num_paths + , ShadeableIntersection * shadeableIntersections + , PathSegment * pathSegments + , Material * materials + , Geom* geoms +) +{ + int idx = blockIdx.x * blockDim.x + threadIdx.x; + if (idx < num_paths) { + ShadeableIntersection intersection = shadeableIntersections[idx]; + PathSegment& path = pathSegments[idx]; + if (intersection.t <= 0.0f) { + path.color = glm::vec3(0.0f); + path.remainingBounces = 0; + return; + } + + // if the intersection exists... + thrust::default_random_engine rng = makeSeededRandomEngine(iter, idx, 0); + thrust::uniform_real_distribution u01(0, 1); + + Material isect_m = materials[intersection.materialId]; + const glm::vec3& wo = -path.ray.direction; + + + if (depth == 0 || path.specularBounce) { + //405: Assumption: light is emitted equally from/to all directions + glm::vec3 Le = isect_m.color * isect_m.emittance; + path.color += path.throughput * Le; + } + + if (isect_m.emittance > 0.f) { + path.color = isect_m.color * isect_m.emittance; + path.remainingBounces = 0; + return; + } - if (geom.type == CUBE) - { - t = boxIntersectionTest(geom, pathSegment.ray, tmp_intersect, tmp_normal, outside); + path.specularBounce = isSpecular(isect_m.bsdf); + + if (!path.specularBounce) { + /*************************************** + **************************************** + **************************************** + ******* Light Importance Sampling ****** + **************************************** + **************************************** + ****************************************/ + //Get f(X) and L(X) + glm::vec3 wi; + float pdf_li = 1.f; + thrust::uniform_real_distribution u02(0, light_count); + int rand_li = u02(rng); + + const Geom light = geoms[rand_li]; + const Material light_m = materials[rand_li]; + glm::vec3 li_x = sample_li(light, light_m, intersection.point, rng, &wi, &pdf_li); //Assuming lights give equal light from anywhere + if (pdf_li < ZeroEpsilon) { + // This is the shadow feeling part of my CPU Code: + // Lines 71-81 + Ray dir_light; + dir_light.origin = intersection.point + intersection.surfaceNormal * EPSILON; + dir_light.direction = wi; + ShadeableIntersection shadow_isect; + getIntersection(dir_light, geoms, geoms_size, shadow_isect); + + // zero out contribution if it doesn't hit anything + bool shadowed = shadow_isect.t > 0.f && (shadow_isect.materialId != light.materialid); + li_x = shadowed ? glm::vec3(0) : li_x; + + const float pdf_bsdf = pdf(isect_m.bsdf, wo, -wi, intersection.surfaceNormal); + + ////This only works because we have one bsdf in each material + const glm::vec3 f_x = f(isect_m, wo, wi) * glm::abs(glm::dot(-wi, intersection.surfaceNormal)); + + glm::vec3 Ld = (f_x * li_x) + / (pdf_li); + + Ld *= light_count; + + path.color += Ld; + path.remainingBounces = 0; + return; } - else if (geom.type == SPHERE) - { - t = sphereIntersectionTest(geom, pathSegment.ray, tmp_intersect, tmp_normal, outside); + } + } +} + + +// A Broken Integrator +__global__ void shadeMaterialMIS( + int iter + , int depth, int depthLimit + , const int light_count + , const int geoms_size + , const int num_paths + , unsigned char* dev_environment, const int map_width, const int map_height, const int map_bpp + , ShadeableIntersection * shadeableIntersections + , PathSegment * pathSegments + , Material * materials + , Geom* geoms +) +{ + int idx = blockIdx.x * blockDim.x + threadIdx.x; + if (idx < num_paths) { + ShadeableIntersection intersection = shadeableIntersections[idx]; + PathSegment& path = pathSegments[idx]; + + if (intersection.t < EPSILON) { + //DO ENVIRONMENT MAPPING + if (dev_environment != NULL) { + path.color = getEnvMapColor(dev_environment, path.ray.direction, map_width, map_height, map_bpp); } - // TODO: add more intersection tests here... triangle? metaball? CSG? - - // Compute the minimum t from the intersection tests to determine what - // scene geometry object was hit first. - if (t > 0.0f && t_min > t) - { - t_min = t; - hit_geom_index = i; - intersect_point = tmp_intersect; - normal = tmp_normal; + else { + path.color = glm::vec3(0.f); } + + path.remainingBounces = 0; + return; } - if (hit_geom_index == -1) - { - intersections[path_index].t = -1.0f; + // if the intersection exists... + thrust::default_random_engine rng = makeSeededRandomEngine(iter, idx, path.remainingBounces); + thrust::uniform_real_distribution u01(0, 1); + + Material isect_m = materials[intersection.materialId]; + const glm::vec3& wo = -path.ray.direction; + + if (depth == 0 || path.specularBounce) { + //405: Assumption: light is emitted equally from/to all directions + glm::vec3 Le = isect_m.color * isect_m.emittance; + path.color += path.throughput * Le; } - else - { - //The ray hits something - intersections[path_index].t = t_min; - intersections[path_index].materialId = geoms[hit_geom_index].materialid; - intersections[path_index].surfaceNormal = normal; + + if (isect_m.emittance > 0.f) { + path.remainingBounces = 0; + return; + } + + //Store a copy. We only add this in the end. + PathSegment gi_Component; + gi_Component.color = glm::vec3(1.f); + gi_Component.ray = path.ray; + float gi_pdf; + // If the material indicates that the object was a light, "light" the ray + scatterRay(gi_Component, intersection, isect_m, rng, gi_pdf); + + __syncthreads(); + + //Random Light Selection + int rand_li = u01(rng)*light_count; + const Geom light = geoms[rand_li]; + const Material light_m = materials[rand_li]; + + //At this point, we've scattered, sampled and gi_component now has a new direction and origin. + // \ ^ + // \ / + // \ / <---- gi_Component + // \/ + // + + if (!path.specularBounce) { + /*************************************** + **************************************** + **************************************** + ******* Light Importance Sampling ****** + **************************************** + **************************************** + ****************************************/ + glm::vec3 Ld = glm::vec3(0.f); + + //Get f(X) and L(X) + glm::vec3 wi; + float pdf_li = 1.f; + + glm::vec3 li_x = sample_li(light, light_m, intersection.point, rng, &wi, &pdf_li); //Assuming lights give equal light from anywhere + if (pdf_li > ZeroEpsilon) { + // This is the shadow feeling part of my CPU Code: + // Lines 71-81 + Ray dir_light_ray; + dir_light_ray.origin = intersection.point + intersection.surfaceNormal * EPSILON; + dir_light_ray.direction = wi; + ShadeableIntersection shadow_isect; + getIntersection(dir_light_ray, geoms, geoms_size, shadow_isect); + + // zero out contribution if it doesn't hit anything + bool shadowed = shadow_isect.t > 0.f && (shadow_isect.materialId != light.materialid); + li_x = shadowed ? glm::vec3(0) : li_x; + + const float pdf_bsdf = pdf(isect_m.bsdf, wo, -wi, intersection.surfaceNormal); + + ////This only works because we have one bsdf in each material + const glm::vec3 f_x = f(isect_m, wo, wi) * glm::abs(glm::dot(-wi, intersection.surfaceNormal)); + + float weight_li = power_heuristic(1, pdf_li, 1, pdf_bsdf); + + Ld = (f_x * li_x * weight_li * path.throughput) + / (pdf_li); + }//END DIRECT LIGHTING + + /*************************************** + **************************************** + **************************************** + ******* BSDF Importance Sampling ******* + **************************************** + **************************************** + ****************************************/ + + //Get f(X) and L(X) + PathSegment indirect_path; + float bsdf_pdf; + indirect_path.ray = path.ray; + indirect_path.color = glm::vec3(1.f); + scatterRay(indirect_path, intersection, isect_m, rng, bsdf_pdf); + wi = indirect_path.ray.direction; + glm::vec3 f_y = indirect_path.color; + + //Only do calculations if bsdfpdf is not zero for efficiency + if (bsdf_pdf > ZeroEpsilon) { + ShadeableIntersection bsdf_direct_isx; + indirect_path.ray.origin = intersection.point + intersection.surfaceNormal * EPSILON; + getIntersection(indirect_path.ray, geoms, geoms_size, bsdf_direct_isx); + + pdf_li = pdfLi(light, intersection, wi); + + //Only add cotribution if object hit is the light + if (bsdf_direct_isx.t > 0 && bsdf_direct_isx.materialId == light.materialid) { + float weight_bsdf = power_heuristic(1, bsdf_pdf, 1, pdf_li); + + glm::vec3 li_y = light_m.emittance * light_m.color; + + Ld += li_y * f_y * weight_bsdf * path.throughput; + } + } + + Ld *= light_count; + + //**************************************** + //**Add Ld to Ray color before GI Stuff*** + //**************************************** + path.color += Ld; + } + + //Update Scene_Ray - This just spawns a new ray for the next loop + path.ray = gi_Component.ray; + path.throughput *= gi_Component.color; + path.specularBounce = isSpecular(isect_m.bsdf); + + if ((isBlack(gi_Component.color)) || intersection.materialId == light.materialid) { + path.remainingBounces = 0; + return; + } + + //Russian Roulette Early Ray Termination + if (depth >= 3) { + float q = glm::max(0.05f, (1 - glm::compMax(path.throughput))); + if (u01(rng) < q) { + path.remainingBounces = 0; + return; + } + path.throughput /= (1 - q); } } } -// LOOK: "fake" shader demonstrating what you might do with the info in -// a ShadeableIntersection, as well as how to use thrust's random number -// generator. Observe that since the thrust random number generator basically -// adds "noise" to the iteration, the image should start off noisy and get -// cleaner as more iterations are computed. -// -// Note that this shader does NOT do a BSDF evaluation! -// Your shaders should handle that - this can allow techniques such as -// bump mapping. -__global__ void shadeFakeMaterial ( - int iter - , int num_paths +// A Backup/Just in Case Integrator +// Who needs source control amirite +__global__ void shadeMaterialMIS_backup( + int iter + , int depth, int depthLimit + , int light_count + , int geoms_size + , int num_paths , ShadeableIntersection * shadeableIntersections , PathSegment * pathSegments , Material * materials - ) + , Geom* geoms +) { - int idx = blockIdx.x * blockDim.x + threadIdx.x; - if (idx < num_paths) - { - ShadeableIntersection intersection = shadeableIntersections[idx]; - if (intersection.t > 0.0f) { // if the intersection exists... - // Set up the RNG - // LOOK: this is how you use thrust's RNG! Please look at - // makeSeededRandomEngine as well. - thrust::default_random_engine rng = makeSeededRandomEngine(iter, idx, 0); - thrust::uniform_real_distribution u01(0, 1); - - Material material = materials[intersection.materialId]; - glm::vec3 materialColor = material.color; - - // If the material indicates that the object was a light, "light" the ray - if (material.emittance > 0.0f) { - pathSegments[idx].color *= (materialColor * material.emittance); - } - // Otherwise, do some pseudo-lighting computation. This is actually more - // like what you would expect from shading in a rasterizer like OpenGL. - // TODO: replace this! you should be able to start with basically a one-liner - else { - float lightTerm = glm::dot(intersection.surfaceNormal, glm::vec3(0.0f, 1.0f, 0.0f)); - pathSegments[idx].color *= (materialColor * lightTerm) * 0.3f + ((1.0f - intersection.t * 0.02f) * materialColor) * 0.7f; - pathSegments[idx].color *= u01(rng); // apply some noise because why not - } - // If there was no intersection, color the ray black. - // Lots of renderers use 4 channel color, RGBA, where A = alpha, often - // used for opacity, in which case they can indicate "no opacity". - // This can be useful for post-processing and image compositing. - } else { - pathSegments[idx].color = glm::vec3(0.0f); - } - } + int idx = blockIdx.x * blockDim.x + threadIdx.x; + if (idx < num_paths) { + ShadeableIntersection intersection = shadeableIntersections[idx]; + PathSegment& path = pathSegments[idx]; + if (intersection.t <= 0.0f) { + path.color = glm::vec3(0.0f); + path.remainingBounces = 0; + return; + } + + // if the intersection exists... + thrust::default_random_engine rng = makeSeededRandomEngine(iter, idx, path.remainingBounces); + thrust::uniform_real_distribution u01(0, 1); + + Material isect_m = materials[intersection.materialId]; + const glm::vec3& wo = -path.ray.direction; + + + if (depth == 0 || path.specularBounce) { + //405: Assumption: light is emitted equally from/to all directions + glm::vec3 Le = isect_m.color * isect_m.emittance; + path.color += path.throughput * Le; + } + + if (isect_m.emittance > 0.f) { + path.color = isect_m.color * isect_m.emittance; + path.remainingBounces = 0; + return; + } + + //Store a copy. We only add this in the end. + PathSegment gi_Component = PathSegment(path); + gi_Component.color = glm::vec3(1.f); + float gi_pdf; + + // If the material indicates that the object was a light, "light" the ray + scatterRay(gi_Component, intersection, isect_m, rng, gi_pdf); + + path.specularBounce = isect_m.bsdf == 1 || isect_m.bsdf == 2; + + thrust::uniform_real_distribution u02(0, light_count); + int rand_li = u02(rng); + const Geom light = geoms[rand_li]; + const Material light_m = materials[rand_li]; + + + //At this point, we've scattered, sampled and gi_component now has a new direction and origin. + // \ ^ + // \ / + // \ / <---- gi_Component + // \/ + // + if (!path.specularBounce) { + /*************************************** + **************************************** + **************************************** + ******* Light Importance Sampling ****** + **************************************** + **************************************** + ****************************************/ + //Get f(X) and L(X) + glm::vec3 wi; + float pdf_li = 1.f; + glm::vec3 li_x = sample_li(light, light_m, intersection.point, rng, &wi, &pdf_li); //Assuming lights give equal light from anywhere + if (pdf_li != 0.f) { + // This is the shadow feeling part of my CPU Code: + // Lines 71-81 + Ray dir_light; + dir_light.origin = intersection.point + intersection.surfaceNormal * EPSILON; + dir_light.direction = wi; + ShadeableIntersection shadow_isect; + getIntersection(dir_light, geoms, geoms_size, shadow_isect); + + // zero out contribution if it doesn't hit anything + bool shadowed = shadow_isect.t > 0.f && (shadow_isect.materialId != light.materialid); + li_x = shadowed ? glm::vec3(0) : li_x; + + const float pdf_bsdf = pdf(isect_m.bsdf, wo, -wi, intersection.surfaceNormal); + + ////This only works because we have one bsdf in each material + const glm::vec3 f_x = f(isect_m, wo, wi) * glm::abs(glm::dot(-wi, intersection.surfaceNormal)); + + float weight_li = power_heuristic(1, pdf_li, 1, pdf_bsdf); + + glm::vec3 Ld = (f_x * li_x * weight_li) + / (pdf_li); + + //Ld *= path.throughput; + //DELET THIS + Ld *= light_count; + + path.color += Ld; + //DELET THIS + //path.color = Ld; + path.remainingBounces = 0; + return; + } + + if (isBlack(gi_Component.color) || intersection.materialId == light.materialid || gi_pdf == 0.f) { + path.remainingBounces = 0; + + if (isBlack(gi_Component.color)) printf("IS HOMIE HOMIE \n"); + if (intersection.materialId == light.materialid) printf("CASE 2: MAT THING \n"); + if (gi_pdf == 0) printf("GI_PDF IS ZERO HOMIE \n"); + return; + } + + } + + //Update Scene_Ray - This just spawns a new ray for the next loop + path.ray = gi_Component.ray; + path.throughput *= gi_Component.color; + path.remainingBounces--; + } } // Add the current iteration's output to the overall image @@ -277,56 +780,45 @@ __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 */ void pathtrace(uchar4 *pbo, int frame, int iter) { - const int traceDepth = hst_scene->state.traceDepth; - const Camera &cam = hst_scene->state.camera; - const int pixelcount = cam.resolution.x * cam.resolution.y; + const int traceDepth = hst_scene->state.traceDepth; + const Camera &cam = hst_scene->state.camera; + const int pixelcount = cam.resolution.x * cam.resolution.y; // 2D block for generating ray from 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); + const dim3 blockSize2d(8, 8); + const dim3 blocksPerGrid2d( + (cam.resolution.x + blockSize2d.x - 1) / blockSize2d.x, + (cam.resolution.y + blockSize2d.y - 1) / blockSize2d.y); // 1D block for path tracing const int blockSize1d = 128; - /////////////////////////////////////////////////////////////////////////// - - // Recap: - // * Initialize array of path rays (using rays that come out of the camera) - // * You can pass the Camera object to that kernel. - // * Each path ray must carry at minimum a (ray, color) pair, - // * where color starts as the multiplicative identity, white = (1, 1, 1). - // * This has already been done for you. - // * For each depth: - // * Compute an intersection in the scene for each path ray. - // A very naive version of this has been implemented for you, but feel - // free to add more primitives and/or a better algorithm. - // Currently, intersection distance is recorded as a parametric distance, - // t, or a "distance along the ray." t = -1.0 indicates no intersection. - // * Color is attenuated (multiplied) by reflections off of any object - // * TODO: Stream compact away all of the terminated paths. - // You may use either your implementation or `thrust::remove_if` or its - // cousins. - // * Note that you can't really use a 2D kernel launch any more - switch - // to 1D. - // * TODO: Shade the rays that intersected something or didn't bottom out. - // That is, color the ray by performing a color computation according - // to the shader, then generate a new ray to continue the ray path. - // We recommend just updating the ray's PathSegment in place. - // Note that this step may come before or after stream compaction, - // since some shaders you write may also cause a path to terminate. - // * Finally, add this iteration's results to the image. This has been done - // for you. - - // TODO: perform one iteration of path tracing - - generateRayFromCamera <<>>(cam, iter, traceDepth, dev_paths); + /////////////////////////////////////////////////////////////////////////// + + // Recap: + // * Initialize array of path rays (using rays that come out of the camera) + // * You can pass the Camera object to that kernel. + // * Each path ray must carry at minimum a (ray, color) pair, + // * where color starts as the multiplicative identity, white = (1, 1, 1). + // * This has already been done for you. + // * For each depth: + // * Compute an intersection in the scene for each path ray. + // A very naive version of this has been implemented for you, but feel + // free to add more primitives and/or a better algorithm. + // Currently, intersection distance is recorded as a parametric distance, + // t, or a "distance along the ray." t = -1.0 indicates no intersection. + // * Color is attenuated (multiplied) by reflections off of any object + // * Finally, add this iteration's results to the image. This has been done + // for you. + + // TODO: perform one iteration of path tracing + generateRayFromCamera << > > (cam, iter, traceDepth, dev_paths); checkCUDAError("generate camera ray"); int depth = 0; @@ -336,58 +828,98 @@ void pathtrace(uchar4 *pbo, int frame, int iter) { // --- PathSegment Tracing Stage --- // Shoot ray into scene, bounce between objects, push shading chunks - bool iterationComplete = false; - while (!iterationComplete) { + bool iterationComplete = false; + while (!iterationComplete && depth < traceDepth) { + + // clean shading chunks + cudaMemset(dev_intersections, 0, pixelcount * sizeof(ShadeableIntersection)); + + // tracing + dim3 numblocksPathSegmentTracing = (num_paths + blockSize1d - 1) / blockSize1d; + +#if CACHE_FIRST + if ((depth == 0 && iter == 1) || depth > 0) { + computeIntersections << > > ( + depth + , num_paths + , dev_paths + , dev_geoms + , hst_scene->geoms.size() + , dev_intersections + ); + checkCUDAError("trace one bounce"); + cudaDeviceSynchronize(); + + if (depth == 0) { + cudaMemcpy(dev_fst_bounce, dev_intersections, sizeof(ShadeableIntersection) * num_paths, cudaMemcpyDeviceToDevice); + } + } + else if (depth == 0) { + cudaMemcpy(dev_intersections, dev_fst_bounce, sizeof(ShadeableIntersection) * num_paths, cudaMemcpyDeviceToDevice); + } +#else + computeIntersections << > > ( + depth + , num_paths + , dev_paths + , dev_geoms + , hst_scene->geoms.size() + , dev_intersections + ); + checkCUDAError("computed intersections"); + cudaDeviceSynchronize(); +#endif + + // --- Shading Stage --- +#if MIS + int lc = hst_scene->light_count; + shadeMaterialMIS << > > ( + iter, + depth, hst_scene->state.traceDepth, + lc, + hst_scene->geoms.size(), + num_paths, + dev_environment, hst_scene->environment_dims[0], hst_scene->environment_dims[1], hst_scene->environment_dims[2], + dev_intersections, + dev_paths, + dev_materials, + dev_geoms + ); + checkCUDAError("shadeMaterialMIS"); +#else + shadeMaterialNaive << > > ( + iter, + num_paths, + dev_intersections, + dev_paths, + dev_materials + ); +#endif + + cudaDeviceSynchronize(); - // clean shading chunks - cudaMemset(dev_intersections, 0, pixelcount * sizeof(ShadeableIntersection)); + // --- Stream Compaction + PathSegment* remaining_end = + thrust::partition(thrust::device, dev_paths, dev_paths + num_paths, hasMoreBounces()); + num_paths = remaining_end - dev_paths; - // 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. + iterationComplete = num_paths == 0; + depth++; } - // Assemble this iteration and apply it to the image - dim3 numBlocksPixels = (pixelcount + blockSize1d - 1) / blockSize1d; - finalGather<<>>(num_paths, dev_image, dev_paths); + // Assemble this iteration and apply it to the image + dim3 numBlocksPixels = (pixelcount + blockSize1d - 1) / blockSize1d; + num_paths = dev_path_end - dev_paths; + finalGather << > > (num_paths, dev_image, dev_paths); - /////////////////////////////////////////////////////////////////////////// + /////////////////////////////////////////////////////////////////////////// - // Send results to OpenGL buffer for rendering - sendImageToPBO<<>>(pbo, cam.resolution, iter, dev_image); + // Send results to OpenGL buffer for rendering + sendImageToPBO << > > (pbo, cam.resolution, iter, dev_image); - // Retrieve image from GPU - cudaMemcpy(hst_scene->state.image.data(), dev_image, - pixelcount * sizeof(glm::vec3), cudaMemcpyDeviceToHost); + // Retrieve image from GPU + cudaMemcpy(hst_scene->state.image.data(), dev_image, + pixelcount * sizeof(glm::vec3), cudaMemcpyDeviceToHost); - checkCUDAError("pathtrace"); + checkCUDAError("pathtrace"); } diff --git a/src/scene.cpp b/src/scene.cpp index cbae043..9dc8fb7 100644 --- a/src/scene.cpp +++ b/src/scene.cpp @@ -1,9 +1,16 @@ #include #include "scene.h" +#include "tiny_obj_loader.h" #include +#include +#include #include #include +//DELET THIS +#include +#define GetCurrentDir _getcwd + Scene::Scene(string filename) { cout << "Reading scene from " << filename << " ..." << endl; cout << " " << endl; @@ -13,6 +20,7 @@ Scene::Scene(string filename) { cout << "Error reading from file - aborting!" << endl; throw; } + while (fp_in.good()) { string line; utilityCore::safeGetline(fp_in, line); @@ -27,7 +35,11 @@ Scene::Scene(string filename) { } else if (strcmp(tokens[0].c_str(), "CAMERA") == 0) { loadCamera(); cout << " " << endl; - } + } else if (strcmp(tokens[0].c_str(), "ENVIRONMENT") == 0) { + loadEnvironment(); + cout << " " << endl; + cout << "Ooh.. environment map, fancy!" << endl; + } } } } @@ -41,6 +53,8 @@ int Scene::loadGeom(string objectid) { cout << "Loading Geom " << id << "..." << endl; Geom newGeom; string line; + bool loadMesh = false; + string meshFile; //load object type utilityCore::safeGetline(fp_in, line); @@ -52,6 +66,17 @@ int Scene::loadGeom(string objectid) { cout << "Creating new cube..." << endl; newGeom.type = CUBE; } + else if (strcmp(line.c_str(), "plane") == 0) { + cout << "Creating new plane..." << endl; + newGeom.type = PLANE; + } + else if (strcmp(line.c_str(), "mesh") == 0) { + cout <<"Found mesh..." << endl; + utilityCore::safeGetline(fp_in, line); + meshFile = line; + loadMesh = true; + newGeom.type = TRIANGLE; + } } //link material @@ -84,11 +109,93 @@ int Scene::loadGeom(string objectid) { newGeom.inverseTransform = glm::inverse(newGeom.transform); newGeom.invTranspose = glm::inverseTranspose(newGeom.transform); - geoms.push_back(newGeom); + if (loadMesh) { + loadOBJ(newGeom, meshFile); + } + else { + geoms.push_back(newGeom); + } + return 1; } } +void Scene::loadOBJ(Geom& base_tri, string& filename) +{ + std::vector shapes; + std::vector materials; // will be discarded + tinyobj::attrib_t attributes; + string errors; + tinyobj::LoadObj(&attributes, &shapes, &materials, &errors, filename.c_str()); + + if (!errors.empty()) { + printf("Error loading obj in loadOBJ because: %s\n", errors.c_str()); + } + + if (errors.size() == 0) { + for (unsigned int i = 0; i < shapes.size(); i++) + { + std::vector &positions = attributes.vertices; + std::vector &normals = attributes.normals; + std::vector &uvs = attributes.texcoords; + const std::vector &indices = shapes[i].mesh.indices; + + for (unsigned int j = 0; j < indices.size() / 3; j ++) + { + Geom tri; + tri.type = TRIANGLE; + tri.materialid = base_tri.materialid; + tri.translation = base_tri.translation; + tri.rotation = base_tri.rotation; + tri.scale = base_tri.scale; + tri.transform = base_tri.transform; + tri.invTranspose = base_tri.invTranspose; + tri.inverseTransform = base_tri.inverseTransform; + + int idx_j0 = indices[3 * j + 0].vertex_index; + int idx_j1 = indices[3 * j + 1].vertex_index; + int idx_j2 = indices[3 * j + 2].vertex_index; + glm::vec3 p0(positions[idx_j0 * 3], positions[idx_j0 * 3 + 1], positions[idx_j0 * 3 + 2]); + glm::vec3 p1(positions[idx_j1 * 3], positions[idx_j1 * 3 + 1], positions[idx_j1 * 3 + 2]); + glm::vec3 p2(positions[idx_j2 * 3], positions[idx_j2 * 3 + 1], positions[idx_j2 * 3 + 2]); + + tri.positions[0] = p0; + tri.positions[1] = p1; + tri.positions[2] = p2; + + idx_j0 = indices[3 * j + 0].normal_index; + idx_j1 = indices[3 * j + 1].normal_index; + idx_j2 = indices[3 * j + 2].normal_index; + //Get Normals Indices + if (normals.size() > 0) //Checking if Normals defined. + { + glm::vec3 n1(normals[idx_j0 * 3], normals[idx_j0 * 3 + 1], normals[idx_j0 * 3 + 2]); + glm::vec3 n2(normals[idx_j1 * 3], normals[idx_j1 * 3 + 1], normals[idx_j1 * 3 + 2]); + glm::vec3 n3(normals[idx_j2 * 3], normals[idx_j2 * 3 + 1], normals[idx_j2 * 3 + 2]); + tri.normals[0] = n1; + tri.normals[1] = n2; + tri.normals[2] = n3; + } + + idx_j0 = indices[3 * j + 0].texcoord_index; + idx_j1 = indices[3 * j + 1].texcoord_index; + idx_j2 = indices[3 * j + 2].texcoord_index; + if (uvs.size() > 0) //Checking if UVs defined. + { + glm::vec2 t1(uvs[idx_j0 * 2], uvs[idx_j0 * 2 + 1]); + glm::vec2 t2(uvs[idx_j1 * 2], uvs[idx_j1 * 2 + 1]); + glm::vec2 t3(uvs[idx_j2 * 2], uvs[idx_j2 * 2 + 1]); + tri.uvs[0] = t1; + tri.uvs[1] = t2; + tri.uvs[2] = t3; + } + + geoms.push_back(tri); + } + } + } +} + int Scene::loadCamera() { cout << "Loading Camera ..." << endl; RenderState &state = this->state; @@ -151,38 +258,89 @@ int Scene::loadCamera() { } int Scene::loadMaterial(string materialid) { - int id = atoi(materialid.c_str()); - if (id != materials.size()) { - cout << "ERROR: MATERIAL ID does not match expected number of materials" << endl; - return -1; - } else { - cout << "Loading Material " << id << "..." << endl; - Material newMaterial; + int id = atoi(materialid.c_str()); + if (id != materials.size()) { + cout << "ERROR: MATERIAL ID does not match expected number of materials" << endl; + return -1; + } + else { + cout << "Loading Material " << id << "..." << endl; + Material newMaterial; - //load static properties - for (int i = 0; i < 7; i++) { - string line; - utilityCore::safeGetline(fp_in, line); - vector tokens = utilityCore::tokenizeString(line); - if (strcmp(tokens[0].c_str(), "RGB") == 0) { - glm::vec3 color( atof(tokens[1].c_str()), atof(tokens[2].c_str()), atof(tokens[3].c_str()) ); - newMaterial.color = color; - } else if (strcmp(tokens[0].c_str(), "SPECEX") == 0) { - newMaterial.specular.exponent = atof(tokens[1].c_str()); - } else if (strcmp(tokens[0].c_str(), "SPECRGB") == 0) { - glm::vec3 specColor(atof(tokens[1].c_str()), atof(tokens[2].c_str()), atof(tokens[3].c_str())); - newMaterial.specular.color = specColor; - } else if (strcmp(tokens[0].c_str(), "REFL") == 0) { - newMaterial.hasReflective = atof(tokens[1].c_str()); - } else if (strcmp(tokens[0].c_str(), "REFR") == 0) { - newMaterial.hasRefractive = atof(tokens[1].c_str()); - } else if (strcmp(tokens[0].c_str(), "REFRIOR") == 0) { - newMaterial.indexOfRefraction = atof(tokens[1].c_str()); - } else if (strcmp(tokens[0].c_str(), "EMITTANCE") == 0) { - newMaterial.emittance = atof(tokens[1].c_str()); - } - } - materials.push_back(newMaterial); - return 1; - } + //load static properties + for (int i = 0; i < 8; i++) { + string line; + utilityCore::safeGetline(fp_in, line); + vector tokens = utilityCore::tokenizeString(line); + if (strcmp(tokens[0].c_str(), "RGB") == 0) { + glm::vec3 color(atof(tokens[1].c_str()), atof(tokens[2].c_str()), atof(tokens[3].c_str())); + newMaterial.color = color; + } + else if (strcmp(tokens[0].c_str(), "SPECEX") == 0) { + newMaterial.specular.exponent = atof(tokens[1].c_str()); + } + else if (strcmp(tokens[0].c_str(), "SPECRGB") == 0) { + glm::vec3 specColor(atof(tokens[1].c_str()), atof(tokens[2].c_str()), atof(tokens[3].c_str())); + newMaterial.specular.color = specColor; + } + else if (strcmp(tokens[0].c_str(), "REFL") == 0) { + newMaterial.hasReflective = atof(tokens[1].c_str()); + } + else if (strcmp(tokens[0].c_str(), "REFR") == 0) { + newMaterial.hasRefractive = atof(tokens[1].c_str()); + } + else if (strcmp(tokens[0].c_str(), "REFRIOR") == 0) { + newMaterial.indexOfRefraction = atof(tokens[1].c_str()); + } + else if (strcmp(tokens[0].c_str(), "EMITTANCE") == 0) { + newMaterial.emittance = atof(tokens[1].c_str()); + } + else if (strcmp(tokens[0].c_str(), "BSDF") == 0) { + newMaterial.bsdf = atoi(tokens[1].c_str()); + if (newMaterial.bsdf == -1) { light_count++; } + } + } + materials.push_back(newMaterial); + return 1; + } +} + +int Scene::loadEnvironment() { + //To Be Filled + char *name = (char*) malloc(sizeof(char) * FILENAME_MAX); + int dim_x; //Width of Image / Scanline Size + int dim_y; //Height of Image / Scanline Count + int dim_bpp; //Bytes Per Pixel / Pixel Depth + + for (int i = 0; i < 2; i++) { + string line; + utilityCore::safeGetline(fp_in, line); + vector tokens = utilityCore::tokenizeString(line); + if (strcmp(tokens[0].c_str(), "FILENAME") == 0) { + strcpy(name, tokens[1].c_str()); + } + else if (strcmp(tokens[0].c_str(), "DIMENSIONS") == 0) { + dim_x = atoi(tokens[1].c_str()); + dim_y = atoi(tokens[2].c_str()); + dim_bpp = atoi(tokens[3].c_str()); + } + } + environment_dims = glm::ivec3(dim_x, dim_y, dim_bpp); + environment = stbi_load(name, &dim_x, &dim_y, &dim_bpp, 0); + + if (environment == NULL) { + printf("STBI Image Loading failed: %s\n", stbi_failure_reason()); + return -1; + } + + /** + printf("First 10 pixels of texture: \n"); + for (int x = 5*dim_bpp; x < 10*dim_bpp; x+= dim_bpp) { + unsigned char* index = (environment + x); + printf("Pixel is: (%d, %d, %d, %d) ", *(index + 0), *(index + 1), *(index + 2), *(index + 3)); + } + **/ + + return 1; } + diff --git a/src/scene.h b/src/scene.h index f29a917..c20eb53 100644 --- a/src/scene.h +++ b/src/scene.h @@ -7,6 +7,9 @@ #include "glm/glm.hpp" #include "utilities.h" #include "sceneStructs.h" +#include +#include +#include using namespace std; @@ -15,8 +18,15 @@ class Scene { ifstream fp_in; int loadMaterial(string materialid); int loadGeom(string objectid); + void loadOBJ(Geom& mesh, string& filename); int loadCamera(); + int loadEnvironment(); + public: + int light_count = 0; + unsigned char* environment = NULL; + glm::ivec3 environment_dims = glm::ivec3(NULL); + Scene(string filename); ~Scene(); diff --git a/src/sceneStructs.h b/src/sceneStructs.h index b38b820..046a9e7 100644 --- a/src/sceneStructs.h +++ b/src/sceneStructs.h @@ -7,9 +7,13 @@ #define BACKGROUND_COLOR (glm::vec3(0.0f)) +#define ETA_REFRACT 1.33f + enum GeomType { SPHERE, CUBE, + PLANE, + TRIANGLE }; struct Ray { @@ -26,6 +30,10 @@ struct Geom { glm::mat4 transform; glm::mat4 inverseTransform; glm::mat4 invTranspose; + //If Is Mesh + glm::vec3 positions[3]; + glm::vec3 normals[3]; + glm::vec2 uvs[3]; }; struct Material { @@ -38,6 +46,7 @@ struct Material { float hasRefractive; float indexOfRefraction; float emittance; + int bsdf; }; struct Camera { @@ -64,6 +73,8 @@ struct PathSegment { glm::vec3 color; int pixelIndex; int remainingBounces; + bool specularBounce; + glm::vec3 throughput; }; // Use with a corresponding PathSegment to do: @@ -72,5 +83,6 @@ struct PathSegment { struct ShadeableIntersection { float t; glm::vec3 surfaceNormal; + glm::vec3 point; int materialId; }; diff --git a/src/tiny_obj_loader.cc b/src/tiny_obj_loader.cc new file mode 100644 index 0000000..e57d044 --- /dev/null +++ b/src/tiny_obj_loader.cc @@ -0,0 +1,2 @@ +#define TINYOBJLOADER_IMPLEMENTATION +#include "tiny_obj_loader.h" diff --git a/src/tiny_obj_loader.h b/src/tiny_obj_loader.h new file mode 100644 index 0000000..6f0515b --- /dev/null +++ b/src/tiny_obj_loader.h @@ -0,0 +1,2063 @@ +/* +The MIT License (MIT) + +Copyright (c) 2012-2017 Syoyo Fujita and many contributors. + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in +all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +THE SOFTWARE. +*/ + +// +// version 1.0.8 : Fix parsing `g` tag just after `usemtl`(#138) +// version 1.0.7 : Support multiple tex options(#126) +// version 1.0.6 : Add TINYOBJLOADER_USE_DOUBLE option(#124) +// version 1.0.5 : Ignore `Tr` when `d` exists in MTL(#43) +// version 1.0.4 : Support multiple filenames for 'mtllib'(#112) +// version 1.0.3 : Support parsing texture options(#85) +// version 1.0.2 : Improve parsing speed by about a factor of 2 for large +// files(#105) +// version 1.0.1 : Fixes a shape is lost if obj ends with a 'usemtl'(#104) +// version 1.0.0 : Change data structure. Change license from BSD to MIT. +// + +// +// Use this in *one* .cc +// #define TINYOBJLOADER_IMPLEMENTATION +// #include "tiny_obj_loader.h" +// + +#ifndef TINY_OBJ_LOADER_H_ +#define TINY_OBJ_LOADER_H_ + +#include +#include +#include + +namespace tinyobj { + +// https://en.wikipedia.org/wiki/Wavefront_.obj_file says ... +// +// -blendu on | off # set horizontal texture blending +// (default on) +// -blendv on | off # set vertical texture blending +// (default on) +// -boost real_value # boost mip-map sharpness +// -mm base_value gain_value # modify texture map values (default +// 0 1) +// # base_value = brightness, +// gain_value = contrast +// -o u [v [w]] # Origin offset (default +// 0 0 0) +// -s u [v [w]] # Scale (default +// 1 1 1) +// -t u [v [w]] # Turbulence (default +// 0 0 0) +// -texres resolution # texture resolution to create +// -clamp on | off # only render texels in the clamped +// 0-1 range (default off) +// # When unclamped, textures are +// repeated across a surface, +// # when clamped, only texels which +// fall within the 0-1 +// # range are rendered. +// -bm mult_value # bump multiplier (for bump maps +// only) +// +// -imfchan r | g | b | m | l | z # specifies which channel of the file +// is used to +// # create a scalar or bump texture. +// r:red, g:green, +// # b:blue, m:matte, l:luminance, +// z:z-depth.. +// # (the default for bump is 'l' and +// for decal is 'm') +// bump -imfchan r bumpmap.tga # says to use the red channel of +// bumpmap.tga as the bumpmap +// +// For reflection maps... +// +// -type sphere # specifies a sphere for a "refl" +// reflection map +// -type cube_top | cube_bottom | # when using a cube map, the texture +// file for each +// cube_front | cube_back | # side of the cube is specified +// separately +// cube_left | cube_right + +#ifdef TINYOBJLOADER_USE_DOUBLE +//#pragma message "using double" +typedef double real_t; +#else +//#pragma message "using float" +typedef float real_t; +#endif + +typedef enum { + TEXTURE_TYPE_NONE, // default + TEXTURE_TYPE_SPHERE, + TEXTURE_TYPE_CUBE_TOP, + TEXTURE_TYPE_CUBE_BOTTOM, + TEXTURE_TYPE_CUBE_FRONT, + TEXTURE_TYPE_CUBE_BACK, + TEXTURE_TYPE_CUBE_LEFT, + TEXTURE_TYPE_CUBE_RIGHT +} texture_type_t; + +typedef struct { + texture_type_t type; // -type (default TEXTURE_TYPE_NONE) + real_t sharpness; // -boost (default 1.0?) + real_t brightness; // base_value in -mm option (default 0) + real_t contrast; // gain_value in -mm option (default 1) + real_t origin_offset[3]; // -o u [v [w]] (default 0 0 0) + real_t scale[3]; // -s u [v [w]] (default 1 1 1) + real_t turbulence[3]; // -t u [v [w]] (default 0 0 0) + // int texture_resolution; // -texres resolution (default = ?) TODO + bool clamp; // -clamp (default false) + char imfchan; // -imfchan (the default for bump is 'l' and for decal is 'm') + bool blendu; // -blendu (default on) + bool blendv; // -blendv (default on) + real_t bump_multiplier; // -bm (for bump maps only, default 1.0) +} texture_option_t; + +typedef struct { + std::string name; + + real_t ambient[3]; + real_t diffuse[3]; + real_t specular[3]; + real_t transmittance[3]; + real_t emission[3]; + real_t shininess; + real_t ior; // index of refraction + real_t dissolve; // 1 == opaque; 0 == fully transparent + // illumination model (see http://www.fileformat.info/format/material/) + int illum; + + int dummy; // Suppress padding warning. + + std::string ambient_texname; // map_Ka + std::string diffuse_texname; // map_Kd + std::string specular_texname; // map_Ks + std::string specular_highlight_texname; // map_Ns + std::string bump_texname; // map_bump, map_Bump, bump + std::string displacement_texname; // disp + std::string alpha_texname; // map_d + std::string reflection_texname; // refl + + texture_option_t ambient_texopt; + texture_option_t diffuse_texopt; + texture_option_t specular_texopt; + texture_option_t specular_highlight_texopt; + texture_option_t bump_texopt; + texture_option_t displacement_texopt; + texture_option_t alpha_texopt; + texture_option_t reflection_texopt; + + // PBR extension + // http://exocortex.com/blog/extending_wavefront_mtl_to_support_pbr + real_t roughness; // [0, 1] default 0 + real_t metallic; // [0, 1] default 0 + real_t sheen; // [0, 1] default 0 + real_t clearcoat_thickness; // [0, 1] default 0 + real_t clearcoat_roughness; // [0, 1] default 0 + real_t anisotropy; // aniso. [0, 1] default 0 + real_t anisotropy_rotation; // anisor. [0, 1] default 0 + real_t pad0; + std::string roughness_texname; // map_Pr + std::string metallic_texname; // map_Pm + std::string sheen_texname; // map_Ps + std::string emissive_texname; // map_Ke + std::string normal_texname; // norm. For normal mapping. + + texture_option_t roughness_texopt; + texture_option_t metallic_texopt; + texture_option_t sheen_texopt; + texture_option_t emissive_texopt; + texture_option_t normal_texopt; + + int pad2; + + std::map unknown_parameter; +} material_t; + +typedef struct { + std::string name; + + std::vector intValues; + std::vector floatValues; + std::vector stringValues; +} tag_t; + +// Index struct to support different indices for vtx/normal/texcoord. +// -1 means not used. +typedef struct { + int vertex_index; + int normal_index; + int texcoord_index; +} index_t; + +typedef struct { + std::vector indices; + std::vector num_face_vertices; // The number of vertices per + // face. 3 = polygon, 4 = quad, + // ... Up to 255. + std::vector material_ids; // per-face material ID + std::vector tags; // SubD tag +} mesh_t; + +typedef struct { + std::string name; + mesh_t mesh; +} shape_t; + +// Vertex attributes +typedef struct { + std::vector vertices; // 'v' + std::vector normals; // 'vn' + std::vector texcoords; // 'vt' +} attrib_t; + +typedef struct callback_t_ { + // W is optional and set to 1 if there is no `w` item in `v` line + void (*vertex_cb)(void *user_data, real_t x, real_t y, real_t z, real_t w); + void (*normal_cb)(void *user_data, real_t x, real_t y, real_t z); + + // y and z are optional and set to 0 if there is no `y` and/or `z` item(s) in + // `vt` line. + void (*texcoord_cb)(void *user_data, real_t x, real_t y, real_t z); + + // called per 'f' line. num_indices is the number of face indices(e.g. 3 for + // triangle, 4 for quad) + // 0 will be passed for undefined index in index_t members. + void (*index_cb)(void *user_data, index_t *indices, int num_indices); + // `name` material name, `material_id` = the array index of material_t[]. -1 + // if + // a material not found in .mtl + void (*usemtl_cb)(void *user_data, const char *name, int material_id); + // `materials` = parsed material data. + void (*mtllib_cb)(void *user_data, const material_t *materials, + int num_materials); + // There may be multiple group names + void (*group_cb)(void *user_data, const char **names, int num_names); + void (*object_cb)(void *user_data, const char *name); + + callback_t_() + : vertex_cb(NULL), + normal_cb(NULL), + texcoord_cb(NULL), + index_cb(NULL), + usemtl_cb(NULL), + mtllib_cb(NULL), + group_cb(NULL), + object_cb(NULL) {} +} callback_t; + +class MaterialReader { + public: + MaterialReader() {} + virtual ~MaterialReader(); + + virtual bool operator()(const std::string &matId, + std::vector *materials, + std::map *matMap, + std::string *err) = 0; +}; + +class MaterialFileReader : public MaterialReader { + public: + explicit MaterialFileReader(const std::string &mtl_basedir) + : m_mtlBaseDir(mtl_basedir) {} + virtual ~MaterialFileReader() {} + virtual bool operator()(const std::string &matId, + std::vector *materials, + std::map *matMap, std::string *err); + + private: + std::string m_mtlBaseDir; +}; + +class MaterialStreamReader : public MaterialReader { + public: + explicit MaterialStreamReader(std::istream &inStream) + : m_inStream(inStream) {} + virtual ~MaterialStreamReader() {} + virtual bool operator()(const std::string &matId, + std::vector *materials, + std::map *matMap, std::string *err); + + private: + std::istream &m_inStream; +}; + +/// Loads .obj from a file. +/// 'attrib', 'shapes' and 'materials' will be filled with parsed shape data +/// 'shapes' will be filled with parsed shape data +/// Returns true when loading .obj become success. +/// Returns warning and error message into `err` +/// 'mtl_basedir' is optional, and used for base directory for .mtl file. +/// In default(`NULL'), .mtl file is searched from an application's working +/// directory. +/// 'triangulate' is optional, and used whether triangulate polygon face in .obj +/// or not. +bool LoadObj(attrib_t *attrib, std::vector *shapes, + std::vector *materials, std::string *err, + const char *filename, const char *mtl_basedir = NULL, + bool triangulate = true); + +/// Loads .obj from a file with custom user callback. +/// .mtl is loaded as usual and parsed material_t data will be passed to +/// `callback.mtllib_cb`. +/// Returns true when loading .obj/.mtl become success. +/// Returns warning and error message into `err` +/// See `examples/callback_api/` for how to use this function. +bool LoadObjWithCallback(std::istream &inStream, const callback_t &callback, + void *user_data = NULL, + MaterialReader *readMatFn = NULL, + std::string *err = NULL); + +/// Loads object from a std::istream, uses GetMtlIStreamFn to retrieve +/// std::istream for materials. +/// Returns true when loading .obj become success. +/// Returns warning and error message into `err` +bool LoadObj(attrib_t *attrib, std::vector *shapes, + std::vector *materials, std::string *err, + std::istream *inStream, MaterialReader *readMatFn = NULL, + bool triangulate = true); + +/// Loads materials into std::map +void LoadMtl(std::map *material_map, + std::vector *materials, std::istream *inStream, + std::string *warning); + +} // namespace tinyobj + +#endif // TINY_OBJ_LOADER_H_ + +#ifdef TINYOBJLOADER_IMPLEMENTATION +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +namespace tinyobj { + +MaterialReader::~MaterialReader() {} + +struct vertex_index { + int v_idx, vt_idx, vn_idx; + vertex_index() : v_idx(-1), vt_idx(-1), vn_idx(-1) {} + explicit vertex_index(int idx) : v_idx(idx), vt_idx(idx), vn_idx(idx) {} + vertex_index(int vidx, int vtidx, int vnidx) + : v_idx(vidx), vt_idx(vtidx), vn_idx(vnidx) {} +}; + +struct tag_sizes { + tag_sizes() : num_ints(0), num_reals(0), num_strings(0) {} + int num_ints; + int num_reals; + int num_strings; +}; + +struct obj_shape { + std::vector v; + std::vector vn; + std::vector vt; +}; + +// See +// http://stackoverflow.com/questions/6089231/getting-std-ifstream-to-handle-lf-cr-and-crlf +static std::istream &safeGetline(std::istream &is, std::string &t) { + t.clear(); + + // The characters in the stream are read one-by-one using a std::streambuf. + // That is faster than reading them one-by-one using the std::istream. + // Code that uses streambuf this way must be guarded by a sentry object. + // The sentry object performs various tasks, + // such as thread synchronization and updating the stream state. + + std::istream::sentry se(is, true); + std::streambuf *sb = is.rdbuf(); + + if (se) { + for (;;) { + int c = sb->sbumpc(); + switch (c) { + case '\n': + return is; + case '\r': + if (sb->sgetc() == '\n') sb->sbumpc(); + return is; + case EOF: + // Also handle the case when the last line has no line ending + if (t.empty()) is.setstate(std::ios::eofbit); + return is; + default: + t += static_cast(c); + } + } + } + + return is; +} + +#define IS_SPACE(x) (((x) == ' ') || ((x) == '\t')) +#define IS_DIGIT(x) \ + (static_cast((x) - '0') < static_cast(10)) +#define IS_NEW_LINE(x) (((x) == '\r') || ((x) == '\n') || ((x) == '\0')) + +// Make index zero-base, and also support relative index. +static inline bool fixIndex(int idx, int n, int *ret) { + if (!ret) { + return false; + } + + if (idx > 0) { + (*ret) = idx - 1; + return true; + } + + if (idx == 0) { + // zero is not allowed according to the spec. + return false; + } + + if (idx < 0) { + (*ret) = n + idx; // negative value = relative + return true; + } + + return false; // never reach here. +} + +static inline std::string parseString(const char **token) { + std::string s; + (*token) += strspn((*token), " \t"); + size_t e = strcspn((*token), " \t\r"); + s = std::string((*token), &(*token)[e]); + (*token) += e; + return s; +} + +static inline int parseInt(const char **token) { + (*token) += strspn((*token), " \t"); + int i = atoi((*token)); + (*token) += strcspn((*token), " \t\r"); + return i; +} + +// Tries to parse a floating point number located at s. +// +// s_end should be a location in the string where reading should absolutely +// stop. For example at the end of the string, to prevent buffer overflows. +// +// Parses the following EBNF grammar: +// sign = "+" | "-" ; +// END = ? anything not in digit ? +// digit = "0" | "1" | "2" | "3" | "4" | "5" | "6" | "7" | "8" | "9" ; +// integer = [sign] , digit , {digit} ; +// decimal = integer , ["." , integer] ; +// float = ( decimal , END ) | ( decimal , ("E" | "e") , integer , END ) ; +// +// Valid strings are for example: +// -0 +3.1417e+2 -0.0E-3 1.0324 -1.41 11e2 +// +// If the parsing is a success, result is set to the parsed value and true +// is returned. +// +// The function is greedy and will parse until any of the following happens: +// - a non-conforming character is encountered. +// - s_end is reached. +// +// The following situations triggers a failure: +// - s >= s_end. +// - parse failure. +// +static bool tryParseDouble(const char *s, const char *s_end, double *result) { + if (s >= s_end) { + return false; + } + + double mantissa = 0.0; + // This exponent is base 2 rather than 10. + // However the exponent we parse is supposed to be one of ten, + // thus we must take care to convert the exponent/and or the + // mantissa to a * 2^E, where a is the mantissa and E is the + // exponent. + // To get the final double we will use ldexp, it requires the + // exponent to be in base 2. + int exponent = 0; + + // NOTE: THESE MUST BE DECLARED HERE SINCE WE ARE NOT ALLOWED + // TO JUMP OVER DEFINITIONS. + char sign = '+'; + char exp_sign = '+'; + char const *curr = s; + + // How many characters were read in a loop. + int read = 0; + // Tells whether a loop terminated due to reaching s_end. + bool end_not_reached = false; + + /* + BEGIN PARSING. + */ + + // Find out what sign we've got. + if (*curr == '+' || *curr == '-') { + sign = *curr; + curr++; + } else if (IS_DIGIT(*curr)) { /* Pass through. */ + } else { + goto fail; + } + + // Read the integer part. + end_not_reached = (curr != s_end); + while (end_not_reached && IS_DIGIT(*curr)) { + mantissa *= 10; + mantissa += static_cast(*curr - 0x30); + curr++; + read++; + end_not_reached = (curr != s_end); + } + + // We must make sure we actually got something. + if (read == 0) goto fail; + // We allow numbers of form "#", "###" etc. + if (!end_not_reached) goto assemble; + + // Read the decimal part. + if (*curr == '.') { + curr++; + read = 1; + end_not_reached = (curr != s_end); + while (end_not_reached && IS_DIGIT(*curr)) { + static const double pow_lut[] = { + 1.0, 0.1, 0.01, 0.001, 0.0001, 0.00001, 0.000001, 0.0000001, + }; + const int lut_entries = sizeof pow_lut / sizeof pow_lut[0]; + + // NOTE: Don't use powf here, it will absolutely murder precision. + mantissa += static_cast(*curr - 0x30) * + (read < lut_entries ? pow_lut[read] : std::pow(10.0, -read)); + read++; + curr++; + end_not_reached = (curr != s_end); + } + } else if (*curr == 'e' || *curr == 'E') { + } else { + goto assemble; + } + + if (!end_not_reached) goto assemble; + + // Read the exponent part. + if (*curr == 'e' || *curr == 'E') { + curr++; + // Figure out if a sign is present and if it is. + end_not_reached = (curr != s_end); + if (end_not_reached && (*curr == '+' || *curr == '-')) { + exp_sign = *curr; + curr++; + } else if (IS_DIGIT(*curr)) { /* Pass through. */ + } else { + // Empty E is not allowed. + goto fail; + } + + read = 0; + end_not_reached = (curr != s_end); + while (end_not_reached && IS_DIGIT(*curr)) { + exponent *= 10; + exponent += static_cast(*curr - 0x30); + curr++; + read++; + end_not_reached = (curr != s_end); + } + exponent *= (exp_sign == '+' ? 1 : -1); + if (read == 0) goto fail; + } + +assemble: + *result = (sign == '+' ? 1 : -1) * + (exponent ? std::ldexp(mantissa * std::pow(5.0, exponent), exponent) + : mantissa); + return true; +fail: + return false; +} + +static inline real_t parseReal(const char **token, double default_value = 0.0) { + (*token) += strspn((*token), " \t"); + const char *end = (*token) + strcspn((*token), " \t\r"); + double val = default_value; + tryParseDouble((*token), end, &val); + real_t f = static_cast(val); + (*token) = end; + return f; +} + +static inline void parseReal2(real_t *x, real_t *y, const char **token, + const double default_x = 0.0, + const double default_y = 0.0) { + (*x) = parseReal(token, default_x); + (*y) = parseReal(token, default_y); +} + +static inline void parseReal3(real_t *x, real_t *y, real_t *z, + const char **token, const double default_x = 0.0, + const double default_y = 0.0, + const double default_z = 0.0) { + (*x) = parseReal(token, default_x); + (*y) = parseReal(token, default_y); + (*z) = parseReal(token, default_z); +} + +static inline void parseV(real_t *x, real_t *y, real_t *z, real_t *w, + const char **token, const double default_x = 0.0, + const double default_y = 0.0, + const double default_z = 0.0, + const double default_w = 1.0) { + (*x) = parseReal(token, default_x); + (*y) = parseReal(token, default_y); + (*z) = parseReal(token, default_z); + (*w) = parseReal(token, default_w); +} + +static inline bool parseOnOff(const char **token, bool default_value = true) { + (*token) += strspn((*token), " \t"); + const char *end = (*token) + strcspn((*token), " \t\r"); + + bool ret = default_value; + if ((0 == strncmp((*token), "on", 2))) { + ret = true; + } else if ((0 == strncmp((*token), "off", 3))) { + ret = false; + } + + (*token) = end; + return ret; +} + +static inline texture_type_t parseTextureType( + const char **token, texture_type_t default_value = TEXTURE_TYPE_NONE) { + (*token) += strspn((*token), " \t"); + const char *end = (*token) + strcspn((*token), " \t\r"); + texture_type_t ty = default_value; + + if ((0 == strncmp((*token), "cube_top", strlen("cube_top")))) { + ty = TEXTURE_TYPE_CUBE_TOP; + } else if ((0 == strncmp((*token), "cube_bottom", strlen("cube_bottom")))) { + ty = TEXTURE_TYPE_CUBE_BOTTOM; + } else if ((0 == strncmp((*token), "cube_left", strlen("cube_left")))) { + ty = TEXTURE_TYPE_CUBE_LEFT; + } else if ((0 == strncmp((*token), "cube_right", strlen("cube_right")))) { + ty = TEXTURE_TYPE_CUBE_RIGHT; + } else if ((0 == strncmp((*token), "cube_front", strlen("cube_front")))) { + ty = TEXTURE_TYPE_CUBE_FRONT; + } else if ((0 == strncmp((*token), "cube_back", strlen("cube_back")))) { + ty = TEXTURE_TYPE_CUBE_BACK; + } else if ((0 == strncmp((*token), "sphere", strlen("sphere")))) { + ty = TEXTURE_TYPE_SPHERE; + } + + (*token) = end; + return ty; +} + +static tag_sizes parseTagTriple(const char **token) { + tag_sizes ts; + + ts.num_ints = atoi((*token)); + (*token) += strcspn((*token), "/ \t\r"); + if ((*token)[0] != '/') { + return ts; + } + (*token)++; + + ts.num_reals = atoi((*token)); + (*token) += strcspn((*token), "/ \t\r"); + if ((*token)[0] != '/') { + return ts; + } + (*token)++; + + ts.num_strings = atoi((*token)); + (*token) += strcspn((*token), "/ \t\r") + 1; + + return ts; +} + +// Parse triples with index offsets: i, i/j/k, i//k, i/j +static bool parseTriple(const char **token, int vsize, int vnsize, int vtsize, + vertex_index *ret) { + if (!ret) { + return false; + } + + vertex_index vi(-1); + + if (!fixIndex(atoi((*token)), vsize, &(vi.v_idx))) { + return false; + } + + (*token) += strcspn((*token), "/ \t\r"); + if ((*token)[0] != '/') { + (*ret) = vi; + return true; + } + (*token)++; + + // i//k + if ((*token)[0] == '/') { + (*token)++; + if (!fixIndex(atoi((*token)), vnsize, &(vi.vn_idx))) { + return false; + } + (*token) += strcspn((*token), "/ \t\r"); + (*ret) = vi; + return true; + } + + // i/j/k or i/j + if (!fixIndex(atoi((*token)), vtsize, &(vi.vt_idx))) { + return false; + } + + (*token) += strcspn((*token), "/ \t\r"); + if ((*token)[0] != '/') { + (*ret) = vi; + return true; + } + + // i/j/k + (*token)++; // skip '/' + if (!fixIndex(atoi((*token)), vnsize, &(vi.vn_idx))) { + return false; + } + (*token) += strcspn((*token), "/ \t\r"); + + (*ret) = vi; + + return true; +} + +// Parse raw triples: i, i/j/k, i//k, i/j +static vertex_index parseRawTriple(const char **token) { + vertex_index vi(static_cast(0)); // 0 is an invalid index in OBJ + + vi.v_idx = atoi((*token)); + (*token) += strcspn((*token), "/ \t\r"); + if ((*token)[0] != '/') { + return vi; + } + (*token)++; + + // i//k + if ((*token)[0] == '/') { + (*token)++; + vi.vn_idx = atoi((*token)); + (*token) += strcspn((*token), "/ \t\r"); + return vi; + } + + // i/j/k or i/j + vi.vt_idx = atoi((*token)); + (*token) += strcspn((*token), "/ \t\r"); + if ((*token)[0] != '/') { + return vi; + } + + // i/j/k + (*token)++; // skip '/' + vi.vn_idx = atoi((*token)); + (*token) += strcspn((*token), "/ \t\r"); + return vi; +} + +static bool ParseTextureNameAndOption(std::string *texname, + texture_option_t *texopt, + const char *linebuf, const bool is_bump) { + // @todo { write more robust lexer and parser. } + bool found_texname = false; + std::string texture_name; + + // Fill with default value for texopt. + if (is_bump) { + texopt->imfchan = 'l'; + } else { + texopt->imfchan = 'm'; + } + texopt->bump_multiplier = 1.0f; + texopt->clamp = false; + texopt->blendu = true; + texopt->blendv = true; + texopt->sharpness = 1.0f; + texopt->brightness = 0.0f; + texopt->contrast = 1.0f; + texopt->origin_offset[0] = 0.0f; + texopt->origin_offset[1] = 0.0f; + texopt->origin_offset[2] = 0.0f; + texopt->scale[0] = 1.0f; + texopt->scale[1] = 1.0f; + texopt->scale[2] = 1.0f; + texopt->turbulence[0] = 0.0f; + texopt->turbulence[1] = 0.0f; + texopt->turbulence[2] = 0.0f; + texopt->type = TEXTURE_TYPE_NONE; + + const char *token = linebuf; // Assume line ends with NULL + + while (!IS_NEW_LINE((*token))) { + token += strspn(token, " \t"); // skip space + if ((0 == strncmp(token, "-blendu", 7)) && IS_SPACE((token[7]))) { + token += 8; + texopt->blendu = parseOnOff(&token, /* default */ true); + } else if ((0 == strncmp(token, "-blendv", 7)) && IS_SPACE((token[7]))) { + token += 8; + texopt->blendv = parseOnOff(&token, /* default */ true); + } else if ((0 == strncmp(token, "-clamp", 6)) && IS_SPACE((token[6]))) { + token += 7; + texopt->clamp = parseOnOff(&token, /* default */ true); + } else if ((0 == strncmp(token, "-boost", 6)) && IS_SPACE((token[6]))) { + token += 7; + texopt->sharpness = parseReal(&token, 1.0); + } else if ((0 == strncmp(token, "-bm", 3)) && IS_SPACE((token[3]))) { + token += 4; + texopt->bump_multiplier = parseReal(&token, 1.0); + } else if ((0 == strncmp(token, "-o", 2)) && IS_SPACE((token[2]))) { + token += 3; + parseReal3(&(texopt->origin_offset[0]), &(texopt->origin_offset[1]), + &(texopt->origin_offset[2]), &token); + } else if ((0 == strncmp(token, "-s", 2)) && IS_SPACE((token[2]))) { + token += 3; + parseReal3(&(texopt->scale[0]), &(texopt->scale[1]), &(texopt->scale[2]), + &token, 1.0, 1.0, 1.0); + } else if ((0 == strncmp(token, "-t", 2)) && IS_SPACE((token[2]))) { + token += 3; + parseReal3(&(texopt->turbulence[0]), &(texopt->turbulence[1]), + &(texopt->turbulence[2]), &token); + } else if ((0 == strncmp(token, "-type", 5)) && IS_SPACE((token[5]))) { + token += 5; + texopt->type = parseTextureType((&token), TEXTURE_TYPE_NONE); + } else if ((0 == strncmp(token, "-imfchan", 8)) && IS_SPACE((token[8]))) { + token += 9; + token += strspn(token, " \t"); + const char *end = token + strcspn(token, " \t\r"); + if ((end - token) == 1) { // Assume one char for -imfchan + texopt->imfchan = (*token); + } + token = end; + } else if ((0 == strncmp(token, "-mm", 3)) && IS_SPACE((token[3]))) { + token += 4; + parseReal2(&(texopt->brightness), &(texopt->contrast), &token, 0.0, 1.0); + } else { + // Assume texture filename + size_t len = strcspn(token, " \t\r"); // untile next space + texture_name = std::string(token, token + len); + token += len; + + token += strspn(token, " \t"); // skip space + + found_texname = true; + } + } + + if (found_texname) { + (*texname) = texture_name; + return true; + } else { + return false; + } +} + +static void InitMaterial(material_t *material) { + material->name = ""; + material->ambient_texname = ""; + material->diffuse_texname = ""; + material->specular_texname = ""; + material->specular_highlight_texname = ""; + material->bump_texname = ""; + material->displacement_texname = ""; + material->reflection_texname = ""; + material->alpha_texname = ""; + for (int i = 0; i < 3; i++) { + material->ambient[i] = 0.f; + material->diffuse[i] = 0.f; + material->specular[i] = 0.f; + material->transmittance[i] = 0.f; + material->emission[i] = 0.f; + } + material->illum = 0; + material->dissolve = 1.f; + material->shininess = 1.f; + material->ior = 1.f; + + material->roughness = 0.f; + material->metallic = 0.f; + material->sheen = 0.f; + material->clearcoat_thickness = 0.f; + material->clearcoat_roughness = 0.f; + material->anisotropy_rotation = 0.f; + material->anisotropy = 0.f; + material->roughness_texname = ""; + material->metallic_texname = ""; + material->sheen_texname = ""; + material->emissive_texname = ""; + material->normal_texname = ""; + + material->unknown_parameter.clear(); +} + +static bool exportFaceGroupToShape( + shape_t *shape, const std::vector > &faceGroup, + const std::vector &tags, const int material_id, + const std::string &name, bool triangulate) { + if (faceGroup.empty()) { + return false; + } + + // Flatten vertices and indices + for (size_t i = 0; i < faceGroup.size(); i++) { + const std::vector &face = faceGroup[i]; + + vertex_index i0 = face[0]; + vertex_index i1(-1); + vertex_index i2 = face[1]; + + size_t npolys = face.size(); + + if (triangulate) { + // Polygon -> triangle fan conversion + for (size_t k = 2; k < npolys; k++) { + i1 = i2; + i2 = face[k]; + + index_t idx0, idx1, idx2; + idx0.vertex_index = i0.v_idx; + idx0.normal_index = i0.vn_idx; + idx0.texcoord_index = i0.vt_idx; + idx1.vertex_index = i1.v_idx; + idx1.normal_index = i1.vn_idx; + idx1.texcoord_index = i1.vt_idx; + idx2.vertex_index = i2.v_idx; + idx2.normal_index = i2.vn_idx; + idx2.texcoord_index = i2.vt_idx; + + shape->mesh.indices.push_back(idx0); + shape->mesh.indices.push_back(idx1); + shape->mesh.indices.push_back(idx2); + + shape->mesh.num_face_vertices.push_back(3); + shape->mesh.material_ids.push_back(material_id); + } + } else { + for (size_t k = 0; k < npolys; k++) { + index_t idx; + idx.vertex_index = face[k].v_idx; + idx.normal_index = face[k].vn_idx; + idx.texcoord_index = face[k].vt_idx; + shape->mesh.indices.push_back(idx); + } + + shape->mesh.num_face_vertices.push_back( + static_cast(npolys)); + shape->mesh.material_ids.push_back(material_id); // per face + } + } + + shape->name = name; + shape->mesh.tags = tags; + + return true; +} + +// Split a string with specified delimiter character. +// http://stackoverflow.com/questions/236129/split-a-string-in-c +static void SplitString(const std::string &s, char delim, + std::vector &elems) { + std::stringstream ss; + ss.str(s); + std::string item; + while (std::getline(ss, item, delim)) { + elems.push_back(item); + } +} + +void LoadMtl(std::map *material_map, + std::vector *materials, std::istream *inStream, + std::string *warning) { + // Create a default material anyway. + material_t material; + InitMaterial(&material); + + // Issue 43. `d` wins against `Tr` since `Tr` is not in the MTL specification. + bool has_d = false; + bool has_tr = false; + + std::stringstream ss; + + std::string linebuf; + while (inStream->peek() != -1) { + safeGetline(*inStream, linebuf); + + // Trim trailing whitespace. + if (linebuf.size() > 0) { + linebuf = linebuf.substr(0, linebuf.find_last_not_of(" \t") + 1); + } + + // Trim newline '\r\n' or '\n' + if (linebuf.size() > 0) { + if (linebuf[linebuf.size() - 1] == '\n') + linebuf.erase(linebuf.size() - 1); + } + if (linebuf.size() > 0) { + if (linebuf[linebuf.size() - 1] == '\r') + linebuf.erase(linebuf.size() - 1); + } + + // Skip if empty line. + if (linebuf.empty()) { + continue; + } + + // Skip leading space. + const char *token = linebuf.c_str(); + token += strspn(token, " \t"); + + assert(token); + if (token[0] == '\0') continue; // empty line + + if (token[0] == '#') continue; // comment line + + // new mtl + if ((0 == strncmp(token, "newmtl", 6)) && IS_SPACE((token[6]))) { + // flush previous material. + if (!material.name.empty()) { + material_map->insert(std::pair( + material.name, static_cast(materials->size()))); + materials->push_back(material); + } + + // initial temporary material + InitMaterial(&material); + + has_d = false; + has_tr = false; + + // set new mtl name + token += 7; + { + std::stringstream sstr; + sstr << token; + material.name = sstr.str(); + } + continue; + } + + // ambient + if (token[0] == 'K' && token[1] == 'a' && IS_SPACE((token[2]))) { + token += 2; + real_t r, g, b; + parseReal3(&r, &g, &b, &token); + material.ambient[0] = r; + material.ambient[1] = g; + material.ambient[2] = b; + continue; + } + + // diffuse + if (token[0] == 'K' && token[1] == 'd' && IS_SPACE((token[2]))) { + token += 2; + real_t r, g, b; + parseReal3(&r, &g, &b, &token); + material.diffuse[0] = r; + material.diffuse[1] = g; + material.diffuse[2] = b; + continue; + } + + // specular + if (token[0] == 'K' && token[1] == 's' && IS_SPACE((token[2]))) { + token += 2; + real_t r, g, b; + parseReal3(&r, &g, &b, &token); + material.specular[0] = r; + material.specular[1] = g; + material.specular[2] = b; + continue; + } + + // transmittance + if ((token[0] == 'K' && token[1] == 't' && IS_SPACE((token[2]))) || + (token[0] == 'T' && token[1] == 'f' && IS_SPACE((token[2])))) { + token += 2; + real_t r, g, b; + parseReal3(&r, &g, &b, &token); + material.transmittance[0] = r; + material.transmittance[1] = g; + material.transmittance[2] = b; + continue; + } + + // ior(index of refraction) + if (token[0] == 'N' && token[1] == 'i' && IS_SPACE((token[2]))) { + token += 2; + material.ior = parseReal(&token); + continue; + } + + // emission + if (token[0] == 'K' && token[1] == 'e' && IS_SPACE(token[2])) { + token += 2; + real_t r, g, b; + parseReal3(&r, &g, &b, &token); + material.emission[0] = r; + material.emission[1] = g; + material.emission[2] = b; + continue; + } + + // shininess + if (token[0] == 'N' && token[1] == 's' && IS_SPACE(token[2])) { + token += 2; + material.shininess = parseReal(&token); + continue; + } + + // illum model + if (0 == strncmp(token, "illum", 5) && IS_SPACE(token[5])) { + token += 6; + material.illum = parseInt(&token); + continue; + } + + // dissolve + if ((token[0] == 'd' && IS_SPACE(token[1]))) { + token += 1; + material.dissolve = parseReal(&token); + + if (has_tr) { + ss << "WARN: Both `d` and `Tr` parameters defined for \"" + << material.name << "\". Use the value of `d` for dissolve." + << std::endl; + } + has_d = true; + continue; + } + if (token[0] == 'T' && token[1] == 'r' && IS_SPACE(token[2])) { + token += 2; + if (has_d) { + // `d` wins. Ignore `Tr` value. + ss << "WARN: Both `d` and `Tr` parameters defined for \"" + << material.name << "\". Use the value of `d` for dissolve." + << std::endl; + } else { + // We invert value of Tr(assume Tr is in range [0, 1]) + // NOTE: Interpretation of Tr is application(exporter) dependent. For + // some application(e.g. 3ds max obj exporter), Tr = d(Issue 43) + material.dissolve = 1.0f - parseReal(&token); + } + has_tr = true; + continue; + } + + // PBR: roughness + if (token[0] == 'P' && token[1] == 'r' && IS_SPACE(token[2])) { + token += 2; + material.roughness = parseReal(&token); + continue; + } + + // PBR: metallic + if (token[0] == 'P' && token[1] == 'm' && IS_SPACE(token[2])) { + token += 2; + material.metallic = parseReal(&token); + continue; + } + + // PBR: sheen + if (token[0] == 'P' && token[1] == 's' && IS_SPACE(token[2])) { + token += 2; + material.sheen = parseReal(&token); + continue; + } + + // PBR: clearcoat thickness + if (token[0] == 'P' && token[1] == 'c' && IS_SPACE(token[2])) { + token += 2; + material.clearcoat_thickness = parseReal(&token); + continue; + } + + // PBR: clearcoat roughness + if ((0 == strncmp(token, "Pcr", 3)) && IS_SPACE(token[3])) { + token += 4; + material.clearcoat_roughness = parseReal(&token); + continue; + } + + // PBR: anisotropy + if ((0 == strncmp(token, "aniso", 5)) && IS_SPACE(token[5])) { + token += 6; + material.anisotropy = parseReal(&token); + continue; + } + + // PBR: anisotropy rotation + if ((0 == strncmp(token, "anisor", 6)) && IS_SPACE(token[6])) { + token += 7; + material.anisotropy_rotation = parseReal(&token); + continue; + } + + // ambient texture + if ((0 == strncmp(token, "map_Ka", 6)) && IS_SPACE(token[6])) { + token += 7; + ParseTextureNameAndOption(&(material.ambient_texname), + &(material.ambient_texopt), token, + /* is_bump */ false); + continue; + } + + // diffuse texture + if ((0 == strncmp(token, "map_Kd", 6)) && IS_SPACE(token[6])) { + token += 7; + ParseTextureNameAndOption(&(material.diffuse_texname), + &(material.diffuse_texopt), token, + /* is_bump */ false); + continue; + } + + // specular texture + if ((0 == strncmp(token, "map_Ks", 6)) && IS_SPACE(token[6])) { + token += 7; + ParseTextureNameAndOption(&(material.specular_texname), + &(material.specular_texopt), token, + /* is_bump */ false); + continue; + } + + // specular highlight texture + if ((0 == strncmp(token, "map_Ns", 6)) && IS_SPACE(token[6])) { + token += 7; + ParseTextureNameAndOption(&(material.specular_highlight_texname), + &(material.specular_highlight_texopt), token, + /* is_bump */ false); + continue; + } + + // bump texture + if ((0 == strncmp(token, "map_bump", 8)) && IS_SPACE(token[8])) { + token += 9; + ParseTextureNameAndOption(&(material.bump_texname), + &(material.bump_texopt), token, + /* is_bump */ true); + continue; + } + + // bump texture + if ((0 == strncmp(token, "map_Bump", 8)) && IS_SPACE(token[8])) { + token += 9; + ParseTextureNameAndOption(&(material.bump_texname), + &(material.bump_texopt), token, + /* is_bump */ true); + continue; + } + + // bump texture + if ((0 == strncmp(token, "bump", 4)) && IS_SPACE(token[4])) { + token += 5; + ParseTextureNameAndOption(&(material.bump_texname), + &(material.bump_texopt), token, + /* is_bump */ true); + continue; + } + + // alpha texture + if ((0 == strncmp(token, "map_d", 5)) && IS_SPACE(token[5])) { + token += 6; + material.alpha_texname = token; + ParseTextureNameAndOption(&(material.alpha_texname), + &(material.alpha_texopt), token, + /* is_bump */ false); + continue; + } + + // displacement texture + if ((0 == strncmp(token, "disp", 4)) && IS_SPACE(token[4])) { + token += 5; + ParseTextureNameAndOption(&(material.displacement_texname), + &(material.displacement_texopt), token, + /* is_bump */ false); + continue; + } + + // reflection map + if ((0 == strncmp(token, "refl", 4)) && IS_SPACE(token[4])) { + token += 5; + ParseTextureNameAndOption(&(material.reflection_texname), + &(material.reflection_texopt), token, + /* is_bump */ false); + continue; + } + + // PBR: roughness texture + if ((0 == strncmp(token, "map_Pr", 6)) && IS_SPACE(token[6])) { + token += 7; + ParseTextureNameAndOption(&(material.roughness_texname), + &(material.roughness_texopt), token, + /* is_bump */ false); + continue; + } + + // PBR: metallic texture + if ((0 == strncmp(token, "map_Pm", 6)) && IS_SPACE(token[6])) { + token += 7; + ParseTextureNameAndOption(&(material.metallic_texname), + &(material.metallic_texopt), token, + /* is_bump */ false); + continue; + } + + // PBR: sheen texture + if ((0 == strncmp(token, "map_Ps", 6)) && IS_SPACE(token[6])) { + token += 7; + ParseTextureNameAndOption(&(material.sheen_texname), + &(material.sheen_texopt), token, + /* is_bump */ false); + continue; + } + + // PBR: emissive texture + if ((0 == strncmp(token, "map_Ke", 6)) && IS_SPACE(token[6])) { + token += 7; + ParseTextureNameAndOption(&(material.emissive_texname), + &(material.emissive_texopt), token, + /* is_bump */ false); + continue; + } + + // PBR: normal map texture + if ((0 == strncmp(token, "norm", 4)) && IS_SPACE(token[4])) { + token += 5; + ParseTextureNameAndOption( + &(material.normal_texname), &(material.normal_texopt), token, + /* is_bump */ false); // @fixme { is_bump will be true? } + continue; + } + + // unknown parameter + const char *_space = strchr(token, ' '); + if (!_space) { + _space = strchr(token, '\t'); + } + if (_space) { + std::ptrdiff_t len = _space - token; + std::string key(token, static_cast(len)); + std::string value = _space + 1; + material.unknown_parameter.insert( + std::pair(key, value)); + } + } + // flush last material. + material_map->insert(std::pair( + material.name, static_cast(materials->size()))); + materials->push_back(material); + + if (warning) { + (*warning) = ss.str(); + } +} + +bool MaterialFileReader::operator()(const std::string &matId, + std::vector *materials, + std::map *matMap, + std::string *err) { + std::string filepath; + + if (!m_mtlBaseDir.empty()) { + filepath = std::string(m_mtlBaseDir) + matId; + } else { + filepath = matId; + } + + std::ifstream matIStream(filepath.c_str()); + if (!matIStream) { + std::stringstream ss; + ss << "WARN: Material file [ " << filepath << " ] not found." << std::endl; + if (err) { + (*err) += ss.str(); + } + return false; + } + + std::string warning; + LoadMtl(matMap, materials, &matIStream, &warning); + + if (!warning.empty()) { + if (err) { + (*err) += warning; + } + } + + return true; +} + +bool MaterialStreamReader::operator()(const std::string &matId, + std::vector *materials, + std::map *matMap, + std::string *err) { + (void)matId; + if (!m_inStream) { + std::stringstream ss; + ss << "WARN: Material stream in error state. " << std::endl; + if (err) { + (*err) += ss.str(); + } + return false; + } + + std::string warning; + LoadMtl(matMap, materials, &m_inStream, &warning); + + if (!warning.empty()) { + if (err) { + (*err) += warning; + } + } + + return true; +} + +bool LoadObj(attrib_t *attrib, std::vector *shapes, + std::vector *materials, std::string *err, + const char *filename, const char *mtl_basedir, bool trianglulate) { + attrib->vertices.clear(); + attrib->normals.clear(); + attrib->texcoords.clear(); + shapes->clear(); + + std::stringstream errss; + + std::ifstream ifs(filename); + if (!ifs) { + errss << "Cannot open file [" << filename << "]" << std::endl; + if (err) { + (*err) = errss.str(); + } + return false; + } + + std::string baseDir; + if (mtl_basedir) { + baseDir = mtl_basedir; + } + MaterialFileReader matFileReader(baseDir); + + return LoadObj(attrib, shapes, materials, err, &ifs, &matFileReader, + trianglulate); +} + +bool LoadObj(attrib_t *attrib, std::vector *shapes, + std::vector *materials, std::string *err, + std::istream *inStream, MaterialReader *readMatFn /*= NULL*/, + bool triangulate) { + std::stringstream errss; + + std::vector v; + std::vector vn; + std::vector vt; + std::vector tags; + std::vector > faceGroup; + std::string name; + + // material + std::map material_map; + int material = -1; + + shape_t shape; + + std::string linebuf; + while (inStream->peek() != -1) { + safeGetline(*inStream, linebuf); + + // Trim newline '\r\n' or '\n' + if (linebuf.size() > 0) { + if (linebuf[linebuf.size() - 1] == '\n') + linebuf.erase(linebuf.size() - 1); + } + if (linebuf.size() > 0) { + if (linebuf[linebuf.size() - 1] == '\r') + linebuf.erase(linebuf.size() - 1); + } + + // Skip if empty line. + if (linebuf.empty()) { + continue; + } + + // Skip leading space. + const char *token = linebuf.c_str(); + token += strspn(token, " \t"); + + assert(token); + if (token[0] == '\0') continue; // empty line + + if (token[0] == '#') continue; // comment line + + // vertex + if (token[0] == 'v' && IS_SPACE((token[1]))) { + token += 2; + real_t x, y, z; + parseReal3(&x, &y, &z, &token); + v.push_back(x); + v.push_back(y); + v.push_back(z); + continue; + } + + // normal + if (token[0] == 'v' && token[1] == 'n' && IS_SPACE((token[2]))) { + token += 3; + real_t x, y, z; + parseReal3(&x, &y, &z, &token); + vn.push_back(x); + vn.push_back(y); + vn.push_back(z); + continue; + } + + // texcoord + if (token[0] == 'v' && token[1] == 't' && IS_SPACE((token[2]))) { + token += 3; + real_t x, y; + parseReal2(&x, &y, &token); + vt.push_back(x); + vt.push_back(y); + continue; + } + + // face + if (token[0] == 'f' && IS_SPACE((token[1]))) { + token += 2; + token += strspn(token, " \t"); + + std::vector face; + face.reserve(3); + + while (!IS_NEW_LINE(token[0])) { + vertex_index vi; + if (!parseTriple(&token, static_cast(v.size() / 3), + static_cast(vn.size() / 3), + static_cast(vt.size() / 2), &vi)) { + if (err) { + (*err) = "Failed parse `f' line(e.g. zero value for face index).\n"; + } + return false; + } + + face.push_back(vi); + size_t n = strspn(token, " \t\r"); + token += n; + } + + // replace with emplace_back + std::move on C++11 + faceGroup.push_back(std::vector()); + faceGroup[faceGroup.size() - 1].swap(face); + + continue; + } + + // use mtl + if ((0 == strncmp(token, "usemtl", 6)) && IS_SPACE((token[6]))) { + token += 7; + std::stringstream ss; + ss << token; + std::string namebuf = ss.str(); + + int newMaterialId = -1; + if (material_map.find(namebuf) != material_map.end()) { + newMaterialId = material_map[namebuf]; + } else { + // { error!! material not found } + } + + if (newMaterialId != material) { + // Create per-face material. Thus we don't add `shape` to `shapes` at + // this time. + // just clear `faceGroup` after `exportFaceGroupToShape()` call. + exportFaceGroupToShape(&shape, faceGroup, tags, material, name, + triangulate); + faceGroup.clear(); + material = newMaterialId; + } + + continue; + } + + // load mtl + if ((0 == strncmp(token, "mtllib", 6)) && IS_SPACE((token[6]))) { + if (readMatFn) { + token += 7; + + std::vector filenames; + SplitString(std::string(token), ' ', filenames); + + if (filenames.empty()) { + if (err) { + (*err) += + "WARN: Looks like empty filename for mtllib. Use default " + "material. \n"; + } + } else { + bool found = false; + for (size_t s = 0; s < filenames.size(); s++) { + std::string err_mtl; + bool ok = (*readMatFn)(filenames[s].c_str(), materials, + &material_map, &err_mtl); + if (err && (!err_mtl.empty())) { + (*err) += err_mtl; // This should be warn message. + } + + if (ok) { + found = true; + break; + } + } + + if (!found) { + if (err) { + (*err) += + "WARN: Failed to load material file(s). Use default " + "material.\n"; + } + } + } + } + + continue; + } + + // group name + if (token[0] == 'g' && IS_SPACE((token[1]))) { + // flush previous face group. + bool ret = exportFaceGroupToShape(&shape, faceGroup, tags, material, name, + triangulate); + (void)ret; // return value not used. + + if (shape.mesh.indices.size() > 0) { + shapes->push_back(shape); + } + + shape = shape_t(); + + // material = -1; + faceGroup.clear(); + + std::vector names; + names.reserve(2); + + while (!IS_NEW_LINE(token[0])) { + std::string str = parseString(&token); + names.push_back(str); + token += strspn(token, " \t\r"); // skip tag + } + + assert(names.size() > 0); + + // names[0] must be 'g', so skip the 0th element. + if (names.size() > 1) { + name = names[1]; + } else { + name = ""; + } + + continue; + } + + // object name + if (token[0] == 'o' && IS_SPACE((token[1]))) { + // flush previous face group. + bool ret = exportFaceGroupToShape(&shape, faceGroup, tags, material, name, + triangulate); + if (ret) { + shapes->push_back(shape); + } + + // material = -1; + faceGroup.clear(); + shape = shape_t(); + + // @todo { multiple object name? } + token += 2; + std::stringstream ss; + ss << token; + name = ss.str(); + + continue; + } + + if (token[0] == 't' && IS_SPACE(token[1])) { + tag_t tag; + + token += 2; + std::stringstream ss; + ss << token; + tag.name = ss.str(); + + token += tag.name.size() + 1; + + tag_sizes ts = parseTagTriple(&token); + + tag.intValues.resize(static_cast(ts.num_ints)); + + for (size_t i = 0; i < static_cast(ts.num_ints); ++i) { + tag.intValues[i] = atoi(token); + token += strcspn(token, "/ \t\r") + 1; + } + + tag.floatValues.resize(static_cast(ts.num_reals)); + for (size_t i = 0; i < static_cast(ts.num_reals); ++i) { + tag.floatValues[i] = parseReal(&token); + token += strcspn(token, "/ \t\r") + 1; + } + + tag.stringValues.resize(static_cast(ts.num_strings)); + for (size_t i = 0; i < static_cast(ts.num_strings); ++i) { + std::stringstream sstr; + sstr << token; + tag.stringValues[i] = sstr.str(); + token += tag.stringValues[i].size() + 1; + } + + tags.push_back(tag); + } + + // Ignore unknown command. + } + + bool ret = exportFaceGroupToShape(&shape, faceGroup, tags, material, name, + triangulate); + // exportFaceGroupToShape return false when `usemtl` is called in the last + // line. + // we also add `shape` to `shapes` when `shape.mesh` has already some + // faces(indices) + if (ret || shape.mesh.indices.size()) { + shapes->push_back(shape); + } + faceGroup.clear(); // for safety + + if (err) { + (*err) += errss.str(); + } + + attrib->vertices.swap(v); + attrib->normals.swap(vn); + attrib->texcoords.swap(vt); + + return true; +} + +bool LoadObjWithCallback(std::istream &inStream, const callback_t &callback, + void *user_data /*= NULL*/, + MaterialReader *readMatFn /*= NULL*/, + std::string *err /*= NULL*/) { + std::stringstream errss; + + // material + std::map material_map; + int material_id = -1; // -1 = invalid + + std::vector indices; + std::vector materials; + std::vector names; + names.reserve(2); + std::string name; + std::vector names_out; + + std::string linebuf; + while (inStream.peek() != -1) { + safeGetline(inStream, linebuf); + + // Trim newline '\r\n' or '\n' + if (linebuf.size() > 0) { + if (linebuf[linebuf.size() - 1] == '\n') + linebuf.erase(linebuf.size() - 1); + } + if (linebuf.size() > 0) { + if (linebuf[linebuf.size() - 1] == '\r') + linebuf.erase(linebuf.size() - 1); + } + + // Skip if empty line. + if (linebuf.empty()) { + continue; + } + + // Skip leading space. + const char *token = linebuf.c_str(); + token += strspn(token, " \t"); + + assert(token); + if (token[0] == '\0') continue; // empty line + + if (token[0] == '#') continue; // comment line + + // vertex + if (token[0] == 'v' && IS_SPACE((token[1]))) { + token += 2; + real_t x, y, z, w; // w is optional. default = 1.0 + parseV(&x, &y, &z, &w, &token); + if (callback.vertex_cb) { + callback.vertex_cb(user_data, x, y, z, w); + } + continue; + } + + // normal + if (token[0] == 'v' && token[1] == 'n' && IS_SPACE((token[2]))) { + token += 3; + real_t x, y, z; + parseReal3(&x, &y, &z, &token); + if (callback.normal_cb) { + callback.normal_cb(user_data, x, y, z); + } + continue; + } + + // texcoord + if (token[0] == 'v' && token[1] == 't' && IS_SPACE((token[2]))) { + token += 3; + real_t x, y, z; // y and z are optional. default = 0.0 + parseReal3(&x, &y, &z, &token); + if (callback.texcoord_cb) { + callback.texcoord_cb(user_data, x, y, z); + } + continue; + } + + // face + if (token[0] == 'f' && IS_SPACE((token[1]))) { + token += 2; + token += strspn(token, " \t"); + + indices.clear(); + while (!IS_NEW_LINE(token[0])) { + vertex_index vi = parseRawTriple(&token); + + index_t idx; + idx.vertex_index = vi.v_idx; + idx.normal_index = vi.vn_idx; + idx.texcoord_index = vi.vt_idx; + + indices.push_back(idx); + size_t n = strspn(token, " \t\r"); + token += n; + } + + if (callback.index_cb && indices.size() > 0) { + callback.index_cb(user_data, &indices.at(0), + static_cast(indices.size())); + } + + continue; + } + + // use mtl + if ((0 == strncmp(token, "usemtl", 6)) && IS_SPACE((token[6]))) { + token += 7; + std::stringstream ss; + ss << token; + std::string namebuf = ss.str(); + + int newMaterialId = -1; + if (material_map.find(namebuf) != material_map.end()) { + newMaterialId = material_map[namebuf]; + } else { + // { error!! material not found } + } + + if (newMaterialId != material_id) { + material_id = newMaterialId; + } + + if (callback.usemtl_cb) { + callback.usemtl_cb(user_data, namebuf.c_str(), material_id); + } + + continue; + } + + // load mtl + if ((0 == strncmp(token, "mtllib", 6)) && IS_SPACE((token[6]))) { + if (readMatFn) { + token += 7; + + std::vector filenames; + SplitString(std::string(token), ' ', filenames); + + if (filenames.empty()) { + if (err) { + (*err) += + "WARN: Looks like empty filename for mtllib. Use default " + "material. \n"; + } + } else { + bool found = false; + for (size_t s = 0; s < filenames.size(); s++) { + std::string err_mtl; + bool ok = (*readMatFn)(filenames[s].c_str(), &materials, + &material_map, &err_mtl); + if (err && (!err_mtl.empty())) { + (*err) += err_mtl; // This should be warn message. + } + + if (ok) { + found = true; + break; + } + } + + if (!found) { + if (err) { + (*err) += + "WARN: Failed to load material file(s). Use default " + "material.\n"; + } + } else { + if (callback.mtllib_cb) { + callback.mtllib_cb(user_data, &materials.at(0), + static_cast(materials.size())); + } + } + } + } + + continue; + } + + // group name + if (token[0] == 'g' && IS_SPACE((token[1]))) { + names.clear(); + + while (!IS_NEW_LINE(token[0])) { + std::string str = parseString(&token); + names.push_back(str); + token += strspn(token, " \t\r"); // skip tag + } + + assert(names.size() > 0); + + // names[0] must be 'g', so skip the 0th element. + if (names.size() > 1) { + name = names[1]; + } else { + name.clear(); + } + + if (callback.group_cb) { + if (names.size() > 1) { + // create const char* array. + names_out.resize(names.size() - 1); + for (size_t j = 0; j < names_out.size(); j++) { + names_out[j] = names[j + 1].c_str(); + } + callback.group_cb(user_data, &names_out.at(0), + static_cast(names_out.size())); + + } else { + callback.group_cb(user_data, NULL, 0); + } + } + + continue; + } + + // object name + if (token[0] == 'o' && IS_SPACE((token[1]))) { + // @todo { multiple object name? } + token += 2; + + std::stringstream ss; + ss << token; + std::string object_name = ss.str(); + + if (callback.object_cb) { + callback.object_cb(user_data, object_name.c_str()); + } + + continue; + } + +#if 0 // @todo + if (token[0] == 't' && IS_SPACE(token[1])) { + tag_t tag; + + token += 2; + std::stringstream ss; + ss << token; + tag.name = ss.str(); + + token += tag.name.size() + 1; + + tag_sizes ts = parseTagTriple(&token); + + tag.intValues.resize(static_cast(ts.num_ints)); + + for (size_t i = 0; i < static_cast(ts.num_ints); ++i) { + tag.intValues[i] = atoi(token); + token += strcspn(token, "/ \t\r") + 1; + } + + tag.floatValues.resize(static_cast(ts.num_reals)); + for (size_t i = 0; i < static_cast(ts.num_reals); ++i) { + tag.floatValues[i] = parseReal(&token); + token += strcspn(token, "/ \t\r") + 1; + } + + tag.stringValues.resize(static_cast(ts.num_strings)); + for (size_t i = 0; i < static_cast(ts.num_strings); ++i) { + std::stringstream ss; + ss << token; + tag.stringValues[i] = ss.str(); + token += tag.stringValues[i].size() + 1; + } + + tags.push_back(tag); + } +#endif + + // Ignore unknown command. + } + + if (err) { + (*err) += errss.str(); + } + + return true; +} +} // namespace tinyobj + +#endif diff --git a/src/utilities.h b/src/utilities.h index abb4f27..58ec9d8 100644 --- a/src/utilities.h +++ b/src/utilities.h @@ -13,6 +13,15 @@ #define TWO_PI 6.2831853071795864769252867665590057683943f #define SQRT_OF_ONE_THIRD 0.5773502691896257645091487805019574556476f #define EPSILON 0.00001f +#define ZeroEpsilon 0.0001f +#define InvPi 0.31830988618379067154f; +#define Inv2Pi 0.15915494309189533577f; +#define Inv4Pi 0.07957747154594766788f; +#define PiOver2 1.57079632679489661923f; +#define PiOver4 0.78539816339744830961f; +#define Sqrt2 1.41421356237309504880f; +#define OneMinusEpsilon 0.99999994f; + namespace utilityCore { extern float clamp(float f, float min, float max); diff --git a/stream_compaction/CMakeLists.txt b/stream_compaction/CMakeLists.txt index ac358c9..bcc484e 100644 --- a/stream_compaction/CMakeLists.txt +++ b/stream_compaction/CMakeLists.txt @@ -1,4 +1,16 @@ set(SOURCE_FILES + "common.h" + "common.cu" + "cpu.h" + "cpu.cu" + "naive.h" + "naive.cu" + "efficient.h" + "efficient.cu" + "thrust.h" + "thrust.cu" + "radix.h" + "radix.cu" ) cuda_add_library(stream_compaction diff --git a/stream_compaction/common.cu b/stream_compaction/common.cu new file mode 100644 index 0000000..b909064 --- /dev/null +++ b/stream_compaction/common.cu @@ -0,0 +1,48 @@ +#include "common.h" + + +void checkCUDAErrorFn(const char *msg, const char *file, int line) { + cudaError_t err = cudaGetLastError(); + if (cudaSuccess == err) { + return; + } + + fprintf(stderr, "CUDA error"); + if (file) { + fprintf(stderr, " (%s:%d)", file, line); + } + fprintf(stderr, ": %s: %s\n", msg, cudaGetErrorString(err)); + exit(EXIT_FAILURE); +} + + +namespace StreamCompaction { + namespace Common { + + /** + * Maps an array to an array of 0s and 1s for stream compaction. Elements + * which map to 0 will be removed, and elements which map to 1 will be kept. + */ + __global__ void kernMapToBoolean(int n, int *bools, const int *idata) { + int i = (threadIdx.x + (blockIdx.x * blockDim.x)); + if (i >= n) { return; } + + bools[i] = (int)(idata[i] != 0); + } + + /** + * Performs scatter on an array. That is, for each element in idata, + * if bools[idx] == 1, it copies idata[idx] to odata[indices[idx]]. + */ + __global__ void kernScatter(int n, int *odata, + const int *idata, const int *bools, const int *indices) { + int i = (threadIdx.x + (blockIdx.x * blockDim.x)); + if (i >= n) { return; } + + if (bools[i] == 1) { + odata[indices[i]] = idata[i]; + } + } + + } +} diff --git a/stream_compaction/common.h b/stream_compaction/common.h new file mode 100644 index 0000000..8b7dfbb --- /dev/null +++ b/stream_compaction/common.h @@ -0,0 +1,213 @@ +#pragma once + +#include +#include + +#include +#include +#include +#include +#include +#include + +#define FILENAME (strrchr(__FILE__, '/') ? strrchr(__FILE__, '/') + 1 : __FILE__) +#define checkCUDAError(msg) checkCUDAErrorFn(msg, FILENAME, __LINE__) + +#define NUM_BANKS 32 +#define LOG_NUM_BANKS 5 +#define CONFLICT_FREE_OFFSET(n) \ + ((n) >> NUM_BANKS + (n) >> (2 * LOG_NUM_BANKS)) + +//All of the following calculations were done for my GTX 1050 and +//they will probably not work with most other non 6.1CC GPU's + +// "x" in my master equation +#define MEM_PER_THREAD 4 + +// "z" in my master equation +#define OPTIMAL_BLOCKS_PER_SM 32 + +// Inequality/RHS in my master equation +#define SHARED_MEM_MAX 96000 + +// A muting term that tones done the optimality +#define TONE_DOWN 0.74f + +//The number of warps on this GPU +#define WARP_SIZE 32 + + + +/** + * Check for CUDA errors; print and exit if there was a problem. + */ +void checkCUDAErrorFn(const char *msg, const char *file = NULL, int line = -1); + +inline int binary(int num) +{ + if (num == 0) + { + return 0; + } + else + { + return (num % 2) + 10 * binary(num / 2); + } +} + +inline void printCPUArrayb(int n, int* arr) { + printf("\n( "); + for (int i = 0; i < n - 1; i++) { + printf("%d, ", binary(arr[i])); + } + printf("%d)", binary(arr[n - 1])); +} + +inline void printGPUArrayb(int n, int* dev_arr) { + int* cpu_arr = (int*)malloc(sizeof(int) * n); + cudaMemcpy(cpu_arr, dev_arr, sizeof(int)* n, cudaMemcpyDeviceToHost); + printCPUArrayb(n, cpu_arr); + free(cpu_arr); +} + + +inline void printCPUArray(int n, int* arr) { + printf("\n( "); + for (int i = 0; i < n-1; i++) { + printf("%d, ", arr[i]); + } + printf("%d)", arr[n-1]); +} + +inline void printGPUArray(int n, int* dev_arr) { + int* cpu_arr = (int*)malloc(sizeof(int) * n); + cudaMemcpy(cpu_arr, dev_arr, sizeof(int)* n, cudaMemcpyDeviceToHost); + printCPUArray(n, cpu_arr); + free(cpu_arr); +} + +inline int ilog2(int x) { + int lg = 0; + while (x >>= 1) { + ++lg; + } + return lg; +} + +inline int getThreadsPerBlock() { + //Get theoretical best y value you can based on GPU specs + //TODO: Find a better way to get the specs + int theoretical_y = SHARED_MEM_MAX * TONE_DOWN / (MEM_PER_THREAD * OPTIMAL_BLOCKS_PER_SM); + + //Find closest multiple of 32 to theoretical_y + for (int i = 0; i < SHARED_MEM_MAX / 32; i++) { + if (i * 32 >= theoretical_y) { + return (i - 1) * 32; + } + } +} + +inline int ilog2ceil(int x) { + return ilog2(x - 1) + 1; +} + +namespace StreamCompaction { + namespace Common { + __global__ void kernMapToBoolean(int n, int *bools, const int *idata); + + __global__ void kernScatter(int n, int *odata, + const int *idata, const int *bools, const int *indices); + + /** + * This class is used for timing the performance + * Uncopyable and unmovable + * + * Adapted from WindyDarian(https://github.com/WindyDarian) + */ + 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; + }; + } +} diff --git a/stream_compaction/cpu.cu b/stream_compaction/cpu.cu new file mode 100644 index 0000000..b49d988 --- /dev/null +++ b/stream_compaction/cpu.cu @@ -0,0 +1,130 @@ +#include +#include "cpu.h" +#include + +#include "common.h" + +namespace StreamCompaction { + namespace CPU { + using StreamCompaction::Common::PerformanceTimer; + PerformanceTimer& timer() + { + static PerformanceTimer timer; + return timer; + } + + /** + * CPU scan (prefix sum). + * For performance analysis, this is supposed to be a simple for loop. + * (Optional) For better understanding before starting moving to GPU, you can simulate your GPU scan in this function first. + */ + void scan(int n, int *odata, const int *idata) { + try { + timer().startCpuTimer(); + } + catch (...){}; + /* FAILED ATTEMPT AT GPU WAY + //Copy idata into odata + for (int i = 0; i < n; i++) { + odata[i] = idata[i]; + } + + // Create Constants + const int logN = ilog2ceil(n); + + printf("LOG N IS : %d \n", logN); + + //Up-Sweep + for (int d = 0; d < logN; d++) { + for (int k = 0; k < n; k += (int)pow(2, d + 1)) { + odata[k + (int)pow(2, d + 1) - 1] += odata[(k + (int)pow(2, d) - 1)]; + } + } + + printf("UPSWEPT: \n ("); + for (int i = 0; i < 15; i++) { + printf("%d, ", odata[i]); + } + printf(".... )\n"); + + //Down-Sweep + odata[n - 1] = 0; + for (int d = logN-1; d >= 0; d--) { + for (int k = 0; k < n; k += (int)pow(2, d + 1)) { + if ((k + (int)pow(2, d+1) - 1) < n) { + int t = odata[(k + (int)pow(2, d) - 1)]; + odata[k + (int)pow(2, d) - 1] = odata[k + (int)pow(2, d + 1) - 1]; + odata[k + (int)pow(2, d + 1) - 1] += t; + } + } + } + + printf("SUMMED: \n ("); + for (int i = 0; i < 15; i++) { + printf("%d, ", odata[i]); + } + printf(".... )\n"); + + */ + odata[0] = 0; + for (int k = 1; k < n; k++) { + odata[k] = idata[k-1] + odata[k-1]; + } + + try { + timer().endCpuTimer(); + } + catch (...) {}; + } + + /** + * CPU stream compaction without using the scan function. + * + * @returns the number of elements remaining after compaction. + */ + int compactWithoutScan(int n, int *odata, const int *idata) { + timer().startCpuTimer(); + int count = 0; + for (int i = 0; i < n; i++) { + if (idata[i] != 0) { + odata[count++] = idata[i]; + } + } + timer().endCpuTimer(); + + return count; + } + + /** + * CPU stream compaction using scan and scatter, like the parallel version. + * + * @returns the number of elements remaining after compaction. + */ + int compactWithScan(int n, int *odata, const int *idata) { + timer().startCpuTimer(); + int count = 0; + int* temp = (int*) malloc(sizeof(int) * n); + int* scanned = (int*) malloc(sizeof(int) * n); + + for (int i = 0; i < n; i++) { + temp[i] = (int) (idata[i] != 0); + count = idata[i] != 0 ? count + 1 : count; + } + + scan(n, scanned, temp); + + for (int i = 0; i < n; i++) { + if (temp[i] == 1) { + odata[scanned[i]] = idata[i]; + } + } + + try { + timer().endCpuTimer(); + } + catch (...) {}; + + return count; + } + } +} diff --git a/stream_compaction/cpu.h b/stream_compaction/cpu.h new file mode 100644 index 0000000..236ce11 --- /dev/null +++ b/stream_compaction/cpu.h @@ -0,0 +1,15 @@ +#pragma once + +#include "common.h" + +namespace StreamCompaction { + namespace CPU { + StreamCompaction::Common::PerformanceTimer& timer(); + + void scan(int n, int *odata, const int *idata); + + int compactWithoutScan(int n, int *odata, const int *idata); + + int compactWithScan(int n, int *odata, const int *idata); + } +} diff --git a/stream_compaction/efficient.cu b/stream_compaction/efficient.cu new file mode 100644 index 0000000..334dc8e --- /dev/null +++ b/stream_compaction/efficient.cu @@ -0,0 +1,684 @@ +#include +#include +#include "common.h" +#include "efficient.h" + +#define BLOCK_SIZE 1024 + +namespace StreamCompaction { + namespace Efficient { + using StreamCompaction::Common::PerformanceTimer; + PerformanceTimer& timer() + { + static PerformanceTimer timer; + return timer; + } + + __global__ void kern_upSweep(int n, int d, int* idata) { + int index = (threadIdx.x + (blockIdx.x * blockDim.x)); + int k = index * (1 << d + 1); + if (index >= n || k >= n) { return; } + + idata[k + (1 << d+1) - 1] += idata[k + (1 << d) - 1]; + } + + __global__ void kern_downSweep(int n, int d, int* idata) { + int k = (threadIdx.x + (blockIdx.x * blockDim.x)) * (1 << d + 1); + if (k >= n) { return; } + + int t = idata[k + (1 << d) - 1]; + idata[k + (1 << d) - 1] = idata[k + (1 << d+1) - 1]; + idata[k + (1 << d+1) - 1] += t; + } + + __global__ void roundN(int n, int nRounded, int* idataRounded, const int* idata) { + int i = (threadIdx.x + (blockIdx.x * blockDim.x)); + if (i >= nRounded) { return; } + + idataRounded[i] = i >= n ? 0 : idata[i]; + } + + __global__ void kern_add_sums(int n, int* sums, int* data) { + // The thread's ID within the entire grid + int threadId = threadIdx.x + (blockIdx.x * blockDim.x); + + if (threadId >= n) { return; } + + int sumIdx = threadId / blockDim.x; + + data[threadId] += sums[sumIdx]; + } + + /** + + __global__ void kern_scan_shared_fixed(int* sums, int *odata, const int *idata) { + //Copy the arary into shared memory + extern __shared__ int temp[]; + + //Halved because we're doing double the loading + const int blockSizeHalf = BLOCK_SIZE * 0.5f; + + //The thread's ID within a block + int threadId = threadIdx.x; + int threadId2 = threadIdx.x + blockDim.x; + + // The thread's ID within the entire grid + int threadId_global = threadIdx.x + (blockIdx.x * blockDim.x); + int threadId2_global = threadId_global + blockDim.x; + + // Offsets for avoiding bank conflicts + int offsetA = CONFLICT_FREE_OFFSET(threadId); + int offsetB = CONFLICT_FREE_OFFSET(threadId2); + + // Load in elts + temp[threadId + offsetA] = idata[threadId_global]; + temp[threadId2 + offsetB] = idata[threadId2_global]; + + // Offset + int stride = 1; + + //Do the UpSweep + //From d = n to 0, log_n times + //#pragma unroll BLOCK_SIZE + for (int d = blockSizeHalf >> 1; d > 0; d >>= 1) { + __syncthreads(); + + if (threadId < d) + { + int ai = stride * (2*threadId + 1) - 1; + int bi = stride * (2*threadId + 2) - 1; + + ai += CONFLICT_FREE_OFFSET(ai); + bi += CONFLICT_FREE_OFFSET(bi); + + temp[bi] += temp[ai]; + } + stride <<= 1; + } + + int idxLast = blockSizeHalf - 1 + CONFLICT_FREE_OFFSET(blockSizeHalf - 1); + + //Set the Sum to the + if (sums != nullptr && threadId == 0) { + sums[blockIdx.x] = temp[idxLast]; + } + + //Temp of n-1 = zero + if (threadId == 0) { temp[idxLast] = 0; } + + //#pragma unroll BLOCK_SIZE + for (int d = 1; d < blockSizeHalf; d <<= 1) { + stride >>= 1; + __syncthreads(); + + if (threadId < d) { + int ai = stride * (2* threadId + 1) - 1; + int bi = stride * (2* threadId + 2) - 1; + + ai += CONFLICT_FREE_OFFSET(ai); + bi += CONFLICT_FREE_OFFSET(bi); + + int t = temp[ai]; //Save Left Child + temp[ai] = temp[bi]; //Store Right Child in left Child's place + temp[bi] += t; //Right child += left child + } + } + + __syncthreads(); + + + //WRITE OUT DA VALUES + odata[threadId_global] = temp[threadId + offsetA]; + odata[threadId_global] = temp[threadId + offsetB]; + } + + + // __global__ void kern_scan_shared(int n, int* sums, int *odata, const int *idata) { + // //Copy the arary into shared memory + // extern __shared__ int temp[]; + // + // //The thread's ID within a block + // int threadId = threadIdx.x; + // + // // The thread's ID within the entire grid + // int threadId_global = threadIdx.x + (blockIdx.x * blockDim.x); + // + // // Load in elts + // temp[threadId] = idata[threadId_global]; + // + // // Offset + // int stride = 1; + // + // //Do the UpSweep + // //From d = n to 0, log_n times + // for (int d = n >> 1; d > 0; d >>= 1) { + // __syncthreads(); + // + // if (threadId < d) + // { + // int ai = stride * (2 * threadId + 1) - 1; + // int bi = stride * (2 * threadId + 2) - 1; + // + // temp[bi] += temp[ai]; + // } + // stride <<= 1; + // } + // + // //Set the Sum to the + // if (sums != nullptr && threadId == 0) { + // sums[blockIdx.x] = temp[n - 1]; + // } + // + // //Temp of n-1 = zero + // if (threadId == 0) { temp[n - 1] = 0; } + // + // //Down Sweep + // for (int d = 1; d < n; d <<= 1) { + // stride >>= 1; + // __syncthreads(); + // + // if (threadId < d) { + // int ai = stride * (2 * threadId + 1) - 1; + // int bi = stride * (2 * threadId + 2) - 1; + // + // int t = temp[ai]; //Save Left Child + // temp[ai] = temp[bi]; //Store Right Child in left Child's place + // temp[bi] += t; //Right child += left child + // } + // } + // + // __syncthreads(); + // + // + // //WRITE OUT DA VALUES + // odata[threadId_global] = temp[threadId]; + // } + + __global__ void kern_scan_shared(int n, int* sums, int *odata, const int *idata) { + //Copy the arary into shared memory + extern __shared__ int temp[]; + + //The thread's ID within a block + int threadId = threadIdx.x; + int threadId2 = threadIdx.x + (n/2); + + // The thread's ID within the entire grid + int threadId_global = threadIdx.x + (blockIdx.x * blockDim.x); + int threadId2_global = threadId_global + (n/2); + + // Offsets for avoiding bank conflicts + int offsetA = CONFLICT_FREE_OFFSET(threadId); + int offsetB = CONFLICT_FREE_OFFSET(threadId2); + + // Load in elts + temp[threadId + offsetA] = idata[threadId_global]; + temp[threadId2 + offsetB] = idata[threadId2_global]; + + // Offset + int stride = 1; + + //Do the UpSweep + //From d = n to 0, log_n times + //#pragma unroll BLOCK_SIZE + for (int d = n >> 1; d > 0; d >>= 1) { + __syncthreads(); + + if (threadId < d) + { + int ai = stride * (2 * threadId + 1) - 1; + int bi = stride * (2 * threadId + 2) - 1; + + ai += CONFLICT_FREE_OFFSET(ai); + bi += CONFLICT_FREE_OFFSET(bi); + + temp[bi] += temp[ai]; + } + stride <<= 1; + } + + int idxLast = n - 1 + CONFLICT_FREE_OFFSET(n - 1); + + //Set the Sum to the + if (sums != nullptr && threadId == 0) { + sums[blockIdx.x] = temp[idxLast]; + } + + //Temp of n-1 = zero + if (threadId == 0) { temp[idxLast] = 0; } + + //#pragma unroll BLOCK_SIZE + for (int d = 1; d < n; d <<= 1) { + stride >>= 1; + __syncthreads(); + + if (threadId < d) { + int ai = stride * (2 * threadId + 1) - 1; + int bi = stride * (2 * threadId + 2) - 1; + + ai += CONFLICT_FREE_OFFSET(ai); + bi += CONFLICT_FREE_OFFSET(bi); + + int t = temp[ai]; //Save Left Child + temp[ai] = temp[bi]; //Store Right Child in left Child's place + temp[bi] += t; //Right child += left child + } + } + + __syncthreads(); + + + //WRITE OUT DA VALUES + odata[threadId_global] = temp[threadId + offsetA]; + odata[threadId_global] = temp[threadId + offsetB]; + } + + + //* + //Performs prefix-sum (aka scan) on idata, storing the result into odata. + // + // *** THIS IS AN EFFICIENT VERSION USING SHARED MEMORY + + void scan_shared(int n, int *odata, const int *idata) { + //First things first, round n to the nearest power of two + int loops = ilog2ceil(n); + int nRounded = 1 << loops; + + // Also round getThreadsPerBlock to nearest power of two to avoid AIOOB errors + // An "if thread >= n" check would also work but that leads to more divergence. + int threads = ilog2ceil(getThreadsPerBlock()); + int threadCount = 1 << threads; + + // Super Hyperthreaded Information Transloading calculation for threads per block + dim3 threadsPerBlock(std::min(BLOCK_SIZE,nRounded)); + dim3 numBlocks(std::ceilf(((float) nRounded / threadsPerBlock.x))); + + int blockCount = numBlocks.x; + + printf("Rounding from: %d to %d because we have %d blocks that are %d ints big\n", n, nRounded, blockCount, threadsPerBlock.x); + + // A copy of idata on the GPU + int *idata_GPU, *idataRounded_GPU, *odata_GPU, *isums_GPU, *osums_GPU; + cudaMalloc((void**)&idata_GPU , sizeof(int) * n); + cudaMalloc((void**)&idataRounded_GPU, sizeof(int) * nRounded); + cudaMalloc((void**)&odata_GPU , sizeof(int) * nRounded); + cudaMemcpy(idata_GPU, idata, sizeof(int) * n, cudaMemcpyHostToDevice); + cudaDeviceSynchronize(); + if (nRounded != n) { + roundN << > > (n, nRounded, idataRounded_GPU, idata_GPU); + } + else { + idataRounded_GPU = idata_GPU; + } + + //printGPUArray(nRounded, idataRounded_GPU); + + try { timer().startGpuTimer(); } + catch (...) {}; + + //If all the data fits within one block + if (blockCount == 1) { + cudaMalloc((void**)&isums_GPU, sizeof(int) * blockCount); + cudaMalloc((void**)&osums_GPU, sizeof(int) * blockCount); + kern_scan_shared<<> >(threadsPerBlock.x, nullptr, odata_GPU, idataRounded_GPU); + } + else { + //If we need to use the sums + //Scan all the blocks separately + cudaMalloc((void**)&isums_GPU, sizeof(int) * blockCount * 2.f); + cudaMalloc((void**)&osums_GPU, sizeof(int) * blockCount * 2.f); + blockCount *= 2; + dim3 numBlocksDouble = dim3(blockCount); + dim3 threadsPerBlockHalf = dim3(threadsPerBlock.x * 0.5f); + kern_scan_shared_fixed << < numBlocksDouble, threadsPerBlockHalf , sizeof(int) * threadsPerBlock.x>> > (isums_GPU, odata_GPU, idataRounded_GPU); + + //printGPUArray(n, odata_GPU); + //Run scan on the sums array + kern_scan_shared <<>> (blockCount, nullptr, osums_GPU, isums_GPU); + //printGPUArray(blockCount, osums_GPU); + + //Add result of scan back to the per-block scans + kern_add_sums << > > (n, osums_GPU, odata_GPU); + } + try { timer().endGpuTimer(); } + catch (...) {}; + + //printGPUArray(n, odata_GPU); + + //printf("\n PRINT THIS BITCH \n\n"); + //printGPUArray(blockCount, osums_GPU); + + cudaMemcpy(odata, odata_GPU, sizeof(int) * n, cudaMemcpyDeviceToHost); + cudaFree(idata_GPU); + cudaFree(idataRounded_GPU); + cudaFree(odata_GPU); + cudaFree(isums_GPU); + cudaFree(osums_GPU); + } + + */ + + __global__ void kern_scan_shared_fixed(int* sums, int *odata, const int *idata) { + //Copy the arary into shared memory + extern __shared__ int temp[]; + + //The thread's ID within a block + int threadId = threadIdx.x; + + // The thread's ID within the entire grid + int threadId_global = threadIdx.x + (blockIdx.x * blockDim.x); + + // Offsets for avoiding bank conflicts + //int offsetA = ; + //int offsetB = ; + + // Load in elts + temp[threadId] = idata[threadId_global]; + + // Offset + int stride = 1; + + //Do the UpSweep + //From d = n to 0, log_n times + #pragma unroll BLOCK_SIZE + for (int d = BLOCK_SIZE >> 1; d > 0; d >>= 1) { + __syncthreads(); + + if (threadId < d) + { + int ai = stride * (2 * threadId + 1) - 1; + int bi = stride * (2 * threadId + 2) - 1; + + temp[bi] += temp[ai]; + } + stride <<= 1; + } + + //Set the Sum to the + if (sums != nullptr && threadId == 0) { + sums[blockIdx.x] = temp[BLOCK_SIZE - 1]; + } + + //Temp of n-1 = zero + if (threadId == 0) { temp[BLOCK_SIZE - 1] = 0; } + + int offset = 0; + + #pragma unroll BLOCK_SIZE + for (int d = 1; d < BLOCK_SIZE; d <<= 1) { + stride >>= 1; + __syncthreads(); + + if (threadId < d) { + int ai = stride * (2 * threadId + 1) - 1; + int bi = stride * (2 * threadId + 2) - 1; + + int t = temp[ai]; //Save Left Child + temp[ai] = temp[bi]; //Store Right Child in left Child's place + temp[bi] += t; //Right child += left child + } + } + + __syncthreads(); + + //WRITE OUT DA VALUES + odata[threadId_global] = temp[threadId]; + } + + __global__ void kern_scan_shared(int n, int* sums, int *odata, const int *idata) { + //Copy the arary into shared memory + extern __shared__ int temp[]; + + //The thread's ID within a block + int threadId = threadIdx.x; + + // The thread's ID within the entire grid + int threadId_global = threadIdx.x + (blockIdx.x * blockDim.x); + + // Offsets for avoiding bank conflicts + //int offsetA = ; + //int offsetB = ; + + // Load in elts + temp[threadId] = idata[threadId_global]; + + // Offset + int stride = 1; + + //Do the UpSweep + //From d = n to 0, log_n times + for (int d = n >> 1; d > 0; d >>= 1) { + __syncthreads(); + + if (threadId < d) + { + int ai = stride * (2 * threadId + 1) - 1; + int bi = stride * (2 * threadId + 2) - 1; + + temp[bi] += temp[ai]; + } + stride <<= 1; + } + + //Set the Sum to the + if (sums != nullptr && threadId == 0) { + sums[blockIdx.x] = temp[n - 1]; + } + + //Temp of n-1 = zero + if (threadId == 0) { temp[n - 1] = 0; } + + int offset = 0; + + //Down Sweep + for (int d = 1; d < n; d <<= 1) { + stride >>= 1; + __syncthreads(); + + if (threadId < d) { + int ai = stride * (2 * threadId + 1) - 1; + int bi = stride * (2 * threadId + 2) - 1; + + int t = temp[ai]; //Save Left Child + temp[ai] = temp[bi]; //Store Right Child in left Child's place + temp[bi] += t; //Right child += left child + } + } + + __syncthreads(); + + + //WRITE OUT DA VALUES + odata[threadId_global] = temp[threadId]; + } + + //* + // Performs prefix-sum (aka scan) on idata, storing the result into odata. + // + //*** THIS IS AN EFFICIENT VERSION USING SHARED MEMORY + // + void scan_shared(int n, int *odata, const int *idata) { + //First things first, round n to the nearest power of two + int loops = ilog2ceil(n); + int nRounded = 1 << loops; + + // Also round getThreadsPerBlock to nearest power of two to avoid AIOOB errors + // An "if thread >= n" check would also work but that leads to more divergence. + int threads = ilog2ceil(getThreadsPerBlock()); + int threadCount = 1 << threads; + + // Super Hyperthreaded Information Transloading calculation for threads per block + dim3 threadsPerBlock(std::min(BLOCK_SIZE, nRounded)); + dim3 numBlocks(std::ceilf(((float)nRounded / threadsPerBlock.x))); + + int blockCount = numBlocks.x; + + printf("Rounding from: %d to %d because we have %d blocks that are %d ints big\n", n, nRounded, blockCount, threadsPerBlock.x); + + // A copy of idata on the GPU + int *idata_GPU, *idataRounded_GPU, *odata_GPU, *isums_GPU, *osums_GPU; + cudaMalloc((void**)&idata_GPU, sizeof(int) * n); + cudaMalloc((void**)&idataRounded_GPU, sizeof(int) * nRounded); + cudaMalloc((void**)&odata_GPU, sizeof(int) * nRounded); + cudaMalloc((void**)&isums_GPU, sizeof(int) * blockCount); + cudaMalloc((void**)&osums_GPU, sizeof(int) * blockCount); + cudaMemcpy(idata_GPU, idata, sizeof(int) * n, cudaMemcpyHostToDevice); + cudaDeviceSynchronize(); + if (nRounded != n) { + roundN << > > (n, nRounded, idataRounded_GPU, idata_GPU); + } + else { + idataRounded_GPU = idata_GPU; + } + + //printGPUArray(nRounded, idataRounded_GPU); + + try { timer().startGpuTimer(); } + catch (...) {}; + + //If all the data fits within one block + if (blockCount == 1) { + kern_scan_shared << > >(threadsPerBlock.x, nullptr, odata_GPU, idataRounded_GPU); + } + else { //If we need to use the sums + //Scan all the blocks separately + kern_scan_shared_fixed << > >(isums_GPU, odata_GPU, idataRounded_GPU); + + //Run scan on the sums array + kern_scan_shared << > > (blockCount, nullptr, osums_GPU, isums_GPU); + + //Add result of scan back to the per-block scans + kern_add_sums << > > (n, osums_GPU, odata_GPU); + } + try { timer().endGpuTimer(); } + catch (...) {}; + + //printGPUArray(n, odata_GPU); + + //printf("\n PRINT THIS BITCH \n\n"); + //printGPUArray(blockCount, osums_GPU); + + cudaMemcpy(odata, odata_GPU, sizeof(int) * n, cudaMemcpyDeviceToHost); + cudaFree(idata_GPU); + cudaFree(idataRounded_GPU); + cudaFree(odata_GPU); + cudaFree(isums_GPU); + cudaFree(osums_GPU); + } + + + /** + * Performs prefix-sum (aka scan) on idata, storing the result into odata. + */ + void scan(int n, int *odata, const int *idata) { + // Super Hyperthreaded Information Transloading calculation for threads per block + dim3 threadsPerBlock(std::min(getThreadsPerBlock(), n)); + dim3 numBlocks(std::ceilf(((float)n / threadsPerBlock.x))); + + //Round Up + int loops = ilog2ceil(n); + int nRounded = 1 << loops; + + // A copy of idata on the GPU + int* idata_GPU; + cudaMalloc((void**)&idata_GPU, sizeof(int) * n); + cudaMemcpy(idata_GPU, idata, sizeof(int) * n, cudaMemcpyHostToDevice); + + try { timer().startGpuTimer(); } + catch (...) {}; + + //Rounded Version of GPU Copy + int* idataRounded_GPU; + cudaMalloc((void**)&idataRounded_GPU, sizeof(int) * nRounded); + //Round the GPU Array + roundN << > > (n, nRounded, idataRounded_GPU, idata_GPU); + + //Up-Sweep: + for (int d = 0; d < loops; d++) { + kern_upSweep<< > >(n, d, idataRounded_GPU); + checkCUDAErrorFn("upSweep failed with error"); + } + + //Set Zero + int zero = 0; + cudaMemcpy(&idataRounded_GPU[nRounded - 1], &zero, sizeof(int), cudaMemcpyHostToDevice); + checkCUDAErrorFn("Zero Copy failed with error"); + + //Down-Sweep: + for (int d = loops - 1; d >= 0; d--) { + kern_downSweep <<> >(nRounded, d, idataRounded_GPU); + checkCUDAErrorFn("downSweep failed with error"); + } + + cudaMemcpy(odata, idataRounded_GPU, sizeof(int) * n, cudaMemcpyDeviceToHost); + + //Free Malloc'd + cudaFree(idataRounded_GPU); + cudaFree(idata_GPU); + /**** PRINTER ****** + printf("After DownSweep: \n ("); + for (int i = nRounded-10; i < nRounded -1; i++) { + printf("%d = %d, ", i, odata[i]); + } + printf("%d = %d) \n\n", nRounded-1, odata[nRounded-1]); + **/ + try { timer().endGpuTimer(); } + catch (...) {}; + } + + /** + * Performs stream compaction on idata, storing the result into odata. + * All zeroes are discarded. + * + * @param n The number of elements in idata. + * @param odata The array into which to store elements. + * @param idata The array of elements to compact. + * @returns The number of elements remaining after compaction. + */ + int compact(int n, int *odata, const int *idata) { + try { timer().startGpuTimer(); } + catch (...) {}; + // Super Hyperthreaded Information Transloading calculation for threads per block + dim3 threadsPerBlock(std::min(getThreadsPerBlock(), n)); + dim3 numBlocks(std::ceilf(((float)n / threadsPerBlock.x))); + + // Create Buffers + int *temp, *scanned, *idata_GPU, *odata_GPU, *count_GPU; + cudaMalloc((void**)&temp , sizeof(int) * n); + cudaMalloc((void**)&scanned , sizeof(int) * n); + cudaMalloc((void**)&idata_GPU, sizeof(int) * n); + cudaMalloc((void**)&odata_GPU, sizeof(int) * n); + + cudaMemcpy(idata_GPU, idata, sizeof(int) * n, cudaMemcpyHostToDevice); + checkCUDAErrorFn("idata memcpy failed with error"); + + //Create Temp Array + Common::kernMapToBoolean << > > (n, temp, idata_GPU); + checkCUDAErrorFn("kern_boolify failed with error"); + + //Scan + scan(n, scanned, temp); + + //Temporarily store "scanned" into odata to get count + cudaMemcpy(odata, scanned, sizeof(int) * n, cudaMemcpyDeviceToHost); + int count = odata[n - 1] + (int)(idata[n - 1] != 0); + + //Compact + Common::kernScatter << > >(n, odata_GPU, idata_GPU, temp, scanned); + checkCUDAErrorFn("kern_compact failed with error"); + + //Bring Back to CPU + cudaMemcpy(odata, odata_GPU, sizeof(int) * count, cudaMemcpyDeviceToHost); + + //Free Up All Malloc'd + cudaFree(temp); + cudaFree(scanned); + cudaFree(idata_GPU); + cudaFree(odata_GPU); + + + try { timer().endGpuTimer(); } + catch (...) {}; + return count; + } + } +} diff --git a/stream_compaction/efficient.h b/stream_compaction/efficient.h new file mode 100644 index 0000000..7ff29eb --- /dev/null +++ b/stream_compaction/efficient.h @@ -0,0 +1,15 @@ +#pragma once + +#include "common.h" + +namespace StreamCompaction { + namespace Efficient { + StreamCompaction::Common::PerformanceTimer& timer(); + + void scan_shared(int n, int *odata, const int *idata); + + void scan(int n, int *odata, const int *idata); + + int compact(int n, int *odata, const int *idata); + } +} diff --git a/stream_compaction/naive.cu b/stream_compaction/naive.cu new file mode 100644 index 0000000..2297edb --- /dev/null +++ b/stream_compaction/naive.cu @@ -0,0 +1,72 @@ +#include +#include +#include "common.h" +#include "naive.h" + +namespace StreamCompaction { + namespace Naive { + using StreamCompaction::Common::PerformanceTimer; + PerformanceTimer& timer() + { + static PerformanceTimer timer; + return timer; + } + + + // TODO: __global__ + __global__ void kern_scan(int n, int d, int* odata, const int* idata) { + int k = threadIdx.x + (blockIdx.x * blockDim.x); + if (k >= n) { return; } + + int two_powd = 1 << (d - 1); + //Using ternary as recommended in lecture + odata[k] = (k >= two_powd) ? idata[k - two_powd] + idata[k] : idata[k]; + } + + __global__ void kern_shiftRight(int n, int* odata, const int* idata) { + int k = threadIdx.x + (blockIdx.x * blockDim.x); + if (k >= n) { return; } + + //Using ternary as recommended in lecture + odata[k] = (k == 0) ? 0 : idata[k - 1]; + } + + /** + * Performs prefix-sum (aka scan) on idata, storing the result into odata. + */ + void scan(int n, int *odata, const int *idata) { + timer().startGpuTimer(); + // Super Hyperthreaded Information Transloading calculation for threads per block + dim3 threadsPerBlock(std::min(getThreadsPerBlock(), n)); + dim3 numBlocks(std::ceilf(((float) n / threadsPerBlock.x) )); + + //New GPU Buffers + int* idata_GPU, *odata_GPU; + cudaMalloc((void**)&idata_GPU, sizeof(int) * n); + cudaMalloc((void**)&odata_GPU, sizeof(int) * n); + cudaMemcpy(idata_GPU, idata, sizeof(int) * n, cudaMemcpyHostToDevice); + + for (int d = 1; d <= ilog2ceil(n); d++) { + kern_scan <<>>(n, d, odata_GPU, idata_GPU); + checkCUDAErrorFn("kern_scan failed with error"); + std::swap(idata_GPU, odata_GPU); + } + kern_shiftRight<< > > (n, odata_GPU, idata_GPU); + + cudaMemcpy(odata, odata_GPU, sizeof(int) * n, cudaMemcpyDeviceToHost); + timer().endGpuTimer(); + + cudaFree(idata_GPU); + cudaFree(odata_GPU); + + //PRINTER + /* + printf("After: \n ( "); + for (int i = 0; i < 15; i++) { + printf("%d, ", odata[i]); + } + printf(" ... ) \n\n"); + */ + } + } +} diff --git a/stream_compaction/naive.h b/stream_compaction/naive.h new file mode 100644 index 0000000..37dcb06 --- /dev/null +++ b/stream_compaction/naive.h @@ -0,0 +1,11 @@ +#pragma once + +#include "common.h" + +namespace StreamCompaction { + namespace Naive { + StreamCompaction::Common::PerformanceTimer& timer(); + + void scan(int n, int *odata, const int *idata); + } +} diff --git a/stream_compaction/radix.cu b/stream_compaction/radix.cu new file mode 100644 index 0000000..d9954a3 --- /dev/null +++ b/stream_compaction/radix.cu @@ -0,0 +1,142 @@ +#include +#include +#include "common.h" +#include "radix.h" +#include "efficient.h" + + +namespace StreamCompaction { + namespace Radix { + using StreamCompaction::Common::PerformanceTimer; + PerformanceTimer& timer() + { + static PerformanceTimer timer; + return timer; + } + + __global__ void kern_initializeHistogram(int n, int* hists) { + int i = threadIdx.x + (blockIdx.x * blockDim.x); + if (i >= n) { return; } + + hists[i] = 0; + } + + + __global__ void kern_generateHistogram(int n, int bit, int desired_bit, const int* idata, int* hists) { + int i = threadIdx.x + (blockIdx.x * blockDim.x); + if (i >= n) { return; } + + int nth_bit = (idata[i] >> bit) & 1; + + //Determine if Zero + hists[i] = (int)(nth_bit == desired_bit); + } + + __global__ void kern_placeItems(int n, int* odata, const int* idata, const int* hists0, const int* hists1, + const int* offset0, const int* offset1) { + int i = threadIdx.x + (blockIdx.x * blockDim.x); + if (i >= n) { return; } + + int zero_count = offset0[n - 1] + hists0[n - 1]; + + //Zero Bit + if (hists0[i] == 1) { + odata[offset0[i]] = idata[i]; + } else { + odata[zero_count + offset1[i]] = idata[i]; + } + } + + /** + * Performs radix sort on a buffer, returns sorted array in odata + * + * @param n The number of elements in idata. + * @param odata The array into which to sort the elements. + * @param idata The array of elements to sort. + */ + void sort(int n, int *odata, const int *idata) { + try { timer().startGpuTimer(); } + catch (...) {}; + //Define Thread and Block Counts + dim3 threadsPerBlock(std::min(getThreadsPerBlock(), n)); + dim3 numBlocks(std::ceilf(((float)n / threadsPerBlock.x))); + + //Allocate GPU Buffers + int* idata_GPU, *odata_GPU; + cudaMalloc((void**)&idata_GPU, sizeof(int) * n); + cudaMalloc((void**)&odata_GPU, sizeof(int) * n); + cudaMemcpy(idata_GPU, idata, sizeof(int) * n, cudaMemcpyHostToDevice); + + int* emptyCounts = (int*) malloc(sizeof(int) * 2 * n); + //Get Maximum Number + int max = -INT_MAX; + for (int i = 0; i < n; i++) { + max = std::max(idata[i], max); + + //Zero out array of counts + *(emptyCounts + 0 * n + i) = 0; + *(emptyCounts + 1 * n + i) = 0; + } + int max_MSB = ilog2ceil(max); + + // ---- For each bit: + for (int bit = 0; bit < max_MSB; bit++) { + //Create a hist[n] Array, with each thread writing to its index + int* histograms_zero; // This is an Array of counts of zero + cudaMalloc((void**)&histograms_zero, sizeof(int) * n); + + //Create a hist[n] Array, with each thread writing to its index + int* histograms_one; // This is an Array of counts of zero + cudaMalloc((void**)&histograms_one, sizeof(int) * n); + + //Create a prefixSum array that lets you do cool stuff to it + int* offset_zero; + cudaMalloc((void**)&offset_zero, sizeof(int) * n); + + //Create a prefixSum array that lets you do cool stuff to it + int* offset_one; + cudaMalloc((void**)&offset_one, sizeof(int) * n); + + // kern: Generate Histogram with Number of zero's and one's for each bit + kern_generateHistogram << > > (n, bit, 0, idata_GPU, histograms_zero); + kern_generateHistogram << > > (n, bit, 1, idata_GPU, histograms_one); + + // kern: Calculate Offsets for each + StreamCompaction::Efficient::scan(n, offset_zero, histograms_zero); + StreamCompaction::Efficient::scan(n, offset_one, histograms_one); + + // kern: Ex Prefix Sum on Histogram + //REMEMBER THE MEMCPY OPTIMIZATION + /** + int last_zero_elt, last_offset_zero; + cudaMemcpy(&last_zero_elt, &histograms_zero[n - 1], sizeof(int), cudaMemcpyDeviceToHost); + cudaMemcpy(&last_offset_zero, &offset_zero[n - 1], sizeof(int), cudaMemcpyDeviceToHost); + int zero_count = last_offset_zero + last_zero_elt; + **/ + + // kern: Place each thing int its location + kern_placeItems << > > (n, odata_GPU, idata_GPU, histograms_zero, histograms_one, + offset_zero, offset_one); + + // Ping Pong Buffers + std::swap(idata_GPU, odata_GPU); + + //Free Old Stuff + cudaFree(histograms_zero); + cudaFree(histograms_one); + cudaFree(offset_zero); + cudaFree(offset_one); + } + + cudaMemcpy(odata, idata_GPU, sizeof(int) * n, cudaMemcpyDeviceToHost); + + cudaFree(idata_GPU); + cudaFree(odata_GPU); + free(emptyCounts); + + try { timer().endGpuTimer(); } + catch (...) {}; + + } + } +} diff --git a/stream_compaction/radix.h b/stream_compaction/radix.h new file mode 100644 index 0000000..aa489d9 --- /dev/null +++ b/stream_compaction/radix.h @@ -0,0 +1,11 @@ +#pragma once + +#include "common.h" + +namespace StreamCompaction { + namespace Radix { + StreamCompaction::Common::PerformanceTimer& timer(); + + void sort(int n, int *odata, const int *idata); + } +} diff --git a/stream_compaction/thrust.cu b/stream_compaction/thrust.cu new file mode 100644 index 0000000..e499bc7 --- /dev/null +++ b/stream_compaction/thrust.cu @@ -0,0 +1,41 @@ +#include +#include +#include +#include +#include +#include "common.h" +#include "thrust.h" + +namespace StreamCompaction { + namespace Thrust { + using StreamCompaction::Common::PerformanceTimer; + PerformanceTimer& timer() + { + static PerformanceTimer timer; + return timer; + } + /** + * Performs prefix-sum (aka scan) on idata, storing the result into odata. + */ + void scan(int n, int *odata, const int *idata) { + int* dv_in, *dv_out; + cudaMalloc((void**)&dv_in, sizeof(int) * n); + cudaMalloc((void**)&dv_out, sizeof(int) * n); + + cudaMemcpy(dv_in, idata, sizeof(int) * n, cudaMemcpyHostToDevice); + + thrust::device_ptr dv_in_thrust(dv_in); + thrust::device_ptr dv_out_thrust(dv_out); + thrust::exclusive_scan(dv_in_thrust, dv_in_thrust + n, dv_out_thrust); + + timer().startGpuTimer(); + thrust::exclusive_scan(dv_in_thrust, dv_in_thrust + n, dv_out_thrust); + timer().endGpuTimer(); + + cudaMemcpy(odata, dv_out, sizeof(int) * n, cudaMemcpyDeviceToHost); + + cudaFree(dv_in); + cudaFree(dv_out); + } + } +} diff --git a/stream_compaction/thrust.h b/stream_compaction/thrust.h new file mode 100644 index 0000000..fe98206 --- /dev/null +++ b/stream_compaction/thrust.h @@ -0,0 +1,11 @@ +#pragma once + +#include "common.h" + +namespace StreamCompaction { + namespace Thrust { + StreamCompaction::Common::PerformanceTimer& timer(); + + void scan(int n, int *odata, const int *idata); + } +}