improved particle tracer

include/mitsuba/render/gatherproc.h

@@ -74,13 +74,13 @@ public:
 	 * of excess photons that had to be discarded. If this is too
 	 * high, the granularity should be decreased.
 	 */
-	inline size_t getExcess() const { return m_excess; }
-	
+	inline size_t getExcessPhotons() const { return m_excess; }
+
 	/**
-	 * \brief Lists the nuber of photons that had to be shot
+	 * \brief Lists the nuber of particles that had to be shot
 	 * in order to fill the photon map.
 	 */
-	inline size_t getShotPhotons() const { return m_numShot; }
+	inline size_t getShotParticles() const { return m_numShot; }
 
 	// ======================================================================
 	/// @{ \name ParallelProcess implementation

include/mitsuba/render/integrator.h

@@ -196,6 +196,17 @@ public:
 		scene = parent.scene;
 		sampler = parent.sampler;
 		depth = parent.depth+1;
+		medium = parent.medium;
+		extra = 0;
+	}
+
+	/// Initialize the query record for a recursive query
+	inline void recursiveQuery(const RadianceQueryRecord &parent) {
+		type = parent.type;
+		scene = parent.scene;
+		sampler = parent.sampler;
+		depth = parent.depth+1;
+		medium = parent.medium;
 		extra = 0;
 	}
 

src/integrators/misc/irrcache.cpp

@@ -27,8 +27,8 @@ MTS_NAMESPACE_BEGIN
  * whenever a Lambertian surface is intersected, an internal irradiance
  * cache is queried for the indirect illumination at the surface position in 
  * question. If this query is successful, the sub-integrator is only 
- * used to compute the remaining types of radiance (direct, in-scatter, 
- * emission) and their sum is returned afterwards.
+ * used to compute the remaining types of illumination (direct, volume, and
+ * emitted radiance) and their sum is returned afterwards.
  * When a query is unsuccessful, a new data point is generated by a final
  * gathering step.
  *
@@ -229,7 +229,7 @@ public:
 		if (rRec.rayIntersect(ray)) {
 			const BSDF *bsdf = its.getBSDF(ray);
 
-			if (bsdf->getType() == BSDF::EDiffuseReflection && 
+			if (bsdf && bsdf->getType() == BSDF::EDiffuseReflection && 
 					(rRec.type & RadianceQueryRecord::EIndirectSurfaceRadiance)) {
 				Spectrum E;
 				if (!m_irrCache->get(its, E)) {

src/integrators/photonmapper/photonmapper.cpp

@@ -165,11 +165,11 @@ public:
 				return false;
 
 			m_globalPhotonMap = proc->getPhotonMap();
-			m_globalPhotonMap->setScale(1 / (Float) proc->getShotPhotons());
+			m_globalPhotonMap->setScale(1 / (Float) proc->getShotParticles());
 			m_globalPhotonMap->setMinPhotons(m_globalMinPhotons);
 			m_globalPhotonMap->balance();
-			Log(EDebug, "Global photon map full. Shot " SIZE_T_FMT " photons, excess due to parallelism: " 
-				SIZE_T_FMT, proc->getShotPhotons(), proc->getExcess());
+			Log(EDebug, "Global photon map full. Shot " SIZE_T_FMT " particles, excess photons due to parallelism: " 
+				SIZE_T_FMT, proc->getShotParticles(), proc->getExcessPhotons());
 			m_globalPhotonMapID = sched->registerResource(m_globalPhotonMap);
 		}
 
@@ -183,6 +183,7 @@ public:
 				m_granularity, 2, m_rrDepth, m_gatherLocally, job);
 
 			proc->bindResource("scene", sceneResID);
+			proc->bindResource("camera", cameraResID);
 			proc->bindResource("sampler", qmcSamplerID);
 
 			m_proc = proc;
@@ -194,11 +195,11 @@ public:
 				return false;
 
 			m_causticPhotonMap = proc->getPhotonMap();
-			m_causticPhotonMap->setScale(1 / (Float) proc->getShotPhotons());
+			m_causticPhotonMap->setScale(1 / (Float) proc->getShotParticles());
 			m_causticPhotonMap->setMinPhotons(m_causticMinPhotons);
 			m_causticPhotonMap->balance();
-			Log(EDebug, "Caustic photon map - excess photons due to parallelism: " 
-				SIZE_T_FMT, proc->getExcess());
+			Log(EDebug, "Caustic photon map full. Shot " SIZE_T_FMT " particles, excess photons due to parallelism: " 
+				SIZE_T_FMT, proc->getShotParticles(), proc->getExcessPhotons());
 			m_causticPhotonMapID = sched->registerResource(m_causticPhotonMap);
 		}
 
@@ -212,6 +213,7 @@ public:
 				m_granularity, m_maxDepth, m_rrDepth, m_gatherLocally, job);
 
 			proc->bindResource("scene", sceneResID);
+			proc->bindResource("camera", cameraResID);
 			proc->bindResource("sampler", qmcSamplerID);
 
 			m_proc = proc;
@@ -223,11 +225,11 @@ public:
 				return false;
 
 			m_volumePhotonMap = proc->getPhotonMap();
-			m_volumePhotonMap->setScale(1 / (Float) proc->getShotPhotons());
+			m_volumePhotonMap->setScale(1 / (Float) proc->getShotParticles());
 			m_volumePhotonMap->setMinPhotons(m_volumeMinPhotons);
 			m_volumePhotonMap->balance();
-			Log(EDebug, "Volumetric photon map - excess photons due to parallelism: " 
-				SIZE_T_FMT, proc->getExcess());
+			Log(EDebug, "Volume photon map full. Shot " SIZE_T_FMT " particles, excess photons due to parallelism: " 
+				SIZE_T_FMT, proc->getShotParticles(), proc->getExcessPhotons());
 			m_volumePhotonMapID = sched->registerResource(m_volumePhotonMap);
 		}
 
@@ -277,6 +279,9 @@ public:
 		   intersection has already been provided). */
 		rRec.rayIntersect(ray);
 
+		if ((rRec.type & RadianceQueryRecord::EVolumeRadiance) && rRec.medium) {
+		}
+
 		if (!its.isValid()) {
 			/* If no intersection could be found, possibly return 
 			   attenuated radiance from a background luminaire */
@@ -294,6 +299,16 @@ public:
 			Li += its.LoSub(rRec.scene, -ray.d);
 
 		const BSDF *bsdf = its.getBSDF(ray);
+				
+		if (its.isMediumTransition())
+			rRec.medium = its.getTargetMedium(ray.d);
+
+		if (bsdf == NULL) {
+			RadianceQueryRecord rRec2;
+			rRec2.recursiveQuery(rRec, rRec.type);
+			return m_parentIntegrator->Li(RayDifferential(its.p, ray.d, ray.time), rRec2);
+		}
+
 		int bsdfType = bsdf->getType();
 
 		Point2 *sampleArray, sample;
@@ -324,11 +339,6 @@ public:
 			}
 		}
 
-//		if (rRec.type & RadianceQueryRecord:EVolumeRadiance) {
-			/* Ray marching */
-//		}
-
-
 		if (bsdfType == BSDF::EDiffuseReflection) {
 			/* Hit a diffuse material - do a direct photon map visualization. */
 			if (rRec.type & RadianceQueryRecord::EIndirectSurfaceRadiance) 

src/integrators/photonmapper/ppm.cpp

@@ -291,11 +291,11 @@ public:
 
 		ref<PhotonMap> photonMap = proc->getPhotonMap();
 		photonMap->balance();
-		Log(EInfo, "Global photon map full. Shot " SIZE_T_FMT " photons, excess due to parallelism: " 
-			SIZE_T_FMT, proc->getShotPhotons(), proc->getExcess());
+		Log(EDebug, "Photon map full. Shot " SIZE_T_FMT " particles, excess photons due to parallelism: " 
+			SIZE_T_FMT, proc->getShotParticles(), proc->getExcessPhotons());
 
 		Log(EInfo, "Gathering ..");
-		m_totalEmitted += proc->getShotPhotons();
+		m_totalEmitted += proc->getShotParticles();
 		film->clear();
 		#pragma omp parallel for schedule(dynamic)
 		for (int blockIdx = 0; blockIdx<(int) m_blocks.size(); ++blockIdx) {

src/integrators/photonmapper/sppm.cpp

@@ -280,11 +280,11 @@ public:
 
 		ref<PhotonMap> photonMap = proc->getPhotonMap();
 		photonMap->balance();
-		Log(EInfo, "Global photon map full. Shot " SIZE_T_FMT " photons, excess due to parallelism: " 
-			SIZE_T_FMT, proc->getShotPhotons(), proc->getExcess());
+		Log(EDebug, "Photon map full. Shot " SIZE_T_FMT " particles, excess photons due to parallelism: " 
+			SIZE_T_FMT, proc->getShotParticles(), proc->getExcessPhotons());
 
 		Log(EInfo, "Gathering ..");
-		m_totalEmitted += proc->getShotPhotons();
+		m_totalEmitted += proc->getShotParticles();
 		film->clear();
 		#pragma omp parallel for schedule(dynamic)
 		for (int blockIdx = 0; blockIdx<(int) m_gatherBlocks.size(); ++blockIdx) {
@@ -309,7 +309,7 @@ public:
 				} else {
 					Float ratio = (N + m_alpha * M) / (N + M);
 					gp.flux = (gp.flux + gp.weight * (flux + 
-						gp.emission * proc->getShotPhotons() * M_PI * gp.radius*gp.radius)) * ratio;
+						gp.emission * proc->getShotParticles() * M_PI * gp.radius*gp.radius)) * ratio;
 					gp.radius = gp.radius * std::sqrt(ratio);
 					gp.N = N + m_alpha * M;
 					contrib = gp.flux / ((Float) m_totalEmitted * gp.radius*gp.radius * M_PI);

src/librender/gatherproc.cpp

@@ -195,12 +195,17 @@ void GatherPhotonProcess::processResult(const WorkResult *wr, bool cancelled) {
 		size_t start = vec.getParticleIndex(i),
 			   end   = vec.getParticleIndex(i+1);
 		++nParticles;
+		bool stop = false;
 		for (size_t j=start; j<end; ++j) {
 			if (!m_photonMap->storePhoton(vec[j])) {
+				cout << "Photon map is full, skipping " << vec.getPhotonCount() -j << " photons" << endl;
 				m_excess += vec.getPhotonCount() - j;
+				stop = true;
 				break;
 			}
 		}
+		if (stop)
+			break;
 	}
 	m_numShot += nParticles;
 	increaseResultCount(vec.getPhotonCount());

src/librender/particleproc.cpp

@@ -32,7 +32,7 @@ ParticleProcess::ParticleProcess(EMode mode, size_t workCount, size_t granularit
 	/* Choose a suitable work unit granularity if none was specified */
 	if (m_granularity == 0)
 		m_granularity = std::max((size_t) 1, workCount /
-			(8 * Scheduler::getInstance()->getWorkerCount()));
+			(16 * Scheduler::getInstance()->getWorkerCount()));
 
 	/* Create a visual progress reporter */
 	m_progress = new ProgressReporter(progressText, workCount,