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; 40 typedef TinyVector<T, N> point_type; 41 typedef TinyVectorView<T, N> point_view_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 */ Polytope()54 Polytope() 55 : graph_() 56 , type_map_(graph_) 57 , vec_map_(graph_) 58 , aligns_map_(graph_) 59 {} 60 61 /** Copy constructor. 62 */ Polytope(const Polytope<N,T> & other)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 */ operator =(const Polytope<N,T> & other)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 */ closed(const node_type n) const90 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 */ closed() const103 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 */ addVertex(const point_view_type & p)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 */ eraseFacet(const node_type u)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 */ getConnected(const node_type u) const159 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 getVertices(const node_type u) const180 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"); 185 ArrayVector<point_view_type> ret; 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 */ litFacets(const point_view_type & p) const196 virtual ArrayVector<node_type> litFacets(const point_view_type & p) const 197 { 198 ArrayVector<node_type> ret; 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 */ tidyUp()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 */ distance(const node_type u,const point_view_type & p) const232 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 */ fill(MultiArrayView<N,unsigned int> & array,const unsigned int label,const point_view_type offset,const point_view_type scale) const247 virtual unsigned int fill( 248 MultiArrayView<N, unsigned int> & array, 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 */ fill(MultiArrayView<N,unsigned int> & array,const unsigned int label,const point_view_type offset) const276 virtual unsigned int fill( 277 MultiArrayView<N, unsigned int> & array, 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 isConnected(const node_type vertex,const node_type facet) const307 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 findNeighbor(const node_type u,const difference_type index) const327 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 assignNeighbors(const node_type u)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 updateInvalidNeighbors(const node_type u)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 openEdge(const node_type u)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 { node_type_iteratorvigra::Polytope::node_type_iterator413 node_type_iterator() 414 {} 415 node_type_iteratorvigra::Polytope::node_type_iterator416 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 operator ++vigra::Polytope::node_type_iterator429 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 operator ==vigra::Polytope::node_type_iterator442 bool operator==(lemon::Invalid i) const 443 { 444 return (static_cast<node_type>(*this) == i); 445 } 446 operator !=vigra::Polytope::node_type_iterator447 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 { out_skip_iteratorvigra::Polytope::out_skip_iterator458 out_skip_iterator() 459 {} 460 out_skip_iteratorvigra::Polytope::out_skip_iterator461 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 operator ++vigra::Polytope::out_skip_iterator476 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 operator ==vigra::Polytope::out_skip_iterator487 bool operator==(lemon::Invalid i) const 488 { 489 return (static_cast<arc_type>(*this) == i); 490 } 491 operator !=vigra::Polytope::out_skip_iterator492 bool operator!=(lemon::Invalid i) const 493 { 494 return (static_cast<arc_type>(*this) != i); 495 } 496 indexvigra::Polytope::out_skip_iterator497 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; 525 typedef typename base_type::point_view_type point_view_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 */ StarPolytope()550 StarPolytope() 551 : base_type() 552 , center_(point_type()) 553 {} 554 555 /** Copy constructor. 556 */ StarPolytope(const point_view_type & center)557 StarPolytope(const point_view_type & center) 558 : base_type() 559 , center_(center) 560 {} 561 562 /** Constructor for the 2-dimensional case taking three vertices and the 563 center. 564 */ StarPolytope(const point_view_type & a,const point_view_type & b,const point_view_type & c,const point_view_type & center)565 StarPolytope( 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 */ 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)587 StarPolytope( 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 */ getCenter() const611 virtual point_type getCenter() const 612 { 613 return center_; 614 } 615 assignNormal(const node_type & u)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 */ addFacet(const node_type & a,const node_type & b)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 */ addFacet(const node_type & a,const node_type & b,const node_type & c)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 close()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 contains(const node_type & n,const point_view_type & p) const760 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 */ contains(const point_view_type & p) const799 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 nVolume(const node_type & n) const811 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 */ nVolume() const830 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 nSurface(const node_type & n) const840 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 */ nSurface() const863 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 886 typedef StarPolytope<N, T> base_type; 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; 890 typedef typename base_type::point_view_type point_view_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 ConvexPolytope()913 ConvexPolytope() 914 : base_type() 915 {} 916 ConvexPolytope(const point_view_type & center)917 ConvexPolytope(const point_view_type & center) 918 : base_type(center) 919 {} 920 ConvexPolytope(const point_view_type & a,const point_view_type & b,const point_view_type & c)921 ConvexPolytope( 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 ConvexPolytope(const point_view_type & a,const point_view_type & b,const point_view_type & c,const point_view_type & d)928 ConvexPolytope( 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 closeFacet(const node_type & vertex,const node_type & facet)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 contains(const node_type & n,const point_view_type & p) const985 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 contains(const point_view_type & p) const997 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 */ addExtremeVertex(const point_view_type & p)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 */ 1050