// 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 "intrinsic_delaunay_triangulation.h"
#include "is_intrinsic_delaunay.h"
#include "tan_half_angle.h"
#include "unique_edge_map.h"
#include "flip_edge.h"
#include <iostream>

template <
  typename Derivedl_in,
  typename DerivedF_in,
  typename Derivedl,
  typename DerivedF>
IGL_INLINE void igl::intrinsic_delaunay_triangulation(
  const Eigen::MatrixBase<Derivedl_in> & l_in,
  const Eigen::MatrixBase<DerivedF_in> & F_in,
  Eigen::PlainObjectBase<Derivedl> & l,
  Eigen::PlainObjectBase<DerivedF> & F)
{
  // We're going to work in place
  l = l_in;
  F = F_in;

  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);
  typedef typename DerivedF::Scalar Index;
  typedef typename Derivedl::Scalar Scalar;
  const Index num_faces = F.rows();

  // Copied from delaunay_triangulation
  bool all_delaunay = false;
  while(!all_delaunay) 
  {
    all_delaunay = true;
    for (size_t uei=0; uei<uE2E.size(); uei++) 
    {
      if (uE2E[uei].size() == 2) 
      {
        if(!is_intrinsic_delaunay(l,F,uE2E,uei)) 
        {
          all_delaunay = false;
          // update l just before flipping edge
          //      .        //
          //     /|\       //
          //   a/ | \d     //
          //   /  e  \     //
          //  /   |   \    //
          // .----|-f--.   //
          //  \   |   /    //
          //   \  |  /     //
          //   b\α|δ/c     //
          //     \|/       //
          //      .        //
          // Compute intrinsic length of oppposite edge
          assert(uE2E[uei].size() == 2 && "edge should have 2 incident faces");
          const Index he1 = uE2E[uei][0];
          const Index he2 = uE2E[uei][1];
          const Index f1 = he1%num_faces;
          const Index c1 = he1/num_faces;
          const Index f2 = he2%num_faces;
          const Index c2 = he2/num_faces;
          assert( std::abs(l(f1,c1)-l(f2,c2) < igl::EPS<Scalar>()) );
          const Scalar e = l(f1,c1);
          const Scalar a = l(f1,(c1+1)%3);
          const Scalar b = l(f1,(c1+2)%3);
          const Scalar c = l(f2,(c2+1)%3);
          const Scalar d = l(f2,(c2+2)%3);
          // tan(α/2)
          const Scalar tan_a_2= tan_half_angle(a,b,e);
          // tan(δ/2)
          const Scalar tan_d_2 = tan_half_angle(d,e,c);
          // tan((α+δ)/2)
          const Scalar tan_a_d_2 = (tan_a_2 + tan_d_2)/(1.0-tan_a_2*tan_d_2);
          // cos(α+δ)
          const Scalar cos_a_d = 
            (1.0 - tan_a_d_2*tan_a_d_2)/(1.0+tan_a_d_2*tan_a_d_2);
          const Scalar f = sqrt(b*b + c*c - 2.0*b*c*cos_a_d);
          // Annotated from flip_edge:
          // Edge to flip [v1,v2] --> [v3,v4]
          // Before:
          // F(f1,:) = [v1,v2,v4] // in some cyclic order
          // F(f2,:) = [v1,v3,v2] // in some cyclic order
          // After: 
          // F(f1,:) = [v1,v3,v4] // in *this* order 
          // F(f2,:) = [v2,v4,v3] // in *this* order
          //
          //          v1                 v1
          //          /|\                / \
          //        c/ | \b            c/f1 \b
          //     v3 /f2|f1\ v4  =>  v3 /__f__\ v4
          //        \  e  /            \ f2  /
          //        d\ | /a            d\   /a
          //          \|/                \ /
          //          v2                 v2
          //
          l(f1,0) = f;
          l(f1,1) = b;
          l(f1,2) = c;
          l(f2,0) = f;
          l(f2,1) = d;
          l(f2,2) = a;
          flip_edge(F, E, uE, EMAP, uE2E, uei);
        }
      }
    }
  }
}

#ifdef IGL_STATIC_LIBRARY
// Explicit template instantiation
// generated by autoexplicit.sh
template void igl::intrinsic_delaunay_triangulation<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -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<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
#endif