8 #include "DGtal/io/readers/VolReader.h"
9 #include "DGtal/io/viewers/Viewer3D.h"
10 #include "DGtal/topology/SurfelAdjacency.h"
11 #include "DGtal/topology/LightImplicitDigitalSurface.h"
12 #include "DGtal/topology/DigitalSurface.h"
13 #include "DGtal/math/linalg/EigenSupport.h"
14 #include "DGtal/dec/DiscreteExteriorCalculus.h"
15 #include "DGtal/dec/DiscreteExteriorCalculusSolver.h"
16 #include "DGtal/dec/DiscreteExteriorCalculusFactory.h"
18 #include "ConfigExamples.h"
23 template <
typename Predicate,
typename Domain>
24 struct FalseOutsideDomain
29 typedef typename Predicate::Point Point;
32 myPredicate(&predicate), myDomain(&adomain)
37 operator()(
const Point& point)
const
39 if (!myDomain->isInside(point))
return false;
40 return (*myPredicate)(point);
43 const Predicate* myPredicate;
44 const Domain* myDomain;
55 const std::string filename = examplesPath +
"samples/Al.100.vol";
65 trace.
info() <<
"domain=" << domain << endl;
67 typedef FalseOutsideDomain<Image, DGtal::Z3i::Domain> ImageExtended;
68 const ImageExtended image_extended(image, domain);
76 const SurfelAdjacency surfel_adjacency(
true);
79 const DigitalSurfaceContainer digital_surface_container(kspace, image_extended, surfel_adjacency, cell_bel);
82 const DigitalSurface digital_surface(digital_surface_container);
85 trace.
info() <<
"digital_surface_size=" << digital_surface.size() << endl;
90 viewer->setWindowTitle(
"alcapone surface");
91 for (DigitalSurface::ConstIterator si=digital_surface.begin(), se=digital_surface.end(); si!=se; si++)
106 const Calculus calculus = CalculusFactory::createFromNSCells<2>(digital_surface.begin(), digital_surface.end(),
false);
108 trace.
info() <<
"calculus=" << calculus << endl;
113 viewer->setWindowTitle(
"alcapone calculus");
126 Calculus::DualForm0 rho(calculus);
127 for (
int index=0; index<rho.length(); index++)
129 const Calculus::SCell cell = rho.getSCell(index);
130 const Calculus::Point point = kspace.
sKCoords(cell);
131 if (point[2]>1)
continue;
132 rho.myContainer(index) = point[1]>100 ? 1 : -1;
134 trace.
info() <<
"rho_range=" << rho.myContainer.minCoeff() <<
"/" << rho.myContainer.maxCoeff() << endl;
140 viewer->setWindowTitle(
"alcapone poisson rho");
146 const Calculus::DualIdentity0 laplace = calculus.laplace<
DUAL>();
150 PoissonSolver solver;
152 const Calculus::DualForm0 phi = solver.
solve(rho);
153 trace.
info() <<
"phi_range=" << phi.myContainer.minCoeff() <<
"/" << phi.myContainer.maxCoeff() << endl;
159 viewer->setWindowTitle(
"alcapone poisson phi");
180 trace.
info() <<
"input_domain=" << input_domain << endl;
183 kspace.
init(input_domain.lowerBound(), input_domain.upperBound(),
true);
193 trace.
info() <<
"input_set_size=" << input_set.size() << endl;
198 viewer->setWindowTitle(
"input set");
199 (*viewer) << input_set;
207 const SurfelAdjacency surfel_adjacency(
true);
210 const DigitalSurfaceContainer digital_surface_container(kspace, input_set, surfel_adjacency, cell_bel);
213 const DigitalSurface digital_surface(digital_surface_container);
216 trace.
info() <<
"surfel_adjacency=" << endl;
217 for (
int ii=0; ii<3; ii++)
219 std::stringstream ss;
220 for (
int jj=0; jj<3; jj++)
221 ss << surfel_adjacency.getAdjacency(ii, jj) <<
" ";
225 trace.
info() <<
"digital_surface_container=" << digital_surface_container << endl;
227 trace.
info() <<
"digital_surface_size=" << digital_surface.size() << endl;
232 viewer->setWindowTitle(
"digital surface");
233 for (DigitalSurface::ConstIterator si=digital_surface.begin(), se=digital_surface.end(); si!=se; si++)
248 const Calculus calculus = CalculusFactory::createFromNSCells<2>(digital_surface.begin(), digital_surface.end());
250 trace.
info() <<
"calculus=" << calculus << endl;
255 viewer->setWindowTitle(
"discrete exterior calculus");
263 const Calculus::PrimalForm2 primal_surfel_area = calculus.hodge<0,
DUAL>() *
264 Calculus::DualForm0(calculus, Eigen::VectorXd::Ones(calculus.kFormLength(0, DUAL)));
266 const Calculus::DualForm2 dual_surfel_area = calculus.hodge<0,
PRIMAL>() *
267 Calculus::PrimalForm0(calculus, Eigen::VectorXd::Ones(calculus.kFormLength(0, PRIMAL)));
269 const double area_th = calculus.kFormLength(2, PRIMAL);
270 const double area_primal = dual_surfel_area.myContainer.sum();
271 const double area_dual = primal_surfel_area.myContainer.sum();
272 trace.
info() <<
"area_th=" << area_th << endl;
273 trace.
info() <<
"area_primal=" << area_primal << endl;
274 trace.
info() <<
"area_dual=" << area_dual << endl;
275 FATAL_ERROR( area_th == area_primal );
276 FATAL_ERROR( area_th == area_dual );
278 const Calculus::DualForm1 dual_edge_length = calculus.hodge<1,
PRIMAL>() *
279 Calculus::PrimalForm1(calculus, Eigen::VectorXd::Ones(calculus.kFormLength(1, PRIMAL)));
281 const bool dual_edge_length_match = dual_edge_length.myContainer == Eigen::VectorXd::Ones(calculus.kFormLength(1, DUAL)) ;
282 trace.
info() <<
"dual_edge_length_match=" << dual_edge_length_match << endl;
283 FATAL_ERROR( dual_edge_length_match );
287 for (Calculus::Index index=0; index<dual_surfel_area.length(); index++)
289 const Calculus::SCell cell = dual_surfel_area.getSCell(index);
290 const Calculus::Scalar area = dual_surfel_area.myContainer(index);
291 trace.
info() << index <<
" " << cell <<
" " << area << endl;
301 int main(
int argc,
char* argv[])
303 QApplication app(argc, argv);
void beginBlock(const std::string &keyword="")
Aim: DiscreteExteriorCalculus represents a calculus in the dec package. This is the main structure in...
Aim: A model of CDigitalSurfaceContainer which defines the digital surface as the boundary of an impl...
static Self diagonal(Component val=1)
SolutionKForm solve(const InputKForm &input_kform) const
const Point & sKCoords(const SCell &c) const
Aim: This class encapsulates its parameter class so that to indicate to the user that the object/poin...
Aim: Define a simple functor using the static cast operator.
Aim: Represent adjacencies between surfel elements, telling if it follows an interior to exterior ord...
Aim: Represents a set of n-1-cells in a nD space, together with adjacency relation between these cell...
Factory for GPL Display3D:
static void draw(Display3D< Space, KSpace > &display, const DGtal::DiscreteExteriorCalculus< dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger > &calculus)
static SCell findABel(const KSpace &K, const PointPredicate &pp, unsigned int nbtries=1000)
DiscreteExteriorCalculusSolver & compute(const Operator &linear_operator)
Aim: This concept represents a digital domain, i.e. a non mutable subset of points of the given digit...
static ImageContainer importVol(const std::string &filename, const Functor &aFunctor=Functor())
Aim: This class provides static members to create DEC structures from various other DGtal structures...
Aim: Defines a predicate on a point.
bool init(const Point &lower, const Point &upper, bool isClosed)
const Point & upperBound() const
Aim: A wrapper class around a STL associative container for storing sets of digital points within som...
Represents a signed cell in a cellular grid space by its Khalimsky coordinates and a boolean value...
Aim: This wraps a linear algebra solver around a discrete exterior calculus.
virtual void show()
Overload QWidget method in order to add a call to updateList() method (to ensure that the lists are w...
const Point & lowerBound() const
Aim: This class is a model of CCellularGridSpaceND. It represents the cubical grid as a cell complex...
Eigen::SparseLU< SparseMatrix > SolverSparseLU