利用拉普拉斯对高频信息的保留，重建出原始图像

操作过程如下

对于两张图像

可以看到两张同样的图片一个脸部模糊一个除了脸其他地方模糊，现在我们需要通过局部能量匹配的融合方法对图像进行增强,即基于能量特征重建

通过一个简单的二阶差分可以看到两幅图边缘差异，即一个脸部模糊一个除了脸其他地方模糊

所谓能量就是二阶差分的边缘像素强度，融合过程中对不同尺度上的两幅图提取能量最强的区域进行融合。

能量匹配过程如下步骤：

若该点匹配度M

若该点匹配度M>e,则根据能量大小分配权重，能量小的权重为

能量大的权重为

融合代码为

E1 = imfilter(F11.^2, f, 'replicate');%计算能量

E2 = imfilter(F12.^2, f, 'replicate');

M = imfilter(F11.*F12, f, 'replicate').^2./(E1.*E2);%计算匹配度

F11(E1

W_min = 0.5*(1 - (1-M)./(1-e));%计算权重

W_max = 1 - W_min;

F21 = W_min.*pyr1{i}+W_max.*pyr2{i};

F22 = W_min.*pyr2{i}+W_max.*pyr1{i};

F21(E1>E2) = F22(E1>E2);%匹配度高时

F11(M>e) = F21(M>e);%按匹配度计算

试验结果：

下面给出完整matlab代码

%局部能量特征

clc

clear

close all

im1 = im2double(imread('lena1.jpg'));

im2 = im2double(imread('lena2.jpg'));

im1=im1(:,:,1);

im2=im2(:,:,1);

level = 4;%金字塔层数

sigma = 1;%金字塔每层提取的细节尺度。小了可能导致细节提取不充分，大了可能会导致融合结果失真。

im1_pyr = MLPD(im1, level, sigma);

im2_pyr = MLPD(im2, level, sigma);

pyr = pyr_fusion(im1_pyr, im2_pyr);

out = pyr{level};

for i = level-1 : -1 : 1

out = pyr{i} + imresize(out,[size(pyr{i},1) size(pyr{i},2)],'bilinear');

end

imshow([im1 im2 out])

%------------------------------函数1 金字塔分解------------------------------

function pyr = MLPD(I, nlev, sigma)

%Laplacian Pyramid Decomposition

%多通道拉普拉斯金字塔分解

% nlev 金字塔层数

% sigma 高斯模糊程度

if ~exist('nlev', 'var')

nlev = 4;

end

if ~exist('sigma', 'var')

sigma = 1;

end

k = fspecial('gaussian', floor(sigma*3)*2+1, sigma);

%构建拉普拉斯金字塔

pyr = cell(nlev,1);

J = I;

for j = 1:nlev-1

J_blur = imfilter(J, k, 'replicate');

J_blur_down = J_blur(1:2:size(J_blur,1)-1,1:2:size(J_blur,2)-1, :); %downsample

J_blur_high = imresize(J_blur_down,[size(J_blur,1) size(J_blur,2)],'bilinear');

pyr{j} = J-J_blur_high;

J=J_blur_down;

end

pyr{nlev}=J_blur_down; %最上一层即为高斯金字塔

end

%------------------------------函数2 金字塔融合------------------------------

function pyr = pyr_fusion(pyr1,pyr2, e)

%pyramid fusion

%局部能量特征的拉普拉斯金字塔融合

if ~exist('e', 'var')

e = 0.6;

end

n = size(pyr1, 1);

pyr = cell(n, 1);

f = fspecial('average', [3 3]);

% 对每一层融合

for i=n:-1:1

F11 = pyr1{i};

F12 = pyr2{i};

E1 = imfilter(F11.^2, f, 'replicate');%计算能量

E2 = imfilter(F12.^2, f, 'replicate');

%figure()

%imshow(E1)

M = imfilter(F11.*F12, f, 'replicate').^2./(E1.*E2);%计算匹配度

F11(E1

W_min = 0.5*(1 - (1-M)./(1-e));%计算权重

W_max = 1 - W_min;

F21 = W_min.*pyr1{i}+W_max.*pyr2{i};

F22 = W_min.*pyr2{i}+W_max.*pyr1{i};

F21(E1>E2) = F22(E1>E2);%匹配度高时

F11(M>e) = F21(M>e);%按匹配度计算

figure()

imshow((F11-min(F11))./(max(F11)-min(F11)))

pyr{i} = F11;

end

end

下面为完整python代码

import cv2,numpy as np

import matplotlib.pyplot as plt

莱娜1=cv2.imread('lena1.jpg')

莱娜2=cv2.imread('lena2.jpg')

莱娜1=cv2.cvtColor(莱娜1,cv2.COLOR_BGR2RGB)

莱娜2=cv2.cvtColor(莱娜2,cv2.COLOR_BGR2RGB)

plt.subplot(121)

plt.imshow(莱娜1)

plt.subplot(122)

plt.imshow(莱娜2)

莱娜1=cv2.resize(莱娜1,(256,256),interpolation=cv2.INTER_CUBIC)

莱娜2=cv2.resize(莱娜2,(256,256),interpolation=cv2.INTER_CUBIC)

层数=8

temp=莱娜1.copy()

for index in range(层数):

dst = cv2.pyrDown(temp)

def 创建拉普拉斯金字塔(path,level):

img=cv2.imread(path)

img=cv2.cvtColor(img,cv2.COLOR_RGB2BGR)

img= cv2.resize(img,(512,512),interpolation=cv2.INTER_CUBIC)

temp = img.copy()#拷贝图像

高斯金字塔 = []

高斯金字塔.append(img)

for i in range(1,level):

下采样 = cv2.pyrDown(temp)

高斯金字塔.append(下采样)

temp = 下采样.copy()

拉普拉斯金字塔 = [高斯金字塔[level-1]]

for i in range(level-1,0,-1):

上采样 = cv2.pyrUp(高斯金字塔[i])

第L层的高斯金字塔 = cv2.subtract(高斯金字塔[i-1],上采样)

拉普拉斯金字塔.append(第L层的高斯金字塔)

print('高斯金字塔层数:',len(高斯金字塔),'拉普拉斯金字塔:',len(拉普拉斯金字塔))

return 拉普拉斯金字塔

def 拉普拉斯卷积(原图,卷积核大小):

卷积核=np.ones([卷积核大小,卷积核大小])*1/(卷积核大小*卷积核大小)

卷积图 = np.zeros([原图.shape[0]-(卷积核大小-卷积核大小//2),原图.shape[1]-(卷积核大小-卷积核大小//2)])

print(原图.shape)

原图行数,原图列数=原图.shape

for y in range(卷积核大小//2,原图行数-(卷积核大小-卷积核大小//2)):

for x in range(卷积核大小//2,原图列数-(卷积核大小-卷积核大小//2)):

半径=卷积核大小//2

拉普拉斯核 = np.array(([-1,-1,-1],[-1,8,-1],[-1,-1,-1]))

窗口=np.sum(原图[x-半径:x+半径+1,y-半径:y+半径+1]*拉普拉斯核)

卷积图[x, y] = 窗口

return 卷积图

莱娜1灰度=cv2.cvtColor(莱娜1,cv2.COLOR_RGB2GRAY)

莱娜1卷积图=拉普拉斯卷积(莱娜1灰度,3)

莱娜2灰度=cv2.cvtColor(莱娜2,cv2.COLOR_RGB2GRAY)

莱娜2卷积图=拉普拉斯卷积(莱娜2灰度,3)

plt.subplot(121)

plt.imshow(莱娜1卷积图,'gray')

plt.subplot(122)

plt.imshow(莱娜2卷积图,'gray')

def 平均值卷积(原图,卷积核大小,通道):

def 单通道卷积(原图,卷积核大小):

原图=cv2.copyMakeBorder(原图,卷积核大小//2,卷积核大小//2,卷积核大小//2,卷积核大小//2,cv2.BORDER_CONSTANT,0)

卷积核=np.ones([卷积核大小,卷积核大小])*1/(卷积核大小*卷积核大小)

卷积图 = np.zeros([原图.shape[0]-(卷积核大小-卷积核大小//2),原图.shape[1]-(卷积核大小-卷积核大小//2)])

原图行数,原图列数=原图.shape

for y in range(卷积核大小//2,原图行数-(卷积核大小-卷积核大小//2)):

for x in range(卷积核大小//2,原图列数-(卷积核大小-卷积核大小//2)):

半径=卷积核大小//2

窗口=原图[x-半径:x+半径+1,y-半径:y+半径+1]

卷积图[x-1, y-1] = np.mean(窗口)

return 卷积图

if 通道==3:

for chanel in range(通道):

原图[:,:,chanel]=单通道卷积(原图[:,:,chanel],卷积核大小)

return 原图

return 单通道卷积(原图,卷积核大小)

层数=3

莱娜=cv2.cvtColor(cv2.imread('lena.jpg'),cv2.COLOR_RGB2BGR)

拉普拉斯金字塔1=创建拉普拉斯金字塔('lena1.jpg',层数)

拉普拉斯金字塔2=创建拉普拉斯金字塔('lena2.jpg',层数)

莱娜灰度1=cv2.cvtColor(cv2.imread('lena1.jpg'),cv2.COLOR_BGR2GRAY)

莱娜灰度2=cv2.cvtColor(cv2.imread('lena2.jpg'),cv2.COLOR_BGR2GRAY)

莱娜1=cv2.cvtColor(cv2.imread('lena1.jpg'),cv2.COLOR_RGB2BGR)

莱娜2=cv2.cvtColor(cv2.imread('lena2.jpg'),cv2.COLOR_RGB2BGR)

莱娜灰度卷积=平均值卷积(莱娜灰度1,3,1)

print(莱娜灰度卷积.shape)

plt.figure('莱娜灰度卷积')

plt.imshow(莱娜灰度卷积,'gray')

能量复原算法

拉普拉斯金字塔4=[]

通道数=3

e=0.6

for i in range(0,层数):

F=np.zeros(拉普拉斯金字塔1[i].shape)

for chanel in range(通道数):

F12=拉普拉斯金字塔1[i][:,:,chanel]

F11=拉普拉斯金字塔2[i][:,:,chanel]

E1 = 平均值卷积(F11,3,1)

E2 = 平均值卷积(F12,3,1)

M = np.power(平均值卷积(np.dot(F11,F12),3,1),2)/np.dot(E1,E2)

F11[E1

W_min = 0.5*(1 - (1-M)/(1-e))

W_max = 1 - W_min

print(拉普拉斯金字塔1[i].shape,拉普拉斯金字塔2[i].shape)

# 两图分别乘以权值

F21 = np.dot(W_min,拉普拉斯金字塔1[i][:,:,chanel])+np.dot(W_max,拉普拉斯金字塔2[i][:,:,chanel])

F22 = np.dot(W_min,拉普拉斯金字塔2[i][:,:,chanel])+np.dot(W_max,拉普拉斯金字塔1[i][:,:,chanel])

#匹配度高时

F21[E1>E2] = F22[E1>E2]

F11[M>e] = F21[M>e] #按匹配度计算

print(F11.shape,F12.shape,E1.shape,E2.shape)

F[:,:,chanel]=F11

F=(F-np.min(F))/(np.max(F)-np.min(F))

plt.figure()

plt.imshow(F)

拉普拉斯金字塔4.append(F)

#融合

ls_ = 拉普拉斯金字塔4[0]

for i in range(1,层数-1):

ls_ = cv2.pyrUp(ls_)

ls_ = cv2.add(ls_, 拉普拉斯金字塔4[i])

python结果

注：这里我用的是jupyter，所以函数定义并没那么集中，以及为了方便明确算法步骤，命名我故意使用的中文，至于为什么融合结果python和matlab有那么大差异我也不太清楚