PBRT
|
00001 00002 /* 00003 pbrt source code Copyright(c) 1998-2012 Matt Pharr and Greg Humphreys. 00004 00005 This file is part of pbrt. 00006 00007 Redistribution and use in source and binary forms, with or without 00008 modification, are permitted provided that the following conditions are 00009 met: 00010 00011 - Redistributions of source code must retain the above copyright 00012 notice, this list of conditions and the following disclaimer. 00013 00014 - Redistributions in binary form must reproduce the above copyright 00015 notice, this list of conditions and the following disclaimer in the 00016 documentation and/or other materials provided with the distribution. 00017 00018 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS 00019 IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED 00020 TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A 00021 PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT 00022 HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, 00023 SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT 00024 LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, 00025 DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY 00026 THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT 00027 (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE 00028 OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 00029 00030 */ 00031 00032 #if defined(_MSC_VER) 00033 #pragma once 00034 #endif 00035 00036 #ifndef PBRT_CORE_MONTECARLO_H 00037 #define PBRT_CORE_MONTECARLO_H 00038 00039 // core/montecarlo.h* 00040 #include "pbrt.h" 00041 #include "geometry.h" 00042 #include "rng.h" 00043 00044 // smallest floating point value less than one; all canonical random samples 00045 // should be <= this. 00046 #ifdef PBRT_IS_WINDOWS 00047 // sadly, MSVC2008 (at least) doesn't support hexidecimal fp constants... 00048 static const float OneMinusEpsilon=0.9999999403953552f; 00049 #else 00050 static const float OneMinusEpsilon=0x1.fffffep-1; 00051 #endif 00052 00053 // Monte Carlo Utility Declarations 00054 struct Distribution1D { 00055 // Distribution1D Public Methods 00056 Distribution1D(const float *f, int n) { 00057 count = n; 00058 func = new float[n]; 00059 memcpy(func, f, n*sizeof(float)); 00060 cdf = new float[n+1]; 00061 // Compute integral of step function at $x_i$ 00062 cdf[0] = 0.; 00063 for (int i = 1; i < count+1; ++i) 00064 cdf[i] = cdf[i-1] + func[i-1] / n; 00065 00066 // Transform step function integral into CDF 00067 funcInt = cdf[count]; 00068 if (funcInt == 0.f) { 00069 for (int i = 1; i < n+1; ++i) 00070 cdf[i] = float(i) / float(n); 00071 } 00072 else { 00073 for (int i = 1; i < n+1; ++i) 00074 cdf[i] /= funcInt; 00075 } 00076 } 00077 ~Distribution1D() { 00078 delete[] func; 00079 delete[] cdf; 00080 } 00081 float SampleContinuous(float u, float *pdf, int *off = NULL) const { 00082 // Find surrounding CDF segments and _offset_ 00083 float *ptr = std::upper_bound(cdf, cdf+count+1, u); 00084 int offset = max(0, int(ptr-cdf-1)); 00085 if (off) *off = offset; 00086 Assert(offset < count); 00087 Assert(u >= cdf[offset] && u < cdf[offset+1]); 00088 00089 // Compute offset along CDF segment 00090 float du = (u - cdf[offset]) / (cdf[offset+1] - cdf[offset]); 00091 Assert(!isnan(du)); 00092 00093 // Compute PDF for sampled offset 00094 if (pdf) *pdf = func[offset] / funcInt; 00095 00096 // Return $x\in{}[0,1)$ corresponding to sample 00097 return (offset + du) / count; 00098 } 00099 int SampleDiscrete(float u, float *pdf) const { 00100 // Find surrounding CDF segments and _offset_ 00101 float *ptr = std::upper_bound(cdf, cdf+count+1, u); 00102 int offset = max(0, int(ptr-cdf-1)); 00103 Assert(offset < count); 00104 Assert(u >= cdf[offset] && u < cdf[offset+1]); 00105 if (pdf) *pdf = func[offset] / (funcInt * count); 00106 return offset; 00107 } 00108 private: 00109 friend struct Distribution2D; 00110 // Distribution1D Private Data 00111 float *func, *cdf; 00112 float funcInt; 00113 int count; 00114 }; 00115 00116 00117 void RejectionSampleDisk(float *x, float *y, RNG &rng); 00118 Vector UniformSampleHemisphere(float u1, float u2); 00119 float UniformHemispherePdf(); 00120 Vector UniformSampleSphere(float u1, float u2); 00121 float UniformSpherePdf(); 00122 Vector UniformSampleCone(float u1, float u2, float thetamax); 00123 Vector UniformSampleCone(float u1, float u2, float thetamax, 00124 const Vector &x, const Vector &y, const Vector &z); 00125 float UniformConePdf(float thetamax); 00126 void UniformSampleDisk(float u1, float u2, float *x, float *y); 00127 void ConcentricSampleDisk(float u1, float u2, float *dx, float *dy); 00128 inline Vector CosineSampleHemisphere(float u1, float u2) { 00129 Vector ret; 00130 ConcentricSampleDisk(u1, u2, &ret.x, &ret.y); 00131 ret.z = sqrtf(max(0.f, 1.f - ret.x*ret.x - ret.y*ret.y)); 00132 return ret; 00133 } 00134 00135 00136 inline float CosineHemispherePdf(float costheta, float phi) { 00137 return costheta * INV_PI; 00138 } 00139 00140 00141 void UniformSampleTriangle(float ud1, float ud2, float *u, float *v); 00142 struct Distribution2D { 00143 // Distribution2D Public Methods 00144 Distribution2D(const float *data, int nu, int nv); 00145 ~Distribution2D(); 00146 void SampleContinuous(float u0, float u1, float uv[2], 00147 float *pdf) const { 00148 float pdfs[2]; 00149 int v; 00150 uv[1] = pMarginal->SampleContinuous(u1, &pdfs[1], &v); 00151 uv[0] = pConditionalV[v]->SampleContinuous(u0, &pdfs[0]); 00152 *pdf = pdfs[0] * pdfs[1]; 00153 } 00154 float Pdf(float u, float v) const { 00155 int iu = Clamp(Float2Int(u * pConditionalV[0]->count), 0, 00156 pConditionalV[0]->count-1); 00157 int iv = Clamp(Float2Int(v * pMarginal->count), 0, 00158 pMarginal->count-1); 00159 if (pConditionalV[iv]->funcInt * pMarginal->funcInt == 0.f) return 0.f; 00160 return (pConditionalV[iv]->func[iu] * pMarginal->func[iv]) / 00161 (pConditionalV[iv]->funcInt * pMarginal->funcInt); 00162 } 00163 private: 00164 // Distribution2D Private Data 00165 vector<Distribution1D *> pConditionalV; 00166 Distribution1D *pMarginal; 00167 }; 00168 00169 00170 void StratifiedSample1D(float *samples, int nsamples, RNG &rng, 00171 bool jitter = true); 00172 void StratifiedSample2D(float *samples, int nx, int ny, RNG &rng, 00173 bool jitter = true); 00174 template <typename T> 00175 void Shuffle(T *samp, uint32_t count, uint32_t dims, RNG &rng) { 00176 for (uint32_t i = 0; i < count; ++i) { 00177 uint32_t other = i + (rng.RandomUInt() % (count - i)); 00178 for (uint32_t j = 0; j < dims; ++j) 00179 swap(samp[dims*i + j], samp[dims*other + j]); 00180 } 00181 } 00182 00183 00184 void LatinHypercube(float *samples, uint32_t nSamples, uint32_t nDim, RNG &rng); 00185 inline double RadicalInverse(int n, int base) { 00186 double val = 0; 00187 double invBase = 1. / base, invBi = invBase; 00188 while (n > 0) { 00189 // Compute next digit of radical inverse 00190 int d_i = (n % base); 00191 val += d_i * invBi; 00192 n *= invBase; 00193 invBi *= invBase; 00194 } 00195 return val; 00196 } 00197 00198 00199 inline void GeneratePermutation(uint32_t *buf, uint32_t b, RNG &rng) { 00200 for (uint32_t i = 0; i < b; ++i) 00201 buf[i] = i; 00202 Shuffle(buf, b, 1, rng); 00203 } 00204 00205 00206 inline double PermutedRadicalInverse(uint32_t n, uint32_t base, 00207 const uint32_t *p) { 00208 double val = 0; 00209 double invBase = 1. / base, invBi = invBase; 00210 00211 while (n > 0) { 00212 uint32_t d_i = p[n % base]; 00213 val += d_i * invBi; 00214 n *= invBase; 00215 invBi *= invBase; 00216 } 00217 return val; 00218 } 00219 00220 00221 class PermutedHalton { 00222 public: 00223 // PermutedHalton Public Methods 00224 PermutedHalton(uint32_t d, RNG &rng); 00225 ~PermutedHalton() { 00226 delete[] b; 00227 delete[] permute; 00228 } 00229 void Sample(uint32_t n, float *out) const { 00230 uint32_t *p = permute; 00231 for (uint32_t i = 0; i < dims; ++i) { 00232 out[i] = min(float(PermutedRadicalInverse(n, b[i], p)), 00233 OneMinusEpsilon); 00234 p += b[i]; 00235 } 00236 } 00237 private: 00238 // PermutedHalton Private Data 00239 uint32_t dims; 00240 uint32_t *b, *permute; 00241 PermutedHalton(const PermutedHalton &); 00242 PermutedHalton &operator=(const PermutedHalton &); 00243 }; 00244 00245 00246 inline float VanDerCorput(uint32_t n, uint32_t scramble = 0); 00247 inline float Sobol2(uint32_t n, uint32_t scramble = 0); 00248 inline float LarcherPillichshammer2(uint32_t n, uint32_t scramble = 0); 00249 inline void Sample02(uint32_t n, const uint32_t scramble[2], float sample[2]); 00250 int LDPixelSampleFloatsNeeded(const Sample *sample, int nPixelSamples); 00251 void LDPixelSample(int xPos, int yPos, float shutterOpen, 00252 float shutterClose, int nPixelSamples, Sample *samples, float *buf, RNG &rng); 00253 Vector SampleHG(const Vector &w, float g, float u1, float u2); 00254 float HGPdf(const Vector &w, const Vector &wp, float g); 00255 00256 // Monte Carlo Inline Functions 00257 inline float BalanceHeuristic(int nf, float fPdf, int ng, float gPdf) { 00258 return (nf * fPdf) / (nf * fPdf + ng * gPdf); 00259 } 00260 00261 00262 inline float PowerHeuristic(int nf, float fPdf, int ng, float gPdf) { 00263 float f = nf * fPdf, g = ng * gPdf; 00264 return (f*f) / (f*f + g*g); 00265 } 00266 00267 00268 00269 // Sampling Inline Functions 00270 inline void Sample02(uint32_t n, const uint32_t scramble[2], 00271 float sample[2]) { 00272 sample[0] = VanDerCorput(n, scramble[0]); 00273 sample[1] = Sobol2(n, scramble[1]); 00274 } 00275 00276 00277 inline float VanDerCorput(uint32_t n, uint32_t scramble) { 00278 // Reverse bits of _n_ 00279 n = (n << 16) | (n >> 16); 00280 n = ((n & 0x00ff00ff) << 8) | ((n & 0xff00ff00) >> 8); 00281 n = ((n & 0x0f0f0f0f) << 4) | ((n & 0xf0f0f0f0) >> 4); 00282 n = ((n & 0x33333333) << 2) | ((n & 0xcccccccc) >> 2); 00283 n = ((n & 0x55555555) << 1) | ((n & 0xaaaaaaaa) >> 1); 00284 n ^= scramble; 00285 return min(((n>>8) & 0xffffff) / float(1 << 24), OneMinusEpsilon); 00286 } 00287 00288 00289 inline float Sobol2(uint32_t n, uint32_t scramble) { 00290 for (uint32_t v = 1 << 31; n != 0; n >>= 1, v ^= v >> 1) 00291 if (n & 0x1) scramble ^= v; 00292 return min(((scramble>>8) & 0xffffff) / float(1 << 24), OneMinusEpsilon); 00293 } 00294 00295 00296 inline float 00297 LarcherPillichshammer2(uint32_t n, uint32_t scramble) { 00298 for (uint32_t v = 1 << 31; n != 0; n >>= 1, v |= v >> 1) 00299 if (n & 0x1) scramble ^= v; 00300 return min(((scramble>>8) & 0xffffff) / float(1 << 24), OneMinusEpsilon); 00301 } 00302 00303 00304 inline void LDShuffleScrambled1D(int nSamples, int nPixel, 00305 float *samples, RNG &rng) { 00306 uint32_t scramble = rng.RandomUInt(); 00307 for (int i = 0; i < nSamples * nPixel; ++i) 00308 samples[i] = VanDerCorput(i, scramble); 00309 for (int i = 0; i < nPixel; ++i) 00310 Shuffle(samples + i * nSamples, nSamples, 1, rng); 00311 Shuffle(samples, nPixel, nSamples, rng); 00312 } 00313 00314 00315 inline void LDShuffleScrambled2D(int nSamples, int nPixel, 00316 float *samples, RNG &rng) { 00317 uint32_t scramble[2] = { rng.RandomUInt(), rng.RandomUInt() }; 00318 for (int i = 0; i < nSamples * nPixel; ++i) 00319 Sample02(i, scramble, &samples[2*i]); 00320 for (int i = 0; i < nPixel; ++i) 00321 Shuffle(samples + 2 * i * nSamples, nSamples, 2, rng); 00322 Shuffle(samples, nPixel, 2 * nSamples, rng); 00323 } 00324 00325 00326 00327 #endif // PBRT_CORE_MONTECARLO_H