DGtal 1.3.0
Loading...
Searching...
No Matches
testIntegralInvariantCovarianceEstimator.cpp
Go to the documentation of this file.
1
32#include <iostream>
33#include <tuple>
34
35#include "DGtal/base/Common.h"
36
38#include "DGtal/shapes/implicit/ImplicitBall.h"
39
41#include "DGtal/shapes/GaussDigitizer.h"
42#include "DGtal/topology/LightImplicitDigitalSurface.h"
43#include "DGtal/topology/DigitalSurface.h"
44#include "DGtal/graph/DepthFirstVisitor.h"
45#include "DGtal/graph/GraphVisitorRange.h"
46
48#include "DGtal/geometry/surfaces/estimation/IIGeometricFunctors.h"
49#include "DGtal/geometry/surfaces/estimation/IntegralInvariantCovarianceEstimator.h"
50
51
53
54
55using namespace DGtal;
56
58// Functions for testing class IntegralInvariantCovarianceEstimator and IIGeometricFunctor.
60
61bool testGaussianCurvature3d( double h, double delta )
62{
68 typedef GraphVisitorRange< Visitor > VisitorRange;
69 typedef VisitorRange::ConstIterator VisitorConstIterator;
70
71 typedef functors::IIGaussianCurvature3DFunctor<Z3i::Space> MyIICurvatureFunctor;
73 typedef MyIICurvatureFunctor::Value Value;
74
75 double re = 5.0;
76 double radius = 5.0;
77 double realValue = 1.0/(radius * radius);
78
79 trace.beginBlock( "Shape initialisation ..." );
80
81 ImplicitShape ishape( Z3i::RealPoint( 0, 0, 0 ), radius );
82 DigitalShape dshape;
83 dshape.attach( ishape );
84 dshape.init( Z3i::RealPoint( -10.0, -10.0, -10.0 ), Z3i::RealPoint( 10.0, 10.0, 10.0 ), h );
85
87 if ( !K.init( dshape.getLowerBound(), dshape.getUpperBound(), true ) )
88 {
89 trace.error() << "Problem with Khalimsky space" << std::endl;
90 return false;
91 }
92
94 Boundary boundary( K, dshape, SurfelAdjacency<Z3i::KSpace::dimension>( true ), bel );
95 MyDigitalSurface surf ( boundary );
96
98
99 trace.beginBlock( "Curvature estimator initialisation ...");
100
101 VisitorRange range( new Visitor( surf, *surf.begin() ));
102 VisitorConstIterator ibegin = range.begin();
103 VisitorConstIterator iend = range.end();
104
105 MyIICurvatureFunctor curvatureFunctor;
106 curvatureFunctor.init( h, re );
107
108 MyIICurvatureEstimator curvatureEstimator( curvatureFunctor );
109 curvatureEstimator.attach( K, dshape );
110 curvatureEstimator.setParams( re/h );
111 curvatureEstimator.init( h, ibegin, iend );
112
113 trace.endBlock();
114
115 trace.beginBlock( "Curvature estimator evaluation ...");
116
117 std::vector< Value > results;
118 std::back_insert_iterator< std::vector< Value > > resultsIt( results );
119 curvatureEstimator.eval( ibegin, iend, resultsIt );
120
121 trace.endBlock();
122
123 trace.beginBlock ( "Comparing results of integral invariant 3D Gaussian curvature ..." );
124
125 double mean = 0.0;
126 unsigned int rsize = static_cast<unsigned int>(results.size());
127
128 if( rsize == 0 )
129 {
130 trace.error() << "ERROR: surface is empty" << std::endl;
131 trace.endBlock();
132 return false;
133 }
134
135 for ( unsigned int i = 0; i < rsize; ++i )
136 {
137 mean += results[ i ];
138 }
139 mean /= rsize;
140
141 if( mean != mean ) //NaN
142 {
143 trace.error() << "ERROR: result is NaN" << std::endl;
144 trace.endBlock();
145 return false;
146 }
147
148 double v = std::abs ( realValue - mean );
149
150 trace.warning() << "True value: " << realValue << std::endl;
151 trace.warning() << "Mean value: " << mean << std::endl;
152 trace.warning() << "Delta: " << delta << " |true - mean|: " << v << std::endl;
153
154 if( v > delta )
155 {
156 trace.endBlock();
157 return false;
158 }
159 trace.endBlock();
160 return true;
161}
162
164{
170 typedef GraphVisitorRange< Visitor > VisitorRange;
171 typedef VisitorRange::ConstIterator VisitorConstIterator;
172
173 typedef functors::IIPrincipalCurvatures3DFunctor<Z3i::Space> MyIICurvatureFunctor;
175 typedef MyIICurvatureFunctor::Value Value;
176
179 typedef MyIICurvatureFunctorTensor::Value ValueTensor;
180
181 double re = 5.0;
182 double radius = 5.0;
183
184 trace.beginBlock( "Shape initialisation ..." );
185
186 ImplicitShape ishape( Z3i::RealPoint( 0, 0, 0 ), radius );
187 DigitalShape dshape;
188 dshape.attach( ishape );
189 dshape.init( Z3i::RealPoint( -10.0, -10.0, -10.0 ), Z3i::RealPoint( 10.0, 10.0, 10.0 ), h );
190
192 if ( !K.init( dshape.getLowerBound(), dshape.getUpperBound(), true ) )
193 {
194 trace.error() << "Problem with Khalimsky space" << std::endl;
195 return false;
196 }
197
199 Boundary boundary( K, dshape, SurfelAdjacency<Z3i::KSpace::dimension>( true ), bel );
200 MyDigitalSurface surf ( boundary );
201
202 trace.endBlock();
203
204 trace.beginBlock( "Curvature estimator initialisation ...");
205
206 VisitorRange range( new Visitor( surf, *surf.begin() ));
207 VisitorConstIterator ibegin = range.begin();
208 VisitorConstIterator iend = range.end();
209
210 MyIICurvatureFunctor curvatureFunctor;
211 curvatureFunctor.init( h, re );
212 MyIICurvatureEstimator curvatureEstimator( curvatureFunctor );
213 curvatureEstimator.attach( K, dshape );
214 curvatureEstimator.setParams( re/h );
215 curvatureEstimator.init( h, ibegin, iend );
216
217 MyIICurvatureFunctorTensor curvatureFunctorTensor;
218 curvatureFunctorTensor.init( h, re );
219 MyIICurvatureEstimatorTensor curvatureEstimatorTensor( curvatureFunctorTensor );
220 curvatureEstimatorTensor.attach( K, dshape );
221 curvatureEstimatorTensor.setParams( re/h );
222 curvatureEstimatorTensor.init( h, ibegin, iend );
223
224 trace.endBlock();
225
226 trace.beginBlock( "Curvature estimator evaluation ...");
227
228 std::vector< Value > results;
229 std::back_insert_iterator< std::vector< Value > > resultsIt( results );
230 curvatureEstimator.eval( ibegin, iend, resultsIt );
231 trace.endBlock();
232
233
234
235 trace.beginBlock( "Checking CurvaturesAndDirections functor");
236 Value val = curvatureEstimator.eval( surf.begin() );
237 ValueTensor valTensor = curvatureEstimatorTensor.eval( surf.begin() );
238 bool ok = (val.first == std::get<0>(valTensor));
239 ok &= (val.second == std::get<1>(valTensor));
240 if (!ok)
241 trace.error()<< "Error comparing principal curvatures between two different functors."<<std::endl;
242
243 trace.endBlock();
244
245
246 trace.beginBlock ( "Comparing results of integral invariant 3D Gaussian curvature ..." );
247
248 unsigned int error_order = 0;
249 unsigned int rsize = static_cast<unsigned int>(results.size());
250
251 if( rsize == 0 )
252 {
253 trace.error() << "ERROR: surface is empty" << std::endl;
254 trace.endBlock();
255 return false;
256 }
257
258 for ( unsigned int i = 0; i < rsize; ++i )
259 {
260 if( std::abs(results[i].first) < std::abs(results[i].second) )
261 {
262 ++error_order;
263 }
264 }
265
266 trace.warning() << "Error order: " << error_order << std::endl;
267 trace.warning() << "If not equals to 0, something is wrong..." << std::endl;
268
269 trace.endBlock();
270
271 if( error_order != 0 )
272 {
273 return false;
274 }
275 return true;
276}
277
279// Standard services - public :
280
281int main( int /*argc*/, char** /*argv*/ )
282{
283 trace.beginBlock ( "Testing class IntegralInvariantCovarianceEstimator and 3d functors" );
284 bool res = testGaussianCurvature3d( 0.6, 0.007 ) && testPrincipalCurvatures3d( 0.6 );
285 trace.emphase() << ( res ? "Passed." : "Error." ) << std::endl;
286 trace.endBlock();
287 return res ? 0 : 1;
288}
289// //
Aim: This class is useful to perform a depth-first exploration of a graph given a starting point or s...
Aim: Represents a set of n-1-cells in a nD space, together with adjacency relation between these cell...
ConstIterator begin() const
Aim: A class for computing the Gauss digitization of some Euclidean shape, i.e. its intersection with...
const Point & getUpperBound() const
void attach(ConstAlias< EuclideanShape > shape)
const Point & getLowerBound() const
void init(const RealPoint &xLow, const RealPoint &xUp, typename RealVector::Component gridStep)
Aim: Transforms a graph visitor into a single pass input range.
Aim: model of CEuclideanOrientedShape and CEuclideanBoundedShape concepts to create a ball in nD....
Definition: ImplicitBall.h:65
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 is a model of CCellularGridSpaceND. It represents the cubical grid as a cell complex,...
bool init(const Point &lower, const Point &upper, bool isClosed)
Specifies the upper and lower bounds for the maximal cells in this space.
Aim: A model of CDigitalSurfaceContainer which defines the digital surface as the boundary of an impl...
Aim: Implements basic operations that will be used in Point and Vector classes.
Definition: PointVector.h:593
static SCell findABel(const KSpace &K, const PointPredicate &pp, unsigned int nbtries=1000)
Aim: Represent adjacencies between surfel elements, telling if it follows an interior to exterior ord...
void beginBlock(const std::string &keyword="")
std::ostream & warning()
std::ostream & emphase()
std::ostream & error()
double endBlock()
Aim: A functor Matrix -> std::pair<RealVector,RealVector> that returns the first and the second princ...
DigitalSurface< MyDigitalSurfaceContainer > MyDigitalSurface
BreadthFirstVisitor< MyDigitalSurface > Visitor
DGtal is the top-level namespace which contains all DGtal functions and types.
Trace trace
Definition: Common.h:154
Represents a signed cell in a cellular grid space by its Khalimsky coordinates and a boolean value.
Aim: A functor Matrix -> Real that returns the Gaussian curvature by diagonalizing the given covarian...
Aim: A functor Matrix -> std::pair<Real,Real> that returns the first and the second principal curvatu...
int main()
Definition: testBits.cpp:56
KSpace K
bool testPrincipalCurvatures3d(double h)
bool testGaussianCurvature3d(double h, double delta)
std::default_random_engine re