23#include "DGtal/io/readers/VolReader.h"
24#include "DGtal/io/viewers/PolyscopeViewer.h"
25#include "DGtal/topology/SurfelAdjacency.h"
26#include "DGtal/topology/LightImplicitDigitalSurface.h"
27#include "DGtal/topology/DigitalSurface.h"
28#include "DGtal/math/linalg/EigenSupport.h"
29#include "DGtal/dec/DiscreteExteriorCalculus.h"
30#include "DGtal/dec/DiscreteExteriorCalculusSolver.h"
31#include "DGtal/dec/DiscreteExteriorCalculusFactory.h"
33#include "ConfigExamples.h"
37template <
typename Predicate,
typename Domain>
43 typedef typename Predicate::Point
Point;
53 if (!
myDomain->isInside(point))
return false;
69 const std::string filename = examplesPath +
"samples/Al.100.vol";
90 const SurfelAdjacency surfel_adjacency(
true);
93 const DigitalSurfaceContainer digital_surface_container(kspace, image_extended, surfel_adjacency, cell_bel);
96 const DigitalSurface digital_surface(digital_surface_container);
99 trace.
info() <<
"digital_surface_size=" << digital_surface.size() << endl;
102 Viewer* viewer =
new Viewer(kspace);
104 for (DigitalSurface::ConstIterator si=digital_surface.begin(), se=digital_surface.end(); si!=se; si++)
120 const Calculus
calculus = CalculusFactory::createFromNSCells<2>(digital_surface.begin(), digital_surface.end(),
false);
125 Viewer* viewer =
new Viewer(kspace);
142 const Calculus::SCell cell = rho.getSCell(
index);
143 const Calculus::Point point = kspace.
sKCoords(cell);
144 if (point[2]>1)
continue;
145 rho.myContainer(
index) = point[1]>100 ? 1 : -1;
147 trace.
info() <<
"rho_range=" << rho.myContainer.minCoeff() <<
"/" << rho.myContainer.maxCoeff() << endl;
151 Viewer* viewer =
new Viewer(kspace);
158 const Calculus::DualIdentity0 laplace =
calculus.laplace<
DUAL>();
162 PoissonSolver solver;
164 const Calculus::DualForm0
phi = solver.solve(rho);
165 trace.
info() <<
"phi_range=" <<
phi.myContainer.minCoeff() <<
"/" <<
phi.myContainer.maxCoeff() << endl;
169 Viewer* viewer =
new Viewer(kspace);
191 trace.
info() <<
"input_domain=" << input_domain << endl;
194 kspace.
init(input_domain.lowerBound(), input_domain.upperBound(),
true);
204 trace.
info() <<
"input_set_size=" << input_set.size() << endl;
207 Viewer* viewer =
new Viewer(kspace);
208 (*viewer) << input_set;
217 const SurfelAdjacency surfel_adjacency(
true);
220 const DigitalSurfaceContainer digital_surface_container(kspace, input_set, surfel_adjacency, cell_bel);
223 const DigitalSurface digital_surface(digital_surface_container);
226 trace.
info() <<
"surfel_adjacency=" << endl;
227 for (
int ii=0; ii<3; ii++)
229 std::stringstream ss;
230 for (
int jj=0; jj<3; jj++)
231 ss << surfel_adjacency.getAdjacency(ii, jj) <<
" ";
235 trace.
info() <<
"digital_surface_container=" << digital_surface_container << endl;
237 trace.
info() <<
"digital_surface_size=" << digital_surface.size() << endl;
240 Viewer* viewer =
new Viewer(kspace);
241 for (DigitalSurface::ConstIterator si=digital_surface.begin(), se=digital_surface.end(); si!=se; si++)
257 const Calculus
calculus = CalculusFactory::createFromNSCells<2>(digital_surface.begin(), digital_surface.end());
262 Viewer* viewer =
new Viewer(kspace);
271 const Calculus::PrimalForm2 primal_surfel_area =
calculus.hodge<0,
DUAL>() *
272 Calculus::DualForm0(
calculus, Eigen::VectorXd::Ones(
calculus.kFormLength(0, DUAL)));
274 const Calculus::DualForm2 dual_surfel_area =
calculus.hodge<0,
PRIMAL>() *
275 Calculus::PrimalForm0(
calculus, Eigen::VectorXd::Ones(
calculus.kFormLength(0, PRIMAL)));
277 const double area_th =
calculus.kFormLength(2, PRIMAL);
278 const double area_primal = dual_surfel_area.myContainer.sum();
279 const double area_dual = primal_surfel_area.myContainer.sum();
280 trace.
info() <<
"area_th=" << area_th << endl;
281 trace.
info() <<
"area_primal=" << area_primal << endl;
282 trace.
info() <<
"area_dual=" << area_dual << endl;
283 FATAL_ERROR( area_th == area_primal );
284 FATAL_ERROR( area_th == area_dual );
286 const Calculus::DualForm1 dual_edge_length =
calculus.hodge<1,
PRIMAL>() *
287 Calculus::PrimalForm1(
calculus, Eigen::VectorXd::Ones(
calculus.kFormLength(1, PRIMAL)));
289 const bool dual_edge_length_match = dual_edge_length.myContainer == Eigen::VectorXd::Ones(
calculus.kFormLength(1, DUAL)) ;
290 trace.
info() <<
"dual_edge_length_match=" << dual_edge_length_match << endl;
291 FATAL_ERROR( dual_edge_length_match );
295 for (Calculus::Index
index=0;
index<dual_surfel_area.length();
index++)
297 const Calculus::SCell cell = dual_surfel_area.getSCell(
index);
298 const Calculus::Scalar area = dual_surfel_area.myContainer(
index);
309int main(
int argc,
char* argv[])
Aim: This class encapsulates its parameter class so that to indicate to the user that the object/poin...
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: This wraps a linear algebra solver around a discrete exterior calculus.
DiscreteExteriorCalculusSolver & compute(const Operator &linear_operator)
Aim: DiscreteExteriorCalculus represents a calculus in the dec package. This is the main structure in...
bool init(const Point &lower, const Point &upper, bool isClosed)
Specifies the upper and lower bounds for the maximal cells in this space.
SignedKhalimskyCell< dim, Integer > SCell
const Point & sKCoords(const SCell &c) const
Return its Khalimsky coordinates.
Aim: A model of CDigitalSurfaceContainer which defines the digital surface as the boundary of an impl...
static Self diagonal(Component val=1)
void show() override
Starts the event loop and display of elements.
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="")
PolygonalCalculus< SH3::RealPoint, SH3::RealVector >::Vector phi(const Face f)
HyperRectDomain< Space > Domain
KhalimskySpaceND< 3, Integer > KSpace
DigitalSetSelector< Domain, BIG_DS+HIGH_BEL_DS >::Type DigitalSet
Eigen::SparseLU< SparseMatrix > SolverSparseLU
static ImageContainer importVol(const std::string &filename, const Functor &aFunctor=Functor())
Aim: This concept represents a digital domain, i.e. a non mutable subset of points of the given digit...
Aim: Defines a predicate on a point.
Aim: Define a simple functor using the static cast operator.
FalseOutsideDomain(DGtal::ConstAlias< Predicate > predicate, DGtal::ConstAlias< Domain > adomain)
BOOST_CONCEPT_ASSERT((DGtal::concepts::CDomain< Domain >))
const Predicate * myPredicate
BOOST_CONCEPT_ASSERT((DGtal::concepts::CPointPredicate< Predicate >))
bool operator()(const Point &point) const
unsigned int index(DGtal::uint32_t n, unsigned int b)
InHalfPlaneBySimple3x3Matrix< Point, double > Functor
ImageContainerBySTLVector< Domain, Value > Image
HyperRectDomain< Space > Domain