Commit 4b83d9c6 authored by Roland Levillain's avatar Roland Levillain
Browse files

mesh-complex-pinv-curv-segm and mesh-complex-pinv-curv-skel apps.

	* apps/mesh-segm-skel/mesh-complex-pinv-curv-segm.cc,
	* apps/mesh-segm-skel/mesh-complex-pinv-curv-skel.cc:
	New.
	* apps/mesh-segm-skel/mesh-complex-max-curv-segm.cc,
	* apps/mesh-segm-skel/mesh-complex-max-curv-skel.cc:
	Remove dead code.
	* apps/mesh-segm-skel/test-mesh-complex-pinv-curv-segm.in,
	* apps/mesh-segm-skel/test-mesh-complex-pinv-curv-skel.in:
	New.
	* apps/mesh-segm-skel/Makefile.am (bin_PROGRAMS): Add
	mesh-complex-pinv-curv-segm and mesh-complex-pinv-curv-skel.
	(mesh_complex_pinv_curv_segm_SOURCES)
	(mesh_complex_pinv_curv_skel_SOURCES):
	New.
	(TESTS): Add test-mesh-complex-pinv-curv-segm and
	test-mesh-complex-pinv-curv-skel.
	(MOSTLYCLEANFILES): Add socket-complex-pinv-curv-segm.off,
	teapot-complex-pinv-curv-segm.off,
	socket-complex-pinv-curv-skel.off and
	teapot-complex-pinv-curv-skel.off.
parent 29bf3b02
2010-07-15 Roland Levillain <roland@lrde.epita.fr>
mesh-complex-pinv-curv-segm and mesh-complex-pinv-curv-skel apps.
* apps/mesh-segm-skel/mesh-complex-pinv-curv-segm.cc,
* apps/mesh-segm-skel/mesh-complex-pinv-curv-skel.cc:
New.
* apps/mesh-segm-skel/mesh-complex-max-curv-segm.cc,
* apps/mesh-segm-skel/mesh-complex-max-curv-skel.cc:
Remove dead code.
* apps/mesh-segm-skel/test-mesh-complex-pinv-curv-segm.in,
* apps/mesh-segm-skel/test-mesh-complex-pinv-curv-skel.in:
New.
* apps/mesh-segm-skel/Makefile.am (bin_PROGRAMS): Add
mesh-complex-pinv-curv-segm and mesh-complex-pinv-curv-skel.
(mesh_complex_pinv_curv_segm_SOURCES)
(mesh_complex_pinv_curv_skel_SOURCES):
New.
(TESTS): Add test-mesh-complex-pinv-curv-segm and
test-mesh-complex-pinv-curv-skel.
(MOSTLYCLEANFILES): Add socket-complex-pinv-curv-segm.off,
teapot-complex-pinv-curv-segm.off,
socket-complex-pinv-curv-skel.off and
teapot-complex-pinv-curv-skel.off.
2010-06-14 Roland Levillain <roland@lrde.epita.fr>
 
Actually use the maximal curvature in mesh-complex-max-curv-segm.
......@@ -102,9 +102,13 @@ MOSTLYCLEANFILES += \
socket-complex-max-curv-segm.off \
teapot-complex-max-curv-segm.off
# FIXME: Implement `mesh-complex-pinv-curv-segm' (factor as much
# code as possible with `mesh-complex-max-curv-segm'.
# ...
# Likewise, but using pseudo-inverse curvature.
bin_PROGRAMS += mesh-complex-pinv-curv-segm
mesh_complex_pinv_curv_segm_SOURCES = mesh-complex-pinv-curv-segm.cc
TESTS += test-mesh-complex-pinv-curv-segm
MOSTLYCLEANFILES += \
socket-complex-pinv-curv-segm.off \
teapot-complex-pinv-curv-segm.off
## Skeletonization.
## ----------------
......@@ -118,8 +122,6 @@ MOSTLYCLEANFILES += \
teapot-max-curv-skel.off \
three-triangles-skel.off
# FIXME: Implement `mesh-complex-{max,pinv}-curv-skel'.
# ...
bin_PROGRAMS += mesh-complex-max-curv-skel
mesh_complex_max_curv_skel_SOURCES = \
mesh-complex-max-curv-skel.cc save_bin_alt.hh
......@@ -127,3 +129,11 @@ TESTS += test-mesh-complex-max-curv-skel
MOSTLYCLEANFILES += \
socket-complex-max-curv-skel.off \
teapot-complex-max-curv-skel.off
bin_PROGRAMS += mesh-complex-pinv-curv-skel
mesh_complex_pinv_curv_skel_SOURCES = \
mesh-complex-pinv-curv-skel.cc save_bin_alt.hh
TESTS += test-mesh-complex-pinv-curv-skel
MOSTLYCLEANFILES += \
socket-complex-pinv-curv-skel.off \
teapot-complex-pinv-curv-skel.off
// Copyright (C) 2008, 2009, 2011, 2013 EPITA Research and Development
// Copyright (C) 2008, 2009, 2010, 2011, 2013 EPITA Research and Development
// Laboratory (LRDE)
//
// This file is part of Olena.
......@@ -29,6 +29,8 @@
/// surface of the (triangle) mesh of a statue, then performing a
/// WST-based segmentation, using a complex-based image.
// Factor with mesh-complex-pinv-curv-segm.cc.
#include <cstdlib>
#include <cmath>
......@@ -87,18 +89,9 @@ int main(int argc, char* argv[])
mln::p_n_faces_fwd_piter<D, G> v(float_ima.domain(), 0);
for_all(v)
{
// FIXME: The program should allow the user to choose the kind
// of measure (curvature).
#if 0
float h = (curv.first(v) + curv.second(v)) / 2;
// Pseudo-inverse curvature.
float h_inv =
1 / float(mln::math::pi) * (atanf(-h) + float(mln::math::pi) / 2);
float_ima(v) = h_inv;
#endif
// Max curvature.
float_ima(v) = mln::math::max(mln::math::sqr(curv.first(v)),
mln::math::sqr(curv.second(v)));
mln::math::sqr(curv.second(v)));
}
// Values on edges.
......
......@@ -29,9 +29,11 @@
/// \file apps/mesh-segm-skel/mesh-complex-max-curv-skel.cc
/// \brief A program computing the maximal curvature values from the
/// surface of the (triangle) mesh of a statue, then computing a
/// skeleton of this surface (by/ thinning), using a complex-based
/// skeleton of this surface (by thinning), using a complex-based
/// image.
// FIXME: Factor with mesh-complex-pinv-curv-skel.cc.
#include <iostream>
#include <mln/core/image/complex_image.hh>
......@@ -92,14 +94,6 @@ main(int argc, char* argv[])
mln::p_n_faces_fwd_piter<D, G> v(float_ima.domain(), 0);
for_all(v)
{
// FIXME: The program should allow the user to choose the kind
// of measure (curvature).
#if 0
float h = (curv.first(v) + curv.second(v)) / 2;
// Pseudo-inverse curvature.
float h_inv = 1 / pi * (atan(-h) + pi / 2);
float_ima(v) = h_inv;
#endif
// Max curvature.
float_ima(v) = mln::math::max(mln::math::sqr(curv.first(v)),
mln::math::sqr(curv.second(v)));
......
// Copyright (C) 2008, 2009, 2010 EPITA Research and Development
// Laboratory (LRDE)
//
// This file is part of the Milena 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/mesh-segm-skel/mesh-complex-pinv-curv-segm.cc
/// \brief A program computing the pseudo-inverse curvature values
/// from the surface of the (triangle) mesh of a statue, then
/// performing a WST-based segmentation, using a complex-based image.
// Factor with mesh-complex-max-curv-segm.cc.
#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 float_ima_t;
// Dimension of the image (and therefore of the complex).
static const unsigned D = float_ima_t::dim;
// Geometry of the image.
typedef mln_geom_(float_ima_t) G;
mln::bin_2complex_image3df bin_input;
mln::io::off::load(bin_input, input_filename);
std::pair<float_ima_t, float_ima_t> curv =
mln::geom::mesh_curvature(bin_input.domain());
// Compute the pseudo_inverse curvature at each vertex.
float_ima_t float_ima(bin_input.domain());
mln::p_n_faces_fwd_piter<D, G> v(float_ima.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);
float_ima(v) = h_inv;
}
// Values on edges.
mln::p_n_faces_fwd_piter<D, G> e(float_ima.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 += float_ima(adj_v);
++n;
}
float_ima(e) = s / n;
// An edge should be adjacent to exactly two vertices.
mln_invariant(n <= 2);
}
// Convert the float image into an unsigned image because some
// algorithms do not handle float images well.
/* FIXME: This is bad: float images should be handled as-is. At
least, we should use a decent conversion, using an optimal affine
transformation (stretching/shrinking) or any other kind of
interpolation. */
typedef mln::unsigned_2complex_image3df ima_t;
ima_t ima(float_ima.domain());
mln_piter_(ima_t) p(float_ima.domain());
for_all(p)
ima(p) = 1000 * float_ima(p);
/*-----------------.
| Simplification. |
`-----------------*/
/// Adjacent 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_ima = mln::morpho::closing::area(ima, adj_edges_nbh, lambda);
/*------.
| 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_ima | mln::p_faces<1, D, G>(closed_ima.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_ima, 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::data::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);
}
// Copyright (C) 2008, 2009, 2010 EPITA Research and Development
// Laboratory (LRDE)
//
// This file is part of the Milena 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/mesh-segm-skel/mesh-complex-pinv-curv-skel.cc
/// \brief A program computing the pseudo-inverse curvature values
/// from the surface of the (triangle) mesh of a statue, then
/// computing a skeleton of this surface (by thinning), using a
/// complex-based image.
// FIXME: Factor with mesh-complex-max-curv-skel.cc.
#include <iostream>
#include <mln/core/image/complex_image.hh>
#include <mln/core/image/complex_neighborhoods.hh>
#include <mln/core/site_set/p_set.hh>
#include <mln/value/label_16.hh>
#include <mln/labeling/regional_minima.hh>
#include <mln/morpho/closing/area.hh>
#include <mln/topo/is_n_face.hh>
#include <mln/topo/is_simple_cell.hh>
#include <mln/topo/detach.hh>
#include <mln/topo/skeleton/breadth_first_thinning.hh>
#include <mln/io/off/load.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"
#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. |
`----------------*/
// Curvature image type.
typedef mln::float_2complex_image3df float_ima_t;
// Dimension of the image (and therefore of the complex).
static const unsigned D = float_ima_t::dim;
// Geometry of the image.
typedef mln_geom_(float_ima_t) G;
mln::bin_2complex_image3df bin_input;
mln::io::off::load(bin_input, input_filename);
std::pair<float_ima_t, float_ima_t> curv =
mln::geom::mesh_curvature(bin_input.domain());
// Compute the pseudo_inverse curvature at each vertex.
float_ima_t float_ima(bin_input.domain());
mln::p_n_faces_fwd_piter<D, G> v(float_ima.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);
float_ima(v) = h_inv;
}
// Values on triangles.
mln::p_n_faces_fwd_piter<D, G> t(float_ima.domain(), 2);
// For each triangle (2-face) T, iterate on the the set of vertices
// (0-faces) transitively adjacent to T.
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: We should be able to pass this value (m) either at the
construction of the neighborhood or at the construction of the
iterator. */
adj_v.iter().set_m(0);
// 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 += float_ima(adj_v);
++n;
}
float_ima(t) = s / n;
// A triangle should be adjacent to exactly two vertices.
mln_invariant(n <= 3);
}
// Convert the float image into an unsigned image because some
// algorithms do not handle float images well.
/* FIXME: This is bad: float images should be handled as-is. At
least, we should use a decent conversion, using an optimal affine
transformation (stretching/shrinking) or any other kind of
interpolation. */
typedef mln::unsigned_2complex_image3df ima_t;
ima_t ima(float_ima.domain());
// Process only triangles, as edges and vertices are set afterwards
// (see FIXME below).
for_all(t)
ima(t) = 1000 * float_ima(t);
/* FIXME: Workaround: Set maximal values on vertices and edges to
exclude them from the set of minimal values. */
for_all(v)
ima(v) = mln_max(mln_value_(ima_t));
mln::p_n_faces_fwd_piter<D, G> e(float_ima.domain(), 1);
for_all(e)
ima(e) = mln_max(mln_value_(ima_t));
/*-----------------.
| Simplification. |
`-----------------*/
/// Adjacent triangles are connected by shared edges.
typedef mln::complex_lower_dim_connected_n_face_neighborhood<D, G> nbh_t;
nbh_t nbh;
ima_t closed_ima = mln::morpho::closing::area(ima, nbh, lambda);
/*---------------.
| Local minima. |
`---------------*/
typedef mln::value::label_16 label_t;
label_t nminima;
/* 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. */
typedef mln_ch_value_(ima_t, label_t) label_ima_t;
label_ima_t minima =
mln::labeling::regional_minima(closed_ima, nbh, nminima);
typedef mln::complex_higher_neighborhood<D, G> higher_nbh_t;
higher_nbh_t higher_nbh;
// Propagate minima values from triangles to edges.
// FIXME: Factor this inside a function.
mln_niter_(higher_nbh_t) adj_t(higher_nbh, e);
for_all(e)
{
label_t ref_adj_minimum = mln::literal::zero;
for_all(adj_t)
if (minima(adj_t) == mln::literal::zero)
{
// If E is adjacent to a non-minimal triangle, then it must
// not belong to a minima.
ref_adj_minimum = mln::literal::zero;
break;
}
else
{
if (ref_adj_minimum == mln::literal::zero)
// If this is the first minimum seen, use it as a reference.
ref_adj_minimum = minima(adj_t);
else
// If this is not the first time a minimum is encountered,
// ensure it is REF_ADJ_MINIMUM.
mln_assertion(minima(adj_t) == ref_adj_minimum);
}
minima(e) = ref_adj_minimum;
}
// Likewise from edges to edges to vertices.
mln_niter_(higher_nbh_t) adj_e(higher_nbh, v);
for_all(v)
{
label_t ref_adj_minimum = mln::literal::zero;
for_all(adj_e)
if (minima(adj_e) == mln::literal::zero)
{
// If V is adjacent to a non-minimal triangle, then it must
// not belong to a minima.
ref_adj_minimum = mln::literal::zero;
break;
}
else
{
if (ref_adj_minimum == mln::literal::zero)
// If this is the first minimum seen, use it as a reference.
ref_adj_minimum = minima(adj_e);
else
// If this is not the first time a minimum is encountered,
// ensure it is REF_ADJ_MINIMUM.
mln_assertion(minima(adj_e) == ref_adj_minimum);
}
minima(v) = ref_adj_minimum;
}
/*-----------------------.
| Initial binary image. |
`-----------------------*/
/* Careful: creating ``holes'' in the surface obviously changes its
topology, but it may also split a single connected component in
two or more components, resulting in a disconnected skeleton. We
may want to improve this step either by forbidding any splitting,
or by incrementally ``digging'' a regional minima as long as no
splitting occurs. */
typedef mln_ch_value_(ima_t, bool) bin_ima_t;
bin_ima_t surface(minima.domain());
mln::data::fill(surface, true);
// Dig ``holes'' in the surface surface by setting minima faces to false.
// FIXME: Use fill with an image_if instead, when available/working.
mln_piter_(bin_ima_t) f(minima.domain());
for_all(f)
if (minima(f) != mln::literal::zero)
surface(f) = false;
/*-----------.
| Skeleton. |
`-----------*/
mln::topo::is_simple_cell<bin_ima_t> is_simple_p;
/* FIXME: Cheat! We'd like to iterate on cells of highest
dimension (2-cells) only, but we cannot constrain the domain of
the input using image_if (yet) like this
breadth_first_thinning(surface | is_n_face<2>, nbh, is_simple_p);
As a workaround, we use the constraint predicate of the
skeleton routine to restrict the iteration to 2-cells. */
mln::topo::is_n_face<bin_ima_t::dim> constraint_p;
bin_ima_t skel =
mln::topo::skeleton::breadth_first_thinning(surface, nbh,
is_simple_p,
mln::topo::detach<D, G>,
constraint_p);
/*---------.
| Output. |
`---------*/