diff --git a/ChangeLog.md b/ChangeLog.md
index cb66f6b21d..fb50abb03e 100644
--- a/ChangeLog.md
+++ b/ChangeLog.md
@@ -245,6 +245,8 @@
(Pablo Hernandez-Cerdan, [#1505](https://github.com/DGtal-team/DGtal/pull/1505))
- Fix fillData in CubicalComplex
(Pablo Hernandez-Cerdan, [#1519](https://github.com/DGtal-team/DGtal/pull/1519))
+ - Add functions to SplitFunctions header to divide a domain into hypercubes.
+ (Pablo Hernandez-Cerdan,[#1521](https://github.com/DGtal-team/DGtal/pull/1521))
- *Shapes package*
- Add a moveTo(const RealPoint& point) method to implicit and star shapes
diff --git a/src/DGtal/topology/CubicalComplexFunctions.h b/src/DGtal/topology/CubicalComplexFunctions.h
index 4fb311893e..9ff7346cbe 100644
--- a/src/DGtal/topology/CubicalComplexFunctions.h
+++ b/src/DGtal/topology/CubicalComplexFunctions.h
@@ -267,7 +267,8 @@ namespace DGtal
const CubicalComplex< TKSpace, TCellContainer >& S2 )
{
typedef CubicalComplex< TKSpace, TCellContainer > CC;
- ASSERT( &(S1.space()) == &(S2.space()) );
+ ASSERT( S1.space().lowerBound() == S2.space().lowerBound() &&
+ S1.space().upperBound() == S2.space().upperBound());
for ( Dimension i = 0; i <= CC::dimension; ++i )
if ( ! functions::isEqual( S1.myCells[ i ], S2.myCells[ i ] ) )
return false;
@@ -294,7 +295,8 @@ namespace DGtal
const CubicalComplex< TKSpace, TCellContainer >& S2 )
{
typedef CubicalComplex< TKSpace, TCellContainer > CC;
- ASSERT( &(S1.space()) == &(S2.space()) );
+ ASSERT( S1.space().lowerBound() == S2.space().lowerBound() &&
+ S1.space().upperBound() == S2.space().upperBound());
for ( Dimension i = 0; i <= CC::dimension; ++i )
if ( ! functions::isEqual( S1.myCells[ i ], S2.myCells[ i ] ) )
return true;
@@ -318,7 +320,8 @@ namespace DGtal
const CubicalComplex< TKSpace, TCellContainer >& S2 )
{
typedef CubicalComplex< TKSpace, TCellContainer > CC;
- ASSERT( &(S1.space()) == &(S2.space()) );
+ ASSERT( S1.space().lowerBound() == S2.space().lowerBound() &&
+ S1.space().upperBound() == S2.space().upperBound());
for ( Dimension i = 0; i <= CC::dimension; ++i )
if ( ! functions::isSubset( S1.myCells[ i ], S2.myCells[ i ] ) )
return false;
@@ -341,7 +344,8 @@ namespace DGtal
const CubicalComplex< TKSpace, TCellContainer >& S2 )
{
typedef CubicalComplex< TKSpace, TCellContainer > CC;
- ASSERT( &(S1.space()) == &(S2.space()) );
+ ASSERT( S1.space().lowerBound() == S2.space().lowerBound() &&
+ S1.space().upperBound() == S2.space().upperBound());
for ( Dimension i = 0; i <= CC::dimension; ++i )
if ( ! functions::isSubset( S2.myCells[ i ], S1.myCells[ i ] ) )
return false;
diff --git a/src/DGtal/topology/SplitFunctions.h b/src/DGtal/topology/SplitFunctions.h
new file mode 100644
index 0000000000..d31fe66467
--- /dev/null
+++ b/src/DGtal/topology/SplitFunctions.h
@@ -0,0 +1,306 @@
+/**
+ * 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 .
+ *
+ **/
+
+#pragma once
+
+/**
+ * @file SplitFunctions.h
+ * @author Pablo Hernandez-Cerdan (\c pablo.hernandez.cerdan@outlook.com)
+ *
+ * @date 2018/01/08
+ *
+ * Defines functions associated to split domains for multi-threading purposes.
+ *
+ * This file is part of the DGtal library.
+ */
+
+#if defined(SplitFunctions_RECURSES)
+#error Recursive header files inclusion detected in SplitFunctions.h
+#else // defined(SplitFunctions_RECURSES)
+/** Prevents recursive inclusion of headers. */
+#define SplitFunctions_RECURSES
+
+#if !defined SplitFunctions_h
+/** Prevents repeated inclusion of headers. */
+#define SplitFunctions_h
+
+//////////////////////////////////////////////////////////////////////////////
+// Inclusions
+#include "DGtal/base/Common.h"
+//////////////////////////////////////////////////////////////////////////////
+namespace DGtal
+{
+ /**
+ * Output struct of @sa getSplit
+ *
+ * @tparam TPoint
+ */
+ template
+ struct SplitBounds {
+ size_t splitIndex;
+ TPoint lowerBound;
+ TPoint upperBound;
+ std::vector splittedRegionIndex;
+ };
+
+ namespace functions {
+ /**
+ * Compute the number of splits given a region defined by lowerBound and upperBound
+ *
+ * @tparam TPoint point type for lower and upper bound
+ * @param[in] requested_number_of_splits number of splits (the actual splits might be less)
+ * @param[in] lowerBound the lower bound point of the domain
+ * @param[in] upperBound the upper bound point of the domain
+ * @param[in,out] splits the splits for each dimension
+ * For example, initialize splits with a std::vector, and use data to get the raw array
+ * needed for this function.
+ * @code
+ * std::vector splits(dimension);
+ * auto number_of_splits = computeSplits(requested_number_of_splits,
+ * lowerBound, upperBound, splits.data());
+ * @endcode
+ *
+ * @return number of splits, lesser or equal than requested_number_of_splits
+ */
+ template
+ size_t computeSplits(
+ const size_t requested_number_of_splits,
+ const TPoint & lowerBound,
+ const TPoint & upperBound,
+ unsigned int splits[]);
+
+ /**
+ * Get the ith split, where i has to be less than the number of splits
+ * computed by @sa computeSplits.
+ *
+ * Returns two points corresponding to the output_lowerBound and output_upperRegion
+ * of the ith split.
+ *
+ * @tparam TPoint point type for lower and upper bound
+ * @param[in] splitIndex the ith split to be computed
+ * @param[in] requested_number_of_splits number of splits (the actual splits might be less)
+ * @param[in] lowerBound the lower bound point of the domain
+ * @param[in] upperBound the upper bound point of the domain
+ *
+ * @return array containing output_lowerBound and output_upperBound
+ */
+ template
+ SplitBounds getSplit(
+ const size_t splitIndex,
+ const size_t requested_number_of_splits,
+ const TPoint & lowerBound,
+ const TPoint & upperBound);
+ } // namespace functions
+
+
+ /**
+ * Return struct of @sa splitComplex
+ *
+ * Holds the splitted sub_complexes and information about the splits.
+ *
+ * @tparam TComplex Complex type
+ */
+ template
+ struct SplittedComplexes {
+ /** Number of splits (might be less than the requested_number_of_splits) */
+ size_t number_of_splits;
+ /** Vector of sub_complexes, with the data copied from the original complex */
+ std::vector sub_complexes;
+ /** number of splits per dimension */
+ std::vector splits;
+ /** Two points defining the domain of the splits */
+ using PointsPair = std::array;
+ /** Vector with the domains of the original splits. */
+ std::vector> splits_bounds;
+ };
+
+ namespace functions {
+
+ /**
+ * Split a CubicalComplex (or VoxelComplex) into sub_complexes.
+ *
+ * @tparam TComplex
+ * @param[in] vc input complex
+ * @param[in] requested_number_of_splits (the actual splits might be less)
+ *
+ * @return SplittedComplexes with the sub_complexes and the splits domain
+ *
+ * @sa computeSplits
+ * @sa getSplit
+ */
+ template < typename TComplex >
+ SplittedComplexes
+ splitComplex(
+ const TComplex & vc ,
+ const size_t requested_number_of_splits
+ );
+
+ /**
+ * Get the number of blocks based on the splits
+ *
+ * @param splits vector with the number of splits per dimension,
+ * i.e [2,2,1] means 1 division in x, and 1 divisions in y, and 0 in z, so 2 blocks.
+ *
+ * @return number of blocks
+ */
+ inline size_t
+ getNumberOfBorderBlocksFromSplits(const std::vector & splits);
+ /**
+ * Get a vector of pair of points (block_lowerBound, block_upperBound)
+ * from a domain represented by lowerBound and upperBound
+ * and the number of splits per dimension.
+ *
+ * @tparam TPoint
+ * @param lowerBound
+ * @param upperBound
+ * @param splits_bounds from SplittedComplexes @sa splitComplex
+ * @param wide_of_block_sub_complex wide of the border, defaults to just
+ * a plane in 3D, or a line in 2D.
+ * If wide is greater than zero, the blocks are hypercubes around the original border.
+ *
+ * @return vector with border blocks defined by its lower and upper bounds
+ * There is one border block per split in each dimension.
+ */
+ template
+ std::vector>
+ getBorderBlocksFromSplits(
+ const TPoint & lowerBound,
+ const TPoint & upperBound,
+ const std::vector> & splits_bounds,
+ const size_t wide_of_block_sub_complex = 0
+ );
+
+
+ /**
+ * Perform an union between out and each complex in complexes.
+ * out |= complex;
+ *
+ * @tparam TComplex complex type
+ * @param[in,out] out complex where the merge is stored
+ * @param[in] complexes vector of complexes
+ */
+ template < typename TComplex >
+ void
+ mergeSubComplexes(TComplex & out,
+ const std::vector &complexes);
+
+ template < typename TComplex >
+ TComplex
+ extractSubComplex(
+ const TComplex & vc ,
+ const typename TComplex::Point & sub_lowerBound,
+ const typename TComplex::Point & sub_upperBound
+ );
+
+ /**
+ * Given a cell and some boundaries, returns if the cell is in the border
+ * or close to the border (if wide_point is used).
+ *
+ * Used internally in:
+ * @sa getBorderVoxels, @sa getBorderVoxelsWithDistanceMap
+ *
+ * @tparam TComplex complex type
+ * @param lowerBound the lowerBound of the hypercube defining a border (usually vc.lowerBound().
+ * @param upperBound the upperBound of the hypercube defining a border (usually vc.upperBound().
+ * @param wide_point the wide of the border (3D), defaults to nullptr, which is equivalent to {1,1,1}
+ * it has to be greater than zero.
+ * @param lowerBound_to_ignore lowerBound defining an hypercube of borders to ignore
+ * @param upperBound_to_ignore upperBound defining an hypercube of borders to ignore
+ *
+ * @return true if cell is in the border, given the input parameters
+ */
+ template < typename TComplex >
+ bool
+ isInBorder(
+ const typename TComplex::KSpace::Cell & cell,
+ const typename TComplex::Point & lowerBound,
+ const typename TComplex::Point & upperBound,
+ const typename TComplex::Point * wide_point = nullptr,
+ const typename TComplex::Point * lowerBound_to_ignore = nullptr,
+ const typename TComplex::Point * upperBound_to_ignore= nullptr
+ );
+
+ /**
+ * Get the voxels from the border of the input vc.
+ * If lowerBound_to_ignore and upperBound_to_ignore are set, the borders
+ * defined by those bounds are ignored.
+ *
+ * lowerBound and upperBound are usually the lower and upper bound
+ * of the input complex vc, but they can also be changed to get the
+ * borders of a smaller hypercube contained in vc.
+ *
+ * @tparam TComplex
+ * @param vc input complex
+ * @param lowerBound the lowerBound of the hypercube defining a border (usually vc.lowerBound().
+ * @param upperBound the upperBound of the hypercube defining a border (usually vc.upperBound().
+ * @param wide_point the wide of the border (3D), defaults to nullptr, which is equivalent to {1,1,1}
+ * it has to be greater than zero.
+ * @param lowerBound_to_ignore lowerBound defining an hypercube of borders to ignore
+ * @param upperBound_to_ignore upperBound defining an hypercube of borders to ignore
+ *
+ * @return vector with iterators of the border
+ */
+ template < typename TComplex >
+ std::vector
+ getBorderVoxels(
+ TComplex & vc,
+ const typename TComplex::Point & lowerBound,
+ const typename TComplex::Point & upperBound,
+ const typename TComplex::Point * wide_point = nullptr,
+ const typename TComplex::Point * lowerBound_to_ignore = nullptr,
+ const typename TComplex::Point * upperBound_to_ignore = nullptr
+ );
+ template < typename TComplex, typename TDistanceTransform >
+ std::vector
+ getBorderVoxelsWithDistanceMap(
+ TComplex & vc,
+ const typename TComplex::Point & lowerBound,
+ const typename TComplex::Point & upperBound,
+ const typename TComplex::Point * wide_point = nullptr,
+ const typename TComplex::Point * lowerBound_to_ignore = nullptr,
+ const typename TComplex::Point * upperBound_to_ignore = nullptr,
+ const TDistanceTransform * dist_map = nullptr
+ );
+
+ /**
+ * Helper function to get the border voxels obtained from
+ * @sa getBorderVoxels to the desired data.
+ *
+ * @tparam TComplex Complex type
+ * @param vc input complex
+ * @param border_iterators iterators from getBorderVoxels
+ * @param data data of the cell to set in the border
+ */
+ template < typename TComplex >
+ void
+ setBorderData(
+ TComplex & vc,
+ const std::vector &border_iterators,
+ const DGtal::uint32_t &data);
+ } // namespace functions
+} // namespace DGtal
+
+///////////////////////////////////////////////////////////////////////////////
+// Includes inline functions.
+#include "DGtal/topology/SplitFunctions.ih"
+
+// //
+///////////////////////////////////////////////////////////////////////////////
+
+#endif // !defined SplitFunctions_h
+
+#undef SplitFunctions_RECURSES
+#endif // else defined(SplitFunctions_RECURSES)
diff --git a/src/DGtal/topology/SplitFunctions.ih b/src/DGtal/topology/SplitFunctions.ih
new file mode 100644
index 0000000000..899d6a4b0e
--- /dev/null
+++ b/src/DGtal/topology/SplitFunctions.ih
@@ -0,0 +1,581 @@
+/**
+ * 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 .
+ *
+ **/
+
+/**
+ * @file SplitFunctions.ih
+ * @author Pablo Hernandez-Cerdan (\c pablo.hernandez.cerdan@outlook.com)
+ *
+ * @date 2019/01/08
+ *
+ * Implementation of inline methods defined in SplitFunctions.h
+ *
+ * This file is part of the DGtal library.
+ */
+
+template
+size_t DGtal::functions::computeSplits(
+ const size_t requested_number_of_splits,
+ const TPoint & lowerBound,
+ const TPoint & upperBound,
+ unsigned int splits[]
+ )
+{
+ const auto dim = lowerBound.dimension;
+ auto region_size = upperBound - lowerBound;
+ for ( size_t i = 0; i < dim; ++i )
+ {
+ if ( region_size[ i ] < 0 )
+ throw std::runtime_error(
+ "upperBound is smaller than lowerBound, region_size is negative" );
+ }
+
+ // Implementation taken from ITK:
+ // itkImageRegionSplitterMultidimensional.cxx ComputeSplitInternal
+ // size of each splited region
+ std::vector splitRegionSize(dim); // Note: stack allocation preferred
+ unsigned int numberOfPieces = 1;
+
+ // initialize arrays
+ for (unsigned int i = 0; i < dim; ++i)
+ {
+ splits[i] = 1;
+ splitRegionSize[i] = region_size[i];
+ }
+
+ while (true)
+ {
+ // find the dimension with the largest size
+ unsigned int maxSplitDim = 0;
+ for (unsigned int i = 1; i < dim; ++i)
+ {
+ if (splitRegionSize[maxSplitDim] < splitRegionSize[i])
+ {
+ maxSplitDim = i;
+ }
+ }
+
+ // calculate the number of additional pieces this split would add
+ unsigned int additionalNumPieces = 1;
+ for (unsigned int i = 0; i < dim; ++i)
+ {
+ if (i != maxSplitDim)
+ {
+ additionalNumPieces *= splits[i];
+ }
+ }
+
+ // check if this would give us too many pieces or
+ // if this would require the subpixel splits
+ if ( numberOfPieces + additionalNumPieces > requested_number_of_splits ||
+ splits[ maxSplitDim ] == static_cast(region_size[ maxSplitDim ]) )
+ {
+ return numberOfPieces;
+ }
+
+ // update the variable with the new split
+ numberOfPieces += additionalNumPieces;
+ ++splits[maxSplitDim];
+ splitRegionSize[maxSplitDim] = region_size[maxSplitDim] / double(splits[maxSplitDim]);
+ }
+
+}
+template
+DGtal::SplitBounds DGtal::functions::getSplit(
+ const size_t splitIndex,
+ const size_t requested_number_of_splits,
+ const TPoint & lowerBound,
+ const TPoint & upperBound
+ )
+{
+ SplitBounds out;
+ out.splitIndex = splitIndex;
+ out.lowerBound = lowerBound;
+ out.upperBound = upperBound;
+ auto & outputLowerBound = out.lowerBound;
+ auto & outputUpperBound = out.upperBound;
+
+ const auto dim = lowerBound.dimension;
+ const auto region_size = upperBound - lowerBound;
+ auto output_region_size = region_size;
+ // Implementation taken from ITK:
+ // itkImageRegionSplitterMultidimensional.cxx GetSplitInternal
+ // number of splits in each dimension
+ std::vector splits(dim); // Note: stack allocation preferred
+
+ // index into splitted regions
+ out.splittedRegionIndex = std::vector(dim); // Note: stack allocation preferred
+ auto & splittedRegionIndex = out.splittedRegionIndex;
+
+ DGtal::functions::computeSplits(requested_number_of_splits,
+ lowerBound, upperBound, &splits[0]);
+
+ // determine which splitted region we are in
+ unsigned int offset = splitIndex;
+ for (unsigned i = dim - 1; i > 0; --i)
+ {
+ unsigned int dimensionOffset = 1;
+ for (unsigned int j = 0; j < i; ++j)
+ {
+ dimensionOffset *= splits[j];
+ }
+
+ splittedRegionIndex[i] = offset / dimensionOffset;
+ offset -= (splittedRegionIndex[i] * dimensionOffset);
+ }
+ splittedRegionIndex[0] = offset;
+
+
+ // Assign the output split region to the input region in-place
+ for (unsigned int i = 0; i < dim; i++)
+ {
+ const auto inputRegionSize = region_size[i];
+ const auto indexOffset =
+ std::floor((splittedRegionIndex[i]) * (inputRegionSize / double(splits[i])));
+
+ outputLowerBound[i] += indexOffset;
+ if (splittedRegionIndex[i] < splits[i] - 1)
+ {
+ // Modified from ITK, minus -1 (index,size in ITK versus lower/upperBounds here)
+ output_region_size[i] =
+ std::floor((splittedRegionIndex[i] + 1) * (inputRegionSize / double(splits[i]))) - indexOffset - 1;
+ }
+ else
+ {
+ // this dimension is falling off the edge of the image
+ output_region_size[i] = output_region_size[i] - indexOffset;
+ }
+ }
+ outputUpperBound = outputLowerBound + output_region_size;
+
+ return out;
+}
+
+template < typename TComplex >
+TComplex
+DGtal::functions::extractSubComplex(
+ const TComplex & vc ,
+ const typename TComplex::Point & sub_lowerBound,
+ const typename TComplex::Point & sub_upperBound
+ ) {
+ const auto & kspace = vc.space();
+ const bool kspace_isClosed = true;
+ typename TComplex::KSpace sub_kspace;
+ sub_kspace.init(sub_lowerBound, sub_upperBound, kspace_isClosed);
+ TComplex sub_complex(sub_kspace);
+
+ const auto sub_lowerCell = sub_kspace.lowerCell();
+ const auto sub_upperCell = sub_kspace.upperCell();
+ const auto & sub_lowerCell_coords = kspace.uKCoords(sub_lowerCell);
+ const auto & sub_upperCell_coords = kspace.uKCoords(sub_upperCell);
+ auto kcell = sub_lowerCell;
+ static_assert(TComplex::KSpace::dimension == 3,
+ "The KSpace::dimension of TComplex isn't 3");
+ for(auto k_x = sub_lowerCell_coords[0]; k_x <= sub_upperCell_coords[0]; k_x++) {
+ for(auto k_y = sub_lowerCell_coords[1]; k_y <= sub_upperCell_coords[1]; k_y++) {
+ for(auto k_z = sub_lowerCell_coords[2]; k_z <= sub_upperCell_coords[2]; k_z++) {
+ kspace.uSetKCoord(kcell, 0, k_x);
+ kspace.uSetKCoord(kcell, 1, k_y);
+ kspace.uSetKCoord(kcell, 2, k_z);
+ const auto dim_cell = kspace.uDim(kcell);
+ const auto iter_cell_map = vc.findCell(dim_cell, kcell);
+ if(iter_cell_map != vc.end(dim_cell)) {
+ sub_complex.insertCell(dim_cell,
+ kcell, iter_cell_map->second);
+ }
+ }
+ }
+ }
+
+ return sub_complex;
+}
+
+template < typename TComplex >
+bool
+DGtal::functions::isInBorder(
+ const typename TComplex::KSpace::Cell &cell,
+ const typename TComplex::Point & lowerBound,
+ const typename TComplex::Point & upperBound,
+ const typename TComplex::Point * wide_point,
+ const typename TComplex::Point * lowerBound_to_ignore,
+ const typename TComplex::Point * upperBound_to_ignore) {
+
+ const typename TComplex::Point wide =
+ wide_point ? *wide_point : typename TComplex::Point(1,1,1);
+
+ const auto cell_coords = TComplex::KSpace::Cell::PreCellularGridSpace::uCoords(cell);
+ // dist to lower and upper
+ const auto distToLower = cell_coords - lowerBound;
+ const auto distToUpper = upperBound - cell_coords;;
+ bool is_close_to_lower_border =
+ distToLower[0] < wide[0]
+ || distToLower[1] < wide[1]
+ || distToLower[2] < wide[2];
+ bool is_close_to_upper_border =
+ distToUpper[0] < wide[0]
+ || distToUpper[1] < wide[1]
+ || distToUpper[2] < wide[2];
+
+ if(lowerBound_to_ignore && upperBound_to_ignore) {
+ if(is_close_to_lower_border) {
+ const auto distToLower_to_ignore = cell_coords - *lowerBound_to_ignore;
+ bool is_close_to_lower_border_to_ignore =
+ distToLower_to_ignore[0] < wide[0]
+ || distToLower_to_ignore[1] < wide[1]
+ || distToLower_to_ignore[2] < wide[2];
+
+ if(is_close_to_lower_border_to_ignore) {
+ // Don't ignore if it is closer to lowerBound than
+ // lowerBound_to_ignore (corners) for any dimension
+ // Compute the border plane normal that results in non-zero
+ // for the dimensions that need checking.
+ const auto lowerBound_border_plane_normal =
+ lowerBound - *lowerBound_to_ignore;
+ for(size_t dim = 0; dim < 3; dim++) {
+ if(lowerBound_border_plane_normal[dim]) {
+ if(distToLower[dim] < wide[dim]
+ && distToLower[dim] <= distToLower_to_ignore[dim]) {
+ is_close_to_lower_border_to_ignore = false;
+ }
+ }
+ }
+ }
+ is_close_to_lower_border &= !is_close_to_lower_border_to_ignore;
+ }
+
+ if(is_close_to_upper_border) {
+ const auto distToUpper_to_ignore = *upperBound_to_ignore - cell_coords;;
+ bool is_close_to_upper_border_to_ignore =
+ distToUpper_to_ignore[0] < wide[0]
+ || distToUpper_to_ignore[1] < wide[1]
+ || distToUpper_to_ignore[2] < wide[2];
+
+ if(is_close_to_upper_border_to_ignore) {
+ // Don't ignore if it is closer to upperBound than
+ // upperBound_to_ignore (corners) for any dimension
+ // Compute the border plane normal that results in non-zero
+ // for the dimensions that need checking.
+ const auto upperBound_border_plane_normal =
+ *upperBound_to_ignore - upperBound;
+ for(size_t dim = 0; dim < 3; dim++) {
+ if(upperBound_border_plane_normal[dim]) {
+ if(distToUpper[dim] < wide[dim]
+ && distToUpper[dim] <= distToUpper_to_ignore[dim]) {
+ is_close_to_upper_border_to_ignore = false;
+ }
+ }
+ }
+ }
+ is_close_to_upper_border &= !is_close_to_upper_border_to_ignore;
+ }
+ }
+
+ const bool isInBorder_not_to_ignore = is_close_to_lower_border
+ || is_close_to_upper_border;
+ return isInBorder_not_to_ignore;
+}
+
+template < typename TComplex >
+std::vector
+DGtal::functions::getBorderVoxels(
+ TComplex & vc,
+ const typename TComplex::Point & lowerBound,
+ const typename TComplex::Point & upperBound,
+ const typename TComplex::Point * wide_point,
+ const typename TComplex::Point * lowerBound_to_ignore,
+ const typename TComplex::Point * upperBound_to_ignore)
+{
+ if(wide_point) {
+ const typename TComplex::Point wide = *wide_point;
+ // wide_point has to be positive
+ ASSERT(wide.min() > 0);
+ // wide point cannot be greater than the domain in any dimension
+ std::array is_close_to_lower_upper_border_to_ignore =
+ {false, false, false};
+ if(lowerBound_to_ignore && upperBound_to_ignore) {
+ const auto dist_lowerToUpper_to_ignore =
+ *upperBound_to_ignore - *lowerBound_to_ignore;
+ for(size_t dim = 0; dim < 3; ++dim) {
+ if(wide[dim] <= dist_lowerToUpper_to_ignore[dim]) {
+ is_close_to_lower_upper_border_to_ignore[dim] = true;
+ }
+ }
+ }
+ const auto dist_lowerToUpper = upperBound - lowerBound;
+ for(size_t dim = 0; dim < 3; ++dim) {
+ if(wide[dim] >= dist_lowerToUpper[dim]
+ && !is_close_to_lower_upper_border_to_ignore[dim]) {
+ std::cerr <<
+ "wide_point in getBorderVoxels is equal or greater than the domain itself\n" <<
+ "wide_point: " << wide << "\n" <<
+ "dist_lowerToUpper: " << dist_lowerToUpper <<
+ std::endl;
+ if(lowerBound_to_ignore && upperBound_to_ignore) {
+ std::cerr << "dist_lowerToUpper_to_ignore: " <<
+ (*upperBound_to_ignore - *lowerBound_to_ignore) << std::endl;
+ }
+ }
+ }
+ }
+
+ using CellIterators = std::vector;
+ CellIterators iterators;
+
+ for(auto it = vc.begin(3); it != vc.end(3); ++it) {
+ const auto & cell = it->first;
+
+ const bool isInBorder_not_to_ignore = isInBorder(cell,
+ lowerBound, upperBound,
+ wide_point, lowerBound_to_ignore, upperBound_to_ignore);
+
+ if(isInBorder_not_to_ignore) {
+ iterators.push_back(it);
+ }
+ }
+ return iterators;
+}
+
+template < typename TComplex, typename TDistanceTransform >
+std::vector
+DGtal::functions::getBorderVoxelsWithDistanceMap(
+ TComplex & vc,
+ const typename TComplex::Point & lowerBound,
+ const typename TComplex::Point & upperBound,
+ const typename TComplex::Point * wide_point,
+ const typename TComplex::Point * lowerBound_to_ignore,
+ const typename TComplex::Point * upperBound_to_ignore,
+ const TDistanceTransform * dist_map)
+{
+ const auto & kspace = vc.space();
+ if(wide_point) {
+ const typename TComplex::Point wide = *wide_point;
+ // wide_point has to be positive
+ ASSERT(wide.min() > 0);
+ // wide point cannot be greater than the domain in any dimension
+ std::array is_close_to_lower_upper_border_to_ignore =
+ {false, false, false};
+ if(lowerBound_to_ignore && upperBound_to_ignore) {
+ const auto dist_lowerToUpper_to_ignore =
+ *upperBound_to_ignore - *lowerBound_to_ignore;
+ for(size_t dim = 0; dim < 3; ++dim) {
+ if(wide[dim] <= dist_lowerToUpper_to_ignore[dim]) {
+ is_close_to_lower_upper_border_to_ignore[dim] = true;
+ }
+ }
+ }
+ const auto dist_lowerToUpper = upperBound - lowerBound;
+ for(size_t dim = 0; dim < 3; ++dim) {
+ if(wide[dim] >= dist_lowerToUpper[dim]
+ && !is_close_to_lower_upper_border_to_ignore[dim]) {
+ std::cerr <<
+ "wide_point in getBorderVoxels is equal or greater than the domain itself\n" <<
+ "wide_point: " << wide << "\n" <<
+ "dist_lowerToUpper: " << dist_lowerToUpper <<
+ std::endl;
+ if(lowerBound_to_ignore && upperBound_to_ignore) {
+ std::cerr << "dist_lowerToUpper_to_ignore: " <<
+ (*upperBound_to_ignore - *lowerBound_to_ignore) << std::endl;
+ }
+ }
+ }
+ }
+
+ using CellIterators = std::vector;
+ CellIterators iterators;
+
+ for(auto it = vc.begin(3); it != vc.end(3); ++it) {
+ const auto & cell = it->first;
+
+ const bool isInBorder_not_to_ignore = isInBorder( cell,
+ lowerBound, upperBound,
+ wide_point, lowerBound_to_ignore, upperBound_to_ignore);
+
+ if(isInBorder_not_to_ignore && dist_map) {
+ const auto cell_coords = kspace.uCoords(cell);
+ const auto distToLower = cell_coords - lowerBound;
+ const auto distToUpper = upperBound - cell_coords;;
+ const auto dist_map_value = (*dist_map)(cell_coords);
+ bool closer_to_block_border_than_object_border = false;
+ for(size_t dim = 0; dim < 3; ++dim) {
+ if( dist_map_value <= distToLower[dim]
+ || dist_map_value <= distToUpper[dim]) {
+ closer_to_block_border_than_object_border = true;
+ }
+ }
+ if(closer_to_block_border_than_object_border) {
+ iterators.push_back(it);
+ }
+ }
+ // With no dmap is exactly the same as getBorderVoxels
+ else if(isInBorder_not_to_ignore) {
+ iterators.push_back(it);
+ }
+ }
+ return iterators;
+}
+
+
+template < typename TComplex >
+void
+DGtal::functions::setBorderData(
+ TComplex & vc,
+ const std::vector &border_iterators,
+ const DGtal::uint32_t &data) {
+ for(auto vit = border_iterators.begin(); vit != border_iterators.end(); ++vit){
+ (*vit)->second = data;
+ // Set the data of all faces as well
+ const auto faces = vc.cellBoundary((*vit)->first);
+ for(const auto & face : faces) {
+ auto it = vc.findCell(face);
+ it->second = data;
+ }
+
+ }
+}
+
+template < typename TComplex >
+DGtal::SplittedComplexes
+DGtal::functions::splitComplex(
+ const TComplex & vc ,
+ const size_t requested_number_of_splits
+ ) {
+ const auto & kspace = vc.space();
+ const auto & lowerBound = kspace.lowerBound();
+ const auto & upperBound = kspace.upperBound();
+ std::vector splits(kspace.dimension);
+ const auto number_of_splits =
+ DGtal::functions::computeSplits(requested_number_of_splits,
+ lowerBound, upperBound, splits.data());
+ std::vector sub_kspaces(number_of_splits);
+ // Initiate output struct
+ SplittedComplexes out;
+ out.number_of_splits = number_of_splits;
+ out.splits = splits;
+ auto & sub_complexes = out.sub_complexes;
+ // Note that when dividing a k-space special care hast to be taken in the
+ // opennes or clossedness of the low dim cells (1-cell and 2-cell in 3D).
+ // To avoid sharing those pointels and linels in the frontier between subspaces,
+ // only one of the colindant sub_kspaces should include them, and the rest omit them.
+ // In case of Cubical and Voxel Complexes, those cells are not that signficant,
+ // at least for thinning, so we make all the sub_kspaces closed.
+ const bool kspace_isClosed = true;
+ for(size_t split = 0; split < number_of_splits; ++split) {
+ const auto split_bounds = DGtal::functions::getSplit(
+ split, requested_number_of_splits, lowerBound, upperBound);
+ // Save the domain of the original splits in the output
+ out.splits_bounds.push_back(split_bounds);
+ const auto & sub_lowerBound = split_bounds.lowerBound;
+ const auto & sub_upperBound = split_bounds.upperBound;
+ sub_kspaces[split].init(sub_lowerBound, sub_upperBound, kspace_isClosed);
+ }
+
+ for(size_t split = 0; split < number_of_splits; ++split) {
+ sub_complexes.emplace_back(TComplex(sub_kspaces[split]));
+ // Note: You don't want to use simplicity table, some spels are going to be FIXED
+ // Also, isSimple from object does not give the same results than isSimpleByThinning
+ // being the latter correct and more general in the case of FIXED data.
+ // sub_complexes[split].copySimplicityTable(vc);
+
+ // Copy cells and cell data from input vc to sub_complexes
+ // First, iterate over the kspace of the sub-domain
+ // Then find cell in original complex, if exists, insert it (copy) into the sub_complex
+ const auto sub_lowerCell = sub_kspaces[split].lowerCell();
+ const auto sub_upperCell = sub_kspaces[split].upperCell();
+ const auto & sub_lowerCell_coords = kspace.uKCoords(sub_lowerCell);
+ const auto & sub_upperCell_coords = kspace.uKCoords(sub_upperCell);
+ auto kcell = sub_lowerCell;
+ static_assert(TComplex::KSpace::dimension == 3,
+ "The KSpace::dimension of TComplex isn't 3");
+ for(auto k_x = sub_lowerCell_coords[0]; k_x <= sub_upperCell_coords[0]; k_x++) {
+ for(auto k_y = sub_lowerCell_coords[1]; k_y <= sub_upperCell_coords[1]; k_y++) {
+ for(auto k_z = sub_lowerCell_coords[2]; k_z <= sub_upperCell_coords[2]; k_z++) {
+ kspace.uSetKCoord(kcell, 0, k_x);
+ kspace.uSetKCoord(kcell, 1, k_y);
+ kspace.uSetKCoord(kcell, 2, k_z);
+ const auto dim_cell = kspace.uDim(kcell);
+ const auto iter_cell_map = vc.findCell(dim_cell, kcell);
+ if(iter_cell_map != vc.end(dim_cell)) {
+ // iter_cell_map->first is equal to kcell...
+ sub_complexes[split].insertCell(dim_cell,
+ iter_cell_map->first, iter_cell_map->second);
+ }
+ }
+ }
+ }
+ }
+ return out;
+}
+
+size_t
+DGtal::functions::getNumberOfBorderBlocksFromSplits(
+ const std::vector & splits)
+{
+ const auto dimension = splits.size();
+ size_t number_of_block_complexes = 0;
+ for(size_t dim = 0; dim < dimension; ++dim) {
+ for(size_t split = 1; split < splits[dim]; ++split) {
+ number_of_block_complexes++;
+ }
+ }
+ return number_of_block_complexes;
+}
+
+template
+std::vector>
+DGtal::functions::getBorderBlocksFromSplits(
+ const TPoint & lowerBound,
+ const TPoint & upperBound,
+ const std::vector> & splits_bounds,
+ const size_t wide_of_block_sub_complex
+ )
+{
+ // The blocks spans the whole domain (large thin layers), and they only differ from lowerBound, upperBound in one dimension where the split is.
+ using PointPair = std::array;
+ std::vector block_lower_upper_bounds;
+ std::vector> visited_index_per_dim(lowerBound.dimension);
+ for(const auto & split_bound : splits_bounds) {
+ const auto & split_index = split_bound.splittedRegionIndex;
+ for(size_t dim = 0; dim < split_bound.lowerBound.dimension; dim++) {
+ if(split_index[dim] > 0 &&
+ visited_index_per_dim[dim].find(split_index[dim]) ==
+ visited_index_per_dim[dim].end()) {
+ visited_index_per_dim[dim].emplace(split_index[dim]);
+ auto block_lowerBound = lowerBound;
+ auto block_upperBound = upperBound;
+ block_lowerBound[dim] = split_bound.lowerBound[dim];
+ block_upperBound[dim] = split_bound.lowerBound[dim];
+ if(wide_of_block_sub_complex) {
+ block_lowerBound[dim] -= wide_of_block_sub_complex;
+ block_upperBound[dim] += wide_of_block_sub_complex - 1;
+ if(block_lowerBound[dim] < lowerBound[dim]) block_lowerBound[dim] = lowerBound[dim];
+ if(block_upperBound[dim] > upperBound[dim]) block_upperBound[dim] = upperBound[dim];
+ }
+ block_lower_upper_bounds.emplace_back(
+ PointPair({block_lowerBound, block_upperBound}));
+ }
+ }
+ }
+ return block_lower_upper_bounds;
+}
+
+template < typename TComplex >
+void
+DGtal::functions::mergeSubComplexes(TComplex & out,
+ const std::vector &complexes)
+{
+ for(auto & sc : complexes) {
+ out |= sc;
+ }
+}
diff --git a/src/DGtal/topology/VoxelComplex.ih b/src/DGtal/topology/VoxelComplex.ih
index 0b8e429853..76eb591197 100644
--- a/src/DGtal/topology/VoxelComplex.ih
+++ b/src/DGtal/topology/VoxelComplex.ih
@@ -828,7 +828,8 @@ DGtal::VoxelComplex::criticalCliquesForD(
#pragma omp single nowait
{for (auto it = cubical.begin(d), itE = cubical.end(d); it != itE; ++it)
#pragma omp task firstprivate(it)
- {auto clique_p = criticalCliquePair(d, it);
+ {
+ auto clique_p = criticalCliquePair(d, it);
const auto &is_critical = clique_p.first;
const auto &clique = clique_p.second;
// Push
diff --git a/src/DGtal/topology/VoxelComplexFunctions.ih b/src/DGtal/topology/VoxelComplexFunctions.ih
index 5c50e418b0..6b42f488f6 100644
--- a/src/DGtal/topology/VoxelComplexFunctions.ih
+++ b/src/DGtal/topology/VoxelComplexFunctions.ih
@@ -81,6 +81,14 @@ asymetricThinningScheme(
typename TComplex::Parent x_y(vc.space());
typename TComplex::Parent x_k(vc.space());
+ // Initialize the constraint set
+ for (auto it = vc.begin(3), itE = vc.end(3) ; it != itE ; ++it ){
+ const auto &constraint_voxel = it->first ;
+ if( Skel(vc, constraint_voxel) == true || it->second.data == TComplex::FIXED){
+ K.insertVoxelCell(constraint_voxel);
+ }
+ }
+
bool stability{false};
uint64_t generation{0};
auto xsize = X.nbCells(3);
@@ -128,7 +136,7 @@ asymetricThinningScheme(
x_k = X - K;
for (auto it = x_k.begin(3), itE = x_k.end(3) ; it != itE ; ++it ){
auto new_voxel = it->first ;
- if( Skel(X, new_voxel) == true){
+ if( Skel(X, new_voxel) == true || it->second.data == TComplex::FIXED){
K.insertVoxelCell(new_voxel);
// K.insertCell(3, new_voxel);
// K.objectSet().insert(K.space().uCoords(new_voxel));
diff --git a/tests/topology/CMakeLists.txt b/tests/topology/CMakeLists.txt
index 2a36f4c15f..66b848c1a5 100644
--- a/tests/topology/CMakeLists.txt
+++ b/tests/topology/CMakeLists.txt
@@ -16,6 +16,7 @@ set(DGTAL_TESTS_SRC
testParDirCollapse
testHalfEdgeDataStructure
testIndexedDigitalSurface
+ testSplitFunctions
)
foreach(FILE ${DGTAL_TESTS_SRC})
diff --git a/tests/topology/testSplitFunctions.cpp b/tests/topology/testSplitFunctions.cpp
new file mode 100644
index 0000000000..421c344d97
--- /dev/null
+++ b/tests/topology/testSplitFunctions.cpp
@@ -0,0 +1,341 @@
+/**
+ * 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 .
+ *
+ **/
+
+/**
+ * @file testSplitFunctions.cpp
+ * @ingroup Tests
+ * @author Pablo Hernandez-Cerdan (\c pablo.hernandez.cerdan@outlook.com)
+ *
+ * @date 2019/01/08
+ *
+ * Testing class for SplitFunctions
+ *
+ * This file is part of the DGtal library.
+ */
+
+///////////////////////////////////////////////////////////////////////////////
+#include "DGtal/base/Common.h"
+#include "DGtal/helpers/StdDefs.h"
+#include "DGtal/topology/SplitFunctions.h"
+#include "DGtalCatch.h"
+#include "DGtal/topology/KhalimskyCellHashFunctions.h"
+#include "DGtal/topology/VoxelComplex.h"
+#include "DGtal/topology/VoxelComplexFunctions.h"
+///////////////////////////////////////////////////////////////////////////////
+using namespace std;
+using namespace DGtal;
+using Point = DGtal::Z3i::Point;
+using DigitalSet = DGtal::DigitalSetByAssociativeContainer<
+DGtal::Z3i::Domain, std::unordered_set>;
+DigitalSet fillSet(const Point & lowerBound, const Point & upperBound) {
+ using Value = typename Point::Component;
+ using Domain = DGtal::Z3i::Domain;
+ Domain domain(lowerBound, upperBound);
+ DigitalSet a_set(domain);
+ for(Value x = lowerBound[0]; x <= upperBound[0]; x++)
+ for(Value y = lowerBound[1]; y <= upperBound[1]; y++)
+ for(Value z = lowerBound[2]; z <= upperBound[2]; z++) {
+ a_set.insert(Point(x,y,z));
+ }
+ return a_set;
+}
+
+TEST_CASE("getBlocksFromSplits", "[blocks]"){
+ using namespace DGtal::functions;
+ using KSpace = DGtal::Z3i::KSpace;
+ Point lowerBound = {0,0,0};
+ Point upperBound = {13,13,13};
+ using Complex = DGtal::VoxelComplex;
+ KSpace kspace;
+ kspace.init(lowerBound, upperBound, true);
+ Complex vc(kspace);
+ vc.construct(fillSet(lowerBound, upperBound));
+ const size_t requested_number_of_splits = 4;
+ auto splitted_complexes = splitComplex(vc, requested_number_of_splits);
+ auto & splits_bounds = splitted_complexes.splits_bounds;
+ auto & splits = splitted_complexes.splits;
+ SECTION("No wide_of_block_sub_complex") {
+ auto block_lower_upper_bounds =
+ getBorderBlocksFromSplits(lowerBound, upperBound, splits_bounds);
+ auto number_of_blocks = getNumberOfBorderBlocksFromSplits(splits);
+ const auto expected_number_of_blocks = 2;
+ CHECK(number_of_blocks == expected_number_of_blocks);
+ REQUIRE(block_lower_upper_bounds.size() == expected_number_of_blocks);
+
+ Point expected_lowerBound_x0 = Point(6, 0, 0);
+ Point expected_upperBound_x0 = Point(6, 13, 13);
+ Point expected_lowerBound_y0 = Point(0, 6, 0);
+ Point expected_upperBound_y0 = Point(13, 6, 13);
+ CHECK(block_lower_upper_bounds[0][0] == expected_lowerBound_x0);
+ CHECK(block_lower_upper_bounds[0][1] == expected_upperBound_x0);
+ CHECK(block_lower_upper_bounds[1][0] == expected_lowerBound_y0);
+ CHECK(block_lower_upper_bounds[1][1] == expected_upperBound_y0);
+ }
+ SECTION("wide_of_block_sub_complex = 2") {
+ size_t wide_of_block_sub_complex = 2;
+ auto block_lower_upper_bounds =
+ getBorderBlocksFromSplits(lowerBound, upperBound, splits_bounds, wide_of_block_sub_complex);
+ Point expected_lowerBound_x0 = Point(4, 0, 0);
+ Point expected_upperBound_x0 = Point(7, 13, 13);
+ Point expected_lowerBound_y0 = Point(0, 4, 0);
+ Point expected_upperBound_y0 = Point(13, 7, 13);
+ CHECK(block_lower_upper_bounds[0][0] == expected_lowerBound_x0);
+ CHECK(block_lower_upper_bounds[0][1] == expected_upperBound_x0);
+ CHECK(block_lower_upper_bounds[1][0] == expected_lowerBound_y0);
+ CHECK(block_lower_upper_bounds[1][1] == expected_upperBound_y0);
+ }
+}
+
+TEST_CASE("computeSplitsEasy", "[computeSplits, getSplit]") {
+ Point lowerBound = {0,0,0};
+ Point upperBound = {4,6,8};
+ size_t requested_number_of_splits = 2;
+ std::vector splits(lowerBound.dimension);
+ trace.beginBlock("computeSplits Easy");
+ trace.info() << "lowerBound: " << lowerBound << std::endl;
+ trace.info() << "upperBound: " << upperBound << std::endl;
+ auto number_of_splits = functions::computeSplits(
+ requested_number_of_splits,
+ lowerBound, upperBound,
+ splits.data());
+ trace.info() << "number_of_splits: " << number_of_splits << std::endl;
+ trace.info() << "splits: " << splits[ 0 ] << ", " << splits[ 1 ] << ", "
+ << splits[ 2 ] << std::endl;
+ CHECK( number_of_splits == requested_number_of_splits );
+ CHECK( splits[0] == 1 );
+ CHECK( splits[1] == 1 );
+ CHECK( splits[2] == 2 );
+ trace.endBlock();
+ trace.beginBlock("getSplit Easy");
+ for (size_t split_index = 0; split_index < number_of_splits; split_index++) {
+ auto split_bounds =
+ DGtal::functions::getSplit(split_index, requested_number_of_splits,
+ lowerBound, upperBound);
+ const auto & outputLowerBound = split_bounds.lowerBound;
+ const auto & outputUpperBound = split_bounds.upperBound;
+ trace.info() << "split_index: " << split_index << std::endl;
+ trace.info() << "outputLowerBound: " << outputLowerBound << std::endl;
+ trace.info() << "outputUpperBound: " << outputUpperBound << std::endl;
+ }
+ size_t split_index = 1;
+ auto split_bounds =
+ DGtal::functions::getSplit(split_index, requested_number_of_splits,
+ lowerBound, upperBound);
+ const auto & outputLowerBound = split_bounds.lowerBound;
+ const auto & outputUpperBound = split_bounds.upperBound;
+ CHECK( outputLowerBound[2] == 4 );
+ CHECK( outputUpperBound[2] == 8 );
+ trace.endBlock();
+}
+////////////////////// splitComplex ///////////////////////////
+// Fixture for a X
+struct Fixture_X_with_tight_bounds {
+ ///////////////////////////////////////////////////////////
+ // type aliases
+ ///////////////////////////////////////////////////////////
+ using Point = DGtal::Z3i::Point;
+ using Domain = DGtal::Z3i::Domain;
+ using KSpace = DGtal::Z3i::KSpace;
+
+ using FixtureDigitalTopology = DGtal::Z3i::DT26_6;
+ using FixtureDigitalSet = DGtal::DigitalSetByAssociativeContainer<
+ DGtal::Z3i::Domain,
+ std::unordered_set>;
+ using FixtureMap = std::unordered_map;
+ using FixtureComplex = DGtal::VoxelComplex;
+
+ ///////////////////////////////////////////////////////////
+ // fixture data
+ FixtureComplex complex_fixture;
+ FixtureDigitalSet set_fixture;
+ KSpace ks_fixture; // needed because ConstAlias in CC constructor.
+ ///////////////////////////////////////////////////////////
+
+ ///////////////////////////////////////////////////////////
+ // Constructor
+ ///////////////////////////////////////////////////////////
+ Fixture_X_with_tight_bounds() : complex_fixture(ks_fixture), set_fixture(create_set()) {
+ create_complex_from_set(set_fixture);
+ }
+
+ ///////////////////////////////////////////////////////////
+ // Function members
+ ///////////////////////////////////////////////////////////
+ FixtureDigitalSet create_set() {
+ using namespace DGtal;
+
+ // with_tight_bounds
+ Point p1(-6, -6, -1);
+ Point p2(6, 6, 1);
+ Domain domain(p1, p2);
+
+ FixtureDigitalSet a_set(domain);
+ std::vector center_set;
+ center_set.reserve(9);
+
+ Point c00(0, 0, 0);
+ center_set.push_back(c00);
+ Point c01x(-1, 0, 0);
+ center_set.push_back(c01x);
+ Point c10x(1, 0, 0);
+ center_set.push_back(c10x);
+ Point c02x(-2, 0, 0);
+ center_set.push_back(c02x);
+ Point c20x(2, 0, 0);
+ center_set.push_back(c20x);
+
+ Point c01y(0, -1, 0);
+ center_set.push_back(c01y);
+ Point c10y(0, 1, 0);
+ center_set.push_back(c10y);
+ Point c02y(0, -2, 0);
+ center_set.push_back(c02y);
+ Point c20y(0, 2, 0);
+ center_set.push_back(c20y);
+
+ Point z_pos(0, 0, 1);
+ int branch_length(4);
+ std::vector diagonals;
+ diagonals.reserve(6);
+ for (const auto &p : center_set) {
+ diagonals.clear();
+ for (int l = 0; l <= branch_length; ++l) {
+ diagonals.push_back({l, l, 0});
+ diagonals.push_back({l, -l, 0});
+ diagonals.push_back({-l, l, 0});
+ diagonals.push_back({-l, -l, 0});
+ for (int z = -1; z <= 1; ++z)
+ for (const auto &d : diagonals)
+ a_set.insert(p + d + (z * z_pos));
+ }
+ }
+
+ return a_set;
+ }
+
+ FixtureComplex &create_complex_from_set(FixtureDigitalSet &input_set) {
+
+ ks_fixture.init(input_set.domain().lowerBound(),
+ input_set.domain().upperBound(), true);
+ complex_fixture = FixtureComplex(ks_fixture);
+ complex_fixture.construct(input_set);
+ return complex_fixture;
+ }
+};
+
+TEST_CASE_METHOD(Fixture_X_with_tight_bounds, "splitComplex", "[parallel]") {
+ using namespace DGtal::functions;
+ auto & vc = complex_fixture;
+ // CHECK(vc.nbCells(0) == 528);
+ // CHECK(vc.nbCells(1) == 1276);
+ // CHECK(vc.nbCells(2) == 1016);
+ // CHECK(vc.nbCells(3) == 267);
+ // Modify data of a cell to check if data is copied to sub complexes
+ Point p{2, 0, 0};
+ const auto modified_cell = vc.space().uSpel(p);
+ const auto modified_cell_dim = vc.space().uDim(modified_cell);
+ auto it_cell= vc.findCell(modified_cell_dim, modified_cell);
+ CHECK(it_cell != vc.end(modified_cell_dim));
+ DGtal::uint32_t magic_number = 15;
+ it_cell->second.data = magic_number;
+ SECTION("SplitComplex") {
+ size_t requested_number_of_splits = 2;
+ trace.beginBlock("SplitComplex");
+ trace.info() << "lowerBound" << vc.space().lowerBound() << std::endl;
+ trace.info() << "upperBound" << vc.space().upperBound() << std::endl;
+ auto out = splitComplex(vc, requested_number_of_splits);
+ auto & sub_complexes = out.sub_complexes;
+ CHECK(out.number_of_splits == requested_number_of_splits);
+ CHECK(sub_complexes.size() == requested_number_of_splits);
+ {
+ size_t sub_index = 0;
+ auto & sc = sub_complexes[sub_index];
+ CHECK(sc.nbCells(0) != 0);
+ CHECK(sc.nbCells(1) != 0);
+ CHECK(sc.nbCells(2) != 0);
+ CHECK(sc.nbCells(3) != 0);
+ trace.info() << "lowerBound S0" << sc.space().lowerBound() << std::endl;
+ trace.info() << "upperBound S0" << sc.space().upperBound() << std::endl;
+ CHECK(sc.space().lowerBound() == typename KSpace::Point(-6,-6,-1));
+ CHECK(sc.space().upperBound() == typename KSpace::Point(-1,6,1));
+ CHECK(sc.space().lowerBound() == out.splits_bounds[sub_index].lowerBound);
+ CHECK(sc.space().upperBound() == out.splits_bounds[sub_index].upperBound);
+ }
+ {
+ size_t sub_index = 1;
+ auto & sc = sub_complexes[sub_index];
+ CHECK(sc.nbCells(0) != 0);
+ CHECK(sc.nbCells(1) != 0);
+ CHECK(sc.nbCells(2) != 0);
+ CHECK(sc.nbCells(3) != 0);
+ trace.info() << "lowerBound S1" << sc.space().lowerBound() << std::endl;
+ trace.info() << "upperBound S1" << sc.space().upperBound() << std::endl;
+ CHECK(sc.space().lowerBound() == typename KSpace::Point(0,-6,-1));
+ CHECK(sc.space().upperBound() == typename KSpace::Point(6,6,1));
+ CHECK(sc.space().lowerBound() == out.splits_bounds[sub_index].lowerBound);
+ CHECK(sc.space().upperBound() == out.splits_bounds[sub_index].upperBound);
+ // Check cell exist in sub_complex
+ auto sc_it_cell = sc.findCell(it_cell->first);
+ CHECK(sc_it_cell != sc.end(modified_cell_dim));
+ // Check data is copied into sub_complexes
+ CHECK(sc[modified_cell].data == magic_number);
+ }
+ const size_t row_voxels_sc0 = 7;
+ CHECK(sub_complexes[0].nbCells(0) + sub_complexes[1].nbCells(0) ==
+ vc.nbCells(0) + (row_voxels_sc0 + 1) * 4);
+ CHECK(sub_complexes[0].nbCells(1) + sub_complexes[1].nbCells(1) ==
+ vc.nbCells(1) + ( (row_voxels_sc0 + 1) * 3 + row_voxels_sc0 * 4));
+ CHECK(sub_complexes[0].nbCells(2) + sub_complexes[1].nbCells(2) ==
+ vc.nbCells(2) + row_voxels_sc0 * 3);
+ CHECK(sub_complexes[0].nbCells(3) + sub_complexes[1].nbCells(3) ==
+ vc.nbCells(3));
+ trace.endBlock();
+ }
+ SECTION("getBorderVoxels") {
+ trace.beginBlock("getBorderVoxels");
+ size_t requested_number_of_splits = 4;
+ auto out = splitComplex(vc, requested_number_of_splits);
+ auto & sub_complexes = out.sub_complexes;
+ const auto lowerBound_to_ignore = vc.space().lowerBound();
+ const auto upperBound_to_ignore = vc.space().upperBound();
+ typename FixtureComplex::Point wide_point = {2,2,2};
+ std::vector> borders;
+ {
+ size_t sub_index = 0;
+ auto & sc = sub_complexes[sub_index];
+ const auto lowerBound = sc.space().lowerBound();
+ const auto upperBound = sc.space().upperBound();
+ auto border_iterators =
+ getBorderVoxels(sc, lowerBound, upperBound, &wide_point,
+ &lowerBound_to_ignore, &upperBound_to_ignore);
+ CHECK(border_iterators.size() == 30);
+ borders.push_back(border_iterators);
+ }
+ {
+ size_t sub_index = 1;
+ auto & sc = sub_complexes[sub_index];
+ const auto lowerBound = sc.space().lowerBound();
+ const auto upperBound = sc.space().upperBound();
+ auto border_iterators =
+ getBorderVoxels(sc, lowerBound, upperBound, &wide_point,
+ &lowerBound_to_ignore, &upperBound_to_ignore);
+ CHECK(border_iterators.size() == 30);
+ borders.push_back(border_iterators);
+ }
+ trace.endBlock();
+ }
+}
+////////////////////// end splitComplex ////////////////////////
diff --git a/tests/topology/testVoxelComplex.cpp b/tests/topology/testVoxelComplex.cpp
index d6add4efb9..2097c6906b 100644
--- a/tests/topology/testVoxelComplex.cpp
+++ b/tests/topology/testVoxelComplex.cpp
@@ -803,7 +803,7 @@ TEST_CASE_METHOD(Fixture_isthmus, "Thin complex", "[isthmus][thin][function]") {
SECTION("with skelIsthmus") {
auto vc_new = asymetricThinningScheme(
vc, selectRandom, skelIsthmus);
- CHECK(vc_new.nbCells(3) == 3);
+ CHECK(vc_new.nbCells(3) == 4);
}
}
//
@@ -1090,30 +1090,9 @@ TEST_CASE_METHOD(Fixture_X, "X DistanceMap", "[x][distance][thin]") {
auto vc_new = persistenceAsymetricThinningScheme(
vc, selectDistMax, oneIsthmusTable,
persistence, verbose);
- // SECTION( "visualize the thining" ){
- // int argc(1);
- // char** argv(nullptr);
- // QApplication app(argc, argv);
- // Viewer3D<> viewer(ks_fixture);
- // viewer.show();
- //
- // viewer.setFillColor(Color(200, 200, 200, 100));
- // for ( auto it = vc_new.begin(3); it!= vc_new.end(3); ++it )
- // viewer << it->first;
- //
- // // All kspace voxels
- // viewer.setFillColor(Color(40, 40, 40, 10));
- // for ( auto it = vc.begin(3); it!= vc.end(3); ++it )
- // viewer << it->first;
- //
- // viewer << Viewer3D<>::updateDisplay;
- // app.exec();
- // }
}
-
+ trace.endBlock();
}
-
- trace.endBlock();
}
// REQUIRE(vc_new.nbCells(3) == 38);