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/graph_image.hh>
00045 #include <mln/core/image/graph_elt_neighborhood.hh>
00046
00047 #include <mln/morpho/closing/area.hh>
00048 #include <mln/labeling/regional_minima.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 std::vector<float> vertex_h_inv(mesh.faces.size(), 0.f);
00088 for (unsigned v = 0; v < mesh.vertices.size(); ++v)
00089 {
00090 float h = (mesh.curv1[v] + mesh.curv2[v]) / 2;
00091 float h_inv = 1 / pi * atan(-h) + pi / 2;
00092 vertex_h_inv[v] = h_inv;
00093 }
00094
00095
00096
00097
00098 typedef int curv_t;
00099 std::vector<curv_t> face_h_inv(mesh.faces.size(), 0.f);
00100 for (unsigned f = 0; f < mesh.faces.size(); ++f)
00101 {
00102 float h_inv =
00103 (vertex_h_inv[mesh.faces[f][0]] +
00104 vertex_h_inv[mesh.faces[f][1]] +
00105 vertex_h_inv[mesh.faces[f][2]])
00106 / 3;
00107
00108
00109
00110
00111 face_h_inv[f] = 1000 * h_inv;
00112 }
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124 mln::util::graph<mln::point3d> g;
00125
00126 for (unsigned i = 0; i < mesh.faces.size(); ++i)
00127 g.add_vertex (mln::point3d(i, i, i));
00128
00129
00130 mesh.need_across_edge();
00131 for (unsigned f = 0; f < mesh.faces.size(); ++f)
00132 for (unsigned e = 0; e < 3; ++e)
00133 {
00134 int f_adj = mesh.across_edge[f][e];
00135 if (f_adj != -1)
00136
00137 g.add_edge(f, f_adj);
00138 }
00139
00140
00141
00142
00143
00144 mln::p_graph<mln::point3d> pg(g);
00145
00146 typedef mln::graph_image<mln::point3d, curv_t> ima_t;
00147 ima_t g_ima(pg, face_h_inv);
00148
00149
00150
00151
00152
00153 typedef mln::graph_elt_neighborhood<mln::point3d> nbh_t;
00154 nbh_t nbh;
00155
00156 ima_t closed_g_ima = mln::morpho::closing::area(g_ima, nbh, lambda);
00157
00158
00159
00160
00161
00162 typedef unsigned label_t;
00163 label_t nlabels;
00164 typedef mln::graph_image<mln::point3d, label_t> label_ima_t;
00165 label_ima_t minima =
00166 mln::labeling::regional_minima(closed_g_ima, nbh, nlabels);
00167 std::cout << "nlabels = " << nlabels << std::endl;
00168
00169
00170
00171
00172
00173
00174
00175
00176
00177
00178
00179
00180
00181
00182
00183 std::vector<bool> face_value (minima.domain().nvertices(), true);
00184 mln_piter_(label_ima_t) pm(minima.domain());
00185 for_all(pm)
00186
00187 if (minima(pm) != 0)
00188 {
00189
00190
00191 mln_psite_(label_ima_t) pp(pm);
00192 face_value[pp.id().to_equiv()] = false;
00193 }
00194
00195
00196 FILE* f_out = fopen(output_filename.c_str(), "wb");
00197 if (!f_out)
00198 {
00199 std::cerr << "Error opening " << output_filename.c_str()
00200 << " for writing." << std::endl;
00201 std::exit(2);
00202 }
00203 write_off_binary(mesh_ptr, face_value, f_out);
00204 fclose(f_out);
00205
00206 delete mesh_ptr;
00207 }