Milena (Olena)
User documentation 2.0a Id
|
00001 // Copyright (C) 2008, 2009 EPITA Research and Development Laboratory (LRDE) 00002 // 00003 // This file is part of Olena. 00004 // 00005 // Olena is free software: you can redistribute it and/or modify it under 00006 // the terms of the GNU General Public License as published by the Free 00007 // Software Foundation, version 2 of the License. 00008 // 00009 // Olena is distributed in the hope that it will be useful, 00010 // but WITHOUT ANY WARRANTY; without even the implied warranty of 00011 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 00012 // General Public License for more details. 00013 // 00014 // You should have received a copy of the GNU General Public License 00015 // along with Olena. If not, see <http://www.gnu.org/licenses/>. 00016 // 00017 // As a special exception, you may use this file as part of a free 00018 // software project without restriction. Specifically, if other files 00019 // instantiate templates or use macros or inline functions from this 00020 // file, or you compile this file and link it with other files to produce 00021 // an executable, this file does not by itself cause the resulting 00022 // executable to be covered by the GNU General Public License. This 00023 // exception does not however invalidate any other reasons why the 00024 // executable file might be covered by the GNU General Public License. 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 // FIXME: Clean up and document. 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); // FIXME 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); // FIXME 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 // Actually, this structure is a tree, despite its confusing name. 00100 template <class I, class N> 00101 struct topo_wst 00102 { 00103 // Typedefs 00104 // -------- 00105 typedef mln_site(I) site; 00106 typedef I image_t; 00107 typedef N neighborhood_t; 00108 00109 // Component tree management 00110 // ------------------------- 00111 private: 00112 struct node { 00113 mln_value(I) level; 00114 int area; 00115 int highest; 00116 typedef mln_site(I) site; 00117 // Can't modify the sites in a p_array 00118 // p_array<mln_site(I)> children; 00119 std::vector<site> children; 00120 00121 void addChildren(const node& n) 00122 { 00123 // typename p_array<mln_site(I)>::fwd_piter it(n.children); 00124 // for (it.start(); 00125 // it.is_valid(); 00126 // it.next()) 00127 // children.append(it.to_site()); 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 // children.append(n); 00135 children.push_back(p); 00136 } 00137 }; // struct node 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 // Ctor 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 // Optimized LCA Algorithm 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 // Ctor 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 // Init: 00274 00275 // Sort the sites in increasing order 00276 p_array<mln_site(I)> S; 00277 S = data::sort_psites_increasing(ima); 00278 00279 // Clear the marker map 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 // if (subtreeRoot.hold(p)) 00287 subtreeRoot(p) = p; 00288 // if (nodes.hold(p)) 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 // FIXME: Should be `n' instead of `q'. 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 // Pour garder une map de correspondance site <-> local_root 00327 // for (int ip = 0; ip < int(S.size()); ++ip) 00328 // { 00329 // site p = S[ip]; 00330 // M(p) = Find_node(p); 00331 // } 00332 00333 mln_piter(I) r(Par_node.domain()); 00334 for_all(r) 00335 Par_node(r) = Find_node(r); 00336 00337 // Find the ``first'' site of ima, according to the forward 00338 // traversal order. 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 // Can't remove the site from the set 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 // typename p_array<mln_site(I)>::fwd_piter it(n.children); 00447 // for_all(it) 00448 // compute_ctree_size(it.to_site()); 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 // typename p_array<mln_site(I)>::fwd_piter it(n.children); 00468 // for_all(it) 00469 // { 00470 // build_euler_tour_rec(it.to_site(), position, d+1); 00471 // depth[position] = d; // Not optimized 00472 // euler[position] = pos[p]; 00473 // repr[position] = p; // Pas necessaire? 00474 // ++position; 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; // Not optimized 00481 euler[position] = pos[p]; 00482 repr[position] = p; // Pas necessaire? 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 // FIXME : free this 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 // compute_tree_size(); // Already done in build_euler_tour 00510 int size = 2 * ctree_size - 1; 00511 int logn = (int)(ceil(log((double)(size))/log(2.0))); 00512 // minim.reserve(size); 00513 minim = new int [logn * size]; // FIXME : Think about freeing this 00514 Minim = new int* [logn]; 00515 00516 Minim[0] = minim; 00517 00518 // Init 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 // Minim[0][size - 1] = size - 1; 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 //Minim[j][i] = size - 1; 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 } // void build_minim () 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) // this node has 1 son, delete it 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 // there is more than one son, recursive call 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 // Remove the nodes with one son 00605 removeOneSonNodes(&Root, newPar_node); 00606 00607 // Update the references on deleted nodes 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 // FIXME: Should be `n' instead of `q'. 00625 mln_niter(N) q(nbh, p); 00626 p_set<mln_site(I)> v; 00627 00628 for_all(q) 00629 // FIXME: Shouldn't it be: `ima.has(q)' instead of 00630 // `ima.domain().has(q)'? 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 // Can't remove the site from the set 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 // Maxima components 00685 mln_ch_value(I, bool) cmax; 00686 initialize(cmax, ima); 00687 00688 // Mark enqueued sites 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 // p_queue < site > m; 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 // Flag C-maxima 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 // if (nodes(Par_node(it.to_site())).children.nsites() == 0) 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 // FIXME: Should be `n' instead of `q'. 00715 mln_niter(N) q(nbh, m.front()); 00716 // FIXME: Shouldn't it be: `cmax.has(q)' instead of 00717 // `cmax.domain().has(q)'? 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 // Main loop 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 // if (nodes(c).children.nsites() == 0) 00746 if (tree.nodes(c).children.size() == 0) 00747 cmax(x) = true; 00748 else 00749 // if (nodes(c).children.nsites() > 1) 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 // FIXME: Should be `n' instead of `q'. 00756 mln_niter(N) q(nbh, x); 00757 // FIXME: Shouldn't it be: `ima.has(q)' instead of 00758 // `ima.domain().has(q)'? 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); // FIXME: 00764 // Just 00765 // invert 00766 // the 00767 // priority. 00768 } 00769 } 00770 } // while(!l.empty()) 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 } // end of namespace mln::morpho::watershed 00782 00783 } // end of namespace mln::morpho 00784 00785 } // end of namespace mln 00786 00787 #endif // ! MLN_MORPHO_WATERSHED_TOPOLOGICAL_HH