00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026 #include "stdafx.h"
00027 #include "montecarlo.h"
00028 #include "geometry.h"
00029 #include "shape.h"
00030 #include "volume.h"
00031
00032
00033 static const int primes[] = {
00034
00035 2, 3, 5, 7, 11, 13, 17, 19, 23, 29,
00036 31, 37, 41, 43, 47, 53, 59, 61, 67, 71,
00037 73, 79, 83, 89, 97, 101, 103, 107, 109, 113,
00038 127, 131, 137, 139, 149, 151, 157, 163, 167, 173,
00039 179, 181, 191, 193, 197, 199, 211, 223, 227, 229,
00040 233, 239, 241, 251, 257, 263, 269, 271, 277, 281,
00041 283, 293, 307, 311, 313, 317, 331, 337, 347, 349,
00042 353, 359, 367, 373, 379, 383, 389, 397, 401, 409,
00043 419, 421, 431, 433, 439, 443, 449, 457, 461, 463,
00044 467, 479, 487, 491, 499, 503, 509, 521, 523, 541,
00045 547, 557, 563, 569, 571, 577, 587, 593, 599, 601,
00046 607, 613, 617, 619, 631, 641, 643, 647, 653, 659,
00047 661, 673, 677, 683, 691, 701, 709, 719, 727, 733,
00048 739, 743, 751, 757, 761, 769, 773, 787, 797, 809,
00049 811, 821, 823, 827, 829, 839, 853, 857, 859, 863,
00050 877, 881, 883, 887, 907, 911, 919, 929, 937, 941,
00051 947, 953, 967, 971, 977, 983, 991, 997, 1009, 1013,
00052 1019, 1021, 1031, 1033, 1039, 1049, 1051, 1061, 1063, 1069,
00053 1087, 1091, 1093, 1097, 1103, 1109, 1117, 1123, 1129, 1151,
00054 1153, 1163, 1171, 1181, 1187, 1193, 1201, 1213, 1217, 1223,
00055 1229, 1231, 1237, 1249, 1259, 1277, 1279, 1283, 1289, 1291,
00056 1297, 1301, 1303, 1307, 1319, 1321, 1327, 1361, 1367, 1373,
00057 1381, 1399, 1409, 1423, 1427, 1429, 1433, 1439, 1447, 1451,
00058 1453, 1459, 1471, 1481, 1483, 1487, 1489, 1493, 1499, 1511,
00059 1523, 1531, 1543, 1549, 1553, 1559, 1567, 1571, 1579, 1583,
00060 1597, 1601, 1607, 1609, 1613, 1619, 1621, 1627, 1637, 1657,
00061 1663, 1667, 1669, 1693, 1697, 1699, 1709, 1721, 1723, 1733,
00062 1741, 1747, 1753, 1759, 1777, 1783, 1787, 1789, 1801, 1811,
00063 1823, 1831, 1847, 1861, 1867, 1871, 1873, 1877, 1879, 1889,
00064 1901, 1907, 1913, 1931, 1933, 1949, 1951, 1973, 1979, 1987,
00065 1993, 1997, 1999, 2003, 2011, 2017, 2027, 2029, 2039, 2053,
00066 2063, 2069, 2081, 2083, 2087, 2089, 2099, 2111, 2113, 2129,
00067 2131, 2137, 2141, 2143, 2153, 2161, 2179, 2203, 2207, 2213,
00068 2221, 2237, 2239, 2243, 2251, 2267, 2269, 2273, 2281, 2287,
00069 2293, 2297, 2309, 2311, 2333, 2339, 2341, 2347, 2351, 2357,
00070 2371, 2377, 2381, 2383, 2389, 2393, 2399, 2411, 2417, 2423,
00071 2437, 2441, 2447, 2459, 2467, 2473, 2477, 2503, 2521, 2531,
00072 2539, 2543, 2549, 2551, 2557, 2579, 2591, 2593, 2609, 2617,
00073 2621, 2633, 2647, 2657, 2659, 2663, 2671, 2677, 2683, 2687,
00074 2689, 2693, 2699, 2707, 2711, 2713, 2719, 2729, 2731, 2741,
00075 2749, 2753, 2767, 2777, 2789, 2791, 2797, 2801, 2803, 2819,
00076 2833, 2837, 2843, 2851, 2857, 2861, 2879, 2887, 2897, 2903,
00077 2909, 2917, 2927, 2939, 2953, 2957, 2963, 2969, 2971, 2999,
00078 3001, 3011, 3019, 3023, 3037, 3041, 3049, 3061, 3067, 3079,
00079 3083, 3089, 3109, 3119, 3121, 3137, 3163, 3167, 3169, 3181,
00080 3187, 3191, 3203, 3209, 3217, 3221, 3229, 3251, 3253, 3257,
00081 3259, 3271, 3299, 3301, 3307, 3313, 3319, 3323, 3329, 3331,
00082 3343, 3347, 3359, 3361, 3371, 3373, 3389, 3391, 3407, 3413,
00083 3433, 3449, 3457, 3461, 3463, 3467, 3469, 3491, 3499, 3511,
00084 3517, 3527, 3529, 3533, 3539, 3541, 3547, 3557, 3559, 3571,
00085 3581, 3583, 3593, 3607, 3613, 3617, 3623, 3631, 3637, 3643,
00086 3659, 3671, 3673, 3677, 3691, 3697, 3701, 3709, 3719, 3727,
00087 3733, 3739, 3761, 3767, 3769, 3779, 3793, 3797, 3803, 3821,
00088 3823, 3833, 3847, 3851, 3853, 3863, 3877, 3881, 3889, 3907,
00089 3911, 3917, 3919, 3923, 3929, 3931, 3943, 3947, 3967, 3989,
00090 4001, 4003, 4007, 4013, 4019, 4021, 4027, 4049, 4051, 4057,
00091 4073, 4079, 4091, 4093, 4099, 4111, 4127, 4129, 4133, 4139,
00092 4153, 4157, 4159, 4177, 4201, 4211, 4217, 4219, 4229, 4231,
00093 4241, 4243, 4253, 4259, 4261, 4271, 4273, 4283, 4289, 4297,
00094 4327, 4337, 4339, 4349, 4357, 4363, 4373, 4391, 4397, 4409,
00095 4421, 4423, 4441, 4447, 4451, 4457, 4463, 4481, 4483, 4493,
00096 4507, 4513, 4517, 4519, 4523, 4547, 4549, 4561, 4567, 4583,
00097 4591, 4597, 4603, 4621, 4637, 4639, 4643, 4649, 4651, 4657,
00098 4663, 4673, 4679, 4691, 4703, 4721, 4723, 4729, 4733, 4751,
00099 4759, 4783, 4787, 4789, 4793, 4799, 4801, 4813, 4817, 4831,
00100 4861, 4871, 4877, 4889, 4903, 4909, 4919, 4931, 4933, 4937,
00101 4943, 4951, 4957, 4967, 4969, 4973, 4987, 4993, 4999, 5003,
00102 5009, 5011, 5021, 5023, 5039, 5051, 5059, 5077, 5081, 5087,
00103 5099, 5101, 5107, 5113, 5119, 5147, 5153, 5167, 5171, 5179,
00104 5189, 5197, 5209, 5227, 5231, 5233, 5237, 5261, 5273, 5279,
00105 5281, 5297, 5303, 5309, 5323, 5333, 5347, 5351, 5381, 5387,
00106 5393, 5399, 5407, 5413, 5417, 5419, 5431, 5437, 5441, 5443,
00107 5449, 5471, 5477, 5479, 5483, 5501, 5503, 5507, 5519, 5521,
00108 5527, 5531, 5557, 5563, 5569, 5573, 5581, 5591, 5623, 5639,
00109 5641, 5647, 5651, 5653, 5657, 5659, 5669, 5683, 5689, 5693,
00110 5701, 5711, 5717, 5737, 5741, 5743, 5749, 5779, 5783, 5791,
00111 5801, 5807, 5813, 5821, 5827, 5839, 5843, 5849, 5851, 5857,
00112 5861, 5867, 5869, 5879, 5881, 5897, 5903, 5923, 5927, 5939,
00113 5953, 5981, 5987, 6007, 6011, 6029, 6037, 6043, 6047, 6053,
00114 6067, 6073, 6079, 6089, 6091, 6101, 6113, 6121, 6131, 6133,
00115 6143, 6151, 6163, 6173, 6197, 6199, 6203, 6211, 6217, 6221,
00116 6229, 6247, 6257, 6263, 6269, 6271, 6277, 6287, 6299, 6301,
00117 6311, 6317, 6323, 6329, 6337, 6343, 6353, 6359, 6361, 6367,
00118 6373, 6379, 6389, 6397, 6421, 6427, 6449, 6451, 6469, 6473,
00119 6481, 6491, 6521, 6529, 6547, 6551, 6553, 6563, 6569, 6571,
00120 6577, 6581, 6599, 6607, 6619, 6637, 6653, 6659, 6661, 6673,
00121 6679, 6689, 6691, 6701, 6703, 6709, 6719, 6733, 6737, 6761,
00122 6763, 6779, 6781, 6791, 6793, 6803, 6823, 6827, 6829, 6833,
00123 6841, 6857, 6863, 6869, 6871, 6883, 6899, 6907, 6911, 6917,
00124 6947, 6949, 6959, 6961, 6967, 6971, 6977, 6983, 6991, 6997,
00125 7001, 7013, 7019, 7027, 7039, 7043, 7057, 7069, 7079, 7103,
00126 7109, 7121, 7127, 7129, 7151, 7159, 7177, 7187, 7193, 7207,
00127 7211, 7213, 7219, 7229, 7237, 7243, 7247, 7253, 7283, 7297,
00128 7307, 7309, 7321, 7331, 7333, 7349, 7351, 7369, 7393, 7411,
00129 7417, 7433, 7451, 7457, 7459, 7477, 7481, 7487, 7489, 7499,
00130 7507, 7517, 7523, 7529, 7537, 7541, 7547, 7549, 7559, 7561,
00131 7573, 7577, 7583, 7589, 7591, 7603, 7607, 7621, 7639, 7643,
00132 7649, 7669, 7673, 7681, 7687, 7691, 7699, 7703, 7717, 7723,
00133 7727, 7741, 7753, 7757, 7759, 7789, 7793, 7817, 7823, 7829,
00134 7841, 7853, 7867, 7873, 7877, 7879, 7883, 7901, 7907, 7919
00135 };
00136
00137
00138
00139
00140 void StratifiedSample1D(float *samp, int nSamples, RNG &rng,
00141 bool jitter) {
00142 float invTot = 1.f / nSamples;
00143 for (int i = 0; i < nSamples; ++i) {
00144 float delta = jitter ? rng.RandomFloat() : 0.5f;
00145 *samp++ = min((i + delta) * invTot, OneMinusEpsilon);
00146 }
00147 }
00148
00149
00150 void StratifiedSample2D(float *samp, int nx, int ny, RNG &rng,
00151 bool jitter) {
00152 float dx = 1.f / nx, dy = 1.f / ny;
00153 for (int y = 0; y < ny; ++y)
00154 for (int x = 0; x < nx; ++x) {
00155 float jx = jitter ? rng.RandomFloat() : 0.5f;
00156 float jy = jitter ? rng.RandomFloat() : 0.5f;
00157 *samp++ = min((x + jx) * dx, OneMinusEpsilon);
00158 *samp++ = min((y + jy) * dy, OneMinusEpsilon);
00159 }
00160 }
00161
00162
00163 void LatinHypercube(float *samples, uint32_t nSamples, uint32_t nDim,
00164 RNG &rng) {
00165
00166 float delta = 1.f / nSamples;
00167 for (uint32_t i = 0; i < nSamples; ++i)
00168 for (uint32_t j = 0; j < nDim; ++j)
00169 samples[nDim * i + j] = min((i + (rng.RandomFloat())) * delta,
00170 OneMinusEpsilon);
00171
00172
00173 for (uint32_t i = 0; i < nDim; ++i) {
00174 for (uint32_t j = 0; j < nSamples; ++j) {
00175 uint32_t other = j + (rng.RandomUInt() % (nSamples - j));
00176 swap(samples[nDim * j + i], samples[nDim * other + i]);
00177 }
00178 }
00179 }
00180
00181
00182 int LDPixelSampleFloatsNeeded(const Sample *sample, int nPixelSamples) {
00183 int n = 5;
00184 for (uint32_t i = 0; i < sample->n1D.size(); ++i)
00185 n += sample->n1D[i];
00186 for (uint32_t i = 0; i < sample->n2D.size(); ++i)
00187 n += 2 * sample->n2D[i];
00188 return nPixelSamples * n;
00189 }
00190
00191
00192 void LDPixelSample(int xPos, int yPos, float shutterOpen,
00193 float shutterClose, int nPixelSamples, Sample *samples,
00194 float *buf, RNG &rng) {
00195
00196 float *imageSamples = buf; buf += 2 * nPixelSamples;
00197 float *lensSamples = buf; buf += 2 * nPixelSamples;
00198 float *timeSamples = buf; buf += nPixelSamples;
00199
00200
00201 uint32_t count1D = samples[0].n1D.size();
00202 uint32_t count2D = samples[0].n2D.size();
00203 const uint32_t *n1D = count1D > 0 ? &samples[0].n1D[0] : NULL;
00204 const uint32_t *n2D = count2D > 0 ? &samples[0].n2D[0] : NULL;
00205 float **oneDSamples = ALLOCA(float *, count1D);
00206 float **twoDSamples = ALLOCA(float *, count2D);
00207 for (uint32_t i = 0; i < count1D; ++i) {
00208 oneDSamples[i] = buf;
00209 buf += n1D[i] * nPixelSamples;
00210 }
00211 for (uint32_t i = 0; i < count2D; ++i) {
00212 twoDSamples[i] = buf;
00213 buf += 2 * n2D[i] * nPixelSamples;
00214 }
00215
00216
00217 LDShuffleScrambled2D(1, nPixelSamples, imageSamples, rng);
00218 LDShuffleScrambled2D(1, nPixelSamples, lensSamples, rng);
00219 LDShuffleScrambled1D(1, nPixelSamples, timeSamples, rng);
00220 for (uint32_t i = 0; i < count1D; ++i)
00221 LDShuffleScrambled1D(n1D[i], nPixelSamples, oneDSamples[i], rng);
00222 for (uint32_t i = 0; i < count2D; ++i)
00223 LDShuffleScrambled2D(n2D[i], nPixelSamples, twoDSamples[i], rng);
00224
00225
00226 for (int i = 0; i < nPixelSamples; ++i) {
00227 samples[i].imageX = xPos + imageSamples[2*i];
00228 samples[i].imageY = yPos + imageSamples[2*i+1];
00229 samples[i].time = Lerp(timeSamples[i], shutterOpen, shutterClose);
00230 samples[i].lensU = lensSamples[2*i];
00231 samples[i].lensV = lensSamples[2*i+1];
00232
00233 for (uint32_t j = 0; j < count1D; ++j) {
00234 int startSamp = n1D[j] * i;
00235 for (uint32_t k = 0; k < n1D[j]; ++k)
00236 samples[i].oneD[j][k] = oneDSamples[j][startSamp+k];
00237 }
00238 for (uint32_t j = 0; j < count2D; ++j) {
00239 int startSamp = 2 * n2D[j] * i;
00240 for (uint32_t k = 0; k < 2*n2D[j]; ++k)
00241 samples[i].twoD[j][k] = twoDSamples[j][startSamp+k];
00242 }
00243 }
00244 }
00245
00246
00247
00248
00249 void RejectionSampleDisk(float *x, float *y, RNG &rng) {
00250 float sx, sy;
00251 do {
00252 sx = 1.f - 2.f * rng.RandomFloat();
00253 sy = 1.f - 2.f * rng.RandomFloat();
00254 } while (sx*sx + sy*sy > 1.f);
00255 *x = sx;
00256 *y = sy;
00257 }
00258
00259
00260 Vector UniformSampleHemisphere(float u1, float u2) {
00261 float z = u1;
00262 float r = sqrtf(max(0.f, 1.f - z*z));
00263 float phi = 2 * M_PI * u2;
00264 float x = r * cosf(phi);
00265 float y = r * sinf(phi);
00266 return Vector(x, y, z);
00267 }
00268
00269
00270 float UniformHemispherePdf() {
00271 return INV_TWOPI;
00272 }
00273
00274
00275 Vector UniformSampleSphere(float u1, float u2) {
00276 float z = 1.f - 2.f * u1;
00277 float r = sqrtf(max(0.f, 1.f - z*z));
00278 float phi = 2.f * M_PI * u2;
00279 float x = r * cosf(phi);
00280 float y = r * sinf(phi);
00281 return Vector(x, y, z);
00282 }
00283
00284
00285 float UniformSpherePdf() {
00286 return 1.f / (4.f * M_PI);
00287 }
00288
00289
00290 void UniformSampleDisk(float u1, float u2, float *x, float *y) {
00291 float r = sqrtf(u1);
00292 float theta = 2.0f * M_PI * u2;
00293 *x = r * cosf(theta);
00294 *y = r * sinf(theta);
00295 }
00296
00297
00298 void ConcentricSampleDisk(float u1, float u2, float *dx, float *dy) {
00299 float r, theta;
00300
00301 float sx = 2 * u1 - 1;
00302 float sy = 2 * u2 - 1;
00303
00304
00305
00306
00307 if (sx == 0.0 && sy == 0.0) {
00308 *dx = 0.0;
00309 *dy = 0.0;
00310 return;
00311 }
00312 if (sx >= -sy) {
00313 if (sx > sy) {
00314
00315 r = sx;
00316 if (sy > 0.0) theta = sy/r;
00317 else theta = 8.0f + sy/r;
00318 }
00319 else {
00320
00321 r = sy;
00322 theta = 2.0f - sx/r;
00323 }
00324 }
00325 else {
00326 if (sx <= sy) {
00327
00328 r = -sx;
00329 theta = 4.0f - sy/r;
00330 }
00331 else {
00332
00333 r = -sy;
00334 theta = 6.0f + sx/r;
00335 }
00336 }
00337 theta *= M_PI / 4.f;
00338 *dx = r * cosf(theta);
00339 *dy = r * sinf(theta);
00340 }
00341
00342
00343 void UniformSampleTriangle(float u1, float u2, float *u, float *v) {
00344 float su1 = sqrtf(u1);
00345 *u = 1.f - su1;
00346 *v = u2 * su1;
00347 }
00348
00349
00350 Distribution2D::Distribution2D(const float *func, int nu, int nv) {
00351 pConditionalV.reserve(nv);
00352 for (int v = 0; v < nv; ++v) {
00353
00354 pConditionalV.push_back(new Distribution1D(&func[v*nu], nu));
00355 }
00356
00357 vector<float> marginalFunc;
00358 marginalFunc.reserve(nv);
00359 for (int v = 0; v < nv; ++v)
00360 marginalFunc.push_back(pConditionalV[v]->funcInt);
00361 pMarginal = new Distribution1D(&marginalFunc[0], nv);
00362 }
00363
00364
00365 Distribution2D::~Distribution2D() {
00366 delete pMarginal;
00367 for (uint32_t i = 0; i < pConditionalV.size(); ++i)
00368 delete pConditionalV[i];
00369 }
00370
00371
00372 PermutedHalton::PermutedHalton(uint32_t d, RNG &rng) {
00373 dims = d;
00374
00375 b = new uint32_t[dims];
00376 uint32_t sumBases = 0;
00377 for (uint32_t i = 0; i < dims; ++i) {
00378 b[i] = primes[i];
00379 sumBases += b[i];
00380 }
00381
00382
00383 permute = new uint32_t[sumBases];
00384 uint32_t *p = permute;
00385 for (uint32_t i = 0; i < dims; ++i) {
00386 GeneratePermutation(p, b[i], rng);
00387 p += b[i];
00388 }
00389 }
00390
00391
00392 float UniformConePdf(float cosThetaMax) {
00393 return 1.f / (2.f * M_PI * (1.f - cosThetaMax));
00394 }
00395
00396
00397 Vector UniformSampleCone(float u1, float u2, float costhetamax) {
00398 float costheta = (1.f - u1) + u1 * costhetamax;
00399 float sintheta = sqrtf(1.f - costheta*costheta);
00400 float phi = u2 * 2.f * M_PI;
00401 return Vector(cosf(phi) * sintheta, sinf(phi) * sintheta, costheta);
00402 }
00403
00404
00405 Vector UniformSampleCone(float u1, float u2, float costhetamax,
00406 const Vector &x, const Vector &y, const Vector &z) {
00407 float costheta = Lerp(u1, costhetamax, 1.f);
00408 float sintheta = sqrtf(1.f - costheta*costheta);
00409 float phi = u2 * 2.f * M_PI;
00410 return cosf(phi) * sintheta * x + sinf(phi) * sintheta * y +
00411 costheta * z;
00412 }
00413
00414
00415 Vector SampleHG(const Vector &w, float g, float u1, float u2) {
00416 float costheta;
00417 if (fabsf(g) < 1e-3)
00418 costheta = 1.f - 2.f * u1;
00419 else {
00420 float sqrTerm = (1.f - g * g) /
00421 (1.f - g + 2.f * g * u1);
00422 costheta = (1.f + g * g - sqrTerm * sqrTerm) / (2.f * g);
00423 }
00424 float sintheta = sqrtf(max(0.f, 1.f-costheta*costheta));
00425 float phi = 2.f * M_PI * u2;
00426 Vector v1, v2;
00427 CoordinateSystem(w, &v1, &v2);
00428 return SphericalDirection(sintheta, costheta, phi, v1, v2, w);
00429 }
00430
00431
00432 float HGPdf(const Vector &w, const Vector &wp, float g) {
00433 return PhaseHG(w, wp, g);
00434 }
00435
00436