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

作者：朱金灿

来源：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. static int iJcobiMatrixCharacterValue(double** pdblCof, long lChannelCount, std::vector<double>& pdblVects, double dblEps,long ljt)
10. {
11.     long i,j,p,q,u,w,t,s,l;
12.     double fm,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) return 1;
31.             if(l > ljt)   return 0;
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.     return 1;
76. }
78. // 根据特征值从大到小排列特征向量矩阵
79. /*
80. pfMatrix [in][out]----- 上一步输出的协方差矩阵
81. nBandNum [in] ------  图像的输入波段数
82. pdblVects  [out] ---- 上一步输出的特征向量矩阵
83. */
84.     static void SortEigenvector(double** pfMatrix,int nBandNum,std::vector<double> &pfVector)
85. {
86.     long p;
87.     double f;
88.     double T;
90.     int count = nBandNum;
91.     for(int i = 0; i < count ; i ++)
92.     {
93.         T = pfMatrix[i][i];
94.         p = i;
95.         for(int j = 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(int j = 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. }