RTK  1.4.0
Reconstruction Toolkit
rtkSchlomka2008NegativeLogLikelihood.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 rtkSchlomka2008NegativeLogLikelihood_h
20 #define rtkSchlomka2008NegativeLogLikelihood_h
21 
23 #include "rtkMacro.h"
24 
25 #include <itkVectorImage.h>
27 #include <itkVariableSizeMatrix.h>
28 
29 namespace rtk
30 {
44 // We have to define the cost function first
46 {
47 public:
48 
53  itkNewMacro( Self );
55 
59 
64 
65  // Constructor
67  {
69  }
70 
71  // Destructor
72  virtual ~Schlomka2008NegativeLogLikelihood() ITK_OVERRIDE
73  {
74  }
75 
76  void Initialize() ITK_OVERRIDE
77  {
78  // This method computes the combined m_IncidentSpectrumAndDetectorResponseProduct
79  // from m_DetectorResponse and m_IncidentSpectrum
80 
81  // In spectral CT, m_DetectorResponse has as many rows as the number of bins,
82  // and m_IncidentSpectrum has only one row (there is only one spectrum illuminating
83  // the object)
85  for (unsigned int i=0; i<m_DetectorResponse.rows(); i++)
86  for (unsigned int j=0; j<m_DetectorResponse.cols(); j++)
88  }
89 
90  // Not used with a simplex optimizer, but may be useful later
91  // for gradient based methods
92  void GetDerivative( const ParametersType & lineIntegrals,
93  DerivativeType & derivatives ) const ITK_OVERRIDE
94  {
95  // Set the size of the derivatives vector
96  derivatives.set_size(m_NumberOfMaterials);
97 
98  // Get some required data
99  vnl_vector<double> attenuationFactors;
100  attenuationFactors.set_size(this->m_NumberOfEnergies);
101  GetAttenuationFactors(lineIntegrals, attenuationFactors);
102  vnl_vector<double> lambdas = ForwardModel(lineIntegrals);
103 
104  // Compute the vector of 1 - m_b / lambda_b
105  vnl_vector<double> weights;
106  weights.set_size(m_NumberOfSpectralBins);
107  for (unsigned int i=0; i<m_NumberOfSpectralBins; i++)
108  weights[i] = 1 - (m_MeasuredData[i] / lambdas[i]);
109 
110  // Prepare intermediate variables
111  vnl_vector<double> intermediate_a;
112  vnl_vector<double> partial_derivative_a;
113 
114  for (unsigned int a=0; a<m_NumberOfMaterials; a++)
115  {
116  // Compute the partial derivatives of lambda_b with respect to the material line integrals
117  intermediate_a = element_product(-attenuationFactors, m_MaterialAttenuations.get_column(a));
118  partial_derivative_a = m_IncidentSpectrumAndDetectorResponseProduct * intermediate_a;
119 
120  // Multiply them together element-wise, then dot product with the weights
121  derivatives[a] = dot_product(partial_derivative_a,weights);
122  }
123  }
124 
125  // Main method
126  MeasureType GetValue( const ParametersType & parameters ) const ITK_OVERRIDE
127  {
128  // Forward model: compute the expected number of counts in each bin
129  vnl_vector<double> forward = ForwardModel(parameters);
130 
131  long double measure = 0;
132  // Compute the negative log likelihood from the lambdas
133  for (unsigned int i=0; i<m_NumberOfSpectralBins; i++)
134  measure += forward[i] - std::log((long double)forward[i]) * m_MeasuredData[i];
135  return measure;
136  }
137 
138  void ComputeFischerMatrix(const ParametersType & lineIntegrals) ITK_OVERRIDE
139  {
140  // Get some required data
141  vnl_vector<double> attenuationFactors;
142  attenuationFactors.set_size(this->m_NumberOfEnergies);
143  GetAttenuationFactors(lineIntegrals, attenuationFactors);
144  vnl_vector<double> lambdas = ForwardModel(lineIntegrals);
145 
146  // Compute the vector of m_b / lambda_b²
147  vnl_vector<double> weights;
148  weights.set_size(m_NumberOfSpectralBins);
149  for (unsigned int i=0; i<m_NumberOfSpectralBins; i++)
150  weights[i] = m_MeasuredData[i] / (lambdas[i] * lambdas[i]);
151 
152  // Prepare intermediate variables
153  vnl_vector<double> intermediate_a;
154  vnl_vector<double> intermediate_a_prime;
155  vnl_vector<double> partial_derivative_a;
156  vnl_vector<double> partial_derivative_a_prime;
157 
158  // Compute the Fischer information matrix
160  for (unsigned int a=0; a<m_NumberOfMaterials; a++)
161  {
162  for (unsigned int a_prime=0; a_prime<m_NumberOfMaterials; a_prime++)
163  {
164  // Compute the partial derivatives of lambda_b with respect to the material line integrals
165  intermediate_a = element_product(-attenuationFactors, m_MaterialAttenuations.get_column(a));
166  intermediate_a_prime = element_product(-attenuationFactors, m_MaterialAttenuations.get_column(a_prime));
167 
168  partial_derivative_a = m_IncidentSpectrumAndDetectorResponseProduct * intermediate_a;
169  partial_derivative_a_prime = m_IncidentSpectrumAndDetectorResponseProduct * intermediate_a_prime;
170 
171  // Multiply them together element-wise, then dot product with the weights
172  partial_derivative_a_prime = element_product(partial_derivative_a, partial_derivative_a_prime);
173  m_Fischer[a][a_prime] = dot_product(partial_derivative_a_prime,weights);
174  }
175  }
176  }
177 
178 private:
179  Schlomka2008NegativeLogLikelihood(const Self &); //purposely not implemented
180  void operator = (const Self &); //purposely not implemented
181 
182 };
183 
184 }// namespace RTK
185 
186 #endif
void ComputeFischerMatrix(const ParametersType &lineIntegrals) override
rtk::ProjectionsDecompositionNegativeLogLikelihood Superclass
void GetAttenuationFactors(const ParametersType &lineIntegrals, vnl_vector< double > &attenuationFactors) const
virtual vnl_vector< double > ForwardModel(const ParametersType &lineIntegrals) const
MeasureType GetValue(const ParametersType &parameters) const override
void GetDerivative(const ParametersType &lineIntegrals, DerivativeType &derivatives) const override
Superclass::MaterialAttenuationsType MaterialAttenuationsType