// This file is part of libigl, a simple c++ geometry processing library.
// 
// Copyright (C) 2018 Alec Jacobson <alecjacobson@gmail.com>
// 
// This Source Code Form is subject to the terms of the Mozilla Public License 
// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
// obtain one at http://mozilla.org/MPL/2.0/.
#include "is_intrinsic_delaunay.h"
#include "unique_edge_map.h"
#include <cassert>
template <
  typename Derivedl,
  typename DerivedF,
  typename DerivedD>
IGL_INLINE void igl::is_intrinsic_delaunay(
  const Eigen::MatrixBase<Derivedl> & l,
  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);
  // loop over all unique edges
  for(int ue = 0;ue < uE2E.size(); ue++)
  {
    const bool ue_is_d = is_intrinsic_delaunay(l,F,uE2E,ue);
    // Set for all instances
    for(int e = 0;e<uE2E[ue].size();e++)
    {
      D( uE2E[ue][e]%num_faces, uE2E[ue][e]/num_faces) = ue_is_d;
    }
  }
}

template <
  typename Derivedl,
  typename DerivedF,
  typename uE2EType,
  typename ueiType>
IGL_INLINE bool igl::is_intrinsic_delaunay(
  const Eigen::MatrixBase<Derivedl> & l,
  const Eigen::MatrixBase<DerivedF> & F,
  const std::vector<std::vector<uE2EType> > & uE2E,
  const ueiType uei)
{
  if(uE2E[uei].size() == 1) return true;
  if(uE2E[uei].size() > 2) return false;
  typedef typename Derivedl::Scalar Scalar;
  typedef typename DerivedF::Scalar Index;

  //      .
  //     /|
  //   c/ |
  //   /  |
  //  /   |
  // .α   | a
  //  \   |
  //   \  |
  //   b\ |  
  //     \| 
  //
  // tan(α/2)
  const auto tan_alpha_over_2 = [](
    const Scalar & a,
    const Scalar & b,
    const Scalar & c)->Scalar
  {
    // Fisher 2007
    return sqrt(((a-b+c)*(a+b-c))/((a+b+c)*(-a+b+c)));
  };

  const auto cot_alpha = [&tan_alpha_over_2](
    const Scalar & a,
    const Scalar & b,
    const Scalar & c)->Scalar
  {
    // Fisher 2007
    const Scalar t = tan_alpha_over_2(a,b,c);
    return (1.0-t*t)/(2*t);
  };

  //      .        //
  //     /|\       //
  //   a/ | \d     //
  //   /  e  \     //
  //  /   |   \    //
  // .α---|-f-β.   //
  //  \   |   /    //
  //   \  |  /     //
  //   b\ | /c     //
  //     \|/       //
  //      .        //

  const Index num_faces = F.rows();
  assert(uE2E[uei].size() == 2 && "edge should have 2 incident faces");
  const Index he_left = uE2E[uei][0];
  const Index he_right = uE2E[uei][1];
  const Index f_left = he_left%num_faces;
  const Index c_left = he_left/num_faces;
  const Index f_right = he_right%num_faces;
  const Index c_right = he_right/num_faces;

  assert( std::abs(l(f_left,c_left)-l(f_right,c_right) < igl::EPS<Scalar>()) );
  const Scalar e = l(f_left,c_left);
  const Scalar a = l(f_left,(c_left+1)%3);
  const Scalar b = l(f_left,(c_left+2)%3);
  const Scalar c = l(f_right,(c_right+1)%3);
  const Scalar d = l(f_right,(c_right+2)%3);

  const Scalar w = cot_alpha(e,a,b) + cot_alpha(e,c,d);
  return w >= 0;
}

#ifdef IGL_STATIC_LIBRARY
// Explicit template instantiation
// generated by autoexplicit.sh
template void igl::is_intrinsic_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> >&);
#endif