From c62e503a390dd5595a94c865951925b74a19a60a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Wed, 23 Nov 2016 00:37:12 +0100 Subject: [PATCH] Added Orbital Tutte parameterizer and an example Only the orbifold type I method with MVC coefs is implented in this commit --- BGL/include/CGAL/boost/graph/Seam_mesh.h | 51 +- .../CMakeLists.txt | 14 +- .../Surface_mesh_parameterization/orbital.cpp | 223 ++++++++ .../Orbital_Tutte_parameterizer_3.h | 535 ++++++++++++++++++ .../internal/shortest_path.h | 167 ++++++ 5 files changed, 970 insertions(+), 20 deletions(-) create mode 100644 Surface_mesh_parameterization/examples/Surface_mesh_parameterization/orbital.cpp create mode 100644 Surface_mesh_parameterization/include/CGAL/Surface_mesh_parameterization/Orbital_Tutte_parameterizer_3.h create mode 100644 Surface_mesh_parameterization/include/CGAL/Surface_mesh_parameterization/internal/shortest_path.h diff --git a/BGL/include/CGAL/boost/graph/Seam_mesh.h b/BGL/include/CGAL/boost/graph/Seam_mesh.h index b03fd7a23d8..0ae90dcc4bc 100644 --- a/BGL/include/CGAL/boost/graph/Seam_mesh.h +++ b/BGL/include/CGAL/boost/graph/Seam_mesh.h @@ -875,6 +875,34 @@ public: /// @} + bool add_seam(TM_vertex_descriptor tmvd_s, TM_vertex_descriptor tmvd_t) + { + std::pair tmed = CGAL::edge(tmvd_s, tmvd_t, tm); + if(!tmed.second) { + std::cout << "was ignored because it is not a valid edge of the mesh (Warning!)" << std::endl; + return false; + } + + if(!is_border(tmed.first, tm)) { // ignore seams that are also a border edge + if(get(sem, tmed.first) == true) { + std::cout << "was ignored because it is already marked as a seam (Warning!)" << std::endl; + return false; + } + + std::cout << "was added to the seam mesh" << std::endl; + + put(sem, tmed.first, true); + put(svm, tmvd_s, true); + put(svm, tmvd_t, true); + ++number_of_seams; + } else { + std::cout << "was ignored because it is on the border of the mesh (Warning!)" << std::endl; + return false; + } + + return true; + } + /// adds the seams to the property maps of the seam mesh. /// /// In input, a seam edge is described by the pair of integers that correspond @@ -926,27 +954,12 @@ public: assert(s < tmvds.size() && t < tmvds.size()); TM_vertex_descriptor tmvd_s = tmvds[s], tmvd_t = tmvds[t]; - std::pair tmed = CGAL::edge(tmvd_s, tmvd_t, tm); - if(!tmed.second) { - std::cout << "Warning: the input seam " << s << " " << t << " was ignored"; - std::cout << " because it is not a valid edge of the mesh" << std::endl; + std::cout << "Seam " << s << " " << t << " "; + if(!add_seam(tmvd_s, tmvd_t)) continue; - } - if(!is_border(tmed.first, tm)) { // ignore seams that are also a border edge - std::cout << "Adding seam: " << s << " " << t << std::endl; - ++number_of_seams; - - put(sem, tmed.first, true); - put(svm, tmvd_s, true); - put(svm, tmvd_t, true); - - if(tmhd == boost::graph_traits::null_halfedge()) { - tmhd = CGAL::halfedge(CGAL::edge(tmvd_s, tmvd_t, tm).first, tm); - } - } else { - std::cout << "Warning: the input seam " << s << " " << t << " was ignored"; - std::cout << " because it is on the border of the mesh" << std::endl; + if(tmhd == boost::graph_traits::null_halfedge()) { + tmhd = CGAL::halfedge(CGAL::edge(tmvd_s, tmvd_t, tm).first, tm); } } diff --git a/Surface_mesh_parameterization/examples/Surface_mesh_parameterization/CMakeLists.txt b/Surface_mesh_parameterization/examples/Surface_mesh_parameterization/CMakeLists.txt index 71ee12d264e..cbc1aca87b8 100644 --- a/Surface_mesh_parameterization/examples/Surface_mesh_parameterization/CMakeLists.txt +++ b/Surface_mesh_parameterization/examples/Surface_mesh_parameterization/CMakeLists.txt @@ -22,12 +22,24 @@ if ( CGAL_FOUND ) # Executables that require Eigen 3.1 include( ${EIGEN3_USE_FILE} ) + set(CMAKE_MODULE_PATH "${EIGEN3_INCLUDE_DIR}/cmake/;${CMAKE_MODULE_PATH}") + + find_package(SPQR) + if (NOT SPQR_FOUND) + message(STATUS "NOTICE: External sparse solver libraries not found") + endif() + + include_directories(${SPQR_INCLUDES}) + create_single_source_cgal_program( "discrete_authalic.cpp" ) create_single_source_cgal_program( "lscm.cpp" ) create_single_source_cgal_program( "seam_Polyhedron_3.cpp" ) create_single_source_cgal_program( "simple_parameterization.cpp" ) create_single_source_cgal_program( "square_border_parameterizer.cpp" ) - + create_single_source_cgal_program( "orbital.cpp" ) + + target_link_libraries(orbital ${SPQR_LIBRARIES}) + else(EIGEN3_FOUND) message(STATUS "NOTICE: The examples require Eigen 3.1 (or greater) and will not be compiled.") endif(EIGEN3_FOUND) diff --git a/Surface_mesh_parameterization/examples/Surface_mesh_parameterization/orbital.cpp b/Surface_mesh_parameterization/examples/Surface_mesh_parameterization/orbital.cpp new file mode 100644 index 00000000000..3f70825a292 --- /dev/null +++ b/Surface_mesh_parameterization/examples/Surface_mesh_parameterization/orbital.cpp @@ -0,0 +1,223 @@ +#include + +#include +#include + +#include +#include + +#include +#include +#include + +#include + +#include +#include +#include + +#include +#include + +#include +#include +#include +#include +#include +#include + +typedef CGAL::Simple_cartesian Kernel; +typedef Kernel::Point_2 Point_2; +typedef Kernel::Point_3 Point_3; +typedef CGAL::Surface_mesh SurfaceMesh; + +typedef boost::graph_traits::vertex_descriptor SM_vertex_descriptor; +typedef boost::graph_traits::halfedge_descriptor SM_halfedge_descriptor; +typedef boost::graph_traits::edge_descriptor SM_edge_descriptor; + +typedef SurfaceMesh::Property_map Seam_edge_pmap; +typedef SurfaceMesh::Property_map Seam_vertex_pmap; + +typedef CGAL::Seam_mesh Mesh; +typedef boost::graph_traits::vertex_descriptor vertex_descriptor; +typedef boost::graph_traits::halfedge_descriptor halfedge_descriptor; + +typedef SurfaceMesh::Property_map UV_pmap; + +namespace SMP = CGAL::Surface_mesh_parameterization; + +/// Read the cones from the input file. +SMP::Error_code read_cones(const SurfaceMesh& sm, const char* filename, + std::vector& cone_vds_in_sm) +{ + std::ifstream in(filename); + std::string vertices_line; + std::getline(in, vertices_line); // read the first line of the file + std::istringstream iss(vertices_line); + int index_1, index_2, index_3; + if(!(iss >> index_1 >> index_2 >> index_3)) { + std::cerr << "Error: problem loading the input cones" << std::endl; + return SMP::ERROR_WRONG_PARAMETER; + } + + std::cout << "Cones: " << index_1 << " " << index_2 << " " << index_3 << std::endl; + + // Locate the cones in the underlying mesh 'sm' + SM_vertex_descriptor vd_1, vd_2, vd_3; + + int counter = 0; + BOOST_FOREACH(SM_vertex_descriptor vd, vertices(sm)) { + if(counter == index_1) { + vd_1 = vd; + } else if(counter == index_2) { + vd_2 = vd; + } else if(counter == index_3) { + vd_3 = vd; + } + ++counter; + } + + CGAL_postcondition(vd_1 != SM_vertex_descriptor() && + vd_2 != SM_vertex_descriptor() && + vd_3 != SM_vertex_descriptor()); + + cone_vds_in_sm[0] = vd_1; + cone_vds_in_sm[1] = vd_2; + cone_vds_in_sm[2] = vd_3; + + return SMP::OK; +} + +/// Locate the cones on the seam mesh (find the corresponding seam mesh +/// vertex_descriptor) and mark them with a tag that indicates whether it is a +/// simple cone or a duplicated cone. +template +void locate_cones(const Mesh& mesh, + const std::vector& cone_vds_in_sm, + ConeMap& cones) +{ + // property map to go from SM_vertex_descriptor to Point_3 + typedef SMP::internal::Kernel_traits::PPM SM_PPM; + const SM_PPM sm_ppmap = get(boost::vertex_point, mesh.mesh()); + + // property map to go from vertex_descriptor to Point_3 + typedef SMP::internal::Kernel_traits::PPM PPM; + const PPM ppmap = get(boost::vertex_point, mesh); + + // to know the ID of the vertex in the seam mesh (debug) + int counter = 0; + + // to check that the duplicated vertex is correctly seen twice + // TMP till a proper check is made + int is_doubled_at_v3 = 0; + + // the cones in the underlying mesh + CGAL_precondition(cone_vds_in_sm.size() == 3); + SM_vertex_descriptor vd_1 = cone_vds_in_sm[0]; + SM_vertex_descriptor vd_2 = cone_vds_in_sm[1]; + SM_vertex_descriptor vd_3 = cone_vds_in_sm[2]; + + BOOST_FOREACH(vertex_descriptor vd, vertices(mesh)) { + if(get(ppmap, vd) == get(sm_ppmap, vd_1)) { + cones.insert(std::make_pair(vd, SMP::Unique_cone)); + std::cout << counter << " [eq to ind1]" << std::endl; + } else if(get(ppmap, vd) == get(sm_ppmap, vd_2)) { + std::cout << counter << " [eq to ind2]" << std::endl; + cones.insert(std::make_pair(vd, SMP::Unique_cone)); + } else if(get(ppmap, vd) == get(sm_ppmap, vd_3)) { + std::cout << counter << " [eq to ind3]" << std::endl; + is_doubled_at_v3++; + cones.insert(std::make_pair(vd, SMP::Duplicated_cone)); + } + ++counter; + } + + std::cout << cone_vds_in_sm.size() << " in sm" << std::endl; + std::cout << cones.size() << " cones" << std::endl; + + if(is_doubled_at_v3 != 2) + exit(0); + + CGAL_postcondition(cones.size() == 4 && is_doubled_at_v3 == 2); +} + +int main(int argc, char * argv[]) +{ + CGAL::Timer task_timer; + task_timer.start(); + + SurfaceMesh sm; // underlying mesh of the seam mesh + + const char* mesh_filename = (argc>1) ? argv[1] : "../data/bunny.off"; + std::ifstream in_mesh(mesh_filename); + if(!in_mesh) { + std::cerr << "Error: problem loading the input data" << std::endl; + return 1; + } + in_mesh >> sm; + + // Selection file that contains the cones and possibly the path between cones + // -- the first line for the cones indices + // -- the second line must be empty + // -- the third line optionally provides the seam edges indices as 'e11 e12 e21 e22 e31 e32' etc. + const char* cone_filename = (argc>2) ? argv[2] : "../data/bunny.selection.txt"; + + // Read the cones and find the corresponding vertex_descriptor in the underlying mesh 'sm' + std::vector cone_vds_in_sm(3); + read_cones(sm, cone_filename, cone_vds_in_sm); + + // Two property maps to store the seam edges and vertices + Seam_edge_pmap seam_edge_pm = sm.add_property_map("e:on_seam", false).first; + Seam_vertex_pmap seam_vertex_pm = sm.add_property_map("v:on_seam",false).first; + + // The seam mesh + Mesh mesh(sm, seam_edge_pm, seam_vertex_pm); + + // Use the path provided between cones to create a seam mesh + SM_halfedge_descriptor smhd = mesh.add_seams(cone_filename); + if(smhd == SM_halfedge_descriptor() ) { + std::cout << "No seams were given in input, computing shortest paths between cones" << std::endl; + std::list seam_edges; + SMP::internal::compute_shortest_paths_between_cones(sm, cone_vds_in_sm, seam_edges); + + // Add the seams to the seam mesh + BOOST_FOREACH(SM_edge_descriptor e, seam_edges) { + std::cout << "Seam " << source(e, sm) << " " << target(e, sm) << " "; + mesh.add_seam(source(e, sm), target(e, sm)); + } + } + + std::cout << mesh.number_of_seam_edges() << " seam edges" << std::endl; + + // Index map of the seam mesh (assuming a single connected component so far) + typedef boost::unordered_map Indices; + Indices indices; + boost::associative_property_map vimap(indices); + int counter = 0; + BOOST_FOREACH(vertex_descriptor vd, vertices(mesh)) { + put(vimap, vd, counter++); + } + + // Mark the cones in the seam mesh + typedef boost::unordered_map Cones; + Cones cmap; + locate_cones(mesh, cone_vds_in_sm, cmap); + + // The 2D points of the uv parametrisation will be written into this map + // Note that this is a halfedge property map, and that uv values + // are only stored for the canonical halfedges representing a vertex + UV_pmap uvmap = sm.add_property_map("h:uv").first; + + // Parameterizer + typedef SMP::Orbital_Tutte_parameterizer_3 Parameterizer; + Parameterizer parameterizer; + + // a halfedge on the (possibly virtual) border + // only used in output (will also be used to handle multiple connected components in the future) + halfedge_descriptor bhd = CGAL::Polygon_mesh_processing::longest_border(mesh, + CGAL::Polygon_mesh_processing::parameters::all_default()).first; + + parameterizer.parameterize(mesh, bhd, cmap, uvmap, vimap); + + std::cout << "Finished in " << task_timer.time() << " seconds" << std::endl; +} diff --git a/Surface_mesh_parameterization/include/CGAL/Surface_mesh_parameterization/Orbital_Tutte_parameterizer_3.h b/Surface_mesh_parameterization/include/CGAL/Surface_mesh_parameterization/Orbital_Tutte_parameterizer_3.h new file mode 100644 index 00000000000..005f23bedca --- /dev/null +++ b/Surface_mesh_parameterization/include/CGAL/Surface_mesh_parameterization/Orbital_Tutte_parameterizer_3.h @@ -0,0 +1,535 @@ +// Copyright (c) 2016 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$ +// $Id$ +// +// +// Author(s) : + +#ifndef CGAL_SURFACE_MESH_PARAMETERIZATION_ORBITAL_TUTTE_PARAMETERIZER_3_H +#define CGAL_SURFACE_MESH_PARAMETERIZATION_ORBITAL_TUTTE_PARAMETERIZER_3_H + +#include +#include +#include + +#include + +#include +#include +#include + +#include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +/// \file Orbital_Tutte_parameterizer_3.h + +// @todo checks that cones are different, are on seams, seam is one connected +// component + +namespace CGAL { + +namespace Surface_mesh_parameterization { + +enum Cone_type +{ + Unique_cone, + Duplicated_cone +}; + +enum Orbifold_type +{ + Square, + Diamond, + Triangle, + Parallelogram +}; + +enum Weight_type +{ + Cotan, + Mvc +}; + +template +< + typename TriangleMesh, + typename SparseLinearAlgebraTraits_d + // ~1s (for bear.off) + = Eigen_solver_traits::EigenType> > + + // ~15s +// = Eigen_solver_traits::EigenType> > + + // ~20s +// = Eigen_solver_traits< > + + // ~35s +// = Eigen_solver_traits::EigenType> > + + // Can't use SuperLU / UmfPackLU, and SparseQR / SparseLU are way too slow +> +class Orbital_Tutte_parameterizer_3 +{ +private: + typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; + typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; + typedef typename boost::graph_traits::face_descriptor face_descriptor; + + typedef typename boost::graph_traits::vertex_iterator vertex_iterator; + typedef typename boost::graph_traits::face_iterator face_iterator; + + // SparseLinearAlgebraTraits_d subtypes: + typedef SparseLinearAlgebraTraits_d Sparse_LA; + typedef typename Sparse_LA::Vector Vector; + typedef typename Sparse_LA::Matrix Matrix; + + // Kernel subtypes + typedef typename internal::Kernel_traits::Kernel Kernel; + typedef typename internal::Kernel_traits::PPM PPM; + typedef typename Kernel::FT NT; + typedef typename Kernel::Vector_2 Vector_2; + typedef typename Kernel::Vector_3 Vector_3; + typedef typename Kernel::Point_2 Point_2; + typedef typename Kernel::Point_3 Point_3; + +public: + // Linear system + /// adds a positional constraint on a vertex x_ind, so that x_ind*w=rhs + void addConstraint(Matrix& A, Vector& B, int ind, double w, Vector_2 rhs) const + { + std::cout << "Constraining " << ind << std::endl; + + A.set_coef(2*ind, 2*ind, w, true /*new_coeff*/); + A.set_coef(2*ind + 1, 2*ind + 1, w, true /*new_coeff*/); + + B[2*ind] = rhs[0]; + B[2*ind + 1] = rhs[1]; + } + + /// adds constraints so that T * x_sinds = x_tinds, where T is a 2x2 + /// matrix, and the Transformation T is modified to affine from + /// linear by requiring that T * x_si - x_ti = T * x_s1 - x_t1 + void addTransConstraints(int s0, int s, int t, + const Eigen::Matrix2d& T, + Matrix& A, Vector& B) const + { + // Matlab lines are commented for comparison. + // Matlab fills together 2*x-1 and 2*x, but C++ fills 2*x and 2*x+1, + // as everything (including loops!) starts at 0 and not 1. + + int t0 = s0; // for clarity + + // iterate on both rows ot the 2x2 matrix T + for(int vert_ind=0; vert_ind<2; ++vert_ind) { + // building up the equations by summing up the terms + + // + // obj.A(end+1, 2*sinds(ind)+[-1,0]) = T(vert_ind,:); + A.set_coef(2*s + vert_ind, 2*s, T(vert_ind, 0), true /*new_coeff*/); + A.set_coef(2*s + vert_ind, 2*s + 1, T(vert_ind, 1), true /*new_coeff*/); + + // - + // obj.A(end, 2*sinds(1)+[-1,0]) = obj.A(end, 2*sinds(1)+[-1,0]) - T(vert_ind,:); + A.add_coef(2*s + vert_ind, 2*s0, - T(vert_ind, 0)); + A.add_coef(2*s + vert_ind, 2*s0 + 1, -T(vert_ind, 1)); + + // - x_ti + // obj.A(end, 2*tinds(ind)+vert_ind-2) = obj.A(end, 2*tinds(ind)+vert_ind-2)-1; + A.add_coef(2*s + vert_ind, 2*t + vert_ind, -1); + + // + x_t1 + // obj.A(end, 2*tinds(1)+vert_ind-2) = obj.A(end, 2*tinds(1)+vert_ind-2)+1; + A.add_coef(2*s + vert_ind, 2*t0 + vert_ind, 1); + + // left hand side is zero + // obj.b=[obj.b; 0]; + B[2*s + vert_ind] = 0; + } + } + + /// Compute the rotational constraint on the border of the mesh. + /// Cone constraints are also added. + template + void AddRotationalConstraint(const TriangleMesh& mesh, + const ConeMap& cmap, + VertexIndexMap vimap, + Matrix& A, Vector& B) const + { + // positions of the cones in the plane TMP + std::vector tcoords(2); + tcoords[0] = Vector_2(-1, -1); + tcoords[1] = Vector_2(1, 1); + + // rotations of the seams TMP + // angles at the cones + std::vector angs(2); + angs[0] = 4; // singularity is [4; 4] for orb1 + angs[1] = 4; // singularity is [4; 4] for orb1 + + // Go through the seam + boost::unordered_set constrained_cones; + for(int i=0; i<2; ++i) { // TMP + // Find a non-duplicated cone that has not yet been constrained + vertex_descriptor start_cone; + int start_cone_index = -1; + + BOOST_FOREACH(vertex_descriptor vd, vertices(mesh)) { + typename ConeMap::const_iterator it = cmap.find(vd); + if(it != cmap.end() && it->second == Unique_cone) { + start_cone_index = get(vimap, vd); + if(constrained_cones.find(start_cone_index) == constrained_cones.end() ) { + start_cone = vd; + constrained_cones.insert(start_cone_index); + break; + } + } + } + CGAL_postcondition(start_cone != vertex_descriptor() && start_cone_index != -1); + + std::cout << "starting from " << start_cone_index << std::endl; + + // Mark 'start_cone' as constraint + addConstraint(A, B, start_cone_index, 1. /*entry in A*/, tcoords[i]); + + // Go through the seam, marking rotation and cone constraints + halfedge_descriptor hd = halfedge(start_cone, mesh); + CGAL_precondition(mesh.has_on_seam(hd)); + halfedge_descriptor bhd = opposite(hd, mesh); + CGAL_precondition(is_border(bhd, mesh)); + + // the rotation angle + double ang = angs[i]; + + // the rotation matrix according to the angle 'ang' + Eigen::Matrix2d R; + R(0,0) = std::cos(2 * CGAL_PI / ang); R(0,1) = - std::sin(2 * CGAL_PI / ang); + R(1,0) = std::sin(2 * CGAL_PI / ang); R(1,1) = std::cos(2 * CGAL_PI / ang); + + // move through the seam from start_cone_index till we reach a duplicated cone + while(true) { + // Get the two halfedges on each side of the stream + halfedge_descriptor hd1 = bhd; // only for clarity + + // the non-border halfedge with same vertices (in the underlying mesh of the seam + // mesh) as bhd is simply bhd with the 'seam' boolean set to false + halfedge_descriptor hd2 = bhd; + hd2.seam = false; + + // Add the rotational constraint between the two halfedges on the seam + vertex_descriptor hd1_target = target(hd1, mesh); + vertex_descriptor hd2_target = target(hd2, mesh); + int hd1t_index = get(vimap, hd1_target); + int hd2t_index = get(vimap, hd2_target); + + std::cout << "hd1/hd2: " << hd1t_index << " " << hd2t_index << std::endl; + + addTransConstraints(start_cone_index, hd1t_index, hd2t_index, R, A, B); + + // Check if we have reached the duplicated cone + vertex_descriptor bhd_target = target(bhd, mesh); // also known as 'hd1_target' + typename ConeMap::const_iterator is_in_map = cmap.find(bhd_target); + if(is_in_map != cmap.end()) { + // starting from a unique cone, the next cone must be the duplicated cone + CGAL_assertion(is_in_map->second == Duplicated_cone); + break; + } + + // move to the next halfedge couple (walking on the border of the seam) + bhd = next(bhd, mesh); + CGAL_postcondition(mesh.has_on_seam(bhd) && is_border(bhd, mesh)); + } + } + } + + // MVC computations + NT compute_w_ij_mvc(const Point_3& pi, const Point_3& pj, const Point_3& pk) const + { + // -> -> + // Compute the angle (pj, pi, pk), the angle between the vectors ij and ik + NT angle = internal::compute_angle_rad(pj, pi, pk); + NT weight = std::tan(0.5 * angle); + + return weight; + } + + /// Compute the coefficients of the mean value Laplacian matrix for the edge + /// `ij` in the face `ijk` + void fill_mvc_matrix(const Point_3& pi, int i, + const Point_3& pj, int j, + const Point_3& pk, int k, Matrix& L) const + { +// std::cout << "compute mvc coef at: " << i << " " << j << " " << k << std::endl; + + // For MVC, the entry of L(i,j) is - [ tan(gamma_ij/2) + tan(delta_ij)/2 ] / |ij| + // where gamma_ij and delta_ij are the angles at i around the edge ij + + // This function computes the angle alpha at i, and add + // -- A(i,j) += tan(alpha / 2) / |ij| + // -- A(i,k) += tan(alpha / 2) / |ik| + // -- A(i,i) -= A(i,j) + A(i,k) + + // The other parts of A(i,j) and A(i,k) will be added when this function + // is called from the neighboring faces of F_ijk that share the vertex i + + // Compute: - tan(alpha / 2) + NT w_i_base = -1.0 * compute_w_ij_mvc(pi, pj, pk); + + // @fixme unefficient: lengths are computed (and inversed!) twice per edge + + // Set w_ij in matrix + Vector_3 edge_ij = pi - pj; + NT len_ij = CGAL::sqrt(edge_ij * edge_ij); + CGAL_assertion(len_ij != 0.0); // two points are identical! + NT w_ij = w_i_base / len_ij; + L.add_coef(2*i, 2*j, w_ij); + L.add_coef(2*i +1, 2*j + 1, w_ij); + + // Set w_ik in matrix + Vector_3 edge_ik = pi - pk; + NT len_ik = CGAL::sqrt(edge_ik * edge_ik); + CGAL_assertion(len_ik != 0.0); // two points are identical! + NT w_ik = w_i_base / len_ik; + L.add_coef(2*i, 2*k, w_ik); + L.add_coef(2*i + 1, 2*k + 1, w_ik); + + // Add to w_ii (w_ii = - sum w_ij) + NT w_ii = - w_ij - w_ik; + L.add_coef(2*i, 2*i, w_ii); + L.add_coef(2*i + 1, 2*i + 1, w_ii); + } + + /// Compute the mean value Laplacian matrix. + template + void mean_value_laplacian(const TriangleMesh& mesh, + VertexIndexMap vimap, + Matrix& L) const + { + const PPM ppmap = get(vertex_point, mesh); + + BOOST_FOREACH(face_descriptor fd, faces(mesh)) { + halfedge_descriptor hd = halfedge(fd, mesh); + + vertex_descriptor vd_i = target(hd, mesh); + vertex_descriptor vd_j = target(next(hd, mesh), mesh); + vertex_descriptor vd_k = source(hd, mesh); + const Point_3& pi = get(ppmap, vd_i); + const Point_3& pj = get(ppmap, vd_j); + const Point_3& pk = get(ppmap, vd_k); + int i = get(vimap, vd_i); + int j = get(vimap, vd_j); + int k = get(vimap, vd_k); + + fill_mvc_matrix(pi, i, pj, j, pk, k, L); + fill_mvc_matrix(pj, j, pk, k, pi, i, L); + fill_mvc_matrix(pk, k, pi, i, pj, j, L); + } + } + + /// Copy the solution into the UV property map. + template + void assign_solution(const TriangleMesh& mesh, + const Vector& X, + VertexUVMap uvmap, + const VertexIndexMap vimap) const + { + std::cout << "size of X: " << X.size() << std::endl; + CGAL_assertion(X.size() == 2*num_vertices(mesh) ); + + BOOST_FOREACH(vertex_descriptor vd, vertices(mesh)) { + int index = get(vimap, vd); + NT u = X(2*index); + NT v = X(2*index + 1); + + put(uvmap, vd, Point_2(u, v)); + } + } + + /// Solve the linear system. + template + Error_code computeFlattening(const TriangleMesh& mesh, + const Matrix& A, const Vector& B, + const Matrix& L, + VertexUVMap uvmap, VertexIndexMap vimap) const + { + std::size_t n = B.size(); + + // Fill the large system M.Xf = Bf: + // ( L A' ) ( Xf ) = ( B ) + // ( A 0 ) ( Xf ) = ( 0 ) + Matrix M(2*n, 2*n); + Vector Bf(2*n); + NT D; + Vector Xf(2*n); + + // full B + for(std::size_t i=0; i::InnerIterator + it(L.eigen_object(), k); it; ++it) { + M.set_coef(it.row(), it.col(), it.value(), true /*new_coeff*/); +// std::cout << it.row() << " " << it.col() << " " << it.value() << '\n'; + } + } + + std::cout << "filled topleft of M" << std::endl; + + for(int k=0; k::InnerIterator + it(A.eigen_object(), k); it; ++it) { + M.set_coef(it.col(), it.row() + n, it.value(), true /*new_coeff*/); // A + M.set_coef(it.row() + n, it.col(), it.value(), true /*new_coeff*/); // A' +// std::cout << it.row() << " " << it.col() << " " << it.value() << '\n'; + } + } + + std::cout << "Filled M and Bf" << std::endl; + +#ifdef SMP_OUTPUT_ORBITAL_MATRICES + std::ofstream outM("linears_M.txt"); + outM << M.eigen_object() << std::endl; +#endif + + CGAL::Timer task_timer; + task_timer.start(); + + SparseLinearAlgebraTraits_d solver; + if(!solver.linear_solver(M, Bf, Xf, D)) { + std::cout << "Could not solve linear system" << std::endl; + return ERROR_CANNOT_SOLVE_LINEAR_SYSTEM; + } + CGAL_assertion(D == 1.0); + std::cout << "Solved the linear system in " << task_timer.time() << " seconds" << std::endl; + + Vector X(n); + for(std::size_t i=0; i + Error_code parameterize(const TriangleMesh& mesh, + halfedge_descriptor bhd, + ConeMap cmap, + VertexUVMap uvmap, + VertexIndexMap vimap) const + { + std::cout << "Flattening" << std::endl; + Error_code status; + + // %%%%%%%%%%%%%%%%%%%%%%% + // Boundary conditions + // %%%%%%%%%%%%%%%%%%%%%%% + std::size_t nbVertices = num_vertices(mesh); + Matrix A(2 * nbVertices, 2 * nbVertices); + Vector B(2 * nbVertices); + + // add rotational constraints + AddRotationalConstraint(mesh, cmap, vimap, A, B); + +#ifdef SMP_OUTPUT_ORBITAL_MATRICES + std::cout << "A and B are filled" << std::endl; + std::ofstream outA("A.txt"), outB("B.txt"); + + for(int k=0; k::InnerIterator it(A.eigen_object(), k); it; ++it) { + outA << "(" << it.row() << ", " << it.col() << ") " << it.value() << std::endl; + } + } + outA << std::endl << A.eigen_object() << std::endl << std::endl; + + outB << B << std::endl; +#endif + + // %%%%%%%%%%%%%%%%%%%% + // Energy (Laplacian) + // %%%%%%%%%%%%%%%%%%%% + + Matrix L(2*nbVertices, 2*nbVertices); + mean_value_laplacian(mesh, vimap, L); + +#ifdef SMP_OUTPUT_ORBITAL_MATRICES + std::ofstream outL("MVCc.txt"); + outL << L.eigen_object() << std::endl; +#endif + + // compute the flattening by solving the boundary conditions + // while satisfying the convex combination property with L + std::cout << "solving system" << std::endl; + status = computeFlattening(mesh, A, B, L, uvmap, vimap); + if(status != OK) + return status; + + std::ofstream out("orbital_result.off"); + IO::output_uvmap_to_off(mesh, bhd, uvmap, out); + + return OK; + } +}; + +} // namespace Surface_mesh_parameterization + +} // namespace CGAL + +#endif // CGAL_SURFACE_MESH_PARAMETERIZATION_ORBITAL_TUTTE_PARAMETERIZER_3_H diff --git a/Surface_mesh_parameterization/include/CGAL/Surface_mesh_parameterization/internal/shortest_path.h b/Surface_mesh_parameterization/include/CGAL/Surface_mesh_parameterization/internal/shortest_path.h new file mode 100644 index 00000000000..e9105501da1 --- /dev/null +++ b/Surface_mesh_parameterization/include/CGAL/Surface_mesh_parameterization/internal/shortest_path.h @@ -0,0 +1,167 @@ +// Copyright (c) 2016 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$ +// $Id$ +// +// +// Author(s) : + +#ifndef CGAL_SURFACE_MESH_PARAMETERIZATION_INTERNAL_SHORTEST_PATH_H +#define CGAL_SURFACE_MESH_PARAMETERIZATION_INTERNAL_SHORTEST_PATH_H + +#include + +#include +#include +#include + +#include +#include +#include + +namespace CGAL { + +namespace Surface_mesh_parameterization { + +namespace internal { + +class Dijkstra_end_exception : public std::exception +{ + const char* what() const throw () + { + return "Reached the target vertex"; + } +}; + +template +void output_shortest_paths_to_selection_file(const TriangleMesh& mesh, + const Seam_container& seams, + std::ofstream& os) +{ + typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; + typedef typename boost::graph_traits::edge_descriptor edge_descriptor; + + boost::unordered_map index_map; + + int counter = 0; + BOOST_FOREACH(vertex_descriptor vd, vertices(mesh)) { + index_map[vd] = counter++; + } + + os << std::endl /* vertices */ << std::endl /* faces */; + + BOOST_FOREACH(edge_descriptor ed, seams) { + // could be made more efficient... + os << index_map[source(ed, mesh)] << " " << index_map[target(ed, mesh)] << " "; + } + + os << std::endl; +} + +/// Visitor to stop Dijkstra when a target turns 'BLACK' (the point has been examined +/// through all its edges) +/// +/// \tparam TriangleMesh +template +class Stop_at_target_Dijkstra_visitor : boost::default_dijkstra_visitor +{ + typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; + typedef typename boost::graph_traits::edge_descriptor edge_descriptor; + + vertex_descriptor destination_vd; + +public: + Stop_at_target_Dijkstra_visitor(vertex_descriptor destination_vd) + : destination_vd(destination_vd) + { } + + void initialize_vertex(const vertex_descriptor& /*s*/, const TriangleMesh& /*mesh*/) const { } + void examine_vertex(const vertex_descriptor& /*s*/, const TriangleMesh& /*mesh*/) const { } + void examine_edge(const edge_descriptor& /*e*/, const TriangleMesh& /*mesh*/) const { } + void edge_relaxed(const edge_descriptor& /*e*/, const TriangleMesh& /*mesh*/) const { } + void discover_vertex(const vertex_descriptor& /*s*/, const TriangleMesh& /*mesh*/) const { } + void edge_not_relaxed(const edge_descriptor& /*e*/, const TriangleMesh& /*mesh*/) const { } + void finish_vertex(const vertex_descriptor &vd, const TriangleMesh& /* mesh*/) const + { + if(vd == destination_vd) + throw Dijkstra_end_exception(); + } +}; + +template +void compute_shortest_paths_between_two_cones(const TriangleMesh& mesh, + typename boost::graph_traits::vertex_descriptor source, + typename boost::graph_traits::vertex_descriptor target, + OutputIterator oi) +{ + if(source == target) { + std::cout << "Warning: the source and target are identical in 'shortest_path' " << std::endl; + return; + } + + typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; + typedef typename boost::graph_traits::edge_descriptor edge_descriptor; + + typedef Stop_at_target_Dijkstra_visitor Stop_visitor; + + typedef boost::unordered_map Pred_umap; + typedef boost::associative_property_map Pred_pmap; + + Pred_umap predecessor; + Pred_pmap pred_pmap(predecessor); + + Stop_visitor vis(target); + + try { + boost::dijkstra_shortest_paths(mesh, source, boost::predecessor_map(pred_pmap).visitor(vis)); + } catch (const std::exception& e) { + std::cout << e.what() << std::endl; + } + + // Draw the path from target to source and collect the edges along the way + vertex_descriptor s, t = target; + do { + s = get(pred_pmap, t); + std::pair e = edge(s, t, mesh); + CGAL_assertion(e.second); + if(e.second) { + *oi++ = e.first; + } + t = s; + } while (s != source); +} + +template +void compute_shortest_paths_between_cones(const TriangleMesh& mesh, + const Cones_vector& cones, + Seam_container& seams) +{ + CGAL_precondition(cones.size() == 3); + compute_shortest_paths_between_two_cones(mesh, cones[0], cones[2], + std::back_inserter(seams)); + compute_shortest_paths_between_two_cones(mesh, cones[1], cones[2], + std::back_inserter(seams)); + + std::ofstream out("shortest_path.selection.txt"); + output_shortest_paths_to_selection_file(mesh, seams, out); +} + +} // namespace internal + +} // namespace Surface_mesh_parameterization + +} // namespace CGAL + +#endif // CGAL_SURFACE_MESH_PARAMETERIZATION_INTERNAL_SHORTEST_PATH_H