DGtal 1.3.0
Loading...
Searching...
No Matches
exampleDECSurface.cpp
1
17
23#include "DGtal/io/readers/VolReader.h"
24#include "DGtal/io/viewers/Viewer3D.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
37
38template <typename Predicate, typename Domain>
39struct FalseOutsideDomain
40{
41 BOOST_CONCEPT_ASSERT(( DGtal::concepts::CPointPredicate<Predicate> ));
42 BOOST_CONCEPT_ASSERT(( DGtal::concepts::CDomain<Domain> ));
43
44 typedef typename Predicate::Point Point;
45
46 FalseOutsideDomain(DGtal::ConstAlias<Predicate> predicate, DGtal::ConstAlias<Domain> adomain) :
47 myPredicate(&predicate), myDomain(&adomain)
48 {
49 }
50
51 bool
52 operator()(const Point& point) const
53 {
54 if (!myDomain->isInside(point)) return false;
55 return (*myPredicate)(point);
56 }
57
58 const Predicate* myPredicate;
59 const Domain* myDomain;
60};
61
62void
63alcapone_3d()
64{
65 using std::endl;
66 using DGtal::trace;
67
68 trace.beginBlock("alcapone");
69
70 const std::string filename = examplesPath + "samples/Al.100.vol";
71
72 trace.beginBlock("digital surface");
73
78 const DGtal::Z3i::Domain domain = image.domain();
79
80 trace.info() << "domain=" << domain << endl;
81
82 typedef FalseOutsideDomain<Image, DGtal::Z3i::Domain> ImageExtended;
83 const ImageExtended image_extended(image, domain);
84
85 DGtal::Z3i::KSpace kspace;
86 kspace.init(domain.lowerBound(), domain.upperBound()+DGtal::Z3i::Point(0,0,1), true);
87
88 const DGtal::Z3i::KSpace::SCell cell_bel = DGtal::Surfaces<DGtal::Z3i::KSpace>::findABel(kspace, image_extended);
89
90 typedef DGtal::SurfelAdjacency<3> SurfelAdjacency;
91 const SurfelAdjacency surfel_adjacency(true);
92
94 const DigitalSurfaceContainer digital_surface_container(kspace, image_extended, surfel_adjacency, cell_bel);
95
97 const DigitalSurface digital_surface(digital_surface_container);
99
100 trace.info() << "digital_surface_size=" << digital_surface.size() << endl;
101
102 {
103 Viewer* viewer = new Viewer(kspace);
104 viewer->show();
105 viewer->setWindowTitle("alcapone surface");
106 for (DigitalSurface::ConstIterator si=digital_surface.begin(), se=digital_surface.end(); si!=se; si++)
107 {
108 const DGtal::Z3i::KSpace::SCell cell = *si;
109 (*viewer) << cell;
110 }
111 (*viewer) << Viewer::updateDisplay;
112 }
113
114 trace.endBlock();
115
116 trace.beginBlock("discrete exterior calculus");
117
121 const Calculus calculus = CalculusFactory::createFromNSCells<2>(digital_surface.begin(), digital_surface.end(), false);
123 trace.info() << "calculus=" << calculus << endl;
124
125 {
126 Viewer* viewer = new Viewer(kspace);
127 viewer->show();
128 viewer->setWindowTitle("alcapone calculus");
129 DisplayFactory::draw(*viewer, calculus);
130 (*viewer) << Viewer::updateDisplay;
131 }
132
133 using DGtal::PRIMAL;
134 using DGtal::DUAL;
135
136 trace.endBlock();
137
138 trace.beginBlock("poisson equation");
139
141 Calculus::DualForm0 rho(calculus);
142 for (int index=0; index<rho.length(); index++)
143 {
144 const Calculus::SCell cell = rho.getSCell(index);
145 const Calculus::Point point = kspace.sKCoords(cell);
146 if (point[2]>1) continue;
147 rho.myContainer(index) = point[1]>100 ? 1 : -1;
148 }
149 trace.info() << "rho_range=" << rho.myContainer.minCoeff() << "/" << rho.myContainer.maxCoeff() << endl;
151
152 {
153 Viewer* viewer = new Viewer(kspace);
154 viewer->show();
155 viewer->setWindowTitle("alcapone poisson rho");
156 DisplayFactory::draw(*viewer, rho);
157 (*viewer) << Viewer::updateDisplay;
158 }
159
161 const Calculus::DualIdentity0 laplace = calculus.laplace<DUAL>();
162
163 typedef DGtal::EigenLinearAlgebraBackend::SolverSparseLU LinearAlgebraSolver;
165 PoissonSolver solver;
166 solver.compute(laplace);
167 const Calculus::DualForm0 phi = solver.solve(rho);
168 trace.info() << "phi_range=" << phi.myContainer.minCoeff() << "/" << phi.myContainer.maxCoeff() << endl;
170
171 {
172 Viewer* viewer = new Viewer(kspace);
173 viewer->show();
174 viewer->setWindowTitle("alcapone poisson phi");
175 DisplayFactory::draw(*viewer, phi);
176 (*viewer) << Viewer::updateDisplay;
177 }
178
179 trace.endBlock();
180
181 trace.endBlock();
182}
183
184void
185pyramid_3d()
186{
187 using std::endl;
188 using DGtal::trace;
189
190 trace.beginBlock("pyramid example");
191
192 trace.beginBlock("digital surface");
193
195 trace.info() << "input_domain=" << input_domain << endl;
196
197 DGtal::Z3i::KSpace kspace;
198 kspace.init(input_domain.lowerBound(), input_domain.upperBound(), true);
199
201 DGtal::Z3i::DigitalSet input_set(input_domain);
202 input_set.insert( DGtal::Z3i::Point(0,0,0) );
203 input_set.insert( DGtal::Z3i::Point(1,0,0) );
204 input_set.insert( DGtal::Z3i::Point(1,1,0) );
205 input_set.insert( DGtal::Z3i::Point(0,1,0) );
206 input_set.insert( DGtal::Z3i::Point(0,0,1) );
208 trace.info() << "input_set_size=" << input_set.size() << endl;
209
210 {
211 Viewer* viewer = new Viewer(kspace);
212 viewer->show();
213 viewer->setWindowTitle("input set");
214 (*viewer) << input_set;
215 (*viewer) << Viewer::updateDisplay;
216 }
217
220
221 typedef DGtal::SurfelAdjacency<3> SurfelAdjacency;
222 const SurfelAdjacency surfel_adjacency(true);
223
225 const DigitalSurfaceContainer digital_surface_container(kspace, input_set, surfel_adjacency, cell_bel);
226
228 const DigitalSurface digital_surface(digital_surface_container);
230
231 trace.info() << "surfel_adjacency=" << endl;
232 for (int ii=0; ii<3; ii++)
233 {
234 std::stringstream ss;
235 for (int jj=0; jj<3; jj++)
236 ss << surfel_adjacency.getAdjacency(ii, jj) << " ";
237 trace.info() << ss.str() << endl;
238 }
239
240 trace.info() << "digital_surface_container=" << digital_surface_container << endl;
241
242 trace.info() << "digital_surface_size=" << digital_surface.size() << endl;
243
244 {
245 Viewer* viewer = new Viewer(kspace);
246 viewer->show();
247 viewer->setWindowTitle("digital surface");
248 for (DigitalSurface::ConstIterator si=digital_surface.begin(), se=digital_surface.end(); si!=se; si++)
249 {
250 const DGtal::Z3i::KSpace::SCell cell = *si;
251 (*viewer) << cell;
252 }
253 (*viewer) << Viewer::updateDisplay;
254 }
255
256 trace.endBlock();
257
258 trace.beginBlock("discrete exterior calculus");
259
263 const Calculus calculus = CalculusFactory::createFromNSCells<2>(digital_surface.begin(), digital_surface.end());
265 trace.info() << "calculus=" << calculus << endl;
266
267 {
268 Viewer* viewer = new Viewer(kspace);
269 viewer->show();
270 viewer->setWindowTitle("discrete exterior calculus");
271 DisplayFactory::draw(*viewer, calculus);
272 (*viewer) << Viewer::updateDisplay;
273 }
274
275 using DGtal::PRIMAL;
276 using DGtal::DUAL;
277
278 const Calculus::PrimalForm2 primal_surfel_area = calculus.hodge<0, DUAL>() *
279 Calculus::DualForm0(calculus, Eigen::VectorXd::Ones(calculus.kFormLength(0, DUAL)));
280
281 const Calculus::DualForm2 dual_surfel_area = calculus.hodge<0, PRIMAL>() *
282 Calculus::PrimalForm0(calculus, Eigen::VectorXd::Ones(calculus.kFormLength(0, PRIMAL)));
283
284 const double area_th = calculus.kFormLength(2, PRIMAL);
285 const double area_primal = dual_surfel_area.myContainer.sum();
286 const double area_dual = primal_surfel_area.myContainer.sum();
287 trace.info() << "area_th=" << area_th << endl;
288 trace.info() << "area_primal=" << area_primal << endl;
289 trace.info() << "area_dual=" << area_dual << endl;
290 FATAL_ERROR( area_th == area_primal );
291 FATAL_ERROR( area_th == area_dual );
292
293 const Calculus::DualForm1 dual_edge_length = calculus.hodge<1, PRIMAL>() *
294 Calculus::PrimalForm1(calculus, Eigen::VectorXd::Ones(calculus.kFormLength(1, PRIMAL)));
295
296 const bool dual_edge_length_match = dual_edge_length.myContainer == Eigen::VectorXd::Ones(calculus.kFormLength(1, DUAL)) ;
297 trace.info() << "dual_edge_length_match=" << dual_edge_length_match << endl;
298 FATAL_ERROR( dual_edge_length_match );
299
300 trace.beginBlock("dual surfel area");
301
302 for (Calculus::Index index=0; index<dual_surfel_area.length(); index++)
303 {
304 const Calculus::SCell cell = dual_surfel_area.getSCell(index);
305 const Calculus::Scalar area = dual_surfel_area.myContainer(index);
306 trace.info() << index << " " << cell << " " << area << endl;
307 }
308
309 trace.endBlock();
310
311 trace.endBlock();
312
313 trace.endBlock();
314}
315
316int main(int argc, char* argv[])
317{
318 QApplication app(argc, argv);
319
320 //pyramid_3d();
321 alcapone_3d();
322
323 return app.exec();
324}
Aim: This class encapsulates its parameter class so that to indicate to the user that the object/poin...
Definition: ConstAlias.h:187
Aim: A wrapper class around a STL associative container for storing sets of digital points within som...
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...
Aim: This class is a model of CCellularGridSpaceND. It represents the cubical grid as a cell complex,...
bool init(const Point &lower, const Point &upper, bool isClosed)
Specifies the upper and lower bounds for the maximal cells in this space.
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)
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()
virtual void show()
Overload QWidget method in order to add a call to updateList() method (to ensure that the lists are w...
Trace trace
Definition: Common.h:154
@ PRIMAL
Definition: Duality.h:61
@ DUAL
Definition: Duality.h:62
Factory for GPL Display3D:
static void draw(Display3D< Space, KSpace > &display, const DGtal::DiscreteExteriorCalculus< dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger > &calculus)
Eigen::SparseLU< SparseMatrix > SolverSparseLU
Definition: EigenSupport.h:110
Represents a signed cell in a cellular grid space by its Khalimsky coordinates and a boolean value.
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.
int main()
Definition: testBits.cpp:56
MyPointD Point
Definition: testClone2.cpp:383
InHalfPlaneBySimple3x3Matrix< Point, double > Functor
Domain domain
ImageContainerBySTLVector< Domain, Value > Image
HyperRectDomain< Space > Domain