Volume and Participating Media
Digital Image Synthesis g g y Yung-Yu Chuang
with slides by Pat Hanrahan and Torsten Moller
Participating media
• We have by far assumed that the scene is in a vacuum Hence radiance is constant along the vacuum. Hence, radiance is constant along the ray. However, some real-world situations such as fog and smoke attenuate and scatter light as fog and smoke attenuate and scatter light.
They participate in rendering.
N t l h
• Natural phenomena
– Fog, smoke, fire – Atmosphere haze – Beam of light through
clouds
– Subsurface scattering
Volume scattering processes
• Absorption (conversion from light to other forms)
E i i
• Emission (contribution from luminous particles)
• Scattering (direction change of particles) – Out-scattering
– In-scattering
– Single scattering v.s. multiple scattering
• Homogeneous v.s. inhomogeneous(heterogeneous) Homogeneous v.s. inhomogeneous(heterogeneous)
emission in-scattering out-scattering absorption
Single scattering and multiple scattering
attenuation attenuation
single scattering multiple scattering
Absorption
The reduction of energy due to conversion of light to another form of energy (e g heat)
to another form of energy (e.g. heat)
L dL
( )
L x ( , ) L dL
L x
ds ) , (
ax
ds
ds x L x x
dL ( x ) ( x ) L ( x ) ds dL ( , )
a( , ) ( , )
Absorption cross-section:
a( x , )
Probability of being absorbed per unit length
Transmittance
ds x L x x
dL ( , )
a( , ) ( , )
dL ( )
xsdL ( ' )
s
sds x x
L x dL
a
( , )
) , (
) ,
(
s
a s
x
x
ds s
x x L
x dL
0
' ) , ' ) (
, ' (
) , '
(
s L x
ax s ds s
x L
0
) ( ' ) , ' ( )
, ( ln ) , (
ln
x Optical distance or depth
x
p p 0
' s
x s '
s ax s ds s
0
' ) , ' ( )
(
Homogenous media: constant ( ) s s
a
s
0
a
( ) s
as
Transmittance and opacity
ds x L x x
dL ( , )
a( , ) ( , )
dL ( )
xsdL ( ' )
s
sds x x
L x dL
a
( , )
) , (
) ,
(
s
a s
x
x
ds s
x x L
x dL
0
' ) , ' ) (
, ' (
) , '
(
) ( L
s L x
ax s ds s
x L
0
) ( ' ) , ' ( )
, ( ln ) , (
ln
) ( ) ( ) ( )
(
( )L T L
L ( x s , ) e
(s)L ( x , ) T
( s ) L ( x , ) L ( x , )
L
s
Transmittance T (s )
)
)
(( s e
sT
Opacity
) , ( x s L )
( s e T
) ( 1 )
( T L ( x s , )
) ( 1 )
( s T
s
Absorption
Emission
• Energy that is added to the environment from luminous particles due to chemical thermal or luminous particles due to chemical, thermal, or nuclear processes that convert energy to visible light
light.
• : emitted radiance added to a ray per it di t t i t i di ti
) , ( x L
veunit distance at a point x in direction
ds x L x
dL ( x ) L ( x ) ds dL ( , )
ve( , )
L dL ( , )
L x ( , ) L dL
L x
ds ) , ( x L
ve ds
Emission
Out-scattering
Light heading in one direction is scattered to other directions due to collisions with particles
directions due to collisions with particles
L dL ( , )
L x
) , (
sx ds
d L dL ( ) ( ) ( ) Scattering cross-section:
ds x L x x
dL ( , )
s( , ) ( , ) Scattering cross-section:
Probability of being scattered per unit length
sy g p g
Extinction
L dL ( , )
L x ( , ) x L dL
ds ) , (
tx ds
ds x L x x
dL ( , )
t( , ) ( , )
t a s
Total cross-section
Alb d
s st a s
W
Albedo
Attenuation due to both absorption and scattering
s
t
x s ds
s ) ( ' , ) '
(
s
tx s ds
0
) , (
)
(
Extinction
• Beam transmittance x '
s
x ) , ( x
L
s
t x s ds
e x x
Tr
0' ) , ' (
) ' (
s: distance between x and x’ x
• Properties of Tr:
• In vaccum Tr ( x x ' ) 1
• In vaccum
• Multiplicative
B ’ l (i h di )
1 ) ( x x Tr
) '' ' ( ) ' (
) ''
( x x Tr x x Tr x x
Tr
• Beer’s law (in homogeneous medium)
ts
e x x
Tr ( x )' x e
Tr ( )
In-scattering
Increased radiance due to scattering from other directions
directions
L dL
( )
L x
'
L dL ( , )
L x
ds ) , (
sx ds
ds d x L x
p x x
dL
s
( , ) ( , ' ) ( , ' ) '
) ,
(
Phase function p ( )
s
p
) , ( ) ,
( ) , ( ) , (
( ) ( )
p p Reciprocity
( ) 1
p d
Energy conserving Phase function p ( )
( ) ( )
p p
2
( ) 1
S
p d
Source term
' ) ' , ( ) ' , ( ) , ( ) , ( ) ,
( x L x x p x L x d
S
ve s
ds x S x
dL ( , ) ( , )
• S is determined by
Volume emission – Volume emission
– Phase function which describes the angular
distribution of scattered radiation (volume analog of distribution of scattered radiation (volume analog of BSDF for surfaces)
Phase functions
Phase angle cos Phase functions
(from the phase of the moon)
( p )
1. Isotropic 1
(cos ) p 4
- simple
2. Rayleigh
2 43 1 cos (cos )
p 4
- Molecules (useful for very small particles whose radii smaller than wavelength of light)
3 Mi tt i 3. Mie scattering
- small spheres (based on Maxwell’s equations; good
model for scattering in the atmosphere due to water
model for scattering in the atmosphere due to water
droplets and fog)
Henyey-Greenstein phase function
Empirical phase function
1 1
2p p
g = -0.3
2
2 32
1 1
(cos )
4 1 2 cos
p g
g g
g = 0.6
0
2 p (cos ) cos d g
0
g : average phase angle
Henyey-Greenstein approximation
• Any phase function can be written in terms of a series of Legendre polynomials (typically n<4) series of Legendre polynomials (typically, n<4)
1 ( 2 1 ) (cos ) )
(cos n b P
p
0
) (cos )
1 2 4 (
) (cos
n
n n
P b n
p
1 ) (
P
) ( )
( P
b
) (
1 ) (
1 0
x x P
x P
cos ) (cos ) (cos
) (cos ), (cos
1
d P
p
P p
b
n n
) 1 3 2 ( ) 1
(
22
x x
P
) (cos ) cos (cos
1
d P
p
n
) 3 5 2 ( ) 1
(
33
x x x
P
...
Schlick approximation
• Approximation to Henyey-Greenstein
2 2
) cos 1 (
1 4 ) 1
(cos
k
p
Schlick k
• K plays a similar role like g
) cos 1 (
4 k
• K plays a similar role like g
– 0: isotropic
-1: back scattering – -1: back scattering
– Could use k 1 . 55 g 0 . 55 g
2Importance sampling for HG
1 1
2(cos ) g
p
2
32(cos )
4 1 2 cos
p
g g
2
2 if 0
1
2
g
otherwise 2
1 1 1
| 2
|
cos 1
2 2 2
g g g g
g
1 2
| 2
| g g g
Pbrt implementation
• core/volume.* volume/*
class VolumeRegion { t0 t1
public:
bool IntersectP(Ray &ray, float *t0, float *t1);
Spectrum sigma a(Point &, Vector &);p g _ ( , ) Spectrum sigma_s(Point &, Vector &);
Spectrum Lve(Point &, Vector &);
// phase functions: pbrt has isotropic, Rayleigh, // phase functions: pbrt has isotropic, Rayleigh, // Mie, HG, Schlick
virtual float p(Point &, Vector &, Vector &);
// attenuation coefficient; s a+s s // attenuation coefficient; s_a+s_s Spectrum sigma_t(Point &, Vector &);
// calculate optical thickness by Monte Carlo or
// l d f l ti
// closed-form solution
Spectrum Tau(Ray &ray, float step=1., float offset=0.5);
t0 step t1
}; t0 step t1
offset
Homogenous volume
• Determined by (constant)
d
– and
– g in phase function
E i i L
s
a– Emission – Spatial extent
L
veSpectrum Tau(Ray &ray, float, float){
float t0 t1;
float t0, t1;
if (!IntersectP(ray,&t0,&t1)) return 0.;
return Distance(ray(t0),ray(t1)) * (sig_a + sig_s);
}
Homogenous volume Varying-density volumes
• Density is varying in the medium and the volume scattering properties at a point is the product of scattering properties at a point is the product of the density at that point and some baseline value
value.
• DensityRegion
– 3D grid, VolumeGrid
– Exponential density, ExponentialDensity
DensityRegion
class DensityRegion : public VolumeRegion { public:
DensityRegion(Spectrum &sig_a, Spectrum &sig_s, float g, Spectrum &Le, Transform &VolumeToWorld);
float Density(Point &Pobj) const = 0;
sigma_a(Point &p, Vector &) {
return Density(WorldToVolume(p)) * sig_a; } Spectrum sigma s(Point &p, Vector &) {p g _ ( p, ) {
return Density(WorldToVolume(p)) * sig_s; } Spectrum sigma_t(Point &p, Vector &) {
return Density(WorldToVolume(p))*(sig a+sig s);}
return Density(WorldToVolume(p)) (sig_a+sig_s);}
Spectrum Lve(Point &p, Vector &) {
return Density(WorldToVolume(p)) * le; } ...
protected:
Transform WorldToVolume;
S i i l
Spectrum sig_a, sig_s, le;
float g;
};
VolumeGrid
• Standard form of given data
T i li i l i f d i
• Tri-linear interpolation of data to give continuous volume
• Often used in volume rendering
VolumeGrid
VolumeGrid(Spectrum &sa, Spectrum &ss, float gg, Spectrum &emit, BBox &e, Transform &v2w,p , , , int nx, int ny, int nz, const float *d);
float VolumeGrid::Density(const Point &Pobj) const { if (!extent.Inside(Pobj)) return 0;
// Compute voxel coordinates and offsets // Compute voxel coordinates and offsets float voxx = (Pobj.x - extent.pMin.x) /
(extent.pMax.x - extent.pMin.x) * nx - .5f;
( p p )
float voxy = (Pobj.y - extent.pMin.y) /
(extent.pMax.y - extent.pMin.y) * ny - .5f;
float voxz = (Pobj.z - extent.pMin.z) /
(extent.pMax.z - extent.pMin.z) * nz - .5f;
VolumeGrid
int vx = Floor2Int(voxx);
int vy = Floor2Int(voxy);
int vz = Floor2Int(voxz);
float dx = voxx - vx, dy = voxy - vy, dz = voxz - vz;
// Trilinearly interpolate density values // Trilinearly interpolate density values
float d00 = Lerp(dx, D(vx, vy, vz), D(vx+1, vy, vz));
float d10 = Lerp(dx, D(vx, vy+1, vz), D(vx+1, vy+1, vz));
float d01 = Lerp(dx, D(vx, vy, vz+1), D(vx+1, vy, vz+1));
float d11 = Lerp(dx, D(vx, vy+1,vz+1),D(vx+1,vy+1,vz+1));
float d0 = Lerp(dy, d00, d10);p( y, , );
float d1 = Lerp(dy, d01, d11);
return Lerp(dz, d0, d1);
}
float D(int x, int y, int z) {
} x = Clamp(x, 0, nx-1);
y = Clamp(y, 0, ny-1);
z = Clamp(z, 0, nz-1);
d i [ * * * ]
return density[z*nx*ny+y*nx+x];
}
Exponential density
• Given by
h
bhd ( )
• Where h is the height
ae
bhh
d ( )
ExponentialDensity in the direction of
the up-vector
ExponentialDensity
class ExponentialDensity : public DensityRegion { public:
ExponentialDensity(Spectrum &sa, Spectrum &ss, float g, Spectrum &emit, BBox &e, Transform &v2w, float aa, float bb, Vector &up)
...
float Density(const Point &Pobj) const {y( j) { if (!extent.Inside(Pobj)) return 0;
float height = Dot(Pobj - extent.pMin, upDir);
return a * expf(-b * height);
return a expf( b height);
} private:
BBox extent; P bj
h BBox extent;
float a, b;
Vector upDir;
}
upDir Pobj
};
extent.pMin
Light transport
• Emission + in-scattering (source term)
' ) ' , ( ) ' , ( ) , ( ) , ( ) ,
( x L x x p x L x d
S
ve s
• Absorption + out-scattering (extinction)
ds x S x
dL ( , ) ( , )
• Absorption + out scattering (extinction)
ds x L x x
dL ( , )
t( , ) ( , )
• Combined
) , ( ) , ( ) , ) (
,
(
x S x L ds x
x dL
t
Infinite length, no surface
• Assume that there is no surface and we have an infinite length we have the solution
infinite length, we have the solution
ds x S x x Tr x
L ( ( , , ) )
( ( ' ) ) ( ( ' , , ) )
0
s t x s ds
T
0' ) , ' (
) ' (
) ' ( S
x e
x
Tr ( ' )
0x '
x L ( x , )
) , ' ( x S
) ' ( x x
Tr L ( x , )
' x s
x '
With surface
• The solution
ds x
S x x Tr x
L x x Tr x
L
d
) , ' ( ) ' ( ) , ( ) (
) , (
0 0
0
from the surface point x
0) , ( x
0 L
x
) , ( x L
0
x
d ) ( x
0x
Tr
) , ( x L
With surface
• The solution
ds x
S x x Tr x
L x x Tr x
L
d
) , ' ( ) ' ( ) , ( ) (
) , (
0 0
0
from the surface point x
0from the participating media )
, ( x
0 L
x
) , ( x L
0
x
d ) ( x
0x
Tr
'
x L ( x , )
' x s x '
Simple atmosphere model
Assumptions
H di
• Homogenous media
• Constant source term (airlight)
( ) ( )
L s L S
( )
t
( )
s L s S
s
( ) 1
s sL ( ) 1
ts S
tsC S C
L s e
S e
C S C
• Fog
• Haze
• Haze
OpenGL fog model
fog
in
f C
fC
C ( 1 )
) (d
GL_EXP
)
)
(( z e
densityzf
GL EXP2
)2
)
(( z e
densityzf
_
z
f end
) (
GL_LINEAR
start z end
f ( )
From http://wiki.delphigl.com/index.php/glFog
VolumeIntegrator
class VolumeIntegrator : public Integrator { public: Beam transmittance for a given public:
virtual Spectrum Transmittance(
const Scene *scene
g ray from mint to maxt const Scene *scene,
const Renderer *renderer, const RayDifferential &ray const RayDifferential &ray,
const Sample *sample, …) const = 0;
} };
Pick up functions Preprocess(), RequestSamples()
d f
and Li() from Integrator.
Emission only
• Solution for the emission-only simplification
) , ' ( ) , '
( x L x
S
evd
ds x
L x x Tr x
L x x Tr x
L
evd
) , ' ( ) ' ( ) , ( ) (
) , (
0 0
0
• Monte Carlo estimator
N
i
i ev i
N
i
i ev
i
Tr x x L x
N t t x
p x L x x Tr
N
10 1 1
) , ( ) ) (
(
) , ( ) (
1
i
i
p x
iN
N
1( )
1Emission only
• Use multiplicativity of Tr
) (
) (
)
( x x Tr x x
1Tr x
1x Tr
i
i
i
i
• Break up integral and compute it incrementally by ray marching y y g
• Tr can get small in a long ray
– Early ray termination – Early ray termination
– Either use Russian Roulette or deterministically
EmissionIntegrator
class EmissionIntegrator : public VolumeIntegrator { public:
EmissionIntegrator(float ss) { stepSize = ss; } void RequestSamples(Sampler *sampler, Sample
*sample Scene *scene);
*sample, Scene *scene);
Spectrum Li(Scene *scene, Renderer *renderer, RayDifferential &ray, Sample *sample, RNG &rng, Spectrum *transmittance, MemoryArena &arena);
Spectrum Transmittance(Scene *scene, Renderer *, RayDifferential &ray, Sample *sample, RNG &rng, RayDifferential &ray, Sample sample, RNG &rng, MemoryArena &arena);
private:
float stepSize;
int tauSampleOffset, scatterSampleOffset;
};
single 1D sample for each
};
single 1D sample for each
EmissionIntegrator::Transmittance
if (!scene->volumeRegion) return Spectrum(1);
float step, offset;
float step, offset;
if (sample) {
step = stepSize;
smaple is non-NULL only for camera rays
p p ;
offset = sample->oneD[tauSampleOffset][0];
} else {
l t f h d d
step = 4.f * stepSize;
offset = rng.RandomFloat();
use larger steps for shadow and indirect rays for efficiency
}
Spectrum tau = scene->volumeRegion->tau(ray, step, offset);
return Exp(-tau);
)
)
((
(s)
( s )
s
a( x s ' , ) ds '
)
( s e
sT
0EmissionIntegrator::Li
VolumeRegion *vr = scene->volumeRegion;
float t0, t1;
if (!vr || !vr->IntersectP(ray, &t0, &t1)
|| (t1-t0) == 0.f) {
*T = Spectrum(1 f);
*T = Spectrum(1.f);
return 0.f;
}
// Do emission-only volume integration in vr Spectrum Lv(0.);
// Prepare for volume integration stepping // Prepare for volume integration stepping int nSamples = Ceil2Int((t1-t0) / stepSize);
float step = (t1 - t0) / nSamples;
Spectrum Tr(1.f);
Point p = ray(t0), pPrev;
Vector w = -ray.d;y ;
t0 += sample->oneD[scatterSampleOffset][0] * step;
EmissionIntegrator::Li
for (int i = 0; i < nSamples; ++i, t0 += step) { pPrev = p; p = ray(t0);
Ray tauRay(pPrev, p - pPrev, 0.f, 1.f, ray.time, ray.depth);
Spectrum stepTau = vr->tau(tauRay,.5f * stepSize, rng.RandomFloat());
Tr *= Exp(-stepTau);
// Possibly terminate if transmittance is small
) (
) (
)
( x x Tr x x
1Tr x
1x Tr
i
i
i
i
// y
if (Tr.y() < 1e-3) {
const float continueProb = .5f;
if (rng RandomFloat() > continueProb) break;
if (rng.RandomFloat() > continueProb) break;
Tr /= continueProb;
}
// Compute emission only source term at p // Compute emission-only source term at _p_
Lv += Tr * vr->Lve(p, w, ray.time);
}
*
t t
N*T = Tr;
return Lv * step;
N
i
i ev
i
x L x
x N Tr
t t
1 0
1
( ) ( , )
Emission only
exponential density
Single scattering
• Consider incidence radiance due to direct illumination
illumination
ds x S x x Tr x
L x x Tr x
L
d
) , ' ( ) ' ( ) , ( ) (
) ,
(
0
0
0
' ) ' , ( ) ' , ( ) , ( ) , ( ) ,
( x L x x p x L x d
S
ve s
d
) ( L
x )
, ( x
0 L
' x x
Single scattering
• Consider incidence radiance due to direct illumination
illumination
ds x S x x Tr x
L x x Tr x
L
d
) , ' ( ) ' ( ) , ( ) (
) ,
(
0
0
' ) ' , ( ) ' , ( ) , ( ) , ( ) ,
( x L x x p x L x d
S
ve s
d
0
) ( L
x )
, ( x
0 L
x
Single scattering
• L
dmay be attenuated by participating media A h i f h i l ld
• At each point of the integral, we could use multiple importance sampling to get
' ) ' , ( ) ' , ( ) ,
(
sx p x L
dx d
But, in practice, we can just pick up light source randomly
source randomly.
Single scattering
Subsurface scattering
• The bidirectional scattering-surface reflectance distribution function (BSSRDF)
distribution function (BSSRDF)
o
ip
op
i) , , ,
( p
o op
i iS ( p
o,
o, p
i,
i)
dA d p
L p
p S p
L ( p , ) S ( p , , p
i,
i) L
i( p
i,
i) cos
id
idA
L
iA H
i i
i i i i o o o
o
o
2
cos ) , ( ) , , , ( )
, (
Subsurface scattering
• Translucent materials have similar mechanism for light scattering as participating media for light scattering as participating media.
Thus, path tracing could be used.
H t l t bj t h
• However, many translucent objects have very high albedo. Taken milk as an example, after 100 tt i t 87 5% f th i id t 100 scattering events, 87.5% of the incident light is still carried by a path, 51% after 500 and 26% ft 1 000
26% after 1,000.
• Efficiently rendering these kinds of translucent scattering media requires a different approach.
Highly translucent materials Main idea
• Assume that the material is homogeneous and the medium is semi infinite we can use
the medium is semi-infinite, we can use diffusion approximation to describe the equlibrium distribution of illumination equlibrium distribution of illumination.
• There is a solution to the diffusion equation by i di l f t li ht t
using a dipole of two light sources to
approximate the overall scattering.
Highly scattering media Principle of similarity
g=0.9 For high albedo objects, an
anisotropically scattering phase anisotropically scattering phase function becomes isotropic after many scattering events.
n=10
n=100 n=1000
Diffusion approximation
• The reduced scattering coefficient
s
s
g
' ( 1 )
• The reduced extinction coefficient
'
' '
'
s a
t
Diffusion approximation
For a point light source
at with power p p L p d
S i
2( , ) )
(
at with power p
fluence
p p
a t
tr
3
') 3 ( 1 ( 3
')
1
tD diffusion coefficient
p
e
pp p
i
tr
( )
The fluence at f i p
ip
coefficient
i
i
D p p
p
p
) 4 ,
from is p (
Dipole diffusion approximation
)
( 1 F
dr
A
z ) z
( 1
) (
dr dr
A F
H r i i id
dr
F n d
F ( )
2( , )
AD z 2
0 d
z
z
z
Dipole diffusion approximation
Light enters the material
at p ( x y 0 ) negative light
at p
i ( x
i, y
i, 0 )
) , , ( x
iy
iz
negative light at
p
ip
opositive light
) , , ( x
iy
iz
at
e
p e p
d
d tr
tr
( )
d x
iy
iz p
o
( , , )
d d
p D p
i o
) 4 , (
o i
i
y z p
x
d
( , ,
)
Dipole diffusion approximation
• Put all together
3 3
) 1 (
) 1 (
4 ) 1 , (
d e d
z d
e d
p z p R
d tr
d tr
o i d
tr
tr
) ( p p R
d( p
i, p
o) R
o
i
p
p
BSSRDF
• BSSRDF based on the diffusion subsurface reflectance approximation
reflectance approximation
) , ( ) (
) , 1 (
) , , ,
( p p
i iF
tR
dp
ip F
t i iS ( p
o,
o, p
i,
i) F
t(
o,
o) R
d( p
i p
o) F
t(
i,
i)
S
i
op
ip
oEvaluating BSSRDF
o
p
ip p
L ( )
p
op
idA d p
L F
p p R F
p L
i i i
i i i i t i
d t
o o o
cos ) , ( ) , ( ) (
) , 1 (
) , (
F R p p F L
ip
i i id
idA
A H
i i t o i d o o
t
( , ) ( ) ( , ) ( , ) cos
2
Evaluating BSSRDF
p L
o o o
)
,
(
F L p d dA
p p R F
A H
i i i i i i i t o
i d o o
t
2
cos ) , ( ) , ( ) (
) , 1 (
A
H
) ( )
(
i idt
E p
F
F h t i l
For homogeneous materials,
o
idA p E p p R F
F p L
o o o
1 ( ) ( ) ( ) ( )
) , (
dA p E p p R F
F
A
i o i d dt
o
t
( , ) ( ) ( ) ( )
Evaluating BSSRDF
op
op
ip L
o o o1 ) ,
(
dA p E p p R F
F
A
i o i d dt
o
t
1 ( , ) ( ) ( ) ( )
j
j j o j d A
i o i d o
o
p R p p E p dA R p p E A
M ( ) ( ) ( ) ( )
Single scattering
Solutions
• Path tracing
Ph i
• Photon mapping
• Hierarchical approach (Jensen 2002)
Three main components
• Sample a large number of random points on the surface and their incident irradiance is
surface and their incident irradiance is computed.
C t hi h f th i t
• Create a hierarchy of these points,
progressively clustering nearby points and i th i i di l
summing their irradiance values.
• At rendering, use a hierarchical integration algorithm to evaluate the subsurface scattering equation.
DipoleSubsurfaceIntegrator
class DipoleSubsurfaceIntegrator : public SurfaceIntegrator {g {
…
int maxSpecularDepth;
float maxError, minSampleDist;
string filename;
vector<IrradiancePoint> irradiancePoints;
vector<IrradiancePoint> irradiancePoints;
BBox octreeBounds;
SubsurfaceOctreeNode *octree;
}
Sample points
• Samples should be reasonably uniformly
distributed A Poisson disk distribution of points distributed. A Poisson disk distribution of points is a good choice.
Th l ith f ti h
• There are some algorithms for generating such distributions. Pbrt uses a 3D version of dart
th i b f i P i h t t
throwing by performing Poisson sphere tests.
The algorithm terminates when a few thousand t t h b j t d i
tests have been rejected in a row.
Bad, but dipole is dad for this anyway
good anyway
Sample points Preprocess
if (scene->lights.size() == 0) return;
<Get SurfacePoints for translucent objects>j
<Compute irradiance values at sample points>
<Create octree of clustered irradiance samples>
<Compute irradiance values at sample points>
for (uint32_t i = 0; i < pts.size(); ++i) { SurfacePoint &sp = pts[i];p p [ ];
Spectrum E(0.f);
for (int j = 0; j < scene->lights.size(); ++j) {
<Add irradiance from light at point>
}
irradiancePoints push back(IrradiancePoint(sp E));
Calculate direct lighting only; it is possible to include indirect lighting but more expensive
irradiancePoints.push_back(IrradiancePoint(sp, E));
}
S
pSurfacePoint
H
L
ip
i id
ip
E ( )
2( , ) | cos | p
L ( ) | cos |
1
S
p
j j
j j
i
p p L p N
E ( )
| cos
| ) , 1 (
)
(
ip
IrradiancePoint
struct SurfacePoint {
…
Point p;
Normal n;
float area, rayEpsilon;
};
struct IrradiancePoint {
…
Point p;
Normal n;
Spectrum E;
float area, rayEpsilon;
};
};
<Create octree of clustered irradiance samples>
octree = octreeArena.Alloc<SubsurfaceOctreeNode>();
for (int i = 0; i < irradiancePoints.size(); ++i)( ; (); ) octreeBounds = Union(octreeBounds,
irradiancePoints[i].p);
for (int i = 0; i < irradiancePoints.size(); ++i) octree->Insert(octreeBounds, &irradiancePoints[i],
octreeArena);
octreeArena);
octree->InitHierarchy();
Computes representative position irradiance and area for each node Computes representative position, irradiance and area for each node.
Positions are weighted by irradiance values to emphasize the points with higher irradiance values.
Li
) ( ) ( ) , 1 ( ) ,
(
o o t o dt o oo
p F F M p
L
if (bssrdf && octree) {
Spectrum sigma_a = bssrdf->sigma_a();
S t i b df > i i ()
Spectrum sigmap_s = bssrdf->sigma_prime_s();
Spectrum sigmap_t = sigmap_s + sigma_a;
if (!sigmap t.IsBlack()) { if (!sigmap_t.IsBlack()) {
DiffusionReflectance Rd(sigma_a, sigmap_s, bssrdf->eta());
Mo = octree->Mo(octreeBounds, p, Rd, …);
FresnelDielectric fresnel(1.f, bssrdf->eta());
Ft S t (1) f l E l t (Ab D t( ))
Ft = Spectrum(1)- fresnel.Evaluate(AbsDot(wo,n));
float Fdt = 1.f - Fdr(bssrdf->eta());
L += (INV PI * Ft) * (Fdt * Mo);( _ ) ( );
} }
Li
L += UniformSampleAllLights(…);
if (ray.depth < maxSpecularDepth) {( y p p p ) { L += SpecularReflect(…);
L += SpecularTransmit(…);
}
return L;
Mo
d j o j jo
o
p R p p E A
M ( ) ( )
j
p
oSubsurfaceOctreeNode::Mo
Spectrum SubsurfaceOctreeNode::Mo(BBox &nodeBound, Point &pt, DiffusionReflectance &Rd, float maxErrorp , , {
If extended solid angle of the node is small enough and the point is not float dw = sumArea / DistanceSquared(pt, p);
if (dw < maxError && !nodeBound Inside(pt))
inside the node, use the representative values of the node to estimate.
if (dw < maxError && !nodeBound.Inside(pt))
return Rd(DistanceSquared(pt, p)) * E * sumArea;
Spectrum Mo = 0.f;
SubsurfaceOctreeNode::Mo
if (isLeaf)
for (int i = 0; i < 8; ++i) {( ; ; ) { if (!ips[i]) break;
Mo += Rd(DistanceSquared(pt, ips[i]->p))
* ips[i]->E * ips[i]->area;
} else { else {
Point pMid=.5f*nodeBound.pMin+.5f*nodeBound.pMax;
for (int child = 0; child < 8; ++child) {( ) { if (!children[child]) continue;
Bbox cBound=octreeChildBound(child,nodeBound,pMid);
Mo+=children[child]->Mo(cBound, pt, Rd, maxError);
} } }
return Mo;
}
