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 #include "pbrt.h"
00026 #include "transport.h"
00027 #include "scene.h"
00028 #include "mc.h"
00029 #include "kdtree.h"
00030 #include "sampling.h"
00031
00032 struct ClosePhoton;
00033
00034
00035 struct Photon {
00036 Photon(const Point &pp, const Spectrum &wt, const Vector &w)
00037 : p(pp), alpha(wt), wi(w) {
00038 }
00039 Photon() { }
00040 Point p;
00041 Spectrum alpha;
00042 Vector wi;
00043 };
00044 struct RadiancePhoton {
00045 RadiancePhoton(const Point &pp, const Normal &nn)
00046 : p(pp), n(nn), Lo(0.f) {
00047 }
00048 RadiancePhoton() { }
00049 Point p;
00050 Normal n;
00051 Spectrum Lo;
00052 };
00053 struct RadiancePhotonProcess {
00054
00055 RadiancePhotonProcess(const Point &pp, const Normal &nn)
00056 : p(pp), n(nn) {
00057 photon = NULL;
00058 }
00059 void operator()(const RadiancePhoton &rp,
00060 float distSquared, float &maxDistSquared) const {
00061 if (Dot(rp.n, n) > 0) {
00062 photon = &rp;
00063 maxDistSquared = distSquared;
00064 }
00065 }
00066 const Point &p;
00067 const Normal &n;
00068 mutable const RadiancePhoton *photon;
00069 };
00070 inline float kernel(const Photon *photon, const Point &p,
00071 float md2) {
00072
00073 float s = (1.f - DistanceSquared(photon->p, p) / md2);
00074 return 3.f / (md2 * M_PI) * s * s;
00075 }
00076
00077 struct PhotonProcess {
00078
00079 PhotonProcess(u_int mp, const Point &p);
00080 void operator()(const Photon &photon, float dist2, float &maxDistSquared) const;
00081 const Point &p;
00082 ClosePhoton *photons;
00083 u_int nLookup;
00084 mutable u_int foundPhotons;
00085 };
00086 struct ClosePhoton {
00087 ClosePhoton(const Photon *p = NULL,
00088 float md2 = INFINITY) {
00089 photon = p;
00090 distanceSquared = md2;
00091 }
00092 bool operator<(const ClosePhoton &p2) const {
00093 return distanceSquared == p2.distanceSquared ? (photon < p2.photon) :
00094 distanceSquared < p2.distanceSquared;
00095 }
00096 const Photon *photon;
00097 float distanceSquared;
00098 };
00099
00100 PhotonProcess::PhotonProcess(u_int mp, const Point &P)
00101 : p(P) {
00102 photons = 0;
00103 nLookup = mp;
00104 foundPhotons = 0;
00105 }
00106
00107 class ExPhotonIntegrator : public SurfaceIntegrator {
00108 public:
00109
00110 ExPhotonIntegrator(int ncaus, int nindir, int nLookup, int mdepth,
00111 float maxdist, bool finalGather, int gatherSamples,
00112 float rrt, float ga);
00113 ~ExPhotonIntegrator();
00114 Spectrum Li(const Scene *scene, const RayDifferential &ray,
00115 const Sample *sample, float *alpha) const;
00116 void RequestSamples(Sample *sample, const Scene *scene);
00117 void Preprocess(const Scene *);
00118 private:
00119 static inline bool unsuccessful(int needed, int found, int shot) {
00120 return (found < needed &&
00121 (found == 0 || found < shot / 1024));
00122 }
00123 static Spectrum LPhoton(KdTree<Photon, PhotonProcess> *map,
00124 int nPaths, int nLookup, BSDF *bsdf, const Intersection &isect,
00125 const Vector &w, float maxDistSquared);
00126
00127 Spectrum estimateE(KdTree<Photon, PhotonProcess> *map, int count,
00128 const Point &p, const Normal &n) const;
00129
00130
00131 int gatherSampleOffset[2], gatherComponentOffset[2];
00132 u_int nCausticPhotons, nIndirectPhotons;
00133 u_int nLookup;
00134 mutable int specularDepth;
00135 int maxSpecularDepth;
00136 float maxDistSquared, rrTreshold;
00137 bool finalGather;
00138 float cosGatherAngle;
00139 int gatherSamples;
00140
00141 int *lightSampleOffset, lightNumOffset;
00142 int *bsdfSampleOffset, *bsdfComponentOffset;
00143 int nCausticPaths, nIndirectPaths;
00144 mutable KdTree<Photon, PhotonProcess> *causticMap;
00145 mutable KdTree<Photon, PhotonProcess> *indirectMap;
00146 mutable KdTree<RadiancePhoton, RadiancePhotonProcess> *radianceMap;
00147 };
00148
00149
00150 Spectrum ExPhotonIntegrator::estimateE(
00151 KdTree<Photon, PhotonProcess> *map, int count,
00152 const Point &p, const Normal &n) const {
00153 if (!map) return 0.f;
00154
00155 PhotonProcess proc(nLookup, p);
00156 proc.photons = (ClosePhoton *)alloca(nLookup *
00157 sizeof(ClosePhoton));
00158 float md2 = maxDistSquared;
00159 map->Lookup(p, proc, md2);
00160
00161 ClosePhoton *photons = proc.photons;
00162 Spectrum E(0.);
00163 for (u_int i = 0; i < proc.foundPhotons; ++i)
00164 if (Dot(n, photons[i].photon->wi) > 0.)
00165 E += photons[i].photon->alpha;
00166 return E / (float(count) * md2 * M_PI);
00167 }
00168 void PhotonProcess::operator()(const Photon &photon,
00169 float distSquared, float &maxDistSquared) const {
00170
00171 static StatsPercentage discarded("Photon Map", "Discarded photons");
00172 discarded.Add(0, 1);
00173 if (foundPhotons < nLookup) {
00174
00175 photons[foundPhotons++] = ClosePhoton(&photon, distSquared);
00176 if (foundPhotons == nLookup) {
00177 std::make_heap(&photons[0], &photons[nLookup]);
00178 maxDistSquared = photons[0].distanceSquared;
00179 }
00180 }
00181 else {
00182
00183 discarded.Add(1, 0);
00184 std::pop_heap(&photons[0], &photons[nLookup]);
00185 photons[nLookup-1] = ClosePhoton(&photon, distSquared);
00186 std::push_heap(&photons[0], &photons[nLookup]);
00187 maxDistSquared = photons[0].distanceSquared;
00188 }
00189 }
00190 Spectrum ExPhotonIntegrator::LPhoton(
00191 KdTree<Photon, PhotonProcess> *map,
00192 int nPaths, int nLookup, BSDF *bsdf,
00193 const Intersection &isect, const Vector &wo,
00194 float maxDistSquared) {
00195 Spectrum L(0.);
00196 if (!map) return L;
00197 BxDFType nonSpecular = BxDFType(BSDF_REFLECTION |
00198 BSDF_TRANSMISSION | BSDF_DIFFUSE | BSDF_GLOSSY);
00199 if (bsdf->NumComponents(nonSpecular) == 0)
00200 return L;
00201 static StatsCounter lookups("Photon Map", "Total lookups");
00202
00203 PhotonProcess proc(nLookup, isect.dg.p);
00204 proc.photons =
00205 (ClosePhoton *)alloca(nLookup * sizeof(ClosePhoton));
00206
00207 ++lookups;
00208 map->Lookup(isect.dg.p, proc, maxDistSquared);
00209
00210 static StatsRatio foundRate("Photon Map", "Photons found per lookup");
00211 foundRate.Add(proc.foundPhotons, 1);
00212
00213 ClosePhoton *photons = proc.photons;
00214 int nFound = proc.foundPhotons;
00215 Normal Nf = Dot(wo, bsdf->dgShading.nn) < 0 ? -bsdf->dgShading.nn :
00216 bsdf->dgShading.nn;
00217
00218 if (bsdf->NumComponents(BxDFType(BSDF_REFLECTION |
00219 BSDF_TRANSMISSION | BSDF_GLOSSY)) > 0) {
00220
00221 for (int i = 0; i < nFound; ++i) {
00222 const Photon *p = photons[i].photon;
00223 BxDFType flag = Dot(Nf, p->wi) > 0.f ?
00224 BSDF_ALL_REFLECTION : BSDF_ALL_TRANSMISSION;
00225 float k = kernel(p, isect.dg.p, maxDistSquared);
00226 L += (k / nPaths) * bsdf->f(wo, p->wi, flag) * p->alpha;
00227 }
00228 }
00229 else {
00230
00231 Spectrum Lr(0.), Lt(0.);
00232 for (int i = 0; i < nFound; ++i) {
00233 if (Dot(Nf, photons[i].photon->wi) > 0.f) {
00234 float k = kernel(photons[i].photon, isect.dg.p,
00235 maxDistSquared);
00236 Lr += (k / nPaths) * photons[i].photon->alpha;
00237 }
00238 else {
00239 float k = kernel(photons[i].photon, isect.dg.p,
00240 maxDistSquared);
00241 Lt += (k / nPaths) * photons[i].photon->alpha;
00242 }
00243 }
00244 L += Lr * bsdf->rho(wo, BSDF_ALL_REFLECTION) * INV_PI +
00245 Lt * bsdf->rho(wo, BSDF_ALL_TRANSMISSION) * INV_PI;
00246 }
00247 return L;
00248 }
00249
00250 ExPhotonIntegrator::ExPhotonIntegrator(int ncaus, int nind,
00251 int nl, int mdepth, float mdist, bool fg,
00252 int gs, float rrt, float ga) {
00253 nCausticPhotons = ncaus;
00254 nIndirectPhotons = nind;
00255 nLookup = nl;
00256 maxDistSquared = mdist * mdist;
00257 maxSpecularDepth = mdepth;
00258 causticMap = indirectMap = NULL;
00259 radianceMap = NULL;
00260 specularDepth = 0;
00261 finalGather = fg;
00262 gatherSamples = gs;
00263 rrTreshold = rrt;
00264 cosGatherAngle = cos(Radians(ga));
00265 }
00266 ExPhotonIntegrator::~ExPhotonIntegrator() {
00267 delete causticMap;
00268 delete indirectMap;
00269 delete radianceMap;
00270 }
00271 void ExPhotonIntegrator::RequestSamples(Sample *sample,
00272 const Scene *scene) {
00273
00274 u_int nLights = scene->lights.size();
00275 lightSampleOffset = new int[nLights];
00276 bsdfSampleOffset = new int[nLights];
00277 bsdfComponentOffset = new int[nLights];
00278 for (u_int i = 0; i < nLights; ++i) {
00279 const Light *light = scene->lights[i];
00280 int lightSamples =
00281 scene->sampler->RoundSize(light->nSamples);
00282 lightSampleOffset[i] = sample->Add2D(lightSamples);
00283 bsdfSampleOffset[i] = sample->Add2D(lightSamples);
00284 bsdfComponentOffset[i] = sample->Add1D(lightSamples);
00285 }
00286 lightNumOffset = -1;
00287
00288 if (finalGather) {
00289 gatherSamples = scene->sampler->RoundSize(max(1,
00290 gatherSamples/2));
00291 gatherSampleOffset[0] = sample->Add2D(gatherSamples);
00292 gatherSampleOffset[1] = sample->Add2D(gatherSamples);
00293 gatherComponentOffset[0] = sample->Add1D(gatherSamples);
00294 gatherComponentOffset[1] = sample->Add1D(gatherSamples);
00295 }
00296 }
00297
00298 void ExPhotonIntegrator::Preprocess(const Scene *scene) {
00299 if (scene->lights.size() == 0) return;
00300 ProgressReporter progress(nCausticPhotons+
00301 nIndirectPhotons, "Shooting photons");
00302 vector<Photon> causticPhotons;
00303 vector<Photon> indirectPhotons;
00304 vector<Photon> directPhotons;
00305 vector<RadiancePhoton> radiancePhotons;
00306 causticPhotons.reserve(nCausticPhotons);
00307 indirectPhotons.reserve(nIndirectPhotons);
00308
00309 static StatsCounter nshot("Photon Map",
00310 "Number of photons shot from lights");
00311 bool causticDone = (nCausticPhotons == 0);
00312 bool indirectDone = (nIndirectPhotons == 0);
00313
00314
00315 int nLights = int(scene->lights.size());
00316 float *lightPower = (float *)alloca(nLights * sizeof(float));
00317 float *lightCDF = (float *)alloca((nLights+1) * sizeof(float));
00318 for (int i = 0; i < nLights; ++i)
00319 lightPower[i] = scene->lights[i]->Power(scene).y();
00320 float totalPower;
00321 ComputeStep1dCDF(lightPower, nLights, &totalPower, lightCDF);
00322
00323 vector<Spectrum> rpReflectances, rpTransmittances;
00324
00325 while (!causticDone || !indirectDone) {
00326 ++nshot;
00327
00328 if (nshot > 500000 &&
00329 (unsuccessful(nCausticPhotons,
00330 causticPhotons.size(),
00331 nshot) ||
00332 unsuccessful(nIndirectPhotons,
00333 indirectPhotons.size(),
00334 nshot))) {
00335 Error("Unable to store enough photons. Giving up.\n");
00336 return;
00337 }
00338
00339
00340 float u[4];
00341 u[0] = RadicalInverse((int)nshot+1, 2);
00342 u[1] = RadicalInverse((int)nshot+1, 3);
00343 u[2] = RadicalInverse((int)nshot+1, 5);
00344 u[3] = RadicalInverse((int)nshot+1, 7);
00345
00346
00347 float lightPdf;
00348 float uln = RadicalInverse((int)nshot+1, 11);
00349 int lightNum = Floor2Int(SampleStep1d(lightPower, lightCDF,
00350 totalPower, nLights, uln, &lightPdf) * nLights);
00351 lightNum = min(lightNum, nLights-1);
00352 const Light *light = scene->lights[lightNum];
00353
00354 RayDifferential photonRay;
00355 float pdf;
00356 Spectrum alpha = light->Sample_L(scene, u[0], u[1], u[2], u[3],
00357 &photonRay, &pdf);
00358 if (pdf == 0.f || alpha.Black()) continue;
00359 alpha /= pdf * lightPdf;
00360
00361 if (!alpha.Black()) {
00362
00363 bool specularPath = false;
00364 Intersection photonIsect;
00365 int nIntersections = 0;
00366 while (scene->Intersect(photonRay, &photonIsect)) {
00367 ++nIntersections;
00368
00369 alpha *= scene->Transmittance(photonRay);
00370 Vector wo = -photonRay.d;
00371 BSDF *photonBSDF = photonIsect.GetBSDF(photonRay);
00372 BxDFType specularType = BxDFType(BSDF_REFLECTION |
00373 BSDF_TRANSMISSION | BSDF_SPECULAR);
00374 bool hasNonSpecular = (photonBSDF->NumComponents() >
00375 photonBSDF->NumComponents(specularType));
00376 if (hasNonSpecular) {
00377
00378 Photon photon(photonIsect.dg.p, alpha, wo);
00379 if (nIntersections == 1) {
00380
00381 directPhotons.push_back(photon);
00382 }
00383 else {
00384
00385 if (specularPath) {
00386
00387 if (!causticDone) {
00388 causticPhotons.push_back(photon);
00389 if (causticPhotons.size() == nCausticPhotons) {
00390 causticDone = true;
00391 nCausticPaths = (int)nshot;
00392 causticMap = new KdTree<Photon, PhotonProcess>(causticPhotons);
00393 }
00394 progress.Update();
00395 }
00396 }
00397 else {
00398
00399 if (!indirectDone) {
00400 indirectPhotons.push_back(photon);
00401 if (indirectPhotons.size() == nIndirectPhotons) {
00402 indirectDone = true;
00403 nIndirectPaths = (int)nshot;
00404 indirectMap = new KdTree<Photon, PhotonProcess>(indirectPhotons);
00405 }
00406 progress.Update();
00407 }
00408 }
00409 }
00410 if (finalGather && RandomFloat() < .125f) {
00411
00412 static StatsCounter rp("Photon Map", "Radiance photons created");
00413 ++rp;
00414 Normal n = photonIsect.dg.nn;
00415 if (Dot(n, photonRay.d) > 0.f) n = -n;
00416 radiancePhotons.push_back(RadiancePhoton(photonIsect.dg.p, n));
00417 Spectrum rho_r = photonBSDF->rho(BSDF_ALL_REFLECTION);
00418 rpReflectances.push_back(rho_r);
00419 Spectrum rho_t = photonBSDF->rho(BSDF_ALL_TRANSMISSION);
00420 rpTransmittances.push_back(rho_t);
00421 }
00422 }
00423
00424 Vector wi;
00425 float pdf;
00426 BxDFType flags;
00427
00428 float u1, u2, u3;
00429 if (nIntersections == 1) {
00430 u1 = RadicalInverse((int)nshot+1, 13);
00431 u2 = RadicalInverse((int)nshot+1, 17);
00432 u3 = RadicalInverse((int)nshot+1, 19);
00433 }
00434 else {
00435 u1 = RandomFloat();
00436 u2 = RandomFloat();
00437 u3 = RandomFloat();
00438 }
00439
00440
00441 Spectrum fr = photonBSDF->Sample_f(wo, &wi, u1, u2, u3,
00442 &pdf, BSDF_ALL, &flags);
00443 if (fr.Black() || pdf == 0.f)
00444 break;
00445 Spectrum anew = alpha * fr *
00446 AbsDot(wi, photonBSDF->dgShading.nn) / pdf;
00447 float continueProb = min(1.f, anew.y() / alpha.y());
00448 if (RandomFloat() > continueProb || nIntersections > 10)
00449 break;
00450 alpha = anew / continueProb;
00451 specularPath = (nIntersections == 1 || specularPath) &&
00452 ((flags & BSDF_SPECULAR) != 0);
00453 photonRay = RayDifferential(photonIsect.dg.p, wi);
00454 }
00455 }
00456 BSDF::FreeAll();
00457 }
00458
00459 progress.Done();
00460
00461
00462 KdTree<Photon, PhotonProcess> directMap(directPhotons);
00463 int nDirectPaths = nshot;
00464 if (finalGather) {
00465 ProgressReporter p2(radiancePhotons.size(), "Computing photon radiances");
00466 for (u_int i = 0; i < radiancePhotons.size(); ++i) {
00467
00468 RadiancePhoton &rp = radiancePhotons[i];
00469 const Spectrum &rho_r = rpReflectances[i];
00470 const Spectrum &rho_t = rpTransmittances[i];
00471 Spectrum E;
00472 Point p = rp.p;
00473 Normal n = rp.n;
00474 if (!rho_r.Black()) {
00475 E = estimateE(&directMap, nDirectPaths, p, n) +
00476 estimateE(indirectMap, nIndirectPaths, p, n) +
00477 estimateE(causticMap, nCausticPaths, p, n);
00478 rp.Lo += E * INV_PI * rho_r;
00479 }
00480 if (!rho_t.Black()) {
00481 E = estimateE(&directMap, nDirectPaths, p, -n) +
00482 estimateE(indirectMap, nIndirectPaths, p, -n) +
00483 estimateE(causticMap, nCausticPaths, p, -n);
00484 rp.Lo += E * INV_PI * rho_t;
00485 }
00486 p2.Update();
00487 }
00488 radianceMap = new KdTree<RadiancePhoton,
00489 RadiancePhotonProcess>(radiancePhotons);
00490 p2.Done();
00491 }
00492 }
00493
00494 Spectrum ExPhotonIntegrator::Li(const Scene *scene,
00495 const RayDifferential &ray, const Sample *sample,
00496 float *alpha) const {
00497
00498 Spectrum L(0.);
00499 Intersection isect;
00500 if (scene->Intersect(ray, &isect)) {
00501 if (alpha) *alpha = 1.;
00502 Vector wo = -ray.d;
00503
00504 L += isect.Le(wo);
00505
00506 BSDF *bsdf = isect.GetBSDF(ray);
00507 const Point &p = bsdf->dgShading.p;
00508 const Normal &n = bsdf->dgShading.nn;
00509 L += UniformSampleAllLights(scene, p, n,
00510 wo, bsdf, sample,
00511 lightSampleOffset, bsdfSampleOffset,
00512 bsdfComponentOffset);
00513
00514
00515 L += LPhoton(causticMap, nCausticPaths, nLookup, bsdf,
00516 isect, wo, maxDistSquared);
00517 if (finalGather) {
00518 #if 1
00519
00520 BxDFType nonSpecular = BxDFType(BSDF_REFLECTION |
00521 BSDF_TRANSMISSION | BSDF_DIFFUSE | BSDF_GLOSSY);
00522 if (bsdf->NumComponents(nonSpecular) > 0) {
00523
00524 u_int nIndirSamplePhotons = 50;
00525 PhotonProcess proc(nIndirSamplePhotons, p);
00526 proc.photons = (ClosePhoton *)alloca(nIndirSamplePhotons *
00527 sizeof(ClosePhoton));
00528 float searchDist2 = maxDistSquared;
00529 while (proc.foundPhotons < nIndirSamplePhotons) {
00530 float md2 = searchDist2;
00531 proc.foundPhotons = 0;
00532 indirectMap->Lookup(p, proc, md2);
00533 searchDist2 *= 2.f;
00534 }
00535
00536 Vector *photonDirs = (Vector *)alloca(nIndirSamplePhotons *
00537 sizeof(Vector));
00538 for (u_int i = 0; i < nIndirSamplePhotons; ++i)
00539 photonDirs[i] = proc.photons[i].photon->wi;
00540
00541 Spectrum Li = 0.;
00542 static StatsCounter gatherRays("Photon Map",
00543 "Final gather rays traced");
00544 for (int i = 0; i < gatherSamples; ++i) {
00545
00546 Vector wi;
00547 float u1 = sample->twoD[gatherSampleOffset[0]][2*i];
00548 float u2 = sample->twoD[gatherSampleOffset[0]][2*i+1];
00549 float u3 = sample->oneD[gatherComponentOffset[0]][i];
00550 float pdf;
00551 Spectrum fr = bsdf->Sample_f(wo, &wi, u1, u2, u3,
00552 &pdf, BxDFType(BSDF_ALL & (~BSDF_SPECULAR)));
00553 if (fr.Black() || pdf == 0.f) continue;
00554
00555 RayDifferential bounceRay(p, wi);
00556 ++gatherRays;
00557 Intersection gatherIsect;
00558 if (scene->Intersect(bounceRay, &gatherIsect)) {
00559
00560 Spectrum Lindir = 0.f;
00561 Normal ng = gatherIsect.dg.nn;
00562 if (Dot(ng, bounceRay.d) > 0) ng = -ng;
00563 RadiancePhotonProcess proc(gatherIsect.dg.p, ng);
00564 float md2 = INFINITY;
00565 radianceMap->Lookup(gatherIsect.dg.p, proc, md2);
00566 if (proc.photon)
00567 Lindir = proc.photon->Lo;
00568 Lindir *= scene->Transmittance(bounceRay);
00569
00570
00571 float photonPdf = 0.f;
00572 float conePdf = UniformConePdf(cosGatherAngle);
00573 for (u_int j = 0; j < nIndirSamplePhotons; ++j)
00574 if (Dot(photonDirs[j], wi) > .999f * cosGatherAngle)
00575 photonPdf += conePdf;
00576 photonPdf /= nIndirSamplePhotons;
00577 float wt = PowerHeuristic(gatherSamples, pdf,
00578 gatherSamples, photonPdf);
00579 Li += fr * Lindir * AbsDot(wi, n) * wt / pdf;
00580 }
00581 }
00582 L += Li / gatherSamples;
00583
00584 Li = 0.;
00585 for (int i = 0; i < gatherSamples; ++i) {
00586
00587 float u1 = sample->oneD[gatherComponentOffset[1]][i];
00588 float u2 = sample->twoD[gatherSampleOffset[1]][2*i];
00589 float u3 = sample->twoD[gatherSampleOffset[1]][2*i+1];
00590 int photonNum = min((int)nIndirSamplePhotons - 1,
00591 Floor2Int(u1 * nIndirSamplePhotons));
00592
00593 Vector vx, vy;
00594 CoordinateSystem(photonDirs[photonNum], &vx, &vy);
00595 Vector wi = UniformSampleCone(u2, u3, cosGatherAngle, vx, vy,
00596 photonDirs[photonNum]);
00597
00598 Spectrum fr = bsdf->f(wo, wi);
00599 if (fr.Black()) continue;
00600
00601 float photonPdf = 0.f;
00602 float conePdf = UniformConePdf(cosGatherAngle);
00603 for (u_int j = 0; j < nIndirSamplePhotons; ++j)
00604 if (Dot(photonDirs[j], wi) > .999f * cosGatherAngle)
00605 photonPdf += conePdf;
00606 photonPdf /= nIndirSamplePhotons;
00607 RayDifferential bounceRay(p, wi);
00608 ++gatherRays;
00609 Intersection gatherIsect;
00610 if (scene->Intersect(bounceRay, &gatherIsect)) {
00611
00612 Spectrum Lindir = 0.f;
00613 Normal ng = gatherIsect.dg.nn;
00614 if (Dot(ng, bounceRay.d) > 0) ng = -ng;
00615 RadiancePhotonProcess proc(gatherIsect.dg.p, ng);
00616 float md2 = INFINITY;
00617 radianceMap->Lookup(gatherIsect.dg.p, proc, md2);
00618 if (proc.photon)
00619 Lindir = proc.photon->Lo;
00620 Lindir *= scene->Transmittance(bounceRay);
00621
00622 float bsdfPdf = bsdf->Pdf(wo, wi);
00623 float wt = PowerHeuristic(gatherSamples, photonPdf,
00624 gatherSamples, bsdfPdf);
00625 Li += fr * Lindir * AbsDot(wi, n) * wt / photonPdf;
00626 }
00627 }
00628 L += Li / gatherSamples;
00629 }
00630 #else
00631
00632 Normal nn = n;
00633 if (Dot(nn, ray.d) > 0.) nn = -n;
00634 RadiancePhotonProcess proc(p, nn);
00635 float md2 = INFINITY;
00636 radianceMap->Lookup(p, proc, md2);
00637 if (proc.photon)
00638 L += proc.photon->Lo;
00639 #endif
00640
00641 }
00642 else {
00643 L += LPhoton(indirectMap, nIndirectPaths, nLookup,
00644 bsdf, isect, wo, maxDistSquared);
00645 }
00646 if (specularDepth++ < maxSpecularDepth) {
00647 Vector wi;
00648
00649 Spectrum f = bsdf->Sample_f(wo, &wi,
00650 BxDFType(BSDF_REFLECTION | BSDF_SPECULAR));
00651 if (!f.Black()) {
00652
00653 RayDifferential rd(p, wi);
00654 rd.hasDifferentials = true;
00655 rd.rx.o = p + isect.dg.dpdx;
00656 rd.ry.o = p + isect.dg.dpdy;
00657
00658 Normal dndx = bsdf->dgShading.dndu * bsdf->dgShading.dudx +
00659 bsdf->dgShading.dndv * bsdf->dgShading.dvdx;
00660 Normal dndy = bsdf->dgShading.dndu * bsdf->dgShading.dudy +
00661 bsdf->dgShading.dndv * bsdf->dgShading.dvdy;
00662 Vector dwodx = -ray.rx.d - wo, dwody = -ray.ry.d - wo;
00663 float dDNdx = Dot(dwodx, n) + Dot(wo, dndx);
00664 float dDNdy = Dot(dwody, n) + Dot(wo, dndy);
00665 rd.rx.d = wi -
00666 dwodx + 2 * Vector(Dot(wo, n) * dndx +
00667 dDNdx * n);
00668 rd.ry.d = wi -
00669 dwody + 2 * Vector(Dot(wo, n) * dndy +
00670 dDNdy * n);
00671 L += scene->Li(rd, sample) * f * AbsDot(wi, n);
00672 }
00673 f = bsdf->Sample_f(wo, &wi,
00674 BxDFType(BSDF_TRANSMISSION | BSDF_SPECULAR));
00675 if (!f.Black()) {
00676
00677 RayDifferential rd(p, wi);
00678 rd.hasDifferentials = true;
00679 rd.rx.o = p + isect.dg.dpdx;
00680 rd.ry.o = p + isect.dg.dpdy;
00681
00682 float eta = bsdf->eta;
00683 Vector w = -wo;
00684 if (Dot(wo, n) < 0) eta = 1.f / eta;
00685
00686 Normal dndx = bsdf->dgShading.dndu * bsdf->dgShading.dudx + bsdf->dgShading.dndv * bsdf->dgShading.dvdx;
00687 Normal dndy = bsdf->dgShading.dndu * bsdf->dgShading.dudy + bsdf->dgShading.dndv * bsdf->dgShading.dvdy;
00688
00689 Vector dwodx = -ray.rx.d - wo, dwody = -ray.ry.d - wo;
00690 float dDNdx = Dot(dwodx, n) + Dot(wo, dndx);
00691 float dDNdy = Dot(dwody, n) + Dot(wo, dndy);
00692
00693 float mu = eta * Dot(w, n) - Dot(wi, n);
00694 float dmudx = (eta - (eta*eta*Dot(w,n))/Dot(wi, n)) * dDNdx;
00695 float dmudy = (eta - (eta*eta*Dot(w,n))/Dot(wi, n)) * dDNdy;
00696
00697 rd.rx.d = wi + eta * dwodx - Vector(mu * dndx + dmudx * n);
00698 rd.ry.d = wi + eta * dwody - Vector(mu * dndy + dmudy * n);
00699 L += scene->Li(rd, sample) * f * AbsDot(wi, n);
00700 }
00701 }
00702 --specularDepth;
00703 }
00704 else {
00705
00706 if (alpha) *alpha = 0.;
00707 for (u_int i = 0; i < scene->lights.size(); ++i)
00708 L += scene->lights[i]->Le(ray);
00709 if (alpha && !L.Black()) *alpha = 1.;
00710 return L;
00711 }
00712 return L;
00713 }
00714 extern "C" DLLEXPORT SurfaceIntegrator *CreateSurfaceIntegrator(const ParamSet ¶ms) {
00715 int nCaustic = params.FindOneInt("causticphotons", 20000);
00716 int nIndirect = params.FindOneInt("indirectphotons", 100000);
00717 int nUsed = params.FindOneInt("nused", 50);
00718 int maxDepth = params.FindOneInt("maxdepth", 5);
00719 bool finalGather = params.FindOneBool("finalgather", true);
00720 int gatherSamples = params.FindOneInt("finalgathersamples", 32);
00721 float maxDist = params.FindOneFloat("maxdist", .1f);
00722 float rrTreshold = params.FindOneFloat("rrthreshold", .05f);
00723 float gatherAngle = params.FindOneFloat("gatherangle", 10.f);
00724 return new ExPhotonIntegrator(nCaustic, nIndirect,
00725 nUsed, maxDepth, maxDist, finalGather, gatherSamples,
00726 rrTreshold, gatherAngle);
00727 }