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 "integrators/photonmap.h"
00028 #include "scene.h"
00029 #include "montecarlo.h"
00030 #include "sampler.h"
00031 #include "progressreporter.h"
00032 #include "intersection.h"
00033 #include "paramset.h"
00034 #include "camera.h"
00035
00036
00037
00038 struct Photon {
00039 Photon(const Point &pp, const Spectrum &wt, const Vector &w)
00040 : p(pp), alpha(wt), wi(w) { }
00041 Photon() { }
00042 Point p;
00043 Spectrum alpha;
00044 Vector wi;
00045 };
00046
00047
00048 class PhotonShootingTask : public Task {
00049 public:
00050 PhotonShootingTask(int tn, float ti, Mutex &m, PhotonIntegrator *in,
00051 ProgressReporter &prog, bool &at, int &ndp,
00052 vector<Photon> &direct, vector<Photon> &indir, vector<Photon> &caustic,
00053 vector<RadiancePhoton> &rps, vector<Spectrum> &rpR, vector<Spectrum> &rpT,
00054 uint32_t &ns, Distribution1D *distrib, const Scene *sc,
00055 const Renderer *sr)
00056 : taskNum(tn), time(ti), mutex(m), integrator(in), progress(prog),
00057 abortTasks(at), nDirectPaths(ndp),
00058 directPhotons(direct), indirectPhotons(indir), causticPhotons(caustic),
00059 radiancePhotons(rps), rpReflectances(rpR), rpTransmittances(rpT),
00060 nshot(ns), lightDistribution(distrib), scene(sc), renderer (sr) { }
00061 void Run();
00062
00063 int taskNum;
00064 float time;
00065 Mutex &mutex;
00066 PhotonIntegrator *integrator;
00067 ProgressReporter &progress;
00068 bool &abortTasks;
00069 int &nDirectPaths;
00070 vector<Photon> &directPhotons, &indirectPhotons, &causticPhotons;
00071 vector<RadiancePhoton> &radiancePhotons;
00072 vector<Spectrum> &rpReflectances, &rpTransmittances;
00073 uint32_t &nshot;
00074 const Distribution1D *lightDistribution;
00075 const Scene *scene;
00076 const Renderer *renderer;
00077 };
00078
00079
00080 struct RadiancePhoton {
00081 RadiancePhoton(const Point &pp, const Normal &nn)
00082 : p(pp), n(nn), Lo(0.f) { }
00083 RadiancePhoton() { }
00084 Point p;
00085 Normal n;
00086 Spectrum Lo;
00087 };
00088
00089
00090 class ComputeRadianceTask : public Task {
00091 public:
00092 ComputeRadianceTask(ProgressReporter &prog, uint32_t tn, uint32_t nt,
00093 vector<RadiancePhoton> &rps, const vector<Spectrum> &rhor,
00094 const vector<Spectrum> &rhot,
00095 uint32_t nlookup, float md2,
00096 int ndirect, KdTree<Photon> *direct,
00097 int nindirect, KdTree<Photon> *indirect,
00098 int ncaus, KdTree<Photon> *caustic)
00099 : progress(prog), taskNum(tn), numTasks(nt), radiancePhotons(rps),
00100 rpReflectances(rhor), rpTransmittances(rhot), nLookup(nlookup),
00101 maxDistSquared(md2),
00102 nDirectPaths(ndirect), nIndirectPaths(nindirect), nCausticPaths(ncaus),
00103 directMap(direct), indirectMap(indirect), causticMap(caustic) { }
00104 void Run();
00105
00106 private:
00107 ProgressReporter &progress;
00108 uint32_t taskNum, numTasks;
00109 vector<RadiancePhoton> &radiancePhotons;
00110 const vector<Spectrum> &rpReflectances, &rpTransmittances;
00111 uint32_t nLookup;
00112 float maxDistSquared;
00113 int nDirectPaths, nIndirectPaths, nCausticPaths;
00114 KdTree<Photon> *directMap, *indirectMap, *causticMap;
00115 };
00116
00117
00118 struct PhotonProcess {
00119
00120 PhotonProcess(uint32_t mp, ClosePhoton *buf);
00121 void operator()(const Point &p, const Photon &photon, float dist2,
00122 float &maxDistSquared);
00123 ClosePhoton *photons;
00124 uint32_t nLookup, nFound;
00125 };
00126
00127
00128 struct ClosePhoton {
00129
00130 ClosePhoton(const Photon *p = NULL, float md2 = INFINITY)
00131 : photon(p), distanceSquared(md2) { }
00132 bool operator<(const ClosePhoton &p2) const {
00133 return distanceSquared == p2.distanceSquared ?
00134 (photon < p2.photon) : (distanceSquared < p2.distanceSquared);
00135 }
00136 const Photon *photon;
00137 float distanceSquared;
00138 };
00139
00140
00141 PhotonProcess::PhotonProcess(uint32_t mp, ClosePhoton *buf) {
00142 photons = buf;
00143 nLookup = mp;
00144 nFound = 0;
00145 }
00146
00147
00148 struct RadiancePhotonProcess {
00149
00150 RadiancePhotonProcess(const Normal &nn)
00151 : n(nn) {
00152 photon = NULL;
00153 }
00154 void operator()(const Point &p, const RadiancePhoton &rp,
00155 float distSquared, float &maxDistSquared) {
00156 if (Dot(rp.n, n) > 0) {
00157 photon = &rp;
00158 maxDistSquared = distSquared;
00159 }
00160 }
00161 const Normal &n;
00162 const RadiancePhoton *photon;
00163 };
00164
00165
00166 inline float kernel(const Photon *photon, const Point &p, float maxDist2);
00167 static Spectrum LPhoton(KdTree<Photon> *map, int nPaths, int nLookup,
00168 ClosePhoton *lookupBuf, BSDF *bsdf, RNG &rng, const Intersection &isect,
00169 const Vector &w, float maxDistSquared);
00170 static Spectrum EPhoton(KdTree<Photon> *map, int count, int nLookup,
00171 ClosePhoton *lookupBuf, float maxDist2, const Point &p, const Normal &n);
00172
00173
00174 inline bool unsuccessful(uint32_t needed, uint32_t found, uint32_t shot) {
00175 return (found < needed && (found == 0 || found < shot / 1024));
00176 }
00177
00178
00179 inline void PhotonProcess::operator()(const Point &p,
00180 const Photon &photon, float distSquared, float &maxDistSquared) {
00181 if (nFound < nLookup) {
00182
00183 photons[nFound++] = ClosePhoton(&photon, distSquared);
00184 if (nFound == nLookup) {
00185 std::make_heap(&photons[0], &photons[nLookup]);
00186 maxDistSquared = photons[0].distanceSquared;
00187 }
00188 }
00189 else {
00190
00191 std::pop_heap(&photons[0], &photons[nLookup]);
00192 photons[nLookup-1] = ClosePhoton(&photon, distSquared);
00193 std::push_heap(&photons[0], &photons[nLookup]);
00194 maxDistSquared = photons[0].distanceSquared;
00195 }
00196 }
00197
00198
00199 inline float kernel(const Photon *photon, const Point &p,
00200 float maxDist2) {
00201 float s = (1.f - DistanceSquared(photon->p, p) / maxDist2);
00202 return 3.f * INV_PI * s * s;
00203 }
00204
00205
00206 Spectrum LPhoton(KdTree<Photon> *map, int nPaths, int nLookup,
00207 ClosePhoton *lookupBuf, BSDF *bsdf, RNG &rng,
00208 const Intersection &isect, const Vector &wo, float maxDist2) {
00209 Spectrum L(0.);
00210 BxDFType nonSpecular = BxDFType(BSDF_REFLECTION |
00211 BSDF_TRANSMISSION | BSDF_DIFFUSE | BSDF_GLOSSY);
00212 if (map && bsdf->NumComponents(nonSpecular) > 0) {
00213 PBRT_PHOTON_MAP_STARTED_LOOKUP(const_cast<DifferentialGeometry *>(&isect.dg));
00214
00215 PhotonProcess proc(nLookup, lookupBuf);
00216 map->Lookup(isect.dg.p, proc, maxDist2);
00217
00218
00219 ClosePhoton *photons = proc.photons;
00220 int nFound = proc.nFound;
00221 Normal Nf = Faceforward(bsdf->dgShading.nn, wo);
00222 if (bsdf->NumComponents(BxDFType(BSDF_REFLECTION |
00223 BSDF_TRANSMISSION | BSDF_GLOSSY)) > 0) {
00224
00225 for (int i = 0; i < nFound; ++i) {
00226 const Photon *p = photons[i].photon;
00227 float k = kernel(p, isect.dg.p, maxDist2);
00228 L += (k / (nPaths * maxDist2)) * bsdf->f(wo, p->wi) *
00229 p->alpha;
00230 }
00231 }
00232 else {
00233
00234 Spectrum Lr(0.), Lt(0.);
00235 for (int i = 0; i < nFound; ++i) {
00236 if (Dot(Nf, photons[i].photon->wi) > 0.f) {
00237 float k = kernel(photons[i].photon, isect.dg.p, maxDist2);
00238 Lr += (k / (nPaths * maxDist2)) * photons[i].photon->alpha;
00239 }
00240 else {
00241 float k = kernel(photons[i].photon, isect.dg.p, maxDist2);
00242 Lt += (k / (nPaths * maxDist2)) * photons[i].photon->alpha;
00243 }
00244 }
00245 L += Lr * bsdf->rho(wo, rng, BSDF_ALL_REFLECTION) * INV_PI +
00246 Lt * bsdf->rho(wo, rng, BSDF_ALL_TRANSMISSION) * INV_PI;
00247 }
00248 PBRT_PHOTON_MAP_FINISHED_LOOKUP(const_cast<DifferentialGeometry *>(&isect.dg),
00249 proc.nFound, proc.nLookup, &L);
00250 }
00251 return L;
00252 }
00253
00254
00255 Spectrum EPhoton(KdTree<Photon> *map, int count, int nLookup,
00256 ClosePhoton *lookupBuf, float maxDist2, const Point &p,
00257 const Normal &n) {
00258 if (!map) return 0.f;
00259
00260 PhotonProcess proc(nLookup, lookupBuf);
00261 float md2 = maxDist2;
00262 map->Lookup(p, proc, md2);
00263 Assert(md2 > 0.f);
00264
00265
00266 if (proc.nFound == 0) return Spectrum(0.f);
00267 ClosePhoton *photons = proc.photons;
00268 Spectrum E(0.);
00269 for (uint32_t i = 0; i < proc.nFound; ++i)
00270 if (Dot(n, photons[i].photon->wi) > 0.)
00271 E += photons[i].photon->alpha;
00272 return E / (count * md2 * M_PI);
00273 }
00274
00275
00276
00277
00278 PhotonIntegrator::PhotonIntegrator(int ncaus, int nind,
00279 int nl, int mdepth, int mphodepth, float mdist, bool fg,
00280 int gs, float ga) {
00281 nCausticPhotonsWanted = ncaus;
00282 nIndirectPhotonsWanted = nind;
00283 nLookup = nl;
00284 maxSpecularDepth = mdepth;
00285 maxPhotonDepth = mphodepth;
00286 maxDistSquared = mdist * mdist;
00287 finalGather = fg;
00288 cosGatherAngle = cos(Radians(ga));
00289 gatherSamples = gs;
00290 nCausticPaths = nIndirectPaths = 0;
00291 causticMap = indirectMap = NULL;
00292 radianceMap = NULL;
00293 lightSampleOffsets = NULL;
00294 bsdfSampleOffsets = NULL;
00295 }
00296
00297
00298 PhotonIntegrator::~PhotonIntegrator() {
00299 delete[] lightSampleOffsets;
00300 delete[] bsdfSampleOffsets;
00301 delete causticMap;
00302 delete indirectMap;
00303 delete radianceMap;
00304 }
00305
00306
00307 void PhotonIntegrator::RequestSamples(Sampler *sampler, Sample *sample,
00308 const Scene *scene) {
00309
00310 uint32_t nLights = scene->lights.size();
00311 lightSampleOffsets = new LightSampleOffsets[nLights];
00312 bsdfSampleOffsets = new BSDFSampleOffsets[nLights];
00313 for (uint32_t i = 0; i < nLights; ++i) {
00314 const Light *light = scene->lights[i];
00315 int nSamples = light->nSamples;
00316 if (sampler) nSamples = sampler->RoundSize(nSamples);
00317 lightSampleOffsets[i] = LightSampleOffsets(nSamples, sample);
00318 bsdfSampleOffsets[i] = BSDFSampleOffsets(nSamples, sample);
00319 }
00320
00321
00322 if (finalGather) {
00323 gatherSamples = max(1, gatherSamples/2);
00324 if (sampler) gatherSamples = sampler->RoundSize(gatherSamples);
00325 bsdfGatherSampleOffsets = BSDFSampleOffsets(gatherSamples, sample);
00326 indirGatherSampleOffsets = BSDFSampleOffsets(gatherSamples, sample);
00327 }
00328 }
00329
00330
00331 void PhotonIntegrator::Preprocess(const Scene *scene,
00332 const Camera *camera, const Renderer *renderer) {
00333 if (scene->lights.size() == 0) return;
00334
00335 Mutex *mutex = Mutex::Create();
00336 int nDirectPaths = 0;
00337 vector<Photon> causticPhotons, directPhotons, indirectPhotons;
00338 vector<RadiancePhoton> radiancePhotons;
00339 bool abortTasks = false;
00340 causticPhotons.reserve(nCausticPhotonsWanted);
00341 indirectPhotons.reserve(nIndirectPhotonsWanted);
00342 uint32_t nshot = 0;
00343 vector<Spectrum> rpReflectances, rpTransmittances;
00344
00345
00346 Distribution1D *lightDistribution = ComputeLightSamplingCDF(scene);
00347
00348
00349 ProgressReporter progress(nCausticPhotonsWanted+nIndirectPhotonsWanted, "Shooting photons");
00350 vector<Task *> photonShootingTasks;
00351 int nTasks = NumSystemCores();
00352 for (int i = 0; i < nTasks; ++i)
00353 photonShootingTasks.push_back(new PhotonShootingTask(
00354 i, camera ? camera->shutterOpen : 0.f, *mutex, this, progress, abortTasks, nDirectPaths,
00355 directPhotons, indirectPhotons, causticPhotons, radiancePhotons,
00356 rpReflectances, rpTransmittances,
00357 nshot, lightDistribution, scene, renderer));
00358 EnqueueTasks(photonShootingTasks);
00359 WaitForAllTasks();
00360 for (uint32_t i = 0; i < photonShootingTasks.size(); ++i)
00361 delete photonShootingTasks[i];
00362 Mutex::Destroy(mutex);
00363 progress.Done();
00364
00365
00366 KdTree<Photon> *directMap = NULL;
00367 if (directPhotons.size() > 0)
00368 directMap = new KdTree<Photon>(directPhotons);
00369 if (causticPhotons.size() > 0)
00370 causticMap = new KdTree<Photon>(causticPhotons);
00371 if (indirectPhotons.size() > 0)
00372 indirectMap = new KdTree<Photon>(indirectPhotons);
00373
00374
00375 if (finalGather && radiancePhotons.size()) {
00376
00377 vector<Task *> radianceTasks;
00378 uint32_t numTasks = 64;
00379 ProgressReporter progRadiance(numTasks, "Computing photon radiances");
00380 for (uint32_t i = 0; i < numTasks; ++i)
00381 radianceTasks.push_back(new ComputeRadianceTask(progRadiance,
00382 i, numTasks, radiancePhotons, rpReflectances, rpTransmittances,
00383 nLookup, maxDistSquared, nDirectPaths, directMap,
00384 nIndirectPaths, indirectMap,
00385 nCausticPaths, causticMap));
00386 EnqueueTasks(radianceTasks);
00387 WaitForAllTasks();
00388 for (uint32_t i = 0; i < radianceTasks.size(); ++i)
00389 delete radianceTasks[i];
00390 progRadiance.Done();
00391 radianceMap = new KdTree<RadiancePhoton>(radiancePhotons);
00392 }
00393 delete directMap;
00394 }
00395
00396
00397 void PhotonShootingTask::Run() {
00398
00399 MemoryArena arena;
00400 RNG rng(31 * taskNum);
00401 vector<Photon> localDirectPhotons, localIndirectPhotons, localCausticPhotons;
00402 vector<RadiancePhoton> localRadiancePhotons;
00403 uint32_t totalPaths = 0;
00404 bool causticDone = (integrator->nCausticPhotonsWanted == 0);
00405 bool indirectDone = (integrator->nIndirectPhotonsWanted == 0);
00406 PermutedHalton halton(6, rng);
00407 vector<Spectrum> localRpReflectances, localRpTransmittances;
00408 while (true) {
00409
00410 const uint32_t blockSize = 4096;
00411 for (uint32_t i = 0; i < blockSize; ++i) {
00412 float u[6];
00413 halton.Sample(++totalPaths, u);
00414
00415 float lightPdf;
00416 int lightNum = lightDistribution->SampleDiscrete(u[0], &lightPdf);
00417 const Light *light = scene->lights[lightNum];
00418
00419
00420 RayDifferential photonRay;
00421 float pdf;
00422 LightSample ls(u[1], u[2], u[3]);
00423 Normal Nl;
00424 Spectrum Le = light->Sample_L(scene, ls, u[4], u[5],
00425 time, &photonRay, &Nl, &pdf);
00426 if (pdf == 0.f || Le.IsBlack()) continue;
00427 Spectrum alpha = (AbsDot(Nl, photonRay.d) * Le) / (pdf * lightPdf);
00428 if (!alpha.IsBlack()) {
00429
00430 PBRT_PHOTON_MAP_STARTED_RAY_PATH(&photonRay, &alpha);
00431 bool specularPath = true;
00432 Intersection photonIsect;
00433 int nIntersections = 0;
00434 while (scene->Intersect(photonRay, &photonIsect)) {
00435 ++nIntersections;
00436
00437 alpha *= renderer->Transmittance(scene, photonRay, NULL, rng, arena);
00438 BSDF *photonBSDF = photonIsect.GetBSDF(photonRay, arena);
00439 BxDFType specularType = BxDFType(BSDF_REFLECTION |
00440 BSDF_TRANSMISSION | BSDF_SPECULAR);
00441 bool hasNonSpecular = (photonBSDF->NumComponents() >
00442 photonBSDF->NumComponents(specularType));
00443 Vector wo = -photonRay.d;
00444 if (hasNonSpecular) {
00445
00446 Photon photon(photonIsect.dg.p, alpha, wo);
00447 bool depositedPhoton = false;
00448 if (specularPath && nIntersections > 1) {
00449 if (!causticDone) {
00450 PBRT_PHOTON_MAP_DEPOSITED_CAUSTIC_PHOTON(&photonIsect.dg, &alpha, &wo);
00451 depositedPhoton = true;
00452 localCausticPhotons.push_back(photon);
00453 }
00454 }
00455 else {
00456
00457
00458
00459
00460 if (nIntersections == 1 && !indirectDone && integrator->finalGather) {
00461 PBRT_PHOTON_MAP_DEPOSITED_DIRECT_PHOTON(&photonIsect.dg, &alpha, &wo);
00462 depositedPhoton = true;
00463 localDirectPhotons.push_back(photon);
00464 }
00465 else if (nIntersections > 1 && !indirectDone) {
00466 PBRT_PHOTON_MAP_DEPOSITED_INDIRECT_PHOTON(&photonIsect.dg, &alpha, &wo);
00467 depositedPhoton = true;
00468 localIndirectPhotons.push_back(photon);
00469 }
00470 }
00471
00472
00473 if (depositedPhoton && integrator->finalGather &&
00474 rng.RandomFloat() < .125f) {
00475 Normal n = photonIsect.dg.nn;
00476 n = Faceforward(n, -photonRay.d);
00477 localRadiancePhotons.push_back(RadiancePhoton(photonIsect.dg.p, n));
00478 Spectrum rho_r = photonBSDF->rho(rng, BSDF_ALL_REFLECTION);
00479 localRpReflectances.push_back(rho_r);
00480 Spectrum rho_t = photonBSDF->rho(rng, BSDF_ALL_TRANSMISSION);
00481 localRpTransmittances.push_back(rho_t);
00482 }
00483 }
00484 if (nIntersections >= integrator->maxPhotonDepth) break;
00485
00486
00487 Vector wi;
00488 float pdf;
00489 BxDFType flags;
00490 Spectrum fr = photonBSDF->Sample_f(wo, &wi, BSDFSample(rng),
00491 &pdf, BSDF_ALL, &flags);
00492 if (fr.IsBlack() || pdf == 0.f) break;
00493 Spectrum anew = alpha * fr *
00494 AbsDot(wi, photonBSDF->dgShading.nn) / pdf;
00495
00496
00497 float continueProb = min(1.f, anew.y() / alpha.y());
00498 if (rng.RandomFloat() > continueProb)
00499 break;
00500 alpha = anew / continueProb;
00501 specularPath &= ((flags & BSDF_SPECULAR) != 0);
00502
00503 if (indirectDone && !specularPath) break;
00504 photonRay = RayDifferential(photonIsect.dg.p, wi, photonRay,
00505 photonIsect.rayEpsilon);
00506 }
00507 PBRT_PHOTON_MAP_FINISHED_RAY_PATH(&photonRay, &alpha);
00508 }
00509 arena.FreeAll();
00510 }
00511
00512
00513 { MutexLock lock(mutex);
00514
00515
00516 if (abortTasks)
00517 return;
00518 if (nshot > 500000 &&
00519 (unsuccessful(integrator->nCausticPhotonsWanted,
00520 causticPhotons.size(), blockSize) ||
00521 unsuccessful(integrator->nIndirectPhotonsWanted,
00522 indirectPhotons.size(), blockSize))) {
00523 Error("Unable to store enough photons. Giving up.\n");
00524 causticPhotons.erase(causticPhotons.begin(), causticPhotons.end());
00525 indirectPhotons.erase(indirectPhotons.begin(), indirectPhotons.end());
00526 radiancePhotons.erase(radiancePhotons.begin(), radiancePhotons.end());
00527 abortTasks = true;
00528 return;
00529 }
00530 progress.Update(localIndirectPhotons.size() + localCausticPhotons.size());
00531 nshot += blockSize;
00532
00533
00534 if (!indirectDone) {
00535 integrator->nIndirectPaths += blockSize;
00536 for (uint32_t i = 0; i < localIndirectPhotons.size(); ++i)
00537 indirectPhotons.push_back(localIndirectPhotons[i]);
00538 localIndirectPhotons.erase(localIndirectPhotons.begin(),
00539 localIndirectPhotons.end());
00540 if (indirectPhotons.size() >= integrator->nIndirectPhotonsWanted)
00541 indirectDone = true;
00542 nDirectPaths += blockSize;
00543 for (uint32_t i = 0; i < localDirectPhotons.size(); ++i)
00544 directPhotons.push_back(localDirectPhotons[i]);
00545 localDirectPhotons.erase(localDirectPhotons.begin(),
00546 localDirectPhotons.end());
00547 }
00548
00549
00550 if (!causticDone) {
00551 integrator->nCausticPaths += blockSize;
00552 for (uint32_t i = 0; i < localCausticPhotons.size(); ++i)
00553 causticPhotons.push_back(localCausticPhotons[i]);
00554 localCausticPhotons.erase(localCausticPhotons.begin(), localCausticPhotons.end());
00555 if (causticPhotons.size() >= integrator->nCausticPhotonsWanted)
00556 causticDone = true;
00557 }
00558
00559 for (uint32_t i = 0; i < localRadiancePhotons.size(); ++i)
00560 radiancePhotons.push_back(localRadiancePhotons[i]);
00561 localRadiancePhotons.erase(localRadiancePhotons.begin(), localRadiancePhotons.end());
00562 for (uint32_t i = 0; i < localRpReflectances.size(); ++i)
00563 rpReflectances.push_back(localRpReflectances[i]);
00564 localRpReflectances.erase(localRpReflectances.begin(), localRpReflectances.end());
00565 for (uint32_t i = 0; i < localRpTransmittances.size(); ++i)
00566 rpTransmittances.push_back(localRpTransmittances[i]);
00567 localRpTransmittances.erase(localRpTransmittances.begin(), localRpTransmittances.end());
00568 }
00569
00570
00571 if (indirectDone && causticDone)
00572 break;
00573 }
00574 }
00575
00576
00577 void ComputeRadianceTask::Run() {
00578
00579 uint32_t taskSize = radiancePhotons.size() / numTasks;
00580 uint32_t excess = radiancePhotons.size() % numTasks;
00581 uint32_t rpStart = min(taskNum, excess) * (taskSize+1) +
00582 max(0, (int)taskNum-(int)excess) * taskSize;
00583 uint32_t rpEnd = rpStart + taskSize + (taskNum < excess ? 1 : 0);
00584 if (taskNum == numTasks-1) Assert(rpEnd == radiancePhotons.size());
00585 ClosePhoton *lookupBuf = new ClosePhoton[nLookup];
00586 for (uint32_t i = rpStart; i < rpEnd; ++i) {
00587
00588 RadiancePhoton &rp = radiancePhotons[i];
00589 const Spectrum &rho_r = rpReflectances[i], &rho_t = rpTransmittances[i];
00590 if (!rho_r.IsBlack()) {
00591
00592 Spectrum E = EPhoton(directMap, nDirectPaths, nLookup, lookupBuf,
00593 maxDistSquared, rp.p, rp.n) +
00594 EPhoton(indirectMap, nIndirectPaths, nLookup, lookupBuf,
00595 maxDistSquared, rp.p, rp.n) +
00596 EPhoton(causticMap, nCausticPaths, nLookup, lookupBuf,
00597 maxDistSquared, rp.p, rp.n);
00598 rp.Lo += INV_PI * rho_r * E;
00599 }
00600 if (!rho_t.IsBlack()) {
00601
00602 Spectrum E = EPhoton(directMap, nDirectPaths, nLookup, lookupBuf,
00603 maxDistSquared, rp.p, -rp.n) +
00604 EPhoton(indirectMap, nIndirectPaths, nLookup, lookupBuf,
00605 maxDistSquared, rp.p, -rp.n) +
00606 EPhoton(causticMap, nCausticPaths, nLookup, lookupBuf,
00607 maxDistSquared, rp.p, -rp.n);
00608 rp.Lo += INV_PI * rho_t * E;
00609 }
00610 }
00611 delete[] lookupBuf;
00612 progress.Update();
00613 }
00614
00615
00616 Spectrum PhotonIntegrator::Li(const Scene *scene, const Renderer *renderer,
00617 const RayDifferential &ray, const Intersection &isect,
00618 const Sample *sample, RNG &rng, MemoryArena &arena) const {
00619 Spectrum L(0.);
00620 Vector wo = -ray.d;
00621
00622 L += isect.Le(wo);
00623
00624
00625 BSDF *bsdf = isect.GetBSDF(ray, arena);
00626 const Point &p = bsdf->dgShading.p;
00627 const Normal &n = bsdf->dgShading.nn;
00628 L += UniformSampleAllLights(scene, renderer, arena, p, n,
00629 wo, isect.rayEpsilon, ray.time, bsdf, sample, rng,
00630 lightSampleOffsets, bsdfSampleOffsets);
00631
00632 ClosePhoton *lookupBuf = arena.Alloc<ClosePhoton>(nLookup);
00633 L += LPhoton(causticMap, nCausticPaths, nLookup, lookupBuf, bsdf,
00634 rng, isect, wo, maxDistSquared);
00635
00636
00637 if (finalGather && indirectMap != NULL) {
00638 #if 1
00639
00640 BxDFType nonSpecular = BxDFType(BSDF_REFLECTION |
00641 BSDF_TRANSMISSION | BSDF_DIFFUSE | BSDF_GLOSSY);
00642 if (bsdf->NumComponents(nonSpecular) > 0) {
00643
00644 const uint32_t nIndirSamplePhotons = 50;
00645 PhotonProcess proc(nIndirSamplePhotons,
00646 arena.Alloc<ClosePhoton>(nIndirSamplePhotons));
00647 float searchDist2 = maxDistSquared;
00648 while (proc.nFound < nIndirSamplePhotons) {
00649 float md2 = searchDist2;
00650 proc.nFound = 0;
00651 indirectMap->Lookup(p, proc, md2);
00652 searchDist2 *= 2.f;
00653 }
00654
00655
00656 Vector *photonDirs = arena.Alloc<Vector>(nIndirSamplePhotons);
00657 for (uint32_t i = 0; i < nIndirSamplePhotons; ++i)
00658 photonDirs[i] = proc.photons[i].photon->wi;
00659
00660
00661 Spectrum Li = 0.;
00662 for (int i = 0; i < gatherSamples; ++i) {
00663
00664 Vector wi;
00665 float pdf;
00666 BSDFSample bsdfSample(sample, bsdfGatherSampleOffsets, i);
00667 Spectrum fr = bsdf->Sample_f(wo, &wi, bsdfSample,
00668 &pdf, BxDFType(BSDF_ALL & ~BSDF_SPECULAR));
00669 if (fr.IsBlack() || pdf == 0.f) continue;
00670 Assert(pdf >= 0.f);
00671
00672
00673 RayDifferential bounceRay(p, wi, ray, isect.rayEpsilon);
00674 Intersection gatherIsect;
00675 if (scene->Intersect(bounceRay, &gatherIsect)) {
00676
00677 Spectrum Lindir = 0.f;
00678 Normal nGather = gatherIsect.dg.nn;
00679 nGather = Faceforward(nGather, -bounceRay.d);
00680 RadiancePhotonProcess proc(nGather);
00681 float md2 = INFINITY;
00682 radianceMap->Lookup(gatherIsect.dg.p, proc, md2);
00683 if (proc.photon != NULL)
00684 Lindir = proc.photon->Lo;
00685 Lindir *= renderer->Transmittance(scene, bounceRay, NULL, rng, arena);
00686
00687
00688
00689
00690 float photonPdf = 0.f;
00691 float conePdf = UniformConePdf(cosGatherAngle);
00692 for (uint32_t j = 0; j < nIndirSamplePhotons; ++j)
00693 if (Dot(photonDirs[j], wi) > .999f * cosGatherAngle)
00694 photonPdf += conePdf;
00695 photonPdf /= nIndirSamplePhotons;
00696 float wt = PowerHeuristic(gatherSamples, pdf, gatherSamples, photonPdf);
00697 Li += fr * Lindir * (AbsDot(wi, n) * wt / pdf);
00698 }
00699 }
00700 L += Li / gatherSamples;
00701
00702
00703 Li = 0.;
00704 for (int i = 0; i < gatherSamples; ++i) {
00705
00706 BSDFSample gatherSample(sample, indirGatherSampleOffsets, i);
00707 int photonNum = min((int)nIndirSamplePhotons - 1,
00708 Floor2Int(gatherSample.uComponent * nIndirSamplePhotons));
00709
00710
00711 Vector vx, vy;
00712 CoordinateSystem(photonDirs[photonNum], &vx, &vy);
00713 Vector wi = UniformSampleCone(gatherSample.uDir[0], gatherSample.uDir[1],
00714 cosGatherAngle, vx, vy, photonDirs[photonNum]);
00715
00716
00717 Spectrum fr = bsdf->f(wo, wi);
00718 if (fr.IsBlack()) continue;
00719 RayDifferential bounceRay(p, wi, ray, isect.rayEpsilon);
00720 Intersection gatherIsect;
00721 PBRT_PHOTON_MAP_STARTED_GATHER_RAY(&bounceRay);
00722 if (scene->Intersect(bounceRay, &gatherIsect)) {
00723
00724 Spectrum Lindir = 0.f;
00725 Normal nGather = gatherIsect.dg.nn;
00726 nGather = Faceforward(nGather, -bounceRay.d);
00727 RadiancePhotonProcess proc(nGather);
00728 float md2 = INFINITY;
00729 radianceMap->Lookup(gatherIsect.dg.p, proc, md2);
00730 if (proc.photon != NULL)
00731 Lindir = proc.photon->Lo;
00732 Lindir *= renderer->Transmittance(scene, bounceRay, NULL, rng, arena);
00733
00734
00735 float photonPdf = 0.f;
00736 float conePdf = UniformConePdf(cosGatherAngle);
00737 for (uint32_t j = 0; j < nIndirSamplePhotons; ++j)
00738 if (Dot(photonDirs[j], wi) > .999f * cosGatherAngle)
00739 photonPdf += conePdf;
00740 photonPdf /= nIndirSamplePhotons;
00741
00742
00743 float bsdfPdf = bsdf->Pdf(wo, wi);
00744 float wt = PowerHeuristic(gatherSamples, photonPdf, gatherSamples, bsdfPdf);
00745 Li += fr * Lindir * AbsDot(wi, n) * wt / photonPdf;
00746 }
00747 PBRT_PHOTON_MAP_FINISHED_GATHER_RAY(&bounceRay);
00748 }
00749 L += Li / gatherSamples;
00750 }
00751 #else
00752
00753 Normal nn = Faceforward(n, -ray.d);
00754 RadiancePhotonProcess proc(nn);
00755 float md2 = INFINITY;
00756 radianceMap->Lookup(p, proc, md2);
00757 if (proc.photon)
00758 L += proc.photon->Lo;
00759 #endif
00760 }
00761 else
00762 L += LPhoton(indirectMap, nIndirectPaths, nLookup, lookupBuf,
00763 bsdf, rng, isect, wo, maxDistSquared);
00764 if (ray.depth+1 < maxSpecularDepth) {
00765 Vector wi;
00766
00767 L += SpecularReflect(ray, bsdf, rng, isect, renderer, scene, sample,
00768 arena);
00769 L += SpecularTransmit(ray, bsdf, rng, isect, renderer, scene, sample,
00770 arena);
00771 }
00772 return L;
00773 }
00774
00775
00776 PhotonIntegrator *CreatePhotonMapSurfaceIntegrator(const ParamSet ¶ms) {
00777 int nCaustic = params.FindOneInt("causticphotons", 20000);
00778 int nIndirect = params.FindOneInt("indirectphotons", 100000);
00779 int nUsed = params.FindOneInt("nused", 50);
00780 if (PbrtOptions.quickRender) nCaustic = nCaustic / 10;
00781 if (PbrtOptions.quickRender) nIndirect = nIndirect / 10;
00782 if (PbrtOptions.quickRender) nUsed = max(1, nUsed / 10);
00783 int maxSpecularDepth = params.FindOneInt("maxspeculardepth", 5);
00784 int maxPhotonDepth = params.FindOneInt("maxphotondepth", 5);
00785 bool finalGather = params.FindOneBool("finalgather", true);
00786 int gatherSamples = params.FindOneInt("finalgathersamples", 32);
00787 if (PbrtOptions.quickRender) gatherSamples = max(1, gatherSamples / 4);
00788 float maxDist = params.FindOneFloat("maxdist", .1f);
00789 float gatherAngle = params.FindOneFloat("gatherangle", 10.f);
00790 return new PhotonIntegrator(nCaustic, nIndirect,
00791 nUsed, maxSpecularDepth, maxPhotonDepth, maxDist, finalGather, gatherSamples,
00792 gatherAngle);
00793 }
00794
00795