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
00029
00030
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
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
00073
00074
00075
00076
00077
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
00084 mesh.need_faces();
00085
00086 mesh.need_curvatures();
00087
00088
00089
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
00097
00098
00099
00100 vertex_h_inv[v] = 1000 * h_inv;
00101 }
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111
00112
00113
00114 std::vector<curv_t> vertex_values (mesh.faces.size(), 0.f);
00115 std::vector<curv_t> edge_values;
00116
00117
00118
00119 mln::util::graph<mln::point3d> g;
00120
00121 for (unsigned i = 0; i < mesh.faces.size(); ++i)
00122 g.add_vertex (mln::point3d(i, i, i));
00123
00124
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
00133 if (g.add_edge(f, f_adj) != mln_max(mln::util::edge_id::equiv))
00134 {
00135
00136
00137
00138
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
00148 edge_values.push_back((vertex_h_inv[adj_vertices[0]] +
00149 vertex_h_inv[adj_vertices[1]])
00150 / 2);
00151 }
00152
00153
00154 mln_assertion(g.edges().size() == edge_values.size());
00155 }
00156 }
00157
00158
00159
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
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
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
00188
00189
00190
00191
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
00204
00205
00206
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
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
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 }