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 "geometry.h"
00026 #include "shape.h"
00027 #include "mc.h"
00028 #include "volume.h"
00029
00030 namespace lux
00031 {
00032
00033
00034 void ComputeStep1dCDF(float *f, int nSteps, float *c,
00035 float *cdf) {
00036
00037 int i;
00038 cdf[0] = 0.;
00039 for (i = 1; i < nSteps+1; ++i)
00040 cdf[i] = cdf[i-1] + f[i-1] / nSteps;
00041
00042 *c = cdf[nSteps];
00043 for (i = 1; i < nSteps+1; ++i)
00044 cdf[i] /= *c;
00045 }
00046 float SampleStep1d(float *f, float *cdf, float c,
00047 int nSteps, float u, float *pdf) {
00048
00049 float *ptr = std::lower_bound(cdf, cdf+nSteps+1, u);
00050 int offset = (int) (ptr-cdf-1);
00051
00052 u = (u - cdf[offset]) / (cdf[offset+1] - cdf[offset]);
00053 *pdf = f[offset] / c;
00054 return (offset + u) / nSteps;
00055 }
00056 void RejectionSampleDisk(float *x, float *y) {
00057 float sx, sy;
00058 do {
00059 sx = 1.f - 2.f * lux::random::floatValue();
00060 sy = 1.f - 2.f * lux::random::floatValue();
00061 } while (sx*sx + sy*sy > 1.f);
00062 *x = sx;
00063 *y = sy;
00064 }
00065 Vector UniformSampleHemisphere(float u1, float u2) {
00066 float z = u1;
00067 float r = sqrtf(max(0.f, 1.f - z*z));
00068 float phi = 2 * M_PI * u2;
00069 float x = r * cosf(phi);
00070 float y = r * sinf(phi);
00071 return Vector(x, y, z);
00072 }
00073 float UniformHemispherePdf(float theta, float phi) {
00074 return INV_TWOPI;
00075 }
00076 Vector UniformSampleSphere(float u1, float u2) {
00077 float z = 1.f - 2.f * u1;
00078 float r = sqrtf(max(0.f, 1.f - z*z));
00079 float phi = 2.f * M_PI * u2;
00080 float x = r * cosf(phi);
00081 float y = r * sinf(phi);
00082 return Vector(x, y, z);
00083 }
00084 float UniformSpherePdf() {
00085 return 1.f / (4.f * M_PI);
00086 }
00087 void UniformSampleDisk(float u1, float u2,
00088 float *x, float *y) {
00089 float r = sqrtf(u1);
00090 float theta = 2.0f * M_PI * u2;
00091 *x = r * cosf(theta);
00092 *y = r * sinf(theta);
00093 }
00094 void ConcentricSampleDisk(float u1, float u2,
00095 float *dx, float *dy) {
00096 float r, theta;
00097
00098 float sx = 2 * u1 - 1;
00099 float sy = 2 * u2 - 1;
00100
00101
00102 if (sx == 0.0 && sy == 0.0) {
00103 *dx = 0.0;
00104 *dy = 0.0;
00105 return;
00106 }
00107 if (sx >= -sy) {
00108 if (sx > sy) {
00109
00110 r = sx;
00111 if (sy > 0.0)
00112 theta = sy/r;
00113 else
00114 theta = 8.0f + sy/r;
00115 }
00116 else {
00117
00118 r = sy;
00119 theta = 2.0f - sx/r;
00120 }
00121 }
00122 else {
00123 if (sx <= sy) {
00124
00125 r = -sx;
00126 theta = 4.0f - sy/r;
00127 }
00128 else {
00129
00130 r = -sy;
00131 theta = 6.0f + sx/r;
00132 }
00133 }
00134 theta *= M_PI / 4.f;
00135 *dx = r*cosf(theta);
00136 *dy = r*sinf(theta);
00137 }
00138 void UniformSampleTriangle(float u1, float u2,
00139 float *u, float *v) {
00140 float su1 = sqrtf(u1);
00141 *u = 1.f - su1;
00142 *v = u2 * su1;
00143 }
00144 float UniformConePdf(float cosThetaMax) {
00145 return 1.f / (2.f * M_PI * (1.f - cosThetaMax));
00146 }
00147 Vector UniformSampleCone(float u1, float u2,
00148 float costhetamax) {
00149 float costheta = Lerp(u1, costhetamax, 1.f);
00150 float sintheta = sqrtf(1.f - costheta*costheta);
00151 float phi = u2 * 2.f * M_PI;
00152 return Vector(cosf(phi) * sintheta,
00153 sinf(phi) * sintheta,
00154 costheta);
00155 }
00156 Vector UniformSampleCone(float u1, float u2, float costhetamax,
00157 const Vector &x, const Vector &y, const Vector &z) {
00158 float costheta = Lerp(u1, costhetamax, 1.f);
00159 float sintheta = sqrtf(1.f - costheta*costheta);
00160 float phi = u2 * 2.f * M_PI;
00161 return cosf(phi) * sintheta * x + sinf(phi) * sintheta * y +
00162 costheta * z;
00163 }
00164 Vector SampleHG(const Vector &w, float g,
00165 float u1, float u2) {
00166 float costheta;
00167 if (fabsf(g) < 1e-3)
00168 costheta = 1.f - 2.f * u1;
00169 else {
00170
00171 float sqrTerm = (1.f - g * g) /
00172 (1.f - g + 2.f * g * u1);
00173 costheta = (1.f + g * g - sqrTerm * sqrTerm) / (2.f * g);
00174 }
00175 float sintheta = sqrtf(max(0.f, 1.f-costheta*costheta));
00176 float phi = 2.f * M_PI * u2;
00177 Vector v1, v2;
00178 CoordinateSystem(w, &v1, &v2);
00179 return SphericalDirection(sintheta, costheta,
00180 phi, v1, v2, w);
00181 }
00182 float HGPdf(const Vector &w, const Vector &wp,
00183 float g) {
00184 return PhaseHG(w, wp, g);
00185 }
00186
00187 }
00188