diff --git a/10020311.PNG b/10020311.PNG new file mode 100755 index 0000000..cec8295 Binary files /dev/null and b/10020311.PNG differ diff --git a/10021224.PNG b/10021224.PNG new file mode 100755 index 0000000..21ff2c9 Binary files /dev/null and b/10021224.PNG differ diff --git a/10021245.PNG b/10021245.PNG new file mode 100755 index 0000000..4430f8f Binary files /dev/null and b/10021245.PNG differ diff --git a/10021259.PNG b/10021259.PNG new file mode 100755 index 0000000..decd0f3 Binary files /dev/null and b/10021259.PNG differ diff --git a/10021534.PNG b/10021534.PNG new file mode 100755 index 0000000..2fbe5ba Binary files /dev/null and b/10021534.PNG differ diff --git a/10021636.PNG b/10021636.PNG new file mode 100755 index 0000000..3b38878 Binary files /dev/null and b/10021636.PNG differ diff --git a/10021740.PNG b/10021740.PNG new file mode 100755 index 0000000..85b72a9 Binary files /dev/null and b/10021740.PNG differ diff --git a/10022047.PNG b/10022047.PNG new file mode 100755 index 0000000..2809988 Binary files /dev/null and b/10022047.PNG differ diff --git a/10022141.PNG b/10022141.PNG new file mode 100755 index 0000000..7165c6b Binary files /dev/null and b/10022141.PNG differ diff --git a/10022256.PNG b/10022256.PNG new file mode 100755 index 0000000..fbaa8bf Binary files /dev/null and b/10022256.PNG differ diff --git a/IMG_1929.MOV b/IMG_1929.MOV new file mode 100755 index 0000000..514f752 Binary files /dev/null and b/IMG_1929.MOV differ diff --git a/README.md b/README.md index 8f246fb..c9dd36f 100755 --- a/README.md +++ b/README.md @@ -5,50 +5,32 @@ Fall 2013 ------------------------------------------------------------------------------- Due Wednesday, 10/02/13 ------------------------------------------------------------------------------- - -------------------------------------------------------------------------------- -NOTE: +Qiong Wang ------------------------------------------------------------------------------- -This project requires an NVIDIA graphics card with CUDA capability! Any card after the Geforce 8xxx series will work. If you do not have an NVIDIA graphics card in the machine you are working on, feel free to use any machine in the SIG Lab or in Moore100 labs. All machines in the SIG Lab and Moore100 are equipped with CUDA capable NVIDIA graphics cards. If this too proves to be a problem, please contact Patrick or Liam as soon as possible. + ------------------------------------------------------------------------------- -INTRODUCTION: +INTRODUCTION ------------------------------------------------------------------------------- -In this project, you will extend your raytracer from Project 1 into a full CUDA based global illumination pathtracer. +This is a fast GPU path tracer initialized by a CUDA ray tracer with ray parallelization. It is one project of CIS 565 GPU Programming course at University of Pennsylvania. -For this project, you may either choose to continue working off of your codebase from Project 1, or you may choose to use the included basecode in this repository. The basecode for Project 2 is the same as the basecode for Project 1, but with some missing components you will need filled in, such as the intersection testing and camera raycasting methods. +The basic algorithm of the path tracing of this project can be summarized as: -How you choose to extend your raytracer into a pathtracer is a fairly open-ended problem; the supplied basecode is meant to serve as one possible set of guidelines for doing so, but you may choose any approach you want in your actual implementation, including completely scrapping the provided basecode in favor of your own from-scratch solution. +1. Build a pool of rays that need to be tested; +2. Construct an image that each color can be accumulated; +3. Launch a kernel to trace one bounce of a ray; +4. Add any new rays to the pool after bounce and check whether the ray hits the light source or not; +5. Remove terminated rays from the pool; +6. Repeat the third step until no ray left in the pool. -------------------------------------------------------------------------------- -CONTENTS: -------------------------------------------------------------------------------- -The Project2 root directory contains the following subdirectories: - -* src/ contains the source code for the project. Both the Windows Visual Studio solution and the OSX makefile reference this folder for all source; the base source code compiles on OSX and Windows without modification. -* scenes/ contains an example scene description file. -* renders/ contains two example renders: the raytraced render from Project 1 (GI_no.bmp), and the same scene rendered with global illumination (GI_yes.bmp). -* PROJ1_WIN/ contains a Windows Visual Studio 2010 project and all dependencies needed for building and running on Windows 7. -* PROJ1_OSX/ contains a OSX makefile, run script, and all dependencies needed for building and running on Mac OSX 10.8. -* PROJ1_NIX/ contains a Linux makefile for building and running on Ubuntu - 12.04 LTS. Note that you will need to set the following environment - variables: - - - PATH=$PATH:/usr/local/cuda-5.5/bin - - LD_LIBRARY_PATH=/usr/local/cuda-5.5/lib64:/lib +Since the number of the ray in the ray pool decreases we need fewer blocks per grid when launching the kernel. Hence, we can have a quite fast execution speed for the path tracer. -The projects build and run exactly the same way as in Project0 and Project1. +Note: Please run the 565Raytracer.sln file rather to run the 565Pathtracer.sln one, you can directly copy the 565Raytracer file from the Project 1 repo. Thank you :) ------------------------------------------------------------------------------- -REQUIREMENTS: +FEATURES IMPLEMENTED ------------------------------------------------------------------------------- -In this project, you are given code for: - -* All of the basecode from Project 1, plus: -* Intersection testing code for spheres and cubes -* Code for raycasting from the camera - -You will need to implement the following features. A number of these required features you may have already implemented in Project 1. If you have, you are ahead of the curve and have less work to do! +Basic features: * Full global illumination (including soft shadows, color bleeding, etc.) by pathtracing rays through the scene. * Properly accumulating emittance and colors to generate a final image @@ -56,92 +38,92 @@ You will need to implement the following features. A number of these required fe * Parallelization by ray instead of by pixel via string compaction * Perfect specular reflection -You are also required to implement at least two of the following features. Some of these features you may have already implemented in Project 1. If you have, you may NOT resubmit those features and instead must pick two new ones to implement. +Optional features: -* Additional BRDF models, such as Cook-Torrance, Ward, etc. Each BRDF model may count as a separate feature. -* Texture mapping -* Bump mapping * Translational motion blur -* Fresnel-based Refraction, i.e. glass -* OBJ Mesh loading and rendering without KD-Tree -* Interactive camera -* Integrate an existing stackless KD-Tree library, such as CUKD (https://github.com/unvirtual/cukd) * Depth of field +* Fresnel-based Refraction, i.e. glass (Still tuning) -Alternatively, implementing just one of the following features can satisfy the "pick two" feature requirement, since these are correspondingly more difficult problems: - -* Physically based subsurface scattering and transmission -* Implement and integrate your own stackless KD-Tree from scratch. -* Displacement mapping -* Deformational motion blur - -As yet another alternative, if you have a feature or features you really want to implement that are not on this list, let us know, and we'll probably say yes! ------------------------------------------------------------------------------- -NOTES ON GLM: +SCREENSHOTS OF THE FEATURES IMPLEMENTED ------------------------------------------------------------------------------- -This project uses GLM, the GL Math library, for linear algebra. You need to know two important points on how GLM is used in this project: +* Global illumination -* In this project, indices in GLM vectors (such as vec3, vec4), are accessed via swizzling. So, instead of v[0], v.x is used, and instead of v[1], v.y is used, and so on and so forth. -* GLM Matrix operations work fine on NVIDIA Fermi cards and later, but pre-Fermi cards do not play nice with GLM matrices. As such, in this project, GLM matrices are replaced with a custom matrix struct, called a cudaMat4, found in cudaMat4.h. A custom function for multiplying glm::vec4s and cudaMat4s is provided as multiplyMV() in intersections.h. +![ScreenShot](https://raw.github.com/GabriellaQiong/Project2-Pathtracer/master/10021534.PNG) -------------------------------------------------------------------------------- -README -------------------------------------------------------------------------------- -All students must replace or augment the contents of this Readme.md in a clear -manner with the following: +* Properly accumulating emittance and colors to generate a final image -* A brief description of the project and the specific features you implemented. -* At least one screenshot of your project running. -* A 30 second or longer video of your project running. To create the video you - can use http://www.microsoft.com/expression/products/Encoder4_Overview.aspx -* A performance evaluation (described in detail below). +![ScreenShot](https://raw.github.com/GabriellaQiong/Project2-Pathtracer/master/10021636.PNG) -------------------------------------------------------------------------------- -PERFORMANCE EVALUATION -------------------------------------------------------------------------------- -The performance evaluation is where you will investigate how to make your CUDA -programs more efficient using the skills you've learned in class. You must have -performed at least one experiment on your code to investigate the positive or -negative effects on performance. +* Antialiasing (Stochastic method) -One such experiment would be to investigate the performance increase involved -with adding a spatial data-structure to your scene data. +![ScreenShot](https://raw.github.com/GabriellaQiong/Project2-Pathtracer/master/anti-aliasing.PNG) -Another idea could be looking at the change in timing between various block -sizes. +* Perfect specular reflection + +![ScreenShot](https://raw.github.com/GabriellaQiong/Project2-Pathtracer/master/10021740.PNG) -A good metric to track would be number of rays per second, or frames per -second, or number of objects displayable at 60fps. +* Translational motion blur -We encourage you to get creative with your tweaks. Consider places in your code -that could be considered bottlenecks and try to improve them. +![ScreenShot](https://raw.github.com/GabriellaQiong/Project2-Pathtracer/master/10022256.PNG) + +* Depth of field + +With small radius for the circle of confusion + +![ScreenShot](https://raw.github.com/GabriellaQiong/Project2-Pathtracer/master/10022047.PNG) + +With big radius for the circle of confusion + +![ScreenShot](https://raw.github.com/GabriellaQiong/Project2-Pathtracer/master/10022141.PNG) -Each student should provide no more than a one page summary of their -optimizations along with tables and or graphs to visually explain any -performance differences. ------------------------------------------------------------------------------- -THIRD PARTY CODE POLICY +VIDEOS OF IMPLEMENTATION ------------------------------------------------------------------------------- -* Use of any third-party code must be approved by asking on the Google group. If it is approved, all students are welcome to use it. Generally, we approve use of third-party code that is not a core part of the project. For example, for the ray tracer, we would approve using a third-party library for loading models, but would not approve copying and pasting a CUDA function for doing refraction. -* Third-party code must be credited in README.md. -* Using third-party code without its approval, including using another student's code, is an academic integrity violation, and will result in you receiving an F for the semester. + +This is the video of the rendering process of the path tracer. + +[![ScreenShot](https://raw.github.com/GabriellaQiong/Project2-Pathtracer/master/videoScreenShot.PNG)](http://www.youtube.com/watch?v=GcbRUaLgz5A) + +The youtube link is here if you cannot open the video in the markdown file: http://www.youtube.com/watch?v=GcbRUaLgz5A + ------------------------------------------------------------------------------- -SELF-GRADING +PERFORMANCE EVALUATION ------------------------------------------------------------------------------- -* On the submission date, email your grade, on a scale of 0 to 100, to Liam, liamboone+cis565@gmail.com, with a one paragraph explanation. Be concise and realistic. Recall that we reserve 30 points as a sanity check to adjust your grade. Your actual grade will be (0.7 * your grade) + (0.3 * our grade). We hope to only use this in extreme cases when your grade does not realistically reflect your work - it is either too high or too low. In most cases, we plan to give you the exact grade you suggest. -* Projects are not weighted evenly, e.g., Project 0 doesn't count as much as the path tracer. We will determine the weighting at the end of the semester based on the size of each project. +Here is the table for the performance evaluation when changing the size of tile and with stream compaction or not. + +| tileSize | with compaction | time for running 10 iteration | approximate fps | +|:---------:|:-----------------:|:---------------------------------------:|:-----------------:| +| 8 | yes | 0 : 22.72 | 0.440 | +| 10 | yes | 0 : 25.97 | 0.385 | +| 8 | no | 0 : 23.39 | 0.428 | +| 10 | no | 0 : 27.47 | 0.364 | + +We can easily find that when the tile size become larger the fps decrease a little somehow. The fps with stream compaction is higher than that without stream compaction. ------------------------------------------------------------------------------- -SUBMISSION +REFERENCES ------------------------------------------------------------------------------- -As with the previous project, you should fork this project and work inside of your fork. Upon completion, commit your finished project back to your fork, and make a pull request to the master repository. -You should include a README.md file in the root directory detailing the following +* Snell's law: http://en.wikipedia.org/wiki/Snell%27s_law + +* Internal Reflection: http://en.wikipedia.org/wiki/Total_internal_reflection + +* Fresnel Equation: http://en.wikipedia.org/wiki/Fresnel_equations + +* Russian Roulette rule from Peter and Karl: https://docs.google.com/file/d/0B72qrSEH6CGfbFV0bGxmLVJiUlU/edit + +* General Phong Specular Model: http://en.wikipedia.org/wiki/Phong_reflection_model + +* Depth of Field: http://http.developer.nvidia.com/GPUGems/gpugems_ch23.html -* A brief description of the project and specific features you implemented -* At least one screenshot of your project running, and at least one screenshot of the final rendered output of your pathtracer -* Instructions for building and running your project if they differ from the base code -* A link to your blog post detailing the project -* A list of all third-party code used +* Stochastic Sampling: http://pages.cpsc.ucalgary.ca/~mario/courses/591-691-W06/PR/3-ray-tracing/3-advanced/readings/Cook_Stochastic_Sampling_TOG86.pdf + +* Usage of thrust in stream compaction: http://stackoverflow.com/questions/12201446/converting-thrustiterators-to-and-from-raw-pointers/12236270#12236270 + +------------------------------------------------------------------------------- +ACKNOWLEDGEMENT +------------------------------------------------------------------------------- +Thanks a lot to Patrick and Liam for the preparation of this project. And special thanks to Liam who helped us debugging so late till 4:30 am in the morning! Thank you all :) diff --git a/anti-aliasing.PNG b/anti-aliasing.PNG new file mode 100755 index 0000000..9ed1a25 Binary files /dev/null and b/anti-aliasing.PNG differ diff --git a/scenes/sampleScene.txt b/scenes/sampleScene.txt index 52d079e..d2ee373 100755 --- a/scenes/sampleScene.txt +++ b/scenes/sampleScene.txt @@ -1,229 +1,239 @@ -MATERIAL 0 //white diffuse -RGB 0.9 0.9 0.9 -SPECEX 0 -SPECRGB 1 1 1 -REFL 0 -REFR 0 -REFRIOR 0 -SCATTER 0 -ABSCOEFF 0 0 0 -RSCTCOEFF 0 -EMITTANCE 0 - -MATERIAL 1 //red diffuse -RGB .63 .06 .04 -SPECEX 0 -SPECRGB 1 1 1 -REFL 0 -REFR 0 -REFRIOR 0 -SCATTER 0 -ABSCOEFF 0 0 0 -RSCTCOEFF 0 -EMITTANCE 0 - -MATERIAL 2 //green diffuse -RGB .15 .48 .09 -SPECEX 0 -SPECRGB 1 1 1 -REFL 0 -REFR 0 -REFRIOR 0 -SCATTER 0 -ABSCOEFF 0 0 0 -RSCTCOEFF 0 -EMITTANCE 0 - -MATERIAL 3 //red glossy -RGB .63 .26 .24 -SPECEX 0 -SPECRGB 1 1 1 -REFL 0 -REFR 0 -REFRIOR 2 -SCATTER 0 -ABSCOEFF 0 0 0 -RSCTCOEFF 0 -EMITTANCE 0 - -MATERIAL 4 //white glossy -RGB 1 1 1 -SPECEX 0 -SPECRGB 1 1 1 -REFL 0 -REFR 0 -REFRIOR 2 -SCATTER 0 -ABSCOEFF 0 0 0 -RSCTCOEFF 0 -EMITTANCE 0 - -MATERIAL 5 //glass -RGB 0 0 0 -SPECEX 0 -SPECRGB 1 1 1 -REFL 0 -REFR 1 -REFRIOR 2.2 -SCATTER 0 -ABSCOEFF .02 5.1 5.7 -RSCTCOEFF 13 -EMITTANCE 0 - -MATERIAL 6 //green glossy -RGB .35 .48 .29 -SPECEX 0 -SPECRGB 1 1 1 -REFL 0 -REFR 0 -REFRIOR 2.6 -SCATTER 0 -ABSCOEFF 0 0 0 -RSCTCOEFF 0 -EMITTANCE 0 - -MATERIAL 7 //light -RGB 1 1 1 -SPECEX 0 -SPECRGB 0 0 0 -REFL 0 -REFR 0 -REFRIOR 0 -SCATTER 0 -ABSCOEFF 0 0 0 -RSCTCOEFF 0 -EMITTANCE 1 - -MATERIAL 8 //light -RGB 1 1 1 -SPECEX 0 -SPECRGB 0 0 0 -REFL 0 -REFR 0 -REFRIOR 0 -SCATTER 0 -ABSCOEFF 0 0 0 -RSCTCOEFF 0 -EMITTANCE 15 - -CAMERA -RES 800 800 -FOVY 25 -ITERATIONS 5000 -FILE test.bmp -frame 0 -EYE 0 4.5 12 -VIEW 0 0 -1 -UP 0 1 0 -frame 1 -EYE 0 4.5 12 -VIEW 0 0 -1 -UP 0 1 0 - -OBJECT 0 -cube -material 0 -frame 0 -TRANS 0 0 0 -ROTAT 0 0 90 -SCALE .01 10 10 -frame 1 -TRANS 0 0 0 -ROTAT 0 0 90 -SCALE .01 10 10 - -OBJECT 1 -cube -material 0 -frame 0 -TRANS 0 5 -5 -ROTAT 0 90 0 -SCALE .01 10 10 -frame 1 -TRANS 0 5 -5 -ROTAT 0 90 0 -SCALE .01 10 10 - -OBJECT 2 -cube -material 0 -frame 0 -TRANS 0 10 0 -ROTAT 0 0 90 -SCALE .01 10 10 -frame 1 -TRANS 0 10 0 -ROTAT 0 0 90 -SCALE .01 10 10 - -OBJECT 3 -cube -material 1 -frame 0 -TRANS -5 5 0 -ROTAT 0 0 0 -SCALE .01 10 10 -frame 1 -TRANS -5 5 0 -ROTAT 0 0 0 -SCALE .01 10 10 - -OBJECT 4 -cube -material 2 -frame 0 -TRANS 5 5 0 -ROTAT 0 0 0 -SCALE .01 10 10 -frame 1 -TRANS 5 5 0 -ROTAT 0 0 0 -SCALE .01 10 10 - -OBJECT 5 -sphere -material 4 -frame 0 -TRANS 0 2 0 -ROTAT 0 180 0 -SCALE 3 3 3 -frame 1 -TRANS 0 2 0 -ROTAT 0 180 0 -SCALE 3 3 3 - -OBJECT 6 -sphere -material 3 -frame 0 -TRANS 2 5 2 -ROTAT 0 180 0 -SCALE 2.5 2.5 2.5 -frame 1 -TRANS 2 5 2 -ROTAT 0 180 0 -SCALE 2.5 2.5 2.5 - -OBJECT 7 -sphere -material 6 -frame 0 -TRANS -2 5 -2 -ROTAT 0 180 0 -SCALE 3 3 3 -frame 1 -TRANS -2 5 -2 -ROTAT 0 180 0 -SCALE 3 3 3 - -OBJECT 8 -cube -material 8 -frame 0 -TRANS 0 10 0 -ROTAT 0 0 90 -SCALE .3 3 3 -frame 1 -TRANS 0 10 0 -ROTAT 0 0 90 -SCALE .3 3 3 +MATERIAL 0 //white diffuse +RGB 1 1 1 +SPECEX 0 +SPECRGB 1 1 1 +REFL 0 +REFR 0 +REFRIOR 0 +SCATTER 0 +ABSCOEFF 0 0 0 +RSCTCOEFF 0 +EMITTANCE 0 + +MATERIAL 1 //red diffuse +RGB .63 .06 .04 +SPECEX 0 +SPECRGB 1 1 1 +REFL 0 +REFR 0 +REFRIOR 0 +SCATTER 0 +ABSCOEFF 0 0 0 +RSCTCOEFF 0 +EMITTANCE 0 + +MATERIAL 2 //green diffuse +RGB .15 .48 .09 +SPECEX 0 +SPECRGB 1 1 1 +REFL 0 +REFR 0 +REFRIOR 0 +SCATTER 0 +ABSCOEFF 0 0 0 +RSCTCOEFF 0 +EMITTANCE 0 + +MATERIAL 3 //red glossy +RGB .63 .06 .04 +SPECEX 500 +SPECRGB 1 1 1 +REFL 0 +REFR 0 +REFRIOR 2 +SCATTER 0 +ABSCOEFF 0 0 0 +RSCTCOEFF 0 +EMITTANCE 0 + +MATERIAL 4 //white glossy +RGB 1 1 1 +SPECEX 500 +SPECRGB 1 1 1 +REFL 0 +REFR 0 +REFRIOR 2 +SCATTER 0 +ABSCOEFF 0 0 0 +RSCTCOEFF 0 +EMITTANCE 0 + +MATERIAL 5 //glass +RGB 0 0 0 +SPECEX 0 +SPECRGB 1 1 1 +REFL 0 +REFR 1 +REFRIOR 2.2 +SCATTER 0 +ABSCOEFF .02 5.1 5.7 +RSCTCOEFF 13 +EMITTANCE 0 + +MATERIAL 6 //green glossy +RGB .15 .48 .09 +SPECEX 500 +SPECRGB 1 1 1 +REFL 0 +REFR 0 +REFRIOR 2.6 +SCATTER 0 +ABSCOEFF 0 0 0 +RSCTCOEFF 0 +EMITTANCE 0 + +MATERIAL 7 //light +RGB 1 1 1 +SPECEX 0 +SPECRGB 0 0 0 +REFL 0 +REFR 0 +REFRIOR 0 +SCATTER 0 +ABSCOEFF 0 0 0 +RSCTCOEFF 0 +EMITTANCE 1 + +MATERIAL 8 //light +RGB 1 1 1 +SPECEX 0 +SPECRGB 0 0 0 +REFL 0 +REFR 0 +REFRIOR 0 +SCATTER 0 +ABSCOEFF 0 0 0 +RSCTCOEFF 0 +EMITTANCE 15 + +MATERIAL 9 //mirror +RGB 0 0 0 +SPECEX 0 +SPECRGB 1 1 1 +REFL 1 +REFR 0 +REFRIOR 1.5 +SCATTER 0 +ABSCOEFF .02 5.1 5.7 +RSCTCOEFF 13 +EMITTANCE 0 + + +MATERIAL 10 //blue diffuse +RGB .13 .16 .84 +SPECEX 500 +SPECRGB 1 1 1 +REFL 0 +REFR 0 +REFRIOR 1.0 +SCATTER 0 +ABSCOEFF 0 0 0 +RSCTCOEFF 0 +EMITTANCE 0 + +MATERIAL 11 //purple diffuse +RGB .43 .16 .54 +SPECEX 0 +SPECRGB 1 1 1 +REFL 0 +REFR 0 +REFRIOR 0 +SCATTER 0 +ABSCOEFF 0 0 0 +RSCTCOEFF 0 +EMITTANCE 0 + +MATERIAL 12 //yellow specular +RGB .5 .5 .05 +SPECEX 500 +SPECRGB 1 1 1 +REFL 0 +REFR 0 +REFRIOR 1.5 +SCATTER 0 +ABSCOEFF 0 0 0 +RSCTCOEFF 0 +EMITTANCE 0 + +CAMERA +RES 800 800 +FOVY 25 +ITERATIONS 5000 +FILE test.bmp +frame 0 +EYE 0 4.5 12 +VIEW 0 0 -1 +UP 0 1 0 + +OBJECT 0 +cube +material 0 +frame 0 +TRANS 0 0 0 +ROTAT 0 0 90 +SCALE .01 10 10 + +OBJECT 1 +cube +material 0 +frame 0 +TRANS 0 5 -5 +ROTAT 0 90 0 +SCALE .01 10 10 + +OBJECT 2 +cube +material 0 +frame 0 +TRANS 0 10 0 +ROTAT 0 0 90 +SCALE .01 10 10 + +OBJECT 3 +cube +material 1 +frame 0 +TRANS -5 5 0 +ROTAT 0 0 0 +SCALE .01 10 10 + +OBJECT 4 +cube +material 2 +frame 0 +TRANS 5 5 0 +ROTAT 0 0 0 +SCALE .01 10 10 + +OBJECT 5 +sphere +material 4 +frame 0 +TRANS 0 2 0 +ROTAT 0 180 0 +SCALE 3 3 3 + +OBJECT 6 +sphere +material 3 +frame 0 +TRANS 2 5 2 +ROTAT 0 180 0 +SCALE 2.5 2.5 2.5 + +OBJECT 7 +sphere +material 10 +frame 0 +TRANS -2 5 -2 +ROTAT 0 180 0 +SCALE 3 3 3 + + +OBJECT 8 +cube +material 8 +frame 0 +TRANS 0 10 0 +ROTAT 0 0 90 +SCALE .3 3 3 \ No newline at end of file diff --git a/src/interactions.h b/src/interactions.h index a09ec95..68f3c05 100755 --- a/src/interactions.h +++ b/src/interactions.h @@ -8,6 +8,16 @@ #include "intersections.h" +struct Fresnel { + float reflectionCoefficient; + float transmissionCoefficient; +}; + +struct AbsorptionAndScatteringProperties{ + glm::vec3 absorptionCoefficient; + float reducedScatteringCoefficient; +}; + //forward declaration __host__ __device__ bool calculateScatterAndAbsorption(ray& r, float& depth, AbsorptionAndScatteringProperties& currentAbsorptionAndScattering, glm::vec3& unabsorbedColor, material m, float randomFloatForScatteringDistance, float randomFloat2, float randomFloat3); __host__ __device__ glm::vec3 getRandomDirectionInSphere(float xi1, float xi2); @@ -17,39 +27,212 @@ __host__ __device__ glm::vec3 calculateReflectionDirection(glm::vec3 normal, glm __host__ __device__ Fresnel calculateFresnel(glm::vec3 normal, glm::vec3 incident, float incidentIOR, float transmittedIOR, glm::vec3 reflectionDirection, glm::vec3 transmissionDirection); __host__ __device__ glm::vec3 calculateRandomDirectionInHemisphere(glm::vec3 normal, float xi1, float xi2); +//TODO (OPTIONAL): IMPLEMENT THIS FUNCTION +__host__ __device__ glm::vec3 calculateTransmission(glm::vec3 absorptionCoefficient, float distance) { + return glm::vec3(0,0,0); +} + +//TODO (OPTIONAL): IMPLEMENT THIS FUNCTION +__host__ __device__ bool calculateScatterAndAbsorption(ray& r, float& depth, AbsorptionAndScatteringProperties& currentAbsorptionAndScattering, + glm::vec3& unabsorbedColor, material m, float randomFloatForScatteringDistance, float randomFloat2, float randomFloat3){ + return false; +} + +//TODO (OPTIONAL): IMPLEMENT THIS FUNCTION -- This is Done but no use all merged with calculateFresnel it's easy using to put together +__host__ __device__ glm::vec3 calculateTransmissionDirection(glm::vec3 normal, glm::vec3 incident, float incidentIOR, float transmittedIOR) { + // According to Snell's law (http://en.wikipedia.org/wiki/Snell%27s_law), calculate the cosine of incident angle and transmitted angle + float ratio = incidentIOR / transmittedIOR; + float cosIncident = glm::dot(normal, - incident); + float cosTransmissionSquare = 1 - ratio * ratio * (1 - cosIncident * cosIncident); + + // Determine whether the square is negative or positive, if negative there is internal reflection + if (cosTransmissionSquare < 0) { + // When the transmitted angle reaches 90 degree the critical angle, there is no light to transmit and all lights reflect internally referring to http://en.wikipedia.org/wiki/Total_internal_reflection + return calculateReflectionDirection(normal, incident); + } + + // Calculate the transmission direction when the cosine of incident angle is negative change the formula + float cosTransmission = glm::sqrt(cosTransmissionSquare); + glm::vec3 transmissionDirection; + if (cosIncident > 0.0f) + transmissionDirection = ratio * incident + (ratio * cosIncident - cosTransmission) * normal; + else + transmissionDirection = ratio * incident + (ratio * cosIncident + cosTransmission) * normal; + return glm::normalize(transmissionDirection); +} + +//TODO (OPTIONAL): IMPLEMENT THIS FUNCTION -- This is Done +__host__ __device__ glm::vec3 calculateReflectionDirection(glm::vec3 normal, glm::vec3 incident) { + // nothing fancy here (simple vector addition) + return (glm::normalize(incident - 2.0f * glm::dot(incident, normal) * normal)); +} + +//TODO (OPTIONAL): IMPLEMENT THIS FUNCTION -- This is Done +__host__ __device__ Fresnel calculateFresnel(glm::vec3 normal, glm::vec3 incident, float incidentIOR, float transmittedIOR, glm::vec3 reflectionDirection, glm::vec3 transmissionDirection) { + + Fresnel fresnel; + float n1 = incidentIOR; + float n2 = transmittedIOR; + + // Return the default value for reflection when n1 or n2 is less than or equal zero + if (n1 <= EPSILON && n2 <= EPSILON) { + fresnel.reflectionCoefficient = 1.0f; + fresnel.transmissionCoefficient = 0.0f; + return fresnel; + } + + // Merged part detailed description in calculateTransmissionDirection + float ratio = n1 / n2; + float cosIncident = glm::dot(normal, - incident); + float cosTransmissionSquare = 1 - ratio * ratio * (1 - cosIncident * cosIncident); + + if (cosTransmissionSquare < 0) { + reflectionDirection = calculateReflectionDirection(normal, incident); + fresnel.reflectionCoefficient = 1.0f; + fresnel.transmissionCoefficient = 0.0f; + return fresnel; + } + + // Calculate the transmission direction when the cosine of incident angle is negative change the formula + float cosTransmission = glm::sqrt(cosTransmissionSquare); + if (cosIncident > 0.0f) + transmissionDirection = ratio * incident + (ratio * cosIncident - cosTransmission) * normal; + else + transmissionDirection = ratio * incident + (ratio * cosIncident + cosTransmission) * normal; + + // Fresnel Equation according to http://en.wikipedia.org/wiki/Fresnel_equations + // Reflection coefficient of s-polarized light + float Rs = 0; + if (!epsilonCheck(n1 * cosIncident + n2 * cosTransmission, 0)) { + Rs = glm::max(pow((n1 * cosIncident - n2 * cosTransmission) / (n1 * cosIncident + n2 * cosTransmission), 2.0f), 0.0f); + } + + // Reflection coefficient of p-polarized light + float Rp = 0; + if (!epsilonCheck(n1 * cosTransmission + n2 * cosIncident, 0)) { + Rp = glm::max(pow((n1 * cosTransmission - n2 * cosIncident) / (n1 * cosTransmission + n2 * cosIncident), 2.0f) , 0.0f); + } + + fresnel.reflectionCoefficient = (Rs + Rp) / 2; + fresnel.transmissionCoefficient = 1 - fresnel.reflectionCoefficient; + reflectionDirection = calculateReflectionDirection(normal, incident); + return fresnel; +} + //LOOK: This function demonstrates cosine weighted random direction generation in a sphere! __host__ __device__ glm::vec3 calculateRandomDirectionInHemisphere(glm::vec3 normal, float xi1, float xi2) { + // PERSONAL NOTES: the cosine weighted random direction is on the up direction normal to the hemisphere + //crucial difference between this and calculateRandomDirectionInSphere: THIS IS COSINE WEIGHTED! - //crucial difference between this and calculateRandomDirectionInSphere: THIS IS COSINE WEIGHTED! + float up = sqrt(xi1); // cos(theta) + float over = sqrt(1 - up * up); // sin(theta) + float around = xi2 * TWO_PI; - float up = sqrt(xi1); // cos(theta) - float over = sqrt(1 - up * up); // sin(theta) - float around = xi2 * 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. - //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 = glm::vec3(0); + 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); + } - 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 ); + //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 ); } //TODO: IMPLEMENT THIS FUNCTION //Now that you know how cosine weighted direction generation works, try implementing non-cosine (uniform) weighted random direction generation. //This should be much easier than if you had to implement calculateRandomDirectionInHemisphere. __host__ __device__ glm::vec3 getRandomDirectionInSphere(float xi1, float xi2) { - return glm::vec3(0,0,0); + // xi1, xi2 in [0, 1) the uniform weighted random direction means each direction has equal probability of being selected + // Reffering to Emmanuel Agu's slides for the difference between cosine-weighted and uniform sampling + float r = sqrt(1.0f - xi1); + float theta = xi2 * TWO_PI; + return glm::vec3(r * cos(theta), r * sin(theta), xi1); } +//TODO (PARTIALLY OPTIONAL): IMPLEMENT THIS FUNCTION +//returns 0 if diffuse scatter, 1 if reflected, 2 if transmitted. +__host__ __device__ int calculateBSDF(ray currentRay, rayPool &nextRay, glm::vec3 intersect, glm::vec3 normal, glm::vec3 emittedColor, + glm::vec3& color, material m, float time, cameraData cam){ + nextRay.ray.origin = intersect; + glm::vec3 reflectionDirection,transmissionDirection; + Fresnel fresnel; + int classifier = -1; // 0 -- diffusion, 1 -- reflection, 2 -- transmission + + // Generate the random numbers + thrust::default_random_engine rng(hash(time)); + thrust::uniform_real_distribution u01(0,1); + float randomNumber = u01(rng); + + // Use the refractive or reflective parameter determining whether the surface is diffuse or reflected or transmitted + if (m.hasReflective && m.hasRefractive) { + // When the light is incidented from medium with IOR larger than 1 to the air there might exist internal reflection + // Using Fresenel equation to find the coefficients both of refraction and reflection + if (glm::dot(normal, currentRay.direction) < EPSILON) + fresnel = calculateFresnel(normal, currentRay.direction, 1.0f, m.indexOfRefraction, reflectionDirection, transmissionDirection); + else + fresnel = calculateFresnel(-normal, currentRay.direction, m.indexOfRefraction, 1.0f, reflectionDirection, transmissionDirection); + + // Using Peter and Karl's method https://docs.google.com/file/d/0B72qrSEH6CGfbFV0bGxmLVJiUlU/edit, do Russian Roulette to choose whether reflection or refraction + if (randomNumber > 0.5) { + classifier = 1; + } else { + classifier = 2; + } + } else if (m.hasReflective) { + classifier = 1; + } else if (m.hasRefractive) { + classifier = 2; + } else { + // Diffision + classifier = 0; + } // if hasReflective && hasReflective + + if (classifier == 0) { + // Diffusion + nextRay.ray.direction = calculateRandomDirectionInHemisphere(normal, randomNumber, 100 * randomNumber); // SO TRICKY here, increase the random angle 100 times the defects disappear!!! + color = m.color * emittedColor; + + // Phong lighting for specular (general one not the Blinn-Phong lighting) http://en.wikipedia.org/wiki/Phong_reflection_model + if( m.specularExponent > EPSILON ){ + float specularCoefficient = 6.0f; + glm::vec3 V = glm::normalize(cam.position-intersect); + glm::vec3 lightReflection=calculateReflectionDirection(normal,-nextRay.ray.direction); + float specular = glm::pow(glm::max(glm::dot(V,lightReflection), 0.0f),m.specularExponent); + color += specularCoefficient * specular * m.specularColor; + } + + return 0; + } else if (classifier == 1) { + // Reflection + nextRay.ray.direction = reflectionDirection; + nextRay.coefficient = fresnel.reflectionCoefficient; + color = emittedColor * nextRay.coefficient; + + // Phong lighting for specular (general one not the Blinn-Phong lighting) http://en.wikipedia.org/wiki/Phong_reflection_model + if( m.specularExponent > EPSILON ){ + float specularCoefficient = 6.0f; + glm::vec3 V = glm::normalize(cam.position-intersect); + glm::vec3 lightReflection=calculateReflectionDirection(normal,-nextRay.ray.direction); + float specular = glm::pow(glm::max(glm::dot(V,lightReflection), 0.0f),m.specularExponent); + color += specularCoefficient * specular * m.specularColor; + } + return 1; + } else if (classifier == 2) { + // Transmission + nextRay.ray.direction = transmissionDirection; + nextRay.coefficient = fresnel.transmissionCoefficient; + color = emittedColor * nextRay.coefficient; + return 2; + } // if classifier + return -1; +}; + #endif diff --git a/src/intersections.h b/src/intersections.h index a6b9469..52cbde9 100755 --- a/src/intersections.h +++ b/src/intersections.h @@ -17,8 +17,7 @@ __host__ __device__ glm::vec3 getPointOnRay(ray r, float t); __host__ __device__ glm::vec3 multiplyMV(cudaMat4 m, glm::vec4 v); __host__ __device__ glm::vec3 getSignOfRay(ray r); __host__ __device__ glm::vec3 getInverseDirectionOfRay(ray r); -__host__ __device__ float boxIntersectionTest(staticGeom sphere, ray r, glm::vec3& intersectionPoint, glm::vec3& normal); -__host__ __device__ float boxIntersectionTest(glm::vec3 boxMin, glm::vec3 boxMax, staticGeom box, ray r, glm::vec3& intersectionPoint, glm::vec3& normal); +__host__ __device__ float boxIntersectionTest(staticGeom box, ray r, glm::vec3& intersectionPoint, glm::vec3& normal); __host__ __device__ float sphereIntersectionTest(staticGeom sphere, ray r, glm::vec3& intersectionPoint, glm::vec3& normal); __host__ __device__ glm::vec3 getRandomPointOnCube(staticGeom cube, float randomSeed); @@ -44,7 +43,8 @@ __host__ __device__ bool epsilonCheck(float a, float b){ //Self explanatory __host__ __device__ glm::vec3 getPointOnRay(ray r, float t){ - return r.origin + float(t-.0001)*glm::normalize(r.direction); + // Add a big epsilon amount .001 forward along the shadow ray (.0001 previously) so as to avoid floating problem + return r.origin + float(t-.001)*glm::normalize(r.direction); } //LOOK: This is a custom function for multiplying cudaMat4 4x4 matrixes with vectors. @@ -69,98 +69,129 @@ __host__ __device__ glm::vec3 getSignOfRay(ray r){ return glm::vec3((int)(inv_direction.x < 0), (int)(inv_direction.y < 0), (int)(inv_direction.z < 0)); } -//Wrapper for cube intersection test for testing against unit cubes -__host__ __device__ float boxIntersectionTest(staticGeom box, ray r, glm::vec3& intersectionPoint, glm::vec3& normal){ - return boxIntersectionTest(glm::vec3(-.5,-.5,-.5), glm::vec3(.5,.5,.5), box, r, intersectionPoint, normal); -} - +//TODO: IMPLEMENT THIS FUNCTION //Cube intersection test, return -1 if no intersection, otherwise, distance to intersection -__host__ __device__ float boxIntersectionTest(glm::vec3 boxMin, glm::vec3 boxMax, staticGeom box, ray r, glm::vec3& intersectionPoint, glm::vec3& normal){ - glm::vec3 currentNormal = glm::vec3(0,0,0); - - ray ro = r; - - glm::vec3 iP0 = multiplyMV(box.inverseTransform,glm::vec4(r.origin, 1.0f)); - glm::vec3 iP1 = multiplyMV(box.inverseTransform,glm::vec4(r.origin+r.direction, 1.0f)); - glm::vec3 iV0 = iP1 - iP0; - - r.origin = iP0; - r.direction = glm::normalize(iV0); - - float tmin, tmax, tymin, tymax, tzmin, tzmax; - - glm::vec3 rsign = getSignOfRay(r); - glm::vec3 rInverseDirection = getInverseDirectionOfRay(r); - - if((int)rsign.x==0){ - tmin = (boxMin.x - r.origin.x) * rInverseDirection.x; - tmax = (boxMax.x - r.origin.x) * rInverseDirection.x; - }else{ - tmin = (boxMax.x - r.origin.x) * rInverseDirection.x; - tmax = (boxMin.x - r.origin.x) * rInverseDirection.x; - } - - if((int)rsign.y==0){ - tymin = (boxMin.y - r.origin.y) * rInverseDirection.y; - tymax = (boxMax.y - r.origin.y) * rInverseDirection.y; - }else{ - tymin = (boxMax.y - r.origin.y) * rInverseDirection.y; - tymax = (boxMin.y - r.origin.y) * rInverseDirection.y; - } - - if ( (tmin > tymax) || (tymin > tmax) ){ - return -1; - } - if (tymin > tmin){ - tmin = tymin; - } - if (tymax < tmax){ - tmax = tymax; - } - - if((int)rsign.z==0){ - tzmin = (boxMin.z - r.origin.z) * rInverseDirection.z; - tzmax = (boxMax.z - r.origin.z) * rInverseDirection.z; - }else{ - tzmin = (boxMax.z - r.origin.z) * rInverseDirection.z; - tzmax = (boxMin.z - r.origin.z) * rInverseDirection.z; - } - - if ( (tmin > tzmax) || (tzmin > tmax) ){ - return -1; - } - if (tzmin > tmin){ - tmin = tzmin; - } - if (tzmax < tmax){ - tmax = tzmax; - } - if(tmin<0){ - return -1; - } - - glm::vec3 osintersect = r.origin + tmin*r.direction; - - if(abs(osintersect.x-abs(boxMax.x))<.001){ - currentNormal = glm::vec3(1,0,0); - }else if(abs(osintersect.y-abs(boxMax.y))<.001){ - currentNormal = glm::vec3(0,1,0); - }else if(abs(osintersect.z-abs(boxMax.z))<.001){ - currentNormal = glm::vec3(0,0,1); - }else if(abs(osintersect.x+abs(boxMin.x))<.001){ - currentNormal = glm::vec3(-1,0,0); - }else if(abs(osintersect.y+abs(boxMin.y))<.001){ - currentNormal = glm::vec3(0,-1,0); - }else if(abs(osintersect.z+abs(boxMin.z))<.001){ - currentNormal = glm::vec3(0,0,-1); - } - - intersectionPoint = multiplyMV(box.transform, glm::vec4(osintersect, 1.0)); - +__host__ __device__ float boxIntersectionTest(staticGeom box, ray r, glm::vec3& intersectionPoint, glm::vec3& normal){ + + // Find the inverse transformed ray in the origin centered coordinates + glm::vec3 ro = multiplyMV(box.inverseTransform, glm::vec4(r.origin, 1.0f)); + glm::vec3 rd = glm::normalize(multiplyMV(box.inverseTransform, glm::vec4(r.direction, 0.0f))); + ray rt; rt.origin = ro; rt.direction = rd; + // Find the inverse direction ray avoiding -/+ 0 determination + glm::vec3 inverseDirection = getInverseDirectionOfRay(rt); + glm::vec3 signInverseDirection = getSignOfRay(rt); + + // Initialize the minimum and maximum distance to intersection to compare + float tmin, tmax, minLimit, maxLimit; + minLimit = -1e7; + maxLimit = 1e7; + tmin = minLimit; + tmax = maxLimit; + + // Find the bounding volumn's maximum and minimum extent t1/t0 + glm::vec3 t0, t1; + t0.x = (0.5 * (2 * signInverseDirection.x - 1) - ro.x) * inverseDirection.x; + t1.x = (0.5 * (1 - 2 * signInverseDirection.x) - ro.x) * inverseDirection.x; + t0.y = (0.5 * (2 * signInverseDirection.y - 1) - ro.y) * inverseDirection.y; + t1.y = (0.5 * (1 - 2 * signInverseDirection.y) - ro.y) * inverseDirection.y; + t0.z = (0.5 * (2 * signInverseDirection.z - 1) - ro.z) * inverseDirection.z; + t1.z = (0.5 * (1 - 2 * signInverseDirection.z) - ro.z) * inverseDirection.z; + + // Surface flags to point out which is the nearest and farthest surface along the intersection direction + int nearFlag, farFlag; // which can be one value in pool of {1, 2, 3}, respectively indicating two surfaces along x || y || z axis + + // Compare maxima and minima to determine the intersections + if(t0.x < 0 && t1.x < 0) + return -1; // whether there is only inverse intersection + + if(t0.x > 0) { + tmin = (tmin > t0.x) ? tmin : t0.x; + if(tmin == t0.x) + nearFlag = 1; // update tmin if t0.x is positive and finite and set x to be the nearest surface axis + } + if(t1.x > 0) { + tmax = (tmax < t1.x) ? tmax : t1.x; + if(tmax == t1.x) + farFlag = 1; // update tmax if t1.x is positive and finite and set x to be the nearest surface axis + } + + // Then compare the current tmin and tmax with y, update + if(t0.y < 0 && t1.y < 0) + return -1; + if(tmin > t1.y || t0.y > tmax) + return -1; // whether the ray passes around the cube + + if(t0.y > 0) { + tmin = (tmin > t0.y) ? tmin : t0.y; + if(tmin == t0.y) + nearFlag = 2; // update tmin if t0.y is positive and larger than tmin and set y to be the nearest surface axis + } + if(t1.y > 0) { + tmax = (tmax < t1.y) ? tmax : t1.y; + if(tmax == t1.y) + farFlag = 2; // update tmax if t1.y is positive and less than tmax and set y to be the nearest surface axis + } + + // Then compare the current tmin and tmax with z, update + if(t0.z < 0 && t1.z < 0) + return -1; + if(tmin > t1.z || t0.z > tmax) + return -1; + + if(t0.z > 0) { + tmin = (tmin > t0.z) ? tmin : t0.z; + if(tmin == t0.z) + nearFlag = 3; // update tmin if t0.z is positive and larger than tmin and set z to be the nearest surface axis + } + if(t1.z > 0 && t1.z < tmax) { + tmax = (tmax < t1.z) ? tmax : t1.z; + if(tmax == t1.z) + farFlag = 3; // update tmax if t1.z is positive and less than tmax and set z to be the nearest surface axis + } + + glm::vec3 intersectionPointTmp, normalTmp; + // Update the normal and intersection point if the tmax or tmax are not all infinity + if(epsilonCheck(tmin, minLimit) == 0) { + intersectionPointTmp = getPointOnRay(rt, tmin); + if(nearFlag == 1) { + normalTmp = glm::vec3(-1, 0, 0); + if(signInverseDirection.x) + normalTmp.x = 1; + } else if(nearFlag == 2) { + normalTmp = glm::vec3(0, -1, 0); + if(signInverseDirection.y) + normalTmp.y = 1; + } else if(nearFlag ==3){ + normalTmp = glm::vec3(0, 0, -1); + if(signInverseDirection.z) + normalTmp.z = 1; + } + } else if(epsilonCheck(tmax, maxLimit) == 0) { + intersectionPointTmp = getPointOnRay(rt, tmax); + if(farFlag == 1) { + normalTmp = glm::vec3(-1, 0, 0); + if(signInverseDirection.x) + normalTmp.x = 1; + } else if(farFlag == 2) { + normalTmp = glm::vec3(0, -1, 0); + if(signInverseDirection.y) + normalTmp.y = 1; + } else { + normalTmp = glm::vec3(0, 0, -1); + if(signInverseDirection.z) + normalTmp.z = 1; + } + } else { + return -1; + } + + // Find the distance between real intersection point and transform back the normal + glm::vec3 realIntersectionPoint = multiplyMV(box.transform, glm::vec4(intersectionPointTmp, 1.0f)); + intersectionPoint = realIntersectionPoint; + normal = glm::normalize(multiplyMV(box.transform, glm::vec4(normalTmp, 0.0f))); + return glm::length(r.origin - realIntersectionPoint); - normal = multiplyMV(box.transform, glm::vec4(currentNormal,0.0)); - return glm::length(intersectionPoint-ro.origin); } //LOOK: Here's an intersection test example from a sphere. Now you just need to figure out cube and, optionally, triangle. @@ -189,9 +220,9 @@ __host__ __device__ float sphereIntersectionTest(staticGeom sphere, ray r, glm:: if (t1 < 0 && t2 < 0) { return -1; } else if (t1 > 0 && t2 > 0) { - t = min(t1, t2); + t = glm::min(t1, t2); // min } else { - t = max(t1, t2); + t = glm::max(t1, t2); // max } glm::vec3 realIntersectionPoint = multiplyMV(sphere.transform, glm::vec4(getPointOnRay(rt, t), 1.0)); @@ -261,20 +292,25 @@ __host__ __device__ glm::vec3 getRandomPointOnCube(staticGeom cube, float random } +//TODO: IMPLEMENT THIS FUNCTION //Generates a random point on a given sphere __host__ __device__ glm::vec3 getRandomPointOnSphere(staticGeom sphere, float randomSeed){ - float radius=.5f; - thrust::default_random_engine rng(hash(randomSeed)); - thrust::uniform_real_distribution u01(-1,1); - thrust::uniform_real_distribution u02(0,TWO_PI); - - float theta = (float)u02(rng); - float cosphi = (float)u01(rng); - float sinphi = sqrt(1 - cosphi*cosphi); - glm::vec3 point = radius*glm::vec3(sinphi*cos(theta),sinphi*sin(theta),cosphi); - glm::vec3 randPoint = multiplyMV(sphere.transform, glm::vec4(point,1.0f)); - - return randPoint; + + // Generate random number + thrust::default_random_engine rng(hash(randomSeed)); + thrust::uniform_real_distribution u01(0,TWO_PI); + thrust::uniform_real_distribution u02(0,PI); + + // Two independent variable expressing the point + float theta, phi; + theta = (float)u01(rng); + phi = (float)u02(rng); + + // Transform back to the real sphere coordinates + glm::vec3 point(cos(theta) * sin(phi), sin(theta) * sin(phi), cos(phi)); + glm::vec3 randPoint = multiplyMV(sphere.transform, glm::vec4(point,1.0f)); + + return randPoint; } #endif diff --git a/src/main.cpp b/src/main.cpp index 81836b1..ed08403 100755 --- a/src/main.cpp +++ b/src/main.cpp @@ -142,8 +142,8 @@ void runCuda(){ gammaSettings gamma; gamma.applyGamma = true; - gamma.gamma = 1.0; - gamma.divisor = 1.0; //renderCam->iterations; + gamma.gamma = 1.0/2.2; + gamma.divisor = 1;//renderCam->iterations; outputImage.setGammaSettings(gamma); string filename = renderCam->imageName; string s; @@ -201,7 +201,7 @@ void runCuda(){ void display(){ runCuda(); - string title = "565Raytracer | " + utilityCore::convertIntToString(iterations) + " Iterations"; + string title = "565Raytracer | Qiong Wang " + utilityCore::convertIntToString(iterations) + " Iterations"; glutSetWindowTitle(title.c_str()); glBindBuffer( GL_PIXEL_UNPACK_BUFFER, pbo); diff --git a/src/raytraceKernel.cu b/src/raytraceKernel.cu index 87a65a6..895d113 100755 --- a/src/raytraceKernel.cu +++ b/src/raytraceKernel.cu @@ -1,28 +1,37 @@ // CIS565 CUDA Raytracer: A parallel raytracer for Patrick Cozzi's CIS565: GPU Computing at the University of Pennsylvania // Written by Yining Karl Li, Copyright (c) 2012 University of Pennsylvania // This file includes code from: -// Rob Farber for CUDA-GL interop, from CUDA Supercomputing For The Masses: http://www.drdobbs.com/architecture-and-design/cuda-supercomputing-for-the-masses-part/222600097 -// Peter Kutz and Yining Karl Li's GPU Pathtracer: http://gpupathtracer.blogspot.com/ -// Yining Karl Li's TAKUA Render, a massively parallel pathtracing renderer: http://www.yiningkarlli.com +// Rob Farber for CUDA-GL interop, from CUDA Supercomputing For The Masses: http://www.drdobbs.com/architecture-and-design/cuda-supercomputing-for-the-masses-part/222600097 +// Peter Kutz and Yining Karl Li's GPU Pathtracer: http://gpupathtracer.blogspot.com/ +// Yining Karl Li's TAKUA Render, a massively parallel pathtracing renderer: http://www.yiningkarlli.com #include #include #include #include "sceneStructs.h" +#include "glm/glm.hpp" #include "utilities.h" #include "raytraceKernel.h" #include "intersections.h" #include "interactions.h" #include -#include "glm/glm.hpp" +#include +#include +#include + +#if CUDA_VERSION >= 5000 + #include +#else + #include +#endif void checkCUDAError(const char *msg) { cudaError_t err = cudaGetLastError(); if( cudaSuccess != err) { - fprintf(stderr, "Cuda error: %s: %s.\n", msg, cudaGetErrorString( err) ); - exit(EXIT_FAILURE); + fprintf(stderr, "Cuda error: %s: %s.\n", msg, cudaGetErrorString( err) ); + exit(EXIT_FAILURE); } -} +} //LOOK: This function demonstrates how to use thrust for random number generation on the GPU! //Function that generates static. @@ -35,41 +44,40 @@ __host__ __device__ glm::vec3 generateRandomNumberFromThread(glm::vec2 resolutio return glm::vec3((float) u01(rng), (float) u01(rng), (float) u01(rng)); } -//Kernel that does the initial raycast from the camera. -__host__ __device__ ray raycastFromCameraKernel(glm::vec2 resolution, float time, int x, int y, glm::vec3 eye, glm::vec3 view, glm::vec3 up, glm::vec2 fov){ - - int index = x + (y * resolution.x); - - thrust::default_random_engine rng(hash(index*time)); - thrust::uniform_real_distribution u01(0,1); - - //standard camera raycast stuff - glm::vec3 E = eye; - glm::vec3 C = view; - glm::vec3 U = up; - float fovx = fov.x; - float fovy = fov.y; - - float CD = glm::length(C); - - glm::vec3 A = glm::cross(C, U); - glm::vec3 B = glm::cross(A, C); - glm::vec3 M = E+C; - glm::vec3 H = (A*float(CD*tan(fovx*(PI/180))))/float(glm::length(A)); - glm::vec3 V = (B*float(CD*tan(-fovy*(PI/180))))/float(glm::length(B)); - - float sx = (x)/(resolution.x-1); - float sy = (y)/(resolution.y-1); - - glm::vec3 P = M + (((2*sx)-1)*H) + (((2*sy)-1)*V); - glm::vec3 PmE = P-E; - glm::vec3 R = E + (float(200)*(PmE))/float(glm::length(PmE)); - - glm::vec3 direction = glm::normalize(R); - //major performance cliff at this point, TODO: find out why! +// Ray cast from the camera with anti-aliasing from stochastic sampling with changes of the x, y data type from integer to float +__host__ __device__ ray raycastFromCameraKernel(glm::vec2 resolution, float time, float x, float y, glm::vec3 eye, glm::vec3 view, glm::vec3 up, glm::vec2 fov, bool DOFflag){ + // Focal length in pixels + float focal = resolution.y / 2.0f / tan(fov.y * (PI / 180)); + + view = glm::normalize(view); + up = glm::normalize(up); + glm::vec3 right = glm::cross(view, up); + + // 3-D raycast vector from the camera in pixels + glm::vec3 rayCast = focal * view + (x - resolution.x / 2.0f) * right - (y - resolution.y / 2.0f) * up; + + // Output the data ray r; r.origin = eye; - r.direction = direction; + r.direction = glm::normalize(rayCast); + + // Depth of field referring to http://http.developer.nvidia.com/GPUGems/gpugems_ch23.html + if (!DOFflag) + return r; + + // Generate random angles + thrust::default_random_engine rng(hash(time)); + thrust::uniform_real_distribution u01(0, PI); + thrust::uniform_real_distribution u02(-1, 1); + + float angle = (float)u01(rng); + // Radius of the circle of confusion (Coc) + float radius = 0.5 * (float)u02(rng); + + glm::vec3 focalPoint = eye + r.direction * focal; + eye = eye + radius * cos(angle) * up + radius * sin(angle) * right; + r.origin = eye; + r.direction = glm::normalize(focalPoint - eye); return r; } @@ -83,7 +91,7 @@ __global__ void clearImage(glm::vec2 resolution, glm::vec3* image){ } } -//Kernel that writes the image to the OpenGL PBO directly. +//Kernel that writes the image to the OpenGL PBO directly. __global__ void sendImageToPBO(uchar4* PBOpos, glm::vec2 resolution, glm::vec3* image){ int x = (blockIdx.x * blockDim.x) + threadIdx.x; @@ -92,7 +100,7 @@ __global__ void sendImageToPBO(uchar4* PBOpos, glm::vec2 resolution, glm::vec3* if(x<=resolution.x && y<=resolution.y){ - glm::vec3 color; + glm::vec3 color; color.x = image[index].x*255.0; color.y = image[index].y*255.0; color.z = image[index].z*255.0; @@ -111,58 +119,159 @@ __global__ void sendImageToPBO(uchar4* PBOpos, glm::vec2 resolution, glm::vec3* // Each thread writes one pixel location in the texture (textel) PBOpos[index].w = 0; - PBOpos[index].x = color.x; + PBOpos[index].x = color.x; PBOpos[index].y = color.y; PBOpos[index].z = color.z; } } +__host__ __device__ int findPrimitiveObject(ray cameraRay, staticGeom* geoms, int numberOfGeoms, glm::vec3& normal, glm::vec3& intersectionPoint, float& distance) { + // Initialize the object index, distance and the intersection point and the normal + int objectIndex = -1; + float updateDistance; + distance = 1e7f; + glm::vec3 updateIntersectionPoint; + glm::vec3 updateNormal; + + // Find the nearest object to the camera + for(int i = 0; i < numberOfGeoms; ++ i) { + glm::vec3 updateIntersectionPoint, updateNormal; + if(geoms[i].type == SPHERE){ + updateDistance = sphereIntersectionTest(geoms[i], cameraRay, updateIntersectionPoint, updateNormal); + if(updateDistance > EPSILON && updateDistance < distance){ + distance = updateDistance; + intersectionPoint = updateIntersectionPoint; + normal = updateNormal; + objectIndex = i; + } + } else if(geoms[i].type == CUBE) { + updateDistance = boxIntersectionTest(geoms[i], cameraRay, updateIntersectionPoint, updateNormal); + if(updateDistance > EPSILON && updateDistance < distance){ + distance = updateDistance; + intersectionPoint = updateIntersectionPoint; + normal = updateNormal; + objectIndex = i; + } + } // if type + } // for i + return objectIndex; +} + +__host__ __device__ void rayPoolUpdate(material materials, glm::vec3 intersectionPoint, glm::vec3 normal, rayPool &rays, float time, cameraData cam, glm::vec3 &colors) { + // Determine if the object itself can be treated as light source + if(materials.emittance > 0) { + // If the object is a light source, then the color in the image is it own + colors = rays.colors * materials.color * materials.emittance; + } else { + ray currentRay; + currentRay.origin = rays.ray.origin; + currentRay.direction = rays.ray.direction; + glm::vec3 currentColor = rays.colors; + calculateBSDF(currentRay, rays, intersectionPoint, normal, currentColor, rays.colors, materials, time, cam); + colors = rays.colors * materials.emittance; + rays.isTerminated = false; + rays.ray.origin = intersectionPoint; + } // if material +} + + //TODO: IMPLEMENT THIS FUNCTION //Core raytracer kernel -__global__ void raytraceRay(glm::vec2 resolution, float time, float bounce, cameraData cam, int rayDepth, glm::vec3* colors, - staticGeom* geoms, int numberOfGeoms, material* materials, int numberOfMaterials){ +__global__ void raytraceRay(glm::vec2 resolution, float time, material* materials, int numberOfMaterials, cameraData cam, int rayDepth, glm::vec3* colors, + staticGeom* geoms, int numberOfGeoms, rayPool* rays){ + + int x = (blockIdx.x * blockDim.x) + threadIdx.x; + int y = (blockIdx.y * blockDim.y) + threadIdx.y; + int index = x + (y * resolution.x); + + // Find the primitive object + glm::vec3 normal, intersectionPoint; + float distance; + int objectIndex = findPrimitiveObject(rays[index].ray, geoms, numberOfGeoms, normal, intersectionPoint, distance); + // If there is no object, return + if(objectIndex == -1) + return; + + // Initialize the material index + int materialIndex = geoms[objectIndex].materialid; + // Check whether the ray is terminated, if not then update the ray pool for the next bounce + if (!rays[index].isTerminated) { + rays[index].isTerminated = true; + if (distance < EPSILON) + return; + rayPoolUpdate(materials[materialIndex], intersectionPoint, normal, rays[index], rayDepth * rays[index].index * time, cam, colors[rays[index].index]); + } +} // end of function + +// Initialize the ray pool +__global__ void rayGenerator (glm::vec2 resolution, float time, material* materials, int numberOfMaterials, cameraData cam, int rayDepth, glm::vec3* colors, + staticGeom* geoms, int numberOfGeoms, rayPool* rays) { int x = (blockIdx.x * blockDim.x) + threadIdx.x; int y = (blockIdx.y * blockDim.y) + threadIdx.y; int index = x + (y * resolution.x); - ray r = raycastFromCameraKernel(resolution, time, x, y, cam.position, cam.view, cam.up, cam.fov); - - if((x<=resolution.x && y<=resolution.y)){ - - float MAX_DEPTH = 100000000000000000; - float depth = MAX_DEPTH; - - for(int i=0; i-EPSILON){ - MAX_DEPTH = depth; - colors[index] = materials[geoms[i].materialid].color; - } - } + // Generate random numbers so as to stochastic sampling each pixel referring http://pages.cpsc.ucalgary.ca/~mario/courses/591-691-W06/PR/3-ray-tracing/3-advanced/readings/Cook_Stochastic_Sampling_TOG86.pdf + thrust::default_random_engine rng(hash(index * time)); + thrust::uniform_real_distribution u01(-0.5,0.5); + thrust::uniform_real_distribution u02(-0.5,0.5); + float xSample = x + u01(rng); + float ySample = y + u02(rng); + // Initialize the first ray from the camera + bool DOFflag = false; + ray cameraRay = raycastFromCameraKernel(resolution, index * time, xSample, ySample, cam.position, cam.view, cam.up, cam.fov, DOFflag); + + // Initialize the pixel color, ray and flags + colors[index] = glm::vec3(1); + rays[index].colors = glm::vec3(1); + rays[index].isTerminated = true; + rays[index].ray.origin = glm::vec3(0); + rays[index].ray.direction = glm::vec3(0); + if (x > resolution.x || y > resolution.y) + return; + + // Record the origin, direction and index for the next bounce + glm::vec3 normal, intersectionPoint; + float distance; + int objectIndex = findPrimitiveObject(cameraRay, geoms, numberOfGeoms, normal, intersectionPoint, distance); + if (distance < EPSILON || objectIndex == -1) + return; + + // Update the ray pool with the first ray from the camera + int materialIndex = geoms[objectIndex].materialid; + rays[index].ray.origin = cameraRay.origin; + rays[index].ray.direction = cameraRay.direction; + rayPoolUpdate(materials[materialIndex], intersectionPoint, normal, rays[index], index*time, cam, colors[index]); + rays[index].index = index; +} - //colors[index] = generateRandomNumberFromThread(resolution, time, x, y); - } +// Filter the noise between iteration and frame +__global__ void frameColorFilter (glm::vec2 resolution, float time, glm::vec3* currentColor, glm::vec3* updateColor) { + int x = (blockIdx.x * blockDim.x) + threadIdx.x; + int y = (blockIdx.y * blockDim.y) + threadIdx.y; + int index = x + (y * resolution.x); + if (x > resolution.x || y > resolution.y) + return; + updateColor[index] = currentColor[index] * (time - 1) / time + updateColor[index] / time; + glm::clamp(updateColor[index], 0, 1); } +// Check whether the bounce is terminated or not +struct IsRayTerminated { + __host__ __device__ bool operator() (const rayPool& rays) { + return rays.isTerminated; + } +}; //TODO: FINISH THIS FUNCTION // Wrapper for the __global__ call that sets up the kernel calls and does a ton of memory management void cudaRaytraceCore(uchar4* PBOpos, camera* renderCam, int frame, int iterations, material* materials, int numberOfMaterials, geom* geoms, int numberOfGeoms){ - int traceDepth = 1; //determines how many bounces the raytracer traces + int traceDepth = 1; //determines how many bounces the raytracer traces + float MAX_DEPTH = 8; // Maximum bounce number + bool motionBlurFlag = true; // set up crucial magic int tileSize = 8; @@ -173,7 +282,7 @@ void cudaRaytraceCore(uchar4* PBOpos, camera* renderCam, int frame, int iteratio glm::vec3* cudaimage = NULL; cudaMalloc((void**)&cudaimage, (int)renderCam->resolution.x*(int)renderCam->resolution.y*sizeof(glm::vec3)); cudaMemcpy( cudaimage, renderCam->image, (int)renderCam->resolution.x*(int)renderCam->resolution.y*sizeof(glm::vec3), cudaMemcpyHostToDevice); - + //package geometry and materials and sent to GPU staticGeom* geomList = new staticGeom[numberOfGeoms]; for(int i=0; iresolution.x*(int)renderCam->resolution.y*sizeof(glm::vec3)); + cudaMemcpy( cudacolor, renderCam->image, (int)renderCam->resolution.x*(int)renderCam->resolution.y*sizeof(glm::vec3), cudaMemcpyHostToDevice); + + // Package rays + rayPool* cudarays = NULL; + cudaMalloc((void**)&cudarays, (int)renderCam->resolution.x*(int)renderCam->resolution.y*sizeof(rayPool)); + //package camera cameraData cam; cam.resolution = renderCam->resolution; @@ -204,11 +337,29 @@ void cudaRaytraceCore(uchar4* PBOpos, camera* renderCam, int frame, int iteratio cam.up = renderCam->ups[frame]; cam.fov = renderCam->fov; - //kernel launches - for(int bounce = 1; bounce <= 1; ++bounce) - { - raytraceRay<<>>(renderCam->resolution, (float)iterations, (float)bounce, cam, traceDepth, cudaimage, cudageoms, numberOfGeoms, cudamaterials, numberOfMaterials); + // Intialize the rays + rayGenerator<<>>(renderCam->resolution, (float)iterations, cudamaterials, numberOfMaterials, cam, traceDepth, cudaimage, cudageoms, numberOfGeoms, cudarays); + checkCUDAError("Kernel failed!"); + + ++ traceDepth; // the ray has been bounced once in the ray generator + dim3 lessBlocksPerGrid; + int rayNumber = (int)renderCam->resolution.x*(int)renderCam->resolution.y; + while (traceDepth <= MAX_DEPTH) { + // Stream Compaction check whether the ray is terminated or not, if it is then move it from the ray pool and decrease the block number per grid for fast computation usage of thrust http://stackoverflow.com/questions/12201446/converting-thrustiterators-to-and-from-raw-pointers/12236270#12236270 + thrust::device_ptr deviceRaysPointer(cudarays); + thrust::device_ptr deviceRaysPointerEnd = thrust::remove_if(deviceRaysPointer, deviceRaysPointer + rayNumber, IsRayTerminated()); + rayNumber = deviceRaysPointerEnd.get() - deviceRaysPointer.get(); + if (rayNumber < EPSILON) + break; + lessBlocksPerGrid = dim3((int)ceil(float(renderCam->resolution.x)/float(tileSize)), (int)ceil(float(rayNumber)/(float)renderCam->resolution.x/float(tileSize))); + //kernel launches for different depth of ray + raytraceRay<<>>(renderCam->resolution, (float)iterations, cudamaterials, numberOfMaterials, cam, traceDepth, cudaimage, cudageoms, numberOfGeoms, cudarays); + ++ traceDepth; } + + // Filter between frames + frameColorFilter<<>>(renderCam->resolution, (float)iterations, cudacolor, cudaimage); + sendImageToPBO<<>>(PBOpos, renderCam->resolution, cudaimage); //retrieve image from GPU @@ -218,10 +369,12 @@ void cudaRaytraceCore(uchar4* PBOpos, camera* renderCam, int frame, int iteratio cudaFree( cudaimage ); cudaFree( cudageoms ); cudaFree( cudamaterials ); - delete [] geomList; + cudaFree( cudacolor ); + cudaFree( cudarays ); + delete geomList; - // make certain the kernel has completed + // make certain the kernel has completed cudaThreadSynchronize(); checkCUDAError("Kernel failed!"); -} +}; \ No newline at end of file diff --git a/src/sceneStructs.h b/src/sceneStructs.h index b10f1cf..fbfd64a 100755 --- a/src/sceneStructs.h +++ b/src/sceneStructs.h @@ -18,6 +18,14 @@ struct ray { glm::vec3 direction; }; +struct rayPool { + ray ray; + glm::vec3 colors; + float coefficient; + int index; + bool isTerminated; +}; + struct geom { enum GEOMTYPE type; int materialid; diff --git a/src/utilities.h b/src/utilities.h index 15b6495..84ec55f 100755 --- a/src/utilities.h +++ b/src/utilities.h @@ -17,13 +17,13 @@ #include #include "cudaMat4.h" -const float PI =3.1415926535897932384626422832795028841971; -const float TWO_PI =6.2831853071795864769252867665590057683943; -const float SQRT_OF_ONE_THIRD =0.5773502691896257645091487805019574556476; -const float E =2.7182818284590452353602874713526624977572; -const float EPSILON =.000000001; -const float ZERO_ABSORPTION_EPSILON =0.00001; -const float RAY_BIAS_AMOUNT =0.0002; +#define PI 3.1415926535897932384626422832795028841971 +#define TWO_PI 6.2831853071795864769252867665590057683943 +#define SQRT_OF_ONE_THIRD 0.5773502691896257645091487805019574556476 +#define NATURAL_E 2.7182818284590452353602874713526624977572 +#define EPSILON .000000001 +#define ZERO_ABSORPTION_EPSILON 0.00001 +#define RAY_BIAS_AMOUNT 0.0002 namespace utilityCore { extern float clamp(float f, float min, float max); diff --git a/videoScreenShot.PNG b/videoScreenShot.PNG new file mode 100755 index 0000000..9681667 Binary files /dev/null and b/videoScreenShot.PNG differ