ENH: make the ImageRegionTileMapSplitter multidimensional

Code/Common/otbImageRegionTileMapSplitter.h

@@ -73,10 +73,10 @@ class ITK_EXPORT ImageRegionTileMapSplitter: public itk::ImageRegionSplitter<VIm
 {
 public:
   /** Standard class typedefs. */
-  typedef ImageRegionTileMapSplitter              Self;
-  typedef itk::ImageRegionSplitter<VImageDimension>  Superclass;
-  typedef itk::SmartPointer<Self>  Pointer;
-  typedef itk::SmartPointer<const Self>  ConstPointer;
+  typedef ImageRegionTileMapSplitter                Self;
+  typedef itk::ImageRegionSplitter<VImageDimension> Superclass;
+  typedef itk::SmartPointer<Self>                   Pointer;
+  typedef itk::SmartPointer<const Self>             ConstPointer;
 
   /** Method for creation through the object factory. */
   itkNewMacro(Self);
@@ -94,12 +94,12 @@ public:
   }
 
   /** Index typedef support. An index is used to access pixel values. */
-  typedef itk::Index<VImageDimension>  IndexType;
-  typedef typename IndexType::IndexValueType  IndexValueType;
+  typedef itk::Index<VImageDimension>        IndexType;
+  typedef typename IndexType::IndexValueType IndexValueType;
 
   /** Size typedef support. A size is used to define region bounds. */
-  typedef itk::Size<VImageDimension>  SizeType;
-  typedef typename SizeType::SizeValueType  SizeValueType;
+  typedef itk::Size<VImageDimension>       SizeType;
+  typedef typename SizeType::SizeValueType SizeValueType;
 
   /** Region typedef support.   */
   typedef itk::ImageRegion<VImageDimension> RegionType;
@@ -120,7 +120,7 @@ public:
                               const RegionType &region);
 
 protected:
-  ImageRegionTileMapSplitter() {}
+  ImageRegionTileMapSplitter(): m_AlignStep(256){}
   virtual ~ImageRegionTileMapSplitter() {}
   void PrintSelf(std::ostream& os, itk::Indent indent) const;
 
@@ -128,6 +128,8 @@ private:
   ImageRegionTileMapSplitter(const ImageRegionTileMapSplitter&); //purposely not implemented
   void operator=(const ImageRegionTileMapSplitter&); //purposely not implemented
 
+  unsigned int m_SplitsPerDimension[VImageDimension];
+  unsigned int m_AlignStep;
 };
 
 
@@ -138,4 +140,3 @@ private:
 #endif
 
 #endif
-

Code/Common/otbImageRegionTileMapSplitter.txx

-
 /*=========================================================================
 
   Program:   ORFEO Toolbox
@@ -21,6 +20,7 @@
 
 #include "otbImageRegionTileMapSplitter.h"
 #include "otbMath.h"
+#include "otbMacro.h"
 
 namespace otb
 {
@@ -35,25 +35,40 @@ ImageRegionTileMapSplitter<VImageDimension>
 {
   int splitAxis;
   const SizeType &regionSize = region.GetSize();
+  const IndexType &regionIndex = region.GetIndex();
+
+  // requested number of splits per dimension
+  unsigned int numPieces = 1;
 
   // split on the outermost dimension available
   splitAxis = VImageDimension - 1;
-  while (regionSize[splitAxis] == 1)
-  {
-    --splitAxis;
-    if (splitAxis < 0)
-    { // cannot split
-      itkDebugMacro("  Cannot Split");
-      return 1;
-    }
-  }
 
   // determine the actual number of pieces that will be generated
-  SizeValueType range = regionSize[splitAxis];
-  int valuesPerPiece = static_cast<int>(vcl_ceil(static_cast<double>(range/((double)requestedNumber*256)) ) ) * 256;
-  int maxPieceUsed = static_cast<int>(vcl_ceil(static_cast<double>(range/(double)valuesPerPiece)));
-
-  return maxPieceUsed + 1;
+  for (unsigned int j = VImageDimension; j > 0; --j)
+    {
+//    otbMsgDevMacro(<< "*** Dimension: " << j-1);
+    SizeValueType range = regionSize[j - 1];
+    unsigned long int remainingToDo = vcl_ceil((double) requestedNumber / numPieces);
+//    unsigned long int remainingToDo = requestedNumber - numPieces;
+    unsigned int maxPieces = (regionIndex[j - 1] + regionSize[j - 1] - 1) / m_AlignStep - regionIndex[j - 1]
+        / m_AlignStep + 1;
+    unsigned int stepPerPiece = 1;
+    if (remainingToDo < maxPieces)
+      {
+      stepPerPiece = static_cast<int> (vcl_floor(static_cast<double> (maxPieces) / remainingToDo));
+      if ((remainingToDo-1)*(stepPerPiece+1) < maxPieces )
+        {
+        stepPerPiece += 1;
+        }
+      }
+    unsigned int maxPieceUsed = static_cast<int> (vcl_ceil(static_cast<double> (maxPieces / (double) stepPerPiece)));
+    m_SplitsPerDimension[j - 1] = maxPieceUsed;
+//    otbMsgDevMacro("*** maxPieces stepPerPiece maxPieceUsed " << maxPieces
+//                      << " " << stepPerPiece << " " << maxPieceUsed);
+    numPieces *= maxPieceUsed;
+    }
+//  otbMsgDevMacro("*** numPieces " << numPieces);
+  return numPieces;
 }
 
 
@@ -67,7 +82,7 @@ ImageRegionTileMapSplitter<VImageDimension>
 {
   int splitAxis;
   RegionType splitRegion;
-  IndexType splitIndex;
+  IndexType splitIndex, regionIndex;
   SizeType splitSize, regionSize;
 
   // Initialize the splitRegion to the requested region
@@ -76,50 +91,49 @@ ImageRegionTileMapSplitter<VImageDimension>
   splitSize = splitRegion.GetSize();
 
   regionSize = region.GetSize();
+  regionIndex = region.GetIndex();
 
-  // split on the outermost dimension available
-  splitAxis = VImageDimension - 1;
-  while (regionSize[splitAxis] == 1)
-  {
-    --splitAxis;
-    if (splitAxis < 0)
-    { // cannot split
-      itkDebugMacro("  Cannot Split");
-      return splitRegion;
+  unsigned int numPieces = GetNumberOfSplits(region, numberOfPieces);
+  if (i > numPieces)
+    {
+    itkDebugMacro("  Cannot Split");
+    return splitRegion;
     }
-  }
-
-  // determine the actual number of pieces that will be generated
-  SizeValueType range = regionSize[splitAxis];
-  int valuesPerPiece = static_cast<int>((int)vcl_ceil(static_cast<double>(range/((double)numberOfPieces*256)))) * 256;
-  int maxPieceUsed = static_cast<int>(vcl_ceil(static_cast<double>(range/(double)valuesPerPiece)));
-
-  // Split the region
-  if ((int) i == maxPieceUsed-1)
-  {
-    splitIndex[splitAxis] += i*valuesPerPiece;
-    // last piece needs to process the "rest" dimension being split
-//     splitSize[splitAxis] = splitSize[splitAxis] - i*valuesPerPiece;
-    splitSize[splitAxis] = regionSize[splitAxis] - splitIndex[splitAxis];
-  }
-  if (((int) i == 0) && ((int) i != maxPieceUsed-1))
-  {
-    //First piece may not contain a whole piece
-    //splitIndex stay at splitRegion.GetIndex();
-    splitSize[splitAxis] = ((static_cast<int>(vcl_floor(static_cast<double>(splitIndex[splitAxis]/valuesPerPiece))))+1)
-                           *valuesPerPiece - splitIndex[splitAxis];
-  }
-  if ((int) i < maxPieceUsed-1)
-  {
-    splitIndex[splitAxis] += i*valuesPerPiece;
-    splitSize[splitAxis] = valuesPerPiece;
-  }
 
+  unsigned int stackSize = 1;
+  for (unsigned int j = 0; j < VImageDimension; ++j)
+    {
+    unsigned int slicePos = (i % (stackSize * m_SplitsPerDimension[j])) / stackSize;
+    stackSize *= m_SplitsPerDimension[j];
+
+    unsigned int generalSplitSize = vcl_ceil(static_cast<double> (regionSize[j]) / (m_SplitsPerDimension[j]
+        * m_AlignStep)) * m_AlignStep;//FIXME check
+    if (slicePos == 0)
+      {
+      splitIndex[j] = regionIndex[j];
+      }
+    else
+      {
+      splitIndex[j] = (regionIndex[j] / generalSplitSize + slicePos) * generalSplitSize;
+      }
+    if (slicePos == 0)
+      {
+      splitSize[j] = generalSplitSize - (regionIndex[j] % generalSplitSize);
+      }
+    else if (slicePos == m_SplitsPerDimension[j] - 1)
+      {
+      splitSize[j] = regionSize[j] - (generalSplitSize - (regionIndex[j] % generalSplitSize))
+          - (m_SplitsPerDimension[j] - 2) * generalSplitSize;
+      }
+    else
+      {
+      splitSize[j] = generalSplitSize;
+      }
+    }
 
   // set the split region ivars
-  splitRegion.SetIndex( splitIndex );
-  splitRegion.SetSize( splitSize );
-
+  splitRegion.SetIndex(splitIndex);
+  splitRegion.SetSize(splitSize);
 
   itkDebugMacro("  Split Piece: " << splitRegion );
 

Testing/Code/Common/CMakeLists.txt

@@ -723,11 +723,7 @@ ADD_TEST(coTvImageRegionTileMapSplitter ${COMMON_TESTS8}
         ${BASELINE_FILES}/coImageRegionTileMapSplitter.txt
                 ${TEMP}/coImageRegionTileMapSplitter.txt
     otbImageRegionTileMapSplitter
-    45 45      # Index
-    1000 1500  # Size
-    10         # Nb split
-    2          # Number of slipt
-    ${TEMP}/coImageRegionTileMapSplitter.txt
+        ${TEMP}/coImageRegionTileMapSplitter.txt
 )
 
 # -------------  otb::ImageOfVectorsToMonoChannelExtractROI ----------------------------

Testing/Code/Common/otbImageRegionTileMapSplitter.cxx

@@ -20,46 +20,128 @@
 
 int otbImageRegionTileMapSplitter( int argc, char * argv[] )
 {
+
+
   const int Dimension = 2;
-  typedef otb::ImageRegionTileMapSplitter< Dimension >  FilterType;
-  typedef FilterType::IndexType IndexType;
-  typedef FilterType::SizeType SizeType;
-  typedef FilterType::RegionType RegionType;
+  typedef otb::ImageRegionTileMapSplitter< Dimension > FilterType;
+  typedef FilterType::IndexType                        IndexType;
+  typedef FilterType::SizeType                         SizeType;
+  typedef FilterType::RegionType                       RegionType;
+
+  std::ofstream outfile(argv[1]);
+  FilterType::Pointer filter = FilterType::New();
 
+  RegionType region, region2;
+  unsigned int nb, nbSplitTheoric, nbAsked;
   IndexType index;
-  index[0] = atoi(argv[1]);
-  index[1] = atoi(argv[2]);
   SizeType size;
-  size[0] = atoi(argv[3]);
-  size[1] = atoi(argv[4]);
-  unsigned int nbSplitTheoric(atoi(argv[5]));
-  unsigned int nbAsked(atoi(argv[6]));
-  const char * outfname(argv[7]);
 
+  //Case 1
+  index[0] = 45;
+  index[1] = 45;
+  size[0] = 1000;
+  size[1] = 1500;
+  nbSplitTheoric = 10;
+  nbAsked = 2;
 
-  RegionType region;
+  region.SetSize(size);
+  region.SetIndex(index);
+
+  nb = filter->GetNumberOfSplits(region,nbSplitTheoric);
+  region2 = filter->GetSplit(nbAsked,nb,region);
+
+  outfile << "\nCase 1 \n";
+  outfile << "Input region: "<<region<<std::endl;
+  outfile << "Input NumberOfSplits: "<<nbSplitTheoric<<std::endl;
+  outfile << "Output GetNumberOfSplits: "<<nb<<std::endl;
+  outfile << "Output GetSplit("<<nbAsked<<","<<nb<<",input region): "<<std::endl;
+  outfile << "Output region: "<< region2<<std::endl;
+
+  //Case 2
+  index[0] = 45;
+  index[1] = 45;
+  size[0] = 1048576;
+  size[1] = 1048576;
+  nbSplitTheoric = 16777216;
+  nbAsked = 2;
 
   region.SetSize(size);
   region.SetIndex(index);
 
-  FilterType::Pointer filter = FilterType::New();
+  nb = filter->GetNumberOfSplits(region,nbSplitTheoric);
+  region2 = filter->GetSplit(nbAsked,nb,region);
+
+  outfile << "\nCase 2 \n";
+  outfile << "Input region: "<<region<<std::endl;
+  outfile << "Input NumberOfSplits: "<<nbSplitTheoric<<std::endl;
+  outfile << "Output GetNumberOfSplits: "<<nb<<std::endl;
+  outfile << "Output GetSplit("<<nbAsked<<","<<nb<<",input region): "<<std::endl;
+  outfile << "Output region: "<< region2<<std::endl;
+
+  //Case 3
+  index[0] = 45;
+  index[1] = 45;
+  size[0] = 1048576;
+  size[1] = 1048576;
+  nbSplitTheoric = 23;
+  nbAsked = 4;
+
+  region.SetSize(size);
+  region.SetIndex(index);
+
+  nb = filter->GetNumberOfSplits(region,nbSplitTheoric);
+  region2 = filter->GetSplit(nbAsked,nb,region);
 
-  unsigned int nb = filter->GetNumberOfSplits(region,nbSplitTheoric);
+  outfile << "\nCase 3 \n";
+  outfile << "Input region: "<<region<<std::endl;
+  outfile << "Input NumberOfSplits: "<<nbSplitTheoric<<std::endl;
+  outfile << "Output GetNumberOfSplits: "<<nb<<std::endl;
+  outfile << "Output GetSplit("<<nbAsked<<","<<nb<<",input region): "<<std::endl;
+  outfile << "Output region: "<< region2<<std::endl;
   
-  RegionType region2 = filter->GetSplit(nbAsked,nb,region);
+  //Case 4
+  index[0] = 45;
+  index[1] = 45;
+  size[0] = 1048576;
+  size[1] = 1024;
+  nbSplitTheoric = 16777216;
+  nbAsked = 16387;
 
-  std::ofstream outfile(outfname);
+  region.SetSize(size);
+  region.SetIndex(index);
 
+  nb = filter->GetNumberOfSplits(region,nbSplitTheoric);
+  region2 = filter->GetSplit(nbAsked,nb,region);
 
+  outfile << "\nCase 4 \n";
   outfile << "Input region: "<<region<<std::endl;
   outfile << "Input NumberOfSplits: "<<nbSplitTheoric<<std::endl;
   outfile << "Output GetNumberOfSplits: "<<nb<<std::endl;
   outfile << "Output GetSplit("<<nbAsked<<","<<nb<<",input region): "<<std::endl;
-  outfile << "    "<<region2<<std::endl;
+  outfile << "Output region: "<< region2<<std::endl;
+
+  //Case 5
+  index[0] = 0;
+  index[1] = 0;
+  size[0] = 513;
+  size[1] = 5376;
+  nbSplitTheoric = 8;
+  nbAsked = 9;
+
+  region.SetSize(size);
+  region.SetIndex(index);
+
+  nb = filter->GetNumberOfSplits(region,nbSplitTheoric);
+  region2 = filter->GetSplit(nbAsked,nb,region);
+
+  outfile << "\nCase 5 \n";
+  outfile << "Input region: "<<region<<std::endl;
+  outfile << "Input NumberOfSplits: "<<nbSplitTheoric<<std::endl;
+  outfile << "Output GetNumberOfSplits: "<<nb<<std::endl;
+  outfile << "Output GetSplit("<<nbAsked<<","<<nb<<",input region): "<<std::endl;
+  outfile << "Output region: "<< region2<<std::endl;
 
   outfile.close();
 
   return EXIT_SUCCESS;
 }
-
-