ICA图像独立主成分分析

Halcom

参考链接：
http://wenku.baidu.com/link?url= ... 0lJ4dJYRj8CSqWHjBae
http://wenku.baidu.com/link?url= ... 0lJ4dJYRj8CSqWHjBae
http://wenku.baidu.com/view/1578be42336c1eb91a375d0e.html

源信号s --> 混合矩阵A --> 中心化、白化-->分离矩阵-->分离信号Y

具体的实现代码如下：

1. clc,clear,close all
   
2. warning off
   
3. %读取和显示图像
   
4. I00= imread('1.bmp');
   
5. I01= imread('2.bmp');
   
6. I02= imread('3.bmp');
   
7. I1=rgb2gray(I00); I1 = imresize(I1,[256,256]);
   
8. I2=rgb2gray(I01); I2 = imresize(I2,[256,256]);
   
9. I3=rgb2gray(I02); I3 = imresize(I3,[256,256]);
   
10. figure(1)
    
11. subplot(321),imshow(I1),
    
12. subplot(322),imhist(I1),
    
13. subplot(323),imshow(I2),
    
14. subplot(324),imhist(I2),
    
15. subplot(325),imshow(I3),
    
16. subplot(326),imhist(I3),
    
17. 
    
18. %对信号进行随机混合得到仿真观测信号
    
19. s1=reshape(I1,[1,256*256]);
    
20. s2=reshape(I2,[1,256*256]);
    
21. s3=reshape(I3,[1,256*256]);
    
22. %将s1、s2、s3变换成一个1*65536的行矩阵
    
23. s1=reshape(I1,[1,256*256]);
    
24. s2=reshape(I2,[1,256*256]);
    
25. s3=reshape(I3,[1,256*256]);     % 将s1、s2、s3变换成一个1*65536的行矩阵
    
26. s=[s1;s2;s3];                       % s为一个3*65536的矩阵
    
27. sig=double(s);                     % sig为一个double型的矩阵
    
28. A=rand(size(sig,1));             % 生成一个大小为3*3，取值在0-1之间的随机矩阵
    
29. mixedsig=A*sig;                  % 混合信号为一个3*65535的矩阵
    
30. ms1=reshape(mixedsig(1,:),[256,256]);
    
31. ms2=reshape(mixedsig(2,:),[256,256]);
    
32. ms3=reshape(mixedsig(3,:),[256,256]);    % 将ms1、ms2、ms3四舍五入转换成无符号整型数
    
33. MI1=uint8(round(ms1));
    
34. MI2=uint8(round(ms2));
    
35. MI3=uint8(round(ms3));
    
36. %显示混合图像?figure(2)?
    
37. subplot(131),imshow(MI1),title('混合图片1');
    
38. subplot(132),imshow(MI2),title('混合图片2');
    
39. subplot(133),imshow(MI3),title('混合图片3');
    
40. mixeds_bak=mixedsig;                                        % 将混合后的数据备份，以便在恢复时直接调用
    
41. mixeds_mean=zeros(3,1);                                   % 标准化
    
42. for i=1:3
    
43.     mixeds_mean(i)=mean(mixedsig(i,:));
    
44. end                                                                      % 计算mixedsig的均值和方差
    
45. for i=1:3
    
46.     for j=1:size(mixedsig,2)
    
47.         mixedsig(i,j)=mixedsig(i,j)-mixeds_mean(i);
    
48.     end
    
49. end  %白化
    
50. mixeds_cov=cov(mixedsig');     % cov为求协方差的函数
    
51. [E,D]=eig(mixeds_cov);            % 对图片矩阵的协方差函数进行特征值分解
    
52. Q=inv(sqrt(D))*(E)';                 % Q为白化矩阵
    
53. mixeds_white=Q*mixedsig;      %  mixeds_white为白化后的图片矩阵
    
54. I=cov(mixeds_white');               % I应为单位阵
    
55. %%%%%%%%%%%%%%%%%%FastICA算法%%%%%%%%%%%%%%%%%%%%%%%%%%?
    
56. X=mixeds_white;%以下对X进行操作
    
57. [variablenum,samplenum]=size(X);
    
58. numofIC=variablenum;%在此应用中，独立元个数等于变量个数
    
59. B=zeros(numofIC,variablenum);%初始化列向量b的寄存矩阵，B=[b1,b2,??,bd]
    
60. for r=1:numofIC
    
61.     i=1;
    
62.     maxiterationsnum=150;%设置最大迭代次数
    
63.     b=2*(rand(numofIC,1)-.5);%随机设置b的初值
    
64.     b=b/norm(b);%对b标准化
    
65.     while i<=maxiterationsnum+1
    
66.         if i==maxiterationsnum %循环结束处理
    
67.             fprintf('\n?第%d分量在%d次迭代内并不收敛.',r,maxiterationsnum);
    
68.             break;
    
69.         end
    
70.         bold=b;%初始化前一步b的寄存器
    
71.         u=1;
    
72.         t=X'*b;
    
73.         g=t.^3;
    
74.         dg=3*t.^2;
    
75.         b=((1-u)*t'*g*b+u*X*g)/samplenum-mean(dg)*b;%核心公式
    
76.         b=b-B*B'*b; % 对b正交化
    
77.         b=b/norm(b);
    
78.         if abs(abs(b'*bold)-1)<1e-9   % 如果收敛，则
    
79.             B(:,r)=b;%保存所得向量b
    
80.             break;
    
81.         end
    
82.         i=i+1;
    
83.     end
    
84. end
    
85. %数据复原?
    
86. ICAeds=B'*Q*mixeds_bak;
    
87. ICAeds_bak=ICAeds;
    
88. ICAeds=abs(55*ICAeds);
    
89. 
    
90. Is1=reshape(ICAeds(1,:),[256,256]);
    
91. Is2=reshape(ICAeds(2,:),[256,256]);
    
92. Is3=reshape(ICAeds(3,:),[256,256]);
    
93. II1=uint8(round(Is1));
    
94. II2=uint8(round(Is2));
    
95. II3=uint8(round(Is3)); %显示分离后的图像
    
96. figure(3)
    
97. subplot(321),imshow(II1),title('分离出的图片1'),
    
98. subplot(322),imhist(II1),title('分离图片1的直方图');
    
99. subplot(323),imshow(II2),title('分离出的图片2'),
    
100. subplot(324),imhist(II2),title('分离图片2的直方图');
     
101. subplot(325),imshow(II3),title('分离出的图片3'),
     
102. subplot(326),imhist(II3),title('分离图片3的直方图');
     
103. III1=imsubtract(I1,II1);
     
104. III2=imsubtract(I2,II2);
     
105. III3=imsubtract(I3,II3);
     
106. %显示分离后的图片与原图的差值图以及对应的直方图%
     
107. figure(4)
     
108. subplot(321),imshow(III1),title('差值图片1'),
     
109. subplot(322),imhist(III1),title('差值图片1的直方图');
     
110. subplot(323),imshow(III2),title('差值图片2'),
     
111. subplot(324),imhist(III2),title('差值图片2的直方图');
     
112. subplot(325),imshow(III3),title('差值图片3'),
     
113. subplot(326),imhist(III3),title('差值图片3的直方图');

复制代码

Halcom

ICA各类程序汇总，打包下载如下：
链接：http://pan.baidu.com/s/1o8Iwy2u 密码：tw0w
具体参见链接如下：
http://halcom.cn/forum.php?mod=v ... 1&extra=#pid229