BUG: fix handling of negative and non-zero region indices in rasterizer

Avoid creation of vector of indices, instead directly set pixel values.
This allows correct 3D index to linear (buffer) index conversion.
It also reduces memory used, and shortens running time.

Closes #2884.

Author: dzenanz

Modules/Core/Mesh/include/itkTriangleMeshToBinaryImageFilter.h

@@ -129,7 +129,6 @@ class ITK_TEMPLATE_EXPORT TriangleMeshToBinaryImageFilter : public ImageSource<T
   using PointVector = std::vector<PointType>;
   using PointArray = std::vector<std::vector<PointType>>;
 
-  using StencilIndexVector = std::vector<int>;
   /** Spacing (size of a pixel) of the output image. The
    * spacing is the geometric distance between image samples.
    * It is stored internally as double, but may be set from
@@ -245,8 +244,6 @@ class ITK_TEMPLATE_EXPORT TriangleMeshToBinaryImageFilter : public ImageSource<T
 
   DirectionType m_Direction;
 
-  StencilIndexVector m_StencilIndex;
-
   void
   PrintSelf(std::ostream & os, Indent indent) const override;
 

Modules/Core/Mesh/include/itkTriangleMeshToBinaryImageFilter.hxx

@@ -185,46 +185,6 @@ TriangleMeshToBinaryImageFilter<TInputMesh, TOutputImage>::GenerateData()
 
   RasterizeTriangles();
 
-  using myIteratorType = itk::ImageRegionIteratorWithIndex<OutputImageType>;
-
-  myIteratorType it(OutputImage, OutputImage->GetLargestPossibleRegion());
-
-  int DataIndex = 0;
-  int StencilId = 0;
-  it.GoToBegin();
-
-  size_t n = m_StencilIndex.size();
-  if (n == 0)
-  {
-    itkWarningMacro(<< "No Image Indices Found.");
-  }
-  else
-  {
-    int StencilMin = m_StencilIndex[0];
-    int StencilMax = m_StencilIndex[n - 1];
-
-    while (!it.IsAtEnd())
-    {
-      if (DataIndex >= StencilMin && DataIndex <= StencilMax)
-      {
-        if (DataIndex == m_StencilIndex[StencilId])
-        {
-          it.Set(m_InsideValue);
-          StencilId++;
-        }
-        else
-        {
-          it.Set(m_OutsideValue);
-        }
-      }
-      else
-      {
-        it.Set(m_OutsideValue);
-      }
-      DataIndex++;
-      ++it;
-    }
-  }
   itkDebugMacro(<< "TriangleMeshToBinaryImageFilter::Update() finished");
 } // end update function
 
@@ -459,8 +419,9 @@ TriangleMeshToBinaryImageFilter<TInputMesh, TOutputImage>::RasterizeTriangles()
     ++cellIt;
   }
 
-  // create the equivalent of vtkStencilData from our zymatrix
-  m_StencilIndex.clear(); // prevent corruption of the filter in later updates
+  OutputImagePointer outputImage = this->GetOutput();
+  outputImage->FillBuffer(m_OutsideValue);
+
   for (int z = extent[4]; z <= extent[5]; ++z)
   {
     for (int y = extent[2]; y <= extent[3]; ++y)
@@ -529,9 +490,14 @@ TriangleMeshToBinaryImageFilter<TInputMesh, TOutputImage>::RasterizeTriangles()
 
         if (x2 >= x1)
         {
+          IndexType ind;
+          ind[1] = y;
+          ind[2] = z;
           for (int idX = x1; idX <= x2; ++idX)
           {
-            m_StencilIndex.push_back(idX + y * m_Size[0] + z * m_Size[0] * m_Size[1]);
+            // TODO: check whether replacing this for loop by ImageScanlineIterator is faster
+            ind[0] = idX;
+            outputImage->SetPixel(ind, m_InsideValue);
           }
         }
         // next x1 value must be at least x2+1