jeremyong / klein

P(R*_{3, 0, 1}) specialized SIMD Geometric Algebra Library
https://jeremyong.com/klein
MIT License
749 stars 57 forks source link

Plane normalization e0 value #31

Open obliviand opened 1 month ago

obliviand commented 1 month ago

Hi, is there a reason that e0 value should be saved rather than normalized with everything else for a plane?

Running the test from issue #23 fails with the e0 having twice the expected value.

For the "plane-normalize" test below I checked the same plane with ganja.js and used those output values. The existing and new tests pass when the _mm_blend_ps and _mm_move_ss are removed.

diff --git a/public/klein/plane.hpp b/public/klein/plane.hpp
index 57a568a..14eef60 100644
--- a/public/klein/plane.hpp
+++ b/public/klein/plane.hpp
@@ -65,11 +65,6 @@ public:
     void normalize() noexcept
     {
         __m128 inv_norm = detail::rsqrt_nr1(detail::hi_dp_bc(p0_, p0_));
-#ifdef KLEIN_SSE_4_1
-        inv_norm = _mm_blend_ps(inv_norm, _mm_set_ss(1.f), 1);
-#else
-        inv_norm = _mm_move_ss(inv_norm, _mm_set_ss(1.f));
-#endif
         p0_ = _mm_mul_ps(inv_norm, p0_);
     }

diff --git a/test/test_metric.cpp b/test/test_metric.cpp
index 0687a4c..3ccc0b5 100644
--- a/test/test_metric.cpp
+++ b/test/test_metric.cpp
@@ -33,6 +33,42 @@ TEST_CASE("measure-point-to-plane")
     CHECK_EQ(std::abs((p1 ^ p2).e0123()), root_two);
 }

+TEST_CASE("plane-normalize")
+{
+    // d*e_0 + a*e_1 + b*e_2 + c*e_3
+    plane p{1.f, 2.f, 3.f, 4.f};
+
+    p.normalize();
+
+    CHECK_EQ(p.e0(), doctest::Approx(1.069f).epsilon(0.001));
+    CHECK_EQ(p.e1(), doctest::Approx(0.267f).epsilon(0.001));
+    CHECK_EQ(p.e2(), doctest::Approx(0.534f).epsilon(0.001));
+    CHECK_EQ(p.e3(), doctest::Approx(0.801f).epsilon(0.001));
+}
+
+TEST_CASE("normalized-plane-normalization")
+{
+    point a{-0.5, 2, -0.5};
+    point b{+0.5, 2, -0.5};
+    point c{+0.5, 2, +0.5};
+
+    // p is already normalized
+    plane p = a & b & c;
+    CHECK_EQ(p.norm(), 1);
+    CHECK_EQ(p.e0(), -2);
+    CHECK_EQ(p.e1(), 0);
+    CHECK_EQ(p.e2(), 1);
+    CHECK_EQ(p.e3(), 0);
+
+    p.normalize();
+    // same checks as before
+    CHECK_EQ(p.norm(), 1);
+    CHECK_EQ(p.e0(), -2);
+    CHECK_EQ(p.e1(), 0);
+    CHECK_EQ(p.e2(), 1);
+    CHECK_EQ(p.e3(), 0);
+}
+
 TEST_CASE("measure-point-to-line")
 {
     line l{0, 1, 0, 1, 0, 0};
jeremyong commented 1 month ago

Thanks feel free to submit a PR for that and I'll accept it!