LDA详解 包括opencv实现代码

voidMyLDA(constMat&_data,constvector<int>&labelVec,int&classNum,intdim,Mat&eigenvalues,Mat&eigenvectors)

{

Matdata=cv::Mat_<double>(_data);

cout<<"/********************判定LDA是否正确******************/"<<endl;

cout<<"data.rows:"<<data.rows<<endl;

cout<<"data.cols:"<<data.cols<<endl;

cout<<"labelVec.size:"<<labelVec.size()<<endl;

cout<<"classNum:"<<classNum<<endl;

cout<<"dim:"<<dim<<endl;

cout<<"判定完毕"<<endl;

assert(dim>0);

intD=data.cols;

intN=data.rows;

if(dim>data.cols)

{

dim=data.cols;

}

/***************************计算类内类间距离*************************/

cout<<"计算类内类间距离..."<<endl;

///类内距离

MatSw;

///类间距离

MatSb;

///每类的均值

vectormeanClass(classNum);

///总均值

MatmeanTotal;

///每类的样本数

vector<int>numClass(classNum);

//initialize

for(inti=0;i<classNum;i++)

{

numClass[i]=0;

meanClass[i]=Mat::zeros(1,D,data.type());//!Dx1imagevector

}

meanTotal=Mat::zeros(1,D,data.type());

//calculatesums

for(inti=0;i

{

Matinstance=data.row(i);

intclassIdx=labelVec[i];

add(meanTotal,instance,meanTotal);

add(meanClass[classIdx],instance,meanClass[classIdx]);

numClass[classIdx]++;

}

//calculatemeans

meanTotal.convertTo(meanTotal,meanTotal.type(),1.0/static_cast<double>(data.rows));

for(inti=0;i<classNum;i++)

{

assert(numClass[i]!=0);

meanClass[i].convertTo(meanClass[i],meanClass[i].type(),

1.0/static_cast<double>(numClass[i]));

}

//subtractclassmeans

for(inti=0;i<N;i++)

{

intclassIdx=labelVec[i];

Matinstance=data.row(i);

subtract(instance,meanClass[classIdx],instance);

}

//calculatewithin-classesscatter

Sw=Mat::zeros(D,D,data.type());

//记录每类的类内距离

introwBegin=0;

for(inti=0;i

{

if(i!=0)

{

rowBegin+=numClass[i-1];

}

Rectr(0,rowBegin,D,numClass[i]);

Mattmp(data,r);

MatSi;

mulTransposed(tmp,Si,TRUE,noArray(),(double)numClass[i]/N);

add(Sw,Si,Sw);

}

//calculatebetween-classesscatter

Sb=Mat::zeros(D,D,data.type());

for(inti=0;i<classNum;i++)

{

Mattmp;

subtract(meanClass[i],meanTotal,tmp);

mulTransposed(tmp,tmp,TRUE,noArray(),(double)numClass[i]/N);

add(Sb,tmp,Sb);

}

/***************************原始LDA*************************/

#if0

///计算Sw^-1*Sb

MatM=Sw.inv()*Sb;

///求Sw^-1*Sb的特征值和特征向量

EigenvalueDecompositiones(M);

eigenvalues=es.eigenvalues();

eigenvectors=es.eigenvectors();

#endif

/***************************正则化ＬＤＡ*************************/

cout<<"计算LDA："<<endl;

//eigendecomposeofSw

cout<<"Sw分解..."<<endl;

MateigValSw;

MateigVectSw;

eigen(Sw,eigValSw,eigVectSw);

eigValSw=eigValSw.reshape(1,1);

eigVectSw=eigVectSw.t();///<特征向量按列排

///特征值的求解过程中有一定的误差，特别是较小的值，因此将这些值置为一个固定值

//regularization,changeeigenvaluesthatarelessthangiventhresholdtoconstantvalue

cout<<endl<<"将小的特征值置为定值..."<<endl;

constfloatthreshold=1.0E-6f*eigValSw.at<double>(0,0);

intindex=dim;

for(inti=0;i<dim;++i)

{

if(eigValSw.at<double>(0,i)<=threshold)

{

index=i;

break;

}

}

//eigenvaluesthatarelessthenthresholdissettotrimValue

floattrimValue=threshold;

if(index>0)

{

trimValue=eigValSw.at<double>(0,index-1)*0.8f*index/dim;

}

for(inti=index;i<dim;++i)

{

eigValSw.at<double>(0,i)=trimValue;

}

cout<<endl<<"计算Sw^(-1/2)=R*I^-1/2"<<endl;

//calculateinversesquarerootofeigenvalues

for(inti=0;i<dim;++i)

{

eigValSw.at<double>(0,i)=1.0f/sqrt(eigValSw.at<double>(0,i));

}

for(inti=0;i<dim;++i)

{

for(intj=0;j<dim;++j)

{

eigVectSw.at<double>(i,j)*=eigValSw.at<double>(0,j);

}

}

cout<<endl<<"计算Sw^(-1/2))^T*Sb*(Sw^(-1/2)"<<endl;

//calculate(Sw^(-1/2))^T*Sb*(Sw^(-1/2))

MatS=eigVectSw.t()*Sb*eigVectSw;

cout<<endl<<"eigendecomposeof(Sw^(-1/2))'*Sb*(Sw^(-1/2))"<<endl;

//eigendecomposeof(Sw^(-1/2))'*Sb*(Sw^(-1/2))

EigenvalueDecompositiones(S);

eigenvalues=es.eigenvalues();

eigenvectors=es.eigenvectors();

cout<<endl<<"eigVectSw*eigenvectors"<<endl;

eigenvectors=eigVectSw*eigenvectors.t();

/***************对特征值特征向量进行排序等处理********************/

cout<<"reshape"<<endl;

eigenvalues=eigenvalues.reshape(1,1);

/////排序：根据特征值的大小将特征值和特征向量排序

cout<<"sorting..."<<endl;

vector<int>sorted_indices;

for(intv=0;v

{

sorted_indices.push_back(v);

}

sortIdxEigenVal(eigenvalues,sorted_indices,CV_SORT_EVERY_ROW+CV_SORT_DESCENDING);

eigenvectors=sortMatrixColumnsByIndices(eigenvectors,sorted_indices);

/////取出需要的维度作为结果

cout<<"取出需要的维度作为结果"<<endl;

eigenvalues=Mat(eigenvalues,Range::all(),Range(0,dim));

eigenvectors=Mat(eigenvectors,Range::all(),Range(0,dim));

cout<<"clear"<<endl;

numClass.clear();

meanClass.clear();

cout<<"eigenvectors.rows:"<<eigenvectors.rows<<endl;

cout<<"eigenvectors.cols:"<<eigenvectors.cols<<endl;

cout<<"/********************MyLDAend******************/"<<endl;

}

///////////////////////////////////////////////////////////////////////////////////////////////

///排序

voidexchange(Mat&data,vector<int>&sorted_indices,intx,inty)

{

doubled=data.at<double>(0,x);

data.at<double>(0,x)=data.at<double>(0,y);

data.at<double>(0,y)=d;

inti=sorted_indices[x];

sorted_indices[x]=sorted_indices[y];

sorted_indices[y]=i;

}

intPartition(Mat&data,vector<int>&sorted_indices,intp,intr)

{

doublekey=data.at<double>(0,r);

inti=p-1;//第一个标记：存放小于key的元素的末位置

for(intj=p;j

{

if(data.at<double>(0,j)>=key)

{

i++;

//exchange(data+i,data+j);

exchange(data,sorted_indices,i,j);

}

}

//exchange(data+i+1,data+r);

exchange(data,sorted_indices,i+1,r);

returni+1;

}

voidQuickSort(Mat&data,vector<int>&sorted_indices,intp,intr)//p:起始位置r：末尾位置

{

if(p

{

intq=Partition(data,sorted_indices,p,r);//核心函数

QuickSort(data,sorted_indices,p,q-1);

QuickSort(data,sorted_indices,q+1,r);

}

}

///inputData是一行数据

voidsortIdxEigenVal(Mat&inputData,vector<int>&sorted_indices,intflag)

{

if(flag!=CV_SORT_EVERY_ROW+CV_SORT_DESCENDING)

{

return;

}

for(inti=0;i

{

Matm=inputData.row(i);

QuickSort(inputData.row(i),sorted_indices,0,m.cols-1);

}

}

////////////////////////////////////////////////////////////////////////////////////////////////

MatsortMatrixColumnsByIndices(Matsrc,vector<int>indices)

{

Matdst(src.rows,src.cols,src.type());

for(size_tidx=0;idx<indices.size();idx++)

{

MatoriginalCol=src.col(indices[idx]);

MatsortedCol=dst.col((int)idx);

originalCol.copyTo(sortedCol);

}

returndst;

};

templatestaticbool

isSymmetric_(InputArraysrc){

Mat_src=src.getMat();

if(_src.cols!=_src.rows)

returnfalse;

for(inti=0;i<_src.rows;i++){

for(intj=0;j<_src.cols;j++){

_Tpa=_src.at<_Tp>(i,j);

_Tpb=_src.at<_Tp>(j,i);

if(a!=b){

returnfalse;

}

}

}

returntrue;

}

templatestaticbool

isSymmetric_(InputArraysrc,doubleeps){

Mat_src=src.getMat();

if(_src.cols!=_src.rows)

returnfalse;

for(inti=0;i<_src.rows;i++){

for(intj=0;j<_src.cols;j++){

_Tpa=_src.at<_Tp>(i,j);

_Tpb=_src.at<_Tp>(j,i);

if(std::abs(a-b)>eps){

returnfalse;

}

}

}

returntrue;

}

staticboolisSymmetric(InputArraysrc,doubleeps=1e-16)

{

Matm=src.getMat();

switch(m.type()){

caseCV_8SC1:returnisSymmetric_<char>(m);break;

caseCV_8UC1:

returnisSymmetric_char>(m);break;

caseCV_16SC1:

returnisSymmetric_<short>(m);break;

caseCV_16UC1:

returnisSymmetric_short>(m);break;

caseCV_32SC1:

returnisSymmetric_<int>(m);break;

caseCV_32FC1:

returnisSymmetric_<float>(m,eps);break;

caseCV_64FC1:

returnisSymmetric_<double>(m,eps);break;

default:

break;

}

returnfalse;

}

/////////////////////////////////////////////////////////////////////////////////////////////////////

classEigenvalueDecomposition{

private:

//Holdsthedatadimension.

intn;

//Storesreal/imagpartofacomplexdivision.

doublecdivr,cdivi;

//Pointertointernalmemory.

double*d,*e,*ort;

double**V,**H;

//Holdsthecomputedeigenvalues.

Mat_eigenvalues;

//Holdsthecomputedeigenvectors.

Mat_eigenvectors;

//Allocatesmemory.

template

_Tp*alloc_1d(intm){

returnnew_Tp[m];

}

//Allocatesmemory.

template

_Tp*alloc_1d(intm,_Tpval){

_Tp*arr=alloc_1d<_Tp>(m);

for(inti=0;i<m;i++)

arr[i]=val;

returnarr;

}

//Allocatesmemory.

template

_Tp**alloc_2d(intm,intn){

_Tp**arr=new_Tp*[m];

for(inti=0;i<m;i++)

arr[i]=new_Tp
;

returnarr;

}

//Allocatesmemory.

template

_Tp**alloc_2d(intm,intn,_Tpval){

_Tp**arr=alloc_2d<_Tp>(m,n);

for(inti=0;i<m;i++){

for(intj=0;j<n;j++){

arr[i][j]=val;

}

}

returnarr;

}

voidcdiv(doublexr,doublexi,doubleyr,doubleyi){

doubler,d;

if(std::abs(yr)>std::abs(yi)){

r=yi/yr;

d=yr+r*yi;

cdivr=(xr+r*xi)/d;

cdivi=(xi-r*xr)/d;

}else{

r=yr/yi;

d=yi+r*yr;

cdivr=(r*xr+xi)/d;

cdivi=(r*xi-xr)/d;

}

}

//NonsymmetricreductionfromHessenbergtorealSchurform.

voidhqr2(){

//ThisisderivedfromtheAlgolprocedurehqr2,

//byMartinandWilkinson,HandbookforAuto.Comp.,

//Vol.ii-LinearAlgebra,andthecorresponding

//FortransubroutineinEISPACK.

//Initialize

intnn=this->n;

intn=nn-1;

intlow=0;

inthigh=nn-1;

doubleeps=pow(2.0,-52.0);

doubleexshift=0.0;

doublep=0,q=0,r=0,s=0,z=0,t,w,x,y;

//Storerootsisolatedbybalancandcomputematrixnorm

doublenorm=0.0;

for(inti=0;i<nn;i++){

if(i<low||i>high){

d[i]=H[i][i];

e[i]=0.0;

}

for(intj=max(i-1,0);j<nn;j++){

norm=norm+std::abs(H[i][j]);

}

}

//Outerloopovereigenvalueindex

intiter=0;

while(n>=low){

//Lookforsinglesmallsub-diagonalelement

intl=n;

while(l>low){

s=std::abs(H[l-1][l-1])+std::abs(H[l][l]);

if(s==0.0){

s=norm;

}

if(std::abs(H[l][l-1])<eps*s){

break;

}

l--;

}

//Checkforconvergence

//Onerootfound

if(l==n){

H

=H

+exshift;

d
=H

;

e
=0.0;

n--;

iter=0;

//Tworootsfound

}elseif(l==n-1){

w=H
[n-1]*H[n-1]
;

p=(H[n-1][n-1]-H

)/2.0;

q=p*p+w;

z=sqrt(std::abs(q));

H

=H

+exshift;

H[n-1][n-1]=H[n-1][n-1]+exshift;

x=H

;

//Realpair

if(q>=0){

if(p>=0){

z=p+z;

}else{

z=p-z;

}

d[n-1]=x+z;

d
=d[n-1];

if(z!=0.0){

d
=x-w/z;

}

e[n-1]=0.0;

e
=0.0;

x=H
[n-1];

s=std::abs(x)+std::abs(z);

p=x/s;

q=z/s;

r=sqrt(p*p+q*q);

p=p/r;

q=q/r;

//Rowmodification

for(intj=n-1;j<nn;j++){

z=H[n-1][j];

H[n-1][j]=q*z+p*H
[j];

H
[j]=q*H
[j]-p*z;

}

//Columnmodification

for(inti=0;i<=n;i++){

z=H[i][n-1];

H[i][n-1]=q*z+p*H[i]
;

H[i]
=q*H[i]
-p*z;

}

//Accumulatetransformations

for(inti=low;i<=high;i++){

z=V[i][n-1];

V[i][n-1]=q*z+p*V[i]
;

V[i]
=q*V[i]
-p*z;

}

//Complexpair

}else{

d[n-1]=x+p;

d
=x+p;

e[n-1]=z;

e
=-z;

}

n=n-2;

iter=0;

//Noconvergenceyet

}else{

//Formshift

x=H

;

y=0.0;

w=0.0;

if(l<n){

y=H[n-1][n-1];

w=H
[n-1]*H[n-1]
;

}

//Wilkinson'soriginaladhocshift

if(iter==10){

exshift+=x;

for(inti=low;i<=n;i++){

H[i][i]-=x;

}

s=std::abs(H
[n-1])+std::abs(H[n-1][n-2]);

x=y=0.75*s;

w=-0.4375*s*s;

}

//MATLAB'snewadhocshift

if(iter==30){

s=(y-x)/2.0;

s=s*s+w;

if(s>0){

s=sqrt(s);

if(y<x){

s=-s;

}

s=x-w/((y-x)/2.0+s);

for(inti=low;i<=n;i++){

H[i][i]-=s;

}

exshift+=s;

x=y=w=0.964;

}

}

iter=iter+1;//(Couldcheckiterationcounthere.)

//Lookfortwoconsecutivesmallsub-diagonalelements

intm=n-2;

while(m>=l){

z=H[m][m];

r=x-z;

s=y-z;

p=(r*s-w)/H[m+1][m]+H[m][m+1];

q=H[m+1][m+1]-z-r-s;

r=H[m+2][m+1];

s=std::abs(p)+std::abs(q)+std::abs(r);

p=p/s;

q=q/s;

r=r/s;

if(m==l){

break;

}

if(std::abs(H[m][m-1])*(std::abs(q)+std::abs(r))<eps*(std::abs(p)

*(std::abs(H[m-1][m-1])+std::abs(z)+std::abs(

H[m+1][m+1])))){

break;

}

m--;

}

for(inti=m+2;i<=n;i++){

H[i][i-2]=0.0;

if(i>m+2){

H[i][i-3]=0.0;

}

}

//DoubleQRstepinvolvingrowsl:nandcolumnsm:n

for(intk=m;k<=n-1;k++){

boolnotlast=(k!=n-1);

if(k!=m){

p=H[k][k-1];

q=H[k+1][k-1];

r=(notlast?H[k+2][k-1]:0.0);

x=std::abs(p)+std::abs(q)+std::abs(r);

if(x!=0.0){

p=p/x;

q=q/x;

r=r/x;

}

}

if(x==0.0){

break;

}

s=sqrt(p*p+q*q+r*r);

if(p<0){

s=-s;

}

if(s!=0){

if(k!=m){

H[k][k-1]=-s*x;

}elseif(l!=m){

H[k][k-1]=-H[k][k-1];

}

p=p+s;

x=p/s;

y=q/s;

z=r/s;

q=q/p;

r=r/p;

//Rowmodification

for(intj=k;j<nn;j++){

p=H[k][j]+q*H[k+1][j];

if(notlast){

p=p+r*H[k+2][j];

H[k+2][j]=H[k+2][j]-p*z;

}

H[k][j]=H[k][j]-p*x;

H[k+1][j]=H[k+1][j]-p*y;

}

//Columnmodification

for(inti=0;i<=min(n,k+3);i++){

p=x*H[i][k]+y*H[i][k+1];

if(notlast){

p=p+z*H[i][k+2];

H[i][k+2]=H[i][k+2]-p*r;

}

H[i][k]=H[i][k]-p;

H[i][k+1]=H[i][k+1]-p*q;

}

//Accumulatetransformations

for(inti=low;i<=high;i++){

p=x*V[i][k]+y*V[i][k+1];

if(notlast){

p=p+z*V[i][k+2];

V[i][k+2]=V[i][k+2]-p*r;

}

V[i][k]=V[i][k]-p;

V[i][k+1]=V[i][k+1]-p*q;

}

}//(s!=0)

}//kloop

}//checkconvergence

}//while(n>=low)

//Backsubstitutetofindvectorsofuppertriangularform

if(norm==0.0){

return;

}

for(n=nn-1;n>=0;n--){

p=d
;

q=e
;

//Realvector

if(q==0){

intl=n;

H

=1.0;

for(inti=n-1;i>=0;i--){

w=H[i][i]-p;

r=0.0;

for(intj=l;j<=n;j++){

r=r+H[i][j]*H[j]
;

}

if(e[i]<0.0){

z=w;

s=r;

}else{

l=i;

if(e[i]==0.0){

if(w!=0.0){

H[i]
=-r/w;

}else{

H[i]
=-r/(eps*norm);

}

//Solverealequations

}else{

x=H[i][i+1];

y=H[i+1][i];

q=(d[i]-p)*(d[i]-p)+e[i]*e[i];

t=(x*s-z*r)/q;

H[i]
=t;

if(std::abs(x)>std::abs(z)){

H[i+1]
=(-r-w*t)/x;

}else{

H[i+1]
=(-s-y*t)/z;

}

}

//Overflowcontrol

t=std::abs(H[i]
);

if((eps*t)*t>1){

for(intj=i;j<=n;j++){

H[j]
=H[j]
/t;

}

}

}

}

//Complexvector

}elseif(q<0){

intl=n-1;

//Lastvectorcomponentimaginarysomatrixistriangular

if(std::abs(H
[n-1])>std::abs(H[n-1]
)){

H[n-1][n-1]=q/H
[n-1];

H[n-1]
=-(H

-p)/H
[n-1];

}else{

cdiv(0.0,-H[n-1]
,H[n-1][n-1]-p,q);

H[n-1][n-1]=cdivr;

H[n-1]
=cdivi;

}

H
[n-1]=0.0;

H

=1.0;

for(inti=n-2;i>=0;i--){

doublera,sa,vr,vi;

ra=0.0;

sa=0.0;

for(intj=l;j<=n;j++){

ra=ra+H[i][j]*H[j][n-1];

sa=sa+H[i][j]*H[j]
;

}

w=H[i][i]-p;

if(e[i]<0.0){

z=w;

r=ra;

s=sa;

}else{

l=i;

if(e[i]==0){

cdiv(-ra,-sa,w,q);

H[i][n-1]=cdivr;

H[i]
=cdivi;

}else{

//Solvecomplexequations

x=H[i][i+1];

y=H[i+1][i];

vr=(d[i]-p)*(d[i]-p)+e[i]*e[i]-q*q;

vi=(d[i]-p)*2.0*q;

if(vr==0.0&&vi==0.0){

vr=eps*norm*(std::abs(w)+std::abs(q)+std::abs(x)

+std::abs(y)+std::abs(z));

}

cdiv(x*r-z*ra+q*sa,

x*s-z*sa-q*ra,vr,vi);

H[i][n-1]=cdivr;

H[i]
=cdivi;

if(std::abs(x)>(std::abs(z)+std::abs(q))){

H[i+1][n-1]=(-ra-w*H[i][n-1]+q

*H[i]
)/x;

H[i+1]
=(-sa-w*H[i]
-q*H[i][n

-1])/x;

}else{

cdiv(-r-y*H[i][n-1],-s-y*H[i]
,z,

q);

H[i+1][n-1]=cdivr;

H[i+1]
=cdivi;

}

}

//Overflowcontrol

t=max(std::abs(H[i][n-1]),std::abs(H[i]
));

if((eps*t)*t>1){

for(intj=i;j<=n;j++){

H[j][n-1]=H[j][n-1]/t;

H[j]
=H[j]
/t;

}

}

}

}

}

}

//Vectorsofisolatedroots

for(inti=0;i<nn;i++){

if(i<low||i>high){

for(intj=i;j<nn;j++){

V[i][j]=H[i][j];

}

}

}

//Backtransformationtogeteigenvectorsoforiginalmatrix

for(intj=nn-1;j>=low;j--){

for(inti=low;i<=high;i++){

z=0.0;

for(intk=low;k<=min(j,high);k++){

z=z+V[i][k]*H[k][j];

}

V[i][j]=z;

}

}

}

//NonsymmetricreductiontoHessenbergform.

voidorthes(){

//ThisisderivedfromtheAlgolproceduresorthesandortran,

//byMartinandWilkinson,HandbookforAuto.Comp.,

//Vol.ii-LinearAlgebra,andthecorresponding

//FortransubroutinesinEISPACK.

intlow=0;

inthigh=n-1;

for(intm=low+1;m<=high-1;m++){

//Scalecolumn.

doublescale=0.0;

for(inti=m;i<=high;i++){

scale=scale+std::abs(H[i][m-1]);

}

if(scale!=0.0){

//ComputeHouseholdertransformation.

doubleh=0.0;

for(inti=high;i>=m;i--){

ort[i]=H[i][m-1]/scale;

h+=ort[i]*ort[i];

}

doubleg=sqrt(h);

if(ort[m]>0){

g=-g;

}

h=h-ort[m]*g;

ort[m]=ort[m]-g;

//ApplyHouseholdersimilaritytransformation

//H=(I-u*u'/h)*H*(I-u*u')/h)

for(intj=m;j<n;j++){

doublef=0.0;

for(inti=high;i>=m;i--){

f+=ort[i]*H[i][j];

}

f=f/h;

for(inti=m;i<=high;i++){

H[i][j]-=f*ort[i];

}

}

for(inti=0;i<=high;i++){

doublef=0.0;

for(intj=high;j>=m;j--){

f+=ort[j]*H[i][j];

}

f=f/h;

for(intj=m;j<=high;j++){

H[i][j]-=f*ort[j];

}

}

ort[m]=scale*ort[m];

H[m][m-1]=scale*g;

}

}

//Accumulatetransformations(Algol'sortran).

for(inti=0;i<n;i++){

for(intj=0;j<n;j++){

V[i][j]=(i==j?1.0:0.0);

}

}

for(intm=high-1;m>=low+1;m--){

if(H[m][m-1]!=0.0){

for(inti=m+1;i<=high;i++){

ort[i]=H[i][m-1];

}

for(intj=m;j<=high;j++){

doubleg=0.0;

for(inti=m;i<=high;i++){

g+=ort[i]*V[i][j];

}

//Doubledivisionavoidspossibleunderflow

g=(g/ort[m])/H[m][m-1];

for(inti=m;i<=high;i++){

V[i][j]+=g*ort[i];

}

}

}

}

}

//Releasesallinternalworkingmemory.

voidrelease(){

//releasestheworkingdata

delete[]d;

delete[]e;

delete[]ort;

for(inti=0;i<n;i++){

delete[]H[i];

delete[]V[i];

}

delete[]H;

delete[]V;

}

//ComputestheEigenvalueDecompositionforamatrixgiveninH.

voidcompute(){

//Allocatememoryfortheworkingdata.

V=alloc_2d<double>(n,n,0.0);

d=alloc_1d<double>(n);

e=alloc_1d<double>(n);

ort=alloc_1d<double>(n);

//ReducetoHessenbergform.

orthes();

//ReduceHessenbergtorealSchurform.

hqr2();

//CopyeigenvaluestoOpenCVMatrix.

_eigenvalues.create(1,n,CV_64FC1);

for(inti=0;i<n;i++){

_eigenvalues.at<double>(0,i)=d[i];

}

//CopyeigenvectorstoOpenCVMatrix.

_eigenvectors.create(n,n,CV_64FC1);

for(inti=0;i<n;i++)

for(intj=0;j<n;j++)

_eigenvectors.at<double>(i,j)=V[i][j];

//Deallocatethememorybyreleasingallinternalworkingdata.

release();

}

public: