PBRT
/home/felix/UBC/projects/AdaptiveLightfieldSampling/pbrt_v2/src/core/montecarlo.h
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