2009年12月29日火曜日

主成分分析法图像处理

1.主成分的一般定义

设有随机变量X1,X2,…,Xp, 其样本均数记为 ,,…, ,样本标准差记为S1,S2,…,Sp。首先作标准化变换

我们有如下的定义:

(1) 若C1=a11x1+a12x2+ … +a1pxp, ,且使 Var(C1)最大,则称C1为第一主成分;

(2) 若C2=a21x1+a22x2+…+a2pxp, ,(a21,a22,…,a2p)垂直于(a11,a12,…,a1p),且使Var(C2)最大,则称C2为第二主成分;

(3) 类似地,可有第三、四、五…主成分,至多有p个。

2.主成分的性质

主成分C1,C2,…,Cp具有如下几个性质:

(1) 主成分间互不相关,即对任意i和j,Ci 和Cj的相关系数

Corr(Ci,Cj)=0 i ¹ j

(2) 组合系数(ai1,ai2,…,aip)构成的向量为单位向量,

(3) 各主成分的方差是依次递减的, 即

Var(C1)≥Var(C2)≥…≥Var(Cp)

(4) 总方差不增不减, 即

Var(C1)+Var(C2)+ … +Var(Cp)

=Var(x1)+Var(x2)+ … +Var(xp)

=p

这一性质说明,主成分是原变量的线性组合,是对原变量信息的一种改组,主成分不增加总信息量,也不减少总信息量。

(5) 主成分和原变量的相关系数 Corr(Ci,xj)=aij =aij

(6) 令X1,X2,…,Xp的相关矩阵为R, (ai1,ai2,…,aip)则是相关矩阵R的第i个特征向量(eigenvector)。而且,特征值li就是第i主成分的方差, 即

Var(Ci)= li

其中li为相关矩阵R的第i个特征值(eigenvalue)

l1≥l2≥…≥lp≥0

3.主成分的数目的选取

前已指出,设有p个随机变量,便有p个主成分。由于总方差不增不减,C1,C2等前几个综合变量的方差较大,而Cp,Cp-1等后几个综合变量的方差较小, 严格说来,只有前几个综合变量才称得上主(要)成份,后几个综合变量实为“次”(要)成份。实践中总是保留前几个,忽略后几个。

保留多少个主成分取决于保留部分的累积方差在方差总和中所占百分比(即累计贡献率),它标志着前几个主成分概括信息之多寡。实践中,粗略规定一个百分比便可决定保留几个主成分;如果多留一个主成分,累积方差增加无几,便不再多留。

4.主成分回归

主成分分析本身往往并不是目的,而是达到目的的一种手段。因此,它多用在大型研究项目的某个中间环节。例如,把它用在多重回归中,便产生了主成分回归。另外,它还可以用于聚类、判别分析等。

在多重回归曾指出,当自变量间高度相关时,某些回归参数的估计值极不稳定,甚至出现有悖常理、难以解释的情形。这时,可先采用主成分分析产生若干主成分,它们必定会将相关性较强的变量综合在同一个主成分中,而不同的主成分又是互相独立的。只要多保留几个主成分,原变量的信息不致过多损失。然后,以这些主成分为自变量进行多重回归就不会再出现共线性的困扰。如果原有p个自变量X1,X2,…,Xp,那么,采用全部p个主成分所作回归完全等价于直接对原变量的回归;采用一部分主成分所作回归虽不完全等价于对原变量的回归,但往往能摆脱某些虚假信息,而出现较合理的结果。

5.主成分分析与图像处理

在研究主成分分析方法基本原理的基础上,用MATLAB软件编写了的主成分分析代码用于神经元结构信息的MOST图片处理的程序,并 MATLB软件对图像进行分析,并对处理结果作比较分析。最后证明图像作主成分分析取得了较好的效果,具有一定的可靠性和实用性。

关于MOST图片的一些信息:MOST图片不适合用主成分分析除噪,去除非神经元的冗余信息,因为主成分分析的精髓在于变换和降维,需要有一个特征向量,变换提取出主成分,而MOST数据为灰度图,一个像素用一个0~255的数字表示,且一个位置只有一幅图像,特征向量为一维,无所谓提取,降维。但当把小块图片区域取出,所有像素点灰度构成一个特征向量,则这个向量维度很高(100*100的图片就是10'000维),但可以考虑用PCA提取主成分,降维,再用SVM,ANN等算法识别小块区域中是否有神经元,实现自动分割。

这里是采用记录CSD过程的图片集为例,应用Salwa Saleh编写的matlab代码,使用PCA方法处理后,得到了表现CSD扩散过程的图片。

6.使用PCA算法处理CSD图片集

(1)、材料:在小鼠发生CSD现象的时候,使用内源信号光学成像(OISI)技术,记录到的时间序列图片集

(2)、使用Salwa Saleh的matlab代码,运行程序,得到如下结果:





图一、PCA方法处理得到的一些成分。



图二、PCA方法处理后,与原始图片相减,得到的CSD扩散过程。

由处理结果可见,PCA方法提取出的主成份,清晰地表现出了CSD扩散过程。

附:Salwa Saleh的matlab代码源码

function [psi,a,ev]=pca(f,N)

% [psi,a,ev]=pca(f,N)

%

% Principal Component Analysis using by the Snapshot Method. (Sirovich, 1987)

% Use SVD instead if (# of images) > (# of pixels/image).

%

% MODEL: f = a * diag(sqrt(ev)) * psi'

%

% INPUT:

% f = matrix of images. (individual images are row vectors)

% N = the number of eigenpairs to compute. [default ALL]

%

% OUTPUT:

% psi = column vectors of spatial principal components.

% a = comumn vectors of temporal principal components.

% ev = PCA eigenvalues in descending order.

if size(f,1) > size(f,2)

disp('ERROR: The number of images exceeds the number of pixels/iamge.')

disp(' Consider using SVD instead.')

warning off; return;

end

disp('Forming the correlation matrix ...');

cor = f*f'; % the pixel-correlation matrix.

if nargin<2 % If N is not specified, compute all eigenpairs.

N=size(f,1);

end;

disp('Diagonalizing ...')

[a,ev]=eig(cor); % compute the temporal eigenvectors/values.

[ev,ind]=sort(diag(ev)); % sort in ascending magnitude.

ev=flipud(ev); % switch from ascending to descending.

a=a(:,flipud(ind)); % order the temporal eigenvectors accordingly.

disp('Calculating PCA eigenvectors ...')

a=a(:,1:N);

psi=f'*a*diag(ev(1:N).^-0.5); % compute the spatial eigenvectors.

参考文献:

[1]Salwa Saleh,Shangbin Chen,Dong Chen,Qingming Luo,Using Principal Component Analysis to Detect the

Wavefront of Cortical Spreading Depression

[2]Lu Chen(2007)“主元分析(PCA)理论分析及应用”

http://www.cad.zju.edu.cn/home/chenlu/pca.htm

[3]Jonathon Shlens.(2005)“A Tutorial on Principal Component Analysis”

http://www.snl.salk.edu/~shlens/pub/notes/pca.pdf

[4]aizaixiyuanqian1985提供(2008)“pca方法实现详细步骤”

http://www.pudn.com/downloads102/sourcecode/math/detail417142.html

[4]“[人脸识别Step2]对人脸图片作PCA降维”

http://www.sw-china.org/forum/simple/index.php?t40.html