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

作者：朱金灿

来源：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

*/

static int iJcobiMatrixCharacterValue(double** pdblCof, long lChannelCount, std::vector<double>& pdblVects, double dblEps,long ljt)

{

long i,j,p,q,u,w,t,s,l;

double fm,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) return 1;

if(l > ljt) return 0;

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;

}

}

return 1;

}

// 根据特征值从大到小排列特征向量矩阵

/*

pfMatrix [in][out]----- 上一步输出的协方差矩阵

nBandNum [in] ------ 图像的输入波段数

pdblVects [out] ---- 上一步输出的特征向量矩阵

*/

static void SortEigenvector(double** pfMatrix,int nBandNum,std::vector<double> &pfVector)

{

long p;

double f;

double T;

int count = nBandNum;

for(int i = 0; i < count ; i ++)

{

T = pfMatrix[i][i];

p = i;

for(int j = 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(int j = 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操作，效果图如下：

第一主成分：



第二主成分：



第三主成分：