diff --git a/include/rtkDisplacedDetectorForOffsetFieldOfViewImageFilter.hxx b/include/rtkDisplacedDetectorForOffsetFieldOfViewImageFilter.hxx index 97a9dd4fd..438098da3 100644 --- a/include/rtkDisplacedDetectorForOffsetFieldOfViewImageFilter.hxx +++ b/include/rtkDisplacedDetectorForOffsetFieldOfViewImageFilter.hxx @@ -67,40 +67,89 @@ DisplacedDetectorForOffsetFieldOfViewImageFilter::Gen << "Consider disabling it by setting m_Disable=true " << "or using the nodisplaced flag of the application you are running"); } - - using FOVFilterType = typename rtk::FieldOfViewImageFilter; - auto fieldofview = FOVFilterType::New(); - fieldofview->SetProjectionsStack(inputPtr.GetPointer()); - fieldofview->SetGeometry(this->GetGeometry()); - bool hasOverlap = fieldofview->ComputeFOVRadius(FOVFilterType::RADIUSBOTH, m_FOVCenterX, m_FOVCenterZ, m_FOVRadius); - double xi = NAN, zi = NAN, ri = NAN; - fieldofview->ComputeFOVRadius(FOVFilterType::RADIUSINF, xi, zi, ri); - double xs = NAN, zs = NAN, rs = NAN; - fieldofview->ComputeFOVRadius(FOVFilterType::RADIUSSUP, xs, zs, rs); - - // 4 cases depending on the position of the two corners - // Case 1: Impossible to account for too large displacements - if (!hasOverlap) - { - itkGenericExceptionMacro(<< "Cannot account for too large detector displacements, a part of" - << " space must be covered by all projections."); - } - // Case 2: Not displaced if less than 10% relative difference between radii - else if (200. * itk::Math::abs(ri - rs) / (ri + rs) < 10.) + else if (this->GetGeometry()->GetSourceToDetectorDistances().front() == 0.) { + itkGenericExceptionMacro(<< "Displaced detector cannot handle parallel geometry. " + << "Consider disabling it by setting m_Disable=true " + << "or using the nodisplaced flag of the application you are running"); } - else if (rs > ri) + + if (!this->GetOffsetsSet()) { - this->SetInPlace(false); - itk::Index<3>::IndexValueType index = - outputLargestPossibleRegion.GetIndex()[0] - outputLargestPossibleRegion.GetSize()[0]; - outputLargestPossibleRegion.SetIndex(0, index); - outputLargestPossibleRegion.SetSize(0, outputLargestPossibleRegion.GetSize()[0] * 2); + using FOVFilterType = typename rtk::FieldOfViewImageFilter; + auto fieldofview = FOVFilterType::New(); + fieldofview->SetProjectionsStack(inputPtr.GetPointer()); + fieldofview->SetGeometry(this->GetGeometry()); + bool hasOverlap = fieldofview->ComputeFOVRadius(FOVFilterType::RADIUSBOTH, m_FOVCenterX, m_FOVCenterZ, m_FOVRadius); + double xi = NAN, zi = NAN, ri = NAN; + fieldofview->ComputeFOVRadius(FOVFilterType::RADIUSINF, xi, zi, ri); + double xs = NAN, zs = NAN, rs = NAN; + fieldofview->ComputeFOVRadius(FOVFilterType::RADIUSSUP, xs, zs, rs); + + // 4 cases depending on the position of the two corners + // Case 1: Impossible to account for too large displacements + if (!hasOverlap) + { + itkGenericExceptionMacro(<< "Cannot account for too large detector displacements, a part of" + << " space must be covered by all projections."); + } + // Case 2: Not displaced if less than 10% relative difference between radii + else if (200. * itk::Math::abs(ri - rs) / (ri + rs) < 10.) + { + } + else if (rs > ri) + { + this->SetInPlace(false); + itk::Index<3>::IndexValueType index = + outputLargestPossibleRegion.GetIndex()[0] - outputLargestPossibleRegion.GetSize()[0]; + outputLargestPossibleRegion.SetIndex(0, index); + outputLargestPossibleRegion.SetSize(0, outputLargestPossibleRegion.GetSize()[0] * 2); + } + else + { + this->SetInPlace(false); + outputLargestPossibleRegion.SetSize(0, outputLargestPossibleRegion.GetSize()[0] * 2); + } } else { - this->SetInPlace(false); - outputLargestPossibleRegion.SetSize(0, outputLargestPossibleRegion.GetSize()[0] * 2); + // If the user manually set offsets, apply the same corner-based logic as the base class + typename Superclass::InputImageType::PointType corner; + inputPtr->TransformIndexToPhysicalPoint(inputPtr->GetLargestPossibleRegion().GetIndex(), corner); + double inferiorCorner = corner[0]; + double superiorCorner = inferiorCorner; + if (inputPtr->GetSpacing()[0] < 0.) + inferiorCorner += inputPtr->GetSpacing()[0] * (outputLargestPossibleRegion.GetSize(0) - 1); + else + superiorCorner += inputPtr->GetSpacing()[0] * (outputLargestPossibleRegion.GetSize(0) - 1); + + inferiorCorner += this->GetMaximumOffset(); + superiorCorner += this->GetMinimumOffset(); + + if (inferiorCorner > 0. || superiorCorner < 0.) + { + itkGenericExceptionMacro(<< "Cannot account for detector displacement larger than 50% of panel size." + << " Corner inf=" << inferiorCorner << " and corner sup=" << superiorCorner); + } + else if ((itk::Math::abs(inferiorCorner + superiorCorner) < + 0.1 * itk::Math::abs(superiorCorner - inferiorCorner)) || + !this->m_PadOnTruncatedSide) + { + // nothing to do, keep region + } + else if (superiorCorner + inferiorCorner > 0.) + { + this->SetInPlace(false); + typename itk::Index<3>::IndexValueType index = + outputLargestPossibleRegion.GetIndex()[0] - outputLargestPossibleRegion.GetSize()[0]; + outputLargestPossibleRegion.SetIndex(0, index); + outputLargestPossibleRegion.SetSize(0, outputLargestPossibleRegion.GetSize()[0] * 2); + } + else + { + this->SetInPlace(false); + outputLargestPossibleRegion.SetSize(0, outputLargestPossibleRegion.GetSize()[0] * 2); + } } outputPtr->SetLargestPossibleRegion(outputLargestPossibleRegion); diff --git a/include/rtkDisplacedDetectorImageFilter.h b/include/rtkDisplacedDetectorImageFilter.h index d331b5dd7..c9b39a96c 100644 --- a/include/rtkDisplacedDetectorImageFilter.h +++ b/include/rtkDisplacedDetectorImageFilter.h @@ -116,6 +116,13 @@ class ITK_TEMPLATE_EXPORT DisplacedDetectorImageFilter : public itk::InPlaceImag itkGetMacro(InferiorCorner, double); itkGetMacro(SuperiorCorner, double); + /** Returns true when offsets were explicitly provided */ + bool + GetOffsetsSet() const + { + return m_OffsetsSet; + } + /** Checks that inputs are correctly set. */ void VerifyPreconditions() const override; diff --git a/include/rtkDisplacedDetectorImageFilter.hxx b/include/rtkDisplacedDetectorImageFilter.hxx index e9444343d..74094628f 100644 --- a/include/rtkDisplacedDetectorImageFilter.hxx +++ b/include/rtkDisplacedDetectorImageFilter.hxx @@ -112,6 +112,12 @@ DisplacedDetectorImageFilter::GenerateOutputInformati << "Consider disabling it by setting m_Disable=true " << "or using the nodisplaced flag of the application you are running"); } + else if (this->GetGeometry()->GetSourceToDetectorDistances().front() == 0.) + { + itkGenericExceptionMacro(<< "Displaced detector cannot handle parallel geometry. " + << "Consider disabling it by setting m_Disable=true " + << "or using the nodisplaced flag of the application you are running"); + } // Compute the X coordinates of the corners of the image typename Superclass::InputImageType::PointType corner; diff --git a/test/rtkdisplaceddetectorcompoffsettest.cxx b/test/rtkdisplaceddetectorcompoffsettest.cxx index 32303f668..cb41cd43e 100644 --- a/test/rtkdisplaceddetectorcompoffsettest.cxx +++ b/test/rtkdisplaceddetectorcompoffsettest.cxx @@ -96,6 +96,30 @@ main(int, char **) CheckImageQuality(streamingCUDA->GetOutput(), streamingCPU->GetOutput(), 1.e-6, 100, 1.); } + { + constexpr double minOffset = -16.; + constexpr double maxOffset = 12.; + + std::cout << "\n\n****** Case 4: manual offsets ******" << std::endl; + + using OffsetDDFType = rtk::DisplacedDetectorForOffsetFieldOfViewImageFilter; + auto cudaddf = OffsetDDFType::New(); + cudaddf->SetInput(slp->GetOutput()); + cudaddf->SetGeometry(geometry); + cudaddf->SetOffsets(minOffset, maxOffset); + cudaddf->InPlaceOff(); + TRY_AND_EXIT_ON_ITK_EXCEPTION(cudaddf->Update()); + + using CPUDDFType = rtk::DisplacedDetectorImageFilter; + auto cpuddf = CPUDDFType::New(); + cpuddf->SetInput(slp->GetOutput()); + cpuddf->SetGeometry(geometry); + cpuddf->SetOffsets(minOffset, maxOffset); + cpuddf->InPlaceOff(); + TRY_AND_EXIT_ON_ITK_EXCEPTION(cpuddf->Update()); + + CheckImageQuality(cudaddf->GetOutput(), cpuddf->GetOutput(), 1.e-6, 100, 1.); + } // If all succeed std::cout << "\n\nTest PASSED! " << std::endl;