fixed HG issues in issue #455

mitsuba-renderer · Dec 18, 2022 · 10d3514 · 10d3514
1 parent 176337c
commit 10d3514
Show file tree

Hide file tree

Showing 3 changed files with 26 additions and 19 deletions.
diff --git a/README.md b/README.md
@@ -134,11 +134,11 @@ Significant features and/or improvements to the code were contributed by
 When using Mitsuba 3 in academic projects, please cite:
 
 ```bibtex
-@software{jakob2022mitsuba3,
+@software{Mitsuba3,
     title = {Mitsuba 3 renderer},
     author = {Wenzel Jakob and Sébastien Speierer and Nicolas Roussel and Merlin Nimier-David and Delio Vicini and Tizian Zeltner and Baptiste Nicolet and Miguel Crespo and Vincent Leroy and Ziyi Zhang},
     note = {https://mitsuba-renderer.org},
-    version = {3.0.1},
-    year = 2022,
+    version = {3.1.1},
+    year = 2022
 }
 ```
diff --git a/src/phase/hg.cpp b/src/phase/hg.cpp
@@ -49,7 +49,7 @@ class HGPhaseFunction final : public PhaseFunction<Float, Spectrum> {
 
     HGPhaseFunction(const Properties &props) : Base(props) {
         ScalarFloat g = props.get<ScalarFloat>("g", 0.8f);
-        if (g >= 1 || g <= -1)
+        if (g >= 1.f || g <= -1.f)
             Log(Error, "The asymmetry parameter must lie in the interval (-1, 1)!");
         m_g = g;
 
@@ -59,12 +59,14 @@ class HGPhaseFunction final : public PhaseFunction<Float, Spectrum> {
     }
 
     void traverse(TraversalCallback *callback) override {
-        callback->put_parameter("g", m_g, ParamFlags::Differentiable | ParamFlags::Discontinuous);
+        callback->put_parameter("g", m_g, ParamFlags::Differentiable |
+                                          ParamFlags::Discontinuous);
     }
 
     MI_INLINE Float eval_hg(Float cos_theta) const {
-        Float temp = 1.0f + m_g * m_g + 2.0f * m_g * cos_theta;
-        return dr::InvFourPi<ScalarFloat> * (1 - m_g * m_g) / (temp * dr::sqrt(temp));
+        Float temp = 1.f + dr::sqr(m_g) + 2.f * m_g * cos_theta;
+        return dr::InvFourPi<ScalarFloat> * (1.f - dr::sqr(m_g)) /
+               (temp * dr::sqrt(temp));
     }
 
     std::pair<Vector3f, Float> sample(const PhaseFunctionContext & /* ctx */,
@@ -74,16 +76,19 @@ class HGPhaseFunction final : public PhaseFunction<Float, Spectrum> {
                                       Mask active) const override {
         MI_MASKED_FUNCTION(ProfilerPhase::PhaseFunctionSample, active);
 
-        Float sqr_term = (1 - m_g * m_g) / (1 - m_g + 2 * m_g * sample2.x());
-        Float cos_theta = (1 + m_g * m_g - sqr_term * sqr_term) / (2 * m_g);
-        dr::masked(cos_theta, m_g < dr::Epsilon<ScalarFloat>) = 1 - 2 * sample2.x();
+        Float sqr_term  = (1.f - dr::sqr(m_g)) / (1.f - m_g + 2.f * m_g * sample2.x()),
+              cos_theta = (1.f + dr::sqr(m_g) - dr::sqr(sqr_term)) / (2.f * m_g);
 
-        Float sin_theta = dr::safe_sqrt(1.0f - cos_theta * cos_theta);
-        auto [sin_phi, cos_phi] = dr::sincos(2 * dr::Pi<ScalarFloat> * sample2.y());
-        auto wo = Vector3f(sin_theta * cos_phi, sin_theta * sin_phi, -cos_theta);
-        wo = mi.to_world(wo);
-        Float pdf = eval_hg(-cos_theta);
-        return { wo, pdf };
+        // Diffuse fallback
+        dr::masked(cos_theta, dr::abs(m_g) < dr::Epsilon<ScalarFloat>) = 1.f - 2.f * sample2.x();
+
+        Float sin_theta = dr::safe_sqrt(1.f - dr::sqr(cos_theta));
+        auto [sin_phi, cos_phi] = dr::sincos(2.f * dr::Pi<ScalarFloat> * sample2.y());
+
+        Vector3f wo = mi.to_world(
+            Vector3f(sin_theta * cos_phi, sin_theta * sin_phi, -cos_theta));
+
+        return { wo, eval_hg(-cos_theta) };
     }
 
     Float eval(const PhaseFunctionContext & /* ctx */, const MediumInteraction3f &mi,
@@ -103,7 +108,6 @@ class HGPhaseFunction final : public PhaseFunction<Float, Spectrum> {
     MI_DECLARE_CLASS()
 private:
     Float m_g;
-
 };
 
 MI_IMPLEMENT_CLASS_VARIANT(HGPhaseFunction, PhaseFunction)

diff --git a/src/phase/tests/test_hg.py b/src/phase/tests/test_hg.py
@@ -1,12 +1,15 @@
 import mitsuba as mi
+import pytest
 
 
 def test01_create(variant_scalar_rgb):
     p = mi.load_dict({"type": "hg", "g": 0.4})
     assert p is not None
 
-def test02_chi2(variants_vec_backends_once_rgb):
-    sample_func, pdf_func = mi.chi2.PhaseFunctionAdapter("hg", '<float name="g" value="0.6"/>')
+
+@pytest.mark.parametrize('g', ['0.6', '-0.6'])
+def test02_chi2(variants_vec_backends_once_rgb, g):
+    sample_func, pdf_func = mi.chi2.PhaseFunctionAdapter("hg", f'<float name="g" value="{g}"/>')
 
     chi2 = mi.chi2.ChiSquareTest(
         domain=mi.chi2.SphericalDomain(),