DGtal  1.3.beta
dec/exampleHeatLaplace.cpp

Example of the Heat Laplace Operator.

See also
Hands on the Operator in DGtal
#include <DGtal/helpers/StdDefs.h>
#include <DGtal/topology/DigitalSurface.h>
#include <DGtal/topology/DigitalSetBoundary.h>
#include <DGtal/topology/SetOfSurfels.h>
#include <DGtal/topology/LightImplicitDigitalSurface.h>
#include <DGtal/math/linalg/EigenSupport.h>
#include <DGtal/dec/DiscreteExteriorCalculus.h>
#include <DGtal/dec/DiscreteExteriorCalculusFactory.h>
#include <DGtal/dec/DiscreteExteriorCalculusSolver.h>
#include <DGtal/dec/VectorField.h>
#include <DGtal/geometry/surfaces/estimation/LocalEstimatorFromSurfelFunctorAdapter.h>
#include <DGtal/geometry/surfaces/estimation/IIGeometricFunctors.h>
#include <DGtal/geometry/surfaces/estimation/IntegralInvariantCovarianceEstimator.h>
#include <DGtal/io/readers/GenericReader.h>
#include <DGtal/io/colormaps/ColorBrightnessColorMap.h>
#include <DGtal/io/viewers/Viewer3D.h>
#include <DGtal/io/colormaps/TickedColorMap.h>
#include "DGtal/io/readers/GenericReader.h"
#include <DGtal/images/IntervalForegroundPredicate.h>
#include <DGtal/images/imagesSetsUtils/SetFromImage.h>
#include <DGtal/shapes/parametric/Ball3D.h>
#include <DGtal/shapes/Shapes.h>
#include <DGtal/shapes/GaussDigitizer.h>
struct Options
{
double h;
double normal_radius;
double convolution_radius;
};
using namespace DGtal;
using namespace Eigen;
typedef Z3i::Space Space;
typedef Z3i::Point Point;
{
return RealPoint( a.norm(), atan2( a[1], a[0] ), acos( a[2] ) );
}
std::function<double(const RealPoint&)> xx_function =
[] ( const RealPoint& p )
{
const RealPoint p_sphere = p / p.norm();
return p_sphere[0] * p_sphere[0];
};
std::function<double(const RealPoint&)> xx_derivative =
[] ( const RealPoint& p )
{
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] ) )
+ 2. * ( sin( p_s[1] ) * sin( p_s[1] ) - cos( p_s[1] ) * cos( p_s[1] ) );
};
template <typename Shape>
void convergence(const Options& options, Shape& shape,
const std::function< double(const RealPoint&) >& input_function,
const std::function< double(const RealPoint&) >& result_function)
{
trace.beginBlock("Laplacian 3D");
trace.beginBlock("Extracting Digital Surface");
Digitizer digitizer;
digitizer.attach(shape);
digitizer.init(shape.getLowerBound() + Z3i::Vector(-1,-1,-1), shape.getUpperBound() + Z3i::Vector(1,1,1), options.h);
Z3i::Domain domain = digitizer.getDomain();
Z3i::KSpace kspace;
bool ok = kspace.init(domain.lowerBound(), domain.upperBound(), true);
if( !ok ) std::cerr << "KSpace init failed" << std::endl;
typedef SetOfSurfels< KSpace, SurfelSet > MySetOfSurfels;
typedef SurfelAdjacency<KSpace::dimension> MySurfelAdjacency;
MySurfelAdjacency surfAdj( true ); // interior in all directions.
MySetOfSurfels theSetOfSurfels( kspace, surfAdj );
Surfaces<KSpace>::sMakeBoundary( theSetOfSurfels.surfelSet(),
kspace, digitizer,
MyDigitalSurface digSurf( theSetOfSurfels );
trace.info() << "Digital Surface has " << digSurf.size() << " surfels." << std::endl;
trace.beginBlock("Initializing Normal Functor");
CanonicSCellEmbedder canonicSCellEmbedder(kspace);
typedef functors::IINormalDirectionFunctor<Space> MyIINormalFunctor;
const double radius = options.normal_radius * pow(options.h, 1. / 3.);
MyIINormalFunctor normalFunctor;
normalFunctor.init(options.h, radius);
MyIINormalEstimator normalEstimator(normalFunctor);
normalEstimator.attach(kspace, digitizer);
normalEstimator.setParams(radius / options.h);
normalEstimator.init(options.h, digSurf.begin(), digSurf.end());
trace.beginBlock("Initializing DEC");
const Calculus calculus = CalculusFactory::createFromNSCells<2>(digSurf.begin(), digSurf.end(), normalEstimator, options.h);
trace.info() << calculus << std::endl;
trace.beginBlock("Computing the input function");
Calculus::DualForm0 input_func(calculus);
for(auto itb = digSurf.begin(), ite = digSurf.end(); itb != ite; itb++)
{
const Calculus::Index i_calc = calculus.getCellIndex( kspace.unsigns( *itb ) );
input_func.myContainer( i_calc ) = input_function( options.h * canonicSCellEmbedder( *itb ) );
}
trace.beginBlock("Computing the Laplace operator");
const double t = options.convolution_radius * pow(options.h, 2. / 3.);
const double K = log( - log1p( t ) ) + 2.;
const Calculus::DualIdentity0 laplace = calculus.heatLaplace<DUAL>(options.h, t, K);
trace.info() << "Matrix has " << ((double)laplace.myContainer.nonZeros() / (double)laplace.myContainer.size() * 100.) << "% of non-zeros elements." << std::endl;
const Eigen::VectorXd laplace_result = (laplace * input_func).myContainer;
Eigen::VectorXd error( digSurf.size() );
Eigen::VectorXd real_laplacian_values( digSurf.size() );
Eigen::VectorXd estimated_laplacian_values( digSurf.size() );
int i = 0;
for(auto itb = digSurf.begin(), ite = digSurf.end(); itb != ite; itb++)
{
const Calculus::Index i_calc = calculus.getCellIndex( kspace.unsigns( *itb ) );
const RealPoint p = options.h * canonicSCellEmbedder( *itb );
const RealPoint p_normalized = p / p.norm();
const double real_laplacian_value = result_function( p_normalized );
const double estimated_laplacian_value = laplace_result( i_calc );
estimated_laplacian_values(i) = estimated_laplacian_value;
real_laplacian_values(i) = real_laplacian_value;
error(i) = estimated_laplacian_value - real_laplacian_value;
++i;
}
trace.info() << "Estimated Laplacian Range : " << estimated_laplacian_values.minCoeff() << " / " << estimated_laplacian_values.maxCoeff() << std::endl;
trace.info() << "Real Laplacian Range : " << real_laplacian_values.minCoeff() << " / " << real_laplacian_values.maxCoeff() << std::endl;
trace.info() << "h = " << options.h << " t = " << t << " K = " << K << std::endl;
trace.info() << "Mean error = " << error.array().abs().mean() << " max error = " << error.array().abs().maxCoeff() << std::endl;
}
int main()
{
Options options;
options.h = 0.1;
options.normal_radius = 2.0;
options.convolution_radius = 0.1;
Ball ball(Point(0.0,0.0,0.0), 1.0);
std::function<double(const RealPoint&, const RealPoint&)> l2_distance =
[](const RealPoint& a, const RealPoint& b) { return (a - b).norm(); };
convergence<Ball>(options, ball, xx_function, xx_derivative);
return 0;
}
Ball
Ball2D< Space > Ball
Definition: testShapeMoveCenter.cpp:58
DGtal::GaussDigitizer
Aim: A class for computing the Gauss digitization of some Euclidean shape, i.e. its intersection with...
Definition: GaussDigitizer.h:79
DGtal::Astroid2D< Space >
DGtal::DiscreteExteriorCalculus
Aim: DiscreteExteriorCalculus represents a calculus in the dec package. This is the main structure in...
Definition: DiscreteExteriorCalculus.h:97
DGtal::HyperRectDomain< Space >
DGtal::Trace::endBlock
double endBlock()
DGtal::KhalimskySpaceND::SurfelSet
std::set< SCell > SurfelSet
Preferred type for defining a set of surfels (always signed cells).
Definition: KhalimskySpaceND.h:450
DGtal::DigitalSurface
Aim: Represents a set of n-1-cells in a nD space, together with adjacency relation between these cell...
Definition: DigitalSurface.h:139
DGtal::SurfelAdjacency< KSpace::dimension >
DGtal::HyperRectDomain::upperBound
const Point & upperBound() const
DGtal::KhalimskySpaceND::init
bool init(const Point &lower, const Point &upper, bool isClosed)
Specifies the upper and lower bounds for the maximal cells in this space.
DGtal::Astroid2D::getLowerBound
RealPoint getLowerBound() const
Definition: Astroid2D.h:122
RealVector
Z3i::RealVector RealVector
Definition: sphereCotangentLaplaceOperator.cpp:71
DGtal::trace
Trace trace
Definition: Common.h:154
DGtal::GaussDigitizer::attach
void attach(ConstAlias< EuclideanShape > shape)
K
KSpace K
Definition: testCubicalComplex.cpp:62
DGtal::Trace::beginBlock
void beginBlock(const std::string &keyword="")
DGtal::IntegralInvariantCovarianceEstimator
Aim: This class implement an Integral Invariant estimator which computes for each surfel the covarian...
Definition: IntegralInvariantCovarianceEstimator.h:115
DGtal::functors::IINormalDirectionFunctor
Aim: A functor Matrix -> RealVector that returns the normal direction by diagonalizing the given cova...
Definition: IIGeometricFunctors.h:67
cartesian_to_spherical
RealPoint cartesian_to_spherical(const RealPoint3D &a)
Definition: sphereCotangentLaplaceOperator.cpp:74
DGtal::Ball2D
Aim: Model of the concept StarShaped represents any circle in the plane.
Definition: Ball2D.h:60
DGtal::SpaceND
Definition: SpaceND.h:95
KSpace
Z3i::KSpace KSpace
Definition: testArithmeticalDSSComputerOnSurfels.cpp:48
DGtal::Trace::info
std::ostream & info()
DGtal::Surfaces::sMakeBoundary
static void sMakeBoundary(SCellSet &aBoundary, const KSpace &aKSpace, const PointPredicate &pp, const Point &aLowerBound, const Point &aUpperBound)
DGtal
DGtal is the top-level namespace which contains all DGtal functions and types.
DGtal::DiscreteExteriorCalculusFactory
Aim: This class provides static members to create DEC structures from various other DGtal structures.
Definition: DiscreteExteriorCalculus.h:70
main
int main(int argc, char **argv)
Definition: testArithmeticDSS-benchmark.cpp:147
LibBoard::error
MessageStream error
DGtal::PointVector::norm
double norm(const NormType type=L_2) const
domain
Domain domain
Definition: testProjection.cpp:88
DGtal::PointVector
Aim: Implements basic operations that will be used in Point and Vector classes.
Definition: PointVector.h:165
DGtal::SetOfSurfels
Aim: A model of CDigitalSurfaceContainer which defines the digital surface as connected surfels....
Definition: SetOfSurfels.h:73
DGtal::DUAL
@ DUAL
Definition: Duality.h:62
Space
SpaceND< 2 > Space
Definition: testSimpleRandomAccessRangeFromPoint.cpp:42
DGtal::KhalimskySpaceND::unsigns
Cell unsigns(const SCell &p) const
Creates an unsigned cell from a signed one.
DGtal::CanonicSCellEmbedder< KSpace >
Point
MyPointD Point
Definition: testClone2.cpp:383
RealPoint
Z2i::RealPoint RealPoint
Definition: testAstroid2D.cpp:46
DGtal::Ball3D
Aim: Model of the concept StarShaped3D represents any Sphere in the space.
Definition: Ball3D.h:60
DGtal::Astroid2D::getUpperBound
RealPoint getUpperBound() const
Definition: Astroid2D.h:131
DGtal::HyperRectDomain::lowerBound
const Point & lowerBound() const
MyDigitalSurface
DigitalSurface< MyDigitalSurfaceContainer > MyDigitalSurface
Definition: greedy-plane-segmentation-ex2.cpp:92
SurfelSet
MyDigitalSurface::SurfelSet SurfelSet
Definition: greedy-plane-segmentation-ex2.cpp:95
DGtal::KhalimskySpaceND
Aim: This class is a model of CCellularGridSpaceND. It represents the cubical grid as a cell complex,...
Definition: KhalimskySpaceND.h:64