本文为上一篇文章的算法实现。
首先，我们来复习一下用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);
事实上，最后的结果并不是很好，不但运行的时间非常长，而且融合后的图像出现了明显的块效应（如下图），有可能是我设计的融合算法有问题，因为分块的时候涉及到了对小数的取整问题，融合后就会丢失一些像素。
