Geometry.h

Go to the documentation of this file.

 
 
#pragma once
 
#include <Debug.h>
#include <array>
#include <cmath>
 
namespace ttk {
 
  namespace Geometry {
 
    template <typename T>
    T angle(const T *vA0, const T *vA1, const T *vB0, const T *vB1);
 
    template <typename T>
    T angle2D(const T *vA0, const T *vA1, const T *vB0, const T *vB1);
 
    // Computes the BAC angle, in the interval [0, 2pi]
    template <typename T>
    inline double angle2DUndirected(const T *vA, const T *vB, const T *vC) {
      double angle = angle2D<T>(vA, vB, vA, vC);
      if(angle < 0) {
        angle += 2 * M_PI;
      }
      return angle;
    }
 
    template <typename T>
    bool areVectorsColinear(const T *vA0,
                            const T *vA1,
                            const T *vB0,
                            const T *vB1,
                            std::array<T, 3> *coefficients = nullptr,
                            const T *tolerance = NULL);
    template <typename T>
    bool isTriangleColinear2D(const T *pptA,
                              const T *pptB,
                              const T *pptC,
                              const T tolerance);
 
    template <typename T>
    int computeBarycentricCoordinates(const T *p0,
                                      const T *p1,
                                      const T *p,
                                      std::array<T, 2> &baryCentrics,
                                      const int &dimension = 3);
 
    template <typename T>
    int computeBarycentricCoordinates(const T *p0,
                                      const T *p1,
                                      const T *p2,
                                      const T *p,
                                      std::array<T, 3> &baryCentrics);
 
    template <typename T>
    bool computeSegmentIntersection(const T &xA,
                                    const T &yA,
                                    const T &xB,
                                    const T &yB,
                                    const T &xC,
                                    const T &yC,
                                    const T &xD,
                                    const T &yD,
                                    T &x,
                                    T &y);
 
    template <typename T>
    int computeTriangleAngles(const T *p0,
                              const T *p1,
                              const T *p2,
                              std::array<T, 3> &angles);
 
    // Get the angle opposite to edge s2 using cosine law
    template <typename T>
    int computeTriangleAngleFromSides(const T s0,
                                      const T s1,
                                      const T s2,
                                      T &angle);
 
    template <typename T>
    int computeTriangleArea(const T *p0, const T *p1, const T *p2, T &area);
 
    template <typename T>
    int
      computeTriangleAreaFromSides(const T s0, const T s1, const T s2, T &area);
 
    template <typename T>
    int crossProduct(const T *vA0,
                     const T *vA1,
                     const T *vB0,
                     const T *vB1,
                     std::array<T, 3> &crossProduct);
 
    template <typename T>
    int crossProduct(const T *vA, const T *vB, T *crossProduct);
 
    template <typename T>
    T distance(const T *p0, const T *p1, const int &dimension = 3);
 
    template <typename T>
    T distance(const std::vector<T> &p0, const std::vector<T> &p1);
 
    template <typename T>
    inline T distance2D(const T *p0, const T *p1) {
      T distance = 0;
      T di0 = (p0[0] - p1[0]);
      T di1 = (p0[1] - p1[1]);
      if(std::is_same<T, float>::value) { // TODO constexpr when c++17
        distance = fmaf(di0, di0, distance);
        distance = fmaf(di1, di1, distance);
      } else {
        distance = fma(di0, di0, distance);
        distance = fma(di1, di1, distance);
      }
 
      return sqrt(distance);
    }
 
    template <typename T>
    T distanceFlatten(const std::vector<std::vector<T>> &p0,
                      const std::vector<std::vector<T>> &p1);
 
    template <typename T>
    T distanceFlatten(const std::vector<std::vector<T>> &p0,
                      const std::vector<std::vector<T>> &p1);
 
    template <typename T>
    T dotProduct(const T *vA0, const T *vA1, const T *vB0, const T *vB1);
 
    template <typename T>
    T dotProduct(const T *vA, const T *vB, const int &dimension = 3);
 
    template <typename T>
    T dotProduct(const std::vector<T> &vA, const std::vector<T> &vB);
 
    template <typename T>
    T dotProductFlatten(const std::vector<std::vector<T>> &vA,
                        const std::vector<std::vector<T>> &vB);
 
    template <typename T, typename Container, size_t dim>
    int getBoundingBox(const Container &points,
                       std::array<std::pair<T, T>, dim> &bBox) {
      if(points.empty()) {
        return -1;
      }
 
      for(size_t i = 0; i < points.size(); i++) {
 
        if(i == 0) {
          for(size_t j = 0; j < dim; j++) {
            bBox[j].first = points[i][j];
            bBox[j].second = points[i][j];
          }
        } else {
          for(size_t j = 0; j < dim; j++) {
            if(points[i][j] < bBox[j].first) {
              bBox[j].first = points[i][j];
            }
            if(points[i][j] > bBox[j].second) {
              bBox[j].second = points[i][j];
            }
          }
        }
      }
 
      return 0;
    }
 
    template <typename T>
    bool isPointInTriangle(const T *p0, const T *p1, const T *p2, const T *p);
 
    template <typename T>
    bool isPointOnSegment(const T &x,
                          const T &y,
                          const T &xA,
                          const T &yA,
                          const T &xB,
                          const T &yB);
 
    template <typename T>
    bool isPointOnSegment(const T *p,
                          const T *pA,
                          const T *pB,
                          const int &dimension = 3);
 
    template <typename T>
    bool isTriangleColinear(const T *p0,
                            const T *p1,
                            const T *p2,
                            const T *tolerance = nullptr);
 
    template <typename T>
    T magnitude(const T *v, const int &dimension = 3);
 
    template <typename T>
    T magnitude(const std::vector<T> &v);
 
    template <typename T>
    T magnitudeFlatten(const std::vector<std::vector<T>> &v);
 
    template <typename T>
    T magnitude(const T *o, const T *d);
 
    template <typename T>
    inline T powInt(const T val, const int n) {
      if(n < 0) {
        return 1.0 / powInt(val, -n);
      } else if(n == 0) {
        return 1;
      } else if(n == 1) {
        return val;
      } else if(n == 2) {
        return val * val;
      } else if(n == 3) {
        return val * val * val;
      } else {
        T ret = val;
        for(int i = 0; i < n - 1; ++i) {
          ret *= val;
        }
        return ret;
      }
    }
 
    template <typename T = double>
    inline T powIntTen(const int n) {
      return powInt(static_cast<T>(10), n);
    }
 
#define TTK_POW_LAMBDA(CALLEXPR, TYPE, EXPN, ...)                      \
  {                                                                    \
    const auto one = [](const TYPE ttkNotUsed(a)) { return TYPE{1}; }; \
    const auto id = [](const TYPE a) { return a; };                    \
    const auto square = [](const TYPE a) { return a * a; };            \
    const auto cube = [](const TYPE a) { return a * a * a; };          \
    const auto powInt                                                  \
      = [EXPN](const TYPE a) { return Geometry::powInt(a, EXPN); };    \
                                                                       \
    if(EXPN == 0) {                                                    \
      CALLEXPR(__VA_ARGS__, one);                                      \
    } else if(EXPN == 1) {                                             \
      CALLEXPR(__VA_ARGS__, id);                                       \
    } else if(EXPN == 2) {                                             \
      CALLEXPR(__VA_ARGS__, square);                                   \
    } else if(EXPN == 3) {                                             \
      CALLEXPR(__VA_ARGS__, cube);                                     \
    } else {                                                           \
      CALLEXPR(__VA_ARGS__, powInt);                                   \
    }                                                                  \
  }
 
    template <typename T1, typename T2>
    inline T1 pow(const T1 val, const T2 n) {
      static_assert(
        std::is_arithmetic<T1>::value && std::is_arithmetic<T2>::value,
        "pow can only be applied on arithmetic values");
 
      if(std::is_integral<T2>::value) {
        return powInt(val, n);
      } else if(std::is_floating_point<T2>::value) {
        return std::pow(val, n);
      }
      // this return statement should be unreachable thanks to the
      // static_assert
      return T1{};
    }
 
    template <typename T>
    int
      subtractVectors(const T *a, const T *b, T *out, const int &dimension = 3);
 
    template <typename T>
    int subtractVectors(const std::vector<T> &a,
                        const std::vector<T> &b,
                        std::vector<T> &out);
 
    template <typename T>
    int addVectors(const T *a, const T *b, T *out, const int &dimension = 3);
 
    template <typename T>
    int addVectors(const std::vector<T> &a,
                   const std::vector<T> &b,
                   std::vector<T> &out);
 
    template <typename T>
    int multiAddVectors(const std::vector<std::vector<T>> &a,
                        const std::vector<std::vector<T>> &b,
                        std::vector<std::vector<T>> &out);
 
    template <typename T>
    int
      multiAddVectorsFlatten(const std::vector<std::vector<std::vector<T>>> &a,
                             const std::vector<std::vector<std::vector<T>>> &b,
                             std::vector<std::vector<T>> &out);
 
    template <typename T>
    int
      scaleVector(const T *a, const T factor, T *out, const int &dimension = 3);
 
    template <typename T>
    int
      scaleVector(const std::vector<T> &a, const T factor, std::vector<T> &out);
 
    template <typename T>
    int vectorProjection(const T *a,
                         const T *b,
                         T *out,
                         const int &dimension = 3);
 
    template <typename T>
    int vectorProjection(const std::vector<T> &a,
                         const std::vector<T> &b,
                         std::vector<T> &out);
 
    template <typename T>
    void addVectorsProjection(const std::vector<T> &a,
                              const std::vector<T> &b,
                              std::vector<T> &a_out,
                              std::vector<T> &b_out);
 
    template <typename T>
    void gramSchmidt(const std::vector<std::vector<T>> &a,
                     std::vector<std::vector<T>> &out);
 
    template <typename T>
    bool isVectorUniform(const std::vector<T> &a);
 
    template <typename T>
    bool isVectorNull(const std::vector<T> &a);
 
    template <typename T>
    bool isVectorNullFlatten(const std::vector<std::vector<T>> &a);
 
    template <typename T>
    int flattenMultiDimensionalVector(const std::vector<std::vector<T>> &a,
                                      std::vector<T> &out);
 
    template <typename T>
    int multiFlattenMultiDimensionalVector(
      const std::vector<std::vector<std::vector<T>>> &a,
      std::vector<std::vector<T>> &out);
 
    template <typename T>
    int unflattenMultiDimensionalVector(const std::vector<T> &a,
                                        std::vector<std::vector<T>> &out,
                                        const int &no_columns = 2);
 
    template <typename T>
    void matrixMultiplication(const std::vector<std::vector<T>> &a,
                              const std::vector<std::vector<T>> &b,
                              std::vector<std::vector<T>> &out);
 
    template <typename T>
    void subtractMatrices(const std::vector<std::vector<T>> &a,
                          const std::vector<std::vector<T>> &b,
                          std::vector<std::vector<T>> &out);
 
    template <typename T>
    void addMatrices(const std::vector<std::vector<T>> &a,
                     const std::vector<std::vector<T>> &b,
                     std::vector<std::vector<T>> &out);
 
    template <typename T>
    void scaleMatrix(const std::vector<std::vector<T>> &a,
                     const T factor,
                     std::vector<std::vector<T>> &out);
 
    template <typename T>
    void transposeMatrix(const std::vector<std::vector<T>> &a,
                         std::vector<std::vector<T>> &out);
 
  } // namespace Geometry
} // namespace ttk