This package allows the estimation of local differential quantities of a surface from a point sample, given either as a mesh or as point cloud.
Note that this package needs the third party libraries Lapack and Blas to be installed to compile the example code.
Consider a sampled smooth surface, and assume we are given a
collection of points P about a given sample p. We aim at
estimating the differential properties up to any fixed order of the
surface at point p from the point set P+ = P { p} - we
denote N=
P+
. More precisely, first order properties
correspond to the normal or the tangent plane; second order properties
provide the principal curvatures and directions, third order
properties provide the directional derivatives of the principal
curvatures along the curvature lines, etc. Most of the time,
estimating first and second order differential quantities is
sufficient. However, some applications involving shape analysis
require estimating third and fourth order differential quantities.
Many different estimators have been proposed in the vast literature of
applied geometry [Pet01] (section 3, page 7), and all
of them need to define a neighborhood around the point at which the
estimation is computed. Our method relies on smooth differential
geometry calculations, carried out on smooth objects fitted from
the sample points.
Datasets amenable to such a processing are naturally unstructured
point clouds, as well as meshes - whose topological information may
be discarded.
Estimating differential properties from discrete date always raises a philosophical issue. On one hand, estimating differential quantities subsumes a smooth surface does exist. In this spirit one wishes to recover its differential properties, so that any estimation method must come with an asymptotic convergence analysis of the results returned. For the method developed in this CGAL package, the interested will find such an analysis in [CP05a], (Theorem 3) - it should be stressed the error bounds proved therein are optimal.
On the other hand, any estimation method may be applied to arbitrarily data - surface unknown, surface piecewise smooth etc. In such a case, no analysis can be carried out, and it is up to the users to check the results match their needs.
Unlike most of the CGAL packages, this package uses approximation methods and is not intended to provide an exact canonical result in any sense. This is why internal computations are performed with a number type possibly different from that of the input data, even if for convenience the results are returned with this original number type. A reasonable choice for this internal number type is for example the double type.
To present the method, we shall need the following notions. Consider a smooth surface. About one of its points, consider a coordinate system whose z-axis does not belong to the tangent space. In such a frame, the surface can locally be written as the graph of a bivariate function. Letting h.o.t. stand for higher order terms, one has :
z(x,y)=JB,d(x,y) + h.o.t. ; JB,d(x,y)= k=0k=d(
i=0i=k
(Bk-i,ixk-iyi)/(i!(k-i)!)).
The degree d polynomial JB,d is the Taylor expansion of the function z, and is called its d-jet. Notice that a d-jet contains Nd=(d+1)(d+2)/2 coefficients.
Recall that an umbilical point of a surface - or umbilic for short, is a point where both principal curvatures are identical. At any point of the surface which is not an umbilic, principal directions d1, d2 are well defined, and these (non oriented) directions together with the normal vector n define two direct orthonormal frames. If v1 is a unit vector of direction d1, there exists a unique unit vector v2 so that (v1,v2,n) is direct; and the other possible frame is (-v1,-v2,n). Both these coordinate systems are known as the Monge coordinate systems. In both these systems, the surface is said to be given in the Monge form and its jet has the following canonical form :
z(x,y) = | (1)/(2)(k1x2 + k2y2)+ (1)/(6)(b0x3+3b1x2y+3b2xy2+b3y3) |
+(1)/(24)(c0x4+4c1x3y+6c2x2y2+4c3xy3+c4y4) + h.o.t. |
The coefficients k1, k2 are the principal curvatures, b0,b3 are the directional derivatives of k1,k2 along their respective curvature line, while b1,b2 are the directional derivatives of k1,k2 along the other curvature lines.
The Monge coordinate system can be computed from any d-jet (d 2), and so are the Monge coefficients. These informations
characterize the local geometry of the surface in a canonical way, and
are the quantities returned by our algorithm.
Further details can be found in section 42.4 and in [CP05a] (section 6).
As usual, the fitting procedure may run into (almost) degenerate cases:
The fitting strategy performed by the class Monge_via_jet_fitting requires the following parameters:
As explained in Section 42.1, the output consists of a coordinate system, the Monge basis, together with the Monge coefficients which are stored in the Monge_form class. In addition, more information on the computational issues are stored in the Monge_via_jet_fitting class.
The Monge_form class provides the following information.
In addition, the class Monge_via_jet_fitting stores
This concept provides the types for the input sample points, together with 3d vectors and a number type. It is used as template for the classe Monge_via_jet_fitting<DataKernel, LocalKernel = Cartesian<double>, SvdTraits = lapack_svd> . Typically, one can use CGAL::Cartesian<double>.
Input points of type DataKernel::Point_3 are converted to LocalKernel::Point_3. For output of the Monge_form class, these types are converted back to Data_Kernel ones. Typically, one can use CGAL::Cartesian<double> which is the default.
This concept provides the number, vector and matrix types for algebra operations required by the fitting method in Monge_via_jet_fitting<DataKernel, LocalKernel = Cartesian<double>, SvdTraits = lapack_svd> . The main method is a linear solver using a singular value decomposition.
To solve the fitting problem, the sample points are first converted from the DataKernel to the LocalKernel (this is done using the CGAL::Cartesian_converter). Then change of coordinate systems and linear algebra operations are performed with this kernel. This implies that the number types LocalKernel::FT and SvdTraits::FT must be identical. Second the Monge basis and coefficients, computed with the LocalKernel, are converted back to the DataKernel (this is done using the CGAL::Cartesian_converter and the CGAL::NT_converter).
The first example illustrates the computation of the local differential quantities from a set of points given in a text file as input. The first point of the list is the one at which the computation is performed. The user has to specify a file for the input points and the degrees d and d'.
File: examples/Jet_fitting_3/Single_estimation.cpp
#include <CGAL/Cartesian.h> #ifndef CGAL_USE_LAPACK int main() { std::cerr << "Skip since LAPACK is not installed" << std::endl; std::cerr << std::endl; return 0; } #else #include <fstream> #include <vector> #include <CGAL/Monge_via_jet_fitting.h> typedef double DFT; typedef CGAL::Cartesian<DFT> Data_Kernel; typedef Data_Kernel::Point_3 DPoint; typedef CGAL::Monge_via_jet_fitting<Data_Kernel> My_Monge_via_jet_fitting; typedef My_Monge_via_jet_fitting::Monge_form My_Monge_form; int main(int argc, char *argv[]) { //check command line if (argc<4) { std::cout << " Usage : " << argv[0] << " <inputPoints.txt> <d_fitting> <d_monge>" << std::endl; return 0; } //open the input file std::ifstream inFile( argv[1]); 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; while (inFile >> x) { inFile >> y >> z; DPoint p(x,y,z); in_points.push_back(p); } inFile.close(); // fct parameters size_t d_fitting = std::atoi(argv[2]); size_t d_monge = std::atoi(argv[3]); My_Monge_form monge_form; My_Monge_via_jet_fitting monge_fit; monge_form = monge_fit(in_points.begin(), in_points.end(), d_fitting, d_monge); //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 << monge_form; std::cout << "condition_number : " << monge_fit.condition_number() << std::endl << "pca_eigen_vals and associated pca_eigen_vecs :" << std::endl; for (int i=0; i<3; i++) std::cout << monge_fit.pca_basis(i).first << std::endl << monge_fit.pca_basis(i).second << std::endl; return 0; } #endif //CGAL_USE_LAPACK
Figs. 42.1 and 42.2 provide illustrations of principal directions of curvature.
Figure 42.2: Principal directions of curvature and normals at vertices of a mesh of the graph of the function f(x,y)=2x2+y2.
![]() | advanced |
![]() |
Figure 42.3: The three bases involved in the estimation.
Input : samples
Output : fitting-basis
Performing a PCA requires diagonalizing a symmetric matrix. This analysis gives an orthonormal basis whose z-axis is provided by the eigenvector associated to the smallest eigenvalue.1 Note one may have to swap the orientation of a vector to get a direct basis.
Let us denote PW F the matrix that changes coordinates from the
world-basis (wx,wy,wz) to the fitting-basis (fx,fy,fz). The
rows of PW
F are the coordinates of the vectors
(fx,fy,fz) in the world-basis. This matrix represents a
orthogonal transformation hence its inverse is its transpose. To obtain
the coordinates of a point in the fitting-basis from the coordinates
in the world-basis, one has to multiply by PW
F.
As mentioned above, the eigenvalues are returned, from which the sampling quality can be assessed. For a good sampling, the eigenvector associated to the smallest eigenvalue should roughly give the normal direction.
Input : samples, fitting-basis
Output : coefficients Ai,j
of the bivariate fitted polynomial in the fitting-basis
Computations are done in the fitting-basis and the origin is the point
p. First, one has to transform coordinates of sample points with a
translation (-p) and multiplication by PW F.
The fitting process consists of finding the coefficients Ai,j of the degree d polynomial
JA,d= k=0k=d(
i=0i=k
(Ak-i,ixk-iyi)/(i!(k-i)!)).
Denote pi=(xi,yi,zi), i=1,..., N the coordinates of the
sample points of P+.
For interpolation the linear equations to solve are A(xi,yi)=zi i=1,...,N, and for approximation one has to minimize i=1N
(A(xi,yi)-zi)2. The linear algebra formulation of the problem is
given by
A = | (A0,0, A1,0,A0,1, ..., A0,d)T |
Z= | (z1, z2,..., zN)T |
M= | (1,xi, yi, (xi2)/(2),..., (xiyid-1)/((d-1)!), (yid)/(d!))i=1,...,N |
The equations for interpolation become MA=Z. For approximation, the system MA=Z is solved in the least square sense, i.e. one seeks the vector A such that A = argminA ||MA-Z||2.
In any case, there is a preconditioning of the matrix M so as to improve the condition number. Assuming the {xi}, {yi} are of order h, the pre-conditioning consists of performing a column scaling by dividing each monomial xikyil by hk+l - refer to Eq. (42.4.2). Practically, the parameter h is chosen as the mean value of the {xi} and {yi}. In other words, the new system is M'Y=(MD-1)(DA)=Z with D the diagonal matrix D=(1,h,h,h2,...,hd,hd), so that the solution A of the original system is A=D-1Y.
There is always a single solution since for under constrained systems we also minimize ||A||2. The method uses a singular value decomposition of the N × Nd matrix M= U S VT, where U is a N × N orthogonal matrix, V is a Nd × Nd orthogonal matrix and S is a N × Nd matrix with the singular values on its diagonal. Denote r the rank of M, we can decompose S= (
Dr | 0r, Nd-r |
0N-r, r | 0N-r, Nd-r |
A= V (
Dr-1 | 0Nd-r, r |
0r, N-r | 0Nd-r, N-r |
One can provide the condition number of the matrix M (after preconditioning) which is the ratio of the maximal and the minimal singular values. It is infinite if the system is under constrained, that is the smallest singular value is zero.
Implementation details. We assume a solve function is provided by the traits SvdTraits. This function 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.
Remark: as an alternative, other methods may be used to solve the
system. A QR decomposition can be substituted to the SVD. One can
also use the normal equation MTMX=MTB and apply methods for square
systems such as LU, QR or Cholesky since MTM is symmetric
definite positive when M has full rank.
The advantages of the SVD
is that it works directly on the rectangular system and gives the
condition number of the system. For more on these alternatives, see
[GvL83] (Chap. 5).
Input : coefficients of the fit Ai,j,
fitting-basis
Output : Monge basis wrt fitting-basis and world-basis
In the fitting basis, we have determined a height function expressed
by Eq. (42.4.2). Computations are done in the fitting-basis.
The partial derivatives, evaluated at (x,y)=(0,0), of the fitted
polynomial JA,d(x,y) are
Ai,j=( i+jJA,d)/(
ix
jy)
Expanding Eq. (42.4.2) yields:
JA,d(x,y) | = A0,0+A1,0x+A0,1y+(1)/(2)(A2,0x2+2A1,1xy+A0,2y2) + (1)/(6)(A3,0x3+3A2,1x2y+...)+ ... |
Input : coefficients of the fit, Monge basis wrt fitting-basis (PF
M)
Output : third and fourth order coefficients of Monge
We use explicit formula. The implicit equation of the fitted polynomial surface in the fitting-basis with origin the point (0,0,A0,0) is Q=0 with
Q=-w-A0,0 + i,j(Ai,juivj)/(i!j!).
The equation in the Monge basis is obtained by substituting (u,v,w)
by PTF M(x,y,z). Denote f(x,y,z)=0 this implicit
equation. By definition of the Monge basis, we have locally (at
(0,0,0))
f(x,y,z)=0 z=g(x,y)
and the Taylor expansion of g at (0,0) are the Monge coefficients
sought.
Let us denote the partial derivatives evaluated at the origin of f
and g by fi,j,k=( i+j+kf)/(
ix
jy
kz) and gi,j=(
i+jg)/(
ix
jy). One has f1,0,0=f0,1,0=f1,1,0=0,
g0,0=g1,0=g0,1=g1,1=0 and g2,0=k1,
g0,2=k2. The partial derivative of order n of f depends on
the matrix PF
M and the partial derivatives of order
at most n of JA,d. The third and fourth order coefficients of are
computed with the implicit function theorem. For instance :
b0=g3,0=-( f3,0,0 f0,0,1 -3 f1,0,1 f2,0,0 )/( f0,0,1 2) | |
b1=g2,1=-(- f0,1,1 f2,0,0 + f2,1,0 f0,0,1 )/( f0,0,1 2) | |
.... |
![]() | advanced |
![]() |
1 | Another possibility is to choose as z-axis the axis of the world-basis with the least angle with the axis determined with the PCA. Then the change of basis reduces to a permutation of axis. |