ENH: Replace ExtractImageFilter with zero-copy sub-region extraction

include/rtkFDKConeBeamReconstructionFilter.h

@@ -23,8 +23,7 @@
 #include "rtkConfiguration.h"
 #include "rtkFDKBackProjectionImageFilter.h"
 #include "rtkFFTRampImageFilter.h"
-
-#include <itkExtractImageFilter.h>
+#include "rtkSubRegionViewImageFilter.h"
 
 namespace rtk
 {
@@ -38,8 +37,8 @@ namespace rtk
  * - rtk::FFTRampImageFilter for ramp filtering,
  * - rtk::FDKBackProjectionImageFilter for backprojection.
  * The input stack of projections is processed piece by piece (the size is
- * controlled with ProjectionSubsetSize) via the use of itk::ExtractImageFilter
- * to extract sub-stacks.
+ * controlled with ProjectionSubsetSize) by extracting sub-stacks directly
+ * from the input buffer pointer (zero-copy).
  *
  * \dot
  * digraph FDKConeBeamReconstructionFilter {
@@ -76,7 +75,7 @@ class ITK_TEMPLATE_EXPORT FDKConeBeamReconstructionFilter : public itk::InPlaceI
   using OutputImageType = TOutputImage;
 
   /** Typedefs of each subfilter of this composite filter */
-  using ExtractFilterType = itk::ExtractImageFilter<InputImageType, OutputImageType>;
+  using ExtractFilterType = rtk::SubRegionViewImageFilter<InputImageType>;
   using WeightFilterType = rtk::FDKWeightProjectionFilter<InputImageType, OutputImageType>;
   using RampFilterType = rtk::FFTRampImageFilter<OutputImageType, OutputImageType, TFFTPrecision>;
   using BackProjectionFilterType = rtk::FDKBackProjectionImageFilter<OutputImageType, OutputImageType>;

include/rtkFDKConeBeamReconstructionFilter.hxx

@@ -41,7 +41,6 @@ FDKConeBeamReconstructionFilter<TInputImage, TOutputImage, TFFTPrecision>::FDKCo
   m_RampFilter->SetInput(m_WeightFilter->GetOutput());
 
   // Default parameters
-  m_ExtractFilter->SetDirectionCollapseToSubmatrix();
   m_WeightFilter->InPlaceOn();
 
   // Default to one projection per subset when FFTW is not available
@@ -92,16 +91,19 @@ FDKConeBeamReconstructionFilter<TInputImage, TOutputImage, TFFTPrecision>::Gener
 
   // We only set the first sub-stack at that point, the rest will be
   // requested in the GenerateData function
-  typename ExtractFilterType::InputImageRegionType projRegion;
+  typename InputImageType::RegionType projRegion;
   projRegion = this->GetInput(1)->GetLargestPossibleRegion();
   unsigned int firstStackSize = std::min(m_ProjectionSubsetSize, (unsigned int)projRegion.GetSize(Dimension - 1));
   projRegion.SetSize(Dimension - 1, firstStackSize);
+
+  m_ExtractFilter->SetInput(this->GetInput(1));
   m_ExtractFilter->SetExtractionRegion(projRegion);
+  if (rtk::IsContiguousSubRegion(this->GetInput(1), projRegion))
+    m_WeightFilter->InPlaceOff();
 
   // Run composite filter update
   m_BackProjectionFilter->SetInput(0, this->GetInput(0));
   m_BackProjectionFilter->SetInPlace(this->GetInPlace());
-  m_ExtractFilter->SetInput(this->GetInput(1));
   m_BackProjectionFilter->UpdateOutputInformation();
 
   // Update output information
@@ -117,13 +119,11 @@ FDKConeBeamReconstructionFilter<TInputImage, TOutputImage, TFFTPrecision>::Gener
 {
   const unsigned int Dimension = this->InputImageDimension;
 
-  // The backprojection works on a small stack of projections, not the full stack
-  typename ExtractFilterType::InputImageRegionType subsetRegion;
+  typename InputImageType::RegionType subsetRegion;
   subsetRegion = this->GetInput(1)->GetLargestPossibleRegion();
   unsigned int nProj = subsetRegion.GetSize(Dimension - 1);
+  unsigned int baseIndex = subsetRegion.GetIndex(Dimension - 1);
 
-  // The progress accumulator tracks the progress of the pipeline
-  // Each filter is equally weighted across all iterations of the stack
   auto progress = itk::ProgressAccumulator::New();
   progress->SetMiniPipelineFilter(this);
   auto frac = (1.0f / 3) / itk::Math::ceil(double(nProj) / m_ProjectionSubsetSize);
@@ -133,22 +133,23 @@ FDKConeBeamReconstructionFilter<TInputImage, TOutputImage, TFFTPrecision>::Gener
 
   for (unsigned int i = 0; i < nProj; i += m_ProjectionSubsetSize)
   {
-    // After the first bp update, we need to use its output as input.
+    subsetRegion.SetIndex(Dimension - 1, baseIndex + i);
+    subsetRegion.SetSize(Dimension - 1, std::min(m_ProjectionSubsetSize, nProj - i));
+    m_ExtractFilter->SetInput(this->GetInput(1));
+    m_ExtractFilter->SetExtractionRegion(subsetRegion);
+    if (rtk::IsContiguousSubRegion(this->GetInput(1), subsetRegion))
+      m_WeightFilter->InPlaceOff();
+
     if (i)
     {
       typename TInputImage::Pointer pimg = m_BackProjectionFilter->GetOutput();
       pimg->DisconnectPipeline();
       m_BackProjectionFilter->SetInput(pimg);
 
-      // Change projection subset
-      subsetRegion.SetIndex(Dimension - 1, i);
-      subsetRegion.SetSize(Dimension - 1, std::min(m_ProjectionSubsetSize, nProj - i));
-      m_ExtractFilter->SetExtractionRegion(subsetRegion);
-
-      // This is required to reset the full pipeline
       m_BackProjectionFilter->GetOutput()->UpdateOutputInformation();
       m_BackProjectionFilter->GetOutput()->PropagateRequestedRegion();
     }
+
     m_BackProjectionFilter->Update();
   }
 

include/rtkFDKVarianceReconstructionFilter.h

@@ -23,8 +23,7 @@
 #include "rtkConfiguration.h"
 #include "rtkFDKBackProjectionImageFilter.h"
 #include "rtkFFTVarianceRampImageFilter.h"
-
-#include <itkExtractImageFilter.h>
+#include "rtkSubRegionViewImageFilter.h"
 
 namespace rtk
 {
@@ -68,7 +67,7 @@ class ITK_TEMPLATE_EXPORT FDKVarianceReconstructionFilter : public itk::InPlaceI
   using OutputImageType = TOutputImage;
 
   /** Typedefs of each subfilter of this composite filter */
-  using ExtractFilterType = itk::ExtractImageFilter<InputImageType, OutputImageType>;
+  using ExtractFilterType = rtk::SubRegionViewImageFilter<InputImageType>;
   using WeightFilterType = rtk::FDKWeightProjectionFilter<InputImageType, OutputImageType>;
   using VarianceRampFilterType = rtk::FFTVarianceRampImageFilter<OutputImageType, OutputImageType, TFFTPrecision>;
   using BackProjectionFilterType = rtk::FDKBackProjectionImageFilter<OutputImageType, OutputImageType>;

include/rtkFDKVarianceReconstructionFilter.hxx

@@ -44,7 +44,6 @@ FDKVarianceReconstructionFilter<TInputImage, TOutputImage, TFFTPrecision>::FDKVa
   m_VarianceRampFilter->SetInput(m_WeightFilter2->GetOutput());
 
   // Default parameters
-  m_ExtractFilter->SetDirectionCollapseToSubmatrix();
   m_WeightFilter1->InPlaceOn();
   m_WeightFilter2->InPlaceOn();
 
@@ -83,6 +82,9 @@ FDKVarianceReconstructionFilter<TInputImage, TOutputImage, TFFTPrecision>::Gener
   m_ExtractFilter->SetInput(this->GetInput(1));
   m_BackProjectionFilter->GetOutput()->SetRequestedRegion(this->GetOutput()->GetRequestedRegion());
   m_BackProjectionFilter->GetOutput()->PropagateRequestedRegion();
+
+  typename Superclass::InputImagePointer inputPtr1 = const_cast<TInputImage *>(this->GetInput(1));
+  inputPtr1->SetRequestedRegion(this->GetInput(1)->GetLargestPossibleRegion());
 }
 
 template <class TInputImage, class TOutputImage, class TFFTPrecision>
@@ -97,16 +99,19 @@ FDKVarianceReconstructionFilter<TInputImage, TOutputImage, TFFTPrecision>::Gener
 
   // We only set the first sub-stack at that point, the rest will be
   // requested in the GenerateData function
-  typename ExtractFilterType::InputImageRegionType projRegion;
+  typename InputImageType::RegionType projRegion;
   projRegion = this->GetInput(1)->GetLargestPossibleRegion();
   unsigned int firstStackSize = std::min(m_ProjectionSubsetSize, (unsigned int)projRegion.GetSize(Dimension - 1));
   projRegion.SetSize(Dimension - 1, firstStackSize);
+
+  m_ExtractFilter->SetInput(this->GetInput(1));
   m_ExtractFilter->SetExtractionRegion(projRegion);
+  if (rtk::IsContiguousSubRegion(this->GetInput(1), projRegion))
+    m_WeightFilter1->InPlaceOff();
 
   // Run composite filter update
   m_BackProjectionFilter->SetInput(0, this->GetInput(0));
   m_BackProjectionFilter->SetInPlace(this->GetInPlace());
-  m_ExtractFilter->SetInput(this->GetInput(1));
   m_BackProjectionFilter->UpdateOutputInformation();
 
   // Update output information
@@ -123,9 +128,10 @@ FDKVarianceReconstructionFilter<TInputImage, TOutputImage, TFFTPrecision>::Gener
   const unsigned int Dimension = this->InputImageDimension;
 
   // The backprojection works on a small stack of projections, not the full stack
-  typename ExtractFilterType::InputImageRegionType subsetRegion;
+  typename InputImageType::RegionType subsetRegion;
   subsetRegion = this->GetInput(1)->GetLargestPossibleRegion();
   unsigned int nProj = subsetRegion.GetSize(Dimension - 1);
+  unsigned int baseIndex = subsetRegion.GetIndex(Dimension - 1);
 
   // The progress accumulator tracks the progress of the pipeline
   // Each filter is equally weighted across all iterations of the stack
@@ -146,15 +152,19 @@ FDKVarianceReconstructionFilter<TInputImage, TOutputImage, TFFTPrecision>::Gener
       pimg->DisconnectPipeline();
       m_BackProjectionFilter->SetInput(pimg);
 
-      // Change projection subset
-      subsetRegion.SetIndex(Dimension - 1, i);
-      subsetRegion.SetSize(Dimension - 1, std::min(m_ProjectionSubsetSize, nProj - i));
-      m_ExtractFilter->SetExtractionRegion(subsetRegion);
-
       // This is required to reset the full pipeline
       m_BackProjectionFilter->GetOutput()->UpdateOutputInformation();
       m_BackProjectionFilter->GetOutput()->PropagateRequestedRegion();
     }
+
+    // Always create the substack for the current subset
+    subsetRegion.SetIndex(Dimension - 1, baseIndex + i);
+    subsetRegion.SetSize(Dimension - 1, std::min(m_ProjectionSubsetSize, nProj - i));
+    m_ExtractFilter->SetInput(this->GetInput(1));
+    m_ExtractFilter->SetExtractionRegion(subsetRegion);
+    if (rtk::IsContiguousSubRegion(this->GetInput(1), subsetRegion))
+      m_WeightFilter1->InPlaceOff();
+
     m_BackProjectionFilter->Update();
   }
 

include/rtkOSEMConeBeamReconstructionFilter.h

@@ -26,10 +26,10 @@
 #include <itkAddImageAdaptor.h>
 #include <itkDivideImageFilter.h>
 #include <itkDivideOrZeroOutImageFilter.h>
-#include <itkExtractImageFilter.h>
 #include <itkMultiplyImageFilter.h>
 
 #include "rtkConstantImageSource.h"
+#include "rtkSubRegionViewImageFilter.h"
 #include "rtkIterativeConeBeamReconstructionFilter.h"
 
 namespace rtk
@@ -45,7 +45,7 @@ namespace rtk
  * - DivideImageFilter,
  * - BackProjectionImageFilter.
  * The input stack of projections is processed piece by piece (the size is
- * controlled with ProjectionSubsetSize) via the use of itk::ExtractImageFilter
+ * controlled with ProjectionSubsetSize) via the use of rtk::SubRegionViewImageFilter
  * to extract sub-stacks.
  *
  * One weighting steps must be applied when processing a given subset:
@@ -69,7 +69,7 @@ namespace rtk
  *
  *  node [shape=box];
  *  ForwardProject [ label="rtk::ForwardProjectionImageFilter" URL="\ref rtk::ForwardProjectionImageFilter"];
- *  Extract [ label="itk::ExtractImageFilter" URL="\ref itk::ExtractImageFilter"];
+ *  Extract [ label="rtk::SubRegionViewImageFilter" URL="\ref rtk::SubRegionViewImageFilter"];
  *  Divide1 [ label="itk::DivideImageFilter" URL="\ref itk::DivideImageFilter"];
  *  Divide [ label="itk::DivideImageFilter" URL="\ref itk::DivideImageFilter"];
  *  ProjectionZero [ label="rtk::ConstantImageSource (full of zero)" URL="\ref rtk::ConstantImageSource"];
@@ -132,7 +132,7 @@ class ITK_TEMPLATE_EXPORT OSEMConeBeamReconstructionFilter
   using ProjectionType = TProjectionImage;
 
   /** Typedefs of each subfilter of this composite filter */
-  using ExtractFilterType = itk::ExtractImageFilter<ProjectionType, ProjectionType>;
+  using ExtractFilterType = rtk::SubRegionViewImageFilter<ProjectionType>;
   using MultiplyFilterType = itk::MultiplyImageFilter<VolumeType, VolumeType, VolumeType>;
   using ForwardProjectionFilterType = rtk::ForwardProjectionImageFilter<ProjectionType, VolumeType>;
   using BackProjectionFilterType = rtk::BackProjectionImageFilter<VolumeType, ProjectionType>;

include/rtkOSEMConeBeamReconstructionFilter.hxx

@@ -53,7 +53,6 @@ OSEMConeBeamReconstructionFilter<TVolumeImage, TProjectionImage>::OSEMConeBeamRe
   m_DivideVolumeFilter->SetInput1(m_MultiplyFilter->GetOutput());
 
   // Default parameters
-  m_ExtractFilter->SetDirectionCollapseToSubmatrix();
 }
 
 template <class TVolumeImage, class TProjectionImage>
@@ -86,7 +85,7 @@ OSEMConeBeamReconstructionFilter<TVolumeImage, TProjectionImage>::GenerateOutput
 
   // We only set the first sub-stack at that point, the rest will be
   // requested in the GenerateData function
-  typename ExtractFilterType::InputImageRegionType projRegion;
+  typename ProjectionType::RegionType projRegion;
 
   // Set forward projection filter
   m_ForwardProjectionFilter = this->InstantiateForwardProjectionFilter(this->m_CurrentForwardProjectionConfiguration);
@@ -97,22 +96,21 @@ OSEMConeBeamReconstructionFilter<TVolumeImage, TProjectionImage>::GenerateOutput
     this->InstantiateBackProjectionFilter(this->m_CurrentBackProjectionConfiguration);
 
   projRegion = this->GetInput(1)->GetLargestPossibleRegion();
-  m_ExtractFilter->SetExtractionRegion(projRegion);
-
-  m_ExtractFilter->SetInput(this->GetInput(1));
-  m_ExtractFilter->UpdateOutputInformation();
 
   // Links with the forward and back projection filters should be set here
   // and not in the constructor, as these filters are set at runtime
   m_ConstantImageSource->SetInformationFromImage(const_cast<TVolumeImage *>(this->GetInput(0)));
   m_ConstantImageSource->SetConstant(0);
 
-  m_OneConstantProjectionStackSource->SetInformationFromImage(
-    const_cast<TProjectionImage *>(m_ExtractFilter->GetOutput()));
+  m_ExtractFilter->SetInput(this->GetInput(1));
+  m_ExtractFilter->SetExtractionRegion(projRegion);
+  m_ExtractFilter->UpdateOutputInformation();
+  if (rtk::IsContiguousSubRegion(this->GetInput(1), projRegion))
+    m_DivideProjectionFilter->InPlaceOff();
+  m_OneConstantProjectionStackSource->SetInformationFromImage(m_ExtractFilter->GetOutput());
   m_OneConstantProjectionStackSource->SetConstant(1);
 
-  m_ZeroConstantProjectionStackSource->SetInformationFromImage(
-    const_cast<TProjectionImage *>(m_ExtractFilter->GetOutput()));
+  m_ZeroConstantProjectionStackSource->SetInformationFromImage(m_ExtractFilter->GetOutput());
   m_ZeroConstantProjectionStackSource->SetConstant(0);
 
   m_BackProjectionFilter->SetInput(0, m_ConstantImageSource->GetOutput());
@@ -163,9 +161,10 @@ OSEMConeBeamReconstructionFilter<TVolumeImage, TProjectionImage>::GenerateData()
   const unsigned int Dimension = this->InputImageDimension;
 
   // The backprojection works on one projection at a time
-  typename ExtractFilterType::InputImageRegionType subsetRegion;
+  typename ProjectionType::RegionType subsetRegion;
   subsetRegion = this->GetInput(1)->GetLargestPossibleRegion();
   unsigned int nProj = subsetRegion.GetSize(Dimension - 1);
+  unsigned int baseIndex = subsetRegion.GetIndex(Dimension - 1);
   subsetRegion.SetSize(Dimension - 1, 1);
 
   // Fill and shuffle randomly the projection order.
@@ -195,12 +194,25 @@ OSEMConeBeamReconstructionFilter<TVolumeImage, TProjectionImage>::GenerateData()
     for (unsigned int i = 0; i < nProj; i++)
     {
       // Change projection subset
-      subsetRegion.SetIndex(Dimension - 1, projOrder[i]);
+      subsetRegion.SetIndex(Dimension - 1, baseIndex + projOrder[i]);
+      m_ExtractFilter->SetInput(this->GetInput(1));
       m_ExtractFilter->SetExtractionRegion(subsetRegion);
       m_ExtractFilter->UpdateOutputInformation();
 
-      m_ZeroConstantProjectionStackSource->SetInformationFromImage(
-        const_cast<TProjectionImage *>(m_ExtractFilter->GetOutput()));
+      // Constant projection for normalization backprojection
+      typename ConstantProjectionSourceType::Pointer oneNormProj = ConstantProjectionSourceType::New();
+      oneNormProj->SetInformationFromImage(m_ExtractFilter->GetOutput());
+      oneNormProj->SetConstant(1);
+
+      // Constant projection for forward projection input
+      typename ConstantProjectionSourceType::Pointer zeroProj = ConstantProjectionSourceType::New();
+      zeroProj->SetInformationFromImage(m_ExtractFilter->GetOutput());
+      zeroProj->SetConstant(0);
+      zeroProj->Update();
+
+      if (rtk::IsContiguousSubRegion(this->GetInput(1), subsetRegion))
+        m_DivideProjectionFilter->InPlaceOff();
+      m_ForwardProjectionFilter->SetInput(0, zeroProj->GetOutput());
 
       // This is required to reset the full pipeline
       m_BackProjectionFilter->GetOutput()->UpdateOutputInformation();
@@ -209,8 +221,8 @@ OSEMConeBeamReconstructionFilter<TVolumeImage, TProjectionImage>::GenerateData()
       m_BackProjectionFilter->Update();
       if (iter == 0 || !m_StoreNormalizationImages)
       {
-        m_OneConstantProjectionStackSource->SetInformationFromImage(
-          const_cast<TProjectionImage *>(m_ExtractFilter->GetOutput()));
+        oneNormProj->SetInformationFromImage(m_ExtractFilter->GetOutput());
+        m_BackProjectionNormalizationFilter->SetInput(1, oneNormProj->GetOutput());
         m_BackProjectionNormalizationFilter->GetOutput()->UpdateOutputInformation();
         m_BackProjectionNormalizationFilter->GetOutput()->PropagateRequestedRegion();
         m_BackProjectionNormalizationFilter->Update();