差分演化算法实现多聚焦图像融合论文算法实现（算法实现）

本文为上一篇文章的算法实现。

首先，我们来复习一下用matlab来显示图像，这个很简单，直接用imread就可以读取图像，然后用imshow显示就好了，如果想在一个figure中显示多个图片，可以用subplot。考虑图像融合，最简单的，直接像素相加，也可以完成两张图像的融合，但是很显然相同位置的像素值变大了，因此融合图像整体会比较亮。如下：

I = imread('clock1.jpg');
J = imread('clock2.jpg');
K = imadd(I,J);
figure;
subplot(131),imshow(I);
subplot(132),imshow(J);
subplot(133),imshow(K);

　　效果如图：

但是根据之前的分析，我们要做的是将源图像分块再融合，选两张原图像块中比较清晰的块放在融合图像中，用DE算法来确定最优块大小。一步一步来吧，我们可以这样分解一下任务：1、求一个图像的分块。2、用文章中推荐的SF算法当作清晰度函数来计算块和全局的清晰度。3、先任意设定块的大小，然后设计融合算法来完成图像的融合。4、找到DE算法（网上应该有，但是没有针对图像融合的），用DE算法求图像的最优块大小。

1、网上可以找到一些不错的图像分块算法，比如：

clc;clear all;close all;
I1 = imread('football.jpg');
rs = size(I1,1);cs = size(I1,2);  %rs：表示图像的行；cs：图像的列
sz = 64;                        %按64个像素进行分块，可自行设置
numr = rs/sz;                   %图像分块的行
numc = cs/sz;                 %图像分快的列
ch = sz; cw = sz;
t1 = (0:numr-1)*ch + 1; t2 = (1:numr)*ch;    %分别求得每一块图像的起始行的像素值
t3 = (0:numc-1)*cw + 1; t4 = (1:numc)*cw;      %分别求得每一块图像的起始列的像素值
 
%figure;
k=0;        %开始分块
for i =1 : numr
for j = 1 : numc
temp = I1(t1(i):t2(i),t3(j):t4(j),:);%暂存分块图像为temp
k = k + 1;
 
subplot(numr,numc,k);%显示分块图像
imshow(temp);
 
pause(0.5);
end
end

　　效果如下图：

这个算法还是不错的，但是不能直接用到我们的程序中，因为我们不需要显示出来（这不是重点），主要是我们分块后需要求每一小块图像的清晰度值，和另外的源图像相应的位置的图像块比较，然后选择清晰度值较大的作为融合图像的分块。现在的问题是怎么存储分块呢，经过不少的尝试，我发现用细胞数组来存储分块是不错的选择。因此在上面的双重for循环内加上如下代码：

1 2	`A{i,j} = temp;` `%将分块存储到细胞数组A中` `BSV{i,j} = spatial_frequencies(A{i,j});%将每一块的SF值也存储到BSV细胞数组中，方便后续的清晰度值比较`

　　ok，基本上没什么问题，把这部分封装成函数吧，命名为spilt.m。输入为一张图片，输出为两个细胞数组。如下：

function [A,BSV] = spilt(img)
A = cell(8,8);
BSV = cell(8,8);
I = imread(img);
rs = size(I,1);cs = size(I,2);  %rs：表示图像的行；cs：图像的列
sz = 64;                        %按64个像素进行分块，可自行设置
numr = rs/sz;                   %图像分块的行
numc = cs/sz;                 %图像分快的列
ch = sz; cw = sz;
t1 = (0:numr-1)*ch + 1; t2 = (1:numr)*ch;    %分别求得每一块图像的起始行的像素值
t3 = (0:numc-1)*cw + 1; t4 = (1:numc)*cw;      %分别求得每一块图像的起始列的像素值
 
%figure;
k=0;        %开始分块
for i =1 : numr
for j = 1 : numc
temp = I(t1(i):t2(i),t3(j):t4(j),:);%暂存分块图像为temp
k = k + 1;
 
%subplot(numr,numc,k);%显示分块图像
%imshow(temp);
   A{i,j} = temp;%将分块存储到细胞数组A中
   BSV{i,j} = spatial_frequencies(A{i,j});%将每一块的SF值也存储到BSV细胞数组中，方便后续的清晰度值比较
pause(0.5);
end
end

　　2、第一步已经完成了(事实上上面的代码已经用到了SF算法)，至于SF算法，网上已经有很多了，找了一个靠谱点的，输入为图像或者图像分块，输出为清晰度值。（就是这个代码本身问题不大，但是后面由它产生了一个bug，找了好久才找出来）如下：

function out_val=spatial_frequencies(img)<br>I =imread(img);
I=double(I);
[m,n]=size(I);
f=0.0;
rf=0.0;
cf=0.0;
for x=1:m
    for y=1:n-1
        rf=rf+(I(x,y+1)-I(x,y))^2;
    end
end
rf=sqrt(rf/m/n);
for x=1:m-1
    for y=1:n
        cf=cf+(I(x+1,y)-I(x,y))^2;
    end
end
cf=sqrt(cf/m/n);
f=sqrt(rf^2+cf^2);
out_val=f;

　　3、现在可以设计融合算法了，其实挺简单的，理解图像就是一个矩阵就可以了。我们先设计分块为8*8的图像融合吧。如下：

clc,clear;
[A,BSVa] = spilt2('clock1.jpg',8,8);
[B,BSVb] = spilt2('clock2.jpg',8,8);
for i =1:8
    for j = 1:8
        if BSVa{i,j} > BSVb{i,j} %比较Ai和Bi的块的清晰度，
            F{i,j} = A{i,j};     %将清晰度比较高的块放到融合块中。
        elseif BSVa{i,j} < BSVb{i,j}
            F{i,j} = B{i,j};
        else
            F{i,j} = (A{i,j} + B{i,j})/2;
        end
    end
end
 
F = cell2mat(F);   %将细胞数组F重新还原成矩阵（即图像）
subplot(131);imshow('clock1.jpg');
subplot(132);imshow('clock2.jpg');
subplot(133);imshow(F);

　　这里要注意一点的是要记得将之前的cell数组还原成矩阵。这个时候就会发现一个问题：error：文件名或 URL 参数必须为字符矢量。

后来找了好久，才找到原因，原来SF函数的imread不应该存在了，因为图片读取在函数外面已经完成了，直接把I =imread(img);注释掉，然后将
I=double(I);改成输入为img就好了。然后再运行融合算法就没问题了。把这个算法封装成函数吧，因为DE算法还要调用融合算法。如下：

function img_fusion = fusion(img1,img2,height,width)
 
[A,BSVa] = spilt2(img1,height,width);
[B,BSVb] = spilt2(img2,height,width);
for i =1:height
    for j = 1:width
        if BSVa{i,j} > BSVb{i,j} %比较Ai和Bi的块的清晰度，
            F{i,j} = A{i,j};     %将清晰度比较高的块放到融合块中。
        elseif BSVa{i,j} < BSVb{i,j}
            F{i,j} = B{i,j};
        else
            F{i,j} = (A{i,j} + B{i,j})/2;
        end
    end
end
 
F = cell2mat(F);   %将细胞数组F重新还原成矩阵（即图像）
img_fusion = F;
%subplot(131);imshow('clock1.jpg');
%subplot(132);imshow('clock2.jpg');
%subplot(133);imshow(F);
end

　　输入为源图像1，2，（这里暂时只考虑两张图片的融合）以及分成m*n块。输出为融合图像。效果如下：

4、下面是改良的用于图像融合的DE算法

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

function out_val=DE(img1,img2)
NP=10;   %种群规模
D=2;    %变量维数
F=0.9;  %交叉率
CR=0.6;  %变异因子
GM=50;   %迭代次数
G=1;     %初始化迭代次数
f_best=0;  %适应度
best_width=32;
best_height=32; 
count=0;
%----------区间划分---------
[row,col]=size(img1);
height_max=row;  %生成种群高的最大值
height_min=1;  %生成种群高的最小值
width_max=col;   %生成种群宽的最大值
width_min=1; %生成种群宽的最小值
Z=cell(1,NP);
for j=1:NP
    Z{j}={[(j-1)*row/NP+1:j*row/NP],[(j-1)*col/NP+1:j*col/NP]};
     
end
%---------初始化种群------
X=cell(1,NP);
for j=1:NP
    Z1=Z{j}{1};
    Z2=Z{j}{2};
    X{j}=[min(Z1)+rand(1)*(max(Z1)-min(Z1)),min(Z2)+rand(1)*(max(Z2)-min(Z2))];
end
%---------迭代------------
while G<=GM
    count=count+1;
      for j=1:NP
        if X{j}(1)>height_max
            height_max=X{j}(1);
        end
        if  X{j}(1)<height_min
            height_min=X{j}(1);
        end
        if X{j}(2)>width_max
            width_max=X{j}(2);
        end
         if  X{j}(2)<width_min
            width_min=X{j}(2);
         end
    end
%---------种群变异---------
    V=cell(1,NP); 
    i=1;
    while i<=NP
        a=randperm(NP);
        b=a(1:3);
        X1=X(b);
        V{i}=X1{1}+F*(X1{2}-X1{3});
        if(((V{i}(1)>0)&(V{i}(2)>0)&V{i}(1)<row&V{i}(2)<col))
            i=i+1;
        end
    end
%-------交叉-------
 U=cell(1,NP);   
    for j=1:NP
        a=randperm(2);
        b=a(1);
        r_rand=rand();
        for i=1:2
            if r_rand>CR&i~=b
                U{j}(i)=X{j}(i);
            else
                U{j}(i)=V{j}(i);
            end
        end
    end
    
%-------选择--------
    G=G+1;
    for j=1:NP
        X{j}=ceil(X{j});
        U{j}=ceil(U{j});
        height_1=X{j}(1);
        width_1=X{j}(2);
        img_fusion_1=fusion(img1,img2,height_1,width_1);
        height_2=U{j}(1);
        width_2=U{j}(2);
        img_fusion_2=fusion(img1,img2,height_2,width_2);
        f_new_1=spatial_frequencies(img_fusion_1);          %适应度函数为SF，计算融合图像整体的清晰度
        f_new_2=spatial_frequencies(img_fusion_2);
        if f_new_1<f_new_2
            f=f_new_2;
            current_width=width_2;
            current_height=height_2;
            X{j}=U{j};
        else
            f=f_new_1;
            current_width=width_1;
            current_height=height_1;
        end
        if f>f_best
            f_best=f;
            best_height=current_height;
            best_width=current_width;
        end
    end
 
%----------判断种群中元素是否相等--------
    val=[0,0];
    for j=1:NP
        val=val+X{j};
    end
    avg_val=val/NP;
    number=0;
    for j=1:NP
        if(X{j}(1)==avg_val(1)&&X{j}(2)==avg_val(2))
            number=number+1;
        end
    end
    if(number==NP)
        H=0;
        w_1=ceil(X{1}(1)/(row/NP));
        w_2=ceil(X{1}(2)/(col/NP));
        Z1=[(w_1-1)*row/NP+1:w_1*row/NP];
        Z2=[(w_2-1)*col/NP+1:w_2*col/NP];
        for j=2:NP
            X{j}=[min(Z1)+rand(1)*(max(Z1)-min(Z1)),min(Z2)+rand(1)*(max(Z2)-min(Z2))];
        end
    end       
end
%------最佳分块----
out_val=[height_min,height_max,width_min,width_max];

　再然后就可以直接调用DE算法计算出图像的最优块大小，

%通过DE算法求出图像的最优分块
img1 = imread('lab_1.tif');
img2 = imread('lab_2.tif');
out = DE(img1,img2);

　最后直接根据最优块大小设计融合函数

%用求得的最优分块来融合图像
img1 = imread('lab_1.tif');
img2 = imread('lab_2.tif');
img_fusion = fusion(img1,img2,60,72);
subplot(131);imshow(img1);
subplot(132);imshow(img2);
subplot(133);imshow(img_fusion);