fixed some inconsistencies in the HG phase function, clarified conventions in the documentation

include/mitsuba/render/medium.h

@@ -74,7 +74,10 @@ public:
 	Float miWeight;
 };
 
-/** \brief Abstract phase function
+/** \brief Abstract phase function.
+ *
+ * The convention used here is that the incident and exitant
+ * direction arguments point away from the scattering event (similar to BSDFs).
  */
 class MTS_EXPORT_RENDER PhaseFunction : public ConfigurableObject {
 public:
@@ -83,7 +86,7 @@ public:
 		EDelta
 	};
 
-	/// Evaluate the phase function for a pair of directions (wi, wo)
+	/// Evaluate the phase function for an outward-pointing pair of directions (wi, wo)
 	virtual Spectrum f(const MediumSamplingRecord &mRec, const Vector &wi, const Vector &wo) const = 0;
 
 	/**

src/medium/heterogeneous-flake.cpp

@@ -596,8 +596,8 @@ Spectrum FlakePhaseFunction::sample(const MediumSamplingRecord &mRec, const Vect
 
 	Vector wo = sphericalDirection(sample.x, sample.y);
 	_wo = frame.toWorld(wo);
-	return Spectrum(std::max((Float) 0, temp.eval(sample.x, sample.y)));
+//	return Spectrum(std::max((Float) 0, temp.eval(sample.x, sample.y)));
-//	return f(mRec, _wi, _wo);
+	return f(mRec, _wi, _wo);
 #else
 	/* Uniform sampling */
 	_wo = squareToSphere(_sample);

src/phase/hg.cpp

@@ -56,19 +56,18 @@ public:
 		if (std::abs(m_g) < Epsilon) {
 			cosTheta = 1 - 2*sample.x;
 		} else {
-			Float sqrTerm = (1 - m_g * m_g) /
+			Float sqrTerm = (1 - m_g * m_g) / (1 - m_g + 2 * m_g * sample.x);
-				(1 - m_g + 2 * m_g * sample.x);
 			cosTheta = (1 + m_g * m_g - sqrTerm * sqrTerm) / (2 * m_g);
 		}
 
-		Float sinTheta = std::sqrt(std::max((Float) 0, 1-cosTheta*cosTheta));
+		Float sinTheta = std::sqrt(std::max((Float) 0, 1.0f-cosTheta*cosTheta));
 		Float phi = 2*M_PI*sample.y, cosPhi = std::cos(phi), sinPhi = std::sin(phi);
 
 		Vector dir(
 			sinTheta * cosPhi,
 			sinTheta * sinPhi,
 			cosTheta);
-		wo = Frame(wi).toWorld(dir);
+		wo = Frame(-wi).toWorld(dir);
 		
 		sampledType = ENormal;
 		return Spectrum(1.0f);
@@ -84,7 +83,7 @@ public:
 
 	Spectrum f(const MediumSamplingRecord &mRec, const Vector &wi, const Vector &wo) const {
 		return Spectrum(1/(4*M_PI) * (1 - m_g*m_g) /
-			std::pow(1.f + m_g*m_g - 2.f * m_g * dot(wi, wo), (Float) 1.5f));
+			std::pow(1.f + m_g*m_g - 2.f * m_g * dot(-wi, wo), (Float) 1.5f));
 	}
 
 	std::string toString() const {

src/phase/kkay.cpp

@@ -111,7 +111,7 @@ public:
 			ESampledType &sampledType, const Point2 &sample) const {
 		wo = squareToSphere(sample);
 		sampledType = ENormal;
-		return Spectrum(1.0f);
+		return f(mRec, wi, wo) / (4 * M_PI);
 	}
 
 	Spectrum sample(const MediumSamplingRecord &mRec, const Vector &wi, Vector &wo, 
@@ -119,17 +119,17 @@ public:
 		wo = squareToSphere(sample);
 		sampledType = ENormal;
 		pdf = 1/(4 * (Float) M_PI);
-		return Spectrum(pdf);
+		return f(mRec, wi, wo);
 	}
 
 	Spectrum f(const MediumSamplingRecord &mRec, const Vector &wi, const Vector &wo) const {
-		Frame frame(mRec.orientation);
 		if (mRec.orientation.length() == 0)
 			return Spectrum(m_kd / (4*M_PI));
 
-		Vector orientation = normalize(mRec.orientation);
+		Frame frame(normalize(mRec.orientation));
 		Vector reflectedLocal = frame.toLocal(wo);
-		reflectedLocal.z = -dot(wi, orientation);
+
+		reflectedLocal.z = -dot(wi, frame.n);
 		Float a = std::sqrt((1-reflectedLocal.z*reflectedLocal.z) / 
 			(reflectedLocal.x*reflectedLocal.x + reflectedLocal.y*reflectedLocal.y));
 		reflectedLocal.y *= a;