ICA图像独立主成分分析
参考链接：http://wenku.baidu.com/link?url= ... 0lJ4dJYRj8CSqWHjBae
http://wenku.baidu.com/link?url= ... 0lJ4dJYRj8CSqWHjBae
http://wenku.baidu.com/view/1578be42336c1eb91a375d0e.html
源信号s --> 混合矩阵A --> 中心化、白化-->分离矩阵-->分离信号Y
具体的实现代码如下：
clc,clear,close all
warning off
%读取和显示图像
I00= imread('1.bmp');
I01= imread('2.bmp');
I02= imread('3.bmp');
I1=rgb2gray(I00); I1 = imresize(I1,);
I2=rgb2gray(I01); I2 = imresize(I2,);
I3=rgb2gray(I02); I3 = imresize(I3,);
figure(1)
subplot(321),imshow(I1),
subplot(322),imhist(I1),
subplot(323),imshow(I2),
subplot(324),imhist(I2),
subplot(325),imshow(I3),
subplot(326),imhist(I3),
%对信号进行随机混合得到仿真观测信号
s1=reshape(I1,);
s2=reshape(I2,);
s3=reshape(I3,);
%将s1、s2、s3变换成一个1*65536的行矩阵
s1=reshape(I1,);
s2=reshape(I2,);
s3=reshape(I3,); % 将s1、s2、s3变换成一个1*65536的行矩阵
s=; % s为一个3*65536的矩阵
sig=double(s); % sig为一个double型的矩阵
A=rand(size(sig,1)); % 生成一个大小为3*3，取值在0-1之间的随机矩阵
mixedsig=A*sig; % 混合信号为一个3*65535的矩阵
ms1=reshape(mixedsig(1,:),);
ms2=reshape(mixedsig(2,:),);
ms3=reshape(mixedsig(3,:),); % 将ms1、ms2、ms3四舍五入转换成无符号整型数
MI1=uint8(round(ms1));
MI2=uint8(round(ms2));
MI3=uint8(round(ms3));
%显示混合图像?figure(2)?
subplot(131),imshow(MI1),title('混合图片1');
subplot(132),imshow(MI2),title('混合图片2');
subplot(133),imshow(MI3),title('混合图片3');
mixeds_bak=mixedsig; % 将混合后的数据备份，以便在恢复时直接调用
mixeds_mean=zeros(3,1); % 标准化
for i=1:3
mixeds_mean(i)=mean(mixedsig(i,:));
end % 计算mixedsig的均值和方差
for i=1:3
for j=1:size(mixedsig,2)
mixedsig(i,j)=mixedsig(i,j)-mixeds_mean(i);
end
end%白化
mixeds_cov=cov(mixedsig'); % cov为求协方差的函数
=eig(mixeds_cov); % 对图片矩阵的协方差函数进行特征值分解
Q=inv(sqrt(D))*(E)'; % Q为白化矩阵
mixeds_white=Q*mixedsig; %mixeds_white为白化后的图片矩阵
I=cov(mixeds_white'); % I应为单位阵
%%%%%%%%%%%%%%%%%%FastICA算法%%%%%%%%%%%%%%%%%%%%%%%%%%?
X=mixeds_white;%以下对X进行操作
=size(X);
numofIC=variablenum;%在此应用中，独立元个数等于变量个数
B=zeros(numofIC,variablenum);%初始化列向量b的寄存矩阵，B=
for r=1:numofIC
i=1;
maxiterationsnum=150;%设置最大迭代次数
b=2*(rand(numofIC,1)-.5);%随机设置b的初值
b=b/norm(b);%对b标准化
while i<=maxiterationsnum+1
if i==maxiterationsnum %循环结束处理
fprintf('\n?第%d分量在%d次迭代内并不收敛.',r,maxiterationsnum);
break;
end
bold=b;%初始化前一步b的寄存器
u=1;
t=X'*b;
g=t.^3;
dg=3*t.^2;
b=((1-u)*t'*g*b+u*X*g)/samplenum-mean(dg)*b;%核心公式
b=b-B*B'*b; % 对b正交化
b=b/norm(b);
if abs(abs(b'*bold)-1)<1e-9 % 如果收敛，则
B(:,r)=b;%保存所得向量b
break;
end
i=i+1;
end
end
%数据复原?
ICAeds=B'*Q*mixeds_bak;
ICAeds_bak=ICAeds;
ICAeds=abs(55*ICAeds);
Is1=reshape(ICAeds(1,:),);
Is2=reshape(ICAeds(2,:),);
Is3=reshape(ICAeds(3,:),);
II1=uint8(round(Is1));
II2=uint8(round(Is2));
II3=uint8(round(Is3)); %显示分离后的图像
figure(3)
subplot(321),imshow(II1),title('分离出的图片1'),
subplot(322),imhist(II1),title('分离图片1的直方图');
subplot(323),imshow(II2),title('分离出的图片2'),
subplot(324),imhist(II2),title('分离图片2的直方图');
subplot(325),imshow(II3),title('分离出的图片3'),
subplot(326),imhist(II3),title('分离图片3的直方图');
III1=imsubtract(I1,II1);
III2=imsubtract(I2,II2);
III3=imsubtract(I3,II3);
%显示分离后的图片与原图的差值图以及对应的直方图%
figure(4)
subplot(321),imshow(III1),title('差值图片1'),
subplot(322),imhist(III1),title('差值图片1的直方图');
subplot(323),imshow(III2),title('差值图片2'),
subplot(324),imhist(III2),title('差值图片2的直方图');
subplot(325),imshow(III3),title('差值图片3'),
subplot(326),imhist(III3),title('差值图片3的直方图');
ICA各类程序汇总，打包下载如下：
链接：http://pan.baidu.com/s/1o8Iwy2u 密码：tw0w
具体参见链接如下：
http://halcom.cn/forum.php?mod=viewthread&tid=204&page=1&extra=#pid229
页:
[1]