mirror of https://github.com/CGAL/cgal
different versions for gsl and clapack
This commit is contained in:
Marc Pouget
parent
e1a108225c
commit
38565feade
|
|
@ -421,6 +421,8 @@ Jet_fitting_3/doc_tex/Jet_fitting_3/jet_fitting_basis.pdf -text svneol=unset#uns
|
|||
Jet_fitting_3/doc_tex/Jet_fitting_3_ref/template_dependence.eps -text
|
||||
Jet_fitting_3/doc_tex/Jet_fitting_3_ref/template_dependence.jpg -text
|
||||
Jet_fitting_3/doc_tex/Jet_fitting_3_ref/template_dependence.pdf -text
|
||||
Jet_fitting_3/examples/Jet_fitting_3/blind_1pt_gsl.C -text
|
||||
Jet_fitting_3/examples/Jet_fitting_3/blind_gsl.C -text
|
||||
Jet_fitting_3/examples/Jet_fitting_3/data/bezierpb3-0.007812.off -text
|
||||
Jet_fitting_3/examples/Jet_fitting_3/data/ellipe0.003.off -text
|
||||
Jet_fitting_3/examples/Jet_fitting_3/data/venus.off -text
|
||||
|
|
|
|||
|
|
@ -0,0 +1,134 @@
|
|||
#include <CGAL/basic.h>
|
||||
#include <CGAL/Cartesian.h>
|
||||
|
||||
#include <stdio.h>
|
||||
#include <iostream>
|
||||
#include <fstream>
|
||||
#include <stdlib.h>
|
||||
|
||||
#include <vector>
|
||||
#include "../../include/CGAL/Monge_via_jet_fitting_gsl.h"
|
||||
#include "GSL.h"
|
||||
|
||||
typedef double DFT;
|
||||
typedef CGAL::Cartesian<DFT> Data_Kernel;
|
||||
typedef Data_Kernel::Point_3 DPoint;
|
||||
typedef CGAL::Monge_rep<Data_Kernel> My_Monge_rep;
|
||||
|
||||
typedef double LFT;
|
||||
typedef CGAL::Cartesian<LFT> Local_Kernel;
|
||||
typedef CGAL::Monge_info<Local_Kernel> My_Monge_info;
|
||||
typedef CGAL::Monge_via_jet_fitting<Data_Kernel, Local_Kernel, GSL> My_Monge_via_jet_fitting;
|
||||
|
||||
|
||||
int main(int argc, char *argv[])
|
||||
{
|
||||
//check command line
|
||||
if (argc<5)
|
||||
{
|
||||
std::cout << "Usage : blind_1pt <inputPoints.txt> <output.txt> <d_fitting>, <d_monge>"
|
||||
<< std::endl;
|
||||
exit(-1);
|
||||
}
|
||||
//open the input file
|
||||
char name_in[20];
|
||||
sprintf(name_in, "%s", argv[1]);
|
||||
std::cout << name_in << '\n';
|
||||
|
||||
std::ifstream inFile( name_in, std::ios::in);
|
||||
if ( !inFile )
|
||||
{
|
||||
std::cerr << "cannot open file for input\n";
|
||||
exit(-1);
|
||||
}
|
||||
//initalize the in_points container
|
||||
double x, y, z;
|
||||
std::vector<DPoint> in_points;
|
||||
char ch[40];
|
||||
while (inFile >> ch) {
|
||||
x = atof(ch);
|
||||
inFile >> ch;
|
||||
y = atof(ch);
|
||||
inFile >> ch;
|
||||
z = atof(ch);
|
||||
DPoint p(x,y,z);
|
||||
in_points.push_back(p);
|
||||
}
|
||||
inFile.close();
|
||||
|
||||
// fct parameters
|
||||
int d_fitting = std::atoi(argv[3]);
|
||||
int d_monge = std::atoi(argv[4]);
|
||||
My_Monge_rep monge_rep;
|
||||
My_Monge_info monge_info;
|
||||
//run the main fct
|
||||
My_Monge_via_jet_fitting do_it(in_points.begin(), in_points.end(),
|
||||
d_fitting, d_monge,
|
||||
monge_rep, monge_info);
|
||||
|
||||
//open a file for output
|
||||
char name_out[20];
|
||||
sprintf(name_out, "%s", argv[2]);
|
||||
std::cout << name_out << '\n';
|
||||
|
||||
std::ofstream outFile( name_out, std::ios::out);
|
||||
if ( !outFile )
|
||||
{
|
||||
std::cerr << "cannot open file for output\n";
|
||||
exit(-1);
|
||||
}
|
||||
|
||||
//OUTPUT on outFile
|
||||
CGAL::set_pretty_mode(outFile);
|
||||
outFile << "vertex : " << in_points[0] << std::endl
|
||||
<< "number of points used : " << in_points.size() << std::endl
|
||||
<< "origin : " << monge_rep.origin_pt() << std::endl
|
||||
<< "d1 : " << monge_rep.d1() << std::endl
|
||||
<< "d2 : " << monge_rep.d2() << std::endl
|
||||
<< "n : " << monge_rep.n() << std::endl
|
||||
<< "cond_nb : " << monge_info.cond_nb() << std::endl
|
||||
<< "pca_eigen_vals " << monge_info.pca_eigen_vals()[0]
|
||||
<< " " << monge_info.pca_eigen_vals()[2]
|
||||
<< " " << monge_info.pca_eigen_vals()[3] << std::endl
|
||||
<< "pca_eigen_vecs : " << std::endl
|
||||
<< monge_info.pca_eigen_vecs()[0] << std::endl
|
||||
<< monge_info.pca_eigen_vecs()[1] << std::endl
|
||||
<< monge_info.pca_eigen_vecs()[2] << std::endl
|
||||
<< std::endl ;
|
||||
if ( d_monge >= 2)
|
||||
outFile << "k1 : " << monge_rep.coefficients()[0] << std::endl
|
||||
<< "k2 : " << monge_rep.coefficients()[1] << std::endl;
|
||||
|
||||
//OUTPUT on std::cout
|
||||
CGAL::set_pretty_mode(std::cout);
|
||||
std::cout << "vertex : " << in_points[0] << std::endl
|
||||
<< "number of points used : " << in_points.size() << std::endl
|
||||
<< "origin : " << monge_rep.origin_pt() << std::endl
|
||||
<< "d1 : " << monge_rep.d1() << std::endl
|
||||
<< "d2 : " << monge_rep.d2() << std::endl
|
||||
<< "n : " << monge_rep.n() << std::endl
|
||||
<< "cond_nb : " << monge_info.cond_nb() << std::endl
|
||||
<< "pca_eigen_vals " << monge_info.pca_eigen_vals()[0]
|
||||
<< " " << monge_info.pca_eigen_vals()[2]
|
||||
<< " " << monge_info.pca_eigen_vals()[3] << std::endl
|
||||
<< "pca_eigen_vecs : " << std::endl
|
||||
<< monge_info.pca_eigen_vecs()[0] << std::endl
|
||||
<< monge_info.pca_eigen_vecs()[1] << std::endl
|
||||
<< monge_info.pca_eigen_vecs()[2] << std::endl
|
||||
<< std::endl ;
|
||||
if ( d_monge >= 2)
|
||||
std::cout << "k1 : " << monge_rep.coefficients()[0] << std::endl
|
||||
<< "k2 : " << monge_rep.coefficients()[1] << std::endl;
|
||||
if ( d_monge >= 3)
|
||||
std::cout << "b0 : " << monge_rep.coefficients()[2] << std::endl
|
||||
<< "b1 : " << monge_rep.coefficients()[3] << std::endl
|
||||
<< "b2 : " << monge_rep.coefficients()[4] << std::endl
|
||||
<< "b3 : " << monge_rep.coefficients()[5] << std::endl;
|
||||
if ( d_monge >= 4)
|
||||
std::cout << "c0 : " << monge_rep.coefficients()[6] << std::endl
|
||||
<< "c1 : " << monge_rep.coefficients()[7] << std::endl
|
||||
<< "c2 : " << monge_rep.coefficients()[8] << std::endl
|
||||
<< "c3 : " << monge_rep.coefficients()[9] << std::endl
|
||||
<< "c4 : " << monge_rep.coefficients()[10] << std::endl;
|
||||
return 1;
|
||||
}
|
||||
|
|
@ -0,0 +1,301 @@
|
|||
#include <CGAL/basic.h>
|
||||
#include <CGAL/Cartesian.h>
|
||||
|
||||
#include <stdio.h>
|
||||
#include <iostream>
|
||||
#include <fstream>
|
||||
#include <stdlib.h>
|
||||
#include <vector>
|
||||
|
||||
#include <boost/property_map.hpp>
|
||||
|
||||
#include "../../include/CGAL/Monge_via_jet_fitting_gsl.h"
|
||||
#include "GSL.h"
|
||||
|
||||
#include "PolyhedralSurf.h"
|
||||
#include "PolyhedralSurf_operations.h"
|
||||
#include "PolyhedralSurf_rings.h"
|
||||
#include "options.h"//parsing command line
|
||||
|
||||
//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;
|
||||
|
||||
|
||||
|
||||
|
||||
//Kernel for local computations
|
||||
typedef double LFT;
|
||||
typedef CGAL::Cartesian<LFT> Local_Kernel;
|
||||
typedef CGAL::Monge_via_jet_fitting<Data_Kernel, Local_Kernel, GSL> My_Monge_via_jet_fitting;
|
||||
typedef CGAL::Monge_rep<Data_Kernel> My_Monge_rep;
|
||||
typedef CGAL::Monge_info<Local_Kernel> My_Monge_info;
|
||||
|
||||
//Syntax requirred by Options
|
||||
static const char *const optv[] = {
|
||||
"?|?",
|
||||
"f:fName <string>", //name of the input off file
|
||||
"d:deg <int>", //degree of the jet
|
||||
"m:mdegree <int>", //degree of the Monge rep
|
||||
"a:nrings <int>", //# rings
|
||||
"p:npoints <int>", //# points
|
||||
"v|",//verbose?
|
||||
NULL
|
||||
};
|
||||
|
||||
// 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());
|
||||
}
|
||||
|
||||
///////////////////////////////////////////////////////////////////////////////////////////////
|
||||
int main(int argc, char *argv[])
|
||||
{
|
||||
|
||||
char *if_name = NULL, //input file name
|
||||
*w_if_name = NULL; //as above, but / replaced by _
|
||||
char* res4openGL_fname;
|
||||
char* verbose_fname;
|
||||
std::ofstream *out_4ogl = NULL, *out_verbose = NULL;
|
||||
|
||||
//parse command line options
|
||||
//--------------------------
|
||||
int optchar;
|
||||
char *optarg;
|
||||
Options opts(*argv, optv);
|
||||
OptArgvIter iter(--argc, ++argv);
|
||||
while ((optchar = opts(iter, (const char *&) optarg))){
|
||||
switch (optchar){
|
||||
case 'f': if_name = optarg; break;
|
||||
case 'd': d_fitting = atoi(optarg); break;
|
||||
case 'm': d_monge = atoi(optarg); break;
|
||||
case 'a': nb_rings = atoi(optarg); break;
|
||||
case 'p': nb_points_to_use = atoi(optarg); break;
|
||||
case 'v': verbose=true; break;
|
||||
default:
|
||||
cerr << "Unknown command line option " << optarg;
|
||||
exit(0);
|
||||
}
|
||||
}
|
||||
//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
|
||||
//--------------------------
|
||||
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;
|
||||
fprintf(stderr, "loadMesh %d Ves %d Facets\n",
|
||||
P.size_of_vertices(), P.size_of_facets());
|
||||
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_rep monge_rep;
|
||||
My_Monge_info monge_info;
|
||||
|
||||
//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;}
|
||||
|
||||
// run the main fct : perform the fitting
|
||||
My_Monge_via_jet_fitting do_it(in_points.begin(), in_points.end(),
|
||||
d_fitting, d_monge,
|
||||
monge_rep, monge_info);
|
||||
|
||||
//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_rep.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_rep.dump_4ogl(*out_4ogl, scale_ppal_dir);
|
||||
|
||||
//verbose txt output
|
||||
if (verbose){
|
||||
|
||||
//debug
|
||||
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_rep.dump_verbose(*out_verbose);
|
||||
monge_info.dump_verbose(*out_verbose);
|
||||
}
|
||||
} //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 1;
|
||||
}
|
||||
|
||||
|
|
@ -6,21 +6,21 @@
|
|||
#---------------------------------------------------------------------#
|
||||
# Choose the right include file from the <cgalroot>/make directory.
|
||||
include $(CGAL_MAKEFILE)
|
||||
#include /home/mpouget/libs/CGAL-3.1/make/makefile_i686_Linux-2.6.12-10-386_g++-4.0.2
|
||||
|
||||
PROFOPT= -g
|
||||
CFLAGS=${PROFOPT}
|
||||
#---------------------------------------------------------------------#
|
||||
# compiler flags
|
||||
#---------------------------------------------------------------------#
|
||||
LAPACK_DIR=/usr/local/CLAPACK
|
||||
LAPACK_INCDIRS = -I$(LAPACK_DIR)/SRC -I$(LAPACK_DIR)
|
||||
GSL_INCDIR = -I$(GSL_DIR)/usr/local/gsl
|
||||
|
||||
CXXFLAGS = \
|
||||
-I../../include \
|
||||
$(CGAL_CXXFLAGS) \
|
||||
$(LONG_NAME_PROBLEM_CXXFLAGS) \
|
||||
$(LAPACK_INCDIRS) \
|
||||
$(GSL_INCDIR) \
|
||||
${CFLAGS}
|
||||
|
||||
#---------------------------------------------------------------------#
|
||||
|
|
@ -30,7 +30,7 @@ F2CDIR = $(LAPACK_DIR)/F2CLIBS
|
|||
LAPACK_LDLIBS = $(LAPACK_DIR)/lapack_LINUX.a \
|
||||
$(LAPACK_DIR)/blas_LINUX.a $(LAPACK_DIR)/tmglib_LINUX.a \
|
||||
$(F2CDIR)/libF77.a $(F2CDIR)/libI77.a -lm -lc
|
||||
#GSL_LIBS=-L${GSL_DIR}/lib -lgsl -lgslcblas -lm
|
||||
GSL_LIBS=-L${GSL_DIR}/lib -lgsl -lgslcblas -lm
|
||||
|
||||
LIBPATH = \
|
||||
$(CGAL_LIBPATH)
|
||||
|
|
@ -38,15 +38,15 @@ LIBPATH = \
|
|||
LDFLAGS = \
|
||||
$(LONG_NAME_PROBLEM_LDFLAGS) \
|
||||
$(CGAL_LDFLAGS) \
|
||||
${LAPACK_LDLIBS}
|
||||
|
||||
${LAPACK_LDLIBS} \
|
||||
${GSL_LIBS}
|
||||
|
||||
#---------------------------------------------------------------------#
|
||||
# target entries
|
||||
#---------------------------------------------------------------------#
|
||||
|
||||
all: \
|
||||
blind blind_1pt
|
||||
blind blind_1pt blind_gsl blind_1pt_gsl
|
||||
|
||||
BLIND_OBJS=options.o PolyhedralSurf.o blind.o
|
||||
blind: $(BLIND_OBJS)
|
||||
|
|
@ -57,11 +57,23 @@ blind_1pt: blind_1pt.o
|
|||
$(CGAL_CXX) $(LIBPATH) -o blind_1pt.exe blind_1pt.o \
|
||||
$(LDFLAGS)
|
||||
|
||||
BLIND_OBJS_GSL=options.o PolyhedralSurf.o blind_gsl.o
|
||||
blind_gsl: $(BLIND_OBJS_GSL)
|
||||
$(CGAL_CXX) $(LIBPATH) -g -o blind_gsl.exe $(BLIND_OBJS_GSL) \
|
||||
$(LDFLAGS)
|
||||
|
||||
blind_1pt_gsl: blind_1pt_gsl.o
|
||||
$(CGAL_CXX) $(LIBPATH) -o blind_1pt_gsl.exe blind_1pt_gsl.o \
|
||||
$(LDFLAGS)
|
||||
|
||||
|
||||
|
||||
rmexe:
|
||||
rm *.exe
|
||||
|
||||
clean: \
|
||||
blind.clean blind_1pt.clean options.clean PolyhedralSurf.clean \
|
||||
blind_gsl.clean blind_1pt_gsl.clean \
|
||||
rmexe
|
||||
|
||||
|
||||
|
|
@ -95,6 +107,80 @@ blind_1pt.o: ../../include/CGAL/jet_fitting_3_assertions.h
|
|||
blind_1pt.o: /usr/include/math.h /usr/include/bits/huge_val.h
|
||||
blind_1pt.o: /usr/include/bits/mathdef.h /usr/include/bits/mathcalls.h
|
||||
blind_1pt.o: LinAlg_lapack.h
|
||||
blind_1pt_gsl.o: /usr/include/stdio.h /usr/include/features.h
|
||||
blind_1pt_gsl.o: /usr/include/sys/cdefs.h /usr/include/gnu/stubs.h
|
||||
blind_1pt_gsl.o: /usr/include/bits/types.h /usr/include/bits/wordsize.h
|
||||
blind_1pt_gsl.o: /usr/include/bits/typesizes.h /usr/include/libio.h
|
||||
blind_1pt_gsl.o: /usr/include/_G_config.h /usr/include/wchar.h
|
||||
blind_1pt_gsl.o: /usr/include/bits/wchar.h /usr/include/gconv.h
|
||||
blind_1pt_gsl.o: /usr/include/bits/stdio_lim.h
|
||||
blind_1pt_gsl.o: /usr/include/bits/sys_errlist.h /usr/include/stdlib.h
|
||||
blind_1pt_gsl.o: /usr/include/sys/types.h /usr/include/time.h
|
||||
blind_1pt_gsl.o: /usr/include/endian.h /usr/include/bits/endian.h
|
||||
blind_1pt_gsl.o: /usr/include/sys/select.h /usr/include/bits/select.h
|
||||
blind_1pt_gsl.o: /usr/include/bits/sigset.h /usr/include/bits/time.h
|
||||
blind_1pt_gsl.o: /usr/include/sys/sysmacros.h
|
||||
blind_1pt_gsl.o: /usr/include/bits/pthreadtypes.h /usr/include/bits/sched.h
|
||||
blind_1pt_gsl.o: /usr/include/alloca.h
|
||||
blind_1pt_gsl.o: ../../include/CGAL/Monge_via_jet_fitting.h
|
||||
blind_1pt_gsl.o: ../../include/CGAL/jet_fitting_3_assertions.h
|
||||
blind_1pt_gsl.o: /usr/include/math.h /usr/include/bits/huge_val.h
|
||||
blind_1pt_gsl.o: /usr/include/bits/mathdef.h /usr/include/bits/mathcalls.h
|
||||
blind_1pt_gsl.o: GSL.h /usr/include/gsl/gsl_linalg.h
|
||||
blind_1pt_gsl.o: /usr/include/gsl/gsl_mode.h
|
||||
blind_1pt_gsl.o: /usr/include/gsl/gsl_permutation.h
|
||||
blind_1pt_gsl.o: /usr/include/gsl/gsl_types.h /usr/include/gsl/gsl_errno.h
|
||||
blind_1pt_gsl.o: /usr/include/errno.h /usr/include/bits/errno.h
|
||||
blind_1pt_gsl.o: /usr/include/linux/errno.h /usr/include/asm/errno.h
|
||||
blind_1pt_gsl.o: /usr/include/asm-i386/errno.h
|
||||
blind_1pt_gsl.o: /usr/include/gsl/gsl_check_range.h
|
||||
blind_1pt_gsl.o: /usr/include/gsl/gsl_vector.h
|
||||
blind_1pt_gsl.o: /usr/include/gsl/gsl_vector_complex_long_double.h
|
||||
blind_1pt_gsl.o: /usr/include/gsl/gsl_complex.h
|
||||
blind_1pt_gsl.o: /usr/include/gsl/gsl_vector_long_double.h
|
||||
blind_1pt_gsl.o: /usr/include/gsl/gsl_block_long_double.h
|
||||
blind_1pt_gsl.o: /usr/include/gsl/gsl_vector_complex.h
|
||||
blind_1pt_gsl.o: /usr/include/gsl/gsl_block_complex_long_double.h
|
||||
blind_1pt_gsl.o: /usr/include/gsl/gsl_vector_complex_double.h
|
||||
blind_1pt_gsl.o: /usr/include/gsl/gsl_vector_double.h
|
||||
blind_1pt_gsl.o: /usr/include/gsl/gsl_block_double.h
|
||||
blind_1pt_gsl.o: /usr/include/gsl/gsl_block_complex_double.h
|
||||
blind_1pt_gsl.o: /usr/include/gsl/gsl_vector_complex_float.h
|
||||
blind_1pt_gsl.o: /usr/include/gsl/gsl_vector_float.h
|
||||
blind_1pt_gsl.o: /usr/include/gsl/gsl_block_float.h
|
||||
blind_1pt_gsl.o: /usr/include/gsl/gsl_block_complex_float.h
|
||||
blind_1pt_gsl.o: /usr/include/gsl/gsl_vector_ulong.h
|
||||
blind_1pt_gsl.o: /usr/include/gsl/gsl_block_ulong.h
|
||||
blind_1pt_gsl.o: /usr/include/gsl/gsl_vector_long.h
|
||||
blind_1pt_gsl.o: /usr/include/gsl/gsl_block_long.h
|
||||
blind_1pt_gsl.o: /usr/include/gsl/gsl_vector_uint.h
|
||||
blind_1pt_gsl.o: /usr/include/gsl/gsl_block_uint.h
|
||||
blind_1pt_gsl.o: /usr/include/gsl/gsl_vector_int.h
|
||||
blind_1pt_gsl.o: /usr/include/gsl/gsl_block_int.h
|
||||
blind_1pt_gsl.o: /usr/include/gsl/gsl_vector_ushort.h
|
||||
blind_1pt_gsl.o: /usr/include/gsl/gsl_block_ushort.h
|
||||
blind_1pt_gsl.o: /usr/include/gsl/gsl_vector_short.h
|
||||
blind_1pt_gsl.o: /usr/include/gsl/gsl_block_short.h
|
||||
blind_1pt_gsl.o: /usr/include/gsl/gsl_vector_uchar.h
|
||||
blind_1pt_gsl.o: /usr/include/gsl/gsl_block_uchar.h
|
||||
blind_1pt_gsl.o: /usr/include/gsl/gsl_vector_char.h
|
||||
blind_1pt_gsl.o: /usr/include/gsl/gsl_block_char.h
|
||||
blind_1pt_gsl.o: /usr/include/gsl/gsl_matrix.h
|
||||
blind_1pt_gsl.o: /usr/include/gsl/gsl_matrix_complex_long_double.h
|
||||
blind_1pt_gsl.o: /usr/include/gsl/gsl_matrix_complex_double.h
|
||||
blind_1pt_gsl.o: /usr/include/gsl/gsl_matrix_complex_float.h
|
||||
blind_1pt_gsl.o: /usr/include/gsl/gsl_matrix_long_double.h
|
||||
blind_1pt_gsl.o: /usr/include/gsl/gsl_matrix_double.h
|
||||
blind_1pt_gsl.o: /usr/include/gsl/gsl_matrix_float.h
|
||||
blind_1pt_gsl.o: /usr/include/gsl/gsl_matrix_ulong.h
|
||||
blind_1pt_gsl.o: /usr/include/gsl/gsl_matrix_long.h
|
||||
blind_1pt_gsl.o: /usr/include/gsl/gsl_matrix_uint.h
|
||||
blind_1pt_gsl.o: /usr/include/gsl/gsl_matrix_int.h
|
||||
blind_1pt_gsl.o: /usr/include/gsl/gsl_matrix_ushort.h
|
||||
blind_1pt_gsl.o: /usr/include/gsl/gsl_matrix_short.h
|
||||
blind_1pt_gsl.o: /usr/include/gsl/gsl_matrix_uchar.h
|
||||
blind_1pt_gsl.o: /usr/include/gsl/gsl_matrix_char.h
|
||||
blind_1pt_gsl.o: /usr/include/gsl/gsl_eigen.h
|
||||
blind.o: /usr/include/stdio.h /usr/include/features.h
|
||||
blind.o: /usr/include/sys/cdefs.h /usr/include/gnu/stubs.h
|
||||
blind.o: /usr/include/bits/types.h /usr/include/bits/wordsize.h
|
||||
|
|
@ -160,6 +246,76 @@ blind.o: /usr/include/gsl/gsl_matrix_uchar.h
|
|||
blind.o: /usr/include/gsl/gsl_matrix_char.h /usr/include/gsl/gsl_eigen.h
|
||||
blind.o: PolyhedralSurf.h PolyhedralSurf_operations.h PolyhedralSurf_rings.h
|
||||
blind.o: options.h
|
||||
blind_gsl.o: /usr/include/stdio.h /usr/include/features.h
|
||||
blind_gsl.o: /usr/include/sys/cdefs.h /usr/include/gnu/stubs.h
|
||||
blind_gsl.o: /usr/include/bits/types.h /usr/include/bits/wordsize.h
|
||||
blind_gsl.o: /usr/include/bits/typesizes.h /usr/include/libio.h
|
||||
blind_gsl.o: /usr/include/_G_config.h /usr/include/wchar.h
|
||||
blind_gsl.o: /usr/include/bits/wchar.h /usr/include/gconv.h
|
||||
blind_gsl.o: /usr/include/bits/stdio_lim.h /usr/include/bits/sys_errlist.h
|
||||
blind_gsl.o: /usr/include/stdlib.h /usr/include/sys/types.h
|
||||
blind_gsl.o: /usr/include/time.h /usr/include/endian.h
|
||||
blind_gsl.o: /usr/include/bits/endian.h /usr/include/sys/select.h
|
||||
blind_gsl.o: /usr/include/bits/select.h /usr/include/bits/sigset.h
|
||||
blind_gsl.o: /usr/include/bits/time.h /usr/include/sys/sysmacros.h
|
||||
blind_gsl.o: /usr/include/bits/pthreadtypes.h /usr/include/bits/sched.h
|
||||
blind_gsl.o: /usr/include/alloca.h
|
||||
blind_gsl.o: ../../include/CGAL/Monge_via_jet_fitting_gsl.h
|
||||
blind_gsl.o: ../../include/CGAL/jet_fitting_3_assertions.h
|
||||
blind_gsl.o: /usr/include/math.h /usr/include/bits/huge_val.h
|
||||
blind_gsl.o: /usr/include/bits/mathdef.h /usr/include/bits/mathcalls.h GSL.h
|
||||
blind_gsl.o: /usr/include/gsl/gsl_linalg.h /usr/include/gsl/gsl_mode.h
|
||||
blind_gsl.o: /usr/include/gsl/gsl_permutation.h /usr/include/gsl/gsl_types.h
|
||||
blind_gsl.o: /usr/include/gsl/gsl_errno.h /usr/include/errno.h
|
||||
blind_gsl.o: /usr/include/bits/errno.h /usr/include/linux/errno.h
|
||||
blind_gsl.o: /usr/include/asm/errno.h /usr/include/asm-i386/errno.h
|
||||
blind_gsl.o: /usr/include/gsl/gsl_check_range.h /usr/include/gsl/gsl_vector.h
|
||||
blind_gsl.o: /usr/include/gsl/gsl_vector_complex_long_double.h
|
||||
blind_gsl.o: /usr/include/gsl/gsl_complex.h
|
||||
blind_gsl.o: /usr/include/gsl/gsl_vector_long_double.h
|
||||
blind_gsl.o: /usr/include/gsl/gsl_block_long_double.h
|
||||
blind_gsl.o: /usr/include/gsl/gsl_vector_complex.h
|
||||
blind_gsl.o: /usr/include/gsl/gsl_block_complex_long_double.h
|
||||
blind_gsl.o: /usr/include/gsl/gsl_vector_complex_double.h
|
||||
blind_gsl.o: /usr/include/gsl/gsl_vector_double.h
|
||||
blind_gsl.o: /usr/include/gsl/gsl_block_double.h
|
||||
blind_gsl.o: /usr/include/gsl/gsl_block_complex_double.h
|
||||
blind_gsl.o: /usr/include/gsl/gsl_vector_complex_float.h
|
||||
blind_gsl.o: /usr/include/gsl/gsl_vector_float.h
|
||||
blind_gsl.o: /usr/include/gsl/gsl_block_float.h
|
||||
blind_gsl.o: /usr/include/gsl/gsl_block_complex_float.h
|
||||
blind_gsl.o: /usr/include/gsl/gsl_vector_ulong.h
|
||||
blind_gsl.o: /usr/include/gsl/gsl_block_ulong.h
|
||||
blind_gsl.o: /usr/include/gsl/gsl_vector_long.h
|
||||
blind_gsl.o: /usr/include/gsl/gsl_block_long.h
|
||||
blind_gsl.o: /usr/include/gsl/gsl_vector_uint.h
|
||||
blind_gsl.o: /usr/include/gsl/gsl_block_uint.h
|
||||
blind_gsl.o: /usr/include/gsl/gsl_vector_int.h
|
||||
blind_gsl.o: /usr/include/gsl/gsl_block_int.h
|
||||
blind_gsl.o: /usr/include/gsl/gsl_vector_ushort.h
|
||||
blind_gsl.o: /usr/include/gsl/gsl_block_ushort.h
|
||||
blind_gsl.o: /usr/include/gsl/gsl_vector_short.h
|
||||
blind_gsl.o: /usr/include/gsl/gsl_block_short.h
|
||||
blind_gsl.o: /usr/include/gsl/gsl_vector_uchar.h
|
||||
blind_gsl.o: /usr/include/gsl/gsl_block_uchar.h
|
||||
blind_gsl.o: /usr/include/gsl/gsl_vector_char.h
|
||||
blind_gsl.o: /usr/include/gsl/gsl_block_char.h /usr/include/gsl/gsl_matrix.h
|
||||
blind_gsl.o: /usr/include/gsl/gsl_matrix_complex_long_double.h
|
||||
blind_gsl.o: /usr/include/gsl/gsl_matrix_complex_double.h
|
||||
blind_gsl.o: /usr/include/gsl/gsl_matrix_complex_float.h
|
||||
blind_gsl.o: /usr/include/gsl/gsl_matrix_long_double.h
|
||||
blind_gsl.o: /usr/include/gsl/gsl_matrix_double.h
|
||||
blind_gsl.o: /usr/include/gsl/gsl_matrix_float.h
|
||||
blind_gsl.o: /usr/include/gsl/gsl_matrix_ulong.h
|
||||
blind_gsl.o: /usr/include/gsl/gsl_matrix_long.h
|
||||
blind_gsl.o: /usr/include/gsl/gsl_matrix_uint.h
|
||||
blind_gsl.o: /usr/include/gsl/gsl_matrix_int.h
|
||||
blind_gsl.o: /usr/include/gsl/gsl_matrix_ushort.h
|
||||
blind_gsl.o: /usr/include/gsl/gsl_matrix_short.h
|
||||
blind_gsl.o: /usr/include/gsl/gsl_matrix_uchar.h
|
||||
blind_gsl.o: /usr/include/gsl/gsl_matrix_char.h /usr/include/gsl/gsl_eigen.h
|
||||
blind_gsl.o: PolyhedralSurf.h PolyhedralSurf_operations.h
|
||||
blind_gsl.o: PolyhedralSurf_rings.h options.h
|
||||
GSL.o: /usr/include/gsl/gsl_linalg.h /usr/include/gsl/gsl_mode.h
|
||||
GSL.o: /usr/include/gsl/gsl_permutation.h /usr/include/stdlib.h
|
||||
GSL.o: /usr/include/features.h /usr/include/sys/cdefs.h
|
||||
|
|
|
|||
|
|
@ -0,0 +1,796 @@
|
|||
#ifndef _MONGE_VIA_JET_FITTING_H_
|
||||
#define _MONGE_VIA_JET_FITTING_H_
|
||||
|
||||
#include <CGAL/basic.h>
|
||||
#include <CGAL/circulator.h>
|
||||
#include <CGAL/Linear_algebraCd.h>
|
||||
#include "jet_fitting_3_assertions.h"
|
||||
|
||||
//#include <CGAL/eigen.h> //for ALTERNATIVE with CGAL eigen code
|
||||
|
||||
#include <math.h>
|
||||
|
||||
CGAL_BEGIN_NAMESPACE
|
||||
|
||||
int fact(int n)
|
||||
{
|
||||
if (n == 0)
|
||||
return(1);
|
||||
return(n * fact(n-1));
|
||||
}
|
||||
////////////////////// CLASS Monge_rep ////////////////////////
|
||||
template <class DataKernel>
|
||||
class Monge_rep {
|
||||
public:
|
||||
typedef typename DataKernel::FT DFT;
|
||||
typedef typename DataKernel::Point_3 DPoint;
|
||||
typedef typename DataKernel::Vector_3 DVector;
|
||||
protected:
|
||||
//point on the fitted surface where diff quantities are computed
|
||||
DPoint m_origin_pt;
|
||||
//the monge trihedron (d1,d2,n) is orthonormal direct
|
||||
DVector m_d1; //maximal ppal dir
|
||||
DVector m_d2; //minimal ppal dir
|
||||
DVector m_n; //normal direction
|
||||
//coeff = (k1, k2, //ppal curv
|
||||
// b0, b1, b2, b3, //third order
|
||||
// c0, c1, c2, c3, c4) //fourth order
|
||||
// if (degree==1) no coeff needed
|
||||
std::vector<DFT> m_coefficients;
|
||||
|
||||
public:
|
||||
//constructor
|
||||
Monge_rep() {
|
||||
m_origin_pt = DPoint(0.,0.,0.);
|
||||
m_d1 = DVector(0.,0.,0.);
|
||||
m_d2 = DVector(0.,0.,0.);
|
||||
m_n = DVector(0.,0.,0.);
|
||||
m_coefficients = std::vector<DFT>();
|
||||
}
|
||||
~Monge_rep() {}
|
||||
//access
|
||||
const DPoint origin_pt() const { return m_origin_pt; }
|
||||
DPoint& origin_pt() { return m_origin_pt; }
|
||||
const DVector d1() const { return m_d1; }
|
||||
DVector& d1() { return m_d1; }
|
||||
const DVector d2() const { return m_d2; }
|
||||
DVector& d2() { return m_d2; }
|
||||
const DVector n() const { return m_n; }
|
||||
DVector& n() { return m_n; }
|
||||
const std::vector<DFT> coefficients() const { return m_coefficients; }
|
||||
std::vector<DFT>& coefficients() { return m_coefficients; }
|
||||
|
||||
//if d>=2, number of coeffs = (d+1)(d+2)/2 -4.
|
||||
//we remove cst, linear and the xy coeff which vanish
|
||||
void set_up(int degree);
|
||||
//switch min-max ppal curv/dir wrt a given normal orientation.
|
||||
// if given_normal.monge_normal < 0 then change the orientation
|
||||
// if z=g(x,y) in the basis (d1,d2,n) then in the basis (d2,d1,-n)
|
||||
// z=h(x,y)=-g(y,x)
|
||||
void comply_wrt_given_normal(const DVector given_normal);
|
||||
void dump_verbose(std::ostream& out_stream);
|
||||
void dump_4ogl(std::ostream& out_stream, const DFT scale);
|
||||
};
|
||||
|
||||
template <class DataKernel>
|
||||
void Monge_rep<DataKernel>::
|
||||
set_up(int degree) {
|
||||
if ( degree >= 2 ) std::fill_n(back_inserter(m_coefficients),
|
||||
(degree+1)*(degree+2)/2-4, 0.);
|
||||
}
|
||||
|
||||
template <class DataKernel>
|
||||
void Monge_rep<DataKernel>::
|
||||
comply_wrt_given_normal(const DVector given_normal)
|
||||
{
|
||||
if ( given_normal*this->n() < 0 )
|
||||
{
|
||||
n() = -n();
|
||||
std::swap(d1(), d2());
|
||||
if ( coefficients().size() >= 2)
|
||||
std::swap(coefficients()[0],coefficients()[1]);
|
||||
if ( coefficients().size() >= 6) {
|
||||
std::swap(coefficients()[2],coefficients()[5]);
|
||||
std::swap(coefficients()[3],coefficients()[4]);}
|
||||
if ( coefficients().size() >= 11) {
|
||||
std::swap(coefficients()[6],coefficients()[10]);
|
||||
std::swap(coefficients()[7],coefficients()[9]);}
|
||||
typename std::vector<DFT>::iterator itb = coefficients().begin(),
|
||||
ite = coefficients().end();
|
||||
for (;itb!=ite;itb++) { *itb = -(*itb); }
|
||||
}
|
||||
}
|
||||
|
||||
template <class DataKernel>
|
||||
void Monge_rep<DataKernel>::
|
||||
dump_verbose(std::ostream& out_stream)
|
||||
{
|
||||
out_stream << "origin : " << origin_pt() << std::endl
|
||||
<< "n : " << n() << std::endl;
|
||||
if ( coefficients().size() >= 2)
|
||||
out_stream << "d1 : " << d1() << std::endl
|
||||
<< "d2 : " << d2() << std::endl
|
||||
<< "k1 : " << coefficients()[0] << std::endl
|
||||
<< "k2 : " << coefficients()[1] << std::endl;
|
||||
if ( coefficients().size() >= 6)
|
||||
out_stream << "b0 : " << coefficients()[2] << std::endl
|
||||
<< "b1 : " << coefficients()[3] << std::endl
|
||||
<< "b2 : " << coefficients()[4] << std::endl
|
||||
<< "b3 : " << coefficients()[5] << std::endl;
|
||||
if ( coefficients().size() >= 11)
|
||||
out_stream << "c0 : " << coefficients()[6] << std::endl
|
||||
<< "c1 : " << coefficients()[7] << std::endl
|
||||
<< "c2 : " << coefficients()[8] << std::endl
|
||||
<< "c3 : " << coefficients()[9] << std::endl
|
||||
<< "c4 : " << coefficients()[10] << std::endl
|
||||
<< "P1 : " <<
|
||||
3*coefficients()[3]*coefficients()[3]
|
||||
+(coefficients()[0]-coefficients()[1])
|
||||
*(coefficients()[6]
|
||||
-3*coefficients()[0]*coefficients()[0]*coefficients()[0]) << std::endl
|
||||
//= 3*b2^2+(k2-k1)(c4-3k2^3)
|
||||
<< "P2 : " <<
|
||||
3*coefficients()[4]*coefficients()[4]
|
||||
-(coefficients()[0]-coefficients()[1])
|
||||
*(coefficients()[10]
|
||||
-3*coefficients()[1]*coefficients()[1]*coefficients()[1] )
|
||||
<< std::endl;
|
||||
|
||||
}
|
||||
|
||||
template <class DataKernel>
|
||||
void Monge_rep<DataKernel>::
|
||||
dump_4ogl(std::ostream& out_stream, const DFT scale)
|
||||
{
|
||||
CGAL_precondition( coefficients().size() >= 2 );
|
||||
out_stream << origin_pt() << " "
|
||||
<< d1() * scale << " "
|
||||
<< d2() * scale << " "
|
||||
<< coefficients()[0] << " "
|
||||
<< coefficients()[1] << " "
|
||||
<< std::endl;
|
||||
}
|
||||
|
||||
////////////////////// CLASS Monge_info ////////////////////////
|
||||
template <class LocalKernel>
|
||||
class Monge_info {
|
||||
public:
|
||||
typedef typename LocalKernel::FT LFT;
|
||||
typedef typename LocalKernel::Vector_3 LVector;
|
||||
protected:
|
||||
LFT m_pca_eigen_vals[3];
|
||||
LVector m_pca_eigen_vecs[3];
|
||||
LFT m_cond_nb;//of the least square system
|
||||
|
||||
public:
|
||||
//constructor
|
||||
Monge_info() {
|
||||
m_cond_nb = 0.;
|
||||
std::fill_n(m_pca_eigen_vals, 3, 0.);
|
||||
std::fill_n(m_pca_eigen_vecs, 3, LVector());
|
||||
}
|
||||
//access
|
||||
const LFT* pca_eigen_vals() const { return m_pca_eigen_vals; }
|
||||
LFT* pca_eigen_vals() { return m_pca_eigen_vals; }
|
||||
const LVector* pca_eigen_vecs() const { return m_pca_eigen_vecs; }
|
||||
LVector* pca_eigen_vecs() { return m_pca_eigen_vecs; }
|
||||
const LFT cond_nb() const { return m_cond_nb; }
|
||||
LFT& cond_nb() { return m_cond_nb; }
|
||||
|
||||
void dump_verbose(std::ostream& out_stream);
|
||||
};
|
||||
|
||||
|
||||
template <class LocalKernel>
|
||||
void Monge_info<LocalKernel>::
|
||||
dump_verbose(std::ostream& out_stream)
|
||||
{
|
||||
out_stream << "cond_nb : " << cond_nb() << std::endl
|
||||
<< "pca_eigen_vals " << pca_eigen_vals()[0]
|
||||
<< " " << pca_eigen_vals()[2]
|
||||
<< " " << pca_eigen_vals()[3] << std::endl
|
||||
<< "pca_eigen_vecs : " << std::endl
|
||||
<< pca_eigen_vecs()[0] << std::endl
|
||||
<< pca_eigen_vecs()[1] << std::endl
|
||||
<< pca_eigen_vecs()[2] << std::endl;
|
||||
}
|
||||
|
||||
|
||||
|
||||
////////////////////// CLASS Monge_via_jet_fitting ////////////////////////
|
||||
template < class DataKernel, class LocalKernel, class LinAlgTraits>
|
||||
class Monge_via_jet_fitting {
|
||||
public:
|
||||
typedef DataKernel Data_Kernel;
|
||||
typedef LocalKernel Local_Kernel;
|
||||
typedef typename std::vector<typename Data_Kernel::Point_3>::iterator Range_Iterator;
|
||||
typedef Monge_rep<Data_Kernel> Monge_rep;
|
||||
typedef Monge_info<Local_Kernel> Monge_info;
|
||||
|
||||
public:
|
||||
Monge_via_jet_fitting(Range_Iterator begin, Range_Iterator end,
|
||||
int d, int dprime,
|
||||
Monge_rep &monge_rep, Monge_info &monge_info);
|
||||
|
||||
protected:
|
||||
typedef typename Local_Kernel::FT LFT;
|
||||
typedef typename Local_Kernel::Point_3 LPoint;
|
||||
typedef typename Local_Kernel::Vector_3 LVector;
|
||||
typedef CGAL::Aff_transformation_3<Local_Kernel> Aff_transformation;
|
||||
|
||||
typedef typename Data_Kernel::FT DFT;
|
||||
typedef typename Data_Kernel::Point_3 DPoint;
|
||||
|
||||
typedef typename LinAlgTraits::Vector LAVector;
|
||||
typedef typename LinAlgTraits::Matrix LAMatrix;
|
||||
|
||||
protected:
|
||||
int deg;
|
||||
int deg_monge;
|
||||
int nb_d_jet_coeff;
|
||||
int nb_input_pts;
|
||||
LFT preconditionning;
|
||||
CGAL::Sqrt<LFT> Lsqrt;
|
||||
|
||||
//translate_p0 changes the origin of the world to p0 the first point
|
||||
// of the input data points
|
||||
//change_world2fitting (coord of a vector in world) = coord of this
|
||||
// vector in fitting. The matrix tranform has as lines the coord of
|
||||
// the basis vectors of fitting in the world coord.
|
||||
//idem for change_fitting2monge
|
||||
Aff_transformation translate_p0, change_world2fitting,
|
||||
change_fitting2monge;
|
||||
|
||||
//eigen val and vect stored in monge_info,
|
||||
// change_world2fitting is computed
|
||||
void compute_PCA(Range_Iterator begin, Range_Iterator end,
|
||||
Monge_info &monge_info);
|
||||
|
||||
//Coordinates of input points are computed in the fitting basis with
|
||||
// p0 as origin.
|
||||
//Preconditionning is computed, M and Z are filled
|
||||
void fill_matrix(Range_Iterator begin, Range_Iterator end,
|
||||
int d, LAMatrix& M, LAVector& Z);
|
||||
|
||||
//A is computed, solving MA=Z in the ls sense
|
||||
//Preconditionning is needed
|
||||
//the condition number of the matrix M is stored in monge_info
|
||||
void solve_linear_system(LAMatrix &M, LAVector &A, const LAVector &Z,
|
||||
Monge_info& monge_info);
|
||||
|
||||
//Classical differential geometric calculus
|
||||
//change_fitting2monge is computed
|
||||
//if deg_monge =1 only 1st order info
|
||||
//if deg_monge >= 2 2nd order info are computed
|
||||
void compute_Monge_basis(const LAVector &A, Monge_rep& monge_rep);
|
||||
|
||||
//if deg_monge >=3 then 3rd (and 4th) order info are computed
|
||||
void compute_Monge_coefficients( LAVector &A, int dprime,
|
||||
Monge_rep& monge_rep);
|
||||
|
||||
//for a trihedron (v1,v2,v3) switches v1 to -v1 if det(v1,v2,v3) < 0
|
||||
void switch_to_direct_orientation(LVector& v1, const LVector& v2,
|
||||
const LVector& v3);
|
||||
};
|
||||
|
||||
//-------------------------------------------------------------------
|
||||
// Implementation
|
||||
//------------------------------------------------------------------
|
||||
|
||||
template < class DataKernel, class LocalKernel, class LinAlgTraits>
|
||||
Monge_via_jet_fitting<DataKernel, LocalKernel, LinAlgTraits>::
|
||||
Monge_via_jet_fitting(Range_Iterator begin, Range_Iterator end,
|
||||
int d, int dprime,
|
||||
Monge_rep& monge_rep,
|
||||
Monge_info& monge_info)
|
||||
{
|
||||
// precondition: on the degrees, jet and monge
|
||||
CGAL_precondition( (d >=1) && (dprime >= 1)
|
||||
&& (dprime <= 4) && (dprime <= d) );
|
||||
this->deg = d;
|
||||
this->deg_monge = dprime;
|
||||
this->nb_d_jet_coeff = (d+1)*(d+2)/2;
|
||||
this->nb_input_pts = end - begin;
|
||||
// precondition: solvable ls system
|
||||
CGAL_precondition( nb_input_pts >= nb_d_jet_coeff );
|
||||
|
||||
//Initialize
|
||||
monge_rep.set_up(dprime);
|
||||
//for the system MA=Z
|
||||
LAMatrix M(nb_input_pts, nb_d_jet_coeff);
|
||||
LAVector A(nb_d_jet_coeff);
|
||||
LAVector Z(nb_input_pts);
|
||||
|
||||
compute_PCA(begin, end, monge_info);
|
||||
fill_matrix(begin, end, d, M, Z);//with precond
|
||||
solve_linear_system(M, A, Z, monge_info); //correct with precond
|
||||
compute_Monge_basis(A, monge_rep);
|
||||
if ( dprime >= 3) compute_Monge_coefficients(A, dprime, monge_rep);
|
||||
}
|
||||
|
||||
template < class DataKernel, class LocalKernel, class LinAlgTraits>
|
||||
void Monge_via_jet_fitting<DataKernel, LocalKernel, LinAlgTraits>::
|
||||
compute_PCA(Range_Iterator begin, Range_Iterator end,
|
||||
Monge_info &monge_info)
|
||||
{
|
||||
LAMatrix Cov(3,3);
|
||||
LAVector eval(3);
|
||||
LAMatrix evec(3,3);
|
||||
|
||||
int n = this->nb_input_pts;
|
||||
LFT x, y, z,
|
||||
sumX = 0., sumY = 0., sumZ = 0.,
|
||||
sumX2 = 0., sumY2 = 0., sumZ2 = 0.,
|
||||
sumXY = 0., sumXZ = 0., sumYZ = 0.,
|
||||
xx, yy, zz, xy, xz, yz;
|
||||
|
||||
for (; begin != end; begin++)
|
||||
{
|
||||
x = (*begin).x();
|
||||
y = (*begin).y();
|
||||
z = (*begin).z();
|
||||
sumX += x / n;
|
||||
sumY += y / n;
|
||||
sumZ += z / n;
|
||||
sumX2 += x * x / n;
|
||||
sumY2 += y * y / n;
|
||||
sumZ2 += z * z / n;
|
||||
sumXY += x * y / n;
|
||||
sumXZ += x * z / n;
|
||||
sumYZ += y * z / n;
|
||||
}
|
||||
xx = sumX2 - sumX * sumX;
|
||||
yy = sumY2 - sumY * sumY;
|
||||
zz = sumZ2 - sumZ * sumZ;
|
||||
xy = sumXY - sumX * sumY;
|
||||
xz = sumXZ - sumX * sumZ;
|
||||
yz = sumYZ - sumY * sumZ;
|
||||
Cov[0][0] = xx;
|
||||
Cov[0][1] = xy;
|
||||
Cov[0][2] = xz;
|
||||
Cov[1][0] = xy;
|
||||
Cov[1][1] = yy;
|
||||
Cov[1][2] = yz;
|
||||
Cov[2][0] = xz;
|
||||
Cov[2][1] = yz;
|
||||
Cov[2][2] = zz;
|
||||
// solve for eigenvalues and eigenvectors.
|
||||
// eigen values are sorted in descending order,
|
||||
// eigen vectors are sorted in accordance.
|
||||
LinAlgTraits::eigen_symm_algo(Cov, eval, evec);
|
||||
|
||||
//store in monge_info
|
||||
for (int i=0; i<3; i++)
|
||||
{
|
||||
monge_info.pca_eigen_vals()[i] = eval[i];//implicit cast LAFT->LFT
|
||||
LVector temp_vect(evec[0][i],evec[1][i],evec[2][i]);
|
||||
monge_info.pca_eigen_vecs()[i] = temp_vect;
|
||||
}
|
||||
|
||||
// //ALTERNATIVE with CGAL eigen code
|
||||
// // assemble covariance matrix as a
|
||||
// // semi-definite matrix.
|
||||
// // Matrix numbering:
|
||||
// // 0
|
||||
// // 1 2
|
||||
// // 3 4 5
|
||||
// LFT covariance[6] = {xx,xy,yy,xz,yz,zz};
|
||||
|
||||
// // solve for eigenvalues and eigenvectors.
|
||||
// // eigen values are sorted in descending order,
|
||||
// // eigen vectors are sorted in accordance.
|
||||
// LFT eigen_values[3];
|
||||
// LFT eigen_vectors[9];
|
||||
// CGAL::CGALi::eigen_symmetric<LFT>(covariance,3,eigen_vectors,eigen_values);
|
||||
// //store in monge_info
|
||||
// for (int i=0; i<3; i++)
|
||||
// {
|
||||
// monge_info.pca_eigen_vals()[i] = eigen_values[i];//implicit cast LAFT->LFT
|
||||
// }
|
||||
// LVector v1(eigen_vectors[0],eigen_vectors[1],eigen_vectors[2]);
|
||||
// monge_info.pca_eigen_vecs()[0] = v1;
|
||||
// LVector v2(eigen_vectors[3],eigen_vectors[4],eigen_vectors[5]);
|
||||
// monge_info.pca_eigen_vecs()[1] = v2;
|
||||
// LVector v3(eigen_vectors[6],eigen_vectors[7],eigen_vectors[8]);
|
||||
// monge_info.pca_eigen_vecs()[2] = v3;
|
||||
// ///end ALTERNATIVE with CGAL eigen code
|
||||
|
||||
switch_to_direct_orientation(monge_info.pca_eigen_vecs()[0],
|
||||
monge_info.pca_eigen_vecs()[1],
|
||||
monge_info.pca_eigen_vecs()[2]);
|
||||
|
||||
//Store the change of basis W->F
|
||||
const LVector* pca_vecs = monge_info.pca_eigen_vecs();
|
||||
Aff_transformation
|
||||
change_basis (pca_vecs[0][0], pca_vecs[0][1], pca_vecs[0][2],
|
||||
pca_vecs[1][0], pca_vecs[1][1], pca_vecs[1][2],
|
||||
pca_vecs[2][0], pca_vecs[2][1], pca_vecs[2][2]);
|
||||
this->change_world2fitting = change_basis;
|
||||
|
||||
/* //debug //test the old method, fitting basis is a permutation of the world basis */
|
||||
/* const LVector* pca_vecs = monge_info.pca_eigen_vecs(); */
|
||||
/* const LVector n_pca = pca_vecs[2]; */
|
||||
/* int index_max=0; */
|
||||
/* x = std::fabs(n_pca[0]); y = std::fabs(n_pca[1]); z = std::fabs(n_pca[2]); */
|
||||
/* if (x>y) if (x>z) index_max = 0; else index_max = 2; */
|
||||
/* else if (y>z) index_max = 1; else index_max = 2; */
|
||||
/* Aff_transformation change_basis; */
|
||||
/* if (index_max == 0) change_basis = Aff_transformation(0,1,0, */
|
||||
/* 0,0,1, */
|
||||
/* 1,0,0); */
|
||||
/* if (index_max == 1) change_basis = Aff_transformation(0,0,1, */
|
||||
/* 1,0,0, */
|
||||
/* 0,1,0); */
|
||||
/* if (index_max == 2) change_basis = Aff_transformation(1,0,0, */
|
||||
/* 0,1,0, */
|
||||
/* 0,0,1); */
|
||||
/* this->change_world2fitting = change_basis; */
|
||||
//test the old method END
|
||||
}
|
||||
|
||||
template < class DataKernel, class LocalKernel, class LinAlgTraits>
|
||||
void Monge_via_jet_fitting<DataKernel, LocalKernel, LinAlgTraits>::
|
||||
fill_matrix(Range_Iterator begin, Range_Iterator end,
|
||||
int d, LAMatrix &M, LAVector &Z)
|
||||
{
|
||||
//origin of fitting coord system = first input data point
|
||||
LPoint point0 = *begin;
|
||||
//transform coordinates of sample points with a
|
||||
//translation ($-p$) and multiplication by $ P_{W\rightarrow F}$.
|
||||
LPoint orig(0.,0.,0.);
|
||||
LVector v_point0_orig(orig - point0);
|
||||
Aff_transformation transl(CGAL::TRANSLATION, v_point0_orig);
|
||||
this->translate_p0 = transl;
|
||||
Aff_transformation transf_points = this->change_world2fitting *
|
||||
this->translate_p0;
|
||||
|
||||
//compute and store transformed points
|
||||
std::vector<LPoint> pts_in_fitting_basis;
|
||||
CGAL_For_all(begin,end){//implicit cast DPoint->LPoint
|
||||
LPoint cur_pt = transf_points(*begin);
|
||||
pts_in_fitting_basis.push_back(cur_pt);
|
||||
}
|
||||
|
||||
//Compute preconditionning
|
||||
LFT precond = 0.;
|
||||
typename std::vector<LPoint>::iterator itb = pts_in_fitting_basis.begin(),
|
||||
ite = pts_in_fitting_basis.end();
|
||||
CGAL_For_all(itb,ite) precond += std::fabs(itb->x()) + std::fabs(itb->y());
|
||||
precond /= 2*this->nb_input_pts;
|
||||
this->preconditionning = precond;
|
||||
//fill matrices M and Z
|
||||
itb = pts_in_fitting_basis.begin();
|
||||
int line_count = 0;
|
||||
LFT x, y;
|
||||
CGAL_For_all(itb,ite) {
|
||||
x = itb->x();
|
||||
y = itb->y();
|
||||
Z[line_count] = itb->z();
|
||||
for (int k=0; k <= d; k++) for (int i=0; i<=k; i++)
|
||||
M[line_count][k*(k+1)/2+i] =
|
||||
std::pow(x,k-i)*std::pow(y,i)
|
||||
/(fact(i)*fact(k-i)*std::pow(this->preconditionning,k));
|
||||
line_count++;
|
||||
}
|
||||
}
|
||||
|
||||
template < class DataKernel, class LocalKernel, class LinAlgTraits>
|
||||
void Monge_via_jet_fitting<DataKernel, LocalKernel, LinAlgTraits>::
|
||||
solve_linear_system(LAMatrix &M, LAVector &A, const LAVector &Z,
|
||||
Monge_info& monge_info)
|
||||
{
|
||||
LinAlgTraits::solve_ls_svd_algo(M, A, Z, monge_info.cond_nb());
|
||||
for (int k=0; k <= this->deg; k++) for (int i=0; i<=k; i++)
|
||||
A[k*(k+1)/2+i] /= std::pow(this->preconditionning,k);
|
||||
}
|
||||
|
||||
template < class DataKernel, class LocalKernel, class LinAlgTraits>
|
||||
void Monge_via_jet_fitting<DataKernel, LocalKernel, LinAlgTraits>::
|
||||
compute_Monge_basis(const LAVector &A, Monge_rep& monge_rep)
|
||||
{
|
||||
// only 1st order info.
|
||||
if ( this->deg_monge == 1 ) {
|
||||
LPoint orig_monge(0., 0., A[0]);
|
||||
LVector normal(-A[1], -A[2], 1.);
|
||||
LFT norm2 = normal * normal;
|
||||
normal = normal / Lsqrt(norm2);
|
||||
monge_rep.origin_pt() =
|
||||
(this->translate_p0.inverse() *
|
||||
this->change_world2fitting.inverse()) (orig_monge );
|
||||
monge_rep.n() = this->change_world2fitting.inverse()(normal);
|
||||
}
|
||||
// else (deg_monge >= 2) then 2nd order info are computed
|
||||
else {
|
||||
//bi-index to uni-index conversion : A(i,j)=A[(i+j)(i+j+1)/2+j]
|
||||
LPoint orig_monge(0., 0., A[0]);
|
||||
//normal = Xu crossprod Xv
|
||||
LVector Xu(1.,0.,A[1]), Xv(0.,1.,A[2]), normal(-A[1], -A[2], 1.);
|
||||
LFT norm2 = normal * normal;
|
||||
normal = normal / Lsqrt(norm2);
|
||||
|
||||
//Surface in fitting_basis : X(u,v)=(u,v,J_A(u,v))
|
||||
//in the basis Xu=(1,0,A[1]), Xv=(0,1,A[2]), Weingarten=-I^{-1}II
|
||||
//first fond form I=(e,f,f,g)
|
||||
// =(Xu.Xu, Xu.Xv, Xu.Xv, Xv.Xv)
|
||||
//second fond form II=(l,m,m,n)/norm2^(1/2)
|
||||
// =(n.Xuu, n.Xuv, n.Xuv, n.Xvv)
|
||||
//ppal curv are the opposite of the eigenvalues of Weingarten or the
|
||||
// eigenvalues of weingarten = -Weingarten = I^{-1}II
|
||||
typedef typename CGAL::Linear_algebraCd<LFT>::Matrix Matrix;
|
||||
|
||||
LFT e = 1+A[1]*A[1], f = A[1]*A[2], g = 1+A[2]*A[2],
|
||||
l = A[3], m = A[4], n = A[5];
|
||||
Matrix weingarten(2,2,0.);
|
||||
weingarten(0,0) = (g*l-f*m)/ (Lsqrt(norm2)*norm2);
|
||||
weingarten(0,1) = (g*m-f*n)/ (Lsqrt(norm2)*norm2);
|
||||
weingarten(1,0) = (e*m-f*l)/ (Lsqrt(norm2)*norm2);
|
||||
weingarten(1,1) = (e*n-f*m)/ (Lsqrt(norm2)*norm2);
|
||||
// Y, Z are normalized GramSchmidt of Xu, Xv
|
||||
// Xu->Y=Xu/||Xu||;
|
||||
// Xv->Z=Xv-(Xu.Xv)Xu/||Xu||^2;
|
||||
// Z-> Z/||Z||
|
||||
LVector Y, Z;
|
||||
LFT normXu = Lsqrt( Xu*Xu );
|
||||
Y = Xu / normXu;
|
||||
LFT XudotXv = Xu * Xv;
|
||||
Z = Xv - XudotXv * Xu / (normXu*normXu);
|
||||
LFT normZ = Lsqrt( Z*Z );
|
||||
Z = Z / normZ;
|
||||
Matrix change_XuXv2YZ(2,2,0.);
|
||||
change_XuXv2YZ(0,0) = 1 / normXu;
|
||||
change_XuXv2YZ(0,1) = -XudotXv / (normXu * normXu * normZ);
|
||||
change_XuXv2YZ(1,0) = 0;
|
||||
change_XuXv2YZ(1,1) = 1 / normZ;
|
||||
LFT det = 0.;
|
||||
Matrix inv = CGAL::Linear_algebraCd<LFT>::inverse ( change_XuXv2YZ, det );
|
||||
//in the new orthonormal basis (Y,Z) of the tangent plane :
|
||||
weingarten = inv *(1/det) * weingarten * change_XuXv2YZ;
|
||||
|
||||
//switch to LinAlgTraits for diagonalization of weingarten
|
||||
LAMatrix W(2,2);
|
||||
for (int i=0; i<=1; i++) for (int j=0; j<=1; j++)
|
||||
W[i][j] = weingarten(i,j);
|
||||
LAVector eval(2);
|
||||
LAMatrix evec(2,2);
|
||||
|
||||
LinAlgTraits::eigen_symm_algo(W, eval, evec);
|
||||
LVector d_max = evec[0][0]*Y + evec[1][0]*Z,
|
||||
d_min = evec[0][1]*Y + evec[1][1]*Z;
|
||||
|
||||
switch_to_direct_orientation(d_max, d_min, normal);
|
||||
Aff_transformation change_basis (d_max[0], d_max[1], d_max[2],
|
||||
d_min[0], d_min[1], d_min[2],
|
||||
normal[0], normal[1], normal[2]);
|
||||
this->change_fitting2monge = change_basis;
|
||||
|
||||
//store the monge basis origin and vectors with their world coord
|
||||
//store ppal curv
|
||||
monge_rep.origin_pt() =
|
||||
(this->translate_p0.inverse() *
|
||||
this->change_world2fitting.inverse()) (orig_monge );
|
||||
monge_rep.d1() = this->change_world2fitting.inverse()(d_max);
|
||||
monge_rep.d2() = this->change_world2fitting.inverse()(d_min);
|
||||
monge_rep.n() = this->change_world2fitting.inverse()(normal);
|
||||
monge_rep.coefficients()[0] = eval[0];
|
||||
monge_rep.coefficients()[1] = eval[1];
|
||||
}
|
||||
//end else
|
||||
}
|
||||
|
||||
template < class DataKernel, class LocalKernel, class LinAlgTraits>
|
||||
void Monge_via_jet_fitting<DataKernel, LocalKernel, LinAlgTraits>::
|
||||
compute_Monge_coefficients( LAVector &A, int dprime,
|
||||
Monge_rep& monge_rep)
|
||||
{
|
||||
//One has the equation w=J_A(u,v) of the fitted surface S
|
||||
// in the fitting_basis
|
||||
//Substituing (u,v,w)=change_fitting2monge^{-1}(x,y,z)
|
||||
//One has the equation f(x,y,z)=0 on this surface S in the monge
|
||||
// basis
|
||||
//The monge form of the surface at the origin is the bivariate fct
|
||||
// g(x,y) s.t. f(x,y,g(x,y))=0
|
||||
//voir les calculs Maple dans monge.mws
|
||||
//Notations are f123= d^3f/dxdydz
|
||||
// g(x,y)=sum (gij x^i y^j/ i!j!) with
|
||||
// g00=g10=g01=g11=0, g20=kmax, g02=kmin
|
||||
//
|
||||
//g(x,y)= 1/2*(k1x^2 +k2y^2)
|
||||
// +1/6*(b0x^3 +3b1x^2y +3b2xy^2 +b3y^3)
|
||||
// +1/24*(c0x^4 +4c1x^3y +6c2x^2y^2 +4c3xy^3 +c4y^4)
|
||||
// +...
|
||||
// p stores change_fitting2monge^{-1}=change_fitting2monge^{T}
|
||||
LFT p[3][3];
|
||||
p[0][0] = this->change_fitting2monge.m(0,0);
|
||||
p[1][0] = this->change_fitting2monge.m(0,1);
|
||||
p[2][0] = this->change_fitting2monge.m(0,2);
|
||||
p[0][1] = this->change_fitting2monge.m(1,0);
|
||||
p[1][1] = this->change_fitting2monge.m(1,1);
|
||||
p[2][1] = this->change_fitting2monge.m(1,2);
|
||||
p[0][2] = this->change_fitting2monge.m(2,0);
|
||||
p[1][2] = this->change_fitting2monge.m(2,1);
|
||||
p[2][2] = this->change_fitting2monge.m(2,2);
|
||||
|
||||
// formula are designed for w=sum( Aij ui vj), but we have J_A = sum( Aij/i!j! ui vj)
|
||||
for (int k=0; k <= this->deg; k++) for (int i=0; i<=k; i++)
|
||||
A[k*(k+1)/2+i] /= fact(k-i)*fact(i);//this is A(k-i;i)
|
||||
|
||||
/* //debug */
|
||||
/* std::cout << "coeff of A" << std::endl */
|
||||
/* << A[0] << " "<< A[1] << " "<< A[2] << std::endl */
|
||||
/* << A[3] << " "<< A[4] << " "<< A[5] << std::endl */
|
||||
/* << A[6] << " "<< A[7] << " "<< A[8] << " "<< A[9]<< std::endl */
|
||||
/* << A[10] << " "<< A[11] << " "<< A[12] << " "<< A[13]<< " " << A[14] << std::endl; */
|
||||
|
||||
|
||||
|
||||
// note f1 = f2 = f12 = 0
|
||||
// LFT f1 = A[1] * p[0][0] + A[2] * p[1][0] - p[2][0];
|
||||
// LFT f2 = A[2] * p[1][1] + A[1] * p[0][1] - p[2][1];
|
||||
// LFT f12 =
|
||||
// 2 * A[3] * p[0][0] * p[0][1]
|
||||
// + 2 * A[5] * p[1][0] * p[1][1]
|
||||
// + A[4] * p[0][1] * p[1][0]
|
||||
// + A[4] * p[0][0] * p[1][1];
|
||||
// -f11 / f3 = kmax
|
||||
// -f22 / f3 = kmin
|
||||
|
||||
LFT f3 = A[1] * p[0][2] + A[2] * p[1][2] - p[2][2];
|
||||
LFT f11 =
|
||||
2 * A[4] * p[0][0] * p[1][0]
|
||||
+ 2 * A[5] * p[1][0] * p[1][0]
|
||||
+ 2 * A[3] * p[0][0] * p[0][0];
|
||||
LFT f13 =
|
||||
A[4] * p[0][0] * p[1][2]
|
||||
+ A[4] * p[0][2] * p[1][0]
|
||||
+ 2 * A[5] * p[1][0] * p[1][2]
|
||||
+ 2 * A[3] * p[0][0] * p[0][2];
|
||||
LFT f22 =
|
||||
2 * A[4] * p[0][1] * p[1][1]
|
||||
+ 2 * A[5] * p[1][1] * p[1][1]
|
||||
+ 2 * A[3] * p[0][1] * p[0][1];
|
||||
LFT f23 =
|
||||
A[4] * p[0][1] * p[1][2]
|
||||
+ 2 * A[5] * p[1][1] * p[1][2]
|
||||
+ A[4] * p[0][2] * p[1][1]
|
||||
+ 2 * A[3] * p[0][1] * p[0][2];
|
||||
LFT f33 =
|
||||
2 * A[5] * p[1][2] * p[1][2]
|
||||
+ 2 * A[3] * p[0][2] * p[0][2]
|
||||
+ 2 * A[4] * p[0][2] * p[1][2];
|
||||
LFT f111 =
|
||||
6 * A[8] * p[0][0] * p[1][0] * p[1][0]
|
||||
+ 6 * A[7] * p[0][0] * p[0][0] * p[1][0]
|
||||
+ 6 * A[6] * p[0][0] * p[0][0] * p[0][0]
|
||||
+ 6 * A[9] * p[1][0] * p[1][0] * p[1][0];
|
||||
LFT f222 =
|
||||
6 * A[7] * p[0][1] * p[0][1] * p[1][1]
|
||||
+ 6 * A[8] * p[0][1] * p[1][1] * p[1][1]
|
||||
+ 6 * A[9] * p[1][1] * p[1][1] * p[1][1]
|
||||
+ 6 * A[6] * p[0][1] * p[0][1] * p[0][1];
|
||||
LFT f112 =
|
||||
2 * A[7] * p[0][0] * p[0][0] * p[1][1]
|
||||
+ 6 * A[6] * p[0][0] * p[0][0] * p[0][1]
|
||||
+ 2 * A[8] * p[0][1] * p[1][0] * p[1][0]
|
||||
+ 4 * A[8] * p[0][0] * p[1][0] * p[1][1]
|
||||
+ 6 * A[9] * p[1][0] * p[1][0] * p[1][1]
|
||||
+ 4 * A[7] * p[0][0] * p[0][1] * p[1][0];
|
||||
LFT f122 =
|
||||
4 * A[8] * p[0][1] * p[1][0] * p[1][1]
|
||||
+ 2 * A[8] * p[0][0] * p[1][1] * p[1][1]
|
||||
+ 6 * A[6] * p[0][0] * p[0][1] * p[0][1]
|
||||
+ 2 * A[7] * p[0][1] * p[0][1] * p[1][0]
|
||||
+ 4 * A[7] * p[0][0] * p[0][1] * p[1][1]
|
||||
+ 6 * A[9] * p[1][0] * p[1][1] * p[1][1];
|
||||
LFT f113 =
|
||||
6*A[6]*p[0][0]*p[0][0]*p[0][2]
|
||||
+6*A[9]*p[1][0]*p[1][0]*p[1][2]
|
||||
+2*A[7]*p[0][0]*p[0][0]*p[1][2]
|
||||
+2*A[8]*p[0][2]*p[1][0]*p[1][0]
|
||||
+4*A[7]*p[0][0]*p[0][2]*p[1][0]
|
||||
+4*A[8]*p[0][0]*p[1][0]*p[1][2];
|
||||
LFT f223 =
|
||||
2*A[8]*p[0][2]*p[1][1]*p[1][1]
|
||||
+6*A[6]*p[0][1]*p[0][1]*p[0][2]
|
||||
+6*A[9]*p[1][1]*p[1][1]*p[1][2]
|
||||
+2*A[7]*p[0][1]*p[0][1]*p[1][2]
|
||||
+4*A[7]*p[0][1]*p[0][2]*p[1][1]
|
||||
+4*A[8]*p[0][1]*p[1][1]*p[1][2];
|
||||
LFT f123 =
|
||||
2*A[8]*p[0][2]*p[1][0]*p[1][1]
|
||||
+2*A[7]*p[0][0]*p[0][1]*p[1][2]
|
||||
+2*A[7]*p[0][0]*p[0][2]*p[1][1]
|
||||
+6*A[9]*p[1][0]*p[1][1]*p[1][2]
|
||||
+2*A[7]*p[0][1]*p[0][2]*p[1][0]
|
||||
+6*A[6]*p[0][0]*p[0][1]*p[0][2]
|
||||
+2*A[8]*p[0][0]*p[1][1]*p[1][2]
|
||||
+2*A[8]*p[0][1]*p[1][0]*p[1][2];
|
||||
|
||||
LFT b0 = 1/(f3*f3)*(-f111*f3+3*f13*f11);
|
||||
LFT b1 = 1/(f3*f3)*(-f112*f3+f23*f11);
|
||||
LFT b2 = 1/(f3*f3)*(-f122*f3+f13*f22);
|
||||
LFT b3 = -1/(f3*f3)*(f222*f3-3*f23*f22);
|
||||
|
||||
monge_rep.coefficients()[2] = b0;
|
||||
monge_rep.coefficients()[3] = b1;
|
||||
monge_rep.coefficients()[4] = b2;
|
||||
monge_rep.coefficients()[5] = b3;
|
||||
|
||||
if ( dprime == 4 )
|
||||
{
|
||||
LFT f1111 =
|
||||
24*A[13]*p[0][0]*p[1][0]*p[1][0]*p[1][0]
|
||||
+24*A[12]*p[0][0]*p[0][0]*p[1][0]*p[1][0]
|
||||
+24*A[11]*p[0][0]*p[0][0]*p[0][0]*p[1][0]
|
||||
+24*A[14]*p[1][0]*p[1][0]*p[1][0]*p[1][0]
|
||||
+24*A[10]*p[0][0]*p[0][0]*p[0][0]*p[0][0];
|
||||
LFT f1112 =
|
||||
6*A[13]*p[0][1]*p[1][0]*p[1][0]*p[1][0]
|
||||
+18*A[13]*p[0][0]*p[1][0]*p[1][0]*p[1][1]
|
||||
+24*A[10]*p[0][0]*p[0][0]*p[0][0]*p[0][1]
|
||||
+12*A[12]*p[0][0]*p[0][1]*p[1][0]*p[1][0]
|
||||
+18*A[11]*p[0][0]*p[0][0]*p[0][1]*p[1][0]
|
||||
+24*A[14]*p[1][0]*p[1][0]*p[1][0]*p[1][1]
|
||||
+6*A[11]*p[0][0]*p[0][0]*p[0][0]*p[1][1]
|
||||
+12*A[12]*p[0][0]*p[0][0]*p[1][0]*p[1][1];
|
||||
LFT f1122 =
|
||||
12*A[11]*p[0][0]*p[0][0]*p[0][1]*p[1][1]
|
||||
+12*A[13]*p[0][0]*p[1][0]*p[1][1]*p[1][1]
|
||||
+12*A[13]*p[0][1]*p[1][0]*p[1][0]*p[1][1]
|
||||
+16*A[12]*p[0][0]*p[0][1]*p[1][0]*p[1][1]
|
||||
+12*A[11]*p[0][0]*p[0][1]*p[0][1]*p[1][0]
|
||||
+24*A[10]*p[0][0]*p[0][0]*p[0][1]*p[0][1]
|
||||
+4*A[12]*p[0][1]*p[0][1]*p[1][0]*p[1][0]
|
||||
+4*A[12]*p[0][0]*p[0][0]*p[1][1]*p[1][1]
|
||||
+24*A[14]*p[1][0]*p[1][0]*p[1][1]*p[1][1];
|
||||
LFT f1222 =
|
||||
6*A[13]*p[0][0]*p[1][1]*p[1][1]*p[1][1]
|
||||
+24*A[10]*p[0][0]*p[0][1]*p[0][1]*p[0][1]
|
||||
+24*A[14]*p[1][0]*p[1][1]*p[1][1]*p[1][1]
|
||||
+6*A[11]*p[0][1]*p[0][1]*p[0][1]*p[1][0]
|
||||
+18*A[11]*p[0][0]*p[0][1]*p[0][1]*p[1][1]
|
||||
+12*A[12]*p[0][0]*p[0][1]*p[1][1]*p[1][1]
|
||||
+12*A[12]*p[0][1]*p[0][1]*p[1][0]*p[1][1]
|
||||
+18*A[13]*p[0][1]*p[1][0]*p[1][1]*p[1][1];
|
||||
LFT f2222 =
|
||||
24*A[13]*p[0][1]*p[1][1]*p[1][1]*p[1][1]
|
||||
+24*A[11]*p[0][1]*p[0][1]*p[0][1]*p[1][1]
|
||||
+24*A[12]*p[0][1]*p[0][1]*p[1][1]*p[1][1]
|
||||
+24*A[10]*p[0][1]*p[0][1]*p[0][1]*p[0][1]
|
||||
+24*A[14]*p[1][1]*p[1][1]*p[1][1]*p[1][1];
|
||||
|
||||
LFT c0 =
|
||||
-1/(f3*f3*f3)*(f1111*(f3*f3)-4*f13*f3*f111+12*f13*f13*f11-6*f113*f3*f11+3*f33*f11*f11);
|
||||
LFT c1 =
|
||||
1/(f3*f3*f3)*(f23*f3*f111+3*f3*f123*f11+3*f13*f3*f112-f1112*(f3*f3)-6*f13*f23*f11);
|
||||
LFT c2 =
|
||||
1/(f3*f3*f3)*(-f33*f22*f11+f113*f3*f22+2*f13*f3*f122-2*f13*f13*f22+f223*f3*f11+2*f23*f3*f112-2*f23*f23*f11-f1122*(f3*f3));
|
||||
LFT c3 =
|
||||
1/(f3*f3*f3)*(-f1222*(f3*f3)-6*f13*f23*f22+3*f123*f3*f22+f13*f3*f222+3*f23*f3*f122);
|
||||
LFT c4 =
|
||||
-1/(f3*f3*f3)*(f2222*(f3*f3)+3*f33*f22*f22-6*f223*f3*f22-4*f23*f3*f222+12*f23*f23*f22) ;
|
||||
|
||||
monge_rep.coefficients()[6] = c0;
|
||||
monge_rep.coefficients()[7] = c1;
|
||||
monge_rep.coefficients()[8] = c2;
|
||||
monge_rep.coefficients()[9] = c3;
|
||||
monge_rep.coefficients()[10] = c4;
|
||||
}
|
||||
}
|
||||
|
||||
template < class DataKernel, class LocalKernel, class LinAlgTraits>
|
||||
void Monge_via_jet_fitting<DataKernel, LocalKernel, LinAlgTraits>::
|
||||
switch_to_direct_orientation(LVector& v1, const LVector& v2,
|
||||
const LVector& v3)
|
||||
{
|
||||
CGAL::Sign orientation = CGAL::sign_of_determinant3x3(v1[0], v2[0], v3[0],
|
||||
v1[1], v2[1], v3[1],
|
||||
v1[2], v2[2], v3[2]);
|
||||
if (orientation == CGAL::NEGATIVE) v1 = -v1;
|
||||
}
|
||||
|
||||
CGAL_END_NAMESPACE
|
||||
|
||||
#endif //_MONGE_VIA_JET_FITTING_H_
|
||||
|
||||
|
||||
|
||||
Loading…
Reference in New Issue