[Olena-patches] 2055: Start a 3-D mesh surface skeletonization tool.

https://svn.lrde.epita.fr/svn/oln/trunk/milena Index: ChangeLog from Roland Levillain &lt;roland(a)lrde.epita.fr&gt; Start a 3-D mesh surface skeletonization tool. * apps/statues/mesh-skel.cc: New. * apps/statues/io.hh: Update header guards. (write_faces_asc_binary, write_off_binary): New functions. * apps/statues/Makefile.am (bin_PROGRAMS): Add mesh-skel. (mesh_skel_SOURCES): New. Add mesh-skel.cc and io.hh. Makefile.am | 3 io.hh | 48 ++++++++++++- mesh-skel.cc | 211 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 258 insertions(+), 4 deletions(-) Index: apps/statues/mesh-skel.cc --- apps/statues/mesh-skel.cc (revision 0) +++ apps/statues/mesh-skel.cc (revision 0) @@ -0,0 +1,211 @@ +// Copyright (C) 2008 EPITA Research and Development Laboratory (LRDE) +// +// This file is part of the Olena Library. This library is free +// software; you can redistribute it and/or modify it under the terms +// of the GNU General Public License version 2 as published by the +// Free Software Foundation. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this library; see the file COPYING. If not, write to +// the Free Software Foundation, 51 Franklin Street, Fifth Floor, +// Boston, MA 02111-1307, USA. +// +// As a special exception, you may use this file as part of a free +// software library without restriction. Specifically, if other files +// instantiate templates or use macros or inline functions from this +// file, or you compile this file and link it with other files to +// produce an executable, this file does not by itself cause the +// resulting executable to be covered by the GNU General Public +// License. This exception does not however invalidate any other +// reasons why the executable file might be covered by the GNU General +// Public License. + +/// \file apps/statues/mesh-segm.cc +/// \brief A program computing a skeleton of the surface of the +/// (triangle) mesh of a statue. + +// FIXME: Factor commons parts between mesh-segm and mesh-skel. + +#include <cstdlib> +#include <cmath> + +#include <utility> +#include <iostream> + +#include <TriMesh.h> + +#include <mln/value/rgb8.hh> + +#include <mln/core/point3d.hh> +#include <mln/make/point3d.hh> + +#include <mln/util/graph.hh> +#include <mln/core/graph_image.hh> +#include <mln/core/graph_elt_neighborhood.hh> + +#include <mln/morpho/closing_area.hh> +#include <mln/labeling/regional_minima.hh> + +#include "io.hh" + + +// Doesn't C++ have a better way to express Pi ? +const float pi = 4 * atanf(1); + + +int main(int argc, char* argv[]) +{ + if (argc != 4) + { + std::cerr << "usage: " << argv[1] << " input.off lambda output.off"; + exit(1); + } + + std::string input_filename = argv[1]; + unsigned lambda = atoi(argv[2]); + std::string output_filename = argv[3]; + + + /*-------. + | Mesh. | + `-------*/ + + // TriMesh is a pain: it systematically allocates on the heap. + // Introduce another name to manipulate the mesh as a (non-pointer) + // object. + TriMesh* mesh_ptr = TriMesh::read(input_filename.c_str()); + if (!mesh_ptr) + exit(2); + TriMesh& mesh = *mesh_ptr; + + // Computes faces (triangles). + mesh.need_faces(); + // Computation of the mean curvature on each vertex of the mesh. + mesh.need_curvatures(); + std::vector<float> vertex_h_inv(mesh.faces.size(), 0.f); + for (unsigned v = 0; v < mesh.vertices.size(); ++v) + { + float h = (mesh.curv1[v] + mesh.curv2[v]) / 2; + float h_inv = 1 / pi * atan(-h) + pi / 2; + vertex_h_inv[v] = h_inv; + } + // Computation of the mean curvature on each face of the mesh. + /* FIXME: Our implementation of the computation of the regional + minima doesn't work well with floats (yet). Convert floating + point values to a proportional integer value for the moment. */ + typedef int curv_t; + std::vector<curv_t> face_h_inv(mesh.faces.size(), 0.f); + for (unsigned f = 0; f < mesh.faces.size(); ++f) + { + float h_inv = + (vertex_h_inv[mesh.faces[f][0]] + + vertex_h_inv[mesh.faces[f][1]] + + vertex_h_inv[mesh.faces[f][2]]) + / 3; + /* FIXME: This coefficient is used to distinguish small + curvature values. We should get rid of it as soon as + labeling::regional_minima works correctly on images of float + values. */ + face_h_inv[f] = 1000 * h_inv; + } + + /*--------. + | Graph. | + `--------*/ + + /* Build a graph whose vertices correspond to the faces of the mesh, + whose edges (between two vertices) correspond to edges (between + two faces) of the mesh. */ + + /* FIXME: We don't have the required site type yet. Simulate with a + dummy type (point3d). */ + mln::util::graph<mln::point3d> g; + // Populate the graph with vertices. + for (unsigned i = 0; i < mesh.faces.size(); ++i) + g.add_vertex (mln::make::point3d(i, i, i)); + + // Populate the graph with edges. + mesh.need_across_edge(); + for (unsigned f = 0; f < mesh.faces.size(); ++f) + for (unsigned e = 0; e < 3; ++e) + { + int f_adj = mesh.across_edge[f][e]; + if (f_adj != -1) + // Add an edge into the graph. + g.add_edge(f, f_adj); + } + + /*--------------. + | Graph image. | + `--------------*/ + + mln::p_graph<mln::point3d> pg(g); + + typedef mln::graph_image<mln::point3d, curv_t> ima_t; + ima_t g_ima(pg, face_h_inv); + + /*-----------------. + | Simplification. | + `-----------------*/ + + typedef mln::graph_elt_neighborhood<mln::point3d> nbh_t; + nbh_t nbh; + + ima_t closed_g_ima (g_ima.domain()); + mln::morpho::closing_area(g_ima, nbh, lambda, closed_g_ima); + + /*------------------. + | Regional minima. | + `------------------*/ + + typedef unsigned label_t; + label_t nlabels; + typedef mln::graph_image<mln::point3d, label_t> label_ima_t; + label_ima_t minima = + mln::labeling::regional_minima(closed_g_ima, nbh, nlabels); + std::cout << "nlabels = " << nlabels << std::endl; + + /*-----------. + | Skeleton. | + `-----------*/ + + // FIXME: To do. + + /*---------. + | Output. | + `---------*/ + + /* FIXME We should created a boolean graph_image instead. But as + we cannot directly save a graph_image as an OFF file now, just + store the values of this would-be image into an array. */ + // Assign a boolean value to graph vertices (mesh faces). + std::vector<bool> face_value (minima.domain().nvertices(), true); + mln_piter_(label_ima_t) pm(minima.domain()); + for_all(pm) + // FIXME: Use literal::zero. + if (minima(pm) != 0) + { + // The face belongs to a regional minima: ``remove'' it from + // the mesh by tagging it as false. + mln_psite_(label_ima_t) pp(pm); + face_value[pp.id().to_equiv()] = false; + } + + // Taken and adapted from TriMesh_io.cc + FILE* f_out = fopen(output_filename.c_str(), "wb"); + if (!f_out) + { + std::cerr << "Error opening " << output_filename.c_str() + << " for writing." << std::endl; + exit(2); + } + write_off_binary(mesh_ptr, face_value, f_out); + fclose(f_out); + + delete mesh_ptr; +} Index: apps/statues/io.hh --- apps/statues/io.hh (revision 2054) +++ apps/statues/io.hh (working copy) @@ -25,11 +25,17 @@ // reasons why the executable file might be covered by the GNU General // Public License. -#ifndef SANDBOX_LEVILLAIN_IO_HH -# define SANDBOX_LEVILLAIN_IO_HH +#ifndef APPS_STATUES_IO_HH +# define APPS_STATUES_IO_HH +/// \file apps/statues/io.hh /// \brief I/O routines adapted from TriMesh's ones. +#include <algorithm> + +#include <TriMesh.h> + + /// Taken from TriMesh_io.cc /// \{ // Convert colors float -> uchar @@ -112,4 +118,40 @@ /// \} -#endif // ! SANDBOX_LEVILLAIN_IO_HH +/// Taken and adapted from TriMesh_io.cc +/// \{ +/// Write a bunch of faces to an ASCII file, removing the ones tagged +/// as ``false'' in \a face_value. +inline void write_faces_asc_binary(TriMesh *mesh, + const std::vector<bool>& face_value, + FILE *f, + const char *before_face, + const char *after_line) +{ + mesh->need_faces(); + for (int i = 0; i < mesh->faces.size(); i++) + if (face_value[i]) + { + fprintf(f, "%s%d %d %d%s\n", + before_face, + mesh->faces[i][0], mesh->faces[i][1], mesh->faces[i][2], + after_line); + } +} + +/// Write an off file, removing faces tagged as ``false'' in \a face_value. +inline void write_off_binary(TriMesh *mesh, + const std::vector<bool>& face_value, + FILE *f) +{ + fprintf(f, "OFF\n"); + mesh->need_faces(); + unsigned long nfaces = + std::count(face_value.begin(), face_value.end(), true); + fprintf(f, "%lu %lu 0\n", (unsigned long) mesh->vertices.size(), nfaces); + write_verts_asc(mesh, f, "", 0, 0, false, 0, ""); + write_faces_asc_binary(mesh, face_value, f, "3 ", ""); +} +/// \} + +#endif // ! APPS_STATUES_IO_HH Index: apps/statues/Makefile.am --- apps/statues/Makefile.am (revision 2054) +++ apps/statues/Makefile.am (working copy) @@ -7,7 +7,8 @@ # Find trimesh archive. AM_LDFLAGS = -L$(top_builddir)/external/trimesh/libsrc -ltrimesh -bin_PROGRAMS = mesh-segm +bin_PROGRAMS = mesh-segm mesh-skel mesh_segm_SOURCES = mesh-segm.cc io.hh +mesh_skel_SOURCES = mesh-skel.cc io.hh TESTS = test