[ VIGRA Homepage | Function Index | Class Index | Namespaces | File List | Main Page ]

polytope.hxx VIGRA

1 #ifndef VIGRA_POLYTOPE_HXX
2 #define VIGRA_POLYTOPE_HXX
3 
4 #ifndef WITH_LEMON
5  #error "Should only be included with flag \"WITH_LEMON\""
6 #endif
7 
8 #include <set>
9 #include <lemon/list_graph.h>
10 #include <lemon/maps.h>
11 
12 #include "config.hxx"
13 #include "error.hxx"
14 #include "tinyvector.hxx"
15 #include "array_vector.hxx"
16 #include "linear_algebra.hxx"
17 #include "numerictraits.hxx"
18 #include "permutation.hxx"
19 
20 namespace vigra {
21 
22 /** \brief Represent an n-dimensional polytope.
23 
24  \tparam N Dimension the polytope.
25  \tparam T Type of the vector components of the polytope vertices.
26 */
27 template <unsigned int N, class T>
28 class Polytope
29 {
30  public:
31 
32  enum Dimension {dimension = N};
33  enum node_enum {INVALID, FACET, VERTEX};
34 
35  template <node_enum NodeType>
36  struct node_type_iterator;
37 
38  typedef T coordinate_type;
39  typedef typename NumericTraits<T>::RealPromote real_type;
42  typedef typename point_type::difference_type difference_type;
43  typedef typename lemon::ListDigraph graph_type;
44  typedef typename graph_type::Node node_type;
45  typedef typename graph_type::Arc arc_type;
46  typedef typename graph_type::NodeIt node_iterator;
47  typedef typename graph_type::OutArcIt out_arc_iterator;
48  typedef typename graph_type::InArcIt in_arc_iterator;
49  typedef node_type_iterator<FACET> facet_iterator;
50  typedef node_type_iterator<VERTEX> vertex_iterator;
51 
52  /** Default constructor creates an empty polytope class.
53  */
55  : graph_()
56  , type_map_(graph_)
57  , vec_map_(graph_)
58  , aligns_map_(graph_)
59  {}
60 
61  /** Copy constructor.
62  */
63  Polytope(const Polytope<N, T> & other)
64  : graph_()
65  , type_map_(graph_)
66  , vec_map_(graph_)
67  , aligns_map_(graph_)
68  {
69  *this = other;
70  }
71 
72  /** Copy from another polytope.
73  */
74  virtual void operator=(const Polytope<N, T> & other)
75  {
76  lemon::digraphCopy(other.graph_, graph_);
77  lemon::mapCopy(other.graph_, other.type_map_, type_map_);
78  lemon::mapCopy(other.graph_, other.vec_map_, vec_map_);
79  lemon::mapCopy(other.graph_, other.aligns_map_, aligns_map_);
80  }
81 
82  virtual bool contains(const point_view_type & p) const = 0;
83 
84  virtual real_type nVolume() const = 0;
85 
86  virtual real_type nSurface() const = 0;
87 
88  /** Check if the facet aligns with other facets at each of its ridges.
89  */
90  virtual bool closed(const node_type n) const
91  {
92  vigra_assert(
93  type_map_[n] == FACET,
94  "Polytope::closed(): Node needs do be a facet");
95  return std::find(
96  aligns_map_[n].begin(),
97  aligns_map_[n].end(),
98  lemon::INVALID) == aligns_map_[n].end();
99  }
100 
101  /** Check if the polytope has a closed surface
102  */
103  virtual bool closed() const
104  {
105  for (facet_iterator n(graph_, type_map_); n != lemon::INVALID; ++n)
106  {
107  if (!(this->closed(n)))
108  {
109  return false;
110  }
111  }
112  return true;
113  }
114 
115 
116  /** Add a vertex to the polytope.
117  */
118  virtual node_type addVertex(const point_view_type & p)
119  {
120  node_type ret = graph_.addNode();
121  type_map_[ret] = VERTEX;
122  vec_map_[ret] = p;
123  for (int i = 0; i < N; i++)
124  {
125  aligns_map_[ret][i] = lemon::INVALID;
126  }
127  return ret;
128  }
129 
130  /** Erase a facet.
131  */
132  virtual void eraseFacet(const node_type u)
133  {
134  vigra_assert(
135  type_map_[u] == FACET,
136  "Polytope::eraseFacet(): Node needs to be a facet");
137  for (auto neighbor : aligns_map_[u])
138  {
139  if (neighbor != lemon::INVALID)
140  {
141  auto it = std::find(
142  aligns_map_[neighbor].begin(),
143  aligns_map_[neighbor].end(),
144  u);
145  vigra_assert(
146  it != aligns_map_[neighbor].end(),
147  "Polytope::eraseFacet(): Inconsistent aligns map");
148  *it = lemon::INVALID;
149  }
150  }
151  graph_.erase(u);
152  }
153 
154  /** Get the connected elements in the graph that represents the polytope.
155  If a facet node is inserted, all of its vertices will be returned, if
156  a vertex node is inserted, all facets having this vertex will be
157  returned.
158  */
159  virtual std::set<node_type> getConnected(const node_type u) const
160  {
161  std::set<node_type> ret;
162  if (type_map_[u] == FACET)
163  {
164  for (out_arc_iterator a(graph_, u); a != lemon::INVALID; ++a)
165  {
166  ret.insert(graph_.target(a));
167  }
168  }
169  else
170  {
171  for (in_arc_iterator a(graph_, u); a != lemon::INVALID; ++a)
172  {
173  ret.insert(graph_.source(a));
174  }
175  }
176  return ret;
177  }
178 
179  // TODO remove
180  virtual ArrayVector<point_view_type> getVertices(const node_type u) const
181  {
182  vigra_assert(
183  type_map_[u] == FACET,
184  "Polytope::getVertices(): Node must be a facet");
186  for (out_arc_iterator a(graph_, u); a != lemon::INVALID; ++a)
187  {
188  ret.push_back(vec_map_[graph_.target(a)]);
189  }
190  return ret;
191  }
192 
193  /** Get all facets whose normal has a positive scalar product with the
194  vector to the given vertex.
195  */
197  {
199  for (facet_iterator n(graph_, type_map_); n != lemon::INVALID; ++n)
200  {
201  if (distance(n, p) > 0)
202  {
203  ret.push_back(n);
204  }
205  }
206  return ret;
207  }
208 
209  /** Remove all vertices that are not part of the polytope mesh.
210  */
211  virtual void tidyUp()
212  {
213  std::set<node_type> to_erase;
214  for (vertex_iterator v(graph_, type_map_); v != lemon::INVALID; ++v)
215  {
216  vigra_assert(
217  type_map_[v] == VERTEX,
218  "Polytope::tidyUp(): vertex not a vertex");
219  in_arc_iterator a(graph_, v);
220  if (a == lemon::INVALID)
221  {
222  to_erase.insert(v);
223  }
224  }
225  for (node_type v : to_erase)
226  {
227  graph_.erase(v);
228  }
229  }
230 
231  /** Get the distance between a facet and a vertex */
232  virtual real_type distance(const node_type u, const point_view_type & p) const
233  {
234  vigra_assert(
235  type_map_[u] == FACET,
236  "Polytope::distance(): Node must be a facet");
237  out_arc_iterator a(graph_, u);
238  vigra_assert(
239  a != lemon::INVALID,
240  "Polytope::distance(): Invalid facet");
241 
242  return dot(p - vec_map_[graph_.target(a)], vec_map_[u]);
243  }
244 
245  /** Label all elements in the array which are inside the polytope.
246  */
247  virtual unsigned int fill(
249  const unsigned int label,
250  const point_view_type offset,
251  const point_view_type scale) const
252  {
253  typedef MultiArrayView<N, unsigned int> array_type;
254 
255  unsigned int ret = 0;
256  typename array_type::iterator it = array.begin();
257  for (it = array.begin(); it != array.end(); it++)
258  {
259  const typename array_type::difference_type coord = it.template get<0>();
260  point_type vec;
261  for (unsigned int i = 0; i < vec.size(); i++)
262  {
263  vec[i] = coord[i]*scale[i] + offset[i];
264  }
265  if (this->contains(vec))
266  {
267  ret++;
268  *it = label;
269  }
270  }
271  return ret;
272  }
273 
274  /** Label all elements in the array which are inside the polytope.
275  */
276  virtual unsigned int fill(
278  const unsigned int label,
279  const point_view_type offset) const
280  {
281  vigra_assert(
282  closed(),
283  "Polytope::fill(): Polytope not closed.");
284  typedef MultiArrayView<N, unsigned int> array_type;
285 
286  unsigned int ret = 0;
287  typename array_type::iterator it = array.begin();
288  for (it = array.begin(); it != array.end(); it++)
289  {
290  const typename array_type::difference_type coord = it.template get<0>();
291  point_type vec;
292  for (unsigned int i = 0; i < vec.size(); i++)
293  {
294  vec[i] = coord[i] + offset[i];
295  }
296  if (this->contains(vec))
297  {
298  ret++;
299  *it = label;
300  }
301  }
302  return ret;
303  }
304 
305  protected:
306 
307  virtual bool isConnected(
308  const node_type vertex,
309  const node_type facet) const
310  {
311  vigra_assert(
312  type_map_[vertex] == VERTEX,
313  "Polytope::isConnected(): First node must be a vertex");
314  vigra_assert(
315  type_map_[facet] == FACET,
316  "Polytope::isConnected(): Second node must be a facet");
317  for (out_arc_iterator a(graph_, facet); a != lemon::INVALID; ++a)
318  {
319  if (graph_.target(a) == vertex)
320  {
321  return true;
322  }
323  }
324  return false;
325  }
326 
327  virtual node_type findNeighbor(
328  const node_type u,
329  const difference_type index) const
330  {
331  vigra_assert(
332  type_map_[u] == FACET,
333  "Polytope::findNeighbor(): Node must be a facet");
334  vigra_assert(
335  index < dimension,
336  "Polytope::findNeighbor(): Invalid index");
337  vigra_assert(
338  countOutArcs(graph_, u) == dimension,
339  "Polytope::findNeighbor(): Bad facet");
340  out_skip_iterator a(graph_, u, index);
341  const node_type first_vertex = graph_.target(a);
342  for (node_type candidate : getConnected(first_vertex))
343  {
344  if (candidate != u)
345  {
346  out_skip_iterator b(a);
347  do
348  {
349  ++b;
350  if (b == lemon::INVALID)
351  {
352  return candidate;
353  }
354  } while (isConnected(graph_.target(b), candidate));
355  }
356  }
357  return lemon::INVALID;
358  }
359 
360  void assignNeighbors(const node_type u)
361  {
362  vigra_assert(
363  type_map_[u] == FACET,
364  "Polytope::assignNeighbors(): Node must be facet");
365  for (int i = 0; i < dimension; i++)
366  {
367  aligns_map_[u][i] = this->findNeighbor(u, i);
368  }
369  }
370 
371  void updateInvalidNeighbors(const node_type u)
372  {
373  vigra_assert(
374  type_map_[u] == FACET,
375  "Polytope::assignNeighbors(): Node must be facet");
376  for (int i = 0; i < dimension; i++)
377  {
378  if (aligns_map_[u][i] == lemon::INVALID)
379  {
380  aligns_map_[u][i] = this->findNeighbor(u, i);
381  }
382  }
383  }
384 
385  ArrayVector<node_type> openEdge(const node_type u)
386  {
387  vigra_assert(
388  type_map_[u] == FACET,
389  "Polytope::openEdge(): Node must be facet");
390  vigra_assert(
391  lemon::countOutArcs(graph_, u) == dimension,
392  "Polytope::openEdge(): Got invalid facet");
393  ArrayVector<node_type> ret;
394  for (int i = 0; i < dimension; i++)
395  {
396  if (aligns_map_[u][i] == lemon::INVALID)
397  {
398  for (out_skip_iterator a(graph_, u, i); a != lemon::INVALID; ++a)
399  {
400  ret.push_back(graph_.target(a));
401  }
402  return ret;
403  }
404  }
405  return ret;
406  }
407 
408  public:
409 
410  template <node_enum NodeType>
411  struct node_type_iterator : public node_type
412  {
413  node_type_iterator()
414  {}
415 
416  explicit node_type_iterator(
417  const graph_type & graph,
418  const typename graph_type::NodeMap<node_enum> & type_map)
419  : graph_(graph)
420  , type_map_(type_map)
421  {
422  graph_.first(static_cast<node_type &>(*this));
423  while (*this != lemon::INVALID && type_map_[*this] != NodeType)
424  {
425  graph_.next(*this);
426  }
427  }
428 
429  node_type_iterator<NodeType> & operator++()
430  {
431  while (*this != lemon::INVALID)
432  {
433  graph_.next(*this);
434  if (type_map_[*this] == NodeType)
435  {
436  return *this;
437  }
438  }
439  return *this;
440  }
441 
442  bool operator==(lemon::Invalid i) const
443  {
444  return (static_cast<node_type>(*this) == i);
445  }
446 
447  bool operator!=(lemon::Invalid i) const
448  {
449  return (static_cast<node_type>(*this) != i);
450  }
451 
452  const graph_type & graph_;
453  const typename graph_type::NodeMap<node_enum> & type_map_;
454  };
455 
456  struct out_skip_iterator : public arc_type
457  {
458  out_skip_iterator()
459  {}
460 
461  explicit out_skip_iterator(
462  const graph_type & graph,
463  const node_type & node,
464  const difference_type skip)
465  : graph_(graph)
466  , skip_(skip)
467  , index_(0)
468  {
469  graph_.firstOut(*this, node);
470  if (skip_ == 0)
471  {
472  graph_.nextOut(*this);
473  }
474  }
475 
476  out_skip_iterator & operator++()
477  {
478  ++index_;
479  graph_.nextOut(*this);
480  if (index_ == skip_)
481  {
482  graph_.nextOut(*this);
483  }
484  return *this;
485  }
486 
487  bool operator==(lemon::Invalid i) const
488  {
489  return (static_cast<arc_type>(*this) == i);
490  }
491 
492  bool operator!=(lemon::Invalid i) const
493  {
494  return (static_cast<arc_type>(*this) != i);
495  }
496 
497  difference_type index() const
498  {
499  return index_;
500  }
501 
502  const graph_type & graph_;
503  const difference_type skip_;
504  difference_type index_;
505  };
506 
507  graph_type graph_;
508  typename graph_type::NodeMap<node_enum> type_map_;
509  typename graph_type::NodeMap<point_type> vec_map_;
510  typename graph_type::NodeMap<TinyVector<node_type, N> > aligns_map_;
511 };
512 
513 /** \brief Specialization of the polytope to polytopes which forms a star
514  domain.
515 */
516 template <unsigned int N, class T>
517 class StarPolytope : public Polytope<N, T>
518 {
519  public:
520 
521  typedef Polytope<N, T> base_type;
522  typedef typename base_type::coordinate_type coordinate_type;
523  typedef typename base_type::real_type real_type;
524  typedef typename base_type::point_type point_type;
526  typedef typename base_type::difference_type difference_type;
527  typedef typename base_type::graph_type graph_type;
528  typedef typename base_type::node_type node_type;
529  typedef typename base_type::arc_type arc_type;
530  typedef typename base_type::node_iterator node_iterator;
531  typedef typename base_type::in_arc_iterator in_arc_iterator;
532  typedef typename base_type::out_arc_iterator out_arc_iterator;
533  typedef typename base_type::out_skip_iterator out_skip_iterator;
534  typedef typename base_type::facet_iterator facet_iterator;
535  typedef typename base_type::vertex_iterator vertex_iterator;
536 
537  using base_type::dimension;
538  using base_type::graph_;
539  using base_type::vec_map_;
540  using base_type::type_map_;
541  using base_type::aligns_map_;
542  using base_type::INVALID;
543  using base_type::FACET;
544  using base_type::VERTEX;
545 
546  public:
547 
548  /** Constructor creates an empty StarPolytope with its center a the orign.
549  */
551  : base_type()
552  , center_(point_type())
553  {}
554 
555  /** Copy constructor.
556  */
558  : base_type()
559  , center_(center)
560  {}
561 
562  /** Constructor for the 2-dimensional case taking three vertices and the
563  center.
564  */
566  const point_view_type & a,
567  const point_view_type & b,
568  const point_view_type & c,
569  const point_view_type & center)
570  : base_type()
571  , center_(center)
572  {
573  vigra_precondition(
574  dimension == 2,
575  "StarPolytope::StarPolytope(): Signature only for use in 2D");
576  node_type na = this->addVertex(a);
577  node_type nb = this->addVertex(b);
578  node_type nc = this->addVertex(c);
579  this->addFacet(nb, nc);
580  this->addFacet(na, nc);
581  this->addFacet(na, nb);
582  }
583 
584  /** Constructor for the 3-dimensional case taking four vertices and the
585  center.
586  */
588  const point_view_type & a,
589  const point_view_type & b,
590  const point_view_type & c,
591  const point_view_type & d,
592  const point_view_type & center)
593  : base_type()
594  , center_(center)
595  {
596  vigra_precondition(
597  dimension == 3,
598  "StarPolytope::StarPolytope(): Signature only for use in 3D");
599  node_type na = this->addVertex(a);
600  node_type nb = this->addVertex(b);
601  node_type nc = this->addVertex(c);
602  node_type nd = this->addVertex(d);
603  this->addFacet(nb, nc, nd);
604  this->addFacet(na, nc, nd);
605  this->addFacet(na, nb, nd);
606  this->addFacet(na, nb, nc);
607  }
608 
609  /** Get the center of the star domain.
610  */
611  virtual point_type getCenter() const
612  {
613  return center_;
614  }
615 
616  virtual void assignNormal(const node_type & u)
617  {
618  vigra_assert(
619  type_map_[u] == FACET,
620  "StarPolytope::assignNormal(): Node needs to be a facet node");
621  MultiArray<2, real_type> mat(dimension-1, dimension);
622  out_arc_iterator a(graph_, u);
623  point_view_type vertex = vec_map_[graph_.target(a)];
624  ++a;
625  for (int i = 0; a != lemon::INVALID; ++a, ++i)
626  {
627  const point_type vec = vec_map_[graph_.target(a)] - vertex;
628  std::copy(vec.begin(), vec.end(), rowVector(mat, i).begin());
629  }
630  point_view_type normal = vec_map_[u];
631  for (int i = 0; i < dimension; i++)
632  {
633  normal[i] = 0;
634  }
635  for (auto permutation : permutations_)
636  {
637  coordinate_type val = 1;
638  for (int i = 0; i < dimension - 1; i++)
639  {
640  val *= mat(i, permutation[i]);
641  }
642  val *= permutation.sign();
643  normal[permutation[dimension - 1]] += val;
644  }
645  if (dot(normal, vertex - center_) < 0)
646  {
647  normal *= -1;
648  }
649  }
650 
651  /** Add a facet to a 2-dimensional polytope.
652  */
653  virtual node_type addFacet(const node_type & a, const node_type & b)
654  {
655  vigra_precondition(
656  dimension == 2,
657  "StarPolytope::addFacet(): Signature only for use in 2D");
658  node_type ret = graph_.addNode();
659  type_map_[ret] = FACET;
660  graph_.addArc(ret, a);
661  graph_.addArc(ret, b);
662  vigra_assert(
663  lemon::countOutArcs(graph_, ret) == dimension,
664  "StarPolytope::addFacet(): Invalid facet created");
665  this->assignNormal(ret);
666  this->assignNeighbors(ret);
667  for (auto facet : aligns_map_[ret])
668  {
669  if (facet != lemon::INVALID)
670  {
671  vigra_assert(
672  type_map_[facet] == FACET,
673  "StarPolytope::addFacet(): Node must be facet");
674  this->updateInvalidNeighbors(facet);
675  }
676  }
677  return ret;
678  }
679 
680  /** Add a facet to a 3-dimensional polytope.
681  */
682  virtual node_type addFacet(
683  const node_type & a,
684  const node_type & b,
685  const node_type & c)
686  {
687  vigra_precondition(
688  dimension == 3,
689  "StarPolytope::addFacet(): Signature only for use in 3D");
690  node_type ret = graph_.addNode();
691  type_map_[ret] = FACET;
692  graph_.addArc(ret, a);
693  graph_.addArc(ret, b);
694  graph_.addArc(ret, c);
695  vigra_assert(
696  lemon::countOutArcs(graph_, ret) == dimension,
697  "StarPolytope::addFacet(): Invalid facet created");
698  this->assignNormal(ret);
699  this->assignNeighbors(ret);
700  for (auto facet : aligns_map_[ret])
701  {
702  if (facet != lemon::INVALID)
703  {
704  vigra_assert(
705  type_map_[facet] == FACET,
706  "StarPolytope::addFacet(): Node must be facet");
707  this->updateInvalidNeighbors(facet);
708  }
709  }
710  return ret;
711  }
712 
713  virtual void close()
714  {
715  vigra_precondition(
716  lemon::countNodes(graph_) == dimension + 1,
717  "StarPolytope::close(): Can only close for dim+1 vertices");
718  // Set center of polytope
719  {
720  vertex_iterator v(graph_, type_map_);
721  center_ = vec_map_[v];
722  for (++v; v != lemon::INVALID; ++v)
723  {
724  center_ += vec_map_[v];
725  }
726  center_ /= static_cast<real_type>(dimension + 1);
727  }
728  // Create facets
729  for (int i = 0; i < dimension + 1; i++)
730  {
731  node_type facet = graph_.addNode();
732  type_map_[facet] = FACET;
733  vertex_iterator v(graph_, type_map_);
734  for (int j = 0; j < dimension; ++j, ++v)
735  {
736  if (i == j)
737  {
738  ++v;
739  }
740  graph_.addArc(facet, v);
741  }
742  vigra_assert(
743  lemon::countOutArcs(graph_, facet) == dimension,
744  "StarPolytope::close(): Invalid facet created");
745  this->assignNormal(facet);
746  this->assignNeighbors(facet);
747  for (auto neighbor : aligns_map_[facet])
748  {
749  if (neighbor != lemon::INVALID)
750  {
751  vigra_assert(
752  type_map_[neighbor] == FACET,
753  "StarPolytope::close(): Node must be facet");
754  this->updateInvalidNeighbors(neighbor);
755  }
756  }
757  }
758  }
759 
760  virtual bool contains(const node_type & n, const point_view_type & p) const
761  {
762  vigra_assert(
763  type_map_[n] == FACET,
764  "StarPolytope::contains(): Node needs do be a facet");
765  ArrayVector<point_view_type> vertices = this->getVertices(n);
766  vertices.push_back(center_);
767  MultiArray<2, coordinate_type> jp_mat(dimension, dimension);
768  MultiArray<2, coordinate_type> jj_mat(dimension, dimension);
769  for (int j = 0; j < dimension + 1; j++)
770  {
771  for (int i = 0, ii = 0; i < dimension; i++, ii++)
772  {
773  if (i == j)
774  {
775  ii++;
776  }
777  {
778  const point_type vec = vertices[ii] - p;
779  std::copy(vec.begin(), vec.end(), rowVector(jp_mat, i).begin());
780  }
781  {
782  const point_type vec = vertices[ii] - vertices[j];
783  std::copy(vec.begin(), vec.end(), rowVector(jj_mat, i).begin());
784  }
785  }
786  const coordinate_type jj_det = linalg::determinant(jj_mat);
787  const coordinate_type jp_det = linalg::determinant(jp_mat);
788  const coordinate_type eps = std::numeric_limits<T>::epsilon() * 2;
789  if (((jj_det > 0) != (jp_det > 0)) && abs(jp_det) > eps)
790  {
791  return false;
792  }
793  }
794  return true;
795  }
796 
797  /** Check if a point is inside the polytope.
798  */
799  virtual bool contains(const point_view_type & p) const
800  {
801  for (facet_iterator n(graph_, type_map_); n != lemon::INVALID; ++n)
802  {
803  if (contains(n, p))
804  {
805  return true;
806  }
807  }
808  return false;
809  }
810 
811  virtual real_type nVolume(const node_type & n) const
812  {
813  vigra_assert(
814  type_map_[n] == FACET,
815  "StarPolytope::nVolume(): Node needs do be a facet");
816  MultiArray<2, coordinate_type> mat(dimension, dimension);
817  real_type fac = 1;
818  out_arc_iterator a(graph_, n);
819  for (int i = 0; i < dimension; ++i, ++a)
820  {
821  fac *= (i+1);
822  const point_type vec = vec_map_[graph_.target(a)] - center_;
823  std::copy(vec.begin(), vec.end(), rowVector(mat, i).begin());
824  }
825  return abs(linalg::determinant(mat) / fac);
826  }
827 
828  /** Calculate the volume of the polytope.
829  */
830  virtual real_type nVolume() const
831  {
832  real_type ret = 0;
833  for (facet_iterator n(graph_, type_map_); n != lemon::INVALID; ++n)
834  {
835  ret += this->nVolume(n);
836  }
837  return ret;
838  }
839 
840  virtual real_type nSurface(const node_type & n) const
841  {
842  vigra_assert(
843  type_map_[n] == FACET,
844  "StarPolytope::nVolume(): Node needs do be a facet");
845  MultiArray<2, coordinate_type> mat(dimension, dimension);
846  real_type factor = vec_map_[n].magnitude();
847  out_arc_iterator a(graph_, n);
848  const point_view_type vec = vec_map_[graph_.target(a)];
849  ++a;
850  for (int i = 1; i < dimension; ++i, ++a)
851  {
852  factor *= i;
853  const point_type tmp = vec_map_[graph_.target(a)] - vec;
854  std::copy(tmp.begin(), tmp.end(), rowVector(mat, i).begin());
855  }
856  const point_type tmp = vec_map_[n];
857  std::copy(tmp.begin(), tmp.end(), rowVector(mat, 0).begin());
858  return abs(linalg::determinant(mat)) / factor;
859  }
860 
861  /** Calculate the surface of the polytope.
862  */
863  virtual real_type nSurface() const
864  {
865  real_type ret = 0;
866  for (facet_iterator n(graph_, type_map_); n != lemon::INVALID; ++n)
867  {
868  ret += this->nSurface(n);
869  }
870  return ret;
871  }
872 
873  protected:
874 
875  PlainChangesPermutations<N> permutations_;
876  point_type center_;
877 };
878 
879 /** Specialization of the StarPolytope to polytopes which have a convex domain.
880 */
881 template <unsigned int N, class T>
882 class ConvexPolytope : public StarPolytope<N, T>
883 {
884  public:
885 
887  typedef typename base_type::coordinate_type coordinate_type;
888  typedef typename base_type::real_type real_type;
889  typedef typename base_type::point_type point_type;
891  typedef typename base_type::difference_type difference_type;
892  typedef typename base_type::graph_type graph_type;
893  typedef typename base_type::node_type node_type;
894  typedef typename base_type::arc_type arc_type;
895  typedef typename base_type::node_iterator node_iterator;
896  typedef typename base_type::in_arc_iterator in_arc_iterator;
897  typedef typename base_type::out_arc_iterator out_arc_iterator;
898  typedef typename base_type::out_skip_iterator out_skip_iterator;
899  typedef typename base_type::facet_iterator facet_iterator;
900  typedef typename base_type::vertex_iterator vertex_iterator;
901 
902  using base_type::dimension;
903  using base_type::graph_;
904  using base_type::vec_map_;
905  using base_type::type_map_;
906  using base_type::aligns_map_;
907  using base_type::INVALID;
908  using base_type::FACET;
909  using base_type::VERTEX;
910 
911  public:
912 
914  : base_type()
915  {}
916 
917  ConvexPolytope(const point_view_type & center)
918  : base_type(center)
919  {}
920 
922  const point_view_type & a,
923  const point_view_type & b,
924  const point_view_type & c)
925  : base_type(a, b, c, (a + b + c) / 3)
926  {}
927 
929  const point_view_type & a,
930  const point_view_type & b,
931  const point_view_type & c,
932  const point_view_type & d)
933  : base_type(a, b, c, d, (a + b + c + d) / 4)
934  {}
935 
936  protected:
937 
938  virtual void closeFacet(
939  const node_type & vertex,
940  const node_type & facet)
941  {
942  vigra_assert(
943  type_map_[vertex] == VERTEX,
944  "ConvexPolytope::closeFacet(): Vertex needs to be a vertex node");
945  vigra_assert(
946  type_map_[facet] == FACET,
947  "ConvexPolytope::closeFacet(): Facet needs to be a facet node");
948  vigra_assert(
949  (this->getConnected(facet)).count(vertex) == 0,
950  "ConvexPolytope::closeFacet(): Cannot close facet with vertex");
951 
952  while (!(this->closed(facet)))
953  {
954  ArrayVector<node_type> vertices = this->openEdge(facet);
955  vigra_assert(
956  vertices.size() == (dimension - 1),
957  "StarPolytope::closeFacet(): Invalid facet");
958  node_type new_facet = graph_.addNode();
959  type_map_[new_facet] = FACET;
960  graph_.addArc(new_facet, vertex);
961  for (auto n : vertices)
962  {
963  graph_.addArc(new_facet, n);
964  }
965  vigra_assert(
966  lemon::countOutArcs(graph_, new_facet) == dimension,
967  "ConvexPolytope::closeFacet(): Invalid facet created");
968  this->assignNormal(new_facet);
969  this->assignNeighbors(new_facet);
970  for (auto neighbor : aligns_map_[new_facet])
971  {
972  if (neighbor != lemon::INVALID)
973  {
974  vigra_assert(
975  type_map_[facet] == FACET,
976  "StarPolytope::addFacet(): Node must be facet");
977  this->updateInvalidNeighbors(neighbor);
978  }
979  }
980  }
981  }
982 
983  public:
984 
985  virtual bool contains(const node_type & n, const point_view_type & p) const
986  {
987  vigra_assert(
988  type_map_[n] == FACET,
989  "ConvexPolytope::contains(): Node needs do be a facet");
990  const out_arc_iterator a(graph_, n);
991  const point_view_type vertex = vec_map_[graph_.target(a)];
992  const point_view_type normal = vec_map_[n];
993  const real_type scalar = dot(p - vertex, normal);
994  return (scalar < std::numeric_limits<T>::epsilon() * 2);
995  }
996 
997  virtual bool contains(const point_view_type & p) const
998  {
999  for (facet_iterator n(graph_, type_map_); n != lemon::INVALID; ++n)
1000  {
1001  if (!contains(n, p))
1002  {
1003  return false;
1004  }
1005  }
1006  return true;
1007  }
1008 
1009  /** Expand the polytope to the given point if it's outside of the current
1010  polytope, such that the new polytope is still convex.
1011  */
1012  virtual void addExtremeVertex(const point_view_type & p)
1013  {
1014  vigra_assert(
1015  this->closed(),
1016  "ConvexPolytope::addExtremeVertex(): Polytope needs to be closed");
1017  ArrayVector<node_type> lit_facets = this->litFacets(p);
1018  if (lit_facets.size() > 0)
1019  {
1020  std::set<node_type> open_facets;
1021  for (node_type lit_facet : lit_facets)
1022  {
1023  for (auto con : aligns_map_[lit_facet])
1024  {
1025  if (con != lemon::INVALID)
1026  {
1027  vigra_assert(
1028  type_map_[con] == FACET,
1029  "ConvexPolytope::addExtremeVertex(): "
1030  "facet not a facet");
1031  open_facets.insert(con);
1032  }
1033  }
1034  open_facets.erase(lit_facet);
1035  this->eraseFacet(lit_facet);
1036  }
1037  this->tidyUp();
1038  node_type new_vertex = this->addVertex(p);
1039  for (auto open_facet : open_facets)
1040  {
1041  this->closeFacet(new_vertex, open_facet);
1042  }
1043  }
1044  }
1045 };
1046 
1047 } /* namespace vigra */
1048 
1049 #endif /* VIGRA_POLYTOPE_HXX */
iterator end()
Definition: multi_array.hxx:1937
PromoteTraits< V1, V2 >::Promote dot(RGBValue< V1, RIDX1, GIDX1, BIDX1 > const &r1, RGBValue< V2, RIDX2, GIDX2, BIDX2 > const &r2)
dot product
Definition: rgbvalue.hxx:906
virtual void eraseFacet(const node_type u)
Definition: polytope.hxx:132
virtual bool contains(const point_view_type &p) const
Definition: polytope.hxx:997
virtual unsigned int fill(MultiArrayView< N, unsigned int > &array, const unsigned int label, const point_view_type offset) const
Definition: polytope.hxx:276
virtual node_type addFacet(const node_type &a, const node_type &b, const node_type &c)
Definition: polytope.hxx:682
StarPolytope()
Definition: polytope.hxx:550
iterator begin()
Definition: multi_array.hxx:1921
Main MultiArray class containing the memory management.
Definition: multi_array.hxx:2474
virtual real_type nVolume() const
Definition: polytope.hxx:830
virtual std::set< node_type > getConnected(const node_type u) const
Definition: polytope.hxx:159
StarPolytope(const point_view_type &a, const point_view_type &b, const point_view_type &c, const point_view_type &d, const point_view_type &center)
Definition: polytope.hxx:587
Polytope(const Polytope< N, T > &other)
Definition: polytope.hxx:63
StarPolytope(const point_view_type &center)
Definition: polytope.hxx:557
Definition: array_vector.hxx:58
virtual ArrayVector< node_type > litFacets(const point_view_type &p) const
Definition: polytope.hxx:196
bool operator!=(FFTWComplex< R > const &a, const FFTWComplex< R > &b)
not equal
Definition: fftw3.hxx:841
virtual real_type nSurface() const
Definition: polytope.hxx:863
std::pair< typename vigra::GridGraph< N, DirectedTag >::vertex_iterator, typename vigra::GridGraph< N, DirectedTag >::vertex_iterator > vertices(vigra::GridGraph< N, DirectedTag > const &g)
Get a (begin, end) iterator pair for the vertices of graph g (API: boost).
Definition: multi_gridgraph.hxx:2840
virtual real_type distance(const node_type u, const point_view_type &p) const
Definition: polytope.hxx:232
bool operator==(FFTWComplex< R > const &a, const FFTWComplex< R > &b)
equal
Definition: fftw3.hxx:825
virtual void addExtremeVertex(const point_view_type &p)
Definition: polytope.hxx:1012
Represent an n-dimensional polytope.
Definition: polytope.hxx:28
Wrapper for fixed size vectors.
Definition: tinyvector.hxx:621
virtual unsigned int fill(MultiArrayView< N, unsigned int > &array, const unsigned int label, const point_view_type offset, const point_view_type scale) const
Definition: polytope.hxx:247
Polytope()
Definition: polytope.hxx:54
virtual bool closed(const node_type n) const
Definition: polytope.hxx:90
MultiArrayView< 2, T, C > rowVector(MultiArrayView< 2, T, C > const &m, MultiArrayIndex d)
Definition: matrix.hxx:697
Definition: polytope.hxx:882
FFTWComplex< R >::NormType abs(const FFTWComplex< R > &a)
absolute value (= magnitude)
Definition: fftw3.hxx:1002
Base class for, and view to, vigra::MultiArray.
Definition: multi_array.hxx:704
virtual node_type addVertex(const point_view_type &p)
Definition: polytope.hxx:118
StarPolytope(const point_view_type &a, const point_view_type &b, const point_view_type &c, const point_view_type &center)
Definition: polytope.hxx:565
size_type size() const
Definition: array_vector.hxx:358
virtual void tidyUp()
Definition: polytope.hxx:211
virtual node_type addFacet(const node_type &a, const node_type &b)
Definition: polytope.hxx:653
virtual point_type getCenter() const
Definition: polytope.hxx:611
virtual bool contains(const point_view_type &p) const
Definition: polytope.hxx:799
Specialization of the polytope to polytopes which forms a star domain.
Definition: polytope.hxx:517
NumericTraits< T >::Promote determinant(MultiArrayView< 2, T, C1 > const &a, std::string method="default")
Definition: linear_solve.hxx:859
virtual void operator=(const Polytope< N, T > &other)
Definition: polytope.hxx:74
virtual bool closed() const
Definition: polytope.hxx:103

© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de)
Heidelberg Collaboratory for Image Processing, University of Heidelberg, Germany

html generated using doxygen and Python
vigra 1.11.0 (Fri May 19 2017)