Loading [MathJax]/extensions/TeX/AMSsymbols.js
DGtal 2.0.0
exampleDECSurface.cpp
1
17
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"
32
33#include "ConfigExamples.h"
34
36
37template <typename Predicate, typename Domain>
39{
42
43 typedef typename Predicate::Point Point;
44
46 myPredicate(&predicate), myDomain(&adomain)
47 {
48 }
49
50 bool
51 operator()(const Point& point) const
52 {
53 if (!myDomain->isInside(point)) return false;
54 return (*myPredicate)(point);
55 }
56
57 const Predicate* myPredicate;
59};
60
61void
62alcapone_3d()
63{
64 using std::endl;
65 using DGtal::trace;
66
67 trace.beginBlock("alcapone");
68
69 const std::string filename = examplesPath + "samples/Al.100.vol";
70
71 trace.beginBlock("digital surface");
72
77 const DGtal::Z3i::Domain domain = image.domain();
78
79 trace.info() << "domain=" << domain << endl;
80
82 const ImageExtended image_extended(image, domain);
83
84 DGtal::Z3i::KSpace kspace;
85 kspace.init(domain.lowerBound(), domain.upperBound()+DGtal::Z3i::Point(0,0,1), true);
86
87 const DGtal::Z3i::KSpace::SCell cell_bel = DGtal::Surfaces<DGtal::Z3i::KSpace>::findABel(kspace, image_extended);
88
89 typedef DGtal::SurfelAdjacency<3> SurfelAdjacency;
90 const SurfelAdjacency surfel_adjacency(true);
91
93 const DigitalSurfaceContainer digital_surface_container(kspace, image_extended, surfel_adjacency, cell_bel);
94
96 const DigitalSurface digital_surface(digital_surface_container);
98
99 trace.info() << "digital_surface_size=" << digital_surface.size() << endl;
100
101 {
102 Viewer* viewer = new Viewer(kspace);
103 viewer->allowReuseList = true;
104 for (DigitalSurface::ConstIterator si=digital_surface.begin(), se=digital_surface.end(); si!=se; si++)
105 {
106 const DGtal::Z3i::KSpace::SCell cell = *si;
107 (*viewer) << cell;
108 }
109 viewer->show();
110 delete viewer;
111 }
112
113 trace.endBlock();
114
115 trace.beginBlock("discrete exterior calculus");
116
120 const Calculus calculus = CalculusFactory::createFromNSCells<2>(digital_surface.begin(), digital_surface.end(), false);
122 trace.info() << "calculus=" << calculus << endl;
123
124 {
125 Viewer* viewer = new Viewer(kspace);
126 (*viewer) << calculus;
127 viewer->show();
128 delete viewer;
129 }
130
131 using DGtal::PRIMAL;
132 using DGtal::DUAL;
133
134 trace.endBlock();
135
136 trace.beginBlock("poisson equation");
137
139 Calculus::DualForm0 rho(calculus);
140 for (int index=0; index<rho.length(); index++)
141 {
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;
146 }
147 trace.info() << "rho_range=" << rho.myContainer.minCoeff() << "/" << rho.myContainer.maxCoeff() << endl;
149
150 {
151 Viewer* viewer = new Viewer(kspace);
152 (*viewer) << rho;
153 viewer->show();
154 delete viewer;
155 }
156
158 const Calculus::DualIdentity0 laplace = calculus.laplace<DUAL>();
159
160 typedef DGtal::EigenLinearAlgebraBackend::SolverSparseLU LinearAlgebraSolver;
162 PoissonSolver solver;
163 solver.compute(laplace);
164 const Calculus::DualForm0 phi = solver.solve(rho);
165 trace.info() << "phi_range=" << phi.myContainer.minCoeff() << "/" << phi.myContainer.maxCoeff() << endl;
167
168 {
169 Viewer* viewer = new Viewer(kspace);
170 (*viewer) << phi;
171 viewer->show();
172 delete viewer;
173 }
174
175 trace.endBlock();
176
177 trace.endBlock();
178}
179
180void
181pyramid_3d()
182{
183 using std::endl;
184 using DGtal::trace;
185
186 trace.beginBlock("pyramid example");
187
188 trace.beginBlock("digital surface");
189
191 trace.info() << "input_domain=" << input_domain << endl;
192
193 DGtal::Z3i::KSpace kspace;
194 kspace.init(input_domain.lowerBound(), input_domain.upperBound(), true);
195
197 DGtal::Z3i::DigitalSet input_set(input_domain);
198 input_set.insert( DGtal::Z3i::Point(0,0,0) );
199 input_set.insert( DGtal::Z3i::Point(1,0,0) );
200 input_set.insert( DGtal::Z3i::Point(1,1,0) );
201 input_set.insert( DGtal::Z3i::Point(0,1,0) );
202 input_set.insert( DGtal::Z3i::Point(0,0,1) );
204 trace.info() << "input_set_size=" << input_set.size() << endl;
205
206 {
207 Viewer* viewer = new Viewer(kspace);
208 (*viewer) << input_set;
209 viewer->show();
210 delete viewer;
211 }
212
215
216 typedef DGtal::SurfelAdjacency<3> SurfelAdjacency;
217 const SurfelAdjacency surfel_adjacency(true);
218
220 const DigitalSurfaceContainer digital_surface_container(kspace, input_set, surfel_adjacency, cell_bel);
221
223 const DigitalSurface digital_surface(digital_surface_container);
225
226 trace.info() << "surfel_adjacency=" << endl;
227 for (int ii=0; ii<3; ii++)
228 {
229 std::stringstream ss;
230 for (int jj=0; jj<3; jj++)
231 ss << surfel_adjacency.getAdjacency(ii, jj) << " ";
232 trace.info() << ss.str() << endl;
233 }
234
235 trace.info() << "digital_surface_container=" << digital_surface_container << endl;
236
237 trace.info() << "digital_surface_size=" << digital_surface.size() << endl;
238
239 {
240 Viewer* viewer = new Viewer(kspace);
241 for (DigitalSurface::ConstIterator si=digital_surface.begin(), se=digital_surface.end(); si!=se; si++)
242 {
243 const DGtal::Z3i::KSpace::SCell cell = *si;
244 (*viewer) << cell;
245 }
246 viewer->show();
247 delete viewer;
248 }
249
250 trace.endBlock();
251
252 trace.beginBlock("discrete exterior calculus");
253
257 const Calculus calculus = CalculusFactory::createFromNSCells<2>(digital_surface.begin(), digital_surface.end());
259 trace.info() << "calculus=" << calculus << endl;
260
261 {
262 Viewer* viewer = new Viewer(kspace);
263 (*viewer) << calculus;
264 viewer->show();
265 delete viewer;
266 }
267
268 using DGtal::PRIMAL;
269 using DGtal::DUAL;
270
271 const Calculus::PrimalForm2 primal_surfel_area = calculus.hodge<0, DUAL>() *
272 Calculus::DualForm0(calculus, Eigen::VectorXd::Ones(calculus.kFormLength(0, DUAL)));
273
274 const Calculus::DualForm2 dual_surfel_area = calculus.hodge<0, PRIMAL>() *
275 Calculus::PrimalForm0(calculus, Eigen::VectorXd::Ones(calculus.kFormLength(0, PRIMAL)));
276
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 );
285
286 const Calculus::DualForm1 dual_edge_length = calculus.hodge<1, PRIMAL>() *
287 Calculus::PrimalForm1(calculus, Eigen::VectorXd::Ones(calculus.kFormLength(1, PRIMAL)));
288
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 );
292
293 trace.beginBlock("dual surfel area");
294
295 for (Calculus::Index index=0; index<dual_surfel_area.length(); index++)
296 {
297 const Calculus::SCell cell = dual_surfel_area.getSCell(index);
298 const Calculus::Scalar area = dual_surfel_area.myContainer(index);
299 trace.info() << index << " " << cell << " " << area << endl;
300 }
301
302 trace.endBlock();
303
304 trace.endBlock();
305
306 trace.endBlock();
307}
308
309int main(int argc, char* argv[])
310{
311 //pyramid_3d();
312 alcapone_3d();
313
314 return 0;
315}
Aim: This class encapsulates its parameter class so that to indicate to the user that the object/poin...
Definition ConstAlias.h:187
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="")
std::ostream & info()
double endBlock()
PolygonalCalculus< SH3::RealPoint, SH3::RealVector >::Vector phi(const Face f)
PolyCalculus * calculus
HyperRectDomain< Space > Domain
Definition StdDefs.h:172
KhalimskySpaceND< 3, Integer > KSpace
Definition StdDefs.h:146
Space::Point Point
Definition StdDefs.h:168
DigitalSetSelector< Domain, BIG_DS+HIGH_BEL_DS >::Type DigitalSet
Definition StdDefs.h:173
Trace trace
@ PRIMAL
Definition Duality.h:61
@ DUAL
Definition Duality.h:62
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...
Definition CDomain.h:130
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 >))
Predicate::Point Point
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)
Definition testBits.cpp:44
int main()
Definition testBits.cpp:56
InHalfPlaneBySimple3x3Matrix< Point, double > Functor
Domain domain
Image image(domain)
ImageContainerBySTLVector< Domain, Value > Image
HyperRectDomain< Space > Domain