00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024 #include "lux.h"
00025 #include "randomgen.h"
00026 #include "shape.h"
00027 #include "mc.h"
00028 #include "volume.h"
00029
00030 #include <limits>
00031
00032 namespace lux
00033 {
00034
00035
00036 void ComputeStep1dCDF(float *f, int nSteps, float *c,
00037 float *cdf) {
00038
00039 int i;
00040 cdf[0] = 0.;
00041 for (i = 1; i < nSteps+1; ++i)
00042 cdf[i] = cdf[i-1] + f[i-1] / nSteps;
00043
00044 *c = cdf[nSteps];
00045 for (i = 1; i < nSteps+1; ++i)
00046 cdf[i] /= *c;
00047 }
00048 float SampleStep1d(float *f, float *cdf, float c,
00049 int nSteps, float u, float *pdf) {
00050
00051 float *ptr = std::lower_bound(cdf, cdf+nSteps+1, u);
00052 int offset = (int) (ptr-cdf-1);
00053
00054 u = (u - cdf[offset]) / (cdf[offset+1] - cdf[offset]);
00055 *pdf = f[offset] / c;
00056 return (offset + u) / nSteps;
00057 }
00058 void RejectionSampleDisk(const TsPack *tspack, float *x, float *y) {
00059 float sx, sy;
00060 do {
00061 sx = 1.f - 2.f * tspack->rng->floatValue();
00062 sy = 1.f - 2.f * tspack->rng->floatValue();
00063 } while (sx*sx + sy*sy > 1.f);
00064 *x = sx;
00065 *y = sy;
00066 }
00067 Vector UniformSampleHemisphere(float u1, float u2) {
00068 float z = u1;
00069 float r = sqrtf(max(0.f, 1.f - z*z));
00070 float phi = 2 * M_PI * u2;
00071 float x = r * cosf(phi);
00072 float y = r * sinf(phi);
00073 return Vector(x, y, z);
00074 }
00075 float UniformHemispherePdf(float theta, float phi) {
00076 return INV_TWOPI;
00077 }
00078 Vector UniformSampleSphere(float u1, float u2) {
00079 float z = 1.f - 2.f * u1;
00080 float r = sqrtf(max(0.f, 1.f - z*z));
00081 float phi = 2.f * M_PI * u2;
00082 float x = r * cosf(phi);
00083 float y = r * sinf(phi);
00084 return Vector(x, y, z);
00085 }
00086 float UniformSpherePdf() {
00087 return 1.f / (4.f * M_PI);
00088 }
00089 void UniformSampleDisk(float u1, float u2,
00090 float *x, float *y) {
00091 float r = sqrtf(u1);
00092 float theta = 2.0f * M_PI * u2;
00093 *x = r * cosf(theta);
00094 *y = r * sinf(theta);
00095 }
00096 void ConcentricSampleDisk(float u1, float u2,
00097 float *dx, float *dy) {
00098 float r, theta;
00099
00100 float sx = 2 * u1 - 1;
00101 float sy = 2 * u2 - 1;
00102
00103
00104 if (sx == 0.0 && sy == 0.0) {
00105 *dx = 0.0;
00106 *dy = 0.0;
00107 return;
00108 }
00109 if (sx >= -sy) {
00110 if (sx > sy) {
00111
00112 r = sx;
00113 if (sy > 0.0)
00114 theta = sy/r;
00115 else
00116 theta = 8.0f + sy/r;
00117 }
00118 else {
00119
00120 r = sy;
00121 theta = 2.0f - sx/r;
00122 }
00123 }
00124 else {
00125 if (sx <= sy) {
00126
00127 r = -sx;
00128 theta = 4.0f - sy/r;
00129 }
00130 else {
00131
00132 r = -sy;
00133 theta = 6.0f + sx/r;
00134 }
00135 }
00136 theta *= M_PI / 4.f;
00137 *dx = r*cosf(theta);
00138 *dy = r*sinf(theta);
00139 }
00140 void UniformSampleTriangle(float u1, float u2,
00141 float *u, float *v) {
00142 float su1 = sqrtf(u1);
00143 *u = 1.f - su1;
00144 *v = u2 * su1;
00145 }
00146 float UniformConePdf(float cosThetaMax) {
00147 return 1.f / (2.f * M_PI * (1.f - cosThetaMax));
00148 }
00149 Vector UniformSampleCone(float u1, float u2,
00150 float costhetamax) {
00151 float costheta = Lerp(u1, costhetamax, 1.f);
00152 float sintheta = sqrtf(1.f - costheta*costheta);
00153 float phi = u2 * 2.f * M_PI;
00154 return Vector(cosf(phi) * sintheta,
00155 sinf(phi) * sintheta,
00156 costheta);
00157 }
00158 Vector UniformSampleCone(float u1, float u2, float costhetamax,
00159 const Vector &x, const Vector &y, const Vector &z) {
00160 float costheta = Lerp(u1, costhetamax, 1.f);
00161 float sintheta = sqrtf(1.f - costheta*costheta);
00162 float phi = u2 * 2.f * M_PI;
00163 return cosf(phi) * sintheta * x + sinf(phi) * sintheta * y +
00164 costheta * z;
00165 }
00166 Vector SampleHG(const Vector &w, float g,
00167 float u1, float u2) {
00168 float costheta;
00169 if (fabsf(g) < 1e-3)
00170 costheta = 1.f - 2.f * u1;
00171 else {
00172
00173 float sqrTerm = (1.f - g * g) /
00174 (1.f - g + 2.f * g * u1);
00175 costheta = (1.f + g * g - sqrTerm * sqrTerm) / (2.f * g);
00176 }
00177 float sintheta = sqrtf(max(0.f, 1.f-costheta*costheta));
00178 float phi = 2.f * M_PI * u2;
00179 Vector v1, v2;
00180 CoordinateSystem(w, &v1, &v2);
00181 return SphericalDirection(sintheta, costheta,
00182 phi, v1, v2, w);
00183 }
00184 float HGPdf(const Vector &w, const Vector &wp,
00185 float g) {
00186 return PhaseHG(w, wp, g);
00187 }
00188
00189 #define A1 (-3.969683028665376e+01)
00190 #define A2 2.209460984245205e+02
00191 #define A3 (-2.759285104469687e+02)
00192 #define A4 1.383577518672690e+02
00193 #define A5 (-3.066479806614716e+01)
00194 #define A6 2.506628277459239e+00
00195
00196 #define B1 (-5.447609879822406e+01)
00197 #define B2 1.615858368580409e+02
00198 #define B3 (-1.556989798598866e+02)
00199 #define B4 6.680131188771972e+01
00200 #define B5 (-1.328068155288572e+01)
00201
00202 #define C1 (-7.784894002430293e-03)
00203 #define C2 (-3.223964580411365e-01)
00204 #define C3 (-2.400758277161838e+00)
00205 #define C4 (-2.549732539343734e+00)
00206 #define C5 4.374664141464968e+00
00207 #define C6 2.938163982698783e+00
00208
00209 #define D1 7.784695709041462e-03
00210 #define D2 3.224671290700398e-01
00211 #define D3 2.445134137142996e+00
00212 #define D4 3.754408661907416e+00
00213
00214 #define P_LOW 0.02425
00215
00216 #define P_HIGH 0.97575
00217
00218 double normsinv(double p)
00219 {
00220 if ((0 < p ) && (p < P_LOW)) {
00221 const double q = sqrt(-2 * log(p));
00222 return (((((C1 * q + C2) * q + C3) * q + C4) * q + C5) * q + C6) /
00223 ((((D1 * q + D2) * q + D3) * q + D4) * q + 1);
00224 } else if ((P_LOW <= p) && (p <= P_HIGH)) {
00225 const double q = p - 0.5;
00226 const double r = q * q;
00227 return (((((A1 * r + A2) * r + A3) * r + A4) * r + A5) * r + A6) * q /
00228 (((((B1 * r + B2) * r + B3) * r + B4) * r + B5) * r + 1);
00229 } else if ((P_HIGH < p) && (p < 1)) {
00230 const double q = sqrt(-2 * log(1 - p));
00231 return -(((((C1 * q + C2) * q + C3) * q + C4) * q + C5) * q + C6) /
00232 ((((D1 * q + D2) * q + D3) * q + D4) * q + 1);
00233 } else {
00234 return INFINITY;
00235 }
00236
00237
00238
00239
00240
00241
00242
00243
00244
00245 }
00246
00247
00248 float GaussianSampleDisk(float u1)
00249 {
00250 return Clamp<float>(normsinv(u1), 0.f, 1.f);
00251 }
00252 float InverseGaussianSampleDisk(float u1)
00253 {
00254 return Clamp<float>(1.f - normsinv(u1), 0.f, 1.f);
00255 }
00256 float ExponentialSampleDisk(float u1, int power)
00257 {
00258 return Clamp(-logf(u1) / power, 0.f, 1.f);
00259 }
00260 float InverseExponentialSampleDisk(float u1, int power)
00261 {
00262 return Clamp(1.f + logf(u1) / power, 0.f, 1.f);
00263 }
00264 float TriangularSampleDisk(float u1)
00265 {
00266 return Clamp(u1 <= .5f ? sqrtf(u1 * .5f) :
00267 1.f - sqrtf((1.f - u1) * .5f), 0.f, 1.f);
00268 }
00269
00270 }
00271