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 00029 00030 // FIXME: Factor commons parts between mesh-segm and mesh-skel. 00031 00032 #include <cstdlib> 00033 #include <cmath> 00034 00035 #include <utility> 00036 #include <iostream> 00037 00038 #include <TriMesh.h> 00039 00040 #include <mln/core/alias/point3d.hh> 00041 #include <mln/core/alias/point3d.hh> 00042 00043 #include <mln/util/graph.hh> 00044 #include <mln/core/image/line_graph_image.hh> 00045 #include <mln/core/image/line_graph_elt_neighborhood.hh> 00046 00047 #include <mln/morpho/closing/area.hh> 00048 #include <mln/morpho/meyer_wst.hh> 00049 00050 #include "io.hh" 00051 00052 00053 // Doesn't C++ have a better way to express Pi ? 00054 const float pi = 4 * atanf(1); 00055 00056 00057 int main(int argc, char* argv[]) 00058 { 00059 if (argc != 4) 00060 { 00061 std::cerr << "usage: " << argv[0] << " input.off lambda output.off" 00062 << std::endl; 00063 std::exit(1); 00064 } 00065 00066 std::string input_filename = argv[1]; 00067 unsigned lambda = atoi(argv[2]); 00068 std::string output_filename = argv[3]; 00069 00070 00071 /*-------. 00072 | Mesh. | 00073 `-------*/ 00074 00075 // TriMesh is a pain: it systematically allocates on the heap. 00076 // Introduce another name to manipulate the mesh as a (non-pointer) 00077 // object. 00078 TriMesh* mesh_ptr = TriMesh::read(input_filename.c_str()); 00079 if (!mesh_ptr) 00080 std::exit(2); 00081 TriMesh& mesh = *mesh_ptr; 00082 00083 // Computes faces (triangles). 00084 mesh.need_faces(); 00085 // Computation of the mean curvature on each vertex of the mesh. 00086 mesh.need_curvatures(); 00087 /* FIXME: Our implementation of the WST doesn't work well with 00088 floats (yet). Convert floating point values to a proportional 00089 integer value for the moment. */ 00090 typedef int curv_t; 00091 std::vector<curv_t> vertex_h_inv(mesh.vertices.size(), 0.f); 00092 for (unsigned v = 0; v < mesh.vertices.size(); ++v) 00093 { 00094 float h = (mesh.curv1[v] + mesh.curv2[v]) / 2; 00095 float h_inv = 1 / pi * atan(-h) + pi / 2; 00096 /* FIXME: This coefficient is used to distinguish small 00097 curvature values. We should get rid of it as soon as 00098 morpho::meyer_wst works correctly on images of float 00099 values. */ 00100 vertex_h_inv[v] = 1000 * h_inv; 00101 } 00102 00103 /*--------. 00104 | Graph. | 00105 `--------*/ 00106 00107 /* Build a graph whose vertices correspond to the faces of the mesh, 00108 whose edges (between two vertices) correspond to edges (between 00109 two faces) of the mesh. */ 00110 00111 // Values associated to the graph sites. 00112 // (Reminder: currently, values associated to vertices in a 00113 // line_graph_image are dummy.) 00114 std::vector<curv_t> vertex_values (mesh.faces.size(), 0.f); 00115 std::vector<curv_t> edge_values; 00116 00117 /* FIXME: We don't have the required site type yet. Simulate with a 00118 dummy type (point3d). */ 00119 mln::util::graph<mln::point3d> g; 00120 // Populate the graph with vertices. 00121 for (unsigned i = 0; i < mesh.faces.size(); ++i) 00122 g.add_vertex (mln::point3d(i, i, i)); 00123 00124 // Populate the graph with edges. 00125 mesh.need_across_edge(); 00126 for (unsigned f = 0; f < mesh.faces.size(); ++f) 00127 for (unsigned e = 0; e < 3; ++e) 00128 { 00129 int f_adj = mesh.across_edge[f][e]; 00130 if (f_adj != -1) 00131 { 00132 // Add an edge into the graph. 00133 if (g.add_edge(f, f_adj) != mln_max(mln::util::edge_id::equiv)) 00134 { 00135 // Find the edge (i.e., the two vertices) common to faces 00136 // F and F_ADJ. 00137 /* FIXME: We lack a proper interface from the TriMesh 00138 structure to do this elegantly. */ 00139 std::vector<int> adj_vertices; 00140 adj_vertices.reserve(2); 00141 for (unsigned i = 0; i < 3; ++i) 00142 for (unsigned j = 0; j < 3; ++j) 00143 if (mesh.faces[f][i] == mesh.faces[f_adj][j]) 00144 adj_vertices.push_back(mesh.faces[f][i]); 00145 mln_assertion(adj_vertices.size() == 2); 00146 00147 // Compute the mean curvature on the edge. 00148 edge_values.push_back((vertex_h_inv[adj_vertices[0]] + 00149 vertex_h_inv[adj_vertices[1]]) 00150 / 2); 00151 } 00152 00153 // Check the consistency of the two arrays. 00154 mln_assertion(g.edges().size() == edge_values.size()); 00155 } 00156 } 00157 00158 /*-------------------. 00159 | Line graph image. | 00160 `-------------------*/ 00161 00162 mln::p_line_graph<mln::point3d> plg(g); 00163 00164 typedef mln::line_graph_image<mln::point3d, curv_t> ima_t; 00165 ima_t lg_ima(plg, vertex_values, edge_values); 00166 00167 /*-----------------. 00168 | Simplification. | 00169 `-----------------*/ 00170 00171 typedef mln::line_graph_elt_neighborhood<mln::point3d> nbh_t; 00172 nbh_t nbh; 00173 00174 ima_t closed_lg_ima = mln::morpho::closing::area(lg_ima, nbh, lambda); 00175 00176 /*------. 00177 | WST. | 00178 `------*/ 00179 00180 typedef unsigned wst_val_t; 00181 wst_val_t nbasins; 00182 typedef mln::line_graph_image<mln::point3d, wst_val_t> wst_ima_t; 00183 wst_ima_t wshed = mln::morpho::meyer_wst(closed_lg_ima, nbh, nbasins); 00184 std::cout << "nbasins = " << nbasins << std::endl; 00185 00186 /*------------------------------------------. 00187 | Label graph vertices (i.e., mesh faces). | 00188 `------------------------------------------*/ 00189 00190 /* FIXME: We should be using wshed.vertex_values_ if 00191 mln::line_graph_image were fully functional... */ 00192 std::vector<unsigned> vertex_label(wshed.domain().nvertices(), 0); 00193 mln_piter_(wst_ima_t) pw(wshed.domain()); 00194 for_all(pw) 00195 if (wshed(pw) != 0) 00196 { 00197 mln_psite_(wst_ima_t) pp(pw); 00198 vertex_label[pp.first_id().to_equiv()] = wshed(pw); 00199 vertex_label[pp.second_id().to_equiv()] = wshed(pw); 00200 } 00201 00202 /*---------. 00203 | Output. | 00204 `---------*/ 00205 00206 // Choose random colors for each basin number. 00207 std::vector<mln::value::rgb8> basin_color (nbasins + 1); 00208 for (unsigned i = 0; i <= nbasins; ++i) 00209 basin_color[i] = mln::value::rgb8(random() % 256, 00210 random() % 256, 00211 random() % 256); 00212 00213 // Assign colors to graph vertices (mesh faces). 00214 std::vector<mln::value::rgb8> face_color (vertex_label.size()); 00215 for (unsigned i = 0; i < vertex_label.size() ; ++i) 00216 face_color[i] = basin_color[vertex_label[i]]; 00217 00218 // Taken and adapted from TriMesh_io.cc 00219 FILE* f_out = fopen(output_filename.c_str(), "wb"); 00220 if (!f_out) 00221 { 00222 std::cerr << "Error opening " << output_filename.c_str() 00223 << " for writing." << std::endl; 00224 std::exit(2); 00225 } 00226 write_off_colored(mesh_ptr, face_color, f_out); 00227 fclose(f_out); 00228 00229 delete mesh_ptr; 00230 }