数字图像处理之【高斯金字塔】与【拉普拉斯金字塔】

1.1 什么是高斯金字塔？

高斯金字塔（Gaussian Pyramid）是一种多分辨率图像表示方法，用于图像处理和计算机视觉领域。它通过对原始图像进行一系列的高斯平滑和下采样操作，生成一组分辨率逐渐降低的图像层次结构。

高斯金字塔的构建过程通常包括以下步骤：

1. 高斯平滑（Gaussian Smoothing）：对原始图像应用高斯滤波器，生成一个平滑后的图像。高斯滤波器是一种低通滤波器，用于减少图像中的高频噪声。

2. 下采样（Downsampling）：将平滑后的图像进行下采样，通常是将图像的宽度和高度各减半，得到较低分辨率的图像。

3. 重复上述步骤：对下采样后的图像重复高斯平滑和下采样过程，直到达到预定的分辨率级别。每一个新生成的图像称为一个金字塔层。

1.2 什么是拉普拉斯金字塔？

拉普拉斯金字塔（Laplacian Pyramid）是基于高斯金字塔的多分辨率图像表示方法，用于图像处理和计算机视觉领域。它通过从高斯金字塔中提取（高频）细节信息，形成一系列细节图像的层次结构。拉普拉斯金字塔可以用于图像压缩、图像增强和图像融合等任务。

构建拉普拉斯金字塔的过程通常包括以下步骤：

1. 构建高斯金字塔：首先，通过高斯平滑和下采样操作构建高斯金字塔。

2. 生成拉普拉斯金字塔层：

   • 从高斯金字塔的每一层生成下一层（低分辨率层）后，将下一层上采样（通常是通过插值方法将图像的宽度和高度各加倍）回到当前层的分辨率。
   • 计算当前层与上采样后的图像的差值，得到当前层的拉普拉斯层。

3. 重复上述步骤：对高斯金字塔的每一层重复生成拉普拉斯层的过程，直到达到最底层。

1.3 示意图

从图中我们可以得到什么信息？拉普拉斯金字塔图像=原图的高频信息。

因此如果我们想回复原图，是不是只要用拉普拉斯金字塔图像加上上采样后的图像就能得到？我想应该是的。但是用文字描述我认为过于繁琐，我们用数学公式来进行清晰的表达。

设${\mathrm{L}}_{\mathrm{n}}$为拉普拉斯金字塔的第$\mathrm{n}$层（自底向上），${\mathrm{G}}_{\mathrm{n}}$为高斯金字塔的第$\mathrm{n}$层（自底向上），$\mathrm{GaussBlur}$为高斯平滑函数，$\mathrm{pyrDown}$为下采样函数，$\mathrm{pyrUp}$为上采样函数，则有
$$\begin{array}{rl}& {\mathrm{G}}_{\mathrm{n}+1}=\mathrm{pyrDown}(\mathrm{GaussBlur}({\mathrm{G}}_{\mathrm{n}}))\\ & {\mathrm{L}}_{\mathrm{n}}={\mathrm{G}}_{\mathrm{n}}-\mathrm{pyrUp}(\mathrm{pyrDown}(\mathrm{GaussBlur}({\mathrm{G}}_{\mathrm{n}})))\\ & {\mathrm{G}}_{\mathrm{n}}={\mathrm{L}}_{\mathrm{n}}+\mathrm{pyrUp}(\mathrm{pyrDown}(\mathrm{GaussBlur}({\mathrm{G}}_{\mathrm{n}})))\end{array}$$

1.4 OpenCV实现

对于以上代码，读者可将图片换成自己喜欢的图片进行测试；函数内部有一部分条件语句是因为如果原图宽或高非偶数，那么下采样后再上采样就会和原图大小不一样。

比如$801\times 801$的图像下采样后尺寸是$401\times 401$（如果原图尺寸是奇数，那么使用向上取整的除法），再上采样就变成了$802\times 802$和原图尺寸不匹配。

为了正确处理这种突发情况，我们需要对图像进行裁剪，刚好裁剪1行和1列。

import cv2
import numpy as np


def correctSize(imgOri, imgDownsample):
    imgUpsample = cv2.pyrUp(imgDownsample)
    if imgOri.shape[0] % 2 == 1:  # 检查原始图像宽度是否为奇数
        imgUpsample = imgUpsample[1:, :, :]  # 如果是，将放大后的图像宽度减少1
        # imgUpsample = imgUpsample[1:, :]  # 如果是，将放大后的图像宽度减少1
    if imgOri.shape[1] % 2 == 1:  # 检查原始图像宽度是否为奇数
        imgUpsample = imgUpsample[:, 1:, :]  # 如果是，将放大后的图像宽度减少1
        # imgUpsample = imgUpsample[1:, :]  # 如果是，将放大后的图像宽度减少1
    return imgOri, imgUpsample


img = cv2.imread('lena.png')
# img = cv2.imread('lena.png',cv2.IMREAD_GRAYSCALE) # 读取灰度图像

print(img.shape)
img12 = cv2.pyrDown(img)
img14 = cv2.pyrDown(img12)
img18 = cv2.pyrDown(img14)
img116 = cv2.pyrDown(img18)
ori, imgUpsample = correctSize(img, img12)
# imgRes = cv2.subtract(img, img12)
# 暂不清楚用subtract有什么隐藏机制，用该函数做减法会损失不少信息，所以就直接用原始减法好了
imgRes = ori - imgUpsample
imgRes = imgRes.clip(0, 255)
# imgRes = cv2.subtract(*correctSize(img12, img14))
# imgRes = cv2.subtract(*correctSize(img14, img18))
# imgRes = cv2.subtract(*correctSize(img18, img116))

cv2.imshow("高斯金字塔第n层图像的高频信息图-拉普拉斯", imgRes)
# cv2.imshow("高斯金字塔第n层图像的高频信息图-拉普拉斯", np.hstack((imgUpsample, imgRes)))

cv2.waitKey()
cv2.destroyAllWindows()

1.5 手搓代码(高斯平滑+下采样+上采样+拉普拉斯图像生成)

首先给出高斯平滑所用的高斯核
$$\frac{1}{256}\left[\begin{array}{ccccc}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\end{array}\right]$$

然后就是手搓代码了，纯手搓，但是运行速度有点慢，如果你图像很大的话

import cv2
import numpy as np

img = cv2.imread('lena.png')

# 定义高斯核，并归一化，用来做卷积
GaussSmoothKernel = np.array([[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]], dtype=np.float32)
GaussSmoothKernel = GaussSmoothKernel / 256

Smoothedimg = np.zeros(img.shape, dtype=np.uint8)

# 下面的单数是封装好的，可用来做图像卷积
# Smoothedimg = cv2.filter2D(img, -1, GaussSmoothKernel)
for i in range(img.shape[0]):
    for j in range(img.shape[1]):
        for k in range(3):
            if i-2 < 0 or i+2 >= img.shape[0] or j-2 < 0 or j+2 >= img.shape[1]:
                # 处理图像边边角角的模糊，在else语句中的操作没法对边角进行平滑，因为边界像素没有相邻像素，无法进行卷积操作
                for m in range(5):
                    for n in range(5):
                        if i+m-2 >= 0 and i+m-2 < img.shape[0] and j+n-2 >= 0 and j+n-2 < img.shape[1]:
                            Smoothedimg[i, j, k] += GaussSmoothKernel[m, n] * img[i+m-2, j+n-2, k]
            else:
                # 图像的边界处的像素无法通过此操作进行卷积操作
                Smoothedimg[i, j, k] = np.sum(GaussSmoothKernel * img[max(i-2, 0):min(i+3, img.shape[0]), max(j-2, 0):min(j+3, img.shape[1]), k], axis=(0, 1))

downsample = Smoothedimg[::2, ::2, :]

# 生成一个0值张量，用来存储上采样后的图像
Upsample = np.zeros((downsample.shape[0]*2, downsample.shape[1]*2, 3), dtype=np.uint8)

for i in range(downsample.shape[0]):
    for j in range(downsample.shape[1]):
        Upsample[i*2, j*2, :] = downsample[i, j, :]
        # 下面为插值操作
        Upsample[i*2+1, j*2, :] = downsample[i, j, :]
        Upsample[i*2, j*2+1, :] = downsample[i, j, :]
        Upsample[i*2+1, j*2+1, :] = downsample[i, j, :]

# 用于处理图像各方向为奇数的情况：这种情况上采样后图像会比原图大一点点，因此需要给它做个裁剪
if img.shape[0] % 2 == 1:
    Upsample = Upsample[1:, :, :]
if img.shape[1] % 2 == 1:
    Upsample = Upsample[:, 1:, :]

# 获取拉普拉斯图像
Laplacian = img - Upsample

# cv2.imshow('img', img)
# cv2.imshow('Smoothedimg', Smoothedimg)
# cv2.imshow('downsample', downsample)
# cv2.imshow('Upsample', Upsample)
# cv2.imshow('Laplacian', Laplacian)
cv2.imshow('ori+smoothedimg+downsample+Upsample+Laplacian', np.vstack((np.hstack((img, Smoothedimg)), np.hstack((Upsample, Laplacian)))))

# 该语句配合上面插值语句使用，把上面的插值语句注释，然后用下面这条语句输出，你将看到下采样丢失了多少信息
# cv2.imshow('DownsampleLossInfo', (Smoothedimg - Upsample).clip(0, 255))

cv2.waitKey(0)
cv2.destroyAllWindows()

看看效果呗，感觉手搓的效果还不赖，当然肯定没法和封装好的比，封装好的函数使用了更为复杂的插值算法：

上图从左到右边，从上到下依次为$原图\to 高斯平滑后的图\to 经过下采样后，再进行上采样恢复了分辨率的图\to 拉普拉斯图$。

1.6 一些后话

其实看了上面的东西我觉得读者也应该理解为什么高斯核被称为低通滤波器了：低频信息总是通过（保留），高频信息却被删除了。这使得图像看起来更加平滑，减少了噪声和细节。

相反，OpenCV的Sobel算子就是一个具有高通滤波器性质的算子，该算子用来做图像锐化，即边缘增强，这意味着它会增强图像的边缘和细节，而使平滑区域变暗或去除。