homogeneous bugfix for channels with sigmaT=0, fixed ppm and sppm

include/mitsuba/core/kdtree.h

@@ -291,7 +291,7 @@ public:
 	//! @{ \name \c stl::vector-like interface
 	// =============================================================
 	/// Clear the kd-tree array
-	inline void clear() { m_nodes.clear(); }
+	inline void clear() { m_nodes.clear(); m_aabb.reset(); }
 	/// Resize the kd-tree array
 	inline void resize(size_t size) { m_nodes.resize(size); }
 	/// Reserve a certain amount of memory for the kd-tree array
@@ -301,7 +301,10 @@ public:
 	/// Return the capacity of the kd-tree
 	inline size_t capacity() const { return m_nodes.capacity(); }
 	/// Append a kd-tree node to the node array
-	inline void push_back(const NodeType &node) { m_nodes.push_back(node); }
+	inline void push_back(const NodeType &node) {
+		m_nodes.push_back(node);
+		m_aabb.expandBy(node.getPosition());
+	}
 	/// Return one of the KD-tree nodes by index
 	inline NodeType &operator[](size_t idx) { return m_nodes[idx]; }
 	/// Return one of the KD-tree nodes by index (const version)
@@ -317,18 +320,19 @@ public:
 	inline size_t getDepth() const { return m_depth; }
 
 	/// Construct the KD-tree hierarchy
-	void build() {
+	void build(bool recomputeAABB = false) {
 		ref<Timer> timer = new Timer();
 
 		SLog(EInfo, "Building a %i-dimensional kd-tree over " SIZE_T_FMT " data points",
 			PointType::dim, m_nodes.size());
-		m_aabb.reset();
 
-		for (size_t i=0; i<m_nodes.size(); ++i)
-			m_aabb.expandBy(m_nodes[i].getPosition());
-		
+		if (recomputeAABB) {
+			m_aabb.reset();
+			for (size_t i=0; i<m_nodes.size(); ++i)
+				m_aabb.expandBy(m_nodes[i].getPosition());
+
+		}
 		int aabbTime = timer->getMilliseconds();
-
 		timer->reset();
 
 		/* Instead of shuffling around the node data itself, only modify
@@ -340,23 +344,29 @@ public:
 			indirection[i] = i;
 
 		m_depth = 0;
+		int constructionTime;
 		if (NodeType::leftBalancedLayout) {
 			std::vector<IndexType> permutation(m_nodes.size());
 			buildLB(0, 1, indirection.begin(), indirection.begin(), 
 				indirection.end(), permutation);
+			constructionTime = timer->getMilliseconds();
+			timer->reset();
 			permute_inplace(&m_nodes[0], permutation);
 		} else {
 			build(1, indirection.begin(), indirection.begin(), indirection.end());
+			constructionTime = timer->getMilliseconds();
+			timer->reset();
 			permute_inplace(&m_nodes[0], indirection);
 		}
 
-		int constructionTime = timer->getMilliseconds();
-		timer->reset();
-
 		int permutationTime = timer->getMilliseconds();
-		
-		SLog(EInfo, "Done after %i ms (breakdown: aabb: %i ms, build: %i ms, permute: %i ms). ",
+
+		if (recomputeAABB)
+			SLog(EInfo, "Done after %i ms (breakdown: aabb: %i ms, build: %i ms, permute: %i ms). ",
 				aabbTime + constructionTime + permutationTime, aabbTime, constructionTime, permutationTime);
+		else
+			SLog(EInfo, "Done after %i ms (breakdown: build: %i ms, permute: %i ms). ",
+				constructionTime + permutationTime, constructionTime, permutationTime);
 	}
 
 	/**

include/mitsuba/core/thread.h

@@ -88,7 +88,9 @@ public:
 #elif MTS_USE_ELF_TLS == 1
 	inline static int getID() { return m_id; }
 #else
-	inline static int getID() { return getThread()->m_id; }
+	inline static int getID() { 
+		return getThread()->m_id;
+	}
 #endif
 
 	/// Return the name of this thread
@@ -116,9 +118,7 @@ public:
 	inline FileResolver *getFileResolver() { return m_fresolver; }
 
 	/// Return the current thread
-	inline static Thread *getThread() {
-		return m_self->get();
-	}
+	inline static Thread *getThread() { return m_self->get(); }
 
 	/// Is this thread still running?
 	inline bool isRunning() const { return m_running; }

include/mitsuba/render/photonmap.h

@@ -165,7 +165,7 @@ public:
 	 * This has to be done once after all photons have been stored, 
 	 * but prior to executing any queries.
 	 */
-	inline void build() { m_kdtree.build(); }
+	inline void build(bool recomputeAABB = false) { m_kdtree.build(recomputeAABB); }
 
 	/// Determine if the photon map is completely filled
 	inline bool isFull() const {

src/integrators/photonmapper/photonmapper.cpp

@@ -346,7 +346,7 @@ public:
 			}
 		}
 
-		if (bsdfType == BSDF::EDiffuseReflection) {
+		if ((bsdfType & BSDF::EAll) == BSDF::EDiffuseReflection) {
 			/* Hit a diffuse material - do a direct photon map visualization. */
 			if (rRec.type & RadianceQueryRecord::EIndirectSurfaceRadiance)
 				LiSurf += m_globalPhotonMap->estimateIrradiance(its.p, 

src/integrators/photonmapper/sppm.cpp

@@ -147,6 +147,8 @@ public:
 			samplers[i] = clonedSampler.get();
 		}
 
+		Thread::initializeOpenMP(Scheduler::getInstance()->getLocalWorkerCount());
+
 		int samplerResID = sched->registerManifoldResource(
 			static_cast<std::vector<SerializableObject*> &>(samplers)); 
 
@@ -219,8 +221,8 @@ public:
 							/* Create hit point if this is a diffuse material or a glossy
 							   one, and there has been a previous interaction with
 							   a glossy material */
-							if (bsdf->getType() == BSDF::EDiffuseReflection || 
-								bsdf->getType() == BSDF::EDiffuseTransmission) {
+							if ((bsdf->getType() & BSDF::EAll) == BSDF::EDiffuseReflection || 
+								(bsdf->getType() & BSDF::EAll) == BSDF::EDiffuseTransmission) {
 								gatherPoint.weight = weight;
 								gatherPoint.depth = depth;
 								if (gatherPoint.its.isLuminaire())
@@ -319,7 +321,6 @@ public:
 		queue->signalRefresh(job, NULL);
 	}
 
-
 	std::string toString() const {
 		return "StochasticProgressivePhotonMapIntegrator[]";
 	}

src/libcore/thread.cpp

@@ -391,7 +391,7 @@ void Thread::initializeOpenMP(size_t threadCount) {
 				#endif
 				counter++;
 			}
-			thread->m_running = true;
+			thread->m_running = false;
 			thread->m_thread = pthread_self();
 			thread->m_joinMutex = new Mutex();
 			thread->m_joined = false;

src/librender/medium.cpp

@@ -67,15 +67,15 @@ void Medium::serialize(Stream *stream, InstanceManager *manager) const {
 
 std::string MediumSamplingRecord::toString() const {
 	std::ostringstream oss;
-	oss << "MediumSamplingRecord[" << std::endl
-		<< "  t = " << t << "," << std::endl
-		<< "  p = " << p.toString() << "," << std::endl
-		<< "  sigmaA = " << sigmaA.toString() << "," << std::endl
-		<< "  sigmaS = " << sigmaS.toString() << "," << std::endl
-		<< "  pdfFailure = " << pdfFailure << "," << std::endl
-		<< "  pdfSuccess = " << pdfSuccess << "," << std::endl
-		<< "  pdfSuccessRev = " << pdfSuccessRev << "," << std::endl
-		<< "  transmittance = " << transmittance.toString()
+	oss << "MediumSamplingRecord[" << endl
+		<< "  t = " << t << "," << endl
+		<< "  p = " << p.toString() << "," << endl
+		<< "  sigmaA = " << sigmaA.toString() << "," << endl
+		<< "  sigmaS = " << sigmaS.toString() << "," << endl
+		<< "  pdfFailure = " << pdfFailure << "," << endl
+		<< "  pdfSuccess = " << pdfSuccess << "," << endl
+		<< "  pdfSuccessRev = " << pdfSuccessRev << "," << endl
+		<< "  transmittance = " << transmittance.toString() << endl
 		<< "]";
 	return oss.str();
 }

src/librender/photon.cpp

@@ -62,7 +62,6 @@ Photon::Photon(const Point &p, const Normal &normal,
 			   uint16_t _depth) {
 	if (!P.isValid()) 
 		SLog(EWarn, "Creating an invalid photon with power: %s", P.toString().c_str());
-
 	/* Possibly convert to single precision floating point
 	   (if Mitsuba is configured to use double precision) */
 	position = p;

src/librender/photonmap.cpp

@@ -143,18 +143,25 @@ struct RawRadianceQuery {
 	inline void operator()(const Photon &photon) {
 		Normal photonNormal(photon.getNormal());
 		Vector wi = -photon.getDirection();
+		Float wiDotGeoN = absDot(photonNormal, wi);
 
 		if (photon.getDepth() > maxDepth
-			|| dot(photonNormal, its.shFrame.n) < 1e-1f 
-			|| dot(photonNormal, wi) < 1e-2f) /// XXX is this latter one really needed?
+			|| dot(photonNormal, its.shFrame.n) < 1e-1f
+			|| wiDotGeoN < 1e-2f)
 			return;
 
+		/* Prevent light leaks due to the use of shading normals -- [Veach, p. 158] */
 		BSDFQueryRecord bRec(its, its.toLocal(wi), its.wi, EImportance);
 
+		Spectrum value = photon.getPower() * bsdf->eval(bRec);
+		if (value.isZero())
+			return;
+
 		/* Account for non-symmetry due to shading normals */
-//		result += photon->getPower() * bsdf->eval(bRec) 
-//				/ dot(photonNormal, wiWorld);
-		result += photon.getPower() * bsdf->eval(bRec);
+		value *= Frame::cosTheta(bRec.wi) / 
+			(wiDotGeoN * Frame::cosTheta(bRec.wo));
+
+		result += value;
 	}
 
 	const Intersection &its;

src/medium/homogeneous.cpp

@@ -152,7 +152,7 @@ public:
 			for (int i=0; i<SPECTRUM_SAMPLES; ++i) {
 				/// Record the highest albedo values across channels
 				Float albedo = m_sigmaS[i] / m_sigmaT[i];
-				if (albedo > m_mediumSamplingWeight)
+				if (albedo > m_mediumSamplingWeight && m_sigmaT[i] != 0)
 					m_mediumSamplingWeight = albedo;
 			}
 			if (m_mediumSamplingWeight > 0) {
@@ -231,7 +231,11 @@ public:
 
 	void configure() {
 		Medium::configure();
-		m_albedo = (m_sigmaS/m_sigmaT).max();
+		m_albedo = 0;
+		for (size_t i=0; i<SPECTRUM_SAMPLES; ++i) {
+			if (m_sigmaT[i] != 0)
+				m_albedo = std::max(m_albedo, m_sigmaS[i]/m_sigmaT[i]);
+		}
 	}
 
 	void serialize(Stream *stream, InstanceManager *manager) const {
@@ -242,7 +246,12 @@ public:
 	}
 
 	Spectrum getTransmittance(const Ray &ray, Sampler *) const {
-		return (m_sigmaT * (ray.mint - ray.maxt)).exp();
+		Float negLength = ray.mint - ray.maxt;
+		Spectrum transmittance;
+		for (size_t i=0; i<SPECTRUM_SAMPLES; ++i)
+			transmittance[i] = m_sigmaT[i] != 0 
+				? std::exp(m_sigmaT[i] * negLength) : (Float) 1.0f;
+		return transmittance;
 	}
 
 	bool sampleDistance(const Ray &ray, MediumSamplingRecord &mRec,