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