[ 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 
48 
49 /*vigra*/
50 #include "graphs.hxx"
51 #include "graph_generalization.hxx"
52 #include "multi_gridgraph.hxx"
53 #include "priority_queue.hxx"
54 #include "union_find.hxx"
55 #include "adjacency_list_graph.hxx"
56 #include "graph_maps.hxx"
57 #include "functorexpression.hxx"
58 #include "array_vector.hxx"
59 
60 namespace vigra{
61 
62 /** \addtogroup GraphDataStructures
63 */
64 //@{
65 
66  namespace detail_graph_algorithms{
67  template <class GRAPH_MAP,class COMPERATOR>
68  struct GraphItemCompare
69  {
70 
71  GraphItemCompare(const GRAPH_MAP & map,const COMPERATOR & comperator)
72  : map_(map),
73  comperator_(comperator){
74 
75  }
76 
77  template<class KEY>
78  bool operator()(const KEY & a, const KEY & b) const{
79  return comperator_(map_[a],map_[b]);
80  }
81 
82  const GRAPH_MAP & map_;
83  const COMPERATOR & comperator_;
84  };
85  } // namespace detail_graph_algorithms
86 
87  /// \brief get a vector of Edge descriptors
88  ///
89  /// Sort the Edge descriptors given weights
90  /// and a comperator
91  template<class GRAPH,class WEIGHTS,class COMPERATOR>
92  void edgeSort(
93  const GRAPH & g,
94  const WEIGHTS & weights,
95  const COMPERATOR & comperator,
96  std::vector<typename GRAPH::Edge> & sortedEdges
97  ){
98  sortedEdges.resize(g.edgeNum());
99  size_t c=0;
100  for(typename GRAPH::EdgeIt e(g);e!=lemon::INVALID;++e){
101  sortedEdges[c]=*e;
102  ++c;
103  }
104  detail_graph_algorithms::GraphItemCompare<WEIGHTS,COMPERATOR> edgeComperator(weights,comperator);
105  std::sort(sortedEdges.begin(),sortedEdges.end(),edgeComperator);
106  }
107 
108 
109  /// \brief copy a lemon node map
110  template<class G,class A,class B>
111  void copyNodeMap(const G & g,const A & a ,B & b){
112  typename G::NodeIt iter(g);
113  while(iter!=lemon::INVALID){
114  b[*iter]=a[*iter];
115  ++iter;
116  }
117 
118  }
119  /// \brief copy a lemon edge map
120  template<class G,class A,class B>
121  void copyEdgeMap(const G & g,const A & a ,B & b){
122  typename G::EdgeIt iter(g);
123  while(iter!=lemon::INVALID){
124  b[*iter]=a[*iter];
125  ++iter;
126  }
127  }
128  /// \brief fill a lemon node map
129  template<class G,class A,class T>
130  void fillNodeMap(const G & g, A & a, const T & value){
131  typename G::NodeIt iter(g);
132  while(iter!=lemon::INVALID){
133  a[*iter]=value;
134  ++iter;
135  }
136  }
137  /// \brief fill a lemon edge map
138  template<class G,class A,class T>
139  void fillEdgeMap(const G & g,A & a ,const T & value){
140  typename G::EdgeIt iter(g);
141  while(iter!=lemon::INVALID){
142  a[*iter]=value;
143  ++iter;
144  }
145  }
146 
147  /// \brief make a region adjacency graph from a graph and labels w.r.t. that graph
148  ///
149  /// \param graphIn : input graph
150  /// \param labels : labels w.r.t. graphIn
151  /// \param[out] rag : region adjacency graph
152  /// \param[out] affiliatedEdges : a vector of edges of graphIn for each edge in rag
153  /// \param ignoreLabel : optional label to ignore (default: -1 means no label will be ignored)
154  ///
155  template<
156  class GRAPH_IN,
157  class GRAPH_IN_NODE_LABEL_MAP
158  >
160  GRAPH_IN graphIn,
161  GRAPH_IN_NODE_LABEL_MAP labels,
162  AdjacencyListGraph & rag,
163  typename AdjacencyListGraph:: template EdgeMap< std::vector<typename GRAPH_IN::Edge> > & affiliatedEdges,
164  const Int64 ignoreLabel=-1
165  ){
166  rag=AdjacencyListGraph();
167  typedef typename GraphMapTypeTraits<GRAPH_IN_NODE_LABEL_MAP>::Value LabelType;
168  typedef GRAPH_IN GraphIn;
169  typedef AdjacencyListGraph GraphOut;
170 
171  typedef typename GraphIn::Edge EdgeGraphIn;
172  typedef typename GraphIn::NodeIt NodeItGraphIn;
173  typedef typename GraphIn::EdgeIt EdgeItGraphIn;
174 
175  typedef typename GraphOut::Edge EdgeGraphOut;
176  for(NodeItGraphIn iter(graphIn);iter!=lemon::INVALID;++iter){
177  const LabelType l=labels[*iter];
178  if(ignoreLabel==-1 || static_cast<Int64>(l)!=ignoreLabel)
179  rag.addNode(l);
180  }
181 
182  // add al edges
183  for(EdgeItGraphIn e(graphIn);e!=lemon::INVALID;++e){
184  const EdgeGraphIn edge(*e);
185  const LabelType lu = labels[graphIn.u(edge)];
186  const LabelType lv = labels[graphIn.v(edge)];
187  if( lu!=lv && ( ignoreLabel==-1 || (static_cast<Int64>(lu)!=ignoreLabel && static_cast<Int64>(lv)!=ignoreLabel) ) ){
188  // if there is an edge between lu and lv no new edge will be added
189  rag.addEdge( rag.nodeFromId(lu),rag.nodeFromId(lv));
190  }
191  }
192  // SET UP HYPEREDGES
193  affiliatedEdges.assign(rag);
194  // add edges
195  for(EdgeItGraphIn e(graphIn);e!=lemon::INVALID;++e){
196  const EdgeGraphIn edge(*e);
197  const LabelType lu = labels[graphIn.u(edge)];
198  const LabelType lv = labels[graphIn.v(edge)];
199  //std::cout<<"edge between ?? "<<lu<<" "<<lv<<"\n";
200  if( lu!=lv && ( ignoreLabel==-1 || (static_cast<Int64>(lu)!=ignoreLabel && static_cast<Int64>(lv)!=ignoreLabel) ) ){
201  //std::cout<<"find edge between "<<lu<<" "<<lv<<"\n";
202  EdgeGraphOut ragEdge= rag.findEdge(rag.nodeFromId(lu),rag.nodeFromId(lv));
203  //std::cout<<"invalid?"<<bool(ragEdge==lemon::INVALID)<<" id "<<rag.id(ragEdge)<<"\n";
204  affiliatedEdges[ragEdge].push_back(edge);
205  //std::cout<<"write done\n";
206  }
207  }
208  }
209 
210  /// \brief shortest path computer
211  template<class GRAPH,class WEIGHT_TYPE>
213  public:
214  typedef GRAPH Graph;
215 
216  typedef typename Graph::Node Node;
217  typedef typename Graph::NodeIt NodeIt;
218  typedef typename Graph::Edge Edge;
219  typedef typename Graph::OutArcIt OutArcIt;
220 
221  typedef WEIGHT_TYPE WeightType;
223  typedef typename Graph:: template NodeMap<Node> PredecessorsMap;
224  typedef typename Graph:: template NodeMap<WeightType> DistanceMap;
226 
227  /// \ brief constructor from graph
228  ShortestPathDijkstra(const Graph & g)
229  : graph_(g),
230  pq_(g.maxNodeId()+1),
231  predMap_(g),
232  distMap_(g)
233  {
234  }
235 
236  /// \brief run shortest path given edge weights
237  ///
238  /// \param weights : edge weights encoding the distance between adjacent nodes (must be non-negative)
239  /// \param source : source node where shortest path should start
240  /// \param target : target node where shortest path should stop. If target is not given
241  /// or <tt>INVALID</tt>, the shortest path from source to all reachable nodes is computed
242  /// \param maxDistance : path search is terminated when the path length exceeds <tt>maxDistance</tt>
243  ///
244  /// When a valid \a target is unreachable from \a source (either because the graph is disconnected
245  /// or \a maxDistance is exceeded), it is set to <tt>lemon::INVALID</tt>. In contrast, if \a target
246  /// was <tt>lemon::INVALID</tt> at the beginning, it will always be set to the last node
247  /// visited in the search.
248  template<class WEIGHTS>
249  void run(const WEIGHTS & weights, const Node & source,
250  const Node & target = lemon::INVALID,
251  WeightType maxDistance=NumericTraits<WeightType>::max())
252  {
253  this->initializeMaps(source);
254  runImpl(weights, target, maxDistance);
255  }
256 
257  /// \brief run shortest path in a region of interest of a \ref GridGraph.
258  ///
259  /// \param start : first point in the desired ROI.
260  /// \param stop : beyond the last point in the desired ROI (i.e. exclusive)
261  /// \param weights : edge weights encoding the distance between adjacent nodes (must be non-negative)
262  /// \param source : source node where shortest path should start
263  /// \param target : target node where shortest path should stop. If target is not given
264  /// or <tt>INVALID</tt>, the shortest path from source to all reachable nodes is computed
265  /// \param maxDistance : path search is terminated when the path length exceeds <tt>maxDistance</tt>
266  ///
267  /// This version of <tt>run()</tt> restricts the path search to the ROI <tt>[start, stop)</tt> and only
268  /// works for instances of \ref GridGraph. Otherwise, it is identical to the standard <tt>run()</tt>
269  /// function.
270  template<class WEIGHTS>
271  void run(Node const & start, Node const & stop,
272  const WEIGHTS & weights, const Node & source,
273  const Node & target = lemon::INVALID,
274  WeightType maxDistance=NumericTraits<WeightType>::max())
275  {
276  vigra_precondition(allLessEqual(start, source) && allLess(source, stop),
277  "ShortestPathDijkstra::run(): source is not within ROI");
278  vigra_precondition(target == lemon::INVALID ||
279  (allLessEqual(start, target) && allLess(target, stop)),
280  "ShortestPathDijkstra::run(): target is not within ROI");
281  this->initializeMaps(source, start, stop);
282  runImpl(weights, target, maxDistance);
283  }
284 
285  /// \brief run shortest path again with given edge weights
286  ///
287  /// This only differs from standard <tt>run()</tt> by initialization: Instead of resetting
288  /// the entire graph, this only resets the nodes that have been visited in the
289  /// previous run, i.e. the contents of the array <tt>discoveryOrder()</tt>.
290  /// This will be much faster if only a small fraction of the nodes has to be reset.
291  template<class WEIGHTS>
292  void reRun(const WEIGHTS & weights, const Node & source,
293  const Node & target = lemon::INVALID,
294  WeightType maxDistance=NumericTraits<WeightType>::max())
295  {
296  this->reInitializeMaps(source);
297  runImpl(weights, target, maxDistance);
298  }
299 
300  /// \brief run shortest path with given edge weights from multiple sources.
301  ///
302  /// This is otherwise identical to standard <tt>run()</tt>, except that
303  /// <tt>source()</tt> returns <tt>lemon::INVALID</tt> after path search finishes.
304  template<class WEIGHTS, class ITER>
305  void
306  runMultiSource(const WEIGHTS & weights, ITER source_begin, ITER source_end,
307  const Node & target = lemon::INVALID,
308  WeightType maxDistance=NumericTraits<WeightType>::max())
309  {
310  this->initializeMapsMultiSource(source_begin, source_end);
311  runImpl(weights, target, maxDistance);
312  }
313 
314  /// \brief get the graph
315  const Graph & graph()const{
316  return graph_;
317  }
318  /// \brief get the source node
319  const Node & source()const{
320  return source_;
321  }
322  /// \brief get the target node
323  const Node & target()const{
324  return target_;
325  }
326 
327  /// \brief check if explicit target is given
328  bool hasTarget()const{
329  return target_!=lemon::INVALID;
330  }
331 
332  /// \brief get an array with all visited nodes, sorted by distance from source
334  return discoveryOrder_;
335  }
336 
337  /// \brief get the predecessors node map (after a call of run)
338  const PredecessorsMap & predecessors()const{
339  return predMap_;
340  }
341 
342  /// \brief get the distances node map (after a call of run)
343  const DistanceMap & distances()const{
344  return distMap_;
345  }
346 
347  /// \brief get the distance to a rarget node (after a call of run)
348  WeightType distance(const Node & target)const{
349  return distMap_[target];
350  }
351 
352 
353  private:
354 
355  template<class WEIGHTS>
356  void runImpl(const WEIGHTS & weights,
357  const Node & target = lemon::INVALID,
358  WeightType maxDistance=NumericTraits<WeightType>::max())
359  {
360  target_ = lemon::INVALID;
361  while(!pq_.empty() ){ //&& !finished){
362  const Node topNode(graph_.nodeFromId(pq_.top()));
363  if(distMap_[topNode] > maxDistance)
364  break; // distance threshold exceeded
365  pq_.pop();
366  discoveryOrder_.push_back(topNode);
367  if(topNode == target)
368  break;
369  // loop over all neigbours
370  for(OutArcIt outArcIt(graph_,topNode);outArcIt!=lemon::INVALID;++outArcIt){
371  const Node otherNode = graph_.target(*outArcIt);
372  const size_t otherNodeId = graph_.id(otherNode);
373 
374  if(pq_.contains(otherNodeId)){
375  const Edge edge(*outArcIt);
376  const WeightType currentDist = distMap_[otherNode];
377  const WeightType alternativeDist = distMap_[topNode]+weights[edge];
378  if(alternativeDist<currentDist){
379  pq_.push(otherNodeId,alternativeDist);
380  distMap_[otherNode]=alternativeDist;
381  predMap_[otherNode]=topNode;
382  }
383  }
384  else if(predMap_[otherNode]==lemon::INVALID){
385  const Edge edge(*outArcIt);
386  const WeightType initialDist = distMap_[topNode]+weights[edge];
387  if(initialDist<=maxDistance)
388  {
389  pq_.push(otherNodeId,initialDist);
390  distMap_[otherNode]=initialDist;
391  predMap_[otherNode]=topNode;
392  }
393  }
394  }
395  }
396  while(!pq_.empty() ){
397  const Node topNode(graph_.nodeFromId(pq_.top()));
398  predMap_[topNode]=lemon::INVALID;
399  pq_.pop();
400  }
401  if(target == lemon::INVALID || discoveryOrder_.back() == target)
402  target_ = discoveryOrder_.back(); // Means that target was reached. If, to the contrary, target
403  // was unreachable within maxDistance, target_ remains INVALID.
404  }
405 
406  void initializeMaps(Node const & source){
407  for(NodeIt n(graph_); n!=lemon::INVALID; ++n){
408  const Node node(*n);
409  predMap_[node]=lemon::INVALID;
410  }
411  distMap_[source]=static_cast<WeightType>(0.0);
412  predMap_[source]=source;
413  discoveryOrder_.clear();
414  pq_.push(graph_.id(source),0.0);
415  source_=source;
416  }
417 
418  void initializeMaps(Node const & source,
419  Node const & start, Node const & stop)
420  {
421  Node left_border = min(start, Node(1)),
422  right_border = min(predMap_.shape()-stop, Node(1)),
423  DONT_TOUCH = Node(lemon::INVALID) - Node(1);
424 
425  initMultiArrayBorder(predMap_.subarray(start-left_border, stop+right_border),
426  left_border, right_border, DONT_TOUCH);
427  predMap_.subarray(start, stop) = lemon::INVALID;
428  predMap_[source]=source;
429 
430  distMap_[source]=static_cast<WeightType>(0.0);
431  discoveryOrder_.clear();
432  pq_.push(graph_.id(source),0.0);
433  source_=source;
434  }
435 
436  template <class ITER>
437  void initializeMapsMultiSource(ITER source, ITER source_end){
438  for(NodeIt n(graph_); n!=lemon::INVALID; ++n){
439  const Node node(*n);
440  predMap_[node]=lemon::INVALID;
441  }
442  discoveryOrder_.clear();
443  for( ; source != source_end; ++source)
444  {
445  distMap_[*source]=static_cast<WeightType>(0.0);
446  predMap_[*source]=*source;
447  pq_.push(graph_.id(*source),0.0);
448  }
449  source_=lemon::INVALID;
450  }
451 
452  void reInitializeMaps(Node const & source){
453  for(unsigned int n=0; n<discoveryOrder_.size(); ++n){
454  predMap_[discoveryOrder_[n]]=lemon::INVALID;
455  }
456  distMap_[source]=static_cast<WeightType>(0.0);
457  predMap_[source]=source;
458  discoveryOrder_.clear();
459  pq_.push(graph_.id(source),0.0);
460  source_=source;
461  }
462 
463  const Graph & graph_;
464  PqType pq_;
465  PredecessorsMap predMap_;
466  DistanceMap distMap_;
467  DiscoveryOrder discoveryOrder_;
468 
469  Node source_;
470  Node target_;
471  };
472 
473  /// \brief get the length in node units of a path
474  template<class NODE,class PREDECESSORS>
475  size_t pathLength(
476  const NODE source,
477  const NODE target,
478  const PREDECESSORS & predecessors
479  ){
480  if(predecessors[target]==lemon::INVALID)
481  return 0;
482  else{
483  NODE currentNode = target;
484  size_t length=1;
485  while(currentNode!=source){
486  currentNode=predecessors[currentNode];
487  length+=1;
488  }
489  return length;
490  }
491  }
492 
493 
494 
495  /// \brief Astar Shortest path search
496  template<class GRAPH,class WEIGHTS,class PREDECESSORS,class DISTANCE,class HEURSTIC>
498  const GRAPH & graph,
499  const typename GRAPH::Node & source,
500  const typename GRAPH::Node & target,
501  const WEIGHTS & weights,
502  PREDECESSORS & predecessors,
503  DISTANCE & distance,
504  const HEURSTIC & heuristic
505  ){
506 
507  typedef GRAPH Graph;
508  typedef typename Graph::Edge Edge;
509  typedef typename Graph::Node Node;
510  typedef typename Graph::NodeIt NodeIt;
511  typedef typename Graph::OutArcIt OutArcIt;
512  typedef typename DISTANCE::value_type DistanceType;
513 
514  typename GRAPH:: template NodeMap<bool> closedSet(graph);
515  vigra::ChangeablePriorityQueue<typename WEIGHTS::value_type> estimatedDistanceOpenSet(graph.maxNodeId()+1);
516  // initialize
517  for(NodeIt n(graph);n!=lemon::INVALID;++n){
518  const Node node(*n);
519  closedSet[node]=false;
520  distance[node]=std::numeric_limits<DistanceType>::infinity();
521  predecessors[node]=lemon::INVALID;
522  }
523  // distance and estimated distance for start node
524  distance[source]=static_cast<DistanceType>(0.0);
525  estimatedDistanceOpenSet.push(graph.id(source),heuristic(source,target));
526 
527  // while any nodes left in openSet
528  while(!estimatedDistanceOpenSet.empty()){
529 
530  // get the node with the lpwest estimated distance in the open set
531  const Node current = graph.nodeFromId(estimatedDistanceOpenSet.top());
532 
533  // reached target?
534  if(current==target)
535  break;
536 
537  // remove current from openSet
538  // add current to closedSet
539  estimatedDistanceOpenSet.pop();
540  closedSet[current]=true;
541 
542  // iterate over neigbours of current
543  for(OutArcIt outArcIt(graph,current);outArcIt!=lemon::INVALID;++outArcIt){
544 
545  // get neigbour node and id
546  const Node neighbour = graph.target(*outArcIt);
547  const size_t neighbourId = graph.id(neighbour);
548 
549  // if neighbour is not yet in closedSet
550  if(!closedSet[neighbour]){
551 
552  // get edge between current and neigbour
553  const Edge edge(*outArcIt);
554 
555  // get tentative score
556  const DistanceType tenativeScore = distance[current] + weights[edge];
557 
558  // neighbour NOT in openSet OR tentative score better than the current distance
559  if(!estimatedDistanceOpenSet.contains(neighbourId) || tenativeScore < distance[neighbour] ){
560  // set predecessors and distance
561  predecessors[neighbour]=current;
562  distance[neighbour]=tenativeScore;
563 
564  // update the estimated cost from neighbour to target
565  // ( and neigbour will be (re)-added to openSet)
566  estimatedDistanceOpenSet.push(neighbourId,distance[neighbour]+heuristic(neighbour,target));
567  }
568  }
569  }
570  }
571  }
572 
573 
574  namespace detail_watersheds_segmentation{
575 
576  struct RawPriorityFunctor{
577  template<class L, class T>
578  T operator()(const L label,const T priority)const{
579  return priority;
580  }
581  };
582 
583  template<class PRIORITY_TYPE,class LABEL_TYPE>
584  struct CarvingFunctor{
585  CarvingFunctor(const LABEL_TYPE backgroundLabel,const PRIORITY_TYPE & factor)
586  : backgroundLabel_(backgroundLabel),
587  factor_(factor){
588  }
589  PRIORITY_TYPE operator()(const LABEL_TYPE label,const PRIORITY_TYPE priority)const{
590  return (label==backgroundLabel_ ? priority*factor_ : priority);
591  }
592  LABEL_TYPE backgroundLabel_;
593  PRIORITY_TYPE factor_;
594  };
595 
596 
597  template<class GRAPH,class EDGE_WEIGHTS,class SEEDS,class PRIORITY_MANIP_FUNCTOR,class LABELS>
598  void edgeWeightedWatershedsSegmentationImpl(
599  const GRAPH & g,
600  const EDGE_WEIGHTS & edgeWeights,
601  const SEEDS & seeds,
602  PRIORITY_MANIP_FUNCTOR & priorManipFunctor,
603  LABELS & labels
604  ){
605  typedef GRAPH Graph;
606  typedef typename Graph::Edge Edge;
607  typedef typename Graph::Node Node;
608  typedef typename Graph::NodeIt NodeIt;
609  typedef typename Graph::OutArcIt OutArcIt;
610 
611  typedef typename EDGE_WEIGHTS::Value WeightType;
612  typedef typename LABELS::Value LabelType;
613  typedef typename Graph:: template NodeMap<bool> NodeBoolMap;
614  typedef PriorityQueue<Node,WeightType,true> PQ;
615 
616  PQ pq;
617  NodeBoolMap inPQ(g);
618  copyNodeMap(g,seeds,labels);
619  fillNodeMap(g,inPQ,false);
620 
621  bool anySeed=false;
622  for(NodeIt n(g);n!=lemon::INVALID;++n){
623  const Node node(*n);
624  if(labels[node]!=static_cast<LabelType>(0)){
625  anySeed=true;
626  for(OutArcIt a(g,node);a!=lemon::INVALID;++a){
627  const Edge edge(*a);
628  const Node neigbour=g.target(*a);
629  //std::cout<<"n- node "<<g.id(neigbour)<<"\n";
630  if(labels[neigbour]==static_cast<LabelType>(0) && !inPQ[neigbour]){
631  const WeightType priority = priorManipFunctor(labels[node],edgeWeights[edge]);
632  pq.push(neigbour,priority);
633  inPQ[neigbour]=true;
634  }
635  }
636  }
637  }
638 
639 
640  if(anySeed){
641 
642  while(!pq.empty()){
643  const Node node = pq.top();
644  const LabelType label = labels[node];
645  //std::cout<<"node "<<g.id(node)<<" with label "<<label<<"\n";
646  if(label!=0){
647  throw std::runtime_error("seems like there are no seeds at all");
648  }
649 
650  pq.pop();
651  bool moreThanOneLabel = false;
652  LabelType labelFound = 0 ;
653  for(OutArcIt a(g,node);a!=lemon::INVALID;++a){
654  const Edge edge(*a);
655  const Node neigbour=g.target(*a);
656  if(labels[neigbour]!=static_cast<LabelType>(0)){
657  if(labelFound==0){
658  labelFound=labels[neigbour];
659  }
660  else{
661  moreThanOneLabel=true;
662  break;
663  }
664  }
665  }
666  if(labelFound!=0 && !moreThanOneLabel ){
667  labels[node]=labelFound;
668  for(OutArcIt a(g,node);a!=lemon::INVALID;++a){
669  const Edge edge(*a);
670  const Node neigbour=g.target(*a);
671  if(labels[neigbour]==static_cast<LabelType>(0)){
672  if(!inPQ[neigbour]){
673  const WeightType priority = priorManipFunctor(labelFound,edgeWeights[edge]);
674  pq.push(neigbour,edgeWeights[edge]);
675  inPQ[neigbour]=true;
676  }
677  }
678  }
679  }
680  }
681 
682 
683  for(NodeIt n(g);n!=lemon::INVALID;++n){
684  const Node node(*n);
685  if(labels[node]==static_cast<LabelType>(0)){
686 
687  WeightType minWeight = std::numeric_limits<WeightType>::infinity();
688  LabelType minWeightLabel = static_cast<LabelType>(0);
689  for(OutArcIt a(g,node);a!=lemon::INVALID;++a){
690  const Edge edge(*a);
691  const Node neigbour=g.target(*a);
692  const WeightType priority = priorManipFunctor(labels[neigbour],edgeWeights[edge]);
693  if(labels[neigbour]!=0 && priority<minWeight){
694  minWeight=priority;
695  minWeightLabel=labels[neigbour];
696  }
697  }
698  labels[node]=minWeightLabel;
699  }
700  }
701  }
702  }
703 
704  } // end namespace detail_watersheds_segmentation
705 
706 
707  /// \brief edge weighted watersheds Segmentataion
708  ///
709  /// \param g: input graph
710  /// \param edgeWeights : edge weights / edge indicator
711  /// \param seeds : seed must be non empty!
712  /// \param[out] labels : resulting nodeLabeling (not necessarily dense)
713  template<class GRAPH,class EDGE_WEIGHTS,class SEEDS,class LABELS>
715  const GRAPH & g,
716  const EDGE_WEIGHTS & edgeWeights,
717  const SEEDS & seeds,
718  LABELS & labels
719  ){
720  detail_watersheds_segmentation::RawPriorityFunctor f;
721  detail_watersheds_segmentation::edgeWeightedWatershedsSegmentationImpl(g,edgeWeights,seeds,f,labels);
722  }
723 
724  /// \brief edge weighted watersheds Segmentataion
725  ///
726  /// \param g: input graph
727  /// \param edgeWeights : edge weights / edge indicator
728  /// \param seeds : seed must be non empty!
729  /// \param backgroundLabel : which label is background
730  /// \param backgroundBias : bias for background
731  /// \param[out] labels : resulting nodeLabeling (not necessarily dense)
732  template<class GRAPH,class EDGE_WEIGHTS,class SEEDS,class LABELS>
734  const GRAPH & g,
735  const EDGE_WEIGHTS & edgeWeights,
736  const SEEDS & seeds,
737  const typename LABELS::Value backgroundLabel,
738  const typename EDGE_WEIGHTS::Value backgroundBias,
739  LABELS & labels
740  ){
741  typedef typename EDGE_WEIGHTS::Value WeightType;
742  typedef typename LABELS::Value LabelType;
743  detail_watersheds_segmentation::CarvingFunctor<WeightType,LabelType> f(backgroundLabel,backgroundBias);
744  detail_watersheds_segmentation::edgeWeightedWatershedsSegmentationImpl(g,edgeWeights,seeds,f,labels);
745  }
746 
747 
748  /// \brief edge weighted watersheds Segmentataion
749  ///
750  /// \param graph: input graph
751  /// \param edgeWeights : edge weights / edge indicator
752  /// \param nodeSizes : size of each node
753  /// \param k : free parameter of felzenszwalb algorithm
754  /// \param[out] nodeLabeling : nodeLabeling (not necessarily dense)
755  /// \param nodeNumStopCond : optional stopping condition
756  template< class GRAPH , class EDGE_WEIGHTS, class NODE_SIZE,class NODE_LABEL_MAP>
758  const GRAPH & graph,
759  const EDGE_WEIGHTS & edgeWeights,
760  const NODE_SIZE & nodeSizes,
761  float k,
762  NODE_LABEL_MAP & nodeLabeling,
763  const int nodeNumStopCond = -1
764  ){
765  typedef GRAPH Graph;
766  typedef typename Graph::Edge Edge;
767  typedef typename Graph::Node Node;
768 
769  typedef typename EDGE_WEIGHTS::Value WeightType;
770  typedef typename EDGE_WEIGHTS::Value NodeSizeType;
771  typedef typename Graph:: template NodeMap<WeightType> NodeIntDiffMap;
772  typedef typename Graph:: template NodeMap<NodeSizeType> NodeSizeAccMap;
773 
774  // initalize node size map and internal diff map
775  NodeIntDiffMap internalDiff(graph);
776  NodeSizeAccMap nodeSizeAcc(graph);
777  copyNodeMap(graph,nodeSizes,nodeSizeAcc);
778  fillNodeMap(graph,internalDiff,static_cast<WeightType>(0.0));
779 
780 
781 
782  // initlaize internal node diff map
783 
784  // sort the edges by their weights
785  std::vector<Edge> sortedEdges;
786  std::less<WeightType> comperator;
787  edgeSort(graph,edgeWeights,comperator,sortedEdges);
788 
789  // make the ufd
790  UnionFindArray<UInt64> ufdArray(graph.maxNodeId()+1);
791 
792 
793  size_t nodeNum = graph.nodeNum();
794 
795 
796  while(true){
797  // iterate over edges is the sorted order
798  for(size_t i=0;i<sortedEdges.size();++i){
799  const Edge e = sortedEdges[i];
800  const size_t rui = ufdArray.findIndex(graph.id(graph.u(e)));
801  const size_t rvi = ufdArray.findIndex(graph.id(graph.v(e)));
802  const Node ru = graph.nodeFromId(rui);
803  const Node rv = graph.nodeFromId(rvi);
804  if(rui!=rvi){
805 
806  //check if to merge or not ?
807  const WeightType w = edgeWeights[e];
808  const NodeSizeType sizeRu = nodeSizeAcc[ru];
809  const NodeSizeType sizeRv = nodeSizeAcc[rv];
810  const WeightType tauRu = static_cast<WeightType>(k)/static_cast<WeightType>(sizeRu);
811  const WeightType tauRv = static_cast<WeightType>(k)/static_cast<WeightType>(sizeRv);
812  const WeightType minIntDiff = std::min(internalDiff[ru]+tauRu,internalDiff[rv]+tauRv);
813  if(w<=minIntDiff){
814  // do merge
815  ufdArray.makeUnion(rui,rvi);
816  --nodeNum;
817  // update size and internal difference
818  const size_t newRepId = ufdArray.findIndex(rui);
819  const Node newRepNode = graph.nodeFromId(newRepId);
820  internalDiff[newRepNode]=w;
821  nodeSizeAcc[newRepNode] = sizeRu+sizeRv;
822  }
823  }
824  if(nodeNum==nodeNumStopCond){
825  break;
826  }
827  }
828  if(nodeNumStopCond==-1){
829  break;
830  }
831  else{
832  if(nodeNum>nodeNumStopCond){
833  k *= 1.2f;
834  }
835  else{
836  break;
837  }
838  }
839  }
840  ufdArray.makeContiguous();
841  for(typename GRAPH::NodeIt n(graph);n!=lemon::INVALID;++n){
842  const Node node(*n);
843  nodeLabeling[node]=ufdArray.findLabel(graph.id(node));
844  }
845  }
846 
847 
848 
849 
850  namespace detail_graph_smoothing{
851 
852  template<
853  class GRAPH,
854  class NODE_FEATURES_IN,
855  class EDGE_WEIGHTS,
856  class WEIGHTS_TO_SMOOTH_FACTOR,
857  class NODE_FEATURES_OUT
858  >
859  void graphSmoothingImpl(
860  const GRAPH & g,
861  const NODE_FEATURES_IN & nodeFeaturesIn,
862  const EDGE_WEIGHTS & edgeWeights,
863  WEIGHTS_TO_SMOOTH_FACTOR & weightsToSmoothFactor,
864  NODE_FEATURES_OUT & nodeFeaturesOut
865 
866  ){
867  typedef GRAPH Graph;
868  typedef typename Graph::Edge Edge;
869  typedef typename Graph::Node Node;
870  typedef typename Graph::NodeIt NodeIt;
871  typedef typename Graph::OutArcIt OutArcIt;
872 
873  typedef typename NODE_FEATURES_IN::Value NodeFeatureInValue;
874  typedef typename NODE_FEATURES_OUT::Reference NodeFeatureOutRef;
875  typedef typename EDGE_WEIGHTS::ConstReference SmoothFactorType;
876 
877 
878  //fillNodeMap(g, nodeFeaturesOut, typename NODE_FEATURES_OUT::value_type(0.0));
879 
880  for(NodeIt n(g);n!=lemon::INVALID;++n){
881 
882  const Node node(*n);
883 
884  NodeFeatureInValue featIn = nodeFeaturesIn[node];
885  NodeFeatureOutRef featOut = nodeFeaturesOut[node];
886 
887  featOut=0;
888  float weightSum = 0.0;
889  size_t degree = 0;
890  for(OutArcIt a(g,node);a!=lemon::INVALID;++a){
891  const Edge edge(*a);
892  const Node neigbour(g.target(*a));
893  SmoothFactorType smoothFactor= weightsToSmoothFactor(edgeWeights[edge]);
894 
895  NodeFeatureInValue neighbourFeat = nodeFeaturesIn[neigbour];
896  neighbourFeat*=smoothFactor;
897  if(degree==0)
898  featOut = neighbourFeat;
899  else
900  featOut += neighbourFeat;
901  weightSum+=smoothFactor;
902  ++degree;
903  }
904  // fixme..set me to right type
905  featIn*=static_cast<float>(degree);
906  weightSum+=static_cast<float>(degree);
907  featOut+=featIn;
908  featOut/=weightSum;
909  }
910  }
911 
912  template<class T>
913  struct ExpSmoothFactor{
914  ExpSmoothFactor(const T lambda,const T edgeThreshold,const T scale)
915  : lambda_(lambda),
916  edgeThreshold_(edgeThreshold),
917  scale_(scale){
918  }
919  T operator()(const T weight){
920  return weight> edgeThreshold_ ? 0 : std::exp(-1.0*lambda_*weight)*scale_;
921  }
922  T lambda_;
923  T edgeThreshold_;
924  T scale_;
925  };
926 
927  } // namespace detail_graph_smoothing
928 
929 
930  /// \brief smooth node features of a graph
931  ///
932  /// \param g : input graph
933  /// \param nodeFeaturesIn : input node features which should be smoothed
934  /// \param edgeIndicator : edge indicator to indicate over which edges one should smooth
935  /// \param lambda : scale edge indicator by lambda bevore taking negative exponent
936  /// \param edgeThreshold : edge threshold
937  /// \param scale : how much smoothing should be applied
938  /// \param[out] nodeFeaturesOut : smoothed node features
939  template<class GRAPH, class NODE_FEATURES_IN,class EDGE_INDICATOR,class NODE_FEATURES_OUT>
941  const GRAPH & g,
942  const NODE_FEATURES_IN & nodeFeaturesIn,
943  const EDGE_INDICATOR & edgeIndicator,
944  const float lambda,
945  const float edgeThreshold,
946  const float scale,
947  NODE_FEATURES_OUT & nodeFeaturesOut
948  ){
949  detail_graph_smoothing::ExpSmoothFactor<float> functor(lambda,edgeThreshold,scale);
950  detail_graph_smoothing::graphSmoothingImpl(g,nodeFeaturesIn,edgeIndicator,functor,nodeFeaturesOut);
951  }
952 
953  /// \brief smooth node features of a graph
954  ///
955  /// \param g : input graph
956  /// \param nodeFeaturesIn : input node features which should be smoothed
957  /// \param edgeIndicator : edge indicator to indicate over which edges one should smooth
958  /// \param lambda : scale edge indicator by lambda bevore taking negative exponent
959  /// \param edgeThreshold : edge threshold
960  /// \param scale : how much smoothing should be applied
961  /// \param iterations : how often should this algorithm be called recursively
962  /// \param[out] nodeFeaturesBuffer : preallocated(!) buffer to store node features temp.
963  /// \param[out] nodeFeaturesOut : smoothed node features
964  template<class GRAPH, class NODE_FEATURES_IN,class EDGE_INDICATOR,class NODE_FEATURES_OUT>
966  const GRAPH & g,
967  const NODE_FEATURES_IN & nodeFeaturesIn,
968  const EDGE_INDICATOR & edgeIndicator,
969  const float lambda,
970  const float edgeThreshold,
971  const float scale,
972  size_t iterations,
973  NODE_FEATURES_OUT & nodeFeaturesBuffer,
974  NODE_FEATURES_OUT & nodeFeaturesOut
975  ){
976 
977  iterations = std::max(size_t(1),iterations);
978  // initial run
979  graphSmoothing(g,nodeFeaturesIn,edgeIndicator,lambda,edgeThreshold,scale,nodeFeaturesOut);
980  iterations -=1;
981 
982  bool outAsIn=true;
983  for(size_t i=0;i<iterations;++i){
984  if(outAsIn){
985  graphSmoothing(g,nodeFeaturesOut,edgeIndicator,lambda,edgeThreshold,scale,nodeFeaturesBuffer);
986  outAsIn=false;
987  }
988  else{
989  graphSmoothing(g,nodeFeaturesBuffer,edgeIndicator,lambda,edgeThreshold,scale,nodeFeaturesOut);
990  outAsIn=true;
991  }
992  }
993  if(!outAsIn){
994  copyNodeMap(g,nodeFeaturesBuffer,nodeFeaturesOut);
995  }
996  }
997 
998  /// project ground truth from a base graph, to a region adjacency graph.
999  ///
1000  ///
1001  ///
1002  ///
1003  template<class RAG, class BASE_GRAPH, class BASE_GRAPH_RAG_LABELS,
1004  class BASE_GRAPH_GT, class RAG_GT, class RAG_GT_QT>
1006  const RAG & rag,
1007  const BASE_GRAPH & baseGraph,
1008  const BASE_GRAPH_RAG_LABELS & baseGraphRagLabels,
1009  const BASE_GRAPH_GT & baseGraphGt,
1010  RAG_GT & ragGt,
1011  RAG_GT_QT & ragGtQt
1012  ){
1013  typedef typename BASE_GRAPH::Node BaseGraphNode;
1014  typedef typename BASE_GRAPH::NodeIt BaseGraphNodeIt;
1015  typedef typename RAG::NodeIt RagNodeIt;
1016  typedef typename RAG::Node RagNode;
1017  typedef typename BASE_GRAPH_RAG_LABELS::Value BaseRagLabelType;
1018  typedef typename BASE_GRAPH_GT::Value GtLabelType;
1019  typedef typename RAG_GT::Value RagGtLabelType;
1020 
1021  // overlap map
1022  typedef std::map<GtLabelType,UInt32> MapType;
1023  typedef typename MapType::const_iterator MapIter;
1024  typedef typename RAG:: template NodeMap<MapType> Overlap;
1025  Overlap overlap(rag);
1026 
1027  for(BaseGraphNodeIt baseNodeIter(baseGraph); baseNodeIter!=lemon::INVALID; ++baseNodeIter ){
1028 
1029  const BaseGraphNode baseNode = *baseNodeIter;
1030 
1031  // get gt label
1032  const GtLabelType gtLabel = baseGraphGt[baseNode];
1033 
1034  // get the label which generated rag
1035  // (node mapping from bg-node to rag-node-id)
1036  const BaseRagLabelType bgRagLabel = baseGraphRagLabels[baseNode];
1037  const RagNode ragNode = rag.nodeFromId(bgRagLabel);
1038 
1039  // fill overlap
1040  overlap[ragNode][gtLabel]+=1;
1041  }
1042 
1043  // select label with max overlap
1044  for(RagNodeIt ragNodeIter(rag); ragNodeIter!=lemon::INVALID; ++ragNodeIter ){
1045  const RagNode ragNode = *ragNodeIter;
1046  const MapType olMap = overlap[ragNode];
1047  UInt32 olSize=0;
1048  RagGtLabelType bestLabel = 0;
1049  for(MapIter olIter = olMap.begin(); olIter!=olMap.end(); ++olIter ){
1050  if(olIter->second > olSize){
1051  olSize = olIter->second;
1052  bestLabel = olIter->first;
1053  }
1054  }
1055  ragGt[ragNode]=bestLabel;
1056  }
1057  }
1058 
1059  /// \brief create edge weights from node weights
1060  ///
1061  /// \param g : input graph
1062  /// \param nodeWeights : node property map holding node weights
1063  /// \param[out] edgeWeights : resulting edge weights
1064  /// \param euclidean : if 'true', multiply the computed weights with the Euclidean
1065  /// distance between the edge's end nodes (default: 'false')
1066  /// \param func : binary function that computes the edge weight from the
1067  /// weights of the edge's end nodes (default: take the average)
1068  template<unsigned int N, class DirectedTag,
1069  class NODEMAP, class EDGEMAP, class FUNCTOR>
1070  void
1072  const GridGraph<N, DirectedTag> & g,
1073  const NODEMAP & nodeWeights,
1074  EDGEMAP & edgeWeights,
1075  bool euclidean,
1076  FUNCTOR const & func)
1077  {
1078  typedef GridGraph<N, DirectedTag> Graph;
1079  typedef typename Graph::Edge Edge;
1080  typedef typename Graph::EdgeIt EdgeIt;
1081  typedef typename MultiArrayShape<N>::type CoordType;
1082 
1083  vigra_precondition(nodeWeights.shape() == g.shape(),
1084  "edgeWeightsFromNodeWeights(): shape mismatch between graph and nodeWeights.");
1085 
1086  for (EdgeIt iter(g); iter!=lemon::INVALID; ++iter)
1087  {
1088  const Edge edge(*iter);
1089  const CoordType uCoord(g.u(edge));
1090  const CoordType vCoord(g.v(edge));
1091  if (euclidean)
1092  {
1093  edgeWeights[edge] = norm(uCoord-vCoord) * func(nodeWeights[uCoord], nodeWeights[vCoord]);
1094  }
1095  else
1096  {
1097  edgeWeights[edge] = func(nodeWeights[uCoord], nodeWeights[vCoord]);
1098  }
1099  }
1100  }
1101 
1102  template<unsigned int N, class DirectedTag,
1103  class NODEMAP, class EDGEMAP>
1104  inline void
1106  const GridGraph<N, DirectedTag> & g,
1107  const NODEMAP & nodeWeights,
1108  EDGEMAP & edgeWeights,
1109  bool euclidean=false)
1110  {
1111  using namespace vigra::functor;
1112  edgeWeightsFromNodeWeights(g, nodeWeights, edgeWeights, euclidean, Param(0.5)*(Arg1()+Arg2()));
1113  }
1114 
1115 
1116  /// \brief create edge weights from an interpolated image
1117  ///
1118  /// \param g : input graph
1119  /// \param interpolatedImage : interpolated image
1120  /// \param[out] edgeWeights : edge weights
1121  /// \param euclidean : if 'true', multiply the weights with the Euclidean
1122  /// distance between the edge's end nodes (default: 'false')
1123  ///
1124  /// For each edge, the function reads the weight from <tt>interpolatedImage[u+v]</tt>,
1125  /// where <tt>u</tt> and <tt>v</tt> are the coordinates of the edge's end points.
1126  template<unsigned int N, class DirectedTag,
1127  class T, class EDGEMAP>
1128  void
1130  const GridGraph<N, DirectedTag> & g,
1131  const MultiArrayView<N, T> & interpolatedImage,
1132  EDGEMAP & edgeWeights,
1133  bool euclidean = false)
1134  {
1135  typedef GridGraph<N, DirectedTag> Graph;
1136  typedef typename Graph::Edge Edge;
1137  typedef typename Graph::EdgeIt EdgeIt;
1138  typedef typename MultiArrayShape<N>::type CoordType;
1139 
1140  vigra_precondition(interpolatedImage.shape() == 2*g.shape()-CoordType(1),
1141  "edgeWeightsFromInterpolatedImage(): interpolated shape must be shape*2-1");
1142 
1143  for (EdgeIt iter(g); iter!=lemon::INVALID; ++iter)
1144  {
1145  const Edge edge(*iter);
1146  const CoordType uCoord(g.u(edge));
1147  const CoordType vCoord(g.v(edge));
1148  if (euclidean)
1149  {
1150  edgeWeights[edge] = norm(uCoord-vCoord) * interpolatedImage[uCoord+vCoord];
1151  }
1152  else
1153  {
1154  edgeWeights[edge] = interpolatedImage[uCoord+vCoord];
1155  }
1156  }
1157  }
1158 
1159 //@}
1160 
1161 } // namespace vigra
1162 
1163 #endif // VIGRA_GRAPH_ALGORITHMS_HXX

© 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.10.0 (Thu Jan 8 2015)