ENH: SampleSurfaceMesh and GeometryMath Speed Optimizations (#1020)

* Add caching routine - approx 15% speedup

* Clear excessive branching and excess copying

* extract costly construction in tight loop

* Consolidate if statement calculations, strip expensive modulo with bitwise calculation, reenable randomness

* strip extraneous coping, further if cleanup

* Initial Ray changes to attempt to preserve pure c++, converted 'DoesRayIntersectBox` if statements

* Add a caching version of Ray to prevent repeated calculations for the same result in loops. Selected to preserve existing API

* extract if statements in favor of integer math for CPU pipeline speed, parameter reordering for clearer distinction in overloads

* clean up remaining ifs and reduce ray misses by appropriately sizing ray

* patch bounding box bug in branchless calculations

* repair overlooked functionality in Ray (cached version)

* Default seed cleans up output and speeds up process

* Modified Ray class

* walk back branchless changes -- processing time is roughly equivalent between branched and branchless, and branched is more readable

* Add type declarations for Clang Compiler

src/Plugins/SimplnxCore/docs/RegularGridSampleSurfaceMesh.md

@@ -10,9 +10,55 @@ This **Filter** "samples" a triangulated surface mesh on a rectilinear grid. The
 
 1. Determine the bounding box and **Triangle** list of each **Feature** by scanning all **Triangles** and noting the **Features** on either side of the **Triangle**
 2. For each **Cell** in the rectilinear grid, determine which bounding box(es) they fall in (*Note:* the bounding box of multiple **Features** can overlap)
-3. For each bounding box a **Cell** falls in, check against that **Feature's** **Triangle** list to determine if the **Cell** falls within that n-sided polyhedra (*Note:* if the surface mesh is conformal, then each **Cell** will only belong to one **Feature**, but if not, the last **Feature** the **Cell** is found to fall inside of will *own* the **Cell**)
+3. For each bounding box a **Cell** falls in, check against that **Feature's** **Triangle** list to determine if the **Cell** falls within that n-sided polyhedra (*Note:* if the surface mesh is conformal, then each **Cell** will only belong to one **Feature**, but if not, the last **Feature** the **Cell** is found to fall inside will *own* the **Cell**)
 4. Assign the **Feature** number that the **Cell** falls within to the *Feature Ids* array in the new rectilinear grid geometry
 
+### Origin, Dimension, and Spacing's Effect on the Output
+
+***Subsequent images in this section will always have the correct example on the left***. The images below were created by thresholding the `Feature Ids` to `1`.
+
+#### Origin
+
+Take a moment to consider the following image:
+
+![Default and Offset Origin](Images/Regular_Grid_Mesh_Offset.png)
+
+In the above picture, the cylinder on the left is a correct output for a specific `.stl` model, the origin used was (-1, -1, -1). On the right, there is a cylinder that had the origin set to (3, 2, -1). The change of origin caused a truncating of the geometry, the white outline is roughly the bounding box for the objects. Thus, when you are defining the origin it is important to consider how the resulting box will intersect the geometry. *For most cases, it is recommended to ensure the origin of the box falls below that of the Geometry to avoid truncation.*
+
+#### Dimension
+
+Take a moment to consider the following image:
+
+![Default and Undersized Voxel](Images/Regular_Grid_Mesh_Undersized.png)
+
+In the above picture, the cylinder on the left is a correct output for a specific `.stl` model, the dimensions used was (440, 440, 640). On the right, there is a cylinder that had the dimensions set to (240, 240, 240). The change in dimensions led to a large truncation of the overall geometry. This parameter is the most difficult to fine tune since setting the dimensions too large can easily blow out your RAM needlessly, while too small will truncate sections potentially leading to output looking deformed (such as a rectangle appearing as a square).
+
+#### Scaling
+
+Take a moment to consider the following image:
+
+![Default and Rescaled](Images/Regular_Grid_Mesh_Rescaled.png)
+
+In the above picture, the cylinder on the left is a correct output for a specific `.stl` model, the scaling used was (0.05, 0.05, 0.05) with dimensions of (440, 440, 640). On the right, there is a cylinder that had the scaling set to (0.15, 0.15, 0.15) with dimensions of (240, 240, 240). You will notice there is no truncation to the geometry despite the smaller box dimensions, thus scaling can be used to reduce the size of the data. **However, there is a caveat, less points means less detail**. The astute may have noticed the cylinder on the right looks more "pixelated" despite having no truncation. Here's a close up image of the edge for easier viewing:
+
+![Default and Rescaled Edge](Images/Regular_Grid_Mesh_Edge.png)
+
+#### Tips
+
+It's easiest to think about these parameters as if you were defining a 3D image to capture. The Dimensions are the pixels (voxels) in your image. They allow you to define the detail (resolution) of the resulting image:
+
+- lower values will save space and run faster, but be pixelated
+- higher values will be detailed, but consume more space and take longer to process.
+
+The spacing allows you to scale the size of the bounding box, which can be thought of as similar to zooming a camera in/out from your focal point.
+
+- larger scaling values can make the capture (bounding box/grid) too large, and your focal point will become pixelated because its "too small" in the overall image. Think about zooming in super close to a picture.
+- smaller scaling values are bring the focal point closer, ensuring more detail is captured. However, if you get too close pieces of the focal point will be missing from the image, but what you do capture will have lots of detail.
+
+The origin can be thought of as moving the camera left and right and up and down. Allowing you to center the focal point in the image or cut out unwanted parts.
+
+All in all, it is best practice to try to fit the bounding box as closely to the actual geometry as possible to ensure you get the most detail out of the used memory.
+
 % Auto generated parameter table will be inserted here
 
 ## License & Copyright

src/simplnx/Common/Ray.hpp

@@ -1,13 +1,13 @@
 #pragma once
 
 #include <stdexcept>
+#include <utility>
 
 #include "simplnx/Common/Array.hpp"
 #include "simplnx/Common/EulerAngle.hpp"
 
 namespace nx::core
 {
-
 /**
  * @class Ray
  * @brief The Ray class describes a ray or line segment in 3D space. Using
@@ -41,9 +41,9 @@ class Ray
    * @param ang
    * @param length
    */
-  Ray(const PointType& origin, const ZXZEulerType& ang, LengthType length)
+  Ray(const PointType& origin, ZXZEulerType ang, LengthType length)
   : m_Origin(origin)
-  , m_Angle(ang)
+  , m_Angle(std::move(ang))
   , m_Length(length)
   {
   }
@@ -81,6 +81,15 @@ class Ray
     return m_Origin;
   }
 
+  /**
+   * @brief Returns the origin point.
+   * @return const PointType&
+   */
+  const PointType& getOriginRef() const
+  {
+    return m_Origin;
+  }
+
   /**
    * @brief Returns the Euler angle describing the Ray's direction.
    * @return ZXZEulerType
@@ -94,7 +103,7 @@ class Ray
    * @brief Sets a new origin point.
    * @param origin
    */
-  void setOrigin(const PointType& origin)
+  virtual void setOrigin(const PointType& origin)
   {
     m_Origin = origin;
   }
@@ -103,7 +112,7 @@ class Ray
    * @brief Sets a new Euler angle.
    * @param ang
    */
-  void setEuler(const ZXZEulerType& ang)
+  virtual void setEuler(const ZXZEulerType& ang)
   {
     m_Angle = ang;
   }
@@ -121,7 +130,7 @@ class Ray
    * @brief Sets a new ray length.
    * @param length
    */
-  void setLength(LengthType length)
+  virtual void setLength(LengthType length)
   {
     m_Length = length;
   }
@@ -131,22 +140,9 @@ class Ray
    * Based on the assumption the point is initially aligned with the global axis' and the that this vector specifically aligned with the local y-axis.
    * @return PointType
    */
-  PointType getEndPoint() const
+  virtual PointType getEndPoint() const
   {
-    const auto sin1 = std::sin(m_Angle[0]);
-    const auto sin2 = std::sin(m_Angle[1]);
-    const auto sin3 = std::sin(m_Angle[2]);
-
-    const auto cos1 = std::cos(m_Angle[0]);
-    const auto cos2 = std::cos(m_Angle[1]);
-    const auto cos3 = std::cos(m_Angle[2]);
-
-    // Reference: https://ntrs.nasa.gov/api/citations/19770019231/downloads/19770019231.pdf Page:23
-    Vec3<T> localXRotationVec((-sin1 * cos2 * sin3) + (cos1 * cos3), (cos1 * cos2 * sin3) + (sin1 * cos3), sin2 * sin3);
-    Vec3<T> localYRotationVec((-sin1 * cos2 * cos3) - (cos1 * cos3), (cos1 * cos2 * cos3) - (sin1 * sin3), sin2 * cos3);
-    Vec3<T> localZRotationVec((sin1 * sin2), -cos1 * sin2, cos2);
-
-    return m_Origin + (localXRotationVec * m_Length) + (localYRotationVec * m_Length) + (localZRotationVec * m_Length);
+    return calculateEndpoint();
   }
 
   /**
@@ -212,9 +208,125 @@ class Ray
   }
 
 protected:
-private:
   PointType m_Origin;
   ZXZEulerType m_Angle;
   LengthType m_Length;
+
+  PointType calculateEndpoint() const
+  {
+    const auto sin1 = std::sin(m_Angle[0]);
+    const auto sin2 = std::sin(m_Angle[1]);
+    const auto sin3 = std::sin(m_Angle[2]);
+
+    const auto cos1 = std::cos(m_Angle[0]);
+    const auto cos2 = std::cos(m_Angle[1]);
+    const auto cos3 = std::cos(m_Angle[2]);
+
+    // Reference: https://ntrs.nasa.gov/api/citations/19770019231/downloads/19770019231.pdf Page:23
+    Vec3<T> localXRotationVec((-sin1 * cos2 * sin3) + (cos1 * cos3), (cos1 * cos2 * sin3) + (sin1 * cos3), sin2 * sin3);
+    Vec3<T> localYRotationVec((-sin1 * cos2 * cos3) - (cos1 * cos3), (cos1 * cos2 * cos3) - (sin1 * sin3), sin2 * cos3);
+    Vec3<T> localZRotationVec((sin1 * sin2), -cos1 * sin2, cos2);
+
+    return m_Origin + (localXRotationVec * m_Length) + (localYRotationVec * m_Length) + (localZRotationVec * m_Length);
+  }
+};
+
+template <typename T>
+class CachedRay : public Ray<T>
+{
+public:
+  /**
+   * @brief Creates a Ray from an origin point, ZXZEuler angle, and target length.
+   * @param origin
+   * @param ang
+   * @param length
+   */
+  CachedRay(const typename Ray<T>::PointType& origin, typename Ray<T>::ZXZEulerType ang, typename Ray<T>::LengthType length)
+  : Ray<T>(origin, std::move(ang), length)
+  {
+    m_Endpoint = this->calculateEndpoint();
+  }
+
+  /**
+   * @brief Copy constructor
+   * @param other
+   */
+  CachedRay(const CachedRay& other)
+  : m_Endpoint(other.m_Endpoint)
+  , Ray<T>(other)
+  {
+  }
+
+  /**
+   * @brief Move constructor
+   * @param other
+   */
+  CachedRay(CachedRay&& other) noexcept
+  : m_Endpoint(std::move(other.m_Endpoint))
+  , Ray<T>(std::move(other))
+  {
+  }
+
+  ~CachedRay() = default;
+
+  /**
+   * @brief Copy assignment operator
+   * @param rhs
+   * @return CachedRay&
+   */
+  CachedRay& operator=(const CachedRay& rhs)
+  {
+    this->m_Origin = rhs.m_Origin;
+    this->m_Angle = rhs.m_Angle;
+    this->m_Length = rhs.m_Length;
+    m_Endpoint = rhs.m_Endpoint;
+    return *this;
+  }
+
+  /**
+   * @brief Move assignment operator
+   * @param rhs
+   * @return Ray&
+   */
+  CachedRay& operator=(CachedRay&& rhs) noexcept
+  {
+    this->m_Origin = std::move(rhs.m_Origin);
+    this->m_Angle = std::move(rhs.m_Angle);
+    this->m_Length = std::move(rhs.m_Length);
+    m_Endpoint = std::move(rhs.m_Endpoint);
+    return *this;
+  }
+
+  void setOrigin(const typename Ray<T>::PointType& origin) override
+  {
+    this->m_Origin = origin;
+    m_Endpoint = this->calculateEndpoint();
+  }
+
+  void setLength(typename Ray<T>::LengthType length) override
+  {
+    this->m_Length = length;
+    m_Endpoint = this->calculateEndpoint();
+  }
+
+  void setEuler(const typename Ray<T>::ZXZEulerType& ang) override
+  {
+    this->m_Angle = ang;
+    m_Endpoint = this->calculateEndpoint();
+  }
+
+  const typename Ray<T>::PointType& getEndPointRef() const
+  {
+    return m_Endpoint;
+  }
+
+  typename Ray<T>::PointType getEndPoint() const override
+  {
+    return m_Endpoint;
+  }
+
+private:
+  // Optional caching for speed
+  typename Ray<T>::PointType m_Endpoint = {};
 };
 } // namespace nx::core

src/simplnx/Utilities/Math/GeometryMath.cpp

@@ -54,10 +54,17 @@ BoundingBox3Df nx::core::GeometryMath::FindBoundingBoxOfVertices(INodeGeometry0D
   return {ll, ur}; // should be valid
 }
 
-BoundingBox3Df nx::core::GeometryMath::FindBoundingBoxOfFace(const TriangleGeom& faces, int32 faceId)
+BoundingBox3Df nx::core::GeometryMath::FindBoundingBoxOfFace(const detail::GeometryStoreCache& cache, const nx::core::TriangleGeom& triangleGeom, int32 faceId, std::vector<usize>& verts)
 {
-  std::array<Point3Df, 3> points;
-  faces.getFaceCoordinates(faceId, points);
+  // Unavoidable branch to verify integrity of subsequent function
+  // specifically because this an exposed function, however, the
+  // CPU should quickly be able to predict once it has run a few
+  // iterations of the function making it relatively lightweight
+  if(verts.size() != cache.NumVertsPerFace)
+  {
+    verts = std::vector<usize>(cache.NumVertsPerFace);
+  }
+  std::array<Point3Df, 3> points = GeometryMath::detail::GetFaceCoordinates<float32>(cache, faceId, verts);
 
   Point3Df ll = points[0];
   Point3Df ur = points[0];
@@ -114,7 +121,7 @@ BoundingBox3Df nx::core::GeometryMath::FindBoundingBoxOfFace(const TriangleGeom&
   return {ll, ur};
 }
 
-BoundingBox3Df nx::core::GeometryMath::FindBoundingBoxOfFaces(const TriangleGeom& faces, const std::vector<int32>& faceIds)
+BoundingBox3Df nx::core::GeometryMath::FindBoundingBoxOfFaces(const nx::core::TriangleGeom& triangleGeom, const std::vector<int32>& faceIds)
 {
   Point3Df ll(0, 0, 0);
   Point3Df ur(0, 0, 0);
@@ -124,9 +131,14 @@ BoundingBox3Df nx::core::GeometryMath::FindBoundingBoxOfFaces(const TriangleGeom
     return {ll, ur};
   }
 
+  detail::GeometryStoreCache cache(triangleGeom.getVertices()->getDataStoreRef(), triangleGeom.getFaces()->getDataStoreRef(), triangleGeom.getNumberOfVerticesPerFace());
+
+  // initialize temp storage 'verts' vector to avoid expensive
+  // calls during tight loops below
+  std::vector<usize> verts(cache.NumVertsPerFace);
   for(const auto& id : faceIds)
   {
-    auto bounds = FindBoundingBoxOfFace(faces, id);
+    auto bounds = FindBoundingBoxOfFace(cache, triangleGeom, id, verts);
     Point3Df min = bounds.getMinPoint();
     Point3Df max = bounds.getMaxPoint();