I'll create an advanced path tracer with a comprehensive feature set. This is a substantial single-file implementation covering the core rendering pipeline plus configurable extras. ## Project Structure ``` pathtracer/ ├── main.cpp ├── config.ini └── Makefile ``` --- ## `main.cpp` ```cpp // ============================================================================ // Advanced Monte Carlo Path Tracer // Single-file C++17 implementation // // Features // -------- // * Unidirectional path tracing with next-event estimation (MIS) // * BVH acceleration structure (SAH split, binned) // * Materials: Lambertian, GGX metal, rough dielectric, plastic (coated), // thin dielectric, emissive // * Microfacet BSDF with Smith height-correlated masking-shadowing // * Multiple importance sampling (BSDF + light sampling, power heuristic) // * Russian roulette with adjustable depth / throughput threshold // * Depth of field (thin lens), configurable bokeh blade count // * Environment lighting: constant, gradient sky, or equirect HDR (.hdr) // * Volumetric homogeneous participating medium (global fog) // * Stratified / low-discrepancy sampling (Sobol + Owen scrambling) // * Blue-noise-ish pixel offsetting via Cranley-Patterson rotation // * Adaptive sampling based on per-pixel variance // * Firefly clamping, tonemapping (Reinhard/ACES/Uncharted2/Linear) // * Multithreaded tile renderer with progress reporting // * OBJ mesh loading with smooth normals, transforms // * Output: PPM and PNG (via tiny embedded PNG writer) // * Full INI configuration file // // Build: g++ -O3 -march=native -std=c++17 -pthread main.cpp -o pt // Run: ./pt config.ini // ============================================================================ #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include // ============================================================================ // SECTION 1 — Math primitives // ============================================================================ namespace ptm { static constexpr float PI = 3.14159265358979323846f; static constexpr float INV_PI = 0.31830988618379067154f; static constexpr float TWO_PI = 6.28318530717958647692f; static constexpr float INV_4PI = 0.07957747154594766788f; static constexpr float EPS = 1e-4f; static constexpr float INF = std::numeric_limits::infinity(); inline float clampf(float v, float a, float b) { return v < a ? a : (v > b ? b : v); } inline float saturate(float v) { return clampf(v, 0.0f, 1.0f); } inline float sqr(float v) { return v * v; } inline float lerpf(float a, float b, float t) { return a + (b - a) * t; } inline float radians(float deg) { return deg * PI / 180.0f; } struct Vec3 { float x = 0, y = 0, z = 0; Vec3() = default; Vec3(float v) : x(v), y(v), z(v) {} Vec3(float a, float b, float c) : x(a), y(b), z(c) {} Vec3 operator-() const { return {-x, -y, -z}; } Vec3 operator+(const Vec3& o) const { return {x + o.x, y + o.y, z + o.z}; } Vec3 operator-(const Vec3& o) const { return {x - o.x, y - o.y, z - o.z}; } Vec3 operator*(const Vec3& o) const { return {x * o.x, y * o.y, z * o.z}; } Vec3 operator/(const Vec3& o) const { return {x / o.x, y / o.y, z / o.z}; } Vec3 operator*(float s) const { return {x * s, y * s, z * s}; } Vec3 operator/(float s) const { float i = 1.0f / s; return {x * i, y * i, z * i}; } Vec3& operator+=(const Vec3& o) { x += o.x; y += o.y; z += o.z; return *this; } Vec3& operator-=(const Vec3& o) { x -= o.x; y -= o.y; z -= o.z; return *this; } Vec3& operator*=(const Vec3& o) { x *= o.x; y *= o.y; z *= o.z; return *this; } Vec3& operator*=(float s) { x *= s; y *= s; z *= s; return *this; } Vec3& operator/=(float s) { float i = 1.0f / s; x *= i; y *= i; z *= i; return *this; } float operator[](int i) const { return (&x)[i]; } float& operator[](int i) { return (&x)[i]; } bool isZero() const { return x == 0 && y == 0 && z == 0; } float maxComp() const { return std::max(x, std::max(y, z)); } float minComp() const { return std::min(x, std::min(y, z)); } float average() const { return (x + y + z) * (1.0f / 3.0f); } // Rec.709 relative luminance float luminance() const { return 0.212671f * x + 0.715160f * y + 0.072169f * z; } bool hasNaN() const { return std::isnan(x) || std::isnan(y) || std::isnan(z); } }; inline Vec3 operator*(float s, const Vec3& v) { return v * s; } inline float dot(const Vec3& a, const Vec3& b) { return a.x * b.x + a.y * b.y + a.z * b.z; } inline float absDot(const Vec3& a, const Vec3& b) { return std::fabs(dot(a, b)); } inline Vec3 cross(const Vec3& a, const Vec3& b) { return {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 lengthSq(const Vec3& v) { return dot(v, v); } inline float length(const Vec3& v) { return std::sqrt(dot(v, v)); } inline Vec3 normalize(const Vec3& v) { return v * (1.0f / length(v)); } inline Vec3 vmin(const Vec3& a, const Vec3& b) { return {std::min(a.x,b.x), std::min(a.y,b.y), std::min(a.z,b.z)}; } inline Vec3 vmax(const Vec3& a, const Vec3& b) { return {std::max(a.x,b.x), std::max(a.y,b.y), std::max(a.z,b.z)}; } inline Vec3 vabs(const Vec3& v) { return {std::fabs(v.x), std::fabs(v.y), std::fabs(v.z)}; } inline Vec3 vexp(const Vec3& v) { return {std::exp(v.x), std::exp(v.y), std::exp(v.z)}; } inline Vec3 vsqrt(const Vec3& v) { return {std::sqrt(v.x), std::sqrt(v.y), std::sqrt(v.z)}; } inline Vec3 vpow(const Vec3& v, float p) { return {std::pow(v.x,p), std::pow(v.y,p), std::pow(v.z,p)}; } inline Vec3 lerp(const Vec3& a, const Vec3& b, float t) { return a + (b - a) * t; } inline Vec3 faceforward(const Vec3& n, const Vec3& v) { return dot(n, v) < 0 ? -n : n; } inline Vec3 clampv(const Vec3& v, float a, float b) { return {clampf(v.x,a,b), clampf(v.y,a,b), clampf(v.z,a,b)}; } inline int maxDim(const Vec3& v) { return (v.x > v.y) ? ((v.x > v.z) ? 0 : 2) : ((v.y > v.z) ? 1 : 2); } struct Vec2 { float x = 0, y = 0; Vec2() = default; Vec2(float a, float b) : x(a), y(b) {} Vec2 operator+(const Vec2& o) const { return {x + o.x, y + o.y}; } Vec2 operator-(const Vec2& o) const { return {x - o.x, y - o.y}; } Vec2 operator*(float s) const { return {x * s, y * s}; } }; // ---- 4x4 matrix (row-major) for object transforms ------------------------- struct Mat4 { float m[4][4]; Mat4() { for (int i=0;i<4;i++) for (int j=0;j<4;j++) m[i][j] = (i==j) ? 1.0f : 0.0f; } static Mat4 translate(const Vec3& t) { Mat4 r; r.m[0][3]=t.x; r.m[1][3]=t.y; r.m[2][3]=t.z; return r; } static Mat4 scale(const Vec3& s) { Mat4 r; r.m[0][0]=s.x; r.m[1][1]=s.y; r.m[2][2]=s.z; return r; } static Mat4 rotate(const Vec3& axisIn, float degrees) { Vec3 a = normalize(axisIn); float s = std::sin(radians(degrees)), c = std::cos(radians(degrees)); Mat4 r; r.m[0][0]=a.x*a.x+(1-a.x*a.x)*c; r.m[0][1]=a.x*a.y*(1-c)-a.z*s; r.m[0][2]=a.x*a.z*(1-c)+a.y*s; r.m[1][0]=a.x*a.y*(1-c)+a.z*s; r.m[1][1]=a.y*a.y+(1-a.y*a.y)*c; r.m[1][2]=a.y*a.z*(1-c)-a.x*s; r.m[2][0]=a.x*a.z*(1-c)-a.y*s; r.m[2][1]=a.y*a.z*(1-c)+a.x*s; r.m[2][2]=a.z*a.z+(1-a.z*a.z)*c; return r; } Mat4 operator*(const Mat4& o) const { Mat4 r; for (int i=0;i<4;i++) for (int j=0;j<4;j++) { float s = 0; for (int k=0;k<4;k++) s += m[i][k]*o.m[k][j]; r.m[i][j] = s; } return r; } Vec3 xformPoint(const Vec3& p) const { float x = m[0][0]*p.x + m[0][1]*p.y + m[0][2]*p.z + m[0][3]; float y = m[1][0]*p.x + m[1][1]*p.y + m[1][2]*p.z + m[1][3]; float z = m[2][0]*p.x + m[2][1]*p.y + m[2][2]*p.z + m[2][3]; float w = m[3][0]*p.x + m[3][1]*p.y + m[3][2]*p.z + m[3][3]; if (w != 1.0f && w != 0.0f) { float iw = 1.0f/w; x*=iw; y*=iw; z*=iw; } return {x, y, z}; } Vec3 xformVector(const Vec3& v) const { return { m[0][0]*v.x + m[0][1]*v.y + m[0][2]*v.z, m[1][0]*v.x + m[1][1]*v.y + m[1][2]*v.z, m[2][0]*v.x + m[2][1]*v.y + m[2][2]*v.z }; } // Inverse-transpose upper 3x3 for normals (Gauss-Jordan on 3x3) Mat4 normalMatrix() const { float a=m[0][0],b=m[0][1],c=m[0][2]; float d=m[1][0],e=m[1][1],f=m[1][2]; float g=m[2][0],h=m[2][1],i=m[2][2]; float det = a*(e*i-f*h) - b*(d*i-f*g) + c*(d*h-e*g); if (std::fabs(det) < 1e-20f) return Mat4(); float id = 1.0f/det; Mat4 r; // inverse then transpose => cofactor matrix / det (no transpose needed twice) r.m[0][0]=(e*i-f*h)*id; r.m[1][0]=(c*h-b*i)*id; r.m[2][0]=(b*f-c*e)*id; r.m[0][1]=(f*g-d*i)*id; r.m[1][1]=(a*i-c*g)*id; r.m[2][1]=(c*d-a*f)*id; r.m[0][2]=(d*h-e*g)*id; r.m[1][2]=(b*g-a*h)*id; r.m[2][2]=(a*e-b*d)*id; r.m[0][3]=r.m[1][3]=r.m[2][3]=0; return r; } }; // ---- Orthonormal basis (Duff et al. branchless) --------------------------- struct Frame { Vec3 s, t, n; Frame() = default; explicit Frame(const Vec3& nn) : n(nn) { float sign = (n.z >= 0.0f) ? 1.0f : -1.0f; float a = -1.0f / (sign + n.z); float b = n.x * n.y * a; s = Vec3(1.0f + sign * n.x * n.x * a, sign * b, -sign * n.x); t = Vec3(b, sign + n.y * n.y * a, -n.y); } Vec3 toLocal(const Vec3& v) const { return {dot(v, s), dot(v, t), dot(v, n)}; } Vec3 toWorld(const Vec3& v) const { return s * v.x + t * v.y + n * v.z; } }; // ---- Sampling helpers ----------------------------------------------------- inline Vec2 sampleConcentricDisk(float u1, float u2) { float sx = 2.0f * u1 - 1.0f, sy = 2.0f * u2 - 1.0f; if (sx == 0.0f && sy == 0.0f) return {0, 0}; float r, theta; if (std::fabs(sx) > std::fabs(sy)) { r = sx; theta = (PI / 4.0f) * (sy / sx); } else { r = sy; theta = (PI / 2.0f) - (PI / 4.0f) * (sx / sy); } return {r * std::cos(theta), r * std::sin(theta)}; } inline Vec3 sampleCosineHemisphere(float u1, float u2) { Vec2 d = sampleConcentricDisk(u1, u2); float z = std::sqrt(std::max(0.0f, 1.0f - d.x * d.x - d.y * d.y)); return {d.x, d.y, z}; } inline Vec3 sampleUniformSphere(float u1, float u2) { float z = 1.0f - 2.0f * u1; float r = std::sqrt(std::max(0.0f, 1.0f - z * z)); float phi = TWO_PI * u2; return {r * std::cos(phi), r * std::sin(phi), z}; } inline Vec3 sampleUniformCone(float u1, float u2, float cosThetaMax) { float cosTheta = (1.0f - u1) + u1 * cosThetaMax; float sinTheta = std::sqrt(std::max(0.0f, 1.0f - cosTheta * cosTheta)); float phi = TWO_PI * u2; return {sinTheta * std::cos(phi), sinTheta * std::sin(phi), cosTheta}; } // Regular n-gon aperture sampling for polygonal bokeh inline Vec2 samplePolygonAperture(float u1, float u2, int blades, float rotation) { if (blades < 3) return sampleConcentricDisk(u1, u2); float nb = float(blades); float t = u1 * nb; int tri = int(t); float ft = t - float(tri); float a0 = TWO_PI * float(tri) / nb + rotation; float a1 = TWO_PI * float(tri + 1) / nb + rotation; Vec2 p0(std::cos(a0), std::sin(a0)); Vec2 p1(std::cos(a1), std::sin(a1)); // Uniform triangle sampling (origin, p0, p1) float su = std::sqrt(ft); float b1 = 1.0f - su, b2 = u2 * su; return { p0.x * b1 + p1.x * b2, p0.y * b1 + p1.y * b2 }; } inline float powerHeuristic(float nf, float fPdf, float ng, float gPdf) { float f = nf * fPdf, g = ng * gPdf; float d = f * f + g * g; return (d <= 0.0f) ? 0.0f : (f * f) / d; } inline float balanceHeuristic(float nf, float fPdf, float ng, float gPdf) { float d = nf * fPdf + ng * gPdf; return (d <= 0.0f) ? 0.0f : (nf * fPdf) / d; } // ============================================================================ // SECTION 2 — Random number generation & QMC // ============================================================================ // PCG32 — small, fast, high quality struct PCG32 { uint64_t state = 0x853c49e6748fea9bULL; uint64_t inc = 0xda3e39cb94b95bdbULL; void seed(uint64_t initState, uint64_t seq = 1) { state = 0; inc = (seq << 1u) | 1u; nextUInt(); state += initState; nextUInt(); } uint32_t nextUInt() { uint64_t old = state; state = old * 6364136223846793005ULL + inc; uint32_t xorshifted = uint32_t(((old >> 18u) ^ old) >> 27u); uint32_t rot = uint32_t(old >> 59u); return (xorshifted >> rot) | (xorshifted << ((~rot + 1u) & 31)); } float nextFloat() { // in [0,1) return std::min(0.99999994f, float(nextUInt()) * 2.3283064365386963e-10f); } Vec2 nextFloat2() { float a = nextFloat(); float b = nextFloat(); return {a, b}; } }; // ---- Sobol (first 5 dimensions via direction numbers) + Owen scrambling --- namespace sobol { // Direction numbers for dimensions 0..4 (Joe & Kuo). 32-bit. static const uint32_t kMatrices[5][32] = { // dim 0 : van der Corput (bit reverse) {0x80000000u,0x40000000u,0x20000000u,0x10000000u,0x08000000u,0x04000000u,0x02000000u,0x01000000u, 0x00800000u,0x00400000u,0x00200000u,0x00100000u,0x00080000u,0x00040000u,0x00020000u,0x00010000u, 0x00008000u,0x00004000u,0x00002000u,0x00001000u,0x00000800u,0x00000400u,0x00000200u,0x00000100u, 0x00000080u,0x00000040u,0x00000020u,0x00000010u,0x00000008u,0x00000004u,0x00000002u,0x00000001u}, // dim 1 {0x80000000u,0xc0000000u,0xa0000000u,0xf0000000u,0x88000000u,0xcc000000u,0xaa000000u,0xff000000u, 0x80800000u,0xc0c00000u,0xa0a00000u,0xf0f00000u,0x88880000u,0xcccc0000u,0xaaaa0000u,0xffff0000u, 0x80008000u,0xc000c000u,0xa000a000u,0xf000f000u,0x88008800u,0xcc00cc00u,0xaa00aa00u,0xff00ff00u, 0x80808080u,0xc0c0c0c0u,0xa0a0a0a0u,0xf0f0f0f0u,0x88888888u,0xccccccccu,0xaaaaaaaau,0xffffffffu}, // dim 2 {0x80000000u,0xc0000000u,0x60000000u,0x90000000u,0xe8000000u,0x5c000000u,0x8e000000u,0xc5000000u, 0x68800000u,0x9cc00000u,0xee600000u,0x55900000u,0x80680000u,0xc09c0000u,0x60ee0000u,0x90550000u, 0xe8808000u,0x5cc0c000u,0x8e606000u,0xc5909000u,0x6868e800u,0x9c9c5c00u,0xeeee8e00u,0x5555c500u, 0x8000e880u,0xc0005cc0u,0x60008e60u,0x9000c590u,0xe8006868u,0x5c009c9cu,0x8e00eeeeu,0xc5005555u}, // dim 3 {0x80000000u,0xc0000000u,0x20000000u,0x50000000u,0xf8000000u,0x74000000u,0xa2000000u,0x93000000u, 0xd8800000u,0x25400000u,0x59e00000u,0xe6d00000u,0x78080000u,0xb40c0000u,0x82020000u,0xc3050000u, 0x208f8000u,0x51474000u,0xfbea2000u,0x75d93000u,0xa0858800u,0x914e5400u,0xdbe79e00u,0x25db6d00u, 0x58800080u,0xe54000c0u,0x79e00020u,0xb6d00050u,0x800800f8u,0xc00c0074u,0x200200a2u,0x50050093u}, // dim 4 {0x80000000u,0x40000000u,0xa0000000u,0x50000000u,0x88000000u,0x24000000u,0x12000000u,0x2d000000u, 0x76800000u,0xb3400000u,0x89a00000u,0x2ed00000u,0x18880000u,0x2c240000u,0x76120000u,0xb32d0000u, 0x88f68000u,0x24b34000u,0x1209a000u,0x2d0ed000u,0x76980800u,0xb3480c00u,0x89b21200u,0x2eed2d00u, 0x18000680u,0x2c000b40u,0x76000ea0u,0xb3000d50u,0x88000208u,0x24000d0cu,0x12000f12u,0x2d000d2du} }; inline uint32_t sample(uint32_t index, int dim) { dim = dim % 5; uint32_t r = 0; for (int b = 0; index; ++b, index >>= 1) if (index & 1) r ^= kMatrices[dim][b]; return r; } // Laine-Karras hash for Owen scrambling inline uint32_t lkHash(uint32_t x, uint32_t seed) { x ^= x * 0x3d20adeau; x += seed; x *= (seed >> 16) | 1u; x ^= x * 0x05526c56u; x ^= x * 0x53a22864u; return x; } inline uint32_t reverseBits(uint32_t x) { x = (x << 16) | (x >> 16); x = ((x & 0x00ff00ffu) << 8) | ((x & 0xff00ff00u) >> 8); x = ((x & 0x0f0f0f0fu) << 4) | ((x & 0xf0f0f0f0u) >> 4); x = ((x & 0x33333333u) << 2) | ((x & 0xccccccccu) >> 2); x = ((x & 0x55555555u) << 1) | ((x & 0xaaaaaaaau) >> 1); return x; } inline uint32_t owenScramble(uint32_t x, uint32_t seed) { return reverseBits(lkHash(reverseBits(x), seed)); } inline uint32_t hashCombine(uint32_t a, uint32_t b) { a ^= b + 0x9e3779b9u + (a << 6) + (a >> 2); return a; } } // namespace sobol enum class SamplerType { Random, Stratified, Sobol }; // Unified sampler: dimension-aware, falls back to PCG for high dimensions. class Sampler { public: void init(SamplerType type, uint32_t seed, int sppTotal) { type_ = type; baseSeed_ = seed; spp_ = std::max(1, sppTotal); strataDim_ = int(std::sqrt(float(spp_))); if (strataDim_ * strataDim_ != spp_) strataDim_ = 0; // not a perfect square } void startPixel(int px, int py, int sampleIndex) { px_ = px; py_ = py; sIdx_ = uint32_t(sampleIndex); dim_ = 0; pixelHash_ = sobol::hashCombine(sobol::hashCombine(uint32_t(px) * 73856093u, uint32_t(py) * 19349663u), baseSeed_); rng_.seed(sobol::hashCombine(pixelHash_, sIdx_ * 0x9e3779b9u), 1); } float get1D() { int d = dim_++; switch (type_) { case SamplerType::Random: return rng_.nextFloat(); case SamplerType::Stratified: if (d == 0 && strataDim_) { int sx = int(sIdx_) % strataDim_; return (float(sx) + rng_.nextFloat()) / float(strataDim_); } if (d == 1 && strataDim_) { int sy = int(sIdx_) / strataDim_; return (float(sy) + rng_.nextFloat()) / float(strataDim_); } return rng_.nextFloat(); case SamplerType::Sobol: { if (d < 5) { uint32_t idx = sobol::owenScramble(sIdx_, sobol::hashCombine(pixelHash_, 0x51633e2du)); uint32_t v = sobol::sample(idx, d); v = sobol::owenScramble(v, sobol::hashCombine(pixelHash_, uint32_t(d) * 0x68bc21ebu)); return std::min(0.99999994f, float(v) * 2.3283064365386963e-10f); } return rng_.nextFloat(); } } return rng_.nextFloat(); } Vec2 get2D() { float a = get1D(); float b = get1D(); return {a, b}; } void setDimension(int d) { dim_ = d; } PCG32& rng() { return rng_; } private: SamplerType type_ = SamplerType::Sobol; uint32_t baseSeed_ = 1u, pixelHash_ = 0u, sIdx_ = 0u; int px_ = 0, py_ = 0, dim_ = 0, spp_ = 1, strataDim_ = 0; PCG32 rng_; }; // ============================================================================ // SECTION 3 — Geometry: rays, AABB, primitives, BVH // ============================================================================ struct Ray { Vec3 o, d; float tMin = EPS; float tMax = INF; Ray() = default; Ray(const Vec3& oo, const Vec3& dd, float t0 = EPS, float t1 = INF) : o(oo), d(dd), tMin(t0), tMax(t1) {} Vec3 at(float t) const { return o + d * t; } }; struct AABB { Vec3 lo{ INF, INF, INF}; Vec3 hi{-INF, -INF, -INF}; void expand(const Vec3& p) { lo = vmin(lo, p); hi = vmax(hi, p); } void expand(const AABB& b) { lo = vmin(lo, b.lo); hi = vmax(hi, b.hi); } Vec3 diagonal() const { return hi - lo; } Vec3 centroid() const { return (lo + hi) * 0.5f; } bool valid() const { return hi.x >= lo.x && hi.y >= lo.y && hi.z >= lo.z; } float surfaceArea() const { if (!valid()) return 0.0f; Vec3 d = diagonal(); return 2.0f * (d.x * d.y + d.y * d.z + d.z * d.x); } // Slab test; returns true if [t0,t1] overlap survives bool hit(const Vec3& o, const Vec3& invD, float tMin, float tMax) const { float t0 = (lo.x - o.x) * invD.x, t1 = (hi.x - o.x) * invD.x; if (invD.x < 0) std::swap(t0, t1); tMin = t0 > tMin ? t0 : tMin; tMax = t1 < tMax ? t1 : tMax; if (tMax < tMin) return false; t0 = (lo.y - o.y) * invD.y; t1 = (hi.y - o.y) * invD.y; if (invD.y < 0) std::swap(t0, t1); tMin = t0 > tMin ? t0 : tMin; tMax = t1 < tMax ? t1 : tMax; if (tMax < tMin) return false; t0 = (lo.z - o.z) * invD.z; t1 = (hi.z - o.z) * invD.z; if (invD.z < 0) std::swap(t0, t1); tMin = t0 > tMin ? t0 : tMin; tMax = t1 < tMax ? t1 : tMax; return tMax >= tMin; } float hitDistance(const Vec3& o, const Vec3& invD, float tMin, float tMax) const { float t0x = (lo.x - o.x) * invD.x, t1x = (hi.x - o.x) * invD.x; float t0y = (lo.y - o.y) * invD.y, t1y = (hi.y - o.y) * invD.y; float t0z = (lo.z - o.z) * invD.z, t1z = (hi.z - o.z) * invD.z; float tn = std::max(std::max(std::min(t0x,t1x), std::min(t0y,t1y)), std::max(std::min(t0z,t1z), tMin)); float tf = std::min(std::min(std::max(t0x,t1x), std::max(t0y,t1y)), std::min(std::max(t0z,t1z), tMax)); return (tf >= tn) ? tn : INF; } }; enum class PrimType : uint8_t { Sphere, Triangle, Quad, Disk }; // Interaction record at a hit point struct Hit { float t = INF; Vec3 p; // position Vec3 ng; // geometric normal (normalized, outward) Vec3 ns; // shading normal (normalized) Vec2 uv{0, 0}; // surface parameterization int matId = 0; int primId = -1; bool backface = false; }; struct Primitive { PrimType type = PrimType::Sphere; int matId = 0; int lightId = -1; // index into light list if emissive area light // Sphere: c = center, r = radius // Triangle: v0,v1,v2 + optional normals n0,n1,n2 // Quad: v0 + edge vectors e1, e2 (parallelogram) // Disk: c = center, ng = normal, r = radius Vec3 v0, v1, v2; Vec3 n0, n1, n2; Vec2 t0, t1, t2; float radius = 1.0f; bool hasVertexNormals = false; bool flipNormal = false; AABB bounds() const { AABB b; switch (type) { case PrimType::Sphere: b.expand(v0 - Vec3(radius)); b.expand(v0 + Vec3(radius)); break; case PrimType::Triangle: b.expand(v0); b.expand(v1); b.expand(v2); break; case PrimType::Quad: b.expand(v0); b.expand(v0 + v1); b.expand(v0 + v2); b.expand(v0 + v1 + v2); break; case PrimType::Disk: { Vec3 n = normalize(n0); Vec3 e = vsqrt(vmax(Vec3(0.0f), Vec3(1.0f) - n * n)) * radius; b.expand(v0 - e); b.expand(v0 + e); break; } } // Guard against zero-thickness boxes Vec3 d = b.diagonal(); for (int i = 0; i < 3; ++i) if (d[i] < 1e-6f) { b.lo[i] -= 1e-5f; b.hi[i] += 1e-5f; } return b; } float area() const { switch (type) { case PrimType::Sphere: return 4.0f * PI * radius * radius; case PrimType::Triangle: return 0.5f * length(cross(v1 - v0, v2 - v0)); case PrimType::Quad: return length(cross(v1, v2)); case PrimType::Disk: return PI * radius * radius; } return 0.0f; } bool intersect(const Ray& r, Hit& h) const { switch (type) { case PrimType::Sphere: { Vec3 oc = r.o - v0; float b = dot(oc, r.d); float c = lengthSq(oc) - radius * radius; float disc = b * b - lengthSq(r.d) * c; if (disc < 0) return false; float sq = std::sqrt(disc); float a = lengthSq(r.d); float t = (-b - sq) / a; if (t < r.tMin || t > r.tMax) { t = (-b + sq) / a; if (t < r.tMin || t > r.tMax) return false; } h.t = t; h.p = r.at(t); Vec3 n = normalize(h.p - v0); // Re-project onto sphere to reduce error h.p = v0 + n * radius; h.backface = (c < 0.0f); h.ng = flipNormal ? -n : n; h.ns = h.ng; float phi = std::atan2(n.z, n.x); if (phi < 0) phi += TWO_PI; h.uv = { phi / TWO_PI, std::acos(clampf(n.y, -1.0f, 1.0f)) * INV_PI }; return true; } case PrimType::Triangle: { // Watertight Möller–Trumbore Vec3 e1 = v1 - v0, e2 = v2 - v0; Vec3 pv = cross(r.d, e2); float det = dot(e1, pv); if (std::fabs(det) < 1e-12f) return false; float invDet = 1.0f / det; Vec3 tv = r.o - v0; float u = dot(tv, pv) * invDet; if (u < -1e-7f || u > 1.0f + 1e-7f) return false; Vec3 qv = cross(tv, e1); float v = dot(r.d, qv) * invDet; if (v < -1e-7f || u + v > 1.0f + 1e-7f) return false; float t = dot(e2, qv) * invDet; if (t < r.tMin || t > r.tMax) return false; h.t = t; h.p = r.at(t); Vec3 gn = normalize(cross(e1, e2)); if (flipNormal) gn = -gn; h.backface = dot(gn, r.d) > 0.0f; h.ng = gn; if (hasVertexNormals) { Vec3 sn = n0 * (1.0f - u - v) + n1 * u + n2 * v; float l = length(sn); h.ns = (l > 1e-8f) ? sn / l : gn; if (flipNormal) h.ns = -h.ns; } else h.ns = gn; h.uv = { t0.x * (1 - u - v) + t1.x * u + t2.x * v, t0.y * (1 - u - v) + t1.y * u + t2.y * v }; return true; } case PrimType::Quad: { Vec3 n = cross(v1, v2); float area2 = lengthSq(n); if (area2 < 1e-16f) return false; Vec3 nn = n / std::sqrt(area2); float denom = dot(nn, r.d); if (std::fabs(denom) < 1e-9f) return false; float t = dot(v0 - r.o, nn) / denom; if (t < r.tMin || t > r.tMax) return false; Vec3 p = r.at(t); Vec3 rel = p - v0; // Solve rel = a*v1 + b*v2 via normal cross-products Vec3 w = n / area2; float a = dot(cross(rel, v2), w) * -1.0f; float bb = dot(cross(v1, rel), w) * -1.0f; // Robust: use direct projection a = dot(cross(v2, rel), w); bb = dot(cross(rel, v1), w); if (a < 0 || a > 1 || bb < 0 || bb > 1) return false; h.t = t; h.p = p; Vec3 gn = flipNormal ? -nn : nn; h.backface = dot(gn, r.d) > 0.0f; h.ng = gn; h.ns = gn; h.uv = {a, bb}; return true; } case PrimType::Disk: { Vec3 nn = normalize(n0); float denom = dot(nn, r.d); if (std::fabs(denom) < 1e-9f) return false; float t = dot(v0 - r.o, nn) / denom; if (t < r.tMin || t > r.tMax) return false; Vec3 p = r.at(t); float d2 = lengthSq(p - v0); if (d2 > radius * radius) return false; h.t = t; h.p = p; Vec3 gn = flipNormal ? -nn : nn; h.backface = dot(gn, r.d) > 0.0f; h.ng = gn; h.ns = gn; float rr = std::sqrt(d2) / radius; h.uv = {rr, 0.0f}; return true; } } return false; } // Sample a point on the primitive (for area light sampling). // Returns position + normal; pdf is w.r.t. area. void sampleArea(const Vec2& u, Vec3& pos, Vec3& nrm, float& pdfArea) const { switch (type) { case PrimType::Sphere: { Vec3 d = sampleUniformSphere(u.x, u.y); pos = v0 + d * radius; nrm = flipNormal ? -d : d; pdfArea = 1.0f / area(); break; } case PrimType::Triangle: { float su = std::sqrt(u.x); float b0 = 1.0f - su, b1 = u.y * su; pos = v0 * b0 + v1 * b1 + v2 * (1.0f - b0 - b1); Vec3 gn = normalize(cross(v1 - v0, v2 - v0)); nrm = flipNormal ? -gn : gn; pdfArea = 1.0f / area(); break; } case PrimType::Quad: { pos = v0 + v1 * u.x + v2 * u.y; Vec3 gn = normalize(cross(v1, v2)); nrm = flipNormal ? -gn : gn; pdfArea = 1.0f / area(); break; } case PrimType::Disk: { Vec2 d = sampleConcentricDisk(u.x, u.y); Frame f(normalize(n0)); pos = v0 + (f.s * d.x + f.t * d.y) * radius; nrm = flipNormal ? -normalize(n0) : normalize(n0); pdfArea = 1.0f / area(); break; } } } }; // ---- BVH with binned SAH ------------------------------------------------- class BVH { public: struct Node { AABB box; uint32_t left = 0; // interior: first child index; leaf: first prim index uint16_t count = 0; // 0 => interior node uint8_t axis = 0; uint8_t pad = 0; }; void build(const std::vector& prims, int maxLeafSize = 4, int nBins = 16) { maxLeaf_ = std::max(1, maxLeafSize); nBins_ = std::max(4, nBins); size_t n = prims.size(); indices_.resize(n); std::iota(indices_.begin(), indices_.end(), 0u); bounds_.resize(n); centroids_.resize(n); for (size_t i = 0; i < n; ++i) { bounds_[i] = prims[i].bounds(); centroids_[i] = bounds_[i].centroid(); } nodes_.clear(); nodes_.reserve(std::max(1, n * 2)); nodes_.push_back(Node{}); if (n > 0) buildRecursive(0, 0, uint32_t(n), 0); else { nodes_[0].box = AABB(); nodes_[0].count = 0; nodes_[0].left = 0; } prims_ = &prims; } bool intersect(Ray& r, Hit& hit) const { if (nodes_.empty() || !prims_) return false; Vec3 invD(1.0f / r.d.x, 1.0f / r.d.y, 1.0f / r.d.z); int dirNeg[3] = {invD.x < 0, invD.y < 0, invD.z < 0}; uint32_t stack[64]; int sp = 0; uint32_t node = 0; bool found = false; while (true) { const Node& nd = nodes_[node]; if (nd.box.hit(r.o, invD, r.tMin, r.tMax)) { if (nd.count > 0) { for (uint16_t i = 0; i < nd.count; ++i) { uint32_t pi = indices_[nd.left + i]; Hit tmp; if ((*prims_)[pi].intersect(r, tmp)) { r.tMax = tmp.t; tmp.primId = int(pi); tmp.matId = (*prims_)[pi].matId; hit = tmp; found = true; } } if (sp == 0) break; node = stack[--sp]; } else { // Traverse nearer child first uint32_t c0 = nd.left, c1 = nd.left + 1; if (dirNeg[nd.axis]) std::swap(c0, c1); stack[sp++] = c1; node = c0; } } else { if (sp == 0) break; node = stack[--sp]; } } return found; } // Any-hit query for shadow rays: early-out, skips primitives whose material // is flagged as non-occluding (handled by caller via callback). bool occluded(const Ray& rayIn, const std::function& skipPrim) const { if (nodes_.empty() || !prims_) return false; Ray r = rayIn; Vec3 invD(1.0f / r.d.x, 1.0f / r.d.y, 1.0f / r.d.z); uint32_t stack[64]; int sp = 0; uint32_t node = 0; while (true) { const Node& nd = nodes_[node]; if (nd.box.hit(r.o, invD, r.tMin, r.tMax)) { if (nd.count > 0) { for (uint16_t i = 0; i < nd.count; ++i) { uint32_t pi = indices_[nd.left + i]; if (skipPrim && skipPrim(int(pi))) continue; Hit tmp; if ((*prims_)[pi].intersect(r, tmp)) return true; } if (sp == 0) break; node = stack[--sp]; } else { stack[sp++] = nd.left + 1; node = nd.left; } } else { if (sp == 0) break; node = stack[--sp]; } } return false; } AABB rootBounds() const { return nodes_.empty() ? AABB() : nodes_[0].box; } size_t nodeCount() const { return nodes_.size(); } int depth() const { return maxDepth_; } private: void buildRecursive(uint32_t nodeIdx, uint32_t first, uint32_t count, int depth) { maxDepth_ = std::max(maxDepth_, depth); AABB nodeBox, centBox; for (uint32_t i = 0; i < count; ++i) { nodeBox.expand(bounds_[indices_[first + i]]); centBox.expand(centroids_[indices_[first + i]]); } nodes_[nodeIdx].box = nodeBox; if (count <= uint32_t(maxLeaf_) || depth > 62) { nodes_[nodeIdx].left = first; nodes_[nodeIdx].count = uint16_t(std::min(count, 65535)); return; } int axis = maxDim(centBox.diagonal()); float extent = centBox.diagonal()[axis]; if (extent < 1e-12f) { nodes_[nodeIdx].left = first; nodes_[nodeIdx].count = uint16_t(std::min(count, 65535)); return; } // Binned SAH struct Bin { AABB box; uint32_t n = 0; }; std::vector bins(nBins_); float k = float(nBins_) / extent; float lo = centBox.lo[axis]; for (uint32_t i = 0; i < count; ++i) { uint32_t pi = indices_[first + i]; int b = std::min(nBins_ - 1, int(k * (centroids_[pi][axis] - lo))); b = std::max(0, b); bins[b].n++; bins[b].box.expand(bounds_[pi]); } std::vector leftArea(nBins_, 0.0f), rightArea(nBins_, 0.0f); std::vector leftCount(nBins_, 0), rightCount(nBins_, 0); AABB acc; uint32_t accN = 0; for (int i = 0; i < nBins_ - 1; ++i) { acc.expand(bins[i].box); accN += bins[i].n; leftArea[i] = acc.surfaceArea(); leftCount[i] = accN; } acc = AABB(); accN = 0; for (int i = nBins_ - 1; i > 0; --i) { acc.expand(bins[i].box); accN += bins[i].n; rightArea[i - 1] = acc.surfaceArea(); rightCount[i - 1] = accN; } float bestCost = INF; int bestSplit = -1; for (int i = 0; i < nBins_ - 1; ++i) { if (leftCount[i] == 0 || rightCount[i] == 0) continue; float c = leftArea[i] * float(leftCount[i]) + rightArea[i] * float(rightCount[i]); if (c < bestCost) { bestCost = c; bestSplit = i; } } float leafCost = nodeBox.surfaceArea() * float(count); if (bestSplit < 0 || (count <= 8 && bestCost >= leafCost)) { nodes_[nodeIdx].left = first; nodes_[nodeIdx].count = uint16_t(std::min(count, 65535)); return; } auto midIt = std::partition(indices_.begin() + first, indices_.begin() + first + count, [&](uint32_t pi) { int b = std::min(nBins_ - 1, int(k * (centroids_[pi][axis] - lo))); return std::max(0, b) <= bestSplit; }); uint32_t leftN = uint32_t(midIt - (indices_.begin() + first)); if (leftN == 0 || leftN == count) leftN = count / 2; uint32_t childBase = uint32_t(nodes_.size()); nodes_.push_back(Node{}); nodes_.push_back(Node{}); nodes_[nodeIdx].left = childBase; nodes_[nodeIdx].count = 0; nodes_[nodeIdx].axis = uint8_t(axis); buildRecursive(childBase, first, leftN, depth + 1); buildRecursive(childBase + 1, first + leftN, count - leftN, depth + 1); } std::vector nodes_; std::vector indices_; std::vector bounds_; std::vector centroids_; const std::vector* prims_ = nullptr; int maxLeaf_ = 4, nBins_ = 16, maxDepth_ = 0; }; // ============================================================================ // SECTION 4 — Textures // ============================================================================ struct Image { int w = 0, h = 0; std::vector px; bool valid() const { return w > 0 && h > 0; } const Vec3& at(int x, int y) const { x = std::max(0, std::min(w - 1, x)); y = std::max(0, std::min(h - 1, y)); return px[size_t(y) * w + x]; } Vec3 bilinear(float u, float v) const { if (!valid()) return Vec3(0.0f); float fx = u * float(w) - 0.5f, fy = v * float(h) - 0.5f; int x0 = int(std::floor(fx)), y0 = int(std::floor(fy)); float tx = fx - float(x0), ty = fy - float(y0); Vec3 a = lerp(at(x0, y0), at(x0 + 1, y0), tx); Vec3 b = lerp(at(x0, y0 + 1), at(x0 + 1, y0 + 1), tx); return lerp(a, b, ty); } }; enum class TexKind { Constant, Checker, Gradient, Image, Noise }; struct Texture { TexKind kind = TexKind::Constant; Vec3 c0{1, 1, 1}, c1{0, 0, 0}; float scaleU = 8.0f, scaleV = 8.0f; std::shared_ptr img; uint32_t noiseSeed = 12345u; Vec3 eval(const Vec2& uv, const Vec3& p) const { switch (kind) { case TexKind::Constant: return c0; case TexKind::Checker: { int iu = int(std::floor(uv.x * scaleU)); int iv = int(std::floor(uv.y * scaleV)); return ((iu + iv) & 1) ? c1 : c0; } case TexKind::Gradient: return lerp(c0, c1, saturate(uv.y)); case TexKind::Image: return img ? img->bilinear(uv.x, 1.0f - uv.y) : c0; case TexKind::Noise: { // Cheap value-noise (3 octaves) modulating between c0 and c1 auto hash3 = [](int x, int y, int z, uint32_t s) { uint32_t h = uint32_t(x) * 374761393u + uint32_t(y) * 668265263u + uint32_t(z) * 2147483647u + s; h = (h ^ (h >> 13)) * 1274126177u; return float((h ^ (h >> 16)) & 0xffffffu) / 16777215.0f; }; auto vnoise = [&](Vec3 q) { Vec3 fl(std::floor(q.x), std::floor(q.y), std::floor(q.z)); Vec3 f = q - fl; f = f * f * (Vec3(3.0f) - f * 2.0f); int X = int(fl.x), Y = int(fl.y), Z = int(fl.z); float n = 0.0f; for (int dz = 0; dz < 2; ++dz) for (int dy = 0; dy < 2; ++dy) for (int dx = 0; dx < 2; ++dx) { float wx = dx ? f.x : 1 - f.x; float wy = dy ? f.y : 1 - f.y; float wz = dz ? f.z : 1 - f.z; n += wx * wy * wz * hash3(X + dx, Y + dy, Z + dz, noiseSeed); } return n; }; float amp = 0.5f, freq = scaleU, sum = 0.0f, norm = 0.0f; for (int o = 0; o < 3; ++o) { sum += amp * vnoise(p * freq); norm += amp; amp *= 0.5f; freq *= 2.0f; } return lerp(c0, c1, saturate(sum / std::max(norm, 1e-6f))); } } return c0; } }; // ============================================================================ // SECTION 5 — BSDFs // ============================================================================ enum class MatType : uint8_t { Diffuse, // Lambertian (optional Oren-Nayar roughness) Conductor, // GGX metal with Fresnel (Schlick / complex IOR approx) Dielectric, // Rough/smooth glass with Beer-Lambert absorption ThinDielectric, // Zero-thickness glass sheet Plastic, // Coated diffuse (specular dielectric layer over Lambert) Emissive, // Area light Mirror // Perfect specular }; struct Material { std::string name; MatType type = MatType::Diffuse; Texture albedo; // base colour / tint Vec3 emission{0, 0, 0}; float emissionScale = 1.0f; bool emitBothSides = false; float roughness = 0.1f; // GGX alpha_x source (perceptual) float anisotropy = 0.0f; // [-1,1]; 0 = isotropic float ior = 1.5f; Vec3 specularTint{1, 1, 1}; Vec3 absorption{0, 0, 0}; // per-unit-distance extinction inside glass float specularWeight = 0.04f; // plastic coat reflectance at normal incidence Vec3 k{0, 0, 0}; // conductor extinction coefficient (optional) bool useComplexIOR = false; // Normal / bump perturbation float bumpStrength = 0.0f; bool isEmissive() const { return (type == MatType::Emissive) || (emission.maxComp() * emissionScale > 0.0f); } Vec3 radiance() const { return emission * emissionScale; } bool isSpecularOnly() const { return type == MatType::Mirror || (type == MatType::Dielectric && roughness < 1e-3f) || type == MatType::ThinDielectric; } }; // GGX / Trowbridge-Reitz microfacet distribution with anisotropy namespace ggx { inline void roughnessToAlpha(float rough, float aniso, float& ax, float& ay) { float r = clampf(rough, 1e-3f, 1.0f); float a = r * r; // squared perceptual mapping float asp = std::sqrt(1.0f - clampf(aniso, -0.99f, 0.99f) * 0.9f); ax = std::max(1e-4f, a / asp); ay = std::max(1e-4f, a * asp); } // D(m) for local-frame microfacet normal m (z = geometric normal) inline float D(const Vec3& m, float ax, float ay) { if (m.z <= 0.0f) return 0.0f; float t = sqr(m.x / ax) + sqr(m.y / ay) + sqr(m.z); return 1.0f / (PI * ax * ay * t * t); } // Smith Lambda term inline float lambda(const Vec3& w, float ax, float ay) { float c2 = w.z * w.z; if (c2 >= 1.0f) return 0.0f; float s2 = 1.0f - c2; float a2 = (sqr(w.x * ax) + sqr(w.y * ay)) / std::max(s2, 1e-9f); float tan2 = s2 / std::max(c2, 1e-9f); return 0.5f * (std::sqrt(1.0f + a2 * tan2 / std::max(a2, 1e-9f) * a2) - 1.0f); } // Simpler, numerically-safe Smith G1 for anisotropic GGX inline float G1(const Vec3& w, float ax, float ay) { if (w.z <= 0.0f) return 0.0f; float ax2 = ax * ax, ay2 = ay * ay; float t = std::sqrt(ax2 * w.x * w.x + ay2 * w.y * w.y + w.z * w.z); return 2.0f * w.z / (w.z + t); } // Height-correlated Smith G2 (Heitz) inline float G2(const Vec3& wo, const Vec3& wi, float ax, float ay) { if (wo.z * wi.z <= 0.0f) return 0.0f; float ax2 = ax * ax, ay2 = ay * ay; float lo = std::sqrt(ax2 * wo.x * wo.x + ay2 * wo.y * wo.y + wo.z * wo.z); float li = std::sqrt(ax2 * wi.x * wi.x + ay2 * wi.y * wi.y + wi.z * wi.z); return 2.0f * std::fabs(wo.z) * std::fabs(wi.z) / (std::fabs(wi.z) * lo + std::fabs(wo.z) * li + 1e-9f); } // Sample visible normal distribution (Heitz 2018) inline Vec3 sampleVNDF(const Vec3& wo, float ax, float ay, float u1, float u2) { Vec3 wo_ = normalize(Vec3(ax * wo.x, ay * wo.y, wo.z)); float lensq = wo_.x * wo_.x + wo_.y * wo_.y; Vec3 T1 = (lensq > 0.0f) ? Vec3(-wo_.y, wo_.x, 0.0f) / std::sqrt(lensq) : Vec3(1, 0, 0); Vec3 T2 = cross(wo_, T1); float r = std::sqrt(u1); float phi = TWO_PI * u2; float t1 = r * std::cos(phi); float t2 = r * std::sin(phi); float s = 0.5f * (1.0f + wo_.z); t2 = (1.0f - s) * std::sqrt(std::max(0.0f, 1.0f - t1 * t1)) + s * t2; Vec3 nh = T1 * t1 + T2 * t2 + wo_ * std::sqrt(std::max(0.0f, 1.0f - t1 * t1 - t2 * t2)); Vec3 m = Vec3(ax * nh.x, ay * nh.y, std::max(0.0f, nh.z)); return normalize(m); } // PDF of sampled microfacet normal (w.r.t. solid angle of m) inline float pdfVNDF(const Vec3& wo, const Vec3& m, float ax, float ay) { if (wo.z <= 0.0f) return 0.0f; return G1(wo, ax, ay) * absDot(wo, m) * D(m, ax, ay) / std::max(wo.z, 1e-9f); } } // namespace ggx inline float fresnelDielectric(float cosThetaI, float eta) { cosThetaI = clampf(cosThetaI, -1.0f, 1.0f); if (cosThetaI < 0.0f) { eta = 1.0f / eta; cosThetaI = -cosThetaI; } float sin2ThetaT = (1.0f - cosThetaI * cosThetaI) / (eta * eta); if (sin2ThetaT >= 1.0f) return 1.0f; // TIR float cosThetaT = std::sqrt(1.0f - sin2ThetaT); float rParl = (eta * cosThetaI - cosThetaT) / (eta * cosThetaI + cosThetaT); float rPerp = (cosThetaI - eta * cosThetaT) / (cosThetaI + eta * cosThetaT); return 0.5f * (rParl * rParl + rPerp * rPerp); } inline Vec3 fresnelSchlick(const Vec3& F0, float cosTheta) { float m = clampf(1.0f - cosTheta, 0.0f, 1.0f); float m5 = m * m * m * m * m; return F0 + (Vec3(1.0f) - F0) * m5; } // Exact conductor Fresnel with complex IOR (per channel) inline Vec3 fresnelConductor(float cosThetaI, const Vec3& eta, const Vec3& k) { cosThetaI = clampf(cosThetaI, 0.0f, 1.0f); float c2 = cosThetaI * cosThetaI; float s2 = 1.0f - c2; Vec3 out; for (int i = 0; i < 3; ++i) { float e2 = eta[i] * eta[i], kk2 = k[i] * k[i]; float t0 = e2 - kk2 - s2; float a2b2 = std::sqrt(std::max(0.0f, t0 * t0 + 4.0f * e2 * kk2)); float t1 = a2b2 + c2; float a = std::sqrt(std::max(0.0f, 0.5f * (a2b2 + t0))); float t2 = 2.0f * a * cosThetaI; float Rs = (t1 - t2) / (t1 + t2 + 1e-12f); float t3 = c2 * a2b2 + s2 * s2; float t4 = t2 * s2; float Rp = Rs * (t3 - t4) / (t3 + t4 + 1e-12f); out[i] = clampf(0.5f * (Rs + Rp), 0.0f, 1.0f); } return out; } inline bool refractVec(const Vec3& wi, const Vec3& n, float eta, Vec3& wt) { float cosThetaI = dot(n, wi); float sin2ThetaI = std::max(0.0f, 1.0f - cosThetaI * cosThetaI); float sin2ThetaT = sin2ThetaI / (eta * eta); if (sin2ThetaT >= 1.0f) return false; float cosThetaT = std::sqrt(1.0f - sin2ThetaT); wt = -wi / eta + n * (cosThetaI / eta - cosThetaT); return true; } // ---- BSDF sample / evaluation interface ---------------------------------- struct BSDFSample { Vec3 wi; // world-space outgoing (incident light) direction Vec3 f{0, 0, 0}; // BSDF value (includes tint, excludes cos term) float pdf = 0.0f; // solid-angle pdf bool specular = false; // delta lobe => skip MIS with lights bool transmission = false; float etaScale = 1.0f; // radiance scaling for refraction (eta^2) }; class BSDF { public: BSDF(const Material& m, const Hit& h, const Vec3& woWorld) : mat_(&m), frame_(h.ns), ng_(h.ng), wo_(frame_.toLocal(woWorld)) { base_ = m.albedo.eval(h.uv, h.p); ggx::roughnessToAlpha(m.roughness, m.anisotropy, ax_, ay_); entering_ = dot(woWorld, h.ns) > 0.0f; } bool valid() const { return std::fabs(wo_.z) > 1e-7f; } // Evaluate f(wo,wi) * |cos| is NOT applied here; caller multiplies cos. Vec3 eval(const Vec3& wiWorld) const { Vec3 wi = frame_.toLocal(wiWorld); if (mat_->isSpecularOnly()) return Vec3(0.0f); switch (mat_->type) { case MatType::Diffuse: return evalDiffuse(wo_, wi); case MatType::Conductor:return evalConductor(wo_, wi); case MatType::Plastic: return evalPlastic(wo_, wi); case MatType::Dielectric: return evalRoughDielectric(wo_, wi); default: return Vec3(0.0f); } } float pdf(const Vec3& wiWorld) const { Vec3 wi = frame_.toLocal(wiWorld); if (mat_->isSpecularOnly()) return 0.0f; switch (mat_->type) { case MatType::Diffuse: return (wi.z * wo_.z > 0.0f) ? std::fabs(wi.z) * INV_PI : 0.0f; case MatType::Conductor: { if (wi.z * wo_.z <= 0.0f) return 0.0f; Vec3 woU = wo_.z < 0 ? -wo_ : wo_; Vec3 wiU = wo_.z < 0 ? -wi : wi; Vec3 m = wiU + woU; if (lengthSq(m) < 1e-12f) return 0.0f; m = normalize(m); return ggx::pdfVNDF(woU, m, ax_, ay_) / (4.0f * absDot(woU, m)); } case MatType::Plastic: { if (wi.z * wo_.z <= 0.0f) return 0.0f; float Fo = fresnelDielectric(std::fabs(wo_.z), mat_->ior); float pSpec = clampf(Fo, 0.1f, 0.9f); Vec3 woU = wo_.z < 0 ? -wo_ : wo_; Vec3 wiU = wo_.z < 0 ? -wi : wi; Vec3 m = normalize(wiU + woU); float pS = ggx::pdfVNDF(woU, m, ax_, ay_) / (4.0f * absDot(woU, m) + 1e-9f); float pD = std::fabs(wi.z) * INV_PI; return pSpec * pS + (1.0f - pSpec) * pD; } case MatType::Dielectric: return pdfRoughDielectric(wo_, wi); default: return 0.0f; } } bool sample(Sampler& s, BSDFSample& bs) const { if (!valid()) return false; switch (mat_->type) { case MatType::Diffuse: return sampleDiffuse(s, bs); case MatType::Conductor: return sampleConductor(s, bs); case MatType::Mirror: return sampleMirror(s, bs); case MatType::Plastic: return samplePlastic(s, bs); case MatType::Dielectric: return sampleDielectric(s, bs); case MatType::ThinDielectric: return sampleThinDielectric(s, bs); case MatType::Emissive: return false; } return false; } bool hasDeltaOnly() const { return mat_->isSpecularOnly(); } const Material& material() const { return *mat_; } Vec3 baseColour() const { return base_; } private: // ---- Diffuse (Lambert + optional Oren-Nayar) -------------------------- Vec3 evalDiffuse(const Vec3& wo, const Vec3& wi) const { if (wi.z * wo.z <= 0.0f) return Vec3(0.0f); if (mat_->roughness <= 1e-3f) return base_ * INV_PI; // Oren-Nayar (Qualitative) float sigma = clampf(mat_->roughness, 0.0f, 1.0f) * (PI * 0.5f); float s2 = sigma * sigma; float A = 1.0f - 0.5f * s2 / (s2 + 0.33f); float B = 0.45f * s2 / (s2 + 0.09f); float cosTi = std::fabs(wi.z), cosTo = std::fabs(wo.z); float sinTi = std::sqrt(std::max(0.0f, 1.0f - cosTi * cosTi)); float sinTo = std::sqrt(std::max(0.0f, 1.0f - cosTo * cosTo)); float maxCos = 0.0f; if (sinTi > 1e-4f && sinTo > 1e-4f) { float dcos = (wi.x * wo.x + wi.y * wo.y) / (sinTi * sinTo); maxCos = std::max(0.0f, dcos); } float sinAlphaTanBeta = (cosTi > cosTo) ? sinTo * (sinTi / std::max(cosTi, 1e-4f)) : sinTi * (sinTo / std::max(cosTo, 1e-4f)); return base_ * INV_PI * (A + B * maxCos * sinAlphaTanBeta); } bool sampleDiffuse(Sampler& s, BSDFSample& bs) const { Vec2 u = s.get2D(); Vec3 wi = sampleCosineHemisphere(u.x, u.y); if (wo_.z < 0.0f) wi.z = -wi.z; bs.pdf = std::fabs(wi.z) * INV_PI; if (bs.pdf < 1e-9f) return false; bs.f = evalDiffuse(wo_, wi); bs.wi = frame_.toWorld(wi); bs.specular = false; return true; } // ---- Conductor -------------------------------------------------------- Vec3 conductorFresnel(float cosTheta) const { if (mat_->useComplexIOR) return fresnelConductor(cosTheta, Vec3(mat_->ior), mat_->k); return fresnelSchlick(base_ * mat_->specularTint, cosTheta); } Vec3 evalConductor(const Vec3& wo, const Vec3& wi) const { if (wi.z * wo.z <= 0.0f) return Vec3(0.0f); Vec3 woU = wo.z < 0 ? -wo : wo; Vec3 wiU = wo.z < 0 ? -wi : wi; Vec3 m = wiU + woU; if (lengthSq(m) < 1e-12f) return Vec3(0.0f); m = normalize(m); float Dv = ggx::D(m, ax_, ay_); float Gv = ggx::G2(woU, wiU, ax_, ay_); Vec3 F = conductorFresnel(absDot(woU, m)); return F * (Dv * Gv / (4.0f * std::fabs(woU.z) * std::fabs(wiU.z) + 1e-9f)); } bool sampleConductor(Sampler& s, BSDFSample& bs) const { Vec3 woU = wo_.z < 0 ? -wo_ : wo_; Vec2 u = s.get2D(); if (mat_->roughness < 1e-3f) { Vec3 wi(-woU.x, -woU.y, woU.z); if (wo_.z < 0) wi = -wi; bs.wi = frame_.toWorld(wi); bs.f = conductorFresnel(std::fabs(wo_.z)) / std::max(std::fabs(wi.z), 1e-6f); bs.pdf = 1.0f; bs.specular = true; return true; } Vec3 m = ggx::sampleVNDF(woU, ax_, ay_, u.x, u.y); Vec3 wiU = m * (2.0f * dot(woU, m)) - woU; if (wiU.z <= 0.0f) return false; float pdfM = ggx::pdfVNDF(woU, m, ax_, ay_); bs.pdf = pdfM / (4.0f * absDot(woU, m) + 1e-9f); if (bs.pdf < 1e-9f) return false; Vec3 wi = (wo_.z < 0) ? -wiU : wiU; bs.f = evalConductor(wo_, wi); bs.wi = frame_.toWorld(wi); bs.specular = false; return true; } bool sampleMirror(Sampler&, BSDFSample& bs) const { Vec3 wi(-wo_.x, -wo_.y, wo_.z); bs.wi = frame_.toWorld(wi); Vec3 F = base_ * mat_->specularTint; bs.f = F / std::max(std::fabs(wi.z), 1e-6f); bs.pdf = 1.0f; bs.specular = true; return true; } // ---- Plastic (coated diffuse) --------------------------------------- Vec3 evalPlastic(const Vec3& wo, const Vec3& wi) const { if (wi.z * wo.z <= 0.0f) return Vec3(0.0f); Vec3 woU = wo.z < 0 ? -wo : wo; Vec3 wiU = wo.z < 0 ? -wi : wi; // Specular coat Vec3 spec(0.0f); Vec3 m = wiU + woU; if (lengthSq(m) > 1e-12f) { m = normalize(m); float Dv = ggx::D(m, ax_, ay_); float Gv = ggx::G2(woU, wiU, ax_, ay_); float F = fresnelDielectric(absDot(woU, m), mat_->ior); spec = mat_->specularTint * (F * Dv * Gv / (4.0f * std::fabs(woU.z) * std::fabs(wiU.z) + 1e-9f)); } // Diffuse base attenuated by coat transmission (approximate) float Fo = fresnelDielectric(std::fabs(woU.z), mat_->ior); float Fi = fresnelDielectric(std::fabs(wiU.z), mat_->ior); Vec3 diff = base_ * INV_PI * (1.0f - Fo) * (1.0f - Fi); return spec + diff; } bool samplePlastic(Sampler& s, BSDFSample& bs) const { float Fo = fresnelDielectric(std::fabs(wo_.z), mat_->ior); float pSpec = clampf(Fo, 0.1f, 0.9f); float uc = s.get1D(); Vec2 u = s.get2D(); Vec3 woU = wo_.z < 0 ? -wo_ : wo_; Vec3 wiU; if (uc < pSpec) { Vec3 m = ggx::sampleVNDF(woU, ax_, ay_, u.x, u.y); wiU = m * (2.0f * dot(woU, m)) - woU; if (wiU.z <= 0.0f) return false; } else { wiU = sampleCosineHemisphere(u.x, u.y); } Vec3 wi = (wo_.z < 0) ? -wiU : wiU; bs.wi = frame_.toWorld(wi); bs.f = evalPlastic(wo_, wi); bs.pdf = pdf(bs.wi); bs.specular = false; return bs.pdf > 1e-9f; } // ---- Dielectric (rough & smooth) ------------------------------------- Vec3 evalRoughDielectric(const Vec3& wo, const Vec3& wi) const { if (mat_->roughness < 1e-3f) return Vec3(0.0f); float eta = (wo.z > 0.0f) ? mat_->ior : 1.0f / mat_->ior; bool reflect = wi.z * wo.z > 0.0f; Vec3 m; if (reflect) m = wi + wo; else m = wi * (1.0f / eta) + wo; // half-vector for refraction if (lengthSq(m) < 1e-12f) return Vec3(0.0f); m = normalize(m); if (m.z < 0.0f) m = -m; if (dot(m, wi) * wi.z < 0.0f || dot(m, wo) * wo.z < 0.0f) return Vec3(0.0f); float F = fresnelDielectric(dot(wo, m), mat_->ior); float Dv = ggx::D(m, ax_, ay_); float Gv = ggx::G2(wo.z > 0 ? wo : -wo, wi.z > 0 ? wi : -wi, ax_, ay_); if (reflect) { return mat_->specularTint * (F * Dv * Gv / std::fabs(4.0f * wo.z * wi.z + 1e-9f)); } else { float denom = sqr(dot(wi, m) / eta + dot(wo, m)); float val = (1.0f - F) * Dv * Gv * std::fabs(dot(wi, m) * dot(wo, m)) / (std::fabs(wi.z * wo.z) * denom * eta * eta + 1e-9f); return base_ * val; } } float pdfRoughDielectric(const Vec3& wo, const Vec3& wi) const { if (mat_->roughness < 1e-3f) return 0.0f; float eta = (wo.z > 0.0f) ? mat_->ior : 1.0f / mat_->ior; bool reflect = wi.z * wo.z > 0.0f; Vec3 m = reflect ? (wi + wo) : (wi * (1.0f / eta) + wo); if (lengthSq(m) < 1e-12f) return 0.0f; m = normalize(m); if (m.z < 0.0f) m = -m; float F = fresnelDielectric(dot(wo, m), mat_->ior); Vec3 woU = wo.z > 0 ? wo : -wo; float pdfM = ggx::pdfVNDF(woU, m, ax_, ay_); if (reflect) { return pdfM * F / (4.0f * absDot(wo, m) + 1e-9f); } else { float denom = sqr(dot(wi, m) / eta + dot(wo, m)); float dwh_dwi = absDot(wi, m) / (denom * eta * eta + 1e-9f); return pdfM * (1.0f - F) * dwh_dwi; } } bool sampleDielectric(Sampler& s, BSDFSample& bs) const { float eta = mat_->ior; if (mat_->roughness < 1e-3f) { // Smooth glass: delta reflection/refraction float cosThetaO = wo_.z; float F = fresnelDielectric(cosThetaO, eta); float uc = s.get1D(); if (uc < F) { Vec3 wi(-wo_.x, -wo_.y, wo_.z); bs.wi = frame_.toWorld(wi); bs.f = mat_->specularTint * (F / std::max(std::fabs(wi.z), 1e-6f)); bs.pdf = F; bs.specular = true; return true; } else { Vec3 n(0, 0, 1); float relEta = (cosThetaO > 0.0f) ? eta : 1.0f / eta; Vec3 wt; if (!refractVec(wo_, faceforward(n, wo_), relEta, wt)) return false; bs.wi = frame_.toWorld(wt); bs.f = base_ * ((1.0f - F) / std::max(std::fabs(wt.z), 1e-6f)); // Radiance transport asymmetry bs.f = bs.f / (relEta * relEta); bs.etaScale = relEta * relEta; bs.pdf = 1.0f - F; bs.specular = true; bs.transmission = true; return true; } } // Rough dielectric via VNDF Vec3 woU = wo_.z > 0 ? wo_ : -wo_; Vec2 u = s.get2D(); Vec3 m = ggx::sampleVNDF(woU, ax_, ay_, u.x, u.y); if (wo_.z < 0) m = m; // m is in upper hemisphere; keep float F = fresnelDielectric(dot(wo_, m), eta); float uc = s.get1D(); Vec3 wi; bool refl = (uc < F); if (refl) { wi = m * (2.0f * dot(wo_, m)) - wo_; if (wi.z * wo_.z <= 0.0f) return false; } else { float relEta = (dot(wo_, m) > 0.0f) ? eta : 1.0f / eta; Vec3 mm = faceforward(m, wo_); if (!refractVec(wo_, mm, relEta, wi)) return false; if (wi.z * wo_.z >= 0.0f) return false; bs.transmission = true; bs.etaScale = relEta * relEta; } bs.wi = frame_.toWorld(wi); bs.f = evalRoughDielectric(wo_, wi); bs.pdf = pdfRoughDielectric(wo_, wi); bs.specular = false; return bs.pdf > 1e-9f; } bool sampleThinDielectric(Sampler& s, BSDFSample& bs) const { float F = fresnelDielectric(std::fabs(wo_.z), mat_->ior); // Account for two interfaces if (F < 1.0f) F = 2.0f * F / (1.0f + F); float uc = s.get1D(); if (uc < F) { Vec3 wi(-wo_.x, -wo_.y, wo_.z); bs.wi = frame_.toWorld(wi); bs.f = mat_->specularTint * (F / std::max(std::fabs(wi.z), 1e-6f)); bs.pdf = F; } else { Vec3 wi = -wo_; // straight through bs.wi = frame_.toWorld(wi); bs.f = base_ * ((1.0f - F) / std::max(std::fabs(wi.z), 1e-6f)); bs.pdf = 1.0f - F; bs.transmission = true; } bs.specular = true; return true; } const Material* mat_; Frame frame_; Vec3 ng_, wo_, base_; float ax_ = 0.1f, ay_ = 0.1f; bool entering_ = true; }; // ============================================================================ // SECTION 6 — Lights (area + environment + distant) // ============================================================================ struct AreaLight { int primIndex = -1; Vec3 radiance{1, 1, 1}; float area = 1.0f; bool twoSided = false; }; struct DistantLight { // directional / sun Vec3 direction{0, -1, 0}; // pointing FROM light Vec3 radiance{1, 1, 1}; float angularRadiusDeg = 0.53f; }; enum class EnvMode { None, Constant, Gradient, HDRI }; // 2D piecewise-constant distribution for env-map importance sampling class Distribution2D { public: void build(const std::vector& f, int nu, int nv) { nu_ = nu; nv_ = nv; condCdf_.assign(size_t(nv) * (nu + 1), 0.0f); margCdf_.assign(nv + 1, 0.0f); condInt_.assign(nv, 0.0f); for (int v = 0; v < nv; ++v) { float* cdf = &condCdf_[size_t(v) * (nu + 1)]; cdf[0] = 0.0f; for (int u = 0; u < nu; ++u) cdf[u + 1] = cdf[u] + f[size_t(v) * nu + u] / float(nu); condInt_[v] = cdf[nu]; if (condInt_[v] > 0.0f) for (int u = 1; u <= nu; ++u) cdf[u] /= condInt_[v]; else for (int u = 1; u <= nu; ++u) cdf[u] = float(u) / float(nu); } margCdf_[0] = 0.0f; for (int v = 0; v < nv; ++v) margCdf_[v + 1] = margCdf_[v] + condInt_[v] / float(nv); total_ = margCdf_[nv]; if (total_ > 0.0f) for (int v = 1; v <= nv; ++v) margCdf_[v] /= total_; else for (int v = 1; v <= nv; ++v) margCdf_[v] = float(v) / float(nv); } bool valid() const { return nu_ > 0 && nv_ > 0 && total_ > 0.0f; } void sample(float u1, float u2, float& u, float& v, float& pdf) const { int vi = findInterval(margCdf_, u2); float dv = margCdf_[vi + 1] - margCdf_[vi]; float fv = (dv > 0.0f) ? (u2 - margCdf_[vi]) / dv : 0.0f; v = (float(vi) + fv) / float(nv_); float pdfV = (total_ > 0.0f) ? condInt_[vi] / total_ : 1.0f; const float* cdf = &condCdf_[size_t(vi) * (nu_ + 1)]; int ui = findIntervalPtr(cdf, nu_ + 1, u1); float du = cdf[ui + 1] - cdf[ui]; float fu = (du > 0.0f) ? (u1 - cdf[ui]) / du : 0.0f; u = (float(ui) + fu) / float(nu_); float pdfU = (condInt_[vi] > 0.0f) ? (cdf[ui + 1] - cdf[ui]) * float(nu_) : 1.0f; pdf = pdfV * pdfU; } float pdf(float u, float v) const { if (!valid()) return 1.0f; int ui = std::min(nu_ - 1, std::max(0, int(u * float(nu_)))); int vi = std::min(nv_ - 1, std::max(0, int(v * float(nv_)))); const float* cdf = &condCdf_[size_t(vi) * (nu_ + 1)]; float pdfU = (condInt_[vi] > 0.0f) ? (cdf[ui + 1] - cdf[ui]) * float(nu_) : 1.0f; float pdfV = (total_ > 0.0f) ? condInt_[vi] / total_ : 1.0f; return pdfU * pdfV; } private: static int findInterval(const std::vector& a, float x) { int lo = 0, hi = int(a.size()) - 1; while (hi - lo > 1) { int mid = (lo + hi) / 2; if (a[mid] <= x) lo = mid; else hi = mid; } return std::min(lo, int(a.size()) - 2); } static int findIntervalPtr(const float* a, int n, float x) { int lo = 0, hi = n - 1; while (hi - lo > 1) { int mid = (lo + hi) / 2; if (a[mid] <= x) lo = mid; else hi = mid; } return std::min(lo, n - 2); } int nu_ = 0, nv_ = 0; float total_ = 0.0f; std::vector condCdf_, margCdf_, condInt_; }; class Environment { public: EnvMode mode = EnvMode::Gradient; Vec3 constColour{0.5f, 0.7f, 1.0f}; Vec3 horizonColour{1.0f, 1.0f, 1.0f}; Vec3 zenithColour{0.35f, 0.55f, 1.0f}; Vec3 groundColour{0.25f, 0.22f, 0.20f}; float intensity = 1.0f; float rotationDeg = 0.0f; std::shared_ptr hdri; bool importanceSample = true; bool visibleToCamera = true; void prepare() { if (mode == EnvMode::HDRI && hdri && hdri->valid() && importanceSample) { int nu = hdri->w, nv = hdri->h; std::vector lum(size_t(nu) * nv); for (int v = 0; v < nv; ++v) { float sinTheta = std::sin(PI * (float(v) + 0.5f) / float(nv)); for (int u = 0; u < nu; ++u) lum[size_t(v) * nu + u] = std::max(0.0f, hdri->px[size_t(v) * nu + u].luminance()) * sinTheta; } dist_.build(lum, nu, nv); } } Vec3 eval(const Vec3& dirIn) const { Vec3 d = rotate(dirIn); switch (mode) { case EnvMode::None: return Vec3(0.0f); case EnvMode::Constant: return constColour * intensity; case EnvMode::Gradient: { if (d.y >= 0.0f) { float t = std::pow(saturate(d.y), 0.5f); return lerp(horizonColour, zenithColour, t) * intensity; } float t = saturate(-d.y * 2.0f); return lerp(horizonColour, groundColour, t) * intensity; } case EnvMode::HDRI: { if (!hdri || !hdri->valid()) return Vec3(0.0f); float u, v; dirToUV(d, u, v); return hdri->bilinear(u, v) * intensity; } } return Vec3(0.0f); } // Sample a direction; returns radiance, sets pdf (solid angle) Vec3 sample(const Vec2& u, Vec3& wi, float& pdf) const { if (mode == EnvMode::None) { pdf = 0; return Vec3(0.0f); } if (mode == EnvMode::HDRI && dist_.valid()) { float su, sv, mapPdf; dist_.sample(u.x, u.y, su, sv, mapPdf); float theta = sv * PI, phi = su * TWO_PI; float sinTheta = std::sin(theta); if (sinTheta <= 0.0f) { pdf = 0; return Vec3(0.0f); } Vec3 dLocal(sinTheta * std::cos(phi), std::cos(theta), sinTheta * std::sin(phi)); wi = unrotate(dLocal); pdf = mapPdf / (2.0f * PI * PI * sinTheta); return hdri->bilinear(su, sv) * intensity; } // Uniform sphere fallback wi = sampleUniformSphere(u.x, u.y); pdf = INV_4PI; return eval(wi); } float pdf(const Vec3& dirIn) const { if (mode == EnvMode::None) return 0.0f; if (mode == EnvMode::HDRI && dist_.valid()) { Vec3 d = rotate(dirIn); float u, v; dirToUV(d, u, v); float theta = v * PI; float sinTheta = std::sin(theta); if (sinTheta <= 0.0f) return 0.0f; return dist_.pdf(u, v) / (2.0f * PI * PI * sinTheta); } return INV_4PI; } bool active() const { return mode != EnvMode::None; } private: static void dirToUV(const Vec3& d, float& u, float& v) { float theta = std::acos(clampf(d.y, -1.0f, 1.0f)); float phi = std::atan2(d.z, d.x); if (phi < 0) phi += TWO_PI; u = phi / TWO_PI; v = theta / PI; } Vec3 rotate(const Vec3& d) const { float a = radians(rotationDeg); float c = std::cos(a), s = std::sin(a); return { c * d.x + s * d.z, d.y, -s * d.x + c * d.z }; } Vec3 unrotate(const Vec3& d) const { float a = -radians(rotationDeg); float c = std::cos(a), s = std::sin(a); return { c * d.x + s * d.z, d.y, -s * d.x + c * d.z }; } Distribution2D dist_; }; // ============================================================================ // SECTION 7 — Camera // ============================================================================ struct Camera { Vec3 position{0, 1, 5}; Vec3 lookAt{0, 1, 0}; Vec3 up{0, 1, 0}; float fovDeg = 40.0f; float aperture = 0.0f; // diameter, world units float focusDistance = -1.0f; // <0 => distance to lookAt int bokehBlades = 0; // 0 => circular float bokehRotationDeg = 0.0f; float shutterOpen = 0.0f, shutterClose = 0.0f; // for future motion blur int width = 800, height = 600; void build() { forward_ = normalize(lookAt - position); right_ = normalize(cross(forward_, up)); upv_ = cross(right_, forward_); aspect_ = float(width) / float(height); halfH_ = std::tan(radians(fovDeg) * 0.5f); halfW_ = halfH_ * aspect_; focus_ = (focusDistance > 0.0f) ? focusDistance : length(lookAt - position); } // px,py in [0,width)x[0,height) with subpixel jitter already applied Ray generateRay(float px, float py, Sampler& s) const { float sx = (2.0f * (px / float(width)) - 1.0f) * halfW_; float sy = (1.0f - 2.0f * (py / float(height))) * halfH_; Vec3 dir = normalize(forward_ + right_ * sx + upv_ * sy); Vec3 origin = position; if (aperture > 1e-6f) { Vec2 ap = (bokehBlades >= 3) ? samplePolygonAperture(s.get1D(), s.get1D(), bokehBlades, radians(bokehRotationDeg)) : sampleConcentricDisk(s.get1D(), s.get1D()); float r = aperture * 0.5f; Vec3 offset = right_ * (ap.x * r) + upv_ * (ap.y * r); float ft = focus_ / std::max(dot(dir, forward_), 1e-6f); Vec3 focalPoint = position + dir * ft; origin = position + offset; dir = normalize(focalPoint - origin); } return Ray(origin, dir); } private: Vec3 forward_, right_, upv_; float aspect_ = 1.0f, halfW_ = 1.0f, halfH_ = 1.0f, focus_ = 1.0f; }; // ============================================================================ // SECTION 8 — Medium (homogeneous global fog) // ============================================================================ struct Medium { bool enabled = false; Vec3 sigmaA{0.0f, 0.0f, 0.0f}; // absorption Vec3 sigmaS{0.0f, 0.0f, 0.0f}; // scattering float g = 0.0f; // Henyey-Greenstein asymmetry int maxScatterEvents = 8; Vec3 sigmaT() const { return sigmaA + sigmaS; } bool isScattering() const { return enabled && sigmaS.maxComp() > 0.0f; } static float hgPhase(float cosTheta, float g) { float d = 1.0f + g * g + 2.0f * g * cosTheta; return INV_4PI * (1.0f - g * g) / (d * std::sqrt(std::max(d, 1e-9f))); } static Vec3 sampleHG(const Vec3& wo, float g, float u1, float u2, float& pdf) { float cosTheta; if (std::fabs(g) < 1e-3f) cosTheta = 1.0f - 2.0f * u1; else { float s = (1.0f - g * g) / (1.0f + g - 2.0f * g * u1); cosTheta = -(1.0f + g * g - s * s) / (2.0f * g); } cosTheta = clampf(cosTheta, -1.0f, 1.0f); float sinTheta = std::sqrt(std::max(0.0f, 1.0f - cosTheta * cosTheta)); float phi = TWO_PI * u2; Frame f(wo); Vec3 wi = f.toWorld(Vec3(sinTheta * std::cos(phi), sinTheta * std::sin(phi), cosTheta)); pdf = hgPhase(cosTheta, g); return wi; } }; // ============================================================================ // SECTION 9 — Scene // ============================================================================ struct RenderSettings { int width = 800, height = 600; int samplesPerPixel = 256; int maxDepth = 24; int minDepthRR = 3; float rrThreshold = 0.05f; bool useNEE = true; bool useMIS = true; int lightSamplesPerBounce = 1; float clampIndirect = 0.0f; // 0 => disabled float clampDirect = 0.0f; int tileSize = 32; int threads = 0; // 0 => hardware_concurrency SamplerType sampler = SamplerType::Sobol; uint32_t seed = 1234u; // Adaptive sampling bool adaptive = false; int adaptiveMinSpp = 32; float adaptiveTolerance = 0.005f; // relative variance target int adaptiveBatch = 16; // Output std::string outputFile = "render.png"; std::string tonemap = "aces"; // linear|reinhard|aces|uncharted2|filmic float exposure = 0.0f; // EV float gamma = 2.2f; bool writeAux = false; // dump normal/albedo/depth buffers bool denoiseSimple = false; // small joint bilateral pass float denoiseSigma = 1.2f; // Debug std::string debugMode = "none"; // none|normal|albedo|depth|uv|bvh|heat bool progressBar = true; int previewEvery = 0; // write partial image every N spp (0=off) }; class Scene { public: std::vector prims; std::vector materials; std::vector areaLights; std::vector distantLights; Environment env; Medium medium; Camera camera; RenderSettings settings; BVH bvh; int addMaterial(const Material& m) { materials.push_back(m); matIndex_[m.name] = int(materials.size()) - 1; return int(materials.size()) - 1; } int findMaterial(const std::string& n) const { auto it = matIndex_.find(n); return (it == matIndex_.end()) ? -1 : it->second; } void finalize() { // Register emissive primitives as area lights areaLights.clear(); for (size_t i = 0; i < prims.size(); ++i) { const Material& m = materials[prims[i].matId]; if (m.isEmissive()) { AreaLight L; L.primIndex = int(i); L.radiance = m.radiance(); L.area = prims[i].area(); L.twoSided = m.emitBothSides; prims[i].lightId = int(areaLights.size()); areaLights.push_back(L); } else prims[i].lightId = -1; } bvh.build(prims, 4, 16); env.prepare(); camera.width = settings.width; camera.height = settings.height; camera.build(); // Precompute light selection weights (uniform + power-proportional blend) lightPower_.clear(); float total = 0.0f; for (const auto& L : areaLights) { float p = L.radiance.luminance() * L.area * PI * (L.twoSided ? 2.0f : 1.0f); lightPower_.push_back(p); total += p; } for (const auto& D : distantLights) { float p = D.radiance.luminance(); lightPower_.push_back(p); total += p; } if (env.active()) { AABB b = bvh.rootBounds(); float r = b.valid() ? std::max(1.0f, length(b.diagonal()) * 0.5f) : 1.0f; float p = env.eval(Vec3(0, 1, 0)).luminance() * 4.0f * PI * r * r * 0.05f; p = std::max(p, 1e-3f); lightPower_.push_back(p); total += p; } lightCdf_.assign(lightPower_.size() + 1, 0.0f); if (total > 0.0f) { for (size_t i = 0; i < lightPower_.size(); ++i) lightCdf_[i + 1] = lightCdf_[i] + lightPower_[i] / total; lightCdf_.back() = 1.0f; } else { for (size_t i = 0; i <= lightPower_.size(); ++i) lightCdf_[i] = float(i) / float(std::max(1, lightPower_.size())); } nAreaLights_ = int(areaLights.size()); nDistant_ = int(distantLights.size()); hasEnvLight_ = env.active(); nLightsTotal_ = nAreaLights_ + nDistant_ + (hasEnvLight_ ? 1 : 0); } int numLights() const { return nLightsTotal_; } bool hasLights() const { return nLightsTotal_ > 0; } // Select a light index proportional to estimated power. int selectLight(float u, float& pmf) const { if (nLightsTotal_ == 0) { pmf = 0; return -1; } int lo = 0, hi = int(lightCdf_.size()) - 1; while (hi - lo > 1) { int mid = (lo + hi) / 2; if (lightCdf_[mid] <= u) lo = mid; else hi = mid; } int idx = std::min(lo, nLightsTotal_ - 1); pmf = lightCdf_[idx + 1] - lightCdf_[idx]; if (pmf <= 0.0f) { pmf = 1.0f / float(nLightsTotal_); } return idx; } float lightPmf(int idx) const { if (idx < 0 || idx >= nLightsTotal_) return 0.0f; float p = lightCdf_[idx + 1] - lightCdf_[idx]; return (p > 0.0f) ? p : 1.0f / float(nLightsTotal_); } int areaLightCount() const { return nAreaLights_; } int distantOffset() const { return nAreaLights_; } int envIndex() const { return hasEnvLight_ ? nAreaLights_ + nDistant_ : -1; } private: std::map matIndex_; std::vector lightPower_, lightCdf_; int nAreaLights_ = 0, nDistant_ = 0, nLightsTotal_ = 0; bool hasEnvLight_ = false; }; // ============================================================================ // SECTION 10 — Integrator // ============================================================================ struct LightSample { Vec3 wi; // direction from shading point toward light Vec3 Li{0, 0, 0}; // incident radiance float pdf = 0.0f; // solid-angle pdf INCLUDING light-selection pmf float distance = INF; int primId = -1; bool isDelta = false; }; class Integrator { public: explicit Integrator(const Scene& sc) : scene_(sc) { sceneRadius_ = 1e4f; AABB b = sc.bvh.rootBounds(); if (b.valid()) sceneRadius_ = std::max(1.0f, length(b.diagonal())); } struct AuxData { Vec3 normal{0, 0, 0}; Vec3 albedo{0, 0, 0}; float depth = 0.0f; int bvhSteps = 0; }; Vec3 radiance(Ray ray, Sampler& sampler, AuxData* aux = nullptr) const { const RenderSettings& S = scene_.settings; Vec3 L(0.0f), beta(1.0f); float etaScaleAccum = 1.0f; bool specularBounce = true; // camera ray behaves like specular float prevBsdfPdf = 0.0f; Vec3 prevPos = ray.o; int scatterEvents = 0; for (int depth = 0; depth < S.maxDepth; ++depth) { Hit hit; bool found = scene_.bvh.intersect(ray, hit); // ---- Participating medium ------------------------------------ if (scene_.medium.enabled) { float tMaxSeg = found ? hit.t : sceneRadius_ * 4.0f; Vec3 st = scene_.medium.sigmaT(); if (st.maxComp() > 0.0f && scatterEvents < scene_.medium.maxScatterEvents) { // Sample a channel for distance sampling int ch = std::min(2, int(sampler.get1D() * 3.0f)); float sigma = st[ch]; if (sigma > 0.0f) { float t = -std::log(1.0f - sampler.get1D()) / sigma; if (t < tMaxSeg) { // Volumetric scattering event Vec3 pMed = ray.at(t); Vec3 Tr = vexp(-st * t); float pdfDist = (Tr[0] * st[0] + Tr[1] * st[1] + Tr[2] * st[2]) / 3.0f; if (pdfDist <= 0.0f) break; beta *= Tr * scene_.medium.sigmaS / pdfDist; if (beta.maxComp() <= 0.0f) break; scatterEvents++; // NEE from medium point if (S.useNEE && scene_.hasLights()) { LightSample ls; if (sampleLightDirect(pMed, Vec3(0.0f), sampler, ls, true)) { float ph = Medium::hgPhase(dot(-ray.d, ls.wi), scene_.medium.g); if (ph > 0.0f && ls.pdf > 0.0f) { Vec3 trs = transmittance(pMed, ls.wi, ls.distance, ls.primId); float w = 1.0f; if (S.useMIS && !ls.isDelta) w = powerHeuristic(1, ls.pdf, 1, ph); L += beta * ls.Li * trs * (ph * w / ls.pdf); } } } // Sample phase function float phPdf; Vec3 wi = Medium::sampleHG(-ray.d, scene_.medium.g, sampler.get1D(), sampler.get1D(), phPdf); prevBsdfPdf = phPdf; specularBounce = false; prevPos = pMed; ray = Ray(pMed, wi); // beta unchanged: phase/pdf == 1 for HG sampling // Russian roulette if (depth >= S.minDepthRR) { float q = std::min(0.95f, beta.maxComp()); if (q < S.rrThreshold || sampler.get1D() > q) break; beta /= q; } continue; } else { // Passed through: attenuate by transmittance Vec3 Tr = vexp(-st * tMaxSeg); float pdfSurv = (Tr[0] + Tr[1] + Tr[2]) / 3.0f; if (pdfSurv <= 0.0f) break; beta *= Tr / pdfSurv; } } } } // ---- Escaped: environment ------------------------------------ if (!found) { Vec3 envL = scene_.env.eval(ray.d); if (!envL.isZero()) { if (specularBounce || !S.useNEE) { L += beta * envL; } else if (S.useMIS) { int ei = scene_.envIndex(); float lightPdf = (ei >= 0) ? scene_.env.pdf(ray.d) * scene_.lightPmf(ei) : 0.0f; float w = powerHeuristic(1, prevBsdfPdf, 1, lightPdf); L += beta * envL * w; } } // Distant lights hit by chance are delta => no contribution break; } if (aux && depth == 0) { aux->normal = hit.ns; aux->depth = hit.t; aux->albedo = scene_.materials[hit.matId].albedo.eval(hit.uv, hit.p); } const Material& mat = scene_.materials[hit.matId]; Vec3 wo = -ray.d; // ---- Emitted radiance --------------------------------------- if (mat.isEmissive()) { bool front = dot(hit.ng, wo) > 0.0f; if (front || mat.emitBothSides) { Vec3 Le = mat.radiance(); if (specularBounce || !S.useNEE) { L += beta * Le; } else if (S.useMIS) { int lid = scene_.prims[hit.primId].lightId; float lightPdf = 0.0f; if (lid >= 0) { float area = scene_.prims[hit.primId].area(); float dist2 = lengthSq(hit.p - prevPos); float cosL = absDot(hit.ng, wo); if (cosL > 1e-7f && area > 0.0f) lightPdf = dist2 / (cosL * area) * scene_.lightPmf(lid); } float w = powerHeuristic(1, prevBsdfPdf, 1, lightPdf); L += beta * Le * w; } } if (mat.type == MatType::Emissive) break; // pure light, absorb } // ---- Build BSDF --------------------------------------------- BSDF bsdf(mat, hit, wo); if (!bsdf.valid()) break; // ---- Next-event estimation ---------------------------------- if (S.useNEE && !bsdf.hasDeltaOnly() && scene_.hasLights()) { int ns = std::max(1, S.lightSamplesPerBounce); Vec3 Ld(0.0f); for (int i = 0; i < ns; ++i) { LightSample ls; if (!sampleLightDirect(hit.p, hit.ns, sampler, ls, false)) continue; if (ls.pdf <= 0.0f || ls.Li.isZero()) continue; float cosSurf = dot(hit.ns, ls.wi); if (cosSurf <= 0.0f && dot(hit.ng, ls.wi) <= 0.0f) { // allow transmission through dielectrics if (mat.type != MatType::Dielectric && mat.type != MatType::ThinDielectric) continue; } Vec3 f = bsdf.eval(ls.wi); if (f.isZero()) continue; float bp = bsdf.pdf(ls.wi); Vec3 tr = transmittance(hit.p + faceforward(hit.ng, ls.wi) * (EPS * shadowBias(hit.t)), ls.wi, ls.distance, ls.primId); if (tr.isZero()) continue; float w = 1.0f; if (S.useMIS && !ls.isDelta) w = powerHeuristic(1, ls.pdf, 1, bp); Ld += f * ls.Li * tr * (std::fabs(cosSurf) * w / ls.pdf); } Ld = Ld / float(ns); if (S.clampDirect > 0.0f) { float m = Ld.maxComp(); if (m > S.clampDirect) Ld *= (S.clampDirect / m); } L += beta * Ld; } // ---- BSDF sampling ------------------------------------------ BSDFSample bs; if (!bsdf.sample(sampler, bs)) break; if (bs.pdf <= 1e-9f || bs.f.isZero()) break; float cosTerm = absDot(hit.ns, bs.wi); if (cosTerm <= 0.0f) break; Vec3 contrib = bs.f * (cosTerm / bs.pdf); if (contrib.hasNaN()) break; beta *= contrib; etaScaleAccum *= bs.etaScale; specularBounce = bs.specular; prevBsdfPdf = bs.pdf; prevPos = hit.p; // Offset origin along the correct side Vec3 offN = faceforward(hit.ng, bs.wi); Vec3 newO = hit.p + offN * (EPS * shadowBias(hit.t)); ray = Ray(newO, bs.wi); // Beer-Lambert absorption when travelling inside a dielectric if (bs.transmission && mat.absorption.maxComp() > 0.0f) { Hit exitHit; Ray probe(newO, bs.wi); if (scene_.bvh.intersect(probe, exitHit) && exitHit.matId == hit.matId) beta *= vexp(-mat.absorption * exitHit.t); } // ---- Russian roulette --------------------------------------- if (depth >= S.minDepthRR) { float rrBeta = (beta * (1.0f / std::max(etaScaleAccum, 1e-6f))).maxComp(); float q = std::min(0.95f, rrBeta); if (q < S.rrThreshold) { if (sampler.get1D() > q) break; beta /= std::max(q, 1e-6f); } else if (rrBeta < 1.0f) { if (sampler.get1D() > q) break; beta /= q; } } if (beta.maxComp() <= 0.0f) break; } if (L.hasNaN()) return Vec3(0.0f); if (S.clampIndirect > 0.0f) { float m = L.maxComp(); if (m > S.clampIndirect) L *= (S.clampIndirect / m); } return L; } // Debug shading modes Vec3 debugShade(Ray ray, Sampler& sampler) const { const std::string& mode = scene_.settings.debugMode; Hit hit; if (!scene_.bvh.intersect(ray, hit)) return (mode == "normal") ? Vec3(0.0f) : scene_.env.eval(ray.d); if (mode == "normal") return (hit.ns + Vec3(1.0f)) * 0.5f; if (mode == "albedo") return scene_.materials[hit.matId].albedo.eval(hit.uv, hit.p); if (mode == "depth") { float d = clampf(hit.t / (sceneRadius_ * 0.5f), 0.0f, 1.0f); return Vec3(1.0f - d); } if (mode == "uv") return Vec3(hit.uv.x, hit.uv.y, 0.0f); if (mode == "heat") { float m = float(hit.primId % 32) / 31.0f; return Vec3(m, 1.0f - m, 0.5f); } return Vec3(0.5f); } private: static float shadowBias(float t) { return std::max(1.0f, t * 0.01f); } Vec3 transmittance(const Vec3& p, const Vec3& wi, float dist, int skipPrim) const { float tmax = (dist == INF) ? sceneRadius_ * 4.0f : dist * (1.0f - 1e-3f); Ray sr(p, wi, EPS, tmax); bool blocked = scene_.bvh.occluded(sr, [&](int pi) { return pi == skipPrim; }); if (blocked) return Vec3(0.0f); if (scene_.medium.enabled) { Vec3 st = scene_.medium.sigmaT(); if (st.maxComp() > 0.0f) return vexp(-st * std::min(tmax, sceneRadius_ * 4.0f)); } return Vec3(1.0f); } bool sampleLightDirect(const Vec3& p, const Vec3& ns, Sampler& sampler, LightSample& out, bool inMedium) const { float pmf; int idx = scene_.selectLight(sampler.get1D(), pmf); if (idx < 0) return false; int nA = scene_.areaLightCount(); int envI = scene_.envIndex(); if (idx < nA) { const AreaLight& L = scene_.areaLights[idx]; const Primitive& pr = scene_.prims[L.primIndex]; Vec2 u = sampler.get2D(); // Cone sampling for spheres is much better than area sampling if (pr.type == PrimType::Sphere) { Vec3 toC = pr.v0 - p; float d2 = lengthSq(toC); float r2 = pr.radius * pr.radius; if (d2 <= r2 * 1.0001f) { // Inside sphere: fall back to area sampling Vec3 pos, nrm; float pdfA; pr.sampleArea(u, pos, nrm, pdfA); Vec3 d = pos - p; float dist = length(d); if (dist < 1e-6f) return false; out.wi = d / dist; float cosL = dot(nrm, -out.wi); if (cosL <= 0.0f && !L.twoSided) return false; out.pdf = pdfA * dist * dist / std::max(std::fabs(cosL), 1e-7f) * pmf; out.Li = L.radiance; out.distance = dist; out.primId = L.primIndex; out.isDelta = false; return out.pdf > 0.0f; } float dc = std::sqrt(d2); float sinThetaMax2 = r2 / d2; float cosThetaMax = std::sqrt(std::max(0.0f, 1.0f - sinThetaMax2)); Frame cf(toC / dc); Vec3 local = sampleUniformCone(u.x, u.y, cosThetaMax); Vec3 wi = cf.toWorld(local); // Compute the actual intersection distance with the sphere float b = dot(wi, -toC); float c = d2 - r2; float disc = b * b - c; float dist = (disc >= 0.0f) ? (-b - std::sqrt(disc)) : (dc - pr.radius); if (dist <= 1e-5f) dist = std::max(1e-5f, dc - pr.radius); out.wi = wi; out.pdf = (1.0f / (TWO_PI * (1.0f - cosThetaMax))) * pmf; out.Li = L.radiance; out.distance = dist; out.primId = L.primIndex; out.isDelta = false; return out.pdf > 0.0f; } Vec3 pos, nrm; float pdfA; pr.sampleArea(u, pos, nrm, pdfA); Vec3 d = pos - p; float dist2 = lengthSq(d); if (dist2 < 1e-10f) return false; float dist = std::sqrt(dist2); out.wi = d / dist; float cosL = dot(nrm, -out.wi); if (cosL <= 0.0f) { if (!L.twoSided) return false; cosL = -cosL; } out.pdf = pdfA * dist2 / std::max(cosL, 1e-7f) * pmf; out.Li = L.radiance; out.distance = dist; out.primId = L.primIndex; out.isDelta = false; return out.pdf > 0.0f; } if (envI >= 0 && idx == envI) { Vec2 u = sampler.get2D(); float pdf; Vec3 wi, Li = scene_.env.sample(u, wi, pdf); if (pdf <= 0.0f || Li.isZero()) return false; out.wi = wi; out.Li = Li; out.pdf = pdf * pmf; out.distance = INF; out.primId = -1; out.isDelta = false; return true; } // Distant / sun light int di = idx - scene_.distantOffset(); if (di >= 0 && di < int(scene_.distantLights.size())) { const DistantLight& D = scene_.distantLights[di]; Vec3 base = normalize(-D.direction); float half = radians(std::max(0.01f, D.angularRadiusDeg)); float cosMax = std::cos(half); Vec2 u = sampler.get2D(); Frame f(base); Vec3 wi = f.toWorld(sampleUniformCone(u.x, u.y, cosMax)); out.wi = wi; out.Li = D.radiance; out.pdf = (1.0f / (TWO_PI * (1.0f - cosMax))) * pmf; out.distance = INF; out.primId = -1; out.isDelta = (D.angularRadiusDeg < 0.02f); return true; } return false; } const Scene& scene_; float sceneRadius_ = 1e4f; }; // ============================================================================ // SECTION 11 — Film, tonemapping, output // ============================================================================ struct Film { int w = 0, h = 0; std::vector sum; // accumulated radiance std::vector sumSq; // for variance std::vector count; std::vector normalBuf, albedoBuf; std::vector depthBuf; void resize(int W, int H) { w = W; h = H; size_t n = size_t(W) * H; sum.assign(n, Vec3(0.0f)); sumSq.assign(n, Vec3(0.0f)); count.assign(n, 0); normalBuf.assign(n, Vec3(0.0f)); albedoBuf.assign(n, Vec3(0.0f)); depthBuf.assign(n, 0.0f); } void add(int x, int y, const Vec3& v) { size_t i = size_t(y) * w + x; sum[i] += v; sumSq[i] += v * v; count[i]++; } Vec3 mean(int x, int y) const { size_t i = size_t(y) * w + x; return count[i] ? sum[i] / float(count[i]) : Vec3(0.0f); } // Relative variance of the mean estimator (per-pixel) float relVariance(int x, int y) const { size_t i = size_t(y) * w + x; int n = count[i]; if (n < 2) return 1e9f; Vec3 m = sum[i] / float(n); Vec3 var = (sumSq[i] / float(n)) - m * m; float v = std::max(0.0f, var.average()) / float(n); float lum = std::max(1e-4f, m.luminance()); return v / (lum * lum); } }; // ---- Tonemap operators --------------------------------------------------- namespace tm { inline Vec3 reinhard(const Vec3& c) { return c / (Vec3(1.0f) + c); } inline Vec3 reinhardExtended(const Vec3& c, float whitePoint) { Vec3 num = c * (Vec3(1.0f) + c / (whitePoint * whitePoint)); return num / (Vec3(1.0f) + c); } inline Vec3 acesFilm(const Vec3& x) { // Narkowicz ACES approximation const float a = 2.51f, b = 0.03f, c = 2.43f, d = 0.59f, e = 0.14f; Vec3 num = x * (x * a + Vec3(b)); Vec3 den = x * (x * c + Vec3(d)) + Vec3(e); return clampv(num / den, 0.0f, 1.0f); } inline Vec3 uncharted2(const Vec3& x) { auto f = [](Vec3 v) { const float A = 0.15f, B = 0.50f, C = 0.10f, D = 0.20f, E = 0.02f, F = 0.30f; return ((v * (v * A + Vec3(C * B)) + Vec3(D * E)) / (v * (v * A + Vec3(B)) + Vec3(D * F))) - Vec3(E / F); }; const float W = 11.2f; Vec3 white = f(Vec3(W)); return f(x * 2.0f) / white; } inline Vec3 filmicHejl(const Vec3& x) { Vec3 c = vmax(Vec3(0.0f), x - Vec3(0.004f)); return (c * (c * 6.2f + Vec3(0.5f))) / (c * (c * 6.2f + Vec3(1.7f)) + Vec3(0.06f)); } inline Vec3 apply(const std::string& name, const Vec3& c, float gamma) { Vec3 o; if (name == "linear") o = clampv(c, 0.0f, 1.0f); else if (name == "reinhard") o = reinhard(c); else if (name == "reinhard2") o = reinhardExtended(c, 4.0f); else if (name == "aces") o = acesFilm(c); else if (name == "uncharted2") o = uncharted2(c); else if (name == "filmic") return clampv(filmicHejl(c), 0.0f, 1.0f); // has gamma baked else o = acesFilm(c); float ig = 1.0f / std::max(gamma, 1e-3f); return { std::pow(saturate(o.x), ig), std::pow(saturate(o.y), ig), std::pow(saturate(o.z), ig) }; } } // namespace tm // ---- Simple joint bilateral denoise using normal + albedo guides --------- inline void denoise(const Film& film, std::vector& colour, float sigmaSpatial) { int W = film.w, H = film.h; std::vector out(colour.size()); int radius = std::max(1, int(std::ceil(sigmaSpatial * 2.0f))); float invSS = 1.0f / (2.0f * sigmaSpatial * sigmaSpatial); const float sigmaN = 0.25f, sigmaA = 0.15f, sigmaC = 0.6f; for (int y = 0; y < H; ++y) { for (int x = 0; x < W; ++x) { size_t ci = size_t(y) * W + x; Vec3 cN = film.normalBuf[ci], cA = film.albedoBuf[ci], cC = colour[ci]; Vec3 acc(0.0f); float wsum = 0.0f; for (int dy = -radius; dy <= radius; ++dy) { int yy = y + dy; if (yy < 0 || yy >= H) continue; for (int dx = -radius; dx <= radius; ++dx) { int xx = x + dx; if (xx < 0 || xx >= W) continue; size_t ni = size_t(yy) * W + xx; float wS = std::exp(-float(dx * dx + dy * dy) * invSS); Vec3 dn = film.normalBuf[ni] - cN; float wN = std::exp(-lengthSq(dn) / (2.0f * sigmaN * sigmaN)); Vec3 da = film.albedoBuf[ni] - cA; float wA = std::exp(-lengthSq(da) / (2.0f * sigmaA * sigmaA)); Vec3 dc = colour[ni] - cC; float wC = std::exp(-lengthSq(dc) / (2.0f * sigmaC * sigmaC)); float w = wS * wN * wA * wC; acc += colour[ni] * w; wsum += w; } } out[ci] = (wsum > 0.0f) ? acc / wsum : cC; } } colour = out; } // ---- PNG writer (zlib "stored" blocks — no external deps) ---------------- namespace png { inline uint32_t crc32(const uint8_t* data, size_t len, uint32_t crc = 0) { static uint32_t table[256]; static bool init = false; if (!init) { for (uint32_t i = 0; i < 256; ++i) { uint32_t c = i; for (int k = 0; k < 8; ++k) c = (c & 1) ? 0xEDB88320u ^ (c >> 1) : (c >> 1); table[i] = c; } init = true; } crc = crc ^ 0xFFFFFFFFu; for (size_t i = 0; i < len; ++i) crc = table[(crc ^ data[i]) & 0xFF] ^ (crc >> 8); return crc ^ 0xFFFFFFFFu; } inline uint32_t adler32(const uint8_t* d, size_t n) { uint32_t a = 1, b = 0; for (size_t i = 0; i < n; ++i) { a = (a + d[i]) % 65521; b = (b + a) % 65521; } return (b << 16) | a; } inline void put32(std::vector& v, uint32_t x) { v.push_back(uint8_t(x >> 24)); v.push_back(uint8_t(x >> 16)); v.push_back(uint8_t(x >> 8)); v.push_back(uint8_t(x)); } inline void chunk(std::ofstream& f, const char tag[4], const std::vector& data) { std::vector hdr; put32(hdr, uint32_t(data.size())); f.write(reinterpret_cast(hdr.data()), hdr.size()); std::vector body(tag, tag + 4); body.insert(body.end(), data.begin(), data.end()); f.write(reinterpret_cast(body.data()), body.size()); std::vector crcv; put32(crcv, crc32(body.data(), body.size())); f.write(reinterpret_cast(crcv.data()), crcv.size()); } inline bool write(const std::string& path, int W, int H, const std::vector& rgb) { std::ofstream f(path, std::ios::binary); if (!f) return false; const uint8_t sig[8] = {0x89, 'P', 'N', 'G', 0x0D, 0x0A, 0x1A, 0x0A}; f.write(reinterpret_cast(sig), 8); std::vector ihdr; put32(ihdr, uint32_t(W)); put32(ihdr, uint32_t(H)); ihdr.push_back(8); // bit depth ihdr.push_back(2); // colour type: truecolour ihdr.push_back(0); ihdr.push_back(0); ihdr.push_back(0); chunk(f, "IHDR", ihdr); // Raw scanlines with filter byte 0 std::vector raw; raw.reserve(size_t(H) * (size_t(W) * 3 + 1)); for (int y = 0; y < H; ++y) { raw.push_back(0); const uint8_t* row = &rgb[size_t(y) * W * 3]; raw.insert(raw.end(), row, row + size_t(W) * 3); } // zlib stream with stored (uncompressed) deflate blocks std::vector z; z.push_back(0x78); z.push_back(0x01); size_t pos = 0; while (pos < raw.size()) { size_t n = std::min(65535, raw.size() - pos); bool last = (pos + n >= raw.size()); z.push_back(last ? 1 : 0); z.push_back(uint8_t(n & 0xFF)); z.push_back(uint8_t((n >> 8) & 0xFF)); uint16_t nn = uint16_t(~n); z.push_back(uint8_t(nn & 0xFF)); z.push_back(uint8_t((nn >> 8) & 0xFF)); z.insert(z.end(), raw.begin() + pos, raw.begin() + pos + n); pos += n; } put32(z, adler32(raw.data(), raw.size())); chunk(f, "IDAT", z); chunk(f, "IEND", {}); return true; } } // namespace png inline bool writePPM(const std::string& path, int W, int H, const std::vector& rgb) { std::ofstream f(path, std::ios::binary); if (!f) return false; f << "P6\n" << W << " " << H << "\n255\n"; f.write(reinterpret_cast(rgb.data()), rgb.size()); return true; } // ---- Radiance .hdr loader (RLE + flat) ----------------------------------- inline std::shared_ptr loadHDR(const std::string& path) { std::ifstream f(path, std::ios::binary); if (!f) { std::cerr << "[hdr] cannot open " << path << "\n"; return nullptr; } std::string line; if (!std::getline(f, line)) return nullptr; if (line.rfind("#?", 0) != 0) { std::cerr << "[hdr] bad magic\n"; return nullptr; } // header while (std::getline(f, line)) { if (line.empty()) break; } if (!std::getline(f, line)) return nullptr; int W = 0, H = 0; if (std::sscanf(line.c_str(), "-Y %d +X %d", &H, &W) != 2) { std::cerr << "[hdr] unsupported resolution line: " << line << "\n"; return nullptr; } auto img = std::make_shared(); img->w = W; img->h = H; img->px.assign(size_t(W) * H, Vec3(0.0f)); std::vector scan(size_t(W) * 4); auto rgbeToVec = [](uint8_t r, uint8_t g, uint8_t b, uint8_t e) -> Vec3 { if (e == 0) return Vec3(0.0f); float s = std::ldexp(1.0f, int(e) - 128 - 8); return Vec3(float(r) * s, float(g) * s, float(b) * s); }; for (int y = 0; y < H; ++y) { uint8_t hdr4[4]; f.read(reinterpret_cast(hdr4), 4); if (!f) break; if (hdr4[0] == 2 && hdr4[1] == 2 && ((hdr4[2] << 8) | hdr4[3]) == W && W >= 8 && W < 32768) { // New RLE for (int c = 0; c < 4; ++c) { int x = 0; while (x < W) { uint8_t cnt; f.read(reinterpret_cast(&cnt), 1); if (!f) break; if (cnt > 128) { uint8_t val; f.read(reinterpret_cast(&val), 1); int n = cnt - 128; for (int i = 0; i < n && x < W; ++i, ++x) scan[size_t(x) * 4 + c] = val; } else { int n = cnt; for (int i = 0; i < n && x < W; ++i, ++x) { uint8_t val; f.read(reinterpret_cast(&val), 1); scan[size_t(x) * 4 + c] = val; } } } } for (int x = 0; x < W; ++x) img->px[size_t(y) * W + x] = rgbeToVec(scan[size_t(x)*4+0], scan[size_t(x)*4+1], scan[size_t(x)*4+2], scan[size_t(x)*4+3]); } else { // Flat img->px[size_t(y) * W + 0] = rgbeToVec(hdr4[0], hdr4[1], hdr4[2], hdr4[3]); for (int x = 1; x < W; ++x) { uint8_t p[4]; f.read(reinterpret_cast(p), 4); img->px[size_t(y) * W + x] = rgbeToVec(p[0], p[1], p[2], p[3]); } } } std::cout << "[hdr] loaded " << path << " (" << W << "x" << H << ")\n"; return img; } // ============================================================================ // SECTION 12 — OBJ loader // ============================================================================ inline bool loadOBJ(const std::string& path, Scene& scene, int matId, const Mat4& xform, bool smoothNormals, bool flipNormals) { std::ifstream f(path); if (!f) { std::cerr << "[obj] cannot open " << path << "\n"; return false; } std::vector V, N; std::vector T; Mat4 nm = xform.normalMatrix(); struct FaceVert { int v = -1, t = -1, n = -1; }; std::vector> tris; std::string line; while (std::getline(f, line)) { if (line.empty() || line[0] == '#') continue; std::istringstream ss(line); std::string tag; ss >> tag; if (tag == "v") { float x, y, z; ss >> x >> y >> z; V.push_back(xform.xformPoint({x, y, z})); } else if (tag == "vn") { float x, y, z; ss >> x >> y >> z; N.push_back(normalize(nm.xformVector({x, y, z}))); } else if (tag == "vt") { float u, v = 0; ss >> u >> v; T.push_back({u, v}); } else if (tag == "f") { std::vector fv; std::string tok; while (ss >> tok) { FaceVert q; // formats: v | v/t | v//n | v/t/n size_t s1 = tok.find('/'); if (s1 == std::string::npos) q.v = std::stoi(tok) - 1; else { q.v = std::stoi(tok.substr(0, s1)) - 1; size_t s2 = tok.find('/', s1 + 1); if (s2 == std::string::npos) q.t = std::stoi(tok.substr(s1 + 1)) - 1; else { if (s2 > s1 + 1) q.t = std::stoi(tok.substr(s1 + 1, s2 - s1 - 1)) - 1; if (s2 + 1 < tok.size()) q.n = std::stoi(tok.substr(s2 + 1)) - 1; } } if (q.v < 0) q.v += int(V.size()) + 1; fv.push_back(q); } for (size_t i = 2; i < fv.size(); ++i) tris.push_back({fv[0], fv[i - 1], fv[i]}); } } // Compute smooth normals if requested and none provided std::vector autoN; if (smoothNormals && N.empty()) { autoN.assign(V.size(), Vec3(0.0f)); for (auto& tr : tris) { const Vec3& a = V[tr[0].v]; const Vec3& b = V[tr[1].v]; const Vec3& c = V[tr[2].v]; Vec3 fn = cross(b - a, c - a); autoN[tr[0].v] += fn; autoN[tr[1].v] += fn; autoN[tr[2].v] += fn; } for (auto& n : autoN) { float l = length(n); if (l > 1e-12f) n /= l; } } size_t before = scene.prims.size(); for (auto& tr : tris) { if (tr[0].v < 0 || tr[1].v < 0 || tr[2].v < 0) continue; if (size_t(tr[0].v) >= V.size() || size_t(tr[1].v) >= V.size() || size_t(tr[2].v) >= V.size()) continue; Primitive p; p.type = PrimType::Triangle; p.matId = matId; p.v0 = V[tr[0].v]; p.v1 = V[tr[1].v]; p.v2 = V[tr[2].v]; p.flipNormal = flipNormals; if (!N.empty() && tr[0].n >= 0 && tr[1].n >= 0 && tr[2].n >= 0 && size_t(tr[0].n) < N.size() && size_t(tr[1].n) < N.size() && size_t(tr[2].n) < N.size()) { p.n0 = N[tr[0].n]; p.n1 = N[tr[1].n]; p.n2 = N[tr[2].n]; p.hasVertexNormals = true; } else if (!autoN.empty()) { p.n0 = autoN[tr[0].v]; p.n1 = autoN[tr[1].v]; p.n2 = autoN[tr[2].v]; p.hasVertexNormals = true; } if (!T.empty()) { auto getT = [&](int i) { return (i >= 0 && size_t(i) < T.size()) ? T[i] : Vec2(0, 0); }; p.t0 = getT(tr[0].t); p.t1 = getT(tr[1].t); p.t2 = getT(tr[2].t); } scene.prims.push_back(p); } std::cout << "[obj] " << path << ": " << (scene.prims.size() - before) << " triangles\n"; return true; } // ============================================================================ // SECTION 13 — INI configuration parser & scene builder // ============================================================================ class Config { public: // sections -> key -> value (repeated sections get numeric suffix) struct Section { std::string name; std::map kv; std::string get(const std::string& k, const std::string& def = "") const { auto it = kv.find(k); return (it == kv.end()) ? def : it->second; } float getf(const std::string& k, float def) const { auto it = kv.find(k); if (it == kv.end()) return def; try { return std::stof(it->second); } catch (...) { return def; } } int geti(const std::string& k, int def) const { auto it = kv.find(k); if (it == kv.end()) return def; try { return std::stoi(it->second); } catch (...) { return def; } } bool getb(const std::string& k, bool def) const { auto it = kv.find(k); if (it == kv.end()) return def; std::string v = it->second; std::transform(v.begin(), v.end(), v.begin(), ::tolower); return v == "1" || v == "true" || v == "yes" || v == "on"; } Vec3 getv(const std::string& k, const Vec3& def) const { auto it = kv.find(k); if (it == kv.end()) return def; std::istringstream ss(it->second); float a, b, c; std::string s = it->second; for (char& ch : s) if (ch == ',') ch = ' '; std::istringstream s2(s); if (s2 >> a) { if (s2 >> b) { if (s2 >> c) return {a, b, c}; return {a, b, b}; } return Vec3(a); } return def; } }; std::vector

sections; bool load(const std::string& path) { std::ifstream f(path); if (!f) { std::cerr << "[cfg] cannot open " << path << "\n"; return false; } Section cur; cur.name = "global"; std::string line; while (std::getline(f, line)) { // strip comments size_t hp = line.find_first_of(";#"); if (hp != std::string::npos) line = line.substr(0, hp); trim(line); if (line.empty()) continue; if (line.front() == '[' && line.back() == ']') { if (!cur.kv.empty() || cur.name != "global") sections.push_back(cur); cur = Section(); cur.name = line.substr(1, line.size() - 2); trim(cur.name); std::transform(cur.name.begin(), cur.name.end(), cur.name.begin(), ::tolower); continue; } size_t eq = line.find('='); if (eq == std::string::npos) continue; std::string k = line.substr(0, eq), v = line.substr(eq + 1); trim(k); trim(v); std::transform(k.begin(), k.end(), k.begin(), ::tolower); cur.kv[k] = v; } sections.push_back(cur); return true; } const Section* first(const std::string& name) const { for (const auto& s : sections) if (s.name == name) return &s; return nullptr; } std::vector all(const std::string& name) const { std::vector out; for (const auto& s : sections) if (s.name == name) out.push_back(&s); return out; } private: static void trim(std::string& s) { size_t a = s.find_first_not_of(" \t\r\n"); size_t b = s.find_last_not_of(" \t\r\n"); s = (a == std::string::npos) ? "" : s.substr(a, b - a + 1); } }; inline MatType parseMatType(const std::string& sIn) { std::string s = sIn; std::transform(s.begin(), s.end(), s.begin(), ::tolower); if (s == "diffuse" || s == "lambert") return MatType::Diffuse; if (s == "conductor" || s == "metal") return MatType::Conductor; if (s == "dielectric" || s == "glass") return MatType::Dielectric; if (s == "thin" || s == "thindielectric") return MatType::ThinDielectric; if (s == "plastic" || s == "coated") return MatType::Plastic; if (s == "emissive" || s == "light") return MatType::Emissive; if (s == "mirror") return MatType::Mirror; return MatType::Diffuse; } inline TexKind parseTexKind(const std::string& sIn) { std::string s = sIn; std::transform(s.begin(), s.end(), s.begin(), ::tolower); if (s == "checker") return TexKind::Checker; if (s == "gradient") return TexKind::Gradient; if (s == "image" || s == "texture") return TexKind::Image; if (s == "noise") return TexKind::Noise; return TexKind::Constant; } inline SamplerType parseSampler(const std::string& sIn) { std::string s = sIn; std::transform(s.begin(), s.end(), s.begin(), ::tolower); if (s == "random") return SamplerType::Random; if (s == "stratified") return SamplerType::Stratified; return SamplerType::Sobol; } inline Mat4 parseTransform(const Config::Section& s) { Mat4 M; Vec3 t = s.getv("translate", Vec3(0.0f)); Vec3 sc = s.getv("scale", Vec3(1.0f)); Vec3 ax = s.getv("rotate_axis", Vec3(0, 1, 0)); float ang = s.getf("rotate_angle", 0.0f); M = Mat4::translate(t) * Mat4::rotate(ax, ang) * Mat4::scale(sc); return M; } inline void buildSceneFromConfig(const Config& cfg, Scene& scene) { // ---------------- Render settings ------------------------------------- RenderSettings& S = scene.settings; if (auto* r = cfg.first("render")) { S.width = r->geti("width", S.width); S.height = r->geti("height", S.height); S.samplesPerPixel = r->geti("spp", S.samplesPerPixel); S.maxDepth = r->geti("max_depth", S.maxDepth); S.minDepthRR = r->geti("rr_start_depth", S.minDepthRR); S.rrThreshold = r->getf("rr_threshold", S.rrThreshold); S.useNEE = r->getb("next_event_estimation", S.useNEE); S.useMIS = r->getb("multiple_importance_sampling", S.useMIS); S.lightSamplesPerBounce = r->geti("light_samples", S.lightSamplesPerBounce); S.clampIndirect = r->getf("clamp_indirect", S.clampIndirect); S.clampDirect = r->getf("clamp_direct", S.clampDirect); S.tileSize = r->geti("tile_size", S.tileSize); S.threads = r->geti("threads", S.threads); S.sampler = parseSampler(r->get("sampler", "sobol")); S.seed = uint32_t(r->geti("seed", int(S.seed))); S.adaptive = r->getb("adaptive", S.adaptive); S.adaptiveMinSpp = r->geti("adaptive_min_spp", S.adaptiveMinSpp); S.adaptiveTolerance = r->getf("adaptive_tolerance", S.adaptiveTolerance); S.adaptiveBatch = r->geti("adaptive_batch", S.adaptiveBatch); S.debugMode = r->get("debug", S.debugMode); S.progressBar = r->getb("progress", S.progressBar); S.previewEvery = r->geti("preview_every", S.previewEvery); } if (auto* o = cfg.first("output")) { S.outputFile = o->get("file", S.outputFile); S.tonemap = o->get("tonemap", S.tonemap); std::transform(S.tonemap.begin(), S.tonemap.end(), S.tonemap.begin(), ::tolower); S.exposure = o->getf("exposure", S.exposure); S.gamma = o->getf("gamma", S.gamma); S.writeAux = o->getb("write_aux", S.writeAux); S.denoiseSimple = o->getb("denoise", S.denoiseSimple); S.denoiseSigma = o->getf("denoise_sigma", S.denoiseSigma); } // ---------------- Camera ---------------------------------------------- if (auto* c = cfg.first("camera")) { Camera& cam = scene.camera; cam.position = c->getv("position", cam.position); cam.lookAt = c->getv("look_at", cam.lookAt); cam.up = c->getv("up", cam.up); cam.fovDeg = c->getf("fov", cam.fovDeg); cam.aperture = c->getf("aperture", cam.aperture); cam.focusDistance = c->getf("focus_distance", cam.focusDistance); cam.bokehBlades = c->geti("bokeh_blades", cam.bokehBlades); cam.bokehRotationDeg = c->getf("bokeh_rotation", cam.bokehRotationDeg); } // ---------------- Environment ----------------------------------------- if (auto* e = cfg.first("environment")) { std::string m = e->get("mode", "gradient"); std::transform(m.begin(), m.end(), m.begin(), ::tolower); if (m == "none") scene.env.mode = EnvMode::None; else if (m == "constant") scene.env.mode = EnvMode::Constant; else if (m == "hdri" || m == "image") scene.env.mode = EnvMode::HDRI; else scene.env.mode = EnvMode::Gradient; scene.env.constColour = e->getv("colour", e->getv("color", scene.env.constColour)); scene.env.horizonColour = e->getv("horizon", scene.env.horizonColour); scene.env.zenithColour = e->getv("zenith", scene.env.zenithColour); scene.env.groundColour = e->getv("ground", scene.env.groundColour); scene.env.intensity = e->getf("intensity", scene.env.intensity); scene.env.rotationDeg = e->getf("rotation", scene.env.rotationDeg); scene.env.importanceSample = e->getb("importance_sample", true); std::string file = e->get("file", ""); if (!file.empty()) { scene.env.hdri = loadHDR(file); if (scene.env.hdri) scene.env.mode = EnvMode::HDRI; else if (scene.env.mode == EnvMode::HDRI) scene.env.mode = EnvMode::Gradient; } } // ---------------- Medium ---------------------------------------------- if (auto* m = cfg.first("medium")) { scene.medium.enabled = m->getb("enabled", false); scene.medium.sigmaA = m->getv("sigma_a", Vec3(0.0f)); scene.medium.sigmaS = m->getv("sigma_s", Vec3(0.0f)); scene.medium.g = m->getf("g", 0.0f); scene.medium.maxScatterEvents = m->geti("max_scatter_events", 8); } // ---------------- Materials ------------------------------------------- for (auto* ms : cfg.all("material")) { Material mt; mt.name = ms->get("name", "mat" + std::to_string(scene.materials.size())); mt.type = parseMatType(ms->get("type", "diffuse")); mt.albedo.kind = parseTexKind(ms->get("texture", "constant")); mt.albedo.c0 = ms->getv("albedo", ms->getv("colour", ms->getv("color", Vec3(0.8f)))); mt.albedo.c1 = ms->getv("albedo2", Vec3(0.2f)); mt.albedo.scaleU = ms->getf("tex_scale_u", ms->getf("tex_scale", 8.0f)); mt.albedo.scaleV = ms->getf("tex_scale_v", mt.albedo.scaleU); mt.albedo.noiseSeed = uint32_t(ms->geti("noise_seed", 12345)); std::string tf = ms->get("texture_file", ""); if (!tf.empty()) { mt.albedo.img = loadHDR(tf); if (mt.albedo.img) mt.albedo.kind = TexKind::Image; } mt.emission = ms->getv("emission", Vec3(0.0f)); mt.emissionScale = ms->getf("emission_scale", 1.0f); mt.emitBothSides = ms->getb("two_sided", false); mt.roughness = ms->getf("roughness", 0.1f); mt.anisotropy = ms->getf("anisotropy", 0.0f); mt.ior = ms->getf("ior", 1.5f); mt.specularTint = ms->getv("specular_tint", Vec3(1.0f)); mt.absorption = ms->getv("absorption", Vec3(0.0f)); mt.k = ms->getv("extinction_k", Vec3(0.0f)); mt.useComplexIOR = ms->getb("complex_ior", false); if (mt.type == MatType::Emissive && mt.emission.isZero()) mt.emission = Vec3(1.0f); scene.addMaterial(mt); } if (scene.materials.empty()) { Material def; def.name = "default"; def.albedo.c0 = Vec3(0.75f); scene.addMaterial(def); } auto resolveMat = [&](const Config::Section* s) -> int { std::string n = s->get("material", ""); int id = n.empty() ? -1 : scene.findMaterial(n); if (id < 0) { // allow inline material via index id = s->geti("material_id", 0); id = std::max(0, std::min(int(scene.materials.size()) - 1, id)); } return id; }; // ---------------- Primitives ------------------------------------------ for (auto* s : cfg.all("sphere")) { Primitive p; p.type = PrimType::Sphere; p.v0 = s->getv("center", Vec3(0.0f)); p.radius = s->getf("radius", 1.0f); p.matId = resolveMat(s); p.flipNormal = s->getb("flip_normal", false); scene.prims.push_back(p); } for (auto* s : cfg.all("quad")) { Primitive p; p.type = PrimType::Quad; p.v0 = s->getv("origin", Vec3(0.0f)); p.v1 = s->getv("edge1", Vec3(1, 0, 0)); p.v2 = s->getv("edge2", Vec3(0, 0, 1)); p.matId = resolveMat(s); p.flipNormal = s->getb("flip_normal", false); scene.prims.push_back(p); } for (auto* s : cfg.all("disk")) { Primitive p; p.type = PrimType::Disk; p.v0 = s->getv("center", Vec3(0.0f)); p.n0 = s->getv("normal", Vec3(0, 1, 0)); p.radius = s->getf("radius", 1.0f); p.matId = resolveMat(s); p.flipNormal = s->getb("flip_normal", false); scene.prims.push_back(p); } for (auto* s : cfg.all("triangle")) { Primitive p; p.type = PrimType::Triangle; p.v0 = s->getv("v0", Vec3(0.0f)); p.v1 = s->getv("v1", Vec3(1, 0, 0)); p.v2 = s->getv("v2", Vec3(0, 1, 0)); p.matId = resolveMat(s); p.flipNormal = s->getb("flip_normal", false); scene.prims.push_back(p); } for (auto* s : cfg.all("box")) { Vec3 lo = s->getv("min", Vec3(-0.5f)); Vec3 hi = s->getv("max", Vec3(0.5f)); int mid = resolveMat(s); bool flip = s->getb("flip_normal", false); Vec3 d = hi - lo; auto addQuad = [&](Vec3 o, Vec3 e1, Vec3 e2) { Primitive p; p.type = PrimType::Quad; p.v0 = o; p.v1 = e1; p.v2 = e2; p.matId = mid; p.flipNormal = flip; scene.prims.push_back(p); }; addQuad(lo, Vec3(d.x,0,0), Vec3(0,d.y,0)); // -Z addQuad(Vec3(lo.x,lo.y,hi.z), Vec3(0,d.y,0), Vec3(d.x,0,0)); // +Z addQuad(lo, Vec3(0,d.y,0), Vec3(0,0,d.z)); // -X addQuad(Vec3(hi.x,lo.y,lo.z), Vec3(0,0,d.z), Vec3(0,d.y,0)); // +X addQuad(lo, Vec3(0,0,d.z), Vec3(d.x,0,0)); // -Y addQuad(Vec3(lo.x,hi.y,lo.z), Vec3(d.x,0,0), Vec3(0,0,d.z)); // +Y } for (auto* s : cfg.all("mesh")) { std::string file = s->get("file", ""); if (file.empty()) continue; loadOBJ(file, scene, resolveMat(s), parseTransform(*s), s->getb("smooth", true), s->getb("flip_normal", false)); } // ---------------- Lights ---------------------------------------------- for (auto* s : cfg.all("distant_light")) { DistantLight D; D.direction = normalize(s->getv("direction", Vec3(0, -1, 0))); D.radiance = s->getv("radiance", Vec3(1.0f)) * s->getf("intensity", 1.0f); D.angularRadiusDeg = s->getf("angular_radius", 0.53f); scene.distantLights.push_back(D); } for (auto* s : cfg.all("sun")) { // alias DistantLight D; Vec3 dir = s->getv("direction", Vec3(0, 0, 0)); if (dir.isZero()) { float el = radians(s->getf("elevation", 45.0f)); float az = radians(s->getf("azimuth", 30.0f)); Vec3 toSun(std::cos(el) * std::cos(az), std::sin(el), std::cos(el) * std::sin(az)); dir = -toSun; } D.direction = normalize(dir); D.radiance = s->getv("radiance", Vec3(1.0f, 0.95f, 0.85f)) * s->getf("intensity", 20.0f); D.angularRadiusDeg = s->getf("angular_radius", 0.53f); scene.distantLights.push_back(D); } } // ---- Built-in default scene (used when no config supplied) --------------- inline void buildDefaultScene(Scene& scene) { scene.settings.width = 900; scene.settings.height = 600; scene.settings.samplesPerPixel = 256; scene.settings.maxDepth = 20; scene.settings.outputFile = "render.png"; scene.settings.tonemap = "aces"; scene.camera.position = {0.0f, 1.6f, 6.5f}; scene.camera.lookAt = {0.0f, 1.0f, 0.0f}; scene.camera.fovDeg = 38.0f; scene.camera.aperture = 0.08f; scene.camera.bokehBlades = 6; scene.env.mode = EnvMode::Gradient; scene.env.intensity = 0.7f; Material floor; floor.name = "floor"; floor.type = MatType::Diffuse; floor.albedo.kind = TexKind::Checker; floor.albedo.c0 = Vec3(0.82f); floor.albedo.c1 = Vec3(0.15f); floor.albedo.scaleU = floor.albedo.scaleV = 12.0f; int mFloor = scene.addMaterial(floor); Material gold; gold.name = "gold"; gold.type = MatType::Conductor; gold.albedo.c0 = Vec3(1.0f, 0.78f, 0.34f); gold.roughness = 0.12f; int mGold = scene.addMaterial(gold); Material glass; glass.name = "glass"; glass.type = MatType::Dielectric; glass.albedo.c0 = Vec3(1.0f); glass.roughness = 0.0f; glass.ior = 1.52f; glass.absorption = Vec3(0.15f, 0.05f, 0.02f); int mGlass = scene.addMaterial(glass); Material plastic; plastic.name = "plastic"; plastic.type = MatType::Plastic; plastic.albedo.c0 = Vec3(0.15f, 0.35f, 0.85f); plastic.roughness = 0.25f; int mPlastic = scene.addMaterial(plastic); Material brushed; brushed.name = "brushed"; brushed.type = MatType::Conductor; brushed.albedo.c0 = Vec3(0.92f, 0.92f, 0.95f); brushed.roughness = 0.3f; brushed.anisotropy = 0.8f; int mBrushed = scene.addMaterial(brushed); Material lightM; lightM.name = "light"; lightM.type = MatType::Emissive; lightM.emission = Vec3(1.0f, 0.92f, 0.8f); lightM.emissionScale = 18.0f; int mLight = scene.addMaterial(lightM); // Ground Primitive ground; ground.type = PrimType::Quad; ground.v0 = {-20, 0, -20}; ground.v1 = {40, 0, 0}; ground.v2 = {0, 0, 40}; ground.matId = mFloor; scene.prims.push_back(ground); auto sphere = [&](Vec3 c, float r, int m) { Primitive p; p.type = PrimType::Sphere; p.v0 = c; p.radius = r; p.matId = m; scene.prims.push_back(p); }; sphere({-2.2f, 1.0f, 0.0f}, 1.0f, mGold); sphere({ 0.0f, 1.0f, 0.0f}, 1.0f, mGlass); sphere({ 2.2f, 1.0f, 0.0f}, 1.0f, mPlastic); sphere({ 0.0f, 0.45f, 2.3f}, 0.45f, mBrushed); // Area light (quad above) Primitive lq; lq.type = PrimType::Quad; lq.v0 = {-1.5f, 4.5f, -1.5f}; lq.v1 = {3.0f, 0, 0}; lq.v2 = {0, 0, 3.0f}; lq.matId = mLight; lq.flipNormal = true; scene.prims.push_back(lq); DistantLight sun; sun.direction = normalize(Vec3(-0.4f, -0.8f, -0.45f)); sun.radiance = Vec3(1.0f, 0.95f, 0.85f) * 3.0f; sun.angularRadiusDeg = 1.0f; scene.distantLights.push_back(sun); } // ============================================================================ // SECTION 14 — Renderer (multithreaded tiles + adaptive sampling) // ============================================================================ class Renderer { public: explicit Renderer(Scene& sc) : scene_(sc), integrator_(sc) {} void render() { const RenderSettings& S = scene_.settings; int W = S.width, H = S.height; film_.resize(W, H); int nThreads = (S.threads > 0) ? S.threads : std::max(1u, std::thread::hardware_concurrency()); // Build tile list ordered from centre outward (nicer previews) struct Tile { int x0, y0, x1, y1; }; std::vector tiles; int ts = std::max(8, S.tileSize); for (int y = 0; y < H; y += ts) for (int x = 0; x < W; x += ts) tiles.push_back({x, y, std::min(x + ts, W), std::min(y + ts, H)}); float cx = W * 0.5f, cy = H * 0.5f; std::sort(tiles.begin(), tiles.end(), [&](const Tile& a, const Tile& b) { float da = sqr((a.x0 + a.x1) * 0.5f - cx) + sqr((a.y0 + a.y1) * 0.5f - cy); float db = sqr((b.x0 + b.x1) * 0.5f - cx) + sqr((b.y0 + b.y1) * 0.5f - cy); return da < db; }); bool debug = (S.debugMode != "none" && !S.debugMode.empty()); int totalSpp = debug ? 1 : S.samplesPerPixel; std::cout << "=====================================================\n"; std::cout << " Resolution : " << W << " x " << H << "\n"; std::cout << " Samples : " << totalSpp << " spp" << (S.adaptive ? " (adaptive)" : "") << "\n"; std::cout << " Max depth : " << S.maxDepth << "\n"; std::cout << " Primitives : " << scene_.prims.size() << "\n"; std::cout << " BVH nodes : " << scene_.bvh.nodeCount() << " depth " << scene_.bvh.depth() << "\n"; std::cout << " Lights : " << scene_.numLights() << " (area " << scene_.areaLightCount() << ", distant " << scene_.distantLights.size() << ", env " << (scene_.env.active() ? 1 : 0) << ")\n"; std::cout << " Threads : " << nThreads << "\n"; std::cout << " Integrator : path tracer" << (S.useNEE ? " + NEE" : "") << (S.useMIS ? " + MIS" : "") << "\n"; if (debug) std::cout << " DEBUG mode : " << S.debugMode << "\n"; std::cout << "=====================================================\n"; auto t0 = std::chrono::steady_clock::now(); std::atomic nextTile{0}; std::atomic pixelsDone{0}; std::atomic raysTraced{0}; long long totalPixels = (long long)W * H; auto worker = [&](int tid) { Sampler sampler; sampler.init(S.sampler, S.seed + uint32_t(tid) * 7919u, totalSpp); while (true) { size_t ti = nextTile.fetch_add(1); if (ti >= tiles.size()) break; const Tile& T = tiles[ti]; for (int y = T.y0; y < T.y1; ++y) { for (int x = T.x0; x < T.x1; ++x) { renderPixel(x, y, sampler, totalSpp, debug, raysTraced); pixelsDone.fetch_add(1); } } } }; std::thread progressThread; std::atomic running{true}; if (S.progressBar) { progressThread = std::thread([&]() { while (running.load()) { long long p = pixelsDone.load(); float frac = float(p) / float(totalPixels); auto now = std::chrono::steady_clock::now(); float el = std::chrono::duration(now - t0).count(); float eta = (frac > 1e-3f) ? el * (1.0f / frac - 1.0f) : 0.0f; int bars = int(frac * 40.0f); std::printf("\r ["); for (int i = 0; i < 40; ++i) std::printf(i < bars ? "=" : " "); std::printf("] %5.1f%% %6.1fs elapsed, ETA %6.1fs ", frac * 100.0f, el, eta); std::fflush(stdout); std::this_thread::sleep_for(std::chrono::milliseconds(250)); } }); } std::vector pool; for (int i = 0; i < nThreads; ++i) pool.emplace_back(worker, i); for (auto& t : pool) t.join(); running = false; if (progressThread.joinable()) progressThread.join(); auto t1 = std::chrono::steady_clock::now(); double secs = std::chrono::duration(t1 - t0).count(); std::printf("\r [========================================] 100.0%% %.2fs total \n", secs); long long rays = raysTraced.load(); std::printf(" Camera rays: %lld (%.2f Mrays/s)\n", rays, double(rays) / secs / 1e6); writeOutput(); } private: void renderPixel(int x, int y, Sampler& sampler, int totalSpp, bool debug, std::atomic& raysTraced) { const RenderSettings& S = scene_.settings; size_t pi = size_t(y) * S.width + x; if (debug) { sampler.startPixel(x, y, 0); Ray r = scene_.camera.generateRay(float(x) + 0.5f, float(y) + 0.5f, sampler); film_.add(x, y, integrator_.debugShade(r, sampler)); raysTraced.fetch_add(1); return; } int done = 0; int batch = S.adaptive ? std::max(1, S.adaptiveBatch) : totalSpp; while (done < totalSpp) { int n = std::min(batch, totalSpp - done); for (int s = 0; s < n; ++s) { sampler.startPixel(x, y, done + s); float jx = sampler.get1D(), jy = sampler.get1D(); Ray ray = scene_.camera.generateRay(float(x) + jx, float(y) + jy, sampler); Integrator::AuxData aux; Vec3 L = integrator_.radiance(ray, sampler, &aux); film_.add(x, y, L); if (done + s == 0) { film_.normalBuf[pi] = aux.normal; film_.albedoBuf[pi] = aux.albedo; film_.depthBuf[pi] = aux.depth; } else { // running average of guides for better denoise quality float w = 1.0f / float(done + s + 1); film_.normalBuf[pi] = lerp(film_.normalBuf[pi], aux.normal, w); film_.albedoBuf[pi] = lerp(film_.albedoBuf[pi], aux.albedo, w); } } done += n; raysTraced.fetch_add(n); if (S.adaptive && done >= S.adaptiveMinSpp) { if (film_.relVariance(x, y) < S.adaptiveTolerance) break; } } } void writeOutput() { const RenderSettings& S = scene_.settings; int W = S.width, H = S.height; std::vector hdr(size_t(W) * H); float expMul = std::pow(2.0f, S.exposure); for (int y = 0; y < H; ++y) for (int x = 0; x < W; ++x) hdr[size_t(y) * W + x] = film_.mean(x, y) * expMul; if (S.denoiseSimple) { std::cout << " Denoising...\n"; denoise(film_, hdr, S.denoiseSigma); } std::vector rgb(size_t(W) * H * 3); bool isDebug = (S.debugMode != "none" && !S.debugMode.empty()); for (size_t i = 0; i < hdr.size(); ++i) { Vec3 c = isDebug ? clampv(hdr[i], 0.0f, 1.0f) : tm::apply(S.tonemap, hdr[i], S.gamma); rgb[i * 3 + 0] = uint8_t(clampf(c.x, 0.0f, 1.0f) * 255.0f + 0.5f); rgb[i * 3 + 1] = uint8_t(clampf(c.y, 0.0f, 1.0f) * 255.0f + 0.5f); rgb[i * 3 + 2] = uint8_t(clampf(c.z, 0.0f, 1.0f) * 255.0f + 0.5f); } std::string out = S.outputFile; bool ok; if (out.size() > 4 && out.substr(out.size() - 4) == ".ppm") ok = writePPM(out, W, H, rgb); else { if (out.size() < 4 || out.substr(out.size() - 4) != ".png") out += ".png"; ok = png::write(out, W, H, rgb); } std::cout << (ok ? " Wrote " : " FAILED writing ") << out << "\n"; if (S.writeAux) { auto dumpVec = [&](const std::vector& buf, const std::string& name, bool remap) { std::vector b(size_t(W) * H * 3); for (size_t i = 0; i < buf.size(); ++i) { Vec3 v = remap ? (buf[i] + Vec3(1.0f)) * 0.5f : buf[i]; v = clampv(v, 0.0f, 1.0f); b[i*3+0] = uint8_t(v.x * 255.0f); b[i*3+1] = uint8_t(v.y * 255.0f); b[i*3+2] = uint8_t(v.z * 255.0f); } png::write(name, W, H, b); std::cout << " Wrote " << name << "\n"; }; dumpVec(film_.normalBuf, "aux_normal.png", true); dumpVec(film_.albedoBuf, "aux_albedo.png", false); float maxD = 1e-6f; for (float d : film_.depthBuf) if (std::isfinite(d)) maxD = std::max(maxD, d); std::vector db(size_t(W) * H * 3); for (size_t i = 0; i < film_.depthBuf.size(); ++i) { float v = clampf(1.0f - film_.depthBuf[i] / maxD, 0.0f, 1.0f); uint8_t g = uint8_t(v * 255.0f); db[i*3+0] = db[i*3+1] = db[i*3+2] = g; } png::write("aux_depth.png", W, H, db); std::cout << " Wrote aux_depth.png\n"; // Variance / sample-count heat map std::vector vb(size_t(W) * H * 3); int maxC = 1; for (int c : film_.count) maxC = std::max(maxC, c); for (int y = 0; y < H; ++y) for (int x = 0; x < W; ++x) { size_t i = size_t(y) * W + x; float t = float(film_.count[i]) / float(maxC); Vec3 c = lerp(Vec3(0.1f, 0.2f, 0.8f), Vec3(1.0f, 0.3f, 0.1f), t); vb[i*3+0] = uint8_t(c.x * 255.0f); vb[i*3+1] = uint8_t(c.y * 255.0f); vb[i*3+2] = uint8_t(c.z * 255.0f); } png::write("aux_samples.png", W, H, vb); std::cout << " Wrote aux_samples.png\n"; } } Scene& scene_; Integrator integrator_; Film film_; }; } // namespace ptm // ============================================================================ // SECTION 15 — main // ============================================================================ int main(int argc, char** argv) { using namespace ptm; Scene scene; if (argc > 1) { Config cfg; if (!cfg.load(argv[1])) { std::cerr << "Falling back to default scene.\n"; buildDefaultScene(scene); } else { buildSceneFromConfig(cfg, scene); if (scene.prims.empty()) { std::cerr << "[warn] config produced no geometry; using default scene.\n"; scene = Scene(); buildDefaultScene(scene); } } } else { std::cout << "No config given — rendering built-in demo scene.\n" << "Usage: " << argv[0] << " scene.ini\n"; buildDefaultScene(scene); } // CLI overrides: key=value after the config path for (int i = 2; i < argc; ++i) { std::string a = argv[i]; size_t eq = a.find('='); if (eq == std::string::npos) continue; std::string k = a.substr(0, eq), v = a.substr(eq + 1); if (k == "spp") scene.settings.samplesPerPixel = std::stoi(v); else if (k == "width") scene.settings.width = std::stoi(v); else if (k == "height") scene.settings.height = std::stoi(v); else if (k == "depth") scene.settings.maxDepth = std::stoi(v); else if (k == "threads")scene.settings.threads = std::stoi(v); else if (k == "out") scene.settings.outputFile = v; else if (k == "debug") scene.settings.debugMode = v; else if (k == "tonemap")scene.settings.tonemap = v; else if (k == "exposure") scene.settings.exposure = std::stof(v); else std::cerr << "[cli] unknown override: " << k << "\n"; } scene.finalize(); Renderer renderer(scene); renderer.render(); return 0; } ``` --- ## `config.ini` ```ini ; ============================================================================ ; Advanced Path Tracer — example scene configuration ; Run: ./pt config.ini ; Override on CLI: ./pt config.ini spp=1024 width=1920 height=1080 ; ============================================================================ [render] width = 1000 height = 640 spp = 512 max_depth = 24 rr_start_depth = 4 rr_threshold = 0.05 next_event_estimation = true multiple_importance_sampling = true light_samples = 1 ; Firefly suppression (0 = off) clamp_indirect = 0 clamp_direct = 0 sampler = sobol ; random | stratified | sobol seed = 20240517 tile_size = 32 threads = 0 ; 0 = all cores ; Adaptive sampling: stop early where the estimate has converged adaptive = true adaptive_min_spp = 64 adaptive_batch = 32 adaptive_tolerance = 0.002 debug = none ; none | normal | albedo | depth | uv | heat progress = true [output] file = render.png ; .png or .ppm tonemap = aces ; linear|reinhard|reinhard2|aces|uncharted2|filmic exposure = 0.0 ; in EV stops gamma = 2.2 write_aux = true ; dump normal/albedo/depth/sample-count buffers denoise = false denoise_sigma = 1.2 [camera] position = 0.0, 1.7, 7.0 look_at = 0.0, 1.0, 0.0 up = 0, 1, 0 fov = 36 aperture = 0.10 ; 0 = pinhole focus_distance = 7.0 ; -1 = auto (distance to look_at) bokeh_blades = 6 ; 0 = circular bokeh_rotation = 15 [environment] mode = gradient ; none | constant | gradient | hdri intensity = 0.6 horizon = 0.85, 0.90, 1.00 zenith = 0.25, 0.45, 0.95 ground = 0.20, 0.18, 0.16 rotation = 0 importance_sample = true ; file = studio.hdr ; uncomment for HDRI lighting [medium] enabled = false sigma_a = 0.002 sigma_s = 0.010 g = 0.3 max_scatter_events = 6 ; ---------------------------------------------------------------- materials -- [material] name = checker_floor type = diffuse texture = checker albedo = 0.85, 0.85, 0.82 albedo2 = 0.08, 0.08, 0.09 tex_scale = 14 [material] name = gold type = conductor albedo = 1.00, 0.77, 0.34 roughness = 0.10 [material] name = brushed_steel type = conductor albedo = 0.93, 0.94, 0.96 roughness = 0.30 anisotropy = 0.85 [material] name = clear_glass type = dielectric albedo = 1, 1, 1 roughness = 0.0 ior = 1.52 absorption = 0.12, 0.04, 0.02 [material] name = frosted_glass type = dielectric albedo = 1, 1, 1 roughness = 0.22 ior = 1.5 [material] name = blue_plastic type = plastic albedo = 0.10, 0.28, 0.80 roughness = 0.20 ior = 1.5 [material] name = red_matte type = diffuse albedo = 0.72, 0.10, 0.08 roughness = 0.6 ; > 0 enables Oren-Nayar [material] name = ceiling_light type = emissive emission = 1.0, 0.93, 0.82 emission_scale = 16.0 [material] name = accent_light type = emissive emission = 0.25, 0.55, 1.0 emission_scale = 30.0 [material] name = mirror type = mirror albedo = 0.96, 0.96, 0.96 ; ----------------------------------------------------------------- geometry -- [quad] material = checker_floor origin = -25, 0, -25 edge1 = 50, 0, 0 edge2 = 0, 0, 50 [sphere] material = gold center = -2.3, 1.0, 0.0 radius = 1.0 [sphere] material = clear_glass center = 0.0, 1.0, 0.0 radius = 1.0 [sphere] material = blue_plastic center = 2.3, 1.0, 0.0 radius = 1.0 [sphere] material = brushed_steel center = -1.1, 0.5, 2.2 radius = 0.5 [sphere] material = frosted_glass center = 1.1, 0.5, 2.2 radius = 0.5 [sphere] material = red_matte center = 0.0, 0.35, 3.4 radius = 0.35 [box] material = mirror min = -5.0, 0.0, -4.2 max = 5.0, 3.2, -4.0 ; --------------------------------------------------------------------- lights [quad] material = ceiling_light origin = -1.8, 5.0, -1.8 edge1 = 3.6, 0, 0 edge2 = 0, 0, 3.6 flip_normal = true [sphere] material = accent_light center = 3.6, 0.30, 2.6 radius = 0.14 [sun] elevation = 42 azimuth = 205 intensity = 3.0 radiance = 1.0, 0.96, 0.88 angular_radius = 1.0 ``` --- ## `Makefile` ```makefile CXX ?= g++ CXXFLAGS := -std=c++17 -O3 -march=native -ffast-math -fno-finite-math-only \ -Wall -Wextra -Wno-unused-parameter -pthread TARGET := pt all: $(TARGET) $(TARGET): main.cpp $(CXX) $(CXXFLAGS) $< -o $@ debug: main.cpp $(CXX) -std=c++17 -O0 -g -fsanitize=address,undefined -pthread $< -o $(TARGET)_dbg run: $(TARGET) ./$(TARGET) config.ini preview: $(TARGET) ./$(TARGET) config.ini spp=16 width=480 height=320 out=preview.png hq: $(TARGET) ./$(TARGET) config.ini spp=4096 width=1920 height=1200 out=hq.png clean: rm -f $(TARGET) $(TARGET)_dbg *.png *.ppm .PHONY: all debug run preview hq clean ``` --- ## Overview of the implementation **Integrator** — Unidirectional path tracing. Each bounce does next-event estimation against a power-weighted light CDF (area lights, distant/sun lights, and the environment as a single "light"), then BSDF importance sampling. The two strategies are combined with the power heuristic, so the returned-emitter path and the light-sampled path don't double-count. Delta lobes skip MIS entirely (`bs.specular` flag). **BSDFs** — All microfacet lobes use anisotropic GGX with visible-normal-distribution sampling (Heitz 2018) and height-correlated Smith masking-shadowing, so energy behaviour is correct at grazing angles. `Diffuse` upgrades to Oren-Nayar when `roughness > 0`. `Dielectric` handles both the smooth delta case and the rough case with the proper refraction Jacobian and `1/η²` radiance scaling. Beer-Lambert absorption fires when a transmission event enters a medium-bearing material. **Acceleration** — Binned SAH BVH, near-child-first traversal using ray direction sign, explicit 64-entry stack, separate `occluded()` path for shadow rays with early exit. **Sampling** — `Sampler` wraps three strategies. The Sobol path uses Owen scrambling via the Laine-Karras hash on both the sample index and the per-dimension output, which decorrelates neighbouring pixels and gives blue-noise-like error distribution at low spp. Dimensions above 4 fall back to PCG32. **Adaptive sampling** — Per-pixel `sumSq` gives the variance of the mean; pixels stop once relative variance drops below `adaptive_tolerance`. The `aux_samples.png` output visualises where budget went. **Notable settings worth trying** | Setting | Effect | |---|---| | `debug=normal\|albedo\|depth\|uv` | instant single-sample diagnostic renders | | `[medium] enabled=true` | global fog with Henyey-Greenstein phase + NEE from scatter points | | `anisotropy=0.85` | brushed-metal streaking | | `bokeh_blades=6` | hexagonal out-of-focus highlights | | `sampler=random` vs `sobol` | compare noise character at 16 spp | | `clamp_indirect=10` | kill fireflies from caustic paths | **Practical caveats.** The PNG writer emits stored (uncompressed) deflate blocks — valid PNG, but files are large; swap in libpng or `stb_image_write` if size matters. The denoiser is a joint bilateral filter, useful for previews but not competitive with OIDN. The medium is homogeneous and global rather than bounded per-object. Spectral rendering, bidirectional/VCM integrators, and subsurface scattering are the natural next additions.