00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024 #if defined(_MSC_VER)
00025 #pragma once
00026 #endif
00027
00028 #ifndef PBRT_CORE_MONTECARLO_H
00029 #define PBRT_CORE_MONTECARLO_H
00030
00031
00032 #include "pbrt.h"
00033 #include "geometry.h"
00034 #include "rng.h"
00035
00036
00037
00038 #ifdef PBRT_IS_WINDOWS
00039
00040 static const float OneMinusEpsilon=0.9999999403953552f;
00041 #else
00042 static const float OneMinusEpsilon=0x1.fffffep-1;
00043 #endif
00044
00045
00046 struct Distribution1D {
00047
00048 Distribution1D(const float *f, int n) {
00049 count = n;
00050 func = new float[n];
00051 memcpy(func, f, n*sizeof(float));
00052 cdf = new float[n+1];
00053
00054 cdf[0] = 0.;
00055 for (int i = 1; i < count+1; ++i)
00056 cdf[i] = cdf[i-1] + func[i-1] / n;
00057
00058
00059 funcInt = cdf[count];
00060 for (int i = 1; i < n+1; ++i)
00061 cdf[i] /= funcInt;
00062 }
00063 ~Distribution1D() {
00064 delete[] func;
00065 delete[] cdf;
00066 }
00067 float SampleContinuous(float u, float *pdf, int *off = NULL) const {
00068
00069 float *ptr = std::upper_bound(cdf, cdf+count+1, u);
00070 int offset = max(0, int(ptr-cdf-1));
00071 if (off) *off = offset;
00072 Assert(offset < count);
00073 Assert(u >= cdf[offset] && u < cdf[offset+1]);
00074
00075
00076 float du = (u - cdf[offset]) / (cdf[offset+1] - cdf[offset]);
00077 Assert(!isnan(du));
00078
00079
00080 if (pdf) *pdf = func[offset] / funcInt;
00081
00082
00083 return (offset + du) / count;
00084 }
00085 int SampleDiscrete(float u, float *pdf) const {
00086
00087 float *ptr = std::upper_bound(cdf, cdf+count+1, u);
00088 int offset = max(0, int(ptr-cdf-1));
00089 Assert(offset < count);
00090 Assert(u >= cdf[offset] && u < cdf[offset+1]);
00091 if (pdf) *pdf = func[offset] / (funcInt * count);
00092 return offset;
00093 }
00094 private:
00095 friend struct Distribution2D;
00096
00097 float *func, *cdf;
00098 float funcInt;
00099 int count;
00100 };
00101
00102
00103 void RejectionSampleDisk(float u1, float u2, float *x, float *y);
00104 Vector UniformSampleHemisphere(float u1, float u2);
00105 float UniformHemispherePdf();
00106 Vector UniformSampleSphere(float u1, float u2);
00107 float UniformSpherePdf();
00108 Vector UniformSampleCone(float u1, float u2, float thetamax);
00109 Vector UniformSampleCone(float u1, float u2, float thetamax,
00110 const Vector &x, const Vector &y, const Vector &z);
00111 float UniformConePdf(float thetamax);
00112 void UniformSampleDisk(float u1, float u2, float *x, float *y);
00113 void ConcentricSampleDisk(float u1, float u2, float *dx, float *dy);
00114 inline Vector CosineSampleHemisphere(float u1, float u2) {
00115 Vector ret;
00116 ConcentricSampleDisk(u1, u2, &ret.x, &ret.y);
00117 ret.z = sqrtf(max(0.f, 1.f - ret.x*ret.x - ret.y*ret.y));
00118 return ret;
00119 }
00120
00121
00122 inline float CosineHemispherePdf(float costheta, float phi) {
00123 return costheta * INV_PI;
00124 }
00125
00126
00127 void UniformSampleTriangle(float ud1, float ud2, float *u, float *v);
00128 struct Distribution2D {
00129
00130 Distribution2D(const float *data, int nu, int nv);
00131 ~Distribution2D();
00132 void SampleContinuous(float u0, float u1, float uv[2],
00133 float *pdf) const {
00134 float pdfs[2];
00135 int v;
00136 uv[1] = pMarginal->SampleContinuous(u1, &pdfs[1], &v);
00137 uv[0] = pConditionalV[v]->SampleContinuous(u0, &pdfs[0]);
00138 *pdf = pdfs[0] * pdfs[1];
00139 }
00140 float Pdf(float u, float v) const {
00141 int iu = Clamp(Float2Int(u * pConditionalV[0]->count), 0,
00142 pConditionalV[0]->count-1);
00143 int iv = Clamp(Float2Int(v * pMarginal->count), 0,
00144 pMarginal->count-1);
00145 if (pConditionalV[iv]->funcInt * pMarginal->funcInt == 0.f) return 0.f;
00146 return (pConditionalV[iv]->func[iu] * pMarginal->func[iv]) /
00147 (pConditionalV[iv]->funcInt * pMarginal->funcInt);
00148 }
00149 private:
00150
00151 vector<Distribution1D *> pConditionalV;
00152 Distribution1D *pMarginal;
00153 };
00154
00155
00156 void StratifiedSample1D(float *samples, int nsamples, RNG &rng,
00157 bool jitter = true);
00158 void StratifiedSample2D(float *samples, int nx, int ny, RNG &rng,
00159 bool jitter = true);
00160 template <typename T>
00161 void Shuffle(T *samp, uint32_t count, uint32_t dims, RNG &rng) {
00162 for (uint32_t i = 0; i < count; ++i) {
00163 uint32_t other = i + (rng.RandomUInt() % (count - i));
00164 for (uint32_t j = 0; j < dims; ++j)
00165 swap(samp[dims*i + j], samp[dims*other + j]);
00166 }
00167 }
00168
00169
00170 void LatinHypercube(float *samples, uint32_t nSamples, uint32_t nDim, RNG &rng);
00171 inline double RadicalInverse(int n, int base) {
00172 double val = 0;
00173 double invBase = 1. / base, invBi = invBase;
00174 while (n > 0) {
00175
00176 int d_i = (n % base);
00177 val += d_i * invBi;
00178 n *= invBase;
00179 invBi *= invBase;
00180 }
00181 return val;
00182 }
00183
00184
00185 inline void GeneratePermutation(uint32_t *buf, uint32_t b, RNG &rng) {
00186 for (uint32_t i = 0; i < b; ++i)
00187 buf[i] = i;
00188 Shuffle(buf, b, 1, rng);
00189 }
00190
00191
00192 inline double PermutedRadicalInverse(uint32_t n, uint32_t base,
00193 const uint32_t *p) {
00194 double val = 0;
00195 double invBase = 1. / base, invBi = invBase;
00196
00197 while (n > 0) {
00198 uint32_t d_i = p[n % base];
00199 val += d_i * invBi;
00200 n *= invBase;
00201 invBi *= invBase;
00202 }
00203 return val;
00204 }
00205
00206
00207 class PermutedHalton {
00208 public:
00209
00210 PermutedHalton(uint32_t d, RNG &rng);
00211 ~PermutedHalton() {
00212 delete[] b;
00213 delete[] permute;
00214 }
00215 void Sample(uint32_t n, float *out) const {
00216 uint32_t *p = permute;
00217 for (uint32_t i = 0; i < dims; ++i) {
00218 out[i] = min(float(PermutedRadicalInverse(n, b[i], p)),
00219 OneMinusEpsilon);
00220 p += b[i];
00221 }
00222 }
00223 private:
00224
00225 uint32_t dims;
00226 uint32_t *b, *permute;
00227 PermutedHalton(const PermutedHalton &);
00228 PermutedHalton &operator=(const PermutedHalton &);
00229 };
00230
00231
00232 inline float VanDerCorput(uint32_t n, uint32_t scramble = 0);
00233 inline float Sobol2(uint32_t n, uint32_t scramble = 0);
00234 inline float LarcherPillichshammer2(uint32_t n, uint32_t scramble = 0);
00235 inline void Sample02(uint32_t n, const uint32_t scramble[2], float sample[2]);
00236 int LDPixelSampleFloatsNeeded(const Sample *sample, int nPixelSamples);
00237 void LDPixelSample(int xPos, int yPos, float shutterOpen,
00238 float shutterClose, int nPixelSamples, Sample *samples, float *buf, RNG &rng);
00239 Vector SampleHG(const Vector &w, float g, float u1, float u2);
00240 float HGPdf(const Vector &w, const Vector &wp, float g);
00241
00242
00243 inline float BalanceHeuristic(int nf, float fPdf, int ng, float gPdf) {
00244 return (nf * fPdf) / (nf * fPdf + ng * gPdf);
00245 }
00246
00247
00248 inline float PowerHeuristic(int nf, float fPdf, int ng, float gPdf) {
00249 float f = nf * fPdf, g = ng * gPdf;
00250 return (f*f) / (f*f + g*g);
00251 }
00252
00253
00254
00255
00256 inline void Sample02(uint32_t n, const uint32_t scramble[2],
00257 float sample[2]) {
00258 sample[0] = VanDerCorput(n, scramble[0]);
00259 sample[1] = Sobol2(n, scramble[1]);
00260 }
00261
00262
00263 inline float VanDerCorput(uint32_t n, uint32_t scramble) {
00264
00265 n = (n << 16) | (n >> 16);
00266 n = ((n & 0x00ff00ff) << 8) | ((n & 0xff00ff00) >> 8);
00267 n = ((n & 0x0f0f0f0f) << 4) | ((n & 0xf0f0f0f0) >> 4);
00268 n = ((n & 0x33333333) << 2) | ((n & 0xcccccccc) >> 2);
00269 n = ((n & 0x55555555) << 1) | ((n & 0xaaaaaaaa) >> 1);
00270 n ^= scramble;
00271 return min(((n>>8) & 0xffffff) / float(1 << 24), OneMinusEpsilon);
00272 }
00273
00274
00275 inline float Sobol2(uint32_t n, uint32_t scramble) {
00276 for (uint32_t v = 1 << 31; n != 0; n >>= 1, v ^= v >> 1)
00277 if (n & 0x1) scramble ^= v;
00278 return min(((scramble>>8) & 0xffffff) / float(1 << 24), OneMinusEpsilon);
00279 }
00280
00281
00282 inline float
00283 LarcherPillichshammer2(uint32_t n, uint32_t scramble) {
00284 for (uint32_t v = 1 << 31; n != 0; n >>= 1, v |= v >> 1)
00285 if (n & 0x1) scramble ^= v;
00286 return min(((scramble>>8) & 0xffffff) / float(1 << 24), OneMinusEpsilon);
00287 }
00288
00289
00290 inline void LDShuffleScrambled1D(int nSamples, int nPixel,
00291 float *samples, RNG &rng) {
00292 uint32_t scramble = rng.RandomUInt();
00293 for (int i = 0; i < nSamples * nPixel; ++i)
00294 samples[i] = VanDerCorput(i, scramble);
00295 for (int i = 0; i < nPixel; ++i)
00296 Shuffle(samples + i * nSamples, nSamples, 1, rng);
00297 Shuffle(samples, nPixel, nSamples, rng);
00298 }
00299
00300
00301 inline void LDShuffleScrambled2D(int nSamples, int nPixel,
00302 float *samples, RNG &rng) {
00303 uint32_t scramble[2] = { rng.RandomUInt(), rng.RandomUInt() };
00304 for (int i = 0; i < nSamples * nPixel; ++i)
00305 Sample02(i, scramble, &samples[2*i]);
00306 for (int i = 0; i < nPixel; ++i)
00307 Shuffle(samples + 2 * i * nSamples, nSamples, 2, rng);
00308 Shuffle(samples, nPixel, 2 * nSamples, rng);
00309 }
00310
00311
00312
00313 #endif // PBRT_CORE_MONTECARLO_H