DGtal  1.0.0
exampleDECSurface.cpp
1 
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"
17 
18 #include "ConfigExamples.h"
19 
22 
23 template <typename Predicate, typename Domain>
24 struct FalseOutsideDomain
25 {
26  BOOST_CONCEPT_ASSERT(( DGtal::concepts::CPointPredicate<Predicate> ));
27  BOOST_CONCEPT_ASSERT(( DGtal::concepts::CDomain<Domain> ));
28 
29  typedef typename Predicate::Point Point;
30 
31  FalseOutsideDomain(DGtal::ConstAlias<Predicate> predicate, DGtal::ConstAlias<Domain> adomain) :
32  myPredicate(&predicate), myDomain(&adomain)
33  {
34  }
35 
36  bool
37  operator()(const Point& point) const
38  {
39  if (!myDomain->isInside(point)) return false;
40  return (*myPredicate)(point);
41  }
42 
43  const Predicate* myPredicate;
44  const Domain* myDomain;
45 };
46 
47 void
48 alcapone_3d()
49 {
50  using std::endl;
51  using DGtal::trace;
52 
53  trace.beginBlock("alcapone");
54 
55  const std::string filename = examplesPath + "samples/Al.100.vol";
56 
57  trace.beginBlock("digital surface");
58 
63  const DGtal::Z3i::Domain domain = image.domain();
64 
65  trace.info() << "domain=" << domain << endl;
66 
67  typedef FalseOutsideDomain<Image, DGtal::Z3i::Domain> ImageExtended;
68  const ImageExtended image_extended(image, domain);
69 
70  DGtal::Z3i::KSpace kspace;
71  kspace.init(domain.lowerBound(), domain.upperBound()+DGtal::Z3i::Point(0,0,1), true);
72 
73  const DGtal::Z3i::KSpace::SCell cell_bel = DGtal::Surfaces<DGtal::Z3i::KSpace>::findABel(kspace, image_extended);
74 
75  typedef DGtal::SurfelAdjacency<3> SurfelAdjacency;
76  const SurfelAdjacency surfel_adjacency(true);
77 
79  const DigitalSurfaceContainer digital_surface_container(kspace, image_extended, surfel_adjacency, cell_bel);
80 
82  const DigitalSurface digital_surface(digital_surface_container);
84 
85  trace.info() << "digital_surface_size=" << digital_surface.size() << endl;
86 
87  {
88  Viewer* viewer = new Viewer(kspace);
89  viewer->show();
90  viewer->setWindowTitle("alcapone surface");
91  for (DigitalSurface::ConstIterator si=digital_surface.begin(), se=digital_surface.end(); si!=se; si++)
92  {
93  const DGtal::Z3i::KSpace::SCell cell = *si;
94  (*viewer) << cell;
95  }
96  (*viewer) << Viewer::updateDisplay;
97  }
98 
99  trace.endBlock();
100 
101  trace.beginBlock("discrete exterior calculus");
102 
106  const Calculus calculus = CalculusFactory::createFromNSCells<2>(digital_surface.begin(), digital_surface.end(), false);
108  trace.info() << "calculus=" << calculus << endl;
109 
110  {
111  Viewer* viewer = new Viewer(kspace);
112  viewer->show();
113  viewer->setWindowTitle("alcapone calculus");
114  DisplayFactory::draw(*viewer, calculus);
115  (*viewer) << Viewer::updateDisplay;
116  }
117 
118  using DGtal::PRIMAL;
119  using DGtal::DUAL;
120 
121  trace.endBlock();
122 
123  trace.beginBlock("poisson equation");
124 
126  Calculus::DualForm0 rho(calculus);
127  for (int index=0; index<rho.length(); index++)
128  {
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;
133  }
134  trace.info() << "rho_range=" << rho.myContainer.minCoeff() << "/" << rho.myContainer.maxCoeff() << endl;
136 
137  {
138  Viewer* viewer = new Viewer(kspace);
139  viewer->show();
140  viewer->setWindowTitle("alcapone poisson rho");
141  DisplayFactory::draw(*viewer, rho);
142  (*viewer) << Viewer::updateDisplay;
143  }
144 
146  const Calculus::DualIdentity0 laplace = calculus.laplace<DUAL>();
147 
148  typedef DGtal::EigenLinearAlgebraBackend::SolverSparseLU LinearAlgebraSolver;
150  PoissonSolver solver;
151  solver.compute(laplace);
152  const Calculus::DualForm0 phi = solver.solve(rho);
153  trace.info() << "phi_range=" << phi.myContainer.minCoeff() << "/" << phi.myContainer.maxCoeff() << endl;
155 
156  {
157  Viewer* viewer = new Viewer(kspace);
158  viewer->show();
159  viewer->setWindowTitle("alcapone poisson phi");
160  DisplayFactory::draw(*viewer, phi);
161  (*viewer) << Viewer::updateDisplay;
162  }
163 
164  trace.endBlock();
165 
166  trace.endBlock();
167 }
168 
169 void
170 pyramid_3d()
171 {
172  using std::endl;
173  using DGtal::trace;
174 
175  trace.beginBlock("pyramid example");
176 
177  trace.beginBlock("digital surface");
178 
180  trace.info() << "input_domain=" << input_domain << endl;
181 
182  DGtal::Z3i::KSpace kspace;
183  kspace.init(input_domain.lowerBound(), input_domain.upperBound(), true);
184 
186  DGtal::Z3i::DigitalSet input_set(input_domain);
187  input_set.insert( DGtal::Z3i::Point(0,0,0) );
188  input_set.insert( DGtal::Z3i::Point(1,0,0) );
189  input_set.insert( DGtal::Z3i::Point(1,1,0) );
190  input_set.insert( DGtal::Z3i::Point(0,1,0) );
191  input_set.insert( DGtal::Z3i::Point(0,0,1) );
193  trace.info() << "input_set_size=" << input_set.size() << endl;
194 
195  {
196  Viewer* viewer = new Viewer(kspace);
197  viewer->show();
198  viewer->setWindowTitle("input set");
199  (*viewer) << input_set;
200  (*viewer) << Viewer::updateDisplay;
201  }
202 
205 
206  typedef DGtal::SurfelAdjacency<3> SurfelAdjacency;
207  const SurfelAdjacency surfel_adjacency(true);
208 
210  const DigitalSurfaceContainer digital_surface_container(kspace, input_set, surfel_adjacency, cell_bel);
211 
212  typedef DGtal::DigitalSurface<DigitalSurfaceContainer> DigitalSurface;
213  const DigitalSurface digital_surface(digital_surface_container);
215 
216  trace.info() << "surfel_adjacency=" << endl;
217  for (int ii=0; ii<3; ii++)
218  {
219  std::stringstream ss;
220  for (int jj=0; jj<3; jj++)
221  ss << surfel_adjacency.getAdjacency(ii, jj) << " ";
222  trace.info() << ss.str() << endl;
223  }
224 
225  trace.info() << "digital_surface_container=" << digital_surface_container << endl;
226 
227  trace.info() << "digital_surface_size=" << digital_surface.size() << endl;
228 
229  {
230  Viewer* viewer = new Viewer(kspace);
231  viewer->show();
232  viewer->setWindowTitle("digital surface");
233  for (DigitalSurface::ConstIterator si=digital_surface.begin(), se=digital_surface.end(); si!=se; si++)
234  {
235  const DGtal::Z3i::KSpace::SCell cell = *si;
236  (*viewer) << cell;
237  }
238  (*viewer) << Viewer::updateDisplay;
239  }
240 
241  trace.endBlock();
242 
243  trace.beginBlock("discrete exterior calculus");
244 
248  const Calculus calculus = CalculusFactory::createFromNSCells<2>(digital_surface.begin(), digital_surface.end());
250  trace.info() << "calculus=" << calculus << endl;
251 
252  {
253  Viewer* viewer = new Viewer(kspace);
254  viewer->show();
255  viewer->setWindowTitle("discrete exterior calculus");
256  DisplayFactory::draw(*viewer, calculus);
257  (*viewer) << Viewer::updateDisplay;
258  }
259 
260  using DGtal::PRIMAL;
261  using DGtal::DUAL;
262 
263  const Calculus::PrimalForm2 primal_surfel_area = calculus.hodge<0, DUAL>() *
264  Calculus::DualForm0(calculus, Eigen::VectorXd::Ones(calculus.kFormLength(0, DUAL)));
265 
266  const Calculus::DualForm2 dual_surfel_area = calculus.hodge<0, PRIMAL>() *
267  Calculus::PrimalForm0(calculus, Eigen::VectorXd::Ones(calculus.kFormLength(0, PRIMAL)));
268 
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 );
277 
278  const Calculus::DualForm1 dual_edge_length = calculus.hodge<1, PRIMAL>() *
279  Calculus::PrimalForm1(calculus, Eigen::VectorXd::Ones(calculus.kFormLength(1, PRIMAL)));
280 
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 );
284 
285  trace.beginBlock("dual surfel area");
286 
287  for (Calculus::Index index=0; index<dual_surfel_area.length(); index++)
288  {
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;
292  }
293 
294  trace.endBlock();
295 
296  trace.endBlock();
297 
298  trace.endBlock();
299 }
300 
301 int main(int argc, char* argv[])
302 {
303  QApplication app(argc, argv);
304 
305  //pyramid_3d();
306  alcapone_3d();
307 
308  return app.exec();
309 }
void beginBlock(const std::string &keyword="")
Aim: DiscreteExteriorCalculus represents a calculus in the dec package. This is the main structure in...
HyperRectDomain< Space > Domain
InHalfPlaneBySimple3x3Matrix< Point, double > Functor
Aim: A model of CDigitalSurfaceContainer which defines the digital surface as the boundary of an impl...
MyDigitalSurface::ConstIterator ConstIterator
Trace trace
Definition: Common.h:144
static Self diagonal(Component val=1)
Aim: This class encapsulates its parameter class so that to indicate to the user that the object/poin...
Definition: ConstAlias.h:186
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...
double endBlock()
Aim: Represents a set of n-1-cells in a nD space, together with adjacency relation between these cell...
Factory for GPL Display3D:
Domain domain
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...
Definition: CDomain.h:129
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)
Specifies the upper and lower bounds for the maximal cells in this space.
const Point & sKCoords(const SCell &c) const
Return its Khalimsky coordinates.
int main(int argc, char **argv)
MyPointD Point
Definition: testClone2.cpp:383
Aim: A wrapper class around a STL associative container for storing sets of digital points within som...
Space::Point Point
Definition: StdDefs.h:168
std::ostream & info()
Represents a signed cell in a cellular grid space by its Khalimsky coordinates and a boolean value.
ImageContainerBySTLVector< Domain, Value > Image
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...
SolutionKForm solve(const InputKForm &input_kform) const
Aim: This class is a model of CCellularGridSpaceND. It represents the cubical grid as a cell complex,...
Eigen::SparseLU< SparseMatrix > SolverSparseLU
Definition: EigenSupport.h:110