Sampling theory

(1)

Sampling and Reconstruction

Digital Image Synthesis g g y Yung-Yu Chuang

with slides by Pat Hanrahan, Torsten Moller and Brian Curless

(2)

Sampling theory

• Sampling theory: the theory of taking discrete sample values (grid of color pixels) from

sample values (grid of color pixels) from functions defined over continuous domains

(incident radiance defined over the film plane) (incident radiance defined over the film plane) and then using those samples to reconstruct new functions that are similar to the original

functions that are similar to the original (reconstruction).

Sampler : selects sample points on the image plane

Filter : blends multiple samples together

(3)

Aliasing

• Reconstruction generates an approximation to the original function Error is called aliasing

the original function. Error is called aliasing.

l l

sampling reconstruction

sample value

l i i

sample position

(4)

Sampling in computer graphics

• Artifacts due to sampling - Aliasing

J i – Jaggies – Moire

Fli k i ll bj t – Flickering small objects – Sparkling highlights

l b h h l ff

– Temporal strobing (such as Wagon-wheel effect)

• Preventing these artifacts - Antialiasing

(5)

Jaggies

Retort sequence by Don Mitchell

Staircase pattern or jaggies

(6)

Moire pattern

• Sampling the equation

equation

)

sin( x

²

 y

²

(7)

Fourier analysis

• Can be used to evaluate the quality between the reconstruction and the original

the reconstruction and the original.

• The concept was introduced to Graphics by

R b t C k i 1986 ( t d d b D Mit h ll) Robert Cook in 1986. (extended by Don Mitchell)

Rob Cook V.P. of Pixar

1981 M.S. Cornell 1981 M.S. Cornell

1987 SIGGRAPH Achievement award

1999 Fellow of ACM

2001 Academic Award with Ed Catmull and Loren Ed Catmull and Loren

Carpenter (for Renderman)

(8)

Fourier transforms

• Most functions can be decomposed into a weighted sum of shifted sinusoids

weighted sum of shifted sinusoids.

• Each function has two representations

– Spatial domain - normal representation

– Frequency domain - spectral representation

• The Fourier transform converts between the spatial and frequency domain

Spatial Frequency

( ) ( )

^{i x}

F  

^

 f x e

^ ^

dx

Spatial Domain

Frequency Domain

( ) 1 ( )

2 f x F  e

i x^

d 







 

) (x

f F (  )

2 

_



(9)

Fourier analysis

spatial domain frequency domain

(10)

Fourier analysis

spatial domain frequency domain

(11)

Fourier analysis

spatial domain frequency domain

(12)

Convolution

• Definition



( ) ( ) ( )

h x    f g  f x g x   x dx  

• Convolution Theorem: Multiplication in the frequency domain is equivalent to convolution in the space

domain domain.

f    g F G

• Symmetric Theorem: Multiplication in the space

domain is equivalent to convolution in the frequency domain is equivalent to convolution in the frequency domain.

f g    F G

f g

(13)

1D convolution theorem example

(14)

2D convolution theorem example

f(x,y) g(x,y) h(x,y)

 

  

F(s

_x

,s

_y

) G(s

_x

,s

_y

) H(s

_x

,s

_y

)

(15)

The delta function

• Dirac delta function, zero width, infinite height and unit area

and unit area

(16)

Sifting and shifting

(17)

Shah/impulse train function

frequency domain spatial domain

,

(18)

Sampling

band limited

(19)

Reconstruction

The reconstructed function is obtained by interpolating

among the samples in some manner

(20)

In math forms

) ( )

III( ) )

(

~ (

F

F ( F ( s ) III(s) )  ) ( s

F    

) ( sinc )

III )

(

~ (

x (x)

x f

f   



^

^



i

i f i x

x

f ( ) sinc ( ) ( )

~



 i

(21)

Reconstruction filters

The sinc filter, while ideal, has two drawbacks:

• It has a large support (slow to compute)

I i d i i i

• It introduces ringing in practice

The box filter is bad because its Fourier transform is a sinc its Fourier transform is a sinc filter which includes high

frequency contribution from th i fi it i f th the infinite series of other copies.

(22)

Aliasing

increase sample decrease sample

increase sample spacing in

spatial domain

p spacing in

frequency domain

(23)

Aliasing

high-frequency d t il l k i t details leak into lower-frequency regions

regions

(24)

Sampling theorem

(25)

Sampling theorem

(26)

Aliasing due to under-sampling

(27)

Sampling theorem

• For band limited functions, we can just increase the sampling rate

increase the sampling rate

• However, few of interesting functions in t hi b d li it d i computer graphics are band limited, in particular, functions with discontinuities.

• It is mostly because the discontinuity always falls between two samples and the samples

provides no information about this discontinuity.

(28)

Aliasing

• Prealiasing: due to sampling under Nyquist rate

P li i d f i f

• Postaliasing: due to use of imperfect

reconstruction filter

(29)

Antialiasing

• Antialiasing = Preventing aliasing

1. Analytically prefilter the signal

– Not solvable in general

2. Uniform supersampling and resample p p g p

3. Nonuniform or stochastic sampling

(30)

Antialiasing (Prefiltering)

(31)

It is blurred, but

better than aliasing

(32)

Uniform supersampling

• Increasing the sampling rate moves each copy of the spectra further apart potentially

of the spectra further apart, potentially reducing the overlap and thus aliasing

R lti l t b l d (filt d)

• Resulting samples must be resampled (filtered) to image sampling rate

s s

s

Pixel   w Sample 

Samples Pixel

(33)

Point vs. Supersampled

Point 4x4 Supersampled

Checkerboard sequence by Tom Duff

(34)

Analytic vs. Supersampled

Exact Area 4x4 Supersampled

(35)

Non-uniform sampling

• Uniform sampling

– The spectrum of uniformly spaced samples is also a set of The spectrum of uniformly spaced samples is also a set of uniformly spaced spikes

– Multiplying the signal by the sampling pattern corresponds to

l i f h h ik (i f )

placing a copy of the spectrum at each spike (in freq. space) – Aliases are coherent (structured), and very noticeable

• Non uniform sampling

• Non-uniform sampling

– Samples at non-uniform locations have a different spectrum; a single spike plus noise

– Sampling a signal in this way converts aliases into broadband noise

Noise is incoherent (structureless) and much less objectionable – Noise is incoherent (structureless), and much less objectionable

• Aliases can’t be removed, but can be made less noticeable

noticeable.

(36)

Antialiasing (nonuniform sampling)

• The impulse train is modified as



^





 



 



 



 iT  

x- 

 2

1 • It turns regular aliasing into noise But random







 

  

i

2 • It turns regular aliasing into noise. But random

noise is less distracting than coherent aliasing.

(37)

Jittered vs. Uniform Supersampling

4x4 Jittered Sampling 4x4 Uniform

(38)

Prefer noise over aliasing

reference aliasing noise

(39)

Jittered sampling

Add uniform random jitter to each sample

(40)

Poisson disk noise (Yellott)

• Blue noise

S h ld b i d l k

• Spectrum should be noisy and lack any concentrated spikes of energy (to avoid

h t li i ) coherent aliasing)

• Spectrum should have deficiency of low-

frequency energy (to hide aliasing in less

noticeable high frequency)

(41)

Distribution of extrafoveal cones

Monkey eye cone distribution

Fourier transform cone distribution

Yellott theory



Aliases replaced by noise p y



Visual system less sensitive to high freq noise

(42)

Example

(43)

Aliasing

frequency domain

function (a) function (b) domain

alias=false

frequency

(44)

Stochastic sampling

(45)

Stochastic sampling

function (a) function (b)

Replace structure

alias by structureless

(high-freq) noise

(46)

Antialiasing (adaptive sampling)

• Take more samples only when necessary.

However in practice it is hard to know where However, in practice, it is hard to know where we need supersampling. Some heuristics could be used

be used.

• It only makes a less aliased image, but may not

b ffi i t th i l li

be more efficient than simple supersampling

particular for complex scenes.

(47)

Application to ray tracing

• Sources of aliasing: object boundary, small objects textures and materials

objects, textures and materials

• Good news: we can do sampling easily

• Bad news: we can’t do prefiltering (because we do not have the whole function)

• Key insight: we can never remove all aliasing,

so we develop techniques to mitigate its impact p q g p

on the quality of the final image.

(48)

pbrt sampling interface

• Creating good sample patterns can substantially improve a ray tracer’s efficiency allowing it to improve a ray tracer s efficiency, allowing it to create a high-quality image with fewer rays.

B l ti di i tl it t

• Because evaluating radiance is costly, it pays to spend time on generating better sampling.

• core/sampling., samplers/

• random.cpp, stratified.cpp, pp pp bestcandidate.cpp,

lowdiscrepancy.cpp, p y pp

(49)

An ineffective sampler

(50)

A more effective sampler

(51)

Main rendering loop for each task

void SamplerRenderer::Run() { Sampler *sampler

= mainSampler >GetSubSampler(taskNum taskCount);

= mainSampler->GetSubSampler(taskNum, taskCount);

...

// Allocate space for samples and intersections int maxSamples = sampler >MaximumSampleCount();

int maxSamples = sampler->MaximumSampleCount();

Sample *samples=origSample->Duplicate(maxSamples);

RayDifferential *rays=new RayDifferential[maxSamples];

Spectrum *Ls = new Spectrum[maxSamples];

Spectrum *Ts = new Spectrum[maxSamples];

Intersection *isects = new Intersection[maxSamples]

...

(52)

Main rendering loop for each task

while ((sCnt=sampler->GetMoreSamples(samples,rng))>0){

for (int i = 0; i < sCnt; ++i) { ...

float rayWeight = camera->GenerateRayDifferential(

samples[i], &rays[i]);

...

if (rayWeight > 0.f)

Ls[i] = rayWeight * renderer->Li(scene, rays[i],

l [i] i [i] [i])

&samples[i], rng, arena, &isects[i], &Ts[i]);

...

} // end for

if ( l >R tR lt ( l L ))

if (sampler->ReportResults(samples, rays, Ls, …)) for (int i = 0; i < sCnt; ++i) {

...

>fil >AddS l ( l [i] L [i]) camera->film->AddSample(samples[i], Ls[i]);

}

} // end while

camera >film >UpdateDisplay(

camera->film->UpdateDisplay(

sampler->xPixelStart, sampler->yPixelStart, sampler->xPixelEnd+1, sampler->yPixelEnd+1);

(53)

Sampler

Generates a good pattern of multidimensional samples.

class Sampler { class Sampler {

…

virtual int **GetMoreSamples(Sample sample virtual int GetMoreSamples(Sample sample,**

RNG &rng) = 0;

virtual int MaximumSampleCount() 0;

random number

generator for pre-allocating memory virtual int MaximumSampleCount() = 0;

virtual bool ReportResults(…);

report radiance for things like adaptive sampling

…

const int xPixelStart, xPixelEnd;

t i t Pi lSt t Pi lE d range of pixels report radiance for things like adaptive sampling

const int yPixelStart, yPixelEnd;

const int samplesPerPixel;

g p

sample per pixel const float shutterOpen, shutterClose;

}

(54)

Sampler

void Sampler::ComputeSubWindow(int num, int count, int *XStart, int *XEnd, int *YStart, int *YEnd) {, , , ) {

int dx=xPixelEnd-xPixelStart,dy=yPixelEnd-yPixelStart;

int nx = count, ny = 1;

while ((nx & 0x1) == 0 && 2 * dx * ny < dy * nx) { nx >>= 1; ny <<= 1;

} }

int xo = num % nx, yo = num / nx;

float tx0=float(xo)/float(nx),tx1=float(xo+1)/float(nx);( ) ( ), ( ) ( ) float ty0=float(yo)/float(ny),ty1=float(yo+1)/float(ny);

*XStart = Floor2Int(Lerp(tx0,xPixelStart,xPixelEnd));

*XEnd = Floor2Int(Lerp(tx1,xPixelStart,xPixelEnd));

*YStart = Floor2Int(Lerp(ty0,yPixelStart,yPixelEnd));

*YEnd = Floor2Int(Lerp(ty1 yPixelStart yPixelEnd));

*YEnd = Floor2Int(Lerp(ty1,yPixelStart,yPixelEnd));

}

(55)

Sample

struct CameraSample { float imageX, imageY;

store required information for generating camera rays

g , g ;

float lensU, lensV;

float time;

};

struct Sample : public CameraSample {

Sample(Sampler *sampler SurfaceIntegrator *surf store required information for one eye ray sample Sample(Sampler *sampler, SurfaceIntegrator *surf,

VolumeIntegrator *vol, const Scene *scene);

uint32 t Add1D(uint32 t num);_ ( _ )

uint32_t Add2D(uint32_t num);

…

Note that it stores all samples

// Sample Public Data

vector<uint32_t> n1D, n2D;

float **oneD **twoD;

Note that it stores all samples required for one eye ray. That is, it may depend on depth.

float **oneD, **twoD;

};

(56)

Sample

• Sample is allocated once in Render(). Sampler is called to fill in the information for each eye ray The

to fill in the information for each eye ray. The

integrator can ask for multiple 1D and/or 2D samples, each with an arbitrary number of entries, e.g.

depending on #lights. For example, WhittedIntegrator does not need samples. DirectLighting needs samples proportional to #lights

proportional to #lights.

• The structure of sample is initiated once and Sampler is responsible for filling in requested sample structure

responsible for filling in requested sample structure

with well-behaved samples.

(57)

Data structure

•Different types of lights require different numbers of samples, usually 2D samples.

samples, usually 2D samples.

•Sampling BRDF requires 2D samples.

•Selection of BRDF components requires 1D samples.

3 1 2

D t D

n1D n2D 2 2 1 1 2 2

sample

allocate together to avoid cache miss filled in by integrators

oneD twoD allocate together to avoid cache miss

mem

lightNumOffset 1 2 3

lightSampleOffset bsdfSampleOffset

1 3 5 2 4 6

PathIntegrator

lightNumOffset lightSampleOffset bsdfSampleOffset

(58)

Sample

Sample::Sample(Sampler *sampler, SurfaceIntegrator

*surf, VolumeIntegrator *vol, Scene *scene) {

if (surf) surf->RequestSamples(sampler,this,scene);

if ( l) l S l ( l hi )

if (vol) vol->RequestSamples(sampler, this, scene);

AllocateSampleMemory();

} }

void Sample::AllocateSampleMemory() { int nPtrs = n1D.size() + n2D.size();

if (!nPtrs) {

D t D NULL t

oneD = twoD = NULL; return;

}

oneD = AllocAligned<float *>(nPtrs);

oneD AllocAligned<float >(nPtrs);

twoD = oneD + n1D.size();

(59)

Sample

int totSamples = 0;

for (uint32 t i = 0; i < n1D.size(); ++i)( _ ; (); ) totSamples += n1D[i];

for (uint32_t i = 0; i < n2D.size(); ++i) totSamples += 2 * n2D[i];

float *mem = AllocAligned<float>(totSamples);

for (uint32_t i = 0; i < n1D.size(); ++i) { oneD[i] = mem; mem += n1D[i];[ ] [ ]

}

for (uint32_t i = 0; i < n2D.size(); ++i) { twoD[i] = mem; mem += 2 * n2D[i];

} } }

(60)

PathIntegrator::RequestSamples

void PathIntegrator::RequestSamples(Sampler *sampler, Sample *sample, const Scene *scene)

{

for (int i = 0; i < SAMPLE_DEPTH; ++i) {

lightSampleOffsets[i]=LightSampleOffsets(1 sample);

lightSampleOffsets[i]=LightSampleOffsets(1,sample);

bsdfSampleOffsets[i]=BSDFSampleOffsets(1,sample);

pathSampleOffsets[i]=BSDFSampleOffsets(1,sample);

} }

(61)

LightSampleOffsets

struct LightSampleOffsets {

LightSampleOffsets(int count, Sample *sample);

int nSamples, componentOffset, posOffset;

};

LightSampleOffsets::LightSampleOffsets(int count, Sample *sample) {

l

nSamples = count;

componentOffset = sample->Add1D(nSamples);

posOffset = sample->Add2D(nSamples);

p p ( p )

}

Li htS l Off t Li htS l LightSampleOffsets LightSample

(62)

DirectLighting::RequestSamples

void DirectLightingIntegrator::RequestSamples(

Sampler *sampler, Sample *sample, Scene *scene) { if ( t t SAMPLE ALL UNIFORM) {

if (strategy == SAMPLE_ALL_UNIFORM) {

uint32_t nLights = scene->lights.size();

lightSampleOffsets=new LightSampleOffsets[nLights];

b dfS l Off t BSDFS l Off t [ Li ht ] bsdfSampleOffsets = new BSDFSampleOffsets[nLights];

for (uint32_t i = 0; i < nLights; ++i) { const Light *light = scene->lights[i];

i S l li h S l

int nSamples = light->nSamples;

if (sampler) nSamples=sampler->RoundSize(nSamples);

lightSampleOffsets[i]

= LightSampleOffsets(nSamples, sample);

bsdfSampleOffsets[i]

= BSDFSampleOffsets(nSamples, sample);

}

lightNumOffset = -1;

}

(63)

DirectLighting::RequestSamples

else {

lightSampleOffsets = new LightSampleOffsets[1];

li htS l Off t [0]

lightSampleOffsets[0]

= LightSampleOffsets(1, sample);

lightNumOffset = sample->Add1D(1);

b dfS l Off t BSDFS l Off t [1]

bsdfSampleOffsets = new BSDFSampleOffsets[1];

bsdfSampleOffsets[0] = BSDFSampleOffsets(1, sample);

} } }

(64)

Random sampler

Just for illustration; does not work well in practice

RandomSampler::RandomSampler(int xstart, int xend, int ystart, int yend, int ns,

float sopen, float sclose) { xPos = xPixelStart;

yPos = yPixelStart;

nSamples = ns;

imageSamples = AllocAligned<float>(5 * nSamples);

lensSamples = imageSamples + 2 * nSamples;p g p p ; timeSamples = lensSamples + 2 * nSamples;

// prepare samples for the first pixel

RNG rng(xstart + ystart * (xend-xstart));

RNG rng(xstart + ystart (xend xstart));

for (int i = 0; i < 5 * nSamples; ++i) imageSamples[i] = rng.RandomFloat();

for (int o = 0; o < 2 * nSamples; o += 2) { for (int o = 0; o < 2 * nSamples; o += 2) {

imageSamples[o] += xPos;

imageSamples[o+1] += yPos; }

l 0

private copy of the current pixel position samplePos = 0;

} #samples consumed for current pixel

(65)

Random sampler

Sampler *RandomSampler::GetSubSampler(

int num, int count) , ) {

int x0, x1, y0, y1;

ComputeSubWindow(num, count, &x0, &x1, &y0, &y1);

if (x0 == x1 || y0 == y1) return NULL;

return new RandomSampler(x0 x1 y0 y1 nSamples return new RandomSampler(x0, x1, y0, y1, nSamples,

shutterOpen, shutterClose);

} }

(66)

Random sampler

int RandomSampler::GetMoreSamples(

Sample *sample, RNG &rng) {p p , g) { if (samplePos==nSamples) {

if(xPixelStart==xPixelEnd || yPixelStart==yPixelEnd) return 0;

if (++xPos == xPixelEnd) {xPos=xPixelStart; ++yPos;}

if (yPos == yPixelEnd) return 0;

for (int i = 0; i < 5 * nSamples; ++i)

generate all samples for one pixel at once

( p )

imageSamples[i] = rng.RandomFloat();

for (int o = 0; o < 2 * nSamples; o += 2) {

imageSamples[o]+=xPos; imageSamples[o+1]+=yPos;

}

samplePos = 0;

}

(67)

Random sampler

sample->imageX = imageSamples[2*samplePos];

sample->imageY = imageSamples[2*samplePos+1];p g g p [ p ];

sample->lensU = lensSamples[2*samplePos];

sample->lensV = lensSamples[2*samplePos+1];

sample->time = Lerp(timeSamples[samplePos],

for (uint32 t i = 0; i < sample >n1D size(); ++i) for (uint32_t i = 0; i < sample->n1D.size(); ++i)

for (uint32_t j = 0; j < sample->n1D[i]; ++j) sample->oneD[i][j] = rng.RandomFloat();p [ ][j] g ()

for (uint32_t i = 0; i < sample->n2D.size(); ++i) for (uint32_t j = 0; j < 2*sample->n2D[i]; ++j)

sample->twoD[i][j] = rng.RandomFloat();

++samplePos;

return 1;

}

(68)

Random sampling

a pixel

completely

random

(69)

Stratified sampling

• Subdivide the sampling domain into non-

overlapping regions (strata) and take a single overlapping regions (strata) and take a single sample from each one so that it is less likely to miss important features

miss important features.

(70)

Stratified sampling

completely random

stratified uniform

stratified jittered

random uniform jittered

turns aliasing into noise

into noise

(71)

Comparison of sampling methods

256 l i l f

256 samples per pixel as reference

1 sample per pixel (no jitter)

(72)

Comparison of sampling methods

1 l i l (ji d) 1 sample per pixel (jittered)

4 samples per pixel (jittered)

(73)

Stratified sampling

reference random stratified

ji d

jittered

(74)

High dimension

• D dimension means N

^D

cells.

S l i k l d i

• Solution: make strata separately and associate

them randomly, also ensuring good distributions.

(75)

StratifiedSampler::GetMoreSamples

if (yPos == yPixelEnd) return 0;

int nSamples = xPixelSamples * yPixelSamples;

// Generate initial stratified samples float *bufp = sampleBuf;

float *imageSamples = bufp; bufp += 2 * nSamples;

float *lensSamples = bufp; bufp += 2 * nSamples;

float *timeSamples = bufp;p p;

StratifiedSample2D(imageSamples, xPixelSamples,

yPixelSamples, rng, jitterSamples);

StratifiedSample2D(lensSamples xPixelSamples StratifiedSample2D(lensSamples, xPixelSamples,

StratifiedSample1D(timeSamples, xPixelSamples *

yPixelSamples rng jitterSamples);

for (int o=0;o<2*xPixelSamples*yPixelSamples;o+=2){

i S l [ ] i S l [ 1]

imageSamples[o]+=xPos; imageSamples[o+1]+=yPos;

}

(76)

StratifiedSampler::GetMoreSamples

Shuffle(lensSamples,xPixelSamples*yPixelSamples,2,rng);

Shuffle(timeSamples,xPixelSamples*yPixelSamples,1,rng);

for (int i = 0; i < nSamples; ++i) {

samples[i].imageX = imageSamples[2*i];

samples[i].imageY = imageSamples[2*i+1];

samples[i].lensU = lensSamples[2*i];

samples[i].lensV = lensSamples[2*i+1];

samples[i].time = Lerp(timeSamples[i], p [ ] p( p [ ],

for (uint32_t j = 0; j < samples[i].n1D.size(); ++j) LatinHypercube(samples[i] oneD[j]

LatinHypercube(samples[i].oneD[j],

samples[i].n1D[j], 1, rng);

for (uint32_t j = 0; j < samples[i].n2D.size(); ++j) LatinHypercube(samples[i] twoD[j]

LatinHypercube(samples[i].twoD[j],

samples[i].n2D[j], 2, rng);

}

if ( i l d) { i lS }

if (++xPos == xPixelEnd) {xPos = xPixelStart; ++yPos;}

return nSamples;

}

(77)

Stratified sampling

void StratifiedSample1D(float *samp, int nSamples, RNG &rng, bool jitter) { /

n stratified samples within [0..1]

float invTot = 1.f / nSamples;

for (int i = 0; i < nSamples; ++i) {

float delta = jitter ? rng.RandomFloat() : 0.5f;

*samp++ = min((i+delta)*invTot, OneMinusEpsilon);

}

} nx*ny stratified samples within [0..1]X[0..1]

}

void StratifiedSample2D(float *samp, int nx, int ny, RNG &rng, bool jitter) {

float dx = 1 f / nx dy = 1 f / ny;

nx ny stratified samples within [0..1]X[0..1]

float dx = 1.f / nx, dy = 1.f / ny;

for (int y = 0; y < ny; ++y)

for (int x = 0; x < nx; ++x) {

float jx = jitter ? rng RandomFloat() : 0 5f;

float jx = jitter ? rng.RandomFloat() : 0.5f;

float jy = jitter ? rng.RandomFloat() : 0.5f;

*samp++ = min((x + jx) * dx, OneMinusEpsilon);

* i (( j ) * d O i il )

*samp++ = min((y + jy) * dy, OneMinusEpsilon);

} }

(78)

Shuffle

template <typename T>

void Shuffle(T *samp, int count, int dims, RNG &rng) ( p, , , g) {

for (int i = 0; i < count; ++i) {

u_int other = i+(rng.RandomUInt()%(count-i));

for (int j = 0; j < dims; ++j)

swap(samp[dims*i + j] samp[dims*other + j]);

swap(samp[dims*i + j], samp[dims*other + j]);

} }

d-dimensional vector swap

}

(79)

Latin hypercube sampling

• Integrators could request an arbitrary n samples.

nx1 or 1xn doesn’t give a good sampling pattern nx1 or 1xn doesn t give a good sampling pattern.

A worst case for stratified sampling

LHS can prevent this to happen

(80)

Latin Hypercube

void LatinHypercube(float *samples,

int nSamples, int nDim, RNG &rng) {

// Generate LHS samples along diagonal float delta = 1.f / nSamples;

for (int i = 0; i < nSamples; ++i) for (int j = 0; j < nDim; ++j)

samples[nDim*i+j] = min((i+(rng.RandomFloat())) p [ j] (( ( g ()))

*delta, OneMinusEpsilon);

// Permute LHS samples in each dimension for (int i = 0; i < nDim; ++i) {

note the difference with shuffle

for (int i = 0; i < nDim; ++i) {

for (int j = 0; j < nSamples; ++j) {

u_int other=j+(rng.RandomUInt() % (nSamples-j));

swap(samples[nDim * j + i]

swap(samples[nDim * j + i],

samples[nDim * other + i]);

} }

(81)

Stratified sampling

(82)

Stratified sampling

1 l d 16 h d l i l

This is better because StratifiedSampler could generate a good LHS pattern for this case

1 camera sample and 16 shadow samples per pixel

16 camera samples and each with 1 shadow sample per pixel

(83)

Low discrepancy sampling

• A possible problem with stratified sampling

• Discrepancy can be used to evaluate the quality of patterns

of patterns

(84)

Low discrepancy sampling

set of N sample points a family of shapes

p p

maximal difference

volume estimated by sample number

real volume

When B is the set of AABBs

with a corner at the origin

with a corner at the origin,

this is called star discrepancy

(85)

1D discrepancy

Uniform is optimal! However, we have learnt that

irregular patterns are perceptually superior to uniform samples. Fortunately, for higher dimension, the low-

discrepancy patterns are less uniform and works discrepancy patterns are less uniform and works reasonably well as sample patterns in practice.

Next, we introduce methods specifically designed for , p y g

generating low-discrepancy sampling patterns.

(86)

Radical inverse

• A positive number n can be expressed in a base b as

• A radical inverse function in base b converts a

nonnegative integer n to a floating-point number in [0,1) g g g p [ , )

inline double RadicalInverse(int n int base) { inline double RadicalInverse(int n, int base) {

double val = 0;

double invBase = 1. / base, invBi = invBase;

while (n > 0) { while (n > 0) {

int d_i = (n % base);

val += d_i * invBi;

/

n /= base;

invBi *= invBase;

}

return val;

}

(87)

van der Corput sequence

• The simplest sequence

R i l li 1D li i h lf l

• Recursively split 1D line in half, sample centers

• Achieve minimal possible discrepancy

(88)

High-dimensional sequence

• Two well-known low-discrepancy sequences

H lt – Halton

– Hammersley

(89)

Halton sequence

• Use relatively prime numbers as bases for each dimension recursively split the dimension

dimension recursively split the dimension into p

_d

parts, sample centers

• Achieve best possible discrepancy for N-D p p y

• Can be used if N is not known in advance

• All prefixes of a sequence are well distributed so as additional samples are added to the

sequence, low discrepancy will be maintained

(90)

Hammersley sequence

• Similar to Halton sequence.

Sli h l b di h H l

• Slightly better discrepancy than Halton.

• Needs to know N in advance.

(91)

Folded radical inverse

• Add the offset i to the ith digit d

_i

and take the modulus b

modulus b.

• It can be used to improve Hammersley and p y

Halton, called Hammersley-Zaremba and

Halton-Zaremba.

(92)

Radial inverse

Halton Hammersley

Better for that there

are fewer clumps.

(93)

Folded radial inverse

Halton Hammersley

The improvement is

more obvious

(94)

Low discrepancy sampling

t tifi d jitt d 1 l / i l stratified jittered, 1 sample/pixel

Hammersley sequence, 1 sample/pixel

(95)

Best candidate sampling

• Stratified sampling doesn’t guarantee good sampling across pixels

sampling across pixels.

• Poisson disk pattern addresses this issue. The

P i di k tt i f i t ith

Poisson disk pattern is a group of points with no two of them closer to each other than some

ifi d di t

specified distance.

• It can be generated by dart throwing. It is time-consuming.

• Best-candidate algorithm by Dan Mitchell. It g y

randomly generates many candidates but only

inserts the one farthest to all previous samples. p p

(96)

Best candidate sampling

stratified jittered best candidate

It avoids holes and clusters

It avoids holes and clusters.

(97)

Best candidate sampling

• Because of it is costly to generate best

candidate pattern pbrt computes a “tilable candidate pattern, pbrt computes a tilable pattern” offline (by treating the square as a rolled torus)

rolled torus).

• tools/samplepat.cpp

→ sampler/bestcandidate o t

→ sampler/bestcandidate.out

(98)

Best candidate sampling

t tifi d jitt d 1 l / i l stratified jittered, 1 sample/pixel

best candidate, 1 sample/pixel

(99)

Best candidate sampling

t tifi d jitt d 4 l / i l stratified jittered, 4 sample/pixel

best candidate, 4 sample/pixel

(100)

Comparisons

reference low-discrepancy best candidate

(101)

Adaptive sampling

• More efficiently generate high-quality images by adding extra samples in parts of the image by adding extra samples in parts of the image that are more complex than others.

b t t t ki d f i l fi t

• pbrt supports two kinds of simple refinement criteria: (1) to check to see if different shapes

i t t d b diff t l i di ti are intersected by different samples, indicating a likely geometric discontinuity and (2) to

h k f i t t b t th l

check for excessive contrast between the colors

of different samples.

(102)

Main rendering loop for each task

while ((sCnt=sampler->GetMoreSamples(samples,rng))>0){

for (int i = 0; i < sCnt; ++i) { ...

float rayWeight = camera->GenerateRayDifferential(

samples[i], &rays[i]);

...

if (rayWeight > 0.f)

Ls[i] = rayWeight * renderer->Li(scene, rays[i],

l [i] i [i] [i])

&samples[i], rng, arena, &isects[i], &Ts[i]);

...

} // end for

if ( l >R tR lt ( l L ))

if (sampler->ReportResults(samples, rays, Ls, …)) for (int i = 0; i < sCnt; ++i) {

...

>fil >AddS l ( l [i] L [i]) camera->film->AddSample(samples[i], Ls[i]);

}

} // end while

camera >film >UpdateDisplay(

camera->film->UpdateDisplay(

sampler->xPixelStart, sampler->yPixelStart, sampler->xPixelEnd+1, sampler->yPixelEnd+1);

(103)

AdaptiveSampler

class AdaptiveSampler : public Sampler {

…

int xPos, yPos;

int minSamples, maxSamples;

current position

least and max number

p , p ;

**float *sampleBuf;** of samples

enum AdaptiveTest { ADAPTIVE_COMPARE_SHAPE_ID, ADAPTIVE_CONTRAST_THRESHOLD };

AdaptiveTest method;

bool supersamplePixel;

which criterion to use

}; whether the current pixel

needs extra samples

(104)

AdaptiveSampler

1. supersamplePixel is set to false initially 2 The initial set of minSamples is generated by 2. The initial set of minSamples is generated by

GetMoreSamples .

3 Rendering loop evaluates these samples and report 3. Rendering loop evaluates these samples and report

them back to ReportResults.

4. If more samples are needed, ReportResults sets to p , p true and leave (xPos, yPos) unchanged.

5. The next call to GetMoreSamples generates a new p g set of maxSamples samples

6. When more samples are not needed or maxSamples

samples have been used, supersamplePixel is set

false and (xPos, yPos) is advanced.

(105)

AdaptiveSampler

int AdaptiveSampler::GetMoreSamples(…) {

…

if (supersamplePixel) {

LDPixelSample(xPos, yPos, shutterOpen, shutterClose, maxSamples, samples,

sampleBuf, rng);

return maxSamples;

} else {

if ( Pos Pi elEnd) ret rn 0 if (yPos == yPixelEnd) return 0;

LDPixelSample(xPos, yPos, shutterOpen, shutterClose, minSamples, samples,

shutterClose, minSamples, samples, sampleBuf, rng);

return minSamples;

} }

(106)

AdaptiveSampler

bool AdaptiveSampler::ReportResults(…) { if (supersamplePixel) {( p p ) {

supersamplePixel = false;

if (++xPos == xPixelEnd) {

xPos = xPixelStart; ++yPos; } return true;

} else if (needsSupersampling( )) { } else if (needsSupersampling(…)) {

supersamplePixel = true;

return false;

} else {

if (++xPos == xPixelEnd) {

xPos = xPixelStart; ++yPos;

}

return true;

} }

(107)

needsSupersampling

bool AdaptiveSampler::needsSupersampling(

Sample *samples, const RayDifferential *rays, p p , y y , const Spectrum *Ls, const Intersection *isects, int count)

{

switch (method) {

case ADAPTIVE COMPARE SHAPE ID:

case ADAPTIVE_COMPARE_SHAPE_ID:

for (int i = 0; i < count-1; ++i)

if (isects[i].shapeId != isects[i+1].shapeId ||( [ ] p [ ] p ||

isects[i].primitiveId != isects[i+1].primitiveId) return true;

Efficient but fails to capture cases like

return false;

(1) coplanar triangles with different ids

but without edges (2) a parametric patch

Co ld fold o er and need more samples

Could fold over and need more samples

(3) shadows, textures …

(108)

needsSupersampling

case ADAPTIVE_CONTRAST_THRESHOLD:

float Lavg = 0.f;g ;

for (int i = 0; i < count; ++i) Lavg += Ls[i].y();

Lavg /= count;

const float maxContrast = 0.5f;

for (int i = 0; i < count; ++i) for (int i = 0; i < count; ++i)

if (fabsf(Ls[i].y() - Lavg) / Lavg > maxContrast) return true;

return false;

}

Not always successful. An example is ImageTexture which has been filtered

return false;

}

For antialiasing. Even if the samples have high contrast, it probably does

N t d l

Not need more samples.

(109)

Adaptive sampling

(110)

Adaptive sampling (geometry)

(111)

Adaptive sampling (contrast)

(112)

Reconstruction filters

• Given the chosen image samples, we can do the following to compute pixel values

following to compute pixel values.

1. reconstruct a continuous function L’ from samples

2 filt L’ t f hi h th

2. prefilter L’ to remove frequency higher than Nyquist limit

3 sample L’ at pixel locations 3. sample L’ at pixel locations

• Because we will only sample L’ at pixel l ti d t d t li itl locations, we do not need to explicitly

reconstruct L’s. Instead, we combine the first

t t

two steps.

(113)

Reconstruction filters

• Ideal reconstruction filters do not exist because of discontinuity in rendering We choose

of discontinuity in rendering. We choose nonuniform sampling, trading off noise for aliasing There is no theory about ideal

aliasing. There is no theory about ideal

reconstruction for nonuniform sampling yet.

I t d id i t l ti bl

• Instead, we consider an interpolation problem

  ^ ^



ⁱ

f x x

ⁱ

y y

ⁱ

L x

ⁱ

y

ⁱ

y

x

I ( , ) ( , )

) ,

( ( )

filter sampled radiance



i

f x ^ x

i

y ^ y

i

y x

I ( , ) ( , ) ( x , y )

) , ( x

_i

y

_i

final value

(114)

Filter

• provides an interface to f(x,y)

il t i t t filt d it t

• Film stores a pointer to a filter and use it to filter the output before writing it to disk.

Filt Filt (fl t fl t )

width, half of support Filter::Filter(float xw, float yw) float Evaluate(float x, float y);

x y is guaranteed to be within the range;

) ( x y f

**• filters/* (box, gaussian, mitchell, sinc,**

x, y is guaranteed to be within the range;

range checking is not necessary

) , ( x y f

filters/ (box, gaussian, mitchell, sinc,

triangle)

(115)

Box filter

• Most commonly used in graphics. It’s just about the worst filter possible incurring postaliasing the worst filter possible, incurring postaliasing by high-frequency leakage.

float BoxFilter::Evaluate(float x, float y)

no need to normalize since the weighted

{

return 1.;

no need to normalize since the weighted sum is divided by the total weight later.

}

(116)

Triangle filter

float TriangleFilter::Evaluate(float x, float y) {

{

return max(0.f, xWidth-fabsf(x)) * max(0.f, yWidth-fabsf(y)); ( , y (y));

}

(117)

Gaussian filter

• Gives reasonably good results in practice

float GaussianFilter::Evaluate(float x float y) float GaussianFilter::Evaluate(float x, float y) {

return Gaussian(x **expX)*Gaussian(y expY);**

**return Gaussian(x, expX)*Gaussian(y, expY);**

}

Gaussian essentially has a infinite support; to compensate this, the value at the end is calculated and subtracted.

this, the value at the end is calculated and subtracted.

(118)

Mitchell filter

• parametric filters, tradeoff between ringing and blurring

and blurring

• Negative lobes improve sharpness; ringing starts t t th i if th b l

to enter the image if they become large.

(119)

Mitchell filter

• Separable filter

• Two parameters, p , B and C, B+2C=1 suggested gg

FFT of a cubic filter.

Mitchell filter is a combination of cubic combination of cubic filters with C⁰ and C¹ Continuity.

(120)

Windowed sinc filter

Lanczos









/ / ) sin

( x

x x

w 

sinc

(121)

Comparisons

box

Mitchell

(122)

Comparisons

windowed sinc

Mitchell

(123)

Comparisons

box Gaussian Mitchell

(124)

Film

• Film class simulates the sensing device in the simulated camera It determines samples’

simulated camera. It determines samples

contributions to the nearby pixels and writes the final floating point image to a file on disk the final floating-point image to a file on disk.

• Tone mapping operations can be used to display th fl ti i t i di l

the floating-point image on a display.

• core/film.*

(125)

Film

class Film { public:

p

Film(int xres, int yres)

: xResolution(xres), yResolution(yres) { } add samples for later reconstruction by weighted average

virtual void AddSample(const CameraSample &sample, const Spectrum &L) = 0;

const Spectrum &L) = 0;

simply sum samples’ contributions, not average them. Pixels with more samples will be brighter. It is used by light transport methods such as MetropolisRender

virtual void Splat(const CameraSample &sample, const Spectrum &L) = 0;

const Spectrum &L) = 0;

(126)

Film

the sample extent could be a bit larger than the pixel extent virtual void GetSampleExtent(int *xstart, p ( ,

int *xend, int *ystart, int *yend) const = 0;

virtual void GetPixelExtent(int *xstart,

int *xend, int *ystart, int *yend) const = 0;

be notified when a region has been recently updated. Do nothing as default.

default.

virtual void UpdateDisplay(int x0, int y0, int x1, int y1, float splatScale = 1.f);

generate the final image for saving to disk or displaying. It accepts a scale factor.

virtual void WriteImage(float splatScale = 1 f)=0;

virtual void WriteImage(float splatScale = 1.f)=0;

const int xResolution, yResolution;y };

(127)

ImageFilm

• film/image.cpp implements the only film

plug in in pbrt It filters samples and writes the plug-in in pbrt. It filters samples and writes the resulting image to disk.

ImageFilm::ImageFilm(int xres, int yres,Filter *filt, float crop[4] string &filename bool openWindow)

float crop[4],string &filename, bool openWindow) {in NDC space. useful for

debugging, or rendering on different computers and

on some system, it can be configured to open a window and show the

image as it’s being rendered ...

i l l k d i l ( i lC

different computers and assembling later

image as it’s being rendered

pixels = new BlockedArray<Pixel>(xPixelCount, yPixelCount);

<precompute filter table>

}

(128)

AddSample

extent

sample  



 

i i i

i i i i i

y y x x f

y x L y y x x y f

x

I ( , )

) , ( ) ,

) ( , (

precomputed

grid of pixels

Filter table

find the nearest neighbor h f l bl

grid of pixels

in the filter table

(129)

Recent progresses on Poisson sampling

• On-the-fly computing

S ll d i [SIGGRAPH 2006]

– Scalloped regions [SIGGRAPH 2006]

• Tile-based

– Recursive Wang tile [SIGGRAPH 2006]

• Parallel

– Li-Yi Wei [SIGGRAPH 2008]

• Show three videos for them Show three videos for them

(130)

Fast Poisson-Disk Sampling

(131)

Fast Poisson-Disk Sampling

(132)