From 07b9e89be58c75c6a46dd274f90dfdbef0672cc0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Thu, 9 Jul 2020 17:37:46 +0200 Subject: [PATCH] Refactor Iterative Authalic parameterization --- .../Fixed_border_iterative_parameterizer_3.h | 743 ----------- .../Iterative_authalic_parameterizer_3.h | 1158 ++++++++++++----- .../Iterative_parameterize.h | 142 -- 3 files changed, 859 insertions(+), 1184 deletions(-) delete mode 100755 Surface_mesh_parameterization/include/CGAL/Surface_mesh_parameterization/Fixed_border_iterative_parameterizer_3.h delete mode 100644 Surface_mesh_parameterization/include/CGAL/Surface_mesh_parameterization/Iterative_parameterize.h diff --git a/Surface_mesh_parameterization/include/CGAL/Surface_mesh_parameterization/Fixed_border_iterative_parameterizer_3.h b/Surface_mesh_parameterization/include/CGAL/Surface_mesh_parameterization/Fixed_border_iterative_parameterizer_3.h deleted file mode 100755 index 4b54e5637fe..00000000000 --- a/Surface_mesh_parameterization/include/CGAL/Surface_mesh_parameterization/Fixed_border_iterative_parameterizer_3.h +++ /dev/null @@ -1,743 +0,0 @@ -// This file has been adopted from CGAL, to integrate the below mentioned paper -// -// Paper : Learning to Reconstruct Symmetric Shapes using Planar Parameterization of 3D Surface -// Author(s) : Hardik Jain, Manuel Wöllhaf, Olaf Hellwich -// Conference : IEEE International Conference on Computer Vision Workshops (ICCVW) 2019 -// - -#ifndef CGAL_SURFACE_MESH_PARAMETERIZATION_FIXED_BORDER_ITERATIVE_PARAMETERIZER_3_H -#define CGAL_SURFACE_MESH_PARAMETERIZATION_FIXED_BORDER_ITERATIVE_PARAMETERIZER_3_H - -#define DEBUG_L0 1 - -#include - -#include -#include -#include -#include -#include - -#include -#include -#include -#include -#include -#include //for 2D functions -#include //for 3D functions - -#if defined(CGAL_EIGEN3_ENABLED) -#include -#endif - -#include - -#include -#include - -/// \file Fixed_border_iterative_parameterizer_3.h - -namespace CGAL { -namespace Surface_mesh_parameterization { - -// ------------------------------------------------------------------------------------ -// Declaration -// ------------------------------------------------------------------------------------ - -/// \ingroup PkgSurfaceMeshParameterizationMethods -/// -/// The class `Fixed_border_iterative_parameterizer_3` -/// is the base class of fixed border parameterization methods (Tutte, Floater, ...). -/// -/// A one-to-one mapping is guaranteed if the border of the surface is mapped onto a convex polygon. -/// -/// This class is a pure virtual class and thus cannot be instantiated. -/// Nevertheless, it implements most of the parameterization algorithm `parameterize()`. -/// Subclasses are *Strategies* \cgalCite{cgal:ghjv-dpero-95} that modify the behavior of this algorithm: -/// - They provide the template parameters `BorderParameterizer_` and `SolverTraits_`. -/// - They implement `compute_w_ij()` to compute w_ij = (i, j), coefficient of matrix A -/// for j neighbor vertex of i. -/// -/// \cgalModels `Parameterizer_3` -/// -/// \tparam TriangleMesh_ must be a model of `FaceGraph`. -/// -/// \tparam BorderParameterizer_ is a Strategy to parameterize the surface border -/// and must be a model of `Parameterizer_3`.
-/// %Default: -/// \code -/// Square_border_arc_length_parameterizer_3 -/// \endcode -/// -/// \tparam SolverTraits_ must be a model of `SparseLinearAlgebraTraits_d`.
-/// Note that the system is *not* symmetric because `Fixed_border_iterative_parameterizer_3` -/// does not remove border vertices from the system.
-/// %Default: If \ref thirdpartyEigen "Eigen" 3.1 (or greater) is available -/// and `CGAL_EIGEN3_ENABLED` is defined, then an overload of `Eigen_solver_traits` -/// is provided as default parameter: -/// \code -/// CGAL::Eigen_solver_traits< -/// Eigen::BiCGSTAB::EigenType, -/// Eigen::IncompleteLUT< double > > > -/// \endcode -/// -/// \sa `CGAL::Surface_mesh_parameterization::Iterative_authalic_parameterizer_3` -/// -template -class Fixed_border_iterative_parameterizer_3 -{ -public: -#ifndef DOXYGEN_RUNNING - typedef typename Default::Get >::type Border_parameterizer; - - typedef typename Default::Get::EigenType, - Eigen::IncompleteLUT > > -#else -#pragma message("Error: You must either provide 'SolverTraits_' or link CGAL with the Eigen library"); - SolverTraits_ // no parameter provided, and Eigen is not enabled: so don't compile! -#endif - >::type Solver_traits; -#else // DOXYGEN_RUNNING - typedef Border_parameterizer_ Border_parameterizer; - typedef SolverTraits_ Solver_traits; -#endif - - typedef TriangleMesh_ TriangleMesh; - - // Private types -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 CGAL::Vertex_around_target_circulator vertex_around_target_circulator; - typedef CGAL::Face_around_target_circulator face_around_target_circulator; - typedef CGAL::Vertex_around_face_circulator vertex_around_face_circulator; - - // Protected types -protected: - // Traits subtypes: - typedef typename internal::Kernel_traits::Kernel Kernel; - typedef typename internal::Kernel_traits::PPM PPM; - typedef typename Kernel::FT NT; - typedef typename Kernel::Point_2 Point_2; - typedef typename Kernel::Point_3 Point_3; - typedef typename Kernel::Vector_3 Vector_3; - - typedef std::unordered_set Vertex_set; - - // Solver traits subtypes: - typedef typename Solver_traits::Vector Vector; - typedef typename Solver_traits::Matrix Matrix; - - typedef CGAL::dynamic_face_property_t Face_double_tag; - typedef typename boost::property_map::type Face_Double_map; - typedef CGAL::dynamic_vertex_property_t Vertex_double_tag; - typedef typename boost::property_map::type Vertex_Double_map; - typedef CGAL::dynamic_vertex_property_t Vertex_int_tag; - typedef typename boost::property_map::type Vertex_Int_map; - typedef CGAL::dynamic_vertex_property_t Vertex_point2_tag; - typedef typename boost::property_map::type Vertex_point2_map; - - // Public operations -public: - /// Constructor - /// - /// \param border_param %Object that maps the surface's border to 2D space - /// \param sparse_la Traits object to access a sparse linear system - Fixed_border_iterative_parameterizer_3(Border_parameterizer border_param = Border_parameterizer(), - Solver_traits sparse_la = Solver_traits()) - : m_borderParameterizer(border_param), m_linearAlgebra(sparse_la), LScounter(0) - { } - - /// Destructor of base class should be virtual. - virtual ~Fixed_border_iterative_parameterizer_3() { } - - // Default copy constructor and operator =() are fine - - /// Compute a one-to-one mapping from a triangular 3D surface mesh - /// to a piece of the 2D space. - /// The mapping is piecewise linear (linear in each triangle). - /// The result is the (u,v) pair image of each vertex of the 3D surface. - /// - /// \tparam VertexUVmap must be a model of `ReadWritePropertyMap` with - /// `boost::graph_traits::%vertex_descriptor` as key type and - /// %Point_2 (type deduced from `TriangleMesh` using `Kernel_traits`) - /// as value type. - /// \tparam VertexIndexMap must be a model of `ReadablePropertyMap` with - /// `boost::graph_traits::%vertex_descriptor` as key type and - /// a unique integer as value type. - /// \tparam VertexParameterizedMap must be a model of `ReadWritePropertyMap` with - /// `boost::graph_traits::%vertex_descriptor` as key type and - /// a Boolean as value type. - /// - /// \param mesh a triangulated surface. - /// \param bhd a halfedge descriptor on the boundary of `mesh`. - /// \param uvmap an instanciation of the class `VertexUVmap`. - /// \param vimap an instanciation of the class `VertexIndexMap`. - /// \param vpmap an instanciation of the class `VertexParameterizedMap`. - /// \param iterations an integer number of iterations to run the parameterization. - /// \param error return error value of the iterative process. - /// - /// \pre `mesh` must be a triangular mesh. - /// \pre The mesh border must be mapped onto a convex polygon. - /// \pre The vertices must be indexed (`vimap` must be initialized) - /// - template - Error_code parameterize(TriangleMesh& mesh, - halfedge_descriptor bhd, - VertexUVmap uvmap1, - VertexIndexMap vimap1, - VertexParameterizedMap vpmap, - int& iterations, - double& error) - { - Error_code status = OK; - - Vertex_set vertices; - - internal::Containers_filler fc(mesh, vertices); - Polygon_mesh_processing::connected_component(face(opposite(bhd, mesh), mesh), mesh, - boost::make_function_output_iterator(fc)); - - int nbVertices = static_cast(vertices.size()); - - if(nbVertices == 0) - return ERROR_EMPTY_MESH; - - // Compute (u,v) for border vertices and mark them as "parameterized" - status = get_border_parameterizer().parameterize(mesh, bhd, uvmap1, vimap1, vpmap); - if (status != OK) - return status; - - // Create two sparse linear systems "A*Xu = Bu" and "A*Xv = Bv" (one line/column per vertex) - Matrix A(nbVertices, nbVertices); - Matrix A_prev(nbVertices, nbVertices); - Vector Xu(nbVertices), Xv(nbVertices), Bu(nbVertices), Bv(nbVertices); - double err[iterations]; - - // Initialize A, Xu, Xv, Bu and Bv after border parameterization - // Fill the border vertices' lines in both linear systems: - // "u = constant" and "v = constant" - // - initialize_system_from_mesh_border(A, Bu, Bv, mesh, bhd, uvmap1, vimap1); - - // Fill the matrix for the inner vertices v_i: compute A's coefficient - // w_ij for each neighbor j; then w_ii = - sum of w_ijs - std::unordered_set main_border; - for(vertex_descriptor v : vertices_around_face(bhd, mesh)) - main_border.insert(v); - - // create new uvmap & vimap - Vertex_point2_map uvmap = get(Vertex_point2_tag(), mesh); - Vertex_Int_map vimap = get(Vertex_int_tag(), mesh); - - for(vertex_descriptor v : vertices) - { - put(uvmap, v, get(uvmap1, v)); - put(vimap, v, get(vimap1, v)); - } - - // update last best uvmap wrt the border - lastBestuvmap = get(Vertex_point2_tag(), mesh); - for(vertex_descriptor v : vertices) - { - // border vertices only - if(main_border.find(v) != main_border.end()) - put(lastBestuvmap, v, get(uvmap, v)); - } - - int lastBesti = 0; - compute_faceArea(mesh); - compute_borderLength_3D(mesh); - double gamma = 1; - bool isChanged = false; - - // iterate it with the new weights - if(DEBUG_L0) - std::cout << std::endl; - - int i = 0; - while(i <= iterations) - { - if(DEBUG_L0) - std::cout << "Iteration " << i << std::flush; - - // update weights for inner vertices - for(vertex_descriptor v : vertices) - { - // inner vertices only - if(main_border.find(v) == main_border.end()) - { - // Compute the line i of matrix A for i inner vertex - if(i == 0) - { - //status = setup_inner_vertex_relations(A, A_prev, Bu, Bv, mesh, v, vimap); - status = setup_inner_vertex_relations_cotangent(A, Bu, Bv, mesh, v, vimap); - if(status != OK) - return status; - } - else - { - status = setup_iter_inner_vertex_relations(A, A_prev, Bu, Bv, mesh, v, vimap, uvmap, gamma); - if(status != OK) - return status; - } - } - } - - // solve linear equations - // Solve "A*Xu = Bu". On success, solution is (1/Du) * Xu. - // Solve "A*Xv = Bv". On success, solution is (1/Dv) * Xv. - NT Du = 0, Dv = 0; - if(!get_linear_algebra_traits().linear_solver(A, Bu, Xu, Du) || - !get_linear_algebra_traits().linear_solver(A, Bv, Xv, Dv)) - { - status = ERROR_CANNOT_SOLVE_LINEAR_SYSTEM; - } - else - { - LScounter = 0; - } - - if(status != OK) - { - if(LScounter < 4) - { - // modify the weights and re-try the linear solver - ++LScounter; - gamma /= 2; - continue; - } - else - { - status = OK; - break; - } - } - - // WARNING: this package does not support homogeneous coordinates! - CGAL_assertion(Du == 1.0); - CGAL_assertion(Dv == 1.0); - - // Copy A to A_prev, it is a computationally inefficient task but neccesary - copy_sparse_matrix(mesh, vertices, vimap, A, A_prev); - - // Copy Xu and Xv coordinates into the (u,v) pair of each vertex - for(vertex_descriptor v : vertices) - { - // inner vertices only - if(main_border.find(v) == main_border.end()) - { - int index = get(vimap,v); - put(uvmap, v, Point_2(Xu[index], Xv[index])); - put(vpmap, v, true); - } - } - - compute_faceWise_L2(mesh, uvmap); - compute_vertexWise_L2(mesh, vertices); - - err[i] = compute_area_distortion(mesh, vertices, main_border, uvmap); - - if(DEBUG_L0) - std::cout << " err " << err[i] << std::flush; - - if(err[i] <= err[lastBesti]) - { - update_uv(vertices, main_border, uvmap); - lastBesti = i; - isChanged = false; - if(DEBUG_L0) - std::cout << " *****" << std::flush; - } - else if(err[i] > 100) - { - break; - } - else - { - if(!isChanged) - { - gamma /= 2; - isChanged = true; - } - } - - ++i; - std::cout << std::endl; - } - - // Check postconditions - if(status != OK) - return status; - - if(i == 0 && i != iterations) - { - // means that the computation terminated for the first iteration may be because system was unsolvable - return ERROR_CANNOT_SOLVE_LINEAR_SYSTEM; - } - - for(vertex_descriptor v : vertices) - { - put(uvmap1, v, get(lastBestuvmap, v)); - put(vimap1, v, get(vimap, v)); - } - - iterations = lastBesti; - error = err[lastBesti]; - return status; - } - - void update_uv(const Vertex_set& vertices, - const std::unordered_set& main_border, - Vertex_point2_map uvmap) - { - for(vertex_descriptor v : vertices) - { - // inner vertices only - if(main_border.find(v) == main_border.end()) - put(lastBestuvmap, v, get(uvmap, v)); - } - } - - /// Initialize A, Bu and Bv after border parameterization. - /// Fill the border vertices' lines in both linear systems: - /// "u = constant" and "v = constant". - /// - /// \tparam VertexUVmap must be a model of `ReadWritePropertyMap` with - /// `boost::graph_traits::%vertex_descriptor` as key type and - /// %Point_2 (type deduced from `TriangleMesh` using `Kernel_traits`) - /// as value type. - /// \tparam VertexIndexMap must be a model of `ReadablePropertyMap` with - /// `boost::graph_traits::%vertex_descriptor` as key type and - /// a unique integer as value type. - /// - /// \param A the matrix in both linear system - /// \param Bu the right hand side vector in the linear system of x coordinates - /// \param Bv the right hand side vector in the linear system of y coordinates - /// \param mesh a triangulated surface. - /// \param bhd a halfedge descriptor on the boundary of `mesh`. - /// \param uvmap an instanciation of the class `VertexUVmap`. - /// \param vimap an instanciation of the class `VertexIndexMap`. - /// - /// \pre Vertices must be indexed (`vimap` must be initialized). - /// \pre A, Bu and Bv must be allocated. - /// \pre Border vertices must be parameterized. - template - void initialize_system_from_mesh_border(Matrix& A, Vector& Bu, Vector& Bv, - const TriangleMesh& mesh, - halfedge_descriptor bhd, - VertexUVmap uvmap, - VertexIndexMap vimap) const - { - for(halfedge_descriptor hd : halfedges_around_face(bhd, mesh)) - { - // Get vertex index in sparse linear system - int index = get(vimap, target(hd, mesh)); - - // Write a diagonal coefficient of A - A.set_coef(index, index, 1, true /*new*/); - - // get the halfedge uv - // Write constant in Bu and Bv - const Point_2& uv = get(uvmap, target(hd, mesh)); - Bu[index] = uv.x(); - Bv[index] = uv.y(); - } - } - - /// Compute w_ij, coefficient of matrix A for j neighbor vertex of i. - /// Implementation note: Subclasses must at least implement compute_w_ij(). - /// - /// \param mesh a triangulated surface. - /// \param main_vertex_v_i the vertex of `mesh` with index `i` - /// \param neighbor_vertex_v_j the vertex of `mesh` with index `j` - virtual NT compute_w_ij(const TriangleMesh& mesh, - vertex_descriptor main_vertex_v_i, - vertex_around_target_circulator neighbor_vertex_v_j) const = 0; - - virtual Error_code setup_inner_vertex_relations(Matrix& A, Matrix& A_prev, - Vector&, Vector&, - const TriangleMesh& mesh, - vertex_descriptor vertex, - Vertex_Int_map& vimap) - { - return OK; // TODO: Need to check - } - - virtual Error_code setup_inner_vertex_relations_cotangent(Matrix& A, - Vector&, Vector&, - const TriangleMesh& mesh, - vertex_descriptor vertex, - Vertex_Int_map& vimap) - { - return OK; // TODO: Need to check - } - - template - Error_code setup_iter_inner_vertex_relations(Matrix& A, Matrix& A_prev, - Vector&, Vector&, - TriangleMesh& mesh, - vertex_descriptor vertex, - VertexIndexMap vimap, - Vertex_point2_map& uvmap, - double gamma) - { - int i = get(vimap, vertex); - - // circulate over vertices around 'vertex' to compute w_ii and w_ijs - NT w_ii = 0; - int vertexIndex = 0; - - //const double sigma = getvL2(mesh, vertex, fL2Map, areaMap); - vertex_around_target_circulator v_j(halfedge(vertex, mesh), mesh), end = v_j; - CGAL_For_all(v_j, end) - { - // Get j index - int j = get(vimap, *v_j); - - // Call to virtual method to do the actual coefficient computation - //NT w_ij = A_prev.get_coef(i,j) / compute_sig_ij(vertex, *v_j) / gamma; - NT w_ij = A_prev.get_coef(i,j) / pow(compute_sig_ij(mesh, uvmap, vertex, *v_j, 1.0),gamma); - - // w_ii = - sum of w_ijs - w_ii -= w_ij; - - // Set w_ij in matrix - A.set_coef(i,j, w_ij, false); - - //A_prev.set_coef(i,j, w_ij, false); - ++vertexIndex; - } - - if(vertexIndex < 2) - return ERROR_NON_TRIANGULAR_MESH; - - // Set w_ii in matrix - A.set_coef(i,i, w_ii, true /*new*/); - return OK; - } - - virtual NT compute_faceArea(TriangleMesh& mesh) = 0; - virtual NT compute_faceWise_L2(TriangleMesh& mesh, Vertex_point2_map &uvmap) = 0; - virtual NT compute_vertexWise_L2(TriangleMesh& mesh, Vertex_set& vertices) = 0; - virtual double compute_sig_ij(TriangleMesh& mesh, Vertex_point2_map &uvmap, vertex_descriptor v_i, vertex_descriptor v_j, double gamma) = 0; - virtual NT compute_borderLength_3D(TriangleMesh& mesh) = 0; - virtual double distError(TriangleMesh& mesh) = 0; - - // Measure L2 stretch - template - double compute_dist_errors(const TriangleMesh& mesh, Vertex_set& vertices, - const std::unordered_set& main_border, - const VertexUVmap uvmap) - { - // iterate fpr all inner vertices and for each vertex - std::vector area_3D; - std::vector area_2D; - std::vector area_dist; - - double A_3D = Polygon_mesh_processing::area(mesh); - double A_2D = 1.; - - for(vertex_descriptor v : vertices) - { - // inner vertices only - if(main_border.find(v) == main_border.end()) - { - double a_2D = 0; - double a_3D = 0; - // find the area of all the adjacent faces to this vertex - face_around_target_circulator f_j(halfedge(v, mesh), mesh), end = f_j; - CGAL_For_all(f_j, end) - { - // get area in original mesh - a_3D += (Polygon_mesh_processing::face_area(*f_j, mesh)/A_3D); - - // get area in parameterised mesh - // iterate for all the vertices of this face and compute area - std::vector uv_points; - for(vertex_descriptor vd : vertices_around_face(halfedge(v, mesh), mesh)) - uv_points.push_back(get(uvmap, vd)); - - a_2D += (abs(CGAL::area(uv_points[0], uv_points[1],uv_points[2])) / A_2D); - } - - area_3D.push_back(a_3D); - area_2D.push_back(a_2D); - area_dist.push_back(square((a_3D/a_2D) - 1)); - } - } - - return sqrt(sum_vector(area_dist)); - } - - template - double compute_area_distortion(TriangleMesh& mesh, - const Vertex_set& vertices, - const std::unordered_set& main_border, - const VertexUVmap uvmap) - { - // iterate for all inner vertices and for each vertex - Face_Double_map area_3DMap = get(Face_double_tag(), mesh); - Face_Double_map area_2DMap = get(Face_double_tag(), mesh);; - std::vector area_dist; - std::vector innerFaces; - double A_3D = 0.0; - double A_2D = 0.0; - - for(face_descriptor fd : faces(mesh)) - { - put(area_3DMap, fd, Polygon_mesh_processing::face_area(fd, mesh)); - - // get area in parameterised mesh - std::vector uv_points; - for(vertex_descriptor vd : vertices_around_face(halfedge(fd, mesh), mesh)) - uv_points.push_back(get(uvmap,vd)); - - put(area_2DMap, fd, abs(CGAL::area(uv_points[0], uv_points[1], uv_points[2]))); - } - - for(face_descriptor fd : faces(mesh)) - { - A_3D += get(area_3DMap, fd); - A_2D += get(area_2DMap, fd); - } - - for(face_descriptor fd : faces(mesh)) - { - double a_3D = get(area_3DMap, fd); - double a_2D = get(area_2DMap, fd); - area_dist.push_back(abs(a_3D/A_3D - a_2D/A_2D)); - } - - return sum_vector(area_dist); - } - - template - T sum_vector(typename std::vector vec) - { - T sum = 0; - for(typename std::vector::iterator it=vec.begin(); it!= vec.end(); ++it) - sum += *it; - - return sum; - } - - template - void copy_sparse_matrix(TriangleMesh& mesh, - Vertex_set& vertices, - VertexIndexMap& vimap, - const Matrix& src, - Matrix& dest) - { - CGAL_assertion(src.row_dimension()==dest.row_dimension()); - CGAL_assertion(src.column_dimension()==dest.column_dimension()); - for(vertex_descriptor vertex : vertices) - { - int i = get(vimap, vertex); - vertex_around_target_circulator v_j(halfedge(vertex, mesh), mesh), end = v_j; - CGAL_For_all(v_j, end) - { - int j = get(vimap, *v_j); - dest.set_coef(i,j, src.get_coef(i,j), false); - } - } - } - - template - void print_matrix(const Vertex_set& vertices, - const VertexIndexMap vimap, - const Matrix &A, - const std::string name) - { - std::cout << "Matrix " << name << "(" << A.row_dimension() << "x" << A.column_dimension() << ")" << std::endl; - - Matrix A1(A.row_dimension(), A.column_dimension()); - int r=0, c=0; - for(vertex_descriptor v1 : vertices) - { - int i = get(vimap, v1); - for(vertex_descriptor v2 : vertices) - { - int j = get(vimap, v2); - A1.set_coef(r, c, A.get_coef(i, j)); - ++c; - } - - ++r; - c = 0; - } - - for(int r = 0; r < A.row_dimension(); ++r) - for(int c = 0; c < A.column_dimension(); ++c) - std::cout << std::setw(10) << A1.get_coef(r,c) << "\t" << std::flush; - std::cout << std::endl; - } - - template - void print_vector(const Vertex_set& vertices, - const VertexIndexMap vimap, - const Vector& A, - const std::string name) - { - std::cout << "Vector " << name << "(" << A.size() << ")" << std::endl; - Vector A1(A.size()); - int r = 0; - for(vertex_descriptor v1 : vertices) - { - int i = get(vimap,v1); - A1.set(r, A(i)); - ++r; - } - - for(int r=0; r #include +#include +#include #include +#include #include -#include +#include +#include #include +#include +#include +#include +#include #include //for 2D functions #include //for 3D functions @@ -25,10 +33,13 @@ #include +#include +#include #include #include #include -#include + +#define DEBUG_L0 1 // @fixme /// \file Iterative_authalic_parameterizer_3.h @@ -42,8 +53,9 @@ namespace Surface_mesh_parameterization { /// /// A one-to-one mapping is guaranteed if the surface's border is mapped onto a convex polygon. /// -/// This class is a strategy called by the main -/// parameterization algorithm `Fixed_border_iterative_parameterizer_3::parameterize()` and it: +/// @fixme wrong description +/// This class is a strategy called by the main parameterization algorithm +/// `Fixed_border_iterative_parameterizer_3::parameterize()` and it: /// - provides the template parameters `BorderParameterizer_` and `SolverTraits_`. /// - implements `compute_w_ij()` to compute w_ij = (i, j), coefficient of the matrix A /// for j neighbor vertex of i, based on Iterative Parameterization algorithm. @@ -77,22 +89,8 @@ template class Iterative_authalic_parameterizer_3 - : public Fixed_border_iterative_parameterizer_3< - TriangleMesh_, - typename Default::Get >::type, - typename Default::Get::EigenType, - Eigen::IncompleteLUT > > >::type > -#else -#pragma message("Error: You must either provide 'SolverTraits_' or link CGAL with the Eigen library") - SolverTraits_>::type > // no parameter provided, and Eigen is not enabled: don't compile -#endif { public: - #ifndef DOXYGEN_RUNNING typedef typename Default::Get >::type Border_parameterizer; @@ -101,7 +99,7 @@ public: #if defined(CGAL_EIGEN3_ENABLED) CGAL::Eigen_solver_traits< Eigen::BiCGSTAB::EigenType, - Eigen::IncompleteLUT > > + Eigen::IncompleteLUT > > #else SolverTraits_ // no parameter provided, and Eigen is not enabled: so don't compile! #endif @@ -113,13 +111,6 @@ public: typedef TriangleMesh_ TriangleMesh; - // Private types -private: - // Superclass - typedef Fixed_border_iterative_parameterizer_3 Base; - // Private types private: typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; @@ -130,145 +121,428 @@ private: typedef CGAL::Face_around_target_circulator face_around_target_circulator; typedef CGAL::Vertex_around_face_circulator vertex_around_face_circulator; - // Traits subtypes: - typedef typename Base::PPM PPM; - typedef typename Base::Kernel Kernel; - typedef typename Base::NT NT; - typedef typename Base::Point_2 Point_2; - typedef typename Base::Point_3 Point_3; - typedef typename Base::Vector_3 Vector_3; + typedef typename internal::Kernel_traits::Kernel Kernel; + + typedef typename Kernel::FT NT; + typedef typename Kernel::Point_2 Point_2; + typedef typename Kernel::Point_3 Point_3; + typedef typename Kernel::Vector_3 Vector_3; + + typedef typename internal::Kernel_traits::PPM PPM; + typedef typename boost::property_traits::reference PPM_ref; + typedef std::unordered_set Vertex_set; // Solver traits subtypes: typedef typename Solver_traits::Vector Vector; typedef typename Solver_traits::Matrix Matrix; - typedef CGAL::dynamic_face_property_t Face_double_tag; - typedef typename boost::property_map::type Face_Double_map; typedef CGAL::dynamic_vertex_property_t Vertex_double_tag; typedef typename boost::property_map::type Vertex_Double_map; + typedef CGAL::dynamic_vertex_property_t Vertex_bool_tag; + typedef typename boost::property_map::type Vertex_bool_map; typedef CGAL::dynamic_vertex_property_t Vertex_int_tag; - typedef typename boost::property_map::type Vertex_Int_map; + typedef typename boost::property_map::type Vertex_int_map; typedef CGAL::dynamic_vertex_property_t Vertex_point2_tag; typedef typename boost::property_map::type Vertex_point2_map; + typedef CGAL::dynamic_face_property_t Face_NT_tag; + typedef typename boost::property_map::type Face_NT_map; + + // Fields +private: + // Object that maps the surface's border onto a 2D space. + Border_parameterizer m_border_parameterizer; + + // Traits object to solve a sparse linear system + Solver_traits m_linear_algebra; + + // Object that keeps the last best UV Map + Vertex_point2_map m_last_best_uv_map; + + // Counter to keep track of failure of Linear solver + int m_linear_solver_failures; + + Vertex_Double_map m_vertex_L2_map; + Face_NT_map m_face_L2_map; + Face_NT_map m_face_areas; + + // Protected accessors +protected: + /// Get the object that maps the surface's border onto a 2D space. + Border_parameterizer& get_border_parameterizer() { return m_border_parameterizer; } + + /// Get the sparse linear algebra (traits object to access the linear system). + Solver_traits& get_linear_algebra_traits() { return m_linear_algebra; } // Public operations public: /// Constructor /// - /// \param border_param %Object that maps the surface's border to 2D space - /// \param sparse_la Traits object to access a sparse linear system. - Iterative_authalic_parameterizer_3(Border_parameterizer border_param = Border_parameterizer(), + /// \param border_parameterizer %Object that maps the surface's border to 2D space + /// \param iterations an integer number of iterations to run the parameterization + /// \param sparse_la Traits object to access a sparse linear system + /// + Iterative_authalic_parameterizer_3(Border_parameterizer border_parameterizer = Border_parameterizer(), Solver_traits sparse_la = Solver_traits()) - : Fixed_border_iterative_parameterizer_3(border_param, sparse_la) + : + m_border_parameterizer(border_parameterizer), + m_linear_algebra(sparse_la) { } - // Default copy constructor and operator =() are fine + // Disable copy constructor and operator =() due to property maps + Iterative_authalic_parameterizer_3(const Iterative_authalic_parameterizer_3&) = delete; + Iterative_authalic_parameterizer_3& operator=(const Iterative_authalic_parameterizer_3&) = delete; -protected: - virtual NT compute_faceArea(TriangleMesh& mesh) + // Distortion functions +public: + // Measure L2 stretch + template + NT compute_area_distortion(const FaceRange& face_range, + TriangleMesh& tmesh, + const VertexUVmap uvmap) { - areaMap = get(Face_double_tag(), mesh); - for(face_descriptor fd : faces(mesh)) - put(areaMap, fd, Polygon_mesh_processing::face_area(fd, mesh)); + Face_NT_map area_2DMap = get(Face_NT_tag(), tmesh); - return 1.0; - } + std::vector area_dist; + NT A_3D = 0; + NT A_2D = 0; - virtual NT compute_borderLength_3D(TriangleMesh& mesh) - { - return 1.0; - } - - virtual NT compute_faceWise_L2(TriangleMesh& mesh, Vertex_point2_map &uvmap) - { - fL2Map = get(Face_double_tag(), mesh); - for(face_descriptor fd : faces(mesh)) + for(face_descriptor f : face_range) { - Point_2 uv_points[3]; - Point_3 mesh_points[3]; - int i = 0; - for(vertex_descriptor vd : CGAL::vertices_around_face(halfedge(fd, mesh), mesh)) - { - uv_points[i] = get(uvmap, vd); - mesh_points[i] = mesh.point(vd); - ++i; - } - put(fL2Map, fd, getfL2(mesh_points, uv_points)); + const NT a_3D = get(m_face_areas, f); + A_3D += a_3D; + + // get area in parameterised mesh + const halfedge_descriptor h = halfedge(f, tmesh); + const NT a_2D = abs(CGAL::area(get(uvmap, source(h, tmesh)), + get(uvmap, target(h, tmesh)), + get(uvmap, target(next(h, tmesh), tmesh)))); + put(area_2DMap, f, a_2D); + A_2D += a_2D; } - return 1.0; + + for(face_descriptor f : face_range) + { + const NT a_3D = get(m_face_areas, f); + const NT a_2D = get(area_2DMap, f); + area_dist.push_back(abs(a_3D/A_3D - a_2D/A_2D)); + } + + return sqrt(std::accumulate(area_dist.begin(), area_dist.end(), NT(0))); } - virtual NT compute_vertexWise_L2(TriangleMesh& mesh, Vertex_set& vertices) + // IO Helpers +public: + template + void print_matrix(const Vertex_set& vertices, + const VertexIndexMap vimap, + const Matrix& A, + const std::string name) { - // update weights for vertices - vL2Map = get(Vertex_double_tag(), mesh); + std::cout << "Matrix " << name << "(" << A.row_dimension() << "x" << A.column_dimension() << ")" << std::endl; + + Matrix A1(A.row_dimension(), A.column_dimension()); + int r=0, c=0; + for(vertex_descriptor v1 : vertices) + { + int i = get(vimap, v1); + for(vertex_descriptor v2 : vertices) + { + int j = get(vimap, v2); + A1.set_coef(r, c, A.get_coef(i, j)); + ++c; + } + + ++r; + c = 0; + } + + for(int r=0; r + void print_vector(const Vertex_set& vertices, + const VertexIndexMap vimap, + const Vector& A, + const std::string name) + { + std::cout << "Vector " << name << "(" << A.size() << ")" << std::endl; + Vector A1(A.size()); + int r = 0; + for(vertex_descriptor v1 : vertices) + { + int i = get(vimap,v1); + A1.set(r, A(i)); + ++r; + } + + for(int r=0; r + void copy_sparse_matrix(const Matrix& src, + Matrix& dest, + const TriangleMesh& tmesh, + const Vertex_set& vertices, + const VertexIndexMap vimap) + { + CGAL_precondition(src.row_dimension() == dest.row_dimension()); + CGAL_precondition(src.column_dimension() == dest.column_dimension()); + + for(vertex_descriptor vertex : vertices) + { + const int i = get(vimap, vertex); + vertex_around_target_circulator v_j(halfedge(vertex, tmesh), tmesh), end = v_j; + CGAL_For_all(v_j, end) + { + const int j = get(vimap, *v_j); + dest.set_coef(i, j, src.get_coef(i, j), false); + } + } + } + +private: + double compute_vertex_L2(const TriangleMesh& tmesh, + const vertex_descriptor v) const + { + NT phi = 0, local_area = 0; + + for(face_descriptor f : CGAL::faces_around_target(halfedge(v, tmesh), tmesh)) + { + if(f == boost::graph_traits::null_face()) + continue; + + phi += CGAL::square(get(m_face_L2_map, f)) * get(m_face_areas, f); + local_area += get(m_face_areas, f); + } + + return sqrt(phi / local_area); + } + + void compute_vertices_L2(const Vertex_set& vertices, + TriangleMesh& tmesh) + { + m_vertex_L2_map = get(Vertex_double_tag(), tmesh); for(vertex_descriptor v : vertices) - put(vL2Map, v, getvL2(mesh, v)); - - return 1.0; + put(m_vertex_L2_map, v, compute_vertex_L2(tmesh, v)); } - virtual Error_code setup_inner_vertex_relations(Matrix& A, Matrix& A_prev, - Vector&, Vector&, - const TriangleMesh& mesh, - vertex_descriptor vertex, - Vertex_Int_map& vimap) + NT get_A(const std::array& uv_points) const { - // circulate over vertices around 'vertex' to compute w_ii and w_ijs - vertex_around_target_circulator v_j(halfedge(vertex, mesh), mesh), end_v_j = v_j; - std::vector NeighborList_; - int neighborsCounter = 0; - double theta_sum = 0.0; + NT A = (((uv_points[1]->x() - uv_points[0]->x()) * (uv_points[2]->y() - uv_points[0]->y())) + - ((uv_points[2]->x() - uv_points[0]->x()) * (uv_points[1]->y() - uv_points[0]->y()))) / NT(2); - // create NeighborList vector with vertex and vector + CGAL_warning(A != NT(0)); // means degenerate face in the param space + if(A == NT(0)) + return NT(1); + + return A; + } + + Point_3 get_Ss(const std::array& mesh_points, + const std::array& uv_points, + const NT den) const + { + const NT dt0 = uv_points[1]->y() - uv_points[2]->y(); + const NT dt1 = uv_points[2]->y() - uv_points[0]->y(); + const NT dt2 = uv_points[0]->y() - uv_points[1]->y(); + Point_3 Ss(den * (mesh_points[0]->x()*dt0 + mesh_points[1]->x()*dt1 + mesh_points[2]->x()*dt2), + den * (mesh_points[0]->y()*dt0 + mesh_points[1]->y()*dt1 + mesh_points[2]->y()*dt2), + den * (mesh_points[0]->z()*dt0 + mesh_points[1]->z()*dt1 + mesh_points[2]->z()*dt2)); + return Ss; + } + + Point_3 get_St(const std::array& mesh_points, + const std::array& uv_points, + const NT den) const + { + const NT ds0 = uv_points[2]->x() - uv_points[1]->x(); + const NT ds1 = uv_points[0]->x() - uv_points[2]->x(); + const NT ds2 = uv_points[1]->x() - uv_points[0]->x(); + Point_3 St(den * (mesh_points[0]->x()*ds0 + mesh_points[1]->x()*ds1 +mesh_points[2]->x()*ds2), + den * (mesh_points[0]->y()*ds0 + mesh_points[1]->y()*ds1 +mesh_points[2]->y()*ds2), + den * (mesh_points[0]->z()*ds0 + mesh_points[1]->z()*ds1 +mesh_points[2]->z()*ds2)); + return St; + } + + NT compute_inner_product(const Point_3& pointA, const Point_3& pointB) const + { + return ((pointA.x())*(pointB.x()) + (pointA.y())*(pointB.y()) + (pointA.z())*(pointB.z())); + } + + template + NT compute_face_L2(const face_descriptor f, + const TriangleMesh& tmesh, + const VertexUVMap uvmap, + const PPM ppmap) const + { + std::array uv_points; + std::array mesh_points; + + int i = 0; + for(vertex_descriptor v : CGAL::vertices_around_face(halfedge(f, tmesh), tmesh)) + { + // just to be safe in case of weird VPM returning temporaries + typename boost::property_traits::reference uvp = get(uvmap, v); + PPM_ref p = get(ppmap, v); + + uv_points[i] = &uvp; + mesh_points[i] = &p; + ++i; + } + + // Formula from Sander et al. 'Texture Mapping Progressive Meshes' + const NT A = get_A(uv_points); + const NT den = NT(1) / (NT(2) * A); + const Point_3 Ss = get_Ss(mesh_points, uv_points, den); + const Point_3 St = get_St(mesh_points, uv_points, den); + + const NT a = compute_inner_product(Ss, Ss); + const NT c = compute_inner_product(St, St); + + return sqrt((a+c) / NT(2)); + } + + template + void compute_faces_L2(const FaceRange& face_range, + TriangleMesh& tmesh, + const VertexUVMap uvmap) + { + m_face_L2_map = get(Face_NT_tag(), tmesh); + const PPM ppmap = get(vertex_point, tmesh); + + for(face_descriptor f : face_range) + put(m_face_L2_map, f, compute_face_L2(f, tmesh, uvmap, ppmap)); + } + + template + void initialize_faces_areas(const FaceRange& face_range, + TriangleMesh& tmesh) + { + m_face_areas = get(Face_NT_tag(), tmesh); + for(face_descriptor f : face_range) + put(m_face_areas, f, Polygon_mesh_processing::face_area(f, tmesh)); + } + + struct Neighbor_list + { + vertex_descriptor vertex; + double angle; + double length; + Vector_3 vector; + Point_2 uv; + double weight; + }; + + NT determinant(Point_2& v0, Point_2& v1) const + { + return (v0.x() * v1.y() - v1.x() * v0.y()); + } + + void get_bary_coords(const Point_2& uv0, const Point_2& uv1, const Point_2& uv2, + NT& tau0, NT& tau1, NT& tau2) const + { + const NT det0 = determinant(uv1, uv2); + const NT det1 = determinant(uv2, uv0); + const NT det2 = determinant(uv0, uv1); + const NT det3 = CGAL::determinant(Vector_2(uv1.x()-uv0.x(), uv1.y()-uv0.y()), + Vector_2(uv2.x()-uv0.x(), uv2.y()-uv0.y())); + CGAL_assertion(det3 > NT(0)); + if(det3 <= NT(0)) + det3 = NT(1); + + tau0 = det0 / det3; + tau1 = det1 / det3; + tau2 = det2 / det3; + } + + double angle(Vector_3& v0, Vector_3& v1) const // @fixme use helper's? + { + return std::acos(v0*v1 / (CGAL::sqrt(v0*v0) * CGAL::sqrt(v1*v1))); + } + + Error_code setup_inner_vertex_relations(Matrix& A, + Matrix& A_prev, + Vector&, + Vector&, + const TriangleMesh& tmesh, + vertex_descriptor v, + Vertex_int_map& vimap) + { + const PPM ppmap = get(vertex_point, tmesh); + + // circulate over vertices around 'vertex' to compute w_ii and w_ijs + std::vector neighbor_list; + int neighborsCounter = 0; + double theta_sum = 0.; + + // create Neighbor_list vector with vertex and vector + vertex_around_target_circulator v_j(halfedge(v, tmesh), tmesh), end_v_j = v_j; CGAL_For_all(v_j, end_v_j) { - NeighborList NL; + Neighbor_list NL; NL.vertex = *v_j; - NL.vector = Vector_3(mesh.point(vertex), mesh.point(*v_j)); + NL.vector = Vector_3(get(ppmap, v), tmesh.point(*v_j)); NL.length = sqrt(NL.vector.squared_length()); - NeighborList_.push_back(NL); + neighbor_list.push_back(NL); ++neighborsCounter; } if(neighborsCounter < 2) return ERROR_NON_TRIANGULAR_MESH; - if(neighborsCounter == 2 && mesh.is_border(vertex)) + if(neighborsCounter == 2 && is_border(v, tmesh)) { - std::cout << "Encountered inner border with valency-2 vertex (" << vertex << "), " + std::cout << "Encountered inner border with valency-2 vertex (" << v << "), " << "initializing with Tutte weights, this can affect optimization" << std::endl; // Tutte weights for(int k=0; k 0) w_ij *= -1.0; w_ii -= w_ij; // Get j index - const int j = get(vimap, NeighborList_[n].vertex); + const int j = get(vimap, neighbor_list[n].vertex); // Set w_ij in matrix - A.set_coef(i,j, w_ij, true /*new*/); - A_prev.set_coef(i,j, w_ij, true); + A.set_coef(i, j, w_ij, true /*new*/); + A_prev.set_coef(i, j, w_ij, true); } // Set w_ii in matrix @@ -330,148 +603,30 @@ protected: return OK; } - virtual Error_code setup_inner_vertex_relations_cotangent(Matrix& A, - Vector&, - Vector&, - const TriangleMesh& mesh, - vertex_descriptor vertex, - Vertex_Int_map& vimap) - { - int i = get(vimap,vertex); - - // circulate over vertices around 'vertex' to compute w_ii and w_ijs - NT w_ii = 0; - int vertexIndex = 0; - int neighborsCounter = 0; - vertex_around_target_circulator v_j(halfedge(vertex, mesh), mesh), end = v_j; - - CGAL_For_all(v_j, end) - { - neighborsCounter++; - } - - bool bcompute_w_ij = true; - if(neighborsCounter < 2) - { - return ERROR_NON_TRIANGULAR_MESH; - } - else if(neighborsCounter==2 && mesh.is_border(vertex)) - { - bcompute_w_ij = false; - std::cout << "Encountered inner border with valency-2 vertex (" << vertex << "), initializing with Tutte weights, this can affect optimization" << std::endl; - } - - CGAL_For_all(v_j, end) - { - // Call to virtual method to do the actual coefficient computation - NT w_ij = -1.0; - if(bcompute_w_ij) - w_ij *= compute_w_ij(mesh, vertex, v_j); - - // w_ii = - sum of w_ijs - w_ii -= w_ij; - - // Get j index - int j = get(vimap, *v_j); - - // Set w_ij in matrix - A.set_coef(i,j, w_ij, true /*new*/); - ++vertexIndex; - } - - // Set w_ii in matrix - A.set_coef(i,i, w_ii, true /*new*/); - return OK; - } - - virtual double compute_sig_ij(TriangleMesh& mesh, - Vertex_point2_map& uvmap, - vertex_descriptor v_i, - vertex_descriptor v_j, - double gamma) - { - double out = (pow(get(vL2Map,v_i),gamma)+pow(get(vL2Map,v_j),gamma)) / 2.0; - if(out <= 0.0) - std::cout << "compute_sig_ij <= 0.0" << std::endl; - return out; - } - - virtual double distError(TriangleMesh& mesh) - { - double varphi = 0.0; - double localArea = 0.0; - for(face_descriptor fd : faces(mesh)) - { - varphi += get(fL2Map,fd) * get(areaMap,fd); - localArea += get(areaMap, fd); - } - - double err = sqrt(varphi / localArea); - return err; - } - -private: - TriangleMesh mesh; - Vertex_Double_map vL2Map; - Face_Double_map areaMap; - Face_Double_map fL2Map; - - double getfL2(Point_3 mesh_points[], Point_2 uv_points[]) const - { - const double A = getA(uv_points); - Point_3 Ss = getSs(mesh_points, uv_points, A); - Point_3 St = getSt(mesh_points, uv_points, A); - - double a = innerProduct(Ss,Ss); - double c = innerProduct(St,St); - return sqrt((a+c)/2.0); - } - - double getvL2(const TriangleMesh& mesh, const vertex_descriptor vertex) const - { - halfedge_descriptor hf = halfedge(vertex, mesh); - face_around_target_circulator f_j(hf, mesh), end_f_j = f_j; - double varphi = 0.0; - double localArea = 0.0; - int i = 0; - CGAL_For_all(f_j, end_f_j) - { - if(*f_j > mesh.number_of_faces()) - continue; - - varphi += get(fL2Map, *f_j) * get(areaMap, *f_j); - localArea += get(areaMap, *f_j); - ++i; - } - - if(mesh.is_border(vertex) && mesh.degree(vertex) != i+1) - std::cerr << std::endl; - else if(!mesh.is_border(vertex) && mesh.degree(vertex) != i) - std::cerr << std::endl; - - return sqrt(varphi / localArea); - } - - /// Compute w_ij = (i, j), coefficient of matrix A for j neighbor vertex of i. - // These weights are shape preserving weights proposed in Floater1997 - NT compute_w_ij(const TriangleMesh& mesh, + /// Compute w_ij, coefficient of matrix A for j neighbor vertex of i. + /// Implementation note: Subclasses must at least implement compute_w_ij(). + /// + /// \param mesh a triangulated surface. + /// \param main_vertex_v_i the vertex of `mesh` with index `i` + /// \param neighbor_vertex_v_j the vertex of `mesh` with index `j` + NT compute_w_ij(const TriangleMesh& tmesh, vertex_descriptor main_vertex_v_i, vertex_around_target_circulator neighbor_vertex_v_j) const { - const PPM ppmap = get(vertex_point, mesh); + const PPM ppmap = get(vertex_point, tmesh); - const Point_3& position_v_i = get(ppmap, main_vertex_v_i); - const Point_3& position_v_j = get(ppmap, *neighbor_vertex_v_j); + const PPM_ref position_v_i = get(ppmap, main_vertex_v_i); + const PPM_ref position_v_j = get(ppmap, *neighbor_vertex_v_j); // Compute the square norm of v_j -> v_i vector Vector_3 edge = position_v_i - position_v_j; - double square_len = edge*edge; + NT square_len = edge*edge; // Compute cotangent of (v_k,v_j,v_i) corner (i.e. cotan of v_j corner) // if v_k is the vertex before v_j when circulating around v_i vertex_around_target_circulator previous_vertex_v_k = neighbor_vertex_v_j; --previous_vertex_v_k; - const Point_3& position_v_k = get(ppmap, *previous_vertex_v_k); + const PPM_ref position_v_k = get(ppmap, *previous_vertex_v_k); NT cotg_psi_ij = internal::cotangent(position_v_k, position_v_j, position_v_i); NT cotg_beta_ij = internal::cotangent(position_v_i, position_v_k, position_v_j); @@ -480,86 +635,491 @@ private: vertex_around_target_circulator next_vertex_v_l = neighbor_vertex_v_j; ++next_vertex_v_l; - const Point_3& position_v_l = get(ppmap, *next_vertex_v_l); + const Point_3 position_v_l = get(ppmap, *next_vertex_v_l); NT cotg_theta_ij = internal::cotangent(position_v_i, position_v_j, position_v_l); NT cotg_alpha_ij = internal::cotangent(position_v_j, position_v_l, position_v_i); - NT weight = 0.0; - CGAL_assertion(square_len != 0.0); // two points are identical! - if(square_len != 0.0) + NT weight = 0; + CGAL_assertion(square_len > NT(0)); // two points are identical! + if(square_len != NT(0)) weight = cotg_beta_ij + cotg_alpha_ij; return weight; } - struct NeighborList + Error_code setup_inner_vertex_relations_cotangent(Matrix& A, + Vector&, + Vector&, + const TriangleMesh& tmesh, + vertex_descriptor v, + Vertex_int_map& vimap) { - vertex_descriptor vertex; - double angle; - double length; - Vector_3 vector; - Point_2 uv; - double weight; - }; + const int i = get(vimap, v); - void baryCoords0(Point_2& uv0, Point_2& uv1, Point_2& uv2, double& tau0, double& tau1, double& tau2) - { - double det0 = determinant(uv1, uv2); - double det1 = determinant(uv2, uv0); - double det2 = determinant(uv0, uv1); - double det3 = CGAL::determinant(Vector_2(uv1.x()-uv0.x(), uv1.y()-uv0.y()), Vector_2(uv2.x()-uv0.x(), uv2.y()-uv0.y())); - if(det3 <= 0.0) - det3 = 1.0; + // circulate over vertices around 'vertex' to compute w_ii and w_ijs + NT w_ii = 0; - tau0 = det0 / det3; - tau1 = det1 / det3; - tau2 = det2 / det3; + bool use_uniform_weights = true; + const int neighborsCounter = degree(v, tmesh); + if(neighborsCounter < 2) + { + return ERROR_NON_TRIANGULAR_MESH; + } + else if(neighborsCounter == 2 && is_border(v, tmesh)) + { + use_uniform_weights = true; + std::cerr << "Encountered inner border with valency-2 vertex (" << get(vimap, v) << ") " + << "initializing with uniform weights, this can affect optimization" << std::endl; + } + + vertex_around_target_circulator vj(halfedge(v, tmesh), tmesh), end = vj; + CGAL_For_all(vj, end) + { + NT w_ij = NT(-1); + if(!use_uniform_weights) + w_ij *= compute_w_ij(tmesh, v, vj); + + // w_ii = - sum of w_ijs + w_ii -= w_ij; + + // Get j index + const int j = get(vimap, *vj); + + // Set w_ij in matrix + A.set_coef(i, j, w_ij, true /*new*/); + } + + // Set w_ii in matrix + A.set_coef(i, i, w_ii, true /*new*/); + + return OK; } - double determinant(Point_2& v0, Point_2& v1) + // Initialize the UV values with a first parameterization of the input. + template + Error_code setup_inner_vertex_relations_MVC(Matrix& A, + Vector&, + Vector&, + TriangleMesh& tmesh, + vertex_descriptor v, + VertexIndexMap& vimap) const { - return (v0.x()*v1.y() - v1.x()*v0.y()); + auto vpm = get_const_property_map(CGAL::vertex_point, tmesh); + CGAL::internal::Mean_value_weight compute_mvc(tmesh, vpm); + + const int i = get(vimap, v); + + // circulate over vertices around 'vertex' to compute w_ii and w_ijs + NT w_ii = 0; + int vertexIndex = 0; + + for(halfedge_descriptor h : CGAL::halfedges_around_target(v, tmesh)) + { + NT w_ij = NT(-1) * compute_mvc(h); + // w_ii = - sum of w_ijs + w_ii -= w_ij; + + // Get j index + int j = get(vimap, source(h, tmesh)); + + // Set w_ij in matrix + A.set_coef(i, j, w_ij, true /*new*/); + vertexIndex++; + } + + if (vertexIndex < 2) + return ERROR_NON_TRIANGULAR_MESH; + + // Set w_ii in matrix + A.set_coef(i,i, w_ii, true /*new*/); + return OK; } - double angle(Vector_3& v0, Vector_3& v1) + template + Error_code setup_iter_inner_vertex_relations(Matrix& A, + Matrix& A_prev, + Vector&, + Vector&, + const TriangleMesh& tmesh, + vertex_descriptor v, + VertexIndexMap vimap, + NT gamma) { - return std::acos(v0*v1 / (CGAL::sqrt(v0*v0) * CGAL::sqrt(v1*v1))); + const int i = get(vimap, v); + + // circulate over vertices around 'v' to compute w_ii and w_ij's + NT w_ii = 0; + int vi = 0; + + vertex_around_target_circulator v_j(halfedge(v, tmesh), tmesh), end = v_j; + CGAL_For_all(v_j, end) + { + // Get j index + const int j = get(vimap, *v_j); + + // NT w_ij = A_prev.get_coef(i, j) / compute_sig_ij(v, *v_j) / gamma; + NT w_ij = A_prev.get_coef(i, j) / pow(compute_sig_ij(v, *v_j, NT(1)), gamma); + + // w_ii = - sum of w_ij's + w_ii -= w_ij; + + // Set w_ij in matrix + A.set_coef(i, j, w_ij, false); + + ++vi; + } + + if(vi < 2) + return ERROR_NON_TRIANGULAR_MESH; + + // Set w_ii in matrix + A.set_coef(i,i, w_ii, true /*new*/); + return OK; } - double getA(Point_2 uv_points[]) const + NT compute_sig_ij(const vertex_descriptor v_i, + const vertex_descriptor v_j, + const NT gamma) { - double A = (((uv_points[1].x() - uv_points[0].x()) * (uv_points[2].y()-uv_points[0].y())) - - ((uv_points[2].x()-uv_points[0].x())*(uv_points[1].y()-uv_points[0].y()))) / 2.; - if(A == 0.0) - return 1.0; - return A; + const NT sig = (pow(get(m_vertex_L2_map, v_i), gamma) + pow(get(m_vertex_L2_map, v_j), gamma)) / 2.; + CGAL_assertion(sig > NT(0)); + return sig; } - Point_3 getSs(Point_3 mesh_points[3], Point_2 uv_points[3], const double &A) const + template + void store_new_best_uv_map(const Vertex_set& vertices, + VertexUVMap uvmap) { - double dt0 = uv_points[1].y() - uv_points[2].y(); - double dt1 = uv_points[2].y() - uv_points[0].y(); - double dt2 = uv_points[0].y() - uv_points[1].y(); - Point_3 Ss((mesh_points[0].x()*dt0 + mesh_points[1].x()*dt1 + mesh_points[2].x()*dt2 )/(2.0*A), - (mesh_points[0].y()*dt0 + mesh_points[1].y()*dt1 + mesh_points[2].y()*dt2 )/(2.0*A), - (mesh_points[0].z()*dt0 + mesh_points[1].z()*dt1 + mesh_points[2].z()*dt2 )/(2.0*A)); - return Ss; + for(vertex_descriptor v : vertices) + put(m_last_best_uv_map, v, get(uvmap, v)); } - Point_3 getSt(Point_3 mesh_points[3], Point_2 uv_points[3], const double &A) const +public: + /// Initialize A, Bu and Bv after border parameterization. + /// Fill the border vertices' lines in both linear systems: + /// "u = constant" and "v = constant". + /// + /// \tparam VertexUVmap must be a model of `ReadWritePropertyMap` with + /// `boost::graph_traits::%vertex_descriptor` as key type and + /// %Point_2 (type deduced from `TriangleMesh` using `Kernel_traits`) + /// as value type. + /// \tparam VertexIndexMap must be a model of `ReadablePropertyMap` with + /// `boost::graph_traits::%vertex_descriptor` as key type and + /// a unique integer as value type. + /// + /// \param A the matrix in both linear system + /// \param Bu the right hand side vector in the linear system of x coordinates + /// \param Bv the right hand side vector in the linear system of y coordinates + /// \param mesh a triangulated surface. + /// \param bhd a halfedge descriptor on the boundary of `mesh`. + /// \param uvmap an instanciation of the class `VertexUVmap`. + /// \param vimap an instanciation of the class `VertexIndexMap`. + /// + /// \pre Vertices must be indexed (`vimap` must be initialized). + /// \pre A, Bu and Bv must be allocated. + /// \pre Border vertices must be parameterized. + template + void initialize_system_from_mesh_border(Matrix& A, Vector& Bu, Vector& Bv, + const TriangleMesh& tmesh, + halfedge_descriptor bhd, + VertexUVmap uvmap, + VertexIndexMap vimap) const { - double ds0 = uv_points[2].x() - uv_points[1].x(); - double ds1 = uv_points[0].x() - uv_points[2].x(); - double ds2 = uv_points[1].x() - uv_points[0].x(); - Point_3 St((mesh_points[0].x()*ds0 + mesh_points[1].x()*ds1 +mesh_points[2].x()*ds2)/(2.0*A), - (mesh_points[0].y()*ds0 + mesh_points[1].y()*ds1 +mesh_points[2].y()*ds2)/(2.0*A), - (mesh_points[0].z()*ds0 + mesh_points[1].z()*ds1 +mesh_points[2].z()*ds2)/(2.0*A)); - return St; + for(halfedge_descriptor hd : halfedges_around_face(bhd, tmesh)) + { + // Get vertex index in sparse linear system + int index = get(vimap, target(hd, tmesh)); + + // Write a diagonal coefficient of A + A.set_coef(index, index, 1, true /*new*/); + + // get the halfedge uv + // Write constant in Bu and Bv + const Point_2& uv = get(uvmap, target(hd, tmesh)); + Bu[index] = uv.x(); + Bv[index] = uv.y(); + } } - double innerProduct(const Point_3& pointA, const Point_3& pointB) const + /// Compute a one-to-one mapping from a triangular 3D surface mesh + /// to a piece of the 2D space. + /// The mapping is piecewise linear (linear in each triangle). + /// The result is the (u,v) pair image of each vertex of the 3D surface. + /// + /// \tparam VertexUVmap must be a model of `ReadWritePropertyMap` with + /// `boost::graph_traits::%vertex_descriptor` as key type and + /// %Point_2 (type deduced from `TriangleMesh` using `Kernel_traits`) + /// as value type. + /// \tparam VertexIndexMap must be a model of `ReadablePropertyMap` with + /// `boost::graph_traits::%vertex_descriptor` as key type and + /// a unique integer as value type. + /// \tparam VertexParameterizedMap must be a model of `ReadWritePropertyMap` with + /// `boost::graph_traits::%vertex_descriptor` as key type and + /// a Boolean as value type. + /// + /// \param tmesh a triangulated surface + /// \param bhd a halfedge descriptor on the boundary of `mesh` + /// \param uvmap an instanciation of the class `VertexUVmap` + /// \param vimap an instanciation of the class `VertexIndexMap` + /// \param vpmap an instanciation of the class `VertexParameterizedMap` + /// \param iterations an integer number of iterations to run the parameterization + /// \param error return error value of the iterative process + /// + /// \pre `tmesh` must be a triangular mesh. + /// \pre The mesh border must be mapped onto a convex polygon. + /// \pre The vertices must be indexed (`vimap` must be initialized). + /// + template + Error_code parameterize(TriangleMesh& tmesh, + halfedge_descriptor bhd, + VertexUVmap uvmap, + VertexIndexMap vimap, + VertexParameterizedMap vpmap, + int& iterations, + double& error) { - return ((pointA.x())*(pointB.x()) + (pointA.y())*(pointB.y()) + (pointA.z())*(pointB.z())); + CGAL_precondition(is_valid_polygon_mesh(tmesh)); + CGAL_precondition(is_border(bhd, tmesh)); + + Error_code status = OK; + + Vertex_set cc_vertices; + std::vector cc_faces; + cc_faces.reserve(num_faces(tmesh)); + + internal::Containers_filler fc(tmesh, cc_vertices, &cc_faces); + Polygon_mesh_processing::connected_component(face(opposite(bhd, tmesh), tmesh), tmesh, + boost::make_function_output_iterator(fc)); + + int nv = static_cast(cc_vertices.size()); + if(nv == 0) + return ERROR_EMPTY_MESH; + + // Compute (u,v) for border vertices and mark them as "parameterized" + status = get_border_parameterizer().parameterize(tmesh, bhd, uvmap, vimap, vpmap); + if (status != OK) + return status; + + // Create two sparse linear systems "A*Xu = Bu" and "A*Xv = Bv" (one line/column per vertex) + Matrix A(nv, nv); + Matrix A_prev(nv, nv); + Vector Xu(nv), Xv(nv), Bu(nv), Bv(nv); + double err[iterations]; + + // Initialize A, Xu, Xv, Bu and Bv after border parameterization + // Fill the border vertices' lines in both linear systems: + // "u = constant" and "v = constant" + initialize_system_from_mesh_border(A, Bu, Bv, tmesh, bhd, uvmap, vimap); + + // Fill the matrix for the inner vertices v_i: compute A's coefficient + // w_ij for each neighbor j; then w_ii = - sum of w_ijs + std::unordered_set main_border; + for(vertex_descriptor v : vertices_around_face(bhd, tmesh)) + main_border.insert(v); // @todo use marks? + + m_last_best_uv_map = get(Vertex_point2_tag(), tmesh); + + initialize_faces_areas(cc_faces, tmesh); + + if(DEBUG_L0) // @fixme clean that stuff + std::cout << std::endl; + + int last_best_i = 0; + NT gamma = 1; // @todo what value should that be + bool is_changed = false; + + // iterate it with the new weights + int i = 0; + while(i < iterations) + { + if(DEBUG_L0) + std::cout << "Iteration " << i << ", gamma = " << gamma << std::flush; + + // update weights for inner vertices + for(vertex_descriptor v : cc_vertices) + { + // inner vertices only + if(main_border.count(v) == 0) + { + // Compute the line i of matrix A for i inner vertex + if(i == 0) + { +#if 0 + status = setup_inner_vertex_relations(A, A_prev, Bu, Bv, tmesh, v, vimap); +#elif 0 + status = setup_inner_vertex_relations_cotangent(A, Bu, Bv, tmesh, v, vimap); +#else + status = setup_inner_vertex_relations_MVC(A, Bu, Bv, tmesh, v, vimap); +#endif + + if(status != OK) + return status; + } + else + { + status = setup_iter_inner_vertex_relations(A, A_prev, Bu, Bv, tmesh, v, vimap, gamma); + if(status != OK) + return status; + } + } + } + + // solve linear equations + // Solve "A*Xu = Bu". On success, solution is (1/Du) * Xu. + // Solve "A*Xv = Bv". On success, solution is (1/Dv) * Xv. + NT Du = 0, Dv = 0; + if(!get_linear_algebra_traits().linear_solver(A, Bu, Xu, Du) || + !get_linear_algebra_traits().linear_solver(A, Bv, Xv, Dv)) + { + if(DEBUG_L0) + std::cout << " Linear solver failure #" << m_linear_solver_failures << std::endl; + + status = ERROR_CANNOT_SOLVE_LINEAR_SYSTEM; + } + else + { + m_linear_solver_failures = 0; + } + + if(status != OK) + { + if(m_linear_solver_failures < 4) + { + // modify the weights and re-try the linear solver + ++m_linear_solver_failures; + gamma /= 2; + continue; + } + else + { + status = OK; + break; + } + } + + // WARNING: this package does not support homogeneous coordinates! + CGAL_postcondition(Du == NT(1)); + CGAL_postcondition(Dv == NT(1)); + + // Copy A to A_prev, it is a computationally inefficient task but neccesary + copy_sparse_matrix(A, A_prev, tmesh, cc_vertices, vimap); + + // Copy Xu and Xv coordinates into the (u,v) pair of each vertex + for(vertex_descriptor v : cc_vertices) + { + // inner vertices only + if(main_border.count(v) == 0) + { + int index = get(vimap,v); + put(uvmap, v, Point_2(Xu[index], Xv[index])); + put(vpmap, v, true); + } + } + + if(DEBUG_L0) + { + std::ofstream out("last_solve.off"); + out.precision(17); + IO::output_uvmap_to_off(tmesh, bhd, uvmap, out); + out.close(); + } + + compute_faces_L2(cc_faces, tmesh, uvmap); + compute_vertices_L2(cc_vertices, tmesh); + + err[i] = compute_area_distortion(cc_faces, tmesh, uvmap); + + if(DEBUG_L0) + std::cout << " err " << err[i] << std::flush; + + if(err[i] <= err[last_best_i]) + { + store_new_best_uv_map(cc_vertices, uvmap); + + last_best_i = i; + is_changed = false; + + if(DEBUG_L0) + std::cout << " *****" << std::flush; + } + else if(err[i] > 100) // @fixme is that reasonnable + { + break; + } + else + { + if(!is_changed) + { + gamma /= 2; + is_changed = true; + } + } + + ++i; + std::cout << std::endl; + } + + // Check postconditions + if(status != OK) + return status; + + if(i == 0 && i != iterations) + { + // means that the computation terminated for the first iteration may be because system was unsolvable + return ERROR_CANNOT_SOLVE_LINEAR_SYSTEM; + } + + for(vertex_descriptor v : cc_vertices) + put(uvmap, v, get(m_last_best_uv_map, v)); + + iterations = last_best_i; + error = err[last_best_i]; + + return status; + } + + template + Error_code parameterize(TriangleMesh& tmesh, + halfedge_descriptor bhd, + VertexUVmap uvmap, + VertexIndexMap vimap, + VertexParameterizedMap vpmap, + int& iterations) // @fixme which default value + { + double unused_error; + return parameterize(tmesh, bhd, uvmap, vimap, vpmap, iterations, unused_error); + } + + template + Error_code parameterize(TriangleMesh& tmesh, + halfedge_descriptor bhd, + VertexUVmap uvmap, + VertexIndexMap vimap, + VertexParameterizedMap vpmap, + const int& iterations = 15) // @fixme which default value + { + int iter = iterations; // need a non-const ref + return parameterize(tmesh, bhd, uvmap, vimap, vpmap, iter); + } + + template + Error_code parameterize(TriangleMesh& tmesh, + halfedge_descriptor bhd, + VertexUVmap uvmap, + const int iterations = 15) // @fixme which default value + { + Vertex_int_map vimap = get(Vertex_int_tag(), tmesh); + internal::fill_index_map_of_cc(bhd, tmesh, vimap); + + Vertex_bool_map vpmap = get(Vertex_bool_tag(), tmesh); + + return parameterize(tmesh, bhd, uvmap, vimap, vpmap, iterations); } }; diff --git a/Surface_mesh_parameterization/include/CGAL/Surface_mesh_parameterization/Iterative_parameterize.h b/Surface_mesh_parameterization/include/CGAL/Surface_mesh_parameterization/Iterative_parameterize.h deleted file mode 100644 index c526ce2f13d..00000000000 --- a/Surface_mesh_parameterization/include/CGAL/Surface_mesh_parameterization/Iterative_parameterize.h +++ /dev/null @@ -1,142 +0,0 @@ -// This file has been adopted from CGAL, to integrate the below mentioned paper -// -// Paper : Learning to Reconstruct Symmetric Shapes using Planar Parameterization of 3D Surface -// Author(s) : Hardik Jain, Manuel Wöllhaf, Olaf Hellwich -// Conference : IEEE International Conference on Computer Vision Workshops (ICCVW) 2019 -// - -#ifndef CGAL_SURFACE_MESH_PARAMETERIZATION_ITERATIVE_PARAMETERIZE_H -#define CGAL_SURFACE_MESH_PARAMETERIZATION_ITERATIVE_PARAMETERIZE_H - -#include - -#include - -#include - -#include -#include - -#include - -#include -#include -#include -#include - -/// \file Iterative_parameterize.h -namespace CGAL { - -namespace Surface_mesh_parameterization { - -/// \ingroup PkgSurfaceMeshParameterizationMainFunction -/// -/// Compute a one-to-one mapping from a 3D triangle surface `mesh` to a -/// simple 2D domain. -/// The mapping is piecewise linear on the triangle mesh. -/// The result is a pair `(u,v)` of parameter coordinates for each vertex of the input mesh. -/// -/// A one-to-one mapping may be guaranteed or not, depending on -/// the chosen Parameterizer algorithm. -/// -/// \tparam TriangleMesh must be a model of `FaceGraph`. -/// \tparam Parameterizer must be a model of `Parameterizer_3`. -/// \tparam HD must be the halfedge_descriptor type corresponding to the graph -/// traits of TriangleMesh. -/// \tparam VertexUVmap must be a model of `ReadWritePropertyMap` with -/// `boost::graph_traits::%vertex_descriptor` as key type and -/// %Point_2 (type deduced from `TriangleMesh` using `Kernel_traits`) -/// as value type. -/// -/// \param mesh a triangulated surface. -/// \param parameterizer a parameterizer. -/// \param bhd a halfedge descriptor on the boundary of `mesh`. -/// \param uvm an instanciation of the class `VertexUVmap`. -/// \param iterations an integer number of iterations to run the parameterization. -/// \param error return error value of the iterative process. -/// -/// \pre `mesh` must be a triangular mesh. -/// \pre The mesh border must be mapped onto a convex polygon -/// (for fixed border parameterizations). -/// \pre The vertices must be indexed (vimap must be initialized). -/// -template -Error_code parameterize(TriangleMesh& mesh, - Parameterizer parameterizer, - HD bhd, - VertexUVmap uvm, - int& iterations, - double& error) -{ - typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; - - typedef boost::unordered_map Indices; - Indices indices; - CGAL::Polygon_mesh_processing::connected_component( - face(opposite(bhd, mesh), mesh), - mesh, - boost::make_function_output_iterator( - internal::Index_map_filler(mesh, indices))); - boost::associative_property_map vipm(indices); - - boost::unordered_set vs; - internal::Bool_property_map > vpm(vs); - - return parameterizer.parameterize(mesh, bhd, uvm, vipm, vpm, iterations, error); -} - -template -Error_code parameterize(TriangleMesh& mesh, - Parameterizer parameterizer, - HD bhd, - VertexUVmap uvm, - int& iterations) -{ - typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; - - typedef boost::unordered_map Indices; - Indices indices; - CGAL::Polygon_mesh_processing::connected_component( - face(opposite(bhd, mesh), mesh), - mesh, - boost::make_function_output_iterator( - internal::Index_map_filler(mesh, indices))); - boost::associative_property_map vipm(indices); - - boost::unordered_set vs; - internal::Bool_property_map > vpm(vs); - double error; - return parameterizer.parameterize(mesh, bhd, uvm, vipm, vpm, iterations, error); -} - -template -Error_code parameterize(TriangleMesh& mesh, - Parameterizer parameterizer, - HD bhd, - VertexUVmap uvm) -{ - typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; - - typedef boost::unordered_map Indices; - Indices indices; - CGAL::Polygon_mesh_processing::connected_component( - face(opposite(bhd, mesh), mesh), - mesh, - boost::make_function_output_iterator( - internal::Index_map_filler(mesh, indices))); - boost::associative_property_map vipm(indices); - - boost::unordered_set vs; - internal::Bool_property_map > vpm(vs); - int iterations = 0; - double error; - return parameterizer.parameterize(mesh, bhd, uvm, vipm, vpm, iterations, error); -} - -} // namespace Surface_mesh_parameterization - -} // namespace CGAL - -#include - -#endif // CGAL_ITERATIVE_PARAMETERIZE_H