DGtal  1.0.0
exampleHeatLaplace.cpp
1 
16 
23 
25 #include <DGtal/helpers/StdDefs.h>
26 
27 #include <DGtal/topology/DigitalSurface.h>
28 #include <DGtal/topology/DigitalSetBoundary.h>
29 #include <DGtal/topology/SetOfSurfels.h>
30 #include <DGtal/topology/LightImplicitDigitalSurface.h>
31 
32 #include <DGtal/math/linalg/EigenSupport.h>
33 
34 #include <DGtal/dec/DiscreteExteriorCalculus.h>
35 #include <DGtal/dec/DiscreteExteriorCalculusFactory.h>
36 #include <DGtal/dec/DiscreteExteriorCalculusSolver.h>
37 #include <DGtal/dec/VectorField.h>
38 
39 #include <DGtal/geometry/surfaces/estimation/LocalEstimatorFromSurfelFunctorAdapter.h>
40 #include <DGtal/geometry/surfaces/estimation/IIGeometricFunctors.h>
41 #include <DGtal/geometry/surfaces/estimation/IntegralInvariantCovarianceEstimator.h>
42 
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"
48 
49 #include <DGtal/images/IntervalForegroundPredicate.h>
50 #include <DGtal/images/imagesSetsUtils/SetFromImage.h>
51 
52 #include <DGtal/shapes/parametric/Ball3D.h>
53 #include <DGtal/shapes/Shapes.h>
54 #include <DGtal/shapes/GaussDigitizer.h>
55 
56 struct Options
57 {
58  double h;
59  double normal_radius;
60  double convolution_radius;
61 };
62 
64 
65 using namespace DGtal;
66 using namespace Eigen;
67 
69 
70 typedef Z3i::Space Space;
73 typedef Z3i::Point Point;
74 
76 
78 {
79  return RealPoint( a.norm(), atan2( a[1], a[0] ), acos( a[2] ) );
80 }
81 
83 std::function<double(const RealPoint&)> xx_function =
84  [] ( const RealPoint& p )
85 {
86  const RealPoint p_sphere = p / p.norm();
87  return p_sphere[0] * p_sphere[0];
88 };
89 
90 std::function<double(const RealPoint&)> xx_derivative =
91  [] ( const RealPoint& p )
92 {
93  const RealPoint p_s = cartesian_to_spherical( p );
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] ) );
96 };
98 
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)
103 {
104  trace.beginBlock("Laplacian 3D");
105 
106  trace.beginBlock("Extracting Digital Surface");
107 
108  typedef Z3i::KSpace KSpace;
109 
111  typedef GaussDigitizer<Z3i::Space, Shape> Digitizer;
112 
113  Digitizer digitizer;
114  digitizer.attach(shape);
115  digitizer.init(shape.getLowerBound() + Z3i::Vector(-1,-1,-1), shape.getUpperBound() + Z3i::Vector(1,1,1), options.h);
116 
117  Z3i::Domain domain = digitizer.getDomain();
118 
119  Z3i::KSpace kspace;
120  bool ok = kspace.init(domain.lowerBound(), domain.upperBound(), true);
121  if( !ok ) std::cerr << "KSpace init failed" << std::endl;
122 
124  typedef SetOfSurfels< KSpace, SurfelSet > MySetOfSurfels;
126  typedef SurfelAdjacency<KSpace::dimension> MySurfelAdjacency;
127 
128  MySurfelAdjacency surfAdj( true ); // interior in all directions.
129  MySetOfSurfels theSetOfSurfels( kspace, surfAdj );
130  Surfaces<KSpace>::sMakeBoundary( theSetOfSurfels.surfelSet(),
131  kspace, digitizer,
132  domain.lowerBound(),
133  domain.upperBound() );
134  MyDigitalSurface digSurf( theSetOfSurfels );
136  trace.info() << "Digital Surface has " << digSurf.size() << " surfels." << std::endl;
137 
138  trace.endBlock();
139 
140  trace.beginBlock("Initializing Normal Functor");
143  CanonicSCellEmbedder canonicSCellEmbedder(kspace);
144 
145  typedef functors::IINormalDirectionFunctor<Space> MyIINormalFunctor;
147 
148  const double radius = options.normal_radius * pow(options.h, 1. / 3.);
149 
150  MyIINormalFunctor normalFunctor;
151  normalFunctor.init(options.h, radius);
152 
153  MyIINormalEstimator normalEstimator(normalFunctor);
154  normalEstimator.attach(kspace, digitizer);
155  normalEstimator.setParams(radius / options.h);
156 
157  normalEstimator.init(options.h, digSurf.begin(), digSurf.end());
159  trace.endBlock();
160 
161  trace.beginBlock("Initializing DEC");
165 
166  const Calculus calculus = CalculusFactory::createFromNSCells<2>(digSurf.begin(), digSurf.end(), normalEstimator, options.h);
168  trace.info() << calculus << std::endl;
169 
170  trace.endBlock();
171  trace.beginBlock("Computing the input function");
173  Calculus::DualForm0 input_func(calculus);
174 
175  for(auto itb = digSurf.begin(), ite = digSurf.end(); itb != ite; itb++)
176  {
177  const Calculus::Index i_calc = calculus.getCellIndex( kspace.unsigns( *itb ) );
178  input_func.myContainer( i_calc ) = input_function( options.h * canonicSCellEmbedder( *itb ) );
179  }
181  trace.endBlock();
182 
183  trace.beginBlock("Computing the Laplace operator");
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;
190  trace.endBlock();
191 
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() );
198 
199  int i = 0;
200  for(auto itb = digSurf.begin(), ite = digSurf.end(); itb != ite; itb++)
201  {
202  const Calculus::Index i_calc = calculus.getCellIndex( kspace.unsigns( *itb ) );
203 
204  const RealPoint p = options.h * canonicSCellEmbedder( *itb );
205  const RealPoint p_normalized = p / p.norm();
206  const RealPoint p_s = cartesian_to_spherical( p );
207 
208  const double real_laplacian_value = result_function( p_normalized );
209  const double estimated_laplacian_value = laplace_result( i_calc );
210 
211  estimated_laplacian_values(i) = estimated_laplacian_value;
212  real_laplacian_values(i) = real_laplacian_value;
213 
214  error(i) = estimated_laplacian_value - real_laplacian_value;
215 
216  ++i;
217  }
218 
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;
221 
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;
224 
225  trace.endBlock();
226 }
227 
228 int main()
229 {
230  Options options;
231 
232  options.h = 0.1;
233  options.normal_radius = 2.0;
234  options.convolution_radius = 0.1;
235 
236  typedef Ball3D<Z3i::Space> Ball;
237  Ball ball(Point(0.0,0.0,0.0), 1.0);
238 
239  std::function<double(const RealPoint&, const RealPoint&)> l2_distance =
240  [](const RealPoint& a, const RealPoint& b) { return (a - b).norm(); };
241 
242  convergence<Ball>(options, ball, xx_function, xx_derivative);
243 
244  return 0;
245 }
void beginBlock(const std::string &keyword="")
Aim: DiscreteExteriorCalculus represents a calculus in the dec package. This is the main structure in...
MyDigitalSurface::SurfelSet SurfelSet
Trace trace
Definition: Common.h:144
std::set< SCell > SurfelSet
Preferred type for defining a set of surfels (always signed cells).
void attach(ConstAlias< EuclideanShape > shape)
Z3i::RealVector RealVector
const Point & upperBound() const
double endBlock()
Aim: Represents a set of n-1-cells in a nD space, together with adjacency relation between these cell...
Domain domain
Aim: This class implement an Integral Invariant estimator which computes for each surfel the covarian...
Aim: Implements basic operations that will be used in Point and Vector classes.
Definition: PointVector.h:165
PointVector< 3, double > RealPoint
Cell unsigns(const SCell &p) const
Creates an unsigned cell from a signed one.
static void sMakeBoundary(SCellSet &aBoundary, const KSpace &aKSpace, const PointPredicate &pp, const Point &aLowerBound, const Point &aUpperBound)
MessageStream error
Aim: This class provides static members to create DEC structures from various other DGtal structures.
RealPoint cartesian_to_spherical(const RealPoint3D &a)
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...
double norm(const NormType type=L_2) const
Z3i::KSpace KSpace
int main(int argc, char **argv)
DGtal is the top-level namespace which contains all DGtal functions and types.
MyPointD Point
Definition: testClone2.cpp:383
std::ostream & info()
const Point & lowerBound() const
Space::Vector Vector
Definition: StdDefs.h:169
Aim: A model of CDigitalSurfaceContainer which defines the digital surface as connected surfels....
Definition: SetOfSurfels.h:73
KSpace K
Aim: This class is a model of CCellularGridSpaceND. It represents the cubical grid as a cell complex,...
DigitalSurface< MyDigitalSurfaceContainer > MyDigitalSurface
Aim: Model of the concept StarShaped3D represents any Sphere in the space.
Definition: Ball3D.h:60