主成分分析实现的一个心得-白红宇

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

阅读量：2430 次

发布时间：2019-05-10

本文共 3098 字，大约阅读时间需要 10 分钟。

作者：朱金灿

来源：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操作，效果图如下：

第一主成分：

第二主成分：

第三主成分：

第四主成分：

第五主成分：

第六主成分：

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

你可能感兴趣的文章