[ 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 
51 namespace vigra{
52 
53 namespace cluster_operators{
54 
55  /// \brief get minimum edge weight from an edge indicator and difference of node features
56  template<
57  class MERGE_GRAPH,
58  class EDGE_INDICATOR_MAP,
59  class EDGE_SIZE_MAP,
60  class NODE_FEATURE_MAP,
61  class NODE_SIZE_MAP,
62  class MIN_WEIGHT_MAP
63  >
65 
66  typedef EdgeWeightNodeFeatures<
67  MERGE_GRAPH,
68  EDGE_INDICATOR_MAP,
69  EDGE_SIZE_MAP,
70  NODE_FEATURE_MAP,
71  NODE_SIZE_MAP,
72  MIN_WEIGHT_MAP
73  > SelfType;
74  public:
75 
76 
77  typedef typename EDGE_INDICATOR_MAP::Value ValueType;
78  typedef ValueType WeightType;
79  typedef MERGE_GRAPH MergeGraph;
80  typedef typename MergeGraph::Graph Graph;
81  typedef typename Graph::Edge BaseGraphEdge;
82  typedef typename Graph::Node BaseGraphNode;
83  typedef typename MergeGraph::Edge Edge;
84  typedef typename MergeGraph::Node Node;
85  typedef typename MergeGraph::EdgeIt EdgeIt;
86  typedef typename MergeGraph::NodeIt NodeIt;
87  typedef typename MergeGraph::IncEdgeIt IncEdgeIt;
88  typedef typename MergeGraph::index_type index_type;
89  typedef MergeGraphItemHelper<MergeGraph,Edge> EdgeHelper;
90  typedef MergeGraphItemHelper<MergeGraph,Node> NodeHelper;
91 
92 
93  typedef typename EDGE_INDICATOR_MAP::Reference EdgeIndicatorReference;
94  typedef typename NODE_FEATURE_MAP::Reference NodeFeatureReference;
95  /// \brief construct cluster operator
97  MergeGraph & mergeGraph,
98  EDGE_INDICATOR_MAP edgeIndicatorMap,
99  EDGE_SIZE_MAP edgeSizeMap,
100  NODE_FEATURE_MAP nodeFeatureMap,
101  NODE_SIZE_MAP nodeSizeMap,
102  MIN_WEIGHT_MAP minWeightEdgeMap,
103  const ValueType beta,
104  const metrics::MetricType metricType,
105  const ValueType wardness=1.0
106  )
107  : mergeGraph_(mergeGraph),
108  edgeIndicatorMap_(edgeIndicatorMap),
109  edgeSizeMap_(edgeSizeMap),
110  nodeFeatureMap_(nodeFeatureMap),
111  nodeSizeMap_(nodeSizeMap),
112  minWeightEdgeMap_(minWeightEdgeMap),
113  pq_(mergeGraph.maxEdgeId()+1),
114  beta_(beta),
115  wardness_(wardness),
116  metric_(metricType)
117  {
118  typedef typename MergeGraph::MergeNodeCallBackType MergeNodeCallBackType;
119  typedef typename MergeGraph::MergeEdgeCallBackType MergeEdgeCallBackType;
120  typedef typename MergeGraph::EraseEdgeCallBackType EraseEdgeCallBackType;
121 
122 
123  MergeNodeCallBackType cbMn(MergeNodeCallBackType:: template from_method<SelfType,&SelfType::mergeNodes>(this));
124  MergeEdgeCallBackType cbMe(MergeEdgeCallBackType:: template from_method<SelfType,&SelfType::mergeEdges>(this));
125  EraseEdgeCallBackType cbEe(EraseEdgeCallBackType:: template from_method<SelfType,&SelfType::eraseEdge>(this));
126 
127  mergeGraph_.registerMergeNodeCallBack(cbMn);
128  mergeGraph_.registerMergeEdgeCallBack(cbMe);
129  mergeGraph_.registerEraseEdgeCallBack(cbEe);
130 
131 
132  for(EdgeIt e(mergeGraph);e!=lemon::INVALID;++e){
133  const Edge edge = *e;
134  const BaseGraphEdge graphEdge=EdgeHelper::itemToGraphItem(mergeGraph_,edge);
135  const index_type edgeId = mergeGraph_.id(edge);
136  const ValueType currentWeight = this->getEdgeWeight(edge);
137  pq_.push(edgeId,currentWeight);
138  minWeightEdgeMap_[graphEdge]=currentWeight;
139  }
140  }
141 
142  /// \brief will be called via callbacks from mergegraph
143  void mergeEdges(const Edge & a,const Edge & b){
144  // update features / weigts etc
145  const BaseGraphEdge aa=EdgeHelper::itemToGraphItem(mergeGraph_,a);
146  const BaseGraphEdge bb=EdgeHelper::itemToGraphItem(mergeGraph_,b);
147  EdgeIndicatorReference va=edgeIndicatorMap_[aa];
148  EdgeIndicatorReference vb=edgeIndicatorMap_[bb];
149  va*=edgeSizeMap_[aa];
150  vb*=edgeSizeMap_[bb];
151  va+=vb;
152  edgeSizeMap_[aa]+=edgeSizeMap_[bb];
153  va/=(edgeSizeMap_[aa]);
154  vb/=edgeSizeMap_[bb];
155  // delete b from pq
156  pq_.deleteItem(b.id());
157  }
158 
159  /// \brief will be called via callbacks from mergegraph
160  void mergeNodes(const Node & a,const Node & b){
161  const BaseGraphNode aa=NodeHelper::itemToGraphItem(mergeGraph_,a);
162  const BaseGraphNode bb=NodeHelper::itemToGraphItem(mergeGraph_,b);
163  NodeFeatureReference va=nodeFeatureMap_[aa];
164  NodeFeatureReference vb=nodeFeatureMap_[bb];
165  va*=nodeSizeMap_[aa];
166  vb*=nodeSizeMap_[bb];
167  va+=vb;
168  nodeSizeMap_[aa]+=nodeSizeMap_[bb];
169  va/=(nodeSizeMap_[aa]);
170  vb/=nodeSizeMap_[bb];
171  }
172 
173  /// \brief will be called via callbacks from mergegraph
174  void eraseEdge(const Edge & edge){
175 
176  //std::cout<<"start to erase edge "<<mergeGraph_.id(edge)<<"\n";
177  // delete edge from pq
178  pq_.deleteItem(edge.id());
179  // get the new region the edge is in
180  // (since the edge is no any more an active edge)
181  //std::cout<<"get the new node \n";
182  const Node newNode = mergeGraph_.inactiveEdgesNode(edge);
183  //std::cout<<"new node "<<mergeGraph_.id(newNode)<<"\n";
184 
185  size_t counter=0;
186  // iterate over all edges of this node
187  for (IncEdgeIt e(mergeGraph_,newNode);e!=lemon::INVALID;++e){
188 
189  //std::cout<<"get inc edge\n";
190  const Edge incEdge(*e);
191 
192  //std::cout<<"get inc graph edge\n";
193  const BaseGraphEdge incGraphEdge = EdgeHelper::itemToGraphItem(mergeGraph_,incEdge);
194 
195  //std::cout<<"get inc edge weight"<<counter<<"\n";
196  // compute the new weight for this edge
197  // (this should involve region differences)
198  const ValueType newWeight = getEdgeWeight(incEdge);
199  // change the weight in pq by repushing
200 
201  //std::cout<<"push\n";
202  pq_.push(incEdge.id(),newWeight);
203  // remember edge weight
204 
205  //std::cout<<"set new\n";
206  minWeightEdgeMap_[incGraphEdge]=newWeight;
207  ++counter;
208  }
209  //std::cout<<"done\n";
210  }
211 
212  /// \brief get the edge which should be contracted next
214  index_type minLabel = pq_.top();
215  while(mergeGraph_.hasEdgeId(minLabel)==false){
216  pq_.deleteItem(minLabel);
217  minLabel = pq_.top();
218  }
219  return Edge(minLabel);
220  }
221 
222  /// \brief get the edge weight of the edge which should be contracted next
223  WeightType contractionWeight(){
224  index_type minLabel = pq_.top();
225  while(mergeGraph_.hasEdgeId(minLabel)==false){
226  pq_.deleteItem(minLabel);
227  minLabel = pq_.top();
228  }
229  return pq_.topPriority();
230 
231  }
232 
233 
234  /// \brief get a reference to the merge
235  MergeGraph & mergeGraph(){
236  return mergeGraph_;
237  }
238  private:
239  ValueType getEdgeWeight(const Edge & e){
240 
241  const Node u = mergeGraph_.u(e);
242  const Node v = mergeGraph_.v(e);
243 
244  const BaseGraphEdge ee=EdgeHelper::itemToGraphItem(mergeGraph_,e);
245  const BaseGraphNode uu=NodeHelper::itemToGraphItem(mergeGraph_,u);
246  const BaseGraphNode vv=NodeHelper::itemToGraphItem(mergeGraph_,v);
247 
248  const ValueType wardFacRaw = 1.0 / ( 1.0/std::log(nodeSizeMap_[uu]) + 1.0/std::log(nodeSizeMap_[vv]) );
249  const ValueType wardFac = (wardFacRaw*wardness_) + (1.0-wardness_);
250 
251  const ValueType fromEdgeIndicator = edgeIndicatorMap_[ee];
252  ValueType fromNodeDist = metric_(nodeFeatureMap_[uu],nodeFeatureMap_[vv]);
253  const ValueType totalWeight = ((1.0-beta_)*fromEdgeIndicator + beta_*fromNodeDist)*wardFac;
254  return totalWeight;
255  }
256 
257 
258  MergeGraph & mergeGraph_;
259  EDGE_INDICATOR_MAP edgeIndicatorMap_;
260  EDGE_SIZE_MAP edgeSizeMap_;
261  NODE_FEATURE_MAP nodeFeatureMap_;
262  NODE_SIZE_MAP nodeSizeMap_;
263  MIN_WEIGHT_MAP minWeightEdgeMap_;
265  ValueType beta_;
266  ValueType wardness_;
267 
268  metrics::Metric<float> metric_;
269  };
270 } // end namespace cluster_operators
271 
272 
273 
274  /// \brief do hierarchical clustering with a given cluster operator
275  template< class CLUSTER_OPERATOR>
277 
278  public:
279  typedef CLUSTER_OPERATOR ClusterOperator;
280  typedef typename ClusterOperator::MergeGraph MergeGraph;
281  typedef typename MergeGraph::Graph Graph;
282  typedef typename Graph::Edge BaseGraphEdge;
283  typedef typename Graph::Node BaseGraphNode;
284  typedef typename MergeGraph::Edge Edge;
285  typedef typename MergeGraph::Node Node;
286  typedef typename CLUSTER_OPERATOR::WeightType ValueType;
287  typedef typename MergeGraph::index_type MergeGraphIndexType;
288 
289  struct Parameter{
290  Parameter(
291  const size_t nodeNumStopCond = 1,
292  const bool buildMergeTree = true,
293  const bool verbose = false
294  )
295  : nodeNumStopCond_ (nodeNumStopCond),
296  buildMergeTreeEncoding_(buildMergeTree),
297  verbose_(verbose){
298  }
299  size_t nodeNumStopCond_;
300  bool buildMergeTreeEncoding_;
301  bool verbose_;
302  };
303 
304  struct MergeItem{
305  MergeItem(
306  const MergeGraphIndexType a,
307  const MergeGraphIndexType b,
308  const MergeGraphIndexType r,
309  const ValueType w
310  ):
311  a_(a),b_(b),r_(r),w_(w){
312  }
313  MergeGraphIndexType a_;
314  MergeGraphIndexType b_;
315  MergeGraphIndexType r_;
316  ValueType w_;
317  };
318 
319  typedef std::vector<MergeItem> MergeTreeEncoding;
320 
321  /// \brief construct HierarchicalClustering from clusterOperator and an optional parameter object
323  ClusterOperator & clusterOperator,
324  const Parameter & parameter = Parameter()
325  )
326  :
327  clusterOperator_(clusterOperator),
328  param_(parameter),
329  mergeGraph_(clusterOperator_.mergeGraph()),
330  graph_(mergeGraph_.graph()),
331  timestamp_(graph_.maxNodeId()+1),
332  toTimeStamp_(graph_.maxNodeId()+1),
333  timeStampIndexToMergeIndex_(graph_.maxNodeId()+1),
334  mergeTreeEndcoding_()
335  {
336  // this can be be made smater since user can pass
337  // stoping condition based on nodeNum
338  mergeTreeEndcoding_.reserve(graph_.nodeNum()*2);
339 
340  for(MergeGraphIndexType nodeId=0;nodeId<=mergeGraph_.maxNodeId();++nodeId){
341  toTimeStamp_[nodeId]=nodeId;
342  }
343  }
344 
345  /// \brief start the clustering
346  void cluster(){
347  if(param_.verbose_)
348  std::cout<<"\n";
349  while(mergeGraph_.nodeNum()>param_.nodeNumStopCond_ && mergeGraph_.edgeNum()>0){
350 
351 
352  const Edge edgeToRemove = clusterOperator_.contractionEdge();
353  if(param_.buildMergeTreeEncoding_){
354  const MergeGraphIndexType uid = mergeGraph_.id(mergeGraph_.u(edgeToRemove));
355  const MergeGraphIndexType vid = mergeGraph_.id(mergeGraph_.v(edgeToRemove));
356  const ValueType w = clusterOperator_.contractionWeight();
357  // do the merge
358  mergeGraph_.contractEdge( edgeToRemove);
359  const MergeGraphIndexType aliveNodeId = mergeGraph_.hasNodeId(uid) ? uid : vid;
360  const MergeGraphIndexType deadNodeId = aliveNodeId==vid ? uid : vid;
361  timeStampIndexToMergeIndex_[timeStampToIndex(timestamp_)]=mergeTreeEndcoding_.size();
362  mergeTreeEndcoding_.push_back(MergeItem( toTimeStamp_[aliveNodeId],toTimeStamp_[deadNodeId],timestamp_,w));
363  toTimeStamp_[aliveNodeId]=timestamp_;
364  timestamp_+=1;
365  }
366  else{
367  // do the merge
368  mergeGraph_.contractEdge( edgeToRemove );
369  }
370  if(param_.verbose_ && mergeGraph_.nodeNum()%10==0)
371  std::cout<<"\rNodes: "<<std::setw(10)<<mergeGraph_.nodeNum()<<std::flush;
372 
373  }
374  if(param_.verbose_)
375  std::cout<<"\n";
376  }
377 
378  /// \brief get the encoding of the merge tree
379  const MergeTreeEncoding & mergeTreeEndcoding()const{
380  return mergeTreeEndcoding_;
381  }
382 
383  /// \brief get the node id's which are the leafes of a treeNodeId
384  template<class OUT_ITER>
385  size_t leafNodeIds(const MergeGraphIndexType treeNodeId, OUT_ITER begin)const{
386  if(treeNodeId<=graph_.maxNodeId()){
387  *begin=treeNodeId;
388  ++begin;
389  return 1;
390  }
391  else{
392  size_t leafNum=0;
393  std::queue<MergeGraphIndexType> queue;
394  queue.push(treeNodeId);
395 
396  while(!queue.empty()){
397 
398  const MergeGraphIndexType id = queue.front();
399  queue.pop();
400  const MergeGraphIndexType mergeIndex = timeStampToMergeIndex(id);
401  const MergeGraphIndexType ab[]= { mergeTreeEndcoding_[mergeIndex].a_, mergeTreeEndcoding_[mergeIndex].b_};
402 
403  for(size_t i=0;i<2;++i){
404  if(ab[i]<=graph_.maxNodeId()){
405  *begin=ab[i];
406  ++begin;
407  ++leafNum;
408  }
409  else{
410  queue.push(ab[i]);
411  }
412  }
413  }
414  return leafNum;
415  }
416  }
417 
418  /// \brief get the graph the merge graph is based on
419  const Graph & graph()const{
420  return graph_;
421  }
422 
423  /// \brief get the merge graph
424  const MergeGraph & mergeGraph()const{
425  return mergeGraph_;
426  }
427 
428  /// \brief get the representative node id
429  const MergeGraphIndexType reprNodeId(const MergeGraphIndexType id)const{
430  return mergeGraph_.reprNodeId(id);
431  }
432  private:
433 
434  MergeGraphIndexType timeStampToIndex(const MergeGraphIndexType timestamp)const{
435  return timestamp- graph_.maxNodeId();
436  }
437 
438 
439  MergeGraphIndexType timeStampToMergeIndex(const MergeGraphIndexType timestamp)const{
440  return timeStampIndexToMergeIndex_[timeStampToIndex(timestamp)];
441  }
442 
443  ClusterOperator & clusterOperator_;
444  Parameter param_;
445  MergeGraph & mergeGraph_;
446  const Graph & graph_;
447  // parameter object
448 
449 
450  // timestamp
451  MergeGraphIndexType timestamp_;
452  std::vector<MergeGraphIndexType> toTimeStamp_;
453  std::vector<MergeGraphIndexType> timeStampIndexToMergeIndex_;
454  // data which can reconstruct the merge tree
455  MergeTreeEncoding mergeTreeEndcoding_;
456 
457 
458  };
459 
460 
461 }
462 
463 #endif // VIGRA_HIERARCHICAL_CLUSTERING_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)