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