ENH: Add computation of the variance for spectral CT.

RTKConsortium · Sep 30, 2024 · cff0080 · cff0080
1 parent 7f5df73
commit cff0080
Show file tree

Hide file tree

Showing 6 changed files with 40 additions and 9 deletions.
diff --git a/applications/rtkdualenergyforwardmodel/rtkdualenergyforwardmodel.cxx b/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;

diff --git a/applications/rtkspectralforwardmodel/rtkspectralforwardmodel.cxx b/applications/rtkspectralforwardmodel/rtkspectralforwardmodel.cxx
@@ -111,18 +111,26 @@ main(int argc, char * argv[])
   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 variance
+  // 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;
 }
diff --git a/applications/rtkspectralforwardmodel/rtkspectralforwardmodel.ggo b/applications/rtkspectralforwardmodel/rtkspectralforwardmodel.ggo
@@ -10,3 +10,4 @@ option "incident"   - "Incident spectrum file"
 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
diff --git a/include/rtkDualEnergyNegativeLogLikelihood.h b/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++)

diff --git a/include/rtkSpectralForwardModelImageFilter.h b/include/rtkSpectralForwardModelImageFilter.h
@@ -112,6 +112,9 @@ class ITK_TEMPLATE_EXPORT SpectralForwardModelImageFilter
   typename DecomposedProjectionsType::ConstPointer
   GetOutputCramerRaoLowerBound();
 
+  typename MeasuredProjectionsType::ConstPointer
+  GetOutputVariances();
+
   itkSetMacro(Thresholds, ThresholdsType);
   itkGetMacro(Thresholds, ThresholdsType);
 

diff --git a/include/rtkSpectralForwardModelImageFilter.hxx b/include/rtkSpectralForwardModelImageFilter.hxx
@@ -256,6 +256,21 @@ SpectralForwardModelImageFilter<DecomposedProjectionsType,
   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,
@@ -494,8 +509,16 @@ 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;
     }