fixed a subtle bug in the simpson integration method in heterogeneous.cpp; added support for querying the maximum value in density data volumes

include/mitsuba/render/volume.h

@@ -56,11 +56,19 @@ public:
 	virtual Vector lookupVector(const Point &p) const;
 
 	/**
-	 * Returns the recommended step size for numerical 
+	 * \brief Return the recommended step size for numerical 
 	 * integration or inifinity if this is not known/applicable
 	 */
 	virtual Float getStepSize() const = 0;
 
+	/**
+	 * \brief Return the maximum floating point value that
+	 * could be returned by \ref lookupFloat.
+	 *
+	 * This is useful when implementing Woodcock-Tracking.
+	 */
+	virtual Float getMaximumFloatValue() const = 0;
+
 	MTS_DECLARE_CLASS()
 protected:
 	/// Virtual destructor

src/medium/heterogeneous.cpp

@@ -133,7 +133,7 @@ public:
 
 		/* Assumes that the density medium does not 
 		   contain values greater than one! */
-		m_maxDensity = m_densityMultiplier;
+		m_maxDensity = m_densityMultiplier * m_density->getMaximumFloatValue();
 		if (m_anisotropicMedium)
 			m_maxDensity *= m_phaseFunction->sigmaDirMax();
 		m_invMaxDensity = 1.0f/m_maxDensity;
@@ -146,8 +146,9 @@ public:
 					m_orientation->getStepSize());
 
 			if (m_stepSize == std::numeric_limits<Float>::infinity()) 
-				Log(EError, "Unable to infer a suitable step size, please specify one "
+				Log(EError, "Unable to infer a suitable step size for deterministic "
-						"manually using the 'stepSize' parameter.");
+						"integration, please specify one manually using the 'stepSize' "
+						"parameter.");
 		}
 		
 		if (m_anisotropicMedium && m_orientation.get() == NULL)
@@ -216,6 +217,7 @@ public:
 
 		/* Compute a suitable step size */
 		uint32_t nSteps = (uint32_t) std::ceil(length / m_stepSize);
+		nSteps += nSteps % 2;
 		const Float stepSize = length/nSteps;
 		const Vector increment = ray.d * stepSize;
 
@@ -442,6 +444,9 @@ public:
 		if (m_method == ESimpsonQuadrature || sampler == NULL) {
 			return Spectrum(std::exp(-integrateDensity(ray)));
 		} else {
+			/* When Woodcock tracking is selected as the sampling method,
+			   we can use this method to get a noisy estimate of 
+			   the transmittance */
 			Float mint, maxt;
 			if (!m_densityAABB.rayIntersect(ray, mint, maxt))
 				return Spectrum(1.0f);
@@ -451,10 +456,10 @@ public:
 			#if defined(HETVOL_STATISTICS)
 				avgRayMarchingStepsTransmittance.incrementBase();
 			#endif
-			int nSteps = 2; /// XXX make configurable
+			int nSamples = 2; /// XXX make configurable
 			Float result = 0;
 
-			for (int i=0; i<nSteps; ++i) {
+			for (int i=0; i<nSamples; ++i) {
 				Float t = mint;
 				while (true) {
 					t -= std::log(1-sampler->next1D()) * m_invMaxDensity;
@@ -474,7 +479,7 @@ public:
 						break;
 				}
 			}
-			return Spectrum(result/nSteps);
+			return Spectrum(result/nSamples);
 		}
 	}
 

src/volume/constvolume.cpp

@@ -95,6 +95,10 @@ public:
 		return std::numeric_limits<Float>::infinity();
 	}
 
+	Float getMaximumFloatValue() const {
+		return m_float;
+	}
+
 	MTS_DECLARE_CLASS()
 protected:
 	int m_type;

src/volume/gridvolume.cpp

@@ -70,6 +70,13 @@ MTS_NAMESPACE_BEGIN
  *
  * Note that Mitsuba expects that entries in direction volumes are either
  * zero or valid unit vectors.
+ *
+ * When using this data source to represent floating point density volumes, 
+ * please ensure that the values are all normalized to lie in the 
+ * range [0, 1] -- otherwise, the Woocock-Tracking integration method in
+ * heterogeneous.cpp will produce incorrect results. You can use
+ * the 'densityMultiplier' parameter of that class to re-scale the
+ * densities if neccessary.
  */
 class GridDataSource : public VolumeDataSource {
 public:
@@ -546,6 +553,10 @@ public:
 	bool supportsVectorLookups() const { return m_channels == 3; }
 	Float getStepSize() const { return m_stepSize; }
 
+	Float getMaximumFloatValue() const {
+		return 1.0f;
+	}
+
 	MTS_DECLARE_CLASS()
 protected: 
 	FINLINE Vector lookupQuantizedDirection(size_t index) const {

src/volume/hgridvolume.cpp

@@ -108,6 +108,7 @@ public:
 					createObject(MTS_CLASS(VolumeDataSource), props));
 			content->configure();
 
+			m_maxFloatValue = contents->getMaximumFloatValue();
 			m_blocks[(m_res.y * block.z + block.y) * m_res.x + block.x] = content;
 			m_stepSize = std::min(m_stepSize, content->getStepSize());
 			m_supportsVectorLookups = m_supportsVectorLookups && content->supportsVectorLookups();
@@ -191,6 +192,9 @@ public:
 			return block->lookupVector(_p);
 	}
 
+	Float getMaximumFloatValue() const {
+		return m_maxFloatValue;
+	}
 
 	MTS_DECLARE_CLASS()
 protected:
@@ -204,7 +208,7 @@ protected:
 	bool m_supportsFloatLookups;
 	bool m_supportsSpectrumLookups;
 	bool m_supportsVectorLookups;
-	Float m_stepSize;
+	Float m_stepSize, m_maxFloatValue;
 };
 
 MTS_IMPLEMENT_CLASS_S(HierarchicalGridDataSource, false, VolumeDataSource);

src/volume/volcache.cpp

@@ -288,6 +288,10 @@ public:
 		delete[] ptr;
 	}
 
+	Float getMaximumFloatValue() const {
+		return m_nested->getMaximumFloatValue();
+	}
+
 	MTS_DECLARE_CLASS()
 protected:
 	ref<VolumeDataSource> m_nested;