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 19, 2024 · fc2ca4b · fc2ca4b
1 parent 6956127
commit fc2ca4b
Show file tree

Hide file tree

Showing 4 changed files with 40 additions and 4 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.variance_given)
+    forward->SetComputeVariances(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.variance_given)
+  {
+    TRY_AND_EXIT_ON_ITK_EXCEPTION(itk::WriteImage(forward->GetOutput(1), args_info.variance_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 "variance"   - "File name for the output variance corresponding to the Cramer Rao Lower Bound"   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.hxx b/include/rtkSpectralForwardModelImageFilter.hxx
@@ -275,10 +275,11 @@ SpectralForwardModelImageFilter<DecomposedProjectionsType,
                                 MaterialAttenuationsImageType>::GenerateOutputInformation()
 {
   Superclass::GenerateOutputInformation();
-
+  this->GetOutput(1)->SetLargestPossibleRegion(this->GetInputDecomposedProjections()->GetLargestPossibleRegion());
   this->m_NumberOfSpectralBins = this->GetInputMeasuredProjections()->GetVectorLength();
   this->m_NumberOfMaterials = this->GetInputDecomposedProjections()->GetVectorLength();
   this->m_NumberOfEnergies = this->GetInputIncidentSpectrum()->GetVectorLength();
+  this->GetOutput(1)->SetVectorLength(this->m_NumberOfMaterials);
 }
 
 template <typename DecomposedProjectionsType,
@@ -470,8 +471,20 @@ SpectralForwardModelImageFilter<DecomposedProjectionsType,
     // If requested, store the variances into output(1)
     if (m_ComputeVariances)
     {
-      output1It.Set(
-        itk::VariableLengthVector<double>(cost->GetVariances(in).data_block(), this->m_NumberOfSpectralBins));
+      if (m_IsSpectralCT)
+      {
+        rtk::ProjectionsDecompositionNegativeLogLikelihood::MeasuredDataType photon_counts;
+        photon_counts.SetSize(forward.size());
+        photon_counts.SetData(forward.data_block());
+        cost->SetMeasuredData(photon_counts);
+        cost->ComputeFischerMatrix(in);
+        output1It.Set(cost->GetCramerRaoLowerBound());
+      }
+      else
+      {
+        output1It.Set(
+          itk::VariableLengthVector<double>(cost->GetVariances(in).data_block(), this->m_NumberOfSpectralBins));
+      }
       ++output1It;
     }