OTB  9.0.0
Orfeo Toolbox
otbMDMDNMFImageFilter.hxx
Go to the documentation of this file.
1 /*
2  * Copyright (C) 2005-2022 Centre National d'Etudes Spatiales (CNES)
3  *
4  * This file is part of Orfeo Toolbox
5  *
6  * https://www.orfeo-toolbox.org/
7  *
8  * Licensed under the Apache License, Version 2.0 (the "License");
9  * you may not use this file except in compliance with the License.
10  * You may obtain a copy of the License at
11  *
12  * http://www.apache.org/licenses/LICENSE-2.0
13  *
14  * Unless required by applicable law or agreed to in writing, software
15  * distributed under the License is distributed on an "AS IS" BASIS,
16  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
17  * See the License for the specific language governing permissions and
18  * limitations under the License.
19  */
20 
21 #ifndef otbMDMDNMFImageFilter_hxx
22 #define otbMDMDNMFImageFilter_hxx
23 
24 #include "otbMDMDNMFImageFilter.h"
25 
26 namespace otb
27 {
28 template <class TInputImage, class TOutputImage>
30 {
31  m_MaxIter = 100;
32  m_CritStopValue = 00.5;
33  m_Delt = 1.;
34  m_LambdD = 0.01;
35  m_LambdS = 0.01;
36 }
37 
41 template <class TInputImage, class TOutputImage>
42 void MDMDNMFImageFilter<TInputImage, TOutputImage>::PrintSelf(std::ostream& os, itk::Indent indent) const
43 {
44  Superclass::PrintSelf(os, indent);
45  os << indent << "Input Endmembers Matrix: " << m_Endmembers << std::endl;
46 }
48 
49 template <class TInputImage, class TOutputImage>
51 {
52  Superclass::GenerateOutputInformation();
53  const unsigned int nbEndmembers = m_Endmembers.columns();
54  if (nbEndmembers != 0)
55  {
56  this->GetOutput()->SetNumberOfComponentsPerPixel(m_Endmembers.columns());
57  }
58  else
59  {
60  throw itk::ExceptionObject(__FILE__, __LINE__, "Endmembers matrix columns size required to know the output size", ITK_LOCATION);
61  }
62 }
63 
64 
65 template <class TInputImage, class TOutputImage>
67 {
68  M.set_size(m.rows() + 1, m.cols());
69 
70  for (unsigned int i = 0; i < M.rows() - 1; ++i)
71  {
72  M.set_row(i, m.get_row(i));
73  }
74  M.set_row(M.rows() - 1, 1.0);
75 }
76 
77 
78 template <class TInputImage, class TOutputImage>
79 double MDMDNMFImageFilter<TInputImage, TOutputImage>::Criterion(const MatrixType& X, const MatrixType& A, const MatrixType& S, const double& itkNotUsed(m_Delt),
80  const double& m_LambdS, const double& m_LambdD)
81 {
82  // This function computes
83  // f = ||Xsu-Asu.S||_2^2 -
84  // m_LambdS * ||1/J*ones - S||_2^2 +
85  // m_LambdD * (trace(transpose(A)*A)-1/L*trace(transpose(A).ones.A));
86  //
87  // = ||E1||_2^2 -
88  // m_LambdS * ||E2||_2^2 +
89  // m_LambdD * (trace(A'*A)-1/L*trace(E3);
90  //
91  // where || ||_2 is the L2 frobenius_norm,
92  // J is the number of endmembers and
93  // L is the number of spectral bands
94 
95  const unsigned int nbEndmembers = A.cols();
96  const unsigned nbBands = A.rows();
97  MatrixType Xsu, Asu, ones, E1, E2; //, E3;
98  double evalf, sumColsOfA, trE3;
99 
100  Xsu.set_size(X.rows() + 1, X.cols());
101  Asu.set_size(A.rows() + 1, A.cols());
102  AddOneRowOfOnes(A, Asu);
103  AddOneRowOfOnes(X, Xsu);
104 
105  //------- Computing the function blocs E1, E2 and E3 -------//
106  // Bloc 1
107  E1 = Xsu - Asu * S;
108 
109  // Bloc 2
110  E2 = S - 1. / ((double)nbEndmembers); // * ones - S;
111 
112  // Computing trace(transpose(A)*A)
113  double trAtA = 0;
114  for (unsigned int i = 0; i < A.columns(); ++i)
115  {
116  trAtA += A.get_column(i).two_norm() * A.get_column(i).two_norm();
117  }
118 
119  // Bloc 3: computing fast trE3 = trace(transpose(A)*ones*A)
120  trE3 = 0;
121  for (unsigned int j = 0; j < nbEndmembers; ++j)
122  {
123  sumColsOfA = A.get_column(j).sum();
124  trE3 += sumColsOfA * sumColsOfA;
125  }
126 
127 
128  /* for (int j=0; j<nbEndmembers; ++j)
129  {
130  sumColsOfA(j) = A.get_column(j).sum();
131  }
132 
133  E3.set_size(nbEndmembers, nbEndmembers);
134  for (int j1=0; j1<nbEndmembers; ++j1)
135  {
136  for (int j2=0; j2<nbEndmembers; ++j2)
137  {
138  E3(j1, j2) = sumColsOfA(j1)*sumColsOfA(j2);
139  }
140  }
141  */
142  //-------------------- Computing f --------------------------//
143  evalf = E1.frobenius_norm() * E1.frobenius_norm() - m_LambdS * E2.frobenius_norm() * E2.frobenius_norm() +
144  m_LambdD * (trAtA - (1. / (static_cast<double>(nbBands)) * trE3));
145  return evalf;
146 }
147 
148 template <class TInputImage, class TOutputImage>
149 void MDMDNMFImageFilter<TInputImage, TOutputImage>::EvalGradS(const MatrixType& X, const MatrixType& A, const MatrixType& S, const double& itkNotUsed(m_Delt),
150  const double& m_LambdS, MatrixType& gradS)
151 {
152  // Calculus of: gradS = 2 * Asu' * (Asu*S-Xsu) - lambd * 2 * (S - 1/J*ones(J, I));
153 
154  MatrixType Xsu, Asu;
155  Xsu.set_size(X.rows() + 1, X.cols());
156  Asu.set_size(A.rows() + 1, A.cols());
157  AddOneRowOfOnes(A, Asu);
158  AddOneRowOfOnes(X, Xsu);
159 
160  gradS = 2. * Asu.transpose() * (Asu * S - Xsu) - m_LambdS * 2. * (S - (1. / (static_cast<double>(S.rows()))));
161 }
162 
163 template <class TInputImage, class TOutputImage>
164 void MDMDNMFImageFilter<TInputImage, TOutputImage>::EvalGradA(const MatrixType& X, const MatrixType& A, const MatrixType& S, const double& itkNotUsed(m_Delt),
165  const double& m_LambdD, MatrixType& gradA)
166 {
167  // Compute gradA
168  // = (A*S-X) * (transpose(S)) + m_LambdD*(A-1/nbBands*ones(L, L)*A)
169  // = (A*S-X) * (transpose(S)) + m_LambdD*A- m_LambdD*/nbBands*ones(L, L)*A)
170 
171  VectorType sumColulmnsOfA;
172  sumColulmnsOfA.set_size(A.cols());
173  unsigned int nbBands = A.rows();
174 
175  // Computing ones*A : all rows are the same,
176  // (ones*A)[i,j] = sum(k=0 to nbBands, A[k,j])
177  for (unsigned int j = 0; j < A.cols(); ++j)
178  {
179  sumColulmnsOfA(j) = A.get_column(j).sum();
180  }
181 
182  // Initialize gradA
183  gradA = (A * S - X) * S.transpose();
184 
185  // 1st update of gradA
186  gradA += A * m_LambdD;
187 
188  // 2nd and last update id gradA, row by row (for performance reasons)
189  for (unsigned int i = 0; i < nbBands; ++i)
190  {
191  gradA.set_row(i, gradA.get_row(i) - (sumColulmnsOfA * m_LambdD / (static_cast<double>(nbBands))));
192  }
193 }
194 
195 
196 template <class TInputImage, class TOutputImage>
198  const double& sig, const double& betinit, const double& m_Delt, const double& m_LambdS,
199  const double& m_LambdD, MatrixType& variMat, double& alph, const bool isDirectEvalDirection)
200 
201 {
202  double evalf, newEvalf, bet;
203  evalf = Call(variMat, fixedMat, X, m_Delt, m_LambdS, m_LambdD, isDirectEvalDirection); // compute evalf
204 
205  MatrixType newVariMat = variMat - alph * gradVariMat;
206  SetNegativeCoefficientsToZero(newVariMat);
207  newEvalf = Call(newVariMat, fixedMat, X, m_Delt, m_LambdS, m_LambdD, isDirectEvalDirection); // compute newEvalf
208  bool bit = ArmijoTest(sig, variMat, newVariMat, evalf, newEvalf, gradVariMat, alph);
209 
210  int count = 1;
211  if (bit == true)
212  {
213  while (bit == true)
214  {
215  bet = pow(betinit, count);
216  alph = alph / bet;
217  newVariMat = variMat - alph * gradVariMat;
218  SetNegativeCoefficientsToZero(newVariMat);
219  newEvalf = Call(newVariMat, fixedMat, X, m_Delt, m_LambdS, m_LambdD, isDirectEvalDirection);
220  bit = ArmijoTest(sig, variMat, newVariMat, evalf, newEvalf, gradVariMat, alph);
221  ++count;
222  }
223  alph = alph * bet;
224  newVariMat = variMat - alph * gradVariMat;
225  SetNegativeCoefficientsToZero(newVariMat);
226  }
227  else
228  {
229  while (bit == false)
230  {
231  bet = pow(betinit, count);
232  alph = alph * bet;
233  newVariMat = variMat - alph * gradVariMat;
234  SetNegativeCoefficientsToZero(newVariMat);
235  newEvalf = Call(newVariMat, fixedMat, X, m_Delt, m_LambdS, m_LambdD, isDirectEvalDirection);
236  bit = ArmijoTest(sig, variMat, newVariMat, evalf, newEvalf, gradVariMat, alph);
237  ++count;
238  }
239  }
240  variMat = newVariMat;
241 }
242 
243 
244 template <class TInputImage, class TOutputImage>
245 double MDMDNMFImageFilter<TInputImage, TOutputImage>::Call(const MatrixType& variMat, const MatrixType& fixedMat, const MatrixType& X, const double& m_Delt,
246  const double& m_LambdS, const double& m_LambdD, const bool isDirectEvalDirection)
247 {
248  double evalf;
249  if (isDirectEvalDirection)
250  {
251  evalf = Criterion(X, variMat, fixedMat, m_Delt, m_LambdS, m_LambdD);
252  }
253  else
254  {
255  evalf = Criterion(X, fixedMat, variMat, m_Delt, m_LambdS, m_LambdD);
256  }
257  return evalf;
258 }
259 
260 
261 template <class TInputImage, class TOutputImage>
263 {
264  for (unsigned int i = 0; i < M.rows(); ++i)
265  {
266  for (unsigned int j = 0; j < M.cols(); ++j)
267  {
268  if (M(i, j) < 0)
269  M(i, j) = 0;
270  }
271  }
272 }
273 
274 template <class TInputImage, class TOutputImage>
276  const MatrixType& M2)
277 {
278  MatrixType M;
279  M.set_size(M1.rows(), M1.cols());
280  for (unsigned int i = 0; i < M.rows(); ++i)
281  {
282  for (unsigned int j = 0; j < M.cols(); ++j)
283  {
284  M(i, j) = M1(i, j) * M2(i, j);
285  }
286  }
287  return M;
288 }
289 
290 template <class TInputImage, class TOutputImage>
292 {
293  double sum = 0;
294  for (unsigned int i = 0; i < M.cols(); ++i)
295  {
296  sum += M.get_column(i).sum();
297  }
298  return sum;
299 }
300 
301 template <class TInputImage, class TOutputImage>
302 bool MDMDNMFImageFilter<TInputImage, TOutputImage>::ArmijoTest(const double& sig, const MatrixType variMat, const MatrixType& newVariMat, const double& evalf,
303  const double& newEvalf, const MatrixType& gradVariMat, const double& alph)
304 {
305  bool bit;
306 
307  // const unsigned int I = variMat.rows();
308  // const unsigned int J = variMat.cols();
309 
310  const MatrixType prod = TermByTermMatrixProduct(gradVariMat, newVariMat - variMat);
311  double sumProd = SumMatrixElements(prod);
312 
313  if (newEvalf - evalf <= sig * alph * sumProd)
314  bit = true;
315  else
316  bit = false;
317 
318  return bit;
319 }
320 
321 
325 template <class TInputImage, class TOutputImage>
327 {
328  this->AllocateOutputs();
329 
330  // Get the input and output pointers
331  InputPointerType inputPtr = const_cast<InputImageType*>(this->GetInput());
332  OutputPointerType outputPtr = this->GetOutput();
333 
334  inputPtr->Update();
335 
336  // Fill the output buffer with black pixels
338  itk::NumericTraits<OutputPixelType>::SetLength(zero, outputPtr->GetNumberOfComponentsPerPixel());
339  outputPtr->FillBuffer(zero);
340 
341  // Adaptation of contribution from A. Huck
342  // Convert input image into matrix
343  typename VectorImageToMatrixImageFilterType::Pointer inputImage2Matrix = VectorImageToMatrixImageFilterType::New();
344  inputImage2Matrix->SetInput(inputPtr);
345  inputImage2Matrix->Update();
346 
347  // useful const variables
348  const unsigned int nbEndmembers = m_Endmembers.columns();
349  const unsigned int nbComponentsPerPixel = inputPtr->GetNumberOfComponentsPerPixel();
350  const unsigned int nbPixels = inputPtr->GetLargestPossibleRegion().GetNumberOfPixels();
351 
352  // otbGenericMsgDebugMacro( << "nbEndmembers " << nbEndmembers );
353  // otbGenericMsgDebugMacro( << "nbComponentsPerPixel " << nbComponentsPerPixel );
354  // otbGenericMsgDebugMacro( << "nbPixels " << nbPixels );
355 
356  // Other declarations
357  double critA, critS, crit = 1;
358  const unsigned int divisorParam = 10;
359 
360  // Tuning the optimized function parameters
361  // const double m_Delt = 1.;
362  // const double m_LambdD = 0.01;
363  m_LambdS *= nbEndmembers;
364 
365  // Tuning the projected gradient parameters
366  double sig = 0.05;
367  double bet = 0.99;
368  double alphA = 1.;
369  double alphS = 1.;
370 
371  MatrixType X = inputImage2Matrix->GetMatrix();
372  // otbGenericMsgDebugMacro( << "X " << X );
373  //------- A and S declaration and initialization --------//
374  //---> These values fit with the ones chosen in the matlab
375  //---"nmf_validationOtb.m" function to validate the OTB code.
376 
377  // A is the endmembers matrix? Output from VCA input of mdmd
378  MatrixType A(this->m_Endmembers);
379 
380  MatrixType S;
381  S.set_size(nbEndmembers, nbPixels);
382  S.fill(1.);
383  // otbGenericMsgDebugMacro( << "S " << S.cols() );
384  //----------- Declaration of useful variables -----------//
385 
386  MatrixType Sold, Sdiff;
387  Sold.set_size(nbEndmembers, nbPixels);
388  Sdiff.set_size(nbEndmembers, nbPixels);
389 
390  MatrixType Aold, Adiff;
391  Aold.set_size(nbComponentsPerPixel, nbEndmembers);
392  Adiff.set_size(nbComponentsPerPixel, nbEndmembers);
393 
394  MatrixType gradS;
395  gradS.set_size(nbEndmembers, nbPixels);
396  gradS.fill(0);
397  MatrixType gradA;
398  gradA.set_size(nbComponentsPerPixel, nbEndmembers);
399  gradS.fill(0);
400  //----------------- Optimization loop -------------------//
401  // FA fA;
402  // FS fS;
403 
404  otbGenericMsgDebugMacro(<< "normX = " << X.fro_norm());
405 
406  unsigned int counter = 0;
407 
408  while ((crit > m_CritStopValue) && (counter < m_MaxIter))
409  {
410 
411  //---------------- Update S -----------------//
412  Sold = S;
413  // otbGenericMsgDebugMacro( << "gradS1 " << gradS );
414  this->EvalGradS(X, A, S, m_Delt, m_LambdS, gradS); // Compute gradS
415  // otbGenericMsgDebugMacro( << "m_LambdS " << m_LambdS );
416  // otbGenericMsgDebugMacro( << "gradS " << gradS );
417  if (counter % divisorParam == 0)
418  {
419 
420  otbGenericMsgDebugMacro(<< "Iteration = " << counter);
421  otbGenericMsgDebugMacro(<< "Criterion = " << Criterion(X, A, S, m_Delt, m_LambdS, m_LambdD));
422  otbGenericMsgDebugMacro(<< "statGradS = " << gradS.fro_norm());
423  otbGenericMsgDebugMacro(<< "gradS(0, 0) = " << gradS(0, 0));
424  otbGenericMsgDebugMacro(<< "alphS = " << alphS);
425  otbGenericMsgDebugMacro(<< "normS = " << S.fro_norm());
426  otbGenericMsgDebugMacro(<< "S(0, 0) = " << S(0, 0));
427  }
428 
429  ProjGradOneStep(X, A, gradS, sig, bet, m_Delt, m_LambdS, m_LambdD, S, alphS, false);
430 
431  if (counter % divisorParam == 0)
432  {
433  otbGenericMsgDebugMacro(<< "alphS = " << alphS);
434  otbGenericMsgDebugMacro(<< "normS = " << S.fro_norm());
435  otbGenericMsgDebugMacro(<< "S(0, 0) = " << S(0, 0));
436  }
437 
438  //---------------- Update A -----------------//
439  Aold = A;
440 
441  this->EvalGradA(X, A, S, m_Delt, m_LambdD, gradA); // Compute gradS
442 
443  if (counter % divisorParam == 0)
444  {
445  otbGenericMsgDebugMacro(<< "gradA(0, 0) = " << gradA(0, 0));
446  }
447 
448 
449  if (counter % divisorParam == 0)
450  {
451  otbGenericMsgDebugMacro(<< "alphA = " << alphA);
452  otbGenericMsgDebugMacro(<< "normA = " << A.fro_norm());
453  otbGenericMsgDebugMacro(<< "A(0, 0) = " << A(0, 0));
454  }
455  ProjGradOneStep(X, S, gradA, sig, bet, m_Delt, m_LambdS, m_LambdD, A, alphA, true);
456 
457  if (counter % divisorParam == 0)
458  {
459  otbGenericMsgDebugMacro(<< "alphA = " << alphA);
460  otbGenericMsgDebugMacro(<< "normA = " << A.fro_norm());
461  otbGenericMsgDebugMacro(<< "A(0, 0) = " << A(0, 0));
462  }
463 
464  //------------ crit evaluation --------------//
465  Adiff = Aold - A;
466  critA = Adiff.absolute_value_max();
467  Sdiff = Sold - S;
468  critS = Sdiff.absolute_value_max();
469  crit = std::max(critA, critS);
470 
471  if (counter % divisorParam == 0)
472  {
473  otbGenericMsgDebugMacro(<< "critA value: " << critA);
474  otbGenericMsgDebugMacro(<< "critS value: " << critS);
475  otbGenericMsgDebugMacro(<< "crit value: " << crit);
476  otbGenericMsgDebugMacro(<< "criterion value: " << Criterion(X, A, S, m_Delt, m_LambdS, m_LambdD));
477  }
478 
479  ++counter;
480  }
481 
482  //--- Putting the rows of in the bands of the output vector image ---//
483  // TODO
484  // Could be improved choosing an imageList for the abundance maps
485  // and a vector list for the endmember spectra (columns of A).
486 
487  itk::ImageRegionIterator<OutputImageType> outputIt(outputPtr, outputPtr->GetRequestedRegion());
488 
489  typename OutputImageType::PixelType vectorPixel;
490  vectorPixel.SetSize(outputPtr->GetNumberOfComponentsPerPixel());
491 
492  unsigned int i = 0;
493  outputIt.GoToBegin();
494  while (!outputIt.IsAtEnd())
495  {
496  for (unsigned int j = 0; j < nbEndmembers; ++j)
497  {
498  vectorPixel.SetElement(j, S(j, i));
499  }
500  outputIt.Set(vectorPixel);
501  ++i;
502  ++outputIt;
503  }
504 }
505 
509 template <class TInputImage, class TOutputImage>
511 {
512  // Call the superclass' implementation of this method
513  Superclass::GenerateInputRequestedRegion();
514 
515  // Get pointers to the input and output
516  InputPointerType inputPtr = const_cast<InputImageType*>(this->GetInput());
517  OutputPointerType outputPtr = this->GetOutput();
518 
519  if (!inputPtr || !outputPtr)
520  {
521  return;
522  }
523 
524  inputPtr->SetRequestedRegion(inputPtr->GetLargestPossibleRegion());
525 }
526 
527 } // end namespace otb
528 
529 #endif
otb::MDMDNMFImageFilter::Criterion
static double Criterion(const MatrixType &X, const MatrixType &A, const MatrixType &S, const double &delt, const double &lambdS, const double &lambdD)
Definition: otbMDMDNMFImageFilter.hxx:79
otb::MDMDNMFImageFilter::TermByTermMatrixProduct
static MatrixType TermByTermMatrixProduct(const MatrixType &M1, const MatrixType &M2)
Definition: otbMDMDNMFImageFilter.hxx:275
otb::MDMDNMFImageFilter::GenerateInputRequestedRegion
void GenerateInputRequestedRegion() override
Definition: otbMDMDNMFImageFilter.hxx:510
otb::MDMDNMFImageFilter::GenerateOutputInformation
void GenerateOutputInformation() override
Definition: otbMDMDNMFImageFilter.hxx:50
otb::MDMDNMFImageFilter::OutputPointerType
OutputImageType::Pointer OutputPointerType
Definition: otbMDMDNMFImageFilter.h:152
otb::VectorImageToMatrixImageFilter::Pointer
itk::SmartPointer< Self > Pointer
Definition: otbVectorImageToMatrixImageFilter.h:158
otb::MDMDNMFImageFilter::Call
static double Call(const MatrixType &variMat, const MatrixType &fixedMat, const MatrixType &X, const double &delt, const double &lambdS, const double &lambdD, const bool isDirectEvalDirection)
Definition: otbMDMDNMFImageFilter.hxx:245
otb::MDMDNMFImageFilter::AddOneRowOfOnes
static void AddOneRowOfOnes(const MatrixType &m, MatrixType &M)
Definition: otbMDMDNMFImageFilter.hxx:66
otb
The "otb" namespace contains all Orfeo Toolbox (OTB) classes.
Definition: otbJoinContainer.h:32
otbMDMDNMFImageFilter.h
otb::MDMDNMFImageFilter::MDMDNMFImageFilter
MDMDNMFImageFilter()
Definition: otbMDMDNMFImageFilter.hxx:29
otbGenericMsgDebugMacro
#define otbGenericMsgDebugMacro(x)
Definition: otbMacro.h:63
otb::MDMDNMFImageFilter::PrintSelf
void PrintSelf(std::ostream &os, itk::Indent indent) const override
Definition: otbMDMDNMFImageFilter.hxx:42
otb::MDMDNMFImageFilter::EvalGradS
static void EvalGradS(const MatrixType &X, const MatrixType &A, const MatrixType &S, const double &delt, const double &lambdS, MatrixType &gradS)
Definition: otbMDMDNMFImageFilter.hxx:149
otb::MDMDNMFImageFilter::OutputPixelType
OutputImageType::PixelType OutputPixelType
Definition: otbMDMDNMFImageFilter.h:157
otb::MDMDNMFImageFilter::SetNegativeCoefficientsToZero
static void SetNegativeCoefficientsToZero(MatrixType &M)
Definition: otbMDMDNMFImageFilter.hxx:262
otb::MDMDNMFImageFilter::ProjGradOneStep
static void ProjGradOneStep(const MatrixType &X, const MatrixType &fixedMat, const MatrixType &gradVariMat, const double &sig, const double &betinit, const double &delt, const double &lambdS, const double &lambdD, MatrixType &variMat, double &alph, const bool isDirectEvalDirection)
Definition: otbMDMDNMFImageFilter.hxx:197
otb::MDMDNMFImageFilter::ArmijoTest
static bool ArmijoTest(const double &sig, const MatrixType variMat, const MatrixType &newVariMat, const double &evalf, const double &newEvalf, const MatrixType &gradVariMat, const double &alph)
Definition: otbMDMDNMFImageFilter.hxx:302
otb::MDMDNMFImageFilter::InputPointerType
InputImageType::Pointer InputPointerType
Definition: otbMDMDNMFImageFilter.h:145
otb::MDMDNMFImageFilter::EvalGradA
static void EvalGradA(const MatrixType &X, const MatrixType &A, const MatrixType &S, const double &delt, const double &lambdD, MatrixType &gradA)
Definition: otbMDMDNMFImageFilter.hxx:164
otb::MDMDNMFImageFilter::InputImageType
TInputImage InputImageType
Definition: otbMDMDNMFImageFilter.h:141
otb::MDMDNMFImageFilter::SumMatrixElements
static double SumMatrixElements(const MatrixType &M)
Definition: otbMDMDNMFImageFilter.hxx:291
otb::MDMDNMFImageFilter::MatrixType
vnl_matrix< PrecisionType > MatrixType
Definition: otbMDMDNMFImageFilter.h:161
otb::MDMDNMFImageFilter::GenerateData
void GenerateData() override
Definition: otbMDMDNMFImageFilter.hxx:326
otb::MDMDNMFImageFilter::VectorType
vnl_vector< double > VectorType
Definition: otbMDMDNMFImageFilter.h:165