DGtal  0.9.3
testVoronoiCovarianceMeasureOnSurface.cpp
1 
30 #include <iostream>
32 #include <vector>
33 #include "DGtal/base/Common.h"
34 #include "DGtal/helpers/StdDefs.h"
35 #include "DGtal/math/Statistic.h"
36 #include "DGtal/math/MPolynomial.h"
37 #include "DGtal/topology/LightImplicitDigitalSurface.h"
38 #include "DGtal/geometry/surfaces/estimation/CSurfelLocalEstimator.h"
39 #include "DGtal/geometry/surfaces/estimation/CDigitalSurfaceLocalEstimator.h"
40 #include "DGtal/geometry/surfaces/estimation/VoronoiCovarianceMeasureOnDigitalSurface.h"
41 #include "DGtal/geometry/surfaces/estimation/VCMDigitalSurfaceLocalEstimator.h"
42 #include "DGtal/geometry/surfaces/estimation/TrueDigitalSurfaceLocalEstimator.h"
43 #include "DGtal/geometry/surfaces/estimation/IIGeometricFunctors.h"
44 #include "DGtal/geometry/surfaces/estimation/IntegralInvariantCovarianceEstimator.h"
45 #include "DGtal/geometry/surfaces/estimation/IntegralInvariantVolumeEstimator.h"
46 #include "DGtal/shapes/GaussDigitizer.h"
47 #include "DGtal/shapes/ShapeGeometricFunctors.h"
48 #include "DGtal/shapes/implicit/ImplicitPolynomial3Shape.h"
49 #include "DGtal/io/readers/MPolynomialReader.h"
50 
52 
54 // Functions for testing class VoronoiCovarianceMeasureOnSurface
56 
57 
62 bool testVoronoiCovarianceMeasureOnSurface()
63 {
64  unsigned int nbok = 0;
65  unsigned int nb = 0;
66 
67  using namespace DGtal;
68  using namespace DGtal::Z3i; // gets Space, Point, Domain
69 
73  typedef GaussDigitizer< Z3i::Space, ImplicitShape > ImplicitDigitalShape;
75  typedef DigitalSurface< SurfaceContainer > Surface;
77  typedef SurfaceContainer::Surfel Surfel;
79  // typedef functors::HatPointFunction<Point,double> KernelFunction;
82  trace.beginBlock("Creating Surface");
83  // std::string poly_str = "1.0-0.16*x^2+0.22*y^2+0.3*z^2";
84  std::string poly_str = "-81.0+x^2+y^2+z^2";
85  // NB (JOL): II is sensitive to orientation !!
86  // std::string poly_str = "81.0-x^2-y^2-z^2";
87  Polynomial3 poly;
88  Polynomial3Reader reader;
89  std::string::const_iterator iter = reader.read( poly, poly_str.begin(), poly_str.end() );
90  if ( iter != poly_str.end() )
91  {
92  std::cerr << "ERROR: I read only <"
93  << poly_str.substr( 0, iter - poly_str.begin() )
94  << ">, and I built P=" << poly << std::endl;
95  return 1;
96  }
97  CountedPtr<ImplicitShape> shape( new ImplicitShape( poly ) );
98 
99  Point p1( -10, -10, -10 );
100  Point p2( 10, 10, 10 );
101  KSpace K;
102  nbok += K.init( p1, p2, true ) ? 1 : 0;
103  nb++;
104  trace.info() << "(" << nbok << "/" << nb << ") "
105  << "K.init() is ok" << std::endl;
106 
107  // Digitizer
108  CountedPtr<ImplicitDigitalShape> dshape( new ImplicitDigitalShape() );
109  dshape->attach( *shape );
110  dshape->init( p1, p2, 1.0 );
111  Surfel bel = Surfaces<KSpace>::findABel( K, *dshape, 10000 );
112  SurfaceContainer* surfaceContainer = new SurfaceContainer
113  ( K, *dshape, SurfelAdjacency<KSpace::dimension>( true ), bel );
114  CountedConstPtrOrConstPtr<Surface> ptrSurface( new Surface( surfaceContainer ) ); // acquired
115  trace.endBlock();
116 
117  trace.beginBlock("Computing VCM on surface." );
118  KernelFunction chi( 1.0, 7.0 );
119  CountedPtr<VCMOnSurface> vcm_surface( new VCMOnSurface( ptrSurface, Pointels,
120  5.0, 7.0, chi, 7.0, Metric(), true ) );
121  trace.endBlock();
122 
123  trace.beginBlock("Wrapping normal estimator." );
124  typedef functors::VCMNormalVectorFunctor<VCMOnSurface> NormalVectorFunctor;
125  typedef VCMDigitalSurfaceLocalEstimator<SurfaceContainer,Metric,
126  KernelFunction, NormalVectorFunctor> VCMNormalEstimator;
127  VCMNormalEstimator estimator( vcm_surface );
128  estimator.init( 1.0, ptrSurface->begin(), ptrSurface->end() );
129  trace.endBlock();
130 
131  trace.beginBlock("Computing II normals on surfels of volume." );
132  typedef functors::IINormalDirectionFunctor<Space> IINormalFunctor;
133  typedef IntegralInvariantCovarianceEstimator<KSpace, ImplicitDigitalShape,
134  IINormalFunctor> IINormalEstimator;
135  IINormalEstimator ii_estimator( K, *dshape );
136  ii_estimator.setParams( 5.0 );
137  ii_estimator.init( 1.0, ptrSurface->begin(), ptrSurface->end() );
138  trace.endBlock();
139 
140  trace.beginBlock("Evaluating normals wrt true normal." );
141  typedef ShapeGeometricFunctors::ShapeNormalVectorFunctor<ImplicitShape> NormalFunctor;
143 
147 
148  TrueNormalEstimator true_estimator;
149  true_estimator.setParams( K, NormalFunctor() );
150  true_estimator.attach( shape );
151  true_estimator.init( 1.0, ptrSurface->begin(), ptrSurface->end() );
152  Statistic<double> error_true;
153  Statistic<double> error_triv_true;
154  Statistic<double> error_ii_true;
155  for ( ConstIterator it = ptrSurface->begin(), itE = ptrSurface->end(); it != itE; ++it )
156  {
157  RealVector n_est = estimator.eval( it );
158  RealVector n_true = true_estimator.eval( it );
159  RealVector n_triv = - vcm_surface->mapSurfel2Normals().find( *it )->second.trivialNormal;
160  RealVector n_ii = ii_estimator.eval( it );
161  error_true.addValue( n_est.dot( n_true ) );
162  error_triv_true.addValue( n_triv.dot( n_true ) );
163  if ( n_ii.dot( n_triv ) < 0 ) n_ii = -n_ii;
164  error_ii_true.addValue( n_ii.dot( n_true ) );
165  }
166  error_true.terminate();
167  error_triv_true.terminate();
168  error_ii_true.terminate();
169  trace.info() << "VCM/true cos angle avg = " << error_true.mean() << std::endl;
170  trace.info() << "VCM/true cos angle dev = " << sqrt( error_true.unbiasedVariance() ) << std::endl;
171  trace.info() << "II/true cos angle avg = " << error_ii_true.mean() << std::endl;
172  trace.info() << "II/true cos angle dev = " << sqrt( error_ii_true.unbiasedVariance() ) << std::endl;
173  trace.info() << "triv/true cos angle avg = " << error_triv_true.mean() << std::endl;
174  trace.info() << "triv/true cos angle dev = " << sqrt( error_triv_true.unbiasedVariance() ) << std::endl;
175  nbok += error_true.mean() > 0.95 ? 1 : 0;
176  nb++;
177  trace.info() << "(" << nbok << "/" << nb << ") "
178  << "VCM/true cos angle avg > 0.95" << std::endl;
179  nbok += sqrt( error_true.unbiasedVariance() ) < 0.05 ? 1 : 0;
180  nb++;
181  trace.info() << "(" << nbok << "/" << nb << ") "
182  << "VCM/true cos angle dev < 0.05" << std::endl;
183  // nbok += error_true.mean() > error_triv_true.mean() ? 1 : 0;
184  // nb++;
185  // trace.info() << "(" << nbok << "/" << nb << ") "
186  // << "VCM/true is closer to 1.0 than triv/true." << std::endl;
187 
188  nbok += error_ii_true.mean() > 0.95 ? 1 : 0;
189  nb++;
190  trace.info() << "(" << nbok << "/" << nb << ") "
191  << "II/true cos angle avg > 0.95" << std::endl;
192  nbok += sqrt( error_ii_true.unbiasedVariance() ) < 0.05 ? 1 : 0;
193  nb++;
194  trace.info() << "(" << nbok << "/" << nb << ") "
195  << "II/true cos angle dev < 0.05" << std::endl;
196  nbok += error_ii_true.mean() > error_triv_true.mean() ? 1 : 0;
197  nb++;
198  trace.info() << "(" << nbok << "/" << nb << ") "
199  << "II/true is closer to 1.0 than triv/true." << std::endl;
200 
201  trace.endBlock();
202 
203  trace.beginBlock("Computing II mean curvatures." );
204  typedef functors::IIMeanCurvature3DFunctor<Space> IICurvatureFunctor;
205  typedef IntegralInvariantVolumeEstimator<KSpace, ImplicitDigitalShape,
206  IICurvatureFunctor> IICurvatureEstimator;
207  IICurvatureEstimator ii_curv_estimator( K, *dshape );
208  ii_curv_estimator.setParams( 5.0 );
209  ii_curv_estimator.init( 1.0, ptrSurface->begin(), ptrSurface->end() );
210  trace.endBlock();
211 
212  trace.beginBlock("Computing VCM mean curvatures." );
214  typedef VCMDigitalSurfaceLocalEstimator<SurfaceContainer,Metric,
215  KernelFunction, VCMCurvatureFunctor> VCMCurvatureEstimator;
216  VCMCurvatureEstimator vcm_curv_estimator( vcm_surface );
217  vcm_curv_estimator.init( 1.0, ptrSurface->begin(), ptrSurface->end() );
218  trace.endBlock();
219 
220  trace.beginBlock("Computing ground truth mean curvatures." );
221  typedef ShapeGeometricFunctors::ShapeMeanCurvatureFunctor<ImplicitShape> CurvatureFunctor;
223 
224  TrueCurvatureEstimator true_curv_estimator;
225  true_curv_estimator.setParams( K, CurvatureFunctor() );
226  true_curv_estimator.attach( shape );
227  true_curv_estimator.init( 1.0, ptrSurface->begin(), ptrSurface->end() );
228  trace.endBlock();
229 
233 
234 
235  trace.beginBlock("Evaluating curvatures." );
236  Statistic<double> stat_vcm_curv;
237  Statistic<double> stat_ii_curv;
238  Statistic<double> stat_true_curv;
239  for ( ConstIterator it = ptrSurface->begin(), itE = ptrSurface->end(); it != itE; ++it )
240  {
241  stat_ii_curv.addValue( ii_curv_estimator.eval( it ) );
242  stat_vcm_curv.addValue( vcm_curv_estimator.eval( it ) );
243  stat_true_curv.addValue( true_curv_estimator.eval( it ) );
244  }
245  stat_ii_curv.terminate();
246  stat_vcm_curv.terminate();
247  stat_true_curv.terminate();
248  trace.info() << "- II curv: E[X]=" << stat_ii_curv.mean()
249  << " min=" << stat_ii_curv.min()
250  << " max=" << stat_ii_curv.max() << std::endl;
251  trace.info() << "- VCM curv: E[X]=" << stat_vcm_curv.mean()
252  << " min=" << stat_vcm_curv.min()
253  << " max=" << stat_vcm_curv.max() << std::endl;
254  trace.info() << "- TRUE curv: E[X]=" << stat_true_curv.mean()
255  << " min=" << stat_true_curv.min()
256  << " max=" << stat_true_curv.max() << std::endl;
257  nbok += ( std::abs( stat_ii_curv.mean() - 0.11 ) < 0.05 ) ? 1 : 0;
258  nb++;
259  trace.info() << "(" << nbok << "/" << nb << ") "
260  << "mean of II curv is around 0.11: "
261  << std::abs( stat_ii_curv.mean() - 0.11 ) << std::endl;
262  nbok += ( std::abs( stat_vcm_curv.mean() - 0.11 ) < 0.05 ) ? 1 : 0;
263  nb++;
264  trace.info() << "(" << nbok << "/" << nb << ") "
265  << "mean of VCM curv is around 0.11: "
266  << std::abs( stat_vcm_curv.mean() - 0.11 ) << std::endl;
267 
268  trace.endBlock();
269  return nbok == nb;
270 }
271 
273 // Standard services - public :
274 
275 int main( int /* argc */, char** /* argv */ )
276 {
277  using namespace std;
278  using namespace DGtal;
279  trace.beginBlock ( "Testing VoronoiCovarianceMeasureOnSurface ..." );
280  bool res = testVoronoiCovarianceMeasureOnSurface();
281  trace.emphase() << ( res ? "Passed." : "Error." ) << endl;
282  trace.endBlock();
283  return res ? 0 : 1;
284 }
285 // //
void beginBlock(const std::string &keyword="")
Iterator read(Polynomial &p, Iterator begin, Iterator end)
Aim: A model of CDigitalSurfaceContainer which defines the digital surface as the boundary of an impl...
Z3i this namespace gathers the standard of types for 3D imagery.
Aim: This class adapts a VoronoiCovarianceMeasureOnDigitalSurface to be a model of CDigitalSurfaceLoc...
MyDigitalSurface::ConstIterator ConstIterator
Aim: Smart pointer based on reference counts.
Definition: CountedPtr.h:79
Trace trace
Definition: Common.h:137
Component dot(const Self &v) const
double mean() const
Aim: Smart or simple const pointer on T. It can be a smart pointer based on reference counts or a sim...
Aim: This class converts a string polynomial expression in a multivariate polynomial.
STL namespace.
double endBlock()
Aim: A functor Surfel -> Quantity that returns the outer normal vector at given surfel.
Aim: Represents a set of n-1-cells in a nD space, together with adjacency relation between these cell...
Aim: Represents a multivariate polynomial, i.e. an element of , where K is some ring or field...
Definition: MPolynomial.h:58
Aim: This class processes a set of sample values for one variable and can then compute different stat...
Definition: Statistic.h:69
Aim: This class implement an Integral Invariant estimator which computes for each surfel the covarian...
Aim: This class implement an Integral Invariant estimator which computes for each surfel the volume o...
static SCell findABel(const KSpace &K, const PointPredicate &pp, unsigned int nbtries=1000)
Aim: implements separable l_p metrics with exact predicates.
void setParams(ConstAlias< KSpace > ks, Clone< GeometricFunctor > fct, const int maxIter=0, const Scalar accuracy=0.1, const Scalar gamma=0.01)
Aim: This class specializes the Voronoi covariance measure for digital surfaces. It adds notably the ...
std::ostream & emphase()
bool init(const Point &lower, const Point &upper, bool isClosed)
Aim: A functor Matrix -> RealVector that returns the normal direction by diagonalizing the given cova...
Aim: A class for computing the Gauss digitization of some Euclidean shape, i.e. its intersection with...
void addValue(Quantity v)
Aim: model of CEuclideanOrientedShape concepts to create a shape from a polynomial.
int main(int argc, char **argv)
DGtal is the top-level namespace which contains all DGtal functions and types.
Aim: This concept describes an object that can process a range of surfels (that are supposed to belon...
Aim: An estimator on digital surfaces that returns the reference local geometric quantity. This is used for comparing estimators.
Aim: A functor Real -> Real that returns the 3d mean curvature by transforming the given volume...
std::ostream & info()
double unbiasedVariance() const
Quantity min() const
KhalimskySpaceND< 3, Integer > KSpace
Definition: StdDefs.h:146
Aim: This concept describes an object that can process a range over some generic digital surface so a...
Aim: A functor Surfel -> Quantity that returns the mean of absolute curvatures at given surfel: (abs(...
KSpace K
Quantity max() const
Aim: This class is a model of CCellularGridSpaceND. It represents the cubical grid as a cell complex...