ENH: Add superior and inferior clip images for Joseph-based forward p…

…rojector

The superior and inferior clip images define for each pixel of the
projections the value of the superior and inferior clip of the ray
emitted from that pixel. Each ray is clipped from
source+inferior_clip*(pixel-source) to
source+superior_clip*(pixel-source) with inferior_clip and superior_clip
equal 0 and 1 by default.

applications/rtkforwardprojections/rtkforwardprojections.cxx

@@ -80,6 +80,23 @@ main(int argc, char * argv[])
     attenuationMap = itk::ReadImage<OutputImageType>(args_info.attenuationmap_arg);
   }
 
+  using ClipImageType = itk::Image<double, Dimension>;
+  ClipImageType::Pointer inferiorClipImage, superiorClipImage;
+  if (args_info.inferiorclipimage_given)
+  {
+    if (args_info.verbose_flag)
+      std::cout << "Reading inferior clip image " << args_info.inferiorclipimage_arg << "..." << std::endl;
+    // Read an existing image to initialize the attenuation map
+    inferiorClipImage = itk::ReadImage<ClipImageType>(args_info.inferiorclipimage_arg);
+  }
+  if (args_info.superiorclipimage_given)
+  {
+    if (args_info.verbose_flag)
+      std::cout << "Reading superior clip image " << args_info.superiorclipimage_arg << "..." << std::endl;
+    // Read an existing image to initialize the attenuation map
+    superiorClipImage = itk::ReadImage<ClipImageType>(args_info.superiorclipimage_arg);
+  }
+
   // Create forward projection image filter
   if (args_info.verbose_flag)
     std::cout << "Projecting volume..." << std::endl;
@@ -116,6 +133,36 @@ main(int argc, char * argv[])
   forwardProjection->SetInput(1, inputVolume);
   if (args_info.attenuationmap_given)
     forwardProjection->SetInput(2, attenuationMap);
+  if (args_info.inferiorclipimage_given)
+  {
+    if (args_info.fp_arg == fp_arg_Joseph)
+    {
+      dynamic_cast<rtk::JosephForwardProjectionImageFilter<OutputImageType, OutputImageType> *>(
+        forwardProjection.GetPointer())
+        ->SetInferiorClipImage(inferiorClipImage);
+    }
+    else if (args_info.fp_arg == fp_arg_JosephAttenuated)
+    {
+      dynamic_cast<rtk::JosephForwardAttenuatedProjectionImageFilter<OutputImageType, OutputImageType> *>(
+        forwardProjection.GetPointer())
+        ->SetInferiorClipImage(inferiorClipImage);
+    }
+  }
+  if (args_info.superiorclipimage_given)
+  {
+    if (args_info.fp_arg == fp_arg_Joseph)
+    {
+      dynamic_cast<rtk::JosephForwardProjectionImageFilter<OutputImageType, OutputImageType> *>(
+        forwardProjection.GetPointer())
+        ->SetSuperiorClipImage(superiorClipImage);
+    }
+    else if (args_info.fp_arg == fp_arg_JosephAttenuated)
+    {
+      dynamic_cast<rtk::JosephForwardAttenuatedProjectionImageFilter<OutputImageType, OutputImageType> *>(
+        forwardProjection.GetPointer())
+        ->SetSuperiorClipImage(superiorClipImage);
+    }
+  }
   if (args_info.sigmazero_given && args_info.fp_arg == fp_arg_Zeng)
     dynamic_cast<rtk::ZengForwardProjectionImageFilter<OutputImageType, OutputImageType> *>(
       forwardProjection.GetPointer())

applications/rtkforwardprojections/rtkforwardprojections.ggo

@@ -14,3 +14,5 @@ option "fp"    f "Forward projection method" values="Joseph","JosephAttenuated",
 option "attenuationmap" - "Attenuation map relative to the volume to perfom the attenuation correction"   string  no
 option "sigmazero" - "PSF value at a distance of 0 meter of the detector"   double  no
 option "alphapsf" - "Slope of the PSF against the detector distance"   double  no 
+option "inferiorclipimage" - "Value of the inferior clip of the ray for each pixel of the projections (only with Joseph-based projector)" string no
+option "superiorclipimage" - "Value of the superior clip of the ray for each pixel of the projections (only with Joseph-based projector)" string no

include/rtkJosephForwardAttenuatedProjectionImageFilter.hxx

@@ -81,6 +81,7 @@ JosephForwardAttenuatedProjectionImageFilter<TInputImage,
                                              TProjectedValueAccumulation,
                                              TComputeAttenuationCorrection>::VerifyInputInformation() const
 {
+  Superclass::VerifyInputInformation();
   using ImageBaseType = const itk::ImageBase<InputImageDimension>;
 
   ImageBaseType * inputPtr1 = nullptr;
@@ -93,7 +94,7 @@ JosephForwardAttenuatedProjectionImageFilter<TInputImage,
     // method since it returns the input as a pointer to a
     // DataObject as opposed to the subclass version which
     // static_casts the input to an TInputImage).
-    if (it.GetName() != "Primary")
+    if (it.GetName() == "_1")
     {
       inputPtr1 = dynamic_cast<ImageBaseType *>(it.GetInput());
     }
@@ -105,7 +106,7 @@ JosephForwardAttenuatedProjectionImageFilter<TInputImage,
 
   for (; !it.IsAtEnd(); ++it)
   {
-    if (it.GetName() != "Primary")
+    if (it.GetName() == "_2")
     {
       auto * inputPtrN = dynamic_cast<ImageBaseType *>(it.GetInput());
       // Physical space computation only matters if we're using two

include/rtkJosephForwardProjectionImageFilter.h

@@ -183,6 +183,8 @@ class ITK_TEMPLATE_EXPORT JosephForwardProjectionImageFilter
   using OutputImageRegionType = typename TOutputImage::RegionType;
   using CoordRepType = double;
   using VectorType = itk::Vector<CoordRepType, TInputImage::ImageDimension>;
+  using TClipImageType = itk::Image<double, TOutputImage::ImageDimension>;
+  using TClipImagePointer = typename TClipImageType::Pointer;
 
   /** Method for creation through the object factory. */
   itkNewMacro(Self);
@@ -254,6 +256,36 @@ class ITK_TEMPLATE_EXPORT JosephForwardProjectionImageFilter
     }
   }
 
+  /** Set/Get the inferior clip image. Each pixel of the image
+   ** corresponds to the value of the inferior clip of the ray
+   ** emitted from that pixel. */
+  void
+  SetInferiorClipImage(const TClipImageType * inferiorClipImage)
+  {
+    // Process object is not const-correct so the const casting is required.
+    this->SetInput("InferiorClipImage", const_cast<TClipImageType *>(inferiorClipImage));
+  }
+  typename TClipImageType::ConstPointer
+  GetInferiorClipImage()
+  {
+    return static_cast<const TClipImageType *>(this->itk::ProcessObject::GetInput("InferiorClipImage"));
+  }
+
+  /** Set/Get the superior clip image. Each pixel of the image
+   ** corresponds to the value of the superior clip of the ray
+   ** emitted from that pixel. */
+  void
+  SetSuperiorClipImage(const TClipImageType * superiorClipImage)
+  {
+    // Process object is not const-correct so the const casting is required.
+    this->SetInput("SuperiorClipImage", const_cast<TClipImageType *>(superiorClipImage));
+  }
+  typename TClipImageType::ConstPointer
+  GetSuperiorClipImage()
+  {
+    return static_cast<const TClipImageType *>(this->itk::ProcessObject::GetInput("SuperiorClipImage"));
+  }
+
   /** Each ray is clipped from source+m_InferiorClip*(pixel-source) to
   ** source+m_SuperiorClip*(pixel-source) with m_InferiorClip and
   ** m_SuperiorClip equal 0 and 1 by default. */
@@ -266,14 +298,17 @@ class ITK_TEMPLATE_EXPORT JosephForwardProjectionImageFilter
   JosephForwardProjectionImageFilter();
   ~JosephForwardProjectionImageFilter() override = default;
 
+  /** Apply changes to the input image requested region. */
+  void
+  GenerateInputRequestedRegion() override;
+
   void
   ThreadedGenerateData(const OutputImageRegionType & outputRegionForThread, ThreadIdType threadId) override;
 
-  /** The two inputs should not be in the same space so there is nothing
-   * to verify. */
+  /** If a third input is given, it should be in the same physical space
+   * than the first one. */
   void
-  VerifyInputInformation() const override
-  {}
+  VerifyInputInformation() const override;
 
   inline OutputPixelType
   BilinearInterpolation(const ThreadIdType     threadId,

include/rtkJosephForwardProjectionImageFilter.hxx

@@ -45,6 +45,130 @@ JosephForwardProjectionImageFilter<TInputImage,
   this->DynamicMultiThreadingOff();
 }
 
+
+template <class TInputImage,
+          class TOutputImage,
+          class TInterpolationWeightMultiplication,
+          class TProjectedValueAccumulation,
+          class TSumAlongRay>
+void
+JosephForwardProjectionImageFilter<TInputImage,
+                                   TOutputImage,
+                                   TInterpolationWeightMultiplication,
+                                   TProjectedValueAccumulation,
+                                   TSumAlongRay>::GenerateInputRequestedRegion()
+{
+  Superclass::GenerateInputRequestedRegion();
+
+  if (this->GetInferiorClipImage().IsNotNull())
+  {
+    TClipImagePointer inputInferiorClipImage = const_cast<TClipImageType *>(this->GetInferiorClipImage().GetPointer());
+    if (!inputInferiorClipImage)
+      return;
+    inputInferiorClipImage->SetRequestedRegion(inputInferiorClipImage->GetRequestedRegion());
+  }
+
+  if (this->GetSuperiorClipImage().IsNotNull())
+  {
+    TClipImagePointer inputSuperiorClipImage = const_cast<TClipImageType *>(this->GetSuperiorClipImage().GetPointer());
+    if (!inputSuperiorClipImage)
+      return;
+    inputSuperiorClipImage->SetRequestedRegion(inputSuperiorClipImage->GetRequestedRegion());
+  }
+}
+
+
+template <class TInputImage,
+          class TOutputImage,
+          class TInterpolationWeightMultiplication,
+          class TProjectedValueAccumulation,
+          class TSumAlongRay>
+void
+JosephForwardProjectionImageFilter<TInputImage,
+                                   TOutputImage,
+                                   TInterpolationWeightMultiplication,
+                                   TProjectedValueAccumulation,
+                                   TSumAlongRay>::VerifyInputInformation() const
+{
+  using ImageBaseType = const itk::ImageBase<TInputImage::ImageDimension>;
+
+  ImageBaseType * inputPtr1 = nullptr;
+
+  itk::InputDataObjectConstIterator it(this);
+  for (; !it.IsAtEnd(); ++it)
+  {
+    // Check whether the output is an image of the appropriate
+    // dimension (use ProcessObject's version of the GetInput()
+    // method since it returns the input as a pointer to a
+    // DataObject as opposed to the subclass version which
+    // static_casts the input to an TInputImage).
+    if (it.GetName() == "Primary")
+    {
+      inputPtr1 = dynamic_cast<ImageBaseType *>(it.GetInput());
+    }
+    if (inputPtr1)
+    {
+      break;
+    }
+  }
+
+  for (; !it.IsAtEnd(); ++it)
+  {
+    if (it.GetName() == "InferiorClipImage" || it.GetName() == "SuperiorClipImage")
+    {
+      auto * inputPtrN = dynamic_cast<ImageBaseType *>(it.GetInput());
+      // Physical space computation only matters if we're using two
+      // images, and not an image and a constant.
+      if (inputPtrN)
+      {
+        // check that the image occupy the same physical space, and that
+        // each index is at the same physical location
+
+        // tolerance for origin and spacing depends on the size of pixel
+        // tolerance for directions a fraction of the unit cube.
+        const double coordinateTol = itk::Math::abs(Self::GetGlobalDefaultCoordinateTolerance() *
+                                                    inputPtr1->GetSpacing()[0]); // use first dimension spacing
+
+        if (!inputPtr1->GetOrigin().GetVnlVector().is_equal(inputPtrN->GetOrigin().GetVnlVector(), coordinateTol) ||
+            !inputPtr1->GetSpacing().GetVnlVector().is_equal(inputPtrN->GetSpacing().GetVnlVector(), coordinateTol) ||
+            !inputPtr1->GetDirection().GetVnlMatrix().as_ref().is_equal(
+              inputPtrN->GetDirection().GetVnlMatrix().as_ref(), Self::GetGlobalDefaultDirectionTolerance()))
+        {
+          std::ostringstream originString, spacingString, directionString;
+          if (!inputPtr1->GetOrigin().GetVnlVector().is_equal(inputPtrN->GetOrigin().GetVnlVector(), coordinateTol))
+          {
+            originString.setf(std::ios::scientific);
+            originString.precision(7);
+            originString << "InputImage Origin: " << inputPtr1->GetOrigin() << ", InputImage" << it.GetName()
+                         << " Origin: " << inputPtrN->GetOrigin() << std::endl;
+            originString << "\tTolerance: " << coordinateTol << std::endl;
+          }
+          if (!inputPtr1->GetSpacing().GetVnlVector().is_equal(inputPtrN->GetSpacing().GetVnlVector(), coordinateTol))
+          {
+            spacingString.setf(std::ios::scientific);
+            spacingString.precision(7);
+            spacingString << "InputImage Spacing: " << inputPtr1->GetSpacing() << ", InputImage" << it.GetName()
+                          << " Spacing: " << inputPtrN->GetSpacing() << std::endl;
+            spacingString << "\tTolerance: " << coordinateTol << std::endl;
+          }
+          if (!inputPtr1->GetDirection().GetVnlMatrix().as_ref().is_equal(
+                inputPtrN->GetDirection().GetVnlMatrix().as_ref(), Self::GetGlobalDefaultDirectionTolerance()))
+          {
+            directionString.setf(std::ios::scientific);
+            directionString.precision(7);
+            directionString << "InputImage Direction: " << inputPtr1->GetDirection() << ", InputImage" << it.GetName()
+                            << " Direction: " << inputPtrN->GetDirection() << std::endl;
+            directionString << "\tTolerance: " << Self::GetGlobalDefaultDirectionTolerance() << std::endl;
+          }
+          itkExceptionMacro(<< "Inputs do not occupy the same physical space! " << std::endl
+                            << originString.str() << spacingString.str() << directionString.str());
+        }
+      }
+    }
+  }
+}
+
+
 template <class TInputImage,
           class TOutputImage,
           class TInterpolationWeightMultiplication,
@@ -85,6 +209,17 @@ JosephForwardProjectionImageFilter<TInputImage,
   itIn = InputRegionIterator::New(this->GetInput(), outputRegionForThread, geometry, volPPToIndex);
   using OutputRegionIterator = itk::ImageRegionIteratorWithIndex<TOutputImage>;
   OutputRegionIterator itOut(this->GetOutput(), outputRegionForThread);
+  using ClipImageIterator = itk::ImageRegionConstIterator<TClipImageType>;
+  ClipImageIterator * itInferiorImage = nullptr;
+  ClipImageIterator * itSuperiorImage = nullptr;
+  if (this->GetInferiorClipImage().IsNotNull())
+  {
+    itInferiorImage = new ClipImageIterator(this->GetInferiorClipImage(), outputRegionForThread);
+  }
+  if (this->GetSuperiorClipImage().IsNotNull())
+  {
+    itSuperiorImage = new ClipImageIterator(this->GetSuperiorClipImage(), outputRegionForThread);
+  }
 
   // Create intersection functions, one for each possible main direction
   typename BoxShape::Pointer    box = BoxShape::New();
@@ -108,8 +243,9 @@ JosephForwardProjectionImageFilter<TInputImage,
   typename BoxShape::VectorType stepMM, np, fp;
   for (unsigned int pix = 0; pix < outputRegionForThread.GetNumberOfPixels(); pix++, itIn->Next(), ++itOut)
   {
-    typename InputRegionIterator::PointType pixelPosition = itIn->GetPixelPosition();
-    typename InputRegionIterator::PointType dirVox = -itIn->GetSourceToPixel();
+    typename InputRegionIterator::PointType  pixelPosition = itIn->GetPixelPosition();
+    typename InputRegionIterator::PointType  dirVox = -itIn->GetSourceToPixel();
+    typename InputRegionIterator ::IndexType pixelIndex = itIn->GetIndex();
 
     // Select main direction
     unsigned int                  mainDir = 0;
@@ -121,16 +257,30 @@ JosephForwardProjectionImageFilter<TInputImage,
         mainDir = i;
     }
 
+    if (this->GetInferiorClipImage().IsNotNull())
+    {
+      double superiorClipImage = 1 - itInferiorImage->Get();
+      superiorClip = std::min(1. - m_InferiorClip, superiorClipImage);
+      ++(*itInferiorImage);
+    }
+    if (this->GetSuperiorClipImage().IsNotNull())
+    {
+      double inferiorClipImage = 1 - itSuperiorImage->Get();
+      inferiorClip = std::max(1. - m_SuperiorClip, inferiorClipImage);
+      ++(*itSuperiorImage);
+    }
+
     // Test if there is an intersection
     BoxShape::ScalarType nearDist = NAN, farDist = NAN;
-    if (box->IsIntersectedByRay(pixelPosition, dirVox, nearDist, farDist) &&
-        farDist >= 0. && // check if detector after the source
-        nearDist <= 1.)  // check if detector after or in the volume
+    bool                 isIntersectByRay = box->IsIntersectedByRay(pixelPosition, dirVox, nearDist, farDist);
+    // Clip the casting between source and pixel of the detector
+    nearDist = std::max(nearDist, inferiorClip);
+    farDist = std::min(farDist, superiorClip);
+
+    if (isIntersectByRay && farDist > 0. && // check if detector after the source
+        nearDist <= 1. &&                   // check if detector after or in the volume
+        farDist > nearDist)
     {
-      // Clip the casting between source and pixel of the detector
-      nearDist = std::max(nearDist, inferiorClip);
-      farDist = std::min(farDist, superiorClip);
-
       // Compute and sort intersections: (n)earest and (f)arthest (p)points
       np = pixelPosition + nearDist * dirVox;
       fp = pixelPosition + farDist * dirVox;