Parallel II Estimator

ChangeLog.md

@@ -29,6 +29,9 @@
 - *Build*
   - Fixing old `DEBUG_VERBOSE` to `DGTAL_DEBUG_VERBOSE` (David Coeurjolly,  [#1834](https://github.com/DGtal-team/DGtal/pull/1834))
 
+- *Geometry*
+  - In DigitalSurfaceConvolver (eg. when using the Integral Invariant estimators), an bug existed when the shapoe domain does not contain the point (0,0,0). This PR fixes it using KCoords to represent the Kernel origin instead of a Spel which requires a valid KSpace (David Coeurjolly, [#1840](https://github.com/DGtal-team/DGtal/pull/1840))
+
 # DGtal 2.1.0
 
 ## New features

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/DigitalSurfaceConvolver.h

@@ -408,7 +408,7 @@ class DigitalSurfaceConvolver
   const DigitalKernel * myKernel; ///< Two choice to iterate over the full kernel. See init() for more information.
   const PairIterators * myKernelMask; ///< Two choice to iterate over the full kernel. See init() for more information.
 
-  Spel myKernelSpelOrigin; ///< Copy of the origin cell of the kernel.
+  KSpace::Point myKernelKCoordsOrigin; ///< Copy of the origin cell of the kernel.
 
   // ------------------------- Hidden services ------------------------------
 
@@ -775,7 +775,7 @@ class DigitalSurfaceConvolver< TFunctor, TKernelFunctor, TKSpace, TDigitalKernel
   const DigitalKernel * myKernel; ///< Two choice to iterate over the full kernel. See init() for more information.
   const PairIterators * myKernelMask; ///< Two choice to iterate over the full kernel. See init() for more information.
 
-  Spel myKernelSpelOrigin; ///< Copy of the origin cell of the kernel.
+  KSpace::Point myKernelKCoordsOrigin; ///< Copy of the origin cell of the kernel.
 
   // ------------------------- Hidden services ------------------------------
 
@@ -1145,7 +1145,7 @@ class DigitalSurfaceConvolver< TFunctor, TKernelFunctor, TKSpace, TDigitalKernel
   const DigitalKernel * myKernel; ///< Two choice to iterate over the full kernel. See init() for more information.
   const PairIterators * myKernelMask; ///< Two choice to iterate over the full kernel. See init() for more information.
 
-  Spel myKernelSpelOrigin; ///< Copy of the origin cell of the kernel.
+  KSpace::Point myKernelKCoordsOrigin; ///< Copy of the origin cell of the kernel.
 
   // ------------------------- Hidden services ------------------------------
 

src/DGtal/geometry/surfaces/DigitalSurfaceConvolver.ih

@@ -77,7 +77,7 @@ DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, d
     myMasks( other.myMasks ),
     myKernel( other.myKernel ),
     myKernelMask( other.myKernelMask ),
-    myKernelSpelOrigin( other.myKernelSpelOrigin )
+    myKernelKCoordsOrigin( other.myKernelKCoordsOrigin )
 {
 }
 
@@ -89,7 +89,8 @@ DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, d
   ConstAlias< PairIterators > fullKernel,
   ConstAlias< std::vector< PairIterators > > masks )
 {
-  myKernelSpelOrigin = myKSpace.sSpel( pOrigin );
+  KhalimskySpaceND<dimension, typename KSpace::Integer> preK;
+  myKernelKCoordsOrigin = preK.sKCoords(preK.sSpel( pOrigin ));
   myKernelMask = &fullKernel;
   myMasks = &masks;
 
@@ -107,7 +108,8 @@ DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, d
   ConstAlias< DigitalKernel > fullKernel,
   ConstAlias< std::vector< PairIterators > > masks )
 {
-  myKernelSpelOrigin = myKSpace.sSpel( pOrigin );
+  KhalimskySpaceND<dimension, typename KSpace::Integer> preK;
+  myKernelKCoordsOrigin = preK.sKCoords(preK.sSpel( pOrigin ));
   myKernel = &fullKernel;
   myMasks = &masks;
 
@@ -589,7 +591,7 @@ DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 2
     myMasks( other.myMasks ),
     myKernel( other.myKernel ),
     myKernelMask( other.myKernelMask ),
-    myKernelSpelOrigin( other.myKernelSpelOrigin )
+    myKernelKCoordsOrigin( other.myKernelKCoordsOrigin )
 {
 }
 
@@ -601,7 +603,8 @@ DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 2
   ConstAlias< PairIterators > fullKernel,
   ConstAlias< std::vector< PairIterators > > masks )
 {
-  myKernelSpelOrigin = myKSpace.sSpel( pOrigin );
+  KhalimskySpaceND<2,typename KSpace::Integer> preK;
+  myKernelKCoordsOrigin = preK.sKCoords(preK.sSpel( pOrigin ));
   myKernelMask = &fullKernel;
   myMasks = &masks;
 
@@ -619,7 +622,8 @@ DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 2
   ConstAlias< DigitalKernel > fullKernel,
   ConstAlias< std::vector< PairIterators > > masks )
 {
-  myKernelSpelOrigin = myKSpace.sSpel( pOrigin );
+  KhalimskySpaceND<2,typename KSpace::Integer> preK;
+  myKernelKCoordsOrigin = preK.sKCoords(preK.sSpel( pOrigin ));
   myKernel = &fullKernel;
   myMasks = &masks;
 
@@ -1048,7 +1052,7 @@ DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 2
   {
     DGtal::Dimension kDim = myKSpace.sOrthDir( *it );
     currentInnerSpel = myKSpace.sDirectIncident( *it, kDim ); /// Spel on the border, but inside the shape
-    shiftInnerSpel = myKSpace.sKCoords( currentInnerSpel ) - myKSpace.sKCoords( myKernelSpelOrigin );
+    shiftInnerSpel = myKSpace.sKCoords( currentInnerSpel ) - myKernelKCoordsOrigin;
 
     /// Check if we can use previous results
     if( useLastResults )
@@ -1213,7 +1217,7 @@ DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 2
   {
     DGtal::Dimension kDim = myKSpace.sOrthDir( *it );
     currentOuterSpel = myKSpace.sIndirectIncident( *it, kDim );
-    shiftOuterSpel = myKSpace.sKCoords( currentOuterSpel ) - myKSpace.sKCoords( myKernelSpelOrigin );
+    shiftOuterSpel = myKSpace.sKCoords( currentOuterSpel ) - myKernelKCoordsOrigin;
 
     ASSERT( currentInnerSpel != currentOuterSpel );
 
@@ -1338,7 +1342,7 @@ DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 2
   {
     DGtal::Dimension kDim = myKSpace.sOrthDir( *it );
     currentInnerSpel = myKSpace.sDirectIncident( *it, kDim ); /// Spel on the border, but inside the shape
-    shiftInnerSpel = myKSpace.sKCoords( currentInnerSpel ) - myKSpace.sKCoords( myKernelSpelOrigin );
+    shiftInnerSpel = myKSpace.sKCoords( currentInnerSpel ) - myKernelKCoordsOrigin;
 
     /// Check if we can use previous results
     if( useLastResults )
@@ -1502,7 +1506,7 @@ DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 2
   {
     DGtal::Dimension kDim = myKSpace.sOrthDir( *it );
     currentOuterSpel = myKSpace.sIndirectIncident( *it, kDim );
-    shiftOuterSpel = myKSpace.sKCoords( currentOuterSpel ) - myKSpace.sKCoords( myKernelSpelOrigin );
+    shiftOuterSpel = myKSpace.sKCoords( currentOuterSpel ) - myKernelKCoordsOrigin;
 
     ASSERT( currentInnerSpel != currentOuterSpel );
 
@@ -1620,7 +1624,7 @@ DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 3
     myMasks( other.myMasks ),
     myKernel( other.myKernel ),
     myKernelMask( other.myKernelMask ),
-    myKernelSpelOrigin( other.myKernelSpelOrigin )
+    myKernelKCoordsOrigin( other.myKernelKCoordsOrigin )
 {
 }
 
@@ -1632,7 +1636,8 @@ DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 3
   ConstAlias< PairIterators > fullKernel,
   ConstAlias< std::vector< PairIterators > > masks )
 {
-  myKernelSpelOrigin = myKSpace.sSpel( pOrigin );
+  KhalimskySpaceND<3,typename KSpace::Integer> preK;
+  myKernelKCoordsOrigin = preK.sKCoords(preK.sSpel( pOrigin ));
   myKernelMask = &fullKernel;
   myMasks = &masks;
 
@@ -1650,10 +1655,10 @@ DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 3
   ConstAlias< DigitalKernel > fullKernel,
   ConstAlias< std::vector< PairIterators > > masks )
 {
-  myKernelSpelOrigin = myKSpace.sSpel( pOrigin );
+  KhalimskySpaceND<3,typename KSpace::Integer> preK;
+  myKernelKCoordsOrigin = preK.sKCoords(preK.sSpel( pOrigin ));
   myKernel = &fullKernel;
   myMasks = &masks;
-
   ASSERT ( myMasks->size () == 27 );
 
   isInitFullMasks = false;
@@ -2099,7 +2104,7 @@ DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 3
   {
     DGtal::Dimension kDim = myKSpace.sOrthDir( *it );
     currentInnerSpel = myKSpace.sDirectIncident( *it, kDim ); /// Spel on the border, but inside the shape
-    shiftInnerSpel = myKSpace.sKCoords( currentInnerSpel ) - myKSpace.sKCoords( myKernelSpelOrigin );
+    shiftInnerSpel = myKSpace.sKCoords( currentInnerSpel ) - myKernelKCoordsOrigin;
 
     /// Check if we can use previous results
     if( useLastResults )
@@ -2269,7 +2274,7 @@ DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 3
   {
     DGtal::Dimension kDim = myKSpace.sOrthDir( *it );
     currentOuterSpel = myKSpace.sIndirectIncident( *it, kDim );
-    shiftOuterSpel = myKSpace.sKCoords( currentOuterSpel ) - myKSpace.sKCoords( myKernelSpelOrigin );
+    shiftOuterSpel = myKSpace.sKCoords( currentOuterSpel ) - myKernelKCoordsOrigin;
 
     ASSERT( currentInnerSpel != currentOuterSpel );
 
@@ -2399,7 +2404,7 @@ DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 3
   {
     DGtal::Dimension kDim = myKSpace.sOrthDir( *it );
     currentInnerSpel = myKSpace.sDirectIncident( *it, kDim ); /// Spel on the border, but inside the shape
-    shiftInnerSpel = myKSpace.sKCoords( currentInnerSpel ) - myKSpace.sKCoords( myKernelSpelOrigin );
+    shiftInnerSpel = myKSpace.sKCoords( currentInnerSpel ) - myKernelKCoordsOrigin;
 
     /// Check if we can use previous results
     if( useLastResults )
@@ -2569,7 +2574,7 @@ DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 3
   {
     DGtal::Dimension kDim = myKSpace.sOrthDir( *it );
     currentOuterSpel = myKSpace.sIndirectIncident( *it, kDim );
-    shiftOuterSpel = myKSpace.sKCoords( currentOuterSpel ) - myKSpace.sKCoords( myKernelSpelOrigin );
+    shiftOuterSpel = myKSpace.sKCoords( currentOuterSpel ) - myKernelKCoordsOrigin;
 
     ASSERT( currentInnerSpel != currentOuterSpel );
 

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;
+    };
+  };
+
+}