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)
+
   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,74 @@
+#pragma once
+
+#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
+  {
+
+      Domain domain; //< The actual split domain
+      uint32_t hintVoxelCount = 0; //< An expected guess for the number of voxel
+  };
+
+  // TODO: CDomainSplitter concept ?
+  /**
+   * @brief Splits a domain evenly along all dimensions
+   *
+   * @tparam The model of domain to split
+   */
+  template<typename Domain>
+  struct EvenDomainSplitter
+  {
+    /**
+     * @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 wanted number of splits
+     */
+    std::vector<SplitInfo<Domain>> 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;
+      std::vector<SplitInfo<Domain>> 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;
+    };
+  };
+}