SLProject  4.3.020
A platform independent 3D computer graphics framework for desktop OS, Android, iOS and online in web browsers
SLPathtracer.cpp File Reference
#include <algorithm>
#include <cmath>
#include <typeinfo>
#include <SLCamera.h>
#include <SLLightRect.h>
#include <SLPathtracer.h>
#include <SLSceneView.h>
#include <GlobalTimer.h>
#include <Profiler.h>
Include dependency graph for SLPathtracer.cpp:

Go to the source code of this file.

Functions

SLfloat rnd01 ()
 
static SLfloat russianRoulette (SLint depth, SLCol4f albedo)
 Russian roulette survival test for the continuation of a path. More...
 
static SLfloat phongLobeWeight (SLfloat exponent, const SLVec3f &sampleDir, const SLVec3f &normal)
 Weight of one sample drawn from the Phong lobe of the given exponent. More...
 
static SLfloat misWeight (SLfloat pdfThis, SLfloat pdfOther)
 Power heuristic (beta = 2) weight of the strategy whose density is pdfThis. More...
 
static SLfloat lightPdfMC (SLRay *ray)
 Solid angle density with which shade() would have sampled this light hit. More...
 

Variables

static const SLfloat PDF_NO_MIS = -1.0f
 Sentinel pdf for SLPathtracer::trace. More...
 
static const SLint RR_START_DEPTH = 3
 Path depth up to which Russian roulette always lets the path continue. More...
 

Detailed Description

Date
July 2014
Authors
Thomas Schneiter, Marcus Hudritsch
Remarks
Please use clangformat to format the code. See more code style on https://github.com/cpvrlab/SLProject4/wiki/SLProject-Coding-Style

Definition in file SLPathtracer.cpp.

Function Documentation

◆ lightPdfMC()

static SLfloat lightPdfMC ( SLRay ray)
static

Solid angle density with which shade() would have sampled this light hit.

Returns 0 if the hit surface is not a light that shade() samples over its area. shade() draws a point uniformly on the rectangle, so the density over the area is 1/area; the conversion to a density over solid angle, as seen from the shaded point, is the usual dist^2 / (cosLight * area).

Definition at line 164 of file SLPathtracer.cpp.

165 {
166  if (!ray->hitNode || typeid(*ray->hitNode) != typeid(SLLightRect))
167  return 0.0f;
168 
169  auto* rect = (SLLightRect*)ray->hitNode;
170  if (!rect->isOn()) return 0.0f;
171 
172  // SLRay::setDir does not normalise, so do not assume a unit direction
173  SLVec3f dir(ray->dir);
174  dir.normalize();
175 
176  // spotDirWS is the normal of the rectangle, as in shade()
177  SLfloat cosLight = -dir.dot(rect->spotDirWS());
178  if (cosLight <= 0.0f) return 0.0f;
179 
180  return (ray->length * ray->length) / (cosLight * rect->area());
181 }
float SLfloat
analog to GLfloat
Definition: SL.h:200
Light node class for a rectangular light source.
Definition: SLLightRect.h:39
SLVec3f dir
Direction vector of ray in WS.
Definition: SLRay.h:79
SLfloat length
length from origin to an intersection
Definition: SLRay.h:80
SLNode * hitNode
Points to the intersected node.
Definition: SLRay.h:108

◆ misWeight()

static SLfloat misWeight ( SLfloat  pdfThis,
SLfloat  pdfOther 
)
static

Power heuristic (beta = 2) weight of the strategy whose density is pdfThis.

Two strategies that can both generate the same path each return an unbiased estimate of it, so simply adding them would count the path twice. Weighting them with w(pdfThis) + w(pdfOther) = 1 counts it exactly once, and the power heuristic puts nearly all of that weight on whichever strategy had the higher density for this particular path, which is the one with the lower variance here. Written as 1/(1+r^2) rather than a^2/(a^2+b^2) so that a very large density cannot overflow.

Definition at line 151 of file SLPathtracer.cpp.

152 {
153  if (pdfThis <= 0.0f) return 0.0f;
154 
155  SLfloat ratio = pdfOther / pdfThis;
156  return 1.0f / (1.0f + ratio * ratio);
157 }

◆ phongLobeWeight()

static SLfloat phongLobeWeight ( SLfloat  exponent,
const SLVec3f sampleDir,
const SLVec3f normal 
)
static

Weight of one sample drawn from the Phong lobe of the given exponent.

SLRay::reflectMC and SLRay::refractMC draw a direction from the Phong lobe

p(w) = (n + 1) / (2 * PI) * cos^n(alpha)

where alpha is the angle to the perfect specular or transmissive direction and n is the shininess or the translucency. The matching normalised Phong BSDF is

f(w) = rho * (n + 2) / (2 * PI) * cos^n(alpha)

and the Monte Carlo estimator of the rendering equation weights the incoming radiance by f * cos(theta) / p, with theta the angle to the SURFACE NORMAL and not to the lobe axis. The cos^n(alpha) and the 2*PI cancel and what is left is

rho * (n + 2) / (n + 1) * cos(theta)

The colour rho is applied by the caller as the material's specular or transmissive colour, so this returns the scalar part.

The cos(theta) is what was missing: the code applied only (n+2)/(n+1), so every glossy sample was too bright by 1/cos(theta). That is a factor of 1 straight along the normal and unbounded at the horizon, which is the direction a wide lobe on a grazing surface samples most often, so the error showed up as a bright rim exactly where a glossy highlight is supposed to fall off.

For a perfect mirror the caller applies no weight at all rather than calling this: (n+2)/(n+1) is the normalisation of the lobe estimator and there is no lobe to normalise when the direction was not drawn from one. At the PERFECT limit of 1000 that factor is 1.001, so it was silently adding a tenth of a percent of energy per specular bounce.

Parameters
exponentshininess for reflection, translucency for transmission
sampleDirsampled direction, already known to be on the right side
normalthe surface normal at the hit point

Definition at line 130 of file SLPathtracer.cpp.

133 {
134  // The absolute value, because SLMesh::preShade does not flip the hit
135  // normal towards the ray: on a back face hit both the normal and the
136  // sampled direction sit on the other side and the cosine comes out
137  // negative although the geometry is the same.
138  SLfloat cosTheta = std::abs(sampleDir.dot(normal));
139 
140  return (exponent + 2.0f) / (exponent + 1.0f) * cosTheta;
141 }
T dot(const SLVec3 &v) const
Definition: SLVec3.h:117
T abs(T a)
Definition: Utils.h:249

◆ rnd01()

SLfloat rnd01 ( )

Uniform random number generator for numbers between 0 and 1 that is used in SLRay, SLLightRect and SLPathtracer.

Remarks
The engine state is thread_local and must stay that way. The ray tracer and the path tracer call rnd01 concurrently from all worker threads (see SLPathtracer::render). A single shared std::mt19937 would be a data race on its 624 word state plus its position index. That is undefined behaviour, and in practice the racing threads hand each other torn and repeated values, so their samples are no longer independent and the noise no longer averages out with 1/sqrt(N). Each thread seeds its own engine from the current time mixed with a shared atomic counter, so that threads created within the same second still get different sequences.

Definition at line 50 of file SLRay.cpp.

51 {
52  static std::atomic<SLuint> seedCounter{0};
53 
54  thread_local std::mt19937 engine((SLuint)std::time(nullptr) * 2654435761u +
55  seedCounter.fetch_add(1u) * 40503u + 1u);
56 
57  thread_local std::uniform_real_distribution<SLfloat> dist(0.0f, 1.0f);
58 
59  return dist(engine);
60 }
unsigned int SLuint
analog to GLuint
Definition: SL.h:198

◆ russianRoulette()

static SLfloat russianRoulette ( SLint  depth,
SLCol4f  albedo 
)
static

Russian roulette survival test for the continuation of a path.

Returns 0 if the path is absorbed and must not be continued, otherwise the survival probability that the caller has to divide its contribution by.

Russian roulette is not the same thing as Monte Carlo, although both are random. Monte Carlo is the estimator itself: sample from a density and average f(x)/p(x). Russian roulette is a technique used inside it, and it answers a different question — how to terminate an unbounded recursion without making the answer wrong. Light bounces between the walls forever, so the recursion has to be stopped somewhere. Stopping it at a fixed depth discards the light that the longer paths would have carried, which is a bias that no number of samples can remove. Stopping it at random and dividing the survivors by their survival probability leaves the mean untouched:

E = q * 0 + (1 - q) * L / (1 - q) = L

The survival probability is the albedo of the surface, so that the 1/survival of a survivor cancels the albedo of that bounce exactly. That is the classic formulation: the colour of the surface decides how likely the path is to continue, rather than how much it is dimmed. A surface with an albedo of 1 absorbs nothing and therefore always survives, which is why the hard depth cap in trace() is still needed as a safety net.

Definition at line 84 of file SLPathtracer.cpp.

85 {
86  if (depth <= RR_START_DEPTH) return 1.0f;
87 
88  SLfloat survival = std::min(albedo.maxXYZ(), 1.0f);
89 
90  if (survival <= 0.0f) return 0.0f; // black surface, nothing to continue
91  if (survival >= 1.0f) return 1.0f; // absorbs nothing, always continues
92 
93  return (rnd01() < survival) ? survival : 0.0f;
94 }
static const SLint RR_START_DEPTH
Path depth up to which Russian roulette always lets the path continue.
SLfloat rnd01()
Definition: SLRay.cpp:50
T maxXYZ()
Definition: SLVec4.h:141

Variable Documentation

◆ PDF_NO_MIS

const SLfloat PDF_NO_MIS = -1.0f
static

Sentinel pdf for SLPathtracer::trace.

Handed down for rays that next event estimation cannot generate: the primary ray, and every specular or transmissive bounce. A light source reached by such a ray contributes its full emission, because shade() never produced a competing estimate of the same path.

Definition at line 29 of file SLPathtracer.cpp.

◆ RR_START_DEPTH

const SLint RR_START_DEPTH = 3
static

Path depth up to which Russian roulette always lets the path continue.

Purely a variance and cost knob: every value here is unbiased, so this trades render time against noise and nothing else.

Roulette is what makes the estimator unbiased, but it is not free. Because the survival probability is the albedo, the 1/survival of a survivor cancels that bounce's attenuation exactly, so a path that survives fifteen bounces arrives carrying full weight instead of 0.75^15. Roulette therefore turns "many paths, each tiny" into "few paths, each full size": the same mean with a much heavier tail. Taking the first bounces deterministically keeps that tail out of the part of the sum that carries most of the energy.

A model of the interreflection series alone (albedo 0.75, exact answer 1/(1-0.75) = 4) says the roulette tail shrinks quickly with this value:

start depth   mean     sd     mean path length
          3   3.998   1.46     7.00
          5   4.000   0.82     9.02
          8   4.000   0.34    12.00

Measured on the real scene, it does not. Raising it from 3 to 8 made the ceiling around the light worse (median absolute residual 3.77 -> 4.05 at 100 spp) for 1.7 times the path length. The model is wrong about what dominates there: the noise in that region is not the roulette tail but the near field of the area light, where a surface centimetres from a 1.0 x 0.65 emitter receives an enormous radiance that neighbouring surfaces then see through a rare bounce. Tracing deeper only finds more of those. So the value is kept low: roulette removes the depth bias, and the near field needs a different fix (solid angle sampling of the rectangle, or the sample clamp).

Definition at line 60 of file SLPathtracer.cpp.