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 "exphotonmap.h"
00025 #include "bxdf.h"
00026 #include "light.h"
00027 #include "paramset.h"
00028 #include "spectrumwavelengths.h"
00029
00030 #include <boost/thread/xtime.hpp>
00031
00032 using namespace lux;
00033
00034
00035 extern boost::thread_specific_ptr<SpectrumWavelengths> thread_wavelengths;
00036
00037
00038 ExPhotonIntegrator* ExPhotonIntegrator::clone() const {
00039 return new ExPhotonIntegrator(*this);
00040 }
00041
00042
00043 EPhotonProcess::EPhotonProcess(u_int mp, const Point &P)
00044 : p(P) {
00045 photons = 0;
00046 nLookup = mp;
00047 foundPhotons = 0;
00048 }
00049
00050 SWCSpectrum ExPhotonIntegrator::estimateE(
00051 KdTree<EPhoton, EPhotonProcess> *map, int count,
00052 const Point &p, const Normal &n) const {
00053 if (!map) return 0.f;
00054
00055 EPhotonProcess proc(nLookup, p);
00056 proc.photons = (EClosePhoton *) alloca(nLookup *
00057 sizeof (EClosePhoton));
00058 float md2 = maxDistSquared;
00059 map->Lookup(p, proc, md2);
00060
00061 EClosePhoton *photons = proc.photons;
00062 SWCSpectrum E(0.);
00063 for (u_int i = 0; i < proc.foundPhotons; ++i) {
00064 if (Dot(n, photons[i].photon->wi) > 0.) {
00065 SWCSpectrum alpha = FromXYZ(photons[i].photon->alpha.c[0],
00066 photons[i].photon->alpha.c[1],
00067 photons[i].photon->alpha.c[2]);
00068 E += alpha;
00069 }
00070 }
00071
00072 return E / (float(count) * md2 * M_PI);
00073 }
00074
00075 void EPhotonProcess::operator()(const EPhoton &photon,
00076 float distSquared, float &maxDistSquared) const {
00077
00078
00079
00080 if (foundPhotons < nLookup) {
00081
00082 photons[foundPhotons++] = EClosePhoton(&photon, distSquared);
00083 if (foundPhotons == nLookup) {
00084 std::make_heap(&photons[0], &photons[nLookup]);
00085 maxDistSquared = photons[0].distanceSquared;
00086 }
00087 } else {
00088
00089
00090 std::pop_heap(&photons[0], &photons[nLookup]);
00091 photons[nLookup - 1] = EClosePhoton(&photon, distSquared);
00092 std::push_heap(&photons[0], &photons[nLookup]);
00093 maxDistSquared = photons[0].distanceSquared;
00094 }
00095 }
00096
00097 SWCSpectrum ExPhotonIntegrator::LPhoton(
00098 KdTree<EPhoton, EPhotonProcess> *map,
00099 int nPaths, int nLookup, BSDF *bsdf,
00100 const Intersection &isect, const Vector &wo,
00101 float maxDistSquared) {
00102 SWCSpectrum L(0.);
00103 if ((nPaths <= 0) || (!map)) return L;
00104
00105 BxDFType nonSpecular = BxDFType(BSDF_REFLECTION |
00106 BSDF_TRANSMISSION | BSDF_DIFFUSE | BSDF_GLOSSY);
00107 if (bsdf->NumComponents(nonSpecular) == 0)
00108 return L;
00109
00110
00111 EPhotonProcess proc(nLookup, isect.dg.p);
00112 proc.photons =
00113 (EClosePhoton *) alloca(nLookup * sizeof (EClosePhoton));
00114
00115
00116 map->Lookup(isect.dg.p, proc, maxDistSquared);
00117
00118
00119
00120
00121 EClosePhoton *photons = proc.photons;
00122 int nFound = proc.foundPhotons;
00123 Normal Nf = Dot(wo, bsdf->dgShading.nn) < 0 ? -bsdf->dgShading.nn :
00124 bsdf->dgShading.nn;
00125
00126 if (bsdf->NumComponents(BxDFType(BSDF_REFLECTION |
00127 BSDF_TRANSMISSION | BSDF_GLOSSY)) > 0) {
00128
00129 for (int i = 0; i < nFound; ++i) {
00130 const EPhoton *p = photons[i].photon;
00131 BxDFType flag = Dot(Nf, p->wi) > 0.f ?
00132 BSDF_ALL_REFLECTION : BSDF_ALL_TRANSMISSION;
00133 float k = Ekernel(p, isect.dg.p, maxDistSquared);
00134
00135 SWCSpectrum alpha = FromXYZ(p->alpha.c[0], p->alpha.c[1], p->alpha.c[2]);
00136 L += (k / nPaths) * bsdf->f(wo, p->wi, flag) * alpha;
00137 }
00138 } else {
00139
00140 SWCSpectrum Lr(0.), Lt(0.);
00141
00142 for (int i = 0; i < nFound; ++i) {
00143 SWCSpectrum alpha = FromXYZ(photons[i].photon->alpha.c[0],
00144 photons[i].photon->alpha.c[1],
00145 photons[i].photon->alpha.c[2]);
00146
00147 if (Dot(Nf, photons[i].photon->wi) > 0.f) {
00148 float k = Ekernel(photons[i].photon, isect.dg.p,
00149 maxDistSquared);
00150 Lr += (k / nPaths) * alpha;
00151 } else {
00152 float k = Ekernel(photons[i].photon, isect.dg.p,
00153 maxDistSquared);
00154 Lt += (k / nPaths) * alpha;
00155 }
00156 }
00157
00158 L += Lr * bsdf->rho(wo, BSDF_ALL_REFLECTION) * INV_PI +
00159 Lt * bsdf->rho(wo, BSDF_ALL_TRANSMISSION) * INV_PI;
00160 }
00161
00162 return L;
00163 }
00164
00165 ExPhotonIntegrator::ExPhotonIntegrator(
00166 LightStrategy st,
00167 int ncaus, int nind, int maxDirPhotons,
00168 int nl, int mdepth, float mdist, bool fg,
00169 int gs, float ga,
00170 RRStrategy grrStrategy, float grrContinueProbability,
00171 bool dbgEnableDirect, bool dbgEnableCaustic,
00172 bool dbgEnableIndirect, bool dbgEnableSpecular) {
00173 lightStrategy = st;
00174
00175 nCausticPhotons = ncaus;
00176 nIndirectPhotons = nind;
00177 maxDirectPhotons = maxDirPhotons;
00178
00179 nLookup = nl;
00180 maxDistSquared = mdist * mdist;
00181 maxSpecularDepth = mdepth;
00182 causticMap = indirectMap = NULL;
00183 radianceMap = NULL;
00184 finalGather = fg;
00185 gatherSamples = gs;
00186 cosGatherAngle = cos(Radians(ga));
00187
00188 gatherRRStrategy = grrStrategy;
00189 gatherRRContinueProbability = grrContinueProbability;
00190
00191 debugEnableDirect = dbgEnableDirect;
00192 debugEnableCaustic = dbgEnableCaustic;
00193 debugEnableIndirect = dbgEnableIndirect;
00194 debugEnableSpecular = dbgEnableSpecular;
00195 }
00196
00197 ExPhotonIntegrator::~ExPhotonIntegrator() {
00198 delete causticMap;
00199 delete indirectMap;
00200 delete radianceMap;
00201 }
00202
00203 void ExPhotonIntegrator::RequestSamples(Sample *sample,
00204 const Scene *scene) {
00205 if (lightStrategy == SAMPLE_AUTOMATIC) {
00206 if (scene->lights.size() > 5)
00207 lightStrategy = SAMPLE_ONE_UNIFORM;
00208 else
00209 lightStrategy = SAMPLE_ALL_UNIFORM;
00210 }
00211
00212
00213 vector<u_int> structure;
00214
00215 structure.push_back(2);
00216 structure.push_back(2);
00217 structure.push_back(1);
00218
00219 structure.push_back(2);
00220 structure.push_back(1);
00221 structure.push_back(2);
00222 structure.push_back(1);
00223
00224 if (lightStrategy == SAMPLE_ONE_UNIFORM)
00225 structure.push_back(1);
00226
00227 sampleOffset = sample->AddxD(structure, maxSpecularDepth + 1);
00228
00229 if (finalGather) {
00230
00231
00232 gatherSamples = gatherSamples / 2;
00233
00234
00235 structure.clear();
00236 structure.push_back(2);
00237 structure.push_back(1);
00238 if (gatherRRStrategy != RR_NONE)
00239 structure.push_back(1);
00240 sampleFinalGather1Offset = sample->AddxD(structure, gatherSamples);
00241
00242 structure.clear();
00243 structure.push_back(2);
00244 structure.push_back(1);
00245 if (gatherRRStrategy != RR_NONE)
00246 structure.push_back(1);
00247 sampleFinalGather2Offset = sample->AddxD(structure, gatherSamples);
00248 }
00249 }
00250
00251 void ExPhotonIntegrator::Preprocess(const Scene *scene) {
00252 if (scene->lights.size() == 0) return;
00253
00254
00255 std::stringstream ss;
00256 ss << "Shooting photons: " << (nCausticPhotons + nIndirectPhotons);
00257 luxError(LUX_NOERROR, LUX_INFO, ss.str().c_str());
00258
00259 vector<EPhoton> causticPhotons;
00260 vector<EPhoton> indirectPhotons;
00261 vector<EPhoton> directPhotons;
00262 vector<ERadiancePhoton> radiancePhotons;
00263 causticPhotons.reserve(nCausticPhotons);
00264 indirectPhotons.reserve(nIndirectPhotons);
00265
00266
00267 SpectrumWavelengths *thr_wl = thread_wavelengths.get();
00268
00269 bool causticDone = (nCausticPhotons == 0);
00270 bool indirectDone = (nIndirectPhotons == 0);
00271 bool directDone = false;
00272
00273
00274 int nLights = int(scene->lights.size());
00275 float *lightPower = (float *)alloca(nLights * sizeof(float));
00276 float *lightCDF = (float *)alloca((nLights+1) * sizeof(float));
00277
00278
00279 const int spectrumSamples = 128;
00280 for (int i = 0; i < nLights; ++i)
00281 lightPower[i] +=0.0f;
00282 for (int j = 0; j < spectrumSamples; j++) {
00283 thr_wl->Sample((float)RadicalInverse(j, 2),
00284 (float)RadicalInverse(j, 3));
00285
00286 for (int i = 0; i < nLights; ++i)
00287 lightPower[i] += scene->lights[i]->Power(scene).y();
00288 }
00289 for (int i = 0; i < nLights; ++i)
00290 lightPower[i] *= 1.0f / spectrumSamples;
00291
00292 float totalPower;
00293 ComputeStep1dCDF(lightPower, nLights, &totalPower, lightCDF);
00294
00295
00296 vector<SWCSpectrum> rpReflectances, rpTransmittances;
00297
00298 boost::xtime lastUpdateTime;
00299 boost::xtime_get(&lastUpdateTime, boost::TIME_UTC);
00300 int nDirectPaths = 0;
00301 int nshot = 0;
00302 while (!causticDone || !indirectDone) {
00303
00304 boost::xtime currentTime;
00305 boost::xtime_get(¤tTime, boost::TIME_UTC);
00306 if (currentTime.sec - lastUpdateTime.sec > 5) {
00307 ss.str("");
00308 ss << "Direct photon count: " << directPhotons.size();
00309 if (maxDirectPhotons > 0)
00310 ss << " (" << (100 * directPhotons.size() / maxDirectPhotons) << "% limit)";
00311 else
00312 ss << " (100% limit)";
00313 ss << " Caustic photonmap size: " << causticPhotons.size();
00314 if (nCausticPhotons > 0)
00315 ss << " (" << (100 * causticPhotons.size() / nCausticPhotons) << "%)";
00316 else
00317 ss << " (100%)";
00318 ss << " Indirect photonmap size: " << indirectPhotons.size();
00319 if (nIndirectPhotons > 0)
00320 ss << " (" << (100 * indirectPhotons.size() / nIndirectPhotons) << "%)";
00321 else
00322 ss << " (100%)";
00323 luxError(LUX_NOERROR, LUX_INFO, ss.str().c_str());
00324
00325 lastUpdateTime = currentTime;
00326 }
00327
00328 ++nshot;
00329
00330
00331 if (nshot > 500000) {
00332 if (indirectDone && unsuccessful(nCausticPhotons, causticPhotons.size(), nshot)) {
00333
00334
00335 luxError(LUX_CONSISTENCY, LUX_WARNING, "Unable to store enough photons in the caustic photonmap. Giving up and disabling the map.");
00336
00337 causticPhotons.clear();
00338 causticDone = true;
00339 nCausticPhotons = 0;
00340 }
00341
00342 if(unsuccessful(nIndirectPhotons, indirectPhotons.size(), nshot)) {
00343 luxError(LUX_CONSISTENCY, LUX_ERROR, "Unable to store enough photons in the indirect photonmap. Giving up. I will unable to reder the image.");
00344 return;
00345 }
00346 }
00347
00348
00349
00350 float u[4];
00351 u[0] = RadicalInverse(nshot, 2);
00352 u[1] = RadicalInverse(nshot, 3);
00353 u[2] = RadicalInverse(nshot, 5);
00354 u[3] = RadicalInverse(nshot, 7);
00355
00356
00357 thr_wl->Sample(RadicalInverse(nshot, 23),
00358 RadicalInverse(nshot, 29));
00359
00360
00361 float lightPdf;
00362 float uln = RadicalInverse(nshot, 11);
00363 int lightNum = Floor2Int(SampleStep1d(lightPower, lightCDF,
00364 totalPower, nLights, uln, &lightPdf) * nLights);
00365 lightNum = min(lightNum, nLights - 1);
00366 const Light *light = scene->lights[lightNum];
00367
00368
00369 RayDifferential photonRay;
00370 float pdf;
00371 SWCSpectrum alpha = light->Sample_L(scene, u[0], u[1], u[2], u[3],
00372 &photonRay, &pdf);
00373 if (pdf == 0.f || alpha.Black()) continue;
00374 alpha /= pdf * lightPdf;
00375
00376 if (!alpha.Black()) {
00377
00378 bool specularPath = false;
00379 Intersection photonIsect;
00380 int nIntersections = 0;
00381 while (scene->Intersect(photonRay, &photonIsect)) {
00382 ++nIntersections;
00383
00384
00385 alpha *= scene->Transmittance(photonRay);
00386 Vector wo = -photonRay.d;
00387
00388 BSDF *photonBSDF = photonIsect.GetBSDF(photonRay, lux::random::floatValue());
00389 BxDFType specularType = BxDFType(BSDF_REFLECTION |
00390 BSDF_TRANSMISSION | BSDF_SPECULAR);
00391 bool hasNonSpecular = (photonBSDF->NumComponents() >
00392 photonBSDF->NumComponents(specularType));
00393 if (hasNonSpecular) {
00394
00395 EPhoton photon(photonIsect.dg.p, alpha, wo);
00396 if (nIntersections == 1) {
00397 if (finalGather && (!directDone)) {
00398
00399 directPhotons.push_back(photon);
00400
00401
00402 if (directPhotons.size() == maxDirectPhotons) {
00403 directDone = true;
00404 nDirectPaths = nshot;
00405 }
00406 }
00407 } else {
00408
00409 if (specularPath) {
00410
00411 if (!causticDone) {
00412 causticPhotons.push_back(photon);
00413
00414 if (causticPhotons.size() == nCausticPhotons) {
00415 causticDone = true;
00416 nCausticPaths = nshot;
00417 causticMap = new KdTree<EPhoton, EPhotonProcess > (causticPhotons);
00418 }
00419 }
00420 } else {
00421
00422 if (!indirectDone) {
00423 indirectPhotons.push_back(photon);
00424
00425 if (indirectPhotons.size() == nIndirectPhotons) {
00426 indirectDone = true;
00427 nIndirectPaths = nshot;
00428 indirectMap = new KdTree<EPhoton, EPhotonProcess > (indirectPhotons);
00429 }
00430 }
00431 }
00432 }
00433
00434 if (finalGather && (!directDone) &&
00435 (lux::random::floatValue() < .125f)) {
00436
00437
00438 Normal n = photonIsect.dg.nn;
00439 if (Dot(n, photonRay.d) > 0.f) n = -n;
00440 radiancePhotons.push_back(ERadiancePhoton(photonIsect.dg.p, n));
00441 SWCSpectrum rho_r = photonBSDF->rho(BSDF_ALL_REFLECTION);
00442 rpReflectances.push_back(rho_r);
00443 SWCSpectrum rho_t = photonBSDF->rho(BSDF_ALL_TRANSMISSION);
00444 rpTransmittances.push_back(rho_t);
00445 }
00446 }
00447
00448
00449 Vector wi;
00450 float pdf;
00451 BxDFType flags;
00452
00453 float u1, u2, u3;
00454 if (nIntersections == 1) {
00455 u1 = RadicalInverse(nshot, 13);
00456 u2 = RadicalInverse(nshot, 17);
00457 u3 = RadicalInverse(nshot, 19);
00458 } else {
00459 u1 = lux::random::floatValue();
00460 u2 = lux::random::floatValue();
00461 u3 = lux::random::floatValue();
00462 }
00463
00464
00465 SWCSpectrum fr = photonBSDF->Sample_f(wo, &wi, u1, u2, u3,
00466 &pdf, BSDF_ALL, &flags);
00467 if (fr.Black() || pdf == 0.f)
00468 break;
00469 SWCSpectrum anew = alpha * fr *
00470 AbsDot(wi, photonBSDF->dgShading.nn) / pdf;
00471 float continueProb = min<float>(1.0f, anew.filter() / alpha.filter());
00472 if (lux::random::floatValue() > continueProb || nIntersections > 10)
00473 break;
00474 alpha = anew / continueProb;
00475 specularPath = (nIntersections == 1 || specularPath) &&
00476 ((flags & BSDF_SPECULAR) != 0);
00477 photonRay = RayDifferential(photonIsect.dg.p, wi);
00478 }
00479 }
00480
00481 BSDF::FreeAll();
00482 }
00483
00484 if (finalGather) {
00485 luxError(LUX_NOERROR, LUX_INFO, "Precompute radiance");
00486
00487
00488 KdTree<EPhoton, EPhotonProcess> directMap(directPhotons);
00489 if (!directDone)
00490 nDirectPaths = nshot;
00491
00492 for (u_int i = 0; i < radiancePhotons.size(); ++i) {
00493
00494 boost::xtime currentTime;
00495 boost::xtime_get(¤tTime, boost::TIME_UTC);
00496 if (currentTime.sec - lastUpdateTime.sec > 5) {
00497 ss.str("");
00498 ss << "Photon: " << i << " (" << (100 * i / radiancePhotons.size()) << "%)" ;
00499 luxError(LUX_NOERROR, LUX_INFO, ss.str().c_str());
00500
00501 lastUpdateTime = currentTime;
00502 }
00503
00504
00505 ERadiancePhoton &rp = radiancePhotons[i];
00506 const SWCSpectrum &rho_r = rpReflectances[i];
00507 const SWCSpectrum &rho_t = rpTransmittances[i];
00508 SWCSpectrum E;
00509 Point p = rp.p;
00510 Normal n = rp.n;
00511
00512 if (!rho_r.Black()) {
00513 E = estimateE(&directMap, nDirectPaths, p, n) +
00514 estimateE(indirectMap, nIndirectPaths, p, n) +
00515 estimateE(causticMap, nCausticPaths, p, n);
00516 rp.Lo += (E * INV_PI * rho_r).ToXYZ();
00517 }
00518
00519 if (!rho_t.Black()) {
00520 E = estimateE(&directMap, nDirectPaths, p, -n) +
00521 estimateE(indirectMap, nIndirectPaths, p, -n) +
00522 estimateE(causticMap, nCausticPaths, p, -n);
00523 rp.Lo += (E * INV_PI * rho_t).ToXYZ();
00524 }
00525 }
00526
00527 radianceMap = new KdTree<ERadiancePhoton,
00528 ERadiancePhotonProcess > (radiancePhotons);
00529 }
00530
00531 luxError(LUX_NOERROR, LUX_INFO, "Photon shooting done");
00532 }
00533
00534 SWCSpectrum ExPhotonIntegrator::Li(const Scene *scene,
00535 const RayDifferential &ray, const Sample *sample,
00536 float *alpha) const {
00537 SampleGuard guard(sample->sampler, sample);
00538
00539 SWCSpectrum L = LiInternal(0, scene, ray, sample, alpha);
00540 sample->AddContribution(sample->imageX, sample->imageY,
00541 L.ToXYZ(), alpha ? *alpha : 1.0f);
00542
00543 return SWCSpectrum(-1.f);
00544 }
00545
00546 SWCSpectrum ExPhotonIntegrator::LiInternal(
00547 const int specularDepth,
00548 const Scene *scene,
00549 const RayDifferential &ray, const Sample *sample,
00550 float *alpha) const {
00551
00552 SWCSpectrum L(0.0f);
00553
00554 if (!indirectMap)
00555 return L;
00556
00557 Intersection isect;
00558 if (scene->Intersect(ray, &isect)) {
00559
00560 float *sampleData = sample->sampler->GetLazyValues(const_cast<Sample *>(sample), sampleOffset, specularDepth);
00561 float *lightSample = &sampleData[0];
00562 float *bsdfSample = &sampleData[2];
00563 float *bsdfComponent = &sampleData[4];
00564
00565 float *reflectionSample = &sampleData[5];
00566 float *reflectionComponent = &sampleData[7];
00567 float *transmissionSample = &sampleData[8];
00568 float *transmissionComponent = &sampleData[10];
00569
00570 float *lightNum;
00571 if (lightStrategy == SAMPLE_ONE_UNIFORM)
00572 lightNum = &sampleData[11];
00573 else
00574 lightNum = NULL;
00575
00576 if (alpha) *alpha = 1.0f;
00577 Vector wo = -ray.d;
00578
00579
00580 L += isect.Le(wo);
00581
00582
00583 BSDF *bsdf = isect.GetBSDF(ray, fabsf(2.f *
00584 bsdfComponent[0]) - 1.0f);
00585 const Point &p = bsdf->dgShading.p;
00586 const Normal &n = bsdf->dgShading.nn;
00587
00588
00589 if (isnan(p.x) || isnan(p.y) || isnan(p.z)) {
00590
00591 std::stringstream ss;
00592 ss << "Internal error in photonmap, received a NaN point in differential shading geometry: point = (" << p << ")";
00593 luxError(LUX_SYSTEM, LUX_ERROR, ss.str().c_str());
00594 return L;
00595 }
00596
00597 if (debugEnableDirect) {
00598
00599 switch (lightStrategy) {
00600 case SAMPLE_ALL_UNIFORM:
00601 L += UniformSampleAllLights(scene, p, n,
00602 wo, bsdf, sample,
00603 lightSample, bsdfSample, bsdfComponent);
00604 break;
00605 case SAMPLE_ONE_UNIFORM:
00606 L += UniformSampleOneLight(scene, p, n,
00607 wo, bsdf, sample,
00608 lightSample, lightNum, bsdfSample, bsdfComponent);
00609 break;
00610 default:
00611 break;
00612 }
00613 }
00614
00615 if (debugEnableCaustic) {
00616
00617 L += LPhoton(causticMap, nCausticPaths, nLookup, bsdf,
00618 isect, wo, maxDistSquared);
00619 }
00620
00621 if (debugEnableIndirect) {
00622 if (finalGather) {
00623
00624 BxDFType nonSpecular = BxDFType(BSDF_REFLECTION |
00625 BSDF_TRANSMISSION | BSDF_DIFFUSE | BSDF_GLOSSY);
00626 if (bsdf->NumComponents(nonSpecular) > 0) {
00627
00628 u_int nIndirSamplePhotons = 50;
00629 EPhotonProcess proc(nIndirSamplePhotons, p);
00630 proc.photons = (EClosePhoton *) alloca(nIndirSamplePhotons *
00631 sizeof (EClosePhoton));
00632 float searchDist2 = maxDistSquared;
00633
00634 int sanityCheckIndex = 0;
00635 while (proc.foundPhotons < nIndirSamplePhotons) {
00636 float md2 = searchDist2;
00637 proc.foundPhotons = 0;
00638 indirectMap->Lookup(p, proc, md2);
00639
00640 searchDist2 *= 2.f;
00641
00642 if(sanityCheckIndex++ > 32) {
00643
00644 std::stringstream ss;
00645 ss << "Internal error in photonmap: point = (" <<
00646 p << ") searchDist2 = " << searchDist2;
00647 luxError(LUX_SYSTEM, LUX_ERROR, ss.str().c_str());
00648 break;
00649 }
00650 }
00651
00652
00653 Vector *photonDirs = (Vector *) alloca(nIndirSamplePhotons *
00654 sizeof (Vector));
00655 for (u_int i = 0; i < nIndirSamplePhotons; ++i)
00656 photonDirs[i] = proc.photons[i].photon->wi;
00657
00658
00659 SWCSpectrum Li = 0.;
00660
00661 for (int i = 0; i < gatherSamples; ++i) {
00662 float *sampleFGData = sample->sampler->GetLazyValues(
00663 const_cast<Sample *>(sample), sampleFinalGather1Offset, i);
00664
00665
00666 Vector wi;
00667 float u1 = sampleFGData[0];
00668 float u2 = sampleFGData[1];
00669 float u3 = sampleFGData[2];
00670 float pdf;
00671 SWCSpectrum fr = bsdf->Sample_f(wo, &wi, u1, u2, u3,
00672 &pdf, BxDFType(BSDF_ALL & (~BSDF_SPECULAR)));
00673 if (fr.Black() || pdf == 0.f) continue;
00674
00675
00676 if (gatherRRStrategy == RR_EFFICIENCY) {
00677 const float dp = AbsDot(wi, n) / pdf;
00678 const float q = min<float>(1.0f, fr.filter() * dp);
00679 if (q < sampleFGData[3])
00680 continue;
00681
00682
00683 fr /= q;
00684 } else if (gatherRRStrategy == RR_PROBABILITY) {
00685 if (gatherRRContinueProbability < sampleFGData[3])
00686 continue;
00687
00688
00689 fr /= gatherRRContinueProbability;
00690 }
00691
00692 RayDifferential bounceRay(p, wi);
00693
00694 Intersection gatherIsect;
00695 if (scene->Intersect(bounceRay, &gatherIsect)) {
00696
00697 SWCSpectrum Lindir = 0.f;
00698 Normal nGather = gatherIsect.dg.nn;
00699 if (Dot(nGather, bounceRay.d) > 0) nGather = -nGather;
00700 ERadiancePhotonProcess proc(gatherIsect.dg.p, nGather);
00701 float md2 = INFINITY;
00702
00703 radianceMap->Lookup(gatherIsect.dg.p, proc, md2);
00704 if (proc.photon) {
00705 Lindir = FromXYZ(proc.photon->Lo.c[0],
00706 proc.photon->Lo.c[1],
00707 proc.photon->Lo.c[2]);
00708 }
00709
00710 Lindir *= scene->Transmittance(bounceRay);
00711
00712
00713 float photonPdf = 0.f;
00714 float conePdf = UniformConePdf(cosGatherAngle);
00715 for (u_int j = 0; j < nIndirSamplePhotons; ++j)
00716 if (Dot(photonDirs[j], wi) > .999f * cosGatherAngle)
00717 photonPdf += conePdf;
00718 photonPdf /= nIndirSamplePhotons;
00719 float wt = PowerHeuristic(gatherSamples, pdf,
00720 gatherSamples, photonPdf);
00721 Li += fr * Lindir * AbsDot(wi, n) * wt / pdf;
00722 }
00723 }
00724 L += Li / gatherSamples;
00725
00726
00727 Li = 0.;
00728 for (int i = 0; i < gatherSamples; ++i) {
00729 float *sampleFGData = sample->sampler->GetLazyValues(
00730 const_cast<Sample *>(sample), sampleFinalGather2Offset, i);
00731
00732
00733 float u1 = sampleFGData[2];
00734 float u2 = sampleFGData[0];
00735 float u3 = sampleFGData[1];
00736 int photonNum = min((int) nIndirSamplePhotons - 1,
00737 Floor2Int(u1 * nIndirSamplePhotons));
00738
00739 Vector vx, vy;
00740 CoordinateSystem(photonDirs[photonNum], &vx, &vy);
00741 Vector wi = UniformSampleCone(u2, u3, cosGatherAngle, vx, vy,
00742 photonDirs[photonNum]);
00743
00744 SWCSpectrum fr = bsdf->f(wo, wi);
00745 if (fr.Black()) continue;
00746
00747
00748 float photonPdf = 0.f;
00749 float conePdf = UniformConePdf(cosGatherAngle);
00750 for (u_int j = 0; j < nIndirSamplePhotons; ++j)
00751 if (Dot(photonDirs[j], wi) > .999f * cosGatherAngle)
00752 photonPdf += conePdf;
00753 photonPdf /= nIndirSamplePhotons;
00754
00755
00756 if (gatherRRStrategy == RR_EFFICIENCY) {
00757 const float dp = 1.0f / photonPdf;
00758 const float q = min<float>(1.0f, fr.filter() * dp);
00759 if (q < sampleFGData[3])
00760 continue;
00761
00762
00763 fr /= q;
00764 } else if (gatherRRStrategy == RR_PROBABILITY) {
00765 if (gatherRRContinueProbability < sampleFGData[3])
00766 continue;
00767
00768
00769 fr /= gatherRRContinueProbability;
00770 }
00771
00772 RayDifferential bounceRay(p, wi);
00773
00774 Intersection gatherIsect;
00775 if (scene->Intersect(bounceRay, &gatherIsect)) {
00776
00777 SWCSpectrum Lindir = 0.f;
00778 Normal nGather = gatherIsect.dg.nn;
00779 if (Dot(nGather, bounceRay.d) > 0) nGather = -nGather;
00780 ERadiancePhotonProcess proc(gatherIsect.dg.p, nGather);
00781 float md2 = INFINITY;
00782
00783 radianceMap->Lookup(gatherIsect.dg.p, proc, md2);
00784 if (proc.photon) {
00785 Lindir = FromXYZ(proc.photon->Lo.c[0],
00786 proc.photon->Lo.c[1],
00787 proc.photon->Lo.c[2]);
00788 }
00789
00790 Lindir *= scene->Transmittance(bounceRay);
00791
00792 float bsdfPdf = bsdf->Pdf(wo, wi);
00793 float wt = PowerHeuristic(gatherSamples, photonPdf,
00794 gatherSamples, bsdfPdf);
00795 Li += fr * Lindir * AbsDot(wi, n) * wt / photonPdf;
00796 }
00797 }
00798
00799 L += Li / gatherSamples;
00800 }
00801 } else {
00802 L += LPhoton(indirectMap, nIndirectPaths, nLookup,
00803 bsdf, isect, wo, maxDistSquared);
00804 }
00805 }
00806
00807 if (debugEnableSpecular && (specularDepth < maxSpecularDepth)) {
00808 float u1 = reflectionSample[0];
00809 float u2 = reflectionSample[1];
00810 float u3 = reflectionComponent[0];
00811
00812 Vector wi;
00813 float pdf;
00814
00815 SWCSpectrum f = bsdf->Sample_f(wo, &wi, u1, u2, u3,
00816 &pdf, BxDFType(BSDF_REFLECTION | BSDF_SPECULAR));
00817 if ((!f.Black()) || (pdf > 0.0f)) {
00818
00819 RayDifferential rd(p, wi);
00820 rd.hasDifferentials = true;
00821 rd.rx.o = p + isect.dg.dpdx;
00822 rd.ry.o = p + isect.dg.dpdy;
00823
00824 Normal dndx = bsdf->dgShading.dndu * bsdf->dgShading.dudx +
00825 bsdf->dgShading.dndv * bsdf->dgShading.dvdx;
00826 Normal dndy = bsdf->dgShading.dndu * bsdf->dgShading.dudy +
00827 bsdf->dgShading.dndv * bsdf->dgShading.dvdy;
00828 Vector dwodx = -ray.rx.d - wo, dwody = -ray.ry.d - wo;
00829 float dDNdx = Dot(dwodx, n) + Dot(wo, dndx);
00830 float dDNdy = Dot(dwody, n) + Dot(wo, dndy);
00831 rd.rx.d = wi -
00832 dwodx + 2 * Vector(Dot(wo, n) * dndx +
00833 dDNdx * n);
00834 rd.ry.d = wi -
00835 dwody + 2 * Vector(Dot(wo, n) * dndy +
00836 dDNdy * n);
00837 L += LiInternal(specularDepth + 1, scene, rd, sample, alpha) *
00838 (f * (1.0f / pdf)) * AbsDot(wi, n);
00839 }
00840
00841 u1 = transmissionSample[0];
00842 u2 = transmissionSample[1];
00843 u3 = transmissionComponent[0];
00844
00845 f = bsdf->Sample_f(wo, &wi, u1, u2, u3,
00846 &pdf, BxDFType(BSDF_TRANSMISSION | BSDF_SPECULAR));
00847 if ((!f.Black()) || (pdf > 0.0f)) {
00848
00849 RayDifferential rd(p, wi);
00850 rd.hasDifferentials = true;
00851 rd.rx.o = p + isect.dg.dpdx;
00852 rd.ry.o = p + isect.dg.dpdy;
00853
00854 float eta = bsdf->eta;
00855 Vector w = -wo;
00856 if (Dot(wo, n) < 0) eta = 1.f / eta;
00857
00858 Normal dndx = bsdf->dgShading.dndu * bsdf->dgShading.dudx + bsdf->dgShading.dndv * bsdf->dgShading.dvdx;
00859 Normal dndy = bsdf->dgShading.dndu * bsdf->dgShading.dudy + bsdf->dgShading.dndv * bsdf->dgShading.dvdy;
00860
00861 Vector dwodx = -ray.rx.d - wo, dwody = -ray.ry.d - wo;
00862 float dDNdx = Dot(dwodx, n) + Dot(wo, dndx);
00863 float dDNdy = Dot(dwody, n) + Dot(wo, dndy);
00864
00865 float mu = eta * Dot(w, n) - Dot(wi, n);
00866 float dmudx = (eta - (eta * eta * Dot(w, n)) / Dot(wi, n)) * dDNdx;
00867 float dmudy = (eta - (eta * eta * Dot(w, n)) / Dot(wi, n)) * dDNdy;
00868
00869 rd.rx.d = wi + eta * dwodx - Vector(mu * dndx + dmudx * n);
00870 rd.ry.d = wi + eta * dwody - Vector(mu * dndy + dmudy * n);
00871 L += LiInternal(specularDepth + 1, scene, rd, sample, alpha) *
00872 (f * (1.0f / pdf)) * AbsDot(wi, n);
00873 }
00874 }
00875 } else {
00876
00877 if (alpha) *alpha = 0.0f;
00878 for (u_int i = 0; i < scene->lights.size(); ++i)
00879 L += scene->lights[i]->Le(ray);
00880 if (alpha && !L.Black()) *alpha = 1.0f;
00881 }
00882
00883 return L * scene->volumeIntegrator->Transmittance(scene, ray, sample, alpha) +
00884 scene->volumeIntegrator->Li(scene, ray, sample, alpha);
00885 }
00886
00887 SurfaceIntegrator* ExPhotonIntegrator::CreateSurfaceIntegrator(const ParamSet ¶ms) {
00888 LightStrategy estrategy;
00889 string st = params.FindOneString("strategy", "auto");
00890 if (st == "one") estrategy = SAMPLE_ONE_UNIFORM;
00891 else if (st == "all") estrategy = SAMPLE_ALL_UNIFORM;
00892 else if (st == "auto") estrategy = SAMPLE_AUTOMATIC;
00893 else {
00894 std::stringstream ss;
00895 ss<<"Strategy '"<<st<<"' for direct lighting unknown. Using \"auto\".";
00896 luxError(LUX_BADTOKEN,LUX_WARNING,ss.str().c_str());
00897 estrategy = SAMPLE_AUTOMATIC;
00898 }
00899
00900 int nCaustic = params.FindOneInt("causticphotons", 20000);
00901 int nIndirect = params.FindOneInt("indirectphotons", 200000);
00902 int maxDirect = params.FindOneInt("maxdirectphotons", 1000000);
00903
00904 int nUsed = params.FindOneInt("nused", 50);
00905 int maxDepth = params.FindOneInt("maxdepth", 4);
00906
00907 bool finalGather = params.FindOneBool("finalgather", true);
00908 int gatherSamples = params.FindOneInt("finalgathersamples", 32);
00909
00910 float maxDist = params.FindOneFloat("maxdist", 0.5f);
00911 float gatherAngle = params.FindOneFloat("gatherangle", 10.0f);
00912
00913 RRStrategy rstrategy;
00914 string rst = params.FindOneString("gatherrrstrategy", "efficiency");
00915 if (rst == "efficiency") rstrategy = RR_EFFICIENCY;
00916 else if (rst == "probability") rstrategy = RR_PROBABILITY;
00917 else if (rst == "none") rstrategy = RR_NONE;
00918 else {
00919 std::stringstream ss;
00920 ss<<"Strategy '"<<st<<"' for russian roulette path termination unknown. Using \"efficiency\".";
00921 luxError(LUX_BADTOKEN,LUX_WARNING,ss.str().c_str());
00922 rstrategy = RR_EFFICIENCY;
00923 }
00924
00925 float rrcontinueProb = params.FindOneFloat("gatherrrcontinueprob", 0.65f);
00926
00927 bool debugEnableDirect = params.FindOneBool("dbg_enabledirect", true);
00928 bool debugEnableCaustic = params.FindOneBool("dbg_enablecaustic", true);
00929 bool debugEnableIndirect = params.FindOneBool("dbg_enableindirect", true);
00930 bool debugEnableSpecular = params.FindOneBool("dbg_enablespecular", true);
00931
00932 return new ExPhotonIntegrator(estrategy, nCaustic, nIndirect, maxDirect,
00933 nUsed, maxDepth, maxDist, finalGather, gatherSamples, gatherAngle,
00934 rstrategy, rrcontinueProb,
00935 debugEnableDirect, debugEnableCaustic, debugEnableIndirect, debugEnableSpecular);
00936 }