小波模极大值用于图像融合

%  小波模极大值用于图像融合

clc;clear
X1=imread('lena64.bmp'); 
 X1=imresize(X1,[184 184]);
I1=(X1);   %把彩色图像转换为灰度图像
figure(1);
subplot(131);imshow(I1); %显示图像1
title('MRI');
X2=imread('21.jpg');
X2=imresize(X2,[184 184]);
I2=rgb2gray(X2/2);   %把彩色图像转换为灰度图像
subplot(132);imshow(I2); %显示图像2
title('CT');

I11=im2double(I1);
I22=im2double(I2);

n=1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%小波系数分解得到各个频率分量系数
[c1,s1]=wavedec2(I11,n,'db4');
[c2,s2]=wavedec2(I22,n,'db4');
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%对第一幅图像（MRI）的处理：
% [h11,v11,d11]=detcoef2('all',c1,s1,1);%提取第一层二维小波分解的细节分量（高频系数）
% [h12,v12,d12]=detcoef2('all',c1,s1,2);%提取第二层二维小波分解的细节分量（高频系数）
% [h13,v13,d13]=detcoef2('all',c1,s1,1);%提取第三层二维小波分解的细节分量（高频系数）
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%第二幅图像（CT）的处理：
% [h21,v21,d21]=detcoef2('all',c2,s2,1);%提取第一层二维小波分解的细节分量（高频系数）
% [h22,v22,d22]=detcoef2('all',c2,s2,2);
[h23,v23,d23]=detcoef2('all',c2,s2,1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%提取最高层（第三层）的低频系数
a31=appcoef2(c1,s1,'db4',1);%提取多层二维小波分解的近似分量（低频系数）
a32=appcoef2(c2,s2,'db4',1);%

A3=(a31+a32);

X=I22;
[SIZEX,SIZEY]=size(X);  %  图像尺寸
    
%  多尺度
m=1.0;
delta=2^m;
%  构造高斯函数的偏导
N=20;  %  滤波器长度（需要调整，必须是偶数）
A=-1/sqrt(2*pi);  %  幅度
for index_x=1:N;
    for index_y=1:N;
        x=index_x-(N+1)/2;
        y=index_y-(N+1)/2;
        phi_x(index_x,index_y)=A*(x/delta^2).*exp(-(x.*x+y.*y)/(2*delta^2));
        phi_y(index_x,index_y)=A*(y/delta^2).*exp(-(x.*x+y.*y)/(2*delta^2));
    end
end;
phi_x=phi_x/norm(phi_x);  %  能量归一化
phi_y=phi_y/norm(phi_y);  %  能量归一化
%  对图象做行列卷积
Gx=conv2(X,phi_x,'same');
Gy=conv2(X,phi_y,'same');
%  求梯度
Grads=sqrt((Gx*Gx)+(Gy*Gy));
%  求幅角（梯度方向）
angle_array=zeros(SIZEX,SIZEY);  %  角度
%  遍历
for i=1:SIZEX;
    for j=1:SIZEY
        if (abs(Gx(i,j))>eps*100)  %  x的绝对值足够大
            p=atan(Gy(i,j)/Gx(i,j))*180/pi;  %  反正切求角度值(1,4象限)
            if (p<0)        %  负的幅角（4象限）
                p=p+360;
            end;
            if (Gx(i,j)<0 & p>180)     %  2象限的特殊处理
                p=p-180;
            elseif (Gx(i,j)<0 & p<180) %  3象限的特殊处理
                p=p+180;
            end
        else  %  90或270度
            p=90;
        end
        angle_array(i,j)=p;  %  幅角赋值
    end
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%对CT边缘图像 梯度幅值 进行小波分解
n=1
[c4,s4]=wavedec2(Grads,n,'db4');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%只需提取高频系数
% [h41,v41,d41]=detcoef2('all',c4,s4,1);%提取第一层二维小波分解的细节分量（高频系数）
% [h42,v42,d42]=detcoef2('all',c4,s4,2);%提取第二层二维小波分解的细节分量（高频系数）
[h43,v43,d43]=detcoef2('all',c4,s4,1);%提取第三层二维小波分解的细节分量（高频系数）

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%对CT边缘图像 角度 进行小波分解
n=1
[c5,s5]=wavedec2(angle_array,n,'db4');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%只需提取高频系数
% [h51,v51,d51]=detcoef2('all',c5,s5,1);%提取第一层二维小波分解的细节分量（高频系数）
% [h52,v52,d52]=detcoef2('all',c5,s5,2);%提取第二层二维小波分解的细节分量（高频系数）
[h53,v53,d53]=detcoef2('all',c5,s5,1);%提取第三层二维小波分解的细节分量（高频系数）

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%角度检测  
%%%%%%%%第三层 水平系数 的检测  其中H3(i,j)为融合后第一层小波系数的水平分量
%%%通过角度值 确定幅度值的方向，再与相邻两点值比较，取最大那个值对应的CT原图小波系数作为融合后该点的小波系数
%%%h53是角度的水平分量，h43是幅度的水平分量
[SIZEX,SIZEY]=size(h53);
H3=h23;
for i=2:SIZEX-1
    for j=2:SIZEY-1
        if ((h53(i,j)>=(-22.5) & h53(i,j)<=22.5) | ...
            (h53(i,j)>=(180-22.5) & h53(i,j)<=(180+22.5)))     %  0/180
            if (h43(i,j)>h43(i+1,j) & h43(i,j)>h43(i-1,j))
               H3(i,j)=h23(i,j);
            elseif h43(i+1,j)>h43(i-1,j)
               H3(i,j)=h23(i+1,j); 
            else
                H3(i,j)=h23(i-1,j);
            end

        elseif ((h53(i,j)>=(90-22.5) & h53(i,j)<=(90+22.5)) | ...
                (h53(i,j)>=(270-22.5) & h53(i,j)<=(270+22.5))) %  90/270
            if (h43(i,j)>h43(i,j+1) & h43(i,j)>h43(i,j-1))
                H3(i,j)=h23(i,j);
            elseif h43(i,j+1)>h43(i,j-1)
                H3(i,j)=h23(i,j+1);
            else
                H3(i,j)=h23(i,j-1);
            end

        elseif ((h53(i,j)>=(45-22.5) & h53(i,j)<=(45+22.5)) | ...
                (h53(i,j)>=(225-22.5) & h53(i,j)<=(225+22.5))) %  45/225
            if (h43(i,j)>h43(i+1,j+1) & h43(i,j)>h43(i-1,j-1))
                H3(i,j)=h23(i,j);
            elseif h43(i+1,j+1)>h43(i-1,j-1)
                H3(i,j)=h23(i+1,j+1);
            else
                 H3(i,j)=h23(i-1,j-1);
            end

        else  %  135/215
            if (h43(i,j)>h43(i+1,j-1) & h43(i,j)>h43(i-1,j+1))
                H3(i,j)=h23(i,j);
            elseif h43(i+1,j-1)>h43(i-1,j+1)
                 H3(i,j)=h23(i+1,j-1);
            else
                H3(i,j)=h23(i-1,j+1);
            end

        end
    end
end

。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。

 

D-20