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
00030
00047 #include <cstdlib>
00048
00049 #include <vector>
00050
00051 #include <mln/value/int_u8.hh>
00052 #include <mln/value/int_u16.hh>
00053
00054 #include <mln/core/routine/duplicate.hh>
00055
00056 #include <mln/core/image/image2d.hh>
00057 #include <mln/core/var.hh>
00058
00059 #include <mln/morpho/line_gradient.hh>
00060 #include <mln/morpho/closing/area_on_vertices.hh>
00061 #include <mln/labeling/regional_minima.hh>
00062 #include <mln/morpho/meyer_wst.hh>
00063
00064 #include <mln/io/pgm/load.hh>
00065 #include <mln/io/pgm/save.hh>
00066
00067 int main(int argc, char* argv[])
00068 {
00069 if (argc != 4)
00070 {
00071 std::cerr
00072 << "usage: " << argv[0] << " max_nregions input.pgm output.pgm"
00073 << std::endl;
00074 std::exit(EXIT_FAILURE);
00075 }
00076
00077
00078
00079
00080
00081 using namespace mln;
00082 using value::int_u8;
00083 using value::int_u16;
00084
00085 typedef int_u8 val_t;
00086 typedef image2d<val_t> orig_ima_t;
00087
00088 orig_ima_t input;
00089 io::pgm::load(input, argv[2]);
00090 if (!input.is_valid())
00091 {
00092 std::cerr << "Error reading input " << argv[2] << std::endl;
00093 std::exit(2);
00094 }
00095
00096
00097
00098
00099
00100
00101 typedef util::site_pair<point2d> P;
00102 typedef fun::i2v::array<P> fedge_site_t;
00103
00104
00105 typedef fun::i2v::array<val_t> fval_t;
00106 fval_t values;
00107 typedef edge_image<P,val_t> lg_ima_t;
00108 lg_ima_t lg_ima = morpho::line_gradient(input);
00109
00110
00111
00112
00113
00114
00115
00116
00117 typedef lg_ima_t::nbh_t nbh_t;
00118 nbh_t nbh;
00119
00120 unsigned area = 0;
00121 unsigned max_area = input.nsites();
00122 unsigned nregions = mln_max(unsigned);
00123 unsigned max_nregions = atoi(argv[1]);
00124
00125 lg_ima_t result = duplicate(lg_ima);
00126 while (area < max_area && nregions > max_nregions)
00127 {
00128 ++area;
00129 std::cerr << "area = " << area << " \t"
00130 << "nregions = " << nregions << std::endl;
00131 lg_ima_t work = duplicate(result);
00132
00133 result = morpho::closing::area_on_vertices(work, nbh, area);
00134
00135
00136 labeling::regional_minima(result, nbh, nregions);
00137 }
00138
00139
00140
00141
00142
00143
00144 typedef int_u16 wst_val_t;
00145 wst_val_t nbasins;
00146 typedef edge_image<P,wst_val_t> wshed_t;
00147
00148 wshed_t wshed = morpho::meyer_wst(result, nbh, nbasins);
00149 std::cout << "nbasins = " << nbasins << std::endl;
00150
00151
00152
00153
00154
00155 const wst_val_t wshed_label = 0;
00156
00157
00158 image2d<wst_val_t> wshed2d(input.domain());
00159
00160
00161
00162
00163
00164
00165 mln_piter_(wshed_t) p(wshed.domain());
00166 for_all(p)
00167 if (wshed(p) != wshed_label)
00168 {
00169 wshed2d(p.first()) = wshed(p);
00170 wshed2d(p.second()) = wshed(p);
00171 }
00172
00173
00174 std::vector<mln_sum_(val_t)> sum(nbasins + 1, 0);
00175 std::vector<unsigned> nsites(nbasins + 1, 0);
00176 mln_piter_(orig_ima_t) q(input.domain());
00177 for_all(q)
00178 {
00179 sum[wshed2d(q)] += input(q);
00180 ++nsites[wshed2d(q)];
00181 }
00182 std::vector<float> average(nbasins + 1);
00183 for (unsigned i = 1; i <= nbasins; ++i)
00184 average[i] = float (sum[i]) / float (nsites[i]);
00185
00186
00187 orig_ima_t output(input.domain());
00188 for_all(q)
00189 output(q) = convert::to<mln_value_(orig_ima_t)>(average[wshed2d(q)]);
00190
00191 std::cout << "area = " << area << " \t"
00192 << "nregions = " << nregions << std::endl;
00193 io::pgm::save(output, argv[3]);
00194 }