### Sampling and Reconstruction

### Digital Image Synthesis g g y *Yung-Yu Chuang*

*with slides by Pat Hanrahan, Torsten Moller and Brian Curless*

**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

**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*

**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

**Jaggies**

**Retort sequence by Don Mitchell**

**Staircase pattern or jaggies**

**Moire pattern**

### • Sampling the equation

### equation

### )

### sin( *x*

^{2}

### *y*

^{2}

**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)

**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

_{}

###

**Fourier analysis**

*spatial domain* *frequency domain*

**Fourier analysis**

*spatial domain* *frequency domain*

**Fourier analysis**

*spatial domain* *frequency domain*

**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*

**1D convolution theorem example**

**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}*)*

**The delta function**

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

### and unit area

**Sifting and shifting**

**Shah/impulse train function**

*frequency domain* *spatial domain*

### ,

**Sampling**

*band limited*

**Reconstruction**

*The reconstructed function is obtained by interpolating *

*The reconstructed function is obtained by interpolating *

*among the samples in some manner*

**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*

**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.

**Aliasing**

*increase sample * *decrease sample *

*increase sample * *spacing in*

*spatial domain*

*p* *spacing in*

*frequency domain*

**Aliasing**

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

*regions*

**Sampling theorem**

**Sampling theorem**

**Aliasing due to under-sampling**

**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.

**Aliasing**

### • Prealiasing: due to sampling under Nyquist rate

### P li i d f i f

### • Postaliasing: due to use of imperfect

### reconstruction filter

**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

**Antialiasing (Prefiltering)**

*It is blurred, but*

*better than aliasing*

**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**

**Point vs. Supersampled**

**Point** **4x4 Supersampled**

**Checkerboard sequence by Tom Duff**

**Analytic vs. Supersampled**

**Exact Area** **4x4 Supersampled**

**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.

**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.

**Jittered vs. Uniform Supersampling**

**4x4 Jittered Sampling** **4x4 Uniform**

**Prefer noise over aliasing**

*reference* *aliasing* *noise*

**Jittered sampling**

**Add uniform random jitter to each sample**

**Add uniform random jitter to each sample**

**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)

**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**

**Example**

**Aliasing**

### frequency domain

### function (a) function (b) domain

### alias=false

### frequency

**Stochastic sampling**

**Stochastic sampling**

### function (a) function (b)

### Replace structure

### alias by structureless

### alias by structureless

### (high-freq) noise

**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.

**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.

**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**

**An ineffective sampler**

**A more effective sampler**

**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 *Ls = new Spectrum[maxSamples];**

**Spectrum *Ts = new Spectrum[maxSamples];**

**Intersection *isects = new Intersection[maxSamples]**

**...**

**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);**

**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;**

**}**

**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));**

**}**

**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;**

**};**

**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.

**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**

**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();**

**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);**

**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];**

**}**
**}**
**} **

**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);**

**}**
**}**

**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**

**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;**

**}**

**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);**

**}**
**}**
**}**

**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*

**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);**

**}**
**}**

**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;**

**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;**

**samplePos = 0;**

**}**

**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], **

**shutterOpen, shutterClose);**

**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;**

**return 1;**

**}**

**Random sampling**

*a pixel*

*completely *

*random*

*random*

**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.

**Stratified sampling**

*completely * *random*

*stratified * *uniform*

*stratified * *jittered*

*random* *uniform* *jittered*

### turns aliasing into noise

### into noise

**Comparison of sampling methods**

*256 l * *i l * *f*

*256 samples per pixel as reference*

*1 sample per pixel (no jitter)*

**Comparison of sampling methods**

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

*4 samples per pixel (jittered)*

**Stratified sampling**

*reference* *random* *stratified *

*ji* *d*

*jittered*

**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.

**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, **

**yPixelSamples, rng, jitterSamples);**

**StratifiedSample1D(timeSamples, xPixelSamples * **

**yPixelSamples** **rng** **jitterSamples);**

**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;**

**}**

**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** **[ ],**

**shutterOpen, shutterClose);**

**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;**

**}**

**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);**

**}**
**}**

**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

**}**

**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

### A worst case for stratified sampling

### LHS can prevent this to happen

**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]);**

**}**
**}**

**}**
**}**

**Stratified sampling**

**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*

**Low discrepancy sampling**

### • A possible problem with stratified sampling

### • Discrepancy can be used to evaluate the quality of patterns

### of patterns

**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

**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.

**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;**

**}**

**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

**High-dimensional sequence**

### • Two well-known low-discrepancy sequences

### H lt – Halton

### – Hammersley

**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

**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.

**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.

**Radial inverse**

*Halton* *Hammersley*

*Better for that there *

*Better for that there *

*are fewer clumps.*

**Folded radial inverse**

*Halton* *Hammersley*

*The improvement is*

*The improvement is*

*more obvious*

**Low discrepancy sampling**

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

*Hammersley sequence, 1 sample/pixel*

**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

**Best candidate sampling**

*stratified jittered* *best candidate*

*It avoids holes and clusters*

*It avoids holes and clusters.*

**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**

**Best candidate sampling**

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

*best candidate, 1 sample/pixel*

**Best candidate sampling**

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

*best candidate, 4 sample/pixel*

**Comparisons**

*reference* *low-discrepancy* *best candidate*

**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.

**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);**

**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

**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. **

**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;**

**}**
**}**

**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;**

**return true;**

**}**
**}**

**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 …

**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.

**Adaptive sampling**

**Adaptive sampling (geometry)**

**Adaptive sampling (contrast)**

**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.

**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

### ^{} ^{}

###

^{i}*f* *x* *x*

^{i}*y* *y*

^{i}*L* *x*

^{i}*y*

^{i}*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

**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)**

**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.

**}**

**Triangle filter**

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

**{**

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

**}**

**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.

**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.

**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^{0} and C^{1}
Continuity.

**Windowed sinc filter**

*Lanczos*

###

###

###

###

### / / ) sin

### ( *x*

*x* *x*

*w*

*sinc*

**Comparisons**

*box*

*Mitchell*

**Comparisons**

*windowed sinc*

*Mitchell*

**Comparisons**

*box* *Gaussian* *Mitchell*

**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.***

**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;**

**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**
**};**

**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>**

**<precompute filter table>**

**}**

**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* ( , )

) , ( ) ,

) ( , (