00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026 #ifndef MLN_MORPHO_WATERSHED_TOPOLOGICAL_HH
00027 # define MLN_MORPHO_WATERSHED_TOPOLOGICAL_HH
00028
00038
00039 # include <vector>
00040 # include <map>
00041 # include <queue>
00042
00043 # include <mln/core/site_set/p_set.hh>
00044 # include <mln/core/site_set/p_priority.hh>
00045 # include <mln/core/site_set/p_queue_fast.hh>
00046
00047 # include <mln/util/ord.hh>
00048
00049 # include <mln/data/sort_psites.hh>
00050 # include <mln/data/fill.hh>
00051
00052
00053
00054
00055 namespace mln
00056 {
00057
00058 namespace morpho
00059 {
00060
00061 namespace watershed
00062 {
00063
00064 template <class I>
00065 mln_site(I) min (const Image<I> &ima_, p_set<mln_site(I)>& components)
00066 {
00067 const I& ima = exact(ima_);
00068
00069 if (components.nsites() == 0)
00070 return mln_site(I)(-1, -1);
00071
00072 typename p_set<mln_site(I)>::fwd_piter it(components);
00073 mln_site(I) min = components[0];
00074
00075 for_all(it)
00076 if (ima(it) < ima(min))
00077 min = it;
00078
00079 return min;
00080 }
00081
00082 template <class I>
00083 mln_site(I) max (p_set<mln_site(I)>& components)
00084 {
00085 if (components.nsites() == 0)
00086 return mln_site(I)(-1, -1);
00087
00088 typename p_set<mln_site(I)>::fwd_piter it(components);
00089 mln_site(I) max = components[0];
00090
00091 for_all(it)
00092 if (ima(it) > ima(max))
00093 max = it;
00094
00095 return max;
00096 }
00097
00098
00099
00100 template <class I, class N>
00101 struct topo_wst
00102 {
00103
00104
00105 typedef mln_site(I) site;
00106 typedef I image_t;
00107 typedef N neighborhood_t;
00108
00109
00110
00111 private:
00112 struct node {
00113 mln_value(I) level;
00114 int area;
00115 int highest;
00116 typedef mln_site(I) site;
00117
00118
00119 std::vector<site> children;
00120
00121 void addChildren(const node& n)
00122 {
00123
00124
00125
00126
00127
00128 for (unsigned i=0; i < n.children.size(); ++i)
00129 children.push_back(n.children[i]);
00130 }
00131
00132 void addChild(const site p)
00133 {
00134
00135 children.push_back(p);
00136 }
00137 };
00138
00139 public:
00140
00141 mln_ch_value(I, site) Par_node;
00142 mln_ch_value(I, site) Par_tree;
00143
00144 mln_ch_value(I, int) Rnk_tree;
00145 mln_ch_value(I, int) Rnk_node;
00146 mln_ch_value(I, site) subtreeRoot;
00147 mln_ch_value(I, node) nodes;
00148 site Root;
00149
00150 private:
00151
00152 void MakeSet_tree(site x);
00153 void MakeSet_node(site x);
00154 site Find_tree(site x);
00155 site Find_node(site x);
00156 void BuildComponentTree();
00157 site MergeNode(site& node1, site& node2);
00158 site Link_tree(site x, site y);
00159 site Link_node(site x, site y);
00160 node MakeNode(int level);
00161
00162 private:
00163 mln_ch_value(I, bool) isproc;
00164
00165
00166 public:
00167 topo_wst(const Image<I>& i,
00168 const Neighborhood<N>& n);
00169
00170 public:
00171 const I &ima;
00172 const N &nbh;
00173
00174 public:
00175 void go();
00176
00177 private:
00178 site min (p_set<site>& components);
00179 site max (p_set<site>& components);
00180 bool highest_fork (p_set<site>& components);
00181 bool highest_fork (p_set<site>& components, site &r);
00182
00183
00184
00185 public:
00186 site lca (site a, site b);
00187
00188 private:
00189 int *euler;
00190 int *depth;
00191 int ctree_size;
00192 std::map<site, int, mln::util::ord<site> > pos;
00193 site *repr;
00194
00195 int *minim;
00196 int **Minim;
00197
00198 void compute_ctree_size (site p);
00199 void build_euler_tour_rec(site p, int &position, int d);
00200 void build_euler_tour();
00201 void build_minim ();
00202 void removeOneSonNodes(site *p, mln_ch_value(I, site) &newPar_node);
00203 void compressTree();
00204 };
00205
00207 template <class T>
00208 typename T::image_t
00209 topological(T &tree);
00210
00211
00212
00213 # ifndef MLN_INCLUDE_ONLY
00214
00215
00216 template <class I, class N>
00217 topo_wst<I, N>::topo_wst(const Image<I>& i,
00218 const Neighborhood<N>& n)
00219 : ima(exact(i)), nbh(exact(n))
00220 {
00221 initialize(Par_node, i);
00222 initialize(Par_tree, i);
00223 initialize(Rnk_tree, i);
00224 initialize(Rnk_node, i);
00225 initialize(subtreeRoot, i);
00226 initialize(nodes, i);
00227 initialize(isproc, i);
00228 }
00229
00230 template <class I, class N>
00231 void topo_wst<I, N>::MakeSet_tree(site x)
00232 {
00233 Par_tree(x) = x;
00234 Rnk_tree(x) = 0;
00235 }
00236
00237 template <class I, class N>
00238 void topo_wst<I, N>::MakeSet_node(site x)
00239 {
00240 Par_node(x) = x;
00241 Rnk_node(x) = 0;
00242 }
00243
00244 template <class I, class N>
00245 typename topo_wst<I, N>::site topo_wst<I, N>::Find_tree(site x)
00246 {
00247 if (Par_tree(x) != x)
00248 Par_tree(x) = Find_tree(Par_tree(x));
00249 return Par_tree(x);
00250 }
00251
00252 template <class I, class N>
00253 typename topo_wst<I, N>::site topo_wst<I, N>::Find_node(site x)
00254 {
00255 if (Par_node(x) != x)
00256 Par_node(x) = Find_node(Par_node(x));
00257 return Par_node(x);
00258 }
00259
00260 template <class I, class N>
00261 void topo_wst<I, N>::go()
00262 {
00263 BuildComponentTree();
00264 compressTree();
00265
00266 build_euler_tour();
00267 build_minim();
00268 }
00269
00270 template <class I, class N>
00271 void topo_wst<I, N>::BuildComponentTree()
00272 {
00273
00274
00275
00276 p_array<mln_site(I)> S;
00277 S = data::sort_psites_increasing(ima);
00278
00279
00280 data::fill(isproc, false);
00281 for (int ip = 0; ip < int(S.nsites()); ++ip)
00282 {
00283 site p = S[ip];
00284 MakeSet_node(p);
00285 MakeSet_tree(p);
00286
00287 subtreeRoot(p) = p;
00288
00289 nodes(p) = MakeNode(ima(p));
00290 }
00291
00292
00293
00294 typename p_array<mln_site(I)>::fwd_piter ip(S);
00295 for_all(ip)
00296 {
00297 site p = ip.to_site();
00298
00299 site curCanonicalElt = Find_tree(p);
00300 site curNode = Find_node(subtreeRoot(curCanonicalElt));
00301
00302
00303 mln_niter(N) q(nbh, ip);
00304 for_all(q)
00305 if (ima.has(q) and isproc(q) and ima(q) <= ima(p))
00306 {
00307 site adjCanonicalElt = Find_tree(q);
00308 site adjNode = Find_node(subtreeRoot(adjCanonicalElt));
00309 if (curNode != adjNode)
00310 {
00311 if (nodes(curNode).level == nodes(adjNode).level)
00312 curNode = MergeNode(adjNode, curNode);
00313 else
00314 {
00315 nodes(curNode).addChild(adjNode);
00316 nodes(curNode).area += nodes(adjNode).area;
00317 nodes(curNode).highest += nodes(adjNode).highest;
00318 }
00319 }
00320
00321 curCanonicalElt = Link_tree(adjCanonicalElt, curCanonicalElt);
00322 subtreeRoot(curCanonicalElt) = curNode;
00323 }
00324 isproc(p) = true;
00325 }
00326
00327
00328
00329
00330
00331
00332
00333 mln_piter(I) r(Par_node.domain());
00334 for_all(r)
00335 Par_node(r) = Find_node(r);
00336
00337
00338
00339 mln_fwd_piter(I) rp(Par_node.domain());;
00340 rp.start();
00341
00342 Root = subtreeRoot(Find_tree(Find_node(rp)));
00343 }
00344
00345
00346 template <class I, class N>
00347 typename topo_wst<I, N>::site topo_wst<I, N>::MergeNode(site& node1,
00348 site& node2)
00349 {
00350 site tmpNode = Link_node(node1, node2);
00351 site tmpNode2;
00352 if (tmpNode == node2)
00353 {
00354 nodes(node2).addChildren(nodes(node1));
00355 tmpNode2 = node1;
00356 }
00357 else
00358 {
00359 nodes(node1).addChildren(nodes(node2));
00360 tmpNode2 = node2;
00361 }
00362 nodes(tmpNode).area += nodes(tmpNode2).area;
00363 nodes(tmpNode).highest += nodes(tmpNode2).highest;
00364 return tmpNode;
00365 }
00366
00367 template <class I, class N>
00368 typename topo_wst<I, N>::site topo_wst<I, N>::Link_tree(site x, site y)
00369 {
00370 if (Rnk_tree(x) > Rnk_tree(y))
00371 std::swap(x, y);
00372 else
00373 if (Rnk_tree(x) == Rnk_tree(y))
00374 Rnk_tree(y) += 1;
00375 Par_tree(x) = y;
00376 return y;
00377 }
00378
00379 template <class I, class N>
00380 typename topo_wst<I, N>::site topo_wst<I, N>::Link_node(site x, site y)
00381 {
00382 if (Rnk_node(x) > Rnk_node(y))
00383 std::swap(x, y);
00384 else
00385 if (Rnk_node(x) == Rnk_node(y))
00386 Rnk_node(y) += 1;
00387 Par_node(x) = y;
00388 return y;
00389 }
00390
00391 template <class I, class N>
00392 typename topo_wst<I, N>::node topo_wst<I, N>::MakeNode(int level)
00393 {
00394 node n;
00395 n.level = level;
00396 n.area = 1;
00397 n.highest = level;
00398 return n;
00399 }
00400
00401
00402 template <class I, class N>
00403 bool topo_wst<I, N>::highest_fork (p_set<site>& components, site &r)
00404 {
00405 if (components.nsites() == 0)
00406 {
00407 std::cerr << "highest fork : empty set" << std::endl;
00408 return false;
00409 }
00410
00411 site
00412 m = min(components),
00413 hfork = m;
00414
00415 typename p_set<site>::fwd_piter it(components);
00416 for_all(it)
00417 {
00418
00419 if (it.to_site() == m)
00420 continue;
00421
00422 site c = lca(hfork, it.to_site());
00423 if (c != it.to_site())
00424 hfork = c;
00425 }
00426
00427 if (nodes(m).level == nodes(hfork).level)
00428 return false;
00429
00430 r = hfork;
00431 return true;
00432 }
00433
00434 template <class I, class N>
00435 bool topo_wst<I, N>::highest_fork (p_set<site>& components) {
00436 site i;
00437 return highest_fork(components, i);
00438 }
00439
00440 template <class I, class N>
00441 void topo_wst<I, N>::compute_ctree_size (site p)
00442 {
00443 ctree_size += 1;
00444 node& n = nodes(p);
00445
00446
00447
00448
00449
00450 for (unsigned i=0; i < n.children.size(); ++i)
00451 compute_ctree_size(n.children[i]);
00452 }
00453
00454
00455 template <class I, class N>
00456 void topo_wst<I, N>::build_euler_tour_rec(site p, int &position, int d)
00457 {
00458 if (pos.find(p) == pos.end())
00459 pos[p] = position;
00460
00461 repr[position] = p;
00462 depth[position] = d;
00463 euler[position] = pos[p];
00464 ++position;
00465 node& n = nodes(p);
00466
00467
00468
00469
00470
00471
00472
00473
00474
00475
00476
00477 for (unsigned i=0; i < n.children.size(); ++i)
00478 {
00479 build_euler_tour_rec(n.children[i], position, d+1);
00480 depth[position] = d;
00481 euler[position] = pos[p];
00482 repr[position] = p;
00483 ++position;
00484 }
00485 }
00486
00487
00488 template <class I, class N>
00489 void topo_wst<I, N>::build_euler_tour ()
00490 {
00491 ctree_size = 0;
00492 compute_ctree_size(Root);
00493
00494 int size = 2 * ctree_size - 1;
00495
00496
00497 euler = new int[size];
00498 depth = new int[size];
00499 repr = new site[size];
00500
00501 int position = 0;
00502 build_euler_tour_rec(Root, position, 0);
00503 }
00504
00505
00506 template <class I, class N>
00507 void topo_wst<I, N>::build_minim ()
00508 {
00509
00510 int size = 2 * ctree_size - 1;
00511 int logn = (int)(ceil(log((double)(size))/log(2.0)));
00512
00513 minim = new int [logn * size];
00514 Minim = new int* [logn];
00515
00516 Minim[0] = minim;
00517
00518
00519 for (int i = 0; i < size - 1; ++i)
00520 if (depth[euler[i]] < depth[euler[i+1]]) {
00521 Minim[0][i] = i;
00522 } else {
00523 Minim[0][i] = i+1;
00524 }
00525
00526
00527
00528 int k;
00529 for (int j = 1; j < logn; j++) {
00530 k = 1 << (j - 1);
00531 Minim[j] = &minim[j * size];
00532 for (int i = 0; i < size; i++) {
00533 if ((i + (k << 1)) >= size) {
00534
00535 }
00536 else {
00537 if (depth[euler[Minim[j - 1][i]]]
00538 <= depth[euler[Minim[j - 1][i + k]]])
00539 Minim[j][i] = Minim[j - 1][i];
00540 else
00541 Minim[j][i] = Minim[j - 1][i + k];
00542 }
00543 }
00544 }
00545
00546 }
00547
00548
00549 template <class I, class N>
00550 typename topo_wst<I, N>::site topo_wst<I, N>::lca (site a, site b)
00551 {
00552 int
00553 m = pos[a],
00554 n = pos[b],
00555 k;
00556
00557 if (m == n)
00558 return repr[m];
00559
00560 if (m > n)
00561 {
00562 k = n;
00563 n = m;
00564 m = k;
00565 }
00566
00567 k = (int)(log((double)(n - m))/log(2.));
00568
00569 if (depth[euler[Minim[k][m]]] < depth[euler[Minim[k][n - (1 << k)]]]) {
00570 return repr[euler[Minim[k][m]]];
00571 } else {
00572 return repr[euler[Minim[k][n - (1 << k)]]];
00573 }
00574 }
00575
00576
00577 template <class I, class N>
00578 void topo_wst<I, N>::removeOneSonNodes(site *p,
00579 mln_ch_value(I, site) &newPar_node)
00580 {
00581 node &n = nodes(*p);
00582
00583 if (n.children.size() == 1)
00584 {
00585 n.area = -1;
00586 newPar_node(*p) = n.children[0];
00587 *p = n.children[0];
00588 removeOneSonNodes(p, newPar_node);
00589 }
00590 else
00591 {
00592 for (unsigned i = 0; i < n.children.size(); ++i)
00593 removeOneSonNodes(&(n.children[i]), newPar_node);
00594 }
00595 }
00596
00597
00598 template <class I, class N>
00599 void topo_wst<I, N>::compressTree()
00600 {
00601 mln_ch_value(I, site) newPar_node;
00602 initialize(newPar_node, Par_node);
00603
00604
00605 removeOneSonNodes(&Root, newPar_node);
00606
00607
00608 mln_piter(I) p(Par_node.domain());
00609 for_all(p)
00610 while (nodes(Par_node(p)).area == -1)
00611 Par_node(p) = newPar_node(Par_node(p));
00612 }
00613
00614 template <class T>
00615 bool w_constructible(T &tree, typename T::site p, typename T::site &r)
00616 {
00617
00618 typedef typename T::image_t I;
00619 typedef typename T::neighborhood_t N;
00620
00621 const I &ima = exact(tree.ima);
00622 const N &nbh = exact(tree.nbh);
00623
00624
00625 mln_niter(N) q(nbh, p);
00626 p_set<mln_site(I)> v;
00627
00628 for_all(q)
00629
00630
00631 if (ima.domain().has(q) && ima(q) < ima(p))
00632 v.insert(tree.Par_node(q));
00633
00634 if (v.nsites() == 0)
00635 return false;
00636
00637 if (v.nsites() == 1)
00638 {
00639 r = v[0];
00640 return true;
00641 }
00642
00643 mln_site(I) c = min(ima, v);
00644 mln_site(I) cmin = c;
00645
00646 typename p_set<mln_site(I)>::fwd_piter it(v);
00647 for_all(it)
00648 {
00649
00650 if (it.to_site() == cmin)
00651 continue;
00652
00653 mln_site(I) c1 = tree.lca(c, it);
00654
00655 if (c1 != it)
00656 c = c1;
00657 }
00658
00659 if (tree.nodes(c).level >= ima(p))
00660 return false;
00661
00662 r = c;
00663 return true;
00664 }
00665
00666 template <class T>
00667 bool w_constructible(T &tree, typename T::site p) {
00668 typename T::site r;
00669 return w_constructible(tree, p, r);
00670 }
00671
00672
00673
00674 template <class T>
00675 typename T::image_t topological(T &tree)
00676 {
00677
00678 typedef typename T::image_t I;
00679 typedef typename T::neighborhood_t N;
00680
00681 I ima = exact(tree.ima);
00682 const N &nbh = exact(tree.nbh);
00683
00684
00685 mln_ch_value(I, bool) cmax;
00686 initialize(cmax, ima);
00687
00688
00689 mln_ch_value(I, bool) enqueued;
00690 initialize(enqueued, ima);
00691
00692 p_priority< mln_value(I), p_queue_fast<mln_site(I)> > l;
00693
00694 std::queue<mln_site(I)> m;
00695 mln_value(I) min = mln_min(mln_value(I));
00696
00697 std::cout << "Init" << std::endl;
00698
00699
00700 data::fill(cmax, false);
00701 data::fill(enqueued, false);
00702
00703 mln_piter(I) it(tree.Par_node.domain());
00704 for_all(it)
00705
00706 if (tree.nodes(tree.Par_node(it)).children.size() == 0)
00707 {
00708 cmax(it) = true;
00709 m.push(it);
00710 }
00711
00712 while (!m.empty())
00713 {
00714
00715 mln_niter(N) q(nbh, m.front());
00716
00717
00718 for_all(q)
00719 if (cmax.domain().has(q) && !cmax(q) && !enqueued(q))
00720 {
00721 enqueued(q) = true;
00722 l.push(mln_max(mln_value(I)) - ima(q), q);
00723 }
00724 m.pop();
00725 }
00726
00727
00728 std::cout << "end" << std::endl;
00729
00730
00731 while (!l.is_empty())
00732 {
00733 mln_site(I) x = l.front();
00734 l.pop();
00735 enqueued(x) = false;
00736
00737 mln_site(I) c;
00738 bool is_w = w_constructible(tree, x, c);
00739
00740 if (is_w)
00741 {
00742 ima(x) = tree.nodes(c).level;
00743 tree.Par_node(x) = c;
00744
00745
00746 if (tree.nodes(c).children.size() == 0)
00747 cmax(x) = true;
00748 else
00749
00750 if (tree.nodes(c).children.size() == 1)
00751 std::cerr << "ERREUR COMPOSANTE BRANCHE "
00752 << tree.nodes(c).children.size() << std::endl;
00753
00754
00755
00756 mln_niter(N) q(nbh, x);
00757
00758
00759 for_all(q)
00760 if (ima.domain().has(q) && !cmax(q) && !enqueued(q))
00761 {
00762 enqueued(q) = true;
00763 l.push(mln_max(mln_value(I)) - ima(q), q);
00764
00765
00766
00767
00768 }
00769 }
00770 }
00771
00772 for_all(it)
00773 ima(it) = tree.nodes(tree.Par_node(it)).level;
00774
00775 return ima;
00776 }
00777
00778 # endif // MLN_INCLUDE_ONLY
00779
00780
00781 }
00782
00783 }
00784
00785 }
00786
00787 #endif // ! MLN_MORPHO_WATERSHED_TOPOLOGICAL_HH