图像金字塔是图像处理和计算机视觉中常用到的概念，常常用于多尺度处理领域(multiscale processing)，尤其早年的图像匹配、识别等算法中都用到了图像金字塔。

高斯金字塔(Gaussian pyramid)

下图为高斯金字塔的示意图，金字塔的底层为原始图像，每向上一层则是通过高斯滤波和1/2采样得到(去掉偶数行和列)。

高斯金字塔

我们可以使用如下Matlab代码来进行得到高斯金字塔：

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

下图就是生成的金字塔图像

金字塔图像

高斯滤波器可以看做一个低通滤波器，那么每经过一次的高斯滤波，图像中仅能够保留某个频率值以下的频率部分，所以高斯金字塔也可以看做一个低通金字塔(每一级只保留某个频率以下的成分)。

拉普拉斯金字塔(Laplacian pyramid)

在进行高斯金字塔运算时，由于不断的进行高斯滤波和下采样，我们丢失了很多高频信号，而拉普拉斯金字塔的目的就是保存这些高频信号，保存这些高频信号所采用的方式就是保存差分图像。比如，拉普拉斯金字塔的第０层，就是原始图像和原始图像下采样(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

%下采样函数

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, :);

％上采样函数

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，注意这里要乘以４！

% 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), :);

％产生拉普拉斯滤波器

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

重建的图像

需要注意的地方

为什么在处理高斯金字塔的时候需要采用滤波呢，直接下采样不可以吗？

如果把图像看做频率信号的话，直接进行下采样则会出现采样不足的情况，消除这种情况的方法是采用低通滤波器(高斯滤波器)对图像进行滤波，将采样不足的高频信号过滤掉，这样在进行下采样的时候就保证了不出现采样不足的情况。

一般在图像处理中，将上面Matlab实现的下采样函数(包括高斯滤波和图像尺寸减半)这一部分叫做Reduce，将上面Matlab实现的上采样函数(包括高斯滤波和图像尺寸增加一倍)这部分叫做Expand，如果用数学方式表达的话，Expand函数如下：

注意前面需要乘以４。

拉普拉斯金字塔可以看做一个带通滤波器，在每一级都保留了图像某个频率值附近的成分。(这一点与高斯金字塔不同，高斯金字塔是低通金字塔)

两个低通滤波器的差值就构成了一个带通滤波器。

两个低通之差，构成带通滤波器

Python实现

下面是高斯金字塔和拉普拉斯金字塔的Opencv-Python实现

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