Merge pull request #1850 from llnl/feature/whitlock/klee_shaping_transform

Klee scaling enhancements

Author: BradWhitlock

RELEASE-NOTES.md

@@ -26,6 +26,7 @@ The Axom project release numbers follow [Semantic Versioning](http://semver.org/
 - Primal: Adds `NURBSPatch::isTriviallyTrimmed()` to check if the trimming curves for a patch lie on the patch boundaries
 - Quest: Adds support for reading mfem files with variable order NURBS curves (requires mfem>4.9).
 - Quest: Adds OMP support for fast GWN methods for STL/Triangulated STEP input and linearized NURBS Curve input.
+- Klee: Adds an optional "center" parameter in scale operators that permits scaling relative to a custom center point.
 
 ### Removed
 

src/axom/core/numerics/transforms.hpp

@@ -8,6 +8,7 @@
 
 #include "axom/config.hpp"
 #include "axom/core/numerics/Matrix.hpp"
+#include "axom/core/ArrayView.hpp"
 
 #include <cassert>
 #include <cmath>
@@ -133,6 +134,33 @@ Matrix<T> axisRotation(double angleRad, double x, double y, double z)
   return M;
 }
 
+/*!
+ * \brief Return translation matrix.
+ *
+ * \param tx The translation in x.
+ * \param ty The translation in y.
+ *
+ * \return A Matrix containing the translation transform.
+ */
+template <typename T = double>
+Matrix<T> translate(T tx, T ty)
+{
+  Matrix<T> M = Matrix<T>::identity(3);
+  M(0, 2) = tx;
+  M(1, 2) = ty;
+  return M;
+}
+
+template <typename T = double>
+Matrix<T> translate(T tx, T ty, T tz)
+{
+  Matrix<T> M = Matrix<T>::identity(4);
+  M(0, 3) = tx;
+  M(1, 3) = ty;
+  M(2, 3) = tz;
+  return M;
+}
+
 /*!
  * \brief Return scaling matrix.
  *
@@ -194,30 +222,66 @@ Matrix<T> scale(T sx, T sy, int ndims = 3)
 }
 
 /*!
- * \brief Return translation matrix.
+ * \brief Return scaling matrix relative to a center point.
  *
- * \param tx The translation in x.
- * \param ty The translation in y.
+ * \param sx The scaling value in x.
+ * \param sy The scaling value in y.
+ * \param center The center point.
  *
- * \return A Matrix containing the translation transform.
+ * \return A 3x3 Matrix containing the scaling transform.
  */
 template <typename T = double>
-Matrix<T> translate(T tx, T ty)
+Matrix<T> scale(T sx, T sy, const axom::ArrayView<T> &center)
 {
-  Matrix<T> M = Matrix<T>::identity(3);
-  M(0, 2) = tx;
-  M(1, 2) = ty;
-  return M;
+  assert(center.size() == 2);
+  const T zero {0};
+  if(axom::utilities::isNearlyEqual(center[0], zero) &&
+     axom::utilities::isNearlyEqual(center[1], zero))
+  {
+    return scale(sx, sy);
+  }
+
+  const auto T0 = translate(-center[0], -center[1]);
+  const auto S = scale(sx, sy, 3);
+  const auto T1 = translate(center[0], center[1]);
+
+  Matrix<T> TS(Matrix<T>::identity(3)), result(Matrix<T>::identity(3));
+  matrix_multiply(T1, S, TS);
+  matrix_multiply(TS, T0, result);
+
+  return result;
 }
 
+/*!
+ * \brief Return scaling matrix relative to a center point.
+ *
+ * \param sx The scaling value in x.
+ * \param sy The scaling value in y.
+ * \param center The center point.
+ *
+ * \return A 4x4 Matrix containing the scaling transform.
+ */
 template <typename T = double>
-Matrix<T> translate(T tx, T ty, T tz)
+Matrix<T> scale(T sx, T sy, T sz, const axom::ArrayView<T> &center)
 {
-  Matrix<T> M = Matrix<T>::identity(4);
-  M(0, 3) = tx;
-  M(1, 3) = ty;
-  M(2, 3) = tz;
-  return M;
+  assert(center.size() == 3);
+  const T zero {0};
+  if(axom::utilities::isNearlyEqual(center[0], zero) &&
+     axom::utilities::isNearlyEqual(center[1], zero) &&
+     axom::utilities::isNearlyEqual(center[2], zero))
+  {
+    return scale(sx, sy, sz, 4);
+  }
+
+  const auto T0 = translate(-center[0], -center[1], -center[2]);
+  const auto S = scale(sx, sy, sz, 4);
+  const auto T1 = translate(center[0], center[1], center[2]);
+
+  Matrix<T> TS(Matrix<T>::identity(4)), result(Matrix<T>::identity(4));
+  matrix_multiply(T1, S, TS);
+  matrix_multiply(TS, T0, result);
+
+  return result;
 }
 
 }  // end namespace transforms

src/axom/core/tests/CMakeLists.txt

@@ -46,6 +46,7 @@ set(core_serial_tests
     numerics_lu.hpp
     numerics_matrix.hpp
     numerics_matvecops.hpp
+    numerics_transforms.hpp
     numerics_polynomial_solvers.hpp
 
     utils_endianness.hpp
@@ -226,4 +227,3 @@ if (ENABLE_BENCHMARKS)
                       COMMAND ${test_name} --benchmark_min_time=0.0001s)
   endforeach()
 endif()
-

src/axom/core/tests/core_serial_main.cpp

@@ -40,6 +40,7 @@
 #include "numerics_lu.hpp"
 #include "numerics_matrix.hpp"
 #include "numerics_matvecops.hpp"
+#include "numerics_transforms.hpp"
 #include "numerics_polynomial_solvers.hpp"
 #include "numerics_quadrature.hpp"
 

src/axom/core/tests/core_static_array.hpp

@@ -25,6 +25,8 @@ struct DevicePair
 
   AXOM_HOST_DEVICE explicit DevicePair(int value) : first(value), second(-value) { }
 
+  AXOM_HOST_DEVICE DevicePair(const DevicePair& obj) : first(obj.first), second(obj.second) { }
+
   AXOM_HOST_DEVICE DevicePair& operator=(const DevicePair& other)
   {
     first = other.first;

src/axom/core/tests/numerics_transforms.hpp

@@ -0,0 +1,95 @@
+// Copyright (c) Lawrence Livermore National Security, LLC and other
+// Axom Project Contributors. See top-level LICENSE and COPYRIGHT
+// files for dates and other details.
+//
+// SPDX-License-Identifier: (BSD-3-Clause)
+
+#include "gtest/gtest.h"
+
+#include "axom/core/ArrayView.hpp"
+#include "axom/core/numerics/matvecops.hpp"
+#include "axom/core/numerics/transforms.hpp"
+
+namespace
+{
+void expect_matrix_near(const axom::numerics::Matrix<double>& actual,
+                        const axom::numerics::Matrix<double>& expected,
+                        double tolerance = 1e-12)
+{
+  ASSERT_EQ(actual.getNumRows(), expected.getNumRows());
+  ASSERT_EQ(actual.getNumColumns(), expected.getNumColumns());
+
+  for(axom::IndexType i = 0; i < actual.getNumRows(); ++i)
+  {
+    for(axom::IndexType j = 0; j < actual.getNumColumns(); ++j)
+    {
+      EXPECT_NEAR(actual(i, j), expected(i, j), tolerance);
+    }
+  }
+}
+
+void expect_vector_near(const double* actual, const double* expected, int size, double tolerance = 1e-12)
+{
+  for(int i = 0; i < size; ++i)
+  {
+    EXPECT_NEAR(actual[i], expected[i], tolerance);
+  }
+}
+}  // namespace
+
+TEST(numerics_transforms, scale_2d_about_center)
+{
+  double centerData[2] = {2.0, 3.0};
+  axom::ArrayView<double> center(centerData, 2);
+
+  auto actual = axom::numerics::transforms::scale(4.0, 5.0, center);
+
+  axom::numerics::Matrix<double> expected = axom::numerics::Matrix<double>::identity(3);
+  expected(0, 0) = 4.0;
+  expected(1, 1) = 5.0;
+  expected(0, 2) = -6.0;
+  expected(1, 2) = -12.0;
+  expect_matrix_near(actual, expected);
+
+  double point[3] = {3.0, 4.0, 1.0};
+  double result[3] = {0.0, 0.0, 0.0};
+  double expectedPoint[3] = {6.0, 8.0, 1.0};
+  axom::numerics::matrix_vector_multiply(actual, point, result);
+  expect_vector_near(result, expectedPoint, 3);
+}
+
+TEST(numerics_transforms, scale_3d_about_center)
+{
+  double centerData[3] = {0.5, 0.5, 0.5};
+  axom::ArrayView<double> center(centerData, 3);
+
+  auto actual = axom::numerics::transforms::scale(2.0, 3.0, 4.0, center);
+
+  axom::numerics::Matrix<double> expected = axom::numerics::Matrix<double>::identity(4);
+  expected(0, 0) = 2.0;
+  expected(1, 1) = 3.0;
+  expected(2, 2) = 4.0;
+  expected(0, 3) = -0.5;
+  expected(1, 3) = -1.0;
+  expected(2, 3) = -1.5;
+  expect_matrix_near(actual, expected);
+
+  double point[4] = {1.5, 1.5, 1.5, 1.0};
+  double result[4] = {0.0, 0.0, 0.0, 0.0};
+  double expectedPoint[4] = {2.5, 3.5, 4.5, 1.0};
+  axom::numerics::matrix_vector_multiply(actual, point, result);
+  expect_vector_near(result, expectedPoint, 4);
+}
+
+TEST(numerics_transforms, scale_about_zero_center_matches_origin_scale)
+{
+  double center2dData[2] = {0.0, 0.0};
+  axom::ArrayView<double> center2d(center2dData, 2);
+  expect_matrix_near(axom::numerics::transforms::scale(2.0, 3.0, center2d),
+                     axom::numerics::transforms::scale(2.0, 3.0));
+
+  double center3dData[3] = {0.0, 0.0, 0.0};
+  axom::ArrayView<double> center3d(center3dData, 3);
+  expect_matrix_near(axom::numerics::transforms::scale(2.0, 3.0, 4.0, center3d),
+                     axom::numerics::transforms::scale(2.0, 3.0, 4.0, 4));
+}

src/axom/klee/AffineMatrixVisitor.cpp

@@ -0,0 +1,49 @@
+// Copyright (c) Lawrence Livermore National Security, LLC and other
+// Axom Project Contributors. See top-level LICENSE and COPYRIGHT
+// files for dates and other details.
+//
+// SPDX-License-Identifier: (BSD-3-Clause)
+#include "axom/klee/AffineMatrixVisitor.hpp"
+
+namespace axom::klee
+{
+
+AffineMatrixVisitor::AffineMatrixVisitor()
+  : GeometryOperatorVisitor()
+  , m_isValid(false)
+  , m_matrix(4, 4)
+{ }
+
+void AffineMatrixVisitor::visit(const klee::Translation& translation)
+{
+  m_matrix = translation.toMatrix();
+  m_isValid = true;
+}
+void AffineMatrixVisitor::visit(const klee::Rotation& rotation)
+{
+  m_matrix = rotation.toMatrix();
+  m_isValid = true;
+}
+void AffineMatrixVisitor::visit(const klee::Scale& scale)
+{
+  m_matrix = scale.toMatrix();
+  m_isValid = true;
+}
+void AffineMatrixVisitor::visit(const klee::UnitConverter& converter)
+{
+  m_matrix = converter.toMatrix();
+  m_isValid = true;
+}
+
+void AffineMatrixVisitor::visit(const klee::CompositeOperator&)
+{
+  SLIC_WARNING_ROOT("CompositeOperator not supported for Shaper query");
+  m_isValid = false;
+}
+void AffineMatrixVisitor::visit(const klee::SliceOperator&)
+{
+  SLIC_WARNING_ROOT("SliceOperator not yet supported for Shaper query");
+  m_isValid = false;
+}
+
+}  // end namespace axom::klee

src/axom/klee/AffineMatrixVisitor.hpp

@@ -0,0 +1,42 @@
+// Copyright (c) Lawrence Livermore National Security, LLC and other
+// Axom Project Contributors. See top-level LICENSE and COPYRIGHT
+// files for dates and other details.
+//
+// SPDX-License-Identifier: (BSD-3-Clause)
+#ifndef AXOM_KLEE_AFFINE_MATRIX_VISITOR_HPP_
+#define AXOM_KLEE_AFFINE_MATRIX_VISITOR_HPP_
+#include "axom/klee/GeometryOperators.hpp"
+
+namespace axom::klee
+{
+
+/*!
+ * \brief Implementation of a GeometryOperatorVisitor for processing klee shape operators
+ *
+ * This class extracts the matrix form of supported operators and marks the operator as unvalid otherwise
+ * To use, check the \a isValid() function after visiting and then call the \a getMatrix() function.
+ */
+class AffineMatrixVisitor : public GeometryOperatorVisitor
+{
+public:
+  AffineMatrixVisitor();
+
+  void visit(const klee::Translation& translation) override;
+  void visit(const klee::Rotation& rotation) override;
+  void visit(const klee::Scale& scale) override;
+  void visit(const klee::UnitConverter& converter) override;
+
+  void visit(const klee::CompositeOperator&) override;
+  void visit(const klee::SliceOperator&) override;
+
+  const numerics::Matrix<double>& getMatrix() const { return m_matrix; }
+  bool isValid() const { return m_isValid; }
+
+private:
+  bool m_isValid;
+  numerics::Matrix<double> m_matrix;
+};
+
+}  // end namespace axom::klee
+
+#endif

src/axom/klee/CMakeLists.txt

@@ -10,6 +10,7 @@ axom_component_requires(NAME       Klee
 # Specify all headers/sources
 #------------------------------------------------------------------------------
 set(klee_headers
+    AffineMatrixVisitor.hpp
     Dimensions.hpp
     Geometry.hpp
     GeometryOperators.hpp
@@ -26,6 +27,7 @@ set(klee_internal_headers
     )
 
 set(klee_sources
+    AffineMatrixVisitor.cpp
     Geometry.cpp
     GeometryOperators.cpp
     KleeError.cpp