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

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

© 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.2