BUG: Mantis-1327: fix the sampling grid used by StreamingShrinkImageFilter

Modules/Filtering/ImageManipulation/include/otbStreamingShrinkImageFilter.h

@@ -78,7 +78,6 @@ public:
                               const RegionType& region) ITK_OVERRIDE;
 
   itkGetMacro(TileSizeAlignment, unsigned int);
-  itkSetMacro(TileSizeAlignment, unsigned int);
 
   itkGetMacro(TileDimension, unsigned int);
 
@@ -246,6 +245,9 @@ private:
 
   /** The shrink factor */
   unsigned int m_ShrinkFactor;
+
+  /** The offset to get the cell center */
+  IndexType m_Offset;
 }; // end of class PersistentStatisticsVectorImageFilter
 
 

Modules/Filtering/ImageManipulation/include/otbStreamingShrinkImageFilter.txx

@@ -43,7 +43,7 @@ StreamingShrinkStreamingManager<TImage>::PrepareStreaming( itk::DataObject * inp
 {
   typedef otb::StreamingShrinkImageRegionSplitter TileSplitterType;
   TileSplitterType::Pointer splitter = TileSplitterType::New();
-  splitter->SetTileSizeAlignment(m_ShrinkFactor);
+  splitter->SetShrinkFactor(m_ShrinkFactor);
   this->m_Splitter = splitter;
 
   unsigned long nbDivisions = this->EstimateOptimalNumberOfDivisions(input, region, 0);
@@ -125,6 +125,9 @@ PersistentShrinkImageFilter<TInputImage, TOutputImage>
                                            inputSpacing = inputPtr->GetSpacing();
   const typename InputImageType::SizeType& inputSize
     = inputPtr->GetLargestPossibleRegion().GetSize();
+  const typename InputImageType::IndexType& inputIndex
+    = inputPtr->GetLargestPossibleRegion().GetIndex();
+  typename InputImageType::IndexType startIndex;
 
   typename OutputImageType::SpacingType shrunkOutputSpacing;
   typename OutputImageType::RegionType  shrunkOutputLargestPossibleRegion;
@@ -134,12 +137,14 @@ PersistentShrinkImageFilter<TInputImage, TOutputImage>
 
   for (unsigned int i = 0; i < OutputImageType::ImageDimension; ++i)
     {
+    startIndex[i] = inputIndex[i] + (m_ShrinkFactor - 1) / 2;
+    if (m_ShrinkFactor > inputSize[i])
+      startIndex[i] = inputIndex[i] + (inputSize[i] - 1) / 2;
+    m_Offset[i] = startIndex[i] % m_ShrinkFactor;
     shrunkOutputSpacing[i] = inputSpacing[i] * static_cast<double>(m_ShrinkFactor);
     shrunkOutputSize[i] = inputSize[i] > m_ShrinkFactor ? inputSize[i] / m_ShrinkFactor : 1;
     
-    shrunkOutputOrigin[i] = inputPtr->GetOrigin()[i] + inputSpacing[i] *
-      (static_cast<double>(inputPtr->GetLargestPossibleRegion().GetIndex(i)) - 0.5)
-      + shrunkOutputSpacing[i] * 0.5;
+    shrunkOutputOrigin[i] = inputPtr->GetOrigin()[i] + inputSpacing[i] * startIndex[i];
 
     // we choose to output a region with a start index [0,0]
     // the origin is set accordingly
@@ -184,12 +189,12 @@ PersistentShrinkImageFilter<TInputImage, TOutputImage>
     {
     const IndexType& inIndex = inIt.GetIndex();
     // TODO the pixel value should be taken near the centre of the cell, not at the corners
-    if (inIndex[0] % m_ShrinkFactor == 0
-        && inIndex[1] % m_ShrinkFactor == 0 )
+    if ((inIndex[0] - m_Offset[0]) % m_ShrinkFactor == 0
+        && (inIndex[1] - m_Offset[1]) % m_ShrinkFactor == 0 )
       {
       IndexType shrunkIndex;
-      shrunkIndex[0] = inIndex[0] / m_ShrinkFactor;
-      shrunkIndex[1] = inIndex[1] / m_ShrinkFactor;
+      shrunkIndex[0] = (inIndex[0] - m_Offset[0]) / m_ShrinkFactor;
+      shrunkIndex[1] = (inIndex[1] - m_Offset[1]) / m_ShrinkFactor;
       if (m_ShrunkOutput->GetLargestPossibleRegion().IsInside(shrunkIndex))
         m_ShrunkOutput->SetPixel(shrunkIndex, inIt.Get());
       }

Modules/Filtering/ImageManipulation/src/otbStreamingShrinkImageFilter.cxx

@@ -27,23 +27,24 @@ StreamingShrinkImageRegionSplitter
   unsigned int theoricalNbPixelPerTile = region.GetNumberOfPixels() / requestedNumber;
   unsigned int theoricalTileDimension = static_cast<unsigned int> (vcl_sqrt(static_cast<double>(theoricalNbPixelPerTile)) );
 
-  // Take the previous multiple of m_TileSizeAlignment (eventually generate more splits than requested)
-  m_TileDimension = theoricalTileDimension / m_TileSizeAlignment * m_TileSizeAlignment;
+  // Take the previous multiple of m_ShrinkFactor (eventually generate more splits than requested)
+  m_TileDimension = theoricalTileDimension / m_ShrinkFactor * m_ShrinkFactor;
 
-  // Minimal tile size is m_TileSizeAlignment * m_TileSizeAlignment
-  if (m_TileDimension < m_TileSizeAlignment)
+  // Minimal tile size is m_ShrinkFactor * m_ShrinkFactor
+  if (m_TileDimension < m_ShrinkFactor)
     {
-    otbMsgDevMacro(<< "Using the minimal tile size : " << m_TileSizeAlignment << " * " << m_TileSizeAlignment);
-    m_TileDimension = m_TileSizeAlignment;
+    otbMsgDevMacro(<< "Using the minimal tile size : " << m_ShrinkFactor << " * " << m_ShrinkFactor);
+    m_TileDimension = m_ShrinkFactor;
     }
 
-  // Use the computed tile size, and generate (m_TileDimension * 1) tiles
   const SizeType&  regionSize = region.GetSize();
-  m_SplitsPerDimension[0] = (regionSize[0] + m_TileDimension - 1) / m_TileDimension;
-  m_SplitsPerDimension[1] = regionSize[1] / m_TileSizeAlignment;
+  // Compute the alignment of the sampling grid
+  m_TileSizeAlignment = (m_ShrinkFactor - 1) / 2;
+  if (m_ShrinkFactor > regionSize[1]) m_TileSizeAlignment = (regionSize[1] - 1)/2;
 
-  if (m_SplitsPerDimension[1] == 0)
-    m_SplitsPerDimension[1] = 1;
+  // Use the computed tile size, and generate (m_TileDimension * 1) tiles
+  m_SplitsPerDimension[0] = (regionSize[0] + m_TileDimension - 1) / m_TileDimension;
+  m_SplitsPerDimension[1] = (regionSize[1] - m_TileSizeAlignment - 1) / m_ShrinkFactor + 1;
 
   unsigned int numPieces = 1;
   for (unsigned int j = 0; j < ImageDimension; ++j)
@@ -82,7 +83,7 @@ StreamingShrinkImageRegionSplitter
 
   // Transform the split index to the actual coordinates
   splitRegion.SetIndex(0, region.GetIndex(0) + m_TileDimension * splitIndex[0]);
-  splitRegion.SetIndex(1, region.GetIndex(1) + m_TileSizeAlignment * splitIndex[1]);
+  splitRegion.SetIndex(1, region.GetIndex(1) + m_ShrinkFactor * splitIndex[1] + m_TileSizeAlignment);
 
   splitRegion.SetSize(0, m_TileDimension);
   splitRegion.SetSize(1, 1);
@@ -98,6 +99,7 @@ StreamingShrinkImageRegionSplitter
 ::PrintSelf(std::ostream& os, itk::Indent indent) const
 {
   Superclass::PrintSelf(os, indent);
+  os << indent << "ShrinkFactor       : " << m_ShrinkFactor << std::endl;
   os << indent << "SplitsPerDimension : " << m_SplitsPerDimension << std::endl;
   os << indent << "TileDimension      : " << m_TileDimension << std::endl;
   os << indent << "TileSizeAlignment  : " << m_TileSizeAlignment << std::endl;