ENH: add projections weights in one step reconstruction

applications/rtkspectralonestep/rtkspectralonestep.cxx

@@ -112,6 +112,14 @@ rtkspectralonestep(const args_info_rtkspectralonestep & args_info)
     spatialRegulWeighsReader->SetFileName(args_info.regul_spatial_weights_arg);
   }
 
+  // Read projections weights if given
+  IncidentSpectrumReaderType::Pointer projectionWeightsReader;
+  if (args_info.projection_weights_given)
+  {
+    projectionWeightsReader = IncidentSpectrumReaderType::New();
+    projectionWeightsReader->SetFileName(args_info.projection_weights_arg);
+  }
+
   // Create input: either an existing volume read from a file or a blank image
   typename itk::ImageSource<MaterialVolumesType>::Pointer inputFilter;
   if (args_info.input_given)
@@ -216,18 +224,33 @@ rtkspectralonestep(const args_info_rtkspectralonestep & args_info)
   if (args_info.subsets_arg != 1)
   {
     using ReorderProjectionsFilterType = rtk::ReorderProjectionsImageFilter<PhotonCountsType>;
-    typename ReorderProjectionsFilterType::Pointer reorder = ReorderProjectionsFilterType::New();
-    reorder->SetInput(photonCountsReader->GetOutput());
-    reorder->SetInputGeometry(geometryReader->GetOutputObject());
-    reorder->SetPermutation(rtk::ReorderProjectionsImageFilter<PhotonCountsType>::SHUFFLE);
-    reorder->Update();
-    mechlemOneStep->SetInputPhotonCounts(reorder->GetOutput());
-    mechlemOneStep->SetGeometry(reorder->GetOutputGeometry());
+    typename ReorderProjectionsFilterType::Pointer reorder_PhotonsCounts = ReorderProjectionsFilterType::New();
+    reorder_PhotonsCounts->SetInput(photonCountsReader->GetOutput());
+    reorder_PhotonsCounts->SetInputGeometry(geometryReader->GetOutputObject());
+    reorder_PhotonsCounts->SetPermutation(rtk::ReorderProjectionsImageFilter<PhotonCountsType>::SHUFFLE);
+    reorder_PhotonsCounts->Update();
+    mechlemOneStep->SetInputPhotonCounts(reorder_PhotonsCounts->GetOutput());
+    mechlemOneStep->SetGeometry(reorder_PhotonsCounts->GetOutputGeometry());
+    if (args_info.projection_weights_given)
+    {
+      // N.B: The permutation of photon counts and weights correspond because we use the SHUFFLE case mode initialized
+      // to zero (cf case (SHUFFLE): in ReorderProjectionsImageFilter.hxx)
+      using ReorderWeightsFilterType = rtk::ReorderProjectionsImageFilter<IncidentSpectrumType>;
+      typename ReorderWeightsFilterType::Pointer reorder_ProjectionWeights = ReorderWeightsFilterType::New();
+      reorder_ProjectionWeights->SetInput(projectionWeightsReader->GetOutput());
+      // Here we get geometry which has not been shuffled
+      reorder_ProjectionWeights->SetInputGeometry(geometryReader->GetOutputObject());
+      reorder_ProjectionWeights->SetPermutation(rtk::ReorderProjectionsImageFilter<IncidentSpectrumType>::SHUFFLE);
+      reorder_ProjectionWeights->Update();
+      mechlemOneStep->SetProjectionWeights(reorder_ProjectionWeights->GetOutput());
+    }
   }
   else
   {
     mechlemOneStep->SetInputPhotonCounts(photonCountsReader->GetOutput());
     mechlemOneStep->SetGeometry(geometryReader->GetOutputObject());
+    if (args_info.projection_weights_given)
+      mechlemOneStep->SetProjectionWeights(projectionWeightsReader->GetOutput());
   }
 
   REPORT_ITERATIONS(mechlemOneStep, MechlemFilterType, MaterialVolumesType);

applications/rtkspectralonestep/rtkspectralonestep.ggo

@@ -13,6 +13,7 @@ option "incident"              - "Incident spectrum file (mhd image)"
 option "attenuations"          a "Material attenuations file"                                        string yes
 option "mask"                  m "Apply a support binary mask: reconstruction kept null outside"     string no
 option "regul_spatial_weights" - "One-component image of spatial regularization weights"             string no
+option "projection_weights"    - "One-component image of projection weights (size of photon counts)" string no
 option "thresholds"            t "Lower threshold of bins, expressed in pulse height"                double yes  multiple
 option "subsets"               - "Number of subsets of projections (should not exceed 6)"            int    no   default="4"
 option "regul_weights"         - "Regularization parameters for each material"                       double multiple no

include/rtkMechlemOneStepSpectralReconstructionFilter.h

@@ -62,6 +62,8 @@ namespace rtk
  * Input3 [shape=Mdiamond];
  * Input4 [label="Input 4 (Spatial Regularization Weights)"];
  * Input4 [shape=Mdiamond];
+ * Input5 [label="Input 5 (Projection Weights)"];
+ * Input5 [shape=Mdiamond];
  * Output [label="Output (Material volumes)"];
  * Output [shape=Mdiamond];
  *
@@ -86,6 +88,7 @@ namespace rtk
  *            URL="\ref rtk::SeparableQuadraticSurrogateRegularizationImageFilter"];
  * MultiplyRegulGradient [ label="itk::MultiplyImageFilter" URL="\ref itk::MultiplyImageFilter" style=dashed];
  * MultiplyRegulHessian [ label="itk::MultiplyImageFilter" URL="\ref itk::MultiplyImageFilter" style=dashed];
+ * MultiplyGradientToBeBackProjected [ label="itk::MultiplyImageFilter" URL="\ref itk::MultiplyImageFilter" style=dashed];
  * AddGradients [ label="itk::AddImageFilter" URL="\ref itk::AddImageFilter"];
  * AddHessians [ label="rtk::AddMatrixAndDiagonalImageFilter" URL="\ref rtk::AddMatrixAndDiagonalImageFilter"];
  * Newton [ label="rtk::GetNewtonUpdateImageFilter" URL="\ref rtk::GetNewtonUpdateImageFilter"];
@@ -112,7 +115,9 @@ namespace rtk
  * Input4 -> MultiplyRegulGradient;
  * SQSRegul -> MultiplyRegulGradient;
  * MultiplyRegulGradient -> AddGradients;
- * BackProjectionGradients -> AddGradients;
+ * Input5 -> MultiplyGradientToBeBackProjected;
+ * BackProjectionGradients -> MultiplyGradientToBeBackProjected
+ * MultiplyGradientToBeBackProjected -> AddGradients;
  * AddGradients -> Newton;
  * Input4 -> MultiplyRegulHessian;
  * SQSRegul -> MultiplyRegulHessian;
@@ -259,6 +264,8 @@ class MechlemOneStepSpectralReconstructionFilter
   SetSupportMask(const SingleComponentImageType * support);
   void
   SetSpatialRegularizationWeights(const SingleComponentImageType * regweights);
+  void
+  SetProjectionWeights(const SingleComponentImageType * weiprojections);
 
   /** Set/Get for the regularization weights */
   itkSetMacro(RegularizationWeights, typename TOutputImage::PixelType);
@@ -309,6 +316,7 @@ class MechlemOneStepSpectralReconstructionFilter
   typename MultiplyFilterType::Pointer                         m_MultiplySupportFilter;
   typename MultiplyGradientFilterType::Pointer                 m_MultiplyRegulGradientsFilter;
   typename MultiplyGradientFilterType::Pointer                 m_MultiplyRegulHessiansFilter;
+  typename MultiplyGradientFilterType::Pointer                 m_MultiplyGradientToBeBackprojectedFilter;
 #endif
 
   /** The inputs of this filter have the same type but not the same meaning
@@ -336,6 +344,8 @@ class MechlemOneStepSpectralReconstructionFilter
   GetSupportMask();
   typename SingleComponentImageType::ConstPointer
   GetSpatialRegularizationWeights();
+  typename SingleComponentImageType::ConstPointer
+  GetProjectionWeights();
 
 #if !defined(ITK_WRAPPING_PARSER)
   /** Functions to instantiate forward and back projection filters with a different

include/rtkMechlemOneStepSpectralReconstructionFilter.hxx

@@ -53,6 +53,7 @@ MechlemOneStepSpectralReconstructionFilter<TOutputImage, TPhotonCounts, TSpectru
   m_SQSRegul = SQSRegularizationType::New();
   m_MultiplyRegulGradientsFilter = MultiplyGradientFilterType::New();
   m_MultiplyRegulHessiansFilter = MultiplyGradientFilterType::New();
+  m_MultiplyGradientToBeBackprojectedFilter = MultiplyGradientFilterType::New();
   m_WeidingerForward = WeidingerForwardModelType::New();
   m_NewtonFilter = NewtonFilterType::New();
   m_NesterovFilter = NesterovFilterType::New();
@@ -111,6 +112,14 @@ MechlemOneStepSpectralReconstructionFilter<TOutputImage, TPhotonCounts, TSpectru
   this->SetNthInput(4, const_cast<SingleComponentImageType *>(regweights));
 }
 
+template <class TOutputImage, class TPhotonCounts, class TSpectrum>
+void
+MechlemOneStepSpectralReconstructionFilter<TOutputImage, TPhotonCounts, TSpectrum>::SetProjectionWeights(
+  const SingleComponentImageType * weiprojections)
+{
+  this->SetNthInput(5, const_cast<SingleComponentImageType *>(weiprojections));
+}
+
 template <class TOutputImage, class TPhotonCounts, class TSpectrum>
 typename TOutputImage::ConstPointer
 MechlemOneStepSpectralReconstructionFilter<TOutputImage, TPhotonCounts, TSpectrum>::GetInputMaterialVolumes()
@@ -148,6 +157,14 @@ typename MechlemOneStepSpectralReconstructionFilter<TOutputImage, TPhotonCounts,
   return static_cast<const SingleComponentImageType *>(this->itk::ProcessObject::GetInput(4));
 }
 
+template <class TOutputImage, class TPhotonCounts, class TSpectrum>
+typename MechlemOneStepSpectralReconstructionFilter<TOutputImage, TPhotonCounts, TSpectrum>::SingleComponentImageType::
+  ConstPointer
+  MechlemOneStepSpectralReconstructionFilter<TOutputImage, TPhotonCounts, TSpectrum>::GetProjectionWeights()
+{
+  return static_cast<const SingleComponentImageType *>(this->itk::ProcessObject::GetInput(5));
+}
+
 template <class TOutputImage, class TPhotonCounts, class TSpectrum>
 void
 MechlemOneStepSpectralReconstructionFilter<TOutputImage, TPhotonCounts, TSpectrum>::SetForwardProjectionFilter(
@@ -293,6 +310,12 @@ MechlemOneStepSpectralReconstructionFilter<TOutputImage, TPhotonCounts, TSpectru
     const_cast<SingleComponentImageType *>(this->GetSpatialRegularizationWeights().GetPointer());
   if (inputPtr4)
     inputPtr4->SetRequestedRegion(inputPtr0->GetRequestedRegion());
+
+  // Input 5 is the image of weights for the gradient to be backprojected (optional)
+  typename SingleComponentImageType::Pointer inputPtr5 =
+    const_cast<SingleComponentImageType *>(this->GetProjectionWeights().GetPointer());
+  if (inputPtr5)
+    inputPtr5->SetRequestedRegion(inputPtr1->GetRequestedRegion());
 }
 
 template <class TOutputImage, class TPhotonCounts, class TSpectrum>
@@ -331,7 +354,16 @@ MechlemOneStepSpectralReconstructionFilter<TOutputImage, TPhotonCounts, TSpectru
   m_WeidingerForward->SetInputProjectionsOfOnes(m_SingleComponentForwardProjectionFilter->GetOutput());
 
   m_GradientsBackProjectionFilter->SetInput(0, m_GradientsSource->GetOutput());
-  m_GradientsBackProjectionFilter->SetInput(1, m_WeidingerForward->GetOutput1());
+  if (this->GetProjectionWeights().GetPointer() != nullptr)
+  {
+    m_MultiplyGradientToBeBackprojectedFilter->SetInput1(m_WeidingerForward->GetOutput1());
+    m_MultiplyGradientToBeBackprojectedFilter->SetInput2(this->GetProjectionWeights());
+    m_GradientsBackProjectionFilter->SetInput(1, m_MultiplyGradientToBeBackprojectedFilter->GetOutput());
+  }
+  else
+  {
+    m_GradientsBackProjectionFilter->SetInput(1, m_WeidingerForward->GetOutput1());
+  }
 
   m_HessiansBackProjectionFilter->SetInput(0, m_HessiansSource->GetOutput());
   m_HessiansBackProjectionFilter->SetInput(1, m_WeidingerForward->GetOutput2());