Milena (Olena)
User documentation 2.0a Id
|
00001 // Copyright (C) 2008, 2009 EPITA Research and Development Laboratory (LRDE) 00002 // 00003 // This file is part of Olena. 00004 // 00005 // Olena is free software: you can redistribute it and/or modify it under 00006 // the terms of the GNU General Public License as published by the Free 00007 // Software Foundation, version 2 of the License. 00008 // 00009 // Olena is distributed in the hope that it will be useful, 00010 // but WITHOUT ANY WARRANTY; without even the implied warranty of 00011 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 00012 // General Public License for more details. 00013 // 00014 // You should have received a copy of the GNU General Public License 00015 // along with Olena. If not, see <http://www.gnu.org/licenses/>. 00016 // 00017 // As a special exception, you may use this file as part of a free 00018 // software project without restriction. Specifically, if other files 00019 // instantiate templates or use macros or inline functions from this 00020 // file, or you compile this file and link it with other files to produce 00021 // an executable, this file does not by itself cause the resulting 00022 // executable to be covered by the GNU General Public License. This 00023 // exception does not however invalidate any other reasons why the 00024 // executable file might be covered by the GNU General Public License. 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 // Doesn't C++ have a better way to express Pi? 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 | Complex image. | 00073 `----------------*/ 00074 00075 // Image type. 00076 typedef mln::float_2complex_image3df ima_t; 00077 // Dimension of the image (and therefore of the complex). 00078 static const unsigned D = ima_t::dim; 00079 // Geometry of the image. 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 // Compute the pseudo_inverse curvature at each vertex. 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 // Pseudo-inverse curvature. 00093 float h_inv = 1 / pi * (atan(-h) + pi / 2); 00094 input(v) = h_inv; 00095 // FIXME: The program should allow the user to choose the kind 00096 // of measure. 00097 // input(v) = mln::math::max(mln::math::sqr(curv.first(v)), 00098 // mln::math::sqr(curv.second(v))); 00099 } 00100 00101 // Values on edges. 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 // Iterate on edges (1-faces). 00107 for_all(e) 00108 { 00109 float s = 0.0f; 00110 unsigned n = 0; 00111 // Iterate on vertices (0-faces). 00112 for_all(adj_v) 00113 { 00114 s += input(adj_v); 00115 ++n; 00116 } 00117 input(e) = s / n; 00118 // An edge should be adjacent to exactly two vertices. 00119 mln_invariant(n <= 2); 00120 } 00121 00122 /*-----------------. 00123 | Simplification. | 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 | WST. | 00136 `------*/ 00137 00138 /* FIXME: Note that the WST is doing too much work, since we have 00139 not constrained the domain of the image to 1-faces only. 00140 It would be great if we could use something like this: 00141 00142 closed_input | mln::p_faces<1, D, G>(closed_input.domain()) 00143 00144 as input of the WST. */ 00145 00146 // Compute the WST on edges. 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 // Label polygons (i.e., propagate labels from edges to polygons). 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 | Output. | 00165 `---------*/ 00166 00167 mln::rgb8_2complex_image3df output(wshed.domain()); 00168 mln::data::fill(output, mln::literal::white); 00169 00170 // FIXME: Use a colorize functor instead. 00171 // Choose random colors for each basin number. 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 }