[Olena-patches] 3106: Add a (non Trimesh-based) curvature computation + segmentation app.

* apps/statues/mesh-complex-max-curv-segm.cc: New. * apps/statues/test-mesh-complex-max-curv-segm.in: New. * apps/statues/Makefile.am (bin_PROGRAMS): Add mesh-complex-max-curv-segm. (mesh-complex-max-curv-segm.cc): New. (TESTS): Add test-mesh-complex-max-curv-segm. --- milena/ChangeLog | 11 ++ milena/apps/statues/Makefile.am | 14 ++ milena/apps/statues/mesh-complex-max-curv-segm.cc | 185 ++++++++++++++++++++ .../statues/test-mesh-complex-max-curv-segm.in | 9 + 4 files changed, 219 insertions(+), 0 deletions(-) create mode 100644 milena/apps/statues/mesh-complex-max-curv-segm.cc create mode 100644 milena/apps/statues/test-mesh-complex-max-curv-segm.in diff --git a/milena/ChangeLog b/milena/ChangeLog index 02fe7c1..f44c118 100644 --- a/milena/ChangeLog +++ b/milena/ChangeLog @@ -1,5 +1,16 @@ 2008-12-30 Roland Levillain &lt;roland(a)lrde.epita.fr&gt; + Add a (non Trimesh-based) curvature computation + segmentation app. + + * apps/statues/mesh-complex-max-curv-segm.cc: New. + * apps/statues/test-mesh-complex-max-curv-segm.in: New. + * apps/statues/Makefile.am + (bin_PROGRAMS): Add mesh-complex-max-curv-segm. + (mesh-complex-max-curv-segm.cc): New. + (TESTS): Add test-mesh-complex-max-curv-segm. + +2008-12-30 Roland Levillain &lt;roland(a)lrde.epita.fr&gt; + Exercise the curvature computation routines ported from Trimesh. * apps/statues/mesh-complex-max-curv.cc: New. diff --git a/milena/apps/statues/Makefile.am b/milena/apps/statues/Makefile.am index d493628..4961960 100644 --- a/milena/apps/statues/Makefile.am +++ b/milena/apps/statues/Makefile.am @@ -77,3 +77,17 @@ TESTS += test-mesh-complex-max-curv bin_PROGRAMS += mesh-complex-segm mesh_complex_segm_SOURCES = mesh-complex-segm.cc TESTS += test-mesh-complex-segm + +# Program computing the max curvature from a mesh, then segmenting it. +# The advantage of this program over `mesh-complex-segm' is that the +# curvature information is much more precise, as data is computed on +# vertices, then average values are computed on edges (in the case of +# `mesh-complex-segm', data are read from polygons, which implies an +# extra averaging step). +bin_PROGRAMS += mesh-complex-max-curv-segm +mesh_complex_max_curv_segm_SOURCES = mesh-complex-max-curv-segm.cc +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'. +# ... diff --git a/milena/apps/statues/mesh-complex-max-curv-segm.cc b/milena/apps/statues/mesh-complex-max-curv-segm.cc new file mode 100644 index 0000000..8a0537e --- /dev/null +++ b/milena/apps/statues/mesh-complex-max-curv-segm.cc @@ -0,0 +1,185 @@ +// 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-complex-max-curv-segm.cc +/// \brief A program computing the maximal curvature values from the +/// surface of the (triangle) mesh of a statue, then performing a +/// WST-based segmentation, using a complex-based image. + +#include <cstdlib> +#include <cmath> + +#include <utility> +#include <iostream> + +#include <mln/core/image/complex_image.hh> +#include <mln/core/image/complex_neighborhoods.hh> + +#include <mln/morpho/closing_area.hh> +#include <mln/morpho/meyer_wst.hh> + +#include <mln/math/max.hh> +#include <mln/math/sqr.hh> + +#include <mln/literal/white.hh> + +#include <mln/io/off/load.hh> +#include <mln/io/off/save.hh> + +#include "trimesh/misc.hh" + + +// Doesn't C++ have a better way to express Pi? +static const float pi = 4 * atanf(1); + + +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]; + unsigned lambda = atoi(argv[2]); + 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; + + mln::bin_2complex_image3df bin_input; + mln::io::off::load(bin_input, input_filename); + std::pair<ima_t, ima_t> curv = mln::geom::mesh_curvature(bin_input.domain()); + + // Compute the pseudo_inverse curvature at each vertex. + ima_t input(bin_input.domain()); + mln::p_n_faces_fwd_piter<D, G> v(input.domain(), 0); + for_all(v) + { + float h = (curv.first(v) + curv.second(v)) / 2; + // Pseudo-inverse curvature. + float h_inv = 1 / pi * (atan(-h) + pi / 2); + input(v) = h_inv; + // FIXME: The program should allow the user to choose the kind + // of measure. +// input(v) = mln::math::max(mln::math::sqr(curv.first(v)), +// mln::math::sqr(curv.second(v))); + } + + // Values on edges. + mln::p_n_faces_fwd_piter<D, G> e(input.domain(), 1); + typedef mln::complex_lower_neighborhood<D, G> adj_vertices_nbh_t; + adj_vertices_nbh_t adj_vertices_nbh; + mln_niter_(adj_vertices_nbh_t) adj_v(adj_vertices_nbh, e); + // Iterate on edges (1-faces). + for_all(e) + { + float s = 0.0f; + unsigned n = 0; + // Iterate on vertices (0-faces). + for_all(adj_v) + { + s += input(adj_v); + ++n; + } + input(e) = s / n; + // An edge should be adjacent to exactly two vertices. + mln_invariant(n <= 2); + } + + /*-----------------. + | Simplification. | + `-----------------*/ + + /// Adjancent edges are connected by shared polygons. + typedef + mln::complex_higher_dim_connected_n_face_neighborhood<D, G> + adj_edges_nbh_t; + adj_edges_nbh_t adj_edges_nbh; + + ima_t closed_input(input.domain()); + mln::morpho::closing_area(input, adj_edges_nbh, lambda, closed_input); + + /*------. + | WST. | + `------*/ + + /* FIXME: Note that the WST is doing too much work, since we have + not constrained the domain of the image to 1-faces only. + It would be great if we could use something like this: + + closed_input | mln::p_faces<1, D, G>(closed_input.domain()) + + as input of the WST. */ + + // Compute the WST on edges. + typedef unsigned wst_val_t; + wst_val_t nbasins; + typedef mln::unsigned_2complex_image3df wst_ima_t; + wst_ima_t wshed = + mln::morpho::meyer_wst(closed_input, adj_edges_nbh, nbasins); + std::cout << "nbasins = " << nbasins << std::endl; + + // Label polygons (i.e., propagate labels from edges to polygons). + typedef mln::complex_higher_neighborhood<D, G> adj_polygons_nbh_t; + adj_polygons_nbh_t adj_polygons_nbh; + mln_niter_(adj_polygons_nbh_t) adj_p(adj_polygons_nbh, e); + for_all(e) + if (wshed(e) != 0) + for_all(adj_p) + wshed(adj_p) = wshed(e); + + /*---------. + | Output. | + `---------*/ + + mln::rgb8_2complex_image3df output(wshed.domain()); + mln::level::fill(output, mln::literal::white); + + // FIXME: Use a colorize functor instead. + // Choose random colors for each basin number. + std::vector<mln::value::rgb8> basin_color (nbasins + 1); + for (unsigned i = 0; i <= nbasins; ++i) + basin_color[i] = mln::value::rgb8(random() % 256, + random() % 256, + random() % 256); + mln_piter_(ima_t) f(wshed.domain()); + for_all(f) + output(f) = basin_color[wshed(f)]; + + mln::io::off::save(output, output_filename); +} diff --git a/milena/apps/statues/test-mesh-complex-max-curv-segm.in b/milena/apps/statues/test-mesh-complex-max-curv-segm.in new file mode 100644 index 0000000..1ec6197 --- /dev/null +++ b/milena/apps/statues/test-mesh-complex-max-curv-segm.in @@ -0,0 +1,9 @@ +#! /bin/sh + +set -ex + +mesh_dir=@top_srcdir@/milena/mesh + +time ./mesh-complex-max-curv-segm $mesh_dir/socket.off 25 socket-complex-max-curv-segm.off +time ./mesh-complex-max-curv-segm $mesh_dir/teapot.off 50 teapot-complex-max-curv-segm.off +time ./mesh-complex-max-curv-segm $mesh_dir/bunny-holefilled.off 500 bunny-holefilled-complex-max-curv-segm.off -- 1.6.0.4