Displaced detector

include/rtkDisplacedDetectorForOffsetFieldOfViewImageFilter.hxx

@@ -67,40 +67,89 @@ DisplacedDetectorForOffsetFieldOfViewImageFilter<TInputImage, TOutputImage>::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<OutputImageType, OutputImageType>;
-  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<OutputImageType, OutputImageType>;
+    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);

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;

include/rtkDisplacedDetectorImageFilter.hxx

@@ -112,6 +112,12 @@ DisplacedDetectorImageFilter<TInputImage, TOutputImage>::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;

test/rtkdisplaceddetectorcompoffsettest.cxx

@@ -96,6 +96,30 @@ main(int, char **)
 
     CheckImageQuality<OutputImageType>(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<OutputImageType>;
+    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<OutputImageType>;
+    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<OutputImageType>(cudaddf->GetOutput(), cpuddf->GetOutput(), 1.e-6, 100, 1.);
+  }
 
   // If all succeed
   std::cout << "\n\nTest PASSED! " << std::endl;