25 #include <DGtal/helpers/StdDefs.h>
27 #include <DGtal/topology/DigitalSurface.h>
28 #include <DGtal/topology/DigitalSetBoundary.h>
29 #include <DGtal/topology/SetOfSurfels.h>
30 #include <DGtal/topology/LightImplicitDigitalSurface.h>
32 #include <DGtal/math/linalg/EigenSupport.h>
34 #include <DGtal/dec/DiscreteExteriorCalculus.h>
35 #include <DGtal/dec/DiscreteExteriorCalculusFactory.h>
36 #include <DGtal/dec/DiscreteExteriorCalculusSolver.h>
37 #include <DGtal/dec/VectorField.h>
39 #include <DGtal/geometry/surfaces/estimation/LocalEstimatorFromSurfelFunctorAdapter.h>
40 #include <DGtal/geometry/surfaces/estimation/IIGeometricFunctors.h>
41 #include <DGtal/geometry/surfaces/estimation/IntegralInvariantCovarianceEstimator.h>
43 #include <DGtal/io/readers/GenericReader.h>
44 #include <DGtal/io/colormaps/ColorBrightnessColorMap.h>
45 #include <DGtal/io/viewers/Viewer3D.h>
46 #include <DGtal/io/colormaps/TickedColorMap.h>
47 #include "DGtal/io/readers/GenericReader.h"
49 #include <DGtal/images/IntervalForegroundPredicate.h>
50 #include <DGtal/images/imagesSetsUtils/SetFromImage.h>
52 #include <DGtal/shapes/parametric/Ball3D.h>
53 #include <DGtal/shapes/Shapes.h>
54 #include <DGtal/shapes/GaussDigitizer.h>
60 double convolution_radius;
65 using namespace DGtal;
66 using namespace Eigen;
79 return RealPoint( a.norm(), atan2( a[1], a[0] ), acos( a[2] ) );
83 std::function<double(
const RealPoint&)> xx_function =
87 return p_sphere[0] * p_sphere[0];
90 std::function<double(
const RealPoint&)> xx_derivative =
94 return 2. * cos( p_s[1] ) * cos( p_s[1] ) * ( 2 * cos( p_s[2] ) * cos( p_s[2] ) - sin( p_s[2] ) * sin( p_s[2] ) )
95 + 2. * ( sin( p_s[1] ) * sin( p_s[1] ) - cos( p_s[1] ) * cos( p_s[1] ) );
99 template <
typename Shape>
100 void convergence(
const Options& options,
Shape& shape,
101 const std::function<
double(
const RealPoint&) >& input_function,
102 const std::function<
double(
const RealPoint&) >& result_function)
121 if( !ok ) std::cerr <<
"KSpace init failed" << std::endl;
128 MySurfelAdjacency surfAdj(
true );
129 MySetOfSurfels theSetOfSurfels( kspace, surfAdj );
136 trace.
info() <<
"Digital Surface has " << digSurf.size() <<
" surfels." << std::endl;
148 const double radius = options.normal_radius * pow(options.h, 1. / 3.);
150 MyIINormalFunctor normalFunctor;
151 normalFunctor.init(options.h, radius);
153 MyIINormalEstimator normalEstimator(normalFunctor);
154 normalEstimator.attach(kspace, digitizer);
155 normalEstimator.setParams(radius / options.h);
157 normalEstimator.init(options.h, digSurf.begin(), digSurf.end());
166 const Calculus calculus = CalculusFactory::createFromNSCells<2>(digSurf.begin(), digSurf.end(), normalEstimator, options.h);
173 Calculus::DualForm0 input_func(calculus);
175 for(
auto itb = digSurf.begin(), ite = digSurf.end(); itb != ite; itb++)
177 const Calculus::Index i_calc = calculus.getCellIndex( kspace.
unsigns( *itb ) );
178 input_func.myContainer( i_calc ) = input_function( options.h * canonicSCellEmbedder( *itb ) );
185 const double t = options.convolution_radius * pow(options.h, 2. / 3.);
186 const double K = log( - log1p( t ) ) + 2.;
187 const Calculus::DualIdentity0 laplace = calculus.heatLaplace<
DUAL>(options.h, t,
K);
189 trace.
info() <<
"Matrix has " << ((double)laplace.myContainer.nonZeros() / (double)laplace.myContainer.size() * 100.) <<
"% of non-zeros elements." << std::endl;
193 const Eigen::VectorXd laplace_result = (laplace * input_func).myContainer;
195 Eigen::VectorXd
error( digSurf.size() );
196 Eigen::VectorXd real_laplacian_values( digSurf.size() );
197 Eigen::VectorXd estimated_laplacian_values( digSurf.size() );
200 for(
auto itb = digSurf.begin(), ite = digSurf.end(); itb != ite; itb++)
202 const Calculus::Index i_calc = calculus.getCellIndex( kspace.
unsigns( *itb ) );
204 const RealPoint p = options.h * canonicSCellEmbedder( *itb );
208 const double real_laplacian_value = result_function( p_normalized );
209 const double estimated_laplacian_value = laplace_result( i_calc );
211 estimated_laplacian_values(i) = estimated_laplacian_value;
212 real_laplacian_values(i) = real_laplacian_value;
214 error(i) = estimated_laplacian_value - real_laplacian_value;
219 trace.
info() <<
"Estimated Laplacian Range : " << estimated_laplacian_values.minCoeff() <<
" / " << estimated_laplacian_values.maxCoeff() << std::endl;
220 trace.
info() <<
"Real Laplacian Range : " << real_laplacian_values.minCoeff() <<
" / " << real_laplacian_values.maxCoeff() << std::endl;
222 trace.
info() <<
"h = " << options.h <<
" t = " << t <<
" K = " <<
K << std::endl;
223 trace.
info() <<
"Mean error = " <<
error.array().abs().mean() <<
" max error = " <<
error.array().abs().maxCoeff() << std::endl;
233 options.normal_radius = 2.0;
234 options.convolution_radius = 0.1;
242 convergence<Ball>(options, ball, xx_function, xx_derivative);
RealPoint getLowerBound() const
RealPoint getUpperBound() const
Aim: Model of the concept StarShaped represents any circle in the plane.
Aim: Model of the concept StarShaped3D represents any Sphere in the space.
Aim: Represents a set of n-1-cells in a nD space, together with adjacency relation between these cell...
Aim: This class provides static members to create DEC structures from various other DGtal structures.
Aim: DiscreteExteriorCalculus represents a calculus in the dec package. This is the main structure in...
Aim: A class for computing the Gauss digitization of some Euclidean shape, i.e. its intersection with...
void attach(ConstAlias< EuclideanShape > shape)
const Point & lowerBound() const
const Point & upperBound() const
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,...
std::set< SCell > SurfelSet
Preferred type for defining a set of surfels (always signed cells).
bool init(const Point &lower, const Point &upper, bool isClosed)
Specifies the upper and lower bounds for the maximal cells in this space.
Cell unsigns(const SCell &p) const
Creates an unsigned cell from a signed one.
Aim: Implements basic operations that will be used in Point and Vector classes.
double norm(const NormType type=L_2) const
Aim: A model of CDigitalSurfaceContainer which defines the digital surface as connected surfels....
static void sMakeBoundary(SCellSet &aBoundary, const KSpace &aKSpace, const PointPredicate &pp, const Point &aLowerBound, const Point &aUpperBound)
void beginBlock(const std::string &keyword="")
Aim: A functor Matrix -> RealVector that returns the normal direction by diagonalizing the given cova...
DigitalSurface< MyDigitalSurfaceContainer > MyDigitalSurface
MyDigitalSurface::SurfelSet SurfelSet
DGtal is the top-level namespace which contains all DGtal functions and types.
RealPoint cartesian_to_spherical(const RealPoint3D &a)
Z3i::RealVector RealVector
int main(int argc, char **argv)