mirror of https://github.com/CGAL/cgal
320 lines
11 KiB
C++
320 lines
11 KiB
C++
#include <CGAL/Cartesian.h>
|
|
#ifndef CGAL_USE_LAPACK
|
|
int main()
|
|
{
|
|
std::cerr << "Skip since LAPACK is not installed" << std::endl;
|
|
std::cerr << std::endl;
|
|
return 0;
|
|
}
|
|
#else
|
|
|
|
#include <CGAL/Monge_via_jet_fitting.h>
|
|
#include <CGAL/Lapack/Linear_algebra_lapack.h>
|
|
|
|
#include <fstream>
|
|
|
|
#include <boost/property_map.hpp>
|
|
|
|
#ifdef CGAL_USE_BOOST_PROGRAM_OPTIONS
|
|
#include <boost/program_options.hpp>
|
|
namespace po = boost::program_options;
|
|
#endif
|
|
|
|
using namespace std;
|
|
|
|
#include "PolyhedralSurf.h"
|
|
#include "PolyhedralSurf_operations.h"
|
|
#include "PolyhedralSurf_rings.h"
|
|
|
|
//Kernel of the PolyhedralSurf
|
|
typedef double DFT;
|
|
typedef CGAL::Cartesian<DFT> Data_Kernel;
|
|
typedef Data_Kernel::Point_3 DPoint;
|
|
typedef Data_Kernel::Vector_3 DVector;
|
|
|
|
//HDS
|
|
typedef PolyhedralSurf::Vertex_handle Vertex_handle;
|
|
typedef PolyhedralSurf::Vertex Vertex;
|
|
typedef PolyhedralSurf::Halfedge_handle Halfedge_handle;
|
|
typedef PolyhedralSurf::Halfedge Halfedge;
|
|
typedef PolyhedralSurf::Vertex_iterator Vertex_iterator;
|
|
typedef PolyhedralSurf::Facet_handle Facet_handle;
|
|
typedef PolyhedralSurf::Facet Facet;
|
|
|
|
struct Hedge_cmp{
|
|
bool operator()(Halfedge_handle a, Halfedge_handle b) const{
|
|
return &*a < &*b;
|
|
}
|
|
};
|
|
|
|
struct Facet_cmp{
|
|
bool operator()(Facet_handle a, Facet_handle b) const{
|
|
return &*a < &*b;
|
|
}
|
|
};
|
|
|
|
//Vertex property map, with std::map
|
|
typedef std::map<Vertex*, int> Vertex2int_map_type;
|
|
typedef boost::associative_property_map< Vertex2int_map_type > Vertex_PM_type;
|
|
typedef T_PolyhedralSurf_rings<PolyhedralSurf, Vertex_PM_type > Poly_rings;
|
|
|
|
//Hedge property map, with enriched Halfedge with its length
|
|
// typedef HEdge_PM<PolyhedralSurf> Hedge_PM_type;
|
|
// typedef T_PolyhedralSurf_hedge_ops<PolyhedralSurf, Hedge_PM_type> Poly_hedge_ops;
|
|
//Hedge property map, with std::map
|
|
typedef std::map<Halfedge_handle, double, Hedge_cmp> Hedge2double_map_type;
|
|
typedef boost::associative_property_map<Hedge2double_map_type> Hedge_PM_type;
|
|
typedef T_PolyhedralSurf_hedge_ops<PolyhedralSurf, Hedge_PM_type> Poly_hedge_ops;
|
|
|
|
// //Facet property map with enriched Facet with its normal
|
|
// typedef Facet_PM<PolyhedralSurf> Facet_PM_type;
|
|
// typedef T_PolyhedralSurf_facet_ops<PolyhedralSurf, Facet_PM_type> Poly_facet_ops;
|
|
//Facet property map, with std::map
|
|
typedef std::map<Facet_handle, Vector_3, Facet_cmp> Facet2normal_map_type;
|
|
typedef boost::associative_property_map<Facet2normal_map_type> Facet_PM_type;
|
|
typedef T_PolyhedralSurf_facet_ops<PolyhedralSurf, Facet_PM_type> Poly_facet_ops;
|
|
|
|
typedef double LFT;
|
|
typedef CGAL::Cartesian<LFT> Local_Kernel;
|
|
typedef CGAL::Monge_via_jet_fitting<Data_Kernel> My_Monge_via_jet_fitting;
|
|
typedef My_Monge_via_jet_fitting::Monge_form My_Monge_form;
|
|
|
|
|
|
// default parameter values and global variables
|
|
unsigned int d_fitting = 2;
|
|
unsigned int d_monge = 2;
|
|
unsigned int nb_rings = 0;//seek min # of rings to get the required #pts
|
|
unsigned int nb_points_to_use = 0;//
|
|
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
|
|
void gather_fitting_points(Vertex* v,
|
|
std::vector<DPoint> &in_points,
|
|
Vertex_PM_type& vpm)
|
|
{
|
|
//container to collect vertices of v on the PolyhedralSurf
|
|
std::vector<Vertex*> 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<Vertex*>::iterator
|
|
itb = gathered.begin(), ite = gathered.end();
|
|
CGAL_For_all(itb,ite) in_points.push_back((*itb)->point());
|
|
}
|
|
|
|
///////////////MAIN///////////////////////////////////////////////////////
|
|
#ifdef CGAL_USE_BOOST_PROGRAM_OPTIONS
|
|
int main(int argc, char *argv[])
|
|
#else
|
|
int main()
|
|
#endif
|
|
{
|
|
string if_name_string;
|
|
const char *if_name = NULL; //input file name
|
|
char *w_if_name = NULL; //as above, but / replaced by _
|
|
char* res4openGL_fname;
|
|
char* verbose_fname = NULL;
|
|
std::ofstream *out_4ogl = NULL, *out_verbose = NULL;
|
|
|
|
try {
|
|
#ifdef CGAL_USE_BOOST_PROGRAM_OPTIONS
|
|
po::options_description desc("Allowed options");
|
|
desc.add_options()
|
|
("help,h", "produce help message.")
|
|
("input-file,f", po::value<string>(&if_name_string)->default_value("data/ellipe0.003.off"),
|
|
"name of the input off file")
|
|
("degree-jet,d", po::value<unsigned int>(&d_fitting)->default_value(2),
|
|
"degree of the jet, 1 <= degre-jet <= 4")
|
|
("degree-monge,m", po::value<unsigned int>(&d_monge)->default_value(2),
|
|
"degree of the Monge rep, 1 <= degree-monge <= degree-jet")
|
|
("nb-rings,a", po::value<unsigned int>(&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<unsigned int>(&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")
|
|
("verbose,v", po::value<bool>(&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;
|
|
}
|
|
#else
|
|
std::cerr << "Command-line options require Boost.ProgramOptions" << std::endl;
|
|
if_name_string = "data/ellipe0.003.off";
|
|
d_fitting = 2;
|
|
d_monge = 2;
|
|
nb_rings = 0;
|
|
nb_points_to_use = 0;
|
|
verbose = false;
|
|
#endif
|
|
}
|
|
catch(exception& e) {
|
|
cerr << "error: " << e.what() << "\n";
|
|
return 1;
|
|
}
|
|
catch(...) {
|
|
cerr << "Exception of unknown type!\n";
|
|
}
|
|
|
|
//modify global variables which are fct of options:
|
|
min_nb_points = (d_fitting + 1) * (d_fitting + 2) / 2;
|
|
if (nb_points_to_use < min_nb_points && nb_points_to_use != 0)
|
|
{std::cerr << "the nb of points asked is not enough to perform the fitting" << std::endl; exit(0);}
|
|
|
|
//prepare output file names
|
|
//--------------------------
|
|
if_name = if_name_string.data();
|
|
assert(if_name != NULL);
|
|
w_if_name = new char[strlen(if_name)+1];
|
|
strcpy(w_if_name, if_name);
|
|
for(unsigned int i=0; i<strlen(w_if_name); i++)
|
|
if (w_if_name[i] == '/') w_if_name[i]='_';
|
|
cerr << if_name << '\n';
|
|
cerr << w_if_name << '\n';
|
|
|
|
res4openGL_fname = new char[strlen(w_if_name) + 10];// append .4ogl.txt
|
|
sprintf(res4openGL_fname, "%s.4ogl.txt", w_if_name);
|
|
out_4ogl = new std::ofstream(res4openGL_fname, std::ios::out);
|
|
assert(out_4ogl!=NULL);
|
|
//if verbose only...
|
|
if(verbose){
|
|
verbose_fname = new char[strlen(w_if_name) + 10];// append .verb.txt
|
|
sprintf(verbose_fname, "%s.verb.txt", w_if_name);
|
|
out_verbose = new std::ofstream( verbose_fname, std::ios::out);
|
|
assert(out_verbose != NULL);
|
|
CGAL::set_pretty_mode(*out_verbose);
|
|
}
|
|
unsigned int nb_vertices_considered = 0;//count vertices for verbose
|
|
|
|
//load the model from <mesh.off>
|
|
//------------------------------
|
|
PolyhedralSurf P;
|
|
std::ifstream stream(if_name);
|
|
stream >> P;
|
|
std::cout << "loadMesh... "<< "Polysurf with " << P.size_of_vertices()
|
|
<< " vertices and " << P.size_of_facets()
|
|
<< " facets. " << std::endl;
|
|
|
|
if(verbose)
|
|
(*out_verbose) << "Polysurf with " << P.size_of_vertices()
|
|
<< " vertices and " << P.size_of_facets()
|
|
<< " facets. " << std::endl;
|
|
//exit if not enough points in the model
|
|
if (min_nb_points > P.size_of_vertices()) exit(0);
|
|
|
|
//create property maps
|
|
//-----------------------------
|
|
//Vertex, using a std::map
|
|
Vertex2int_map_type vertex2props;
|
|
Vertex_PM_type vpm(vertex2props);
|
|
|
|
//Hedge, with enriched hedge
|
|
//HEdgePM_type hepm = get_hepm(boost::edge_weight_t(), P);
|
|
//Hedge, using a std::map
|
|
Hedge2double_map_type hedge2props;
|
|
Hedge_PM_type hepm(hedge2props);
|
|
|
|
//Facet PM, with enriched Facet
|
|
//FacetPM_type fpm = get_fpm(boost::vertex_attribute_t(), P);
|
|
//Facet PM, with std::map
|
|
Facet2normal_map_type facet2props;
|
|
Facet_PM_type fpm(facet2props);
|
|
|
|
//initialize Polyhedral data : length of edges, normal of facets
|
|
Poly_hedge_ops::compute_edges_length(P, hepm);
|
|
Poly_facet_ops::compute_facets_normals(P, fpm);
|
|
|
|
//MAIN LOOP: perform calculation for each vertex
|
|
//----------------------------------------------
|
|
std::vector<DPoint> in_points; //container for data points
|
|
Vertex_iterator vitb, vite;
|
|
|
|
//initialize the tag of all vertices to -1
|
|
vitb = P.vertices_begin(); vite = P.vertices_end();
|
|
CGAL_For_all(vitb,vite) put(vpm, &(*vitb), -1);
|
|
|
|
vitb = P.vertices_begin(); vite = P.vertices_end();
|
|
for (; vitb != vite; vitb++) {
|
|
//initialize
|
|
Vertex* v = &(*vitb);
|
|
in_points.clear();
|
|
My_Monge_form monge_form;
|
|
|
|
//gather points around the vertex using rings
|
|
gather_fitting_points(v, in_points, vpm);
|
|
|
|
//skip if the nb of points is to small
|
|
if ( in_points.size() < min_nb_points )
|
|
{std::cerr << "not enough pts for fitting this vertex" << in_points.size() << std::endl;
|
|
continue;}
|
|
|
|
// perform the fitting
|
|
My_Monge_via_jet_fitting monge_fit;
|
|
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 DVector normal_mesh = Poly_facet_ops::compute_vertex_average_unit_normal(v, fpm);
|
|
monge_form.comply_wrt_given_normal(normal_mesh);
|
|
|
|
//OpenGL output. Scaling for ppal dir, may be optimized with a
|
|
//global mean edges length computed only once on all edges of P
|
|
DFT scale_ppal_dir = Poly_hedge_ops::compute_mean_edges_length_around_vertex(v, hepm)/2;
|
|
|
|
(*out_4ogl) << v->point() << " ";
|
|
monge_form.dump_4ogl(*out_4ogl, scale_ppal_dir);
|
|
|
|
//verbose txt output
|
|
if (verbose) {
|
|
std::vector<DPoint>::iterator itbp = in_points.begin(), itep = in_points.end();
|
|
(*out_verbose) << "in_points list : " << std::endl ;
|
|
for (;itbp!=itep;itbp++) (*out_verbose) << *itbp << std::endl ;
|
|
|
|
(*out_verbose) << "--- vertex " << ++nb_vertices_considered
|
|
<< " : " << v->point() << std::endl
|
|
<< "number of points used : " << in_points.size() << std::endl
|
|
;// << monge_form;
|
|
}
|
|
} //all vertices processed
|
|
|
|
//cleanup filenames
|
|
//------------------
|
|
delete res4openGL_fname;
|
|
out_4ogl->close();
|
|
delete out_4ogl;
|
|
if(verbose) {
|
|
delete verbose_fname;
|
|
out_verbose->close();
|
|
delete out_verbose;
|
|
}
|
|
return 0;
|
|
}
|
|
#endif // CGAL_USE_LAPACK
|