21 #ifndef otbMDMDNMFImageFilter_hxx
22 #define otbMDMDNMFImageFilter_hxx
28 template <
class TInputImage,
class TOutputImage>
32 m_CritStopValue = 00.5;
41 template <
class TInputImage,
class TOutputImage>
44 Superclass::PrintSelf(os, indent);
45 os << indent <<
"Input Endmembers Matrix: " << m_Endmembers << std::endl;
49 template <
class TInputImage,
class TOutputImage>
52 Superclass::GenerateOutputInformation();
53 const unsigned int nbEndmembers = m_Endmembers.columns();
54 if (nbEndmembers != 0)
56 this->GetOutput()->SetNumberOfComponentsPerPixel(m_Endmembers.columns());
60 throw itk::ExceptionObject(__FILE__, __LINE__,
"Endmembers matrix columns size required to know the output size", ITK_LOCATION);
65 template <
class TInputImage,
class TOutputImage>
68 M.set_size(m.rows() + 1, m.cols());
70 for (
unsigned int i = 0; i < M.rows() - 1; ++i)
72 M.set_row(i, m.get_row(i));
74 M.set_row(M.rows() - 1, 1.0);
78 template <
class TInputImage,
class TOutputImage>
80 const double& m_LambdS,
const double& m_LambdD)
95 const unsigned int nbEndmembers = A.cols();
96 const unsigned nbBands = A.rows();
98 double evalf, sumColsOfA, trE3;
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);
110 E2 = S - 1. / ((double)nbEndmembers);
114 for (
unsigned int i = 0; i < A.columns(); ++i)
116 trAtA += A.get_column(i).two_norm() * A.get_column(i).two_norm();
121 for (
unsigned int j = 0; j < nbEndmembers; ++j)
123 sumColsOfA = A.get_column(j).sum();
124 trE3 += sumColsOfA * sumColsOfA;
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));
148 template <
class TInputImage,
class TOutputImage>
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);
160 gradS = 2. * Asu.transpose() * (Asu * S - Xsu) - m_LambdS * 2. * (S - (1. / (
static_cast<double>(S.rows()))));
163 template <
class TInputImage,
class TOutputImage>
172 sumColulmnsOfA.set_size(A.cols());
173 unsigned int nbBands = A.rows();
177 for (
unsigned int j = 0; j < A.cols(); ++j)
179 sumColulmnsOfA(j) = A.get_column(j).sum();
183 gradA = (A * S - X) * S.transpose();
186 gradA += A * m_LambdD;
189 for (
unsigned int i = 0; i < nbBands; ++i)
191 gradA.set_row(i, gradA.get_row(i) - (sumColulmnsOfA * m_LambdD / (
static_cast<double>(nbBands))));
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)
202 double evalf, newEvalf, bet;
203 evalf = Call(variMat, fixedMat, X, m_Delt, m_LambdS, m_LambdD, isDirectEvalDirection);
205 MatrixType newVariMat = variMat - alph * gradVariMat;
206 SetNegativeCoefficientsToZero(newVariMat);
207 newEvalf = Call(newVariMat, fixedMat, X, m_Delt, m_LambdS, m_LambdD, isDirectEvalDirection);
208 bool bit = ArmijoTest(sig, variMat, newVariMat, evalf, newEvalf, gradVariMat, alph);
215 bet = pow(betinit, count);
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);
224 newVariMat = variMat - alph * gradVariMat;
225 SetNegativeCoefficientsToZero(newVariMat);
231 bet = pow(betinit, count);
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);
240 variMat = newVariMat;
244 template <
class TInputImage,
class TOutputImage>
246 const double& m_LambdS,
const double& m_LambdD,
const bool isDirectEvalDirection)
249 if (isDirectEvalDirection)
251 evalf = Criterion(X, variMat, fixedMat, m_Delt, m_LambdS, m_LambdD);
255 evalf = Criterion(X, fixedMat, variMat, m_Delt, m_LambdS, m_LambdD);
261 template <
class TInputImage,
class TOutputImage>
264 for (
unsigned int i = 0; i < M.rows(); ++i)
266 for (
unsigned int j = 0; j < M.cols(); ++j)
274 template <
class TInputImage,
class TOutputImage>
279 M.set_size(M1.rows(), M1.cols());
280 for (
unsigned int i = 0; i < M.rows(); ++i)
282 for (
unsigned int j = 0; j < M.cols(); ++j)
284 M(i, j) = M1(i, j) * M2(i, j);
290 template <
class TInputImage,
class TOutputImage>
294 for (
unsigned int i = 0; i < M.cols(); ++i)
296 sum += M.get_column(i).sum();
301 template <
class TInputImage,
class TOutputImage>
303 const double& newEvalf,
const MatrixType& gradVariMat,
const double& alph)
310 const MatrixType prod = TermByTermMatrixProduct(gradVariMat, newVariMat - variMat);
311 double sumProd = SumMatrixElements(prod);
313 if (newEvalf - evalf <= sig * alph * sumProd)
325 template <
class TInputImage,
class TOutputImage>
328 this->AllocateOutputs();
338 itk::NumericTraits<OutputPixelType>::SetLength(zero, outputPtr->GetNumberOfComponentsPerPixel());
339 outputPtr->FillBuffer(zero);
344 inputImage2Matrix->SetInput(inputPtr);
345 inputImage2Matrix->Update();
348 const unsigned int nbEndmembers = m_Endmembers.columns();
349 const unsigned int nbComponentsPerPixel = inputPtr->GetNumberOfComponentsPerPixel();
350 const unsigned int nbPixels = inputPtr->GetLargestPossibleRegion().GetNumberOfPixels();
357 double critA, critS, crit = 1;
358 const unsigned int divisorParam = 10;
363 m_LambdS *= nbEndmembers;
371 MatrixType X = inputImage2Matrix->GetMatrix();
381 S.set_size(nbEndmembers, nbPixels);
387 Sold.set_size(nbEndmembers, nbPixels);
388 Sdiff.set_size(nbEndmembers, nbPixels);
391 Aold.set_size(nbComponentsPerPixel, nbEndmembers);
392 Adiff.set_size(nbComponentsPerPixel, nbEndmembers);
395 gradS.set_size(nbEndmembers, nbPixels);
398 gradA.set_size(nbComponentsPerPixel, nbEndmembers);
406 unsigned int counter = 0;
408 while ((crit > m_CritStopValue) && (counter < m_MaxIter))
414 this->EvalGradS(X, A, S, m_Delt, m_LambdS, gradS);
417 if (counter % divisorParam == 0)
429 ProjGradOneStep(X, A, gradS, sig, bet, m_Delt, m_LambdS, m_LambdD, S, alphS,
false);
431 if (counter % divisorParam == 0)
441 this->EvalGradA(X, A, S, m_Delt, m_LambdD, gradA);
443 if (counter % divisorParam == 0)
449 if (counter % divisorParam == 0)
455 ProjGradOneStep(X, S, gradA, sig, bet, m_Delt, m_LambdS, m_LambdD, A, alphA,
true);
457 if (counter % divisorParam == 0)
466 critA = Adiff.absolute_value_max();
468 critS = Sdiff.absolute_value_max();
469 crit = std::max(critA, critS);
471 if (counter % divisorParam == 0)
487 itk::ImageRegionIterator<OutputImageType> outputIt(outputPtr, outputPtr->GetRequestedRegion());
489 typename OutputImageType::PixelType vectorPixel;
490 vectorPixel.SetSize(outputPtr->GetNumberOfComponentsPerPixel());
493 outputIt.GoToBegin();
494 while (!outputIt.IsAtEnd())
496 for (
unsigned int j = 0; j < nbEndmembers; ++j)
498 vectorPixel.SetElement(j, S(j, i));
500 outputIt.Set(vectorPixel);
509 template <
class TInputImage,
class TOutputImage>
513 Superclass::GenerateInputRequestedRegion();
519 if (!inputPtr || !outputPtr)
524 inputPtr->SetRequestedRegion(inputPtr->GetLargestPossibleRegion());