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 #include <iostream>
00027 #include <iomanip>
00028
00029 #include <vector>
00030
00031 #include <mln/core/image/image2d.hh>
00032 #include <mln/core/routine/duplicate.hh>
00033 #include <mln/core/alias/neighb2d.hh>
00034 #include <mln/core/site_set/p_queue_fast.hh>
00035 #include <mln/labeling/blobs.hh>
00036 #include <mln/io/pgm/load.hh>
00037 #include <mln/debug/println.hh>
00038 #include <mln/draw/line.hh>
00039 #include <mln/pw/all.hh>
00040 #include <mln/binarization/threshold.hh>
00041
00042 #include <mln/value/int_u8.hh>
00043 #include <mln/value/label_8.hh>
00044 #include <mln/core/alias/point2d.hh>
00045
00046 #include <mln/core/site_set/p_faces.hh>
00047 #include <mln/core/image/complex_image.hh>
00048 #include <mln/core/alias/complex_geometry.hh>
00049 #include <mln/core/alias/complex_image.hh>
00050
00051
00052 #include <mln/core/image/complex_neighborhoods.hh>
00053 #include <mln/core/image/complex_neighborhood_piter.hh>
00054 #include <mln/core/image/complex_windows.hh>
00055
00056 #include <mln/data/fill.hh>
00057 #include <mln/morpho/gradient.hh>
00058 #include <mln/labeling/colorize.hh>
00059
00060 #include "chain.hh"
00061
00062
00063
00064
00065 namespace mln
00066 {
00067
00068 template <typename I, typename N>
00069 mln_concrete(I)
00070 influence_zones(const I& input, const N& nbh)
00071 {
00072 mln_concrete(I) output = duplicate(input);
00073
00074 p_queue_fast<mln_site(I)> q;
00075
00076 {
00077
00078 mln_piter(I) p(input.domain());
00079 mln_niter(N) n(nbh, p);
00080 for_all(p)
00081 if (input(p) != 0)
00082 for_all(n) if (input.has(n))
00083 if (input(n) == 0)
00084 {
00085 q.push(p);
00086 break;
00087 }
00088 }
00089 {
00090
00091 mln_site(I) p;
00092 mln_niter(N) n(nbh, p);
00093 while (! q.is_empty())
00094 {
00095 p = q.pop_front();
00096 mln_invariant(output(p) != 0);
00097 for_all(n) if (input.has(n))
00098 if (output(n) == 0)
00099 {
00100 output(n) = output(p);
00101 q.push(n);
00102 }
00103 }
00104 }
00105
00106 return output;
00107
00108 }
00109
00110
00111 namespace io
00112 {
00113
00114 namespace neato
00115 {
00116
00117
00118
00119 void
00120 save(const complex_image<1, discrete_plane_1complex_geometry,
00121 value::int_u8>& ima,
00122 const std::string& filename,
00123 const std::string& bgcolor = "#0000C0",
00124 const std::string& fontcolor = "#0000C0",
00125 bool empty_vertex_label = true)
00126 {
00127 typedef value::int_u8 V;
00128 typedef complex_image<1, discrete_plane_1complex_geometry, V> I;
00129 const unsigned D = 1;
00130 typedef discrete_plane_1complex_geometry G;
00131
00132 std::ofstream g(filename.c_str());
00133 g << "graph wst" << std::endl
00134 << "{" << std::endl
00135 << " graph [bgcolor = \"" << bgcolor << "\"]" << std::endl
00136 << " edge [color = \"#FFFFFF\"]" << std::endl
00137 << " node [color = \"#FFFFFF\", height=\"5\", width=\"5\","
00138 << " fontsize=\"100\", fontcolor = \"" << fontcolor << "\"]"
00139 << std::endl;
00140
00141
00142 p_n_faces_fwd_piter<D, G> v(ima.domain(), 0);
00143 typedef complex_higher_neighborhood<D, G> e_nbh_t;
00144 e_nbh_t e_nbh;
00145 for_all(v)
00146 {
00147 V vertex_color = ima(v);
00148 std::ostringstream vertex_color_str;
00149
00150 vertex_color_str << '#'
00151 << std::hex
00152 << std::setfill('0')
00153 << std::setw(2) << vertex_color
00154 << std::setw(2) << vertex_color
00155 << std::setw(2) << vertex_color
00156 << std::dec;
00157
00158 g << " v" << v.unproxy_().face_id()
00159 << " [pos = \""
00160 << std::fixed << std::setprecision(1)
00161 << (float)v.to_site().front()[1] << ", "
00162 << -(float)v.to_site().front()[0]
00163 << "\", color = \"" << vertex_color_str.str()
00164 << "\", fillcolor = \"" << vertex_color_str.str()
00165 << "\", pin = \"true\", style=\"filled,setlinewidth(3)\"";
00166 if (empty_vertex_label)
00167 g << ", label = \"\"";
00168 g << "];"
00169 << std::endl;
00170 }
00171
00172
00173 p_n_faces_fwd_piter<D, G> e(ima.domain(), 1);
00174 typedef complex_lower_neighborhood<D, G> v_nbh_t;
00175 v_nbh_t v_nbh;
00176 mln_niter_(v_nbh_t) adj_v(v_nbh, e);
00177 for_all(e)
00178 {
00179 V edge_color = ima(e);
00180 std::ostringstream edge_color_str;
00181 edge_color_str << '#'
00182 << std::hex
00183 << std::setfill('0')
00184 << std::setw(2) << edge_color
00185 << std::setw(2) << edge_color
00186 << std::setw(2) << edge_color
00187 << std::dec;
00188
00189
00190 adj_v.start();
00191 topo::face<1> v1 = adj_v.unproxy_().face();
00192 point2d p1 = adj_v.to_site().front();
00193 adj_v.next();
00194 topo::face<1> v2 = adj_v.unproxy_().face();
00195 point2d p2 = adj_v.to_site().front();
00196 adj_v.next();
00197 mln_invariant(!adj_v.is_valid());
00198
00199 g << " v" << v1.face_id() << " -- v" << v2.face_id() << " ";
00200 g << "[color = \"" << edge_color_str.str()
00201 << "\", style=\"setlinewidth(10)\"];" << std::endl;
00202 }
00203
00204 g << "}" << std::endl;
00205 g.close();
00206 }
00207
00208
00209 void
00210 save(const complex_image<1, discrete_plane_1complex_geometry,
00211 value::rgb8>& ima,
00212 const std::string& filename,
00213 const std::string& bgcolor = "#0000C0",
00214 const std::string& fontcolor = "#0000C0",
00215 bool empty_vertex_label = true)
00216 {
00217 typedef value::rgb8 V;
00218 typedef complex_image<1, discrete_plane_1complex_geometry, V> I;
00219 const unsigned D = 1;
00220 typedef discrete_plane_1complex_geometry G;
00221
00222 std::ofstream g(filename.c_str());
00223 g << "graph wst" << std::endl
00224 << "{" << std::endl
00225 << " graph [bgcolor = \"" << bgcolor << "\"]" << std::endl
00226 << " edge [color = \"#FFFFFF\"]" << std::endl
00227 << " node [color = \"#FFFFFF\", height=\"5\", width=\"5\","
00228 << " fontsize=\"100\", fontcolor = \"" << fontcolor << "\"]"
00229 << std::endl;
00230
00231
00232 p_n_faces_fwd_piter<D, G> v(ima.domain(), 0);
00233 typedef complex_higher_neighborhood<D, G> e_nbh_t;
00234 e_nbh_t e_nbh;
00235 for_all(v)
00236 {
00237 V vertex_color = ima(v);
00238 std::ostringstream vertex_color_str;
00239
00240 vertex_color_str << '#'
00241 << std::hex
00242 << std::setfill('0')
00243 << std::setw(2) << vertex_color.red()
00244 << std::setw(2) << vertex_color.green()
00245 << std::setw(2) << vertex_color.blue()
00246 << std::dec;
00247
00248 g << " v" << v.unproxy_().face_id()
00249 << " [pos = \""
00250 << std::fixed << std::setprecision(1)
00251 << (float)v.to_site().front()[1] << ", "
00252 << -(float)v.to_site().front()[0]
00253 << "\", color = \"" << vertex_color_str.str()
00254 << "\", fillcolor = \"" << vertex_color_str.str()
00255 << "\", pin = \"true\", style=\"filled,setlinewidth(3)\"";
00256 if (empty_vertex_label)
00257 g << ", label = \"\"";
00258 g << "];"
00259 << std::endl;
00260 }
00261
00262
00263 p_n_faces_fwd_piter<D, G> e(ima.domain(), 1);
00264 typedef complex_lower_neighborhood<D, G> v_nbh_t;
00265 v_nbh_t v_nbh;
00266 mln_niter_(v_nbh_t) adj_v(v_nbh, e);
00267 for_all(e)
00268 {
00269 V edge_color = ima(e);
00270 std::ostringstream edge_color_str;
00271 edge_color_str << '#'
00272 << std::hex
00273 << std::setfill('0')
00274 << std::setw(2) << edge_color.red()
00275 << std::setw(2) << edge_color.green()
00276 << std::setw(2) << edge_color.blue()
00277 << std::dec;
00278
00279
00280 adj_v.start();
00281 topo::face<1> v1 = adj_v.unproxy_().face();
00282 point2d p1 = adj_v.to_site().front();
00283 adj_v.next();
00284 topo::face<1> v2 = adj_v.unproxy_().face();
00285 point2d p2 = adj_v.to_site().front();
00286 adj_v.next();
00287 mln_invariant(!adj_v.is_valid());
00288
00289 g << " v" << v1.face_id() << " -- v" << v2.face_id() << " ";
00290 g << "[color = \"" << edge_color_str.str()
00291 << "\", style=\"setlinewidth(10)\"];" << std::endl;
00292 }
00293
00294 g << "}" << std::endl;
00295 g.close();
00296 }
00297
00298 }
00299
00300 }
00301
00302 }
00303
00304
00305
00306 mln::int_u8_1complex_image2d
00307 make_complex_image(const mln::image2d<mln::value::int_u8>& input)
00308 {
00309 using namespace mln;
00310 using mln::value::int_u8;
00311
00312
00313
00314
00315
00316 border::thickness = 0;
00317
00318 unsigned nlabels;
00319 image2d<unsigned> label =
00320 labeling::blobs(mln::binarization::threshold(input, 1), c4(), nlabels);
00321
00322 std::cout << "n seeds = " << nlabels << std::endl;
00323 {
00324 image2d<int_u8> lab(label.domain());
00325 data::paste(label, lab);
00326 }
00327
00328 image2d<unsigned> iz = influence_zones(label, c4());
00329 {
00330 image2d<int_u8> IZ(iz.domain());
00331 data::paste(iz, IZ);
00332 }
00333
00334
00335
00336
00337 std::vector< std::vector<bool> > adj(nlabels + 1);
00338 for (unsigned l = 1; l <= nlabels; ++l)
00339 adj[l].resize(nlabels + 1, false);
00340
00341 {
00342 box2d::piter p(iz.domain());
00343 for_all(p)
00344 {
00345 point2d r = p + right, b = p + down;
00346 if (iz.has(r) && iz(p) != iz(r))
00347 {
00348 if (iz(p) <= iz(r))
00349 adj[iz(p)][iz(r)] = true;
00350 else
00351 adj[iz(r)][iz(p)] = true;
00352 }
00353 if (iz.has(b) && iz(p) != iz(b))
00354 {
00355 if (iz(p) <= iz(b))
00356 adj[iz(p)][iz(b)] = true;
00357 else
00358 adj[iz(b)][iz(p)] = true;
00359 }
00360 }
00361 }
00362
00363
00364
00365
00366 const unsigned D = 1;
00367
00368 topo::complex<D> c;
00369
00370 typedef point2d P;
00371 typedef geom::complex_geometry<D, P> G;
00372 G geom;
00373
00374
00375 typedef topo::n_face<0, D> vertex;
00376 typedef topo::n_face<1, D> edge;
00377
00378 {
00379
00380
00381 std::vector<vertex> v;
00382 {
00383 box2d::piter p(label.domain());
00384 for_all(p)
00385 if (label(p) != 0)
00386 {
00387 geom.add_location(p);
00388 v.push_back(c.add_face());
00389 }
00390 }
00391
00392 std::cout << "v size = " << v.size() << std::endl;
00393
00394
00395 std::vector<edge> e;
00396 {
00397 for (unsigned l = 1; l <= nlabels; ++l)
00398 for (unsigned ll = l + 1; ll <= nlabels; ++ll)
00399 if (adj[l][ll])
00400 e.push_back( c.add_face(-v[l-1] + v[ll-1]) );
00401 }
00402
00403 std::cout << "e size = " << e.size() << std::endl;
00404
00405 }
00406
00407
00408
00409
00410
00411
00412 p_complex<D, G> pc(c, geom);
00413
00414
00415
00416
00417
00418
00419
00420 typedef complex_image<D, G, int_u8> output_t;
00421
00422
00423 output_t output(pc);
00424
00425
00426 p_n_faces_fwd_piter<D, G> v(output.domain(), 0);
00427 for_all(v)
00428 output(v) = input(v.to_site().front());
00429
00430
00431 p_n_faces_fwd_piter<D, G> e(output.domain(), 1);
00432 for_all(e)
00433 output(e) = 128;
00434
00435 return output;
00436 }
00437
00438
00439
00440 template <unsigned D, typename G, typename V>
00441 mln::complex_higher_dim_connected_n_face_window_p<D, G>
00442 make_elt_win(const mln::complex_image<D, G, V>& )
00443 {
00444 return mln::complex_higher_dim_connected_n_face_window_p<D, G>();
00445 }
00446
00447
00448 template <unsigned D, typename G, typename V>
00449 mln::complex_higher_dim_connected_n_face_neighborhood<D, G>
00450 make_elt_nbh(const mln::complex_image<D, G, V>& )
00451 {
00452 return mln::complex_higher_dim_connected_n_face_neighborhood<D, G>();
00453 }
00454
00455
00456 int main(int argc, char* argv[])
00457 {
00458 if (argc != 4)
00459 {
00460 std::cerr << "usage: " << argv[0] << " seeds.pgm lambda output.neato"
00461 << std::endl;
00462 std::exit(1);
00463 }
00464 std::string input_filename = argv[1];
00465 unsigned lambda = atoi(argv[2]);
00466 std::string output_filename = argv[3];
00467
00468 using namespace mln;
00469 using mln::value::int_u8;
00470
00471 typedef int_u8_1complex_image2d int_u8_graph_image2d;
00472
00473 typedef int_u8_graph_image2d input;
00474 typedef value::label_8 label;
00475 typedef mln_ch_value_(input, label) output;
00476 label nbasins;
00477
00478
00479 image2d<int_u8> seeds = io::pgm::load<int_u8>(input_filename);
00480
00481 typedef int_u8_graph_image2d ima_t;
00482 ima_t ima = make_complex_image(seeds);
00483 io::neato::save(ima, "apps/graph.neato");
00484
00485
00486 input g = morpho::gradient(ima, make_elt_win(ima));
00487
00488 #if 0
00489
00490 io::neato::save(g, "apps/graph-g.neato");
00491 #endif
00492
00493
00494 output s = chain(g, make_elt_nbh(g), lambda, nbasins);
00495 io::neato::save(labeling::colorize(value::rgb8(), s, nbasins),
00496 output_filename);
00497 }