RTK  2.6.0
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);
53 #ifdef itkOverrideGetNameOfClassMacro
54  itkOverrideGetNameOfClassMacro(DualEnergyNegativeLogLikelihood);
55 #else
57 #endif
58 
62 
67 
68  // Constructor
70 
71  // Destructor
72  ~DualEnergyNegativeLogLikelihood() override = default;
73 
74  void
75  Initialize() override
76  {
77  // This method computes the combined m_IncidentSpectrumAndDetectorResponseProduct
78  // from m_DetectorResponse and m_IncidentSpectrum
79  m_Thresholds.SetSize(2);
80  m_Thresholds[0] = 1;
82 
83  // In dual energy CT, one possible design is to illuminate the object with
84  // either a low energy or a high energy spectrum, alternating between the two. In that case
85  // m_DetectorResponse has only one row (there is a single detector) and m_IncidentSpectrum
86  // has two rows (one for high energy, the other for low)
88  for (unsigned int i = 0; i < 2; i++)
89  for (unsigned int j = 0; j < m_DetectorResponse.cols(); j++)
91  }
92 
93  // Not used with a simplex optimizer, but may be useful later
94  // for gradient based methods
95  void
96  GetDerivative(const ParametersType & itkNotUsed(lineIntegrals),
97  DerivativeType & itkNotUsed(derivatives)) const override
98  {
99  itkExceptionMacro(<< "Not implemented");
100  }
101 
102  // Main method
104  GetValue(const ParametersType & parameters) const override
105  {
106  // Forward model: compute the expected total energy measured by the detector for each spectrum
107  vnl_vector<double> forward = ForwardModel(parameters);
108  vnl_vector<double> variances = GetVariances(parameters);
109 
110  long double measure = 0;
111  // TODO: Improve this estimation
112  // We assume that the variance of the integrated energy is equal to the mean
113  // From equation (5) of "Cramer-Rao lower bound of basis image noise in multiple-energy x-ray imaging",
114  // PMB 2009, Roessl et al, we replace the variance with the mean
115 
116  // Compute the negative log likelihood from the expectedEnergies
117  for (unsigned int i = 0; i < this->m_NumberOfMaterials; i++)
118  measure += std::log((long double)variances[i]) +
119  (forward[i] - this->m_MeasuredData[i]) * (forward[i] - this->m_MeasuredData[i]) / variances[i];
120  measure *= 0.5;
121 
122  return measure;
123  }
124 
125  vnl_vector<double>
126  GetVariances(const ParametersType & lineIntegrals) const override
127  {
128  vnl_vector<double> attenuationFactors;
129  attenuationFactors.set_size(m_NumberOfEnergies);
130  GetAttenuationFactors(lineIntegrals, attenuationFactors);
131 
132  // Apply detector response, getting the lambdas
133  vnl_vector<double> intermediate;
134  intermediate.set_size(m_NumberOfEnergies);
135  for (unsigned int i = 0; i < m_NumberOfEnergies; i++)
136  intermediate[i] = i + 1;
137  intermediate = element_product(attenuationFactors, intermediate);
138  return (m_IncidentSpectrumAndDetectorResponseProduct * intermediate);
139  }
140 
141 protected:
143 };
144 
145 } // namespace rtk
146 
147 #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