Geometry: Added spherical geometry methods. Includes unit tests.

git-svn-id: http://www.openflipper.org/svnrepo/OpenFlipper/branches/Free@15498 383ad7c9-94d9-4d36-a494-682f7c89f535

ACG/CMakeLists.txt

@@ -89,6 +89,8 @@ if (GTEST_FOUND)
   link_directories ( ${GTEST_LIBRARY_DIR})
 
   add_executable (ACG_tests ${TEST_SOURCES})
+  set_target_properties(ACG_tests PROPERTIES
+  		COMPILE_FLAGS "-UNDEBUG")
   target_link_libraries(ACG_tests
         ${GTEST_LIBRARIES} ${OPENMESH_LIBRARY}
   )

ACG/Geometry/Spherical.hh

+/*
+ * Spherical.hh
+ *
+ *  Created on: Sep 18, 2012
+ *      Author: ebke
+ */
+
+#ifndef ACG_GEOMETRY_SPHERICAL_HH_
+#define ACG_GEOMETRY_SPHERICAL_HH_
+
+#include <ACG/Math/GLMatrixT.hh>
+#include <cmath>
+#include <stdexcept>
+
+namespace ACG {
+namespace Geometry {
+
+static const double epsilon = 1e-6;
+
+namespace Spherical_Private {
+
+template<class Vec>
+static inline typename Vec::value_type angleBetween(const Vec &v1, const Vec &v2, const Vec &normal) {
+    typedef typename Vec::value_type Scalar;
+
+    /*
+     * Handle unstable 0 degree special case.
+     * (0 degrees can easily become 360 degrees.)
+     */
+    const Scalar dotProd = v1 | v2;
+    if (std::abs(dotProd - 1.0) < epsilon) return 0;
+
+    /*
+     * General case: Use determinant to check >/< 180 degree.
+     */
+    const Scalar det = GLMatrixT<Scalar>(normal, v1, v2).determinant();
+
+    /*
+     * Interpret arc cos accrodingly.
+     */
+    const Scalar arcos_angle = std::acos(std::max(-1.0, std::min(1.0, dotProd)));
+
+    if (det >= -1e-6)
+        return arcos_angle;
+    else
+        return 2 * M_PI - arcos_angle;
+}
+
+template<class Vec>
+static inline typename Vec::value_type angleBetween(const Vec &v1, const Vec &v2) {
+    const Vec normal = (v1 % v2).normalized();
+    return angleBetween(v1, v2, normal);
+}
+
+} /* namespace Spherical_Private */
+
+/**
+ * Compute inner angle sum of the spherical triangle specified by
+ * the three supplied unit vectors.
+ *
+ * @param n0 Must be unit length.
+ * @param n1 Must be unit length.
+ * @param n2 Must be unit length.
+ *
+ * @return Inner angle sum of the specified spherical triangle.
+ */
+template<class Vec>
+static inline typename Vec::value_type sphericalInnerAngleSum(const Vec &n0, const Vec &n1, const Vec &n2) {
+    typedef typename Vec::value_type Scalar;
+
+    const Scalar a = Spherical_Private::angleBetween(n0, n1);
+    const Scalar b = Spherical_Private::angleBetween(n1, n2);
+    const Scalar c = Spherical_Private::angleBetween(n2, n0);
+    if (a < epsilon || b < epsilon || c < epsilon) return M_PI;
+
+    const Scalar s = .5 * (a + b + c);
+    const Scalar sin_s = std::sin(s);
+    const Scalar sin_a = std::sin(a);
+    const Scalar sin_b = std::sin(b);
+    const Scalar sin_c = std::sin(c);
+
+#ifndef NDEBUG
+    if (std::sin(s - a) < -1e-4) throw std::logic_error("ACG::Geometry::sphericalInnerAngleSum(): 0 > s-a");
+    if (std::sin(s - b) < -1e-4) throw std::logic_error("ACG::Geometry::sphericalInnerAngleSum(): 0 > s-b");
+    if (std::sin(s - c) < -1e-4) throw std::logic_error("ACG::Geometry::sphericalInnerAngleSum(): 0 > s-c");
+#endif
+
+    const Scalar alpha_2 = std::acos(std::min(1.0, std::sqrt(sin_s * std::max(.0, std::sin(s - a)) / (sin_b * sin_c))));
+    const Scalar beta_2  = std::acos(std::min(1.0, std::sqrt(sin_s * std::max(.0, std::sin(s - b)) / (sin_c * sin_a))));
+    const Scalar gamma_2 = std::acos(std::min(1.0, std::sqrt(sin_s * std::max(.0, std::sin(s - c)) / (sin_a * sin_b))));
+
+    return 2 * (alpha_2 + beta_2 + gamma_2);
+}
+
+/**
+ * Compute gauss curvature of spherical polyhedral spanned by the
+ * supplied unit length normals. Normals have to be supplied in
+ * CCW order.
+ *
+ * @return Gauss curvature of the supplied spherical polyhedral.
+ */
+template<class Vec, class INPUT_ITERATOR>
+static inline typename Vec::value_type sphericalPolyhedralGaussCurv(INPUT_ITERATOR normals_begin, INPUT_ITERATOR normals_end) {
+    typedef typename Vec::value_type Scalar;
+
+    if (normals_begin == normals_end) return 0;
+    Vec n0 = *(normals_begin++);
+
+    if (normals_begin == normals_end) return 0;
+    Vec n2 = *(normals_begin++);
+
+    Scalar result = 0;
+    while (normals_begin != normals_end) {
+        /*
+         * Next triangle.
+         */
+        const Vec n1 = n2;
+        n2 = *(normals_begin++);
+
+        const Scalar sign = GLMatrixT<Scalar>(n0, n1, n2).determinant() >= 0 ? 1 : -1;
+        result += sign * (sphericalInnerAngleSum(n0, n1, n2) - M_PI);
+    }
+
+    return result;
+}
+
+} /* namespace Geometry */
+} /* namespace ACG */
+
+#endif /* ACG_GEOMETRY_SPHERICAL_HH_ */

ACG/tests/Geometry/Spherical_test.cc

+/*
+ * Spherical_test.cc
+ *
+ *  Created on: Sep 18, 2012
+ *      Author: ebke
+ */
+
+#include <gtest/gtest.h>
+
+#include <ACG/Math/VectorT.hh>
+#include <ACG/Math/GLMatrixT.hh>
+#include <ACG/Geometry/Spherical.hh>
+#include <boost/assign.hpp>
+
+namespace {
+
+using ACG::Vec3d;
+
+inline Vec3d rot(const Vec3d &ref, const Vec3d &normal, double angle) {
+    ACG::GLMatrixT<double> rmat;
+    rmat.identity();
+    rmat.rotate(angle * M_1_PI * 180.0, normal, ACG::MULT_FROM_LEFT);
+    return rmat.transform_vector(ref);
+}
+
+class Spherical : public testing::Test {
+
+    protected:
+
+        virtual void SetUp() {
+        }
+
+        virtual void TearDown() {
+        }
+
+};
+
+TEST_F(Spherical, sphericalInnerAngleSum_zeroTriangle) {
+    {
+        const Vec3d n1(1, 0, 0);
+        EXPECT_NEAR(M_PI, ACG::Geometry::sphericalInnerAngleSum(n1, n1, n1), 1e-6);
+    }
+
+    /*
+     * Jitter along one axis.
+     */
+    {
+        const Vec3d n1(1, 0, 0);
+        const Vec3d axis = (n1 % Vec3d(3, 1, 2).normalized()).normalized();
+        EXPECT_NEAR(M_PI, ACG::Geometry::sphericalInnerAngleSum(
+                rot(n1, axis, .1), n1, n1), 1e-6);
+        EXPECT_NEAR(M_PI, ACG::Geometry::sphericalInnerAngleSum(
+                n1, rot(n1, axis, .1), n1), 1e-6);
+        EXPECT_NEAR(M_PI, ACG::Geometry::sphericalInnerAngleSum(
+                n1, n1, rot(n1, axis, .1)), 1e-6);
+        EXPECT_NEAR(M_PI, ACG::Geometry::sphericalInnerAngleSum(
+                rot(n1, axis, .1), rot(n1, axis, -.05), rot(n1, axis, .07)), 1e-6);
+    }
+
+    {
+        const Vec3d n1 = Vec3d(4, 5, 6).normalized();
+        const Vec3d axis = (n1 % Vec3d(3, 1, 2).normalized()).normalized();
+        EXPECT_NEAR(M_PI, ACG::Geometry::sphericalInnerAngleSum(
+                rot(n1, axis, .1), rot(n1, axis, -.05), rot(n1, axis, .07)), 1e-6);
+    }
+}
+
+TEST_F(Spherical, sphericalPolyhedralGaussCurv_pointPolyhedral) {
+    std::vector<Vec3d> normals = boost::assign::list_of<Vec3d>
+    (Vec3d(1, 0, 0))
+    (Vec3d(1, 0, 0))
+    (Vec3d(1, 0, 0))
+    (Vec3d(1, 0, 0));
+
+    EXPECT_NEAR(0, ACG::Geometry::sphericalPolyhedralGaussCurv<Vec3d>(normals.begin(), normals.end()), 1e-6);
+
+    const Vec3d v = Vec3d(3, 2, 7).normalized();
+    normals.clear();
+    for (int i = 0; i < 7; ++i) normals.push_back(v);
+
+    EXPECT_NEAR(0, ACG::Geometry::sphericalPolyhedralGaussCurv<Vec3d>(normals.begin(), normals.end()), 1e-6);
+}
+
+TEST_F(Spherical, sphericalPolyhedralGaussCurv_linePolyhedral) {
+    /*
+     * Jitter along one axis.
+     */
+    const Vec3d n1 = Vec3d(4, 5, 6).normalized();
+    const Vec3d axis = (n1 % Vec3d(3, 1, 2)).normalized();
+
+    const std::vector<Vec3d> normals = boost::assign::list_of<Vec3d>
+    (rot(n1, axis, .1))
+    (rot(n1, axis, .2))
+    (rot(n1, axis, .05))
+    (rot(n1, axis, .09))
+    (rot(n1, axis, -.2))
+    (rot(n1, axis, .01))
+    (rot(n1, axis, -.1))
+    (rot(n1, axis, -.2));
+
+    EXPECT_NEAR(0, ACG::Geometry::sphericalPolyhedralGaussCurv<Vec3d>(normals.begin(), normals.end()), 1e-6);
+}
+
+TEST_F(Spherical, sphericalPolyhedralGaussCurv_cubeCorner) {
+    {
+        std::vector<Vec3d> normals = boost::assign::list_of<Vec3d>
+        (Vec3d(1, 0, 0))
+        (Vec3d(0, 1, 0))
+        (Vec3d(0, 0, 1));
+
+        EXPECT_NEAR(M_PI_2, ACG::Geometry::sphericalPolyhedralGaussCurv<Vec3d>(normals.begin(), normals.end()), 1e-6);
+    }
+}
+
+
+TEST_F(Spherical, sphericalPolyhedralGaussCurv_houseCorner) {
+    std::vector<Vec3d> normals = boost::assign::list_of<Vec3d>
+    (Vec3d(0, 0, 1))
+    (Vec3d(0, 1, 0))
+    (Vec3d(0, 1, 0))
+    (Vec3d(0, 1, 0))
+    (Vec3d(1, 0, 0));
+
+    EXPECT_NEAR(-M_PI_2, ACG::Geometry::sphericalPolyhedralGaussCurv<Vec3d>(normals.begin(), normals.end()), 1e-6);
+}
+
+} /* namespace */