1 /************************************************************************/ 2 /* */ 3 /* Copyright 2014 by Thorsten Beier and Ullrich Koethe */ 4 /* */ 5 /* This file is part of the VIGRA computer vision library. */ 6 /* The VIGRA Website is */ 7 /* http://hci.iwr.uni-heidelberg.de/vigra/ */ 8 /* Please direct questions, bug reports, and contributions to */ 9 /* ullrich.koethe@iwr.uni-heidelberg.de or */ 10 /* vigra@informatik.uni-hamburg.de */ 11 /* */ 12 /* Permission is hereby granted, free of charge, to any person */ 13 /* obtaining a copy of this software and associated documentation */ 14 /* files (the "Software"), to deal in the Software without */ 15 /* restriction, including without limitation the rights to use, */ 16 /* copy, modify, merge, publish, distribute, sublicense, and/or */ 17 /* sell copies of the Software, and to permit persons to whom the */ 18 /* Software is furnished to do so, subject to the following */ 19 /* conditions: */ 20 /* */ 21 /* The above copyright notice and this permission notice shall be */ 22 /* included in all copies or substantial portions of the */ 23 /* Software. */ 24 /* */ 25 /* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND */ 26 /* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES */ 27 /* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND */ 28 /* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT */ 29 /* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, */ 30 /* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING */ 31 /* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR */ 32 /* OTHER DEALINGS IN THE SOFTWARE. */ 33 /* */ 34 /************************************************************************/ 35 36 /** 37 * This header provides definitions of graph-related algorithms 38 */ 39 40 #ifndef VIGRA_GRAPH_ALGORITHMS_HXX 41 #define VIGRA_GRAPH_ALGORITHMS_HXX 42 43 /*std*/ 44 #include <algorithm> 45 #include <vector> 46 #include <functional> 47 #include <set> 48 #include <iomanip> 49 50 /*vigra*/ 51 #include "graphs.hxx" 52 #include "graph_generalization.hxx" 53 #include "multi_gridgraph.hxx" 54 #include "priority_queue.hxx" 55 #include "union_find.hxx" 56 #include "adjacency_list_graph.hxx" 57 #include "graph_maps.hxx" 58 59 #include "timing.hxx" 60 //#include "openmp_helper.hxx" 61 62 63 #include "functorexpression.hxx" 64 #include "array_vector.hxx" 65 66 namespace vigra{ 67 68 /** \addtogroup GraphDataStructures 69 */ 70 //@{ 71 72 namespace detail_graph_algorithms{ 73 template <class GRAPH_MAP,class COMPERATOR> 74 struct GraphItemCompare 75 { 76 GraphItemComparevigra::detail_graph_algorithms::GraphItemCompare77 GraphItemCompare(const GRAPH_MAP & map,const COMPERATOR & comperator) 78 : map_(map), 79 comperator_(comperator){ 80 81 } 82 83 template<class KEY> operator ()vigra::detail_graph_algorithms::GraphItemCompare84 bool operator()(const KEY & a, const KEY & b) const{ 85 return comperator_(map_[a],map_[b]); 86 } 87 88 const GRAPH_MAP & map_; 89 const COMPERATOR & comperator_; 90 }; 91 } // namespace detail_graph_algorithms 92 93 /// \brief get a vector of Edge descriptors 94 /// 95 /// Sort the Edge descriptors given weights 96 /// and a comperator 97 template<class GRAPH,class WEIGHTS,class COMPERATOR> edgeSort(const GRAPH & g,const WEIGHTS & weights,const COMPERATOR & comperator,std::vector<typename GRAPH::Edge> & sortedEdges)98 void edgeSort( 99 const GRAPH & g, 100 const WEIGHTS & weights, 101 const COMPERATOR & comperator, 102 std::vector<typename GRAPH::Edge> & sortedEdges 103 ){ 104 sortedEdges.resize(g.edgeNum()); 105 size_t c=0; 106 for(typename GRAPH::EdgeIt e(g);e!=lemon::INVALID;++e){ 107 sortedEdges[c]=*e; 108 ++c; 109 } 110 detail_graph_algorithms::GraphItemCompare<WEIGHTS,COMPERATOR> edgeComperator(weights,comperator); 111 std::sort(sortedEdges.begin(),sortedEdges.end(),edgeComperator); 112 } 113 114 115 /// \brief copy a lemon node map 116 template<class G,class A,class B> copyNodeMap(const G & g,const A & a,B & b)117 void copyNodeMap(const G & g,const A & a ,B & b){ 118 typename G::NodeIt iter(g); 119 while(iter!=lemon::INVALID){ 120 b[*iter]=a[*iter]; 121 ++iter; 122 } 123 124 } 125 /// \brief copy a lemon edge map 126 template<class G,class A,class B> copyEdgeMap(const G & g,const A & a,B & b)127 void copyEdgeMap(const G & g,const A & a ,B & b){ 128 typename G::EdgeIt iter(g); 129 while(iter!=lemon::INVALID){ 130 b[*iter]=a[*iter]; 131 ++iter; 132 } 133 } 134 /// \brief fill a lemon node map 135 template<class G,class A,class T> fillNodeMap(const G & g,A & a,const T & value)136 void fillNodeMap(const G & g, A & a, const T & value){ 137 typename G::NodeIt iter(g); 138 while(iter!=lemon::INVALID){ 139 a[*iter]=value; 140 ++iter; 141 } 142 } 143 /// \brief fill a lemon edge map 144 template<class G,class A,class T> fillEdgeMap(const G & g,A & a,const T & value)145 void fillEdgeMap(const G & g,A & a ,const T & value){ 146 typename G::EdgeIt iter(g); 147 while(iter!=lemon::INVALID){ 148 a[*iter]=value; 149 ++iter; 150 } 151 } 152 153 /// \brief make a region adjacency graph from a graph and labels w.r.t. that graph 154 /// 155 /// \param graphIn : input graph 156 /// \param labels : labels w.r.t. graphIn 157 /// \param[out] rag : region adjacency graph 158 /// \param[out] affiliatedEdges : a vector of edges of graphIn for each edge in rag 159 /// \param ignoreLabel : optional label to ignore (default: -1 means no label will be ignored) 160 /// 161 template< 162 class GRAPH_IN, 163 class GRAPH_IN_NODE_LABEL_MAP 164 > makeRegionAdjacencyGraph(GRAPH_IN graphIn,GRAPH_IN_NODE_LABEL_MAP labels,AdjacencyListGraph & rag,typename AdjacencyListGraph::template EdgeMap<std::vector<typename GRAPH_IN::Edge>> & affiliatedEdges,const Int64 ignoreLabel=-1)165 void makeRegionAdjacencyGraph( 166 GRAPH_IN graphIn, 167 GRAPH_IN_NODE_LABEL_MAP labels, 168 AdjacencyListGraph & rag, 169 typename AdjacencyListGraph:: template EdgeMap< std::vector<typename GRAPH_IN::Edge> > & affiliatedEdges, 170 const Int64 ignoreLabel=-1 171 ){ 172 rag=AdjacencyListGraph(); 173 typedef typename GraphMapTypeTraits<GRAPH_IN_NODE_LABEL_MAP>::Value LabelType; 174 typedef GRAPH_IN GraphIn; 175 typedef AdjacencyListGraph GraphOut; 176 177 typedef typename GraphIn::Edge EdgeGraphIn; 178 typedef typename GraphIn::NodeIt NodeItGraphIn; 179 typedef typename GraphIn::EdgeIt EdgeItGraphIn; 180 typedef typename GraphOut::Edge EdgeGraphOut; 181 182 183 for(NodeItGraphIn iter(graphIn);iter!=lemon::INVALID;++iter){ 184 const LabelType l=labels[*iter]; 185 if(ignoreLabel==-1 || static_cast<Int64>(l)!=ignoreLabel) 186 rag.addNode(l); 187 } 188 189 for(EdgeItGraphIn e(graphIn);e!=lemon::INVALID;++e){ 190 const EdgeGraphIn edge(*e); 191 const LabelType lu = labels[graphIn.u(edge)]; 192 const LabelType lv = labels[graphIn.v(edge)]; 193 if( lu!=lv && ( ignoreLabel==-1 || (static_cast<Int64>(lu)!=ignoreLabel && static_cast<Int64>(lv)!=ignoreLabel) ) ){ 194 // if there is an edge between lu and lv no new edge will be added 195 rag.addEdge( rag.nodeFromId(lu),rag.nodeFromId(lv)); 196 } 197 } 198 199 //SET UP HYPEREDGES 200 affiliatedEdges.assign(rag); 201 for(EdgeItGraphIn e(graphIn);e!=lemon::INVALID;++e){ 202 const EdgeGraphIn edge(*e); 203 const LabelType lu = labels[graphIn.u(edge)]; 204 const LabelType lv = labels[graphIn.v(edge)]; 205 //std::cout<<"edge between ?? "<<lu<<" "<<lv<<"\n"; 206 if( lu!=lv && ( ignoreLabel==-1 || (static_cast<Int64>(lu)!=ignoreLabel && static_cast<Int64>(lv)!=ignoreLabel) ) ){ 207 //std::cout<<"find edge between "<<lu<<" "<<lv<<"\n"; 208 EdgeGraphOut ragEdge= rag.findEdge(rag.nodeFromId(lu),rag.nodeFromId(lv)); 209 //std::cout<<"invalid?"<<bool(ragEdge==lemon::INVALID)<<" id "<<rag.id(ragEdge)<<"\n"; 210 affiliatedEdges[ragEdge].push_back(edge); 211 //std::cout<<"write done\n"; 212 } 213 } 214 } 215 216 template<unsigned int DIM, class DTAG, class AFF_EDGES> affiliatedEdgesSerializationSize(const GridGraph<DIM,DTAG> &,const AdjacencyListGraph & rag,const AFF_EDGES & affEdges)217 size_t affiliatedEdgesSerializationSize( 218 const GridGraph<DIM,DTAG> &, 219 const AdjacencyListGraph & rag, 220 const AFF_EDGES & affEdges 221 ){ 222 size_t size = 0; 223 224 typedef typename AdjacencyListGraph::EdgeIt EdgeIt; 225 typedef typename AdjacencyListGraph::Edge Edge; 226 227 for(EdgeIt iter(rag); iter!=lemon::INVALID; ++iter){ 228 const Edge e(*iter); 229 size+=1; 230 size+=affEdges[e].size()*(DIM+1); 231 } 232 return size; 233 } 234 235 template<class OUT_ITER, unsigned int DIM, class DTAG, class AFF_EDGES> serializeAffiliatedEdges(const GridGraph<DIM,DTAG> &,const AdjacencyListGraph & rag,const AFF_EDGES & affEdges,OUT_ITER outIter)236 void serializeAffiliatedEdges( 237 const GridGraph<DIM,DTAG> &, 238 const AdjacencyListGraph & rag, 239 const AFF_EDGES & affEdges, 240 OUT_ITER outIter 241 ){ 242 243 typedef typename AdjacencyListGraph::EdgeIt EdgeIt; 244 typedef typename AdjacencyListGraph::Edge Edge; 245 typedef typename GridGraph<DIM,DTAG>::Edge GEdge; 246 247 for(EdgeIt iter(rag); iter!=lemon::INVALID; ++iter){ 248 249 const Edge edge = *iter; 250 const size_t numAffEdge = affEdges[edge].size(); 251 *outIter = numAffEdge; ++outIter; 252 253 for(size_t i=0; i<numAffEdge; ++i){ 254 const GEdge gEdge = affEdges[edge][i]; 255 for(size_t d=0; d<DIM+1; ++d){ 256 *outIter = gEdge[d]; ++outIter; 257 } 258 } 259 } 260 } 261 262 template<class IN_ITER, unsigned int DIM, class DTAG, class AFF_EDGES> deserializeAffiliatedEdges(const GridGraph<DIM,DTAG> &,const AdjacencyListGraph & rag,AFF_EDGES & affEdges,IN_ITER begin,IN_ITER)263 void deserializeAffiliatedEdges( 264 const GridGraph<DIM,DTAG> &, 265 const AdjacencyListGraph & rag, 266 AFF_EDGES & affEdges, 267 IN_ITER begin, 268 IN_ITER 269 ){ 270 271 typedef typename AdjacencyListGraph::EdgeIt EdgeIt; 272 typedef typename AdjacencyListGraph::Edge Edge; 273 typedef typename GridGraph<DIM,DTAG>::Edge GEdge; 274 275 affEdges.assign(rag); 276 277 for(EdgeIt iter(rag); iter!=lemon::INVALID; ++iter){ 278 279 const Edge edge = *iter; 280 const size_t numAffEdge = *begin; ++begin; 281 282 for(size_t i=0; i<numAffEdge; ++i){ 283 GEdge gEdge; 284 for(size_t d=0; d<DIM+1; ++d){ 285 gEdge[d]=*begin; ++begin; 286 } 287 affEdges[edge].push_back(gEdge); 288 } 289 } 290 } 291 292 293 294 295 /// \brief shortest path computer 296 template<class GRAPH,class WEIGHT_TYPE> 297 class ShortestPathDijkstra{ 298 public: 299 typedef GRAPH Graph; 300 301 typedef typename Graph::Node Node; 302 typedef typename Graph::NodeIt NodeIt; 303 typedef typename Graph::Edge Edge; 304 typedef typename Graph::OutArcIt OutArcIt; 305 306 typedef WEIGHT_TYPE WeightType; 307 typedef ChangeablePriorityQueue<WeightType> PqType; 308 typedef typename Graph:: template NodeMap<Node> PredecessorsMap; 309 typedef typename Graph:: template NodeMap<WeightType> DistanceMap; 310 typedef ArrayVector<Node> DiscoveryOrder; 311 312 /// \ brief constructor from graph ShortestPathDijkstra(const Graph & g)313 ShortestPathDijkstra(const Graph & g) 314 : graph_(g), 315 pq_(g.maxNodeId()+1), 316 predMap_(g), 317 distMap_(g) 318 { 319 } 320 321 /// \brief run shortest path given edge weights 322 /// 323 /// \param weights : edge weights encoding the distance between adjacent nodes (must be non-negative) 324 /// \param source : source node where shortest path should start 325 /// \param target : target node where shortest path should stop. If target is not given 326 /// or <tt>INVALID</tt>, the shortest path from source to all reachable nodes is computed 327 /// \param maxDistance : path search is terminated when the path length exceeds <tt>maxDistance</tt> 328 /// 329 /// When a valid \a target is unreachable from \a source (either because the graph is disconnected 330 /// or \a maxDistance is exceeded), it is set to <tt>lemon::INVALID</tt>. In contrast, if \a target 331 /// was <tt>lemon::INVALID</tt> at the beginning, it will always be set to the last node 332 /// visited in the search. 333 template<class WEIGHTS> run(const WEIGHTS & weights,const Node & source,const Node & target=lemon::INVALID,WeightType maxDistance=NumericTraits<WeightType>::max ())334 void run(const WEIGHTS & weights, const Node & source, 335 const Node & target = lemon::INVALID, 336 WeightType maxDistance=NumericTraits<WeightType>::max()) 337 { 338 this->initializeMaps(source); 339 runImpl(weights, target, maxDistance); 340 } 341 342 /// \brief run shortest path in a region of interest of a \ref GridGraph. 343 /// 344 /// \param start : first point in the desired ROI. 345 /// \param stop : beyond the last point in the desired ROI (i.e. exclusive) 346 /// \param weights : edge weights encoding the distance between adjacent nodes (must be non-negative) 347 /// \param source : source node where shortest path should start 348 /// \param target : target node where shortest path should stop. If target is not given 349 /// or <tt>INVALID</tt>, the shortest path from source to all reachable nodes is computed 350 /// \param maxDistance : path search is terminated when the path length exceeds <tt>maxDistance</tt> 351 /// 352 /// This version of <tt>run()</tt> restricts the path search to the ROI <tt>[start, stop)</tt> and only 353 /// works for instances of \ref GridGraph. Otherwise, it is identical to the standard <tt>run()</tt> 354 /// function. 355 template<class WEIGHTS> run(Node const & start,Node const & stop,const WEIGHTS & weights,const Node & source,const Node & target=lemon::INVALID,WeightType maxDistance=NumericTraits<WeightType>::max ())356 void run(Node const & start, Node const & stop, 357 const WEIGHTS & weights, const Node & source, 358 const Node & target = lemon::INVALID, 359 WeightType maxDistance=NumericTraits<WeightType>::max()) 360 { 361 vigra_precondition(allLessEqual(start, source) && allLess(source, stop), 362 "ShortestPathDijkstra::run(): source is not within ROI"); 363 vigra_precondition(target == lemon::INVALID || 364 (allLessEqual(start, target) && allLess(target, stop)), 365 "ShortestPathDijkstra::run(): target is not within ROI"); 366 this->initializeMaps(source, start, stop); 367 runImpl(weights, target, maxDistance); 368 } 369 370 /// \brief run shortest path again with given edge weights 371 /// 372 /// This only differs from standard <tt>run()</tt> by initialization: Instead of resetting 373 /// the entire graph, this only resets the nodes that have been visited in the 374 /// previous run, i.e. the contents of the array <tt>discoveryOrder()</tt>. 375 /// This will be much faster if only a small fraction of the nodes has to be reset. 376 template<class WEIGHTS> reRun(const WEIGHTS & weights,const Node & source,const Node & target=lemon::INVALID,WeightType maxDistance=NumericTraits<WeightType>::max ())377 void reRun(const WEIGHTS & weights, const Node & source, 378 const Node & target = lemon::INVALID, 379 WeightType maxDistance=NumericTraits<WeightType>::max()) 380 { 381 this->reInitializeMaps(source); 382 runImpl(weights, target, maxDistance); 383 } 384 385 /// \brief run shortest path with given edge weights from multiple sources. 386 /// 387 /// This is otherwise identical to standard <tt>run()</tt>, except that 388 /// <tt>source()</tt> returns <tt>lemon::INVALID</tt> after path search finishes. 389 template<class WEIGHTS, class ITER> 390 void runMultiSource(const WEIGHTS & weights,ITER source_begin,ITER source_end,const Node & target=lemon::INVALID,WeightType maxDistance=NumericTraits<WeightType>::max ())391 runMultiSource(const WEIGHTS & weights, ITER source_begin, ITER source_end, 392 const Node & target = lemon::INVALID, 393 WeightType maxDistance=NumericTraits<WeightType>::max()) 394 { 395 this->initializeMapsMultiSource(source_begin, source_end); 396 runImpl(weights, target, maxDistance); 397 } 398 399 400 /// \brief run shortest path with given edge weights from multiple sources. 401 /// 402 /// This is otherwise identical to standard <tt>run()</tt>, except that 403 /// <tt>source()</tt> returns <tt>lemon::INVALID</tt> after path search finishes. 404 template<class EFGE_WEIGHTS,class NODE_WEIGHTS, class ITER> 405 void runMultiSource(const EFGE_WEIGHTS & edgeWeights,const NODE_WEIGHTS & nodeWeights,ITER source_begin,ITER source_end,const Node & target=lemon::INVALID,WeightType maxDistance=NumericTraits<WeightType>::max ())406 runMultiSource( 407 const EFGE_WEIGHTS & edgeWeights, 408 const NODE_WEIGHTS & nodeWeights, 409 ITER source_begin, 410 ITER source_end, 411 const Node & target = lemon::INVALID, 412 WeightType maxDistance = NumericTraits<WeightType>::max()) 413 { 414 this->initializeMapsMultiSource(source_begin, source_end); 415 runImplWithNodeWeights(edgeWeights, nodeWeights, target, maxDistance); 416 } 417 418 /// \brief get the graph graph() const419 const Graph & graph()const{ 420 return graph_; 421 } 422 /// \brief get the source node source() const423 const Node & source()const{ 424 return source_; 425 } 426 /// \brief get the target node target() const427 const Node & target()const{ 428 return target_; 429 } 430 431 /// \brief check if explicit target is given hasTarget() const432 bool hasTarget()const{ 433 return target_!=lemon::INVALID; 434 } 435 436 /// \brief get an array with all visited nodes, sorted by distance from source discoveryOrder() const437 const DiscoveryOrder & discoveryOrder() const{ 438 return discoveryOrder_; 439 } 440 441 /// \brief get the predecessors node map (after a call of run) predecessors() const442 const PredecessorsMap & predecessors()const{ 443 return predMap_; 444 } 445 446 /// \brief get the distances node map (after a call of run) distances() const447 const DistanceMap & distances()const{ 448 return distMap_; 449 } 450 451 /// \brief get the distance to a rarget node (after a call of run) distance(const Node & target) const452 WeightType distance(const Node & target)const{ 453 return distMap_[target]; 454 } 455 456 457 private: 458 459 template<class WEIGHTS> runImpl(const WEIGHTS & weights,const Node & target=lemon::INVALID,WeightType maxDistance=NumericTraits<WeightType>::max ())460 void runImpl(const WEIGHTS & weights, 461 const Node & target = lemon::INVALID, 462 WeightType maxDistance=NumericTraits<WeightType>::max()) 463 { 464 ZeroNodeMap<Graph, WEIGHT_TYPE> zeroNodeMap; 465 this->runImplWithNodeWeights(weights,zeroNodeMap, target, maxDistance); 466 } 467 468 469 template<class EDGE_WEIGHTS, class NODE_WEIGHTS> runImplWithNodeWeights(const EDGE_WEIGHTS & edgeWeights,const NODE_WEIGHTS & nodeWeights,const Node & target=lemon::INVALID,WeightType maxDistance=NumericTraits<WeightType>::max ())470 void runImplWithNodeWeights( 471 const EDGE_WEIGHTS & edgeWeights, 472 const NODE_WEIGHTS & nodeWeights, 473 const Node & target = lemon::INVALID, 474 WeightType maxDistance=NumericTraits<WeightType>::max()) 475 { 476 target_ = lemon::INVALID; 477 while(!pq_.empty() ){ //&& !finished){ 478 const Node topNode(graph_.nodeFromId(pq_.top())); 479 if(distMap_[topNode] > maxDistance) 480 break; // distance threshold exceeded 481 pq_.pop(); 482 discoveryOrder_.push_back(topNode); 483 if(topNode == target) 484 break; 485 // loop over all neigbours 486 for(OutArcIt outArcIt(graph_,topNode);outArcIt!=lemon::INVALID;++outArcIt){ 487 const Node otherNode = graph_.target(*outArcIt); 488 const size_t otherNodeId = graph_.id(otherNode); 489 const WeightType otherNodeWeight = nodeWeights[otherNode]; 490 if(pq_.contains(otherNodeId)){ 491 const Edge edge(*outArcIt); 492 const WeightType currentDist = distMap_[otherNode]; 493 const WeightType alternativeDist = distMap_[topNode]+edgeWeights[edge]+otherNodeWeight; 494 if(alternativeDist<currentDist){ 495 pq_.push(otherNodeId,alternativeDist); 496 distMap_[otherNode]=alternativeDist; 497 predMap_[otherNode]=topNode; 498 } 499 } 500 else if(predMap_[otherNode]==lemon::INVALID){ 501 const Edge edge(*outArcIt); 502 const WeightType initialDist = distMap_[topNode]+edgeWeights[edge]+otherNodeWeight; 503 if(initialDist<=maxDistance) 504 { 505 pq_.push(otherNodeId,initialDist); 506 distMap_[otherNode]=initialDist; 507 predMap_[otherNode]=topNode; 508 } 509 } 510 } 511 } 512 while(!pq_.empty() ){ 513 const Node topNode(graph_.nodeFromId(pq_.top())); 514 predMap_[topNode]=lemon::INVALID; 515 pq_.pop(); 516 } 517 if(target == lemon::INVALID || discoveryOrder_.back() == target) 518 target_ = discoveryOrder_.back(); // Means that target was reached. If, to the contrary, target 519 // was unreachable within maxDistance, target_ remains INVALID. 520 } 521 initializeMaps(Node const & source)522 void initializeMaps(Node const & source){ 523 for(NodeIt n(graph_); n!=lemon::INVALID; ++n){ 524 const Node node(*n); 525 predMap_[node]=lemon::INVALID; 526 } 527 distMap_[source]=static_cast<WeightType>(0.0); 528 predMap_[source]=source; 529 discoveryOrder_.clear(); 530 pq_.push(graph_.id(source),0.0); 531 source_=source; 532 } 533 initializeMaps(Node const & source,Node const & start,Node const & stop)534 void initializeMaps(Node const & source, 535 Node const & start, Node const & stop) 536 { 537 Node left_border = min(start, Node(1)), 538 right_border = min(predMap_.shape()-stop, Node(1)), 539 DONT_TOUCH = Node(lemon::INVALID) - Node(1); 540 541 initMultiArrayBorder(predMap_.subarray(start-left_border, stop+right_border), 542 left_border, right_border, DONT_TOUCH); 543 predMap_.subarray(start, stop) = lemon::INVALID; 544 predMap_[source]=source; 545 546 distMap_[source]=static_cast<WeightType>(0.0); 547 discoveryOrder_.clear(); 548 pq_.push(graph_.id(source),0.0); 549 source_=source; 550 } 551 552 template <class ITER> initializeMapsMultiSource(ITER source,ITER source_end)553 void initializeMapsMultiSource(ITER source, ITER source_end){ 554 for(NodeIt n(graph_); n!=lemon::INVALID; ++n){ 555 const Node node(*n); 556 predMap_[node]=lemon::INVALID; 557 } 558 discoveryOrder_.clear(); 559 for( ; source != source_end; ++source) 560 { 561 distMap_[*source]=static_cast<WeightType>(0.0); 562 predMap_[*source]=*source; 563 pq_.push(graph_.id(*source),0.0); 564 } 565 source_=lemon::INVALID; 566 } 567 reInitializeMaps(Node const & source)568 void reInitializeMaps(Node const & source){ 569 for(unsigned int n=0; n<discoveryOrder_.size(); ++n){ 570 predMap_[discoveryOrder_[n]]=lemon::INVALID; 571 } 572 distMap_[source]=static_cast<WeightType>(0.0); 573 predMap_[source]=source; 574 discoveryOrder_.clear(); 575 pq_.push(graph_.id(source),0.0); 576 source_=source; 577 } 578 579 const Graph & graph_; 580 PqType pq_; 581 PredecessorsMap predMap_; 582 DistanceMap distMap_; 583 DiscoveryOrder discoveryOrder_; 584 585 Node source_; 586 Node target_; 587 }; 588 589 /// \brief get the length in node units of a path 590 template<class NODE,class PREDECESSORS> pathLength(const NODE source,const NODE target,const PREDECESSORS & predecessors)591 size_t pathLength( 592 const NODE source, 593 const NODE target, 594 const PREDECESSORS & predecessors 595 ){ 596 if(predecessors[target]==lemon::INVALID) 597 return 0; 598 else{ 599 NODE currentNode = target; 600 size_t length=1; 601 while(currentNode!=source){ 602 currentNode=predecessors[currentNode]; 603 length+=1; 604 } 605 return length; 606 } 607 } 608 609 /// \brief Astar Shortest path search 610 template<class GRAPH,class WEIGHTS,class PREDECESSORS,class DISTANCE,class HEURSTIC> shortestPathAStar(const GRAPH & graph,const typename GRAPH::Node & source,const typename GRAPH::Node & target,const WEIGHTS & weights,PREDECESSORS & predecessors,DISTANCE & distance,const HEURSTIC & heuristic)611 void shortestPathAStar( 612 const GRAPH & graph, 613 const typename GRAPH::Node & source, 614 const typename GRAPH::Node & target, 615 const WEIGHTS & weights, 616 PREDECESSORS & predecessors, 617 DISTANCE & distance, 618 const HEURSTIC & heuristic 619 ){ 620 621 typedef GRAPH Graph; 622 typedef typename Graph::Edge Edge; 623 typedef typename Graph::Node Node; 624 typedef typename Graph::NodeIt NodeIt; 625 typedef typename Graph::OutArcIt OutArcIt; 626 typedef typename DISTANCE::value_type DistanceType; 627 628 typename GRAPH:: template NodeMap<bool> closedSet(graph); 629 vigra::ChangeablePriorityQueue<typename WEIGHTS::value_type> estimatedDistanceOpenSet(graph.maxNodeId()+1); 630 // initialize 631 for(NodeIt n(graph);n!=lemon::INVALID;++n){ 632 const Node node(*n); 633 closedSet[node]=false; 634 distance[node]=std::numeric_limits<DistanceType>::infinity(); 635 predecessors[node]=lemon::INVALID; 636 } 637 // distance and estimated distance for start node 638 distance[source]=static_cast<DistanceType>(0.0); 639 estimatedDistanceOpenSet.push(graph.id(source),heuristic(source,target)); 640 641 // while any nodes left in openSet 642 while(!estimatedDistanceOpenSet.empty()){ 643 644 // get the node with the lpwest estimated distance in the open set 645 const Node current = graph.nodeFromId(estimatedDistanceOpenSet.top()); 646 647 // reached target? 648 if(current==target) 649 break; 650 651 // remove current from openSet 652 // add current to closedSet 653 estimatedDistanceOpenSet.pop(); 654 closedSet[current]=true; 655 656 // iterate over neigbours of current 657 for(OutArcIt outArcIt(graph,current);outArcIt!=lemon::INVALID;++outArcIt){ 658 659 // get neigbour node and id 660 const Node neighbour = graph.target(*outArcIt); 661 const size_t neighbourId = graph.id(neighbour); 662 663 // if neighbour is not yet in closedSet 664 if(!closedSet[neighbour]){ 665 666 // get edge between current and neigbour 667 const Edge edge(*outArcIt); 668 669 // get tentative score 670 const DistanceType tenativeScore = distance[current] + weights[edge]; 671 672 // neighbour NOT in openSet OR tentative score better than the current distance 673 if(!estimatedDistanceOpenSet.contains(neighbourId) || tenativeScore < distance[neighbour] ){ 674 // set predecessors and distance 675 predecessors[neighbour]=current; 676 distance[neighbour]=tenativeScore; 677 678 // update the estimated cost from neighbour to target 679 // ( and neigbour will be (re)-added to openSet) 680 estimatedDistanceOpenSet.push(neighbourId,distance[neighbour]+heuristic(neighbour,target)); 681 } 682 } 683 } 684 } 685 } 686 687 688 template< 689 class GRAPH, 690 class EDGE_WEIGHTS, 691 class NODE_WEIGHTS, 692 class SEED_NODE_MAP, 693 class WEIGHT_TYPE 694 > shortestPathSegmentation(const GRAPH & graph,const EDGE_WEIGHTS & edgeWeights,const NODE_WEIGHTS & nodeWeights,SEED_NODE_MAP & seeds)695 void shortestPathSegmentation( 696 const GRAPH & graph, 697 const EDGE_WEIGHTS & edgeWeights, 698 const NODE_WEIGHTS & nodeWeights, 699 SEED_NODE_MAP & seeds 700 ){ 701 702 typedef GRAPH Graph; 703 typedef typename Graph::Node Node; 704 typedef typename Graph::NodeIt NodeIt; 705 typedef WEIGHT_TYPE WeightType; 706 707 // find seeds 708 std::vector<Node> seededNodes; 709 for(NodeIt n(graph);n!=lemon::INVALID;++n){ 710 const Node node(*n); 711 // not a seed 712 if(seeds[node]!=0){ 713 seededNodes.push_back(node); 714 } 715 } 716 717 // do shortest path 718 typedef ShortestPathDijkstra<Graph, WeightType> Sp; 719 typedef typename Sp::PredecessorsMap PredecessorsMap; 720 Sp sp(graph); 721 sp.runMultiSource(edgeWeights, nodeWeights, seededNodes.begin(), seededNodes.end()); 722 const PredecessorsMap & predMap = sp.predecessors(); 723 // do the labeling 724 for(NodeIt n(graph);n!=lemon::INVALID;++n){ 725 Node node(*n); 726 if(seeds[node]==0){ 727 Node pred=predMap[node]; 728 while(seeds[pred]==0){ 729 pred=predMap[pred]; 730 } 731 seeds[node]=seeds[pred]; 732 } 733 } 734 } 735 736 namespace detail_watersheds_segmentation{ 737 738 struct RawPriorityFunctor{ 739 template<class LabelType, class T> operator ()vigra::detail_watersheds_segmentation::RawPriorityFunctor740 T operator()(const LabelType /*label*/,const T priority)const{ 741 return priority; 742 } 743 }; 744 745 746 747 template<class PRIORITY_TYPE,class LABEL_TYPE> 748 struct CarvingFunctor{ CarvingFunctorvigra::detail_watersheds_segmentation::CarvingFunctor749 CarvingFunctor(const LABEL_TYPE backgroundLabel, 750 const PRIORITY_TYPE & factor, 751 const PRIORITY_TYPE & noPriorBelow 752 ) 753 : backgroundLabel_(backgroundLabel), 754 factor_(factor), 755 noPriorBelow_(noPriorBelow){ 756 } operator ()vigra::detail_watersheds_segmentation::CarvingFunctor757 PRIORITY_TYPE operator()(const LABEL_TYPE label,const PRIORITY_TYPE priority)const{ 758 if(priority>=noPriorBelow_) 759 return (label==backgroundLabel_ ? priority*factor_ : priority); 760 else{ 761 return priority; 762 } 763 } 764 LABEL_TYPE backgroundLabel_; 765 PRIORITY_TYPE factor_; 766 PRIORITY_TYPE noPriorBelow_; 767 }; 768 769 770 template< 771 class GRAPH, 772 class EDGE_WEIGHTS, 773 class SEEDS, 774 class PRIORITY_MANIP_FUNCTOR, 775 class LABELS 776 > edgeWeightedWatershedsSegmentationImpl(const GRAPH & g,const EDGE_WEIGHTS & edgeWeights,const SEEDS & seeds,PRIORITY_MANIP_FUNCTOR & priorManipFunctor,LABELS & labels)777 void edgeWeightedWatershedsSegmentationImpl( 778 const GRAPH & g, 779 const EDGE_WEIGHTS & edgeWeights, 780 const SEEDS & seeds, 781 PRIORITY_MANIP_FUNCTOR & priorManipFunctor, 782 LABELS & labels 783 ){ 784 typedef GRAPH Graph; 785 typedef typename Graph::Edge Edge; 786 typedef typename Graph::Node Node; 787 typedef typename Graph::NodeIt NodeIt; 788 typedef typename Graph::OutArcIt OutArcIt; 789 790 typedef typename EDGE_WEIGHTS::Value WeightType; 791 typedef typename LABELS::Value LabelType; 792 //typedef typename Graph:: template EdgeMap<bool> EdgeBoolMap; 793 typedef PriorityQueue<Edge,WeightType,true> PQ; 794 795 PQ pq; 796 //EdgeBoolMap inPQ(g); 797 copyNodeMap(g,seeds,labels); 798 //fillEdgeMap(g,inPQ,false); 799 800 801 // put edges from nodes with seed on pq 802 for(NodeIt n(g);n!=lemon::INVALID;++n){ 803 const Node node(*n); 804 if(labels[node]!=static_cast<LabelType>(0)){ 805 for(OutArcIt a(g,node);a!=lemon::INVALID;++a){ 806 const Edge edge(*a); 807 const Node neigbour=g.target(*a); 808 //std::cout<<"n- node "<<g.id(neigbour)<<"\n"; 809 if(labels[neigbour]==static_cast<LabelType>(0)){ 810 const WeightType priority = priorManipFunctor(labels[node],edgeWeights[edge]); 811 pq.push(edge,priority); 812 //inPQ[edge]=true; 813 } 814 } 815 } 816 } 817 818 819 while(!pq.empty()){ 820 821 const Edge edge = pq.top(); 822 pq.pop(); 823 824 const Node u = g.u(edge); 825 const Node v = g.v(edge); 826 const LabelType lU = labels[u]; 827 const LabelType lV = labels[v]; 828 829 830 if(lU==0 && lV==0){ 831 throw std::runtime_error("both have no labels"); 832 } 833 else if(lU!=0 && lV!=0){ 834 // nothing to do 835 } 836 else{ 837 838 const Node unlabeledNode = lU==0 ? u : v; 839 const LabelType label = lU==0 ? lV : lU; 840 841 // assign label to unlabeled node 842 labels[unlabeledNode] = label; 843 844 // iterate over the nodes edges 845 for(OutArcIt a(g,unlabeledNode);a!=lemon::INVALID;++a){ 846 const Edge otherEdge(*a); 847 const Node targetNode=g.target(*a); 848 if(labels[targetNode] == 0){ 849 //if(inPQ[otherEdge] == false && labels[targetNode] == 0){ 850 const WeightType priority = priorManipFunctor(label,edgeWeights[otherEdge]); 851 pq.push(otherEdge,priority); 852 // inPQ[otherEdge]=true; 853 } 854 } 855 } 856 857 } 858 859 } 860 861 } // end namespace detail_watersheds_segmentation 862 863 864 /// \brief edge weighted watersheds Segmentataion 865 /// 866 /// \param g: input graph 867 /// \param edgeWeights : edge weights / edge indicator 868 /// \param seeds : seed must be non empty! 869 /// \param[out] labels : resulting nodeLabeling (not necessarily dense) 870 template<class GRAPH,class EDGE_WEIGHTS,class SEEDS,class LABELS> edgeWeightedWatershedsSegmentation(const GRAPH & g,const EDGE_WEIGHTS & edgeWeights,const SEEDS & seeds,LABELS & labels)871 void edgeWeightedWatershedsSegmentation( 872 const GRAPH & g, 873 const EDGE_WEIGHTS & edgeWeights, 874 const SEEDS & seeds, 875 LABELS & labels 876 ){ 877 detail_watersheds_segmentation::RawPriorityFunctor fPriority; 878 detail_watersheds_segmentation::edgeWeightedWatershedsSegmentationImpl(g,edgeWeights,seeds,fPriority,labels); 879 } 880 881 882 /// \brief edge weighted watersheds Segmentataion 883 /// 884 /// \param g: input graph 885 /// \param edgeWeights : edge weights / edge indicator 886 /// \param seeds : seed must be non empty! 887 /// \param backgroundLabel : which label is background 888 /// \param backgroundBias : bias for background 889 /// \param noPriorBelow : don't bias the background if edge indicator is below this value 890 /// \param[out] labels : resulting nodeLabeling (not necessarily dense) 891 template<class GRAPH,class EDGE_WEIGHTS,class SEEDS,class LABELS> carvingSegmentation(const GRAPH & g,const EDGE_WEIGHTS & edgeWeights,const SEEDS & seeds,const typename LABELS::Value backgroundLabel,const typename EDGE_WEIGHTS::Value backgroundBias,const typename EDGE_WEIGHTS::Value noPriorBelow,LABELS & labels)892 void carvingSegmentation( 893 const GRAPH & g, 894 const EDGE_WEIGHTS & edgeWeights, 895 const SEEDS & seeds, 896 const typename LABELS::Value backgroundLabel, 897 const typename EDGE_WEIGHTS::Value backgroundBias, 898 const typename EDGE_WEIGHTS::Value noPriorBelow, 899 LABELS & labels 900 ){ 901 typedef typename EDGE_WEIGHTS::Value WeightType; 902 typedef typename LABELS::Value LabelType; 903 detail_watersheds_segmentation::CarvingFunctor<WeightType,LabelType> fPriority(backgroundLabel,backgroundBias, noPriorBelow); 904 detail_watersheds_segmentation::edgeWeightedWatershedsSegmentationImpl(g,edgeWeights,seeds,fPriority,labels); 905 } 906 907 /// \brief edge weighted watersheds Segmentataion 908 /// 909 /// \param graph: input graph 910 /// \param edgeWeights : edge weights / edge indicator 911 /// \param nodeSizes : size of each node 912 /// \param k : free parameter of felzenszwalb algorithm 913 /// \param[out] nodeLabeling : nodeLabeling (not necessarily dense) 914 /// \param nodeNumStopCond : optional stopping condition 915 template< class GRAPH , class EDGE_WEIGHTS, class NODE_SIZE,class NODE_LABEL_MAP> felzenszwalbSegmentation(const GRAPH & graph,const EDGE_WEIGHTS & edgeWeights,const NODE_SIZE & nodeSizes,float k,NODE_LABEL_MAP & nodeLabeling,const int nodeNumStopCond=-1)916 void felzenszwalbSegmentation( 917 const GRAPH & graph, 918 const EDGE_WEIGHTS & edgeWeights, 919 const NODE_SIZE & nodeSizes, 920 float k, 921 NODE_LABEL_MAP & nodeLabeling, 922 const int nodeNumStopCond = -1 923 ){ 924 typedef GRAPH Graph; 925 typedef typename Graph::Edge Edge; 926 typedef typename Graph::Node Node; 927 928 typedef typename EDGE_WEIGHTS::Value WeightType; 929 typedef typename EDGE_WEIGHTS::Value NodeSizeType; 930 typedef typename Graph:: template NodeMap<WeightType> NodeIntDiffMap; 931 typedef typename Graph:: template NodeMap<NodeSizeType> NodeSizeAccMap; 932 933 // initalize node size map and internal diff map 934 NodeIntDiffMap internalDiff(graph); 935 NodeSizeAccMap nodeSizeAcc(graph); 936 copyNodeMap(graph,nodeSizes,nodeSizeAcc); 937 fillNodeMap(graph,internalDiff,static_cast<WeightType>(0.0)); 938 939 940 941 // initlaize internal node diff map 942 943 // sort the edges by their weights 944 std::vector<Edge> sortedEdges; 945 std::less<WeightType> comperator; 946 edgeSort(graph,edgeWeights,comperator,sortedEdges); 947 948 // make the ufd 949 UnionFindArray<UInt64> ufdArray(graph.maxNodeId()+1); 950 951 952 size_t nodeNum = graph.nodeNum(); 953 954 955 while(true){ 956 // iterate over edges is the sorted order 957 for(size_t i=0;i<sortedEdges.size();++i){ 958 const Edge e = sortedEdges[i]; 959 const size_t rui = ufdArray.findIndex(graph.id(graph.u(e))); 960 const size_t rvi = ufdArray.findIndex(graph.id(graph.v(e))); 961 const Node ru = graph.nodeFromId(rui); 962 const Node rv = graph.nodeFromId(rvi); 963 if(rui!=rvi){ 964 965 //check if to merge or not ? 966 const WeightType w = edgeWeights[e]; 967 const NodeSizeType sizeRu = nodeSizeAcc[ru]; 968 const NodeSizeType sizeRv = nodeSizeAcc[rv]; 969 const WeightType tauRu = static_cast<WeightType>(k)/static_cast<WeightType>(sizeRu); 970 const WeightType tauRv = static_cast<WeightType>(k)/static_cast<WeightType>(sizeRv); 971 const WeightType minIntDiff = std::min(internalDiff[ru]+tauRu,internalDiff[rv]+tauRv); 972 if(w<=minIntDiff){ 973 // do merge 974 ufdArray.makeUnion(rui,rvi); 975 --nodeNum; 976 // update size and internal difference 977 const size_t newRepId = ufdArray.findIndex(rui); 978 const Node newRepNode = graph.nodeFromId(newRepId); 979 internalDiff[newRepNode]=w; 980 nodeSizeAcc[newRepNode] = sizeRu+sizeRv; 981 } 982 } 983 if(nodeNumStopCond >= 0 && nodeNum==static_cast<size_t>(nodeNumStopCond)){ 984 break; 985 } 986 } 987 if(nodeNumStopCond==-1){ 988 break; 989 } 990 else{ 991 if(nodeNumStopCond >= 0 && nodeNum>static_cast<size_t>(nodeNumStopCond)){ 992 k *= 1.2f; 993 } 994 else{ 995 break; 996 } 997 } 998 } 999 ufdArray.makeContiguous(); 1000 for(typename GRAPH::NodeIt n(graph);n!=lemon::INVALID;++n){ 1001 const Node node(*n); 1002 nodeLabeling[node]=ufdArray.findLabel(graph.id(node)); 1003 } 1004 } 1005 1006 1007 1008 1009 namespace detail_graph_smoothing{ 1010 1011 template< 1012 class GRAPH, 1013 class NODE_FEATURES_IN, 1014 class EDGE_WEIGHTS, 1015 class WEIGHTS_TO_SMOOTH_FACTOR, 1016 class NODE_FEATURES_OUT 1017 > graphSmoothingImpl(const GRAPH & g,const NODE_FEATURES_IN & nodeFeaturesIn,const EDGE_WEIGHTS & edgeWeights,WEIGHTS_TO_SMOOTH_FACTOR & weightsToSmoothFactor,NODE_FEATURES_OUT & nodeFeaturesOut)1018 void graphSmoothingImpl( 1019 const GRAPH & g, 1020 const NODE_FEATURES_IN & nodeFeaturesIn, 1021 const EDGE_WEIGHTS & edgeWeights, 1022 WEIGHTS_TO_SMOOTH_FACTOR & weightsToSmoothFactor, 1023 NODE_FEATURES_OUT & nodeFeaturesOut 1024 1025 ){ 1026 typedef GRAPH Graph; 1027 typedef typename Graph::Edge Edge; 1028 typedef typename Graph::Node Node; 1029 typedef typename Graph::NodeIt NodeIt; 1030 typedef typename Graph::OutArcIt OutArcIt; 1031 1032 typedef typename NODE_FEATURES_IN::Value NodeFeatureInValue; 1033 typedef typename NODE_FEATURES_OUT::Reference NodeFeatureOutRef; 1034 typedef typename EDGE_WEIGHTS::ConstReference SmoothFactorType; 1035 1036 1037 //fillNodeMap(g, nodeFeaturesOut, typename NODE_FEATURES_OUT::value_type(0.0)); 1038 1039 for(NodeIt n(g);n!=lemon::INVALID;++n){ 1040 1041 const Node node(*n); 1042 1043 NodeFeatureInValue featIn = nodeFeaturesIn[node]; 1044 NodeFeatureOutRef featOut = nodeFeaturesOut[node]; 1045 1046 featOut=0; 1047 float weightSum = 0.0; 1048 size_t degree = 0; 1049 for(OutArcIt a(g,node);a!=lemon::INVALID;++a){ 1050 const Edge edge(*a); 1051 const Node neigbour(g.target(*a)); 1052 SmoothFactorType smoothFactor= weightsToSmoothFactor(edgeWeights[edge]); 1053 1054 NodeFeatureInValue neighbourFeat = nodeFeaturesIn[neigbour]; 1055 neighbourFeat*=smoothFactor; 1056 if(degree==0) 1057 featOut = neighbourFeat; 1058 else 1059 featOut += neighbourFeat; 1060 weightSum+=smoothFactor; 1061 ++degree; 1062 } 1063 // fixme..set me to right type 1064 featIn*=static_cast<float>(degree); 1065 weightSum+=static_cast<float>(degree); 1066 featOut+=featIn; 1067 featOut/=weightSum; 1068 } 1069 } 1070 1071 template<class T> 1072 struct ExpSmoothFactor{ ExpSmoothFactorvigra::detail_graph_smoothing::ExpSmoothFactor1073 ExpSmoothFactor(const T lambda,const T edgeThreshold,const T scale) 1074 : lambda_(lambda), 1075 edgeThreshold_(edgeThreshold), 1076 scale_(scale){ 1077 } operator ()vigra::detail_graph_smoothing::ExpSmoothFactor1078 T operator()(const T weight){ 1079 return weight> edgeThreshold_ ? 0 : std::exp(-1.0*lambda_*weight)*scale_; 1080 } 1081 T lambda_; 1082 T edgeThreshold_; 1083 T scale_; 1084 }; 1085 1086 } // namespace detail_graph_smoothing 1087 1088 1089 /// \brief smooth node features of a graph 1090 /// 1091 /// \param g : input graph 1092 /// \param nodeFeaturesIn : input node features which should be smoothed 1093 /// \param edgeIndicator : edge indicator to indicate over which edges one should smooth 1094 /// \param lambda : scale edge indicator by lambda before taking negative exponent 1095 /// \param edgeThreshold : edge threshold 1096 /// \param scale : how much smoothing should be applied 1097 /// \param[out] nodeFeaturesOut : smoothed node features 1098 template<class GRAPH, class NODE_FEATURES_IN,class EDGE_INDICATOR,class NODE_FEATURES_OUT> graphSmoothing(const GRAPH & g,const NODE_FEATURES_IN & nodeFeaturesIn,const EDGE_INDICATOR & edgeIndicator,const float lambda,const float edgeThreshold,const float scale,NODE_FEATURES_OUT & nodeFeaturesOut)1099 void graphSmoothing( 1100 const GRAPH & g, 1101 const NODE_FEATURES_IN & nodeFeaturesIn, 1102 const EDGE_INDICATOR & edgeIndicator, 1103 const float lambda, 1104 const float edgeThreshold, 1105 const float scale, 1106 NODE_FEATURES_OUT & nodeFeaturesOut 1107 ){ 1108 detail_graph_smoothing::ExpSmoothFactor<float> functor(lambda,edgeThreshold,scale); 1109 detail_graph_smoothing::graphSmoothingImpl(g,nodeFeaturesIn,edgeIndicator,functor,nodeFeaturesOut); 1110 } 1111 1112 /// \brief smooth node features of a graph 1113 /// 1114 /// \param g : input graph 1115 /// \param nodeFeaturesIn : input node features which should be smoothed 1116 /// \param edgeIndicator : edge indicator to indicate over which edges one should smooth 1117 /// \param lambda : scale edge indicator by lambda before taking negative exponent 1118 /// \param edgeThreshold : edge threshold 1119 /// \param scale : how much smoothing should be applied 1120 /// \param iterations : how often should this algorithm be called recursively 1121 /// \param[out] nodeFeaturesBuffer : preallocated(!) buffer to store node features temp. 1122 /// \param[out] nodeFeaturesOut : smoothed node features 1123 template<class GRAPH, class NODE_FEATURES_IN,class EDGE_INDICATOR,class NODE_FEATURES_OUT> recursiveGraphSmoothing(const GRAPH & g,const NODE_FEATURES_IN & nodeFeaturesIn,const EDGE_INDICATOR & edgeIndicator,const float lambda,const float edgeThreshold,const float scale,size_t iterations,NODE_FEATURES_OUT & nodeFeaturesBuffer,NODE_FEATURES_OUT & nodeFeaturesOut)1124 void recursiveGraphSmoothing( 1125 const GRAPH & g, 1126 const NODE_FEATURES_IN & nodeFeaturesIn, 1127 const EDGE_INDICATOR & edgeIndicator, 1128 const float lambda, 1129 const float edgeThreshold, 1130 const float scale, 1131 size_t iterations, 1132 NODE_FEATURES_OUT & nodeFeaturesBuffer, 1133 NODE_FEATURES_OUT & nodeFeaturesOut 1134 ){ 1135 1136 iterations = std::max(size_t(1),iterations); 1137 // initial run 1138 graphSmoothing(g,nodeFeaturesIn,edgeIndicator,lambda,edgeThreshold,scale,nodeFeaturesOut); 1139 iterations -=1; 1140 1141 bool outAsIn=true; 1142 for(size_t i=0;i<iterations;++i){ 1143 if(outAsIn){ 1144 graphSmoothing(g,nodeFeaturesOut,edgeIndicator,lambda,edgeThreshold,scale,nodeFeaturesBuffer); 1145 outAsIn=false; 1146 } 1147 else{ 1148 graphSmoothing(g,nodeFeaturesBuffer,edgeIndicator,lambda,edgeThreshold,scale,nodeFeaturesOut); 1149 outAsIn=true; 1150 } 1151 } 1152 if(!outAsIn){ 1153 copyNodeMap(g,nodeFeaturesBuffer,nodeFeaturesOut); 1154 } 1155 } 1156 1157 1158 template<class GRAPH, class NODE_MAP, class EDGE_MAP> nodeGtToEdgeGt(const GRAPH & g,const NODE_MAP & nodeGt,const Int64 ignoreLabel,EDGE_MAP & edgeGt)1159 void nodeGtToEdgeGt( 1160 const GRAPH & g, 1161 const NODE_MAP & nodeGt, 1162 const Int64 ignoreLabel, 1163 EDGE_MAP & edgeGt 1164 ){ 1165 typedef typename GRAPH::Node Node; 1166 typedef typename GRAPH::EdgeIt EdgeIt; 1167 typedef typename GRAPH::Edge Edge; 1168 1169 for(EdgeIt edgeIt(g); edgeIt!=lemon::INVALID; ++edgeIt){ 1170 const Edge edge(*edgeIt); 1171 const Node u = g.u(edge); 1172 const Node v = g.v(edge); 1173 1174 const Int64 lU = static_cast<Int64>(nodeGt[u]); 1175 const Int64 lV = static_cast<Int64>(nodeGt[v]); 1176 1177 if(ignoreLabel==-1 || (lU!=ignoreLabel || lV!=ignoreLabel)){ 1178 edgeGt[edge] = lU == lV ? 0 : 1; 1179 } 1180 else{ 1181 edgeGt[edge] = 2; 1182 } 1183 } 1184 } 1185 1186 1187 1188 1189 1190 1191 /// project ground truth from a base graph, to a region adjacency graph. 1192 /// 1193 /// 1194 /// 1195 /// 1196 template<class RAG, class BASE_GRAPH, class BASE_GRAPH_RAG_LABELS, 1197 class BASE_GRAPH_GT, class RAG_GT, class RAG_GT_QT> projectGroundTruth(const RAG & rag,const BASE_GRAPH & baseGraph,const BASE_GRAPH_RAG_LABELS & baseGraphRagLabels,const BASE_GRAPH_GT & baseGraphGt,RAG_GT & ragGt,RAG_GT_QT &)1198 void projectGroundTruth( 1199 const RAG & rag, 1200 const BASE_GRAPH & baseGraph, 1201 const BASE_GRAPH_RAG_LABELS & baseGraphRagLabels, 1202 const BASE_GRAPH_GT & baseGraphGt, 1203 RAG_GT & ragGt, 1204 RAG_GT_QT & 1205 ){ 1206 typedef typename BASE_GRAPH::Node BaseGraphNode; 1207 typedef typename BASE_GRAPH::NodeIt BaseGraphNodeIt; 1208 typedef typename RAG::NodeIt RagNodeIt; 1209 typedef typename RAG::Node RagNode; 1210 typedef typename BASE_GRAPH_RAG_LABELS::Value BaseRagLabelType; 1211 typedef typename BASE_GRAPH_GT::Value GtLabelType; 1212 typedef typename RAG_GT::Value RagGtLabelType; 1213 1214 // overlap map 1215 typedef std::map<GtLabelType,UInt32> MapType; 1216 typedef typename MapType::const_iterator MapIter; 1217 typedef typename RAG:: template NodeMap<MapType> Overlap; 1218 Overlap overlap(rag); 1219 1220 size_t i=0; 1221 //::cout<<"\n"; 1222 for(BaseGraphNodeIt baseNodeIter(baseGraph); baseNodeIter!=lemon::INVALID; ++baseNodeIter , ++i ){ 1223 1224 //if (i%2000 == 0){ 1225 // std::cout<<"\r"<<std::setw(10)<< float(i)/float(baseGraph.nodeNum())<<std::flush; 1226 //} 1227 1228 const BaseGraphNode baseNode = *baseNodeIter; 1229 1230 // get gt label 1231 const GtLabelType gtLabel = baseGraphGt[baseNode]; 1232 1233 // get the label which generated rag 1234 // (node mapping from bg-node to rag-node-id) 1235 const BaseRagLabelType bgRagLabel = baseGraphRagLabels[baseNode]; 1236 const RagNode ragNode = rag.nodeFromId(bgRagLabel); 1237 1238 // fill overlap 1239 overlap[ragNode][gtLabel]+=1; 1240 } 1241 //std::cout<<"\n"; 1242 // select label with max overlap 1243 for(RagNodeIt ragNodeIter(rag); ragNodeIter!=lemon::INVALID; ++ragNodeIter ){ 1244 const RagNode ragNode = *ragNodeIter; 1245 const MapType olMap = overlap[ragNode]; 1246 UInt32 olSize=0; 1247 RagGtLabelType bestLabel = 0; 1248 for(MapIter olIter = olMap.begin(); olIter!=olMap.end(); ++olIter ){ 1249 if(olIter->second > olSize){ 1250 olSize = olIter->second; 1251 bestLabel = olIter->first; 1252 } 1253 } 1254 ragGt[ragNode]=bestLabel; 1255 } 1256 } 1257 1258 1259 1260 /// \brief Find indices of points on the edges 1261 /// 1262 /// \param rag : Region adjacency graph of the labels array 1263 /// \param graph : Graph of labels array 1264 /// \param affiliatedEdges : The affiliated edges of the region adjacency graph 1265 /// \param labelsArray : The label image 1266 /// \param node : The node (of the region adjacency graph), whose edges shall be found 1267 template<class RAGGRAPH, class GRAPH, class RAGEDGES, unsigned int N, class T> ragFindEdges(const RAGGRAPH & rag,const GRAPH & graph,const RAGEDGES & affiliatedEdges,MultiArrayView<N,T> labelsArray,const typename RAGGRAPH::Node & node)1268 MultiArray<2, MultiArrayIndex> ragFindEdges( 1269 const RAGGRAPH & rag, 1270 const GRAPH & graph, 1271 const RAGEDGES & affiliatedEdges, 1272 MultiArrayView<N, T> labelsArray, 1273 const typename RAGGRAPH::Node & node 1274 ){ 1275 typedef typename GRAPH::Node Node; 1276 typedef typename GRAPH::Edge Edge; 1277 typedef typename RAGGRAPH::OutArcIt RagOutArcIt; 1278 typedef typename RAGGRAPH::Edge RagEdge; 1279 typedef typename GraphDescriptorToMultiArrayIndex<GRAPH>::IntrinsicNodeMapShape NodeCoordinate; 1280 1281 T nodeLabel = rag.id(node); 1282 1283 // Find edges and write them into a set. 1284 std::set< NodeCoordinate > edgeCoordinates; 1285 for (RagOutArcIt iter(rag, node); iter != lemon::INVALID; ++iter) 1286 { 1287 const RagEdge ragEdge(*iter); 1288 const std::vector<Edge> & affEdges = affiliatedEdges[ragEdge]; 1289 for (int i=0; i<affEdges.size(); ++i) 1290 { 1291 Node u = graph.u(affEdges[i]); 1292 Node v = graph.v(affEdges[i]); 1293 T uLabel = labelsArray[u]; 1294 T vLabel = labelsArray[v]; 1295 1296 NodeCoordinate coords; 1297 if (uLabel == nodeLabel) 1298 { 1299 coords = GraphDescriptorToMultiArrayIndex<GRAPH>::intrinsicNodeCoordinate(graph, u); 1300 } 1301 else if (vLabel == nodeLabel) 1302 { 1303 coords = GraphDescriptorToMultiArrayIndex<GRAPH>::intrinsicNodeCoordinate(graph, v); 1304 } 1305 else 1306 { 1307 vigra_precondition(false, "You should not come to this part of the code."); 1308 } 1309 edgeCoordinates.insert(coords); 1310 } 1311 } 1312 1313 // Fill the return array. 1314 MultiArray<2, MultiArrayIndex> edgePoints(Shape2(edgeCoordinates.size(), N)); 1315 edgePoints.init(0); 1316 int next = 0; 1317 typedef typename std::set< NodeCoordinate >::iterator setIter; 1318 for (setIter iter = edgeCoordinates.begin(); iter!=edgeCoordinates.end(); ++iter) 1319 { 1320 NodeCoordinate coords = *iter; 1321 for (int k=0; k<coords.size(); ++k) 1322 { 1323 edgePoints(next, k) = coords[k]; 1324 } 1325 next++; 1326 } 1327 return edgePoints; 1328 } 1329 /// \brief create edge weights from node weights 1330 /// 1331 /// \param g : input graph 1332 /// \param nodeWeights : node property map holding node weights 1333 /// \param[out] edgeWeights : resulting edge weights 1334 /// \param euclidean : if 'true', multiply the computed weights with the Euclidean 1335 /// distance between the edge's end nodes (default: 'false') 1336 /// \param func : binary function that computes the edge weight from the 1337 /// weights of the edge's end nodes (default: take the average) 1338 template<unsigned int N, class DirectedTag, 1339 class NODEMAP, class EDGEMAP, class FUNCTOR> 1340 void edgeWeightsFromNodeWeights(const GridGraph<N,DirectedTag> & g,const NODEMAP & nodeWeights,EDGEMAP & edgeWeights,bool euclidean,FUNCTOR const & func)1341 edgeWeightsFromNodeWeights( 1342 const GridGraph<N, DirectedTag> & g, 1343 const NODEMAP & nodeWeights, 1344 EDGEMAP & edgeWeights, 1345 bool euclidean, 1346 FUNCTOR const & func) 1347 { 1348 typedef GridGraph<N, DirectedTag> Graph; 1349 typedef typename Graph::Edge Edge; 1350 typedef typename Graph::EdgeIt EdgeIt; 1351 typedef typename MultiArrayShape<N>::type CoordType; 1352 1353 vigra_precondition(nodeWeights.shape() == g.shape(), 1354 "edgeWeightsFromNodeWeights(): shape mismatch between graph and nodeWeights."); 1355 1356 for (EdgeIt iter(g); iter!=lemon::INVALID; ++iter) 1357 { 1358 const Edge edge(*iter); 1359 const CoordType uCoord(g.u(edge)); 1360 const CoordType vCoord(g.v(edge)); 1361 if (euclidean) 1362 { 1363 edgeWeights[edge] = norm(uCoord-vCoord) * func(nodeWeights[uCoord], nodeWeights[vCoord]); 1364 } 1365 else 1366 { 1367 edgeWeights[edge] = func(nodeWeights[uCoord], nodeWeights[vCoord]); 1368 } 1369 } 1370 } 1371 1372 template<unsigned int N, class DirectedTag, 1373 class NODEMAP, class EDGEMAP> 1374 inline void edgeWeightsFromNodeWeights(const GridGraph<N,DirectedTag> & g,const NODEMAP & nodeWeights,EDGEMAP & edgeWeights,bool euclidean=false)1375 edgeWeightsFromNodeWeights( 1376 const GridGraph<N, DirectedTag> & g, 1377 const NODEMAP & nodeWeights, 1378 EDGEMAP & edgeWeights, 1379 bool euclidean=false) 1380 { 1381 using namespace vigra::functor; 1382 edgeWeightsFromNodeWeights(g, nodeWeights, edgeWeights, euclidean, Param(0.5)*(Arg1()+Arg2())); 1383 } 1384 1385 1386 /// \brief create edge weights from an interpolated image 1387 /// 1388 /// \param g : input graph 1389 /// \param interpolatedImage : interpolated image 1390 /// \param[out] edgeWeights : edge weights 1391 /// \param euclidean : if 'true', multiply the weights with the Euclidean 1392 /// distance between the edge's end nodes (default: 'false') 1393 /// 1394 /// For each edge, the function reads the weight from <tt>interpolatedImage[u+v]</tt>, 1395 /// where <tt>u</tt> and <tt>v</tt> are the coordinates of the edge's end points. 1396 template<unsigned int N, class DirectedTag, 1397 class T, class EDGEMAP> 1398 void edgeWeightsFromInterpolatedImage(const GridGraph<N,DirectedTag> & g,const MultiArrayView<N,T> & interpolatedImage,EDGEMAP & edgeWeights,bool euclidean=false)1399 edgeWeightsFromInterpolatedImage( 1400 const GridGraph<N, DirectedTag> & g, 1401 const MultiArrayView<N, T> & interpolatedImage, 1402 EDGEMAP & edgeWeights, 1403 bool euclidean = false) 1404 { 1405 typedef GridGraph<N, DirectedTag> Graph; 1406 typedef typename Graph::Edge Edge; 1407 typedef typename Graph::EdgeIt EdgeIt; 1408 typedef typename MultiArrayShape<N>::type CoordType; 1409 1410 vigra_precondition(interpolatedImage.shape() == 2*g.shape()-CoordType(1), 1411 "edgeWeightsFromInterpolatedImage(): interpolated shape must be shape*2-1"); 1412 1413 for (EdgeIt iter(g); iter!=lemon::INVALID; ++iter) 1414 { 1415 const Edge edge(*iter); 1416 const CoordType uCoord(g.u(edge)); 1417 const CoordType vCoord(g.v(edge)); 1418 if (euclidean) 1419 { 1420 edgeWeights[edge] = norm(uCoord-vCoord) * interpolatedImage[uCoord+vCoord]; 1421 } 1422 else 1423 { 1424 edgeWeights[edge] = interpolatedImage[uCoord+vCoord]; 1425 } 1426 } 1427 } 1428 1429 template<class GRAPH> 1430 struct ThreeCycle{ 1431 1432 typedef typename GRAPH::Node Node; 1433 ThreeCyclevigra::ThreeCycle1434 ThreeCycle(const Node & a, const Node & b, const Node c){ 1435 nodes_[0] = a; 1436 nodes_[1] = b; 1437 nodes_[2] = c; 1438 std::sort(nodes_, nodes_+3); 1439 } operator <vigra::ThreeCycle1440 bool operator < (const ThreeCycle & other)const{ 1441 if(nodes_[0] < other.nodes_[0]){ 1442 return true; 1443 } 1444 else if(nodes_[0] == other.nodes_[0]){ 1445 if(nodes_[1] < other.nodes_[1]){ 1446 return true; 1447 } 1448 else if(nodes_[1] == other.nodes_[1]){ 1449 if(nodes_[2] < other.nodes_[2]){ 1450 return true; 1451 } 1452 else{ 1453 return false; 1454 } 1455 } 1456 else{ 1457 return false; 1458 } 1459 } 1460 else{ 1461 return false; 1462 } 1463 } 1464 1465 Node nodes_[3]; 1466 1467 1468 }; 1469 1470 1471 template<class GRAPH> find3Cycles(const GRAPH & g,MultiArray<1,TinyVector<Int32,3>> & cyclesArray)1472 void find3Cycles( 1473 const GRAPH & g, 1474 MultiArray<1, TinyVector<Int32, 3> > & cyclesArray 1475 ){ 1476 typedef typename GRAPH::Node Node; 1477 typedef typename GRAPH::Edge Edge; 1478 typedef typename GRAPH::EdgeIt EdgeIt; 1479 typedef typename GRAPH::OutArcIt OutArcIt; 1480 1481 typedef ThreeCycle<GRAPH> Cycle; 1482 1483 std::set< Cycle > cycles; 1484 typedef typename std::set<Cycle>::const_iterator SetIter; 1485 for (EdgeIt iter(g); iter!=lemon::INVALID; ++iter){ 1486 const Edge edge(*iter); 1487 const Node u = g.u(edge); 1488 const Node v = g.v(edge); 1489 1490 // find a node n which is connected to u and v 1491 for(OutArcIt outArcIt(g,u); outArcIt!=lemon::INVALID;++outArcIt){ 1492 const Node w = g.target(*outArcIt); 1493 if(w != v){ 1494 const Edge e = g.findEdge(w,v); 1495 if(e != lemon::INVALID){ 1496 // found cycle 1497 cycles.insert(Cycle(u, v, w)); 1498 } 1499 } 1500 } 1501 } 1502 cyclesArray.reshape(TinyVector<UInt32,1>(cycles.size())); 1503 UInt32 i=0; 1504 for(SetIter iter=cycles.begin(); iter!=cycles.end(); ++iter){ 1505 1506 const Cycle & c = *iter; 1507 for(size_t j=0;j<3; ++j){ 1508 cyclesArray(i)[j] = g.id(c.nodes_[j]); 1509 } 1510 ++i; 1511 } 1512 } 1513 1514 template<class GRAPH> find3CyclesEdges(const GRAPH & g,MultiArray<1,TinyVector<Int32,3>> & cyclesArray)1515 void find3CyclesEdges( 1516 const GRAPH & g, 1517 MultiArray<1, TinyVector<Int32, 3> > & cyclesArray 1518 ){ 1519 typedef typename GRAPH::Node Node; 1520 typedef typename GRAPH::Edge Edge; 1521 typedef typename GRAPH::EdgeIt EdgeIt; 1522 typedef typename GRAPH::OutArcIt OutArcIt; 1523 1524 typedef ThreeCycle<GRAPH> Cycle; 1525 1526 std::set< Cycle > cycles; 1527 typedef typename std::set<Cycle>::const_iterator SetIter; 1528 for (EdgeIt iter(g); iter!=lemon::INVALID; ++iter){ 1529 const Edge edge(*iter); 1530 const Node u = g.u(edge); 1531 const Node v = g.v(edge); 1532 1533 // find a node n which is connected to u and v 1534 for(OutArcIt outArcIt(g,u); outArcIt!=lemon::INVALID;++outArcIt){ 1535 const Node w = g.target(*outArcIt); 1536 if(w != v){ 1537 const Edge e = g.findEdge(w,v); 1538 if(e != lemon::INVALID){ 1539 // found cycle 1540 cycles.insert(Cycle(u, v, w)); 1541 } 1542 } 1543 } 1544 } 1545 cyclesArray.reshape(TinyVector<UInt32,1>(cycles.size())); 1546 UInt32 i=0; 1547 for(SetIter iter=cycles.begin(); iter!=cycles.end(); ++iter){ 1548 1549 const Cycle & c = *iter; 1550 const Node u = c.nodes_[0]; 1551 const Node v = c.nodes_[1]; 1552 const Node w = c.nodes_[2]; 1553 1554 cyclesArray(i)[0] = g.id(g.findEdge(u, v)); 1555 cyclesArray(i)[1] = g.id(g.findEdge(u, w)); 1556 cyclesArray(i)[2] = g.id(g.findEdge(v, w)); 1557 ++i; 1558 } 1559 } 1560 1561 //@} 1562 1563 } // namespace vigra 1564 1565 #endif // VIGRA_GRAPH_ALGORITHMS_HXX 1566