图像分离

I1=imread ('input1.jpg');          
I2=imread ('input2.jpg'); 
I3=imread ('input3.jpg'); 
I1=rgb2gray(I1);
I2=rgb2gray(I2);
I3=rgb2gray(I3);
subplot(4,3,1),imshow(I1),title('输入图像1'), 
subplot(4,3,2),imshow(I2),title('输入图像2'), 
subplot(4,3,3),imshow(I3),title('输入图像3'), 
 
% 计算图片数据的维数 
[x1,y1,z1]=size(I1); 
[x2,y2,z2]=size(I2); 
[x3,y3,z3]=size(I3); 
 
% 将其重新排列为一维行向量并组成矩阵 
s1=[reshape(I1,[1,x1*y1*z1])]; 
s2=[reshape(I2,[1,x2*y2*z2])]; 
s3=[reshape(I3,[1,x3*y3*z3])]; 
S_all=[s1;s2;s3];   % 图片个数即为变量数,图片的像素数即为采样数 
                              % 因此S_all是一个变量个数*采样个数的矩阵 
S=double(S_all);      % 将图像数据转换为双精度格式 ,图片处理前要将数据格式用double转化为双精度格式(以保证计算精度)
Sweight=rand(size(S_all,1));     % 取一随机矩阵,作为信号混合的权矩阵 ,注意size(S_all,1)产生三个数
MixedS=Sweight*S;           % 得到三个图像的混合信号矩阵 
 
% 将混合矩阵重新排列为原始的图片矩阵并输出 
ms1=reshape(MixedS(1,:),[x1,y1,z1]); 
ms2=reshape(MixedS(2,:),[x2,y2,z2]); 
ms3=reshape(MixedS(3,:),[x3,y3,z3]); 
 
I1_mixed=uint8(round(ms1)); %显示图像前用uint8转换数据格式
I2_mixed=uint8(round(ms2)); 
I3_mixed=uint8(round(ms3)); 
 
subplot(4,3,4),imshow(I1_mixed),title('混合图像1'), 
subplot(4,3,5),imshow(I2_mixed),title('混合图像2'), 
subplot(4,3,6),imshow(I3_mixed),title('混合图像3'), 
 
MixedS_bak=MixedS;      % 将混合后的数据备份,以便在恢复时直接调用 
 
%%%%%%%%%%%%%%%%%%%%%%%%%%  标准化  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
MixedS_mean=zeros(3,1); 
for i=1:3 
    MixedS_mean(i)=mean(MixedS(i,:)); 
end                                         % 计算MixedS的均值与方差 


for i=1:3 
    for j=1:size(MixedS,2) 
        MixedS(i,j)=MixedS(i,j)-MixedS_mean(i); 
    end 
end 
%%%%%%%%%%%%%%%%%%%%%%%%%%%  白化  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
 
MixedS_cov=cov(MixedS');                    % cov为求协方差的函数 
[E,D]=eig(MixedS_cov);                     % 对图片矩阵的协方差函数进行特征值分解 
Q=inv(sqrt(D))*(E)';                        % Q为白化矩阵 
MixedS_white=Q*MixedS;                      % MixedS_white为白化后的图片矩阵 
IsI=cov(MixedS_white');                     % IsI应为单位阵             
 
%%%%%%%%%%%%%%%%%%%%%%%% FASTICA算法  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
X=MixedS_white;                             % 以下算法将对X进行操作 
[VariableNum,SampleNum]=size(X); 
numofIC=VariableNum;                        % 在此应用中,独立元个数等于变量个数 
B=zeros(numofIC,VariableNum);               % 初始化列向量b的寄存矩阵,B=[b1,b2,...,bd] 
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;     
                                            % 核心公式,参见理论部分公式2.52 
        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;                     % 参见理论部分公式2.55 
ICAedS_bak=ICAedS; 
ICAedS=abs(55*ICAedS); 
 
% t1=ICAedS(1,:); 
% t2=ICAedS(2,:); 
% t3=ICAedS(3,:); 
%  
% for i=1:size(t1,2) 
%     t1(1,i)=t1(1,i)+abs(min(t1)); 
%     t1(1,i)=t1(1,i)*256/(max(t1)-min(t1)); 
%         t2(1,i)=t2(1,i)+abs(min(t2)); 
%     t2(1,i)=t2(1,i)*256/(max(t2)-min(t2)); 
%         t3(1,i)=t3(1,i)+abs(min(t3)); 
%     t3(1,i)=t3(1,i)*256/(max(t3)-min(t3)); 
%          
% end 
%  
% ICAedS(1,:)=t1; 
% ICAedS(2,:)=t2; 
% ICAedS(3,:)=t3; 
 
 
% 将计算后的混合矩阵重新排列为图片矩阵并输出 
is1=reshape(ICAedS(1,:),[x1,y1,z1]); 
is2=reshape(ICAedS(2,:),[x2,y2,z2]); 
is3=reshape(ICAedS(3,:),[x3,y3,z3]); 
 
I1_icaed =uint8 (round(is1)); 
I2_icaed =uint8 (round(is2)); 
I3_icaed =uint8 (round(is3)); 
 
subplot(4,3,7),imshow(I1_icaed),title('ICA解混图像1'), 
subplot(4,3,8),imshow(I2_icaed),title('ICA解混图像2'), 
subplot(4,3,9),imshow(I3_icaed),title('ICA解混图像3'), 
 
%%%%%%%%%%%%%%%%%%%%%%%%%%  PCA计算并构图  %%%%%%%%%%%%%%%%%%%%%%%%%% 
[V,D]=eig(MixedS_cov); 
Vtmp=zeros(size(V,1),1); 
for j=1:2 
    for i=1:2 
        if D(i,i)<D(i+1,i+1) 
            tmp=D(i,i);Vtmp=V(:,i); 
            D(i,i)=D(i+1,i+1);V(:,i)=V(:,i+1); 
            D(i+1,i+1)=tmp;V(:,i+1)=Vtmp; 
        end 
    end 
end 
 
t1=(MixedS'*V(:,1))'; 
t2=(MixedS'*V(:,2))';    
t3=(MixedS'*V(:,3))'; 
 
for i=1:size(t1,2)                    % 亮度调整                        
    t1(1,i)=t1(1,i)+abs(min(t1)); 
    t1(1,i)=t1(1,i)*256/(max(t1)-min(t1)); 
    t2(1,i)=t2(1,i)+abs(min(t2)); 
    t2(1,i)=t2(1,i)*256/(max(t2)-min(t2)); 
    t3(1,i)=t3(1,i)+abs(min(t3)); 
    t3(1,i)=t3(1,i)*256/(max(t3)-min(t3)); 
end 
%%%%%%%%%%%%%%%%%%%%%%%   构图    %%%%%%%%%%%%%%%%%%%%%% 
ts1=reshape(t1,[x1,y1,z1]); 
ts2=reshape(t2,[x2,y2,z2]); 
ts3=reshape(t3,[x3,y3,z3]); 
 
T1_icaed =uint8 (round(ts1)); 
T2_icaed =uint8 (round(ts2)); 
T3_icaed =uint8 (round(ts3)); 
 
subplot(4,3,10),imshow(T1_icaed),title('PCA解混图像1'), 
subplot(4,3,11),imshow(T2_icaed),title('PCA解混图像2'), 
subplot(4,3,12),imshow(T3_icaed),title('PCA解混图像3'), 
 
%%%%%%%%%%%%%%%%%%%%%%%%%%  输出一维矢量图  %%%%%%%%%%%%%%%%%%%%%%%%%% 
% subplot(3,3,1),plot(S(1,:)),title('输入图像1一维矢量图') 
% subplot(3,3,2),plot(S(2,:)),title('输入图像2一维矢量图') 
% subplot(3,3,3),plot(S(3,:)),title('输入图像3一维矢量图') 
% subplot(3,3,4),plot(ICAedS(1,:)),title('ICA解混图像1一维矢量图') 
% subplot(3,3,5),plot(ICAedS(2,:)),title('ICA解混图像2一维矢量图') 
% subplot(3,3,6),plot(ICAedS(3,:)),title('ICA解混图像3一维矢量图') 
% subplot(3,3,7),plot(t1),title('PCA解混图像1一维矢量图') 
% subplot(3,3,8),plot(t2),title('PCA解混图像2一维矢量图') 
% subplot(3,3,9),plot(t3),title('PCA解混图像3一维矢量图') 




im = imread(imgname);
im=rgb2gray(im);
imnoise_1=imnoise(im,'speckle',1);
imnoise_1=uint8(imnoise_1);
g2=size(imnoise_1);
figure(2)imshow(imnoise_1);
title('Noise image');
imnoise_2=double(imnoise_1);
imnoise_3=imnoise_2+var;
imnoise_3=log(imnoise_3);
immoise_3=double(imnoise_3);
imnoise_4=exp(imnoise_3);
imnoise_4=imnoise_4-var;
imnoise_5=uint8(imnoise_4);
figure(5)imshow(imnoise_5)title('recovery speckle image');


end
阅读更多
文章标签: 图像
个人分类: 理论
下一篇信噪比
想对作者说点什么? 我来说一句

没有更多推荐了,返回首页

关闭
关闭