diff --git a/.gitattributes b/.gitattributes index b98356860a3..e3ead25149a 100644 --- a/.gitattributes +++ b/.gitattributes @@ -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 diff --git a/Jet_fitting_3/demo/Jet_fitting_3/makefile b/Jet_fitting_3/demo/Jet_fitting_3/makefile index 06b499854be..6cf59e4050c 100644 --- a/Jet_fitting_3/demo/Jet_fitting_3/makefile +++ b/Jet_fitting_3/demo/Jet_fitting_3/makefile @@ -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 diff --git a/Jet_fitting_3/doc_tex/Jet_fitting_3/main.tex b/Jet_fitting_3/doc_tex/Jet_fitting_3/main.tex index f12b35c8545..5fd56d6f0c6 100644 --- a/Jet_fitting_3/doc_tex/Jet_fitting_3/main.tex +++ b/Jet_fitting_3/doc_tex/Jet_fitting_3/main.tex @@ -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 diff --git a/Jet_fitting_3/doc_tex/Jet_fitting_3_ref/LinAlgTraits.tex b/Jet_fitting_3/doc_tex/Jet_fitting_3_ref/LinAlgTraits.tex index 67d6e6c140a..203caca6e83 100644 --- a/Jet_fitting_3/doc_tex/Jet_fitting_3_ref/LinAlgTraits.tex +++ b/Jet_fitting_3/doc_tex/Jet_fitting_3_ref/LinAlgTraits.tex @@ -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 . diff --git a/Jet_fitting_3/examples/Jet_fitting_3/blind_1pt.C b/Jet_fitting_3/examples/Jet_fitting_3/blind_1pt.C index 123c09e5da6..1f4ef29d12f 100644 --- a/Jet_fitting_3/examples/Jet_fitting_3/blind_1pt.C +++ b/Jet_fitting_3/examples/Jet_fitting_3/blind_1pt.C @@ -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; } diff --git a/Jet_fitting_3/test/Jet_fitting_3/blind_1pt.C b/Jet_fitting_3/test/Jet_fitting_3/blind_1pt.C new file mode 100644 index 00000000000..3d4615d763c --- /dev/null +++ b/Jet_fitting_3/test/Jet_fitting_3/blind_1pt.C @@ -0,0 +1,65 @@ +#include +#include + +#include +#include +#include +#include + +#include +#include "../../include/CGAL/Monge_via_jet_fitting.h" +#include "include/LinAlg_lapack.h" + +typedef double DFT; +typedef CGAL::Cartesian Data_Kernel; +typedef Data_Kernel::Point_3 DPoint; +typedef Data_Kernel::Vector_3 DVector; +typedef CGAL::Monge_rep My_Monge_rep; + +typedef double LFT; +typedef CGAL::Cartesian Local_Kernel; +typedef CGAL::Monge_info My_Monge_info; +typedef CGAL::Monge_via_jet_fitting 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 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"; +} diff --git a/Jet_fitting_3/test/Jet_fitting_3/data/in_points_d4.txt b/Jet_fitting_3/test/Jet_fitting_3/data/in_points_d4.txt new file mode 100644 index 00000000000..e294b4abf51 --- /dev/null +++ b/Jet_fitting_3/test/Jet_fitting_3/data/in_points_d4.txt @@ -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 \ No newline at end of file diff --git a/Jet_fitting_3/test/Jet_fitting_3/include/LinAlg_lapack.h b/Jet_fitting_3/test/Jet_fitting_3/include/LinAlg_lapack.h new file mode 100644 index 00000000000..1f9d6cc9dc3 --- /dev/null +++ b/Jet_fitting_3/test/Jet_fitting_3/include/LinAlg_lapack.h @@ -0,0 +1,110 @@ +#ifndef _LAPACK_H_ +#define _LAPACK_H_ + +#include +#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 diff --git a/Jet_fitting_3/test/Jet_fitting_3/makefile b/Jet_fitting_3/test/Jet_fitting_3/makefile new file mode 100644 index 00000000000..de4bc9e4e7e --- /dev/null +++ b/Jet_fitting_3/test/Jet_fitting_3/makefile @@ -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 /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