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
00031 #include <cstdlib>
00032 #include <cmath>
00033
00034 #include <utility>
00035 #include <iostream>
00036
00037 #include <mln/core/image/complex_image.hh>
00038 #include <mln/core/image/complex_neighborhoods.hh>
00039
00040 #include <mln/morpho/closing/area.hh>
00041 #include <mln/morpho/meyer_wst.hh>
00042
00043 #include <mln/math/max.hh>
00044 #include <mln/math/sqr.hh>
00045
00046 #include <mln/literal/white.hh>
00047
00048 #include <mln/io/off/load.hh>
00049 #include <mln/io/off/save.hh>
00050
00051 #include "trimesh/misc.hh"
00052
00053
00054
00055 static const float pi = 4 * atanf(1);
00056
00057
00058 int main(int argc, char* argv[])
00059 {
00060 if (argc != 4)
00061 {
00062 std::cerr << "usage: " << argv[0] << " input.off lambda output.off"
00063 << std::endl;
00064 std::exit(1);
00065 }
00066
00067 std::string input_filename = argv[1];
00068 unsigned lambda = atoi(argv[2]);
00069 std::string output_filename = argv[3];
00070
00071
00072
00073
00074
00075
00076 typedef mln::float_2complex_image3df ima_t;
00077
00078 static const unsigned D = ima_t::dim;
00079
00080 typedef mln_geom_(ima_t) G;
00081
00082 mln::bin_2complex_image3df bin_input;
00083 mln::io::off::load(bin_input, input_filename);
00084 std::pair<ima_t, ima_t> curv = mln::geom::mesh_curvature(bin_input.domain());
00085
00086
00087 ima_t input(bin_input.domain());
00088 mln::p_n_faces_fwd_piter<D, G> v(input.domain(), 0);
00089 for_all(v)
00090 {
00091 float h = (curv.first(v) + curv.second(v)) / 2;
00092
00093 float h_inv = 1 / pi * (atan(-h) + pi / 2);
00094 input(v) = h_inv;
00095
00096
00097
00098
00099 }
00100
00101
00102 mln::p_n_faces_fwd_piter<D, G> e(input.domain(), 1);
00103 typedef mln::complex_lower_neighborhood<D, G> adj_vertices_nbh_t;
00104 adj_vertices_nbh_t adj_vertices_nbh;
00105 mln_niter_(adj_vertices_nbh_t) adj_v(adj_vertices_nbh, e);
00106
00107 for_all(e)
00108 {
00109 float s = 0.0f;
00110 unsigned n = 0;
00111
00112 for_all(adj_v)
00113 {
00114 s += input(adj_v);
00115 ++n;
00116 }
00117 input(e) = s / n;
00118
00119 mln_invariant(n <= 2);
00120 }
00121
00122
00123
00124
00125
00127 typedef
00128 mln::complex_higher_dim_connected_n_face_neighborhood<D, G>
00129 adj_edges_nbh_t;
00130 adj_edges_nbh_t adj_edges_nbh;
00131
00132 ima_t closed_input = mln::morpho::closing::area(input, adj_edges_nbh, lambda);
00133
00134
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146
00147 typedef unsigned wst_val_t;
00148 wst_val_t nbasins;
00149 typedef mln::unsigned_2complex_image3df wst_ima_t;
00150 wst_ima_t wshed =
00151 mln::morpho::meyer_wst(closed_input, adj_edges_nbh, nbasins);
00152 std::cout << "nbasins = " << nbasins << std::endl;
00153
00154
00155 typedef mln::complex_higher_neighborhood<D, G> adj_polygons_nbh_t;
00156 adj_polygons_nbh_t adj_polygons_nbh;
00157 mln_niter_(adj_polygons_nbh_t) adj_p(adj_polygons_nbh, e);
00158 for_all(e)
00159 if (wshed(e) != 0)
00160 for_all(adj_p)
00161 wshed(adj_p) = wshed(e);
00162
00163
00164
00165
00166
00167 mln::rgb8_2complex_image3df output(wshed.domain());
00168 mln::data::fill(output, mln::literal::white);
00169
00170
00171
00172 std::vector<mln::value::rgb8> basin_color (nbasins + 1);
00173 for (unsigned i = 0; i <= nbasins; ++i)
00174 basin_color[i] = mln::value::rgb8(random() % 256,
00175 random() % 256,
00176 random() % 256);
00177 mln_piter_(ima_t) f(wshed.domain());
00178 for_all(f)
00179 output(f) = basin_color[wshed(f)];
00180
00181 mln::io::off::save(output, output_filename);
00182 }