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

graph_algorithms.hxx VIGRA

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 
77  GraphItemCompare(const GRAPH_MAP & map,const COMPERATOR & comperator)
78  : map_(map),
79  comperator_(comperator){
80 
81  }
82 
83  template<class KEY>
84  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>
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>
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>
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>
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>
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  >
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>
217  size_t affiliatedEdgesSerializationSize(
218  const GridGraph<DIM,DTAG> & gridGraph,
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>
236  void serializeAffiliatedEdges(
237  const GridGraph<DIM,DTAG> & gridGraph,
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>
263  void deserializeAffiliatedEdges(
264  const GridGraph<DIM,DTAG> & gridGraph,
265  const AdjacencyListGraph & rag,
266  AFF_EDGES & affEdges,
267  IN_ITER begin,
268  IN_ITER end
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
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>
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>
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>
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
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
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
419  const Graph & graph()const{
420  return graph_;
421  }
422  /// \brief get the source node
423  const Node & source()const{
424  return source_;
425  }
426  /// \brief get the target node
427  const Node & target()const{
428  return target_;
429  }
430 
431  /// \brief check if explicit target is given
432  bool hasTarget()const{
433  return target_!=lemon::INVALID;
434  }
435 
436  /// \brief get an array with all visited nodes, sorted by distance from source
437  const DiscoveryOrder & discoveryOrder() const{
438  return discoveryOrder_;
439  }
440 
441  /// \brief get the predecessors node map (after a call of run)
442  const PredecessorsMap & predecessors()const{
443  return predMap_;
444  }
445 
446  /// \brief get the distances node map (after a call of run)
447  const DistanceMap & distances()const{
448  return distMap_;
449  }
450 
451  /// \brief get the distance to a rarget node (after a call of run)
452  WeightType distance(const Node & target)const{
453  return distMap_[target];
454  }
455 
456 
457  private:
458 
459  template<class WEIGHTS>
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>
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 
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 
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>
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 
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>
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>
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  >
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  int label = 0 ;
728  Node pred=predMap[node];
729  while(seeds[pred]==0){
730  pred=predMap[pred];
731  }
732  seeds[node]=seeds[pred];
733  }
734  }
735  }
736 
737  namespace detail_watersheds_segmentation{
738 
739  struct RawPriorityFunctor{
740  template<class L, class T>
741  T operator()(const L label,const T priority)const{
742  return priority;
743  }
744  };
745 
746 
747 
748  template<class PRIORITY_TYPE,class LABEL_TYPE>
749  struct CarvingFunctor{
750  CarvingFunctor(const LABEL_TYPE backgroundLabel,
751  const PRIORITY_TYPE & factor,
752  const PRIORITY_TYPE & noPriorBelow
753  )
754  : backgroundLabel_(backgroundLabel),
755  factor_(factor),
756  noPriorBelow_(noPriorBelow){
757  }
758  PRIORITY_TYPE operator()(const LABEL_TYPE label,const PRIORITY_TYPE priority)const{
759  if(priority>=noPriorBelow_)
760  return (label==backgroundLabel_ ? priority*factor_ : priority);
761  else{
762  return priority;
763  }
764  }
765  LABEL_TYPE backgroundLabel_;
766  PRIORITY_TYPE factor_;
767  PRIORITY_TYPE noPriorBelow_;
768  };
769 
770 
771  template<
772  class GRAPH,
773  class EDGE_WEIGHTS,
774  class SEEDS,
775  class PRIORITY_MANIP_FUNCTOR,
776  class LABELS
777  >
778  void edgeWeightedWatershedsSegmentationImpl(
779  const GRAPH & g,
780  const EDGE_WEIGHTS & edgeWeights,
781  const SEEDS & seeds,
782  PRIORITY_MANIP_FUNCTOR & priorManipFunctor,
783  LABELS & labels
784  ){
785  typedef GRAPH Graph;
786  typedef typename Graph::Edge Edge;
787  typedef typename Graph::Node Node;
788  typedef typename Graph::NodeIt NodeIt;
789  typedef typename Graph::OutArcIt OutArcIt;
790 
791  typedef typename EDGE_WEIGHTS::Value WeightType;
792  typedef typename LABELS::Value LabelType;
793  //typedef typename Graph:: template EdgeMap<bool> EdgeBoolMap;
794  typedef PriorityQueue<Edge,WeightType,true> PQ;
795 
796  PQ pq;
797  //EdgeBoolMap inPQ(g);
798  copyNodeMap(g,seeds,labels);
799  //fillEdgeMap(g,inPQ,false);
800 
801 
802  // put edges from nodes with seed on pq
803  for(NodeIt n(g);n!=lemon::INVALID;++n){
804  const Node node(*n);
805  if(labels[node]!=static_cast<LabelType>(0)){
806  for(OutArcIt a(g,node);a!=lemon::INVALID;++a){
807  const Edge edge(*a);
808  const Node neigbour=g.target(*a);
809  //std::cout<<"n- node "<<g.id(neigbour)<<"\n";
810  if(labels[neigbour]==static_cast<LabelType>(0)){
811  const WeightType priority = priorManipFunctor(labels[node],edgeWeights[edge]);
812  pq.push(edge,priority);
813  //inPQ[edge]=true;
814  }
815  }
816  }
817  }
818 
819 
820  while(!pq.empty()){
821 
822  const Edge edge = pq.top();
823  pq.pop();
824 
825  const Node u = g.u(edge);
826  const Node v = g.v(edge);
827  const LabelType lU = labels[u];
828  const LabelType lV = labels[v];
829 
830 
831  if(lU==0 && lV==0){
832  throw std::runtime_error("both have no labels");
833  }
834  else if(lU!=0 && lV!=0){
835  // nothing to do
836  }
837  else{
838 
839  const Node unlabeledNode = lU==0 ? u : v;
840  const LabelType label = lU==0 ? lV : lU;
841 
842  // assign label to unlabeled node
843  labels[unlabeledNode] = label;
844 
845  // iterate over the nodes edges
846  for(OutArcIt a(g,unlabeledNode);a!=lemon::INVALID;++a){
847  const Edge otherEdge(*a);
848  const Node targetNode=g.target(*a);
849  if(labels[targetNode] == 0){
850  //if(inPQ[otherEdge] == false && labels[targetNode] == 0){
851  const WeightType priority = priorManipFunctor(label,edgeWeights[otherEdge]);
852  pq.push(otherEdge,priority);
853  // inPQ[otherEdge]=true;
854  }
855  }
856  }
857 
858  }
859 
860  }
861 
862  } // end namespace detail_watersheds_segmentation
863 
864 
865  /// \brief edge weighted watersheds Segmentataion
866  ///
867  /// \param g: input graph
868  /// \param edgeWeights : edge weights / edge indicator
869  /// \param seeds : seed must be non empty!
870  /// \param[out] labels : resulting nodeLabeling (not necessarily dense)
871  template<class GRAPH,class EDGE_WEIGHTS,class SEEDS,class LABELS>
872  void edgeWeightedWatershedsSegmentation(
873  const GRAPH & g,
874  const EDGE_WEIGHTS & edgeWeights,
875  const SEEDS & seeds,
876  LABELS & labels
877  ){
878  detail_watersheds_segmentation::RawPriorityFunctor fPriority;
879  detail_watersheds_segmentation::edgeWeightedWatershedsSegmentationImpl(g,edgeWeights,seeds,fPriority,labels);
880  }
881 
882 
883  /// \brief edge weighted watersheds Segmentataion
884  ///
885  /// \param g: input graph
886  /// \param edgeWeights : edge weights / edge indicator
887  /// \param seeds : seed must be non empty!
888  /// \param backgroundLabel : which label is background
889  /// \param backgroundBias : bias for background
890  /// \param noPriorBelow : don't bias the background if edge indicator is below this value
891  /// \param[out] labels : resulting nodeLabeling (not necessarily dense)
892  template<class GRAPH,class EDGE_WEIGHTS,class SEEDS,class LABELS>
893  void carvingSegmentation(
894  const GRAPH & g,
895  const EDGE_WEIGHTS & edgeWeights,
896  const SEEDS & seeds,
897  const typename LABELS::Value backgroundLabel,
898  const typename EDGE_WEIGHTS::Value backgroundBias,
899  const typename EDGE_WEIGHTS::Value noPriorBelow,
900  LABELS & labels
901  ){
902  typedef typename EDGE_WEIGHTS::Value WeightType;
903  typedef typename LABELS::Value LabelType;
904  detail_watersheds_segmentation::CarvingFunctor<WeightType,LabelType> fPriority(backgroundLabel,backgroundBias, noPriorBelow);
905  detail_watersheds_segmentation::edgeWeightedWatershedsSegmentationImpl(g,edgeWeights,seeds,fPriority,labels);
906  }
907 
908  /// \brief edge weighted watersheds Segmentataion
909  ///
910  /// \param graph: input graph
911  /// \param edgeWeights : edge weights / edge indicator
912  /// \param nodeSizes : size of each node
913  /// \param k : free parameter of felzenszwalb algorithm
914  /// \param[out] nodeLabeling : nodeLabeling (not necessarily dense)
915  /// \param nodeNumStopCond : optional stopping condition
916  template< class GRAPH , class EDGE_WEIGHTS, class NODE_SIZE,class NODE_LABEL_MAP>
917  void felzenszwalbSegmentation(
918  const GRAPH & graph,
919  const EDGE_WEIGHTS & edgeWeights,
920  const NODE_SIZE & nodeSizes,
921  float k,
922  NODE_LABEL_MAP & nodeLabeling,
923  const int nodeNumStopCond = -1
924  ){
925  typedef GRAPH Graph;
926  typedef typename Graph::Edge Edge;
927  typedef typename Graph::Node Node;
928 
929  typedef typename EDGE_WEIGHTS::Value WeightType;
930  typedef typename EDGE_WEIGHTS::Value NodeSizeType;
931  typedef typename Graph:: template NodeMap<WeightType> NodeIntDiffMap;
932  typedef typename Graph:: template NodeMap<NodeSizeType> NodeSizeAccMap;
933 
934  // initalize node size map and internal diff map
935  NodeIntDiffMap internalDiff(graph);
936  NodeSizeAccMap nodeSizeAcc(graph);
937  copyNodeMap(graph,nodeSizes,nodeSizeAcc);
938  fillNodeMap(graph,internalDiff,static_cast<WeightType>(0.0));
939 
940 
941 
942  // initlaize internal node diff map
943 
944  // sort the edges by their weights
945  std::vector<Edge> sortedEdges;
946  std::less<WeightType> comperator;
947  edgeSort(graph,edgeWeights,comperator,sortedEdges);
948 
949  // make the ufd
950  UnionFindArray<UInt64> ufdArray(graph.maxNodeId()+1);
951 
952 
953  size_t nodeNum = graph.nodeNum();
954 
955 
956  while(true){
957  // iterate over edges is the sorted order
958  for(size_t i=0;i<sortedEdges.size();++i){
959  const Edge e = sortedEdges[i];
960  const size_t rui = ufdArray.findIndex(graph.id(graph.u(e)));
961  const size_t rvi = ufdArray.findIndex(graph.id(graph.v(e)));
962  const Node ru = graph.nodeFromId(rui);
963  const Node rv = graph.nodeFromId(rvi);
964  if(rui!=rvi){
965 
966  //check if to merge or not ?
967  const WeightType w = edgeWeights[e];
968  const NodeSizeType sizeRu = nodeSizeAcc[ru];
969  const NodeSizeType sizeRv = nodeSizeAcc[rv];
970  const WeightType tauRu = static_cast<WeightType>(k)/static_cast<WeightType>(sizeRu);
971  const WeightType tauRv = static_cast<WeightType>(k)/static_cast<WeightType>(sizeRv);
972  const WeightType minIntDiff = std::min(internalDiff[ru]+tauRu,internalDiff[rv]+tauRv);
973  if(w<=minIntDiff){
974  // do merge
975  ufdArray.makeUnion(rui,rvi);
976  --nodeNum;
977  // update size and internal difference
978  const size_t newRepId = ufdArray.findIndex(rui);
979  const Node newRepNode = graph.nodeFromId(newRepId);
980  internalDiff[newRepNode]=w;
981  nodeSizeAcc[newRepNode] = sizeRu+sizeRv;
982  }
983  }
984  if(nodeNum==nodeNumStopCond){
985  break;
986  }
987  }
988  if(nodeNumStopCond==-1){
989  break;
990  }
991  else{
992  if(nodeNum>nodeNumStopCond){
993  k *= 1.2f;
994  }
995  else{
996  break;
997  }
998  }
999  }
1000  ufdArray.makeContiguous();
1001  for(typename GRAPH::NodeIt n(graph);n!=lemon::INVALID;++n){
1002  const Node node(*n);
1003  nodeLabeling[node]=ufdArray.findLabel(graph.id(node));
1004  }
1005  }
1006 
1007 
1008 
1009 
1010  namespace detail_graph_smoothing{
1011 
1012  template<
1013  class GRAPH,
1014  class NODE_FEATURES_IN,
1015  class EDGE_WEIGHTS,
1016  class WEIGHTS_TO_SMOOTH_FACTOR,
1017  class NODE_FEATURES_OUT
1018  >
1019  void graphSmoothingImpl(
1020  const GRAPH & g,
1021  const NODE_FEATURES_IN & nodeFeaturesIn,
1022  const EDGE_WEIGHTS & edgeWeights,
1023  WEIGHTS_TO_SMOOTH_FACTOR & weightsToSmoothFactor,
1024  NODE_FEATURES_OUT & nodeFeaturesOut
1025 
1026  ){
1027  typedef GRAPH Graph;
1028  typedef typename Graph::Edge Edge;
1029  typedef typename Graph::Node Node;
1030  typedef typename Graph::NodeIt NodeIt;
1031  typedef typename Graph::OutArcIt OutArcIt;
1032 
1033  typedef typename NODE_FEATURES_IN::Value NodeFeatureInValue;
1034  typedef typename NODE_FEATURES_OUT::Reference NodeFeatureOutRef;
1035  typedef typename EDGE_WEIGHTS::ConstReference SmoothFactorType;
1036 
1037 
1038  //fillNodeMap(g, nodeFeaturesOut, typename NODE_FEATURES_OUT::value_type(0.0));
1039 
1040  for(NodeIt n(g);n!=lemon::INVALID;++n){
1041 
1042  const Node node(*n);
1043 
1044  NodeFeatureInValue featIn = nodeFeaturesIn[node];
1045  NodeFeatureOutRef featOut = nodeFeaturesOut[node];
1046 
1047  featOut=0;
1048  float weightSum = 0.0;
1049  size_t degree = 0;
1050  for(OutArcIt a(g,node);a!=lemon::INVALID;++a){
1051  const Edge edge(*a);
1052  const Node neigbour(g.target(*a));
1053  SmoothFactorType smoothFactor= weightsToSmoothFactor(edgeWeights[edge]);
1054 
1055  NodeFeatureInValue neighbourFeat = nodeFeaturesIn[neigbour];
1056  neighbourFeat*=smoothFactor;
1057  if(degree==0)
1058  featOut = neighbourFeat;
1059  else
1060  featOut += neighbourFeat;
1061  weightSum+=smoothFactor;
1062  ++degree;
1063  }
1064  // fixme..set me to right type
1065  featIn*=static_cast<float>(degree);
1066  weightSum+=static_cast<float>(degree);
1067  featOut+=featIn;
1068  featOut/=weightSum;
1069  }
1070  }
1071 
1072  template<class T>
1073  struct ExpSmoothFactor{
1074  ExpSmoothFactor(const T lambda,const T edgeThreshold,const T scale)
1075  : lambda_(lambda),
1076  edgeThreshold_(edgeThreshold),
1077  scale_(scale){
1078  }
1079