00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023 #include "photonmap.h"
00024 #include "light.h"
00025 #include "mc.h"
00026 #include "spectrumwavelengths.h"
00027 #include "error.h"
00028 #include "osfunc.h"
00029 #include "mc.h"
00030
00031 #include <fstream>
00032 #include <boost/thread/xtime.hpp>
00033
00034 using namespace lux;
00035
00036 namespace lux
00037 {
00038
00039 SWCSpectrum BasicColorPhoton::GetSWCSpectrum(const TsPack* tspack, u_int nb) const
00040 {
00041 SWCSpectrum result(0.f);
00042 for (u_int i = 0; i < WAVELENGTH_SAMPLES; ++i) {
00043 if (i > 0) {
00044 const float delta0 = 1.f - (fabsf(tspack->swl->w[i] - w[i - 1]) * WAVELENGTH_SAMPLES) /
00045 ((WAVELENGTH_END - WAVELENGTH_START) * nb);
00046 if (delta0 > 0.f)
00047 result.c[i] += delta0 * alpha.c[i - 1];
00048 }
00049 const float delta1 = 1.f - (fabsf(tspack->swl->w[i] - w[i]) * WAVELENGTH_SAMPLES) /
00050 ((WAVELENGTH_END - WAVELENGTH_START) * nb);
00051 if (delta1 > 0.f)
00052 result.c[i] += delta1 * alpha.c[i];
00053 if (i < WAVELENGTH_SAMPLES - 1) {
00054 const float delta2 = 1.f - (fabsf(tspack->swl->w[i] - w[i + 1]) * WAVELENGTH_SAMPLES) /
00055 ((WAVELENGTH_END - WAVELENGTH_START) * nb);
00056 if (delta2 > 0.f)
00057 result.c[i] += delta2 * alpha.c[i + 1];
00058 }
00059 }
00060 return result;
00061 }
00062
00063 void BasicColorPhoton::save(bool isLittleEndian, std::basic_ostream<char> &stream) const
00064 {
00065
00066 for (u_int i = 0; i < 3; ++i)
00067 osWriteLittleEndianFloat(isLittleEndian, stream, p[i]);
00068
00069
00070 for (u_int i = 0; i < WAVELENGTH_SAMPLES; ++i)
00071 osWriteLittleEndianFloat(isLittleEndian, stream, alpha.c[i]);
00072
00073
00074 for (u_int i = 0; i < WAVELENGTH_SAMPLES; ++i)
00075 osWriteLittleEndianFloat(isLittleEndian, stream, w[i]);
00076 }
00077
00078 void BasicColorPhoton::load(bool isLittleEndian, std::basic_istream<char> &stream)
00079 {
00080
00081 for (u_int i = 0; i < 3; ++i)
00082 p[i] = osReadLittleEndianFloat(isLittleEndian, stream);
00083
00084
00085 for (u_int i = 0; i < WAVELENGTH_SAMPLES; ++i)
00086 alpha.c[i] = osReadLittleEndianFloat(isLittleEndian, stream);
00087
00088
00089 for (u_int i = 0; i < WAVELENGTH_SAMPLES; ++i)
00090 w[i] = osReadLittleEndianFloat(isLittleEndian, stream);
00091 }
00092
00093 void LightPhoton::save(bool isLittleEndian, std::basic_ostream<char> &stream) const
00094 {
00095
00096 BasicColorPhoton::save(isLittleEndian, stream);
00097
00098
00099 for (u_int i = 0; i < 3; ++i)
00100 osWriteLittleEndianFloat(isLittleEndian, stream, wi[i]);
00101 }
00102
00103 void LightPhoton::load(bool isLittleEndian, std::basic_istream<char> &stream)
00104 {
00105
00106 BasicColorPhoton::load(isLittleEndian, stream);
00107
00108
00109 for (u_int i = 0; i < 3; ++i)
00110 w[i] = osReadLittleEndianFloat(isLittleEndian, stream);
00111 }
00112
00113 void RadiancePhoton::save(bool isLittleEndian, std::basic_ostream<char> &stream) const
00114 {
00115
00116 BasicColorPhoton::save(isLittleEndian, stream);
00117
00118
00119 for (u_int i = 0; i < 3; ++i)
00120 osWriteLittleEndianFloat(isLittleEndian, stream, n[i]);
00121 }
00122
00123 void RadiancePhoton::load(bool isLittleEndian, std::basic_istream<char> &stream)
00124 {
00125
00126 BasicColorPhoton::load(isLittleEndian, stream);
00127
00128
00129 for (u_int i = 0; i < 3; ++i)
00130 n[i] = osReadLittleEndianFloat(isLittleEndian, stream);
00131 }
00132
00133 SWCSpectrum RadiancePhotonMap::LPhoton(const TsPack *tspack,
00134 const Intersection& isect, const Vector& wo,
00135 const BxDFType bxdfType) const
00136 {
00137 SWCSpectrum L(0.f);
00138 if (!photonmap)
00139 return L;
00140
00141 Normal ng = isect.dg.nn;
00142 if (Dot(wo, ng) < 0.f)
00143 ng = -ng;
00144
00145 const Point& p = isect.dg.p;
00146
00147 if ((bxdfType & BSDF_REFLECTION) != 0) {
00148
00149 NearPhotonProcess<RadiancePhoton> procRefl(p, ng);
00150 float md2Refl = maxDistSquared;
00151 lookup(p, procRefl, md2Refl);
00152 if (procRefl.photon)
00153 L += procRefl.photon->GetSWCSpectrum(tspack, 1);
00154 }
00155
00156 if ((bxdfType & BSDF_TRANSMISSION ) != 0) {
00157
00158 NearPhotonProcess<RadiancePhoton> procTransm(p, -ng);
00159 float md2Transm = maxDistSquared;
00160 lookup(p, procTransm, md2Transm);
00161 if (procTransm.photon)
00162 L += procTransm.photon->GetSWCSpectrum(tspack, 1);
00163 }
00164
00165 return L;
00166 }
00167
00168 SWCSpectrum LightPhotonMap::EPhoton(const TsPack *tspack, const Point &p,
00169 const Normal &n) const
00170 {
00171 SWCSpectrum E(0.f);
00172 if ((nPaths <= 0) || (!photonmap))
00173 return E;
00174
00175
00176 NearSetPhotonProcess<LightPhoton> proc(nLookup, p);
00177 proc.photons =
00178 (ClosePhoton<LightPhoton> *) alloca(nLookup * sizeof (ClosePhoton<LightPhoton>));
00179 float md2 = maxDistSquared;
00180 lookup(p, proc, md2);
00181
00182
00183 ClosePhoton<LightPhoton> *photons = proc.photons;
00184 for (u_int i = 0; i < proc.foundPhotons; ++i) {
00185 if (Dot(n, photons[i].photon->wi) > 0.f)
00186 E += photons[i].photon->GetSWCSpectrum(tspack, proc.foundPhotons);
00187 }
00188
00189 return E / (nPaths * md2 * M_PI);
00190 }
00191
00192 SWCSpectrum LightPhotonMap::LPhoton(const TsPack* tspack, const BSDF *bsdf,
00193 const Intersection &isect, const Vector &wo,
00194 const BxDFType bxdfType) const
00195 {
00196 SWCSpectrum L(0.f);
00197 if ((nPaths <= 0) || (!photonmap))
00198 return L;
00199
00200 if (bsdf->NumComponents(bxdfType) == 0)
00201 return L;
00202
00203
00204 NearSetPhotonProcess<LightPhoton> proc(nLookup, isect.dg.p);
00205 proc.photons =
00206 (ClosePhoton<LightPhoton> *) alloca(nLookup * sizeof (ClosePhoton<LightPhoton>));
00207
00208 float md2 = maxDistSquared;
00209 lookup(isect.dg.p, proc, md2);
00210
00211
00212 const ClosePhoton<LightPhoton> *photons = proc.photons;
00213 const u_int nFound = proc.foundPhotons;
00214 const Normal Nf = Dot(wo, isect.dg.nn) < 0 ? -isect.dg.nn : isect.dg.nn;
00215 const float distSquared = md2;
00216
00217
00218 for (u_int i = 0; i < nFound; ++i) {
00219 const LightPhoton *p = photons[i].photon;
00220 BxDFType flag = BxDFType(bxdfType &
00221 (Dot(Nf, p->wi) > 0.f ?
00222 BSDF_ALL_REFLECTION : BSDF_ALL_TRANSMISSION));
00223 float k = Ekernel(p, isect.dg.p, distSquared);
00224 const SWCSpectrum alpha = p->GetSWCSpectrum(tspack, nFound);
00225
00226 L += (k / nPaths) * bsdf->f(tspack, p->wi, wo, flag) * alpha;
00227 }
00228
00229 return L;
00230 }
00231
00232 SWCSpectrum LightPhotonMap::LPhotonDiffuseApprox(const TsPack* tspack,
00233 const BSDF *bsdf, const Intersection &isect, const Vector &wo,
00234 const BxDFType bxdfType) const
00235 {
00236 SWCSpectrum L(0.f);
00237 if ((nPaths <= 0) || (!photonmap))
00238 return L;
00239
00240 if (bsdf->NumComponents(bxdfType) == 0)
00241 return L;
00242
00243
00244 NearSetPhotonProcess<LightPhoton> proc(nLookup, isect.dg.p);
00245 proc.photons =
00246 (ClosePhoton<LightPhoton> *) alloca(nLookup * sizeof (ClosePhoton<LightPhoton>));
00247
00248 float md2 = maxDistSquared;
00249 lookup(isect.dg.p, proc, md2);
00250
00251
00252 const ClosePhoton<LightPhoton> *photons = proc.photons;
00253 const u_int nFound = proc.foundPhotons;
00254 const Normal Nf = Dot(wo, isect.dg.nn) < 0 ? -isect.dg.nn : isect.dg.nn;
00255 const float distSquared = md2;
00256
00257
00258 SWCSpectrum Lr(0.f), Lt(0.f);
00259
00260 for (u_int i = 0; i < nFound; ++i) {
00261 const LightPhoton *p = photons[i].photon;
00262 const SWCSpectrum alpha = p->GetSWCSpectrum(tspack, nFound);
00263
00264 float k = Ekernel(p, isect.dg.p, distSquared);
00265 if (Dot(Nf, photons[i].photon->wi) > 0.f)
00266 Lr += (k / nPaths) * alpha;
00267 else
00268 Lt += (k / nPaths) * alpha;
00269 }
00270
00271 if ((bxdfType & BSDF_REFLECTION) != 0)
00272 L += Lr * bsdf->rho(tspack, wo, BxDFType(bxdfType & BSDF_ALL_REFLECTION)) * INV_PI;
00273 if ((bxdfType & BSDF_TRANSMISSION) != 0)
00274 L += Lt * bsdf->rho(tspack, wo, BxDFType(bxdfType & BSDF_ALL_TRANSMISSION)) * INV_PI;
00275
00276 return L;
00277 }
00278
00279 SWCSpectrum LightPhotonMap::LDiffusePhoton(const TsPack* tspack,
00280 const BSDF *bsdf, const Intersection &isect, const Vector &wo) const
00281 {
00282 SWCSpectrum L(0.f);
00283 if ((nPaths <= 0) || (!photonmap))
00284 return L;
00285
00286 const BxDFType diffuseType = BxDFType(BSDF_DIFFUSE | BSDF_REFLECTION | BSDF_TRANSMISSION);
00287
00288 if (bsdf->NumComponents(diffuseType) == 0)
00289 return L;
00290
00291
00292 NearSetPhotonProcess<LightPhoton> proc(nLookup, isect.dg.p);
00293 proc.photons =
00294 (ClosePhoton<LightPhoton> *) alloca(nLookup * sizeof (ClosePhoton<LightPhoton>));
00295
00296 float md2 = maxDistSquared;
00297 lookup(isect.dg.p, proc, md2);
00298
00299
00300 const ClosePhoton<LightPhoton> *photons = proc.photons;
00301 const u_int nFound = proc.foundPhotons;
00302 const Normal Nf = Dot(wo, isect.dg.nn) < 0 ? -isect.dg.nn : isect.dg.nn;
00303 const float distSquared = md2;
00304
00305
00306 SWCSpectrum Lr(0.f), Lt(0.f);
00307
00308 for (u_int i = 0; i < nFound; ++i) {
00309 const LightPhoton *p = photons[i].photon;
00310 const SWCSpectrum alpha = p->GetSWCSpectrum(tspack, nFound);
00311
00312 float k = Ekernel(p, isect.dg.p, distSquared);
00313 if (Dot(Nf, photons[i].photon->wi) > 0.f)
00314 Lr += (k / nPaths) * alpha;
00315 else
00316 Lt += (k / nPaths) * alpha;
00317 }
00318
00319 L = Lr * bsdf->rho(tspack, wo, BxDFType(BSDF_DIFFUSE | BSDF_REFLECTION)) * INV_PI +
00320 Lt * bsdf->rho(tspack, wo, BxDFType(BSDF_DIFFUSE | BSDF_TRANSMISSION)) * INV_PI;
00321
00322 return L;
00323 }
00324
00325 static bool unsuccessful(int needed, int found, int shot)
00326 {
00327 return (found < needed && (found == 0 || found < shot / 1024));
00328 }
00329
00330 void PhotonMapPreprocess(const TsPack *tspack, const Scene *scene,
00331 const string *mapFileName, const BxDFType photonBxdfType,
00332 const BxDFType radianceBxdfType, u_int nDirectPhotons,
00333 u_int nRadiancePhotons, RadiancePhotonMap *radianceMap,
00334 u_int nIndirectPhotons, LightPhotonMap *indirectMap,
00335 u_int nCausticPhotons, LightPhotonMap *causticMap,
00336 u_int maxDepth)
00337 {
00338 if (scene->lights.size() == 0)
00339 return;
00340
00341 std::stringstream ss;
00342
00343
00344 if (mapFileName) {
00345
00346 std::ifstream ifs(mapFileName->c_str(), std::ios_base::in | std::ios_base::binary);
00347
00348 if (ifs.good()) {
00349 ss.str("");
00350 ss << "Found photon maps file: " << *mapFileName;
00351 luxError(LUX_NOERROR, LUX_INFO, ss.str().c_str());
00352
00353 bool isLittleEndian = osIsLittleEndian();
00354
00355 bool ok = true;
00356
00357
00358 int storedPhotonBxdfType;
00359 int storedRadianceBxdfType;
00360 u_int storedNDirectPhotons;
00361 u_int storedNRadiancePhotons;
00362 u_int storedNIndirectPhotons;
00363 u_int storedNCausticPhotons;
00364 storedPhotonBxdfType = osReadLittleEndianInt(isLittleEndian, ifs);
00365 storedRadianceBxdfType = osReadLittleEndianInt(isLittleEndian, ifs);
00366 storedNDirectPhotons = osReadLittleEndianUInt(isLittleEndian, ifs);
00367 storedNRadiancePhotons = osReadLittleEndianUInt(isLittleEndian, ifs);
00368 storedNIndirectPhotons = osReadLittleEndianUInt(isLittleEndian, ifs);
00369 storedNCausticPhotons = osReadLittleEndianUInt(isLittleEndian, ifs);
00370 if (storedPhotonBxdfType != photonBxdfType ||
00371 storedRadianceBxdfType != radianceBxdfType ||
00372 storedNDirectPhotons != nDirectPhotons ||
00373 storedNRadiancePhotons != nRadiancePhotons ||
00374 storedNIndirectPhotons != nIndirectPhotons ||
00375 storedNCausticPhotons != nCausticPhotons) {
00376 luxError(LUX_NOERROR, LUX_INFO, "Some photon map settings changed, rebuilding photon maps...");
00377 ok = false;
00378 }
00379 if (ok) {
00380
00381 u_int storedNLights;
00382 storedNLights = osReadLittleEndianUInt(isLittleEndian, ifs);
00383 if (storedNLights != scene->lights.size() ) {
00384 luxError(LUX_NOERROR, LUX_INFO, "Scene changed, rebuilding photon maps...");
00385 ok = false;
00386 }
00387 }
00388
00389
00390 if (ok) {
00391 luxError(LUX_NOERROR, LUX_INFO, "Reading radiance photon map...");
00392 RadiancePhotonMap::load(ifs, radianceMap);
00393 ss.str("");
00394 ss << "Read " << radianceMap->getPhotonCount() << " radiance photons";
00395 luxError(LUX_NOERROR, LUX_INFO, ss.str().c_str());
00396
00397 if (!ifs.good()) {
00398 luxError(LUX_NOERROR, LUX_INFO, "Failed to read all photon maps");
00399 ok = false;
00400 }
00401 }
00402 if (ok) {
00403 luxError(LUX_NOERROR, LUX_INFO, "Reading indirect photon map...");
00404 LightPhotonMap::load(ifs, indirectMap);
00405 ss.str("");
00406 ss << "Read " << indirectMap->getPhotonCount() << " light photons";
00407 luxError(LUX_NOERROR, LUX_INFO, ss.str().c_str());
00408
00409 if (!ifs.good()) {
00410 luxError(LUX_NOERROR, LUX_INFO, "Failed to read all photon maps");
00411 ok = false;
00412 }
00413 }
00414 if (ok) {
00415 luxError(LUX_NOERROR, LUX_INFO, "Reading caustic photon map...");
00416 LightPhotonMap::load(ifs, causticMap);
00417 ss.str("");
00418 ss << "Read " << causticMap->getPhotonCount() << " light photons";
00419 luxError(LUX_NOERROR, LUX_INFO, ss.str().c_str());
00420
00421 if (!ifs.good()) {
00422 luxError(LUX_NOERROR, LUX_INFO, "Failed to read all photon maps");
00423 ok = false;
00424 }
00425 }
00426
00427
00428 ifs.close();
00429
00430
00431 if (ok)
00432 return;
00433 } else {
00434 luxError(LUX_NOERROR, LUX_INFO, "Photon maps file doesn't exist yet");
00435 ifs.close();
00436 }
00437 }
00438
00439
00440 bool computeRadianceMap = (nRadiancePhotons > 0);
00441
00442
00443 ss.str("");
00444 ss << "Shooting photons (target: " << (nCausticPhotons + nIndirectPhotons) << ")...";
00445 luxError(LUX_NOERROR, LUX_INFO, ss.str().c_str());
00446
00447 vector<LightPhoton> directPhotons;
00448 directPhotons.reserve(nDirectPhotons);
00449 bool directDone = (nDirectPhotons == 0);
00450
00451 vector<LightPhoton> causticPhotons;
00452 causticPhotons.reserve(nCausticPhotons);
00453 bool causticDone = (nCausticPhotons == 0);
00454
00455 vector<LightPhoton> indirectPhotons;
00456 indirectPhotons.reserve(nIndirectPhotons);
00457 bool indirectDone = (nIndirectPhotons == 0);
00458
00459 vector<RadiancePhoton> radiancePhotons;
00460 radiancePhotons.reserve(nRadiancePhotons);
00461 bool radianceDone = (nRadiancePhotons == 0);
00462
00463
00464 SpectrumWavelengths *thr_wl = tspack->swl;
00465
00466
00467 u_int nLights = int(scene->lights.size());
00468 float *lightPower = (float *)alloca(nLights * sizeof(float));
00469 float *lightCDF = (float *)alloca((nLights + 1) * sizeof(float));
00470
00471
00472 const u_int spectrumSamples = 128;
00473 for (u_int i = 0; i < nLights; ++i)
00474 lightPower[i] = 0.f;
00475 for (u_int j = 0; j < spectrumSamples; ++j) {
00476 thr_wl->Sample(RadicalInverse(j, 2), RadicalInverse(j, 3));
00477
00478 for (u_int i = 0; i < nLights; ++i)
00479 lightPower[i] += scene->lights[i]->Power(tspack, scene).Y(tspack);
00480 }
00481 for (u_int i = 0; i < nLights; ++i)
00482 lightPower[i] /= spectrumSamples;
00483
00484 float totalPower;
00485 ComputeStep1dCDF(lightPower, nLights, &totalPower, lightCDF);
00486
00487
00488 vector<SWCSpectrum> rpReflectances;
00489 rpReflectances.reserve(nRadiancePhotons);
00490 vector<SWCSpectrum> rpTransmittances;
00491 rpTransmittances.reserve(nRadiancePhotons);
00492
00493 boost::xtime photonShootingStartTime;
00494 boost::xtime lastUpdateTime;
00495 boost::xtime_get(&photonShootingStartTime, boost::TIME_UTC);
00496 boost::xtime_get(&lastUpdateTime, boost::TIME_UTC);
00497 u_int nshot = 0;
00498 while (!radianceDone || !directDone || !causticDone || !indirectDone) {
00499
00500 boost::xtime currentTime;
00501 boost::xtime_get(¤tTime, boost::TIME_UTC);
00502 if (currentTime.sec - lastUpdateTime.sec > 5) {
00503 ss.str("");
00504 ss << "Photon shooting progress: Direct[" << directPhotons.size();
00505 if (nDirectPhotons > 0)
00506 ss << " (" << (100 * directPhotons.size() / nDirectPhotons) << "% limit)";
00507 else
00508 ss << " (100% limit)";
00509 ss << "] Caustic[" << causticPhotons.size();
00510 if (nCausticPhotons > 0)
00511 ss << " (" << (100 * causticPhotons.size() / nCausticPhotons) << "%)";
00512 else
00513 ss << " (100%)";
00514 ss << "] Indirect[" << indirectPhotons.size();
00515 if (nIndirectPhotons > 0)
00516 ss << " (" << (100 * indirectPhotons.size() / nIndirectPhotons) << "%)";
00517 else
00518 ss << " (100%)";
00519 ss << "] Radiance[" << radiancePhotons.size();
00520 if (nRadiancePhotons > 0)
00521 ss << " (" << (100 * radiancePhotons.size() / nRadiancePhotons) << "% limit)";
00522 else
00523 ss << " (100% limit)";
00524 ss << "]";
00525 luxError(LUX_NOERROR, LUX_INFO, ss.str().c_str());
00526
00527 lastUpdateTime = currentTime;
00528 }
00529
00530 ++nshot;
00531
00532
00533 if (nshot > 500000) {
00534 if (indirectDone && unsuccessful(nCausticPhotons, causticPhotons.size(), nshot)) {
00535
00536
00537 luxError(LUX_CONSISTENCY, LUX_WARNING, "Unable to store enough photons in the caustic photonmap. Giving up and disabling the map.");
00538
00539 causticPhotons.clear();
00540 causticDone = true;
00541 nCausticPhotons = 0;
00542 }
00543
00544 if (unsuccessful(nIndirectPhotons, indirectPhotons.size(), nshot)) {
00545 luxError(LUX_CONSISTENCY, LUX_ERROR, "Unable to store enough photons in the indirect photonmap. Unable to render the image.");
00546 return;
00547 }
00548 }
00549
00550 thr_wl->Sample(RadicalInverse(nshot, 2), RadicalInverse(nshot, 3));
00551
00552
00553
00554 float u[4];
00555 u[0] = RadicalInverse(nshot, 5);
00556 u[1] = RadicalInverse(nshot, 7);
00557 u[2] = RadicalInverse(nshot, 11);
00558 u[3] = RadicalInverse(nshot, 13);
00559
00560
00561 float lightPdf;
00562 float uln = RadicalInverse(nshot, 17);
00563 u_int lightNum = Floor2Int(SampleStep1d(lightPower, lightCDF,
00564 totalPower, nLights, uln, &lightPdf) * nLights);
00565 lightNum = min(lightNum, nLights - 1);
00566 const Light *light = scene->lights[lightNum];
00567
00568
00569 RayDifferential photonRay;
00570 float pdf;
00571 SWCSpectrum alpha = light->Sample_L(tspack, scene,
00572 u[0], u[1], u[2], u[3], &photonRay, &pdf);
00573 if (pdf == 0.f || alpha.Black())
00574 continue;
00575 alpha /= pdf * lightPdf;
00576
00577 if (!alpha.Black()) {
00578
00579 bool specularPath = false;
00580 Intersection photonIsect;
00581 u_int nIntersections = 0;
00582 while (scene->Intersect(photonRay, &photonIsect)) {
00583 ++nIntersections;
00584
00585
00586
00587
00588 Vector wo = -photonRay.d;
00589
00590 BSDF *photonBSDF = photonIsect.GetBSDF(tspack, photonRay);
00591 if (photonBSDF->NumComponents(photonBxdfType) > 0) {
00592
00593 LightPhoton photon(tspack, photonIsect.dg.p, alpha, wo);
00594
00595 if (nIntersections == 1) {
00596 if (computeRadianceMap && (!directDone)) {
00597
00598 directPhotons.push_back(photon);
00599
00600
00601 if (directPhotons.size() == nDirectPhotons)
00602 directDone = true;
00603 }
00604 } else {
00605
00606 if (specularPath) {
00607
00608 if (!causticDone) {
00609 causticPhotons.push_back(photon);
00610
00611 if (causticPhotons.size() == nCausticPhotons) {
00612 causticDone = true;
00613 causticMap->init(nshot, causticPhotons);
00614 }
00615 }
00616 } else {
00617
00618 if (!indirectDone) {
00619 indirectPhotons.push_back(photon);
00620
00621 if (indirectPhotons.size() == nIndirectPhotons) {
00622 indirectDone = true;
00623 indirectMap->init(nshot, indirectPhotons);
00624 }
00625 }
00626 }
00627 }
00628
00629 if (computeRadianceMap &&
00630 (!radianceDone) &&
00631 (photonBSDF->NumComponents(radianceBxdfType) > 0) &&
00632 (tspack->rng->floatValue() < 0.125f)) {
00633 SWCSpectrum rho_t =
00634 photonBSDF->rho(tspack, BxDFType(radianceBxdfType & BSDF_ALL_TRANSMISSION));
00635 SWCSpectrum rho_r =
00636 photonBSDF->rho(tspack, BxDFType(radianceBxdfType & BSDF_ALL_REFLECTION));
00637
00638 if(!rho_t.Black() || !rho_r.Black()) {
00639
00640 Normal n = photonIsect.dg.nn;
00641 if (Dot(n, photonRay.d) > 0.f)
00642 n = -n;
00643 radiancePhotons.push_back(RadiancePhoton(tspack, photonIsect.dg.p, n));
00644
00645 rpReflectances.push_back(rho_r);
00646 rpTransmittances.push_back(rho_t);
00647
00648 if (radiancePhotons.size() == nRadiancePhotons)
00649 radianceDone = true;
00650 }
00651 }
00652 }
00653
00654
00655 Vector wi;
00656 float pdfo;
00657 BxDFType flags;
00658
00659 float u1, u2, u3;
00660 if (nIntersections == 1) {
00661 u1 = RadicalInverse(nshot, 19);
00662 u2 = RadicalInverse(nshot, 23);
00663 u3 = RadicalInverse(nshot, 29);
00664 } else {
00665 u1 = tspack->rng->floatValue();
00666 u2 = tspack->rng->floatValue();
00667 u3 = tspack->rng->floatValue();
00668 }
00669
00670
00671 SWCSpectrum fr;
00672 if (!photonBSDF->Sample_f(tspack, wo, &wi, u1, u2, u3, &fr, &pdfo, BSDF_ALL, &flags))
00673 break;
00674 SWCSpectrum anew = alpha * fr * AbsDot(wi, photonIsect.dg.nn) * AbsDot(wo, photonBSDF->dgShading.nn) / (AbsDot(wo, photonIsect.dg.nn) * pdfo);
00675 float continueProb = min(1.f, anew.Filter(tspack) / alpha.Filter(tspack));
00676 if (tspack->rng->floatValue() > continueProb || nIntersections > maxDepth)
00677 break;
00678 alpha = anew / continueProb;
00679 specularPath = (nIntersections == 1 || specularPath) &&
00680 ((flags & BSDF_SPECULAR) != 0 || pdfo > 100.f);
00681 photonRay = RayDifferential(photonIsect.dg.p, wi);
00682 }
00683 }
00684
00685 BSDF::FreeAll(tspack);
00686 }
00687
00688 boost::xtime photonShootingEndTime;
00689 boost::xtime_get(&photonShootingEndTime, boost::TIME_UTC);
00690 ss.str("");
00691 ss << "Photon shooting done (" << ( photonShootingEndTime.sec - photonShootingStartTime.sec ) << "s)";
00692 luxError(LUX_NOERROR, LUX_INFO, ss.str().c_str());
00693
00694 if (computeRadianceMap) {
00695 luxError(LUX_NOERROR, LUX_INFO, "Computing radiance photon map...");
00696
00697
00698 LightPhotonMap directMap(radianceMap->nLookup, radianceMap->maxDistSquared);
00699 directMap.init(nDirectPhotons, directPhotons);
00700
00701 for (u_int i = 0; i < radiancePhotons.size(); ++i) {
00702
00703 boost::xtime currentTime;
00704 boost::xtime_get(¤tTime, boost::TIME_UTC);
00705 if (currentTime.sec - lastUpdateTime.sec > 5) {
00706 ss.str("");
00707 ss << "Radiance photon map computation progress: " << i << " (" << (100 * i / radiancePhotons.size()) << "%)";
00708 luxError(LUX_NOERROR, LUX_INFO, ss.str().c_str());
00709
00710 lastUpdateTime = currentTime;
00711 }
00712
00713
00714 RadiancePhoton &rp = radiancePhotons[i];
00715 const SWCSpectrum &rho_r = rpReflectances[i];
00716 const SWCSpectrum &rho_t = rpTransmittances[i];
00717 const Point& p = rp.p;
00718 const Normal& n = rp.n;
00719 SWCSpectrum alpha(0.f);
00720 for (u_int j = 0; j < WAVELENGTH_SAMPLES; ++j)
00721 tspack->swl->w[j] = rp.w[j];
00722
00723 if (!rho_r.Black()) {
00724 SWCSpectrum E = directMap.EPhoton(tspack, p, n);
00725 E += indirectMap->EPhoton(tspack, p, n);
00726 E += causticMap->EPhoton(tspack, p, n);
00727
00728 alpha += E * INV_PI * rho_r;
00729 }
00730
00731 if (!rho_t.Black()) {
00732 SWCSpectrum E = directMap.EPhoton(tspack, p, -n);
00733 E += indirectMap->EPhoton(tspack, p, -n);
00734 E += causticMap->EPhoton(tspack, p, -n);
00735
00736 alpha += E * INV_PI * rho_t;
00737 }
00738
00739 rp.alpha = alpha;
00740 }
00741
00742 radianceMap->init(radiancePhotons);
00743
00744
00745 boost::xtime radianceComputeEndTime;
00746 boost::xtime_get(&radianceComputeEndTime, boost::TIME_UTC);
00747 ss.str("");
00748 ss << "Radiance photon map computed (" << ( radianceComputeEndTime.sec - photonShootingEndTime.sec ) << "s)";
00749 luxError(LUX_NOERROR, LUX_INFO, ss.str().c_str());
00750 }
00751
00752
00753 if (mapFileName) {
00754 luxError(LUX_NOERROR, LUX_INFO, "Saving photon maps to file" );
00755
00756 std::ofstream ofs(mapFileName->c_str(), std::ios_base::out | std::ios_base::binary);
00757 if(ofs.good()) {
00758 ss.str("");
00759 ss << "Writing photon maps to '" << (*mapFileName) << "'...";
00760 luxError(LUX_NOERROR, LUX_INFO, ss.str().c_str());
00761
00762 bool isLittleEndian = osIsLittleEndian();
00763
00764
00765 osWriteLittleEndianUInt(isLittleEndian, ofs, photonBxdfType);
00766 osWriteLittleEndianUInt(isLittleEndian, ofs, radianceBxdfType);
00767 osWriteLittleEndianUInt(isLittleEndian, ofs, nDirectPhotons);
00768 osWriteLittleEndianUInt(isLittleEndian, ofs, nRadiancePhotons);
00769 osWriteLittleEndianUInt(isLittleEndian, ofs, nIndirectPhotons);
00770 osWriteLittleEndianUInt(isLittleEndian, ofs, nCausticPhotons);
00771 osWriteLittleEndianUInt(isLittleEndian, ofs, scene->lights.size());
00772
00773
00774
00775 if (radianceMap) {
00776 radianceMap->save(ofs);
00777
00778 ss.str("");
00779 ss << "Written " << radianceMap->getPhotonCount() << " radiance photons";
00780 luxError(LUX_NOERROR, LUX_INFO, ss.str().c_str());
00781 } else
00782 osWriteLittleEndianInt(isLittleEndian, ofs, 0);
00783
00784
00785 if (indirectMap) {
00786 indirectMap->save(ofs);
00787
00788 ss.str("");
00789 ss << "Written " << indirectMap->getPhotonCount() << " indirect photons";
00790 luxError(LUX_NOERROR, LUX_INFO, ss.str().c_str());
00791 } else
00792 osWriteLittleEndianInt(isLittleEndian, ofs, 0);
00793
00794
00795 if (causticMap) {
00796 causticMap->save(ofs);
00797
00798 ss.str("");
00799 ss << "Written " << causticMap->getPhotonCount() << " caustic photons";
00800 luxError(LUX_NOERROR, LUX_INFO, ss.str().c_str());
00801 } else
00802 osWriteLittleEndianInt(isLittleEndian, ofs, 0);
00803
00804 if(!ofs.good()) {
00805 ss.str("");
00806 ss << "Error while writing photon maps to file";
00807 luxError(LUX_SYSTEM, LUX_SEVERE, ss.str().c_str());
00808 }
00809
00810 ofs.close();
00811 } else {
00812 ss.str("");
00813 ss << "Cannot open file '" << (*mapFileName) << "' for writing photon maps";
00814 luxError(LUX_SYSTEM, LUX_SEVERE, ss.str().c_str());
00815 }
00816 }
00817 }
00818
00819 SWCSpectrum PhotonMapFinalGatherWithImportaceSampling(const TsPack* tspack,
00820 const Scene *scene, const Sample *sample, int sampleFinalGather1Offset,
00821 int sampleFinalGather2Offset, int gatherSamples, float cosGatherAngle,
00822 PhotonMapRRStrategy rrStrategy, float rrContinueProbability,
00823 const LightPhotonMap *indirectMap, const RadiancePhotonMap *radianceMap,
00824 const Vector &wo, const BSDF *bsdf, const BxDFType bxdfType)
00825 {
00826 SWCSpectrum L(0.f);
00827
00828
00829 if (bsdf->NumComponents(bxdfType) > 0) {
00830 const Point &p = bsdf->dgShading.p;
00831 const Normal &n = bsdf->dgShading.nn;
00832
00833
00834 u_int nIndirSamplePhotons = indirectMap->nLookup;
00835 NearSetPhotonProcess<LightPhoton> proc(nIndirSamplePhotons, p);
00836 proc.photons =
00837 (ClosePhoton<LightPhoton> *) alloca(nIndirSamplePhotons * sizeof(ClosePhoton<LightPhoton>));
00838 float searchDist2 = indirectMap->maxDistSquared;
00839
00840 int sanityCheckIndex = 0;
00841 while (proc.foundPhotons < nIndirSamplePhotons) {
00842 float md2 = searchDist2;
00843 proc.foundPhotons = 0;
00844 indirectMap->lookup(p, proc, md2);
00845
00846 searchDist2 *= 2.f;
00847
00848 if (sanityCheckIndex++ > 32) {
00849
00850 std::stringstream ss;
00851 ss << "Internal error in photonmap: point = (" <<
00852 p << ") searchDist2 = " << searchDist2;
00853 luxError(LUX_SYSTEM, LUX_ERROR, ss.str().c_str());
00854 break;
00855 }
00856 }
00857
00858
00859 Vector *photonDirs = (Vector *) alloca(nIndirSamplePhotons * sizeof(Vector));
00860 for (u_int i = 0; i < nIndirSamplePhotons; ++i)
00861 photonDirs[i] = proc.photons[i].photon->wi;
00862
00863 const float scaledCosGatherAngle = 0.999f * cosGatherAngle;
00864
00865 SWCSpectrum Li(0.f);
00866 for (int i = 0; i < gatherSamples ; ++i) {
00867 float *sampleFGData = sample->sampler->GetLazyValues(
00868 const_cast<Sample *>(sample), sampleFinalGather1Offset, i);
00869
00870
00871 Vector wi;
00872 float u1 = sampleFGData[0];
00873 float u2 = sampleFGData[1];
00874 float u3 = sampleFGData[2];
00875 float pdf;
00876 SWCSpectrum fr;
00877 if (!bsdf->Sample_f(tspack, wo, &wi, u1, u2, u3, &fr, &pdf, bxdfType, NULL, NULL, true))
00878 continue;
00879
00880
00881 if (rrStrategy == RR_EFFICIENCY) {
00882 const float dp = AbsDot(wi, n) / pdf;
00883 const float q = min(1.f, fr.Filter(tspack) * dp);
00884 if (q < sampleFGData[3])
00885 continue;
00886
00887
00888 fr /= q;
00889 } else if (rrStrategy == RR_PROBABILITY) {
00890 if (rrContinueProbability < sampleFGData[3])
00891 continue;
00892
00893
00894 fr /= rrContinueProbability;
00895 }
00896
00897
00898 RayDifferential bounceRay(p, wi);
00899 Intersection gatherIsect;
00900 if (scene->Intersect(bounceRay, &gatherIsect)) {
00901
00902 Normal nGather = gatherIsect.dg.nn;
00903 if (Dot(nGather, bounceRay.d) > 0)
00904 nGather = -nGather;
00905 NearPhotonProcess<RadiancePhoton> procir(gatherIsect.dg.p, nGather);
00906 float md2 = radianceMap->maxDistSquared;
00907
00908 radianceMap->lookup(gatherIsect.dg.p, procir, md2);
00909 if (procir.photon) {
00910 SWCSpectrum Lindir = procir.photon->GetSWCSpectrum(tspack, 1);
00911
00912 scene->Transmittance(tspack, bounceRay, sample, &Lindir);
00913
00914
00915 float photonPdf = 0.f;
00916 float conePdf = UniformConePdf(cosGatherAngle);
00917 for (u_int j = 0; j < nIndirSamplePhotons; ++j) {
00918 if (Dot(photonDirs[j], wi) > scaledCosGatherAngle)
00919 photonPdf += conePdf;
00920 }
00921 photonPdf /= nIndirSamplePhotons;
00922 float wt = PowerHeuristic(gatherSamples, pdf, gatherSamples, photonPdf);
00923
00924 if (bounceRay.maxt < sqrtf(md2))
00925 wt *= (1.f + bounceRay.maxt / sqrtf(md2)) / 2.f;
00926 Li += fr * Lindir * (AbsDot(wi, n) * wt / pdf);
00927 }
00928 }
00929 }
00930 L += Li / gatherSamples;
00931
00932
00933 Li = 0.f;
00934 for (int i = 0; i < gatherSamples; ++i) {
00935 float *sampleFGData = sample->sampler->GetLazyValues(
00936 const_cast<Sample *>(sample), sampleFinalGather2Offset, i);
00937
00938
00939 float u1 = sampleFGData[2];
00940 float u2 = sampleFGData[0];
00941 float u3 = sampleFGData[1];
00942 int photonNum = min<u_int>(nIndirSamplePhotons - 1,
00943 Floor2Int(u1 * nIndirSamplePhotons));
00944
00945 Vector vx, vy;
00946 CoordinateSystem(photonDirs[photonNum], &vx, &vy);
00947 Vector wi = UniformSampleCone(u2, u3, cosGatherAngle, vx, vy,
00948 photonDirs[photonNum]);
00949
00950 SWCSpectrum fr = bsdf->f(tspack, wi, wo, bxdfType);
00951 if (fr.Black())
00952 continue;
00953
00954
00955 float photonPdf = 0.f;
00956 float conePdf = UniformConePdf(cosGatherAngle);
00957 for (u_int j = 0; j < nIndirSamplePhotons; ++j) {
00958 if (Dot(photonDirs[j], wi) > scaledCosGatherAngle)
00959 photonPdf += conePdf;
00960 }
00961 photonPdf /= nIndirSamplePhotons;
00962
00963
00964 if (rrStrategy == RR_EFFICIENCY) {
00965 const float dp = 1.f / photonPdf;
00966 const float q = min(1.f, fr.Filter(tspack) * dp);
00967 if (q < sampleFGData[3])
00968 continue;
00969
00970
00971 fr /= q;
00972 } else if (rrStrategy == RR_PROBABILITY) {
00973 if (rrContinueProbability < sampleFGData[3])
00974 continue;
00975
00976
00977 fr /= rrContinueProbability;
00978 }
00979
00980 RayDifferential bounceRay(p, wi);
00981 Intersection gatherIsect;
00982 if (scene->Intersect(bounceRay, &gatherIsect)) {
00983
00984 Normal nGather = gatherIsect.dg.nn;
00985 if (Dot(nGather, bounceRay.d) > 0)
00986 nGather = -nGather;
00987 NearPhotonProcess<RadiancePhoton> procir(gatherIsect.dg.p, nGather);
00988 float md2 = radianceMap->maxDistSquared;
00989
00990 radianceMap->lookup(gatherIsect.dg.p, procir, md2);
00991 if (procir.photon) {
00992 SWCSpectrum Lindir = procir.photon->GetSWCSpectrum(tspack, 1);
00993
00994 scene->Transmittance(tspack, bounceRay, sample, &Lindir);
00995
00996 float bsdfPdf = bsdf->Pdf(tspack, wi, wo, bxdfType);
00997 float wt = PowerHeuristic(gatherSamples, photonPdf, gatherSamples, bsdfPdf);
00998
00999 if (bounceRay.maxt < sqrtf(md2))
01000 wt *= (1.f + bounceRay.maxt / sqrtf(md2)) / 2.f;
01001 Li += fr * Lindir * (AbsDot(wi, n) * wt / photonPdf);
01002 }
01003 }
01004 }
01005
01006 L += Li / gatherSamples;
01007 }
01008
01009 return L;
01010 }
01011
01012 SWCSpectrum PhotonMapFinalGather(const TsPack *tspack, const Scene *scene,
01013 const Sample *sample, int sampleFinalGatherOffset, int gatherSamples,
01014 PhotonMapRRStrategy rrStrategy, float rrContinueProbability,
01015 const LightPhotonMap *indirectMap, const RadiancePhotonMap *radianceMap,
01016 const Vector &wo, const BSDF *bsdf, const BxDFType bxdfType)
01017 {
01018 SWCSpectrum L(0.f);
01019
01020
01021 if (bsdf->NumComponents(bxdfType) > 0) {
01022 const Point &p = bsdf->dgShading.p;
01023 const Normal &n = bsdf->dgShading.nn;
01024
01025
01026 SWCSpectrum Li(0.f);
01027 for (int i = 0; i < gatherSamples ; ++i) {
01028 float *sampleFGData = sample->sampler->GetLazyValues(
01029 const_cast<Sample *>(sample), sampleFinalGatherOffset, i);
01030
01031
01032 Vector wi;
01033 float u1 = sampleFGData[0];
01034 float u2 = sampleFGData[1];
01035 float u3 = sampleFGData[2];
01036 float pdf;
01037 SWCSpectrum fr;
01038 if (!bsdf->Sample_f(tspack, wo, &wi, u1, u2, u3, &fr, &pdf, bxdfType, NULL, NULL, true))
01039 continue;
01040
01041
01042 if (rrStrategy == RR_EFFICIENCY) {
01043 const float dp = AbsDot(wi, n) / pdf;
01044 const float q = min(1.f, fr.Filter(tspack) * dp);
01045 if (q < sampleFGData[3])
01046 continue;
01047
01048
01049 fr /= q;
01050 } else if (rrStrategy == RR_PROBABILITY) {
01051 if (rrContinueProbability < sampleFGData[3])
01052 continue;
01053
01054
01055 fr /= rrContinueProbability;
01056 }
01057
01058
01059 RayDifferential bounceRay(p, wi);
01060 Intersection gatherIsect;
01061 if (scene->Intersect(bounceRay, &gatherIsect)) {
01062
01063 Normal nGather = gatherIsect.dg.nn;
01064 if (Dot(nGather, bounceRay.d) > 0)
01065 nGather = -nGather;
01066 NearPhotonProcess<RadiancePhoton> proc(gatherIsect.dg.p, nGather);
01067 float md2 = radianceMap->maxDistSquared;
01068
01069 radianceMap->lookup(gatherIsect.dg.p, proc, md2);
01070 if (proc.photon) {
01071 SWCSpectrum Lindir = proc.photon->GetSWCSpectrum(tspack, 1);
01072 scene->Transmittance(tspack, bounceRay, sample, &Lindir);
01073
01074 Li += fr * Lindir * (AbsDot(wi, n) / pdf);
01075 }
01076 }
01077 }
01078 L += Li / gatherSamples;
01079 }
01080
01081 return L;
01082 }
01083
01084 void LightPhotonMap::load(std::basic_istream<char> &stream, LightPhotonMap *map)
01085 {
01086 if (!map)
01087 return;
01088
01089 bool isLittleEndian = osIsLittleEndian();
01090
01091
01092 int count;
01093 count = osReadLittleEndianInt(isLittleEndian, stream);
01094
01095 int npaths;
01096 npaths = osReadLittleEndianInt(isLittleEndian, stream);
01097
01098 vector<LightPhoton> photons(count);
01099 for (int i = 0; i < count; ++i)
01100 photons[i].load(isLittleEndian, stream);
01101
01102 if (count > 0)
01103 map->init(npaths, photons);
01104 }
01105
01106 void LightPhotonMap::save(std::basic_ostream<char> &stream) const
01107 {
01108 bool isLittleEndian = osIsLittleEndian();
01109
01110
01111 osWriteLittleEndianInt(isLittleEndian, stream, photonCount);
01112 osWriteLittleEndianInt(isLittleEndian, stream, nPaths);
01113
01114 if (photonmap != NULL) {
01115 LightPhoton *photons = photonmap->getNodeData();
01116 for (int i = 0; i < photonCount; i++)
01117 photons[i].save(isLittleEndian, stream);
01118 }
01119 }
01120
01121 void RadiancePhotonMap::load(std::basic_istream<char> &stream, RadiancePhotonMap *map)
01122 {
01123 if(!map)
01124 return;
01125
01126 bool isLittleEndian = osIsLittleEndian();
01127
01128
01129 int count;
01130 count = osReadLittleEndianInt(isLittleEndian, stream);
01131
01132 vector<RadiancePhoton> photons(count);
01133 for (int i = 0; i < count; ++i)
01134 photons[i].load(isLittleEndian, stream);
01135
01136 if (count > 0)
01137 map->init(photons);
01138 }
01139
01140 void RadiancePhotonMap::save(std::basic_ostream<char> &stream) const
01141 {
01142 bool isLittleEndian = osIsLittleEndian();
01143
01144
01145 osWriteLittleEndianInt(isLittleEndian, stream, photonCount);
01146
01147 if (photonmap != NULL) {
01148 RadiancePhoton *photons = photonmap->getNodeData();
01149 for (int i = 0; i < photonCount; ++i)
01150 photons[i].save(isLittleEndian, stream);
01151 }
01152 }
01153
01154 }