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

polytope.hxx
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
20namespace 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*/
27template <unsigned int N, class T>
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*/
516template <unsigned int N, class T>
517class StarPolytope : public Polytope<N, T>
518{
519 public:
520
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*/
881template <unsigned int N, class T>
882class 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 */
size_type size() const
Definition array_vector.hxx:358
Definition array_vector.hxx:514
Definition polytope.hxx:883
virtual bool contains(const point_view_type &p) const
Definition polytope.hxx:997
virtual void addExtremeVertex(const point_view_type &p)
Definition polytope.hxx:1012
Base class for, and view to, vigra::MultiArray.
Definition multi_fwd.hxx:127
iterator end()
Definition multi_array.hxx:1939
iterator begin()
Definition multi_array.hxx:1923
Main MultiArray class containing the memory management.
Definition multi_fwd.hxx:131
Represent an n-dimensional polytope.
Definition polytope.hxx:29
Polytope()
Definition polytope.hxx:54
virtual real_type distance(const node_type u, const point_view_type &p) const
Definition polytope.hxx:232
virtual void eraseFacet(const node_type u)
Definition polytope.hxx:132
virtual void operator=(const Polytope< N, T > &other)
Definition polytope.hxx:74
virtual bool closed(const node_type n) const
Definition polytope.hxx:90
virtual void tidyUp()
Definition polytope.hxx:211
virtual node_type addVertex(const point_view_type &p)
Definition polytope.hxx:118
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
virtual ArrayVector< node_type > litFacets(const point_view_type &p) const
Definition polytope.hxx:196
virtual std::set< node_type > getConnected(const node_type u) const
Definition polytope.hxx:159
Polytope(const Polytope< N, T > &other)
Definition polytope.hxx:63
virtual bool closed() const
Definition polytope.hxx:103
virtual unsigned int fill(MultiArrayView< N, unsigned int > &array, const unsigned int label, const point_view_type offset) const
Definition polytope.hxx:276
Specialization of the polytope to polytopes which forms a star domain.
Definition polytope.hxx:518
StarPolytope()
Definition polytope.hxx:550
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
virtual point_type getCenter() const
Definition polytope.hxx:611
StarPolytope(const point_view_type &center)
Definition polytope.hxx:557
virtual bool contains(const point_view_type &p) const
Definition polytope.hxx:799
virtual real_type nSurface() const
Definition polytope.hxx:863
virtual real_type nVolume() const
Definition polytope.hxx:830
virtual node_type addFacet(const node_type &a, const node_type &b, const node_type &c)
Definition polytope.hxx:682
virtual node_type addFacet(const node_type &a, const node_type &b)
Definition polytope.hxx:653
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
size_type size() const
Definition tinyvector.hxx:913
Wrapper for fixed size vectors.
Definition tinyvector.hxx:1254
Class for fixed size vectors.
Definition tinyvector.hxx:1008
NumericTraits< T >::Promote determinant(MultiArrayView< 2, T, C1 > const &a, std::string method="default")
Definition linear_solve.hxx:859
FFTWComplex< R >::NormType abs(const FFTWComplex< R > &a)
absolute value (= magnitude)
Definition fftw3.hxx:1002
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

© 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.12.2