ENH: Check rotation axis is along Y in FDK reconstruction

include/rtkFDKConeBeamReconstructionFilter.h

@@ -111,6 +111,12 @@ class ITK_TEMPLATE_EXPORT FDKConeBeamReconstructionFilter : public itk::InPlaceI
   itkGetMacro(ProjectionSubsetSize, unsigned int);
   itkSetMacro(ProjectionSubsetSize, unsigned int);
 
+  /** Get / Set the maximum angular deviation, in degrees, allowed between the rotation
+   * axis and the Y axis. FDK reconstruction assumes a rotation around Y and will throw
+   * if the deviation exceeds this tolerance. Default is 30 degrees. */
+  itkGetConstMacro(AngularTolerance, double);
+  itkSetMacro(AngularTolerance, double);
+
   /** Get / Set and init the backprojection filter. The set function takes care
    * of initializing the mini-pipeline and the ramp filter must therefore be
    * created before calling this set function. */
@@ -151,6 +157,9 @@ class ITK_TEMPLATE_EXPORT FDKConeBeamReconstructionFilter : public itk::InPlaceI
   /** Number of projections processed at a time. */
   unsigned int m_ProjectionSubsetSize{ 16 };
 
+  /** Maximum angular deviation, in degrees, between the rotation axis and the Y axis. */
+  double m_AngularTolerance{ 30. };
+
   /** Geometry propagated to subfilters of the mini-pipeline. */
   ThreeDCircularProjectionGeometry::Pointer m_Geometry;
 }; // end of class

include/rtkFDKConeBeamReconstructionFilter.hxx

@@ -22,6 +22,8 @@
 
 #include <itkProgressAccumulator.h>
 
+#include <cmath>
+
 namespace rtk
 {
 
@@ -63,6 +65,23 @@ FDKConeBeamReconstructionFilter<TInputImage, TOutputImage, TFFTPrecision>::Verif
 
   if (this->m_Geometry.IsNull())
     itkExceptionMacro(<< "Geometry has not been set.");
+
+  // Check that the rotation axis is along Y (index 1), as required by the ramp filter
+  // which always applies along the first direction of the projection images. Some
+  // misalignment is allowed, up to the configurable angular tolerance.
+  const double cosTolerance = std::cos(this->m_AngularTolerance * itk::Math::pi / 180.);
+  const auto & matrices = this->m_Geometry->GetRotationMatrices();
+  for (unsigned int i = 0; i < matrices.size(); i++)
+  {
+    const double dot = matrices[i][1][1];
+    if (dot < cosTolerance)
+    {
+      itkExceptionMacro(<< "The rotation axis is not along the Y direction for projection #" << i
+                        << ": it is tilted by " << std::acos(dot) * 180. / itk::Math::pi
+                        << " degrees while the angular tolerance is " << this->m_AngularTolerance
+                        << " degrees. FDK reconstruction assumes a rotation around Y.");
+    }
+  }
 }
 
 template <class TInputImage, class TOutputImage, class TFFTPrecision>

test/rtkfdktest.cxx

@@ -1,6 +1,7 @@
 #include <itkImageRegionConstIterator.h>
 #include <itkImageRegionSplitterDirection.h>
 #include <itkStreamingImageFilter.h>
+#include <itkTestingMacros.h>
 
 #include "rtkConstantImageSource.h"
 #include "rtkDrawSheppLoganFilter.h"
@@ -180,5 +181,21 @@ rtkfdktest(int, char *[])
   TRY_AND_EXIT_ON_ITK_EXCEPTION(dsl->UpdateLargestPossibleRegion())
   CheckImageQuality<OutputImageType>(fov->GetOutput(), dsl->GetOutput(), 0.03, 26, 2.0);
   std::cout << "Test PASSED! " << std::endl;
+
+  std::cout << "\n\n****** Case 6: rotation axis tilted beyond tolerance should throw ******" << std::endl;
+
+  // Geometry with an out-of-plane angle larger than the default tolerance (30 degrees)
+  auto badGeometry = rtk::ThreeDCircularProjectionGeometry::New();
+  for (unsigned int noProj = 0; noProj < NumberOfProjectionImages; noProj++)
+    badGeometry->AddProjection(600., 1200., noProj * 360. / NumberOfProjectionImages, 0, 0, 60);
+
+  auto badFeldkamp = FDKType::New();
+  badFeldkamp->SetInput(0, tomographySource->GetOutput());
+  badFeldkamp->SetInput(1, slp->GetOutput());
+  badFeldkamp->SetGeometry(badGeometry);
+
+  ITK_TRY_EXPECT_EXCEPTION(badFeldkamp->Update());
+  std::cout << "Test PASSED! " << std::endl;
+
   return EXIT_SUCCESS;
 }