主成分分析实现的一个心得

作者：朱金灿<?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为特征向量矩阵。由此我怀疑我的特征矩阵求取有问题。后来从网上找了一种求特征矩阵的办法，进行主成分分析的效果。下面是具体的实现代码：

1. //计算特征向量
2. /*
3. pdblCof[in][out]-----协方差矩阵
4. lChannelCount[in]------图像的输入波段数
5. pdblVects[out]----特征向量矩阵
6. dblEps[in]----误差范围，我取为0.0000001
7. Ljt[in]-----循环次数，我取为1000000
8. */
9. staticintiJcobiMatrixCharacterValue(double**pdblCof,longlChannelCount,std::vector<double>&pdblVects,doubledblEps,longljt)
10. {
11. longi,j,p,q,u,w,t,s,l;
12. doublefm,cn,sn,omega,x,y,d;
13. l=1;
14. for(i=0;i<lChannelCount;i++)
15. {
16. pdblVects[i*lChannelCount+i]=1.0;
17. for(j=0;j<lChannelCount;j++)
18. if(i!=j)pdblVects[i*lChannelCount+j]=0.0;
19. }
21. while(1){
22. fm=0.0;
23. for(i=0;i<lChannelCount;i++)
24. for(j=0;j<lChannelCount;j++)
25. {
26. d=fabs(pdblCof[i][j]);
27. if((i!=j)&&(d>fm))
28. {fm=d;p=i;q=j;}
29. }
30. if(fm<dblEps)return1;
31. if(l>ljt)return0;
33. l+=1;
34. u=p*lChannelCount+q;w=p*lChannelCount+p;t=q*lChannelCount+p;s=q*lChannelCount+q;
36. x=-pdblCof[p][q];
37. y=(pdblCof[q][q]-pdblCof[p][p])/2.0;
39. omega=x/sqrt(x*x+y*y);
40. if(y<0)omega=-omega;
42. sn=1.0+sqrt(1.0-omega*omega);
43. sn=omega/sqrt(2.0*sn);
44. cn=sqrt(1.0-sn*sn);
45. fm=pdblCof[p][p];
47. pdblCof[p][p]=fm*cn*cn+pdblCof[q][q]*sn*sn+pdblCof[p][q]*omega;
48. pdblCof[q][q]=fm*sn*sn+pdblCof[q][q]*cn*cn-pdblCof[p][q]*omega;
49. pdblCof[p][q]=0.0;
50. pdblCof[q][p]=0.0;
52. for(j=0;j<lChannelCount;j++)
53. if((j!=p)&&(j!=q))
54. {
55. fm=pdblCof[p][j];
56. pdblCof[p][j]=fm*cn+pdblCof[q][j]*sn;
57. pdblCof[q][j]=-fm*sn+pdblCof[q][j]*cn;
58. }
60. for(i=0;i<lChannelCount;i++)
61. if((i!=p)&&(i!=q)){
63. fm=pdblCof[i][p];
64. pdblCof[i][p]=fm*cn+pdblCof[i][q]*sn;
65. pdblCof[i][q]=-fm*sn+pdblCof[i][q]*cn;
66. }
68. for(i=0;i<lChannelCount;i++)
69. {
70. fm=pdblVects[i*lChannelCount+p];
71. pdblVects[i*lChannelCount+p]=fm*cn+pdblVects[i*lChannelCount+q]*sn;
72. pdblVects[i*lChannelCount+q]=-fm*sn+pdblVects[i*lChannelCount+q]*cn;
73. }
74. }
75. return1;
76. }
78. //根据特征值从大到小排列特征向量矩阵
79. /*
80. pfMatrix[in][out]-----上一步输出的协方差矩阵
81. nBandNum[in]------图像的输入波段数
82. pdblVects[out]----上一步输出的特征向量矩阵
83. */
84. staticvoidSortEigenvector(double**pfMatrix,intnBandNum,std::vector<double>&pfVector)
85. {
86. longp;
87. doublef;
88. doubleT;
90. intcount=nBandNum;
91. for(inti=0;i<count;i++)
92. {
93. T=pfMatrix[i][i];
94. p=i;
95. for(intj=i;j<count;j++)
96. if(T<pfMatrix[j][j])
97. {
98. T=pfMatrix[j][j];
99. p=j;
100. }
101. if(p!=i)
102. {
103. f=pfMatrix[p][p];
104. pfMatrix[p][p]=pfMatrix[i][i];
105. pfMatrix[i][i]=f;
107. for(intj=0;j<count;j++)
108. {
109. f=pfVector[j*count+p];
110. pfVector[j*count+p]=pfVector[j*count+i];
111. pfVector[j*count+i]=f;
112. }
113. }
114. }
115. }

执行上面两步之后，所得到的特征矩阵为用于和原图像相乘的矩阵A.

对一幅TM图像的第1，2，3，4，5，7通道执行PCA操作，效果图如下：

第一主成分：

第二主成分：

第三主成分：

第四主成分：

第五主成分：

第六主成分：

从上面可以看出，正确的图像看起来视觉非常和谐，毫无刺目的感觉。