[ 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  >
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>
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;
308  typedef typename Graph:: template NodeMap<Node> PredecessorsMap;
309  typedef typename Graph:: template NodeMap<WeightType> DistanceMap;
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
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>
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
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;
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>
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>
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>
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  T operator()(const T weight){
1080  return weight> edgeThreshold_ ? 0 : std::exp(-1.0*lambda_*weight)*scale_;
1081  }
1082  T lambda_;
1083  T edgeThreshold_;
1084  T scale_;
1085  };
1086 
1087  } // namespace detail_graph_smoothing
1088 
1089 
1090  /// \brief smooth node features of a graph
1091  ///
1092  /// \param g : input graph
1093  /// \param nodeFeaturesIn : input node features which should be smoothed
1094  /// \param edgeIndicator : edge indicator to indicate over which edges one should smooth
1095  /// \param lambda : scale edge indicator by lambda before taking negative exponent
1096  /// \param edgeThreshold : edge threshold
1097  /// \param scale : how much smoothing should be applied
1098  /// \param[out] nodeFeaturesOut : smoothed node features
1099  template<class GRAPH, class NODE_FEATURES_IN,class EDGE_INDICATOR,class NODE_FEATURES_OUT>
1101  const GRAPH & g,
1102  const NODE_FEATURES_IN & nodeFeaturesIn,
1103  const EDGE_INDICATOR & edgeIndicator,
1104  const float lambda,
1105  const float edgeThreshold,
1106  const float scale,
1107  NODE_FEATURES_OUT & nodeFeaturesOut
1108  ){
1109  detail_graph_smoothing::ExpSmoothFactor<float> functor(lambda,edgeThreshold,scale);
1110  detail_graph_smoothing::graphSmoothingImpl(g,nodeFeaturesIn,edgeIndicator,functor,nodeFeaturesOut);
1111  }
1112 
1113  /// \brief smooth node features of a graph
1114  ///
1115  /// \param g : input graph
1116  /// \param nodeFeaturesIn : input node features which should be smoothed
1117  /// \param edgeIndicator : edge indicator to indicate over which edges one should smooth
1118  /// \param lambda : scale edge indicator by lambda before taking negative exponent
1119  /// \param edgeThreshold : edge threshold
1120  /// \param scale : how much smoothing should be applied
1121  /// \param iterations : how often should this algorithm be called recursively
1122  /// \param[out] nodeFeaturesBuffer : preallocated(!) buffer to store node features temp.
1123  /// \param[out] nodeFeaturesOut : smoothed node features
1124  template<class GRAPH, class NODE_FEATURES_IN,class EDGE_INDICATOR,class NODE_FEATURES_OUT>
1126  const GRAPH & g,
1127  const NODE_FEATURES_IN & nodeFeaturesIn,
1128  const EDGE_INDICATOR & edgeIndicator,
1129  const float lambda,
1130  const float edgeThreshold,
1131  const float scale,
1132  size_t iterations,
1133  NODE_FEATURES_OUT & nodeFeaturesBuffer,
1134  NODE_FEATURES_OUT & nodeFeaturesOut
1135  ){
1136 
1137  iterations = std::max(size_t(1),iterations);
1138  // initial run
1139  graphSmoothing(g,nodeFeaturesIn,edgeIndicator,lambda,edgeThreshold,scale,nodeFeaturesOut);
1140  iterations -=1;
1141 
1142  bool outAsIn=true;
1143  for(size_t i=0;i<iterations;++i){
1144  if(outAsIn){
1145  graphSmoothing(g,nodeFeaturesOut,edgeIndicator,lambda,edgeThreshold,scale,nodeFeaturesBuffer);
1146  outAsIn=false;
1147  }
1148  else{
1149  graphSmoothing(g,nodeFeaturesBuffer,edgeIndicator,lambda,edgeThreshold,scale,nodeFeaturesOut);
1150  outAsIn=true;
1151  }
1152  }
1153  if(!outAsIn){
1154  copyNodeMap(g,nodeFeaturesBuffer,nodeFeaturesOut);
1155  }
1156  }
1157 
1158 
1159  template<class GRAPH, class NODE_MAP, class EDGE_MAP>
1160  void nodeGtToEdgeGt(
1161  const GRAPH & g,
1162  const NODE_MAP & nodeGt,
1163  const Int64 ignoreLabel,
1164  EDGE_MAP & edgeGt
1165  ){
1166  typedef typename GRAPH::Node Node;
1167  typedef typename GRAPH::EdgeIt EdgeIt;
1168  typedef typename GRAPH::Edge Edge;
1169 
1170  for(EdgeIt edgeIt(g); edgeIt!=lemon::INVALID; ++edgeIt){
1171  const Edge edge(*edgeIt);
1172  const Node u = g.u(edge);
1173  const Node v = g.v(edge);
1174 
1175  const Int64 lU = static_cast<Int64>(nodeGt[u]);
1176  const Int64 lV = static_cast<Int64>(nodeGt[v]);
1177 
1178  if(ignoreLabel==-1 || (lU!=ignoreLabel || lV!=ignoreLabel)){
1179  edgeGt[edge] = lU == lV ? 0 : 1;
1180  }
1181  else{
1182  edgeGt[edge] = 2;
1183  }
1184  }
1185  }
1186 
1187 
1188 
1189 
1190 
1191 
1192  /// project ground truth from a base graph, to a region adjacency graph.
1193  ///
1194  ///
1195  ///
1196  ///
1197  template<class RAG, class BASE_GRAPH, class BASE_GRAPH_RAG_LABELS,
1198  class BASE_GRAPH_GT, class RAG_GT, class RAG_GT_QT>
1200  const RAG & rag,
1201  const BASE_GRAPH & baseGraph,
1202  const BASE_GRAPH_RAG_LABELS & baseGraphRagLabels,
1203  const BASE_GRAPH_GT & baseGraphGt,
1204  RAG_GT & ragGt,
1205  RAG_GT_QT & ragGtQt
1206  ){
1207  typedef typename BASE_GRAPH::Node BaseGraphNode;
1208  typedef typename BASE_GRAPH::NodeIt BaseGraphNodeIt;
1209  typedef typename RAG::NodeIt RagNodeIt;
1210  typedef typename RAG::Node RagNode;
1211  typedef typename BASE_GRAPH_RAG_LABELS::Value BaseRagLabelType;
1212  typedef typename BASE_GRAPH_GT::Value GtLabelType;
1213  typedef typename RAG_GT::Value RagGtLabelType;
1214 
1215  // overlap map
1216  typedef std::map<GtLabelType,UInt32> MapType;
1217  typedef typename MapType::const_iterator MapIter;
1218  typedef typename RAG:: template NodeMap<MapType> Overlap;
1219  Overlap overlap(rag);
1220 
1221  size_t i=0;
1222  //::cout<<"\n";
1223  for(BaseGraphNodeIt baseNodeIter(baseGraph); baseNodeIter!=lemon::INVALID; ++baseNodeIter , ++i ){
1224 
1225  //if (i%2000 == 0){
1226  // std::cout<<"\r"<<std::setw(10)<< float(i)/float(baseGraph.nodeNum())<<std::flush;
1227  //}
1228 
1229  const BaseGraphNode baseNode = *baseNodeIter;
1230 
1231  // get gt label
1232  const GtLabelType gtLabel = baseGraphGt[baseNode];
1233 
1234  // get the label which generated rag
1235  // (node mapping from bg-node to rag-node-id)
1236  const BaseRagLabelType bgRagLabel = baseGraphRagLabels[baseNode];
1237  const RagNode ragNode = rag.nodeFromId(bgRagLabel);
1238 
1239  // fill overlap
1240  overlap[ragNode][gtLabel]+=1;
1241  }
1242  //std::cout<<"\n";
1243  // select label with max overlap
1244  for(RagNodeIt ragNodeIter(rag); ragNodeIter!=lemon::INVALID; ++ragNodeIter ){
1245  const RagNode ragNode = *ragNodeIter;
1246  const MapType olMap = overlap[ragNode];
1247  UInt32 olSize=0;
1248  RagGtLabelType bestLabel = 0;
1249  for(MapIter olIter = olMap.begin(); olIter!=olMap.end(); ++olIter ){
1250  if(olIter->second > olSize){
1251  olSize = olIter->second;
1252  bestLabel = olIter->first;
1253  }
1254  }
1255  ragGt[ragNode]=bestLabel;
1256  }
1257  }
1258 
1259 
1260 
1261  /// \brief Find indices of points on the edges
1262  ///
1263  /// \param rag : Region adjacency graph of the labels array
1264  /// \param graph : Graph of labels array
1265  /// \param affiliatedEdges : The affiliated edges of the region adjacency graph
1266  /// \param labelsArray : The label image
1267  /// \param node : The node (of the region adjacency graph), whose edges shall be found
1268  template<class RAGGRAPH, class GRAPH, class RAGEDGES, unsigned int N, class T>
1270  const RAGGRAPH & rag,
1271  const GRAPH & graph,
1272  const RAGEDGES & affiliatedEdges,
1273  MultiArrayView<N, T> labelsArray,
1274  const typename RAGGRAPH::Node & node
1275  ){
1276  typedef typename GRAPH::Node Node;
1277  typedef typename GRAPH::Edge Edge;
1278  typedef typename RAGGRAPH::OutArcIt RagOutArcIt;
1279  typedef typename RAGGRAPH::Edge RagEdge;
1280  typedef typename GraphDescriptorToMultiArrayIndex<GRAPH>::IntrinsicNodeMapShape NodeCoordinate;
1281 
1282  T nodeLabel = rag.id(node);
1283 
1284  // Find edges and write them into a set.
1285  std::set< NodeCoordinate > edgeCoordinates;
1286  for (RagOutArcIt iter(rag, node); iter != lemon::INVALID; ++iter)
1287  {
1288  const RagEdge ragEdge(*iter);
1289  const std::vector<Edge> & affEdges = affiliatedEdges[ragEdge];
1290  for (int i=0; i<affEdges.size(); ++i)
1291  {
1292  Node u = graph.u(affEdges[i]);
1293  Node v = graph.v(affEdges[i]);
1294  T uLabel = labelsArray[u];
1295  T vLabel = labelsArray[v];
1296 
1297  NodeCoordinate coords;
1298  if (uLabel == nodeLabel)
1299  {
1300  coords = GraphDescriptorToMultiArrayIndex<GRAPH>::intrinsicNodeCoordinate(graph, u);
1301  }
1302  else if (vLabel == nodeLabel)
1303  {
1304  coords = GraphDescriptorToMultiArrayIndex<GRAPH>::intrinsicNodeCoordinate(graph, v);
1305  }
1306  else
1307  {
1308  vigra_precondition(false, "You should not come to this part of the code.");
1309  }
1310  edgeCoordinates.insert(coords);
1311  }
1312  }
1313 
1314  // Fill the return array.
1315  MultiArray<2, MultiArrayIndex> edgePoints(Shape2(edgeCoordinates.size(), N));
1316  edgePoints.init(0);
1317  int next = 0;
1318  typedef typename std::set< NodeCoordinate >::iterator setIter;
1319  for (setIter iter = edgeCoordinates.begin(); iter!=edgeCoordinates.end(); ++iter)
1320  {
1321  NodeCoordinate coords = *iter;
1322  for (int k=0; k<coords.size(); ++k)
1323  {
1324  edgePoints(next, k) = coords[k];
1325  }
1326  next++;
1327  }
1328  return edgePoints;
1329  }
1330  /// \brief create edge weights from node weights
1331  ///
1332  /// \param g : input graph
1333  /// \param nodeWeights : node property map holding node weights
1334  /// \param[out] edgeWeights : resulting edge weights
1335  /// \param euclidean : if 'true', multiply the computed weights with the Euclidean
1336  /// distance between the edge's end nodes (default: 'false')
1337  /// \param func : binary function that computes the edge weight from the
1338  /// weights of the edge's end nodes (default: take the average)
1339  template<unsigned int N, class DirectedTag,
1340  class NODEMAP, class EDGEMAP, class FUNCTOR>
1341  void
1343  const GridGraph<N, DirectedTag> & g,
1344  const NODEMAP & nodeWeights,
1345  EDGEMAP & edgeWeights,
1346  bool euclidean,
1347  FUNCTOR const & func)
1348  {
1349  typedef GridGraph<N, DirectedTag> Graph;
1350  typedef typename Graph::Edge Edge;
1351  typedef typename Graph::EdgeIt EdgeIt;
1352  typedef typename MultiArrayShape<N>::type CoordType;
1353 
1354  vigra_precondition(nodeWeights.shape() == g.shape(),
1355  "edgeWeightsFromNodeWeights(): shape mismatch between graph and nodeWeights.");
1356 
1357  for (EdgeIt iter(g); iter!=lemon::INVALID; ++iter)
1358  {
1359  const Edge edge(*iter);
1360  const CoordType uCoord(g.u(edge));
1361  const CoordType vCoord(g.v(edge));
1362  if (euclidean)
1363  {
1364  edgeWeights[edge] = norm(uCoord-vCoord) * func(nodeWeights[uCoord], nodeWeights[vCoord]);
1365  }
1366  else
1367  {
1368  edgeWeights[edge] = func(nodeWeights[uCoord], nodeWeights[vCoord]);
1369  }
1370  }
1371  }
1372 
1373  template<unsigned int N, class DirectedTag,
1374  class NODEMAP, class EDGEMAP>
1375  inline void
1377  const GridGraph<N, DirectedTag> & g,
1378  const NODEMAP & nodeWeights,
1379  EDGEMAP & edgeWeights,
1380  bool euclidean=false)
1381  {
1382  using namespace vigra::functor;
1383  edgeWeightsFromNodeWeights(g, nodeWeights, edgeWeights, euclidean, Param(0.5)*(Arg1()+Arg2()));
1384  }
1385 
1386 
1387  /// \brief create edge weights from an interpolated image
1388  ///
1389  /// \param g : input graph
1390  /// \param interpolatedImage : interpolated image
1391  /// \param[out] edgeWeights : edge weights
1392  /// \param euclidean : if 'true', multiply the weights with the Euclidean
1393  /// distance between the edge's end nodes (default: 'false')
1394  ///
1395  /// For each edge, the function reads the weight from <tt>interpolatedImage[u+v]</tt>,
1396  /// where <tt>u</tt> and <tt>v</tt> are the coordinates of the edge's end points.
1397  template<unsigned int N, class DirectedTag,
1398  class T, class EDGEMAP>
1399  void
1401  const GridGraph<N, DirectedTag> & g,
1402  const MultiArrayView<N, T> & interpolatedImage,
1403  EDGEMAP & edgeWeights,
1404  bool euclidean = false)
1405  {
1406  typedef GridGraph<N, DirectedTag> Graph;
1407  typedef typename Graph::Edge Edge;
1408  typedef typename Graph::EdgeIt EdgeIt;
1409  typedef typename MultiArrayShape<N>::type CoordType;
1410 
1411  vigra_precondition(interpolatedImage.shape() == 2*g.shape()-CoordType(1),
1412  "edgeWeightsFromInterpolatedImage(): interpolated shape must be shape*2-1");
1413 
1414  for (EdgeIt iter(g); iter!=lemon::INVALID; ++iter)
1415  {
1416  const Edge edge(*iter);
1417  const CoordType uCoord(g.u(edge));
1418  const CoordType vCoord(g.v(edge));
1419  if (euclidean)
1420  {
1421  edgeWeights[edge] = norm(uCoord-vCoord) * interpolatedImage[uCoord+vCoord];
1422  }
1423  else
1424  {
1425  edgeWeights[edge] = interpolatedImage[uCoord+vCoord];
1426  }
1427  }
1428  }
1429 
1430  template<class GRAPH>
1431  struct ThreeCycle{
1432 
1433  typedef typename GRAPH::Node Node;
1434 
1435  ThreeCycle(const Node & a, const Node & b, const Node c){
1436  nodes_[0] = a;
1437  nodes_[1] = b;
1438  nodes_[2] = c;
1439  std::sort(nodes_, nodes_+3);
1440  }
1441  bool operator < (const ThreeCycle & other)const{
1442  if(nodes_[0] < other.nodes_[0]){
1443  return true;
1444  }
1445  else if(nodes_[0] == other.nodes_[0]){
1446  if(nodes_[1] < other.nodes_[1]){
1447  return true;
1448  }
1449  else if(nodes_[1] == other.nodes_[1]){
1450  if(nodes_[2] < other.nodes_[2]){
1451  return true;
1452  }
1453  else{
1454  return false;
1455  }
1456  }
1457  else{
1458  return false;
1459  }
1460  }
1461  else{
1462  return false;
1463  }
1464  }
1465 
1466  Node nodes_[3];
1467 
1468 
1469  };
1470 
1471 
1472  template<class GRAPH>
1473  void find3Cycles(
1474  const GRAPH & g,
1475  MultiArray<1, TinyVector<Int32, 3> > & cyclesArray
1476  ){
1477  typedef typename GRAPH::Node Node;
1478  typedef typename GRAPH::Edge Edge;
1479  typedef typename GRAPH::EdgeIt EdgeIt;
1480  typedef typename GRAPH::OutArcIt OutArcIt;
1481 
1482  typedef ThreeCycle<GRAPH> Cycle;
1483 
1484  std::set< Cycle > cycles;
1485  typedef typename std::set<Cycle>::const_iterator SetIter;
1486  for (EdgeIt iter(g); iter!=lemon::INVALID; ++iter){
1487  const Edge edge(*iter);
1488  const Node u = g.u(edge);
1489  const Node v = g.v(edge);
1490 
1491  // find a node n which is connected to u and v
1492  for(OutArcIt outArcIt(g,u); outArcIt!=lemon::INVALID;++outArcIt){
1493  const Node w = g.target(*outArcIt);
1494  if(w != v){
1495  const Edge e = g.findEdge(w,v);
1496  if(e != lemon::INVALID){
1497  // found cycle
1498  cycles.insert(Cycle(u, v, w));
1499  }
1500  }
1501  }
1502  }
1503  cyclesArray.reshape(TinyVector<UInt32,1>(cycles.size()));
1504  UInt32 i=0;
1505  for(SetIter iter=cycles.begin(); iter!=cycles.end(); ++iter){
1506 
1507  const Cycle & c = *iter;
1508  for(size_t j=0;j<3; ++j){
1509  cyclesArray(i)[j] = g.id(c.nodes_[j]);
1510  }
1511  ++i;
1512  }
1513  }
1514 
1515  template<class GRAPH>
1516  void find3CyclesEdges(
1517  const GRAPH & g,
1518  MultiArray<1, TinyVector<Int32, 3> > & cyclesArray
1519  ){
1520  typedef typename GRAPH::Node Node;
1521  typedef typename GRAPH::Edge Edge;
1522  typedef typename GRAPH::EdgeIt EdgeIt;
1523  typedef typename GRAPH::OutArcIt OutArcIt;
1524 
1525  typedef ThreeCycle<GRAPH> Cycle;
1526 
1527  std::set< Cycle > cycles;
1528  typedef typename std::set<Cycle>::const_iterator SetIter;
1529  for (EdgeIt iter(g); iter!=lemon::INVALID; ++iter){
1530  const Edge edge(*iter);
1531  const Node u = g.u(edge);
1532  const Node v = g.v(edge);
1533 
1534  // find a node n which is connected to u and v
1535  for(OutArcIt outArcIt(g,u); outArcIt!=lemon::INVALID;++outArcIt){
1536  const Node w = g.target(*outArcIt);
1537  if(w != v){
1538  const Edge e = g.findEdge(w,v);
1539  if(e != lemon::INVALID){
1540  // found cycle
1541  cycles.insert(Cycle(u, v, w));
1542  }
1543  }
1544  }
1545  }
1546  cyclesArray.reshape(TinyVector<UInt32,1>(cycles.size()));
1547  UInt32 i=0;
1548  for(SetIter iter=cycles.begin(); iter!=cycles.end(); ++iter){
1549 
1550  const Cycle & c = *iter;
1551  const Node u = c.nodes_[0];
1552  const Node v = c.nodes_[1];
1553  const Node w = c.nodes_[2];
1554 
1555  cyclesArray(i)[0] = g.id(g.findEdge(u, v));
1556  cyclesArray(i)[1] = g.id(g.findEdge(u, w));
1557  cyclesArray(i)[2] = g.id(g.findEdge(v, w));
1558  ++i;
1559  }
1560  }
1561 
1562 //@}
1563 
1564 } // namespace vigra
1565 
1566 #endif // VIGRA_GRAPH_ALGORITHMS_HXX
vigra::GridGraph< N, DirectedTag >::degree_size_type degree(typename vigra::GridGraph< N, DirectedTag >::vertex_descriptor const &v, vigra::GridGraph< N, DirectedTag > const &g)
Return total number of edges (incoming and outgoing) of vertex v (API: boost).
Definition: multi_gridgraph.hxx:2828
Define a grid graph in arbitrary dimensions.
Definition: multi_fwd.hxx:217
const PredecessorsMap & predecessors() const
get the predecessors node map (after a call of run)
Definition: graph_algorithms.hxx:442
void initMultiArrayBorder(...)
Write values to the specified border values in the array.
vigra::GridGraph< N, DirectedTag >::vertex_descriptor target(typename vigra::GridGraph< N, DirectedTag >::edge_descriptor const &e, vigra::GridGraph< N, DirectedTag > const &g)
Get a vertex descriptor for the end vertex of edge e in graph g (API: boost).
Definition: multi_gridgraph.hxx:2954
void run(Node const &start, Node const &stop, const WEIGHTS &weights, const Node &source, const Node &target=lemon::INVALID, WeightType maxDistance=NumericTraits< WeightType >::max())
run shortest path in a region of interest of a GridGraph.
Definition: graph_algorithms.hxx:356
void fillNodeMap(const G &g, A &a, const T &value)
fill a lemon node map
Definition: graph_algorithms.hxx:136
WeightType distance(const Node &target) const
get the distance to a rarget node (after a call of run)
Definition: graph_algorithms.hxx:452
linalg::TemporaryMatrix< T > exp(MultiArrayView< 2, T, C > const &v)
MultiArray< 2, MultiArrayIndex > ragFindEdges(const RAGGRAPH &rag, const GRAPH &graph, const RAGEDGES &affiliatedEdges, MultiArrayView< N, T > labelsArray, const typename RAGGRAPH::Node &node)
Find indices of points on the edges.
Definition: graph_algorithms.hxx:1269
Heap-based priority queue compatible to BucketQueue.
Definition: priority_queue.hxx:290
void shortestPathAStar(const GRAPH &graph, const typename GRAPH::Node &source, const typename GRAPH::Node &target, const WEIGHTS &weights, PREDECESSORS &predecessors, DISTANCE &distance, const HEURSTIC &heuristic)
Astar Shortest path search.
Definition: graph_algorithms.hxx:611
detail::GenericEdge< index_type > Edge
edge descriptor
Definition: adjacency_list_graph.hxx:249
Main MultiArray class containing the memory management.
Definition: multi_array.hxx:2422
const DistanceMap & distances() const
get the distances node map (after a call of run)
Definition: graph_algorithms.hxx:447
void felzenszwalbSegmentation(const GRAPH &graph, const EDGE_WEIGHTS &edgeWeights, const NODE_SIZE &nodeSizes, float k, NODE_LABEL_MAP &nodeLabeling, const int nodeNumStopCond=-1)
edge weighted watersheds Segmentataion
Definition: graph_algorithms.hxx:917
void edgeWeightsFromNodeWeights(const GridGraph< N, DirectedTag > &g, const NODEMAP &nodeWeights, EDGEMAP &edgeWeights, bool euclidean, FUNCTOR const &func)
create edge weights from node weights
Definition: graph_algorithms.hxx:1342
bool allLess(TinyVectorBase< V1, SIZE, D1, D2 > const &l, TinyVectorBase< V2, SIZE, D3, D4 > const &r)
pointwise less-than
Definition: tinyvector.hxx:1375
undirected adjacency list graph in the LEMON API
Definition: adjacency_list_graph.hxx:227
Definition: accessor.hxx:43
Node nodeFromId(const index_type id) const
Get node descriptor for given node ID i (API: LEMON). Return Node(lemon::INVALID) when the ID does no...
detail::SelectIntegerType< 64, detail::SignedIntTypes >::type Int64
64-bit signed int
Definition: sized_int.hxx:177
const DiscoveryOrder & discoveryOrder() const
get an array with all visited nodes, sorted by distance from source
Definition: graph_algorithms.hxx:437
void 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)
smooth node features of a graph
Definition: graph_algorithms.hxx:1125
FFTWComplex< R >::NormType norm(const FFTWComplex< R > &a)
norm (= magnitude)
Definition: fftw3.hxx:1037
const Node & source() const
get the source node
Definition: graph_algorithms.hxx:423
shortest path computer
Definition: graph_algorithms.hxx:297
const Graph & graph() const
get the graph
Definition: graph_algorithms.hxx:419
Definition: inspectimage.hxx:255
void copyNodeMap(const G &g, const A &a, B &b)
copy a lemon node map
Definition: graph_algorithms.hxx:117
const difference_type & shape() const
Definition: multi_array.hxx:1596
MultiArray & init(const U &init)
Definition: multi_array.hxx:2799
void edgeWeightsFromInterpolatedImage(const GridGraph< N, DirectedTag > &g, const MultiArrayView< N, T > &interpolatedImage, EDGEMAP &edgeWeights, bool euclidean=false)
create edge weights from an interpolated image
Definition: graph_algorithms.hxx:1400
Edge findEdge(const Node &a, const Node &b) const
Get a descriptor for the edge connecting vertices u and v, or lemon::INVALID if no such edge exists (...
vigra::GridGraph< N, DirectedTag >::vertex_descriptor source(typename vigra::GridGraph< N, DirectedTag >::edge_descriptor const &e, vigra::GridGraph< N, DirectedTag > const &g)
Get a vertex descriptor for the start vertex of edge e in graph g (API: boost).
Definition: multi_gridgraph.hxx:2943
const Node & target() const
get the target node
Definition: graph_algorithms.hxx:427
void runMultiSource(const WEIGHTS &weights, ITER source_begin, ITER source_end, const Node &target=lemon::INVALID, WeightType maxDistance=NumericTraits< WeightType >::max())
run shortest path with given edge weights from multiple sources.
Definition: graph_algorithms.hxx:391
void copyEdgeMap(const G &g, const A &a, B &b)
copy a lemon edge map
Definition: graph_algorithms.hxx:127
size_t pathLength(const NODE source, const NODE target, const PREDECESSORS &predecessors)
get the length in node units of a path
Definition: graph_algorithms.hxx:591
Class for fixed size vectors.This class contains an array of size SIZE of the specified VALUETYPE...
Definition: accessor.hxx:940
void fillEdgeMap(const G &g, A &a, const T &value)
fill a lemon edge map
Definition: graph_algorithms.hxx:145
void 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)
make a region adjacency graph from a graph and labels w.r.t. that graph
Definition: graph_algorithms.hxx:165
void edgeSort(const GRAPH &g, const WEIGHTS &weights, const COMPERATOR &comperator, std::vector< typename GRAPH::Edge > &sortedEdges)
get a vector of Edge descriptors
Definition: graph_algorithms.hxx:98
std::pair< typename vigra::GridGraph< N, DirectedTag >::edge_descriptor, bool > edge(typename vigra::GridGraph< N, DirectedTag >::vertex_descriptor const &u, typename vigra::GridGraph< N, DirectedTag >::vertex_descriptor const &v, vigra::GridGraph< N, DirectedTag > const &g)
Return the pair (edge_descriptor, true) when an edge between vertices u and v exists, or (lemon::INVALID, false) otherwise (API: boost).
Definition: multi_gridgraph.hxx:2990
void 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 &ragGtQt)
Definition: graph_algorithms.hxx:1199
MultiArrayShape< 2 >::type Shape2
shape type for MultiArray<2, T>
Definition: multi_shape.hxx:254
void 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)
edge weighted watersheds Segmentataion
Definition: graph_algorithms.hxx:893
bool operator<(FixedPoint< IntBits1, FracBits1 > l, FixedPoint< IntBits2, FracBits2 > r)
less than
Definition: fixedpoint.hxx:512
void edgeWeightedWatershedsSegmentation(const GRAPH &g, const EDGE_WEIGHTS &edgeWeights, const SEEDS &seeds, LABELS &labels)
edge weighted watersheds Segmentataion
Definition: graph_algorithms.hxx:872
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())
run shortest path with given edge weights from multiple sources.
Definition: graph_algorithms.hxx:406
detail::SelectIntegerType< 32, detail::UnsignedIntTypes >::type UInt32
32-bit unsigned int
Definition: sized_int.hxx:183
Base class for, and view to, vigra::MultiArray.
Definition: multi_array.hxx:652
Node v(Edge const &e) const
Get the end node of the given edge e (API: LEMON, the boost::graph API provides the free function bo...
Definition: multi_gridgraph.hxx:2390
void 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)
smooth node features of a graph
Definition: graph_algorithms.hxx:1100
bool allLessEqual(TinyVectorBase< V1, SIZE, D1, D2 > const &l, TinyVectorBase< V2, SIZE, D3, D4 > const &r)
pointwise less-equal
Definition: tinyvector.hxx:1399
void run(const WEIGHTS &weights, const Node &source, const Node &target=lemon::INVALID, WeightType maxDistance=NumericTraits< WeightType >::max())
run shortest path given edge weights
Definition: graph_algorithms.hxx:334
Node u(Edge const &e) const
Get the start node of the given edge e (API: LEMON, the boost::graph API provides the free function ...
Definition: multi_gridgraph.hxx:2382
void reRun(const WEIGHTS &weights, const Node &source, const Node &target=lemon::INVALID, WeightType maxDistance=NumericTraits< WeightType >::max())
run shortest path again with given edge weights
Definition: graph_algorithms.hxx:377
detail_adjacency_list_graph::ItemIter< GraphType, Edge > EdgeIt
edge iterator
Definition: adjacency_list_graph.hxx:253
ShortestPathDijkstra(const Graph &g)
\ brief constructor from graph
Definition: graph_algorithms.hxx:313
bool hasTarget() const
check if explicit target is given
Definition: graph_algorithms.hxx:432

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

html generated using doxygen and Python
vigra 1.11.0