instrospect-> introspect, minor changes

This commit is contained in:
Marc Pouget 2006-05-26 14:32:20 +00:00
parent e7d3c22cb9
commit 44264f1d9f
9 changed files with 360 additions and 70 deletions

2
.gitattributes vendored
View File

@ -449,6 +449,8 @@ 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
Jet_fitting_3/examples/Jet_fitting_3/test.off -text
Jet_fitting_3/test/Jet_fitting_3/blind_1pt.C -text
Jet_fitting_3/test/Jet_fitting_3/makefile -text
Kernel_23/doc_tex/Kernel_23_ref/fig/IsoCuboid.eps -text
Kernel_23/doc_tex/Kernel_23_ref/fig/IsoCuboid.gif -text svneol=unset#unset
Kernel_23/doc_tex/Kernel_23_ref/fig/IsoCuboid.pdf -text svneol=unset#unset

View File

@ -46,7 +46,7 @@ QTLDFLAGS = \
$(CGAL_QT_LDFLAGS) -lGL -lglut
INTRO_LD_FLAGS=-L$(INTROSPEC_DIR) -lInstrospect
INTRO_LD_FLAGS=-L$(INTROSPEC_DIR) -lIntrospect
#---------------------------------------------------------------------#
# target entries

View File

@ -1,7 +1,7 @@
\ccUserChapter{Estimation of local differential properties of sampled
surfaces via polynomial fitting}
\label{chap:Jet_fitting_3}
\ccChapterAuthor{Marc Pouget and Frédéric Cazals}
\ccChapterAuthor{Marc Pouget and Frederic Cazals}
\minitoc

View File

@ -35,8 +35,6 @@ the \ccc{LocalKernel} concept~: \ccc{LocalKernel::FT}.
% +------------------------------------------------------------------
\ccNestedType{FT}{The scalar type.}
\ccGlue
\ccNestedType{Vector}{The Vector type.}
\ccGlue
\ccNestedType{Matrix }{The Matrix type.}
%\ccCreation
@ -45,22 +43,28 @@ the \ccc{LocalKernel} concept~: \ccc{LocalKernel::FT}.
%\ccConstructor{LinAlgTraits();}{default constructor.}
\ccOperations
The Matrix has classical access to its elements.
Usual brackets operator for Vector and Matrix types.
\ccMethod{void set_elt(size_t i, size_t j, const FT value)}{}
\ccGlue
\ccMethod{FT get_elt(size_t i, size_t j)}{}
\ccMethod{void eigen_symm_algo(Matrix& S, Vector& eval, Matrix&
The LinAlgTraits has an eigenanalysis and a singular value
decomposition algorithm.
\ccMethod{void eigen_symm_algo(Matrix& S, FT* eval, Matrix&
evec);} {Performs an eigenanalysis of a real symmetric matrix. Eigen
values are sorted in descending order, eigen vectors are sorted in
values are sorted in ascending order, eigen vectors are sorted in
accordance. }
\ccMethod{void solve_ls_svd_algo(Matrix& M, Vector& X, const Vector&
B, double &cond_nb);}{ Solves the system MX=B (in the least square
sense if M is not square) using a Singular Value Decomposition and
gives the condition number of M.}
\ccMethod{void solve_ls_svd_algo(Matrix& M, FT* B, FT* cond_nb);}
{ Solves the system MX=B (in the least square sense if M is not
square) using a Singular Value Decomposition and gives the condition
number of M. The solution is stored in B.}
\ccHasModels
% +------------------------------------------------------------------
\ccc{GSL}. See the implementation in the exemple directory
exemples/Jet\_fitting\_3/GSL.h .
\ccc{Lapack}. See the implementation in the exemple directory
exemples/Jet\_fitting\_3/LinAlg\_lapack.h .

View File

@ -81,65 +81,16 @@ int main(int argc, char *argv[])
//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;
if ( d_monge >= 3)
outFile << "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)
outFile << "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;
<< "number of points used : " << in_points.size() << std::endl;
monge_rep.dump_verbose(outFile);
monge_info.dump_verbose(outFile);
//OUTPUT on std::cout
//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;
<< "number of points used : " << in_points.size() << std::endl;
monge_rep.dump_verbose(std::cout);
monge_info.dump_verbose(std::cout);
return 1;
}

View File

@ -0,0 +1,65 @@
#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.h"
#include "include/LinAlg_lapack.h"
typedef double DFT;
typedef CGAL::Cartesian<DFT> Data_Kernel;
typedef Data_Kernel::Point_3 DPoint;
typedef Data_Kernel::Vector_3 DVector;
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, Lapack> My_Monge_via_jet_fitting;
int main()
{
//open the input file
std::ifstream inFile( "data/in_points_d4.txt", 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 = 4;
int d_monge = 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);
monge_rep.comply_wrt_given_normal( -monge_rep.n() );
//OUTPUT on std::cout
monge_rep.dump_verbose( std::cout );
monge_rep.dump_4ogl( std::cout, 1 );
monge_info.dump_verbose( std::cout );
std::cout << "success\n";
}

View File

@ -0,0 +1,17 @@
0 0 0
0 1 .2
0 -1 0.2
1 0 0.1
-1 0 0.1
1 1 0.3
-1 1 .3
-1 -1 .3
1 -1 .3
2 0 .4
0 2 .8
-2 0 .4
0 -2 .8
2 2 1.2
-2 -2 1.2
2 -2 1.2
-2 2 1.2

View File

@ -0,0 +1,110 @@
#ifndef _LAPACK_H_
#define _LAPACK_H_
#include <stdlib.h>
#include "blaswrap.h"
#include "f2c.h"
extern "C" {
#include "clapack.h"
}
////////////////////////class Lapack_matrix/////////////////////
//in Lapack, matrices are one-dimensional arrays
// and elements are column-major ordered
// this class is a wrapper defining set and get in the usual way
class Lapack_matrix{
protected:
double* m_matrix;
public:
size_t nb_rows;
size_t nb_columns;
//contructor
// initializes all the elements of the matrix to zero.
Lapack_matrix(size_t n1, size_t n2) {
m_matrix = (double*) calloc (n1*n2, sizeof(double));
nb_rows = n1;
nb_columns = n2;
}
//destructor
// ~Lapack_matrix();// { free m_matrix; }
//access
const double* matrix() const { return m_matrix; }
double* matrix() { return m_matrix; }
void set_elt(size_t i, size_t j, const double value) { m_matrix[j*nb_rows+i] = value; }
double get_elt(size_t i, size_t j) { return m_matrix[j*nb_rows+i]; }
};
//Lapack_matrix::~Lapack_matrix() { free m_matrix; } bug!
////////////////////////class Lapack/////////////////////
class Lapack{
public:
typedef Lapack_matrix Matrix;
// solve for eigenvalues and eigenvectors.
// eigen values are sorted in ascending order,
// eigen vectors are sorted in accordance.
static
void eigen_symm_algo(Matrix& S, double* eval, Matrix& evec);
//solve MX=B using SVD and give the condition number of M
static
void solve_ls_svd_algo(Matrix& M, double* B, double &cond_nb);
};
void Lapack::eigen_symm_algo(Matrix& S, double* eval, Matrix& evec)
{
char jobz = 'V';
char range = 'A';
char uplo = 'U';
integer n = S.nb_rows,
lda = S.nb_rows,
il = 0, iu = 0,
m = n,
lwork = 26*n,
liwork = 10*n,
info;
double vl=0, vu=0, abstol = 0;
integer* isuppz = (integer*) malloc(2*n*sizeof(integer));
double* work = (double*) malloc(lwork*sizeof(double));
integer* iwork = (integer*) malloc(liwork*sizeof(integer));
dsyevr_(&jobz, &range, &uplo, &n, S.matrix(), &lda, &vl, &vu, &il, &iu, &abstol, &m,
eval, evec.matrix(), &n, isuppz, work, &lwork, iwork, &liwork, &info);
//clean up
free(isuppz);
free(work);
free(iwork);
}
void Lapack::solve_ls_svd_algo(Matrix& M, double* B, double &cond_nb)
{
integer m = M.nb_rows,
n = M.nb_columns,
nrhs = 1,
lda = m,
ldb = m,
rank,
lwork = 5*m,
info;
//c style
double* sing_values = (double*) malloc(n*sizeof(double));
double* work = (double*) malloc(lwork*sizeof(double));
double rcond = -1;
dgelss_(&m, &n, &nrhs, M.matrix(), &lda, B, &ldb, sing_values,
&rcond, &rank, work, &lwork, &info);
assert(info==0);
cond_nb = sing_values[0]/sing_values[n-1];
//clean up
free(sing_values);
free(work);
}
#endif

View File

@ -0,0 +1,141 @@
# Created by the script create_makefile
# This is the makefile for compiling a CGAL application.
#---------------------------------------------------------------------#
# include platform specific settings
#---------------------------------------------------------------------#
# Choose the right include file from the <cgalroot>/make directory.
include $(CGAL_MAKEFILE)
PROFOPT= -g
CFLAGS=${PROFOPT}
#---------------------------------------------------------------------#
# compiler flags
#---------------------------------------------------------------------#
LAPACK_INCDIRS = -I$(LAPACK_DIR)/SRC -I$(LAPACK_DIR)
CXXFLAGS = \
-I../../include \
$(CGAL_CXXFLAGS) \
$(LONG_NAME_PROBLEM_CXXFLAGS) \
$(LAPACK_INCDIRS) \
${CFLAGS}
#---------------------------------------------------------------------#
# linker flags
#---------------------------------------------------------------------#
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
LIBPATH = \
$(CGAL_LIBPATH)
LDFLAGS = \
$(LONG_NAME_PROBLEM_LDFLAGS) \
${LAPACK_LDLIBS} \
$(CGAL_LDFLAGS) \
#---------------------------------------------------------------------#
# target entries
#---------------------------------------------------------------------#
all: \
blind_1pt
blind_1pt: blind_1pt.o
$(CGAL_CXX) $(LIBPATH) -o blind_1pt.exe blind_1pt.o \
$(LDFLAGS)
rmexe:
rm *.exe
clean: \
blind_1pt.clean
rmexe
depend:
makedepend *.[Ch]
#---------------------------------------------------------------------#
# suffix rules
#---------------------------------------------------------------------#
.C$(OBJ_EXT):
$(CGAL_CXX) $(CXXFLAGS) $(OBJ_OPT) $<
# DO NOT DELETE
blind_1pt.o: /usr/include/stdio.h /usr/include/features.h
blind_1pt.o: /usr/include/sys/cdefs.h /usr/include/gnu/stubs.h
blind_1pt.o: /usr/include/bits/types.h /usr/include/bits/wordsize.h
blind_1pt.o: /usr/include/bits/typesizes.h /usr/include/libio.h
blind_1pt.o: /usr/include/_G_config.h /usr/include/wchar.h
blind_1pt.o: /usr/include/bits/wchar.h /usr/include/gconv.h
blind_1pt.o: /usr/include/bits/stdio_lim.h /usr/include/bits/sys_errlist.h
blind_1pt.o: /usr/include/stdlib.h /usr/include/sys/types.h
blind_1pt.o: /usr/include/time.h /usr/include/endian.h
blind_1pt.o: /usr/include/bits/endian.h /usr/include/sys/select.h
blind_1pt.o: /usr/include/bits/select.h /usr/include/bits/sigset.h
blind_1pt.o: /usr/include/bits/time.h /usr/include/sys/sysmacros.h
blind_1pt.o: /usr/include/bits/pthreadtypes.h /usr/include/bits/sched.h
blind_1pt.o: /usr/include/alloca.h ../../include/CGAL/Monge_via_jet_fitting.h
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: include/GSL.h /usr/include/gsl/gsl_linalg.h
blind_1pt.o: /usr/include/gsl/gsl_mode.h /usr/include/gsl/gsl_permutation.h
blind_1pt.o: /usr/include/gsl/gsl_types.h /usr/include/gsl/gsl_errno.h
blind_1pt.o: /usr/include/errno.h /usr/include/bits/errno.h
blind_1pt.o: /usr/include/linux/errno.h /usr/include/asm/errno.h
blind_1pt.o: /usr/include/asm-i386/errno.h /usr/include/gsl/gsl_check_range.h
blind_1pt.o: /usr/include/gsl/gsl_vector.h
blind_1pt.o: /usr/include/gsl/gsl_vector_complex_long_double.h
blind_1pt.o: /usr/include/gsl/gsl_complex.h
blind_1pt.o: /usr/include/gsl/gsl_vector_long_double.h
blind_1pt.o: /usr/include/gsl/gsl_block_long_double.h
blind_1pt.o: /usr/include/gsl/gsl_vector_complex.h
blind_1pt.o: /usr/include/gsl/gsl_block_complex_long_double.h
blind_1pt.o: /usr/include/gsl/gsl_vector_complex_double.h
blind_1pt.o: /usr/include/gsl/gsl_vector_double.h
blind_1pt.o: /usr/include/gsl/gsl_block_double.h
blind_1pt.o: /usr/include/gsl/gsl_block_complex_double.h
blind_1pt.o: /usr/include/gsl/gsl_vector_complex_float.h
blind_1pt.o: /usr/include/gsl/gsl_vector_float.h
blind_1pt.o: /usr/include/gsl/gsl_block_float.h
blind_1pt.o: /usr/include/gsl/gsl_block_complex_float.h
blind_1pt.o: /usr/include/gsl/gsl_vector_ulong.h
blind_1pt.o: /usr/include/gsl/gsl_block_ulong.h
blind_1pt.o: /usr/include/gsl/gsl_vector_long.h
blind_1pt.o: /usr/include/gsl/gsl_block_long.h
blind_1pt.o: /usr/include/gsl/gsl_vector_uint.h
blind_1pt.o: /usr/include/gsl/gsl_block_uint.h
blind_1pt.o: /usr/include/gsl/gsl_vector_int.h
blind_1pt.o: /usr/include/gsl/gsl_block_int.h
blind_1pt.o: /usr/include/gsl/gsl_vector_ushort.h
blind_1pt.o: /usr/include/gsl/gsl_block_ushort.h
blind_1pt.o: /usr/include/gsl/gsl_vector_short.h
blind_1pt.o: /usr/include/gsl/gsl_block_short.h
blind_1pt.o: /usr/include/gsl/gsl_vector_uchar.h
blind_1pt.o: /usr/include/gsl/gsl_block_uchar.h
blind_1pt.o: /usr/include/gsl/gsl_vector_char.h
blind_1pt.o: /usr/include/gsl/gsl_block_char.h /usr/include/gsl/gsl_matrix.h
blind_1pt.o: /usr/include/gsl/gsl_matrix_complex_long_double.h
blind_1pt.o: /usr/include/gsl/gsl_matrix_complex_double.h
blind_1pt.o: /usr/include/gsl/gsl_matrix_complex_float.h
blind_1pt.o: /usr/include/gsl/gsl_matrix_long_double.h
blind_1pt.o: /usr/include/gsl/gsl_matrix_double.h
blind_1pt.o: /usr/include/gsl/gsl_matrix_float.h
blind_1pt.o: /usr/include/gsl/gsl_matrix_ulong.h
blind_1pt.o: /usr/include/gsl/gsl_matrix_long.h
blind_1pt.o: /usr/include/gsl/gsl_matrix_uint.h
blind_1pt.o: /usr/include/gsl/gsl_matrix_int.h
blind_1pt.o: /usr/include/gsl/gsl_matrix_ushort.h
blind_1pt.o: /usr/include/gsl/gsl_matrix_short.h
blind_1pt.o: /usr/include/gsl/gsl_matrix_uchar.h
blind_1pt.o: /usr/include/gsl/gsl_matrix_char.h /usr/include/gsl/gsl_eigen.h