DGtal 1.3.0
Loading...
Searching...
No Matches
testVoronoiCovarianceMeasureOnSurface.cpp
1
31#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
62bool 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
70 typedef MPolynomial< 3, double > Polynomial3;
71 typedef MPolynomialReader<3, double> Polynomial3Reader;
72 typedef ImplicitPolynomial3Shape<Z3i::Space> ImplicitShape;
73 typedef GaussDigitizer< Z3i::Space, ImplicitShape > ImplicitDigitalShape;
77 typedef SurfaceContainer::Surfel Surfel;
79 typedef functors::HatPointFunction<Point,double> KernelFunction;
80 //typedef functors::BallConstantPointFunction<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*4.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( -20, -20, -20 );
100 Point p2( 20, 20, 20 );
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, 5.0 );
119 CountedPtr<VCMOnSurface> vcm_surface( new VCMOnSurface( ptrSurface, Pointels,
120 10.0, 5.0, chi, 1.5, 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." );
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." );
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.055555 ) < 0.05 ) ? 1 : 0;
258 nb++;
259 trace.info() << "(" << nbok << "/" << nb << ") "
260 << "mean of II curv is around 0.055555: "
261 << std::abs( stat_ii_curv.mean() - 0.055555 ) << std::endl;
262 nbok += ( std::abs( stat_vcm_curv.mean() - 0.055555 ) < 0.05 ) ? 1 : 0;
263 nb++;
264 trace.info() << "(" << nbok << "/" << nb << ") "
265 << "mean of VCM curv is around 0.055555: "
266 << std::abs( stat_vcm_curv.mean() - 0.055555 ) << std::endl;
267
268 trace.endBlock();
269 return nbok == nb;
270}
271
273// Standard services - public :
274
275int 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// //
Aim: Smart or simple const pointer on T. It can be a smart pointer based on reference counts or a sim...
Aim: Smart pointer based on reference counts.
Definition: CountedPtr.h:80
Aim: Represents a set of n-1-cells in a nD space, together with adjacency relation between these cell...
DigitalSurfaceContainer::SurfelConstIterator ConstIterator
Aim: implements separable l_p metrics with exact predicates.
Aim: A class for computing the Gauss digitization of some Euclidean shape, i.e. its intersection with...
Aim: model of CEuclideanOrientedShape concepts to create a shape from a polynomial.
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...
Aim: A model of CDigitalSurfaceContainer which defines the digital surface as the boundary of an impl...
Aim: This class converts a string polynomial expression in a multivariate polynomial.
Iterator read(Polynomial &p, Iterator begin, Iterator end)
Aim: Represents a multivariate polynomial, i.e. an element of , where K is some ring or field.
Definition: MPolynomial.h:955
auto dot(const PointVector< dim, OtherComponent, OtherStorage > &v) const -> decltype(DGtal::dotProduct(*this, v))
Dot product with a PointVector.
Aim: This class processes a set of sample values for one variable and can then compute different stat...
Definition: Statistic.h:70
double mean() const
double unbiasedVariance() const
Quantity max() const
Quantity min() const
void addValue(Quantity v)
Aim: A utility class for constructing surfaces (i.e. set of (n-1)-cells).
Definition: Surfaces.h:79
Aim: Represent adjacencies between surfel elements, telling if it follows an interior to exterior ord...
void beginBlock(const std::string &keyword="")
std::ostream & emphase()
std::ostream & info()
double endBlock()
Aim: An estimator on digital surfaces that returns the reference local geometric quantity....
void setParams(ConstAlias< KSpace > ks, Clone< GeometricFunctor > fct, const int maxIter=20, const Scalar accuracy=0.0001, const Scalar gamma=0.5)
Aim: This class adapts a VoronoiCovarianceMeasureOnDigitalSurface to be a model of CDigitalSurfaceLoc...
Aim: This class specializes the Voronoi covariance measure for digital surfaces. It adds notably the ...
Aim: A functor Matrix -> RealVector that returns the normal direction by diagonalizing the given cova...
SH3::DigitalSurface Surface
MyDigitalSurface::ConstIterator ConstIterator
Z3i this namespace gathers the standard of types for 3D imagery.
DGtal is the top-level namespace which contains all DGtal functions and types.
Trace trace
Definition: Common.h:154
STL namespace.
Aim: This concept describes an object that can process a range over some generic digital surface so a...
Aim: This concept describes an object that can process a range of surfels (that are supposed to belon...
Aim: A functor Real -> Real that returns the 3d mean curvature by transforming the given volume....
Aim: A functor RealPoint -> Quantity that returns the mean curvature at given point.
Aim: A functor RealPoint -> Quantity that returns the normal vector at given point.
Aim: A functor Surfel -> Quantity that returns the mean of absolute curvatures at given surfel: (abs(...
Aim: A functor Surfel -> Quantity that returns the outer normal vector at given surfel.
int main()
Definition: testBits.cpp:56
MyPointD Point
Definition: testClone2.cpp:383
KSpace K