圖像金字塔是圖像處理和計(jì)算機(jī)視覺中常用到的概念嵌巷,常常用于多尺度處理領(lǐng)域(multiscale processing)重窟,尤其早年的圖像匹配、識(shí)別等算法中都用到了圖像金字塔享郊。
高斯金字塔(Gaussian pyramid)
下圖為高斯金字塔的示意圖览祖,金字塔的底層為原始圖像,每向上一層則是通過高斯濾波和1/2采樣得到(去掉偶數(shù)行和列)炊琉。
我們可以使用如下Matlab代碼來進(jìn)行得到高斯金字塔:
function [ pyr ] = gaussian_pyramid( I,nlev )
%GAUSSIAN_PYRAMID Summary of this function goes here
% Detailed explanation goes here
pyr = cell(nlev,1);
pyr{1} = I;
filter = fspecial('gaussian');
for i=2:nlev
% gaussian filter
I = imfilter(I,filter,'symmetric');
% downsample
I = I(1:2:end,1:2:end);
pyr{i} = I;
end
end
下圖就是生成的金字塔圖像
高斯濾波器可以看做一個(gè)低通濾波器展蒂,那么每經(jīng)過一次的高斯濾波,圖像中僅能夠保留某個(gè)頻率值以下的頻率部分苔咪,所以高斯金字塔也可以看做一個(gè)低通金字塔(每一級只保留某個(gè)頻率以下的成分)锰悼。
拉普拉斯金字塔(Laplacian pyramid)
在進(jìn)行高斯金字塔運(yùn)算時(shí),由于不斷的進(jìn)行高斯濾波和下采樣团赏,我們丟失了很多高頻信號箕般,而拉普拉斯金字塔的目的就是保存這些高頻信號,保存這些高頻信號所采用的方式就是保存差分圖像舔清。比如丝里,拉普拉斯金字塔的第0層,就是原始圖像和原始圖像下采樣(Reduce)后再次上采樣(Expand)的圖像的差值鸠踪。
function pyr = laplacian_pyramid(I,nlev)
r = size(I,1);
c = size(I,2);
if ~exist('nlev')
% compute the highest possible pyramid
nlev = floor(log(min(r,c)) / log(2));
end
% recursively build pyramid
pyr = cell(nlev,1);
filter = pyramid_filter;
J = I;
for l = 1:nlev - 1
% apply low pass filter, and downsample
I = downsample(J,filter);
odd = 2*size(I) - size(J); % for each dimension, check if the upsampled version has to be odd
% in each level, store difference between image and upsampled low pass version
pyr{l} = J - upsample(I,odd,filter);
J = I; % continue with low pass image
end
pyr{nlev} = J; % the coarest level contains the residual low pass image
%下采樣函數(shù)
function R = downsample(I, filter)
border_mode = 'symmetric';
% low pass, convolve with separable filter
R = imfilter(I,filter,border_mode); %horizontal
R = imfilter(R,filter',border_mode); %vertical
% decimate
r = size(I,1);
c = size(I,2);
R = R(1:2:r, 1:2:c, :);
%上采樣函數(shù)
function R = upsample(I,odd,filter)
% increase resolution
I = padarray(I,[1 1 0],'replicate'); % pad the image with a 1-pixel border
r = 2*size(I,1);
c = 2*size(I,2);
k = size(I,3);
R = zeros(r,c,k);
R(1:2:r, 1:2:c, :) = 4*double(I); % increase size 2 times; the padding is now 2 pixels wide丙者,注意這里要乘以4!
% interpolate, convolve with separable filter
R = imfilter(R,filter); %horizontal
R = imfilter(R,filter'); %vertical
% remove the border
R = R(3:r - 2 - odd(1), 3:c - 2 - odd(2), :);
%產(chǎn)生拉普拉斯濾波器
function f = pyramid_filter()
f = [.05, .25, .4, .25, .05]; % original [Burt and Adelson, 1983]
%f = [.0625, .25, .375, .25, .0625]; % binom-5
f = f'*f;
end
通過上面的代碼营密,我們可以得到拉普拉斯金字塔如下所示械媒。
由于拉普拉斯金字塔保留了高頻信號,那么我們可以用它來重建原始圖像。
function R = reconstruct_laplacian_pyramid(pyr)
r = size(pyr{1},1);
c = size(pyr{1},2);
nlev = length(pyr);
% start with low pass residual
R = pyr{nlev};
filter = pyramid_filter;
for l = nlev - 1 : -1 : 1
% upsample, and add to current level
odd = 2*size(R) - size(pyr{l});
R = pyr{l} + upsample(R,odd,filter);
%figure
%imshow(R,[]);
%imwrite(mat2gray(R),[num2str(l),'.jpg']);
end
需要注意的地方
為什么在處理高斯金字塔的時(shí)候需要采用濾波呢纷捞,直接下采樣不可以嗎痢虹?
如果把圖像看做頻率信號的話,直接進(jìn)行下采樣則會(huì)出現(xiàn)采樣不足的情況主儡,消除這種情況的方法是采用低通濾波器(高斯濾波器)對圖像進(jìn)行濾波奖唯,將采樣不足的高頻信號過濾掉,這樣在進(jìn)行下采樣的時(shí)候就保證了不出現(xiàn)采樣不足的情況糜值。-
一般在圖像處理中丰捷,將上面Matlab實(shí)現(xiàn)的下采樣函數(shù)(包括高斯濾波和圖像尺寸減半)這一部分叫做Reduce,將上面Matlab實(shí)現(xiàn)的上采樣函數(shù)(包括高斯濾波和圖像尺寸增加一倍)這部分叫做Expand寂汇,如果用數(shù)學(xué)方式表達(dá)的話病往,Expand函數(shù)如下:
注意前面需要乘以4。
拉普拉斯金字塔可以看做一個(gè)帶通濾波器骄瓣,在每一級都保留了圖像某個(gè)頻率值附近的成分停巷。(這一點(diǎn)與高斯金字塔不同,高斯金字塔是低通金字塔)
兩個(gè)低通濾波器的差值就構(gòu)成了一個(gè)帶通濾波器榕栏。
Python實(shí)現(xiàn)
下面是高斯金字塔和拉普拉斯金字塔的Opencv-Python實(shí)現(xiàn)
import cv2
import numpy as np
def gaussian_pyr(img,lev):
img = img.astype(np.float)
g_pyr = [img]
cur_g = img;
for index in range(lev):
cur_g = cv2.pyrDown(cur_g)
g_pyr.append(cur_g)
return g_pyr
def laplacian_pyr(img,lev):
img = img.astype(np.float)
g_pyr = gaussian_pyr(img,lev)
l_pyr = []
for index in range(lev):
cur_g = g_pyr[index]
next_g = cv2.pyrUp(g_pyr[index+1])
cur_l = cv2.subtract(cur_g,next_g)
l_pyr.append(cur_l)
l_pyr.append(g_pyr[-1])
return l_pyr
def lpyr_recons(l_pyr):
lev = len(l_pyr)
cur_l = l_pyr[-1]
for index in range(lev-2,-1,-1):
print(index)
cur_l = cv2.pyrUp(cur_l)
next_l = l_pyr[index]
cur_l = cur_l + next_l
return cur_l