Skip to content
Snippets Groups Projects
Commit 864ab5e8 authored by Julien Malik's avatar Julien Malik
Browse files

ENH: remove MDMD code

parent 7af536a9
No related branches found
No related tags found
No related merge requests found
#include "otbVectorImage.h"
#include "vnl/vnl_matrix.h"
#include "vnl/vnl_cross.h"
#include "vnl/vnl_vector.h"
#include "itkImageRegionConstIterator.h"
#include "itkImageRegionIterator.h"
#include "otbImageFileReader.h"
#include "otbImageFileWriter.h"
typedef double PixelType;
typedef otb::VectorImage<PixelType, 2> VectorImageType;
typedef vnl_matrix<double> VnlMatrixType;
typedef vnl_vector<double> VectorType;
typedef itk::VariableLengthVector<double> ItkVectorType;
typedef otb::ImageFileReader<VectorImageType> ReaderType;
typedef otb::ImageFileWriter<VectorImageType> WriterType;
#define active_debug 0
#define active_debug_loop 0
#define divisor 1
#define divisorParam 10
#define maxIter 100
void AddOneRowOfOnes(const vnl_matrix<double> & m,
vnl_matrix<double> & M)
{
M.set_size(m.rows()+1, m.cols());
for (int i=0; i<M.rows()-1; ++i)
{
M.set_row(i, m.get_row(i));
}
M.set_row(M.rows()-1, 1.0);
}
double Criterion(const VnlMatrixType & X,
const VnlMatrixType & A,
const VnlMatrixType & S,
const double &delt,
const double &lambdS,
const double &lambdD)
{
// This function computes
// f = ||Xsu-Asu.S||_2^2 -
// lambdS * ||1/J*ones - S||_2^2 +
// lambdD * (trace(transpose(A)*A)-1/L*trace(transpose(A).ones.A));
//
// = ||E1||_2^2 -
// lambdS * ||E2||_2^2 +
// lambdD * (trace(A'*A)-1/L*trace(E3);
//
// where || ||_2 is the L2 frobenius_norm,
// J is the number of endmembers and
// L is the number of spectral bands
unsigned int nbEndmembers = A.cols();
unsigned nbBands = A.rows();
VnlMatrixType Xsu, Asu, ones, E1, E2; //, E3;
double evalf, sumColsOfA, trE3;
Xsu.set_size(X.rows()+1, X.cols());
Asu.set_size(A.rows()+1, A.cols());
AddOneRowOfOnes(A, Asu);
AddOneRowOfOnes(X, Xsu);
//------- Computing the function blocs E1, E2 and E3 -------//
// Bloc 1
E1 = Xsu - Asu*S;
// Bloc 2
E2 = S - 1./ ((double) nbEndmembers);// * ones - S;
// Computing trace(transpose(A)*A)
double trAtA = 0;
for (int i=0; i<A.columns(); ++i)
{
trAtA += A.get_column(i).two_norm() * A.get_column(i).two_norm();
}
// Bloc 3: computing fast trE3 = trace(transpose(A)*ones*A)
trE3 = 0;
for (int j=0; j<nbEndmembers; ++j)
{
sumColsOfA = A.get_column(j).sum();
trE3 += sumColsOfA * sumColsOfA;
}
/* for (int j=0; j<nbEndmembers; ++j)
{
sumColsOfA(j) = A.get_column(j).sum();
}
E3.set_size(nbEndmembers, nbEndmembers);
for (int j1=0; j1<nbEndmembers; ++j1)
{
for (int j2=0; j2<nbEndmembers; ++j2)
{
E3(j1,j2) = sumColsOfA(j1)*sumColsOfA(j2);
}
}
*/
//-------------------- Computing f --------------------------//
evalf = E1.frobenius_norm() * E1.frobenius_norm()
- lambdS * E2.frobenius_norm() * E2.frobenius_norm()
+ lambdD * (trAtA - (1./ ((double) nbBands) * trE3));
return evalf;
}
class F
{
public:
virtual double Call(const VnlMatrixType & variMat,
const VnlMatrixType & fixedMat,
const VnlMatrixType & X,
const double & delt,
const double & lambdS,
const double & lambdD) = 0;
};
class FA: public F
{
public:
double Call(const VnlMatrixType & variMat,
const VnlMatrixType & fixedMat,
const VnlMatrixType & X,
const double & delt,
const double & lambdS,
const double & lambdD)
{
double evalf;
evalf = Criterion(X, variMat, fixedMat, delt, lambdS, lambdD);
return evalf;
}
};
class FS: public F
{
public:
double Call(const VnlMatrixType & variMat,
const VnlMatrixType & fixedMat,
const VnlMatrixType & X,
const double & delt,
const double & lambdS,
const double & lambdD)
{
double evalf;
evalf = Criterion(X, fixedMat, variMat, delt, lambdS, lambdD);
return evalf;
}
};
void evalGradS(const VnlMatrixType &X,
const VnlMatrixType &A,
const VnlMatrixType &S,
const double &delt,
const double &lambdS,
VnlMatrixType & gradS)
{
// Calculus of: gradS = 2 * Asu' * (Asu*S-Xsu) - lambd * 2 * (S - 1/J*ones(J,I));
VnlMatrixType Xsu, Asu, ones;
Xsu.set_size(X.rows()+1, X.cols());
Asu.set_size(A.rows()+1, A.cols());
ones.set_size(S.rows(), S.cols());
ones.fill(1.);
AddOneRowOfOnes(A, Asu);
AddOneRowOfOnes(X, Xsu);
gradS = 2. * Asu.transpose() * (Asu*S-Xsu) - lambdS * 2. * (S - (1./((double) S.rows()))); //))*ones);
}
void evalGradA(const VnlMatrixType &X,
const VnlMatrixType &A,
const VnlMatrixType &S,
const double &delt,
const double &lambdD,
VnlMatrixType &gradA)
{
// Compute gradA
// = (A*S-X) * (transpose(S)) + lambdD*(A-1/nbBands*ones(L,L)*A)
// = (A*S-X) * (transpose(S)) + lambdD*A- lambdD*/nbBands*ones(L,L)*A)
VnlMatrixType onesA;
VectorType sumColulmnsOfA;
sumColulmnsOfA.set_size(A.cols());
unsigned int nbBands = A.rows();
// Computing vector onesA
for (int j=0; j<onesA.size(); ++j)
{
sumColulmnsOfA(j) = A.get_column(j).sum();
}
// Initialize gradA
gradA = (A*S-X) * S.transpose();
// 1st update of gradA
gradA += A*lambdD;
// 2nd and last update id gradA, row by row (for performance reasons)
for (int i=0; i<nbBands; ++i)
{
gradA.set_row(i, gradA.get_row(i) - ( sumColulmnsOfA*lambdD/((double)nbBands)));
}
}
void SetNegativeCoefficientsToZero(VnlMatrixType & M)
{
for (int i=0; i<M.rows(); ++i)
{
for (int j=0; j<M.cols(); ++j)
{
if (M(i,j)<0)
M(i,j) = 0;
}
}
}
VnlMatrixType TermByTermMatrixProduct(const VnlMatrixType & M1, const VnlMatrixType & M2)
{
VnlMatrixType M;
M.set_size(M1.rows(), M1.cols());
for (int i=0; i<M.rows(); ++i)
{
for (int j=0; j<M.cols(); ++j)
{
M(i,j) = M1(i,j) * M2(i,j);
}
}
return M;
}
double SumMatrixElements(const VnlMatrixType & M)
{
double sum = 0;
for (int i = 0; i<M.cols(); ++i)
{
sum += M.get_column(i).sum();
}
return sum;
}
bool ArmijoTest(const double & sig,
const VnlMatrixType variMat,
const VnlMatrixType & newVariMat,
const double & evalf,
const double & newEvalf,
const VnlMatrixType & gradVariMat,
const double & alph)
{
bool bit;
//const unsigned int I = variMat.rows();
//const unsigned int J = variMat.cols();
const VnlMatrixType prod = TermByTermMatrixProduct(gradVariMat, newVariMat-variMat);
double sumProd = SumMatrixElements(prod);
if (newEvalf-evalf <= sig*alph*sumProd)
bit = true;
else
bit = false;
return bit;
}
void huckProjGradOneStep(
const VnlMatrixType & X,
const VnlMatrixType & fixedMat,
const VnlMatrixType & gradVariMat,
F & f,
const double & sig,
const double & betinit,
const double & delt,
const double & lambdS,
const double & lambdD,
VnlMatrixType & variMat,
double & alph)
{
double evalf, newEvalf, bet;
evalf = f.Call(variMat, fixedMat, X, delt, lambdS, lambdD); // compute evalf
VnlMatrixType newVariMat = variMat - alph*gradVariMat;
SetNegativeCoefficientsToZero(newVariMat);
newEvalf = f.Call(newVariMat, fixedMat, X, delt, lambdS, lambdD); // compute newEvalf
bool bit = ArmijoTest(sig, variMat, newVariMat, evalf, newEvalf, gradVariMat, alph);
int count = 1;
if (bit == true)
{
while (bit == true)
{
bet = pow(betinit, count);
alph = alph/bet;
newVariMat = variMat - alph*gradVariMat;
SetNegativeCoefficientsToZero(newVariMat);
newEvalf = f.Call(newVariMat, fixedMat, X, delt, lambdS, lambdD);
bit = ArmijoTest(sig, variMat, newVariMat, evalf, newEvalf, gradVariMat, alph);
++count;
}
alph = alph*bet;
newVariMat = variMat - alph*gradVariMat;
SetNegativeCoefficientsToZero(newVariMat);
}
else
{
while (bit == false)
{
bet = pow(betinit, count);
alph = alph*bet;
newVariMat = variMat - alph*gradVariMat;
SetNegativeCoefficientsToZero(newVariMat);
newEvalf = f.Call(newVariMat, fixedMat, X, delt, lambdS, lambdD);
bit = ArmijoTest(sig, variMat, newVariMat, evalf, newEvalf, gradVariMat, alph);
++count;
}
}
variMat = newVariMat;
}
//-------------------------------------------------------------------//
//--------------------------- MAIN ------------------------------//
//-------------------------------------------------------------------//
int main(int argc, char * argv[])
{
typedef itk::ImageRegionConstIterator<VectorImageType> ConstIteratorType;
typedef itk::ImageRegionIterator<VectorImageType> IteratorType;
//--------------------- Input parameters --------------------//
// Reading input arguments
const char * inputFilename = argv[1];
const char * outputFilename = argv[2];
const int nbEndmembers = atoi(argv[3]);
if (active_debug)
std::cout << "Loading params: OK" << std::endl;
// Stop value for the criterion
const double critStopValue = 0.005;
// Tunning the optimized function parameters
const double delt = 1.;
const double lambdD = 0.01;
const double lambdS = 0.01 * (double) nbEndmembers;
// Tunning the projected gradient parameters
double sig = 0.05;
double bet = 0.99;
double alphA = 1.;
double alphS = 1.;
// Other declarations
double critA, critS, crit = 1;
int i; // a counter
if (active_debug)
std::cout << "Loading all parameters: OK" << std::endl;
// Definition of the region to unmix
VectorImageType::RegionType::IndexType index;
index[0] = 0;
index[1] = 0;
VectorImageType::RegionType::SizeType size;
size[0] = 60;
size[1] = 75;
VectorImageType::RegionType region;
region.SetIndex(index);
region.SetSize(size);
if (active_debug)
std::cout << "Building region: OK" << std::endl;
//----- Setting the hsi spectral pixel into matrix columns ----//
ReaderType::Pointer reader = ReaderType::New();
reader->SetFileName(inputFilename);
reader->GetOutput()->SetRegions(region);
reader->GetOutput();
reader->Update();
ConstIteratorType inputIt(reader->GetOutput(), region);
VnlMatrixType X;
X.set_size(reader->GetOutput()->GetNumberOfComponentsPerPixel(), region.GetNumberOfPixels());
i = 0;
for (inputIt.GoToBegin(); !inputIt.IsAtEnd(); ++inputIt)
{
for (int j = 0; j<reader->GetOutput()->GetNumberOfComponentsPerPixel(); ++j)
{
X(j,i) = inputIt.Get().GetElement(j); //inputIt.Get()[j];
}
++i;
}
//------- A and S declaration and initialization --------//
//---> These values fit with the ones chosen in the matlab "nmf_validationOtb.m" function to validate the OTB code.
VnlMatrixType A;
A.set_size(reader->GetOutput()->GetNumberOfComponentsPerPixel(), nbEndmembers);
A.fill(100.);
A.set_column(1,200.);
A.set_column(2,300.);
A.set_column(3,400.);
A.set_column(4,500.);
VnlMatrixType S;
S.set_size(nbEndmembers, region.GetNumberOfPixels());
S.fill(1.);
//----------- Declaration of useful variables -----------//
VnlMatrixType Sold, Sdiff;
Sold.set_size(nbEndmembers, region.GetNumberOfPixels());
Sdiff.set_size(nbEndmembers, region.GetNumberOfPixels());
VnlMatrixType Aold, Adiff;
Aold.set_size(reader->GetOutput()->GetNumberOfComponentsPerPixel(), nbEndmembers);
Adiff.set_size(reader->GetOutput()->GetNumberOfComponentsPerPixel(), nbEndmembers);
VnlMatrixType gradS;
gradS.set_size(nbEndmembers, region.GetNumberOfPixels());
VnlMatrixType gradA;
gradA.set_size(reader->GetOutput()->GetNumberOfComponentsPerPixel(), nbEndmembers);
//----------------- Optimization loop -------------------//
FA fA;
FS fS;
std::cout << "normX = " << X.fro_norm() << std::endl;
unsigned int counter = 0;
while ((crit > critStopValue) && (counter < maxIter))
{
//---------------- Update S -----------------//
Sold = S;
evalGradS(X, A, S, delt, lambdS, gradS); // Compute gradS
if (counter%divisorParam == 0)
{
std::cout << "Iteration = " << counter << std::endl;
std::cout << "Criterion = " << Criterion(X, A, S, delt, lambdS, lambdD) << std::endl;
std::cout << "statGradS = " << gradS.fro_norm() << std::endl;
std::cout << "gradS(0,0) = " << gradS(0,0) << std::endl;
std::cout << "alphS = " << alphS << std::endl;
std::cout << "normS = " << S.fro_norm() << std::endl;
std::cout << "S(0,0) = " << S(0,0) << std::endl;
}
huckProjGradOneStep(X, A, gradS, fS, sig, bet, delt,lambdS, lambdD, S, alphS);
if (counter%divisorParam == 0)
{
std::cout << "alphS = " << alphS << std::endl;
std::cout << "normS = " << S.fro_norm() << std::endl;
std::cout << "S(0,0) = " << S(0,0) << std::endl;
}
//---------------- Update A -----------------//
Aold = A;
evalGradA(X, A, S, delt, lambdD, gradA); // Compute gradS
if (counter%divisorParam == 0)
{
std::cout << "gradA(0,0) = " << gradA(0,0) << std::endl;
}
if (counter%divisorParam == 0)
{
std::cout << "alphA = " << alphA << std::endl;
std::cout << "normA = " << A.fro_norm() << std::endl;
std::cout << "A(0,0) = " << A(0,0) << std::endl;
}
huckProjGradOneStep(X, S, gradA, fA, sig, bet, delt, lambdS, lambdD, A, alphA);
if (counter%divisorParam == 0)
{
std::cout << "alphA = " << alphA << std::endl;
std::cout << "normA = " << A.fro_norm() << std::endl;
std::cout << "A(0,0) = " << A(0,0) << std::endl;
}
//------------ crit evaluation --------------//
Adiff = Aold-A;
critA = Adiff.absolute_value_max();
Sdiff = Sold - S;
critS = Sdiff.absolute_value_max();
crit = max(critA,critS);
if (counter%divisorParam == 0)
{
std::cout << "critA value: " << critA << std::endl;
std::cout << "critS value: " << critS << std::endl;
std::cout << "crit value: " << crit << std::endl;
std::cout << "criterion value: " << Criterion(X, A, S, delt, lambdS, lambdD) << std::endl;
std::cout << std::endl;
}
++counter;
}
//--- Putting the rows of in the bands of the output vector image ---//
//---> Could be impoved choosing an imageList for the abundance maps and a vector list for the endmember spectra (columns of A).
VectorImageType::RegionType::IndexType outputindex;
outputindex[0] = 0;
outputindex[1] = 0;
VectorImageType::RegionType outputRegion;
outputRegion.SetSize(size);
outputRegion.SetIndex(outputindex);
VectorImageType::Pointer outputImage = VectorImageType::New();
outputImage->SetNumberOfComponentsPerPixel(nbEndmembers);
outputImage->SetRegions(outputRegion);
outputImage->Allocate();
IteratorType outputIt(outputImage, outputRegion);
VectorImageType::PixelType vectorPixel;
vectorPixel.SetSize(outputImage->GetNumberOfComponentsPerPixel());
i = 0;
for (outputIt.GoToBegin(); !outputIt.IsAtEnd(); ++outputIt)
{
for (int j=0; j<nbEndmembers; ++j)
{
vectorPixel.SetElement(j, S(j,i));
}
outputIt.Set(vectorPixel);
++i;
}
//------ Write the abundance map vector output image -------//
WriterType::Pointer writer = WriterType::New();
writer->SetFileName(outputFilename);
writer->SetInput(outputImage);
writer->Update();
std::cout << "Number of iterations: " << counter << std::endl;
std::cout << A << std::endl;
return EXIT_SUCCESS;
}
En pièce jointe, je t'envoie :
· Le .cxx à faire évoluer en filtre. Il reste encore à choisir un
filtre duquel hériter, et cela dépend:
o De l’input du filtre mdmd qui doit être compatible avec l’output du
filtre VCA ; à mon avis, le mieux est de partir sur un filtre
composite qui initialise par défaut avec VCA, et peut optionnellement
initialiser différemment si d’autre algorithmes de démixage sont
implémentés pas la suite ;
o De l’output : en effet, deux informations sont pertinentes : la
matrices des spectres et la matrice des abondances. D’après moi, le
plus intuitif est probablement de leur associer un vector list et un
image list, respectivement. Reste à trancher sur ce que doit retourner
la méthode GetOuput()…
· Des éléments sur les tests de validation. En effet, j’ai observé que
les sorties de mdmd_matlab et mdmd_otb ne sont pas
identiques. Dans les deux fichiers .txt (en pj) j’ai comparé
l’évolution de certains paramètres intéressants sur les 100
premières itérations, avec un pas de 10. Note que
l’exécution de mdmd requiert plusieurs milliers voir
dizaines de milliers d’itérations. En gros, il n’y a pas de
bug apparent, ni à la compilation ni à l’exécution, et les
résultats concordent sur les premières itérations et
s’écartent progressivement. Il est fortement possible que
cela soit lié :
o A l’optimisation algorithmique de certains calculs matriciels que
j’ai remplacé par des boucles dans le souci de réduire la complexité
algorithmique ; donc approximations numériques en chaîne ;
o A la non-convexité de la fonctionnelle minimisée : une mauvaise
initialisation peut nous amener à un point ou « zone » stationnaire.
o A une erreur de code. J’en doute, mais cette possibilité n’est pas à
exclure… Il faudra vérifier la performance de démixage lorsque VCA OTB
sera fonctionnel.
En ce qui concerne l’obtention des fichiers .txt en pj :
o Pour la version OTB, ces informations peuvent être observées en
jouant avec les valeurs des constantes de pré compilation «
divisorParam » et « maxIter ». Lancer avec la ligne de commande :
./Mdmd cupriteSubHsi.tif mdmdTest.tif 5
o Pour la version matlab, il faut modifier la fonction enc_mdmd (ligne
1.17) : commenter l’appel de nmf et décommenter celui de
nmf_validationOtb. Dans la fonction nmf_validationOtb, tu as aussi la
possibilité de jouer avec les paramètres maxIter et divisorParam. Il
faut ensuite lancer le script buildToolsForMdmdTest.m. Je t’enverrai
l’ensemble des codes matlab dans un prochain mail.
normX = 603654.214
Iteration = 0
Criterion = 90094702329
statGradS = 1285321873.221
gradS(1,1) = -1706192.08
alphS = 1
normS = 150
S(1,1) = 1
alphS = 1.0809e-007
normS = 281.4423
S(1,1) = 1.1844
gradA(1,1) = 9957501.9272
alphA = 1
normA = 2345.2079
A(1,1) = 100
alphA = 3.0369e-005
normA = 2153.7625
A(1,1) = 0
critA = 500
critS = 2.3101
crit = 500
Criterion = 15309260462.5938
Iteration = 10
Criterion = 2686889068.8731
statGradS = 783553420.4516
gradS(1,1) = -1960100.0349
alphS = 1.4925e-008
normS = 73.6853
S(1,1) = 0.47935
alphS = 1.2836e-008
normS = 68.5848
S(1,1) = 0.50451
gradA(1,1) = 167401.0128
alphA = 0.00036352
normA = 8648.1097
A(1,1) = 501.0179
alphA = 0.00042267
normA = 9204.2799
A(1,1) = 430.2618
critA = 159.619
critS = 0.26956
crit = 159.619
Criterion = 2337128817.8505
Iteration = 20
Criterion = 1595837997.2004
statGradS = 514347986.8684
gradS(1,1) = -2395985.1198
alphS = 1.082e-008
normS = 64.5575
S(1,1) = 0.43852
alphS = 1.1151e-008
normS = 66.3822
S(1,1) = 0.46524
gradA(1,1) = 72680.8909
alphA = 0.00047208
normA = 9356.9452
A(1,1) = 493.7102
alphA = 0.00044446
normA = 9098.6778
A(1,1) = 461.4067
critA = 150.2003
critS = 0.23721
crit = 150.2003
Criterion = 1462527572.6131
Iteration = 30
Criterion = 1053584906.9984
statGradS = 175434152.3404
gradS(1,1) = 2707489.6621
alphS = 1.5228e-008
normS = 74.8763
S(1,1) = 0.46052
alphS = 1.5075e-008
normS = 74.1781
S(1,1) = 0.4197
gradA(1,1) = -42084.0382
alphA = 0.00035629
normA = 8218.0184
A(1,1) = 383.2538
alphA = 0.00035989
normA = 8306.6489
A(1,1) = 398.3994
critA = 173.1547
critS = 0.20005
crit = 173.1547
Criterion = 1012161036.8095
Iteration = 40
Criterion = 655198203.4228
statGradS = 298221005.5584
gradS(1,1) = 8029255.2735
alphS = 1.0187e-008
normS = 61.4219
S(1,1) = 0.34414
alphS = 9.5907e-009
normS = 60.0239
S(1,1) = 0.26714
gradA(1,1) = -29385.3752
alphA = 0.00053259
normA = 10142.5983
A(1,1) = 458.078
alphA = 0.0005657
normA = 10374.7819
A(1,1) = 474.7012
critA = 144.9396
critS = 0.18088
crit = 144.9396
Criterion = 626331229.6128
Iteration = 50
Criterion = 434848221.4427
statGradS = 368651394.7997
gradS(1,1) = 11447511.8921
alphS = 6.9531e-009
normS = 51.2871
S(1,1) = 0.22499
alphS = 6.7466e-009
normS = 50.6527
S(1,1) = 0.14776
gradA(1,1) = 7364.3605
alphA = 0.0008123
normA = 12203.0817
A(1,1) = 515.2094
alphA = 0.00082051
normA = 12352.9388
A(1,1) = 509.1669
critA = 109.5425
critS = 0.16252
crit = 109.5425
Criterion = 418793776.9749
Iteration = 60
Criterion = 343547786.5383
statGradS = 384313500.5553
gradS(1,1) = 11776869.6192
alphS = 6.0405e-009
normS = 48.4072
S(1,1) = 0.14876
alphS = 6.1015e-009
normS = 48.5162
S(1,1) = 0.076905
gradA(1,1) = 40985.0319
alphA = 0.00090726
normA = 12940.2071
A(1,1) = 499.0921
alphA = 0.00089819
normA = 12894.9461
A(1,1) = 462.28
critA = 70.0152
critS = 0.17264
crit = 70.0152
Criterion = 342185614.9388
Iteration = 70
Criterion = 256484071.5246
statGradS = 306288599.1876
gradS(1,1) = 7590872.3464
alphS = 7.3114e-009
normS = 53.5817
S(1,1) = 0.12257
alphS = 7.5352e-009
normS = 54.4638
S(1,1) = 0.065368
gradA(1,1) = 71501.9484
alphA = 0.00072001
normA = 11625.3676
A(1,1) = 374.5915
alphA = 0.00069863
normA = 11426.3698
A(1,1) = 324.6383
critA = 60.0957
critS = 0.17431
crit = 60.0957
Criterion = 248933011.4906
Iteration = 80
Criterion = 154988443.11
statGradS = 184269367.5074
gradS(1,1) = 3611337.3885
alphS = 1.0499e-008
normS = 63.7444
S(1,1) = 0.12572
alphS = 1.082e-008
normS = 64.8281
S(1,1) = 0.08664
gradA(1,1) = 84088.8308
alphA = 0.00050143
normA = 9784.1624
A(1,1) = 230.7511
alphA = 0.00048653
normA = 9619.8164
A(1,1) = 189.8391
critA = 42.3723
critS = 0.14848
crit = 42.3723
Criterion = 152720365.3
Iteration = 90
Criterion = 98162265.2132
statGradS = 104617796.6329
gradS(1,1) = 1734026.3648
alphS = 1.4193e-008
normS = 74.2287
S(1,1) = 0.14269
alphS = 1.4628e-008
normS = 75.136
S(1,1) = 0.11733
gradA(1,1) = 81157.397
alphA = 0.00036352
normA = 8453.58
A(1,1) = 130.5066
alphA = 0.00035989
normA = 8354.8058
A(1,1) = 101.2989
critA = 30.1449
critS = 0.11667
crit = 30.1449
Criterion = 96507189.1139
Number of iterations: 100
A =
1.0e+003 *
0.0687 0.2304 0.4551 0.6798 0.9044
0.1705 0.4777 0.8509 1.2241 1.5974
0.1820 0.5584 1.0059 1.4535 1.9011
0.3656 0.6996 1.0869 1.4742 1.8614
0.8737 1.0188 1.1714 1.3240 1.4765
1.2724 1.2597 1.2061 1.1525 1.0989
1.3722 1.3273 1.2295 1.1317 1.0339
1.0266 1.0285 0.9930 0.9574 0.9219
1.0197 1.0677 1.0689 1.0701 1.0713
0.8494 0.9260 0.9700 1.0141 1.0581
\ No newline at end of file
ahk@ubuntu:~/OTB/Examples/ProjetCnes/AnomalyDetection$ ./Mdmd cupriteSubHsi.tif mdmdTest.tif 5
normX = 603654
Iteration = 0
Criterion = 9.00947e+10
statGradS = 1.28532e+09
gradS(0,0) = -1.70619e+06
alphS = 1
normS = 150
S(0,0) = 1
alphS = 1.08086e-07
normS = 281.442
S(0,0) = 1.18442
gradA(0,0) = 9.9575e+06
alphA = 1
normA = 2345.21
A(0,0) = 100
alphA = 3.03687e-05
normA = 2153.76
A(0,0) = 0
critA value: 500
critS value: 2.31013
crit value: 500
criterion value: 1.53093e+10
Iteration = 10
Criterion = 2.68652e+09
statGradS = 7.83454e+08
gradS(0,0) = -1.95985e+06
alphS = 1.49246e-08
normS = 73.6866
S(0,0) = 0.479355
alphS = 1.28361e-08
normS = 68.5867
S(0,0) = 0.504512
gradA(0,0) = 167473
alphA = 0.000363524
normA = 8647.92
A(0,0) = 501.015
alphA = 0.000422674
normA = 9203.98
A(0,0) = 430.229
critA value: 159.607
critS value: 0.269527
crit value: 159.607
criterion value: 2.33655e+09
Iteration = 20
Criterion = 1.59508e+09
statGradS = 5.146e+08
gradS(0,0) = -2.39745e+06
alphS = 1.08201e-08
normS = 64.548
S(0,0) = 0.438478
alphS = 1.11513e-08
normS = 66.372
S(0,0) = 0.465212
gradA(0,0) = 72818.6
alphA = 0.000472083
normA = 9358.43
A(0,0) = 493.849
alphA = 0.000444457
normA = 9100.01
A(0,0) = 461.484
critA value: 150.078
critS value: 0.237329
crit value: 150.078
criterion value: 1.46217e+09
Iteration = 30
Criterion = 1.06721e+09
statGradS = 1.73285e+08
gradS(0,0) = 2.69726e+06
alphS = 1.52277e-08
normS = 74.949
S(0,0) = 0.460691
alphS = 1.52277e-08
normS = 74.2361
S(0,0) = 0.419618
gradA(0,0) = -42883.7
alphA = 0.00035629
normA = 8210.26
A(0,0) = 382.699
alphA = 0.000359889
normA = 8301.62
A(0,0) = 398.132
critA value: 175.94
critS value: 0.19938
crit value: 175.94
criterion value: 1.03154e+09
Iteration = 40
Criterion = 7.03605e+08
statGradS = 3.30188e+08
gradS(0,0) = 8.8162e+06
alphS = 9.98416e-09
normS = 60.47
S(0,0) = 0.342045
alphS = 9.39989e-09
normS = 58.9103
S(0,0) = 0.259174
gradA(0,0) = -32395.9
alphA = 0.000548896
normA = 10317.8
A(0,0) = 468.965
alphA = 0.000583014
normA = 10586.9
A(0,0) = 487.853
critA value: 154.074
critS value: 0.192659
crit value: 154.074
criterion value: 6.79176e+08
Iteration = 50
Criterion = 4.98772e+08
statGradS = 4.30776e+08
gradS(0,0) = 1.33197e+07
alphS = 6.54621e-09
normS = 49.9772
S(0,0) = 0.22549
alphS = 6.48075e-09
normS = 49.4899
S(0,0) = 0.139169
gradA(0,0) = 12714.6
alphA = 0.000845623
normA = 12509.6
A(0,0) = 534.193
alphA = 0.000871508
normA = 12626
A(0,0) = 523.113
critA value: 105.851
critS value: 0.190364
crit value: 105.851
criterion value: 4.93914e+08
Iteration = 60
Criterion = 3.80354e+08
statGradS = 3.93095e+08
gradS(0,0) = 1.0988e+07
alphS = 6.54621e-09
normS = 50.6818
S(0,0) = 0.162198
alphS = 6.74659e-09
normS = 51.3901
S(0,0) = 0.0880664
gradA(0,0) = 61111.5
alphA = 0.000812302
normA = 12271.6
A(0,0) = 452.376
alphA = 0.000788176
normA = 12084.9
A(0,0) = 404.209
critA value: 70.2173
critS value: 0.196986
crit value: 70.2173
criterion value: 3.70205e+08
Iteration = 70
Criterion = 2.30588e+08
statGradS = 2.47836e+08
gradS(0,0) = 5.30463e+06
alphS = 9.68762e-09
normS = 61.4907
S(0,0) = 0.148233
alphS = 9.98416e-09
normS = 62.7932
S(0,0) = 0.095271
gradA(0,0) = 90569.3
alphA = 0.000532594
normA = 10102.7
A(0,0) = 270.97
alphA = 0.000516775
normA = 9893.02
A(0,0) = 224.166
critA value: 55.3946
critS value: 0.180005
crit value: 55.3946
criterion value: 2.17339e+08
Iteration = 80
Criterion = 1.34165e+08
statGradS = 1.39846e+08
gradS(0,0) = 2.40736e+06
alphS = 1.43366e-08
normS = 74.8298
S(0,0) = 0.158344
alphS = 1.47754e-08
normS = 76.0876
S(0,0) = 0.122774
gradA(0,0) = 95683.9
alphA = 0.000359889
normA = 8361.74
A(0,0) = 136.831
alphA = 0.0003492
normA = 8229.55
A(0,0) = 103.418
critA value: 33.8398
critS value: 0.152228
crit value: 33.8398
criterion value: 1.2713e+08
Iteration = 90
Criterion = 8.73534e+07
statGradS = 7.57751e+07
gradS(0,0) = 1.06525e+06
alphS = 1.93816e-08
normS = 86.0671
S(0,0) = 0.178105
alphS = 1.95773e-08
normS = 87.0371
S(0,0) = 0.15725
gradA(0,0) = 95819.8
alphA = 0.000271615
normA = 7334.9
A(0,0) = 60.3666
alphA = 0.000263548
normA = 7259.52
A(0,0) = 35.1135
critA value: 27.7495
critS value: 0.111613
crit value: 27.7495
criterion value: 8.41154e+07
Number of iterations: 100
14.7873 173.733 391.879 610.024 828.17
67.4916 368.579 733.588 1098.6 1463.61
60.5014 429.145 866.819 1304.49 1742.17
235.566 561.825 939.885 1317.94 1696
737.104 876.644 1022.73 1168.81 1314.89
1134.32 1117.62 1060.56 1003.51 946.456
1231.9 1183.13 1082.53 981.923 881.319
912.844 911.043 872.836 834.629 796.422
896.43 939.267 937.997 936.726 935.456
737.206 808.711 849.732 890.753 931.774
ahk@ubuntu:~/OTB/Examples/ProjetCnes/AnomalyDetection$
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment