RTK  2.4.1
Reconstruction Toolkit
rtkDualEnergyNegativeLogLikelihood.h
Go to the documentation of this file.
1 /*=========================================================================
2  *
3  * Copyright RTK Consortium
4  *
5  * Licensed under the Apache License, Version 2.0 (the "License");
6  * you may not use this file except in compliance with the License.
7  * You may obtain a copy of the License at
8  *
9  * https://www.apache.org/licenses/LICENSE-2.0.txt
10  *
11  * Unless required by applicable law or agreed to in writing, software
12  * distributed under the License is distributed on an "AS IS" BASIS,
13  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14  * See the License for the specific language governing permissions and
15  * limitations under the License.
16  *
17  *=========================================================================*/
18 
19 #ifndef rtkDualEnergyNegativeLogLikelihood_h
20 #define rtkDualEnergyNegativeLogLikelihood_h
21 
23 #include "rtkMacro.h"
24 
25 #include <itkVectorImage.h>
27 #include <itkVariableSizeMatrix.h>
28 
29 namespace rtk
30 {
42 // We have to define the cost function first
44 {
45 public:
46  ITK_DISALLOW_COPY_AND_MOVE(DualEnergyNegativeLogLikelihood);
47 
52  itkNewMacro(Self);
54 
58 
63 
64  // Constructor
66 
67  // Destructor
68  ~DualEnergyNegativeLogLikelihood() override = default;
69 
70  void
71  Initialize() override
72  {
73  // This method computes the combined m_IncidentSpectrumAndDetectorResponseProduct
74  // from m_DetectorResponse and m_IncidentSpectrum
75  m_Thresholds.SetSize(2);
76  m_Thresholds[0] = 1;
78 
79  // In dual energy CT, one possible design is to illuminate the object with
80  // either a low energy or a high energy spectrum, alternating between the two. In that case
81  // m_DetectorResponse has only one row (there is a single detector) and m_IncidentSpectrum
82  // has two rows (one for high energy, the other for low)
84  for (unsigned int i = 0; i < 2; i++)
85  for (unsigned int j = 0; j < m_DetectorResponse.cols(); j++)
87  }
88 
89  // Not used with a simplex optimizer, but may be useful later
90  // for gradient based methods
91  void
92  GetDerivative(const ParametersType & itkNotUsed(lineIntegrals),
93  DerivativeType & itkNotUsed(derivatives)) const override
94  {
95  itkExceptionMacro(<< "Not implemented");
96  }
97 
98  // Main method
100  GetValue(const ParametersType & parameters) const override
101  {
102  // Forward model: compute the expected total energy measured by the detector for each spectrum
103  vnl_vector<double> forward = ForwardModel(parameters);
104  vnl_vector<double> variances = GetVariances(parameters);
105 
106  long double measure = 0;
107  // TODO: Improve this estimation
108  // We assume that the variance of the integrated energy is equal to the mean
109  // From equation (5) of "Cramer-Rao lower bound of basis image noise in multiple-energy x-ray imaging",
110  // PMB 2009, Roessl et al, we replace the variance with the mean
111 
112  // Compute the negative log likelihood from the expectedEnergies
113  for (unsigned int i = 0; i < this->m_NumberOfMaterials; i++)
114  measure += std::log((long double)variances[i]) +
115  (forward[i] - this->m_MeasuredData[i]) * (forward[i] - this->m_MeasuredData[i]) / variances[i];
116  measure *= 0.5;
117 
118  return measure;
119  }
120 
121  vnl_vector<double>
122  GetVariances(const ParametersType & lineIntegrals) const override
123  {
124  vnl_vector<double> attenuationFactors;
125  attenuationFactors.set_size(m_NumberOfEnergies);
126  GetAttenuationFactors(lineIntegrals, attenuationFactors);
127 
128  // Apply detector response, getting the lambdas
129  vnl_vector<double> intermediate;
130  intermediate.set_size(m_NumberOfEnergies);
131  for (unsigned int i = 0; i < m_NumberOfEnergies; i++)
132  intermediate[i] = i + 1;
133  intermediate = element_product(attenuationFactors, intermediate);
134  return (m_IncidentSpectrumAndDetectorResponseProduct * intermediate);
135  }
136 
137 protected:
139 };
140 
141 } // namespace rtk
142 
143 #endif
virtual vnl_vector< double > ForwardModel(const ParametersType &lineIntegrals) const
MeasureType GetValue(const ParametersType &parameters) const override
Superclass::MaterialAttenuationsType MaterialAttenuationsType
void GetAttenuationFactors(const ParametersType &lineIntegrals, vnl_vector< double > &attenuationFactors) const
Superclass::DetectorResponseType DetectorResponseType
Superclass::IncidentSpectrumType IncidentSpectrumType
void GetDerivative(const ParametersType &, DerivativeType &) const override
vnl_vector< double > GetVariances(const ParametersType &lineIntegrals) const override
~DualEnergyNegativeLogLikelihood() override=default