RTK  2.0.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  * http://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_ASSIGN(DualEnergyNegativeLogLikelihood);
47 
52  itkNewMacro( Self );
54 
58 
63 
64  // Constructor
66  {
68  }
69 
70  // Destructor
71  ~DualEnergyNegativeLogLikelihood() override = default;
72 
73  void Initialize() override
74  {
75  // This method computes the combined m_IncidentSpectrumAndDetectorResponseProduct
76  // from m_DetectorResponse and m_IncidentSpectrum
77  m_Thresholds.SetSize(2);
78  m_Thresholds[0]=1;
80 
81  // In dual energy CT, one possible design is to illuminate the object with
82  // either a low energy or a high energy spectrum, alternating between the two. In that case
83  // m_DetectorResponse has only one row (there is a single detector) and m_IncidentSpectrum
84  // has two rows (one for high energy, the other for low)
86  for (unsigned int i=0; i<2; i++)
87  for (unsigned int j=0; j<m_DetectorResponse.cols(); j++)
89  }
90 
91  // Not used with a simplex optimizer, but may be useful later
92  // for gradient based methods
93  void GetDerivative( const ParametersType & itkNotUsed(lineIntegrals),
94  DerivativeType & itkNotUsed(derivatives)) const override
95  {
96  itkExceptionMacro(<< "Not implemented");
97  }
98 
99  // Main method
100  MeasureType 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]) + (forward[i] - this->m_MeasuredData[i]) * (forward[i] - this->m_MeasuredData[i]) / variances[i];
115  measure *= 0.5;
116 
117  return measure;
118  }
119 
120  vnl_vector<double> GetVariances( const ParametersType & lineIntegrals ) const override
121  {
122  vnl_vector<double> attenuationFactors;
123  attenuationFactors.set_size(m_NumberOfEnergies);
124  GetAttenuationFactors(lineIntegrals, attenuationFactors);
125 
126  // Apply detector response, getting the lambdas
127  vnl_vector<double> intermediate;
128  intermediate.set_size(m_NumberOfEnergies);
129  for (unsigned int i=0; i<m_NumberOfEnergies; i++)
130  intermediate[i]=i+1;
131  intermediate = element_product(attenuationFactors, intermediate);
132  return (m_IncidentSpectrumAndDetectorResponseProduct * intermediate);
133  }
134 
135 protected:
137 
138 };
139 
140 }// namespace RTK
141 
142 #endif
Superclass::MaterialAttenuationsType MaterialAttenuationsType
MeasureType GetValue(const ParametersType &parameters) const override
Superclass::DetectorResponseType DetectorResponseType
void GetAttenuationFactors(const ParametersType &lineIntegrals, vnl_vector< double > &attenuationFactors) const
Superclass::IncidentSpectrumType IncidentSpectrumType
virtual vnl_vector< double > ForwardModel(const ParametersType &lineIntegrals) const
~DualEnergyNegativeLogLikelihood() override=default
vnl_vector< double > GetVariances(const ParametersType &lineIntegrals) const override
void GetDerivative(const ParametersType &, DerivativeType &) const override