add postprocess function to the API (code related changes)

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] =
+        sdf_values[facet_begin] = -1.0;
-          0.0;  // no intersection is found, for now we are returning 0 not a special indicator
       }
     }
   }

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 Polyhedron>
+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<Polyhedron>           Default_filter;
+
+public:
+
+  template <class Filter, class SDFPropertyMap>
+  std::pair<double, double>
+  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 <class SDFPropertyMap>
+  std::pair<double, double>
+  postprocess(const Polyhedron& mesh, SDFPropertyMap sdf_pmap) {
+    return postprocess<Default_filter>(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<class SDFPropertyMap>
+  void check_missing_sdf_values(const Polyhedron& mesh,
+                                SDFPropertyMap sdf_values) {
+    std::vector<Facet_const_handle> still_missing_facets;
+    double min_sdf = (std::numeric_limits<double>::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<Facet_const_handle>::iterator it =
+          still_missing_facets.begin();
+        it != still_missing_facets.end(); ++it) {
+      sdf_values[*it] = min_sdf;
+    }
+  }
+
+  template<class SDFPropertyMap>
+  std::pair<double, double> min_max_value(const Polyhedron& mesh,
+                                          SDFPropertyMap sdf_values) {
+    double max_sdf = -(std::numeric_limits<double>::max)();
+    double min_sdf = (std::numeric_limits<double>::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<class SDFPropertyMap>
+  std::pair<double, double> 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<int>(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 <class SDFPropertyMap>
   std::pair<double, double>
   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<double, double> min_max_sdf_values = linear_normalize_sdf_values(
-          sdf_pmap);
 
-    // return minimum and maximum sdf values before normalization
+    Postprocess_sdf_values<Polyhedron> p;
-    return min_max_sdf_values;
+    return postprocess_req ? p.template postprocess<Filter>(mesh,
+           sdf_pmap) : // return minimum and maximum sdf values before normalization
+           p.min_max_value(mesh, sdf_pmap);
   }
 
   template <class FacetSegmentMap, class SDFPropertyMap>
@@ -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<class SDFPropertyMap>
-  std::pair<double, double> linear_normalize_sdf_values(SDFPropertyMap
-      sdf_values) {
-    double max_sdf = -(std::numeric_limits<double>::max)();
-    double min_sdf = (std::numeric_limits<double>::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<int>(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<class SDFPropertyMap>
-  void check_zero_sdf_values(SDFPropertyMap sdf_values) {
-    std::vector<Facet_const_handle> still_zero_facets;
-    double min_sdf = (std::numeric_limits<double>::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<Facet_const_handle>::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.

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:
+ * It is possible to compute raw SDF values (see \ref Surface_mesh_segmentationRawSDF) and apply post-processing steps (see \ref Surface_mesh_segmentationPostprocessing).
- *   - Facets with no SDF values (i.e. zero) are assigned the average SDF value of their one-ring edge-adjacent neighbors.
+ * For raw SDF values, -1.0 is used as an indicator for no SDF value.
- *     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].
  *
  * @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 <bool Fast_sdf_calculation_mode, class Polyhedron,
          class SDFPropertyMap, class GeomTraits
@@ -48,11 +45,13 @@ 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())
 {
   internal::Surface_mesh_segmentation<Polyhedron, GeomTraits, Fast_sdf_calculation_mode>
   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<true, Polyhedron, SDFPropertyMap, GeomTraits>
-         (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<class Polyhedron, class SDFPropertyMap>
+std::pair<double, double>
+postprocess_sdf_values(const Polyhedron& polyhedron, SDFPropertyMap sdf_values)
+{
+  CGAL_precondition(polyhedron.is_pure_triangle());
+  return internal::Postprocess_sdf_values<Polyhedron>().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<Fast_sdf_calculation_mode, Polyhedron, SDFPropertyMap, typename Polyhedron::Traits>
-         (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<true, Polyhedron, SDFPropertyMap, typename Polyhedron::Traits>
-         (polyhedron, sdf_values, cone_angle, number_of_rays, traits);
+         (polyhedron, sdf_values, cone_angle, number_of_rays, postprocess, traits);
 }
 
 template <class Polyhedron, class SDFPropertyMap, class SegmentPropertyMap>