Regularize quadric surfaces to deduplicate opposing surfaces (#1237)

* Add exterior bounding box for gentrap
* Add example showing failure to deduplicate
* Un-inline full GQ/SQ constructors
* Regularize signs on quadrics
* Fix function decorator and add comments with regularization test

scripts/cmake-presets/goldfinger.json

@@ -24,7 +24,7 @@
         "CMAKE_CXX_STANDARD":   {"type": "STRING",   "value": "17"},
         "CMAKE_CXX_EXTENSIONS": {"type": "BOOL",   "value": "OFF"},
         "CMAKE_FIND_FRAMEWORK": {"type": "STRING", "value": "LAST"},
-        "CMAKE_INSTALL_PREFIX": "${sourceDir}/install-${presetName}"
+        "CMAKE_INSTALL_PREFIX": "${sourceDir}/install-${presetName}",
         "CMAKE_CXX_FLAGS_RELWITHDEBINFO": "-O2 -g -DNDEBUG -fno-inline -fno-omit-frame-pointer",
         "CMAKE_CXX_FLAGS": "-Wall -Wextra -Werror -Wno-error=deprecated -pedantic -fdiagnostics-color=always"
       }

src/orange/surf/GeneralQuadric.cc

@@ -7,10 +7,39 @@
 //---------------------------------------------------------------------------//
 #include "GeneralQuadric.hh"
 
+#include "corecel/Assert.hh"
+
 #include "SimpleQuadric.hh"
 
 namespace celeritas
 {
+//---------------------------------------------------------------------------//
+/*!
+ * Construct with all coefficients.
+ *
+ * TODO: normalize so that largest eigenvalue is unity? Or what? (It would be
+ * nice to have "slightly twisted planes" have order-epsilon cross terms as
+ * opposed to order 1/eps linear terms.)
+ */
+GeneralQuadric::GeneralQuadric(Real3 const& abc,
+                               Real3 const& def,
+                               Real3 const& ghi,
+                               real_type j)
+    : a_(abc[0])
+    , b_(abc[1])
+    , c_(abc[2])
+    , d_(def[0])
+    , e_(def[1])
+    , f_(def[2])
+    , g_(ghi[0])
+    , h_(ghi[1])
+    , i_(ghi[2])
+    , j_(j)
+{
+    CELER_EXPECT(a_ != 0 || b_ != 0 || c_ != 0 || d_ != 0 || e_ != 0 || f_ != 0
+                 || g_ != 0 || h_ != 0 || i_ != 0);
+}
+
 //---------------------------------------------------------------------------//
 /*!
  * Promote from a simple quadric.

src/orange/surf/GeneralQuadric.hh

@@ -59,11 +59,11 @@ class GeneralQuadric
   public:
     //// CONSTRUCTORS ////
 
-    // Construct with radius
-    explicit inline CELER_FUNCTION GeneralQuadric(Real3 const& abc,
-                                                  Real3 const& def,
-                                                  Real3 const& ghi,
-                                                  real_type j);
+    // Construct from coefficients
+    explicit GeneralQuadric(Real3 const& abc,
+                            Real3 const& def,
+                            Real3 const& ghi,
+                            real_type j);
 
     // Construct from raw data
     template<class R>
@@ -114,31 +114,6 @@ class GeneralQuadric
 
 //---------------------------------------------------------------------------//
 // INLINE DEFINITIONS
-//---------------------------------------------------------------------------//
-/*!
- * Construct with all coefficients.
- *
- * TODO: normalize?
- */
-CELER_FUNCTION GeneralQuadric::GeneralQuadric(Real3 const& abc,
-                                              Real3 const& def,
-                                              Real3 const& ghi,
-                                              real_type j)
-    : a_(abc[0])
-    , b_(abc[1])
-    , c_(abc[2])
-    , d_(def[0])
-    , e_(def[1])
-    , f_(def[2])
-    , g_(ghi[0])
-    , h_(ghi[1])
-    , i_(ghi[2])
-    , j_(j)
-{
-    CELER_EXPECT(a_ != 0 || b_ != 0 || c_ != 0 || d_ != 0 || e_ != 0 || f_ != 0
-                 || g_ != 0 || h_ != 0 || i_ != 0);
-}
-
 //---------------------------------------------------------------------------//
 /*!
  * Construct from raw data.

src/orange/surf/SimpleQuadric.cc

@@ -7,6 +7,7 @@
 //---------------------------------------------------------------------------//
 #include "SimpleQuadric.hh"
 
+#include "corecel/Assert.hh"
 #include "corecel/cont/Range.hh"
 #include "corecel/math/Algorithms.hh"
 #include "corecel/math/ArrayOperators.hh"
@@ -18,6 +19,27 @@
 
 namespace celeritas
 {
+//---------------------------------------------------------------------------//
+/*!
+ * Construct with coefficients.
+ *
+ * The quadric is ill-defined if all non-constants are zero.
+ *
+ * TODO: normalize so that largest eigenvalue is unity?
+ */
+SimpleQuadric::SimpleQuadric(Real3 const& abc, Real3 const& def, real_type g)
+    : a_(abc[0])
+    , b_(abc[1])
+    , c_(abc[2])
+    , d_(def[0])
+    , e_(def[1])
+    , f_(def[2])
+    , g_(g)
+{
+    CELER_EXPECT(a_ != 0 || b_ != 0 || c_ != 0 || d_ != 0 || e_ != 0
+                 || f_ != 0);
+}
+
 //---------------------------------------------------------------------------//
 /*!
  * Promote from a plane.

src/orange/surf/SimpleQuadric.hh

@@ -62,7 +62,6 @@ class SimpleQuadric
     //// CONSTRUCTORS ////
 
     // Construct from coefficients
-    inline CELER_FUNCTION
     SimpleQuadric(Real3 const& abc, Real3 const& def, real_type g);
 
     // Construct from raw data
@@ -120,28 +119,6 @@ class SimpleQuadric
 
 //---------------------------------------------------------------------------//
 // INLINE DEFINITIONS
-//---------------------------------------------------------------------------//
-/*!
- * Construct with coefficients.
- *
- * The quadric is ill-defined if all non-constants are zero.
- *
- * TODO: normalize?
- */
-CELER_FUNCTION
-SimpleQuadric::SimpleQuadric(Real3 const& abc, Real3 const& def, real_type g)
-    : a_(abc[0])
-    , b_(abc[1])
-    , c_(abc[2])
-    , d_(def[0])
-    , e_(def[1])
-    , f_(def[2])
-    , g_(g)
-{
-    CELER_EXPECT(a_ != 0 || b_ != 0 || c_ != 0 || d_ != 0 || e_ != 0
-                 || f_ != 0);
-}
-
 //---------------------------------------------------------------------------//
 /*!
  * Construct from raw data.

src/orange/surf/SurfaceSimplifier.cc

@@ -58,6 +58,40 @@ class ZeroSnapper
     SoftZero<> soft_zero_;
 };
 
+struct SignCount
+{
+    int pos{0};
+    int neg{0};
+};
+
+SignCount count_signs(Span<real_type const, 3> arr, real_type tol)
+{
+    SignCount result;
+    for (auto v : arr)
+    {
+        if (v < -tol)
+            ++result.neg;
+        else if (v > tol)
+            ++result.pos;
+    }
+    return result;
+}
+
+template<class S>
+S negate_coefficients(S const& orig)
+{
+    auto arr = make_array(orig.data());
+    for (real_type& v : arr)
+    {
+        v = negate(v);
+    }
+
+    // TODO: make_span doesn't use the correct overload and creates a
+    // dynamic extent span
+    using SpanT = decltype(orig.data());
+    return S{SpanT{arr.data(), arr.size()}};
+}
+
 //---------------------------------------------------------------------------//
 }  // namespace
 
@@ -236,9 +270,6 @@ auto SurfaceSimplifier::operator()(Sphere const& s) const
 //---------------------------------------------------------------------------//
 /*!
  * Simple quadric with near-zero terms can be another second-order surface.
- *
- * TODO: renormalize so that second-order terms are O(1) (and simplifying
- * quadrics that are scaled by a constant)?
  */
 auto SurfaceSimplifier::operator()(SimpleQuadric const& sq)
     -> Optional<Plane,
@@ -252,48 +283,28 @@ auto SurfaceSimplifier::operator()(SimpleQuadric const& sq)
                 SimpleQuadric>
 {
     // Determine possible simplifications by calculating number of zeros
-    int num_pos{0};
-    int num_neg{0};
-    for (auto v : sq.second())
-    {
-        if (v < -tol_)
-            ++num_neg;
-        else if (v > tol_)
-            ++num_pos;
-    }
+    auto signs = count_signs(sq.second(), tol_);
 
-    if (num_pos == 0 && num_neg == 0)
+    if (signs.pos == 0 && signs.neg == 0)
     {
         // It's a plane
         return detail::QuadricPlaneConverter{tol_}(sq);
     }
-    else if (num_neg > num_pos)
+    else if (signs.neg > signs.pos)
     {
-        // Normalize sign so that it has more positive signs than negative
-        auto arr = make_array(sq.data());
-        for (auto& v : arr)
-        {
-            v = negate(v);
-        }
-
-        // Flip sense
+        // More positive signs than negative:
+        // flip the sense and reverse the values
         *sense_ = flip_sense(*sense_);
-
-        // Construct reversed-sign quadric
-        // Todo: make_span doesn't use the correct overload and creates a
-        // dynamic extent span
-        return SimpleQuadric{
-            Span<decltype(arr)::value_type, SimpleQuadric::StorageSpan::extent>{
-                make_span(arr)}};
+        return negate_coefficients(sq);
     }
-    else if (num_pos == 3)
+    else if (signs.pos == 3)
     {
         // Could be a sphere
         detail::QuadricSphereConverter to_sphere{tol_};
         if (auto s = to_sphere(sq))
             return *s;
     }
-    else if (num_pos == 2 && num_neg == 1)
+    else if (signs.pos == 2 && signs.neg == 1)
     {
         // Cone: one second-order term less than zero, others equal
         detail::QuadricConeConverter to_cone{tol_};
@@ -304,7 +315,7 @@ auto SurfaceSimplifier::operator()(SimpleQuadric const& sq)
         if (auto c = to_cone(AxisTag<Axis::z>{}, sq))
             return *c;
     }
-    else if (num_pos == 2 && num_neg == 0)
+    else if (signs.pos == 2 && signs.neg == 0)
     {
         // Cyl: one second-order term is zero, others are equal
         detail::QuadricCylConverter to_cyl{tol_};
@@ -322,21 +333,35 @@ auto SurfaceSimplifier::operator()(SimpleQuadric const& sq)
 
 //---------------------------------------------------------------------------//
 /*!
- * Quadric with no cross terms is "simple".
+ * Quadric can be regularized or simplified.
  *
- * TODO: guard against different-signed GQs?
+ * - When no cross terms are present, it's "simple".
+ * - When the higher-order terms are negative, the signs need to be flipped.
  */
 auto SurfaceSimplifier::operator()(GeneralQuadric const& gq)
-    -> Optional<SimpleQuadric>
+    -> Optional<SimpleQuadric, GeneralQuadric>
 {
-    auto cross = gq.cross();
-    if (std::all_of(cross.begin(), cross.end(), SoftZero{tol_}))
+    // Cross term signs
+    auto csigns = count_signs(gq.cross(), tol_);
+    if (csigns.pos == 0 && csigns.neg == 0)
     {
         // No cross terms
         return SimpleQuadric{
             make_array(gq.second()), make_array(gq.first()), gq.zeroth()};
     }
 
+    // Second-order term signs
+    auto ssigns = count_signs(gq.second(), tol_);
+    if (ssigns.neg > ssigns.pos
+        || (ssigns.neg == 0 && ssigns.pos == 0 && csigns.neg > csigns.pos))
+    {
+        // More positive signs than negative:
+        // flip the sense and reverse the values
+        *sense_ = flip_sense(*sense_);
+        return negate_coefficients(gq);
+    }
+
+    // No simplification
     return {};
 }
 

src/orange/surf/SurfaceSimplifier.hh

@@ -79,8 +79,8 @@ class SurfaceSimplifier
              SimpleQuadric>
     operator()(SimpleQuadric const&);
 
-    // Quadric with no cross terms is simple
-    Optional<SimpleQuadric> operator()(GeneralQuadric const&);
+    // Quadric can be normalized or simplified
+    Optional<SimpleQuadric, GeneralQuadric> operator()(GeneralQuadric const&);
 
     //! Default: no simplifcation
     template<class S>