Scratchapixel 2.0
Sign in
This project contains the following files (right-click files you'd like to download):
whitted.cpp
//[header] // A simple program to demonstrate how to implement Whitted-style ray-tracing //[/header] //[compile] // Download the whitted.cpp file to a folder. // Open a shell/terminal, and run the following command where the files is saved: // // c++ -o whitted whitted.cpp -std=c++11 -O3 // // Run with: ./whitted. Open the file ./out.png in Photoshop or any program // reading PPM files. //[/compile] //[ignore] // Copyright (C) 2012 www.scratchapixel.com // // This program is free software: you can redistribute it and/or modify // it under the terms of the GNU General Public License as published by // the Free Software Foundation, either version 3 of the License, or // (at your option) any later version. // // This program is distributed in the hope that it will be useful, // but WITHOUT ANY WARRANTY; without even the implied warranty of // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the // GNU General Public License for more details. // // You should have received a copy of the GNU General Public License // along with this program. If not, see <http://www.gnu.org/licenses/>. //[/ignore] #include <cstdio> #include <cstdlib> #include <memory> #include <vector> #include <utility> #include <cstdint> #include <iostream> #include <fstream> #include <cmath> const float kInfinity = std::numeric_limits<float>::max(); class Vec3f { public: Vec3f() : x(0), y(0), z(0) {} Vec3f(float xx) : x(xx), y(xx), z(xx) {} Vec3f(float xx, float yy, float zz) : x(xx), y(yy), z(zz) {} Vec3f operator * (const float &r) const { return Vec3f(x * r, y * r, z * r); } Vec3f operator * (const Vec3f &v) const { return Vec3f(x * v.x, y * v.y, z * v.z); } Vec3f operator - (const Vec3f &v) const { return Vec3f(x - v.x, y - v.y, z - v.z); } Vec3f operator + (const Vec3f &v) const { return Vec3f(x + v.x, y + v.y, z + v.z); } Vec3f operator - () const { return Vec3f(-x, -y, -z); } Vec3f& operator += (const Vec3f &v) { x += v.x, y += v.y, z += v.z; return *this; } friend Vec3f operator * (const float &r, const Vec3f &v) { return Vec3f(v.x * r, v.y * r, v.z * r); } friend std::ostream & operator << (std::ostream &os, const Vec3f &v) { return os << v.x << ", " << v.y << ", " << v.z; } float x, y, z; }; class Vec2f { public: Vec2f() : x(0), y(0) {} Vec2f(float xx) : x(xx), y(xx) {} Vec2f(float xx, float yy) : x(xx), y(yy) {} Vec2f operator * (const float &r) const { return Vec2f(x * r, y * r); } Vec2f operator + (const Vec2f &v) const { return Vec2f(x + v.x, y + v.y); } float x, y; }; Vec3f normalize(const Vec3f &v) { float mag2 = v.x * v.x + v.y * v.y + v.z * v.z; if (mag2 > 0) { float invMag = 1 / sqrtf(mag2); return Vec3f(v.x * invMag, v.y * invMag, v.z * invMag); } return v; } inline float dotProduct(const Vec3f &a, const Vec3f &b) { return a.x * b.x + a.y * b.y + a.z * b.z; } Vec3f crossProduct(const Vec3f &a, const Vec3f &b) { return Vec3f( a.y * b.z - a.z * b.y, a.z * b.x - a.x * b.z, a.x * b.y - a.y * b.x ); } inline float clamp(const float &lo, const float &hi, const float &v) { return std::max(lo, std::min(hi, v)); } inline float deg2rad(const float &deg) { return deg * M_PI / 180; } inline Vec3f mix(const Vec3f &a, const Vec3f& b, const float &mixValue) { return a * (1 - mixValue) + b * mixValue; } struct Options { uint32_t width; uint32_t height; float fov; float imageAspectRatio; uint8_t maxDepth; Vec3f backgroundColor; float bias; }; class Light { public: Light(const Vec3f &p, const Vec3f &i) : position(p), intensity(i) {} Vec3f position; Vec3f intensity; }; enum MaterialType { DIFFUSE_AND_GLOSSY, REFLECTION_AND_REFRACTION, REFLECTION }; class Object { public: Object() : materialType(DIFFUSE_AND_GLOSSY), ior(1.3), Kd(0.8), Ks(0.2), diffuseColor(0.2), specularExponent(25) {} virtual ~Object() {} virtual bool intersect(const Vec3f &, const Vec3f &, float &, uint32_t &, Vec2f &) const = 0; virtual void getSurfaceProperties(const Vec3f &, const Vec3f &, const uint32_t &, const Vec2f &, Vec3f &, Vec2f &) const = 0; virtual Vec3f evalDiffuseColor(const Vec2f &) const { return diffuseColor; } // material properties MaterialType materialType; float ior; float Kd, Ks; Vec3f diffuseColor; float specularExponent; }; bool solveQuadratic(const float &a, const float &b, const float &c, float &x0, float &x1) { float discr = b * b - 4 * a * c; if (discr < 0) return false; else if (discr == 0) x0 = x1 = - 0.5 * b / a; else { float q = (b > 0) ? -0.5 * (b + sqrt(discr)) : -0.5 * (b - sqrt(discr)); x0 = q / a; x1 = c / q; } if (x0 > x1) std::swap(x0, x1); return true; } class Sphere : public Object{ public: Sphere(const Vec3f &c, const float &r) : center(c), radius(r), radius2(r * r) {} bool intersect(const Vec3f &orig, const Vec3f &dir, float &tnear, uint32_t &index, Vec2f &uv) const { // analytic solution Vec3f L = orig - center; float a = dotProduct(dir, dir); float b = 2 * dotProduct(dir, L); float c = dotProduct(L, L) - radius2; float t0, t1; if (!solveQuadratic(a, b, c, t0, t1)) return false; if (t0 < 0) t0 = t1; if (t0 < 0) return false; tnear = t0; return true; } void getSurfaceProperties(const Vec3f &P, const Vec3f &I, const uint32_t &index, const Vec2f &uv, Vec3f &N, Vec2f &st) const { N = normalize(P - center); } Vec3f center; float radius, radius2; }; bool rayTriangleIntersect( const Vec3f &v0, const Vec3f &v1, const Vec3f &v2, const Vec3f &orig, const Vec3f &dir, float &tnear, float &u, float &v) { Vec3f edge1 = v1 - v0; Vec3f edge2 = v2 - v0; Vec3f pvec = crossProduct(dir, edge2); float det = dotProduct(edge1, pvec); if (det == 0 || det < 0) return false; Vec3f tvec = orig - v0; u = dotProduct(tvec, pvec); if (u < 0 || u > det) return false; Vec3f qvec = crossProduct(tvec, edge1); v = dotProduct(dir, qvec); if (v < 0 || u + v > det) return false; float invDet = 1 / det; tnear = dotProduct(edge2, qvec) * invDet; u *= invDet; v *= invDet; return true; } class MeshTriangle : public Object { public: MeshTriangle( const Vec3f *verts, const uint32_t *vertsIndex, const uint32_t &numTris, const Vec2f *st) { uint32_t maxIndex = 0; for (uint32_t i = 0; i < numTris * 3; ++i) if (vertsIndex[i] > maxIndex) maxIndex = vertsIndex[i]; maxIndex += 1; vertices = std::unique_ptr<Vec3f[]>(new Vec3f[maxIndex]); memcpy(vertices.get(), verts, sizeof(Vec3f) * maxIndex); vertexIndex = std::unique_ptr<uint32_t[]>(new uint32_t[numTris * 3]); memcpy(vertexIndex.get(), vertsIndex, sizeof(uint32_t) * numTris * 3); numTriangles = numTris; stCoordinates = std::unique_ptr<Vec2f[]>(new Vec2f[maxIndex]); memcpy(stCoordinates.get(), st, sizeof(Vec2f) * maxIndex); } bool intersect(const Vec3f &orig, const Vec3f &dir, float &tnear, uint32_t &index, Vec2f &uv) const { bool intersect = false; for (uint32_t k = 0; k < numTriangles; ++k) { const Vec3f & v0 = vertices[vertexIndex[k * 3]]; const Vec3f & v1 = vertices[vertexIndex[k * 3 + 1]]; const Vec3f & v2 = vertices[vertexIndex[k * 3 + 2]]; float t, u, v; if (rayTriangleIntersect(v0, v1, v2, orig, dir, t, u, v) && t < tnear) { tnear = t; uv.x = u; uv.y = v; index = k; intersect |= true; } } return intersect; } void getSurfaceProperties(const Vec3f &P, const Vec3f &I, const uint32_t &index, const Vec2f &uv, Vec3f &N, Vec2f &st) const { const Vec3f &v0 = vertices[vertexIndex[index * 3]]; const Vec3f &v1 = vertices[vertexIndex[index * 3 + 1]]; const Vec3f &v2 = vertices[vertexIndex[index * 3 + 2]]; Vec3f e0 = normalize(v1 - v0); Vec3f e1 = normalize(v2 - v1); N = normalize(crossProduct(e0, e1)); const Vec2f &st0 = stCoordinates[vertexIndex[index * 3]]; const Vec2f &st1 = stCoordinates[vertexIndex[index * 3 + 1]]; const Vec2f &st2 = stCoordinates[vertexIndex[index * 3 + 2]]; st = st0 * (1 - uv.x - uv.y) + st1 * uv.x + st2 * uv.y; } Vec3f evalDiffuseColor(const Vec2f &st) const { float scale = 5; float pattern = (fmodf(st.x * scale, 1) > 0.5) ^ (fmodf(st.y * scale, 1) > 0.5); return mix(Vec3f(0.815, 0.235, 0.031), Vec3f(0.937, 0.937, 0.231), pattern); } std::unique_ptr<Vec3f[]> vertices; uint32_t numTriangles; std::unique_ptr<uint32_t[]> vertexIndex; std::unique_ptr<Vec2f[]> stCoordinates; }; // [comment] // Compute reflection direction // [/comment] Vec3f reflect(const Vec3f &I, const Vec3f &N) { return I - 2 * dotProduct(I, N) * N; } // [comment] // Compute refraction direction using Snell's law // // We need to handle with care the two possible situations: // // - When the ray is inside the object // // - When the ray is outside. // // If the ray is outside, you need to make cosi positive cosi = -N.I // // If the ray is inside, you need to invert the refractive indices and negate the normal N // [/comment] Vec3f refract(const Vec3f &I, const Vec3f &N, const float &ior) { float cosi = clamp(-1, 1, dotProduct(I, N)); float etai = 1, etat = ior; Vec3f n = N; if (cosi < 0) { cosi = -cosi; } else { std::swap(etai, etat); n= -N; } float eta = etai / etat; float k = 1 - eta * eta * (1 - cosi * cosi); return k < 0 ? 0 : eta * I + (eta * cosi - sqrtf(k)) * n; } // [comment] // Compute Fresnel equation // // \param I is the incident view direction // // \param N is the normal at the intersection point // // \param ior is the mateural refractive index // // \param[out] kr is the amount of light reflected // [/comment] void fresnel(const Vec3f &I, const Vec3f &N, const float &ior, float &kr) { float cosi = clamp(-1, 1, dotProduct(I, N)); float etai = 1, etat = ior; if (cosi > 0) { std::swap(etai, etat); } // Compute sini using Snell's law float sint = etai / etat * sqrtf(std::max(0.f, 1 - cosi * cosi)); // Total internal reflection if (sint >= 1) { kr = 1; } else { float cost = sqrtf(std::max(0.f, 1 - sint * sint)); cosi = fabsf(cosi); float Rs = ((etat * cosi) - (etai * cost)) / ((etat * cosi) + (etai * cost)); float Rp = ((etai * cosi) - (etat * cost)) / ((etai * cosi) + (etat * cost)); kr = (Rs * Rs + Rp * Rp) / 2; } // As a consequence of the conservation of energy, transmittance is given by: // kt = 1 - kr; } // [comment] // Returns true if the ray intersects an object, false otherwise. // // \param orig is the ray origin // // \param dir is the ray direction // // \param objects is the list of objects the scene contains // // \param[out] tNear contains the distance to the cloesest intersected object. // // \param[out] index stores the index of the intersect triangle if the interesected object is a mesh. // // \param[out] uv stores the u and v barycentric coordinates of the intersected point // // \param[out] *hitObject stores the pointer to the intersected object (used to retrieve material information, etc.) // // \param isShadowRay is it a shadow ray. We can return from the function sooner as soon as we have found a hit. // [/comment] bool trace( const Vec3f &orig, const Vec3f &dir, const std::vector<std::unique_ptr<Object>> &objects, float &tNear, uint32_t &index, Vec2f &uv, Object **hitObject) { *hitObject = nullptr; for (uint32_t k = 0; k < objects.size(); ++k) { float tNearK = kInfinity; uint32_t indexK; Vec2f uvK; if (objects[k]->intersect(orig, dir, tNearK, indexK, uvK) && tNearK < tNear) { *hitObject = objects[k].get(); tNear = tNearK; index = indexK; uv = uvK; } } return (*hitObject != nullptr); } // [comment] // Implementation of the Whitted-syle light transport algorithm (E [S*] (D|G) L) // // This function is the function that compute the color at the intersection point // of a ray defined by a position and a direction. Note that thus function is recursive (it calls itself). // // If the material of the intersected object is either reflective or reflective and refractive, // then we compute the reflection/refracton direction and cast two new rays into the scene // by calling the castRay() function recursively. When the surface is transparent, we mix // the reflection and refraction color using the result of the fresnel equations (it computes // the amount of reflection and refractin depending on the surface normal, incident view direction // and surface refractive index). // // If the surface is duffuse/glossy we use the Phong illumation model to compute the color // at the intersection point. // [/comment] Vec3f castRay( const Vec3f &orig, const Vec3f &dir, const std::vector<std::unique_ptr<Object>> &objects, const std::vector<std::unique_ptr<Light>> &lights, const Options &options, uint32_t depth, bool test = false) { if (depth > options.maxDepth) { return options.backgroundColor; } Vec3f hitColor = options.backgroundColor; float tnear = kInfinity; Vec2f uv; uint32_t index = 0; Object *hitObject = nullptr; if (trace(orig, dir, objects, tnear, index, uv, &hitObject)) { Vec3f hitPoint = orig + dir * tnear; Vec3f N; // normal Vec2f st; // st coordinates hitObject->getSurfaceProperties(hitPoint, dir, index, uv, N, st); Vec3f tmp = hitPoint; switch (hitObject->materialType) { case REFLECTION_AND_REFRACTION: { Vec3f reflectionDirection = normalize(reflect(dir, N)); Vec3f refractionDirection = normalize(refract(dir, N, hitObject->ior)); Vec3f reflectionRayOrig = (dotProduct(reflectionDirection, N) < 0) ? hitPoint - N * options.bias : hitPoint + N * options.bias; Vec3f refractionRayOrig = (dotProduct(refractionDirection, N) < 0) ? hitPoint - N * options.bias : hitPoint + N * options.bias; Vec3f reflectionColor = castRay(reflectionRayOrig, reflectionDirection, objects, lights, options, depth + 1, 1); Vec3f refractionColor = castRay(refractionRayOrig, refractionDirection, objects, lights, options, depth + 1, 1); float kr; fresnel(dir, N, hitObject->ior, kr); hitColor = reflectionColor * kr + refractionColor * (1 - kr); break; } case REFLECTION: { float kr; fresnel(dir, N, hitObject->ior, kr); Vec3f reflectionDirection = reflect(dir, N); Vec3f reflectionRayOrig = (dotProduct(reflectionDirection, N) < 0) ? hitPoint + N * options.bias : hitPoint - N * options.bias; hitColor = castRay(reflectionRayOrig, reflectionDirection, objects, lights, options, depth + 1) * kr; break; } default: { // [comment] // We use the Phong illumation model int the default case. The phong model // is composed of a diffuse and a specular reflection component. // [/comment] Vec3f lightAmt = 0, specularColor = 0; Vec3f shadowPointOrig = (dotProduct(dir, N) < 0) ? hitPoint + N * options.bias : hitPoint - N * options.bias; // [comment] // Loop over all lights in the scene and sum their contribution up // We also apply the lambert cosine law here though we haven't explained yet what this means. // [/comment] for (uint32_t i = 0; i < lights.size(); ++i) { Vec3f lightDir = lights[i]->position - hitPoint; // square of the distance between hitPoint and the light float lightDistance2 = dotProduct(lightDir, lightDir); lightDir = normalize(lightDir); float LdotN = std::max(0.f, dotProduct(lightDir, N)); Object *shadowHitObject = nullptr; float tNearShadow = kInfinity; // is the point in shadow, and is the nearest occluding object closer to the object than the light itself? bool inShadow = trace(shadowPointOrig, lightDir, objects, tNearShadow, index, uv, &shadowHitObject) && tNearShadow * tNearShadow < lightDistance2; lightAmt += (1 - inShadow) * lights[i]->intensity * LdotN; Vec3f reflectionDirection = reflect(-lightDir, N); specularColor += powf(std::max(0.f, -dotProduct(reflectionDirection, dir)), hitObject->specularExponent) * lights[i]->intensity; } hitColor = lightAmt * hitObject->evalDiffuseColor(st) * hitObject->Kd + specularColor * hitObject->Ks; break; } } } return hitColor; } // [comment] // The main render function. This where we iterate over all pixels in the image, generate // primary rays and cast these rays into the scene. The content of the framebuffer is // saved to a file. // [/comment] void render( const Options &options, const std::vector<std::unique_ptr<Object>> &objects, const std::vector<std::unique_ptr<Light>> &lights) { Vec3f *framebuffer = new Vec3f[options.width * options.height]; Vec3f *pix = framebuffer; float scale = tan(deg2rad(options.fov * 0.5)); float imageAspectRatio = options.width / (float)options.height; Vec3f orig(0); for (uint32_t j = 0; j < options.height; ++j) { for (uint32_t i = 0; i < options.width; ++i) { // generate primary ray direction float x = (2 * (i + 0.5) / (float)options.width - 1) * imageAspectRatio * scale; float y = (1 - 2 * (j + 0.5) / (float)options.height) * scale; Vec3f dir = normalize(Vec3f(x, y, -1)); *(pix++) = castRay(orig, dir, objects, lights, options, 0); } } // save framebuffer to file std::ofstream ofs; ofs.open("./out.ppm"); ofs << "P6\n" << options.width << " " << options.height << "\n255\n"; for (uint32_t i = 0; i < options.height * options.width; ++i) { char r = (char)(255 * clamp(0, 1, framebuffer[i].x)); char g = (char)(255 * clamp(0, 1, framebuffer[i].y)); char b = (char)(255 * clamp(0, 1, framebuffer[i].z)); ofs << r << g << b; } ofs.close(); delete [] framebuffer; } // [comment] // In the main function of the program, we create the scene (create objects and lights) // as well as set the options for the render (image widht and height, maximum recursion // depth, field-of-view, etc.). We then call the render function(). // [/comment] int main(int argc, char **argv) { // creating the scene (adding objects and lights) std::vector<std::unique_ptr<Object>> objects; std::vector<std::unique_ptr<Light>> lights; Sphere *sph1 = new Sphere(Vec3f(-1, 0, -12), 2); sph1->materialType = DIFFUSE_AND_GLOSSY; sph1->diffuseColor = Vec3f(0.6, 0.7, 0.8); Sphere *sph2 = new Sphere(Vec3f(0.5, -0.5, -8), 1.5); sph2->ior = 1.5; sph2->materialType = REFLECTION_AND_REFRACTION; objects.push_back(std::unique_ptr<Sphere>(sph1)); objects.push_back(std::unique_ptr<Sphere>(sph2)); Vec3f verts[4] = {{-5,-3,-6}, {5,-3,-6}, {5,-3,-16}, {-5,-3,-16}}; uint32_t vertIndex[6] = {0, 1, 3, 1, 2, 3}; Vec2f st[4] = {{0, 0}, {1, 0}, {1, 1}, {0, 1}}; MeshTriangle *mesh = new MeshTriangle(verts, vertIndex, 2, st); mesh->materialType = DIFFUSE_AND_GLOSSY; objects.push_back(std::unique_ptr<MeshTriangle>(mesh)); lights.push_back(std::unique_ptr<Light>(new Light(Vec3f(-20, 70, 20), 0.5))); lights.push_back(std::unique_ptr<Light>(new Light(Vec3f(30, 50, -12), 1))); // setting up options Options options; options.width = 640; options.height = 480; options.fov = 90; options.backgroundColor = Vec3f(0.235294, 0.67451, 0.843137); options.maxDepth = 5; options.bias = 0.00001; // finally, render render(options, objects, lights); return 0; }