DGtal  1.5.beta
IIGeometricFunctors.h
1 
17 #pragma once
18 
31 #if defined(IIGeometricFunctors_RECURSES)
32 #error Recursive header files inclusion detected in IIGeometricFunctors.h
33 #else // defined(IIGeometricFunctors_RECURSES)
35 #define IIGeometricFunctors_RECURSES
36 
37 #if !defined IIGeometricFunctors_h
39 #define IIGeometricFunctors_h
40 
42 // Inclusions
43 #include <iostream>
44 #include <tuple>
45 #include "DGtal/base/Common.h"
46 #include "DGtal/math/linalg/EigenDecomposition.h"
47 #include "DGtal/math/linalg/CMatrix.h"
49 
50 // @since 0.8 In DGtal::functors
51 namespace DGtal {
52  namespace functors {
53 
55  // template class IINormalDirectionFunctor
66  template <typename TSpace, typename TMatrix=SimpleMatrix< typename TSpace::RealVector::Component, TSpace::dimension, TSpace::dimension> >
68  {
69  // ----------------------- Standard services ------------------------------
70  public:
72  typedef TSpace Space;
73  typedef typename Space::RealVector RealVector;
74  typedef typename RealVector::Component Component;
75  typedef TMatrix Matrix;
76  typedef Matrix Argument;
78  typedef Quantity Value;
79 
82 
86  IINormalDirectionFunctor( const Self& /* other */ ) {}
89  Self& operator=( const Self& /* other */ ) { return *this; }
98  Value operator()( const Argument& arg ) const
99  {
102 
103  ASSERT ( !std::isnan(eigenValues[0]) ); // NaN
104 #ifdef DEBUG
105  for( Dimension i_dim = 1; i_dim < Space::dimension; ++i_dim )
106  {
107  ASSERT ( std::abs(eigenValues[i_dim - 1]) <= std::abs(eigenValues[i_dim]) );
108  }
109 #endif
110 
111  return eigenVectors.column( 0 ); // normal vector is associated to smallest eigenvalue.
112  }
113 
118  void init( Component /* h */, Component /* r */ ) {}
119 
120  private:
125  }; // end of class IINormalDirectionFunctor
126 
127 
129  // template class IITangentDirectionFunctor
140  template <typename TSpace, typename TMatrix=SimpleMatrix< typename TSpace::RealVector::Component, TSpace::dimension, TSpace::dimension> >
142  {
143  // ----------------------- Standard services ------------------------------
144  public:
146  typedef TSpace Space;
147  typedef typename Space::RealVector RealVector;
149  typedef TMatrix Matrix;
150  typedef Matrix Argument;
152  typedef Quantity Value;
153 
156  BOOST_STATIC_ASSERT(( Space::dimension == 2 ));
157 
161  IITangentDirectionFunctor( const Self& /* other */ ) {}
164  Self& operator=( const Self& /* other */ ) { return *this; }
173  Value operator()( const Argument& arg ) const
174  {
177 
178  ASSERT ( !std::isnan(eigenValues[0]) ); // NaN
179 #ifdef DEBUG
180  for( Dimension i_dim = 1; i_dim < Space::dimension; ++i_dim )
181  {
182  ASSERT ( std::abs(eigenValues[i_dim - 1]) <= std::abs(eigenValues[i_dim]) );
183  }
184 #endif
185  return eigenVectors.column( 1 ); // tangent vector is associated to greatest eigenvalue.
186  }
187 
192  void init( Component /* h */, Component /* r */ ) {}
193 
194  private:
199  }; // end of class IITangentDirectionFunctor
200 
201 
203  // template class IIFirstPrincipalDirectionFunctor
218  template <typename TSpace, typename TMatrix=SimpleMatrix< typename TSpace::RealVector::Component, TSpace::dimension, TSpace::dimension> >
220  {
221  // ----------------------- Standard services ------------------------------
222  public:
224  typedef TSpace Space;
225  typedef typename Space::RealVector RealVector;
227  typedef TMatrix Matrix;
228  typedef Matrix Argument;
230  typedef Quantity Value;
231 
234  BOOST_STATIC_ASSERT(( Space::dimension >= 2 ));
235 
239  IIFirstPrincipalDirectionFunctor( const Self& /* other */ ) {}
242  Self& operator=( const Self& /* other */ ) { return *this; }
251  Value operator()( const Argument& arg ) const
252  {
255 
256  ASSERT ( !std::isnan(eigenValues[0]) ); // NaN
257 #ifdef DEBUG
258  for( Dimension i_dim = 1; i_dim < Space::dimension; ++i_dim )
259  {
260  ASSERT ( std::abs(eigenValues[i_dim - 1]) <= std::abs(eigenValues[i_dim]) );
261  }
262 #endif
263 
264  return eigenVectors.column( Space::dimension - 1 ); // first principal curvature direction is associated to greatest eigenvalue.
265  }
266 
271  void init( Component /* h */, Component /* r */ ) {}
272 
273  private:
278  }; // end of class IIFirstPrincipalDirectionFunctor
279 
280 
282  // template class IISecondPrincipalDirectionFunctor
297  template <typename TSpace, typename TMatrix=SimpleMatrix< typename TSpace::RealVector::Component, TSpace::dimension, TSpace::dimension> >
299  {
300  // ----------------------- Standard services ------------------------------
301  public:
303  typedef TSpace Space;
304  typedef typename Space::RealVector RealVector;
306  typedef TMatrix Matrix;
307  typedef Matrix Argument;
309  typedef Quantity Value;
310 
313  BOOST_STATIC_ASSERT(( Space::dimension >= 3 ));
314 
318  IISecondPrincipalDirectionFunctor( const Self& /* other */ ) {}
321  Self& operator=( const Self& /* other */ ) { return *this; }
330  Value operator()( const Argument& arg ) const
331  {
334 
335  ASSERT ( !std::isnan(eigenValues[0]) ); // NaN
336 #ifdef DEBUG
337  for( Dimension i_dim = 1; i_dim < Space::dimension; ++i_dim )
338  {
339  ASSERT ( std::abs(eigenValues[i_dim - 1]) <= std::abs(eigenValues[i_dim]) );
340  }
341 #endif
342 
343  return eigenVectors.column( Space::dimension - 2 ); // second principal curvature direction is associated to greatest eigenvalue.
344  }
345 
350  void init( Component /* h */, Component /* r */ ) {}
351 
352  private:
357  }; // end of class IISecondPrincipalDirectionFunctor
358 
360  // template class IIPrincipalDirectionsFunctor
375  template <typename TSpace, typename TMatrix=SimpleMatrix< typename TSpace::RealVector::Component, TSpace::dimension, TSpace::dimension> >
377  {
378  // ----------------------- Standard services ------------------------------
379  public:
381  typedef TSpace Space;
382  typedef typename Space::RealVector RealVector;
384  typedef TMatrix Matrix;
385  typedef Matrix Argument;
386  typedef std::pair<RealVector,RealVector> Quantity;
387  typedef Quantity Value;
388 
391  BOOST_STATIC_ASSERT(( Space::dimension >= 3 ));
392 
396  IIPrincipalDirectionsFunctor( const Self& /* other */ ) {}
399  Self& operator=( const Self& /* other */ ) { return *this; }
408  Value operator()( const Argument& arg ) const
409  {
412 
413  ASSERT ( !std::isnan(eigenValues[0]) ); // NaN
414 #ifdef DEBUG
415  for( Dimension i_dim = 1; i_dim < Space::dimension; ++i_dim )
416  {
417  ASSERT ( std::abs(eigenValues[i_dim - 1]) <= std::abs(eigenValues[i_dim]) );
418  }
419 #endif
420 
421  return Value(
422  eigenVectors.column( Space::dimension - 1 ),
423  eigenVectors.column( Space::dimension - 2 )
424  );
425  }
426 
431  void init( Component /* h */, Component /* r */ ) {}
432 
433  private:
438  }; // end of class IIPrincipalDirectionsFunctor
439 
440 
442  // template class IIPrincipalCurvaturesAndDirectionsFunctor
457  template <typename TSpace,
458  typename TMatrix=SimpleMatrix< typename TSpace::RealVector::Component, TSpace::dimension,
459  TSpace::dimension> >
461  {
462  // ----------------------- Standard services ------------------------------
463  public:
465  typedef TSpace Space;
466  typedef typename Space::RealVector RealVector;
468  typedef TMatrix Matrix;
469  typedef Matrix Argument;
470  typedef std::tuple<double, double, RealVector, RealVector> Quantity;
471  typedef Quantity Value;
472 
475  BOOST_STATIC_ASSERT(( Space::dimension == 3 ));
476 
484  Value operator()( const Argument& arg ) const
485  {
486  Argument cp_arg = arg;
487  cp_arg *= dh5;
490 
491  ASSERT ( !std::isnan(eigenValues[0]) ); // NaN
492  ASSERT ( (std::abs(eigenValues[0]) <= std::abs(eigenValues[1]))
493  && (std::abs(eigenValues[1]) <= std::abs(eigenValues[2])) );
494 
495  Quantity res(d6_PIr6 * ( eigenValues[2] - ( 3.0 * eigenValues[1] )) + d8_5r,
496  d6_PIr6 * ( eigenValues[1] - ( 3.0 * eigenValues[2] )) + d8_5r,
497  eigenVectors.column( 1 ),
498  eigenVectors.column( 2 ));
499  return res;
500  }
501 
506  void init( Component h , Component r )
507  {
508  double r3 = r * r * r;
509  double r6 = r3 * r3;
510  d6_PIr6 = 6.0 / ( M_PI * r6 );
511  d8_5r = 8.0 / ( 5.0 * r );
512  double h2 = h * h;
513  dh5 = h2 * h2 * h;
514  }
515 
516  private:
517  double dh5;
518  double d6_PIr6;
519  double d8_5r;
524  }; // end of class IIPrincipalCurvaturesAndDirectionsFunctor
525 
527  // template class IICurvatureFunctor
538  template <typename TSpace>
540  {
541  // ----------------------- Standard services ------------------------------
542  public:
544  typedef TSpace Space;
545  typedef typename Space::RealVector RealVector;
549  typedef Quantity Value;
550 
552  BOOST_STATIC_ASSERT(( Space::dimension == 2 ));
553 
560  Value operator()( const Argument& arg ) const
561  {
562  Quantity cp_quantity = arg;
563  cp_quantity *= dh2;
564  return d3_r * ( dPI_2 - d1_r2 * cp_quantity );
565  }
566 
573  void init( Component h, Component r )
574  {
575  d1_r2 = 1.0 / ( r * r );
576  dPI_2 = M_PI / 2.0;
577  d3_r = 3.0 / r;
578  dh2 = h * h;
579  }
580 
581  private:
586  }; // end of class IICurvatureFunctor
587 
588 
590  // template class IIMeanCurvature3DFunctor
601  template <typename TSpace>
603  {
604  // ----------------------- Standard services ------------------------------
605  public:
607  typedef TSpace Space;
608  typedef typename Space::RealVector RealVector;
612  typedef Quantity Value;
613 
615  BOOST_STATIC_ASSERT(( Space::dimension == 3 ));
616 
623  Value operator()( const Argument& arg ) const
624  {
625  Quantity cp_quantity = arg;
626  cp_quantity *= dh3;
627  return d8_3r - d_4_PIr4 * cp_quantity;
628  }
629 
636  void init( Component h, Component r )
637  {
638  d8_3r = 8.0 / ( 3.0 * r );
639  double r2 = r * r;
640  d_4_PIr4 = 4.0 / ( M_PI * r2 * r2 );
641  dh3 = h * h * h;
642  }
643 
644  private:
648  }; // end of class IIMeanCurvature3DFunctor
649 
651  // template class IIGaussianCurvature3DFunctor
665  template <typename TSpace, typename TMatrix=SimpleMatrix< typename TSpace::RealVector::Component, TSpace::dimension, TSpace::dimension> >
667  {
668  // ----------------------- Standard services ------------------------------
669  public:
671  typedef TSpace Space;
672  typedef typename Space::RealVector RealVector;
674  typedef TMatrix Matrix;
675  typedef Matrix Argument;
677  typedef Quantity Value;
678 
681  BOOST_STATIC_ASSERT(( Space::dimension == 3 ));
682 
690  Value operator()( const Argument& arg ) const
691  {
692  Argument cp_arg = arg;
693  cp_arg *= dh5;
696 
697  ASSERT ( !std::isnan(eigenValues[0]) ); // NaN
698  ASSERT ( (std::abs(eigenValues[0]) <= std::abs(eigenValues[1]))
699  && (std::abs(eigenValues[1]) <= std::abs(eigenValues[2])) );
700 
701  Value k1 = d6_PIr6 * ( eigenValues[2] - ( 3.0 * eigenValues[1] )) + d8_5r;
702  Value k2 = d6_PIr6 * ( eigenValues[1] - ( 3.0 * eigenValues[2] )) + d8_5r;
703  return k1 * k2;
704  }
705 
712  void init( Component h, Component r )
713  {
714  double r3 = r * r * r;
715  double r6 = r3 * r3;
716  d6_PIr6 = 6.0 / ( M_PI * r6 );
717  d8_5r = 8.0 / ( 5.0 * r );
718  double h2 = h * h;
719  dh5 = h2 * h2 * h;
720  }
721 
722  private:
726 
731  }; // end of class IIGaussianCurvature3DFunctor
732 
734  // template class IIFirstPrincipalCurvature3DFunctor
748  template <typename TSpace, typename TMatrix=SimpleMatrix< typename TSpace::RealVector::Component, TSpace::dimension, TSpace::dimension> >
750  {
751  // ----------------------- Standard services ------------------------------
752  public:
754  typedef TSpace Space;
755  typedef typename Space::RealVector RealVector;
757  typedef TMatrix Matrix;
758  typedef Matrix Argument;
760  typedef Quantity Value;
761 
764  BOOST_STATIC_ASSERT(( Space::dimension == 3 ));
765 
773  Value operator()( const Argument& arg ) const
774  {
775  Argument cp_arg = arg;
776  cp_arg *= dh5;
779 
780  ASSERT ( !std::isnan(eigenValues[0]) ); // NaN
781  ASSERT ( (std::abs(eigenValues[0]) <= std::abs(eigenValues[1]))
782  && (std::abs(eigenValues[1]) <= std::abs(eigenValues[2])) );
783 
784 
785  return d6_PIr6 * ( eigenValues[2] - ( 3.0 * eigenValues[1] )) + d8_5r;
786  }
787 
794  void init( Component h, Component r )
795  {
796  double r3 = r * r * r;
797  double r6 = r3 * r3;
798  d6_PIr6 = 6.0 / ( M_PI * r6 );
799  d8_5r = 8.0 / ( 5.0 * r );
800  double h2 = h * h;
801  dh5 = h2 * h2 * h;
802  }
803 
804  private:
808 
813  }; // end of class IIFirstPrincipalCurvature3DFunctor
814 
816  // template class IISecondPrincipalCurvature3DFunctor
830  template <typename TSpace, typename TMatrix=SimpleMatrix< typename TSpace::RealVector::Component, TSpace::dimension, TSpace::dimension> >
832  {
833  // ----------------------- Standard services ------------------------------
834  public:
836  typedef TSpace Space;
837  typedef typename Space::RealVector RealVector;
839  typedef TMatrix Matrix;
840  typedef Matrix Argument;
842  typedef Quantity Value;
843 
846  BOOST_STATIC_ASSERT(( Space::dimension == 3 ));
847 
855  Value operator()( const Argument& arg ) const
856  {
857  Argument cp_arg = arg;
858  cp_arg *= dh5;
861 
862  ASSERT ( !std::isnan(eigenValues[0]) ); // NaN
863  ASSERT ( (std::abs(eigenValues[0]) <= std::abs(eigenValues[1]))
864  && (std::abs(eigenValues[1]) <= std::abs(eigenValues[2])) );
865 
866  return d6_PIr6 * ( eigenValues[1] - ( 3.0 * eigenValues[2] )) + d8_5r;
867  }
868 
875  void init( Component h, Component r )
876  {
877  double r3 = r * r * r;
878  double r6 = r3 * r3;
879  d6_PIr6 = 6.0 / ( M_PI * r6 );
880  d8_5r = 8.0 / ( 5.0 * r );
881  double h2 = h * h;
882  dh5 = h2 * h2 * h;
883  }
884 
885  private:
889 
894  }; // end of class IISecondPrincipalCurvature3DFunctor
895 
896 
898  // template class IIPrincipalCurvatures3DFunctor
911  template <typename TSpace, typename TMatrix=SimpleMatrix< typename TSpace::RealVector::Component, TSpace::dimension, TSpace::dimension> >
913  {
914  // ----------------------- Standard services ------------------------------
915  public:
917  typedef TSpace Space;
918  typedef typename Space::RealVector RealVector;
920  typedef TMatrix Matrix;
921  typedef Matrix Argument;
922  typedef std::pair<Component, Component> Quantity;
923  typedef Quantity Value;
924 
927  BOOST_STATIC_ASSERT(( Space::dimension == 3 ));
928 
937  Value operator()( const Argument& arg ) const
938  {
939  Argument cp_arg = arg;
940  cp_arg *= dh5;
943 
944  ASSERT ( !std::isnan(eigenValues[0]) ); // NaN
945  ASSERT ( (std::abs(eigenValues[0]) <= std::abs(eigenValues[1]))
946  && (std::abs(eigenValues[1]) <= std::abs(eigenValues[2])) );
947 
948  return Value(
949  d6_PIr6 * ( eigenValues[2] - ( 3.0 * eigenValues[1] )) + d8_5r,
950  d6_PIr6 * ( eigenValues[1] - ( 3.0 * eigenValues[2] )) + d8_5r
951  );
952  }
953 
960  void init( Component h, Component r )
961  {
962  double r3 = r * r * r;
963  double r6 = r3 * r3;
964  d6_PIr6 = 6.0 / ( M_PI * r6 );
965  d8_5r = 8.0 / ( 5.0 * r );
966  double h2 = h * h;
967  dh5 = h2 * h2 * h;
968  }
969 
970  private:
971  double dh5;
972  double d6_PIr6;
973  double d8_5r;
974 
979  }; // end of class IIPrincipalCurvatures3DFunctor
980 
981 } // namespace functors
982 
983 } // namespace DGtal
984 
985 
986 // //
988 
989 #endif // !defined IIGeometricFunctors_h
990 
991 #undef IIGeometricFunctors_RECURSES
992 #endif // else defined(IIGeometricFunctors_RECURSES)
static void getEigenDecomposition(const Matrix &matrix, Matrix &eigenVectors, Vector &eigenValues)
Compute both eigen vectors and eigen values from an input matrix.
Aim: Implements basic operations that will be used in Point and Vector classes.
Definition: PointVector.h:593
TEuclideanRing Component
Type for Vector elements.
Definition: PointVector.h:614
Aim: implements basic MxN Matrix services (M,N>=1).
Definition: SimpleMatrix.h:76
Aim: A functor Matrix -> RealVector that returns the first principal curvature direction by diagonali...
BOOST_CONCEPT_ASSERT((concepts::CMatrix< Matrix >))
BOOST_CONCEPT_ASSERT((concepts::CSpace< TSpace >))
RealVector eigenValues
A data member only used for temporary calculations.
IIFirstPrincipalDirectionFunctor(const Self &)
Copy constructor. Nothing to do.
Matrix eigenVectors
A data member only used for temporary calculations.
IIFirstPrincipalDirectionFunctor< TSpace > Self
Aim: A functor Matrix -> RealVector that returns the normal direction by diagonalizing the given cova...
IINormalDirectionFunctor(const Self &)
Copy constructor. Nothing to do.
BOOST_CONCEPT_ASSERT((concepts::CMatrix< Matrix >))
Value operator()(const Argument &arg) const
RealVector eigenValues
A data member only used for temporary calculations.
BOOST_CONCEPT_ASSERT((concepts::CSpace< TSpace >))
Matrix eigenVectors
A data member only used for temporary calculations.
IINormalDirectionFunctor< TSpace > Self
Aim: A functor Matrix -> std::pair<RealVector,RealVector> that returns the first and the second princ...
Matrix eigenVectors
A data member only used for temporary calculations.
IIPrincipalCurvaturesAndDirectionsFunctor< TSpace > Self
std::tuple< double, double, RealVector, RealVector > Quantity
RealVector eigenValues
A data member only used for temporary calculations.
Aim: A functor Matrix -> std::pair<RealVector,RealVector> that returns the first and the second princ...
Value operator()(const Argument &arg) const
RealVector eigenValues
A data member only used for temporary calculations.
std::pair< RealVector, RealVector > Quantity
IIPrincipalDirectionsFunctor(const Self &)
Copy constructor. Nothing to do.
BOOST_CONCEPT_ASSERT((concepts::CMatrix< Matrix >))
IIPrincipalDirectionsFunctor< TSpace > Self
BOOST_CONCEPT_ASSERT((concepts::CSpace< TSpace >))
Matrix eigenVectors
A data member only used for temporary calculations.
Aim: A functor Matrix -> RealVector that returns the second principal curvature direction by diagonal...
Matrix eigenVectors
A data member only used for temporary calculations.
RealVector eigenValues
A data member only used for temporary calculations.
IISecondPrincipalDirectionFunctor(const Self &)
Copy constructor. Nothing to do.
BOOST_CONCEPT_ASSERT((concepts::CSpace< TSpace >))
IISecondPrincipalDirectionFunctor< TSpace > Self
BOOST_CONCEPT_ASSERT((concepts::CMatrix< Matrix >))
Aim: A functor Matrix -> RealVector that returns the tangent direction by diagonalizing the given cov...
Value operator()(const Argument &arg) const
IITangentDirectionFunctor< TSpace > Self
BOOST_CONCEPT_ASSERT((concepts::CSpace< TSpace >))
Matrix eigenVectors
A data member only used for temporary calculations.
RealVector eigenValues
A data member only used for temporary calculations.
IITangentDirectionFunctor(const Self &)
Copy constructor. Nothing to do.
BOOST_CONCEPT_ASSERT((concepts::CMatrix< Matrix >))
BOOST_STATIC_ASSERT((Space::dimension==2))
DGtal is the top-level namespace which contains all DGtal functions and types.
DGtal::uint32_t Dimension
Definition: Common.h:136
Aim: Represent any static or dynamic sized matrix having sparse or dense representation.
Definition: CMatrix.h:91
Aim: Defines the concept describing a digital space, ie a cartesian product of integer lines.
Definition: CSpace.h:106
Aim: A functor Real -> Real that returns the 2d curvature by transforming the given volume....
IICurvatureFunctor< TSpace > Self
BOOST_STATIC_ASSERT((Space::dimension==2))
void init(Component h, Component r)
Value operator()(const Argument &arg) const
BOOST_CONCEPT_ASSERT((concepts::CSpace< TSpace >))
Aim: A functor Matrix -> Real that returns the first principal curvature value by diagonalizing the g...
BOOST_CONCEPT_ASSERT((concepts::CMatrix< Matrix >))
IIFirstPrincipalCurvature3DFunctor< TSpace > Self
RealVector eigenValues
A data member only used for temporary calculations.
Matrix eigenVectors
A data member only used for temporary calculations.
BOOST_CONCEPT_ASSERT((concepts::CSpace< TSpace >))
Aim: A functor Matrix -> Real that returns the Gaussian curvature by diagonalizing the given covarian...
IIGaussianCurvature3DFunctor< TSpace > Self
BOOST_CONCEPT_ASSERT((concepts::CSpace< TSpace >))
BOOST_CONCEPT_ASSERT((concepts::CMatrix< Matrix >))
Value operator()(const Argument &arg) const
Matrix eigenVectors
A data member only used for temporary calculations.
RealVector eigenValues
A data member only used for temporary calculations.
Aim: A functor Real -> Real that returns the 3d mean curvature by transforming the given volume....
IIMeanCurvature3DFunctor< TSpace > Self
BOOST_STATIC_ASSERT((Space::dimension==3))
Value operator()(const Argument &arg) const
BOOST_CONCEPT_ASSERT((concepts::CSpace< TSpace >))
Aim: A functor Matrix -> std::pair<Real,Real> that returns the first and the second principal curvatu...
BOOST_CONCEPT_ASSERT((concepts::CMatrix< Matrix >))
Matrix eigenVectors
A data member only used for temporary calculations.
IIPrincipalCurvatures3DFunctor< TSpace > Self
RealVector eigenValues
A data member only used for temporary calculations.
BOOST_CONCEPT_ASSERT((concepts::CSpace< TSpace >))
Aim: A functor Matrix -> Real that returns the second principal curvature value by diagonalizing the ...
RealVector eigenValues
A data member only used for temporary calculations.
IISecondPrincipalCurvature3DFunctor< TSpace > Self
BOOST_CONCEPT_ASSERT((concepts::CSpace< TSpace >))
Matrix eigenVectors
A data member only used for temporary calculations.
BOOST_CONCEPT_ASSERT((concepts::CMatrix< Matrix >))