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

applications/CMakeLists.txt

@@ -89,6 +89,7 @@ add_subdirectory(rtkmaskcollimation)
 add_subdirectory(rtkspectraldenoiseprojections)
 add_subdirectory(rtkprojectionmatrix)
 add_subdirectory(rtkvectorconjugategradient)
+add_subdirectory(rtkdualenergyforwardmodel)
 
 add_subdirectory(rtkcheckimagequality)
 

applications/rtkdualenergyforwardmodel/rtkdualenergyforwardmodel.cxx

@@ -134,9 +134,7 @@ main(int argc, char * argv[])
   // If requested, write the variances
   if (args_info.variances_given)
   {
-    writer->SetInput(forward->GetOutput(1));
-    writer->SetFileName(args_info.variances_arg);
-    writer->Update();
+    TRY_AND_EXIT_ON_ITK_EXCEPTION(itk::WriteImage(forward->GetOutputVariances(), args_info.variances_arg))
   }
 
   return EXIT_SUCCESS;

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,28 @@ main(int argc, char * argv[])
   forward->SetMaterialAttenuations(materialAttenuations);
   forward->SetThresholds(thresholds);
   forward->SetIsSpectralCT(true);
+  if (args_info.cramer_rao_given)
+    forward->SetComputeCramerRaoLowerBound(true);
+  if (args_info.variances_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 Cramer-Rao lower bound
+  if (args_info.cramer_rao_given)
+  {
+    TRY_AND_EXIT_ON_ITK_EXCEPTION(itk::WriteImage(forward->GetOutputCramerRaoLowerBound(), args_info.cramer_rao_arg))
+  }
+
+  // If requested, write the variance
+  if (args_info.variances_given)
+  {
+    TRY_AND_EXIT_ON_ITK_EXCEPTION(itk::WriteImage(forward->GetOutputVariances(), args_info.variances_arg))
+  }
+
+
   return EXIT_SUCCESS;
 }

applications/rtkspectralforwardmodel/rtkspectralforwardmodel.ggo

@@ -9,3 +9,5 @@ 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
+option "variances"  - "Output variances of photon counts, file name"                      string                       no

include/rtkDualEnergyNegativeLogLikelihood.h

@@ -108,10 +108,8 @@ class DualEnergyNegativeLogLikelihood : public rtk::ProjectionsDecompositionNega
     vnl_vector<double> variances = GetVariances(parameters);
 
     long double measure = 0;
-    // TODO: Improve this estimation
-    // We assume that the variance of the integrated energy is equal to the mean
     // From equation (5) of "Cramer-Rao lower bound of basis image noise in multiple-energy x-ray imaging",
-    // PMB 2009, Roessl et al, we replace the variance with the mean
+    // PMB 2009, Roessl et al.
 
     // Compute the negative log likelihood from the expectedEnergies
     for (unsigned int i = 0; i < this->m_NumberOfMaterials; i++)

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()
   {

include/rtkSpectralForwardModelImageFilter.h

@@ -109,6 +109,12 @@ class ITK_TEMPLATE_EXPORT SpectralForwardModelImageFilter
   typename MaterialAttenuationsImageType::ConstPointer
   GetMaterialAttenuations();
 
+  typename DecomposedProjectionsType::ConstPointer
+  GetOutputCramerRaoLowerBound();
+
+  typename MeasuredProjectionsType::ConstPointer
+  GetOutputVariances();
+
   itkSetMacro(Thresholds, ThresholdsType);
   itkGetMacro(Thresholds, ThresholdsType);
 
@@ -127,6 +133,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 +174,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
 

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,36 @@ 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,
+          typename DetectorResponseImageType,
+          typename MaterialAttenuationsImageType>
+typename MeasuredProjectionsType::ConstPointer
+SpectralForwardModelImageFilter<DecomposedProjectionsType,
+                                MeasuredProjectionsType,
+                                IncidentSpectrumImageType,
+                                DetectorResponseImageType,
+                                MaterialAttenuationsImageType>::GetOutputVariances()
+{
+  return static_cast<const MeasuredProjectionsType *>(this->GetOutput(1));
+}
+
 template <typename DecomposedProjectionsType,
           typename MeasuredProjectionsType,
           typename IncidentSpectrumImageType,
@@ -253,11 +288,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 +312,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 +438,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);
 
@@ -470,11 +509,34 @@ 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)
+      {
+        // For spectral CT, Poisson noise is assumed therefore the variances is equal to the photon counts.
+        output1It.Set(outputPixel);
+      }
+      else
+      {
+        output1It.Set(
+          itk::VariableLengthVector<double>(cost->GetVariances(in).data_block(), this->m_NumberOfSpectralBins));
+      }
       ++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;