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

hierarchical_clustering.hxx VIGRA

1 /************************************************************************/
2 /* */
3 /* Copyright 2011 by 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 
38 #ifndef VIGRA_HIERARCHICAL_CLUSTERING_HXX
39 #define VIGRA_HIERARCHICAL_CLUSTERING_HXX
40 
41 
42 
43 /*std*/
44 #include <queue>
45 #include <iomanip>
46 
47 /*vigra*/
48 #include "priority_queue.hxx"
49 #include "metrics.hxx"
50 #include "merge_graph_adaptor.hxx"
51 
52 namespace vigra{
53 
54 /** \addtogroup GraphDataStructures
55 */
56 //@{
57 
58 namespace cluster_operators{
59 
60 template<
61  class MERGE_GRAPH,
62  class EDGE_INDICATOR_MAP,
63  class EDGE_SIZE_MAP,
64  class NODE_SIZE_MAP,
65  class MIN_WEIGHT_MAP
66  >
67 class EdgeWeightedUcm
68 {
69  typedef EdgeWeightedUcm<
70  MERGE_GRAPH,
71  EDGE_INDICATOR_MAP,
72  EDGE_SIZE_MAP,
73  NODE_SIZE_MAP,
74  MIN_WEIGHT_MAP
75  > SelfType;
76 
77  public:
78 
79  typedef typename EDGE_INDICATOR_MAP::Value ValueType;
80  typedef ValueType WeightType;
81  typedef MERGE_GRAPH MergeGraph;
82  typedef typename MergeGraph::Graph Graph;
83  typedef typename Graph::Edge BaseGraphEdge;
84  typedef typename Graph::Node BaseGraphNode;
85  typedef typename MergeGraph::Edge Edge;
86  typedef typename MergeGraph::Node Node;
87  typedef typename MergeGraph::EdgeIt EdgeIt;
88  typedef typename MergeGraph::NodeIt NodeIt;
89  typedef typename MergeGraph::IncEdgeIt IncEdgeIt;
90  typedef typename MergeGraph::index_type index_type;
91  typedef MergeGraphItemHelper<MergeGraph,Edge> EdgeHelper;
92  typedef MergeGraphItemHelper<MergeGraph,Node> NodeHelper;
93 
94  typedef typename MergeGraph::MergeNodeCallBackType MergeNodeCallBackType;
95  typedef typename MergeGraph::MergeEdgeCallBackType MergeEdgeCallBackType;
96  typedef typename MergeGraph::EraseEdgeCallBackType EraseEdgeCallBackType;
97 
98  typedef typename EDGE_INDICATOR_MAP::Reference EdgeIndicatorReference;
99  /// \brief construct cluster operator
100  EdgeWeightedUcm(
101  MergeGraph & mergeGraph,
102  EDGE_INDICATOR_MAP edgeIndicatorMap,
103  EDGE_SIZE_MAP edgeSizeMap,
104  NODE_SIZE_MAP nodeSizeMap,
105  MIN_WEIGHT_MAP minWeightEdgeMap,
106  const ValueType wardness=1.0
107  )
108  : mergeGraph_(mergeGraph),
109  edgeIndicatorMap_(edgeIndicatorMap),
110  edgeSizeMap_(edgeSizeMap),
111  nodeSizeMap_(nodeSizeMap),
112  minWeightEdgeMap_(minWeightEdgeMap),
113  pq_(mergeGraph.maxEdgeId()+1),
114  wardness_(wardness)
115  {
116  MergeNodeCallBackType cbMn(MergeNodeCallBackType:: template from_method<SelfType,&SelfType::mergeNodes>(this));
117  MergeEdgeCallBackType cbMe(MergeEdgeCallBackType:: template from_method<SelfType,&SelfType::mergeEdges>(this));
118  EraseEdgeCallBackType cbEe(EraseEdgeCallBackType:: template from_method<SelfType,&SelfType::eraseEdge>(this));
119 
120  mergeGraph_.registerMergeNodeCallBack(cbMn);
121  mergeGraph_.registerMergeEdgeCallBack(cbMe);
122  mergeGraph_.registerEraseEdgeCallBack(cbEe);
123 
124  for(EdgeIt e(mergeGraph_);e!=lemon::INVALID;++e){
125  const Edge edge = *e;
126  const BaseGraphEdge graphEdge=EdgeHelper::itemToGraphItem(mergeGraph_,edge);
127  const index_type edgeId = mergeGraph_.id(edge);
128  const ValueType currentWeight = this->getEdgeWeight(edge);
129  pq_.push(edgeId,currentWeight);
130  minWeightEdgeMap_[graphEdge]=currentWeight;
131  }
132  }
133 
134  void setWardness(const float w){
135  wardness_ = w;
136  }
137 
138  void resetMgAndPq(){
139  mergeGraph_.reset();
140 
141  MergeNodeCallBackType cbMn(MergeNodeCallBackType:: template from_method<SelfType,&SelfType::mergeNodes>(this));
142  MergeEdgeCallBackType cbMe(MergeEdgeCallBackType:: template from_method<SelfType,&SelfType::mergeEdges>(this));
143  EraseEdgeCallBackType cbEe(EraseEdgeCallBackType:: template from_method<SelfType,&SelfType::eraseEdge>(this));
144 
145  mergeGraph_.registerMergeNodeCallBack(cbMn);
146  mergeGraph_.registerMergeEdgeCallBack(cbMe);
147  mergeGraph_.registerEraseEdgeCallBack(cbEe);
148 
149  pq_.reset();
150  for(EdgeIt e(mergeGraph_);e!=lemon::INVALID;++e){
151  const Edge edge = *e;
152  const BaseGraphEdge graphEdge=EdgeHelper::itemToGraphItem(mergeGraph_,edge);
153  const index_type edgeId = mergeGraph_.id(edge);
154  const ValueType currentWeight = this->getEdgeWeight(edge);
155  pq_.push(edgeId,currentWeight);
156  minWeightEdgeMap_[graphEdge]=currentWeight;
157  }
158  }
159 
160  /// \brief will be called via callbacks from mergegraph
161  void mergeEdges(const Edge & a,const Edge & b){
162  // update features / weigts etc
163  const BaseGraphEdge aa=EdgeHelper::itemToGraphItem(mergeGraph_,a);
164  const BaseGraphEdge bb=EdgeHelper::itemToGraphItem(mergeGraph_,b);
165  EdgeIndicatorReference va=edgeIndicatorMap_[aa];
166  EdgeIndicatorReference vb=edgeIndicatorMap_[bb];
167  va*=edgeSizeMap_[aa];
168  vb*=edgeSizeMap_[bb];
169 
170  va+=vb;
171  edgeSizeMap_[aa]+=edgeSizeMap_[bb];
172  va/=(edgeSizeMap_[aa]);
173  vb/=edgeSizeMap_[bb];
174  // delete b from pq
175  pq_.deleteItem(b.id());
176  }
177 
178  /// \brief will be called via callbacks from mergegraph
179  void mergeNodes(const Node & a,const Node & b){
180  const BaseGraphNode aa=NodeHelper::itemToGraphItem(mergeGraph_,a);
181  const BaseGraphNode bb=NodeHelper::itemToGraphItem(mergeGraph_,b);
182  nodeSizeMap_[aa]+=nodeSizeMap_[bb];
183  }
184 
185  /// \brief will be called via callbacks from mergegraph
186  void eraseEdge(const Edge & edge){
187 
188  //std::cout<<"start to erase edge "<<mergeGraph_.id(edge)<<"\n";
189  // delete edge from pq
190  pq_.deleteItem(edge.id());
191  // get the new region the edge is in
192  // (since the edge is no any more an active edge)
193  //std::cout<<"get the new node \n";
194  const Node newNode = mergeGraph_.inactiveEdgesNode(edge);
195  //std::cout<<"new node "<<mergeGraph_.id(newNode)<<"\n";
196 
197  // iterate over all edges of this node
198  for (IncEdgeIt e(mergeGraph_,newNode);e!=lemon::INVALID;++e)
199  {
200  //std::cout<<"get inc edge\n";
201  const Edge incEdge(*e);
202 
203  //std::cout<<"get inc graph edge\n";
204  const BaseGraphEdge incGraphEdge = EdgeHelper::itemToGraphItem(mergeGraph_,incEdge);
205 
206  //std::cout<<"get inc edge weight"<<counter<<"\n";
207  // compute the new weight for this edge
208  // (this should involve region differences)
209  const ValueType newWeight = getEdgeWeight(incEdge);
210 
211  // change the weight in pq by repushing
212  pq_.push(incEdge.id(),newWeight);
213  minWeightEdgeMap_[incGraphEdge]=newWeight;
214 
215  }
216  //std::cout<<"done\n";
217  }
218 
219  /// \brief get the edge which should be contracted next
220  Edge contractionEdge(){
221  index_type minLabel = pq_.top();
222  while(mergeGraph_.hasEdgeId(minLabel)==false){
223  eraseEdge(Edge(minLabel));
224  minLabel = pq_.top();
225  }
226  return Edge(minLabel);
227  }
228 
229  /// \brief get the edge weight of the edge which should be contracted next
230  WeightType contractionWeight(){
231  index_type minLabel = pq_.top();
232  while(mergeGraph_.hasEdgeId(minLabel)==false){
233  eraseEdge(Edge(minLabel));
234  minLabel = pq_.top();
235  }
236  return pq_.topPriority();
237 
238  }
239 
240 
241  /// \brief get a reference to the merge
242  MergeGraph & mergeGraph(){
243  return mergeGraph_;
244  }
245 
246  bool done(){
247 
248  index_type minLabel = pq_.top();
249  while(mergeGraph_.hasEdgeId(minLabel)==false){
250  eraseEdge(Edge(minLabel));
251  minLabel = pq_.top();
252  }
253  return mergeGraph_.edgeNum()==0 || mergeGraph_.nodeNum()==1;
254  }
255 
256  private:
257  ValueType getEdgeWeight(const Edge & e){
258 
259  const Node u = mergeGraph_.u(e);
260  const Node v = mergeGraph_.v(e);
261 
262  const BaseGraphEdge ee=EdgeHelper::itemToGraphItem(mergeGraph_,e);
263  const BaseGraphNode uu=NodeHelper::itemToGraphItem(mergeGraph_,u);
264  const BaseGraphNode vv=NodeHelper::itemToGraphItem(mergeGraph_,v);
265 
266  const float sizeU = nodeSizeMap_[uu] ;
267  const float sizeV = nodeSizeMap_[vv] ;
268 
269  const ValueType wardFac = 2.0 / ( 1.0/std::pow(sizeU,wardness_) + 1/std::pow(sizeV,wardness_) );
270  //const ValueType wardFac = (wardFacRaw*wardness_) + (1.0-wardness_);
271 
272  const ValueType fromEdgeIndicator = edgeIndicatorMap_[ee];
273  const ValueType totalWeight = fromEdgeIndicator*wardFac;
274  return totalWeight;
275  }
276 
277 
278  MergeGraph & mergeGraph_;
279  EDGE_INDICATOR_MAP edgeIndicatorMap_;
280  EDGE_SIZE_MAP edgeSizeMap_;
281  NODE_SIZE_MAP nodeSizeMap_;
282  MIN_WEIGHT_MAP minWeightEdgeMap_;
284  ValueType wardness_;;
285 };
286 
287  /// \brief This Cluster Operator is a MONSTER.
288  /// It can really do a lot.
289  ///
290  /// Each edge has a single scalar weight w_e.
291  /// Each node has a feature vector f_n.
292  /// (all f_n's have the same length).
293  /// Edges and nodes have a length / size
294  ///
295  /// The total edge weight is computed via a complicated formula
296  ///
297  /// The main idea is the following.
298  /// Use a mixture between the edge weights w_e,
299  /// and node based edge weights which are computed
300  /// via a metric which measures the 'difference' between
301  /// the u/v feature vectors f_n.
302  ///
303  /// Furthermore a 'Ward'-like regularization can be applied.
304  /// This is useful if one have clusters with sizes
305  /// in the same magnitude (or 'similar' sizes).
306  /// The amount of 'ward'-regularization is controlled
307  /// with the 'wardness' parameter.
308  ///
309  /// Also labels (in the sense of seeds) can be attached to get a 'watershed-ish'
310  /// behavior (nodes with different labels will never be merged)
311  /// The '0'-Label is used to indicate that there is no label at all.
312  /// If certain connected regions share the same seed/label
313  /// it is not guaranteed that they will merge. But a certain prior / multiplier
314  /// must be specified. The total weight of an edge where the u/v node have
315  /// the same label is multiplied with this very multiplier.
316  template<
317  class MERGE_GRAPH,
318  class EDGE_INDICATOR_MAP,
319  class EDGE_SIZE_MAP,
320  class NODE_FEATURE_MAP,
321  class NODE_SIZE_MAP,
322  class MIN_WEIGHT_MAP,
323  class NODE_LABEL_MAP
324  >
326 
327  typedef EdgeWeightNodeFeatures<
328  MERGE_GRAPH,
329  EDGE_INDICATOR_MAP,
330  EDGE_SIZE_MAP,
331  NODE_FEATURE_MAP,
332  NODE_SIZE_MAP,
333  MIN_WEIGHT_MAP,
334  NODE_LABEL_MAP
335  > SelfType;
336  public:
337 
338 
339  typedef typename EDGE_INDICATOR_MAP::Value ValueType;
340  typedef ValueType WeightType;
341  typedef MERGE_GRAPH MergeGraph;
342  typedef typename MergeGraph::Graph Graph;
343  typedef typename Graph::Edge BaseGraphEdge;
344  typedef typename Graph::Node BaseGraphNode;
345  typedef typename MergeGraph::Edge Edge;
346  typedef typename MergeGraph::Node Node;
347  typedef typename MergeGraph::EdgeIt EdgeIt;
348  typedef typename MergeGraph::NodeIt NodeIt;
349  typedef typename MergeGraph::IncEdgeIt IncEdgeIt;
350  typedef typename MergeGraph::index_type index_type;
351  typedef MergeGraphItemHelper<MergeGraph,Edge> EdgeHelper;
352  typedef MergeGraphItemHelper<MergeGraph,Node> NodeHelper;
353 
354 
355  typedef typename EDGE_INDICATOR_MAP::Reference EdgeIndicatorReference;
356  typedef typename NODE_FEATURE_MAP::Reference NodeFeatureReference;
357  /// \brief construct cluster operator
359  MergeGraph & mergeGraph,
360  EDGE_INDICATOR_MAP edgeIndicatorMap,
361  EDGE_SIZE_MAP edgeSizeMap,
362  NODE_FEATURE_MAP nodeFeatureMap,
363  NODE_SIZE_MAP nodeSizeMap,
364  MIN_WEIGHT_MAP minWeightEdgeMap,
365  NODE_LABEL_MAP nodeLabelMap,
366  const ValueType beta,
367  const metrics::MetricType metricType,
368  const ValueType wardness=static_cast<ValueType>(1.0),
369  const ValueType gamma = static_cast<ValueType>(10000000.0),
370  const ValueType sameLabelMultiplier = static_cast<ValueType>(0.8)
371  )
372  : mergeGraph_(mergeGraph),
373  edgeIndicatorMap_(edgeIndicatorMap),
374  edgeSizeMap_(edgeSizeMap),
375  nodeFeatureMap_(nodeFeatureMap),
376  nodeSizeMap_(nodeSizeMap),
377  minWeightEdgeMap_(minWeightEdgeMap),
378  nodeLabelMap_(nodeLabelMap),
379  pq_(mergeGraph.maxEdgeId()+1),
380  beta_(beta),
381  wardness_(wardness),
382  gamma_(gamma),
383  sameLabelMultiplier_(sameLabelMultiplier),
384  metric_(metricType)
385  {
386  typedef typename MergeGraph::MergeNodeCallBackType MergeNodeCallBackType;
387  typedef typename MergeGraph::MergeEdgeCallBackType MergeEdgeCallBackType;
388  typedef typename MergeGraph::EraseEdgeCallBackType EraseEdgeCallBackType;
389 
390 
391  MergeNodeCallBackType cbMn(MergeNodeCallBackType:: template from_method<SelfType,&SelfType::mergeNodes>(this));
392  MergeEdgeCallBackType cbMe(MergeEdgeCallBackType:: template from_method<SelfType,&SelfType::mergeEdges>(this));
393  EraseEdgeCallBackType cbEe(EraseEdgeCallBackType:: template from_method<SelfType,&SelfType::eraseEdge>(this));
394 
395  mergeGraph_.registerMergeNodeCallBack(cbMn);
396  mergeGraph_.registerMergeEdgeCallBack(cbMe);
397  mergeGraph_.registerEraseEdgeCallBack(cbEe);
398 
399 
400 
401  for(EdgeIt e(mergeGraph);e!=lemon::INVALID;++e){
402  const Edge edge = *e;
403  const BaseGraphEdge graphEdge=EdgeHelper::itemToGraphItem(mergeGraph_,edge);
404  const index_type edgeId = mergeGraph_.id(edge);
405  const ValueType currentWeight = this->getEdgeWeight(edge);
406  pq_.push(edgeId,currentWeight);
407  minWeightEdgeMap_[graphEdge]=currentWeight;
408  }
409 
410  }
411 
412  /// \brief will be called via callbacks from mergegraph
413  void mergeEdges(const Edge & a,const Edge & b){
414  // update features / weigts etc
415  const BaseGraphEdge aa=EdgeHelper::itemToGraphItem(mergeGraph_,a);
416  const BaseGraphEdge bb=EdgeHelper::itemToGraphItem(mergeGraph_,b);
417  EdgeIndicatorReference va=edgeIndicatorMap_[aa];
418  EdgeIndicatorReference vb=edgeIndicatorMap_[bb];
419  va*=edgeSizeMap_[aa];
420  vb*=edgeSizeMap_[bb];
421 
422 
423  va+=vb;
424  edgeSizeMap_[aa]+=edgeSizeMap_[bb];
425  va/=(edgeSizeMap_[aa]);
426  vb/=edgeSizeMap_[bb];
427  // delete b from pq
428  pq_.deleteItem(b.id());
429  }
430 
431  /// \brief will be called via callbacks from mergegraph
432  void mergeNodes(const Node & a,const Node & b){
433  const BaseGraphNode aa=NodeHelper::itemToGraphItem(mergeGraph_,a);
434  const BaseGraphNode bb=NodeHelper::itemToGraphItem(mergeGraph_,b);
435  NodeFeatureReference va=nodeFeatureMap_[aa];
436  NodeFeatureReference vb=nodeFeatureMap_[bb];
437  va*=nodeSizeMap_[aa];
438  vb*=nodeSizeMap_[bb];
439  va+=vb;
440  nodeSizeMap_[aa]+=nodeSizeMap_[bb];
441  va/=(nodeSizeMap_[aa]);
442  vb/=nodeSizeMap_[bb];
443 
444 
445  // update labels
446  const UInt32 labelA = nodeLabelMap_[aa];
447  const UInt32 labelB = nodeLabelMap_[bb];
448 
449  if(labelA!=0 && labelB!=0 && labelA!=labelB){
450  throw std::runtime_error("both nodes have labels");
451  }
452  else{
453  const UInt32 newLabel = std::max(labelA, labelB);
454  nodeLabelMap_[aa] = newLabel;
455  }
456  }
457 
458  /// \brief will be called via callbacks from mergegraph
459  void eraseEdge(const Edge & edge){
460 
461  //std::cout<<"start to erase edge "<<mergeGraph_.id(edge)<<"\n";
462  // delete edge from pq
463  pq_.deleteItem(edge.id());
464  // get the new region the edge is in
465  // (since the edge is no any more an active edge)
466  //std::cout<<"get the new node \n";
467  const Node newNode = mergeGraph_.inactiveEdgesNode(edge);
468  //std::cout<<"new node "<<mergeGraph_.id(newNode)<<"\n";
469 
470 
471  // iterate over all edges of this node
472  for (IncEdgeIt e(mergeGraph_,newNode);e!=lemon::INVALID;++e){
473 
474  //std::cout<<"get inc edge\n";
475  const Edge incEdge(*e);
476 
477  //std::cout<<"get inc graph edge\n";
478  const BaseGraphEdge incGraphEdge = EdgeHelper::itemToGraphItem(mergeGraph_,incEdge);
479 
480  //std::cout<<"get inc edge weight"<<counter<<"\n";
481  // compute the new weight for this edge
482  // (this should involve region differences)
483  const ValueType newWeight = getEdgeWeight(incEdge);
484 
485 
486  // change the weight in pq by repushing
487 
488  //std::cout<<"push\n";
489  pq_.push(incEdge.id(),newWeight);
490  minWeightEdgeMap_[incGraphEdge]=newWeight;
491 
492  }
493  //std::cout<<"done\n";
494  }
495 
496  /// \brief get the edge which should be contracted next
498  index_type minLabel = pq_.top();
499  while(mergeGraph_.hasEdgeId(minLabel)==false){
500  pq_.deleteItem(minLabel);
501  minLabel = pq_.top();
502  }
503  return Edge(minLabel);
504  }
505 
506  /// \brief get the edge weight of the edge which should be contracted next
507  WeightType contractionWeight(){
508  index_type minLabel = pq_.top();
509  while(mergeGraph_.hasEdgeId(minLabel)==false){
510  pq_.deleteItem(minLabel);
511  minLabel = pq_.top();
512  }
513  return pq_.topPriority();
514 
515  }
516 
517 
518  /// \brief get a reference to the merge
519  MergeGraph & mergeGraph(){
520  return mergeGraph_;
521  }
522 
523  bool done(){
524 
525  index_type minLabel = pq_.top();
526  while(mergeGraph_.hasEdgeId(minLabel)==false){
527  pq_.deleteItem(minLabel);
528  minLabel = pq_.top();
529  }
530  const ValueType p = pq_.topPriority();
531 
532  return p>= gamma_;
533  }
534 
535  private:
536  ValueType getEdgeWeight(const Edge & e){
537 
538  const Node u = mergeGraph_.u(e);
539  const Node v = mergeGraph_.v(e);
540 
541  const BaseGraphEdge ee=EdgeHelper::itemToGraphItem(mergeGraph_,e);
542  const BaseGraphNode uu=NodeHelper::itemToGraphItem(mergeGraph_,u);
543  const BaseGraphNode vv=NodeHelper::itemToGraphItem(mergeGraph_,v);
544 
545  const float sizeU = nodeSizeMap_[uu];
546  const float sizeV = nodeSizeMap_[vv];
547 
548 
549  const ValueType wardFac = 2.0 / ( 1.0/std::pow(sizeU,wardness_) + 1/std::pow(sizeV,wardness_) );
550 
551  const ValueType fromEdgeIndicator = edgeIndicatorMap_[ee];
552  ValueType fromNodeDist = metric_(nodeFeatureMap_[uu],nodeFeatureMap_[vv]);
553  ValueType totalWeight = ((1.0-beta_)*fromEdgeIndicator + beta_*fromNodeDist)*wardFac;
554 
555 
556  const UInt32 labelA = nodeLabelMap_[uu];
557  const UInt32 labelB = nodeLabelMap_[vv];
558 
559  if(labelA!=0 && labelB!=0){
560  if(labelA == labelB){
561  totalWeight*=sameLabelMultiplier_;
562  }
563  else{
564  totalWeight += gamma_;
565  }
566  }
567  return totalWeight;
568  }
569 
570 
571  MergeGraph & mergeGraph_;
572  EDGE_INDICATOR_MAP edgeIndicatorMap_;
573  EDGE_SIZE_MAP edgeSizeMap_;
574  NODE_FEATURE_MAP nodeFeatureMap_;
575  NODE_SIZE_MAP nodeSizeMap_;
576  MIN_WEIGHT_MAP minWeightEdgeMap_;
577  NODE_LABEL_MAP nodeLabelMap_;
579  ValueType beta_;
580  ValueType wardness_;
581  ValueType gamma_;
582  ValueType sameLabelMultiplier_;
583  metrics::Metric<float> metric_;
584  };
585 
586 
587 
588 } // end namespace cluster_operators
589 
590 /** \brief Options object for hierarchical clustering.
591 
592  <b>\#include</b> <vigra/hierarchical_clustering.hxx><br/>
593  Namespace: vigra
594 
595  This class allows to set various parameters of \ref hierarchicalClustering().
596  See there for usage examples.
597 */
599 {
600  public:
601 
603  const size_t nodeNumStopCond = 1,
604  const bool buildMergeTree = false,
605  const bool verbose = false)
606  : nodeNumStopCond_ (nodeNumStopCond)
607  , maxMergeWeight_(NumericTraits<double>::max())
608  , nodeFeatureImportance_(0.5)
609  , sizeImportance_(1.0)
610  , nodeFeatureMetric_(metrics::ManhattanMetric)
611  , buildMergeTreeEncoding_(buildMergeTree)
612  , verbose_(verbose)
613  {}
614 
615  /** Stop merging when the number of clusters reaches this threshold.
616 
617  Default: 1 (don't stop early)
618  */
620  {
621  nodeNumStopCond_ = count;
622  return *this;
623  }
624 
625  /** Stop merging when the weight of the cheapest edge exceeds this threshold.
626 
627  Default: infinity (don't stop early)
628  */
630  {
631  maxMergeWeight_ = val;
632  return *this;
633  }
634 
635  /** Importance of node features relative to edge weights.
636 
637  Must be between 0 and 1, with 0 meaning that node features
638  are ignored, and 1 meaning that edge weights are ignored.
639 
640  Default: 0.5 (equal importance)
641  */
643  {
644  vigra_precondition(0.0 <= val && val <= 1.0,
645  "ClusteringOptions::nodePropertyImportance(val): 0 <= val <= 1 required.");
646  nodeFeatureImportance_ = val;
647  return *this;
648  }
649 
650  /** Importance of node size.
651 
652  Must be between 0 and 1, with 0 meaning that size is ignored,
653  and 1 meaning that the algorithm prefers to keep cluster sizes
654  balanced.
655 
656  Default: 1.0 (prefer like-sized clusters)
657  */
659  {
660  vigra_precondition(0.0 <= val && val <= 1.0,
661  "ClusteringOptions::sizeImportance(val): 0 <= val <= 1 required.");
662  sizeImportance_ = val;
663  return *this;
664  }
665 
666  /** Metric to be used when transforming node features into cluster distances.
667 
668  The cluster (= node) distance is the respective norm of the difference
669  vector between the corresponding node feature vectors.
670 
671  Default: metrics::ManhattanMetric (L1-norm of the feature difference)
672  */
674  {
675  nodeFeatureMetric_ = metric;
676  return *this;
677  }
678 
679  ClusteringOptions & buildMergeTreeEncoding(bool val=true)
680  {
681  buildMergeTreeEncoding_ = val;
682  return *this;
683  }
684 
685  /** Display progress information.
686 
687  Default: false
688  */
689  ClusteringOptions & verbose(bool val=true)
690  {
691  verbose_ = val;
692  return *this;
693  }
694 
695  size_t nodeNumStopCond_;
696  double maxMergeWeight_;
697  double nodeFeatureImportance_;
698  double sizeImportance_;
699  metrics::MetricType nodeFeatureMetric_;
700  bool buildMergeTreeEncoding_;
701  bool verbose_;
702 };
703 
704 // \brief do hierarchical clustering with a given cluster operator
705 template< class CLUSTER_OPERATOR>
706 class HierarchicalClusteringImpl
707 {
708  public:
709  typedef CLUSTER_OPERATOR ClusterOperator;
710  typedef typename ClusterOperator::MergeGraph MergeGraph;
711  typedef typename MergeGraph::Graph Graph;
712  typedef typename Graph::Edge BaseGraphEdge;
713  typedef typename Graph::Node BaseGraphNode;
714  typedef typename MergeGraph::Edge Edge;
715  typedef typename MergeGraph::Node Node;
716  typedef typename CLUSTER_OPERATOR::WeightType ValueType;
717  typedef typename MergeGraph::index_type MergeGraphIndexType;
718 
719  typedef ClusteringOptions Parameter;
720 
721  struct MergeItem{
722  MergeItem(
723  const MergeGraphIndexType a,
724  const MergeGraphIndexType b,
725  const MergeGraphIndexType r,
726  const ValueType w
727  ):
728  a_(a),b_(b),r_(r),w_(w){
729  }
730  MergeGraphIndexType a_;
731  MergeGraphIndexType b_;
732  MergeGraphIndexType r_;
733  ValueType w_;
734  };
735 
736  typedef std::vector<MergeItem> MergeTreeEncoding;
737 
738  /// \brief construct HierarchicalClusteringImpl from clusterOperator and an optional parameter object
739  HierarchicalClusteringImpl(
740  ClusterOperator & clusterOperator,
741  const Parameter & parameter = Parameter()
742  )
743  :
744  clusterOperator_(clusterOperator),
745  param_(parameter),
746  mergeGraph_(clusterOperator_.mergeGraph()),
747  graph_(mergeGraph_.graph()),
748  timestamp_(graph_.maxNodeId()+1),
749  toTimeStamp_(),
750  timeStampIndexToMergeIndex_(),
751  mergeTreeEndcoding_()
752  {
753  if(param_.buildMergeTreeEncoding_){
754  // this can be be made smater since user can pass
755  // stoping condition based on nodeNum
756  mergeTreeEndcoding_.reserve(graph_.nodeNum()*2);
757  toTimeStamp_.resize(graph_.maxNodeId()+1);
758  timeStampIndexToMergeIndex_.resize(graph_.maxNodeId()+1);
759  for(MergeGraphIndexType nodeId=0;nodeId<=mergeGraph_.maxNodeId();++nodeId){
760  toTimeStamp_[nodeId]=nodeId;
761  }
762  }
763 
764 
765 
766  }
767 
768  /// \brief start the clustering
769  void cluster(){
770  if(param_.verbose_)
771  std::cout<<"\n";
772  while(mergeGraph_.nodeNum()>param_.nodeNumStopCond_ && mergeGraph_.edgeNum()>0 && !clusterOperator_.done()){
773 
774  const Edge edgeToRemove = clusterOperator_.contractionEdge();
775  if(param_.buildMergeTreeEncoding_){
776  const MergeGraphIndexType uid = mergeGraph_.id(mergeGraph_.u(edgeToRemove));
777  const MergeGraphIndexType vid = mergeGraph_.id(mergeGraph_.v(edgeToRemove));
778  const ValueType w = clusterOperator_.contractionWeight();
779  // do the merge
780  mergeGraph_.contractEdge( edgeToRemove);
781  const MergeGraphIndexType aliveNodeId = mergeGraph_.hasNodeId(uid) ? uid : vid;
782  const MergeGraphIndexType deadNodeId = aliveNodeId==vid ? uid : vid;
783  timeStampIndexToMergeIndex_[timeStampToIndex(timestamp_)]=mergeTreeEndcoding_.size();
784  mergeTreeEndcoding_.push_back(MergeItem( toTimeStamp_[aliveNodeId],toTimeStamp_[deadNodeId],timestamp_,w));
785  toTimeStamp_[aliveNodeId]=timestamp_;
786  timestamp_+=1;
787  }
788  else{
789  //std::cout<<"constract\n";
790  // do the merge
791  mergeGraph_.contractEdge( edgeToRemove );
792  }
793  if(param_.verbose_ && mergeGraph_.nodeNum()%1==0){
794  std::cout<<"\rNodes: "<<std::setw(10)<<mergeGraph_.nodeNum()<<std::flush;
795  }
796 
797  }
798  if(param_.verbose_)
799  std::cout<<"\n";
800  }
801 
802  /// \brief get the encoding of the merge tree
803  const MergeTreeEncoding & mergeTreeEndcoding()const{
804  return mergeTreeEndcoding_;
805  }
806 
807  template<class EDGE_MAP>
808  void ucmTransform(EDGE_MAP & edgeMap)const{
809  typedef typename Graph::EdgeIt BaseGraphEdgeIt;
810 
811  for(BaseGraphEdgeIt iter(graph()); iter!=lemon::INVALID; ++iter ){
812  const BaseGraphEdge edge=*iter;
813  edgeMap[edge] = edgeMap[mergeGraph().reprGraphEdge(edge)];
814  }
815  }
816 
817  /// \brief get the node id's which are the leafes of a treeNodeId
818  template<class OUT_ITER>
819  size_t leafNodeIds(const MergeGraphIndexType treeNodeId, OUT_ITER begin)const{
820  if(treeNodeId<=graph_.maxNodeId()){
821  *begin=treeNodeId;
822  ++begin;
823  return 1;
824  }
825  else{
826  size_t leafNum=0;
827  std::queue<MergeGraphIndexType> queue;
828  queue.push(treeNodeId);
829 
830  while(!queue.empty()){
831 
832  const MergeGraphIndexType id = queue.front();
833  queue.pop();
834  const MergeGraphIndexType mergeIndex = timeStampToMergeIndex(id);
835  const MergeGraphIndexType ab[]= { mergeTreeEndcoding_[mergeIndex].a_, mergeTreeEndcoding_[mergeIndex].b_};
836 
837  for(size_t i=0;i<2;++i){
838  if(ab[i]<=graph_.maxNodeId()){
839  *begin=ab[i];
840  ++begin;
841  ++leafNum;
842  }
843  else{
844  queue.push(ab[i]);
845  }
846  }
847  }
848  return leafNum;
849  }
850  }
851 
852  /// \brief get the graph the merge graph is based on
853  const Graph & graph()const{
854  return graph_;
855  }
856 
857  /// \brief get the merge graph
858  const MergeGraph & mergeGraph()const{
859  return mergeGraph_;
860  }
861 
862  /// \brief get the representative node id
863  const MergeGraphIndexType reprNodeId(const MergeGraphIndexType id)const{
864  return mergeGraph_.reprNodeId(id);
865  }
866 private:
867 
868  MergeGraphIndexType timeStampToIndex(const MergeGraphIndexType timestamp)const{
869  return timestamp- graph_.maxNodeId();
870  }
871 
872 
873  MergeGraphIndexType timeStampToMergeIndex(const MergeGraphIndexType timestamp)const{
874  return timeStampIndexToMergeIndex_[timeStampToIndex(timestamp)];
875  }
876 
877  ClusterOperator & clusterOperator_;
878  Parameter param_;
879  MergeGraph & mergeGraph_;
880  const Graph & graph_;
881  // parameter object
882 
883 
884  // timestamp
885  MergeGraphIndexType timestamp_;
886  std::vector<MergeGraphIndexType> toTimeStamp_;
887  std::vector<MergeGraphIndexType> timeStampIndexToMergeIndex_;
888  // data which can reconstruct the merge tree
889  MergeTreeEncoding mergeTreeEndcoding_;
890 
891 
892 };
893 
894 /********************************************************/
895 /* */
896 /* hierarchicalClustering */
897 /* */
898 /********************************************************/
899 
900 /** \brief Reduce the number of nodes in a graph by iteratively contracting
901  the cheapest edge.
902 
903  <b> Declarations:</b>
904 
905  \code
906  namespace vigra {
907  template <class GRAPH,
908  class EDGE_WEIGHT_MAP, class EDGE_LENGTH_MAP,
909  class NODE_FEATURE_MAP, class NOSE_SIZE_MAP,
910  class NODE_LABEL_MAP>
911  void
912  hierarchicalClustering(GRAPH const & graph,
913  EDGE_WEIGHT_MAP const & edgeWeights, EDGE_LENGTH_MAP const & edgeLengths,
914  NODE_FEATURE_MAP const & nodeFeatures, NOSE_SIZE_MAP const & nodeSizes,
915  NODE_LABEL_MAP & labelMap,
916  ClusteringOptions options = ClusteringOptions());
917  }
918  \endcode
919 
920  Hierarchical clustering is a simple and versatile image segmentation
921  algorithm that typically operates either directly on the pixels (e.g. on
922  a \ref vigra::GridGraph) or on a region adjacency graph over suitable
923  superpixels (e.g. on an \ref vigra::AdjacencyListGraph). The graph is
924  passed to the function in its first argument. After clustering is completed,
925  the parameter \a labelMap contains a mapping from original node IDs to
926  the ID of the cluster each node belongs to. Cluster IDs correspond to
927  the ID of an arbitrarily chosen representative node within each cluster,
928  i.e. they form a sparse subset of the original IDs.
929 
930  Properties of the graph's edges and nodes are provided in the property maps
931  \a edgeWeights, \a edgeLengths, \a nodeFeatures, and \a nodeSizes. These maps
932  are indexed by edge or node ID and return the corresponding feature. Features
933  must by arithmetic scalars or, in case of node features, scalars or vectors
934  of scalars (precisely: objects that provide <tt>begin()</tt> and <tt>end()</tt>
935  to create an STL range). Edge weights are typically derived from an edge
936  indicator such as the gradient magnitude, and node features are either the
937  responses of a filter family (when clustering on the pixel grid), or region
938  statistics as computed by \ref FeatureAccumulators (when clustering on
939  superpixels).
940 
941  In each step, the algorithm merges the two nodes \f$u\f$ and \f$v\f$ whose
942  cluster distance is smallest, where the cluster distance is defined as
943 
944  \f[
945  d_{uv} = \left( (1-\beta) w_{uv} + \beta || f_u - f_v ||_M \right)
946  \cdot \frac{2}{s_u^{-\omega} + s_v^{-\omega}}
947  \f]
948 
949  with \f$ w_{uv} \f$ denoting the weight of edge \f$uv\f$, \f$f_u\f$ and \f$f_v\f$
950  being the node features (possibly vectors to be compared with metric \f$M\f$),
951  and \f$s_u\f$ and \f$s_v\f$ the corresponding node sizes. The metric is defined
952  in the option object by calling \ref vigra::ClusteringOptions::nodeFeatureMetric()
953  and must be selected from the tags defined in \ref vigra::metrics::MetricType.
954 
955  The parameters \f$0 \le \beta \le 1\f$ and \f$0 \le \omega \le 1\f$ control the
956  relative influence of the inputs: With \f$\beta = 0\f$, the node features are
957  ignored, whereas with \f$\beta = 1\f$ the edge weights are ignored. Similarly,
958  with \f$\omega = 0\f$, the node size is ignored, whereas with \f$\omega = 1\f$,
959  cluster distances are scaled by the harmonic mean of the cluster sizes, making
960  the merging of small clusters more favorable. The parameters are defined in the
961  option object by calling \ref vigra::ClusteringOptions::nodeFeatureImportance() and
962  \ref vigra::ClusteringOptions::sizeImportance() respectively.
963 
964  After each merging step, the features of the resulting cluster \f$z\f$ and the weights
965  of its outgoing edges are updated by mean of the corresponding properties of the original
966  clusters \f$u\f$ and \f$v\f$, weighted by the respective node sizes \f$s_z\f$ and
967  edge lengths \f$l_{zy}\f$:
968 
969  \f{eqnarray*}{
970  s_z & = & s_u + s_v \\
971  f_z & = & \frac{s_u f_u + s_v f_v}{s_z} \\
972  l_{zy} & = & l_{uy} + l_{vy} \textrm{ for all nodes }y\textrm{ connected to }u\textrm{ or }v \\
973  w_{zy} & = & \frac{l_{uy} w_{uy} + l_{vy} w_{vy}}{l_{zy}}
974  \f}
975 
976  Clustering normally stops when only one cluster remains. This default can be overridden
977  by the option object parameters \ref vigra::ClusteringOptions::minRegionCount()
978  and \ref vigra::ClusteringOptions::maxMergeDistance() to stop at a particular number of
979  clusters or a particular cluster distance respectively.
980 
981  <b> Usage:</b>
982 
983  <b>\#include</b> <vigra/hierarchical_clustering.hxx><br>
984  Namespace: vigra
985 
986  A fully worked example can be found in <a href="graph_agglomerative_clustering_8cxx-example.html">graph_agglomerative_clustering.cxx</a>
987 */
989 
990 template <class GRAPH,
991  class EDGE_WEIGHT_MAP, class EDGE_LENGTH_MAP,
992  class NODE_FEATURE_MAP, class NOSE_SIZE_MAP,
993  class NODE_LABEL_MAP>
994 void
995 hierarchicalClustering(GRAPH const & graph,
996  EDGE_WEIGHT_MAP const & edgeWeights, EDGE_LENGTH_MAP const & edgeLengths,
997  NODE_FEATURE_MAP const & nodeFeatures, NOSE_SIZE_MAP const & nodeSizes,
998  NODE_LABEL_MAP & labelMap,
999  ClusteringOptions options = ClusteringOptions())
1000 {
1001  typedef typename NODE_LABEL_MAP::Value LabelType;
1002  typedef MergeGraphAdaptor<GRAPH> MergeGraph;
1003  typedef typename GRAPH::template EdgeMap<float> EdgeUltrametric;
1004  typedef typename GRAPH::template NodeMap<LabelType> NodeSeeds;
1005 
1006  MergeGraph mergeGraph(graph);
1007 
1008  // create property maps to store the computed ultrametric and
1009  // to provide optional cannot-link constraints;
1010  // we don't use these options here and therefore leave the maps empty
1011  EdgeUltrametric edgeUltrametric(graph);
1012  NodeSeeds nodeSeeds(graph);
1013 
1014  // create an operator that stores all property maps needed for
1015  // hierarchical clustering and updates them after every merge step
1016  typedef cluster_operators::EdgeWeightNodeFeatures<
1017  MergeGraph,
1018  EDGE_WEIGHT_MAP,
1019  EDGE_LENGTH_MAP,
1020  NODE_FEATURE_MAP,
1021  NOSE_SIZE_MAP,
1022  EdgeUltrametric,
1023  NodeSeeds>
1024  MergeOperator;
1025 
1026  MergeOperator mergeOperator(mergeGraph,
1027  edgeWeights, edgeLengths,
1028  nodeFeatures, nodeSizes,
1029  edgeUltrametric, nodeSeeds,
1030  options.nodeFeatureImportance_,
1031  options.nodeFeatureMetric_,
1032  options.sizeImportance_,
1033  options.maxMergeWeight_);
1034 
1035  typedef HierarchicalClusteringImpl<MergeOperator> Clustering;
1036 
1037  Clustering clustering(mergeOperator, options);
1038  clustering.cluster();
1039 
1040  for(typename GRAPH::NodeIt node(graph); node != lemon::INVALID; ++node)
1041  {
1042  labelMap[*node] = mergeGraph.reprNodeId(graph.id(*node));
1043  }
1044 }
1045 
1046 //@}
1047 
1048 } // namespace vigra
1049 
1050 #endif // VIGRA_HIERARCHICAL_CLUSTERING_HXX
This Cluster Operator is a MONSTER. It can really do a lot.
Definition: hierarchical_clustering.hxx:325
MergeGraph & mergeGraph()
get a reference to the merge
Definition: hierarchical_clustering.hxx:519
ClusteringOptions & minRegionCount(size_t count)
Definition: hierarchical_clustering.hxx:619
priority_type topPriority() const
get top priority
Definition: priority_queue.hxx:496
ClusteringOptions & sizeImportance(double val)
Definition: hierarchical_clustering.hxx:658
ClusteringOptions & maxMergeDistance(double val)
Definition: hierarchical_clustering.hxx:629
ClusteringOptions & verbose(bool val=true)
Definition: hierarchical_clustering.hxx:689
Manhattan distance (L1 norm)
Definition: metrics.hxx:236
Options object for hierarchical clustering.
Definition: hierarchical_clustering.hxx:598
ClusteringOptions & nodeFeatureImportance(double val)
Definition: hierarchical_clustering.hxx:642
void push(const value_type i, const priority_type p)
Insert a index with a given priority.
Definition: priority_queue.hxx:475
MetricType
Tags to select a metric for vector distance computation.
Definition: metrics.hxx:229
doxygen_overloaded_function(template<...> void separableConvolveBlockwise) template< unsigned int N
Separated convolution on ChunkedArrays.
WeightType contractionWeight()
get the edge weight of the edge which should be contracted next
Definition: hierarchical_clustering.hxx:507
void hierarchicalClustering(...)
Reduce the number of nodes in a graph by iteratively contracting the cheapest edge.
EdgeWeightNodeFeatures(MergeGraph &mergeGraph, EDGE_INDICATOR_MAP edgeIndicatorMap, EDGE_SIZE_MAP edgeSizeMap, NODE_FEATURE_MAP nodeFeatureMap, NODE_SIZE_MAP nodeSizeMap, MIN_WEIGHT_MAP minWeightEdgeMap, NODE_LABEL_MAP nodeLabelMap, const ValueType beta, const metrics::MetricType metricType, const ValueType wardness=static_cast< ValueType >(1.0), const ValueType gamma=static_cast< ValueType >(10000000.0), const ValueType sameLabelMultiplier=static_cast< ValueType >(0.8))
construct cluster operator
Definition: hierarchical_clustering.hxx:358
void mergeEdges(const Edge &a, const Edge &b)
will be called via callbacks from mergegraph
Definition: hierarchical_clustering.hxx:413
void eraseEdge(const Edge &edge)
will be called via callbacks from mergegraph
Definition: hierarchical_clustering.hxx:459
double gamma(double x)
The gamma function.
Definition: mathutil.hxx:1577
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 deleteItem(const value_type i)
deleqte the priority associated with index i
Definition: priority_queue.hxx:516
ClusteringOptions & nodeFeatureMetric(metrics::MetricType metric)
Definition: hierarchical_clustering.hxx:673
void mergeNodes(const Node &a, const Node &b)
will be called via callbacks from mergegraph
Definition: hierarchical_clustering.hxx:432
detail::SelectIntegerType< 32, detail::UnsignedIntTypes >::type UInt32
32-bit unsigned int
Definition: sized_int.hxx:183
const_reference top() const
get index with top priority
Definition: priority_queue.hxx:490
Edge contractionEdge()
get the edge which should be contracted next
Definition: hierarchical_clustering.hxx:497

© 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 (Thu Mar 17 2016)