From dd5bf0ec3ed26ecef2d853adfdffa37a7218141d Mon Sep 17 00:00:00 2001 From: iyaz Date: Thu, 8 Aug 2013 02:33:20 +0300 Subject: [PATCH] add postprocess function to the API (code related changes) --- .../SDF_calculation.h | 3 +- .../Surface_mesh_segmentation.h | 251 ++++++++++-------- .../include/CGAL/mesh_segmentation.h | 57 +++- 3 files changed, 186 insertions(+), 125 deletions(-) diff --git a/Surface_mesh_segmentation/include/CGAL/internal/Surface_mesh_segmentation/SDF_calculation.h b/Surface_mesh_segmentation/include/CGAL/internal/Surface_mesh_segmentation/SDF_calculation.h index 76ee5eedba7..79959d5ebe0 100644 --- a/Surface_mesh_segmentation/include/CGAL/internal/Surface_mesh_segmentation/SDF_calculation.h +++ b/Surface_mesh_segmentation/include/CGAL/internal/Surface_mesh_segmentation/SDF_calculation.h @@ -197,8 +197,7 @@ public: if(sdf_value) { sdf_values[facet_begin] = *sdf_value; } else { - sdf_values[facet_begin] = - 0.0; // no intersection is found, for now we are returning 0 not a special indicator + sdf_values[facet_begin] = -1.0; } } } diff --git a/Surface_mesh_segmentation/include/CGAL/internal/Surface_mesh_segmentation/Surface_mesh_segmentation.h b/Surface_mesh_segmentation/include/CGAL/internal/Surface_mesh_segmentation/Surface_mesh_segmentation.h index 4c606598589..1839d347b58 100644 --- a/Surface_mesh_segmentation/include/CGAL/internal/Surface_mesh_segmentation/Surface_mesh_segmentation.h +++ b/Surface_mesh_segmentation/include/CGAL/internal/Surface_mesh_segmentation/Surface_mesh_segmentation.h @@ -25,6 +25,139 @@ namespace CGAL /// @cond CGAL_DOCUMENT_INTERNAL namespace internal { + +// Postprocess functions for sdf values +template +class Postprocess_sdf_values +{ + + typedef typename Polyhedron::Facet Facet; + typedef typename Polyhedron::Facet_const_handle Facet_const_handle; + typedef typename Polyhedron::Facet_const_iterator Facet_const_iterator; + + typedef Bilateral_filtering Default_filter; + +public: + + template + std::pair + postprocess(const Polyhedron& mesh, SDFPropertyMap sdf_pmap) { + check_missing_sdf_values(mesh, sdf_pmap); + Filter()(mesh, get_window_size(mesh), sdf_pmap); + return linear_normalize_sdf_values(mesh, sdf_pmap); + } + + template + std::pair + postprocess(const Polyhedron& mesh, SDFPropertyMap sdf_pmap) { + return postprocess(mesh, sdf_pmap); + } + + /** + * Finds facets which have missing (indicated by -1) sdf values. + * Sdf values on these facets are assigned to average sdf value of its neighbors. + * If still there is any facet which has no sdf value, assigns minimum sdf value to it. + * This is meaningful since (being an outlier) zero sdf values might effect normalization & log extremely. + * @param[in, out] sdf_values `ReadWritePropertyMap` with `Polyhedron::Facet_const_handle` as key and `double` as value type + */ + template + void check_missing_sdf_values(const Polyhedron& mesh, + SDFPropertyMap sdf_values) { + std::vector still_missing_facets; + double min_sdf = (std::numeric_limits::max)(); + // If there is any facet which has no sdf value, assign average sdf value of its neighbors + for(Facet_const_iterator facet_it = mesh.facets_begin(); + facet_it != mesh.facets_end(); ++facet_it) { + double sdf_value = sdf_values[facet_it]; + CGAL_assertion(sdf_value == -1 || sdf_value >= 0); // validity check + if(sdf_value != -1.0) { + min_sdf = (std::min)(sdf_value, min_sdf); + continue; + } + + typename Facet::Halfedge_around_facet_const_circulator facet_circulator = + facet_it->facet_begin(); + double total_neighbor_sdf = 0.0; + int nb_valid_neighbors = 0; + do { + if(!facet_circulator->opposite()->is_border()) { + double neighbor_sdf = sdf_values[facet_circulator->opposite()->facet()]; + if(neighbor_sdf != -1) { + total_neighbor_sdf += neighbor_sdf; + ++nb_valid_neighbors; + } + } + } while( ++facet_circulator != facet_it->facet_begin()); + + if(nb_valid_neighbors == 0) { + still_missing_facets.push_back(facet_it); + } else { + sdf_value = total_neighbor_sdf / nb_valid_neighbors; + sdf_values[facet_it] = sdf_value; + // trying to update min_sdf is pointless, since it is interpolated one of the neighbors sdf will be smaller than it + } + } + // If still there is any facet which has no sdf value, assign minimum sdf value. + // This is meaningful since (being an outlier) 0 sdf values might effect normalization & log extremely. + for(typename std::vector::iterator it = + still_missing_facets.begin(); + it != still_missing_facets.end(); ++it) { + sdf_values[*it] = min_sdf; + } + } + + template + std::pair min_max_value(const Polyhedron& mesh, + SDFPropertyMap sdf_values) { + double max_sdf = -(std::numeric_limits::max)(); + double min_sdf = (std::numeric_limits::max)(); + for(Facet_const_iterator facet_it = mesh.facets_begin(); + facet_it != mesh.facets_end(); ++facet_it) { + double sdf_value = sdf_values[facet_it]; + max_sdf = (std::max)(sdf_value, max_sdf); + min_sdf = (std::min)(sdf_value, min_sdf); + } + return std::make_pair(min_sdf, max_sdf); + } + /** + * Normalize sdf values between [0-1]. + * @param sdf_values `ReadWritePropertyMap` with `Polyhedron::Facet_const_handle` as key and `double` as value type + * @return minimum and maximum SDF values before normalization + */ + template + std::pair linear_normalize_sdf_values(const Polyhedron& mesh, + SDFPropertyMap sdf_values) { + double min_sdf, max_sdf; + boost::tie(min_sdf, max_sdf) = min_max_value(mesh, sdf_values); + + if(min_sdf == max_sdf) { + CGAL_warning(min_sdf == max_sdf && !"Linear normalization is not applicable!"); + return std::make_pair(min_sdf, max_sdf); + } + + const double max_min_dif = max_sdf - min_sdf; + for(Facet_const_iterator facet_it = mesh.facets_begin(); + facet_it != mesh.facets_end(); ++facet_it) { + sdf_values[facet_it] = (sdf_values[facet_it] - min_sdf) / max_min_dif; + } + return std::make_pair(min_sdf, max_sdf); + } + + /** + * Simple window-size determination function for smoothing. + * It is proportional to square root of size of facets in polyhedron. + * @return size of the window + * - 0-2000 -> 1 + * - 2000-8000 -> 2 + * - 8000-18000 -> 3 + * - ... + */ + int get_window_size(const Polyhedron& mesh) { + double facet_sqrt = std::sqrt(mesh.size_of_facets() / 2000.0); + return static_cast(facet_sqrt) + 1; + } +}; + /** * @brief Main entry point for mesh segmentation algorithm. * @@ -80,7 +213,6 @@ private: private: const Polyhedron& mesh; GeomTraits traits; - int m_window_size; // member functions public: /** @@ -88,7 +220,7 @@ public: * @param mesh `CGAL Polyhedron` on which other functions operate. */ Surface_mesh_segmentation(const Polyhedron& mesh, GeomTraits traits) - : mesh(mesh), traits(traits), m_window_size( get_default_window_size() ) { + : mesh(mesh), traits(traits) { CGAL_precondition(mesh.is_pure_triangle()); } @@ -96,19 +228,16 @@ public: template std::pair calculate_sdf_values(double cone_angle, int number_of_rays, - SDFPropertyMap sdf_pmap) { + SDFPropertyMap sdf_pmap, bool postprocess_req) { // calculate sdf values SDF_calculation_class sdf_calculator(mesh, false, true, traits); sdf_calculator.calculate_sdf_values(mesh.facets_begin(), mesh.facets_end(), cone_angle, number_of_rays, sdf_pmap); - // apply post-processing steps - check_zero_sdf_values(sdf_pmap); - Filter()(mesh, get_window_size(), sdf_pmap); - std::pair min_max_sdf_values = linear_normalize_sdf_values( - sdf_pmap); - // return minimum and maximum sdf values before normalization - return min_max_sdf_values; + Postprocess_sdf_values p; + return postprocess_req ? p.template postprocess(mesh, + sdf_pmap) : // return minimum and maximum sdf values before normalization + p.min_max_value(mesh, sdf_pmap); } template @@ -157,21 +286,6 @@ public: return number_of_centers; } - /** - * Return the current window size. It is used when calling Filter::operator(). - * By default it is the value returned by `get_default_window_size()`. - */ - int get_window_size() const { - return m_window_size; - } - - /** - * Set the window size used to call `Filter::operator()` to `s` - */ - void set_window_size(int s) { - m_window_size=s; - } - private: /** @@ -220,93 +334,6 @@ private: } } - /** - * Normalize sdf values between [0-1]. - * @param sdf_values `ReadWritePropertyMap` with `Polyhedron::Facet_const_handle` as key and `double` as value type - * @return minimum and maximum SDF values before normalization - */ - template - std::pair linear_normalize_sdf_values(SDFPropertyMap - sdf_values) { - double max_sdf = -(std::numeric_limits::max)(); - double min_sdf = (std::numeric_limits::max)(); - for(Facet_const_iterator facet_it = mesh.facets_begin(); - facet_it != mesh.facets_end(); ++facet_it) { - double sdf_value = sdf_values[facet_it]; - max_sdf = (std::max)(sdf_value, max_sdf); - min_sdf = (std::min)(sdf_value, min_sdf); - } - const double max_min_dif = max_sdf - min_sdf; - for(Facet_const_iterator facet_it = mesh.facets_begin(); - facet_it != mesh.facets_end(); ++facet_it) { - sdf_values[facet_it] = (sdf_values[facet_it] - min_sdf) / max_min_dif; - } - return std::make_pair(min_sdf, max_sdf); - } - - - /** - * Simple window-size determination function for smoothing. - * It is proportional to square root of size of facets in polyhedron. - * @return size of the window - * - 0-2000 -> 1 - * - 2000-8000 -> 2 - * - 8000-18000 -> 3 - * - ... - */ - int get_default_window_size() { - double facet_sqrt = std::sqrt(mesh.size_of_facets() / 2000.0); - return static_cast(facet_sqrt) + 1; - } - - /** - * Finds facets which have zero sdf values. - * Sdf values on these facets are assigned to average sdf value of its neighbors. - * If still there is any facet which has no sdf value, assigns minimum sdf value to it. - * This is meaningful since (being an outlier) zero sdf values might effect normalization & log extremely. - * @param[in, out] sdf_values `ReadWritePropertyMap` with `Polyhedron::Facet_const_handle` as key and `double` as value type - */ - template - void check_zero_sdf_values(SDFPropertyMap sdf_values) { - std::vector still_zero_facets; - double min_sdf = (std::numeric_limits::max)(); - // If there is any facet which has no sdf value, assign average sdf value of its neighbors - for(Facet_const_iterator facet_it = mesh.facets_begin(); - facet_it != mesh.facets_end(); ++facet_it) { - double sdf_value = sdf_values[facet_it]; - if(sdf_value == 0.0) { - typename Facet::Halfedge_around_facet_const_circulator facet_circulator = - facet_it->facet_begin(); - double total_neighbor_sdf = 0.0; - do { - if(!facet_circulator->opposite()->is_border()) { - total_neighbor_sdf += sdf_values[facet_circulator->opposite()->facet()]; - } - } while( ++facet_circulator != facet_it->facet_begin()); - - sdf_value = total_neighbor_sdf / 3.0; - sdf_values[facet_it] = sdf_value; - - if(sdf_value == 0.0) { - still_zero_facets.push_back( - facet_it); // if sdf_value is still zero, put facet to zero-list - } else { - min_sdf = (std::min)(sdf_value, - min_sdf); // find min_sdf other than zero - } - } else { - min_sdf = (std::min)(sdf_value, min_sdf); // find min_sdf other than zero - } - } - // If still there is any facet which has no sdf value, assign minimum sdf value. - // This is meaningful since (being an outlier) 0 sdf values might effect normalization & log extremely. - for(typename std::vector::iterator it = - still_zero_facets.begin(); - it != still_zero_facets.end(); ++it) { - sdf_values[*it] = min_sdf; - } - } - /** * Receives probability-matrix with probabilities betwen [0-1], and returns log-normalized probabilities * which are suitable to use in graph-cut. diff --git a/Surface_mesh_segmentation/include/CGAL/mesh_segmentation.h b/Surface_mesh_segmentation/include/CGAL/mesh_segmentation.h index 037677d41fc..1460be0c963 100644 --- a/Surface_mesh_segmentation/include/CGAL/mesh_segmentation.h +++ b/Surface_mesh_segmentation/include/CGAL/mesh_segmentation.h @@ -17,12 +17,8 @@ namespace CGAL * @brief Function computing the Shape Diameter Function over a surface mesh. * * This function implements the Shape Diameter Function (SDF) as described in \cite shapira2008consistent. - * After the computation of SDF values for each facet, the following post-processing steps are applied: - * - Facets with no SDF values (i.e. zero) are assigned the average SDF value of their one-ring edge-adjacent neighbors. - * If there is still a facet having a zero SDF value, the minimum SDF value greater than zero is assigned to it. Note that this step is not inherited from the paper. - * The main reason for avoiding zero SDF values is that it can obstruct log-normalization process which takes place at the beginning of segment_from_sdf_values()`. - * - SDF values are smoothed with bilateral filtering. - * - SDF values are linearly normalized between [0,1]. + * It is possible to compute raw SDF values (see \ref Surface_mesh_segmentationRawSDF) and apply post-processing steps (see \ref Surface_mesh_segmentationPostprocessing). + * For raw SDF values, -1.0 is used as an indicator for no SDF value. * * @pre @a polyhedron.is_pure_triangle() * @tparam Fast_sdf_calculation_mode regardless of `GeomTraits`, use inexact predicates while traversing AABB tree nodes. @@ -34,8 +30,9 @@ namespace CGAL * @param[out] sdf_values the SDF value of each facet * @param cone_angle opening angle for cone, expressed in radians * @param number_of_rays number of rays picked from cone for each facet. In general, increasing the number of rays beyond the default value has little influence upon the resulting segmentation. + * @param postprocess if true apply post-processing steps in `CGAL::postprocess_sdf_values`, otherwise return raw SDF values. * @param traits traits object - * @return minimum and maximum SDF values before linear normalization + * @return minimum and maximum SDF values if @a postprocess is true, otherwise minimum and maximum SDF values before linear normalization */ template algorithm(polyhedron, traits); - return algorithm.calculate_sdf_values(cone_angle, number_of_rays, sdf_values); + return algorithm.calculate_sdf_values(cone_angle, number_of_rays, sdf_values, + postprocess); } /// @cond SKIP_IN_MANUAL @@ -66,13 +65,47 @@ compute_sdf_values(const Polyhedron& polyhedron, SDFPropertyMap sdf_values, double cone_angle = 2.0 / 3.0 * CGAL_PI, int number_of_rays = 25, + bool postprocess = true, GeomTraits traits = GeomTraits()) { return compute_sdf_values - (polyhedron, sdf_values, cone_angle, number_of_rays, traits); + (polyhedron, sdf_values, cone_angle, number_of_rays, postprocess, traits); } /// @endcond +/*! + * \ingroup PkgSurfaceSegmentation + * @brief Function post-processing raw SDF values computed per facet. + * + * Post-processing steps applied are: + * - Facets with -1.0 SDF values are assigned the average SDF value of their one-ring edge-adjacent neighbors. + * If there is still a facet having -1.0 SDF value, the minimum valid SDF value assigned to it. Note that this step is not inherited from the paper. + * The main reason for not assigning 0 to facets with no SDF values (i.e. -1.0) is that it can obstruct log-normalization process which takes place at the beginning of `CGAL::segment_from_sdf_values`. + * - SDF values are smoothed with bilateral filtering. + * - SDF values are linearly normalized between [0,1]. + * + * See the section \ref Surface_mesh_segmentationPostprocessing for more details. + * + * @pre @a polyhedron.is_pure_triangle() + * @pre Raw values should be either equal to -1.0 or, greater or equal than 0.0 + * @tparam Polyhedron a %CGAL polyhedron + * @tparam SDFPropertyMap a `ReadWritePropertyMap` with `Polyhedron::Facet_const_handle` as key and `double` as value type + * + * @param polyhedron surface mesh on which SDF values are computed + * @param[in, out] sdf_values the SDF value of each facet + * + * @return minimum and maximum SDF values before linear normalization + */ +template +std::pair +postprocess_sdf_values(const Polyhedron& polyhedron, SDFPropertyMap sdf_values) +{ + CGAL_precondition(polyhedron.is_pure_triangle()); + return internal::Postprocess_sdf_values().postprocess(polyhedron, + sdf_values); +} + + /*! * \ingroup PkgSurfaceSegmentation * @brief Function computing the segmentation of a surface mesh given an SDF value per facet. @@ -206,10 +239,11 @@ compute_sdf_values(const Polyhedron& polyhedron, SDFPropertyMap sdf_values, double cone_angle = 2.0 / 3.0 * CGAL_PI, int number_of_rays = 25, + bool postprocess = true, typename Polyhedron::Traits traits = typename Polyhedron::Traits()) { return compute_sdf_values - (polyhedron, sdf_values, cone_angle, number_of_rays, traits); + (polyhedron, sdf_values, cone_angle, number_of_rays, postprocess, traits); } template < class Polyhedron, class SDFPropertyMap> @@ -218,10 +252,11 @@ compute_sdf_values(const Polyhedron& polyhedron, SDFPropertyMap sdf_values, double cone_angle = 2.0 / 3.0 * CGAL_PI, int number_of_rays = 25, + bool postprocess = true, typename Polyhedron::Traits traits = typename Polyhedron::Traits()) { return compute_sdf_values - (polyhedron, sdf_values, cone_angle, number_of_rays, traits); + (polyhedron, sdf_values, cone_angle, number_of_rays, postprocess, traits); } template