// Copyright (c) 2018  Carnegie Mellon University (USA), GeometryFactory (France)
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org).
// You can redistribute it and/or modify it under the terms of the GNU
// General Public License as published by the Free Software Foundation,
// either version 3 of the License, or (at your option) any later version.
//
// Licensees holding a valid commercial license may use this file in
// accordance with the commercial license agreement provided with the software.
//
// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
//
// $URL: https://github.com/CGAL/cgal/blob/releases/CGAL-4.14.3/Heat_method_3/include/CGAL/Heat_method_3/Surface_mesh_geodesic_distances_3.h $
// $Id: Surface_mesh_geodesic_distances_3.h 96ef7fa 2019-01-07T11:45:44+01:00 Andreas Fabri
// SPDX-License-Identifier: GPL-3.0+
//
//
// Author(s) : Christina Vaz, Keenan Crane, Andreas Fabri

#ifndef CGAL_HEAT_METHOD_3_SURFACE_MESH_GEODESIC_DISTANCES_3_H
#define CGAL_HEAT_METHOD_3_SURFACE_MESH_GEODESIC_DISTANCES_3_H

#include <CGAL/license/Heat_method_3.h>
#include <CGAL/Heat_method_3/internal/Intrinsic_Delaunay_triangulation_3.h>
#include <CGAL/Heat_method_3/internal/V2V.h>
#include <CGAL/disable_warnings.h>

#include <CGAL/property_map.h>
#include <CGAL/double.h>
#include <CGAL/boost/graph/properties.h>
#include <CGAL/Dynamic_property_map.h>
#include <CGAL/squared_distance_3.h>
#include <CGAL/Polygon_mesh_processing/measure.h>
#include <CGAL/number_utils.h>
#ifdef CGAL_EIGEN3_ENABLED
#include <CGAL/Eigen_solver_traits.h>
#endif

#include <boost/foreach.hpp>
#include <boost/range/has_range_iterator.hpp>

#include <vector>
#include <set>

namespace CGAL {
namespace Heat_method_3 {

/**
 * A tag class used to specify that the heat method is applied to the input mesh.
 */
struct Direct
{};

/**
 * A tag class used to specify that the heat method applies the intrinsic Delaunay triangulation to the input mesh.
 */
struct Intrinsic_Delaunay
{};

namespace internal {
template <typename TriangleMesh,
          typename Traits,
          typename LA,
          typename VertexPointMap>
class Surface_mesh_geodesic_distances_3
{
protected:
  /// Polygon_mesh typedefs
  typedef typename boost::graph_traits<TriangleMesh>            graph_traits;
  typedef typename graph_traits::vertex_descriptor              vertex_descriptor;
  typedef typename graph_traits::edge_descriptor                edge_descriptor;
  typedef typename graph_traits::halfedge_descriptor            halfedge_descriptor;
  typedef typename graph_traits::face_descriptor                face_descriptor;
  typedef typename std::set<vertex_descriptor>                  Vertex_const_range;
  /// Geometric typedefs
  typedef typename Traits::Point_3                              Point_3;
  typedef typename Traits::FT                                   FT;
  typedef typename Traits::Vector_3                             Vector_3;

  // Property map typedefs
  typedef typename boost::property_traits<VertexPointMap>::reference VertexPointMap_reference;

  typedef typename LA::Matrix Matrix;
  typedef typename LA::Vector Vector;
  typedef typename LA::Index Index;

  // The Vertex_id_map is a property map where you can associate an index to a vertex_descriptor
  typedef typename boost::graph_traits<TriangleMesh>::vertices_size_type vertices_size_type;

  typedef CGAL::dynamic_vertex_property_t<Index> Vertex_property_tag;
  typedef typename boost::property_map<TriangleMesh, Vertex_property_tag >::const_type Vertex_id_map;
  Vertex_id_map vertex_id_map;

  typedef CGAL::dynamic_face_property_t<Index> Face_property_tag;
  typedef typename boost::property_map<TriangleMesh, Face_property_tag >::const_type Face_id_map;
  Face_id_map face_id_map;

public:

  /*!
    \brief Constructor
  */
  Surface_mesh_geodesic_distances_3(const TriangleMesh& tm)
    : vertex_id_map(get(Vertex_property_tag(),tm)), face_id_map(get(Face_property_tag(),tm)), v2v(tm), tm(tm), vpm(get(vertex_point,tm))
  {
    build();
  }


  /*!
    \brief Constructor
  */
  Surface_mesh_geodesic_distances_3(const TriangleMesh& tm, VertexPointMap vpm)
    : v2v(tm), tm(tm), vpm(vpm)
  {
    build();
  }


  /**
   * returns the triangle mesh the algorithm is running on.
   */
  const TriangleMesh& triangle_mesh() const{
    return tm;
  }


private:

  const Matrix&
  mass_matrix() const
  {
    return m_mass_matrix;
  }


  const Matrix&
  cotan_matrix() const
  {
    return m_cotan_matrix;
  }


  const VertexPointMap&
  vertex_point_map() const
  {
    return vpm;
  }


  const Vertex_id_map&
  get_vertex_id_map() const
  {
    return vertex_id_map;
  }

public:

  /**
   * adds `vd` to the source set, returning `false` iff `vd` is already in the set.
   */
  template <typename VD>
  bool
  add_source(VD vd)
  {
    m_source_change_flag = true;
    return m_sources.insert(v2v(vd)).second;
  }

  /**
   * adds vertices in the `vrange` to the source set'.
   * \tparam VertexConstRange a model of the concept `ConstRange` with value type `boost::graph_traits<TriangleMesh>::%vertex_descriptor`
   */
  template <typename VertexConstRange>
  void
  add_sources(const VertexConstRange& vrange)
  {
    typedef typename std::iterator_traits<typename VertexConstRange::const_iterator>::value_type value_type;
    m_source_change_flag = true;
    BOOST_FOREACH(value_type vd, vrange){
      m_sources.insert(v2v(vd));
    }
  }


  /**
   * removes `vd` from the source set, returning `true` iff `vd` was in the set.
   */
  template <typename VD>
  bool
  remove_source(VD vd)
  {
    m_source_change_flag = true;
    return (m_sources.erase(v2v(vd)) == 1);
  }


  /**
   * clears the current source set.
   */
  void
  clear_sources()
  {
    m_source_change_flag = true;
    m_sources.clear();
    return;
  }

  /**
   * returns the set of  source vertices.
   */

  const Vertex_const_range&
  sources() const
  {
    return m_sources;
  }

private:

  double
  summation_of_edges() const
  {
    typename Traits::Compute_squared_distance_3 squared_distance = Traits().compute_squared_distance_3_object();
    double edge_sum = 0;
    CGAL::Vertex_around_face_iterator<TriangleMesh> vbegin, vend, vmiddle;
    BOOST_FOREACH(face_descriptor f, faces(tm)) {
      halfedge_descriptor hd = halfedge(f,tm);

      vertex_descriptor neighbor_one = target(hd,tm);
      halfedge_descriptor hd2 = next(hd,tm);
      vertex_descriptor neighbor_two = target(hd2,tm);
      halfedge_descriptor hd3 = next(hd2,tm);
      vertex_descriptor current = target(hd3,tm);

      VertexPointMap_reference pi = get(vpm,current);
      VertexPointMap_reference pj = get(vpm, neighbor_one);
      VertexPointMap_reference pk = get(vpm, neighbor_two);
      edge_sum += (CGAL::is_border(opposite(hd,tm),tm)?1.0:0.5) * std::sqrt(to_double(squared_distance(pi,pj)));
      edge_sum += (CGAL::is_border(opposite(hd2,tm),tm)?1.0:0.5) * std::sqrt(to_double(squared_distance(pj,pk))) ;
      edge_sum += (CGAL::is_border(opposite(hd3,tm),tm)?1.0:0.5) * std::sqrt(to_double(squared_distance(pk,pi))) ;
    }
    return edge_sum;
  }


  double
  time_step() const
  {
    return m_time_step;
  }


  void
  update_kronecker_delta()
  {
    //currently just working with a single vertex in source set, add the first one for now
    Index i;
    Matrix K(static_cast<int>(num_vertices(tm)), 1);
    if(sources().empty()) {
      i = 0;
      K.set_coef(i,0, 1, true);
    } else {
      BOOST_FOREACH(vertex_descriptor vd, sources()){
        i = get(vertex_id_map, vd);
        K.set_coef(i,0, 1, true);
      }
    }
    m_kronecker.swap(K);
  }


  const Matrix&
  kronecker_delta() const
  {
    return m_kronecker;
  }


  void factor_cotan_laplace()
  {
    Matrix A, A0;
    A0 = m_time_step * m_cotan_matrix;
    A = m_mass_matrix + A0;

    double d=0;

    if(! la.factor(A,d)) {
      // decomposition failed
      CGAL_error_msg("Eigen Decomposition in cotan failed");
    }
  }

  void solve_cotan_laplace()
  {
    if(! la.linear_solver(m_kronecker, m_solved_u)) {
      // solving failed
      CGAL_error_msg("Eigen Solving in cotan failed");
    }
  }


  void
  compute_unit_gradient()
  {
    typename Traits::Construct_vector_3 construct_vector = Traits().construct_vector_3_object();
    typename Traits::Construct_sum_of_vectors_3 sum = Traits().construct_sum_of_vectors_3_object();
    typename Traits::Compute_scalar_product_3 scalar_product = Traits().compute_scalar_product_3_object();
    typename Traits::Construct_cross_product_vector_3 cross_product = Traits().construct_cross_product_vector_3_object();
    typename Traits::Construct_scaled_vector_3 scale = Traits().construct_scaled_vector_3_object();
    if(m_X.empty()){
      m_X.resize(num_faces(tm));
    }
    CGAL::Vertex_around_face_iterator<TriangleMesh> vbegin, vend, vmiddle;
    BOOST_FOREACH(face_descriptor f, faces(tm)) {
      boost::tie(vbegin, vend) = vertices_around_face(halfedge(f,tm),tm);
      vertex_descriptor current = *(vbegin);
      vertex_descriptor neighbor_one = *(++vbegin);
      vertex_descriptor neighbor_two = *(++vbegin);
      VertexPointMap_reference p_i = get(vpm,current);
      VertexPointMap_reference p_j = get(vpm, neighbor_one);
      VertexPointMap_reference p_k = get(vpm, neighbor_two);
      Index face_i = get(face_id_map, f);

      Index i = get(vertex_id_map, current);
      Index j = get(vertex_id_map, neighbor_one);
      Index k = get(vertex_id_map, neighbor_two);

      //get area of face_i
      //get outward unit normal
      //cross that with eij, ejk, eki
      //so (Ncross eij) *uk and so on
      //sum all of those then multiply by 1./(2a)

      Vector_3 v_ij = construct_vector(p_i,p_j);
      Vector_3 v_ik = construct_vector(p_i,p_k);

      Vector_3 cross = cross_product(v_ij, v_ik);
      double N_cross = (CGAL::sqrt(to_double(scalar_product(cross,cross))));
      Vector_3 unit_cross = scale(cross, 1./N_cross);
      double area_face = N_cross * (1./2);
      double u_i = CGAL::abs(m_solved_u(i));
      double u_j = CGAL::abs(m_solved_u(j));
      double u_k = CGAL::abs(m_solved_u(k));
      double r_Mag = 1./(std::max)((std::max)(u_i, u_j),u_k);
      /* normalize heat values so that they have roughly unit magnitude */
      if(!std::isinf(r_Mag)) {
        u_i = u_i * r_Mag;
        u_j = u_j * r_Mag;
        u_k = u_k * r_Mag;
      }
      Vector_3 edge_sums = scale(cross_product(unit_cross,v_ij), u_k);
      edge_sums = sum(edge_sums, scale(cross_product(unit_cross, construct_vector(p_j,p_k)), u_i));
      edge_sums = sum(edge_sums, scale(cross_product(unit_cross, construct_vector(p_k,p_i)), u_j));
      edge_sums = scale(edge_sums, (1./area_face));
      double e_magnitude = CGAL::sqrt(to_double(scalar_product(edge_sums,edge_sums)));
      m_X[face_i] = scale(edge_sums,(1./e_magnitude));
    }
  }


  void
  compute_divergence()
  {
    typename Traits::Construct_cross_product_vector_3 cross_product = Traits().construct_cross_product_vector_3_object();
    typename Traits::Compute_scalar_product_3 scalar_product = Traits().compute_scalar_product_3_object();
    typename Traits::Construct_vector_3 construct_vector = Traits().construct_vector_3_object();
    Matrix indexD(dimension,1);
    CGAL::Vertex_around_face_iterator<TriangleMesh> vbegin, vend, vmiddle;
    BOOST_FOREACH(face_descriptor f, faces(tm)) {
      boost::tie(vbegin, vend) = vertices_around_face(halfedge(f,tm),tm);
      vertex_descriptor current = *(vbegin);
      vertex_descriptor neighbor_one = *(++vbegin);
      vertex_descriptor neighbor_two = *(++vbegin);
      Index i = get(vertex_id_map, current);
      Index j = get(vertex_id_map, neighbor_one);
      Index k = get(vertex_id_map, neighbor_two);
      VertexPointMap_reference p_i = get(vpm,current);
      VertexPointMap_reference p_j = get(vpm, neighbor_one);
      VertexPointMap_reference p_k = get(vpm, neighbor_two);
      Index face_i = get(face_id_map, f);

      Vector_3 v_ij = construct_vector(p_i,p_j);
      Vector_3 v_ik = construct_vector(p_i,p_k);
      Vector_3 cross = cross_product(v_ij, v_ik);
      double norm_cross = CGAL::sqrt(to_double(scalar_product(cross,cross)));
      double dot = to_double(scalar_product(v_ij, v_ik));
      double cotan_i = dot/norm_cross;

      Vector_3 v_ji = construct_vector(p_j, p_i);
      Vector_3 v_jk = construct_vector(p_j, p_k);

      cross = cross_product(v_ji, v_jk);
      dot = to_double(scalar_product(v_ji, v_jk));
      double cotan_j = dot/norm_cross;

      Vector_3 v_ki = construct_vector(p_k,p_i);
      Vector_3 v_kj = construct_vector(p_k,p_j);

      cross = cross_product(v_ki, v_kj);
      dot = to_double(scalar_product(v_ki,v_kj));
      double cotan_k = dot/norm_cross;

      const Vector_3& a  = m_X[face_i];
      double i_entry = (to_double(scalar_product(a,v_ij)) * cotan_k) + (to_double(scalar_product(a,v_ik)) *  cotan_j);
      double j_entry = (to_double(scalar_product(a,v_jk)) * cotan_i) + (to_double(scalar_product(a,v_ji)) * cotan_k);
      double k_entry = (to_double(scalar_product(a,v_ki)) * cotan_j) + (to_double(scalar_product(a,v_kj)) * cotan_i);

      indexD.add_coef(i, 0, (1./2)*i_entry);
      indexD.add_coef(j, 0, (1./2)*j_entry);
      indexD.add_coef(k, 0, (1./2)*k_entry);
    }
    indexD.swap(m_index_divergence);
  }


  // modifies m_solved_phi
  void
  value_at_source_set(const Vector& phi)
  {
    Vector source_set_val(dimension);
    if(sources().empty()) {
      for(int k = 0; k<dimension; k++) {
        source_set_val(k,0) = phi.coeff(0,0);
      }
      source_set_val = phi-source_set_val;
    } else {
      for(int i = 0; i<dimension; i++) {
        double min_val = (std::numeric_limits<double>::max)();
        Index vd_index;
        //go through the distances to the sources and leave the minimum distance;
        BOOST_FOREACH(vertex_descriptor vd, sources()){
          vd_index = get(vertex_id_map, vd);
          double new_d = CGAL::abs(-phi.coeff(vd_index,0)+phi.coeff(i,0));
          if(phi.coeff(vd_index,0)==phi.coeff(i,0)) {
            min_val = 0.;
          }
          if(new_d < min_val) {
            min_val = new_d;
          }
        }
        source_set_val(i,0) = min_val;
      }
    }
    m_solved_phi.swap(source_set_val);
  }


  void
  factor_phi()
  {
    double d = 1;
    if(! la_cotan.factor(m_cotan_matrix, d)) {
      // decomposition failed
      CGAL_error_msg("Eigen Decomposition in solve_phi() failed");
    }
  }

  
  void
  solve_phi()
  {
    Vector phi;
  
    if(! la_cotan.linear_solver(m_index_divergence, phi)) {
      // solving failed
      CGAL_error_msg("Eigen Solving in solve_phi() failed");
    }
    value_at_source_set(phi);
  }


  // this function returns a (number of vertices)x1 vector where
  // the ith index has the distance from the first vertex to the ith vertex
  const Vector&
  distances() const
  {
    return m_solved_phi;
  }

public:

  /**
   *  Updates the distance property map after changes in the source set.
   **/
  template<class VertexDistanceMap>
  void estimate_geodesic_distances(VertexDistanceMap vdm)
  {
    if(is_empty(tm)){
      return;
    }
    double d=0;
    if(m_source_change_flag) {
      //don't need to recompute Mass matrix, cotan matrix or timestep reflect that in this function
      update_kronecker_delta(); //                                changes m_kronecker
      solve_cotan_laplace();    // depends on m_kronecker         changes m_solved_u
      compute_unit_gradient();  // depends on m_solved_u          changes m_X
      compute_divergence();     // depends on m_X                 changes m_index_divergence
      solve_phi();              // depends on m_index_divergence
      m_source_change_flag = false;
    }
    BOOST_FOREACH(vertex_descriptor vd, vertices(tm)){
      Index i_d = get(vertex_id_map, vd);
      d = m_solved_phi(i_d,0);
      put(vdm,vd,d);
    }
  }

private:

  void
  build()
  {
    typename Traits::Construct_cross_product_vector_3 cross_product = Traits().construct_cross_product_vector_3_object();
    typename Traits::Compute_scalar_product_3 scalar_product = Traits().compute_scalar_product_3_object();
    typename Traits::Construct_vector_3 construct_vector = Traits().construct_vector_3_object();

    if(is_empty(tm)){
      return;
    }
    m_source_change_flag = false;

    CGAL_precondition(is_triangle_mesh(tm));
    Index i = 0;
    BOOST_FOREACH(vertex_descriptor vd, vertices(tm)){
      put(vertex_id_map, vd, i++);
    }
    Index face_i = 0;
    BOOST_FOREACH(face_descriptor fd, faces(tm)){
      put(face_id_map, fd, face_i++);
    }
    dimension = static_cast<int>(num_vertices(tm));
    {
      Matrix M(dimension,dimension);
      m_mass_matrix.swap(M);
    }
    {
      Matrix M(dimension,dimension);
      m_cotan_matrix.swap(M);
    }
    CGAL::Vertex_around_face_iterator<TriangleMesh> vbegin, vend, vmiddle;

    BOOST_FOREACH(face_descriptor f, faces(tm)) {
      boost::tie(vbegin, vend) = vertices_around_face(halfedge(f,tm),tm);
      vertex_descriptor current = *(vbegin);
      vertex_descriptor neighbor_one = *(++vbegin);
      vertex_descriptor neighbor_two = *(++vbegin);
      Index i = get(vertex_id_map, current);
      Index j = get(vertex_id_map, neighbor_one);
      Index k = get(vertex_id_map, neighbor_two);
      Point_3 pi, pj, pk;

      VertexPointMap_reference p_i = get(vpm,current);
      VertexPointMap_reference p_j = get(vpm, neighbor_one);
      VertexPointMap_reference p_k = get(vpm, neighbor_two);
      pi = p_i;
      pj = p_j;
      pk = p_k;

      Vector_3 v_ij = construct_vector(p_i,p_j);
      Vector_3 v_ik = construct_vector(p_i,p_k);

      Vector_3 cross = cross_product(v_ij, v_ik);
      double dot = to_double(scalar_product(v_ij,v_ik));

      double norm_cross = (CGAL::sqrt(to_double(scalar_product(cross,cross))));

      double cotan_i = dot/norm_cross;
      m_cotan_matrix.add_coef(j,k ,-(1./2)*cotan_i);
      m_cotan_matrix.add_coef(k,j,-(1./2)* cotan_i);
      m_cotan_matrix.add_coef(j,j,(1./2)*cotan_i);
      m_cotan_matrix.add_coef(k,k,(1./2)* cotan_i);

      Vector_3 v_ji = construct_vector(p_j,p_i);
      Vector_3 v_jk = construct_vector(p_j,p_k);

      cross = cross_product(v_ji, v_jk);
      dot = to_double(scalar_product(v_ji, v_jk));
      double cotan_j = dot/norm_cross;
      m_cotan_matrix.add_coef(i,k ,-(1./2)*cotan_j);
      m_cotan_matrix.add_coef(k,i,-(1./2)* cotan_j);
      m_cotan_matrix.add_coef(i,i,(1./2)* cotan_j);
      m_cotan_matrix.add_coef(k,k,(1./2)* cotan_j);

      Vector_3 v_ki = construct_vector(p_k,p_i);
      Vector_3 v_kj = construct_vector(p_k,p_j);
      cross = cross_product(v_ki, v_kj);
      dot = to_double(scalar_product(v_ki,v_kj));
      double cotan_k = dot/norm_cross;
      m_cotan_matrix.add_coef(i,j,-(1./2)*cotan_k);
      m_cotan_matrix.add_coef(j,i,-(1./2)* cotan_k);
      m_cotan_matrix.add_coef(i,i,(1./2)* cotan_k);
      m_cotan_matrix.add_coef(j,j,(1./2)* cotan_k);

      //double area_face = CGAL::Polygon_mesh_processing::face_area(f,tm);
      //cross is 2*area
      m_mass_matrix.add_coef(i,i, (1./6.)*norm_cross);
      m_mass_matrix.add_coef(j,j, (1./6.)*norm_cross);
      m_mass_matrix.add_coef(k,k, (1./6.)*norm_cross);
      m_cotan_matrix.add_coef(i,i, 1e-8);
      m_cotan_matrix.add_coef(j,j, 1e-8);
      m_cotan_matrix.add_coef(k,k, 1e-8);
    }
    m_cotan_matrix.assemble_matrix();
    m_mass_matrix.assemble_matrix();

    m_time_step = 1./(num_edges(tm));
    m_time_step = m_time_step*summation_of_edges();
    m_time_step = m_time_step*m_time_step;
    update_kronecker_delta();
    factor_cotan_laplace();
    solve_cotan_laplace();
    compute_unit_gradient();
    compute_divergence();
    factor_phi();
  }

  LA la;
  LA la_cotan;
  int dimension;
  V2V<TriangleMesh> v2v;
  const TriangleMesh& tm;
  VertexPointMap vpm;
  std::set<vertex_descriptor> m_sources;
  double m_time_step;
  Matrix m_kronecker;
  Matrix m_mass_matrix, m_cotan_matrix;
  Vector m_solved_u;
  std::vector<Vector_3> m_X;
  Matrix m_index_divergence;
  Vector m_solved_phi;
  bool m_source_change_flag;

};

template <typename TriangleMesh,
          typename Traits,
          typename Mode,
          typename LA,
          typename VertexPointMap>
struct Base_helper
  : public Surface_mesh_geodesic_distances_3<TriangleMesh, Traits, LA, VertexPointMap>
{
  typedef Surface_mesh_geodesic_distances_3<TriangleMesh, Traits, LA, VertexPointMap> type;

  Base_helper(const TriangleMesh& tm, VertexPointMap vpm)
    : type(tm, vpm)
  {}

  Base_helper(const TriangleMesh& tm)
    : type(tm)
  {}

  type& base()
  {
    return static_cast<type&>(*this);
  }

  const type& base() const
  {
    return static_cast<const type&>(*this);
  }

  template <class VertexDistanceMap>
  void estimate_geodesic_distances(VertexDistanceMap vdm)
  {
    base().estimate_geodesic_distances(vdm);
  }
};

template<class TriangleMesh,
         class Traits,
         class VertexPointMap>
struct Idt_storage
{
  Intrinsic_Delaunay_triangulation_3<TriangleMesh, Traits> m_idt;

  Idt_storage(const TriangleMesh& tm, VertexPointMap vpm)
    : m_idt(tm, vpm)
  {}

  Idt_storage(const TriangleMesh& tm)
    : m_idt(tm)
  {}
};

template <typename TriangleMesh,
          typename Traits,
          typename LA,
          typename VertexPointMap>
struct Base_helper<TriangleMesh, Traits, Intrinsic_Delaunay, LA, VertexPointMap>
  : public Idt_storage<TriangleMesh, Traits, VertexPointMap>
  , public Surface_mesh_geodesic_distances_3<Intrinsic_Delaunay_triangulation_3<TriangleMesh, Traits>,
                                             Traits,
                                             LA,
                                             typename Intrinsic_Delaunay_triangulation_3<TriangleMesh, Traits>::Vertex_point_map>
{
  typedef CGAL::Heat_method_3::Intrinsic_Delaunay_triangulation_3<TriangleMesh, Traits> Idt;
  typedef Idt_storage<TriangleMesh, Traits, VertexPointMap> Idt_wrapper;

  typedef Surface_mesh_geodesic_distances_3<Idt, Traits, LA, typename Idt::Vertex_point_map> type;

  Base_helper(const TriangleMesh& tm, VertexPointMap vpm)
    : Idt_wrapper(tm, vpm)
    , type(this->m_idt)
  {}

  Base_helper(const TriangleMesh& tm)
    :  Idt_wrapper(tm)
    , type(this->m_idt)
  {}

  type& base()
  {
    return static_cast<type&>(*this);
  }

  const type& base() const
  {
    return static_cast<const type&>(*this);
  }

  template <class VertexDistanceMap>
  void estimate_geodesic_distances(VertexDistanceMap vdm)
  {
    base().estimate_geodesic_distances(this->m_idt.vertex_distance_map(vdm));
  }
};

} // namespace internal


/**
 * \ingroup PkgHeatMethod
 *
 * Class `Surface_mesh_geodesic_distances_3` computes estimated geodesic distances for a set of source vertices where sources can be added and removed.
 * The class performs a preprocessing step that only depends on the mesh, so that the distance computation takes less
 * time after changes to the set of sources.
 *
 * \tparam TriangleMesh a triangulated surface mesh, model of `FaceListGraph` and `HalfedgeListGraph`
 * \tparam Mode must be `Intrinsic_Delaunay` to indicate that an intrinsic Delaunay triangulation is internally constructed
 *                              or `Direct` to indicate that the input mesh should be used as is.
 *                              If `Intrinsic_Delaunay`, then the type `TriangleMesh` must have an internal property for `vertex_point`
 *                              and its value type must be the same as the value type of `VertexPointMap`.
 * \tparam VertexPointMap a model of `ReadablePropertyMap` with
 *         `boost::graph_traits<TriangleMesh>::%vertex_descriptor` as key and
 *         `Traits::Point_3` as value type.
 *         The default is `typename boost::property_map< TriangleMesh, vertex_point_t>::%const_type`.
 * \tparam LA a model of `SparseLinearAlgebraWithFactorTraits_d`. If `CGAL_EIGEN3_ENABLED` is defined,
 *         then `Eigen_solver_traits<Eigen::SimplicialLDLT<typename Eigen_sparse_matrix<double>::%EigenType > >`
 *         is used as default
 * \tparam Traits a model of `HeatMethodTraits_3`. The default is the Kernel of the value type
 *         of the vertex point map (extracted using `Kernel_traits`).
 *
 */
template <typename TriangleMesh,
          typename Mode = Direct,
          typename VertexPointMap = Default,
          typename LA = Default,
          typename Traits = Default >
class Surface_mesh_geodesic_distances_3
#ifndef DOXYGEN_RUNNING
  : public
    internal::Base_helper<
      TriangleMesh,
      typename Default::Get<Traits,
        typename Kernel_traits<
          typename boost::property_traits<
            typename Default::Get<
              VertexPointMap,
              typename boost::property_map< TriangleMesh, vertex_point_t>::const_type
            >::type
          >::value_type
        >::Kernel
      >::type,
      Mode,
      #ifdef CGAL_EIGEN3_ENABLED
      typename Default::Get<
        LA,
        Eigen_solver_traits<Eigen::SimplicialLDLT<typename Eigen_sparse_matrix<double>::EigenType > >
      >::type,
      #else
      LA,
      #endif
      typename Default::Get<
        VertexPointMap,
        typename boost::property_map< TriangleMesh, vertex_point_t>::const_type
      >::type
    >
#endif
{
  // extract real types from Default
  #ifdef CGAL_EIGEN3_ENABLED
  typedef typename Default::Get<
    LA,
    Eigen_solver_traits<Eigen::SimplicialLDLT<typename Eigen_sparse_matrix<double>::EigenType > >
  >::type LA_type;
  #else
  typedef LA LA_type;
  #endif

  typedef typename Default::Get<
    VertexPointMap,
    typename boost::property_map< TriangleMesh, vertex_point_t>::const_type>::type Vertex_point_map;
  typedef
    typename Default::Get<Traits,
                          typename Kernel_traits<
                            typename boost::property_traits<Vertex_point_map>::value_type
                          >::Kernel
    >::type Traits_type;

  typedef internal::Base_helper<TriangleMesh, Traits_type, Mode, LA_type, Vertex_point_map> Base_helper;

  const typename Base_helper::type& base() const
  {
    return Base_helper::base();
  }

  typename Base_helper::type& base()
  {
    return Base_helper::base();
  }

public:

  /// Vertex descriptor type
  typedef typename boost::graph_traits<TriangleMesh>::vertex_descriptor vertex_descriptor;
  #ifndef DOXYGEN_RUNNING
  /// Source vertex iterator type
  typedef typename Base_helper::type::Vertex_const_range Vertex_const_range;
  #else
  /// a model of `ConstRange` with an iterator that is model of `ForwardIterator` with `vertex_descriptor` as value type.
  typedef unspecified_type Vertex_const_range;
  /// Vertex point map type derived from `VertexPointMap`.
  typedef unspecified_type Vertex_point_map;
  #endif

  /*!
    \brief Constructor
  */
  Surface_mesh_geodesic_distances_3(const TriangleMesh& tm)
    : Base_helper(tm)
  {}

  /*!
    \brief Constructor
  */
  Surface_mesh_geodesic_distances_3(const TriangleMesh& tm, Vertex_point_map vpm)
    : Base_helper(tm, vpm)
  {}

  /**
   * returns the triangle mesh the algorithm is running on.
   */
  const TriangleMesh& triangle_mesh() const
  {
    return base().triangle_mesh();
  }

  /**
   * adds `vd` to the source set, returning `false` if `vd` is already in the set.
   */
  bool
  add_source(vertex_descriptor vd)
  {
    return base().add_source(vd);
  }

  /**
   * adds the range of vertices to the source set.
   * \tparam VertexConstRange a model of the concept `ConstRange` with value type `boost::graph_traits<TriangleMesh>::%vertex_descriptor`
   */
  template <typename VertexConstRange>
  void
  add_sources(const VertexConstRange& vrange)
  {
    base().add_sources(vrange);
  }

  /**
   * removes `vd` from the source set, returning `true` if `vd` was in the set.
   */
  bool
  remove_source(vertex_descriptor vd)
  {
    return base().remove_source(vd);
  }

  /**
   * clears the current source set.
   */
  void
  clear_sources()
  {
    base().clear_sources();
  }

  /**
   * get estimated distance from the current source set to a vertex `vd`.
   * \warning The return type is `double` even when used with an exact kernel.
   */
  double
  estimate_geodesic_distance(vertex_descriptor vd) const
  {
    return base().estimate_geodesic_distance(vd);
  }

  /**
   * returns the source set.
   */
  const Vertex_const_range&
  sources() const
  {
    return base().sources();
  }


  /**
   * fills the distance property map with the estimated geodesic distance of each vertex to the closest source vertex.
   * \tparam VertexDistanceMap a property map model of `WritablePropertyMap`
   * with `vertex_descriptor` as key type and `double` as value type.
   * \param vdm the vertex distance map to be filled
   * \warning The key type is `double` even when used with an exact kernel.
   **/
  template <class VertexDistanceMap>
  void estimate_geodesic_distances(VertexDistanceMap vdm)
  {
    Base_helper::estimate_geodesic_distances(vdm);
  }
};

#if defined(DOXYGEN_RUNNING) || defined(CGAL_EIGEN3_ENABLED)
/// \ingroup PkgHeatMethod
/// computes for each vertex of the triangle mesh `tm` the estimated geodesic distance to a given source vertex.
/// This function is provided only if \ref thirdpartyEigen "Eigen" 3.3 (or greater) is available and `CGAL_EIGEN3_ENABLED` is defined.
/// \tparam TriangleMesh a triangulated surface mesh, model of `FaceListGraph` and `HalfedgeListGraph`.
///         It must have an internal vertex point property map with the value type being a 3D point from a cgal Kernel model
/// \tparam VertexDistanceMap a property map model of `WritablePropertyMap`
///         with `boost::graph_traits<TriangleMesh>::%vertex_descriptor` as key type and `double` as value type.
/// \tparam Mode either the tag `Direct` or `Intrinsic_Delaunay`, which determines if the geodesic distance
///              is computed directly on the mesh or if the intrinsic Delaunay triangulation is applied first.
///              The default is `Intrinsic_Delaunay`.
/// \warning The return type is `double` even when used with an exact kernel.
///
/// \sa `CGAL::Heat_method_3::Surface_mesh_geodesic_distances_3`
template <typename TriangleMesh, typename VertexDistanceMap, typename Mode>
void
estimate_geodesic_distances(const TriangleMesh& tm,
                            VertexDistanceMap vdm,
                            typename boost::graph_traits<TriangleMesh>::vertex_descriptor source,
                            Mode)
{
  CGAL::Heat_method_3::Surface_mesh_geodesic_distances_3<TriangleMesh, Mode> hm(tm);
  hm.add_source(source);
  hm.estimate_geodesic_distances(vdm);
}


#ifndef DOXYGEN_RUNNING
template <typename TriangleMesh, typename VertexDistanceMap>
void
estimate_geodesic_distances(const TriangleMesh& tm,
                            VertexDistanceMap vdm,
                            typename boost::graph_traits<TriangleMesh>::vertex_descriptor source)
{
  CGAL::Heat_method_3::Surface_mesh_geodesic_distances_3<TriangleMesh, Intrinsic_Delaunay> hm(tm);
  hm.add_source(source);
  hm.estimate_geodesic_distances(vdm);
}
#endif


/// \ingroup PkgHeatMethod
/// computes for each vertex  of the triangle mesh `tm` the estimated geodesic distance to a given set of source vertices.
/// This function is provided only if \ref thirdpartyEigen "Eigen" 3.3 (or greater) is available and `CGAL_EIGEN3_ENABLED` is defined.
/// \tparam TriangleMesh a triangulated surface mesh, model of `FaceListGraph` and `HalfedgeListGraph`
///         It must have an internal vertex point property map with the value type being a 3D point from a cgal Kernel model
/// \tparam VertexDistanceMap a property map model of `WritablePropertyMap`
///         with `boost::graph_traits<TriangleMesh>::%vertex_descriptor` as key type and `double` as value type.
/// \tparam VertexConstRange a model of the concept `ConstRange` with value type `boost::graph_traits<TriangleMesh>::%vertex_descriptor`
/// \tparam Mode either the tag `Direct` or `Intrinsic_Delaunay`, which determines if the geodesic distance
///              is computed directly on the mesh or if the intrinsic Delaunay triangulation is applied first.
///              The default is `Intrinsic_Delaunay`.
/// \warning The return type is `double` even when used with an exact kernel.
/// \sa `CGAL::Heat_method_3::Surface_mesh_geodesic_distances_3`
template <typename TriangleMesh, typename VertexDistanceMap, typename VertexConstRange, typename Mode>
void
estimate_geodesic_distances(const TriangleMesh& tm,
                            VertexDistanceMap vdm,
                            const VertexConstRange& sources,
                            Mode
#ifndef DOXYGEN_RUNNING
                            , typename boost::enable_if<
                              typename boost::has_range_const_iterator<VertexConstRange>
                                     >::type* = 0
#endif
)
{
  CGAL::Heat_method_3::Surface_mesh_geodesic_distances_3<TriangleMesh, Mode> hm(tm);
  hm.add_sources(sources);
  hm.estimate_geodesic_distances(vdm);
}

#ifndef DOXYGEN_RUNNING
template <typename TriangleMesh, typename VertexDistanceMap, typename VertexConstRange>
void
estimate_geodesic_distances(const TriangleMesh& tm,
                            VertexDistanceMap vdm,
                            const VertexConstRange& sources,
                            typename boost::enable_if<
                              typename boost::has_range_const_iterator<VertexConstRange>
                                     >::type* = 0)
{
  CGAL::Heat_method_3::Surface_mesh_geodesic_distances_3<TriangleMesh, Intrinsic_Delaunay> hm(tm);
  hm.add_sources(sources);
  hm.estimate_geodesic_distances(vdm);
}
#endif

#endif

} // namespace Heat_method_3
} // namespace CGAL

#include <CGAL/enable_warnings.h>
#endif // CGAL_HEAT_METHOD_3_SURFACE_MESH_GEODESIC_DISTANCES_3_H
