DGtal  0.9.2
exampleDECSurface.cpp
1 #include "DGtal/io/readers/VolReader.h"
3 #include "DGtal/io/viewers/Viewer3D.h"
4 #include "DGtal/topology/SurfelAdjacency.h"
5 #include "DGtal/topology/LightImplicitDigitalSurface.h"
6 #include "DGtal/topology/DigitalSurface.h"
7 #include "DGtal/math/linalg/EigenSupport.h"
8 #include "DGtal/dec/DiscreteExteriorCalculus.h"
9 #include "DGtal/dec/DiscreteExteriorCalculusSolver.h"
10 #include "DGtal/dec/DiscreteExteriorCalculusFactory.h"
11 
12 #include "ConfigExamples.h"
13 
16 
17 template <typename Predicate, typename Domain>
18 struct FalseOutsideDomain
19 {
20  BOOST_CONCEPT_ASSERT(( DGtal::concepts::CPointPredicate<Predicate> ));
21  BOOST_CONCEPT_ASSERT(( DGtal::concepts::CDomain<Domain> ));
22 
23  typedef typename Predicate::Point Point;
24 
25  FalseOutsideDomain(DGtal::ConstAlias<Predicate> predicate, DGtal::ConstAlias<Domain> domain) :
26  predicate(&predicate), domain(&domain)
27  {
28  }
29 
30  bool
31  operator()(const Point& point) const
32  {
33  if (!domain->isInside(point)) return false;
34  return (*predicate)(point);
35  }
36 
37  const Predicate* predicate;
38  const Domain* domain;
39 };
40 
41 void
42 alcapone_3d()
43 {
44  using std::endl;
45  using DGtal::trace;
46 
47  trace.beginBlock("alcapone");
48 
49  const std::string filename = examplesPath + "samples/Al.100.vol";
50 
51  trace.beginBlock("digital surface");
52 
55  typedef DGtal::functors::Cast<bool> Functor;
56  const Image image = DGtal::VolReader<Image, Functor>::importVol(filename, Functor());
57  const DGtal::Z3i::Domain domain = image.domain();
58 
59  trace.info() << "domain=" << domain << endl;
60 
61  typedef FalseOutsideDomain<Image, DGtal::Z3i::Domain> ImageExtended;
62  const ImageExtended image_extended(image, domain);
63 
64  DGtal::Z3i::KSpace kspace;
65  kspace.init(domain.lowerBound(), domain.upperBound()+DGtal::Z3i::Point(0,0,1), true);
66 
67  const DGtal::Z3i::KSpace::SCell cell_bel = DGtal::Surfaces<DGtal::Z3i::KSpace>::findABel(kspace, image_extended);
68 
69  typedef DGtal::SurfelAdjacency<3> SurfelAdjacency;
70  const SurfelAdjacency surfel_adjacency(true);
71 
73  const DigitalSurfaceContainer digital_surface_container(kspace, image_extended, surfel_adjacency, cell_bel);
74 
76  const DigitalSurface digital_surface(digital_surface_container);
78 
79  trace.info() << "digital_surface_size=" << digital_surface.size() << endl;
80 
81  {
82  Viewer* viewer = new Viewer(kspace);
83  viewer->show();
84  viewer->setWindowTitle("alcapone surface");
85  for (DigitalSurface::ConstIterator si=digital_surface.begin(), se=digital_surface.end(); si!=se; si++)
86  {
87  const DGtal::Z3i::KSpace::SCell cell = *si;
88  (*viewer) << cell;
89  }
90  (*viewer) << Viewer::updateDisplay;
91  }
92 
93  trace.endBlock();
94 
95  trace.beginBlock("discrete exterior calculus");
96 
100  const Calculus calculus = CalculusFactory::createFromNSCells<2>(digital_surface.begin(), digital_surface.end(), false);
102  trace.info() << "calculus=" << calculus << endl;
103 
104  {
105  Viewer* viewer = new Viewer(kspace);
106  viewer->show();
107  viewer->setWindowTitle("alcapone calculus");
108  DisplayFactory::draw(*viewer, calculus);
109  (*viewer) << Viewer::updateDisplay;
110  }
111 
112  using DGtal::PRIMAL;
113  using DGtal::DUAL;
114 
115  trace.endBlock();
116 
117  trace.beginBlock("poisson equation");
118 
120  Calculus::DualForm0 rho(calculus);
121  for (int index=0; index<rho.length(); index++)
122  {
123  const Calculus::SCell cell = rho.getSCell(index);
124  const Calculus::Point point = kspace.sKCoords(cell);
125  if (point[2]>1) continue;
126  rho.myContainer(index) = point[1]>100 ? 1 : -1;
127  }
128  trace.info() << "rho_range=" << rho.myContainer.minCoeff() << "/" << rho.myContainer.maxCoeff() << endl;
130 
131  {
132  Viewer* viewer = new Viewer(kspace);
133  viewer->show();
134  viewer->setWindowTitle("alcapone poisson rho");
135  DisplayFactory::draw(*viewer, rho);
136  (*viewer) << Viewer::updateDisplay;
137  }
138 
140  const Calculus::DualIdentity0 laplace = calculus.laplace<DUAL>();
141 
142  typedef DGtal::EigenLinearAlgebraBackend::SolverSparseLU LinearAlgebraSolver;
144  PoissonSolver solver;
145  solver.compute(laplace);
146  const Calculus::DualForm0 phi = solver.solve(rho);
147  trace.info() << "phi_range=" << phi.myContainer.minCoeff() << "/" << phi.myContainer.maxCoeff() << endl;
149 
150  {
151  Viewer* viewer = new Viewer(kspace);
152  viewer->show();
153  viewer->setWindowTitle("alcapone poisson phi");
154  DisplayFactory::draw(*viewer, phi);
155  (*viewer) << Viewer::updateDisplay;
156  }
157 
158  trace.endBlock();
159 
160  trace.endBlock();
161 }
162 
163 void
164 pyramid_3d()
165 {
166  using std::endl;
167  using DGtal::trace;
168 
169  trace.beginBlock("pyramid example");
170 
171  trace.beginBlock("digital surface");
172 
174  trace.info() << "input_domain=" << input_domain << endl;
175 
176  DGtal::Z3i::KSpace kspace;
177  kspace.init(input_domain.lowerBound(), input_domain.upperBound(), true);
178 
180  DGtal::Z3i::DigitalSet input_set(input_domain);
181  input_set.insert( DGtal::Z3i::Point(0,0,0) );
182  input_set.insert( DGtal::Z3i::Point(1,0,0) );
183  input_set.insert( DGtal::Z3i::Point(1,1,0) );
184  input_set.insert( DGtal::Z3i::Point(0,1,0) );
185  input_set.insert( DGtal::Z3i::Point(0,0,1) );
187  trace.info() << "input_set_size=" << input_set.size() << endl;
188 
189  {
190  Viewer* viewer = new Viewer(kspace);
191  viewer->show();
192  viewer->setWindowTitle("input set");
193  (*viewer) << input_set;
194  (*viewer) << Viewer::updateDisplay;
195  }
196 
199 
200  typedef DGtal::SurfelAdjacency<3> SurfelAdjacency;
201  const SurfelAdjacency surfel_adjacency(true);
202 
204  const DigitalSurfaceContainer digital_surface_container(kspace, input_set, surfel_adjacency, cell_bel);
205 
206  typedef DGtal::DigitalSurface<DigitalSurfaceContainer> DigitalSurface;
207  const DigitalSurface digital_surface(digital_surface_container);
209 
210  trace.info() << "surfel_adjacency=" << endl;
211  for (int ii=0; ii<3; ii++)
212  {
213  std::stringstream ss;
214  for (int jj=0; jj<3; jj++)
215  ss << surfel_adjacency.getAdjacency(ii, jj) << " ";
216  trace.info() << ss.str() << endl;
217  }
218 
219  trace.info() << "digital_surface_container=" << digital_surface_container << endl;
220 
221  trace.info() << "digital_surface_size=" << digital_surface.size() << endl;
222 
223  {
224  Viewer* viewer = new Viewer(kspace);
225  viewer->show();
226  viewer->setWindowTitle("digital surface");
227  for (DigitalSurface::ConstIterator si=digital_surface.begin(), se=digital_surface.end(); si!=se; si++)
228  {
229  const DGtal::Z3i::KSpace::SCell cell = *si;
230  (*viewer) << cell;
231  }
232  (*viewer) << Viewer::updateDisplay;
233  }
234 
235  trace.endBlock();
236 
237  trace.beginBlock("discrete exterior calculus");
238 
242  const Calculus calculus = CalculusFactory::createFromNSCells<2>(digital_surface.begin(), digital_surface.end());
244  trace.info() << "calculus=" << calculus << endl;
245 
246  {
247  Viewer* viewer = new Viewer(kspace);
248  viewer->show();
249  viewer->setWindowTitle("discrete exterior calculus");
250  DisplayFactory::draw(*viewer, calculus);
251  (*viewer) << Viewer::updateDisplay;
252  }
253 
254  using DGtal::PRIMAL;
255  using DGtal::DUAL;
256 
257  const Calculus::PrimalForm2 primal_surfel_area = calculus.hodge<0, DUAL>() *
258  Calculus::DualForm0(calculus, Eigen::VectorXd::Ones(calculus.kFormLength(0, DUAL)));
259 
260  const Calculus::DualForm2 dual_surfel_area = calculus.hodge<0, PRIMAL>() *
261  Calculus::PrimalForm0(calculus, Eigen::VectorXd::Ones(calculus.kFormLength(0, PRIMAL)));
262 
263  const double area_th = calculus.kFormLength(2, PRIMAL);
264  const double area_primal = dual_surfel_area.myContainer.sum();
265  const double area_dual = primal_surfel_area.myContainer.sum();
266  trace.info() << "area_th=" << area_th << endl;
267  trace.info() << "area_primal=" << area_primal << endl;
268  trace.info() << "area_dual=" << area_dual << endl;
269  FATAL_ERROR( area_th == area_primal );
270  FATAL_ERROR( area_th == area_dual );
271 
272  const Calculus::DualForm1 dual_edge_length = calculus.hodge<1, PRIMAL>() *
273  Calculus::PrimalForm1(calculus, Eigen::VectorXd::Ones(calculus.kFormLength(1, PRIMAL)));
274 
275  const bool dual_edge_length_match = dual_edge_length.myContainer == Eigen::VectorXd::Ones(calculus.kFormLength(1, DUAL)) ;
276  trace.info() << "dual_edge_length_match=" << dual_edge_length_match << endl;
277  FATAL_ERROR( dual_edge_length_match );
278 
279  trace.beginBlock("dual surfel area");
280 
281  for (Calculus::Index index=0; index<dual_surfel_area.length(); index++)
282  {
283  const Calculus::SCell cell = dual_surfel_area.getSCell(index);
284  const Calculus::Scalar area = dual_surfel_area.myContainer(index);
285  trace.info() << index << " " << cell << " " << area << endl;
286  }
287 
288  trace.endBlock();
289 
290  trace.endBlock();
291 
292  trace.endBlock();
293 }
294 
295 int main(int argc, char* argv[])
296 {
297  QApplication app(argc, argv);
298 
299  //pyramid_3d();
300  alcapone_3d();
301 
302  return app.exec();
303 }
304 
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...
Trace trace
Definition: Common.h:130
static Self diagonal(Component val=1)
SolutionKForm solve(const InputKForm &input_kform) const
const Point & sKCoords(const SCell &c) const
virtual void show()
Overload QWidget method in order to add a call to updateList() method (to ensure that the lists are w...
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:
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)
const Point & upperBound() const
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...
Aim: This wraps a linear algebra solver around a discrete exterior calculus.
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
Definition: EigenSupport.h:105