ENH: ProgressReporter for FDK reconstruction

applications/rtkfdk/rtkfdk.cxx

@@ -32,6 +32,7 @@
 #endif
 #include "rtkFDKWarpBackProjectionImageFilter.h"
 #include "rtkCyclicDeformationImageFilter.h"
+#include "rtkProgressCommands.h"
 
 #include <itkStreamingImageFilter.h>
 #include <itkImageRegionSplitterDirection.h>
@@ -146,7 +147,11 @@ main(int argc, char * argv[])
   f->GetRampFilter()->SetTruncationCorrection(args_info.pad_arg);                                                      \
   f->GetRampFilter()->SetHannCutFrequency(args_info.hann_arg);                                                         \
   f->GetRampFilter()->SetHannCutFrequencyY(args_info.hannY_arg);                                                       \
-  f->SetProjectionSubsetSize(args_info.subsetsize_arg)
+  f->SetProjectionSubsetSize(args_info.subsetsize_arg); \
+  if(args_info.verbose_flag) \
+  { \
+    f->AddObserver(itk::ProgressEvent(), progressCommand); \
+  }
 
   // FDK reconstruction filtering
   using FDKCPUType = rtk::FDKConeBeamReconstructionFilter<OutputImageType>;
@@ -155,9 +160,14 @@ main(int argc, char * argv[])
   using FDKCUDAType = rtk::CudaFDKConeBeamReconstructionFilter;
   FDKCUDAType::Pointer feldkampCUDA;
 #endif
+
   itk::Image<OutputPixelType, Dimension> * pfeldkamp = nullptr;
   if (!strcmp(args_info.hardware_arg, "cpu"))
   {
+    // Progress reporting
+    using PercentageProgressCommandType = rtk::PercentageProgressCommand<FDKCPUType>;
+    PercentageProgressCommandType::Pointer progressCommand = PercentageProgressCommandType::New();
+
     feldkamp = FDKCPUType::New();
     SET_FELDKAMP_OPTIONS(feldkamp);
 
@@ -180,6 +190,10 @@ main(int argc, char * argv[])
       return EXIT_FAILURE;
     }
 
+    // Progress reporting
+    using PercentageProgressCommandType = rtk::PercentageProgressCommand<FDKCUDAType>;
+    PercentageProgressCommandType::Pointer progressCommand = PercentageProgressCommandType::New();
+
     feldkampCUDA = FDKCUDAType::New();
     SET_FELDKAMP_OPTIONS(feldkampCUDA);
     pfeldkamp = feldkampCUDA->GetOutput();

include/rtkCudaFFTProjectionsConvolutionImageFilter.hxx

@@ -135,6 +135,11 @@ CudaFFTProjectionsConvolutionImageFilter<TParentImageFilter>::GPUGenerateData()
                        *(float **)(cuPadImgP->GetCudaDataManager()->GetGPUBufferPointer()),
                        *(float2 **)(this->m_KernelFFTCUDA->GetCudaDataManager()->GetGPUBufferPointer()));
 
+  /* Impossible to get a pixel-wise progress reporting with CUDA
+   * so this filter only reports 0 and 100% completion. */
+  itk::ProgressReporter progress(this, 0, 1);
+  progress.CompletedPixel();
+
   // CUDA Cropping and Graft Output
   using CropFilter = CudaCropImageFilter;
   CropFilter::Pointer                            cf = CropFilter::New();

include/rtkFDKConeBeamReconstructionFilter.hxx

@@ -21,6 +21,8 @@
 
 #include "rtkFDKConeBeamReconstructionFilter.h"
 
+#include <itkProgressAccumulator.h>
+
 namespace rtk
 {
 
@@ -111,6 +113,14 @@ FDKConeBeamReconstructionFilter<TInputImage, TOutputImage, TFFTPrecision>::Gener
   subsetRegion = this->GetInput(1)->GetLargestPossibleRegion();
   unsigned int nProj = subsetRegion.GetSize(Dimension - 1);
 
+  // The progress accumulator tracks the progress of the pipeline
+  // Each filter is equally weighted across all iterations of the stack
+  itk::ProgressAccumulator::Pointer progress = itk::ProgressAccumulator::New();
+  progress->SetMiniPipelineFilter(this);
+  progress->RegisterInternalFilter(m_WeightFilter, (1.0f/3) / itk::Math::ceil(double(nProj) / m_ProjectionSubsetSize));
+  progress->RegisterInternalFilter(m_RampFilter, (1.0f/3) / itk::Math::ceil(double(nProj) / m_ProjectionSubsetSize));
+  progress->RegisterInternalFilter(m_BackProjectionFilter, (1.0f/3) / itk::Math::ceil(double(nProj) / m_ProjectionSubsetSize));
+
   for (unsigned int i = 0; i < nProj; i += m_ProjectionSubsetSize)
   {
     // After the first bp update, we need to use its output as input.

include/rtkProgressCommands.h

@@ -0,0 +1,105 @@
+/*=========================================================================
+ *
+ *  Copyright RTK Consortium
+ *
+ *  Licensed under the Apache License, Version 2.0 (the "License");
+ *  you may not use this file except in compliance with the License.
+ *  You may obtain a copy of the License at
+ *
+ *         http://www.apache.org/licenses/LICENSE-2.0.txt
+ *
+ *  Unless required by applicable law or agreed to in writing, software
+ *  distributed under the License is distributed on an "AS IS" BASIS,
+ *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ *  See the License for the specific language governing permissions and
+ *  limitations under the License.
+ *
+ *=========================================================================*/
+
+#ifndef rtkProgressCommands_h
+#define rtkIterationCommands_h
+
+namespace rtk
+{
+
+/** \class ProgressCommand
+ * \brief Abstract superclass to all progress callbacks.
+ * Derived classes must implement the Run() method.
+ *
+ * \author Aurélien Coussat
+ *
+ * \ingroup RTK
+ *
+ */
+template <typename TCaller>
+class ProgressCommand : public itk::Command
+{
+protected:
+  /** Callback function to be redefined by derived classes. */
+  virtual void
+  Run(const TCaller * caller) = 0;
+
+public:
+  /** Standard class typedefs. */
+  typedef ProgressCommand         Self;
+  typedef itk::Command            Superclass;
+  typedef itk::SmartPointer<Self> Pointer;
+
+  void
+  Execute(itk::Object * caller, const itk::EventObject & event) override
+  {
+    Execute((const itk::Object *)caller, event);
+  }
+
+  void
+  Execute(const itk::Object * caller, const itk::EventObject & event) override
+  {
+    if (!itk::ProgressEvent().CheckEvent(&event))
+    {
+      return;
+    }
+    const auto * cCaller = dynamic_cast<const TCaller *>(caller);
+    if (cCaller)
+    {
+      Run(cCaller);
+    } // TODO fail when cast fails
+  }
+};
+
+/** \class PercentageProgressCommand
+ * \brief Outputs every time a filter progresses by at least one percent.
+ *
+ * \author Aurélien Coussat
+ *
+ * \ingroup RTK
+ *
+ */
+template <typename TCaller>
+class PercentageProgressCommand : public ProgressCommand<TCaller>
+{
+public:
+  /** Standard class typedefs. */
+  typedef PercentageProgressCommand Self;
+  typedef ProgressCommand<TCaller>  Superclass;
+  typedef itk::SmartPointer<Self>   Pointer;
+  itkNewMacro(Self);
+
+  int percentage = -1;
+
+protected:
+  void
+  Run(const TCaller * caller) override
+  {
+    int newPercentage = (int)(caller->GetProgress() * 100.);
+    if (newPercentage > percentage)
+    {
+      std::cout << caller->GetNameOfClass() << " " << newPercentage << "% completed."
+                << std::endl; // TODO allow string personnalization
+      percentage = newPercentage;
+    }
+  }
+};
+
+} // end namespace rtk
+
+#endif

src/rtkCudaFDKBackProjectionImageFilter.cxx

@@ -24,6 +24,7 @@
 #include <itkImageRegionIteratorWithIndex.h>
 #include <itkLinearInterpolateImageFunction.h>
 #include <itkMacro.h>
+#include <itkProgressReporter.h>
 
 namespace rtk
 {
@@ -103,6 +104,9 @@ CudaFDKBackProjectionImageFilter ::GPUGenerateData()
       fMatrix[j + (iProj - iFirstProj) * 12] = matrix[j / 4][j % 4];
   }
 
+  // Slab-wise progress reporting
+  itk::ProgressReporter progress(this, 0, itk::Math::ceil(double(nProj) / SLAB_SIZE));
+
   for (unsigned int i = 0; i < nProj; i += SLAB_SIZE)
   {
     // If nProj is not a multiple of SLAB_SIZE, the last slab will contain less than SLAB_SIZE projections
@@ -113,6 +117,8 @@ CudaFDKBackProjectionImageFilter ::GPUGenerateData()
 
     // Re-use the output as input
     pin = pout;
+
+    progress.CompletedPixel();
   }
 
   delete[] fMatrix;

src/rtkCudaFDKWeightProjectionFilter.cxx

@@ -122,6 +122,11 @@ CudaFDKWeightProjectionFilter ::GPUGenerateData()
                          proj_row,
                          proj_col);
 
+  /* Impossible to get a pixel-wise progress reporting with CUDA
+   * so this filter only reports 0 and 100% completion. */
+  itk::ProgressReporter progress(this, 0, 1);
+  progress.CompletedPixel();
+
   delete[] geomMatrix;
 }