Towards adding an example to TinyAD that uses CGAL::Surface_mesh

TinyAD4CGAL.cpp

#include <Eigen/Core>

#include <CGAL/Simple_cartesian.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/Surface_mesh_parameterization/IO/File_off.h>
#include <CGAL/Surface_mesh_parameterization/parameterize.h>
#include <CGAL/Surface_mesh_parameterization/Barycentric_mapping_parameterizer_3.h>
#include <CGAL/Polygon_mesh_processing/measure.h>

namespace TinyAD {
  template <typename T>
  inline Eigen::Index idx_from_handle(const CGAL::SM_Index<T>& _h)
  {
    return _h.idx();
  }

}

#include <TinyAD/ScalarFunction.hh>
#include <TinyAD/Utils/NewtonDirection.hh>
#include <TinyAD/Utils/NewtonDecrement.hh>
#include <TinyAD/Utils/LineSearch.hh>

typedef CGAL::Simple_cartesian<double> K;
typedef K::Point_3 Point_3;
typedef K::Point_2 Point_2;

typedef CGAL::Surface_mesh<K::Point_3> Mesh;
typedef boost::graph_traits<Mesh>::vertex_descriptor vertex_descriptor;
typedef boost::graph_traits<Mesh>::halfedge_descriptor halfedge_descriptor;
typedef boost::graph_traits<Mesh>::face_descriptor face_descriptor;

namespace SMP = CGAL::Surface_mesh_parameterization;


Eigen::Vector3d toeigen(const Point_3& p)
{
  return Eigen::Vector3d(p.x(),p.y(),p.z());
}

Eigen::Vector2d toeigen(const Point_2& p)
{
  return Eigen::Vector2d(p.x(),p.y());
}


int main(int argc, char* argv[])
{
  const std::string filename = (argc>1) ? argv[1] : CGAL::data_file_path("meshes/nefertiti.off");
  Mesh mesh;
  CGAL::IO::read_polygon_mesh(filename, mesh);

  // compute initial embedding

  // a halfedge on the border
  halfedge_descriptor bhd = CGAL::Polygon_mesh_processing::longest_border(mesh).first;

  // The UV property map that holds the parameterized values
  typedef Mesh::Property_map<vertex_descriptor, Point_2>  UV_pmap;
  UV_pmap uv_map = mesh.add_property_map<vertex_descriptor, Point_2>("h:uv").first;

  SMP::Barycentric_mapping_parameterizer_3<Mesh> tutte;

  SMP::parameterize(mesh, tutte, bhd, uv_map);

  // todo: Use an adapter so that the put of a 2D CGAL point
  // converts into a 2D Eigen point
  typedef Mesh::Property_map<vertex_descriptor, Eigen::Vector2d>  UVe_pmap;
  UVe_pmap ph_param = mesh.add_property_map<vertex_descriptor, Eigen::Vector2d>("h:uve").first;
  for(vertex_descriptor vd : vertices(mesh)){
    ph_param[vd] = toeigen(uv_map[vd]);
  }

  Mesh::Property_map<face_descriptor,Eigen::Matrix2d> ph_rest_shapes;
  bool created;
  boost::tie(ph_rest_shapes, created) = mesh.add_property_map<face_descriptor,Eigen::Matrix2d>("f:rest_shapes");
  assert(created);
  for (auto t : faces(mesh)){
    // Get 3D vertex positions
    Eigen::Vector3d ar_3d = toeigen(mesh.point(target(halfedge(t,mesh), mesh)));
    Eigen::Vector3d br_3d = toeigen(mesh.point(target(next(halfedge(t,mesh), mesh), mesh)));
    Eigen::Vector3d cr_3d = toeigen(mesh.point(source(halfedge(t,mesh), mesh)));

    // Set up local 2D coordinate system
    Eigen::Vector3d n = (br_3d - ar_3d).cross(cr_3d - ar_3d);
    Eigen::Vector3d b1 = (br_3d - ar_3d).normalized();
    Eigen::Vector3d b2 = n.cross(b1).normalized();

    // Express a, b, c in local 2D coordiante system
    Eigen::Vector2d ar_2d(0.0, 0.0);
    Eigen::Vector2d br_2d((br_3d - ar_3d).dot(b1), 0.0);
    Eigen::Vector2d cr_2d((cr_3d - ar_3d).dot(b1), (cr_3d - ar_3d).dot(b2));

    // Save 2-by-2 matrix with edge vectors as colums
    ph_rest_shapes[t]= TinyAD::col_mat(br_2d - ar_2d, cr_2d - ar_2d);
  }

  // Set up function with 2D vertex positions as variables.
  auto func = TinyAD::scalar_function<2>(mesh.vertices());

  // Add objective term per triangle. Each connecting 3 vertices.
  func.add_elements<3>(mesh.faces(), [&] (auto& element) -> TINYAD_SCALAR_TYPE(element)
  {
    // Element is evaluated with either double or TinyAD::Double<6>
    using T = TINYAD_SCALAR_TYPE(element);

    // Get variable 2D vertex positions of triangle t
    face_descriptor t = element.handle;
    Eigen::Vector2<T> a = element.variables(target(halfedge(t,mesh),mesh));
    Eigen::Vector2<T> b = element.variables(target(next(halfedge(t,mesh),mesh), mesh));
    Eigen::Vector2<T> c = element.variables(source(halfedge(t,mesh),mesh));

    // Triangle flipped?
    Eigen::Matrix2<T> M = TinyAD::col_mat(b - a, c - a);
    if (M.determinant() <= 0.0)
      return (T)INFINITY;

    // Get constant 2D rest shape and area of triangle t
    Eigen::Matrix2d Mr = ph_rest_shapes[t];
    double A = 0.5 * Mr.determinant();

    // Compute symmetric Dirichlet energy
    Eigen::Matrix2<T> J = M * Mr.inverse();
    return A * (J.squaredNorm() + J.inverse().squaredNorm());
  });

    // Assemble inital x vector from parametrization property.
    // x_from_data(...) takes a lambda function that maps
    // each variable handle (OpenMesh::SmartVertexHandle) to its initial 2D value (Eigen::Vector2d).
    Eigen::VectorXd x = func.x_from_data([&] (vertex_descriptor v) {
        return ph_param[v];
    });

    // Projected Newton
    TinyAD::LinearSolver solver;
    int max_iters = 1000;
    double convergence_eps = 1e-2;
    for (int i = 0; i < max_iters; ++i)
    {
        auto [f, g, H_proj] = func.eval_with_hessian_proj(x);
        TINYAD_DEBUG_OUT("Energy in iteration " << i << ": " << f);
        Eigen::VectorXd d = TinyAD::newton_direction(g, H_proj, solver);
        if (TinyAD::newton_decrement(d, g) < convergence_eps)
            break;
        x = TinyAD::line_search(x, d, f, g, func);
    }
    TINYAD_DEBUG_OUT("Final energy: " << func.eval(x));

    // Write final x vector to parametrization property.
    // x_to_data(...) takes a lambda function that writes the final value
    // of each variable (Eigen::Vector2d) back our parametrization property.
    func.x_to_data(x, [&] (vertex_descriptor v, const Eigen::Vector2d& _p) {
        ph_param[v] = _p;
    });
  return 0;
}