* apps/mesh-segm-skel/mesh-complex-segm.cc, * apps/mesh-segm-skel/mesh-complex-skel.cc, * apps/mesh-segm-skel/mesh-complex-max-curv-segm.cc, * apps/mesh-segm-skel/mesh-complex-max-curv-skel.cc: Convert the floating point value input to an unsigned integer image to avoid precision errors due to float comparisons. --- milena/ChangeLog | 11 ++++ .../mesh-segm-skel/mesh-complex-max-curv-segm.cc | 39 ++++++++---- .../mesh-segm-skel/mesh-complex-max-curv-skel.cc | 61 +++++++++++++------- milena/apps/mesh-segm-skel/mesh-complex-segm.cc | 24 ++++++-- milena/apps/mesh-segm-skel/mesh-complex-skel.cc | 41 ++++++++++--- 5 files changed, 127 insertions(+), 49 deletions(-) diff --git a/milena/ChangeLog b/milena/ChangeLog index 3de4207..2ce43e0 100644 --- a/milena/ChangeLog +++ b/milena/ChangeLog @@ -1,5 +1,16 @@ 2010-06-14 Roland Levillain &lt;roland(a)lrde.epita.fr&gt; + Improve the precision of some mesh-segm-skel apps. + + * apps/mesh-segm-skel/mesh-complex-segm.cc, + * apps/mesh-segm-skel/mesh-complex-skel.cc, + * apps/mesh-segm-skel/mesh-complex-max-curv-segm.cc, + * apps/mesh-segm-skel/mesh-complex-max-curv-skel.cc: + Convert the floating point value input to an unsigned integer + image to avoid precision errors due to float comparisons. + +2010-06-14 Roland Levillain &lt;roland(a)lrde.epita.fr&gt; + New application: mesh-complex-max-curv-skel. * apps/mesh-segm-skel/mesh-complex-max-curv-skel.cc: New. diff --git a/milena/apps/mesh-segm-skel/mesh-complex-max-curv-segm.cc b/milena/apps/mesh-segm-skel/mesh-complex-max-curv-segm.cc index 477ccdf..1f43381 100644 --- a/milena/apps/mesh-segm-skel/mesh-complex-max-curv-segm.cc +++ b/milena/apps/mesh-segm-skel/mesh-complex-max-curv-segm.cc @@ -75,25 +75,26 @@ int main(int argc, char* argv[]) `----------------*/ // Image type. - typedef mln::float_2complex_image3df ima_t; + typedef mln::float_2complex_image3df float_ima_t; // Dimension of the image (and therefore of the complex). - static const unsigned D = ima_t::dim; + static const unsigned D = float_ima_t::dim; // Geometry of the image. - typedef mln_geom_(ima_t) G; + typedef mln_geom_(float_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()); + std::pair<float_ima_t, float_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); + 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); - input(v) = h_inv; + float_ima(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)), @@ -101,7 +102,7 @@ int main(int argc, char* argv[]) } // Values on edges. - mln::p_n_faces_fwd_piter<D, G> e(input.domain(), 1); + 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); @@ -113,14 +114,26 @@ int main(int argc, char* argv[]) // Iterate on vertices (0-faces). for_all(adj_v) { - s += input(adj_v); + s += float_ima(adj_v); ++n; } - input(e) = s / 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. | `-----------------*/ @@ -131,7 +144,7 @@ int main(int argc, char* argv[]) adj_edges_nbh_t; adj_edges_nbh_t adj_edges_nbh; - ima_t closed_input = mln::morpho::closing::area(input, adj_edges_nbh, lambda); + ima_t closed_ima = mln::morpho::closing::area(ima, adj_edges_nbh, lambda); /*------. | WST. | @@ -141,7 +154,7 @@ int main(int argc, char* argv[]) 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()) + closed_ima | mln::p_faces<1, D, G>(closed_ima.domain()) as input of the WST. */ @@ -150,7 +163,7 @@ int main(int argc, char* argv[]) 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); + 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). diff --git a/milena/apps/mesh-segm-skel/mesh-complex-max-curv-skel.cc b/milena/apps/mesh-segm-skel/mesh-complex-max-curv-skel.cc index 6265a25..a174dd6 100644 --- a/milena/apps/mesh-segm-skel/mesh-complex-max-curv-skel.cc +++ b/milena/apps/mesh-segm-skel/mesh-complex-max-curv-skel.cc @@ -75,37 +75,38 @@ main(int argc, char* argv[]) | Complex image. | `----------------*/ - // Image type. - typedef mln::float_2complex_image3df ima_t; + // Curvature image type. + typedef mln::float_2complex_image3df float_ima_t; // Dimension of the image (and therefore of the complex). - static const unsigned D = ima_t::dim; + static const unsigned D = float_ima_t::dim; // Geometry of the image. - typedef mln_geom_(ima_t) G; + typedef mln_geom_(float_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()); + std::pair<float_ima_t, float_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); + float_ima_t float_ima(bin_input.domain()); + 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. + // 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); - input(v) = h_inv; + float_ima(v) = h_inv; #endif // Max curvature. - input(v) = mln::math::max(mln::math::sqr(curv.first(v)), - mln::math::sqr(curv.second(v))); + float_ima(v) = mln::math::max(mln::math::sqr(curv.first(v)), + mln::math::sqr(curv.second(v))); } // Values on triangles. - mln::p_n_faces_fwd_piter<D, G> t(input.domain(), 2); + 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; @@ -123,21 +124,34 @@ main(int argc, char* argv[]) // Iterate on vertices (0-faces). for_all(adj_v) { - s += input(adj_v); + s += float_ima(adj_v); ++n; } - input(t) = s / 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) - input(v) = mln_max(float); - mln::p_n_faces_fwd_piter<D, G> e(input.domain(), 1); + ima(v) = mln_max(mln_value_(ima_t)); + mln::p_n_faces_fwd_piter<D, G> e(float_ima.domain(), 1); for_all(e) - input(e) = mln_max(float); + ima(e) = mln_max(mln_value_(ima_t)); /*-----------------. | Simplification. | @@ -147,7 +161,7 @@ main(int argc, char* argv[]) typedef mln::complex_lower_dim_connected_n_face_neighborhood<D, G> nbh_t; nbh_t nbh; - ima_t closed_input = mln::morpho::closing::area(input, nbh, lambda); + ima_t closed_ima = mln::morpho::closing::area(ima, nbh, lambda); /*---------------. | Local minima. | @@ -161,7 +175,7 @@ main(int argc, char* argv[]) to 2-faces. */ typedef mln_ch_value_(ima_t, label_t) label_ima_t; label_ima_t minima = - mln::labeling::regional_minima(closed_input, nbh, nminima); + mln::labeling::regional_minima(closed_ima, nbh, nminima); typedef mln::complex_higher_neighborhood<D, G> higher_nbh_t; higher_nbh_t higher_nbh; @@ -201,7 +215,7 @@ main(int argc, char* argv[]) for_all(adj_e) if (minima(adj_e) == mln::literal::zero) { - // If V is adjcent to a non-minimal triangle, then it must + // If V is adjacent to a non-minimal triangle, then it must // not belong to a minima. ref_adj_minimum = mln::literal::zero; break; @@ -223,6 +237,13 @@ main(int argc, char* argv[]) | 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); diff --git a/milena/apps/mesh-segm-skel/mesh-complex-segm.cc b/milena/apps/mesh-segm-skel/mesh-complex-segm.cc index 8d349fe..655a46c 100644 --- a/milena/apps/mesh-segm-skel/mesh-complex-segm.cc +++ b/milena/apps/mesh-segm-skel/mesh-complex-segm.cc @@ -66,15 +66,27 @@ int main(int argc, char* argv[]) | Complex image. | `----------------*/ - // Image type. - typedef mln::float_2complex_image3df ima_t; + // Input type. + typedef mln::float_2complex_image3df float_input_t; // Dimension of the image (and therefore of the complex). - static const unsigned D = ima_t::dim; + static const unsigned D = float_input_t::dim; // Geometry of the image. - typedef mln_geom_(ima_t) G; - - ima_t input; - mln::io::off::load(input, input_filename); + typedef mln_geom_(float_input_t) G; + + float_input_t float_input; + mln::io::off::load(float_input, input_filename); + + // 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 input(float_input.domain()); + mln_piter_(ima_t) p(input.domain()); + for_all(p) + input(p) = 1000 * float_input(p); // Values on edges. mln::p_n_faces_fwd_piter<D, G> e(input.domain(), 1); diff --git a/milena/apps/mesh-segm-skel/mesh-complex-skel.cc b/milena/apps/mesh-segm-skel/mesh-complex-skel.cc index b7e1f1f..5cd97b4 100644 --- a/milena/apps/mesh-segm-skel/mesh-complex-skel.cc +++ b/milena/apps/mesh-segm-skel/mesh-complex-skel.cc @@ -72,23 +72,37 @@ main(int argc, char* argv[]) `----------------*/ // Image type. - typedef mln::float_2complex_image3df ima_t; + typedef mln::float_2complex_image3df float_input_t; // Dimension of the image (and therefore of the complex). - static const unsigned D = ima_t::dim; + static const unsigned D = float_input_t::dim; // Geometry of the image. - typedef mln_geom_(ima_t) G; - - ima_t input; - mln::io::off::load(input, input_filename); + typedef mln_geom_(float_input_t) G; + + float_input_t float_input; + mln::io::off::load(float_input, input_filename); + + // 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 input(float_input.domain()); + // Process only triangles, as edges and vertices are set afterwards + // (see FIXME below). + mln::p_n_faces_fwd_piter<D, G> t(input.domain(), 2); + for_all(t) + input(t) = 1000 * float_input(t); /* 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); + input(v) = mln_max(mln_value_(ima_t)); mln::p_n_faces_fwd_piter<D, G> e(input.domain(), 1); for_all(e) - input(e) = mln_max(float); + input(e) = mln_max(mln_value_(ima_t)); /*-----------------. | Simplification. | @@ -126,7 +140,7 @@ main(int argc, char* argv[]) for_all(adj_t) if (minima(adj_t) == mln::literal::zero) { - // If E is adjcent to a non-minimal triangle, then it must + // If E is adjacent to a non-minimal triangle, then it must // not belong to a minima. ref_adj_minimum = mln::literal::zero; break; @@ -152,7 +166,7 @@ main(int argc, char* argv[]) for_all(adj_e) if (minima(adj_e) == mln::literal::zero) { - // If V is adjcent to a non-minimal triangle, then it must + // If V is adjacent to a non-minimal triangle, then it must // not belong to a minima. ref_adj_minimum = mln::literal::zero; break; @@ -174,6 +188,13 @@ main(int argc, char* argv[]) | 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); -- 1.5.6.5