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 "sh.h"
00028 #include "scene.h"
00029 #include "integrator.h"
00030 #include "intersection.h"
00031 #include "montecarlo.h"
00032 #include "imageio.h"
00033 #include <float.h>
00034
00035
00036 static void legendrep(float x, int lmax, float *out) {
00037 #define P(l,m) out[SHIndex(l,m)]
00038
00039 P(0,0) = 1.f;
00040 P(1,0) = x;
00041 for (int l = 2; l <= lmax; ++l)
00042 {
00043 P(l, 0) = ((2*l-1)*x*P(l-1,0) - (l-1)*P(l-2,0)) / l;
00044 Assert(!isnan(P(l, 0)));
00045 Assert(!isinf(P(l, 0)));
00046 }
00047
00048
00049 float neg = -1.f;
00050 float dfact = 1.f;
00051 float xroot = sqrtf(max(0.f, 1.f - x*x));
00052 float xpow = xroot;
00053 for (int l = 1; l <= lmax; ++l) {
00054 P(l, l) = neg * dfact * xpow;
00055 Assert(!isnan(P(l, l)));
00056 Assert(!isinf(P(l, l)));
00057 neg *= -1.f;
00058 dfact *= 2*l + 1;
00059 xpow *= xroot;
00060 }
00061
00062
00063 for (int l = 2; l <= lmax; ++l)
00064 {
00065 P(l, l-1) = x * (2*l-1) * P(l-1, l-1);
00066 Assert(!isnan(P(l, l-1)));
00067 Assert(!isinf(P(l, l-1)));
00068 }
00069
00070
00071 for (int l = 3; l <= lmax; ++l)
00072 for (int m = 1; m <= l-2; ++m)
00073 {
00074 P(l, m) = ((2 * (l-1) + 1) * x * P(l-1,m) -
00075 (l-1+m) * P(l-2,m)) / (l - m);
00076 Assert(!isnan(P(l, m)));
00077 Assert(!isinf(P(l, m)));
00078 }
00079 #if 0
00080
00081
00082 for (int l = 1; l <= lmax; ++l) {
00083 float fa = 1.f, fb = fact(2*l);
00084
00085 for (int m = -l; m < 0; ++m) {
00086 float neg = ((-m) & 0x1) ? -1.f : 1.f;
00087 P(l,m) = neg * fa/fb * P(l,-m);
00088 fb /= l-m;
00089 fa *= (l+m+1) > 1 ? (l+m+1) : 1.;
00090 }
00091 }
00092 #endif
00093 #undef P
00094 }
00095
00096
00097 static inline float fact(float v);
00098 static inline float divfact(int a, int b);
00099 static inline float K(int l, int m) {
00100 return sqrtf((2.f * l + 1.f) * INV_FOURPI * divfact(l, m));
00101 }
00102
00103
00104 static inline float divfact(int a, int b) {
00105 if (b == 0) return 1.f;
00106 float fa = a, fb = fabsf(b);
00107 float v = 1.f;
00108 for (float x = fa-fb+1.f; x <= fa+fb; x += 1.f)
00109 v *= x;
00110 return 1.f / v;
00111 }
00112
00113
00114
00115 static float dfact(float v) {
00116 if (v <= 1.f) return 1.f;
00117 return v * dfact(v - 2.f);
00118 }
00119
00120
00121 static inline float fact(float v) {
00122 if (v <= 1.f) return 1.f;
00123 return v * fact(v - 1.f);
00124 }
00125
00126
00127 static void sinCosIndexed(float s, float c, int n,
00128 float *sout, float *cout) {
00129 float si = 0, ci = 1;
00130 for (int i = 0; i < n; ++i) {
00131
00132 *sout++ = si;
00133 *cout++ = ci;
00134 float oldsi = si;
00135 si = si * c + ci * s;
00136 ci = ci * c - oldsi * s;
00137 }
00138 }
00139
00140
00141 static void toZYZ(const Matrix4x4 &m, float *alpha, float *beta, float *gamma) {
00142 #define M(a, b) (m.m[a][b])
00143
00144 float sy = sqrtf(M(2,1)*M(2,1) + M(2,0)*M(2,0));
00145 if (sy > 16*FLT_EPSILON) {
00146 *gamma = -atan2f(M(1,2), -M(0,2));
00147 *beta = -atan2f(sy, M(2,2));
00148 *alpha = -atan2f(M(2,1), M(2,0));
00149 } else {
00150 *gamma = 0;
00151 *beta = -atan2f(sy, M(2,2));
00152 *alpha = -atan2f(-M(1,0), M(1,1));
00153 }
00154 #undef M
00155 }
00156
00157
00158 static inline float lambda(float l) {
00159 return sqrtf((4.f * M_PI) / (2.f * l + 1.));
00160 }
00161
00162
00163
00164
00165 void SHEvaluate(const Vector &w, int lmax, float *out) {
00166 if (lmax > 28)
00167 Severe("SHEvaluate() runs out of numerical precision for lmax > 28. "
00168 "If you need more bands, try recompiling using doubles.");
00169
00170 Assert(w.Length() > .995f && w.Length() < 1.005f);
00171 legendrep(w.z, lmax, out);
00172
00173
00174 float *Klm = ALLOCA(float, SHTerms(lmax));
00175 for (int l = 0; l <= lmax; ++l)
00176 for (int m = -l; m <= l; ++m)
00177 Klm[SHIndex(l, m)] = K(l, m);
00178
00179
00180 float *sins = ALLOCA(float, lmax+1), *coss = ALLOCA(float, lmax+1);
00181 float xyLen = sqrtf(max(0.f, 1.f - w.z*w.z));
00182 if (xyLen == 0.f) {
00183 for (int i = 0; i <= lmax; ++i) sins[i] = 0.f;
00184 for (int i = 0; i <= lmax; ++i) coss[i] = 1.f;
00185 }
00186 else
00187 sinCosIndexed(w.y / xyLen, w.x / xyLen, lmax+1, sins, coss);
00188
00189
00190 static const float sqrt2 = sqrtf(2.f);
00191 for (int l = 0; l <= lmax; ++l) {
00192 for (int m = -l; m < 0; ++m)
00193 {
00194 out[SHIndex(l, m)] = sqrt2 * Klm[SHIndex(l, m)] *
00195 out[SHIndex(l, -m)] * sins[-m];
00196 Assert(!isnan(out[SHIndex(l,m)]));
00197 Assert(!isinf(out[SHIndex(l,m)]));
00198 }
00199 out[SHIndex(l, 0)] *= Klm[SHIndex(l, 0)];
00200 for (int m = 1; m <= l; ++m)
00201 {
00202 out[SHIndex(l, m)] *= sqrt2 * Klm[SHIndex(l, m)] * coss[m];
00203 Assert(!isnan(out[SHIndex(l,m)]));
00204 Assert(!isinf(out[SHIndex(l,m)]));
00205 }
00206 }
00207 }
00208
00209
00210 #if 0
00211
00212 void SHEvaluate(float costheta, float cosphi, float sinphi, int lmax, float *out) {
00213 legendrep(costheta, lmax, out);
00214
00215 float *Klm = ALLOCA(float, SHTerms(lmax));
00216 klm(lmax, Klm);
00217
00218 float sqrt2 = sqrtf(2.f);
00219 float sins[(lmax+1)], coss[(lmax+1)];
00220 sinCosIndexed(sinphi, cosphi, lmax+1, sins, coss);
00221
00222 for (int l = 0; l <= lmax; ++l) {
00223 for (int m = -l; m < 0; ++m)
00224
00225 out[SHIndex(l, m)] = sqrt2 * Klm[SHIndex(l, m)] *
00226 out[SHIndex(l, -m)] * sins[-m];
00227
00228 out[SHIndex(l, 0)] *= Klm[SHIndex(l, 0)];
00229
00230 for (int m = 1; m <= l; ++m)
00231 out[SHIndex(l, m)] *= sqrt2 * Klm[SHIndex(l, m)] * coss[m];
00232 }
00233 }
00234
00235
00236 #endif
00237 void SHWriteImage(const char *filename, const Spectrum *c, int lmax, int yres) {
00238 int xres = 2 * yres;
00239 float *rgb = new float[xres * yres * 3];
00240 float *Ylm = ALLOCA(float, SHTerms(lmax));
00241 for (int y = 0; y < yres; ++y) {
00242 float theta = (float(y) + 0.5f) / float(yres) * M_PI;
00243 for (int x = 0; x < xres; ++x) {
00244 float phi = (float(x) + 0.5f) / float(xres) * 2.f * M_PI;
00245
00246 Vector w = SphericalDirection(sinf(theta), cosf(theta), phi);
00247 SHEvaluate(w, lmax, Ylm);
00248 Spectrum val = 0.f;
00249 for (int i = 0; i < SHTerms(lmax); ++i)
00250 val += Ylm[i] * c[i];
00251 int offset = xres * y + x;
00252 val.ToRGB(&rgb[3*offset]);
00253 }
00254 }
00255
00256 WriteImage(filename, rgb, NULL, xres, yres, xres, yres, 0, 0);
00257 delete[] rgb;
00258 }
00259
00260
00261 void SHProjectIncidentDirectRadiance(const Point &p, float pEpsilon,
00262 float time, MemoryArena &arena, const Scene *scene,
00263 bool computeLightVis, int lmax, RNG &rng, Spectrum *c_d) {
00264
00265 Spectrum *c = arena.Alloc<Spectrum>(SHTerms(lmax));
00266 for (uint32_t i = 0; i < scene->lights.size(); ++i) {
00267 Light *light = scene->lights[i];
00268 light->SHProject(p, pEpsilon, lmax, scene, computeLightVis, time,
00269 rng, c);
00270 for (int j = 0; j < SHTerms(lmax); ++j)
00271 c_d[j] += c[j];
00272 }
00273 SHReduceRinging(c_d, lmax);
00274 }
00275
00276
00277 void SHProjectIncidentIndirectRadiance(const Point &p, float pEpsilon,
00278 float time, const Renderer *renderer, Sample *origSample,
00279 const Scene *scene, int lmax, RNG &rng, int ns, Spectrum *c_i) {
00280 Sample *sample = origSample->Duplicate(1);
00281 MemoryArena arena;
00282 uint32_t scramble[2] = { rng.RandomUInt(), rng.RandomUInt() };
00283 int nSamples = RoundUpPow2(ns);
00284 float *Ylm = ALLOCA(float, SHTerms(lmax));
00285 for (int i = 0; i < nSamples; ++i) {
00286
00287 float u[2];
00288 Sample02(i, scramble, u);
00289 Vector wi = UniformSampleSphere(u[0], u[1]);
00290 float pdf = UniformSpherePdf();
00291
00292
00293 Spectrum Li = 0.f;
00294 RayDifferential ray(p, wi, pEpsilon, INFINITY, time);
00295
00296
00297 sample->time = time;
00298 for (uint32_t j = 0; j < sample->n1D.size(); ++j)
00299 for (uint32_t k = 0; k < sample->n1D[j]; ++k)
00300 sample->oneD[j][k] = rng.RandomFloat();
00301 for (uint32_t j = 0; j < sample->n2D.size(); ++j)
00302 for (uint32_t k = 0; k < 2 * sample->n2D[j]; ++k)
00303 sample->twoD[j][k] = rng.RandomFloat();
00304 Li = renderer->Li(scene, ray, sample, rng, arena);
00305
00306
00307 SHEvaluate(wi, lmax, Ylm);
00308 for (int j = 0; j < SHTerms(lmax); ++j)
00309 c_i[j] += Ylm[j] * Li / (pdf * nSamples);
00310 arena.FreeAll();
00311 }
00312 delete[] sample;
00313 }
00314
00315
00316 void SHReduceRinging(Spectrum *c, int lmax, float lambda) {
00317 for (int l = 0; l <= lmax; ++l) {
00318 float scale = 1.f / (1.f + lambda * l * l * (l + 1) * (l + 1));
00319 for (int m = -l; m <= l; ++m)
00320 c[SHIndex(l, m)] *= scale;
00321 }
00322 }
00323
00324
00325 void SHRotate(const Spectrum *c_in, Spectrum *c_out, const Matrix4x4 &m,
00326 int lmax, MemoryArena &arena) {
00327 float alpha, beta, gamma;
00328 toZYZ(m, &alpha, &beta, &gamma);
00329 Spectrum *work = arena.Alloc<Spectrum>(SHTerms(lmax));
00330 SHRotateZ(c_in, c_out, gamma, lmax);
00331 SHRotateXPlus(c_out, work, lmax);
00332 SHRotateZ(work, c_out, beta, lmax);
00333 SHRotateXMinus(c_out, work, lmax);
00334 SHRotateZ(work, c_out, alpha, lmax);
00335 }
00336
00337
00338 void SHRotateZ(const Spectrum *c_in, Spectrum *c_out, float alpha,
00339 int lmax) {
00340 Assert(c_in != c_out);
00341 c_out[0] = c_in[0];
00342 if (lmax == 0) return;
00343
00344 float *ct = ALLOCA(float, lmax+1);
00345 float *st = ALLOCA(float, lmax+1);
00346 sinCosIndexed(sinf(alpha), cosf(alpha), lmax+1, st, ct);
00347 for (int l = 1; l <= lmax; ++l) {
00348
00349 for (int m = -l; m < 0; ++m)
00350 c_out[SHIndex(l, m)] =
00351 ( ct[-m] * c_in[SHIndex(l, m)] +
00352 -st[-m] * c_in[SHIndex(l, -m)]);
00353 c_out[SHIndex(l, 0)] = c_in[SHIndex(l, 0)];
00354 for (int m = 1; m <= l; ++m)
00355 c_out[SHIndex(l, m)] =
00356 (ct[m] * c_in[SHIndex(l, m)] +
00357 st[m] * c_in[SHIndex(l, -m)]);
00358 }
00359 }
00360
00361
00362 void SHConvolveCosTheta(int lmax, const Spectrum *c_in,
00363 Spectrum *c_out) {
00364 static const float c_costheta[18] = { 0.8862268925, 1.0233267546,
00365 0.4954159260, 0.0000000000, -0.1107783690, 0.0000000000,
00366 0.0499271341, 0.0000000000, -0.0285469331, 0.0000000000,
00367 0.0185080823, 0.0000000000, -0.0129818395, 0.0000000000,
00368 0.0096125342, 0.0000000000, -0.0074057109, 0.0000000000 };
00369 for (int l = 0; l <= lmax; ++l)
00370 for (int m = -l; m <= l; ++m) {
00371 int o = SHIndex(l, m);
00372 if (l < 18) c_out[o] = lambda(l) * c_in[o] * c_costheta[l];
00373 else c_out[o] = 0.f;
00374 }
00375 }
00376
00377
00378 void SHConvolvePhong(int lmax, float n, const Spectrum *c_in,
00379 Spectrum *c_out) {
00380 for (int l = 0; l <= lmax; ++l) {
00381 float c_phong = expf(-(l*l) / (2.f * n));
00382 for (int m = -l; m <= l; ++m) {
00383 int o = SHIndex(l, m);
00384 c_out[o] = lambda(l) * c_in[o] * c_phong;
00385 }
00386 }
00387 }
00388
00389
00390 void SHComputeDiffuseTransfer(const Point &p, const Normal &n,
00391 float rayEpsilon, const Scene *scene, RNG &rng, int nSamples,
00392 int lmax, Spectrum *c_transfer) {
00393 for (int i = 0; i < SHTerms(lmax); ++i)
00394 c_transfer[i] = 0.f;
00395 uint32_t scramble[2] = { rng.RandomUInt(), rng.RandomUInt() };
00396 float *Ylm = ALLOCA(float, SHTerms(lmax));
00397 for (int i = 0; i < nSamples; ++i) {
00398
00399 float u[2];
00400 Sample02(i, scramble, u);
00401 Vector w = UniformSampleSphere(u[0], u[1]);
00402 float pdf = UniformSpherePdf();
00403 if (Dot(w, n) > 0.f && !scene->IntersectP(Ray(p, w, rayEpsilon))) {
00404
00405 SHEvaluate(w, lmax, Ylm);
00406 for (int j = 0; j < SHTerms(lmax); ++j)
00407 c_transfer[j] += (Ylm[j] * AbsDot(w, n)) / (pdf * nSamples);
00408 }
00409 }
00410 }
00411
00412
00413 void SHComputeTransferMatrix(const Point &p, float rayEpsilon,
00414 const Scene *scene, RNG &rng, int nSamples, int lmax,
00415 Spectrum *T) {
00416 for (int i = 0; i < SHTerms(lmax)*SHTerms(lmax); ++i)
00417 T[i] = 0.f;
00418 uint32_t scramble[2] = { rng.RandomUInt(), rng.RandomUInt() };
00419 float *Ylm = ALLOCA(float, SHTerms(lmax));
00420 for (int i = 0; i < nSamples; ++i) {
00421
00422 float u[2];
00423 Sample02(i, scramble, u);
00424 Vector w = UniformSampleSphere(u[0], u[1]);
00425 float pdf = UniformSpherePdf();
00426 if (!scene->IntersectP(Ray(p, w, rayEpsilon))) {
00427
00428 SHEvaluate(w, lmax, Ylm);
00429 for (int j = 0; j < SHTerms(lmax); ++j)
00430 for (int k = 0; k < SHTerms(lmax); ++k)
00431 T[j*SHTerms(lmax)+k] += (Ylm[j] * Ylm[k]) / (pdf * nSamples);
00432 }
00433 }
00434 }
00435
00436
00437 void SHComputeBSDFMatrix(const Spectrum &Kd, const Spectrum &Ks,
00438 float roughness, RNG &rng, int nSamples, int lmax, Spectrum *B) {
00439 for (int i = 0; i < SHTerms(lmax)*SHTerms(lmax); ++i)
00440 B[i] = 0.f;
00441
00442 MemoryArena arena;
00443 DifferentialGeometry dg(Point(0,0,0), Vector(1,0,0), Vector(0,1,0),
00444 Normal(0,0,0), Normal(0,0,0), 0, 0, NULL);
00445 BSDF *bsdf = BSDF_ALLOC(arena, BSDF)(dg, Normal(0,0,1));
00446 bsdf->Add(BSDF_ALLOC(arena, Lambertian)(Spectrum(Kd)));
00447 Fresnel *fresnel = BSDF_ALLOC(arena, FresnelDielectric)(1.5f, 1.f);
00448 bsdf->Add(BSDF_ALLOC(arena, Microfacet)(Ks, fresnel,
00449 BSDF_ALLOC(arena, Blinn)(1.f / roughness)));
00450
00451
00452 float *Ylm = new float[SHTerms(lmax) * nSamples];
00453 Vector *w = new Vector[nSamples];
00454 uint32_t scramble[2] = { rng.RandomUInt(), rng.RandomUInt() };
00455 for (int i = 0; i < nSamples; ++i) {
00456 float u[2];
00457 Sample02(i, scramble, u);
00458 w[i] = UniformSampleSphere(u[0], u[1]);
00459 SHEvaluate(w[i], lmax, &Ylm[SHTerms(lmax)*i]);
00460 }
00461
00462
00463 for (int osamp = 0; osamp < nSamples; ++osamp) {
00464 const Vector &wo = w[osamp];
00465 for (int isamp = 0; isamp < nSamples; ++isamp) {
00466 const Vector &wi = w[isamp];
00467
00468 Spectrum f = bsdf->f(wo, wi);
00469 if (!f.IsBlack()) {
00470 float pdf = UniformSpherePdf() * UniformSpherePdf();
00471 f *= fabsf(CosTheta(wi)) / (pdf * nSamples * nSamples);
00472 for (int i = 0; i < SHTerms(lmax); ++i)
00473 for (int j = 0; j < SHTerms(lmax); ++j)
00474 B[i*SHTerms(lmax)+j] += f * Ylm[isamp*SHTerms(lmax)+j] *
00475 Ylm[osamp*SHTerms(lmax)+i];
00476 }
00477 }
00478 }
00479
00480
00481 delete[] w;
00482 delete[] Ylm;
00483 }
00484
00485
00486 void SHMatrixVectorMultiply(const Spectrum *M, const Spectrum *v,
00487 Spectrum *vout, int lmax) {
00488 for (int i = 0; i < SHTerms(lmax); ++i) {
00489 vout[i] = 0.f;
00490 for (int j = 0; j < SHTerms(lmax); ++j)
00491 vout[i] += M[SHTerms(lmax) * i + j] * v[j];
00492 }
00493 }
00494
00495