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 Photon;
00033 struct ClosePhoton;
00034 struct PhotonProcess;
00035 class PhotonIntegrator : public SurfaceIntegrator {
00036 public:
00037
00038 PhotonIntegrator(int ncaus, int ndir, int nindir, int nLookup, int mdepth,
00039 float maxdist, bool finalGather, int gatherSamples,
00040 bool directWithPhotons);
00041 ~PhotonIntegrator();
00042 Spectrum Li(const Scene *scene, const RayDifferential &ray,
00043 const Sample *sample, float *alpha) const;
00044 void RequestSamples(Sample *sample, const Scene *scene);
00045 void Preprocess(const Scene *);
00046 private:
00047
00048 static inline bool unsuccessful(int needed, int found, int shot) {
00049 return (found < needed &&
00050 (found == 0 || found < shot / 1024));
00051 }
00052 static Spectrum LPhoton(KdTree<Photon, PhotonProcess> *map,
00053 int nPaths, int nLookup, BSDF *bsdf, const Intersection &isect,
00054 const Vector &w, float maxDistSquared);
00055
00056 u_int nCausticPhotons, nIndirectPhotons, nDirectPhotons;
00057 u_int nLookup;
00058 mutable int specularDepth;
00059 int maxSpecularDepth;
00060 float maxDistSquared;
00061 bool directWithPhotons, finalGather;
00062 int gatherSamples;
00063
00064 int *lightSampleOffset, lightNumOffset;
00065 int *bsdfSampleOffset, *bsdfComponentOffset;
00066 int gatherSampleOffset, gatherComponentOffset;
00067 int nCausticPaths, nDirectPaths, nIndirectPaths;
00068 mutable KdTree<Photon, PhotonProcess> *causticMap;
00069 mutable KdTree<Photon, PhotonProcess> *directMap;
00070 mutable KdTree<Photon, PhotonProcess> *indirectMap;
00071 };
00072 struct Photon {
00073
00074 Photon(const Point &pp, const Spectrum &wt, const Vector &w)
00075 : p(pp), alpha(wt), wi(w) {
00076 }
00077 Photon() { }
00078 Point p;
00079 Spectrum alpha;
00080 Vector wi;
00081 };
00082 struct PhotonProcess {
00083
00084 PhotonProcess(u_int mp, const Point &p);
00085 void operator()(const Photon &photon, float dist2, float &maxDistSquared) const;
00086 const Point &p;
00087 ClosePhoton *photons;
00088 u_int nLookup;
00089 mutable u_int foundPhotons;
00090 };
00091 struct ClosePhoton {
00092 ClosePhoton(const Photon *p = NULL,
00093 float md2 = INFINITY) {
00094 photon = p;
00095 distanceSquared = md2;
00096 }
00097 bool operator<(const ClosePhoton &p2) const {
00098 return distanceSquared == p2.distanceSquared ? (photon < p2.photon) :
00099 distanceSquared < p2.distanceSquared;
00100 }
00101 const Photon *photon;
00102 float distanceSquared;
00103 };
00104
00105 PhotonIntegrator::PhotonIntegrator(int ncaus, int ndir, int nind,
00106 int nl, int mdepth, float mdist, bool fg,
00107 int gs, bool dp) {
00108 nCausticPhotons = ncaus;
00109 nIndirectPhotons = nind;
00110 nDirectPhotons = ndir;
00111 nLookup = nl;
00112 maxDistSquared = mdist * mdist;
00113 maxSpecularDepth = mdepth;
00114 causticMap = directMap = indirectMap = NULL;
00115 specularDepth = 0;
00116 finalGather = fg;
00117 gatherSamples = gs;
00118 directWithPhotons = dp;
00119 }
00120 PhotonIntegrator::~PhotonIntegrator() {
00121 delete causticMap;
00122 delete directMap;
00123 delete indirectMap;
00124 }
00125 void PhotonIntegrator::RequestSamples(Sample *sample,
00126 const Scene *scene) {
00127
00128 u_int nLights = scene->lights.size();
00129 lightSampleOffset = new int[nLights];
00130 bsdfSampleOffset = new int[nLights];
00131 bsdfComponentOffset = new int[nLights];
00132 for (u_int i = 0; i < nLights; ++i) {
00133 const Light *light = scene->lights[i];
00134 int lightSamples =
00135 scene->sampler->RoundSize(light->nSamples);
00136 lightSampleOffset[i] = sample->Add2D(lightSamples);
00137 bsdfSampleOffset[i] = sample->Add2D(lightSamples);
00138 bsdfComponentOffset[i] = sample->Add1D(lightSamples);
00139 }
00140 lightNumOffset = -1;
00141 if (finalGather) {
00142 gatherSamples = scene->sampler->RoundSize(gatherSamples);
00143 gatherSampleOffset = sample->Add2D(gatherSamples);
00144 gatherComponentOffset = sample->Add1D(gatherSamples);
00145 }
00146 }
00147 void PhotonIntegrator::Preprocess(const Scene *scene) {
00148 if (scene->lights.size() == 0) return;
00149 ProgressReporter progress(nCausticPhotons+nDirectPhotons+
00150 nIndirectPhotons, "Shooting photons");
00151 vector<Photon> causticPhotons;
00152 vector<Photon> directPhotons;
00153 vector<Photon> indirectPhotons;
00154 causticPhotons.reserve(nCausticPhotons);
00155 directPhotons.reserve(nDirectPhotons);
00156 indirectPhotons.reserve(nIndirectPhotons);
00157
00158 static StatsCounter nshot("Photon Map",
00159 "Number of photons shot from lights");
00160 bool causticDone = (nCausticPhotons == 0);
00161 bool directDone = (nDirectPhotons == 0);
00162 bool indirectDone = (nIndirectPhotons == 0);
00163 while (!causticDone || !directDone || !indirectDone) {
00164 ++nshot;
00165
00166 if (nshot > 500000 &&
00167 (unsuccessful(nCausticPhotons,
00168 causticPhotons.size(),
00169 nshot) ||
00170 unsuccessful(nDirectPhotons,
00171 directPhotons.size(),
00172 nshot) ||
00173 unsuccessful(nIndirectPhotons,
00174 indirectPhotons.size(),
00175 nshot))) {
00176 Error("Unable to store enough photons. Giving up.\n");
00177 return;
00178 }
00179
00180
00181 float u[4];
00182 u[0] = (float)RadicalInverse((int)nshot+1, 2);
00183 u[1] = (float)RadicalInverse((int)nshot+1, 3);
00184 u[2] = (float)RadicalInverse((int)nshot+1, 5);
00185 u[3] = (float)RadicalInverse((int)nshot+1, 7);
00186
00187 int nLights = int(scene->lights.size());
00188 int lightNum =
00189 min(Floor2Int(nLights * (float)RadicalInverse((int)nshot+1, 11)),
00190 nLights-1);
00191 Light *light = scene->lights[lightNum];
00192 float lightPdf = 1.f / nLights;
00193
00194 RayDifferential photonRay;
00195 float pdf;
00196 Spectrum alpha =
00197 light->Sample_L(scene, u[0], u[1], u[2], u[3],
00198 &photonRay, &pdf);
00199 if (pdf == 0.f || alpha.Black()) continue;
00200 alpha /= pdf * lightPdf;
00201 if (!alpha.Black()) {
00202
00203 bool specularPath = false;
00204 Intersection photonIsect;
00205 int nIntersections = 0;
00206 while (scene->Intersect(photonRay, &photonIsect)) {
00207 ++nIntersections;
00208
00209 alpha *= scene->Transmittance(photonRay);
00210 Vector wo = -photonRay.d;
00211 BSDF *photonBSDF = photonIsect.GetBSDF(photonRay);
00212 BxDFType specularType = BxDFType(BSDF_REFLECTION |
00213 BSDF_TRANSMISSION | BSDF_SPECULAR);
00214 bool hasNonSpecular = (photonBSDF->NumComponents() >
00215 photonBSDF->NumComponents(specularType));
00216 if (hasNonSpecular) {
00217
00218 Photon photon(photonIsect.dg.p, alpha, wo);
00219 if (nIntersections == 1) {
00220
00221 if (!directDone) {
00222 directPhotons.push_back(photon);
00223 if (directPhotons.size() == nDirectPhotons) {
00224 directDone = true;
00225 nDirectPaths = (int)nshot;
00226 directMap =
00227 new KdTree<Photon,
00228 PhotonProcess>(directPhotons);
00229 }
00230 progress.Update();
00231 }
00232 }
00233 else if (specularPath) {
00234
00235 if (!causticDone) {
00236 causticPhotons.push_back(photon);
00237 if (causticPhotons.size() == nCausticPhotons) {
00238 causticDone = true;
00239 nCausticPaths = (int)nshot;
00240 causticMap =
00241 new KdTree<Photon,
00242 PhotonProcess>(causticPhotons);
00243 }
00244 progress.Update();
00245 }
00246 }
00247 else {
00248
00249 if (!indirectDone) {
00250 indirectPhotons.push_back(photon);
00251 if (indirectPhotons.size() == nIndirectPhotons) {
00252 indirectDone = true;
00253 nIndirectPaths = (int)nshot;
00254 indirectMap =
00255 new KdTree<Photon,
00256 PhotonProcess>(indirectPhotons);
00257 }
00258 progress.Update();
00259 }
00260 }
00261 }
00262
00263 Vector wi;
00264 float pdf;
00265 BxDFType flags;
00266
00267 float u1, u2, u3;
00268 if (nIntersections == 1) {
00269 u1 = (float)RadicalInverse((int)nshot+1, 13);
00270 u2 = (float)RadicalInverse((int)nshot+1, 17);
00271 u3 = (float)RadicalInverse((int)nshot+1, 19);
00272 }
00273 else {
00274 u1 = RandomFloat();
00275 u2 = RandomFloat();
00276 u3 = RandomFloat();
00277 }
00278 Spectrum fr = photonBSDF->Sample_f(wo, &wi, u1, u2, u3,
00279 &pdf, BSDF_ALL, &flags);
00280 if (fr.Black() || pdf == 0.f)
00281 break;
00282 specularPath = (nIntersections == 1 || specularPath) &&
00283 ((flags & BSDF_SPECULAR) != 0);
00284 alpha *= fr * AbsDot(wi, photonBSDF->dgShading.nn) / pdf;
00285 photonRay = RayDifferential(photonIsect.dg.p, wi);
00286
00287 if (nIntersections > 3) {
00288 float continueProbability = .5f;
00289 if (RandomFloat() > continueProbability)
00290 break;
00291 alpha /= continueProbability;
00292 }
00293 }
00294 }
00295 BSDF::FreeAll();
00296 }
00297 progress.Done();
00298 }
00299 Spectrum PhotonIntegrator::Li(const Scene *scene,
00300 const RayDifferential &ray, const Sample *sample,
00301 float *alpha) const {
00302
00303 Spectrum L(0.);
00304 Intersection isect;
00305 if (scene->Intersect(ray, &isect)) {
00306 if (alpha) *alpha = 1.;
00307 Vector wo = -ray.d;
00308
00309 L += isect.Le(wo);
00310
00311 BSDF *bsdf = isect.GetBSDF(ray);
00312 const Point &p = bsdf->dgShading.p;
00313 const Normal &n = bsdf->dgShading.nn;
00314
00315 if (directWithPhotons)
00316 L += LPhoton(directMap, nDirectPaths, nLookup,
00317 bsdf, isect, wo, maxDistSquared);
00318 else
00319 L += UniformSampleAllLights(scene, p, n,
00320 wo, bsdf, sample,
00321 lightSampleOffset, bsdfSampleOffset,
00322 bsdfComponentOffset);
00323
00324
00325 L += LPhoton(causticMap, nCausticPaths, nLookup, bsdf,
00326 isect, wo, maxDistSquared);
00327 if (finalGather) {
00328
00329 Spectrum Li(0.);
00330 for (int i = 0; i < gatherSamples; ++i) {
00331
00332 Vector wi;
00333 float u1 = sample->twoD[gatherSampleOffset][2*i];
00334 float u2 = sample->twoD[gatherSampleOffset][2*i+1];
00335 float u3 = sample->oneD[gatherComponentOffset][i];
00336 float pdf;
00337 Spectrum fr = bsdf->Sample_f(wo, &wi, u1, u2, u3,
00338 &pdf, BxDFType(BSDF_ALL & (~BSDF_SPECULAR)));
00339 if (fr.Black() || pdf == 0.f) continue;
00340 RayDifferential bounceRay(p, wi);
00341 static StatsCounter gatherRays("Photon Map",
00342 "Final gather rays traced");
00343 ++gatherRays;
00344 Intersection gatherIsect;
00345 if (scene->Intersect(bounceRay, &gatherIsect)) {
00346
00347 BSDF *gatherBSDF = gatherIsect.GetBSDF(bounceRay);
00348 Vector bounceWo = -bounceRay.d;
00349 Spectrum Lindir =
00350 LPhoton(directMap, nDirectPaths, nLookup,
00351 gatherBSDF, gatherIsect, bounceWo, maxDistSquared) +
00352 LPhoton(indirectMap, nIndirectPaths, nLookup,
00353 gatherBSDF, gatherIsect, bounceWo, maxDistSquared) +
00354 LPhoton(causticMap, nCausticPaths, nLookup,
00355 gatherBSDF, gatherIsect, bounceWo, maxDistSquared);
00356 Lindir *= scene->Transmittance(bounceRay);
00357 Li += fr * Lindir * AbsDot(wi, n) / pdf;
00358 }
00359 }
00360 L += Li / float(gatherSamples);
00361 }
00362 else
00363 L += LPhoton(indirectMap, nIndirectPaths, nLookup,
00364 bsdf, isect, wo, maxDistSquared);
00365 if (specularDepth++ < maxSpecularDepth) {
00366 Vector wi;
00367
00368 Spectrum f = bsdf->Sample_f(wo, &wi,
00369 BxDFType(BSDF_REFLECTION | BSDF_SPECULAR));
00370 if (!f.Black()) {
00371
00372 RayDifferential rd(p, wi);
00373 rd.hasDifferentials = true;
00374 rd.rx.o = p + isect.dg.dpdx;
00375 rd.ry.o = p + isect.dg.dpdy;
00376
00377 Normal dndx = bsdf->dgShading.dndu * bsdf->dgShading.dudx +
00378 bsdf->dgShading.dndv * bsdf->dgShading.dvdx;
00379 Normal dndy = bsdf->dgShading.dndu * bsdf->dgShading.dudy +
00380 bsdf->dgShading.dndv * bsdf->dgShading.dvdy;
00381 Vector dwodx = -ray.rx.d - wo, dwody = -ray.ry.d - wo;
00382 float dDNdx = Dot(dwodx, n) + Dot(wo, dndx);
00383 float dDNdy = Dot(dwody, n) + Dot(wo, dndy);
00384 rd.rx.d = wi -
00385 dwodx + 2 * Vector(Dot(wo, n) * dndx +
00386 dDNdx * n);
00387 rd.ry.d = wi -
00388 dwody + 2 * Vector(Dot(wo, n) * dndy +
00389 dDNdy * n);
00390 L += scene->Li(rd, sample) * f * AbsDot(wi, n);
00391 }
00392 f = bsdf->Sample_f(wo, &wi,
00393 BxDFType(BSDF_TRANSMISSION | BSDF_SPECULAR));
00394 if (!f.Black()) {
00395
00396 RayDifferential rd(p, wi);
00397 rd.hasDifferentials = true;
00398 rd.rx.o = p + isect.dg.dpdx;
00399 rd.ry.o = p + isect.dg.dpdy;
00400
00401 float eta = bsdf->eta;
00402 Vector w = -wo;
00403 if (Dot(wo, n) < 0) eta = 1.f / eta;
00404
00405 Normal dndx = bsdf->dgShading.dndu * bsdf->dgShading.dudx + bsdf->dgShading.dndv * bsdf->dgShading.dvdx;
00406 Normal dndy = bsdf->dgShading.dndu * bsdf->dgShading.dudy + bsdf->dgShading.dndv * bsdf->dgShading.dvdy;
00407
00408 Vector dwodx = -ray.rx.d - wo, dwody = -ray.ry.d - wo;
00409 float dDNdx = Dot(dwodx, n) + Dot(wo, dndx);
00410 float dDNdy = Dot(dwody, n) + Dot(wo, dndy);
00411
00412 float mu = eta * Dot(w, n) - Dot(wi, n);
00413 float dmudx = (eta - (eta*eta*Dot(w,n))/Dot(wi, n)) * dDNdx;
00414 float dmudy = (eta - (eta*eta*Dot(w,n))/Dot(wi, n)) * dDNdy;
00415
00416 rd.rx.d = wi + eta * dwodx - Vector(mu * dndx + dmudx * n);
00417 rd.ry.d = wi + eta * dwody - Vector(mu * dndy + dmudy * n);
00418 L += scene->Li(rd, sample) * f * AbsDot(wi, n);
00419 }
00420 }
00421 --specularDepth;
00422 }
00423 else {
00424
00425 if (alpha) *alpha = 0.;
00426 for (u_int i = 0; i < scene->lights.size(); ++i)
00427 L += scene->lights[i]->Le(ray);
00428 if (alpha && !L.Black()) *alpha = 1.;
00429 return L;
00430 }
00431 return L;
00432 }
00433 Spectrum PhotonIntegrator::LPhoton(
00434 KdTree<Photon, PhotonProcess> *map,
00435 int nPaths, int nLookup, BSDF *bsdf,
00436 const Intersection &isect, const Vector &wo,
00437 float maxDistSquared) {
00438 Spectrum L(0.);
00439 if (!map) return L;
00440 BxDFType nonSpecular = BxDFType(BSDF_REFLECTION |
00441 BSDF_TRANSMISSION | BSDF_DIFFUSE | BSDF_GLOSSY);
00442 if (bsdf->NumComponents(nonSpecular) == 0)
00443 return L;
00444 static StatsCounter lookups("Photon Map", "Total lookups");
00445
00446 PhotonProcess proc(nLookup, isect.dg.p);
00447 proc.photons =
00448 (ClosePhoton *)alloca(nLookup * sizeof(ClosePhoton));
00449
00450 ++lookups;
00451 map->Lookup(isect.dg.p, proc, maxDistSquared);
00452
00453 static StatsRatio foundRate("Photon Map", "Photons found per lookup");
00454 foundRate.Add(proc.foundPhotons, 1);
00455 float scale = 1.f / (float(nPaths) * maxDistSquared * M_PI);
00456
00457 ClosePhoton *photons = proc.photons;
00458 int nFound = proc.foundPhotons;
00459 Normal Nf = Dot(wo, bsdf->dgShading.nn) < 0 ? -bsdf->dgShading.nn :
00460 bsdf->dgShading.nn;
00461 if (bsdf->NumComponents(BxDFType(BSDF_REFLECTION |
00462 BSDF_TRANSMISSION | BSDF_GLOSSY)) > 0) {
00463
00464 for (int i = 0; i < nFound; ++i) {
00465 BxDFType flag = Dot(Nf, photons[i].photon->wi) > 0.f ?
00466 BSDF_ALL_REFLECTION : BSDF_ALL_TRANSMISSION;
00467 L += bsdf->f(wo, photons[i].photon->wi, flag) *
00468 (scale * photons[i].photon->alpha);
00469 }
00470 }
00471 else {
00472
00473 Spectrum Lr(0.), Lt(0.);
00474 for (int i = 0; i < nFound; ++i)
00475 if (Dot(Nf, photons[i].photon->wi) > 0.f)
00476 Lr += photons[i].photon->alpha;
00477 else
00478 Lt += photons[i].photon->alpha;
00479 L += (scale * INV_PI) * (Lr * bsdf->rho(wo, BSDF_ALL_REFLECTION) +
00480 Lt * bsdf->rho(wo, BSDF_ALL_TRANSMISSION));
00481 }
00482 return L;
00483 }
00484 PhotonProcess::PhotonProcess(u_int mp, const Point &P)
00485 : p(P) {
00486 photons = 0;
00487 nLookup = mp;
00488 foundPhotons = 0;
00489 }
00490 void PhotonProcess::operator()(const Photon &photon,
00491 float distSquared, float &maxDistSquared) const {
00492 static StatsPercentage discarded("Photon Map", "Discarded photons");
00493 discarded.Add(0, 1);
00494 if (foundPhotons < nLookup) {
00495
00496 photons[foundPhotons++] = ClosePhoton(&photon, distSquared);
00497 if (foundPhotons == nLookup) {
00498 std::make_heap(&photons[0], &photons[nLookup]);
00499 maxDistSquared = photons[0].distanceSquared;
00500 }
00501 }
00502 else {
00503
00504 discarded.Add(1, 0);
00505 std::pop_heap(&photons[0], &photons[nLookup]);
00506 photons[nLookup-1] = ClosePhoton(&photon, distSquared);
00507 std::push_heap(&photons[0], &photons[nLookup]);
00508 maxDistSquared = photons[0].distanceSquared;
00509 }
00510 }
00511 extern "C" DLLEXPORT SurfaceIntegrator *CreateSurfaceIntegrator(const ParamSet ¶ms) {
00512 int nCaustic = params.FindOneInt("causticphotons", 20000);
00513 int nDirect = params.FindOneInt("directphotons", 100000);
00514 int nIndirect = params.FindOneInt("indirectphotons", 100000);
00515 int nUsed = params.FindOneInt("nused", 50);
00516 int maxDepth = params.FindOneInt("maxdepth", 5);
00517 bool finalGather = params.FindOneBool("finalgather", true);
00518 bool directPhotons = params.FindOneBool("directwithphotons", false);
00519 int gatherSamples = params.FindOneInt("finalgathersamples", 32);
00520 float maxDist = params.FindOneFloat("maxdist", .1f);
00521 return new PhotonIntegrator(nCaustic, nDirect, nIndirect,
00522 nUsed, maxDepth, maxDist, finalGather, gatherSamples,
00523 directPhotons);
00524 }