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_TREE_IMLP_DUAL_UNION_FIND_HH
00027 # define MLN_MORPHO_TREE_IMLP_DUAL_UNION_FIND_HH
00028
00040
00041 # include <mln/data/fill.hh>
00042
00043 # include <mln/geom/nsites.hh>
00044 # include <mln/geom/translate.hh>
00045
00046 # include <mln/morpho/tree/data.hh>
00047
00048 # include <mln/util/timer.hh>
00049
00050 namespace mln
00051 {
00052
00053 namespace morpho
00054 {
00055
00056 namespace tree
00057 {
00058
00059 namespace impl
00060 {
00061
00062 namespace generic
00063 {
00064
00075 template <typename I, typename S, typename N>
00076 data< I, p_array<mln_psite(I)> >
00077 dual_union_find(const Image<I>& f,
00078 const Image<I>& m,
00079 const Site_Set<S>& s_f,
00080 const Site_Set<S>& s_m,
00081 const Neighborhood<N>& nbh);
00082
00083 }
00084
00085 }
00086
00087 # ifndef MLN_INCLUDE_ONLY
00088
00089 namespace internal
00090 {
00091
00092
00093 static util::timer t_prop;
00094
00095
00096 template <typename I>
00097 inline
00098 mln_psite(I) find_root(I& zpar,
00099 const mln_psite(I)& p)
00100 {
00101 if (zpar(p) == p)
00102 return p;
00103
00104 t_prop.resume();
00105 mln_psite(I) x = zpar(p);
00106 while (zpar(x) != x)
00107 x = zpar(x);
00108
00109 mln_psite(I) tmp;
00110 for (mln_psite(I) y = p; y != x; y = tmp)
00111 {
00112 tmp = zpar(y);
00113 zpar(y) = x;
00114 }
00115 t_prop.stop();
00116
00117 return x;
00118 }
00119
00120 template <typename I>
00121 inline
00122 void
00123 update_m_parent(const I& f,
00124 mln_ch_value(I, mln_psite(I))& parent,
00125 mln_ch_value(I, bool)& deja_vu,
00126 p_array< mln_psite(I) >& sset,
00127 const mln_domain(I)& d_ext,
00128 const mln_psite(I)& p)
00129 {
00130 typedef mln_psite(I) P;
00131
00132 P q = parent(p);
00133 P x = parent(q);
00134
00135 mln_assertion(d_ext.has(q));
00136
00137 while (d_ext.has(x) && f(q) == f(x) && q != x)
00138 {
00139 q = x;
00140 x = parent(q);
00141 }
00142
00143 if (!d_ext.has(x))
00144 {
00145 if (f(x) == f(parent(x)))
00146 x = (parent(q) = parent(x));
00147 if (f(q) != f(x))
00148 {
00149 x = q;
00150 if (!deja_vu(q))
00151 {
00152 deja_vu(q) = true;
00153 sset.append(q);
00154 }
00155 }
00156
00157 }
00158 else
00159 {
00160 if (x != q)
00161 {
00162 update_m_parent(f, parent, deja_vu, sset, d_ext, q);
00163 x = q;
00164 }
00165 if (!deja_vu(q))
00166 {
00167 deja_vu(q) = true;
00168 sset.append(q);
00169 }
00170 }
00171
00172 for (P i = p, tmp = parent(i); i != q; i = tmp, tmp = parent(i))
00173 parent(i) = x;
00174 }
00175
00176 }
00177
00178 namespace impl
00179 {
00180
00181 namespace generic
00182 {
00183
00184
00185 template <typename I, typename S, typename N>
00186 data< I, p_array<mln_psite(I)> >
00187 dual_union_find(const Image<I>& f_,
00188 const Image<I>& m_,
00189 const Site_Set<S>& s_f_,
00190 const Site_Set<S>& s_m_,
00191 const Neighborhood<N>& nbh_)
00192 {
00193 trace::entering("morpho::tree::impl::generic::dual_union_find");
00194
00195 util::timer tm;
00196 tm.start();
00197 internal::t_prop.reset();
00198
00199 typedef mln_psite(I) P;
00200 typedef unsigned L;
00201
00202 const I& f = exact(f_);
00203 const I& m = exact(m_);
00204 const S& s_f = exact(s_f_);
00205 const S& s_m = exact(s_m_);
00206 const N& nbh = exact(nbh_);
00207
00208
00209 mln_psite(I)::delta dp(literal::zero);
00210 mln_domain(I) d_f = f.domain();
00211 mln_domain(I) d_ext = f.domain();
00212 mln_domain(I) d = f.domain();
00213
00214
00215 dp[0] = d.pmax()[0] - d.pmin()[0] + 1;
00216 d.pmax() += dp;
00217 d_ext.pmin() += dp;
00218 d_ext.pmax() += dp;
00219
00220
00221 mln_concrete(I) fext;
00222 mln_ch_value(I, P) parent;
00223 p_array<mln_psite(I)> s;
00224
00225
00226 fext = geom::translate(m, dp.to_vec(), f, d);
00227 initialize(parent, fext);
00228 s.reserve(geom::nsites(fext));
00229
00230
00231 {
00232
00233 mln_ch_value(I, bool) deja_vu;
00234 mln_ch_value(I, P) zpar;
00235
00236 initialize(deja_vu, fext);
00237 initialize(zpar, fext);
00238 mln::data::fill(deja_vu, false);
00239
00240 mln_bkd_piter(S) p_f(s_f);
00241 mln_bkd_piter(S) p_m(s_m);
00242 p_f.start();
00243 p_m.start();
00244
00245
00246 while (p_m.is_valid() || p_f.is_valid())
00247 {
00248 mln_bkd_piter(S)& it = (!p_f.is_valid() || (p_m.is_valid() && f(p_f) <= m(p_m))) ? p_m : p_f;
00249
00250 P p = it;
00251 P ext = p + dp;
00252
00253
00254
00255
00256
00257
00258
00259 mln_assertion(!(deja_vu(p) && deja_vu(ext)));
00260 if (deja_vu(ext))
00261 {
00262 mln_assertion(m(p) >= f(p));
00263
00264 parent(p) = p;
00265 zpar(p) = p;
00266
00267 P r = internal::find_root(zpar, ext);
00268 zpar(r) = p;
00269 parent(r) = p;
00270
00271 deja_vu(p) = true;
00272 }
00273 else if (deja_vu(p))
00274 {
00275 mln_assertion(f(p) > m(p));
00276 parent(p) = ext;
00277 zpar(p) = ext;
00278 parent(ext) = ext;
00279 zpar(ext) = ext;
00280
00281 mln_niter(N) n(nbh, ext);
00282 for_all(n)
00283 if (d_ext.has(n) && deja_vu(n))
00284 {
00285 P r = internal::find_root(zpar, n);
00286
00287 if (r != ext)
00288 {
00289 parent(r) = ext;
00290 zpar(r) = ext;
00291 }
00292 }
00293 deja_vu(ext) = true;
00294 }
00295 else if (f(p) <= m(p))
00296 {
00297 zpar(ext) = ext;
00298 mln_niter(N) n(nbh, ext);
00299 for_all(n)
00300 if (d_ext.has(n) && deja_vu(n))
00301 {
00302 P r = internal::find_root(zpar, n);
00303 if (r != ext)
00304 {
00305 zpar(r) = ext;
00306 parent(r) = ext;
00307 }
00308 }
00309 deja_vu(ext) = true;
00310 }
00311 else
00312 {
00313 deja_vu(p) = true;
00314 }
00315 it.next();
00316 }
00317 }
00318 std::cout << ">> MAJ zpar: " << internal::t_prop << " s" << std::endl;
00319 std::cout << "Parent construction: " << tm << " s" << std::endl;
00320 tm.restart();
00321
00322
00323 {
00324 mln_ch_value(I, bool) deja_vu(d_ext);
00325 mln::data::fill(deja_vu, false);
00326 mln_fwd_piter(S) p(s_f);
00327 for_all(p)
00328 {
00329 P q = parent(p);
00330 if (!f.domain().has(q))
00331 internal::update_m_parent(fext, parent, deja_vu, s, d_ext, p);
00332 else if (fext(parent(q)) == f(q))
00333 parent(p) = parent(q);
00334 s.append(p);
00335
00336 mln_assertion((q = parent(p)) == parent(q) || fext(q) != fext(parent(q)));
00337 }
00338 }
00339 std::cout << "Canonization: " << tm << " s" << std::endl;
00340
00341
00342
00343 tree::data<I, p_array<mln_psite(I)> > tree(fext, parent, s);
00344 trace::exiting("morpho::tree::impl::generic::dual_union_find");
00345
00346 return tree;
00347 }
00348
00349 }
00350
00351 }
00352
00353 # endif // ! MLN_INCLUDE_ONLY
00354
00355 }
00356
00357 }
00358
00359 }
00360
00361 #endif // !MLN_MORPHO_TREE_IMLP_DUAL_UNION_FIND_HH