diff --git a/include/rtkFDKConeBeamReconstructionFilter.h b/include/rtkFDKConeBeamReconstructionFilter.h index 85bc1104f..00eea83bb 100644 --- a/include/rtkFDKConeBeamReconstructionFilter.h +++ b/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 diff --git a/include/rtkFDKConeBeamReconstructionFilter.hxx b/include/rtkFDKConeBeamReconstructionFilter.hxx index aac8d6d25..05bb61cc3 100644 --- a/include/rtkFDKConeBeamReconstructionFilter.hxx +++ b/include/rtkFDKConeBeamReconstructionFilter.hxx @@ -22,6 +22,8 @@ #include +#include + namespace rtk { @@ -63,6 +65,23 @@ FDKConeBeamReconstructionFilter::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 diff --git a/test/rtkfdktest.cxx b/test/rtkfdktest.cxx index b1d8197fd..5906bee2a 100644 --- a/test/rtkfdktest.cxx +++ b/test/rtkfdktest.cxx @@ -1,6 +1,7 @@ #include #include #include +#include #include "rtkConstantImageSource.h" #include "rtkDrawSheppLoganFilter.h" @@ -180,5 +181,21 @@ rtkfdktest(int, char *[]) TRY_AND_EXIT_ON_ITK_EXCEPTION(dsl->UpdateLargestPossibleRegion()) CheckImageQuality(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; }