作者:朱金灿<?xml:namespace prefix = o ns = "urn:schemas-microsoft-com:office:office" />
来源:blog.csdn.net/clever101
主成份分析(Principal Component Analysis,PCA)也叫做主成份变换、主分量分析或 —L(Karhunen—Loeve)变换,是建立在统计特征基础上的多维(如多波段)正交线
性变换。它是遥感图像处理中最常用也是最有用的变换算法之一。
这次我要实现一个主成分分析算法,图是做出来了,但是和著名的遥感软件PCI和ENVI的效果比起来很差。如第一主成分的图如下:
上面噪音极多,而且看起来不合谐。我知道自己的算法有问题,在排除了自己的读取图像的问题后。我考虑到是不是求取特征矩阵时出了问题,因为主成分的输出数据是Y=X*A。
其中X为原图像,Y为目标图像,A为特征向量矩阵。由此我怀疑我的特征矩阵求取有问题。后来从网上找了一种求特征矩阵的办法,进行主成分分析的效果。下面是具体的实现代码:
- //计算特征向量
- /*
- pdblCof[in][out]-----协方差矩阵
- lChannelCount[in]------图像的输入波段数
- pdblVects[out]----特征向量矩阵
- dblEps[in]----误差范围,我取为0.0000001
- Ljt[in]-----循环次数,我取为1000000
- */
- staticintiJcobiMatrixCharacterValue(double**pdblCof,longlChannelCount,std::vector<double>&pdblVects,doubledblEps,longljt)
- {
- longi,j,p,q,u,w,t,s,l;
- doublefm,cn,sn,omega,x,y,d;
- l=1;
- for(i=0;i<lChannelCount;i++)
- {
- pdblVects[i*lChannelCount+i]=1.0;
- for(j=0;j<lChannelCount;j++)
- if(i!=j)pdblVects[i*lChannelCount+j]=0.0;
- }
- while(1){
- fm=0.0;
- for(i=0;i<lChannelCount;i++)
- for(j=0;j<lChannelCount;j++)
- {
- d=fabs(pdblCof[i][j]);
- if((i!=j)&&(d>fm))
- {fm=d;p=i;q=j;}
- }
- if(fm<dblEps)return1;
- if(l>ljt)return0;
- l+=1;
- u=p*lChannelCount+q;w=p*lChannelCount+p;t=q*lChannelCount+p;s=q*lChannelCount+q;
- x=-pdblCof[p][q];
- y=(pdblCof[q][q]-pdblCof[p][p])/2.0;
- omega=x/sqrt(x*x+y*y);
- if(y<0)omega=-omega;
- sn=1.0+sqrt(1.0-omega*omega);
- sn=omega/sqrt(2.0*sn);
- cn=sqrt(1.0-sn*sn);
- fm=pdblCof[p][p];
- pdblCof[p][p]=fm*cn*cn+pdblCof[q][q]*sn*sn+pdblCof[p][q]*omega;
- pdblCof[q][q]=fm*sn*sn+pdblCof[q][q]*cn*cn-pdblCof[p][q]*omega;
- pdblCof[p][q]=0.0;
- pdblCof[q][p]=0.0;
- for(j=0;j<lChannelCount;j++)
- if((j!=p)&&(j!=q))
- {
- fm=pdblCof[p][j];
- pdblCof[p][j]=fm*cn+pdblCof[q][j]*sn;
- pdblCof[q][j]=-fm*sn+pdblCof[q][j]*cn;
- }
- for(i=0;i<lChannelCount;i++)
- if((i!=p)&&(i!=q)){
- fm=pdblCof[i][p];
- pdblCof[i][p]=fm*cn+pdblCof[i][q]*sn;
- pdblCof[i][q]=-fm*sn+pdblCof[i][q]*cn;
- }
- for(i=0;i<lChannelCount;i++)
- {
- fm=pdblVects[i*lChannelCount+p];
- pdblVects[i*lChannelCount+p]=fm*cn+pdblVects[i*lChannelCount+q]*sn;
- pdblVects[i*lChannelCount+q]=-fm*sn+pdblVects[i*lChannelCount+q]*cn;
- }
- }
- return1;
- }
- //根据特征值从大到小排列特征向量矩阵
- /*
- pfMatrix[in][out]-----上一步输出的协方差矩阵
- nBandNum[in]------图像的输入波段数
- pdblVects[out]----上一步输出的特征向量矩阵
- */
- staticvoidSortEigenvector(double**pfMatrix,intnBandNum,std::vector<double>&pfVector)
- {
- longp;
- doublef;
- doubleT;
- intcount=nBandNum;
- for(inti=0;i<count;i++)
- {
- T=pfMatrix[i][i];
- p=i;
- for(intj=i;j<count;j++)
- if(T<pfMatrix[j][j])
- {
- T=pfMatrix[j][j];
- p=j;
- }
- if(p!=i)
- {
- f=pfMatrix[p][p];
- pfMatrix[p][p]=pfMatrix[i][i];
- pfMatrix[i][i]=f;
- for(intj=0;j<count;j++)
- {
- f=pfVector[j*count+p];
- pfVector[j*count+p]=pfVector[j*count+i];
- pfVector[j*count+i]=f;
- }
- }
- }
- }
执行上面两步之后,所得到的特征矩阵为用于和原图像相乘的矩阵A.
对一幅TM图像的第1,2,3,4,5,7通道执行PCA操作,效果图如下:
第一主成分:
第二主成分:
第三主成分:
第四主成分:
第五主成分:
第六主成分:
从上面可以看出,正确的图像看起来视觉非常和谐,毫无刺目的感觉。