二、拉普拉斯金字塔簡介
1 拉普拉斯金字塔
在高斯金字塔的運算過程中，圖像經過卷積和下采樣操作會丢失部分高頻細節資訊。為描述這些高頻資訊，人們定義了拉普拉斯金字塔(Laplacian Pyramid， LP)。用高斯金字塔的每一層圖像減去其上一層圖像上采樣并高斯卷積之後的預測圖像，得到一系列的內插補點圖像即為 LP 分解圖像。
将Gl内插方法得到放大圖像Gl，使Gl的尺寸與Gl-1的尺寸相同，即放大算子Expand。
該式子實作兩個步驟：在偶數行和列插入0，然後使用下采樣中的高斯核進行濾波處理，得到和l-1層一樣大小的圖像。
N為拉普拉斯金字塔頂層的層号LPl是拉普拉斯金字塔分解的第L層圖像。由LP0，LP1、LP2…LPN構成的金字塔即為拉普拉斯金字塔。它的每一層L0圖像是高斯金字塔本層G0圖像與其高一層圖像G1經内插放大後圖像G1的差，此過程相當于帶通濾波，是以拉普拉斯金字塔又稱為帶通金字塔分解。
内插方法：opencv中有實作的函數pyrup。可以得到*G1。然後在兩個函數作差，相減就可以得到拉普拉斯金字塔。
求得每個圖像的拉普拉斯金字塔後需要對相應層次的圖像進行融合，具體的融合規則有，取大、取小，等等。
首先對原始圖像Gi進行向下取樣，然後向上采用，最後與原始圖像相減，得到拉普拉斯金字塔圖像。
下面這張圖也是比較常見的：
使用原始圖像 套入公式得到拉普拉斯金字塔第0層。
使用原始圖像向下采樣Od 代入公式 得到 拉普拉色金字塔第1層。
這張圖在看一些資料的時候也經常看到：
2 融合應用
圖像拉普拉斯金字塔分解的目的是将源圖像分别分解到不同的空間頻帶上，融合過程是在各空間頻率層上分别進行的，這樣就可以針對不同分解層的不同頻帶上的特征與細節，采用不同的融合算子以達到突出特定頻帶上特征與細節的目的。即有可能将來自不同圖像的特征與細節融合在一起。
（1）頂層處理
設LAl和LBl分别為源圖像A,B經過拉普拉斯金字塔分解後得到的第l層圖像，融合後的結果為LFl。當l=N時，LAN和LBN分别為源圖像A，B經過拉普拉斯金字塔分解後得到的頂層圖像。對于頂層圖像的融合，首先計算以其各個像素為中心的區域大小為M*N(M、N取奇數且M >= 3、N >= 3)的區域平均梯度：
其中，Ix與Iy分别為像素f(x,y)在x與y方向上的一階差分，定義如下：
是以對于頂層圖像中的每一個像素LAN(i, j)和LBN(i, j)都可以得到與之相對應的區域平均梯度GA(i, j)和GB(i, j)。由于平均梯度反映了圖像中的微小細節反差和紋理變化特征，同時也反映出圖像的清晰度。一般來說平均梯度越大，圖像層次也豐富，則圖像越清晰。是以頂層圖像的融合結果為：
（2）各層次處理
當0<l<N時，則對于經過拉普拉斯金字塔分解的第l層圖像，首先計算其區域能量：
則其他層次圖像的融合結果為：
在得到金字塔各個層次的融合圖像LF1、LF2、LFN後。通過前面的重構，便可得到最終的融合圖像。
第二種融合規則：
采用最高層系數取平均，其餘各層系數絕對值取大的融合政策進行融合。融合後圖像的系數（灰階值）越接近較清晰圖像的灰階值就說明融合效果好。
三、部分源代碼
clear all;
close all;
clc;
im1 = (imread('紅外.jpg'));
im2 = imread('可見光.jpg');
figure(1);
imshow(im1);
title('紅外圖像');
figure(2);
imshow(im2);
title('可見光圖像');
if length(size(im1)) ==3%判斷輸入圖像是否是真彩色圖像
im1 = rgb2gray(im1);
end
if length(size(im2)) ==3
im2 = rgb2gray(im2);
end
im1 = double(im1);
im2 = double(im2);
%%
%拉普拉斯濾波器
w = [1 4 6 4 1; 4 16 24 16 4; 6 24 36 24 6; 4 16 24 16 4; 1 4 6 4 1]/256;
G = cell(1,5);
H = cell(1,5);
I = cell(1,5);
G{1} = im1;%第一層為原圖像
H{1} = im2;
for k=2:5
G{k} = imfilter(G{k-1},w,'replicate');%濾波
H{k} = imfilter(H{k-1},w,'replicate');%濾波
G{k} = G{k}(1:2:size(G{k},1),1:2:size(G{k},2)); %下采樣
H{k} = H{k}(1:2:size(H{k},1),1:2:size(H{k},2)); %下采樣
end
for k=5:-1:2
G{k-1} = G{k-1}(1:2*size(G{k},1),1:2*size(G{k},2));%調整圖像尺寸
H{k-1} = H{k-1}(1:2*size(H{k},1),1:2*size(H{k},2));%調整圖像尺寸
end
for k=1:4
G{k} = G{k} - us(G{k+1},w);%獲得拉普拉斯金字塔
H{k} = H{k} - us(H{k+1},w);
a = abs(G{k}) > abs(H{k});%兩幅圖像的1-4級拉普拉斯金字塔取絕對值大的
I{k} = a.*G{k} + ~a.*H{k};
end
I{5} = H{5};%定義融合後頂層圖像尺寸，并為圖像邊緣指派，圖像中間的值會被比較平均梯度後所取的值替換
gg = g(G{5});
gh = g(H{5});
for i = 2:size(G{5},1)-2
for j = 2:size(G{5},2)-2
if gg(i,j) > gh(i,j)%比較某點鄰域平均梯度大小，取大的
I{5}(i,j) = G{5}(i,j);
else
I{5}(i,j) = H{5}(i,j);
end
end
end
for k=4:-1:1
I{k} = I{k}+us(I{k+1},w);%重構圖像
end
function f = g(b)%區域平均梯度函數（3*3區域）
[Y,X] = size(b);
g = zeros(Y,X);%儲存梯度
f = g;%儲存平均梯度
for i = 1:Y-1%求梯度，最後一行、列求不到，後面會
for j = 1:X-1
g(i,j) = sqrt((((b(i+1,j)-b(i,j))^2+(b(i,j+1)-b(i,j))^2))/2);%某點梯度
end
end
function f = us(im,w)
w=w*4;
[i,j] = size(im);
[m,n] = size(w);
w11 = w(1:2:m,1:2:n);%插入濾波器
w12 = w(1:2:m,2:2:n);
w21 = w(2:2:m,1:2:n);
w22 = w(2:2:m,2:2:n);
四、運作結果
五、matlab版本及參考文獻
1 matlab版本
2014a
2 參考文獻
[1] 蔡利梅.MATLAB圖像處理——理論、算法與執行個體分析[M].清華大學出版社，2020.
[2]楊丹,趙海濱,龍哲.MATLAB圖像處理執行個體詳解[M].清華大學出版社，2013.