mc.cpp

Go to the documentation of this file.
00001 
00002 /*
00003  * pbrt source code Copyright(c) 1998-2007 Matt Pharr and Greg Humphreys
00004  *
00005  * All Rights Reserved.
00006  * For educational use only; commercial use expressly forbidden.
00007  * NO WARRANTY, express or implied, for this software.
00008  * (See file License.txt for complete license)
00009  */
00010 
00011 // mc.cpp*
00012 #include "pbrt.h"
00013 #include "geometry.h"
00014 #include "shape.h"
00015 #include "mc.h"
00016 #include "volume.h"
00017 // MC Function Definitions
00018 void ComputeStep1dCDF(float *f, int nSteps, float *c,
00019                 float *cdf) {
00020         // Compute integral of step function at $x_i$
00021         int i;
00022         cdf[0] = 0.;
00023         for (i = 1; i < nSteps+1; ++i)
00024                 cdf[i] = cdf[i-1] + f[i-1] / nSteps;
00025         // Transform step function integral into cdf
00026         *c = cdf[nSteps];
00027         for (i = 1; i < nSteps+1; ++i)
00028                 cdf[i] /= *c;
00029 }
00030 float SampleStep1d(float *f, float *cdf, float c,
00031                 int nSteps, float u, float *pdf) {
00032         // Find surrounding cdf segments
00033         float *ptr = std::lower_bound(cdf, cdf+nSteps+1, u);
00034         int offset = (int) (ptr-cdf-1);
00035         // Return offset along current cdf segment
00036         u = (u - cdf[offset]) / (cdf[offset+1] - cdf[offset]);
00037         *pdf = f[offset] / c;
00038         return (offset + u) / nSteps;
00039 }
00040 void RejectionSampleDisk(float *x, float *y) {
00041         float sx, sy;
00042         do {
00043                 sx = 1.f - 2.f * RandomFloat();
00044                 sy = 1.f - 2.f * RandomFloat();
00045         } while (sx*sx + sy*sy > 1.f);
00046         *x = sx;
00047         *y = sy;
00048 }
00049 COREDLL Vector UniformSampleHemisphere(float u1, float u2) {
00050         float z = u1;
00051         float r = sqrtf(max(0.f, 1.f - z*z));
00052         float phi = 2 * M_PI * u2;
00053         float x = r * cosf(phi);
00054         float y = r * sinf(phi);
00055         return Vector(x, y, z);
00056 }
00057 COREDLL float UniformHemispherePdf(float theta, float phi) {
00058         return INV_TWOPI;
00059 }
00060 COREDLL Vector UniformSampleSphere(float u1, float u2) {
00061         float z = 1.f - 2.f * u1;
00062         float r = sqrtf(max(0.f, 1.f - z*z));
00063         float phi = 2.f * M_PI * u2;
00064         float x = r * cosf(phi);
00065         float y = r * sinf(phi);
00066         return Vector(x, y, z);
00067 }
00068 COREDLL float UniformSpherePdf() {
00069         return 1.f / (4.f * M_PI);
00070 }
00071 COREDLL void UniformSampleDisk(float u1, float u2,
00072                 float *x, float *y) {
00073         float r = sqrtf(u1);
00074         float theta = 2.0f * M_PI * u2;
00075         *x = r * cosf(theta);
00076         *y = r * sinf(theta);
00077 }
00078 COREDLL void ConcentricSampleDisk(float u1, float u2,
00079                 float *dx, float *dy) {
00080         float r, theta;
00081         // Map uniform random numbers to $[-1,1]^2$
00082         float sx = 2 * u1 - 1;
00083         float sy = 2 * u2 - 1;
00084         // Map square to $(r,\theta)$
00085         // Handle degeneracy at the origin
00086         if (sx == 0.0 && sy == 0.0) {
00087                 *dx = 0.0;
00088                 *dy = 0.0;
00089                 return;
00090         }
00091         if (sx >= -sy) {
00092                 if (sx > sy) {
00093                         // Handle first region of disk
00094                         r = sx;
00095                         if (sy > 0.0)
00096                                 theta = sy/r;
00097                         else
00098                                 theta = 8.0f + sy/r;
00099                 }
00100                 else {
00101                         // Handle second region of disk
00102                         r = sy;
00103                         theta = 2.0f - sx/r;
00104                 }
00105         }
00106         else {
00107                 if (sx <= sy) {
00108                         // Handle third region of disk
00109                         r = -sx;
00110                         theta = 4.0f - sy/r;
00111                 }
00112                 else {
00113                         // Handle fourth region of disk
00114                         r = -sy;
00115                         theta = 6.0f + sx/r;
00116                 }
00117         }
00118         theta *= M_PI / 4.f;
00119         *dx = r*cosf(theta);
00120         *dy = r*sinf(theta);
00121 }
00122 COREDLL void UniformSampleTriangle(float u1, float u2,
00123                 float *u, float *v) {
00124         float su1 = sqrtf(u1);
00125         *u = 1.f - su1;
00126         *v = u2 * su1;
00127 }
00128 COREDLL float UniformConePdf(float cosThetaMax) {
00129         return 1.f / (2.f * M_PI * (1.f - cosThetaMax));
00130 }
00131 Vector UniformSampleCone(float u1, float u2,
00132                 float costhetamax) {
00133         float costheta = Lerp(u1, costhetamax, 1.f);
00134         float sintheta = sqrtf(1.f - costheta*costheta);
00135         float phi = u2 * 2.f * M_PI;
00136         return Vector(cosf(phi) * sintheta,
00137                       sinf(phi) * sintheta,
00138                           costheta);
00139 }
00140 COREDLL Vector UniformSampleCone(float u1, float u2, float costhetamax,
00141                 const Vector &x, const Vector &y, const Vector &z) {
00142         float costheta = Lerp(u1, costhetamax, 1.f);
00143         float sintheta = sqrtf(1.f - costheta*costheta);
00144         float phi = u2 * 2.f * M_PI;
00145         return cosf(phi) * sintheta * x + sinf(phi) * sintheta * y +
00146                 costheta * z;
00147 }
00148 COREDLL Vector SampleHG(const Vector &w, float g,
00149                 float u1, float u2) {
00150         float costheta;
00151         if (fabsf(g) < 1e-3)
00152                 costheta = 1.f - 2.f * u1;
00153         else
00154                 costheta = -1.f / (2.f * g) *
00155                         (1.f + g*g - ((1.f-g*g) * (1.f-g+2.f*g*u1)));
00156         float sintheta = sqrtf(max(0.f, 1.f-costheta*costheta));
00157         float phi = 2.f * M_PI * u2;
00158         Vector v1, v2;
00159         CoordinateSystem(w, &v1, &v2);
00160         return SphericalDirection(sintheta, costheta,
00161                 phi, v1, v2, w);
00162 }
00163 COREDLL float HGPdf(const Vector &w, const Vector &wp,
00164                 float g) {
00165         return PhaseHG(w, wp, g);
00166 }

Generated on Wed Sep 26 14:01:20 2007 for pbrt by  doxygen 1.5.1