libigl_Martin/include/igl/is_delaunay.cpp

#include "is_delaunay.h"
#include "unique_edge_map.h"
#include <cassert>

template <
  typename DerivedV,
  typename DerivedF,
  typename DerivedD>
IGL_INLINE void igl::is_delaunay(
  const Eigen::MatrixBase<DerivedV> & V,
  const Eigen::MatrixBase<DerivedF> & F,
  Eigen::PlainObjectBase<DerivedD> & D)
{
  typedef Eigen::Matrix<typename DerivedF::Scalar,Eigen::Dynamic,2> MatrixX2I;
  typedef Eigen::Matrix<typename DerivedF::Scalar,Eigen::Dynamic,1> VectorXI;
  MatrixX2I E,uE;
  VectorXI EMAP;
  std::vector<std::vector<typename DerivedF::Scalar> > uE2E;
  igl::unique_edge_map(F, E, uE, EMAP, uE2E);
  const int num_faces = F.rows();
  D.setConstant(F.rows(),F.cols(),false);
  const auto D_at = [&D,&num_faces](const int he)->typename DerivedD::Scalar&
  {
    const int f = he%num_faces;
    const int c = he/num_faces;
    return D(f,c);
  };
  typedef typename DerivedV::Scalar Scalar;
  // Should use Shewchuk's predicates instead.
  const auto float_incircle = [](
    const Scalar pa[2],
    const Scalar pb[2],
    const Scalar pc[2],
    const Scalar pd[2])->short
  {
    const Eigen::Matrix3d A = (Eigen::Matrix3d(3,3)<<
      pa[0]-pd[0], pa[1]-pd[1],(pa[0]-pd[0])*(pa[0]-pd[0])+(pa[1]-pd[1])*(pa[1]-pd[1]),
      pb[0]-pd[0], pb[1]-pd[1],(pb[0]-pd[0])*(pb[0]-pd[0])+(pb[1]-pd[1])*(pb[1]-pd[1]),
      pc[0]-pd[0], pc[1]-pd[1],(pc[0]-pd[0])*(pc[0]-pd[0])+(pc[1]-pd[1])*(pc[1]-pd[1])
      ).finished();
    const Scalar detA = A.determinant();
    return (Scalar(0) < detA) - (detA < Scalar(0));
  };
  // loop over all unique edges
  for(int ue = 0;ue < uE2E.size(); ue++)
  {
    bool ue_is_d = false;
    // Is boundary?
    switch(uE2E[ue].size())
    {
      case 1:
        ue_is_d = true;
        break;
      case 2:
      {
        ue_is_d = is_delaunay(V,F,uE2E,float_incircle,ue);
        break;
      }
      default:
        ue_is_d = false;
        break;
    }
    for(int e = 0;e<uE2E[ue].size();e++)
    {
      D_at(uE2E[ue][e]) = ue_is_d;
    }
  }
}

template <
  typename DerivedV,
  typename DerivedF,
  typename uE2EType,
  typename InCircle,
  typename ueiType>
IGL_INLINE bool igl::is_delaunay(
  const Eigen::MatrixBase<DerivedV> & V,
  const Eigen::MatrixBase<DerivedF> & F,
  const std::vector<std::vector<uE2EType> > & uE2E,
  const InCircle incircle,
  const ueiType uei)
{
  const int num_faces = F.rows();
  typedef typename DerivedV::Scalar Scalar;
  const auto& half_edges = uE2E[uei];
  assert((half_edges.size() == 2)  && "uE2E[uei].size() should be 2");
  const size_t f1 = half_edges[0] % num_faces;
  const size_t f2 = half_edges[1] % num_faces;
  const size_t c1 = half_edges[0] / num_faces;
  const size_t c2 = half_edges[1] / num_faces;
  assert(c1 < 3);
  assert(c2 < 3);
  assert(f1 != f2);
  const size_t v1 = F(f1, (c1+1)%3);
  const size_t v2 = F(f1, (c1+2)%3);
  const size_t v4 = F(f1, c1);
  const size_t v3 = F(f2, c2);
  const Scalar p1[] = {V(v1, 0), V(v1, 1)};
  const Scalar p2[] = {V(v2, 0), V(v2, 1)};
  const Scalar p3[] = {V(v3, 0), V(v3, 1)};
  const Scalar p4[] = {V(v4, 0), V(v4, 1)};
  auto orientation = incircle(p1, p2, p4, p3);
  return orientation <= 0;
}

#ifdef IGL_STATIC_LIBRARY
// Explicit template instantiation
// generated by autoexplicit.sh
template void igl::is_delaunay<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<bool, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<bool, -1, -1, 0, -1, -1> >&);
// generated by autoexplicit.sh
template bool igl::is_delaunay<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, int, short (*)(double const*, double const*, double const*, double const*), unsigned long>(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > > const&, short (*)(double const*, double const*, double const*, double const*), unsigned long);
#endif