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