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

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

    首先，我们来复习一下用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循环内加上如下代码：

   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)
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算法

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);

　事实上，最后的结果并不是很好，不但运行的时间非常长，而且融合后的图像出现了明显的块效应（如下图），有可能是我设计的融合算法有问题，因为分块的时候涉及到了对小数的取整问题，融合后就会丢失一些像素。

 

转载于:https://www.cnblogs.com/libai123456/p/9157370.html