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" 62 bool testVoronoiCovarianceMeasureOnSurface()
64 unsigned int nbok = 0;
67 using namespace DGtal;
77 typedef SurfaceContainer::Surfel Surfel;
84 std::string poly_str =
"-81.0+x^2+y^2+z^2";
89 std::string::const_iterator iter = reader.
read( poly, poly_str.begin(), poly_str.end() );
90 if ( iter != poly_str.end() )
92 std::cerr <<
"ERROR: I read only <" 93 << poly_str.substr( 0, iter - poly_str.begin() )
94 <<
">, and I built P=" << poly << std::endl;
99 Point p1( -10, -10, -10 );
100 Point p2( 10, 10, 10 );
102 nbok +=
K.
init( p1, p2,
true ) ? 1 : 0;
104 trace.
info() <<
"(" << nbok <<
"/" << nb <<
") " 105 <<
"K.init() is ok" << std::endl;
109 dshape->attach( *shape );
110 dshape->init( p1, p2, 1.0 );
112 SurfaceContainer* surfaceContainer =
new SurfaceContainer
118 KernelFunction chi( 1.0, 7.0 );
120 5.0, 7.0, chi, 7.0, Metric(),
true ) );
126 KernelFunction, NormalVectorFunctor> VCMNormalEstimator;
127 VCMNormalEstimator estimator( vcm_surface );
128 estimator.init( 1.0, ptrSurface->begin(), ptrSurface->end() );
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() );
141 typedef ShapeGeometricFunctors::ShapeNormalVectorFunctor<ImplicitShape> NormalFunctor;
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() );
155 for (
ConstIterator it = ptrSurface->begin(), itE = ptrSurface->end(); it != itE; ++it )
158 RealVector n_true = true_estimator.eval( it );
159 RealVector n_triv = - vcm_surface->mapSurfel2Normals().find( *it )->second.trivialNormal;
163 if ( n_ii.dot( n_triv ) < 0 ) n_ii = -n_ii;
164 error_ii_true.
addValue( n_ii.dot( n_true ) );
169 trace.
info() <<
"VCM/true cos angle avg = " << error_true.
mean() << std::endl;
171 trace.
info() <<
"II/true cos angle avg = " << error_ii_true.
mean() << std::endl;
173 trace.
info() <<
"triv/true cos angle avg = " << error_triv_true.
mean() << std::endl;
175 nbok += error_true.
mean() > 0.95 ? 1 : 0;
177 trace.
info() <<
"(" << nbok <<
"/" << nb <<
") " 178 <<
"VCM/true cos angle avg > 0.95" << std::endl;
181 trace.
info() <<
"(" << nbok <<
"/" << nb <<
") " 182 <<
"VCM/true cos angle dev < 0.05" << std::endl;
188 nbok += error_ii_true.
mean() > 0.95 ? 1 : 0;
190 trace.
info() <<
"(" << nbok <<
"/" << nb <<
") " 191 <<
"II/true cos angle avg > 0.95" << std::endl;
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;
198 trace.
info() <<
"(" << nbok <<
"/" << nb <<
") " 199 <<
"II/true is closer to 1.0 than triv/true." << std::endl;
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() );
215 KernelFunction, VCMCurvatureFunctor> VCMCurvatureEstimator;
216 VCMCurvatureEstimator vcm_curv_estimator( vcm_surface );
217 vcm_curv_estimator.init( 1.0, ptrSurface->begin(), ptrSurface->end() );
221 typedef ShapeGeometricFunctors::ShapeMeanCurvatureFunctor<ImplicitShape> CurvatureFunctor;
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() );
239 for (
ConstIterator it = ptrSurface->begin(), itE = ptrSurface->end(); it != itE; ++it )
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 ) );
249 <<
" min=" << stat_ii_curv.
min()
250 <<
" max=" << stat_ii_curv.
max() << std::endl;
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;
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;
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;
275 int main(
int ,
char** )
278 using namespace DGtal;
280 bool res = testVoronoiCovarianceMeasureOnSurface();
281 trace.
emphase() << ( res ?
"Passed." :
"Error." ) << endl;
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.
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.
Aim: A functor Surfel -> Quantity that returns the outer normal vector at given surfel.
auto dot(const PointVector< dim, OtherComponent, OtherStorage > &v) const -> decltype(DGtal::dotProduct(*this, v))
Dot product with a PointVector.
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.
Aim: This class processes a set of sample values for one variable and can then compute different stat...
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.
Aim: This class specializes the Voronoi covariance measure for digital surfaces. It adds notably the ...
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 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)
void setParams(ConstAlias< KSpace > ks, Clone< GeometricFunctor > fct, const int maxIter=20, const Scalar accuracy=0.0001, const Scalar gamma=0.5)
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....
Aim: A functor Real -> Real that returns the 3d mean curvature by transforming the given volume....
double unbiasedVariance() const
KhalimskySpaceND< 3, Integer > KSpace
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(...
Aim: This class is a model of CCellularGridSpaceND. It represents the cubical grid as a cell complex,...