RTK  1.4.0
Reconstruction Toolkit
rtkProjectionsDecompositionNegativeLogLikelihood.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 rtkProjectionsDecompositionNegativeLogLikelihood_h
20 #define rtkProjectionsDecompositionNegativeLogLikelihood_h
21 
23 #include <itkVectorImage.h>
25 #include <itkVariableSizeMatrix.h>
26 #include "rtkMacro.h"
27 
28 namespace rtk
29 {
38 // We have to define the cost function first
40 {
41 public:
42 
47  itkNewMacro( Self );
49 
50 // enum { SpaceDimension=m_NumberOfMaterials };
51 
52  typedef Superclass::ParametersType ParametersType;
53  typedef Superclass::DerivativeType DerivativeType;
54  typedef Superclass::MeasureType MeasureType;
55 
56  typedef vnl_matrix<double> DetectorResponseType;
57  typedef vnl_matrix<double> MaterialAttenuationsType;
58  typedef vnl_matrix<float> IncidentSpectrumType;
62 
63  // Constructor
65  {
68  m_Initialized = false;
69  }
70 
71  // Destructor
73  {
74  }
75 
76  MeasureType GetValue( const ParametersType & itkNotUsed(parameters)) const ITK_OVERRIDE {
77  long double measure = 0;
78  return measure;
79  }
80  void GetDerivative( const ParametersType & itkNotUsed(lineIntegrals),
81  DerivativeType & itkNotUsed(derivatives)) const ITK_OVERRIDE {itkExceptionMacro(<< "Not implemented");}
82  virtual void Initialize() {}
83 
85  {
86  // Return the inverses of the diagonal components (i.e. the inverse variances, to be used directly in WLS reconstruction)
88  diag.SetSize(m_NumberOfMaterials);
89  diag.Fill(0);
90 
91  for (unsigned int mat=0; mat<m_NumberOfMaterials; mat++)
92  diag[mat] = 1./m_Fischer.GetInverse()[mat][mat];
93  return diag;
94  }
95 
97  {
98  // Return the whole Fischer information matrix
100  fischer.SetSize(m_NumberOfMaterials * m_NumberOfMaterials);
101  fischer.Fill(0);
102 
103  for (unsigned int i=0; i<m_NumberOfMaterials; i++)
104  for (unsigned int j=0; j<m_NumberOfMaterials; j++)
105  fischer[i * m_NumberOfMaterials + j] = m_Fischer[i][j];
106  return fischer;
107  }
108 
109  virtual void ComputeFischerMatrix(const ParametersType & itkNotUsed(lineIntegrals)) {}
110 
111  unsigned int GetNumberOfParameters(void) const ITK_OVERRIDE
112  {
113  return m_NumberOfMaterials;
114  }
115 
116  virtual vnl_vector<double> ForwardModel(const ParametersType & lineIntegrals) const
117  {
118  vnl_vector<double> attenuationFactors;
119  attenuationFactors.set_size(m_NumberOfEnergies);
120  GetAttenuationFactors(lineIntegrals, attenuationFactors);
121 
122  // Apply detector response, getting the lambdas
123  return (m_IncidentSpectrumAndDetectorResponseProduct * attenuationFactors);
124  }
125 
126  void GetAttenuationFactors(const ParametersType & lineIntegrals, vnl_vector<double> & attenuationFactors) const
127  {
128  // Apply attenuation at each energy
129  vnl_vector<double> vnlLineIntegrals;
130 
131  // Initialize the line integrals vnl vector
132  vnlLineIntegrals.set_size(m_NumberOfMaterials);
133  for (unsigned int m=0; m<m_NumberOfMaterials; m++)
134  vnlLineIntegrals[m] = lineIntegrals[m];
135 
136  // Apply the material attenuations matrix
137  attenuationFactors = this->m_MaterialAttenuations * vnlLineIntegrals;
138 
139  // Compute the negative exponential
140  for (unsigned int energy = 0; energy<m_NumberOfEnergies; energy++)
141  {
142  attenuationFactors[energy] = std::exp(-attenuationFactors[energy]);
143  }
144  }
145 
147  {
149  initialGuess.SetSize(m_NumberOfMaterials);
150 
151  // Compute the mean attenuation in each bin, weighted by the input spectrum
152  // Needs to be done for each pixel, since the input spectrum is variable
153  MeanAttenuationInBinType MeanAttenuationInBin;
154  MeanAttenuationInBin.SetSize(this->m_NumberOfMaterials, this->m_NumberOfSpectralBins);
155  MeanAttenuationInBin.Fill(0);
156 
157  for (unsigned int mat = 0; mat<this->m_NumberOfMaterials; mat++)
158  {
159  for (unsigned int bin=0; bin<m_NumberOfSpectralBins; bin++)
160  {
161  double accumulate = 0;
162  double accumulateWeights = 0;
163  for (int energy=m_Thresholds[bin]-1; (energy<m_Thresholds[bin+1]) && (energy < (int)(this->m_MaterialAttenuations.rows())); energy++)
164  {
165  accumulate += this->m_MaterialAttenuations[energy][mat] * this->m_IncidentSpectrum[0][energy];
166  accumulateWeights += this->m_IncidentSpectrum[0][energy];
167  }
168  MeanAttenuationInBin[mat][bin] = accumulate / accumulateWeights;
169  }
170  }
171 
172  for (unsigned int mat = 0; mat<m_NumberOfMaterials; mat++)
173  {
174  // Initialise to a very high value
175  initialGuess[mat] = 1e10;
176  for (unsigned int bin = 0; bin<m_NumberOfSpectralBins; bin++)
177  {
178  // Compute the length of current material required to obtain the attenuation
179  // observed in current bin. Keep only the minimum among all bins
180  double requiredLength = this->BinwiseLogTransform()[bin] / MeanAttenuationInBin[mat][bin];
181  if (initialGuess[mat] > requiredLength)
182  initialGuess[mat] = requiredLength;
183  }
184  }
185 
186  return initialGuess;
187  }
188 
190  {
191  itk::VariableLengthVector<double> logTransforms;
192  logTransforms.SetSize(m_NumberOfSpectralBins);
193 
194  vnl_vector<double> ones, nonAttenuated;
195  ones.set_size(m_NumberOfEnergies);
196  ones.fill(1.0);
197 
198  // The way m_IncidentSpectrumAndDetectorResponseProduct works is
199  // it is mutliplied by the vector of attenuations factors (here
200  // filled with ones, since we want the non-attenuated signal)
201  nonAttenuated = m_IncidentSpectrumAndDetectorResponseProduct * ones;
202 
203  for (unsigned int i=0; i<m_MeasuredData.GetSize(); i++)
204  {
205  // Divide by the actually measured photon counts and apply log
206  if (m_MeasuredData[i] > 0)
207  logTransforms[i] = log(nonAttenuated[i] / m_MeasuredData[i]);
208  }
209 
210  return logTransforms;
211  }
212 
213  virtual vnl_vector<double> GetVariances( const ParametersType & itkNotUsed(lineIntegrals) ) const
214  {
215  vnl_vector<double> meaninglessResult;
216  meaninglessResult.set_size(m_NumberOfSpectralBins);
217  meaninglessResult.fill(0.);
218  return(meaninglessResult);
219  }
220 
221  itkSetMacro(MeasuredData, MeasuredDataType)
222  itkGetMacro(MeasuredData, MeasuredDataType)
223 
224  itkSetMacro(DetectorResponse, DetectorResponseType)
225  itkGetMacro(DetectorResponse, DetectorResponseType)
226 
227  itkSetMacro(MaterialAttenuations, MaterialAttenuationsType)
228  itkGetMacro(MaterialAttenuations, MaterialAttenuationsType)
229 
230  itkSetMacro(NumberOfEnergies, unsigned int)
231  itkGetMacro(NumberOfEnergies, unsigned int)
232 
233  itkSetMacro(NumberOfMaterials, unsigned int)
234  itkGetMacro(NumberOfMaterials, unsigned int)
235 
236  itkSetMacro(IncidentSpectrum, IncidentSpectrumType)
237  itkGetMacro(IncidentSpectrum, IncidentSpectrumType)
238 
239  itkSetMacro(NumberOfSpectralBins, unsigned int)
240  itkGetMacro(NumberOfSpectralBins, unsigned int)
241 
242  itkSetMacro(Thresholds, ThresholdsType)
243  itkGetMacro(Thresholds, ThresholdsType)
244 
245 protected:
246  MaterialAttenuationsType m_MaterialAttenuations;
247  DetectorResponseType m_DetectorResponse;
248  MeasuredDataType m_MeasuredData;
249  ThresholdsType m_Thresholds;
250  IncidentSpectrumType m_IncidentSpectrum;
252  unsigned int m_NumberOfEnergies;
253  unsigned int m_NumberOfMaterials;
257 
258 private:
259  ProjectionsDecompositionNegativeLogLikelihood(const Self &); //purposely not implemented
260  void operator = (const Self &); //purposely not implemented
261 
262 };
263 
264 }// namespace RTK
265 
266 #endif
void GetDerivative(const ParametersType &, DerivativeType &) const override
void GetAttenuationFactors(const ParametersType &lineIntegrals, vnl_vector< double > &attenuationFactors) const
virtual vnl_vector< double > ForwardModel(const ParametersType &lineIntegrals) const
#define itkSetMacro(name, type)
virtual vnl_vector< double > GetVariances(const ParametersType &) const