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_RD_HH
00027 # define MLN_MORPHO_RD_HH
00028
00029
00030
00031
00035
00036 # include <vector>
00037 # include <mln/core/concept/image.hh>
00038 # include <mln/core/concept/neighborhood.hh>
00039 # include <mln/trait/value_.hh>
00040 # include <mln/data/fill.hh>
00041 # include <mln/data/compare.hh>
00042 # include <mln/util/ord.hh>
00043
00044
00045 namespace mln
00046 {
00047
00048 namespace morpho
00049 {
00050
00051 template <typename I, typename N>
00052 I
00053 Rd(const Image<I>& f, const Image<I>& g, const Neighborhood<N>& nbh);
00054
00055
00056 # ifndef MLN_INCLUDE_ONLY
00057
00058 namespace impl
00059 {
00060
00061 template <typename I>
00062 inline
00063 std::vector<unsigned> compute_histo(const I& ima)
00064 {
00065 std::vector<unsigned> h(256, 0);
00066 mln_piter(I) p(ima.domain());
00067 for_all(p)
00068 ++h[ima(p)];
00069 return h;
00070 }
00071
00072 template <typename I>
00073 inline
00074 std::vector<mln_psite(I)> histo_reverse_sort(const I& ima)
00075 {
00076 std::vector<unsigned> h = compute_histo(ima);
00077
00078 std::vector<int> loc(256);
00079 loc[255] = 0;
00080 for (int l = 254; l >= 0; --l)
00081 loc[l] = loc[l+1] + h[l+1];
00082 std::vector<mln_psite(I)> vec(ima.domain().nsites());
00083
00084 mln_piter(I) p(ima.domain());
00085 for_all(p)
00086 vec[loc[ima(p)]++] = p;
00087 return vec;
00088 }
00089
00090
00091 template <typename I, typename N>
00092 struct Rd
00093 {
00094 typedef mln_psite(I) point;
00095 typedef mln_value(I) value;
00096
00097
00098 const I& f;
00099 const I& g;
00100 const N& nbh;
00101
00102
00103 I o;
00104
00105
00106 mln_ch_value(I, bool) is_proc;
00107 mln_ch_value(I, point) parent;
00108 std::vector<point> S;
00109
00110 inline
00111 Rd(const I& f, const I& g, const N& nbh)
00112 : f(f), g(g), nbh(nbh),
00113 o(f.domain()),
00114 is_proc(f.domain()),
00115 parent(f.domain())
00116 {
00117
00118 data::fill(o, f);
00119 S = histo_reverse_sort(g);
00120 data::fill(is_proc, false);
00121
00122
00123 for (unsigned i = 0; i < S.size(); ++i)
00124 {
00125 point p = S[i];
00126
00127 make_set(p);
00128 mln_niter(N) n(nbh, p);
00129 for_all(n)
00130 {
00131 if (f.domain().has(n))
00132 assert(is_proc(n) == is_proc__(n, p));
00133 if (f.has(n) && is_proc(n))
00134 do_union(n, p);
00135 }
00136 is_proc(p) = true;
00137 }
00138
00139
00140 for (int i = S.size() - 1; i >= 0; --i)
00141 {
00142 point p = S[i];
00143 if (parent(p) == p)
00144 {
00145 if (o(p) == mln_max(value))
00146 o(p) = g(p);
00147 else
00148 o(p) = o(parent(p));
00149 }
00150 }
00151
00152 }
00153
00154 inline
00155 bool is_proc__(const point& n, const point& p) const
00156 {
00157 return g(n) > g(p) || (g(n) == g(p) && util::ord_strict(point(n), p));
00158 }
00159
00160 inline
00161 void make_set(const point& p)
00162 {
00163 parent(p) = p;
00164 }
00165
00166 inline
00167 point find_root(const point& x)
00168 {
00169 if (parent(x) == x)
00170 return x;
00171 else
00172 return parent(x) = find_root(parent(x));
00173 }
00174
00175 inline
00176 bool equiv(const point& r, const point& p)
00177 {
00178 return g(r) == g(p) || g(p) >= o(r);
00179 }
00180
00181 inline
00182 void do_union(const point& n, const point& p)
00183 {
00184 point r = find_root(n);
00185 if (r != p)
00186 {
00187 if (equiv(r, p))
00188 {
00189 parent(r) = p;
00190 if (o(r) > o(p))
00191 o(p) = o(r);
00192 }
00193 else
00194 o(p) = mln_max(value);
00195 }
00196 }
00197
00198 };
00199
00200 }
00201
00202
00203
00204
00205 template <typename I, typename N>
00206 inline
00207 I
00208 Rd(const Image<I>& f, const Image<I>& g, const Neighborhood<N>& nbh)
00209 {
00210 assert(f <= g);
00211 impl::Rd<I,N> run(exact(f), exact(g), exact(nbh));
00212 return run.o;
00213 }
00214
00215 # endif // ! MLN_INCLUDE_ONLY
00216
00217 }
00218
00219 }
00220
00221
00222 #endif // ! MLN_MORPHO_RD_HH