Parallel II Estimator

cmake/CheckDGtalOptionalDependencies.cmake

@@ -230,6 +230,8 @@ endif()
 set(OPENMP_FOUND_DGTAL 0)
 if(DGTAL_WITH_OPENMP)
   include(openmp)
+  target_link_libraries(DGtal PUBLIC OpenMP::OpenMP_CXX)
+  target_compile_definitions(DGtal PUBLIC -DDGTAL_WITH_OPENMP)
   set(DGtalLibDependencies ${DGtalLibDependencies} OpenMP::OpenMP_CXX)
   set(OPENMP_FOUND_DGTAL 1)
   set(DGTAL_WITH_OPENMP 1)

src/DGtal/geometry/surfaces/estimation/DomainSplitter.h

@@ -0,0 +1,135 @@
+#pragma once
+
+#include "DGtal/kernel/domains/CDomain.h"
+#include <vector>
+
+namespace DGtal
+{
+  /**
+   * @brief Data structure returned by Domain splitters
+   *
+   * This class offers the ability for the splitter to give information
+   * about the surfaces or object within the splits.
+   * For example, the hintVoxelCount field may allow to pre-allocate buffers.
+   */
+  template<typename Domain>
+  struct SplitInfo
+  {
+    BOOST_CONCEPT_ASSERT(( concepts::CDomain< Domain > ));
+
+    Domain domain; //< The actual split domain
+    uint32_t hintVoxelCount = 0; //< An expected guess for the number of voxel
+  };
+
+  /**
+   * @brief Splits a domain evenly along all dimensions
+   *
+   * @tparam The model of domain to split
+   */
+  template<typename Domain>
+  struct RegularDomainSplitter
+  {
+    BOOST_CONCEPT_ASSERT(( concepts::CDomain< Domain > ));
+
+    //Output spllitted domain type
+    typedef std::vector<SplitInfo<Domain>> SplitDomainsInfo;
+
+    /**
+     * @brief Splits a domain
+     *
+     * This functions may split the domain in fewer part than wanted to ensure
+     * even subdomains.
+     *
+     * @param d The domain to split
+     * @param splitHint The targeted number of splits
+     */
+    SplitDomainsInfo operator()(const Domain& d, uint32_t splitHint) const
+    {
+      // Find best match possible for even splitting
+      const uint32_t splitCount = std::floor(std::log(splitHint) / std::log(Domain::dimension));
+      const uint32_t totalSplits = std::pow(splitCount, Domain::dimension);
+
+      if (splitCount == 0)
+        return { SplitInfo{d, 0} };
+
+      auto splitSize = (d.upperBound() - d.lowerBound()) / (int32_t)splitCount;
+      SplitDomainsInfo result;
+      result.reserve(totalSplits);
+
+      for (uint32_t i = 0; i < totalSplits; ++i)
+      {
+        auto start = d.lowerBound();
+        auto idx = i;
+        for (uint32_t j = 0; j < Domain::dimension; ++j)
+        {
+          auto k = idx % splitCount;
+          start[j] += k * splitSize[j] + k; // +k ensure no overlap between domains
+          idx /= splitCount;
+        }
+
+        // Make correction to ensure it remains within the domain
+        auto end = start + splitSize;
+        for (uint32_t j = 0; j < Domain::dimension; ++j)
+          end[j] = std::clamp(end[j], d.lowerBound()[j], d.upperBound()[j]);
+
+        result.emplace_back(Domain(start, end), 0);
+      }
+
+      return result;
+    };
+  };
+
+
+  /**
+   * @brief Splits a domain along one of the domain grid axis.
+   *
+   * @tparam The model of domain to split
+   */
+  template<typename Domain>
+  struct AxisDomainSplitter
+  {
+
+    BOOST_CONCEPT_ASSERT(( concepts::CDomain< Domain > ));
+
+    //Output spllitted domain type
+    typedef std::vector<SplitInfo<Domain>> SplitDomainsInfo;
+
+    /**
+     * @brief Regularly splits a domain along one axis
+     *
+     * @param d The domain to split
+     * @param splitHint The targeted number of splits (clamped to the width of the domain)
+     * @param dim the split axis
+     */
+    SplitDomainsInfo operator()(const Domain& d, uint32_t splitHint, typename Domain::Dimension dim) const
+    {
+      SplitDomainsInfo result;
+      if (splitHint == 0)
+        return result;
+
+      auto lower = d.lowerBound();
+      auto upper = d.upperBound();
+      auto length = upper[dim] - lower[dim] + 1;
+      uint32_t splitCount = splitHint;
+      if (splitCount > length)
+        splitCount = length;
+
+      result.reserve(splitCount);
+      auto base = length / splitCount;
+      auto rem = length % splitCount;
+
+      auto start = lower;
+      for (uint32_t i = 0; i < splitCount; ++i)
+      {
+        auto size = base + (i < rem ? 1 : 0);
+        auto end = upper;
+        end[dim] = start[dim] + size - 1;
+        result.emplace_back(Domain(start, end), 0);
+        start[dim] = end[dim] + 1;
+      }
+
+      return result;
+    };
+  };
+
+}

src/DGtal/geometry/surfaces/estimation/ParallelIIEstimator.h

@@ -0,0 +1,176 @@
+/**
+ *  This program is free software: you can redistribute it and/or modify
+ *  it under the terms of the GNU Lesser General Public License as
+ *  published by the Free Software Foundation, either version 3 of the
+ *  License, or  (at your option) any later version.
+ *
+ *  This program 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 program.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ **/
+
+#include <DGtal/geometry/surfaces/estimation/DomainSplitter.h>
+#include <DGtal/topology/helpers/Surfaces.h>
+
+#ifndef DGTAL_WITH_OPENMP
+#error You need to have activated OpenMP (DGTAL_WITH_OPENMP) to include this file.
+#endif
+
+// We require openmp for this estimator, even though it could run without
+#include <omp.h>
+
+namespace DGtal
+{
+  /**
+   * @brief Run an Integral Invariant estimator in parallel
+   *
+   * This class is meant as an almost perfect replacement of
+   * other IIEstimator. The only difference is the constructor
+   * which needs the number of threads.
+   *
+   * TODO: Value to control output (ie. value + loc, just value or ordered value)
+   * TODO: Surface construction on subdomain fonctor ?
+   *
+   * @tparam TEstimator The model of Estimator
+   * @tparam TSplitter The function to split the domain
+   */
+  template<class TEstimator, typename TSplitter>
+  class ParallelIIEstimator
+  {
+  public:
+    using Estimator = TEstimator;
+    using Splitter = TSplitter;
+    using Domain = typename TEstimator::Domain;
+    using Scalar = typename TEstimator::Scalar;
+
+    using KSpace = typename TEstimator::KSpace;
+    using PointPredicate = typename TEstimator::PointPredicate;
+    using Surfel = typename KSpace::Surfel;
+    using SurfelSet = typename KSpace::SurfelSet;
+    using EstimatorQuantity = typename TEstimator::Quantity;
+    using Quantity = EstimatorQuantity;
+
+    // Building the surface
+    using Boundary = LightImplicitDigitalSurface<KSpace, PointPredicate>;
+    using Surface = DigitalSurface<Boundary>;
+    using Visitor = DepthFirstVisitor<Surface>;
+    using VisitorRange = GraphVisitorRange<Visitor>;
+
+
+    /**
+     * @brief Constructor
+     *
+     * @tparam Args The estimator constructor arguments
+     *
+     * @param nbThread The number of thread to run in parallel. -1 means as many as possible
+     * @param args Constructor arguments to underlying estimators
+     */
+    template<typename... Args>
+    ParallelIIEstimator(int32_t nbThread, Args&&... args);
+
+    /**
+     * Clears the object. It is now invalid.
+     */
+    void clear();
+
+    /// @return the grid step.
+    Scalar h() const;
+
+    /**
+     * Attach a shape, defined as a functor spel -> boolean
+     *
+     * @param[in] K the cellular grid space in which the shape is defined.
+     * @param aPointPredicate the shape of interest. The alias can be secured
+     * if a some counted pointer is handed.
+    */
+    void attach(ConstAlias<KSpace> K,
+                ConstAlias<PointPredicate> pp);
+
+    /**
+     * Set specific parameters: the radius of the ball.
+     *
+     * @param[in] dRadius the "digital" radius of the kernel (buy may be non integer).
+     */
+    void setParams(double dRadius);
+
+
+    /**
+     * Checks the validity/consistency of the object.
+     * @return 'true' if the object is valid, 'false' otherwise.
+     */
+     bool isValid() const;
+
+    /**
+     * Model of CDigitalSurfaceLocalEstimator. Initialisation.
+     *
+     * @tparam SurfelConstIterator any model of forward readable iterator on Surfel.
+     * @param[in] _h grid size (must be >0).
+     * @param[in] ite iterator on the first surfel of the surface.
+     * @param[in] itb iterator after the last surfel of the surface.
+     */
+    template<typename ItA, typename ItB>
+    void init(double h_, ItA, ItB);
+
+    /**
+     * -- Estimation --
+     *
+     * Compute the integral invariant volume at surfel *it of
+     * a shape, then apply the VolumeFunctor to extract some
+     * geometric information.
+     *
+     * @tparam SurfelConstIterator type of Iterator on a Surfel
+     *
+     * @param[in] it iterator pointing on the surfel of the shape where
+     * we wish to evaluate some geometric information.
+     *
+     * @return quantity (normal vector) at surfel *it
+     */
+    template<typename It>
+    Quantity eval(It it);
+    /**
+     * -- Estimation --
+     *
+     * Compute the integral invariant volume for a range of
+     * surfels [itb,ite) on a shape, then apply the
+     * VolumeFunctor to extract some geometric information.
+     * Return the result on an OutputIterator (param).
+     *
+     * @tparam OutputIterator type of Iterator of an array of Quantity
+     * @tparam SurfelConstIterator type of Iterator on a Surfel
+     *
+     * @param[in] itb iterator defining the start of the range of surfels
+     * where we wish to compute some geometric information.
+     *
+     * @param[in] ite iterator defining the end of the range of surfels
+     * where we wish to compute some geometric information.
+     *
+     * @param[in] result output iterator of results of the computation.
+     * @return the updated output iterator after all outputs.
+     */
+    template<typename It, typename Oit>
+    Oit eval(It itb, It ite, Oit rslt);
+
+    /**
+     * Writes/Displays the object on an output stream.
+     * @param out the output stream where the object is written.
+     */
+    void selfDisplay ( std::ostream & out ) const;
+
+  private:
+    std::vector<Estimator> myEstimators;
+    Splitter mySplitter;
+
+    CountedConstPtrOrConstPtr<PointPredicate> myPointPredicate;
+    CountedConstPtrOrConstPtr<KSpace> myKSpace;
+
+    double myH;
+    double myRadius;
+  };
+}
+
+#include "ParallelIIEstimator.ih"