ENH: Add option to compute the Cramer-Rao Lower Bound.

Add option to obtain the Cramer-Rao Lower Bound (CRLB) when computing the photon counts using rtkSpectralForwardModelImageFilter. The CRLB corresponds to the lower bound of the variance of any unbiased estimator. Therefore, it could be used to estimate the variance that would be obtained after material decomposition.
RTKConsortium · Sep 30, 2024 · 7f5df73 · 7f5df73
1 parent 93ac0c3
commit 7f5df73
Show file tree

Hide file tree

Showing 5 changed files with 78 additions and 9 deletions.
diff --git a/applications/rtkspectralforwardmodel/rtkspectralforwardmodel.cxx b/applications/rtkspectralforwardmodel/rtkspectralforwardmodel.cxx
@@ -32,7 +32,6 @@ main(int argc, char * argv[])
 
   using PixelValueType = float;
   constexpr unsigned int Dimension = 3;
-
   using DecomposedProjectionType = itk::VectorImage<PixelValueType, Dimension>;
   using SpectralProjectionsType = itk::VectorImage<PixelValueType, Dimension>;
   using IncidentSpectrumImageType = itk::VectorImage<PixelValueType, Dimension - 1>;
@@ -110,11 +109,20 @@ main(int argc, char * argv[])
   forward->SetMaterialAttenuations(materialAttenuations);
   forward->SetThresholds(thresholds);
   forward->SetIsSpectralCT(true);
+  if (args_info.cramer_rao_given)
+    forward->SetComputeCramerRaoLowerBound(true);
 
   TRY_AND_EXIT_ON_ITK_EXCEPTION(forward->Update())
 
   // Write output
   TRY_AND_EXIT_ON_ITK_EXCEPTION(itk::WriteImage(forward->GetOutput(), args_info.output_arg))
 
+  // If requested, write the variance
+  if (args_info.cramer_rao_given)
+  {
+    TRY_AND_EXIT_ON_ITK_EXCEPTION(itk::WriteImage(forward->GetOutputCramerRaoLowerBound(), args_info.cramer_rao_arg))
+  }
+
+
   return EXIT_SUCCESS;
 }
diff --git a/applications/rtkspectralforwardmodel/rtkspectralforwardmodel.ggo b/applications/rtkspectralforwardmodel/rtkspectralforwardmodel.ggo
@@ -9,3 +9,4 @@ option "detector"   d "Detector response file"
 option "incident"   - "Incident spectrum file"                                            string                       yes
 option "attenuations" a "Material attenuations file"                                      string                       yes
 option "thresholds" t "Lower threshold of bins, expressed in pulse height"                double                       yes multiple
+option "cramer_rao" - "File name for the output Cramer Rao Lower Bound (estimation of the variance in the material decomposed images)"  string   no
diff --git a/include/rtkProjectionsDecompositionNegativeLogLikelihood.h b/include/rtkProjectionsDecompositionNegativeLogLikelihood.h
@@ -106,6 +106,20 @@ class ProjectionsDecompositionNegativeLogLikelihood : public itk::SingleValuedCo
     return diag;
   }
 
+  virtual itk::VariableLengthVector<float>
+  GetCramerRaoLowerBound()
+  {
+    // Return the inverses of the diagonal components (i.e. the inverse variances, to be used directly in WLS
+    // reconstruction)
+    itk::VariableLengthVector<double> diag;
+    diag.SetSize(m_NumberOfMaterials);
+    diag.Fill(0);
+
+    for (unsigned int mat = 0; mat < m_NumberOfMaterials; mat++)
+      diag[mat] = m_Fischer.GetInverse()[mat][mat];
+    return diag;
+  }
+
   virtual itk::VariableLengthVector<float>
   GetFischerMatrix()
   {

diff --git a/include/rtkSpectralForwardModelImageFilter.h b/include/rtkSpectralForwardModelImageFilter.h
@@ -109,6 +109,9 @@ class ITK_TEMPLATE_EXPORT SpectralForwardModelImageFilter
   typename MaterialAttenuationsImageType::ConstPointer
   GetMaterialAttenuations();
 
+  typename DecomposedProjectionsType::ConstPointer
+  GetOutputCramerRaoLowerBound();
+
   itkSetMacro(Thresholds, ThresholdsType);
   itkGetMacro(Thresholds, ThresholdsType);
 
@@ -127,6 +130,9 @@ class ITK_TEMPLATE_EXPORT SpectralForwardModelImageFilter
   itkSetMacro(ComputeVariances, bool);
   itkGetMacro(ComputeVariances, bool);
 
+  itkSetMacro(ComputeCramerRaoLowerBound, bool);
+  itkGetMacro(ComputeCramerRaoLowerBound, bool);
+
 protected:
   SpectralForwardModelImageFilter();
   ~SpectralForwardModelImageFilter() override = default;
@@ -165,8 +171,9 @@ class ITK_TEMPLATE_EXPORT SpectralForwardModelImageFilter
   unsigned int m_NumberOfIterations;
   unsigned int m_NumberOfMaterials;
   bool         m_OptimizeWithRestarts;
-  bool         m_IsSpectralCT;     // If not, it is dual energy CT
-  bool         m_ComputeVariances; // Only implemented for dual energy CT
+  bool         m_IsSpectralCT;               // If not, it is dual energy CT
+  bool         m_ComputeVariances;           // Only implemented for dual energy CT
+  bool         m_ComputeCramerRaoLowerBound; // Only implemented for spectral CT
 
 }; // end of class
 

diff --git a/include/rtkSpectralForwardModelImageFilter.hxx b/include/rtkSpectralForwardModelImageFilter.hxx
@@ -40,14 +40,19 @@ SpectralForwardModelImageFilter<DecomposedProjectionsType,
   m_NumberOfSpectralBins = 8;
   m_IsSpectralCT = true;
   m_ComputeVariances = false;
+  m_ComputeCramerRaoLowerBound = false;
 
-  this->SetNumberOfIndexedOutputs(2); // decomposed projections, inverse variance of decomposition noise
+  this->SetNumberOfIndexedOutputs(
+    3); // dual energy/ photon counts, variance of the distribution, Cramer-Rao lower bound.
 
-  // Dual energy projections (mean of the distribution)
+  // Dual energy projections (mean of the distribution) or photon counts in case of spectral CT.
   this->SetNthOutput(0, this->MakeOutput(0));
 
   // Variance of the distribution
   this->SetNthOutput(1, this->MakeOutput(1));
+
+  // Cramer-rao lower bound
+  this->SetNthOutput(2, this->MakeOutput(2));
 }
 
 template <typename DecomposedProjectionsType,
@@ -236,6 +241,21 @@ SpectralForwardModelImageFilter<DecomposedProjectionsType,
   return static_cast<const MaterialAttenuationsImageType *>(this->itk::ProcessObject::GetInput("MaterialAttenuations"));
 }
 
+template <typename DecomposedProjectionsType,
+          typename MeasuredProjectionsType,
+          typename IncidentSpectrumImageType,
+          typename DetectorResponseImageType,
+          typename MaterialAttenuationsImageType>
+typename DecomposedProjectionsType::ConstPointer
+SpectralForwardModelImageFilter<DecomposedProjectionsType,
+                                MeasuredProjectionsType,
+                                IncidentSpectrumImageType,
+                                DetectorResponseImageType,
+                                MaterialAttenuationsImageType>::GetOutputCramerRaoLowerBound()
+{
+  return static_cast<const DecomposedProjectionsType *>(this->GetOutput(2));
+}
+
 template <typename DecomposedProjectionsType,
           typename MeasuredProjectionsType,
           typename IncidentSpectrumImageType,
@@ -253,11 +273,13 @@ SpectralForwardModelImageFilter<DecomposedProjectionsType,
   switch (idx)
   {
     case 0:
-      output = (DecomposedProjectionsType::New()).GetPointer();
+      output = (MeasuredProjectionsType::New()).GetPointer();
       break;
     case 1:
-      output = (DecomposedProjectionsType::New()).GetPointer();
+      output = (MeasuredProjectionsType::New()).GetPointer();
       break;
+    case 2:
+      output = (DecomposedProjectionsType::New()).GetPointer();
   }
   return output.GetPointer();
 }
@@ -275,10 +297,11 @@ SpectralForwardModelImageFilter<DecomposedProjectionsType,
                                 MaterialAttenuationsImageType>::GenerateOutputInformation()
 {
   Superclass::GenerateOutputInformation();
-
   this->m_NumberOfSpectralBins = this->GetInputMeasuredProjections()->GetVectorLength();
   this->m_NumberOfMaterials = this->GetInputDecomposedProjections()->GetVectorLength();
   this->m_NumberOfEnergies = this->GetInputIncidentSpectrum()->GetVectorLength();
+  this->GetOutput(2)->SetLargestPossibleRegion(this->GetInputDecomposedProjections()->GetLargestPossibleRegion());
+  this->GetOutput(2)->SetVectorLength(this->m_NumberOfMaterials);
 }
 
 template <typename DecomposedProjectionsType,
@@ -400,7 +423,8 @@ SpectralForwardModelImageFilter<DecomposedProjectionsType,
   ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////
   // Walk the output projection stack. For each pixel, set the cost function's member variables and run the optimizer.
   itk::ImageRegionIterator<MeasuredProjectionsType>        output0It(this->GetOutput(0), outputRegionForThread);
-  itk::ImageRegionIterator<DecomposedProjectionsType>      output1It(this->GetOutput(1), outputRegionForThread);
+  itk::ImageRegionIterator<MeasuredProjectionsType>        output1It(this->GetOutput(1), outputRegionForThread);
+  itk::ImageRegionIterator<DecomposedProjectionsType>      output2It(this->GetOutput(2), outputRegionForThread);
   itk::ImageRegionConstIterator<DecomposedProjectionsType> inIt(this->GetInputDecomposedProjections(),
                                                                 outputRegionForThread);
 
@@ -475,6 +499,21 @@ SpectralForwardModelImageFilter<DecomposedProjectionsType,
       ++output1It;
     }
 
+    // If requested, compute the Cramer Rao Lower bound (only for spectral CT).
+    if (m_IsSpectralCT)
+    {
+      if (m_ComputeCramerRaoLowerBound)
+      {
+        rtk::ProjectionsDecompositionNegativeLogLikelihood::MeasuredDataType photon_counts;
+        photon_counts.SetSize(forward.size());
+        photon_counts.SetData(forward.data_block());
+        cost->SetMeasuredData(photon_counts);
+        cost->ComputeFischerMatrix(in);
+        output2It.Set(cost->GetCramerRaoLowerBound());
+        ++output2It;
+      }
+    }
+
     // Move forward
     ++output0It;
     ++inIt;