From 9483280b5b07c0c5eb8a5fd0ff060cb0f4ff1b92 Mon Sep 17 00:00:00 2001 From: Andreas Fabri Date: Wed, 11 Mar 2015 22:58:55 +0100 Subject: [PATCH] Add example using Surface_mesh --- .../examples/Ridges_3/Ridges_Umbilics_SM.cpp | 405 ++++++++++++++++++ 1 file changed, 405 insertions(+) create mode 100644 Ridges_3/examples/Ridges_3/Ridges_Umbilics_SM.cpp diff --git a/Ridges_3/examples/Ridges_3/Ridges_Umbilics_SM.cpp b/Ridges_3/examples/Ridges_3/Ridges_Umbilics_SM.cpp new file mode 100644 index 00000000000..9de794c67ec --- /dev/null +++ b/Ridges_3/examples/Ridges_3/Ridges_Umbilics_SM.cpp @@ -0,0 +1,405 @@ + +//this is an enriched Polyhedron with facet normals +#include +#include +#include +#include "PolyhedralSurf_rings.h" +#include "compute_normals.h" +#include +#include +#include +#include +#include + +#ifdef CGAL_USE_BOOST_PROGRAM_OPTIONS +#include +namespace po = boost::program_options; +#endif + +using namespace std; + +typedef CGAL::Simple_cartesian Kernel; +typedef Kernel::FT FT; +typedef Kernel::Point_3 Point_3; +typedef Kernel::Vector_3 Vector_3; + +typedef CGAL::Surface_mesh PolyhedralSurf; + +typedef boost::graph_traits::vertex_descriptor vertex_descriptor; +typedef boost::graph_traits::vertex_iterator vertex_iterator; +typedef boost::graph_traits::face_descriptor face_descriptor; +/* + struct Face_cmp{ + bool operator()(face_descriptor a, face_descriptor b) const{ + return &*a < &*b; + } + }; +*/ +typedef T_PolyhedralSurf_rings Poly_rings; +typedef CGAL::Monge_via_jet_fitting Monge_via_jet_fitting; +typedef Monge_via_jet_fitting::Monge_form Monge_form; + +typedef PolyhedralSurf::Property_map Vertex2FT_property_map; +typedef PolyhedralSurf::Property_map Vertex2Vector_property_map; + + +typedef std::map Face2Vector_map; +typedef boost::associative_property_map< Face2Vector_map > Face2Vector_property_map; + +//RIDGES +typedef CGAL::Ridge_line Ridge_line; +typedef CGAL::Ridge_approximation < PolyhedralSurf, + Vertex2FT_property_map, + Vertex2Vector_property_map > Ridge_approximation; +//UMBILICS +typedef CGAL::Umbilic Umbilic; +typedef CGAL::Umbilic_approximation < PolyhedralSurf, + Vertex2FT_property_map, + Vertex2Vector_property_map > Umbilic_approximation; + +//create property maps + + +PolyhedralSurf::Property_map +vertex2k1_pm, vertex2k2_pm, + vertex2b0_pm, vertex2b3_pm, + vertex2P1_pm, vertex2P2_pm, vertex2P2_p; + +PolyhedralSurf::Property_map vertex2d1_pm, vertex2d2_pm; + +PolyhedralSurf::Property_map face2normal_pm; + +// default fct parameter values and global variables +unsigned int d_fitting = 3; +unsigned int d_monge = 3; +unsigned int nb_rings = 0;//seek min # of rings to get the required #pts +unsigned int nb_points_to_use = 0;// +CGAL::Ridge_order tag_order = CGAL::Ridge_order_3; +double umb_size = 2; +bool verbose = false; +unsigned int min_nb_points = (d_fitting + 1) * (d_fitting + 2) / 2; + +/* gather points around the vertex v using rings on the + polyhedralsurf. the collection of points resorts to 3 alternatives: + 1. the exact number of points to be used + 2. the exact number of rings to be used + 3. nothing is specified +*/ +template +void gather_fitting_points(vertex_descriptor v, + std::vector &in_points, + Poly_rings& poly_rings, + VertexPointMap vpm) +{ + //container to collect vertices of v on the PolyhedralSurf + std::vector gathered; + //initialize + in_points.clear(); + + //OPTION -p nb_points_to_use, with nb_points_to_use != 0. Collect + //enough rings and discard some points of the last collected ring to + //get the exact "nb_points_to_use" + if ( nb_points_to_use != 0 ) { + poly_rings.collect_enough_rings(v, nb_points_to_use, gathered);//, vpm); + if ( gathered.size() > nb_points_to_use ) gathered.resize(nb_points_to_use); + } + else { // nb_points_to_use=0, this is the default and the option -p is not considered; + // then option -a nb_rings is checked. If nb_rings=0, collect + // enough rings to get the min_nb_points required for the fitting + // else collect the nb_rings required + if ( nb_rings == 0 ) + poly_rings.collect_enough_rings(v, min_nb_points, gathered);//, vpm); + else poly_rings.collect_i_rings(v, nb_rings, gathered);//, vpm); + } + + //store the gathered points + std::vector::const_iterator + itb = gathered.begin(), ite = gathered.end(); + CGAL_For_all(itb,ite) in_points.push_back(get(vpm,*itb)); +} + +/* Use the jet_fitting package and the class Poly_rings to compute + diff quantities. +*/ +void compute_differential_quantities(PolyhedralSurf& P, Poly_rings& poly_rings) +{ + //container for approximation points + std::vector in_points; + + typedef boost::property_map::type VPM; + VPM vpm = get(CGAL::vertex_point,P); + + //MAIN LOOP + vertex_iterator vitb = P.vertices_begin(), vite = P.vertices_end(); + for (; vitb != vite; vitb++) { + //initialize + vertex_descriptor v = * vitb; + in_points.clear(); + Monge_form monge_form; + Monge_via_jet_fitting monge_fit; + + //gather points around the vertex using rings + gather_fitting_points(v, in_points, poly_rings, vpm); + + //exit if the nb of points is too small + if ( in_points.size() < min_nb_points ) + {std::cerr << "Too few points to perform the fitting" << std::endl; exit(1);} + + //For Ridges we need at least 3rd order info + assert( d_monge >= 3); + // run the main fct : perform the fitting + monge_form = monge_fit(in_points.begin(), in_points.end(), + d_fitting, d_monge); + + //switch min-max ppal curv/dir wrt the mesh orientation + const Vector_3 normal_mesh = computeFacetsAverageUnitNormal(P,v, face2normal_pm, Kernel()); + monge_form.comply_wrt_given_normal(normal_mesh); + + //Store monge data needed for ridge computations in property maps + vertex2d1_pm[v] = monge_form.maximal_principal_direction(); + vertex2d2_pm[v] = monge_form.minimal_principal_direction(); + vertex2k1_pm[v] = monge_form.coefficients()[0]; + vertex2k2_pm[v] = monge_form.coefficients()[1]; + vertex2b0_pm[v] = monge_form.coefficients()[2]; + vertex2b3_pm[v] = monge_form.coefficients()[5]; + if ( d_monge >= 4) { + //= 3*b1^2+(k1-k2)(c0-3k1^3) + vertex2P1_pm[v] = + 3*monge_form.coefficients()[3]*monge_form.coefficients()[3] + +(monge_form.coefficients()[0]-monge_form.coefficients()[1]) + *(monge_form.coefficients()[6] + -3*monge_form.coefficients()[0]*monge_form.coefficients()[0] + *monge_form.coefficients()[0]); + //= 3*b2^2+(k2-k1)(c4-3k2^3) + vertex2P2_pm[v] = + 3*monge_form.coefficients()[4]*monge_form.coefficients()[4] + +(-monge_form.coefficients()[0]+monge_form.coefficients()[1]) + *(monge_form.coefficients()[10] + -3*monge_form.coefficients()[1]*monge_form.coefficients()[1] + *monge_form.coefficients()[1]); + } + } //END FOR LOOP +} + + +///////////////MAIN/////////////////////////////////////////////////////// +#ifdef CGAL_USE_BOOST_PROGRAM_OPTIONS +int main(int argc, char *argv[]) +#else +int main() +#endif +{ + std::string if_name, of_name;// of_name same as if_name with '/' -> '_' + + try { +#ifdef CGAL_USE_BOOST_PROGRAM_OPTIONS + unsigned int int_tag; + po::options_description desc("Allowed options"); + desc.add_options() + ("help,h", "produce help message.") + ("input-file,f", po::value(&if_name)->default_value("data/poly2x^2+y^2-0.062500.off"), + "name of the input off file") + ("degree-jet,d", po::value(&d_fitting)->default_value(3), + "degree of the jet, 3 <= degre-jet <= 4") + ("degree-monge,m", po::value(&d_monge)->default_value(3), + "degree of the Monge rep, 3<= degree-monge <= degree-jet") + ("nb-rings,a", po::value(&nb_rings)->default_value(0), + "number of rings to collect neighbors. 0 means collect enough rings to make appro possible a>=1 fixes the nb of rings to be collected") + ("nb-points,p", po::value(&nb_points_to_use)->default_value(0), + "number of neighbors to use. 0 means this option is not considered, this is the default p>=1 fixes the nb of points to be used") + ("ridge_order,t", po::value(&int_tag)->default_value(3), + "Order of differential quantities used, must be 3 or 4") + ("umbilic-patch-size,u", po::value(&umb_size)->default_value(2), + "size of umbilic patches (as multiple of 1ring size)") + ("verbose,v", po::value(&verbose)->default_value(false), + "verbose output on text file") + ; + + po::variables_map vm; + po::store(po::parse_command_line(argc, argv, desc), vm); + po::notify(vm); + + if (vm.count("help")) { + cout << desc << "\n"; + return 1; + } + + + if (vm.count("ridge_order")){ + if ( int_tag == 3 ) tag_order = CGAL::Ridge_order_3; + if ( int_tag == 4 ) tag_order = CGAL::Ridge_order_4; + if ( int_tag != 3 && int_tag != 4 ) + {cerr << "ridge_order must be CGAL::Ridge_order_3 or CGAL::Ridge_order_4"; + return 1;} + } +#else + std::cerr << "Command-line options require Boost.ProgramOptions" << std::endl; + if_name = "data/poly2x^2+y^2-0.062500.off"; + d_fitting = 3; + d_monge = 3; + nb_rings = 0; + nb_points_to_use = 0; + umb_size = 2; + verbose = false; +#endif + } + + catch(exception& e) { + cerr << "error: " << e.what() << "\n"; + return 1; + } + catch(...) { + cerr << "Exception of unknown type!\n"; + } + + //modify global variables + min_nb_points = (d_fitting + 1) * (d_fitting + 2) / 2; + + //prepare output file names + assert(!if_name.empty()); + of_name = if_name; + for(unsigned int i=0; i + PolyhedralSurf P; + std::ifstream stream(if_name.c_str()); + stream >> P; + fprintf(stderr, "loadMesh %d Ves %d Facets\n", + (int)num_vertices(P), (int)num_faces(P)); + if(verbose) + out_verb << "Polysurf with " << num_vertices(P) + << " vertices and " << num_faces(P) + << " facets. " << std::endl; + + +vertex2k1_pm = P.add_property_map("v:k1",0).first; +vertex2k2_pm = P.add_property_map("v:k2",0).first; +vertex2b0_pm = P.add_property_map("v:b0",0).first; +vertex2b3_pm = P.add_property_map("v:b3",0).first; +vertex2P1_pm = P.add_property_map("v:P1",0).first; +vertex2P2_pm = P.add_property_map("v:P2",0).first; +vertex2P2_p = P.add_property_map("v:P2_p",0).first; + +vertex2d1_pm = P.add_property_map("v:d1",Vector_3(0,0,0)).first; +vertex2d2_pm = P.add_property_map("v:d2",Vector_3(0,0,0)).first; + +face2normal_pm = P.add_property_map("f:n",Vector_3(0,0,0)).first; + + //exit if not enough points in the model + if (min_nb_points > num_vertices(P)) + {std::cerr << "not enough points in the model" << std::endl; exit(1);} + + //initialize Polyhedral data : normal of facets + compute_facets_normals(P,face2normal_pm, Kernel()); + + //create a Poly_rings object + Poly_rings poly_rings(P); + + std::cout << "Compute differential quantities via jet fitting..." << std::endl; + //initialize the diff quantities property maps + compute_differential_quantities(P, poly_rings); + + //--------------------------------------------------------------------------- + //Ridges + //-------------------------------------------------------------------------- + std::cout << "Compute ridges..." << std::endl; + Ridge_approximation ridge_approximation(P, + vertex2k1_pm, vertex2k2_pm, + vertex2b0_pm, vertex2b3_pm, + vertex2d1_pm, vertex2d2_pm, + vertex2P1_pm, vertex2P2_pm ); + std::vector ridge_lines; + back_insert_iterator > ii(ridge_lines); + + //Find MAX_RIDGE, MIN_RIDGE, CREST_RIDGES + // ridge_approximation.compute_max_ridges(ii, tag_order); + // ridge_approximation.compute_min_ridges(ii, tag_order); + ridge_approximation.compute_crest_ridges(ii, tag_order); + + // or with the global function + CGAL::compute_max_ridges(P, + vertex2k1_pm, vertex2k2_pm, + vertex2b0_pm, vertex2b3_pm, + vertex2d1_pm, vertex2d2_pm, + vertex2P1_pm, vertex2P2_pm, + ii, tag_order); + + std::vector::iterator iter_lines = ridge_lines.begin(), + iter_end = ridge_lines.end(); + //OpenGL output + + typedef boost::property_map::type VPM; + VPM vpm = get(CGAL::vertex_point,P); + + for (;iter_lines!=iter_end;iter_lines++) (*iter_lines)->dump_4ogl(out_4ogl, vpm); + + + for (iter_lines = ridge_lines.begin();iter_lines!=iter_end;iter_lines++){ + //verbose txt output + if (verbose){ + (*iter_lines)->dump_verbose(out_verb,vpm); + } + delete *iter_lines; + } + + //--------------------------------------------------------------------------- + // UMBILICS + //-------------------------------------------------------------------------- + std::cout << "Compute umbilics..." << std::endl; + std::vector umbilics; + back_insert_iterator > umb_it(umbilics); + + //explicit construction of the class + // Umbilic_approximation umbilic_approximation(P, +// vertex2k1_pm, vertex2k2_pm, +// vertex2d1_pm, vertex2d2_pm); +// umbilic_approximation.compute(umb_it, umb_size); + //or global function call + CGAL::compute_umbilics(P, + vertex2k1_pm, vertex2k2_pm, + vertex2d1_pm, vertex2d2_pm, + umb_it, umb_size); + + std::vector::iterator iter_umb = umbilics.begin(), + iter_umb_end = umbilics.end(); + // output + std::cout << "nb of umbilics " << umbilics.size() << std::endl; + for (;iter_umb!=iter_umb_end;iter_umb++) std::cout << **iter_umb; + + //verbose txt output + if (verbose) { + out_verb << "nb of umbilics " << umbilics.size() << std::endl; + } + for ( iter_umb = umbilics.begin();iter_umb!=iter_umb_end;iter_umb++){ + if (verbose) { + out_verb << **iter_umb; + } + delete *iter_umb; + } + return 0; +}