Scratchapixel 2.0
Sign in
This project contains the following files (right-click files you'd like to download):
skycolor.cpp
//[header] // Simulating the color of the sky (Nishita model). // // See "Display of The Earth Taking into Account Atmospheric Scattering" for more information. //[/header] //[compile] // Download the acceleration.cpp and teapotdata.h file to a folder. // Open a shell/terminal, and run the following command where the files is saved: // // clang++ -std=c++11 -o skycolor skycolor.cpp -O3 // // You can use c++ if you don't use clang++ // // Run with: ./skycolor. Open the resulting image (ppm) in Photoshop or any program // reading PPM files. //[/compile] //[ignore] // Copyright (C) 2016 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] #if defined(WIN32) || defined(_WIN32) #include "stdafx.h" #endif #include <cassert> #include <iostream> #include <fstream> #include <algorithm> #include <cmath> #include <chrono> #include <random> #include <limits> #ifndef M_PI #define M_PI (3.14159265358979323846f) #endif const float kInfinity = std::numeric_limits<float>::max(); template<typename T> class Vec3 { public: Vec3() : x(0), y(0), z(0) {} Vec3(T xx) : x(xx), y(xx), z(xx) {} Vec3(T xx, T yy, T zz) : x(xx), y(yy), z(zz) {} Vec3 operator * (const T& r) const { return Vec3(x * r, y * r, z * r); } Vec3 operator * (const Vec3<T> &v) const { return Vec3(x * v.x, y * v.y, z * v.z); } Vec3 operator + (const Vec3& v) const { return Vec3(x + v.x, y + v.y, z + v.z); } Vec3 operator - (const Vec3& v) const { return Vec3(x - v.x, y - v.y, z - v.z); } template<typename U> Vec3 operator / (const Vec3<U>& v) const { return Vec3(x / v.x, y / v.y, z / v.z); } friend Vec3 operator / (const T r, const Vec3& v) { return Vec3(r / v.x, r / v.y, r / v.z); } const T& operator [] (size_t i) const { return (&x)[i]; } T& operator [] (size_t i) { return (&x)[i]; } T length2() const { return x * x + y * y + z * z; } T length() const { return std::sqrt(length2()); } Vec3& operator += (const Vec3<T>& v) { x += v.x, y += v.y, z += v.z; return *this; } Vec3& operator *= (const float& r) { x *= r, y *= r, z *= r; return *this; } friend Vec3 operator * (const float&r, const Vec3& v) { return Vec3(v.x * r, v.y * r, v.z * r); } friend std::ostream& operator << (std::ostream& os, const Vec3<T>& v) { os << v.x << " " << v.y << " " << v.z << std::endl; return os; } T x, y, z; }; template<typename T> void normalize(Vec3<T>& vec) { T len2 = vec.length2(); if (len2 > 0) { T invLen = 1 / std::sqrt(len2); vec.x *= invLen, vec.y *= invLen, vec.z *= invLen; } } template<typename T> T dot(const Vec3<T>& va, const Vec3<T>& vb) { return va.x * vb.x + va.y * vb.y + va.z * vb.z; } using Vec3f = Vec3<float>; // [comment] // The atmosphere class. Stores data about the planetory body (its radius), the atmosphere itself // (thickness) and things such as the Mie and Rayleigh coefficients, the sun direction, etc. // [/comment] class Atmosphere { public: Atmosphere( Vec3f sd = Vec3f(0, 1, 0), float er = 6360e3, float ar = 6420e3, float hr = 7994, float hm = 1200) : sunDirection(sd), earthRadius(er), atmosphereRadius(ar), Hr(hr), Hm(hm) {} Vec3f computeIncidentLight(const Vec3f& orig, const Vec3f& dir, float tmin, float tmax) const; Vec3f sunDirection; // The sun direction (normalized) float earthRadius; // In the paper this is usually Rg or Re (radius ground, eart) float atmosphereRadius; // In the paper this is usually R or Ra (radius atmosphere) float Hr; // Thickness of the atmosphere if density was uniform (Hr) float Hm; // Same as above but for Mie scattering (Hm) static const Vec3f betaR; static const Vec3f betaM; }; const Vec3f Atmosphere::betaR(3.8e-6f, 13.5e-6f, 33.1e-6f); const Vec3f Atmosphere::betaM(21e-6f); bool solveQuadratic(float a, float b, float c, float& x1, float& x2) { if (b == 0) { // Handle special case where the the two vector ray.dir and V are perpendicular // with V = ray.orig - sphere.centre if (a == 0) return false; x1 = 0; x2 = std::sqrtf(-c / a); return true; } float discr = b * b - 4 * a * c; if (discr < 0) return false; float q = (b < 0.f) ? -0.5f * (b - std::sqrtf(discr)) : -0.5f * (b + std::sqrtf(discr)); x1 = q / a; x2 = c / q; return true; } // [comment] // A simple routine to compute the intersection of a ray with a sphere // [/comment] bool raySphereIntersect(const Vec3f& orig, const Vec3f& dir, const float& radius, float& t0, float& t1) { // They ray dir is normalized so A = 1 float A = dir.x * dir.x + dir.y * dir.y + dir.z * dir.z; float B = 2 * (dir.x * orig.x + dir.y * orig.y + dir.z * orig.z); float C = orig.x * orig.x + orig.y * orig.y + orig.z * orig.z - radius * radius; if (!solveQuadratic(A, B, C, t0, t1)) return false; if (t0 > t1) std::swap(t0, t1); return true; } // [comment] // This is where all the magic happens. We first raymarch along the primary ray (from the camera origin // to the point where the ray exits the atmosphere or intersect with the planetory body). For each // sample along the primary ray, we then "cast" a light ray and raymarch along that ray as well. // We basically shoot a ray in the direction of the sun. // [/comment] Vec3f Atmosphere::computeIncidentLight(const Vec3f& orig, const Vec3f& dir, float tmin, float tmax) const { float t0, t1; if (!raySphereIntersect(orig, dir, atmosphereRadius, t0, t1) || t1 < 0) return 0; if (t0 > tmin && t0 > 0) tmin = t0; if (t1 < tmax) tmax = t1; uint32_t numSamples = 16; uint32_t numSamplesLight = 8; float segmentLength = (tmax - tmin) / numSamples; float tCurrent = tmin; Vec3f sumR(0), sumM(0); // mie and rayleigh contribution float opticalDepthR = 0, opticalDepthM = 0; float mu = dot(dir, sunDirection); // mu in the paper which is the cosine of the angle between the sun direction and the ray direction float phaseR = 3.f / (16.f * M_PI) * (1 + mu * mu); float g = 0.76f; float phaseM = 3.f / (8.f * M_PI) * ((1.f - g * g) * (1.f + mu * mu)) / ((2.f + g * g) * pow(1.f + g * g - 2.f * g * mu, 1.5f)); for (uint32_t i = 0; i < numSamples; ++i) { Vec3f samplePosition = orig + (tCurrent + segmentLength * 0.5f) * dir; float height = samplePosition.length() - earthRadius; // compute optical depth for light float hr = exp(-height / Hr) * segmentLength; float hm = exp(-height / Hm) * segmentLength; opticalDepthR += hr; opticalDepthM += hm; // light optical depth float t0Light, t1Light; raySphereIntersect(samplePosition, sunDirection, atmosphereRadius, t0Light, t1Light); float segmentLengthLight = t1Light / numSamplesLight, tCurrentLight = 0; float opticalDepthLightR = 0, opticalDepthLightM = 0; uint32_t j; for (j = 0; j < numSamplesLight; ++j) { Vec3f samplePositionLight = samplePosition + (tCurrentLight + segmentLengthLight * 0.5f) * sunDirection; float heightLight = samplePositionLight.length() - earthRadius; if (heightLight < 0) break; opticalDepthLightR += exp(-heightLight / Hr) * segmentLengthLight; opticalDepthLightM += exp(-heightLight / Hm) * segmentLengthLight; tCurrentLight += segmentLengthLight; } if (j == numSamplesLight) { Vec3f tau = betaR * (opticalDepthR + opticalDepthLightR) + betaM * 1.1f * (opticalDepthM + opticalDepthLightM); Vec3f attenuation(exp(-tau.x), exp(-tau.y), exp(-tau.z)); sumR += attenuation * hr; sumM += attenuation * hm; } tCurrent += segmentLength; } // [comment] // We use a magic number here for the intensity of the sun (20). We will make it more // scientific in a future revision of this lesson/code // [/comment] return (sumR * betaR * phaseR + sumM * betaM * phaseM) * 20; } void renderSkydome(const Vec3f& sunDir, const char *filename) { Atmosphere atmosphere(sunDir); auto t0 = std::chrono::high_resolution_clock::now(); #if 1 // [comment] // Render fisheye // [/comment] const unsigned width = 512, height = 512; Vec3f *image = new Vec3f[width * height], *p = image; memset(image, 0x0, sizeof(Vec3f) * width * height); for (unsigned j = 0; j < height; ++j) { float y = 2.f * (j + 0.5f) / float(height - 1) - 1.f; for (unsigned i = 0; i < width; ++i, ++p) { float x = 2.f * (i + 0.5f) / float(width - 1) - 1.f; float z2 = x * x + y * y; if (z2 <= 1) { float phi = std::atan2(y, x); float theta = std::acos(1 - z2); Vec3f dir(sin(theta) * cos(phi), cos(theta), sin(theta) * sin(phi)); // 1 meter above sea level *p = atmosphere.computeIncidentLight(Vec3f(0, atmosphere.earthRadius + 1, 0), dir, 0, kInfinity); } } fprintf(stderr, "\b\b\b\b\%3d%c", (int)(100 * j / (width - 1)), '%'); } #else // [comment] // Render from a normal camera // [/comment] const unsigned width = 640, height = 480; Vec3f *image = new Vec3f[width * height], *p = image; memset(image, 0x0, sizeof(Vec3f) * width * height); float aspectRatio = width / float(height); float fov = 65; float angle = std::tan(fov * M_PI / 180 * 0.5f); unsigned numPixelSamples = 4; Vec3f orig(0, atmosphere.earthRadius + 1000, 30000); // camera position std::default_random_engine generator; std::uniform_real_distribution<float> distribution(0, 1); // to generate random floats in the range [0:1] for (unsigned y = 0; y < height; ++y) { for (unsigned x = 0; x < width; ++x, ++p) { for (unsigned m = 0; m < numPixelSamples; ++m) { for (unsigned n = 0; n < numPixelSamples; ++n) { float rayx = (2 * (x + (m + distribution(generator)) / numPixelSamples) / float(width) - 1) * aspectRatio * angle; float rayy = (1 - (y + (n + distribution(generator)) / numPixelSamples) / float(height) * 2) * angle; Vec3f dir(rayx, rayy, -1); normalize(dir); // [comment] // Does the ray intersect the planetory body? (the intersection test is against the Earth here // not against the atmosphere). If the ray intersects the Earth body and that the intersection // is ahead of us, then the ray intersects the planet in 2 points, t0 and t1. But we // only want to comupute the atmosphere between t=0 and t=t0 (where the ray hits // the Earth first). If the viewing ray doesn't hit the Earth, or course the ray // is then bounded to the range [0:INF]. In the method computeIncidentLight() we then // compute where this primary ray intersects the atmosphere and we limit the max t range // of the ray to the point where it leaves the atmosphere. // [/comment] float t0, t1, tMax = kInfinity; if (raySphereIntersect(orig, dir, atmosphere.earthRadius, t0, t1) && t1 > 0) tMax = std::max(0.f, t0); // [comment] // The *viewing or camera ray* is bounded to the range [0:tMax] // [/comment] *p += atmosphere.computeIncidentLight(orig, dir, 0, tMax); } } *p *= 1.f / (numPixelSamples * numPixelSamples); } fprintf(stderr, "\b\b\b\b%3d%c", (int)(100 * y / (width - 1)), '%'); } #endif std::cout << "\b\b\b\b" << ((std::chrono::duration<float>)(std::chrono::high_resolution_clock::now() - t0)).count() << " seconds" << std::endl; // Save result to a PPM image (keep these flags if you compile under Windows) std::ofstream ofs(filename, std::ios::out | std::ios::binary); ofs << "P6\n" << width << " " << height << "\n255\n"; p = image; for (unsigned j = 0; j < height; ++j) { for (unsigned i = 0; i < width; ++i, ++p) { #if 1 // Apply tone mapping function (*p)[0] = (*p)[0] < 1.413f ? pow((*p)[0] * 0.38317f, 1.0f / 2.2f) : 1.0f - exp(-(*p)[0]); (*p)[1] = (*p)[1] < 1.413f ? pow((*p)[1] * 0.38317f, 1.0f / 2.2f) : 1.0f - exp(-(*p)[1]); (*p)[2] = (*p)[2] < 1.413f ? pow((*p)[2] * 0.38317f, 1.0f / 2.2f) : 1.0f - exp(-(*p)[2]); #endif ofs << (unsigned char)(std::min(1.f, (*p)[0]) * 255) << (unsigned char)(std::min(1.f, (*p)[1]) * 255) << (unsigned char)(std::min(1.f, (*p)[2]) * 255); } } ofs.close(); delete[] image; } int main() { #if 1 // [comment] // Render a sequence of images (sunrise to sunset) // [/comment] unsigned nangles = 128; for (unsigned i = 0; i < nangles; ++i) { char filename[1024]; sprintf(filename, "./skydome.%04d.ppm", i); float angle = i / float(nangles - 1) * M_PI * 0.6; fprintf(stderr, "Rendering image %d, angle = %0.2f\n", i, angle * 180 / M_PI); renderSkydome(Vec3f(0, cos(angle), -sin(angle)), filename); } #else // [comment] // Render one single image // [/comment] float angle = M_PI * 0; Vec3f sunDir(0, std::cos(angle), -std::sin(angle)); std::cerr << "Sun direction: " << sunDir << std::endl; renderSkydome(sunDir, "./skydome.ppm"); #endif return 0; }