- Olena-patches

3110: Start to factor code among tests on complex-based images.
by Roland Levillain 30 Dec '08

30 Dec '08

* tests/core/image/complex_image.hh: New. * tests/core/image/Makefile.am (noinst_HEADERS): New. Set to complex_image.hh. --- milena/ChangeLog | 8 ++ milena/tests/core/image/Makefile.am | 2 + milena/tests/core/image/complex_image.hh | 165 ++++++++++++++++++++++++++++++ 3 files changed, 175 insertions(+), 0 deletions(-) create mode 100644 milena/tests/core/image/complex_image.hh diff --git a/milena/ChangeLog b/milena/ChangeLog index 623b97b..e218828 100644 --- a/milena/ChangeLog +++ b/milena/ChangeLog @@ -1,5 +1,13 @@ 2008-12-30 Roland Levillain <roland(a)lrde.epita.fr> + Start to factor code among tests on complex-based images. + + * tests/core/image/complex_image.hh: New. + * tests/core/image/Makefile.am (noinst_HEADERS): New. + Set to complex_image.hh. + +2008-12-30 Roland Levillain <roland(a)lrde.epita.fr> + Add a test for apps/statues/test-mesh-max-curv. * apps/statues/test-mesh-max-curv.in: New. diff --git a/milena/tests/core/image/Makefile.am b/milena/tests/core/image/Makefile.am index cfbba60..72f0c93 100644 --- a/milena/tests/core/image/Makefile.am +++ b/milena/tests/core/image/Makefile.am @@ -30,6 +30,8 @@ check_PROGRAMS = \ ## tr_image ## value_enc_image +noinst_HEADERS = complex_image.hh + ##bgraph_image_SOURCES = bgraph_image.cc cast_image_SOURCES = cast_image.cc complex_image_SOURCES = complex_image.cc diff --git a/milena/tests/core/image/complex_image.hh b/milena/tests/core/image/complex_image.hh new file mode 100644 index 0000000..5487d46 --- /dev/null +++ b/milena/tests/core/image/complex_image.hh @@ -0,0 +1,165 @@ +// 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 tests/core/image/complex_image.hh +/// \brief Shared code for tests on complex-based images. + +#ifndef TESTS_CORE_IMAGE_COMPLEX_IMAGE_HH +# define TESTS_CORE_IMAGE_COMPLEX_IMAGE_HH + +# include <mln/core/image/complex_image.hh> +# include <mln/core/alias/complex_image.hh> +# include <mln/value/int_u8.hh> + +// FIXME: Keep this? (See below.) +# include <mln/core/site_set/p_faces.hh> + + +inline +mln::int_u8_2complex_image2d +make_test_complex_image() +{ + /*----------. + | Complex. | + `----------*/ + + /* A 2-d (simplicial) complex and its adjacency graph. + + c 0 1 2 3 + r .------------------------ + | v0 e3 v3 + 0 | o-----------o v0----e3----v3 + | / \ ,-----. / / \ | / + | / . \ \ t1/ / / \ t1 / + 1 | e0 / / \ e1\ / / e4 e0. ,e1Ž `e4 + | / /t0 \ \ ' / / t0 \ / + | / `-----' \ / / | \ / + 2 | o-----------o v1----e2----v2 + | v1 e2 v2 + + v = vertex + e = edge + t = triangle + */ + + + const unsigned D = 2; + + mln::topo::complex<D> c; + + // 0-faces (points). + mln::topo::n_face<0, D> v0 = c.add_face(); + mln::topo::n_face<0, D> v1 = c.add_face(); + mln::topo::n_face<0, D> v2 = c.add_face(); + mln::topo::n_face<0, D> v3 = c.add_face(); + + // 1-faces (segments). + mln::topo::n_face<1, D> e0 = c.add_face(v0 + v1); + mln::topo::n_face<1, D> e1 = c.add_face(v0 + v2); + mln::topo::n_face<1, D> e2 = c.add_face(v1 + v2); + mln::topo::n_face<1, D> e3 = c.add_face(v0 + v3); + mln::topo::n_face<1, D> e4 = c.add_face(v2 + v3); + + // 2-faces (triangles). + mln::topo::n_face<2, D> t0 = c.add_face(e0 + e1 + e2); + mln::topo::n_face<2, D> t1 = c.add_face(e1 + e3 + e4); + + + /*------------------------------. + | Complex geometry (location). | + `------------------------------*/ + + typedef mln::discrete_plane_2complex_geometry G; + G geom; + geom.add_location(mln::point2d(0,1)); // 0-face #0. + geom.add_location(mln::point2d(2,0)); // 0-face #1. + geom.add_location(mln::point2d(2,2)); // 0-face #2. + geom.add_location(mln::point2d(0,3)); // 0-face #3. + + + /*---------------------. + | Complex-based pset. | + `---------------------*/ + + // A pset. + mln::p_complex<D, G> pc(c, geom); + mln::topo::face<D> af(e0); + // An associated psite. + mln::complex_psite<D, G> cs(pc, af); + + + /*--------------------. + | Faces-based psets. | + `--------------------*/ + + /* FIXME: Not that p_faces have become less interesting since the + introduction of p_complex_faces_{fwd,bkd}_piter_. Keep this part + of the test? */ + + /* FIXME: Move this parts out of this function? (into an existing + test?) */ + + // Pset of 0-faces. + mln::p_faces<0, D, G> pf0(c); + // Pset of 1-faces. + mln::p_faces<1, D, G> pf1(c); + // Pset of 2-faces. + mln::p_faces<2, D, G> pf2(c); + + // Some psites on faces. + mln::faces_psite<0, D, G> fs0(pf0, v0); + mln::faces_psite<1, D, G> fs1(pf1, e0); + mln::faces_psite<2, D, G> fs2(pf2, t0); + + + /*----------------------. + | Complex-based image. | + `----------------------*/ + + using mln::value::int_u8; + + // An image type built on a 2-complex with mln::int_u8 values on + // each face. + typedef mln::complex_image<D, G, int_u8> ima_t; + + // Values. + mln::metal::vec<D + 1, std::vector< int_u8 > > values; + // Assign 0 to 0-faces, 1 to 1-faces and 2 to 2-faces. + for (unsigned d = 0; d <= D; ++d) + for (unsigned n = 0; n < pc.cplx().nfaces(d); ++n) + values[d].push_back(d); + + // Create and init an image based on PC. + ima_t ima(pc, values); + + // Check the value associated to edge E0 (through complex psite CS). + mln_postcondition(ima(cs) == 1u); + + return ima; +} + +#endif // ! TESTS_CORE_IMAGE_COMPLEX_IMAGE_HH -- 1.6.0.4

1 0

3109: configure.ac: Configure milena/apps/statues/test-mesh-max-curv.
by Roland Levillain 30 Dec '08

30 Dec '08

--- ChangeLog | 4 ++++ configure.ac | 2 ++ 2 files changed, 6 insertions(+), 0 deletions(-) diff --git a/ChangeLog b/ChangeLog index 61de06b..41080a8 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,5 +1,9 @@ 2008-12-30 Roland Levillain <roland(a)lrde.epita.fr> + * configure.ac: Configure milena/apps/statues/test-mesh-max-curv. + +2008-12-30 Roland Levillain <roland(a)lrde.epita.fr> + * configure.ac: Configure milena/apps/statues/test-mesh-complex-max-curv-segm. diff --git a/configure.ac b/configure.ac index b156eb7..96dbeeb 100644 --- a/configure.ac +++ b/configure.ac @@ -267,6 +267,8 @@ AC_CONFIG_FILES([ # Configure tests. # FIXME: Consider using `sed' instead of `configure' to create these tests. +AC_CONFIG_FILES([milena/apps/statues/test-mesh-max-curv], + [chmod +x milena/apps/statues/test-mesh-max-curv]) AC_CONFIG_FILES([milena/apps/statues/test-mesh-complex-max-curv], [chmod +x milena/apps/statues/test-mesh-complex-max-curv]) AC_CONFIG_FILES([milena/apps/statues/test-mesh-segm], -- 1.6.0.4

1 0

3108: Add a test for apps/statues/test-mesh-max-curv.
by Roland Levillain 30 Dec '08

30 Dec '08

* apps/statues/test-mesh-max-curv.in: New. * apps/statues/Makefile.am (TESTS): Add test-mesh-max-curv. --- milena/ChangeLog | 7 +++++++ milena/apps/statues/Makefile.am | 1 + milena/apps/statues/test-mesh-max-curv.in | 7 +++++++ 3 files changed, 15 insertions(+), 0 deletions(-) create mode 100644 milena/apps/statues/test-mesh-max-curv.in diff --git a/milena/ChangeLog b/milena/ChangeLog index f44c118..623b97b 100644 --- a/milena/ChangeLog +++ b/milena/ChangeLog @@ -1,5 +1,12 @@ 2008-12-30 Roland Levillain <roland(a)lrde.epita.fr> + Add a test for apps/statues/test-mesh-max-curv. + + * apps/statues/test-mesh-max-curv.in: New. + * apps/statues/Makefile.am (TESTS): Add test-mesh-max-curv. + +2008-12-30 Roland Levillain <roland(a)lrde.epita.fr> + Add a (non Trimesh-based) curvature computation + segmentation app. * apps/statues/mesh-complex-max-curv-segm.cc: New. diff --git a/milena/apps/statues/Makefile.am b/milena/apps/statues/Makefile.am index 4961960..b0225c8 100644 --- a/milena/apps/statues/Makefile.am +++ b/milena/apps/statues/Makefile.am @@ -59,6 +59,7 @@ mesh_max_curv_SOURCES = mesh-max-curv.cc io.hh mesh_max_curv_CPPFLAGS = $(AM_CPPFLAGS) $(CPPFLAGS_trimesh) mesh_max_curv_LDFLAGS = $(LDFLAGS_trimesh) mesh_max_curv_LDADD = $(LDADD_trimesh) +TESTS += test-mesh-max-curv ## ---------------------------------------------------- ## ## Complex-based applications, independent of Trimesh. ## diff --git a/milena/apps/statues/test-mesh-max-curv.in b/milena/apps/statues/test-mesh-max-curv.in new file mode 100644 index 0000000..9e9b2ce --- /dev/null +++ b/milena/apps/statues/test-mesh-max-curv.in @@ -0,0 +1,7 @@ +#! /bin/sh + +set -ex + +mesh_dir=@top_srcdir@/milena/mesh + +time ./mesh-max-curv $mesh_dir/teapot.off teapot-max-curv.off -- 1.6.0.4

1 0

3107: configure.ac: Configure milena/apps/statues/test-mesh-complex-max-curv-segm.
by Roland Levillain 30 Dec '08

30 Dec '08

--- ChangeLog | 5 +++++ configure.ac | 2 ++ 2 files changed, 7 insertions(+), 0 deletions(-) diff --git a/ChangeLog b/ChangeLog index e141d92..61de06b 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,6 +1,11 @@ 2008-12-30 Roland Levillain <roland(a)lrde.epita.fr> * configure.ac: Configure + milena/apps/statues/test-mesh-complex-max-curv-segm. + +2008-12-30 Roland Levillain <roland(a)lrde.epita.fr> + + * configure.ac: Configure milena/apps/statues/test-mesh-complex-max-curv. 2008-12-30 Roland Levillain <roland(a)lrde.epita.fr> diff --git a/configure.ac b/configure.ac index aa5c962..b156eb7 100644 --- a/configure.ac +++ b/configure.ac @@ -273,6 +273,8 @@ AC_CONFIG_FILES([milena/apps/statues/test-mesh-segm], [chmod +x milena/apps/statues/test-mesh-segm]) AC_CONFIG_FILES([milena/apps/statues/test-mesh-complex-segm], [chmod +x milena/apps/statues/test-mesh-complex-segm]) +AC_CONFIG_FILES([milena/apps/statues/test-mesh-complex-max-curv-segm], + [chmod +x milena/apps/statues/test-mesh-complex-max-curv-segm]) AC_OUTPUT -- 1.6.0.4

1 0

3106: Add a (non Trimesh-based) curvature computation + segmentation app.
by Roland Levillain 30 Dec '08

30 Dec '08

* 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 <roland(a)lrde.epita.fr> + 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 <roland(a)lrde.epita.fr> + 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

1 0

3105: configure.ac: Configure milena/apps/statues/test-mesh-complex-max-curv.
by Roland Levillain 30 Dec '08

30 Dec '08

--- ChangeLog | 5 +++++ configure.ac | 2 ++ 2 files changed, 7 insertions(+), 0 deletions(-) diff --git a/ChangeLog b/ChangeLog index 61c54b7..e141d92 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,5 +1,10 @@ 2008-12-30 Roland Levillain <roland(a)lrde.epita.fr> + * configure.ac: Configure + milena/apps/statues/test-mesh-complex-max-curv. + +2008-12-30 Roland Levillain <roland(a)lrde.epita.fr> + Configure segmentation tests in milena/apps/statues/. * configure.ac: Configure milena/apps/statues/test-mesh-segm and diff --git a/configure.ac b/configure.ac index 6a4c5a2..aa5c962 100644 --- a/configure.ac +++ b/configure.ac @@ -267,6 +267,8 @@ AC_CONFIG_FILES([ # Configure tests. # FIXME: Consider using `sed' instead of `configure' to create these tests. +AC_CONFIG_FILES([milena/apps/statues/test-mesh-complex-max-curv], + [chmod +x milena/apps/statues/test-mesh-complex-max-curv]) AC_CONFIG_FILES([milena/apps/statues/test-mesh-segm], [chmod +x milena/apps/statues/test-mesh-segm]) AC_CONFIG_FILES([milena/apps/statues/test-mesh-complex-segm], -- 1.6.0.4

1 0

3104: Exercise the curvature computation routines ported from Trimesh.
by Roland Levillain 30 Dec '08

30 Dec '08

* apps/statues/mesh-complex-max-curv.cc: New. * test-mesh-complex-max-curv.in: New test. * apps/statues/Makefile.am (bin_PROGRAMS): Add mesh-complex-max-curv. (mesh_complex_max_curv_SOURCES): New. --- milena/ChangeLog | 10 ++ milena/apps/statues/Makefile.am | 18 ++- milena/apps/statues/mesh-complex-max-curv.cc | 165 +++++++++++++++++++++ milena/apps/statues/test-mesh-complex-max-curv.in | 7 + 4 files changed, 197 insertions(+), 3 deletions(-) create mode 100644 milena/apps/statues/mesh-complex-max-curv.cc create mode 100644 milena/apps/statues/test-mesh-complex-max-curv.in diff --git a/milena/ChangeLog b/milena/ChangeLog index f08a8d9..02fe7c1 100644 --- a/milena/ChangeLog +++ b/milena/ChangeLog @@ -1,5 +1,15 @@ 2008-12-30 Roland Levillain <roland(a)lrde.epita.fr> + Exercise the curvature computation routines ported from Trimesh. + + * apps/statues/mesh-complex-max-curv.cc: New. + * test-mesh-complex-max-curv.in: New test. + * apps/statues/Makefile.am (bin_PROGRAMS): Add + mesh-complex-max-curv. + (mesh_complex_max_curv_SOURCES): New. + +2008-12-30 Roland Levillain <roland(a)lrde.epita.fr> + Port Trimesh's curvature computation routines to Milena. * apps/statues/trimesh/: New directory. diff --git a/milena/apps/statues/Makefile.am b/milena/apps/statues/Makefile.am index 8dd7199..d493628 100644 --- a/milena/apps/statues/Makefile.am +++ b/milena/apps/statues/Makefile.am @@ -21,6 +21,11 @@ LDADD_trimesh = -L$(top_builddir)/external/trimesh/libsrc -ltrimesh $(GLLIBS) # FIXME: Have program names embed the type of structure used (either # graph or complex). +bin_PROGRAMS = +TESTS = + +## FIXME: Programs' names should reflect a dependency on Trimesh, if any. + ## ------------------------------------------------ ## ## Graph-based applications, depending on Trimesh. ## ## ------------------------------------------------ ## @@ -31,6 +36,8 @@ LDADD_trimesh = -L$(top_builddir)/external/trimesh/libsrc -ltrimesh $(GLLIBS) # mesh_segm_CPPFLAGS = $(AM_CPPFLAGS) $(CPPFLAGS_trimesh) # mesh_segm_LDFLAGS = $(LDFLAGS_trimesh) # mesh_segm_LDADD = $(LDADD_trimesh) +# +# TESTS += test-mesh-segm # mesh_skel_SOURCES = mesh-skel.cc io.hh # mesh_skel_CPPFLAGS = $(AM_CPPFLAGS) $(CPPFLAGS_trimesh) @@ -60,8 +67,13 @@ mesh_max_curv_LDADD = $(LDADD_trimesh) noinst_HEADERS = trimesh/misc.hh EXTRA_DIST = trimesh/README +# A small test exercising the curvature computation routines ported +# from Trimesh to Milena. +bin_PROGRAMS += mesh-complex-max-curv +mesh_complex_max_curv_SOURCES = mesh-complex-max-curv.cc +TESTS += test-mesh-complex-max-curv + +# Segmentation program working on precomputed meshes with curvatures data. bin_PROGRAMS += mesh-complex-segm mesh_complex_segm_SOURCES = mesh-complex-segm.cc - - -TESTS = test-mesh-segm test-mesh-complex-segm +TESTS += test-mesh-complex-segm diff --git a/milena/apps/statues/mesh-complex-max-curv.cc b/milena/apps/statues/mesh-complex-max-curv.cc new file mode 100644 index 0000000..0a75da5 --- /dev/null +++ b/milena/apps/statues/mesh-complex-max-curv.cc @@ -0,0 +1,165 @@ +// 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.cc +/// \brief A program computing the max curvature at each (2-)face of +/// the surface of the (triangle) mesh of a statue, 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/level/fill.hh> +#include <mln/literal/zero.hh> + +#include <mln/math/max.hh> +#include <mln/math/sqr.hh> +#include <mln/accu/min_max.hh> +#include <mln/fun/v2v/linear.hh> +#include <mln/level/transform.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 != 3) + { + std::cerr << "usage: " << argv[0] << " input.off output.off" + << std::endl; + std::exit(1); + } + + std::string input_filename = argv[1]; + std::string output_filename = argv[2]; + + /*----------------. + | 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 input; + mln::io::off::load(input, input_filename); + +// // ------------------------------------------------------------ +// // FIXME: TEST. +// mln::complex_image< 2, mln::space_2complex_geometry, mln::algebra::vec<3, float> > +// normal = mesh_normal(input.domain()); +// mln::p_n_faces_fwd_piter<D, G> v_(normal.domain(), 0); +// for_all(v_) +// std::cout << normal(v_) << std::endl; +// std::exit(0); +// // ------------------------------------------------------------ + + std::pair<ima_t, ima_t> curv = mln::geom::mesh_curvature(input.domain()); + + // Compute the max curvature at each vertex. + ima_t max_curv(input.domain()); + mln::level::fill(max_curv, mln::literal::zero); + mln::p_n_faces_fwd_piter<D, G> v(max_curv.domain(), 0); + for_all(v) + max_curv(v) = mln::math::max(mln::math::sqr(curv.first(v)), + mln::math::sqr(curv.second(v))); + + // Propagate these values on triangles. + mln::p_n_faces_fwd_piter<D, G> t(max_curv.domain(), 2); + typedef mln::complex_m_face_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, t); + /* FIXME: Not really user friendly! The `m' value should pass at + the construction of ADJ_V. */ + adj_v.iter().set_m(0); + mln::accu::min_max<float> acc; + // Iterate on triangles (2-faces). + for_all(t) + { + float s = 0.0f; + unsigned n = 0; + // Iterate on vertices (0-faces). + for_all(adj_v) + { + s += max_curv(adj_v); + ++n; + } + float m = s / n; + max_curv(t) = m; + acc.take(m); + // A triangle should be adjacent to exactly three vertices. + mln_invariant(n <= 3); + } + + // Normalize values between 0 and 1. + /* Shrink the values of FACE_M into the range 0..1, as these are + the only values accepted a an RGB floating-point component in the + OFF file format. */ + ima_t output(max_curv.domain()); + mln::level::fill(output, mln::literal::zero); + std::pair<float, float> min_max(acc); + // FIXME: Taken from mln/level/stretch.hh (this should be factored). + float min = min_max.first; + float max = min_max.second; + std::cout << min << std::endl; + std::cout << max << std::endl; + // Don't normalize actually if the curvature is constant (i.e., + // if min == max). + if (min != max) + { + float m = 0.0f; + float M = 1.0f; + float a = (M - m) / (max - min); + float b = (m * max - M * min) / (max - min); + mln::fun::v2v::linear<float, float, float> f(a, b); + output = mln::level::transform(max_curv, f); + } + + // Output. + mln::io::off::save(output, output_filename); +} diff --git a/milena/apps/statues/test-mesh-complex-max-curv.in b/milena/apps/statues/test-mesh-complex-max-curv.in new file mode 100644 index 0000000..e138a12 --- /dev/null +++ b/milena/apps/statues/test-mesh-complex-max-curv.in @@ -0,0 +1,7 @@ +#! /bin/sh + +set -ex + +mesh_dir=@top_srcdir@/milena/mesh + +time ./mesh-complex-max-curv $mesh_dir/teapot.off teapot-complex-max-curv.off -- 1.6.0.4

1 0

3103: Port Trimesh's curvature computation routines to Milena.
by Roland Levillain 30 Dec '08

30 Dec '08

* apps/statues/trimesh/: New directory. * apps/statues/trimesh/misc.hh: New. * apps/statues/trimesh/README: New. * apps/statues/Makefile.am (noinst_HEADERS): New. Add trimesh/misc.hh. (EXTRA_DIST): New. Add trimesh/README. --- milena/ChangeLog | 12 + milena/apps/statues/Makefile.am | 3 + milena/apps/statues/trimesh/README | 4 + milena/apps/statues/trimesh/misc.hh | 657 +++++++++++++++++++++++++++++++++++ 4 files changed, 676 insertions(+), 0 deletions(-) create mode 100644 milena/apps/statues/trimesh/README create mode 100644 milena/apps/statues/trimesh/misc.hh diff --git a/milena/ChangeLog b/milena/ChangeLog index fae6012..f08a8d9 100644 --- a/milena/ChangeLog +++ b/milena/ChangeLog @@ -1,5 +1,17 @@ 2008-12-30 Roland Levillain <roland(a)lrde.epita.fr> + Port Trimesh's curvature computation routines to Milena. + + * apps/statues/trimesh/: New directory. + * apps/statues/trimesh/misc.hh: New. + * apps/statues/trimesh/README: New. + * apps/statues/Makefile.am (noinst_HEADERS): New. + Add trimesh/misc.hh. + (EXTRA_DIST): New. + Add trimesh/README. + +2008-12-30 Roland Levillain <roland(a)lrde.epita.fr> + Simplify the normalization of vectors. * mln/algebra/vec.hh (mln::algebra::vec<n, T>::normalize): diff --git a/milena/apps/statues/Makefile.am b/milena/apps/statues/Makefile.am index 0d1bfc1..8dd7199 100644 --- a/milena/apps/statues/Makefile.am +++ b/milena/apps/statues/Makefile.am @@ -57,6 +57,9 @@ mesh_max_curv_LDADD = $(LDADD_trimesh) ## Complex-based applications, independent of Trimesh. ## ## ---------------------------------------------------- ## +noinst_HEADERS = trimesh/misc.hh +EXTRA_DIST = trimesh/README + bin_PROGRAMS += mesh-complex-segm mesh_complex_segm_SOURCES = mesh-complex-segm.cc diff --git a/milena/apps/statues/trimesh/README b/milena/apps/statues/trimesh/README new file mode 100644 index 0000000..20fe1e7 --- /dev/null +++ b/milena/apps/statues/trimesh/README @@ -0,0 +1,4 @@ +Files from this directory are adaptation of some Trimesh code into +Milena. We might want to push them into Milena someday. + +Don't forget to add authors / copyright owners where necessary. diff --git a/milena/apps/statues/trimesh/misc.hh b/milena/apps/statues/trimesh/misc.hh new file mode 100644 index 0000000..78af913 --- /dev/null +++ b/milena/apps/statues/trimesh/misc.hh @@ -0,0 +1,657 @@ +// 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. + +// FIXME: Split this file. +// FIXME: Address license-related issues? + +#ifndef MILENA_APPS_STATUES_MISC_HH +# define MILENA_APPS_STATUES_MISC_HH + +# include <algorithm> // For std::swap. +# include <utility> // For std::pair. + +# include <mln/algebra/vec.hh> +# include <mln/algebra/mat.hh> + +# include <mln/norm/l2.hh> + + +/** See http://gcc.gnu.org/onlinedocs/gcc-4.3.2/gcc/Other-Builtins.html + for a definition of __builtin_expect. + + Macros likely and unlikely are for optimization purpose only (the + developer can use them to help the compiler improve the + basic-block sewing pass. */ +// FIXME: We should turn `likely' and `unlikely' into functions. +#ifndef likely +# if !defined(__GNUC__) || (__GNUC__ == 2 && __GNUC_MINOR__ < 96) +# define likely(x) (x) +# define unlikely(x) (x) +# else +# define likely(x) (__builtin_expect((x), 1)) +# define unlikely(x) (__builtin_expect((x), 0)) +# endif +#endif + +namespace mln +{ + + namespace algebra + { + + // From Trimesh. + // FIXME: Doc (input, return value, etc.). + // FIXME: Idea: Add reference to Numerical Recipes. + + /** \brief Perform LDL^T decomposition of a symmetric positive + definite matrix. Like Cholesky, but no square roots. + Overwrites lower triangle of matrix. + + From Trimesh's ldltdc routine. */ + template <unsigned N, typename T> + inline + bool + ldlt_decomp(mat<N, N, T>& A, vec<N, T>& rdiag) + { + vec<N - 1, T> v; + for (int i = 0; i < N; ++i) + { + for (int k = 0; k < i; ++k) + v[k] = A(i, k) * rdiag[k]; + for (int j = i; j < N; ++j) + { + T sum = A(i, j); + for (int k = 0; k < i; k++) + sum -= v[k] * A(j, k); + if (i == j) + { + if (unlikely(sum <= T(0))) + return false; + rdiag[i] = T(1) / sum; + } + else + A(j, i) = sum; + } + } + return true; + } + + // From Trimesh. + // FIXME: Doc (input, return value, etc.). + // FIXME: Idea: Add reference to Numerical Recipes. + + /// \brief Solve A x = B after mln::algebra::ldlt_decomp. + template <unsigned N, typename T> + inline + void + ldlt_solve(const mat<N, N, T>& A, const vec<N, T>& rdiag, + const vec<N, T>& B, vec<N, T>& x) + { + for (unsigned i = 0; i < N; ++i) + { + T sum = B[i]; + for (int k = 0; k < i; ++k) + sum -= A(i, k) * x[k]; + x[i] = sum * rdiag[i]; + } + for (int i = N - 1; i >= 0; --i) + { + T sum = 0; + for (int k = i + 1; k < N; ++k) + sum += A(k, i) * x[k]; + x[i] -= sum * rdiag[i]; + } + } + + } // end of namespace mln::algebra + +} // end of namespace mln + + +// FIXME: To be moved elsewhere (in mln/geom)? + +// Inspired from Trimesh. + +namespace mln +{ + + namespace geom + { + + // FIXME: More doc (see Trimesh's). + + /* FIXME: We should restrict attached data to vertices in return + value. */ + /** \brief Compute normals at vertices. + + Inspired from the method Trimesh::need_normals of the Trimesh + library. + \see http://www.cs.princeton.edu/gfx/proj/trimesh2/ + + For simplicity purpose, and contrary to Trimesh, this routine + only compute normals from a mesh, not from a cloud of + points. */ + inline + complex_image< 2, mln::space_2complex_geometry, algebra::vec<3, float> > + mesh_normal(const p_complex<2, space_2complex_geometry>& mesh) + { + // Shortcuts. + static const unsigned D = 2; + typedef space_2complex_geometry G; + typedef algebra::vec<3, float> vec3f; + typedef complex_image< D, G, vec3f > normal_t; + + normal_t normal(mesh); + level::fill(normal, literal::zero); + + mln::p_n_faces_fwd_piter<D, G> f(mesh, 2); + // A neighborhood where neighbors are the set of 0-faces + // transitively adjacent to the reference point. + typedef mln::complex_m_face_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, 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. Alas, we can't inherit ctors (yet), so we + can't rely on a simple mixin approach. */ + adj_v.iter().set_m(0); + + // Iterate on 2-faces (triangles). + for_all(f) + { + /// Psites of 0-faces transitively adjacent to F. + std::vector<mln_psite_(normal_t)> p; + p.reserve(3); + for_all(adj_v) + p.push_back(adj_v); + /// The should be exactly three vertices adjacent to F. + mln_assertion(p.size() == 3); + + /* FIXME: Is the explicit conversion to_site() really needed? + If not, we could provide a local shortcut like + + vec(p[0]) // Instead of p[0].to_site().front().to_vec() + + to shorten these lines. */ + // FIXME: Factor as well? + vec3f a = + p[0].to_site().front().to_vec() - p[1].to_site().front().to_vec(); + vec3f b = + p[1].to_site().front().to_vec() - p[2].to_site().front().to_vec(); + vec3f c = + p[2].to_site().front().to_vec() - p[0].to_site().front().to_vec(); + + // FIXME: Factor as well? + float l2a = norm::sqr_l2(a); + float l2b = norm::sqr_l2(b); + float l2c = norm::sqr_l2(c); + vec3f face_normal = algebra::vprod(a, b); + + normal(p[0]) += face_normal * (1.0f / (l2a * l2c)); + normal(p[1]) += face_normal * (1.0f / (l2b * l2a)); + normal(p[2]) += face_normal * (1.0f / (l2c * l2b)); + } + + // Normalize normals. We don't need to iterate on all faces, just + // 0-faces. + mln::p_n_faces_fwd_piter<D, G> v(mesh, 0); + for_all(v) + normal(v).normalize(); + + return normal; + } + + + /* FIXME: We should restrict attached data to vertices in return + value. */ + /** \brief Compute the area ``belonging'' to normals at vertices. + + Inspired from the method Trimesh::need_pointareas of the + Trimesh library. + \see http://www.cs.princeton.edu/gfx/proj/trimesh2/ + + From the documentation of Trimesh: + + ``Compute the area "belonging" to each vertex or each corner + of a triangle (defined as Voronoi area restricted to the + 1-ring of a vertex, or to the triangle).'' + */ + inline + std::pair< + complex_image< 2, mln::space_2complex_geometry, algebra::vec<3, float> >, + complex_image< 2, mln::space_2complex_geometry, float > + > + mesh_corner_point_area(const p_complex<2, space_2complex_geometry>& mesh) + { + // Shortcuts. + static const unsigned D = 2; + typedef space_2complex_geometry G; + typedef algebra::vec<3, float> vec3f; + + // Hold data for 2-faces only. + typedef complex_image< D, G, vec3f > corner_area_t; + // Hold data for 0-faces only. + typedef complex_image< D, G, float > point_area_t; + // Packed output. + typedef std::pair<corner_area_t, point_area_t> output_t; + + // Initialize output. + output_t output(mesh, mesh); + corner_area_t& corner_area = output.first; + point_area_t& point_area = output.second; + level::fill(corner_area, literal::zero); + level::fill(point_area, literal::zero); + + mln::p_n_faces_fwd_piter<D, G> f(mesh, 2); + // A neighborhood where neighbors are the set of 0-faces + // transitively adjacent to the reference point. + typedef mln::complex_m_face_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, 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. Alas, we can't inherit ctors (yet), so we + can't rely on a simple mixin approach. */ + adj_v.iter().set_m(0); + + for_all(f) + { + /// Psites of 0-faces transitively adjacent to F. + std::vector<mln_psite_(corner_area_t)> p; + p.reserve(3); + for_all(adj_v) + p.push_back(adj_v); + /// The should be exactly three vertices adjacent to F. + mln_assertion(p.size() == 3); + + // (Opposite) edge vectors. + algebra::vec<3, vec3f> e; + // FIXME: See above. + e.set + (p[2].to_site().front().to_vec() - p[1].to_site().front().to_vec(), + p[0].to_site().front().to_vec() - p[2].to_site().front().to_vec(), + p[1].to_site().front().to_vec() - p[0].to_site().front().to_vec()); + + // Compute corner weights. + float area = 0.5f * norm::l2(algebra::vprod(e[0], e[1])); + vec3f sqr_norm; + sqr_norm.set(norm::sqr_l2(e[0]), + norm::sqr_l2(e[1]), + norm::sqr_l2(e[2])); + vec3f edge_weight; + edge_weight.set + (sqr_norm[0] * (sqr_norm[1] + sqr_norm[2] - sqr_norm[0]), + sqr_norm[1] * (sqr_norm[2] + sqr_norm[0] - sqr_norm[1]), + sqr_norm[2] * (sqr_norm[0] + sqr_norm[1] - sqr_norm[2])); + + if (edge_weight[0] <= 0.0f) + { + corner_area(f)[1] = -0.25f * sqr_norm[2] * area / (e[0] * e[2]); + corner_area(f)[2] = -0.25f * sqr_norm[1] * area / (e[0] * e[1]); + corner_area(f)[0] = area - corner_area(f)[1] - corner_area(f)[2]; + } + else if (edge_weight[1] <= 0.0f) + { + corner_area(f)[2] = -0.25f * sqr_norm[0] * area / (e[1] * e[0]); + corner_area(f)[0] = -0.25f * sqr_norm[2] * area / (e[1] * e[2]); + corner_area(f)[1] = area - corner_area(f)[2] - corner_area(f)[0]; + } + else if (edge_weight[2] <= 0.0f) + { + corner_area(f)[0] = -0.25f * sqr_norm[1] * area / (e[2] * e[1]); + corner_area(f)[1] = -0.25f * sqr_norm[0] * area / (e[2] * e[0]); + corner_area(f)[2] = area - corner_area(f)[0] - corner_area(f)[1]; + } + else + { + float ewscale = + 0.5f * area / (edge_weight[0] + edge_weight[1] + edge_weight[2]); + for (int i = 0; i < 3; ++i) + corner_area(f)[i] = + ewscale * (edge_weight[(i+1)%3] + edge_weight[(i+2)%3]); + } + + for (int i = 0; i < 3; ++i) + point_area(p[i]) += corner_area(f)[i]; + } + + return output; + } + + + namespace internal + { + + /** \brief \a i + 1 and \a i - 1 modulo 3 + + Defined as macros in Trimesh. + + #define NEXT(i) ((i)<2 ? (i)+1 : (i)-2) + #define PREV(i) ((i)>0 ? (i)-1 : (i)+2) + + According to Trimesh doc, + + ``This way of computing it tends to be faster than using %'' + + but I (Roland) just can't believe it... We should profile + some code first. + + \{ */ + static inline unsigned next(unsigned i) { return (i + 1) % 3; } + static inline unsigned prev(unsigned i) { return (i - 1) % 3; } + /** \} */ + + + /// Rotate a coordinate system to be perpendicular to the given normal + inline + void + rot_coord_sys(const algebra::vec<3, float> &old_u, + const algebra::vec<3, float> &old_v, + const algebra::vec<3, float> &new_norm, + algebra::vec<3, float> &new_u, + algebra::vec<3, float> &new_v) + { + new_u = old_u; + new_v = old_v; + algebra::vec<3, float> old_norm = vprod(old_u, old_v); + float ndot = old_norm * new_norm; + if (unlikely(ndot <= -1.0f)) + { + new_u = -new_u; + new_v = -new_v; + return; + } + algebra::vec<3, float> perp_old = new_norm - ndot * old_norm; + algebra::vec<3, float> dperp = + 1.0f / (1 + ndot) * (old_norm + new_norm); + new_u -= dperp * (new_u * perp_old); + new_v -= dperp * (new_v * perp_old); + } + + + // FIXME: Add const to formals whenever we can. + + // Reproject a curvature tensor from the basis spanned by old_u and old_v + // (which are assumed to be unit-length and perpendicular) to the + // new_u, new_v basis. + inline + void + proj_curv(const algebra::vec<3, float>& old_u, + const algebra::vec<3, float>& old_v, + float old_ku, float old_kuv, float old_kv, + const algebra::vec<3, float>& new_u, + const algebra::vec<3, float>& new_v, + float& new_ku, float& new_kuv, float& new_kv) + { + algebra::vec<3, float> r_new_u, r_new_v; + rot_coord_sys(new_u, new_v, vprod(old_u, old_v), r_new_u, r_new_v); + + float u1 = r_new_u * old_u; + float v1 = r_new_u * old_v; + float u2 = r_new_v * old_u; + float v2 = r_new_v * old_v; + new_ku = old_ku * u1*u1 + old_kuv * (2.0f * u1*v1) + old_kv * v1*v1; + new_kuv = old_ku * u1*u2 + old_kuv * (u1*v2 + u2*v1) + old_kv * v1*v2; + new_kv = old_ku * u2*u2 + old_kuv * (2.0f * u2*v2) + old_kv * v2*v2; + } + + /** Given a curvature tensor, find principal directions and + curvatures. Makes sure that pdir1 and pdir2 are + perpendicular to normal. */ + inline + void + diagonalize_curv(const algebra::vec<3, float>& old_u, + const algebra::vec<3, float>& old_v, + float ku, float kuv, float kv, + const algebra::vec<3, float>& new_norm, + algebra::vec<3, float>& pdir1, + algebra::vec<3, float>& pdir2, + float& k1, float& k2) + { + algebra::vec<3, float> r_old_u, r_old_v; + rot_coord_sys(old_u, old_v, new_norm, r_old_u, r_old_v); + + float c = 1, s = 0, tt = 0; + if (likely(kuv != 0.0f)) + { + // Jacobi rotation to diagonalize. + float h = 0.5f * (kv - ku) / kuv; + tt = (h < 0.0f) ? + 1.0f / (h - sqrt(1.0f + h*h)) : + 1.0f / (h + sqrt(1.0f + h*h)); + c = 1.0f / sqrt(1.0f + tt*tt); + s = tt * c; + } + + k1 = ku - tt * kuv; + k2 = kv + tt * kuv; + + if (fabs(k1) >= fabs(k2)) + pdir1 = c*r_old_u - s*r_old_v; + else + { + std::swap(k1, k2); + pdir1 = s*r_old_u + c*r_old_v; + } + pdir2 = vprod(new_norm, pdir1); + } + + } // end of namespace mln::geom::internal + + + /** \brief Compute the principal curvatures of a surface at + vertices. + + These princpal curvatures are names kappa_1 and kappa_2 in + + Sylvie Philipp-Foliguet, Michel Jordan Laurent Najman and + Jean Cousty. Artwork 3D Model Database Indexing and + Classification. + + \param[in] ima The surface (triangle mesh) on which the + curvature is to be computed. */ + /* FIXME: We should restrict attached data to vertices in return + value. */ + /* FIXME: Have a typedef for this return type? */ + inline + std::pair< + complex_image< 2, mln::space_2complex_geometry, float >, + complex_image< 2, mln::space_2complex_geometry, float > + > + mesh_curvature(const p_complex<2, space_2complex_geometry>& mesh) + { + // Shortcuts. + static const unsigned D = 2; + typedef space_2complex_geometry G; + typedef algebra::vec<3, float> vec3f; + typedef algebra::mat<3, 3, float> mat3f; + + typedef complex_image< D, G, vec3f > corner_area_t; + typedef complex_image< D, G, float > point_area_t; + + // Normals at vertices. + typedef complex_image< 2, mln::space_2complex_geometry, vec3f > normal_t; + normal_t normal = mesh_normal(mesh); + + // Areas ``belonging'' to a normal. + typedef complex_image< D, G, vec3f > corner_area_t; + typedef complex_image< D, G, float > point_area_t; + typedef std::pair<corner_area_t, point_area_t> corner_point_area_t; + corner_point_area_t corner_point_area = mesh_corner_point_area(mesh); + // Shortcuts. + corner_area_t& corner_area = corner_point_area.first; + point_area_t& point_area = corner_point_area.second; + + // Curvatures at each vertex. + typedef complex_image< 2, mln::space_2complex_geometry, float > curv_t; + typedef std::pair<curv_t, curv_t> output_t; + output_t output(mesh, mesh); + curv_t& curv1 = output.first; + curv_t& curv2 = output.second; + // ?? + complex_image< 2, mln::space_2complex_geometry, float > curv12(mesh); + // Principal directions at each vertex. + complex_image< 2, mln::space_2complex_geometry, vec3f > pdir1(mesh); + complex_image< 2, mln::space_2complex_geometry, vec3f > pdir2(mesh); + + /*-------------------------------------------------. + | Set up an initial coordinate system per vertex. | + `-------------------------------------------------*/ + + mln::p_n_faces_fwd_piter<D, G> f(mesh, 2); + // A neighborhood where neighbors are the set of 0-faces + // transitively adjacent to the reference point. + typedef mln::complex_m_face_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, 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. Alas, we can't inherit ctors (yet), so we + can't rely on a simple mixin approach. */ + adj_v.iter().set_m(0); + + for_all(f) + { + /// Psites of 0-faces transitively adjacent to F. + std::vector<mln_psite_(curv_t)> p; + p.reserve(3); + for_all(adj_v) + p.push_back(adj_v); + /// The should be exactly three vertices adjacent to F. + mln_assertion(p.size() == 3); + + // FIXME: See above. + pdir1(p[0]) = + p[1].to_site().front().to_vec() - p[0].to_site().front().to_vec(); + pdir1(p[1]) = + p[2].to_site().front().to_vec() - p[1].to_site().front().to_vec(); + pdir1(p[2]) = + p[0].to_site().front().to_vec() - p[2].to_site().front().to_vec(); + } + + mln::p_n_faces_fwd_piter<D, G> v(mesh, 0); + for_all(v) + { + pdir1(v) = algebra::vprod(pdir1(v), normal(v)); + pdir1(v).normalize(); + pdir2(v) = algebra::vprod(normal(v), pdir1(v)); + } + + /*-----------------------------. + | Compute curvature per-face. | + `-----------------------------*/ + + for_all(f) + { + std::vector<mln_psite_(curv_t)> p; + p.reserve(3); + for_all(adj_v) + p.push_back(adj_v); + mln_assertion(p.size() == 3); + + // (Opposite) edge vectors. + algebra::vec<3, vec3f> e; + // FIXME: See above. + e.set + (p[2].to_site().front().to_vec() - p[1].to_site().front().to_vec(), + p[0].to_site().front().to_vec() - p[2].to_site().front().to_vec(), + p[1].to_site().front().to_vec() - p[0].to_site().front().to_vec()); + + // N-T-B coordinate system per face. + vec3f t = e[0]; + t.normalize(); + vec3f n = algebra::vprod(e[0], e[1]); + /* FIXME: Why don't we normalize N? Is it a normal vector by + construction? Or maybe we don't need it to be normal? */ + vec3f b = algebra::vprod(n, t); + b.normalize(); + + // Estimate curvature based on variation of normals along edges. + vec3f m = literal::zero; + mat3f w = literal::zero; + for (int j = 0; j < 3; ++j) + { + float u = e[j] * t; + float v = e[j] * b; + w(0, 0) += u * u; + w(0, 1) += u * v; + /* FIXME: Those two lines were commented in Trimesh's + original code. + + //w(1, 1) += v*v + u*u; + //w(1, 2) += u*v; + + */ + w(2, 2) += v * v; + vec3f dn = + normal(p[internal::prev(j)]) - normal(p[internal::next(j)]); + float dnu = dn * t; + float dnv = dn * b; + m[0] += dnu*u; + m[1] += dnu*v + dnv*u; + m[2] += dnv*v; + } + w(1, 1) = w(0, 0) + w(2, 2); + w(1, 2) = w(0, 1); + + // Least squares solution. + vec3f diag; + bool ldlt_decomp_sucess_p = algebra::ldlt_decomp(w, diag); + mln_assertion(ldlt_decomp_sucess_p); + algebra::ldlt_solve(w, diag, m, m); + + // Push it back out to the vertices + for (int j = 0; j < 3; ++j) + { + float c1, c12, c2; + internal::proj_curv(t, b, m[0], m[1], m[2], + pdir1(p[j]), pdir2(p[j]), c1, c12, c2); + float wt = corner_area(f)[j] / point_area(p[j]); + curv1(p[j]) += wt * c1; + curv12(p[j]) += wt * c12; + curv2(p[j]) += wt * c2; + } + } + + /*-------------------------------------------------------------. + | Compute principal directions and curvatures at each vertex. | + `-------------------------------------------------------------*/ + + for_all(v) + internal::diagonalize_curv(pdir1(v), pdir2(v), + curv1(v), curv12(v), curv2(v), + normal(v), pdir1(v), pdir2(v), + curv1(v), curv2(v)); + + return output; + } + + } // end of namespace mln::geom + +} // end of namespace mln + +#endif // ! MILENA_APPS_STATUES_MISC_HH -- 1.6.0.4

1 0

3102: Simplify the normalization of vectors.
by Roland Levillain 30 Dec '08

30 Dec '08

* mln/algebra/vec.hh (mln::algebra::vec<n, T>::normalize): Use mln::norm::l2 instead of ad hoc code. --- milena/ChangeLog | 7 +++++++ milena/mln/algebra/vec.hh | 30 +++++++++++++++++++++++------- 2 files changed, 30 insertions(+), 7 deletions(-) diff --git a/milena/ChangeLog b/milena/ChangeLog index 0f3bc29..fae6012 100644 --- a/milena/ChangeLog +++ b/milena/ChangeLog @@ -1,5 +1,12 @@ 2008-12-30 Roland Levillain <roland(a)lrde.epita.fr> + Simplify the normalization of vectors. + + * mln/algebra/vec.hh (mln::algebra::vec<n, T>::normalize): + Use mln::norm::l2 instead of ad hoc code. + +2008-12-30 Roland Levillain <roland(a)lrde.epita.fr> + More complex-related aliases. * mln/core/alias/complex_geometry.hh diff --git a/milena/mln/algebra/vec.hh b/milena/mln/algebra/vec.hh index 2a61693..1434930 100644 --- a/milena/mln/algebra/vec.hh +++ b/milena/mln/algebra/vec.hh @@ -39,6 +39,7 @@ # include <mln/core/concept/object.hh> # include <mln/literal/zero.hh> +# include <mln/norm/l2.hh> # include <mln/trait/all.hh> # include <mln/trait/value_.hh> # include <mln/fun/i2v/all_to.hh> @@ -64,7 +65,10 @@ namespace mln struct zero_t; } - + namespace norm { + template <unsigned n, typename C> + mln_sum(C) l2(const algebra::vec<n,C>& vec); + } namespace trait @@ -211,6 +215,18 @@ namespace mln unsigned size() const; + /* FIXME: What if the vector is null? Even if we choose not to + handle this case, we should *state* in the documentation the + behavior of this method. + + I (Roland) have added an assertion to detect ``erroneous'' + cases, but we might want something different. + ``Implementation defined'' or ``undefined behavior'' is fine + by me, as long as the documentation mentions it. + + FWIW, Trimesh's developers chose to set all the coordinates + of a vector being normalized to 1, when it norm is equal or + lower (sic) to zero. */ const vec<n, T>& normalize(); /// Constructor; coordinates are set by function \p f. @@ -306,7 +322,8 @@ namespace mln operator-(const vec<n,T>& lhs, const vec<n,U>& rhs); // vec * vec - + + /// Scalar product (dot product). template <unsigned n, typename T, typename U> mln_sum_x(T,U) operator*(const vec<n,T>& lhs, const vec<n,U>& rhs); @@ -331,6 +348,7 @@ namespace mln // vprod // FIXME: Generalize... + /// Vectorial product (cross product). template <typename T, typename U> vec<3, mln_trait_op_times(T,U)> // FIXME: Sum of product... vprod(const vec<3, T>& lhs, const vec<3, U>& rhs); @@ -431,12 +449,10 @@ namespace mln inline const vec<n, T>& vec<n, T>::normalize() { - float n_l2 = 0; - for (unsigned i = 0; i < n; ++i) - n_l2 += data_[i] * data_[i]; - n_l2 = float(std::sqrt(n_l2)); + float l2_norm = float(norm::l2(*this)); + mln_assertion(l2_norm > 0.f); for (unsigned i = 0; i < n; ++i) - data_[i] = static_cast<T>(data_[i] / n_l2); + data_[i] = static_cast<T>(data_[i] / l2_norm); return *this; } -- 1.6.0.4

1 0

3101: More complex-related aliases.
by Roland Levillain 30 Dec '08

30 Dec '08

* mln/core/alias/complex_geometry.hh (mln::discrete_plane_2complex_geometry) * mln/core/alias/complex_image.hh (mln::int_u8_2complex_image2d): New typedefs. Typos & aesthetic changes. --- milena/ChangeLog | 11 +++++++++++ milena/mln/core/alias/complex_geometry.hh | 7 ++++++- milena/mln/core/alias/complex_image.hh | 20 ++++++++++++++++++-- 3 files changed, 35 insertions(+), 3 deletions(-) diff --git a/milena/ChangeLog b/milena/ChangeLog index aa168d3..0f3bc29 100644 --- a/milena/ChangeLog +++ b/milena/ChangeLog @@ -1,5 +1,16 @@ 2008-12-30 Roland Levillain <roland(a)lrde.epita.fr> + More complex-related aliases. + + * mln/core/alias/complex_geometry.hh + (mln::discrete_plane_2complex_geometry) + * mln/core/alias/complex_image.hh + (mln::int_u8_2complex_image2d): + New typedefs. + Typos & aesthetic changes. + +2008-12-30 Roland Levillain <roland(a)lrde.epita.fr> + Clean up apps/statues/ a bit. * apps/statues/mesh-max-curv.cc, apps/statues/mesh-pinv-curv.cc, diff --git a/milena/mln/core/alias/complex_geometry.hh b/milena/mln/core/alias/complex_geometry.hh index d4e3a75..44fa899 100644 --- a/milena/mln/core/alias/complex_geometry.hh +++ b/milena/mln/core/alias/complex_geometry.hh @@ -40,8 +40,13 @@ namespace mln { /// \brief Type alias for the geometry of a 2-complex located in a + /// discrete 2-dimensional plane (with integer coordinates). + typedef mln::geom::complex_geometry<2, point2d> + discrete_plane_2complex_geometry; + + /// \brief Type alias for the geometry of a 2-complex located in a /// 3-dimensional space (with floating-point coordinates). - typedef mln::geom::complex_geometry<2,point3df> space_2complex_geometry; + typedef mln::geom::complex_geometry<2, point3df> space_2complex_geometry; } // end of namespace mln diff --git a/milena/mln/core/alias/complex_image.hh b/milena/mln/core/alias/complex_image.hh index 6336d6b..0d8e379 100644 --- a/milena/mln/core/alias/complex_image.hh +++ b/milena/mln/core/alias/complex_image.hh @@ -40,6 +40,22 @@ namespace mln { + /*------------------------------. + | 2-d plane 2-complex aliases. | + `------------------------------*/ + + /// \brief Type alias for an 8-bit gray-level image based on a + /// 2-complex, where 0-faces are located at discrete (integer) + /// 2-dimensional points. + typedef + mln::complex_image<2, mln::discrete_plane_2complex_geometry, + mln::value::int_u8> + int_u8_2complex_image2d; + + + /*------------------------------. + | 3-d space 2-complex aliases. | + `------------------------------*/ /// \brief Type alias for a binary image based on a 2-complex, where /// 0-faces are located at floating-point 3-dimensional points. @@ -47,14 +63,14 @@ namespace mln mln::complex_image<2, mln::space_2complex_geometry, bool> bin_2complex_image3df; - /// \brief Type alias for an 8-bit grey-level image based on a + /// \brief Type alias for an 8-bit gray-level image based on a /// 2-complex, where 0-faces are located at floating-point /// 3-dimensional points. typedef mln::complex_image<2, mln::space_2complex_geometry, mln::value::int_u8> int_u8_2complex_image3df; - /// \brief Type alias for a grey-level image based on a 2-complex, + /// \brief Type alias for a gray-level image based on a 2-complex, /// where 0-faces are located at floating-point 3-dimensional /// points. typedef -- 1.6.0.4

1 0