[Olena-patches] [PATCH 05/10] Write and application computing skeletons on complex-based images.

* apps/statues/mesh-complex-skel.cc: New. * apps/statues/save_bin_alt.hh: New. Workaround over mln::io::off::save's limitations. * mln/topo/complex.hh [NDEBUG]: Remove a warning. * apps/statues/Makefile.am (AM_CXXFLAGS): Add -Wall -Wextra. (bin_PROGRAMS): Add mesh-complex-skel. (mesh_complex_skel_SOURCES): New. --- milena/ChangeLog | 12 + milena/apps/statues/Makefile.am | 16 +- milena/apps/statues/mesh-complex-skel.cc | 704 ++++++++++++++++++++++++++++++ milena/apps/statues/save_bin_alt.hh | 188 ++++++++ milena/mln/topo/complex.hh | 2 + 5 files changed, 921 insertions(+), 1 deletions(-) create mode 100644 milena/apps/statues/mesh-complex-skel.cc create mode 100644 milena/apps/statues/save_bin_alt.hh diff --git a/milena/ChangeLog b/milena/ChangeLog index ce5e952..ab602ac 100644 --- a/milena/ChangeLog +++ b/milena/ChangeLog @@ -1,5 +1,17 @@ 2009-04-06 Roland Levillain <roland@lrde.epita.fr> + Write and application computing skeletons on complex-based images. + + * apps/statues/mesh-complex-skel.cc: New. + * apps/statues/save_bin_alt.hh: New. + Workaround over mln::io::off::save's limitations. + * mln/topo/complex.hh [NDEBUG]: Remove a warning. + * apps/statues/Makefile.am (AM_CXXFLAGS): Add -Wall -Wextra. + (bin_PROGRAMS): Add mesh-complex-skel. + (mesh_complex_skel_SOURCES): New. + +2009-04-06 Roland Levillain <roland@lrde.epita.fr> + Allow the comparison of faces of different dimensions using `<'. * mln/topo/face.hh (operator< (const face<D>&, const face<D>&)): diff --git a/milena/apps/statues/Makefile.am b/milena/apps/statues/Makefile.am index e132895..9c7c9ca 100644 --- a/milena/apps/statues/Makefile.am +++ b/milena/apps/statues/Makefile.am @@ -8,7 +8,8 @@ include $(top_srcdir)/external/trimesh/gluit/gluit.mk AM_CPPFLAGS = -I$(top_srcdir)/milena CPPFLAGS_trimesh = -I$(top_srcdir)/external/trimesh/include # Produce fast code. -AM_CXXFLAGS = -O3 -DNDEBUG -ggdb +# FIXME: This is GCC-dependent; use a variable computed by configure instead. +AM_CXXFLAGS = -Wall -Wextra -ggdb -O3 -DNDEBUG # Find the trimesh library and its dependencies. # # Don't use TRIMESH_LDFLAGS, since it looks like the LDFLAGS of the @@ -39,6 +40,9 @@ TESTS = # # TESTS += test-mesh-segm +# FIXME: mesh_skel is unfinished. Anyway, it should be superseded by +# another program, using no Trimesh code. + # mesh_skel_SOURCES = mesh-skel.cc io.hh # mesh_skel_CPPFLAGS = $(AM_CPPFLAGS) $(CPPFLAGS_trimesh) # mesh_skel_LDFLAGS = $(LDFLAGS_trimesh) @@ -68,6 +72,9 @@ TESTS += test-mesh-max-curv noinst_HEADERS = trimesh/misc.hh EXTRA_DIST = trimesh/README +## Segmentation. +## ------------ + # A small program exercising the curvature computation routines ported # from Trimesh to Milena. bin_PROGRAMS += mesh-complex-max-curv @@ -96,3 +103,10 @@ TESTS += test-mesh-complex-max-curv-segm # FIXME: Implement `mesh-complex-pinv-curv-segm' (factor as much # code as possible with `mesh-complex-max-curv-segm'. # ... + +## Skeletonization. +## ---------------- + +# Skeletonization program working on precomputed meshes with curvatures data. +bin_PROGRAMS += mesh-complex-skel +mesh_complex_skel_SOURCES = mesh-complex-skel.cc diff --git a/milena/apps/statues/mesh-complex-skel.cc b/milena/apps/statues/mesh-complex-skel.cc new file mode 100644 index 0000000..e58c5cd --- /dev/null +++ b/milena/apps/statues/mesh-complex-skel.cc @@ -0,0 +1,704 @@ +// Copyright (C) 2008, 2009 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-complex-skel.cc +/// \brief A program computing a skeleton of the surface of the +/// (triangle) mesh of a statue (by thinning), using a complex-based +/// image. + +#include <cstdlib> +#include <cmath> + +#include <algorithm> +#include <utility> +#include <iostream> + +#include <mln/core/image/complex_image.hh> +#include <mln/core/image/complex_neighborhoods.hh> +#include <mln/core/image/image_if.hh> + +#include <mln/core/site_set/p_set.hh> + +#include <mln/pw/value.hh> +#include <mln/pw/cst.hh> + +#include <mln/value/label_16.hh> +#include <mln/literal/white.hh> + +#include <mln/core/routine/duplicate.hh> +#include <mln/labeling/regional_minima.hh> + +#include <mln/io/off/load.hh> +#include <mln/io/off/save.hh> +/* FIXME: Remove as soon as mln::io::off::save is able to save a + morphed mln::complex_image (i.e., seen through image_if). */ +#include "save_bin_alt.hh" + + +// FIXME: This is bad. Remove as soon as all routines have been moved +// to the library. +using namespace mln; + +// Fwd. decl. +template <typename I> struct is_simple_triangle; + +// // A very simple and limited predicate for the simplicity of a point. +// template <typename I> +// class is_simple_triangle +// : public mln::Function_p2b< is_simple_triangle<I> > +// { +// public: +// /// Dimension of the image (and therefore of the complex). +// static const unsigned D = I::dim; +// /// Geometry of the image. +// typedef mln_geom(I) G; +// /// Psite type. +// typedef mln::complex_psite<D, G> psite; + +// /// Result type of the functor. +// typedef bool result; + +// is_simple_triangle() +// : ima_(0), nbh_() +// { +// } + +// is_simple_triangle(const mln::Image<I>& ima) +// : ima_(mln::exact(&ima)), nbh_() +// { +// } + +// /* FIXME: Rename as init() or something like this? See how other +// functors are written. */ +// /// Set the underlying image. +// void set_image(const mln::Image<I>& ima) +// { +// ima_ = mln::exact(&ima); +// } + +// bool operator()(const psite& p) const +// { +// mln_precondition(ima_); +// unsigned nneighbs = 0; +// mln_niter(nbh_t) n(nbh_, p); +// for_all(n) +// if ((*ima_)(n)) +// ++nneighbs; +// mln_invariant(n <= 3); +// /* FIXME: Dummy: A triangle is considered simple if it is not +// completely attached nor completely detached (i.e, if it has +// more than 0 adjacent triangles but less than three). This test +// is obiously simplistic and wrong, but we'll use it as a first +// approximation. We should use the notion of free pair later. +// And provide efficient, mask-based version as well. */ +// return 0 < nneighbs && nneighbs < 3; +// } + +// private: +// const I* ima_; + +// // FIXME: Make this type a parameter of the functor? +// typedef mln::complex_lower_dim_connected_n_face_neighborhood<D, G> nbh_t; +// nbh_t nbh_; +// }; + + +// FIXME: Doc. +template <unsigned D, typename G> +bool +is_facet(const mln::complex_psite<D, G>& f) +{ + typedef mln::complex_higher_neighborhood<D, G> higher_adj_nbh_t; + higher_adj_nbh_t higher_adj_nbh; + mln_niter(higher_adj_nbh_t) n(higher_adj_nbh, f); + for_all(n) + // If the neighborhood is not empty, then F is included in a face + // of higher dimension. + return false; + // Otherwise, F is a facet. + return true; +} + + +/** \brief Compute the set of faces of the cell corresponding to the + facet \a f. + + \pre \a f is a facet (it does not belong to any face of higher + dimension) + + \return a set of faces containing the attachment. */ +template <unsigned D, typename G> +mln::p_set< complex_psite<D, G> > +cell(const complex_psite<D, G>& f) +{ + // Ensure the face F is a facet. + mln_precondition(is_facet(f)); + + typedef complex_psite<D, G> psite; + typedef p_set<psite> faces_t; + + // Compute the cell F^HAT. + faces_t f_hat; + /* FIXME: We need a cell-iterator here + (see https://trac.lrde.org/olena/ticket/162). */ + typedef complex_m_face_neighborhood<D, G> m_faces_nbh_t; + m_faces_nbh_t m_faces_nbh; + mln_niter(m_faces_nbh_t) g(m_faces_nbh, f); + for (unsigned m = 0; m < f.n(); ++m) + { + g.iter().set_m(m); + for_all(g) + { + f_hat.insert(g); +// // FIXME: Debug. +// std::cerr << g << std::endl; + } + } + f_hat.insert(f); +// // FIXME: Debug. +// std::cerr << f << std::endl; + return f_hat; +} + + +/** \brief Compute the attachment of the cell corresponding to the + facet \a f to the image \a ima. + + \pre ima is an image of Boolean values. + + \return a set of faces containing the attachment. + + We do not use the fomal definition of the attachment here (see + couprie.08.pami). We use the following (equivalent) definition: + an N-face F in CELL is in the attachment of CELL to IMA if it is + adjacent to at least an (N-1)-face or an (N+1)-face that does not + belong to CELL. */ +template <unsigned D, typename G, typename V> +mln::p_set< complex_psite<D, G> > +attachment(const mln::complex_psite<D, G>& f, + const mln::complex_image<D, G, V>& ima) +{ + mlc_equal(V, bool)::check(); + + typedef complex_psite<D, G> psite; + typedef p_set<psite> faces_t; + + faces_t f_hat = cell(f); +// // FIXME: Debug. +// std::cerr << "f_hat.nsites() = " << f_hat.nsites() << std::endl; + faces_t att_f; + + typedef complex_lower_higher_neighborhood<D, G> adj_nbh_t; + adj_nbh_t adj_nbh; + mln_piter(faces_t) g(f_hat); + mln_niter(adj_nbh_t) n(adj_nbh, g); + for_all(g) + for_all(n) + if (ima(n) && !f_hat.has(n)) + { + att_f.insert(g); + break; + } +// // FIXME: Debug. +// std::cerr << "att_f.nsites() = " << att_f.nsites() << std::endl; + return att_f; +} + + +/** \brief Compute the detachment of the cell corresponding to the + facet \a f to the image \a ima. + + \pre ima is an image of Boolean values. + + \return a set of faces containing the detachment. + + We do not use the fomal definition of the detachment here (see + couprie.08.pami). We use the following (equivalent) definition: + an N-face F in CELL is not in the detachment of CELL from IMA if + it is adjacent to at least an (N-1)-face or an (N+1)-face that + does not belong to CELL. */ +template <unsigned D, typename G, typename V> +mln::p_set< complex_psite<D, G> > +detachment(const mln::complex_psite<D, G>& f, + const mln::complex_image<D, G, V>& ima) +{ + mlc_equal(V, bool)::check(); + + typedef complex_psite<D, G> psite; + typedef p_set<psite> faces_t; + + faces_t f_hat = cell(f); + // Initialize DETACH_F to F_HAT. + faces_t detach_f = f_hat; + + typedef complex_lower_higher_neighborhood<D, G> adj_nbh_t; + adj_nbh_t adj_nbh; + mln_piter(faces_t) g(f_hat); + mln_niter(adj_nbh_t) n(adj_nbh, g); + for_all(g) + for_all(n) + if (ima(n) && !f_hat.has(n)) + { + detach_f.remove(g); + break; + } + return detach_f; +} + + +/** \brief A predicate for the simplicity of a point based on the + collapse property of the attachment. + + The functor does not actually take a cell as input, but a face + that is expected to be a 2-facet. */ +template <typename I> +class is_simple_cell + : public mln::Function_p2b< is_simple_cell<I> > +{ +public: + /// Dimension of the image (and therefore of the complex). + static const unsigned D = I::dim; + /// Geometry of the image. + typedef mln_geom(I) G; + /// Psite type. + typedef mln::complex_psite<D, G> psite; + + /// Result type of the functor. + typedef bool result; + + is_simple_cell() + : ima_(0) + { + } + + is_simple_cell(const mln::Image<I>& ima) + : ima_(mln::exact(&ima)) + { + } + + /* FIXME: Rename as init() or something like this? See how other + functors are written. */ + /// Set the underlying image. + void set_image(const mln::Image<I>& ima) + { + ima_ = mln::exact(&ima); + } + + // Based on the algorithm A2 from couprie.08.pami. + bool operator()(const psite& p) const + { +// // FIXME: Debug. +// std::cerr << "is_simple_cell(" << p << ")" << std::endl; + + mln_precondition(ima_); + + typedef p_set<psite> faces_t; + + // Compute the attachment of the cell corresponding to P to he + // domain of *IMA_. + faces_t att = attachment(p, *ima_); + + // A cell with an empty attachment is not simple. + /* FIXME: Why does p_set not provide an empty() predicate? */ + if (att.nsites() == 0) + return false; + + /* FIXME: This part could be moved to its own function/method + (`is_collapsible'). Moreover, the code could be split: looking + up for a free pair (g, h) could be performed in another + routine. */ + // Try to collapse the attachment (to a single point). + { + typedef complex_lower_neighborhood<D, G> lower_adj_nbh_t; + typedef complex_higher_neighborhood<D, G> higher_adj_nbh_t; + lower_adj_nbh_t lower_adj_nbh; + higher_adj_nbh_t higher_adj_nbh; + while (att.nsites() > 1) + { + + bool simple_pair_collapsed = false; + /* FIXME: The selection of G and H is probably suboptimal + (i.e., selecting G has a face of highest avalaible + dimension in CELL is probably smarter). */ + mln_piter(faces_t) g(att); + for_all(g) + /* G cannot have dimension 0, since we later look for an + adjacent face H of lower dimension. */ + if (static_cast<psite>(g).n() > 0) + { + // Check whether G is a facet within ATT. + bool g_is_facet = true; + mln_niter(higher_adj_nbh_t) f(higher_adj_nbh, g); + for_all(f) + if (att.has(f)) + { + g_is_facet = false; + break; + } + if (!g_is_facet) + break; + + // Look for a face H stricly included in G. + bool gh_is_simple_pair = false; + mln_niter(lower_adj_nbh_t) h(lower_adj_nbh, g); + for_all(h) + { + bool h_strictly_in_g = true; + if (att.has(h)) + { + mln_niter(higher_adj_nbh_t) i(higher_adj_nbh, h); + for_all(i) + if (i != g && att.has(i)) + { + h_strictly_in_g = false; + break; + } + } + if (h_strictly_in_g) + { + gh_is_simple_pair = true; + att.remove(g); + att.remove(h); + mln_invariant(att.nsites() > 0); + break; + } + } // for_all(h) + + // If a free pair (G, H) has been found and removed, + // restart the free pair look up from the beginning. + if (gh_is_simple_pair) + { + simple_pair_collapsed = true; + break; + } + } // for_all(g) + + if (!simple_pair_collapsed) + // If no free pair (G, H) was found, then the attachment is + // not collapsible. + return false; + + } // while (att.nsites() > 1) + + mln_postcondition(att.nsites() == 1); + mln_postcondition(att[0].n() == 0); + // If the attachment is collapsible to a 0-face, then the cell + // corresponding to the (face) P is simple. + return true; + } + } + +private: + const I* ima_; +}; + + +// FIXME: Doc. +template <unsigned D, typename G, typename V> +void +detach (const mln::complex_psite<D, G>f, mln::complex_image<D, G, V>& ima) +{ + mlc_equal(V, bool)::check(); + + typedef mln::complex_psite<D, G> psite; + typedef p_set<psite> faces_t; + + // Compute the detachment of P from IMA. + faces_t detach = detachment(f, ima); + // Detach all its faces from IMA. +#if 0 + // FIXME: Does not work yet? Check again. + data::fill(ima | detach, false); +#endif + mln_piter(faces_t) p(detach); + for_all(p) + ima(p) = false; +} + + +// FIXME: Move this elsewhere. +// FIXME: Use a generic `I' instead of `mln::bin_2complex_image3df'. +// FIXME: Add a constraint/priority argument (K/P). +/** \brief Breadth-First Thinning. + + A non generic implementation of a binary skeleton on a triangle + surface (mesh). + + \a N is the adjacency relation between triangles, but it probably + does not make sense to use something other than an instance of + mln::complex_lower_dim_connected_n_face_neighborhood if the image type is + hard-coded as mln::bin_2complex_image3df. */ +template <typename N, typename F> +mln::bin_2complex_image3df +skeleton(const mln::bin_2complex_image3df& input, + const mln::Neighborhood<N>& nbh_, + mln::Function_p2b<F>& is_simple_) +{ + const N& nbh = exact(nbh_); + F& is_simple = exact(is_simple_); + + typedef mln::bin_2complex_image3df I; + typedef mln_psite(I) psite; + + I output = mln::duplicate(input); + // Attach the work image to IS_SIMPLE. + is_simple.set_image(output); + + typedef mln::p_set<psite> set_t; + set_t set; + +// // FIXME: Debug. +// unsigned count = 0; + + // Populate the set with candiate simple points. + mln_piter(I) p(output.domain()); + for_all(p) + { + /* FIXME: Cheat! We'd like to iterate on 2-cells only, but + we cannot constrain the domain of the input using + image_if (yet). Hence this filter to restrict the + dimension of the considered cells. A less dirty + approaach would be to use a constraint predicate (see + comment above). */ + if (p.unproxy_().n() != I::dim) + continue; + +// // FIXME: Debug. +// std::cerr << "skeleton: Iteration (for_all(p)) : " +// << count++ << std::endl; + +// // FIXME: Debug. +// std::cerr << " output(p) = " << output(p) << std::endl; +// std::cerr << " is_simple(p) = " << is_simple(p) << std::endl; + + if (output(p) && is_simple(p)) + set.insert(p); + } + +// // FIXME: Debug. +// std::cerr << "skeleton: set.nsites() = " << set.nsites() << std::endl; + +// // FIXME: Debug. +// count = 0; + + while (!set.is_empty()) + { +// // FIXME: Debug. +// mln_piter(set_t) ps(set); +// std::cerr << "set = " << std::endl; +// for_all(ps) +// std::cerr << " " << static_cast<psite>(ps) << std::endl; +// std::cerr << std::endl; + +// // FIXME: Debug. +// std::cerr << "skeleton: Iteration (while (!set.is_empty())) : " +// << count++ << std::endl; + + set_t next_set; + // FIXME: Use a piter on SET instead of this hand-made iteration. + for (unsigned i = 0; i < set.nsites(); ++i) + { + psite p = set[i]; + + /* FIXME: Cheat! We'd like to iterate on 2-cells only, but + we cannot constrain the domain of the input using + image_if (yet). Hence this filter to restrict the + dimension of the considered cells. A less dirty + approaach would be to use a constraint predicate (see + comment above). */ + if (p.n() != I::dim) + continue; + + /* FIXME: We compute the cell and attachment of P twice: + within is_simple and within detach. How could we reuse + this elegantly, without breaking the genericity of the + skeleton algorithm? */ + if (is_simple(p)) + { + // FIXME: `detach' could be a functor, as is `is_simple'. + detach(p, output); + mln_niter(N) n(nbh, p); + for_all(n) + { + /* FIXME: Cheat! We'd like to iterate on 2-cells only, but + we cannot constrain the domain of the input using + image_if (yet). Hence this filter to restrict the + dimension of the considered cells. A less dirty + approaach would be to use a constraint predicate (see + comment above). */ + if (n.unproxy_().n() != I::dim) + continue; + + if (output.domain().has(n) && output(n) && is_simple(n)) + next_set.insert(n); + } + } + } + set.clear(); + std::swap(set, next_set); + } + + return output; +} + + +// Fwd. decl. +template <unsigned N> struct is_n_face; + +/// A functor testing wheter a mln::complex_psite is an \N-face. +template <unsigned N> +struct is_n_face : public mln::Function_p2b< is_n_face<N> > +{ + typedef bool result; + + template <unsigned D, typename G> + bool operator()(const mln::complex_psite<D, G>& p) + { + return p.n() == N; + } +}; + + +int main(int argc, char* argv[]) +{ + if (argc != 4) + { + std::cerr << "usage: " << argv[0] << " input.off lambda output.off" + << std::endl; + std::exit(1); + } + + std::string input_filename = argv[1]; + // FIXME: Use lambda to filter the input. +#if 0 + unsigned lambda = atoi(argv[2]); +#endif + std::string output_filename = argv[3]; + + /*----------------. + | Complex image. | + `----------------*/ + + // Image type. + typedef mln::float_2complex_image3df ima_t; + // Dimension of the image (and therefore of the complex). + static const unsigned D = ima_t::dim; + // Geometry of the image. + typedef mln_geom_(ima_t) G; + + ima_t input; + mln::io::off::load(input, input_filename); + + /* FIXME: Workaround: Set maximal values on vertices and edges to + exclude them from the set of minimal values. */ + mln::p_n_faces_fwd_piter<D, G> v(input.domain(), 0); + for_all(v) + input(v) = mln_max(float); + mln::p_n_faces_fwd_piter<D, G> e(input.domain(), 1); + for_all(e) + input(e) = mln_max(float); + + /*---------------. + | Local minima. | + `---------------*/ + + typedef mln::value::label_16 label_t; + label_t nminima; + + typedef mln::complex_lower_dim_connected_n_face_neighborhood<D, G> nbh_t; + nbh_t nbh; + + /* FIXME: We should use something like `ima_t | p_n_faces(2)' instead + of `ima_t' here. Or better: `input' should only associate data + to 2-faces. */ + mln_ch_value_(ima_t, label_t) minima = + mln::labeling::regional_minima(input, nbh, nminima); + + /*-----------------------. + | Initial binary image. | + `-----------------------*/ + + typedef mln_ch_value_(ima_t, bool) bin_ima_t; + bin_ima_t surface(input.domain()); + mln::data::fill(surface, true); + // Dig ``holes'' in the surface surface by setting minima faces to false. +#if 0 + /* FIXME: The `is_n_face<2>' predicate is required because the + domain of SURFACE also comprises 0-faces and 1-faces... */ + mln::data::fill(surface + | is_n_face<2>() + | mln::pw::value(minima) != mln::pw::cst(mln::literal::zero), + false); +#endif + /* FIXME: The code above does not compile (yet). Probably because + of an incompatibility between complex_image and image_if. Use + this hand-made filling for the moment. */ + mln::p_n_faces_fwd_piter<D, G> t(input.domain(), 2); + for_all(t) + if (minima(t) != mln::literal::zero) + surface(t) = false; + +// #if 0 +// /* FIXME: Debug. But alas, this does not work (yet). +// Use workaround mln::io::off::save_bin_salt instead (bad!) */ +// # if 0 +// mln::io::off::save(surface | mln::pw::value(surface) == mln::pw::cst(true), +// "surface.off"); +// # endif +// mln::io::off::save_bin_alt(surface, "surface.off"); +// #endif + +// // FIXME: Debug. +// typedef mln_ch_value_(ima_t, bool) bin_ima_t; +// bin_ima_t surface; +// mln::io::off::load(surface, +// "/Users/roland/src/olena/milena/mesh/two-triangles.off"); +// // Set one of the two triangle to false. +// p_n_faces_fwd_piter<D, G> t(surface.domain(), 2); +// t.start(); +// mln_precondition(t.is_valid()); +// surface(t) = false; + + /*-----------. + | Skeleton. | + `-----------*/ + + is_simple_cell<bin_ima_t> is_simple_p; + bin_ima_t skel = skeleton(surface, nbh, is_simple_p); + + /*---------. + | Output. | + `---------*/ + + /* FIXME: This does not work (yet). + Use workaround mln::io::off::save_bin_salt instead (bad!) */ +#if 0 + mln::io::off::save(skel | mln::pw::value(skel) == mln::pw::cst(true), + output_filename); +#endif + mln::io::off::save_bin_alt(skel, output_filename); +} diff --git a/milena/apps/statues/save_bin_alt.hh b/milena/apps/statues/save_bin_alt.hh new file mode 100644 index 0000000..b2def9b --- /dev/null +++ b/milena/apps/statues/save_bin_alt.hh @@ -0,0 +1,188 @@ +// Copyright (C) 2009 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. + +#ifndef APPS_STATUES_SAVE_BIN_ALT_HH +# define APPS_STATUES_SAVE_BIN_ALT_HH + +/*--------------------------------------------------------------------. +| FIXME: Copied and adjusted (in a hurry) from mln/io/off/save.hh, | +| because fixing image_if + complex_image was much too long. Sorry. | +`--------------------------------------------------------------------*/ + +# include <cstdlib> + +# include <iostream> +# include <fstream> +# include <sstream> + +# include <string> + +# include <mln/core/alias/complex_image.hh> +# include <mln/core/image/complex_neighborhoods.hh> +# include <mln/core/image/complex_neighborhood_piter.hh> + +namespace mln +{ + + namespace io + { + + namespace off + { + + /** FIXME: Similar to + mln::io::off::save(const bin_2complex_image3df&, const std::string&), + but does not save faces whose value is `false'. */ + template <typename I> + void save_bin_alt(const I& ima, + const std::string& filename) + { + const std::string me = "mln::io::off::save"; + + std::ofstream ostr(filename.c_str()); + if (!ostr) + { + std::cerr << me << ": `" << filename << "' invalid file." + << std::endl; + /* FIXME: Too violent. We should allow the use of + exceptions, at least to have Milena's code behave + correctly in interpreted environments (std::exit() or + std::abort() causes the termination of a Python + interpreter, for instance!). */ + std::exit(1); + } + + /*---------. + | Header. | + `---------*/ + + static const unsigned D = I::dim; + typedef mln_geom(I) G; + + /* ``The .off files in the Princeton Shape Benchmark conform + to the following standard.'' */ + + /* ``OFF files are all ASCII files beginning with the + keyword OFF. '' */ + ostr << "OFF" << std::endl; + + // A comment. + ostr << "# Generated by Milena 1.0 http://olena.lrde.epita.fr\n" + << "# EPITA Research and Development Laboratory (LRDE)" + << std::endl; + + // Count the number of 2-faces set to `true'. + unsigned n2faces = 0; + p_n_faces_fwd_piter<D, G> f(ima.domain(), 2); + for_all(f) + if (ima(f)) + ++n2faces; + + /* ``The next line states the number of vertices, the number + of faces, and the number of edges. The number of edges can + be safely ignored.'' */ + /* FIXME: This is too long. We shall be able to write + + ima.domain().template nfaces<0>() + + or even + + ima.template nfaces<0>(). + */ + ostr << ima.domain().cplx().template nfaces<0>() << ' ' + << n2faces << ' ' + << ima.domain().cplx().template nfaces<1>() << std::endl; + + /*-------. + | Data. | + `-------*/ + + // ------------------------------------------ // + // Vertices & geometry (vertices locations). // + // ------------------------------------------ // + + /* ``The vertices are listed with x, y, z coordinates, written + one per line.'' */ + + // Traverse the 0-faces (vertices). + p_n_faces_fwd_piter<D, G> v(ima.domain(), 0); + for_all(v) + { + mln_invariant(v.to_site().size() == 1); + ostr << v.to_site().front()[0] << ' ' + << v.to_site().front()[1] << ' ' + << v.to_site().front()[2] << std::endl; + } + + // --------------- // + // Faces & edges. // + // --------------- // + + /* ``After the list of vertices, the faces are listed, with one + face per line. For each face, the number of vertices is + specified, followed by indices into the list of + vertices.'' */ + + // Traverse the 2-faces (polygons). + typedef complex_m_face_neighborhood<D, G> nbh_t; + // A neighborhood where neighbors are the set of 0-faces + // transitively adjacent to the reference point. + nbh_t nbh; + mln_fwd_niter(nbh_t) u(nbh, f); + /* FIXME: We should be able to pas this value (m) either at + the construction of the neighborhood or at the construction + of the iterator. */ + u.iter().set_m(0); + + // For each (2-)face, iterate on (transitively) ajacent + // vertices (0-faces). + for_all(f) + if (ima(f)) + { + unsigned nvertices = 0; + std::ostringstream vertices; + for_all(u) + { + // FIXME: Likewise, this is a bit too long... + vertices << ' ' << u.unproxy_().face().face_id(); + ++nvertices; + } + ostr << nvertices << vertices.str(); + ostr << std::endl; + } + + ostr.close(); + } + + } // end of namespace mln::io::off + + } // end of namespace mln::io + +} // end of namespace mln + + +#endif // ! APPS_STATUES_SAVE_BIN_ALT_HH diff --git a/milena/mln/topo/complex.hh b/milena/mln/topo/complex.hh index 2b392ac..c182743 100644 --- a/milena/mln/topo/complex.hh +++ b/milena/mln/topo/complex.hh @@ -951,6 +951,8 @@ namespace mln { // If we reached this method, then N should be 0. mln_precondition(n == 0); + // Prevent ``unused variable'' warnings when NDEBUG is defined. + (void) n; return f(faces_); } -- 1.6.1.2