Scratchapixel 2.0
Sign in
This project contains the following files (right-click files you'd like to download):
distfields.cppgeometry.h
//[header] // This program generate and render the Utah teapot //[/header] //[compile] // Download the distfields.cpp and geometry.h files to a folder. // Open a shell/terminal, and run the following command where the files are saved: // // clang++ -o distfields distfields.cpp -std=c++14 -O3 // // Run with: ./distfileds. Open the file ./spheretrace.ppm 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 <cstdlib> #include <cstdio> #include <fstream> #include "geometry.h" #include <cmath> #include <vector> #include <limits> #include <memory> constexpr float kInfinity = std::numeric_limits<float>::max(); Matrix44f lookAt(const Vec3f& from, const Vec3f& to, Vec3f tmp = Vec3f(0, 1, 0)) { // because the forward axis in a right hand coordinate system points backward we compute -(to - from) Vec3f z = from - to; z.normalize(); Vec3f x = tmp.normalize().crossProduct(z); x.normalize(); // this is just in case Up is not normalized Vec3f y = z.crossProduct(x); Matrix44f camToWorld; // set x-axis camToWorld[0][0] = x[0]; camToWorld[0][1] = x[1]; camToWorld[0][2] = x[2]; // set y-axis camToWorld[1][0] = y[0]; camToWorld[1][1] = y[1]; camToWorld[1][2] = y[2]; // set z-axis camToWorld[2][0] = z[0]; camToWorld[2][1] = z[1]; camToWorld[2][2] = z[2]; // set position camToWorld[3][0] = from[0]; camToWorld[3][1] = from[1]; camToWorld[3][2] = from[2]; return camToWorld; } float degToRad(const float& angle) { return angle / 180.f * M_PI; } // [comment] // ImplicitShape base class // [/comment] class ImplicitShape { public: virtual float getDistance(const Vec3f& from) const = 0; virtual ~ImplicitShape() {} }; // [comment] // Implicit sphere surface // [/comment] class Sphere : public ImplicitShape { public: Sphere(const Vec3f& c, const float& r) : center(c), radius(r) {} float getDistance(const Vec3f& from) const { return (from - center).length() - radius; } float radius; Vec3f center; }; // [comment] // Implicit plane surface // [/comment] class Plane : public ImplicitShape { public: Plane(const Vec3f& nn = Vec3f(0, 1, 0), const Vec3f& pp = Vec3f(0)) : n(nn), pointOnPlane(pp) {} float getDistance(const Vec3f& from) const { return n.dotProduct(from - pointOnPlane); } Vec3f n, pointOnPlane; }; // [comment] // Implicit torus surface // [/comment] class Torus : public ImplicitShape { public: Torus(const float& _r0, const float& _r1) : r0(_r0), r1(_r1) {} float getDistance(const Vec3f& from) const { // reduce 3D point to 2D point float tmpx = sqrtf(from.x * from.x + from.z * from.z) - r0; float tmpy = from.y; // distance to cicle return sqrtf(tmpx * tmpx + tmpy * tmpy) - r1; } float r0, r1; }; // [comment] // Implicit cube surface // [/comment] class Cube : public ImplicitShape { public: Cube(const Vec3f &_corner) : corner(_corner) {} float getDistance(const Vec3f& from) const { #if 0 // first transform the input point into the object's "object-space". float scale = 2.f; // this matrix doesn't scale the object Matrix44f objectToWorld(0.542903, -0.545887, 0.638172, 0, 0.778733, 0.611711, -0.139228, 0, -0.314374, 0.572553, 0.7572, 0, 0, 1.459974, 0, 1); Matrix44f worldToObject = objectToWorld.inverse(); Vec3f fromObjectSpace = from; worldToObject.multVecMatrix(from, fromObjectSpace); #else Vec3f fromObjectSpace = from; float scale = 1; #endif fromObjectSpace *= 1.f / scale; fromObjectSpace.x = std::fabs(fromObjectSpace.x); fromObjectSpace.y = std::fabs(fromObjectSpace.y); fromObjectSpace.z = std::fabs(fromObjectSpace.z); // now compute the distance from the point to the neares point on the surface of the object Vec3f d = fromObjectSpace - corner; Vec3f dmax = d; dmax.x = std::max(dmax.x, 0.f); dmax.y = std::max(dmax.y, 0.f); dmax.z = std::max(dmax.z, 0.f); // don't forget to apply the scale back return scale * (std::min(std::max(d.x, std::max(d.y, d.z)), 0.f) + dmax.length()); } Vec3f corner{0.25, 0.25, 0.25}; }; struct unionFunc { float operator() (float a, float b) const { return std::min(a, b); } }; struct subtractFunc { float operator() (float a, float b) const { return std::max(-a, b); } }; struct intersectionFunc { float operator() (float a, float b) const { return std::max(a, b); } }; struct blendFunc { blendFunc(const float &_k) : k(_k) {} float operator() (float a, float b) const { float res = exp(-k * a) + exp(-k * b); return -log(std::max(0.0001f, res)) / k; } float k; }; struct mixFunc { mixFunc(const float &_t) : t(_t) {} float operator() (float a, float b) const { return a * (1 -t) + b * t; } float t; }; // [comment] // Combining implict shapes using CSG // [/comment] template<typename Op, typename ... Args> class CSG : public ImplicitShape { public: CSG( const std::shared_ptr<ImplicitShape> s1, const std::shared_ptr<ImplicitShape> s2, Args&& ... args) : op(std::forward<Args>(args) ...), shape1(s1), shape2(s2) {} float getDistance(const Vec3f& from) const { return op(shape1->getDistance(from), shape2->getDistance(from)); } Op op; const std::shared_ptr<ImplicitShape> shape1, shape2; }; // [comment] // Blobbies // [/comment] class SoftObject : public ImplicitShape { struct Blob { float R; // radius Vec3f c; // blob center }; public: SoftObject() { #if 1 blobbies.push_back({2.0, Vec3f(-1, 0, 0)}); blobbies.push_back({1.5, Vec3f( 1, 0, 0)}); #else for (size_t i = 0; i < 20; ++i) { float radius = (0.3 + drand48() * 1.3); Vec3f c((0.5 - drand48()) * 3, (0.5 - drand48()) * 3, (0.5 - drand48()) * 3); blobbies.push_back({radius, c}); } #endif } float getDistance(const Vec3f& from) const { float sumDensity = 0; float sumRi = 0; float minDistance = kInfinity; for (const auto& blob: blobbies) { float r = (blob.c - from).length(); if (r <= blob.R) { // this can be factored for speed if you want sumDensity += 2 * (r * r * r) / (blob.R * blob.R * blob.R) - 3 * (r * r) / (blob.R * blob.R) + 1; } minDistance = std::min(minDistance, r - blob.R); sumRi += blob.R; } return std::max(minDistance, (magic - sumDensity) / ( 3 / 2.f * sumRi)); } float magic{ 0.2 }; std::vector<Blob> blobbies; }; using Union = CSG<unionFunc>; using Subtract = CSG<subtractFunc>; using Intersect = CSG<intersectionFunc>; using Blend = CSG<blendFunc, float>; using Mix = CSG<mixFunc, float>; class PointLight { public: PointLight(const Vec3f& p, const Vec3f& c, const float& i) : pos(p), col(c), intensity(i) {} Vec3f pos; Vec3f col; float intensity; }; std::vector<std::shared_ptr<ImplicitShape>> makeScene() { std::vector<std::shared_ptr<ImplicitShape>> shapes; #if 0 shapes.push_back(std::make_shared<Plane>(Vec3f(0, 1, 0), Vec3f(0, -2, 0))); shapes.push_back(std::make_shared<Cube>()); shapes.push_back(std::make_shared<Torus>(2, 0.65)); #elif 0 shapes.push_back(std::make_shared<Plane>(Vec3f(0, 1, 0), Vec3f(0, -2, 0))); shapes.push_back(std::make_shared<Blend>( std::make_shared<Cube>(Vec3f(1.5)), std::make_shared<Torus>(2, 0.65), 5)); #elif 0 shapes.push_back(std::make_shared<Blend>( std::make_shared<Plane>(Vec3f(0, 1, 0), Vec3f(0, 0, 0)), std::make_shared<Torus>(2, 0.65), 5)); #elif 0 shapes.push_back(std::make_shared<Plane>(Vec3f(0, 1, 0), Vec3f(0, -2, 0))); shapes.push_back(std::make_shared<Mix>( std::make_shared<Cube>(Vec3f(1)), std::make_shared<Sphere>(Vec3f(0), 1), 0.5)); #else shapes.push_back(std::make_shared<Plane>(Vec3f(0, 1, 0), Vec3f(0, -2, 0))); shapes.push_back(std::make_shared<SoftObject>()); #endif return shapes; } bool sphereTraceShadow( const Vec3f& rayOrigin, const Vec3f& rayDirection, const float& maxDistance, const std::vector<std::shared_ptr<ImplicitShape>> scene) { constexpr float threshold = 10e-5; float t = 0; while (t < maxDistance) { float minDistance = kInfinity; Vec3f from = rayOrigin + t * rayDirection; for (auto shape : scene) { float d = shape->getDistance(from); if (d < minDistance) minDistance = d; // did we find an intersection? if (minDistance <= threshold * t) return true; } // no intersection, move along the ray by minDistance t += minDistance; } return false; } Vec3f shade( const Vec3f& rayOrigin, const Vec3f& rayDirection, const float& t, const ImplicitShape *shape, const std::vector<std::shared_ptr<ImplicitShape>> scene, const std::vector<std::unique_ptr<PointLight>> &lights) { constexpr float delta = 10e-5; Vec3f p = rayOrigin + t * rayDirection; Vec3f n = Vec3f( shape->getDistance(p + Vec3f(delta, 0, 0)) - shape->getDistance(p + Vec3f(-delta, 0, 0)), shape->getDistance(p + Vec3f(0, delta, 0)) - shape->getDistance(p + Vec3f(0, -delta, 0)), shape->getDistance(p + Vec3f(0, 0, delta)) - shape->getDistance(p + Vec3f(0, 0, -delta))); n.normalize(); Vec3f R = 0; // loop over all lights in the scene and add their contribution to P's brightness for (const auto& light: lights) { Vec3f lightDir = light->pos - p; if (lightDir.dotProduct(n) > 0) { float dist2 = lightDir.norm(); lightDir.normalize(); bool shadow = 1 - sphereTraceShadow(p, lightDir, sqrtf(dist2), scene); R += shadow * lightDir.dotProduct(n) * light->col * light->intensity / (4 * M_PI * dist2); } } return R; } Vec3f sphereTrace( const Vec3f& rayOrigin, const Vec3f& rayDirection, const std::vector<std::shared_ptr<ImplicitShape>>& scene, const std::vector<std::unique_ptr<PointLight>>& lights) { constexpr float maxDistance = 100; float t = 0; uint32_t numSteps = 0; const ImplicitShape *isectShape = nullptr; constexpr float threshold = 10e-6; while (t < maxDistance) { float minDistance = kInfinity; Vec3f from = rayOrigin + t * rayDirection; for (const auto& shape : scene) { float d = shape->getDistance(from); if (d < minDistance) { minDistance = d; isectShape = shape.get(); } } if (minDistance <= threshold * t) { return shade(rayOrigin, rayDirection, t, isectShape, scene, lights); } t += minDistance; numSteps++; } return 0; } int main(int argc, char **argv) { srand48(13); Matrix44f camToWorld = lookAt(Vec3f(0, 1, 9), 0); std::vector<std::shared_ptr<ImplicitShape>> scene = makeScene(); std::vector<std::unique_ptr<PointLight>> lights; lights.push_back(std::make_unique<PointLight>(Vec3f( 20, 30, 20), Vec3f(1.0, 0.9, 0.7), 4000)); lights.push_back(std::make_unique<PointLight>(Vec3f(-20, 30, -20), Vec3f(0.8, 0.9, 1.0), 4000)); lights.push_back(std::make_unique<PointLight>(Vec3f( -5, 10, 20), Vec3f(1.0, 1.0, 1.0), 3000)); constexpr uint32_t width = 640, height = 480; constexpr float ratio = width / static_cast<float>(height); constexpr float fov = 60; float angle = tan(degToRad(fov * 0.5)); Vec3f *buffer = new Vec3f[width * height]; Vec3f rayOrigin; camToWorld.multVecMatrix(Vec3f(0), rayOrigin); for (uint32_t j = 0; j < height; ++j) { for (uint32_t i = 0; i < width; ++i) { float x = (2 * i / static_cast<float>(width) - 1) * ratio * angle; float y = (1 - j / static_cast<float>(height) * 2) * angle; Vec3f rayDirection; camToWorld.multDirMatrix(Vec3f(x, y, -1).normalize(), rayDirection); ImplicitShape *tmp; Vec3f pixelColor = sphereTrace(rayOrigin, rayDirection, scene, lights); buffer[width * j + i] = pixelColor; } } std::ofstream ofs; ofs.open("./spheretrace.ppm"); ofs << "P6\n" << width << " " << height << "\n255\n"; for (uint32_t i = 0; i < width * height; ++i) { unsigned char r = static_cast<unsigned char>(std::min(1.0f, buffer[i][0]) * 255); unsigned char g = static_cast<unsigned char>(std::min(1.0f, buffer[i][1]) * 255); unsigned char b = static_cast<unsigned char>(std::min(1.0f, buffer[i][2]) * 255); ofs << r << g << b; } ofs.close(); delete [] buffer; return 0; }