DGtal  1.1.0
exampleDECSurface.cpp
1 
16 
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 
35 typedef DGtal::Viewer3D<DGtal::Z3i::Space, DGtal::Z3i::KSpace> Viewer;
36 typedef DGtal::Display3DFactory<DGtal::Z3i::Space, DGtal::Z3i::KSpace> DisplayFactory;
37 
38 template <typename Predicate, typename Domain>
39 struct 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 
62 void
63 alcapone_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 
75  typedef DGtal::ImageContainerBySTLVector<DGtal::Z3i::Domain, bool> Image;
76  typedef DGtal::functors::Cast<bool> Functor;
77  const Image image = DGtal::VolReader<Image, Functor>::importVol(filename, Functor());
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 
93  typedef DGtal::LightImplicitDigitalSurface<DGtal::Z3i::KSpace, ImageExtended> DigitalSurfaceContainer;
94  const DigitalSurfaceContainer digital_surface_container(kspace, image_extended, surfel_adjacency, cell_bel);
95 
96  typedef DGtal::DigitalSurface<DigitalSurfaceContainer> DigitalSurface;
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 
119  typedef DGtal::DiscreteExteriorCalculusFactory<DGtal::EigenLinearAlgebraBackend> CalculusFactory;
120  typedef DGtal::DiscreteExteriorCalculus<2, 3, DGtal::EigenLinearAlgebraBackend> Calculus;
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;
164  typedef DGtal::DiscreteExteriorCalculusSolver<Calculus, LinearAlgebraSolver, 0, DUAL, 0, DUAL> PoissonSolver;
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 
184 void
185 pyramid_3d()
186 {
187  using std::endl;
188  using DGtal::trace;
189 
190  trace.beginBlock("pyramid example");
191 
192  trace.beginBlock("digital surface");
193 
194  const DGtal::Z3i::Domain input_domain(DGtal::Z3i::Point::diagonal(-1), DGtal::Z3i::Point::diagonal(2));
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 
219  const DGtal::Z3i::KSpace::SCell cell_bel = DGtal::Surfaces<DGtal::Z3i::KSpace>::findABel(kspace, input_set);
220 
221  typedef DGtal::SurfelAdjacency<3> SurfelAdjacency;
222  const SurfelAdjacency surfel_adjacency(true);
223 
224  typedef DGtal::LightImplicitDigitalSurface<DGtal::Z3i::KSpace, DGtal::Z3i::DigitalSet> DigitalSurfaceContainer;
225  const DigitalSurfaceContainer digital_surface_container(kspace, input_set, surfel_adjacency, cell_bel);
226 
227  typedef DGtal::DigitalSurface<DigitalSurfaceContainer> DigitalSurface;
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 
261  typedef DGtal::DiscreteExteriorCalculusFactory<DGtal::EigenLinearAlgebraBackend> CalculusFactory;
262  typedef DGtal::DiscreteExteriorCalculus<2, 3, DGtal::EigenLinearAlgebraBackend> Calculus;
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 
316 int main(int argc, char* argv[])
317 {
318  QApplication app(argc, argv);
319 
320  //pyramid_3d();
321  alcapone_3d();
322 
323  return app.exec();
324 }
ConstIterator
MyDigitalSurface::ConstIterator ConstIterator
Definition: greedy-plane-segmentation-ex2.cpp:93
DGtal::DiscreteExteriorCalculus
Aim: DiscreteExteriorCalculus represents a calculus in the dec package. This is the main structure in...
Definition: DiscreteExteriorCalculus.h:98
DGtal::ConstAlias
Aim: This class encapsulates its parameter class so that to indicate to the user that the object/poin...
Definition: ConstAlias.h:187
DGtal::HyperRectDomain< Space >
DGtal::Trace::endBlock
double endBlock()
DGtal::VolReader::importVol
static ImageContainer importVol(const std::string &filename, const Functor &aFunctor=Functor())
DGtal::ImageContainerBySTLVector
Definition: ImageContainerBySTLVector.h:127
DGtal::DigitalSurface
Aim: Represents a set of n-1-cells in a nD space, together with adjacency relation between these cell...
Definition: DigitalSurface.h:140
DGtal::SurfelAdjacency
Aim: Represent adjacencies between surfel elements, telling if it follows an interior to exterior ord...
Definition: SurfelAdjacency.h:66
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::Display3DFactory::draw
static void draw(Display3D< Space, KSpace > &display, const DGtal::DiscreteExteriorCalculus< dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger > &calculus)
DGtal::trace
Trace trace
Definition: Common.h:150
DGtal::DiscreteExteriorCalculusSolver::compute
DiscreteExteriorCalculusSolver & compute(const Operator &linear_operator)
DGtal::Display3D< SpaceND< 3 >, KhalimskySpaceND< 3 > >::updateDisplay
@ updateDisplay
Definition: Display3D.h:249
DGtal::Trace::beginBlock
void beginBlock(const std::string &keyword="")
Functor
InHalfPlaneBySimple3x3Matrix< Point, double > Functor
Definition: testConvexHull2DReverse.cpp:51
DGtal::PointVector< dim, Integer >::diagonal
static Self diagonal(Component val=1)
DGtal::functors::Cast
Aim: Define a simple functor using the static cast operator.
Definition: BasicFunctors.h:394
DGtal::SignedKhalimskyCell
Represents a signed cell in a cellular grid space by its Khalimsky coordinates and a boolean value.
Definition: KhalimskySpaceND.h:209
DGtal::PRIMAL
@ PRIMAL
Definition: Duality.h:61
DGtal::Trace::info
std::ostream & info()
DGtal::LightImplicitDigitalSurface
Aim: A model of CDigitalSurfaceContainer which defines the digital surface as the boundary of an impl...
Definition: LightImplicitDigitalSurface.h:116
DGtal::Viewer3D
Definition: Viewer3D.h:133
DGtal::concepts::CDomain
Aim: This concept represents a digital domain, i.e. a non mutable subset of points of the given digit...
Definition: CDomain.h:130
Image
ImageContainerBySTLVector< Domain, Value > Image
Definition: testSimpleRandomAccessRangeFromPoint.cpp:45
DGtal::Surfaces::findABel
static SCell findABel(const KSpace &K, const PointPredicate &pp, unsigned int nbtries=1000)
Domain
HyperRectDomain< Space > Domain
Definition: testSimpleRandomAccessRangeFromPoint.cpp:44
DGtal::KhalimskySpaceND::sKCoords
const Point & sKCoords(const SCell &c) const
Return its Khalimsky coordinates.
DGtal::DiscreteExteriorCalculusFactory
Aim: This class provides static members to create DEC structures from various other DGtal structures.
Definition: DiscreteExteriorCalculusFactory.h:62
DGtal::DiscreteExteriorCalculusSolver::solve
SolutionKForm solve(const InputKForm &input_kform) const
main
int main(int argc, char **argv)
Definition: testArithmeticDSS-benchmark.cpp:147
domain
Domain domain
Definition: testProjection.cpp:88
DGtal::PointVector< dim, Integer >
DGtal::DiscreteExteriorCalculusSolver
Aim: This wraps a linear algebra solver around a discrete exterior calculus.
Definition: DiscreteExteriorCalculusSolver.h:70
DGtal::DUAL
@ DUAL
Definition: Duality.h:62
DGtal::concepts::CPointPredicate
Aim: Defines a predicate on a point.
Definition: CPointPredicate.h:81
DGtal::Viewer3D::show
virtual void show()
Overload QWidget method in order to add a call to updateList() method (to ensure that the lists are w...
Point
MyPointD Point
Definition: testClone2.cpp:383
DGtal::EigenLinearAlgebraBackend::SolverSparseLU
Eigen::SparseLU< SparseMatrix > SolverSparseLU
Definition: EigenSupport.h:110
DGtal::Display3DFactory
Factory for GPL Display3D:
Definition: Display3DFactory.h:87
DGtal::DigitalSetByAssociativeContainer
Aim: A wrapper class around a STL associative container for storing sets of digital points within som...
Definition: DigitalSetByAssociativeContainer.h:90
DGtal::KhalimskySpaceND
Aim: This class is a model of CCellularGridSpaceND. It represents the cubical grid as a cell complex,...
Definition: KhalimskySpaceND.h:394
Generated on Fri Oct 9 2020 08:57:53 for DGtal by doxygen 1.8.20