差分演化算法實現多聚焦圖像融合論文算法實現（算法實現）

本文為上一篇文章的算法實現。

    首先，我們來復習一下用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);

　事實上，最后的結果並不是很好，不但運行的時間非常長，而且融合后的圖像出現了明顯的塊效應（如下圖），有可能是我設計的融合算法有問題，因為分塊的時候涉及到了對小數的取整問題，融合后就會丟失一些像素。