OpenCV算法学习笔记之平滑算法
此系列的其他文章:
OpenCV算法学习笔记之初识OpenCV
OpenCV算法学习笔记之几何变换
OpenCV算法学习笔记之对比度增强
更多文章可以访问我的博客Aengus | Blog
平滑技术主要用于去除图像的噪音,常用的图像平滑算法有高斯平滑,均值平滑,中值平滑,具有保持边缘作用的平滑算法的双边滤波、导向滤波等。
二维离散卷积运算
卷积定义与full卷积
离散卷积运算,针对矩阵与卷积核矩阵
,有以下计算步骤:
第一步:将逆时针翻转180°得到
;
第二步:沿着
按照先行后列的顺序移动,每移动到一个固定的位置,
与
重合的部分就相乘求和,直到
每个位置都进行了运算。
如下,其中与
为:
第一步:
第二步:
二维卷积运算
将以上运算结果存到矩阵中,即
,该矩阵称为
与
的“full卷积”结果,用符号★表示,记
,卷积核
通常也称为卷积掩码或卷积算子。
显然高为,宽为
的矩阵
与高为
,宽为
的矩阵
full卷积后的结果是一个高为
,宽为
的矩阵。
* 矩阵计算
同样full卷积可以利用矩阵进行计算,计算过程较为复杂。有以下几个步骤:
第一步:在和
的右侧与下侧填充为零,将其尺寸扩展到
,其中
,
,扩展后矩阵记为
和
;
第二步:把按行堆叠,重构成一个
的列向量
。步骤是将
的第一行转置,使之成为
最上面的
个元素,然后将其他行转置,依次放在下面;
第三步:基于的每一行构造一个
的循环矩阵。具体步骤为:以构造基于
的第
行的循环矩阵
为例,首先将该行转置为
的第一列,然后其后的每一列都是由前一行向下移位得到的,即
构造循环矩阵
这样就构造出了个
个循环矩阵;
第四步:以第三步得到的个矩阵为块矩阵,构造块循环矩阵
,即
与第三步基于行向量构造循环矩阵类似,为一个
的矩阵;
第五步:令矩阵,这里的乘法是矩阵的乘法,最终所得
为一个
的矩阵;
第六步:将第五步得到的矩阵重新排列成一个的矩阵,即取前
个排为第一行,接着再取
个作第二行,以此类推直到所有的都取完,所得结果即为full卷积结果
。
valid卷积
由full卷积计算过程可知如果靠近
的边界。就会有部分延伸到未定义的值,忽略边界,只考虑
能完全覆盖
的情况,该过程就称为valid卷积。就最上面的例子而言,满足情况的只有
vaild卷积
因此。
高为,宽为
的矩阵
与高为
,宽为
的矩阵
valid卷积后的结果是一个高为
,宽为
的矩阵,只有
且
时才会存在valid卷积。如果存在valid卷积,显然是full卷积的一部分,用Python语法表示两者的关系为:
用OpenCV表示则为:
same卷积
以上两种卷积的结果要么比原图像大,要么比原图像小,可以利用same卷积得到与原图相同大小的结果。
same卷积的不同之处在于先选定一个“锚点”,然后用这个“锚点”分别去乘原图像的每一个位置并求和,还是用最上面的例子,选定最左上角为“锚点”,注意第一步还都是将卷积核翻转,same卷积的过程为:
same卷积过程
最终卷积结果。same卷积结果也是full卷积的一部分,假设
的锚点位置为第
行第
列(位置索引从0开始),则用Python语法表示same卷积结果与full卷积结果的关系为:
OpenCV语法表示为:
大部分情况下,为了更方便地指定卷积核的锚点,通常选定卷积核的宽高为奇数,令卷积核的中心点为锚点。
对于full卷积和same卷积,矩阵的边界处的值由于缺乏完整的邻接值,因此卷积运算常常在这些区域进行特殊处理,常用的方法为边界扩充,有以下几种方式:
- 在矩阵
的边界外填充常数,常常填充0;
 - 重复矩阵
边界处的行和列;
 - 卷绕输入矩阵,即矩阵平铺;
 - 以矩阵边界为中心,令矩阵外某位置未定义的灰度值等于图像内其镜像位置的灰度值。
 
最常用的是第4种方式,此种方式可以产生最小程度的干扰。
代码实现
边界填充
OpenCV提供函数void copyMakeBorder(InputArray src, OutputArray dst, int top, int buttom, int left, int right, int borderType, const Scalar& value=Scalar())实现对边界的扩充,其中top、bottom、left、right分别代表上方、下方、左侧、右侧填充的行数,borderType代表填充的类型,参数如下:
BORDER_REPLCATE  // 边界复制
BORDER_CONSTANT  // 常数扩充
BORDER_REFLECT   // 反射扩充
BORDER_REFLECT_101 // 以边界为中心的反射扩充
BORDER_WRAP      // 平铺扩充
value代表borderType=BORDER_CONSTANT的时候的常数值。默认的填充方式为BORDER_REFLECT_101。
二维卷积运算
Python的科学计算包Scipy提供函数convolve2d(in1, in2, mode='full', boundary='fill', fillvalue=0)实现了卷积功能,参数解释如下表所示:
| 参数 | 解释 | 
|---|---|
| in1 | 输入的二维矩阵 | 
| in2 | 输入的二维矩阵,代表卷积核 | 
| mode | 卷积类型,'full','valid','same' | 
| boundary | 边界扩充方式,'fill','wrap','symm' | 
| fillvalue | 当boundary='fill'时,边界填充常数的值,默认为0 | 
当卷积类型为same时,卷积核in2的锚点位置因其尺寸不同而不同,假设它的高宽分别为、
:
- 
和
都为奇数时,锚点位置默认为中心点
;
 - 
为偶数、
为奇数时,锚点的位置默认为
;
 - 
为奇数、
为偶数时,锚点的位置默认为
;
 - 当
和
均为偶数时,锚点位置默认在右下角
;
 
对于边界扩充,值“fill”等价于函数copyMakeBorder的BORDER_CONSTANT,值“symm”等价于BORDER_REFLECT,值“wrap”等价于BORDER_WRAP。
如果想利用convolve2d()函数实现same卷积,可以利用我们之前提到的full卷积结果和same卷积结果的关系,以下为Python示例代码:
import  numpy as np
from scipy import signal
src = np.array([[1, 2], [3, 4]], np.float32)
# src的高和宽
h1, w1 = src.shape[:2]
# 卷积核
kernel = np.array([[-1, -2], [2, 1]], np.float32)
# 卷积核的高和宽
h2, w2 = kernel.shape[:2]
# 计算 full 卷积
c_full= signal.convolve2d(src, kernel, mode='full')
# 指定锚点位置
kr, kc = 0, 0
# 根据锚点位置,从full卷积结果中截取得到same卷积
c_same = c_full[h2-kr-1:h1+h2-kr-1, w2-kc-1:w1+w2-kc-1]
C++实现
OpenCV中没有直接给出卷积运算的函数,但是我们可以利用两个函数实现。对于卷积运算的第一步,使用函数flip使得输入的卷积核逆时针逆转180°;然后通过函数void filter2D(InputArray src, OutputArray dst, int ddepth, InputArray kernel, Point anchor=Point(-1, -1), double delta=0, int borderType=BORDER_DEFAULT)实现第二步,参数解释如下表所示:
| 参数 | 解释 | 
|---|---|
| src | 输入矩阵 | 
| dst | 输出矩阵 | 
| ddepth | 输出矩阵的数据类型(位深) | 
| kernel | 卷积核,数据类型为CV_32F/CV_64F | 
| anchor | 锚点位置 | 
| delta | 默认值为0 | 
| borderType | 边界扩充类型 | 
需要注意的是输入矩阵和输出矩阵的数据类型的对应:
当src.depth()=CV_8U时,ddepth=-1/CV_16S/CV_32F/CV_64F;
当src.depth()=CV_16U/CV_16S时,ddepth=-1/CV_32F/CV_64F;
当src.depth()=CV_32F时,ddepth=-1/CV_32F/CV_64F;
当src.depth()=CV_64F时,ddepth=-1/CV_64F;
当参数ddept=-1时,代表输出矩阵和输入矩阵的数据类型一致,而卷积核kernel的数据类型必须为CV_32F/CV_64F,否则即使程序不报错,结果也可能出错;
// 翻转后的卷积核
Mat kernelFlip;
flip(kernel, kernelFlip, -1);
// 卷积运算的第二步
filter2D(src, dst, ddepth, kernelFlip, anchor, 0.0, borderType);
参数anchor指定锚点位置,Point(-1, -1)代表默认值,也就是我们上面提到的Python的科学计算包Scipy提供函数convolve2d的默认值。
可分离的卷积核
如果一个卷积核可以由至少两个尺寸比它小的卷积核full卷积而成,并且在计算过程中在所有边界处均进行扩充零的操作,且满足
其中 的尺寸均比
小,则称该卷积核是可分离的。
在图像处理中使用的卷积核通常可以分离为一维水平方向和一维垂直方向上的卷积核,如:
主要注意的是full卷积是不满足交换律的,但是一维水平方向和一维垂直方向的卷积核的full卷积是满足交换律的。
离散卷积的性质
full卷积的性质
如果卷积核是可分离的,且
,则有:
same卷积的性质
我们这里用符号表示same卷积运算,设卷积核
是full卷积可分离的,即
,且
的行和列都为奇数,这时我们取锚点为中心点,对边界进行扩充,则same卷积有类似性质:
如果对边界扩充的常数不为0,则得到的卷积结果可能是不同的,但是只有上下左右的边界处不相同;一般情况下,用非零常数填充边界,若的尺寸是
,那么结果只有矩阵的上侧的前
行,左侧的前
,下侧的最后
行,右侧的最后
行不同,其他中间部分相同。不同的填充方式最终结果也只有边界处不同。
利用same卷积的结合律可以大大减小运算量,这也是图像处理中利用可分离卷积核进行计算的重要原因。
高斯平滑
高斯卷积核的构建及分离性
构造宽为,高为
的高斯卷积核算子
,其中高宽均为奇数,锚点位置在
,步骤如下:
第一步:计算高斯矩阵:
其中
代表位置索引,其中
,且均为整数;
第二步:计算高斯矩阵的和:
第三步:高斯矩阵除以其本身的和,即归一化,得到的便是高斯卷积算子:
我们用Python实现以上步骤:
def get_gauss_kernel(sigma, h, w):
    # 第一步:构建高斯矩阵
    gauss_matrix = np.zeros([h, w], np.float32)
    # 得到中心点位置
    c_h = (h - 1)/2
    c_w = (w - 1)/2
    # 计算gauss(sigma, r, c)
    for r in range(h):
        for c in range(w):
            norm2 = math.pow(r - c_h, 2) + math.pow(c - c_w, 2)
            gauss_matrix[r][c] = math.exp(-norm2/(2*math.pow(sigma, 2)))
    # 第二步:计算高斯矩阵和
    gauss_sum = np.sum(gauss_matrix)
    # 第三步:归一化
    gauss_kernel = gauss_matrix/gauss_sum
    return gauss_kernel
高斯卷积算子翻转180度后和原来的是相同的。
可分离性
因为,所以高斯卷积核可分离成一维水平方向上的高斯核和一维垂直方向上的高斯核。基于这种可分离性,OpenCV只给出了构建一维垂直方向上的高斯卷积核函数
Mat getGaussianKernel(int ksize, double sigma, int ktype=CV_64F),其中ksize代表高斯核的行数,为正奇数;sigma为标准差;ktype决定输出的数据类型为CV_32F或CV_64F。
如果我们需要水平方向的高斯核,只需将结果进行转置就可以了,转置函数为transpose(src)。
高斯核的二项式近似
二项式函数
一维高斯函数为
其中期望,方差
。
时,高斯卷积算子和二项式函数有如下关系:
因为一维高斯卷积核的尺寸是奇数,那么对应二项式的指数为偶数,常用的二项式近似如下:
| n | 窗口大小 | |||
|---|---|---|---|---|
| 2 | 1 | 0.5 | 3×3 | 1 2 1 | 
| 4 | 2 | 1 | 5×5 | 1 4 6 4 1 | 
| 6 | 3 | 1.5 | 7×7 | 1 6 15 20 15 6 1 | 
| 8 | 4 | 2 | 9×9 | 1 8 28 56 70 56 28 8 1 | 
实现
Python
定义gauss_blur实现高斯平滑,首先得到水平方向上的高斯卷积,然后进行垂直方向的高斯卷积,其中sigma、h、w分别代表高斯核的标准差,高,宽,后面两个参数同signal.convolve2d,代表矩阵的边界扩充方式,常用的方式为symm,具体代码如下:
def gauss-blur(image, sigma, h, w, _boundary='fill', _fillvalue=0):
    # 构建水平方向上的卷积核
    gauss_kernel_x = cv2.getGaussianKernel(sigma, w, cv2.CV_64F)
    # 转置
    gauss_kernel_x = np.transpoose(gauss_kernel_x)
    # 图像与水平高斯卷积核卷积
    gauss_blur_x = signal.convlove2d(image, gauss_kernel_x, mode='same', 
                   boundary=_boundary, fillvalue=_fillvalue)
    # 进行垂直方向上的卷积操作
    gauss_kernel_y = cv2.getGaussianKernel(sigma, h, cv2.CV_64F)
    gauss_blur_xy = signal.convlove2d(gauss_blur_x, gauss_kernel_y, mode='same', 
                    boundary=_boundary, fillvalue=_fillvalue)
    
    return gauss_blur_xy 
上面的代码返回的结果是浮点类型,需要用命令astype(numpy.uint8)进行数据类型转换,从而进行灰度显示,否则用imshow会显示为黑色。
随着高斯核尺寸和标准差的增大,图像平滑效果越来越明显,图像越来越模糊。
C++
由于高斯卷积核是可分离的,所以可以通过定义可分离卷积函数sepConv2D_X_Y实现图像高斯平滑。
Mat gaussBlur(const Mat &image, Size winSize, float sigma, int ddepth=CV_64F, Point anchor=Point(-1, -1), int borderType=BORDER_DEFAULT){
    // 卷积核的宽高,均为奇数
    CV_Assert(winSize.width%2 == 1 && winSize.height%2 == 1);  // 断言
    // 构建垂直方向上的卷积核算子
    Mat gaussKernelY = getGaussianKernel(sigma, winSize.height, CV_64F);
    // 构建水平方向上的卷积核算子
    Mat gaussKernelX = getGaussianKernel(sigma, winSize.width, CV_64F);
    gaussKernelX = gaussKernelX.t();  // 转置
    // 分离的高斯卷积核
    Mat blurImage;
    sepConv2D_X_Y(image, blurImage, ddepth, gaussKernelY, gaussKernelX, Point(-1, -1));
    
    return blurImage;
}
OpenCV提供函数void GaussianBlur(InputArray src, OutputArray dst, Size ksize, double sigmaX, double sigmaY=0, int borderType=BORDER_DEFAULT),该函数可以处理灰度或彩色图像,对彩色图像的处理是对每个通道分别进行处理;ksize是高斯卷积核尺寸,高宽均为奇数且可以不同;sigmaX, sigmaY是水平和垂直方向上的标准差,sigmaY默认为0,代表和sigmaX相同。
可以看出上述函数也是通过卷积核分离实现的。当卷积核较小时结果对标准差的变化不是很敏感;当卷积核尺寸较大时结果对标准差的变化就比较敏感。
均值平滑
均值平滑顾名思义就是图像的每一个位置的邻域的平均值作为该位置的输出值。代码实现也比较简单,将高斯算子代码换为均值算子代码即可。
均值卷积核的构建
高为,宽为
的均值卷积算子的构建方法很简单,令所有元素均为
即可:
其中高和宽均为奇数,锚点位置在。
均值卷积算子也是可分离的。
快速均值平滑
由于卷积运算是比较耗时的计算,所以我们可以根据均值平滑的定义利用图像的积分进行快速均值平滑操作。
行
列图像
的积分
定义为
也就是任意一个位置的积分等于最左上角到该位置的矩阵的和。举例如下:
图像积分
利用矩阵积分,可以计算出矩阵中任意矩形区域的和:
均值平滑的本质就是计算任意一个点的领域的平均值,而平均值是由该领域的和除以邻域的面积得到的。
图像积分实现
为了在快速均值平滑中省去判断边界的问题,所以对积分后的图像矩阵的上方和左侧进行补零,尺寸为。代码如下:
def integral(image):
    rows, cols = image.shape
    # 行积分运算
    inte_image_c = np.zeros((rows, cols), np.float32)
    for r in range(rows):
        for c in range(cols):
            if c == 0:
                inte_image_c[r][c] = image[r][c]
            else:
                inte_image_r[r][c] = inte_image_c[r][c-1] + image[r][c]
    # 列积分运算
    inte_image = np.zeros(image.shape, np.float32)
    for c in range(cols):
        for r in range(rows):
            if r == 0:
                inte_image[r][c] = inte_image_c[r][c]
            else:
                inte_image[r][c] = inte_image[r-1][c] + inte_image_c[r][c]
    # 补零
    inte_image0 = np.zeros((rows+1, cols+1), np.float32)
    inte_image0[1:rows+1, 1:cols+1] = imte_image
    
    return inte_image0
OpenCV提供函数void integral(InputArray src, OutputArray sum, int sdept=-1)实现了图像积分,其中sum是输出矩阵,大小为,
sdepth代表输出矩阵的数据类型,有CV_32S/CV_32F/CV_64F。该函数也会在最上面一行和最左列进行补零操作。
Python实现
下面我们定义函数fast_mean_blur实现均值平滑:
def fast_mean_blur(image, size, border_type=cv2.BORDER_DEFAULT):
    """
    :param image: 图像矩阵
    :param size: 平滑算子尺寸
    :param border_type: 边界扩充类型
    :return: 均值平滑后的图像
    """
    half_h = (size[0]-1)/2
    half_w = (size[1]-1)/2
    # 算子矩阵面积的倒数
    ratio = 1.0/(size[0]*size[1])
    # 边界扩充
    padd_image = cv2.copyMakeBorder(image, half_h, half_h, half_w, half_w, border_type)
    # 图像积分
    padd_integral = integral(padd_image)
    # 图像的高宽
    rows, cols = image.shape
    # 均值平滑后的结果
    mean_blur_image = np.zeros(image.shape, np.float32)
    r, c = 0, 0
    for h in range(half_h, half_h+rows, 1):
        for w in range(half_w, half_w+cols, 1):
            mean_blur_image[r][c] = (padd_integral[h+half_h+1][w+half_w+1] + 
                                     padd_integral[h-half_h][w-half_w] - 
                                     padd_integral[h+half_h+1][w-half_w] - 
                                     padd_integral[h-half_h][w+half_w+1]) * ratio
            c += 1
        r += 1
        c = 0
    
    return mean_blur_image
函数返回的是浮点型,如果输入的是8位图,还需要利用命令astype(numpy.uint8)进行转换。
C++实现
C++实现代码和Python实现的思路类似,这里就不给出C++代码了。
OpenCV提供函数boxFilter和blur均实现了快速均值平滑,且均可以处理多通道图像。
void boxFilter( InputArray src, OutputArray dst, int ddepth,
                             Size ksize, Point anchor = Point(-1,-1),
                             bool normalize = true,
                             int borderType = BORDER_DEFAULT );
void blur( InputArray src, OutputArray dst,
                        Size ksize, Point anchor = Point(-1,-1),
                        int borderType = BORDER_DEFAULT );
对于blur函数中的anchor,若输入图像矩阵高宽均为奇数,则Point(-1, -1)代表中心点。boxFilter(src, dst, src.type(), ksize, Point(-1, -1), true, BORDER_DEFAULT)的功能和blur是一样的。
中值平滑
中值平滑原理就是把某个位置的邻域的灰度值进行排序,然后取中位数作为这个位置的输出值。对于中值平滑较为合适的边界处理方式为镜像补充。
中值平滑最重要的作用是去除椒盐噪声。椒盐噪声是指在图像传输的过程中由于解码错误等原因出现孤立白点或黑点的现象。除此之外,中值平滑还有一个特别重要的特性就是不会降低边缘的锐利程度,也就是说具有一定的保边作用。
一般来说,如果图像中出现较亮或者较暗的物体,若其大小小于中值平滑的窗口半径,那么它们基本会被过滤掉,而较大的目标则几乎被原封不动的保存下来。
实现
Python
首先利用命令ndarray[r1:r2+1, c1:c2+1]得到ndarray从左上角(r1, c1)到右下角(r2, c2)的矩形区域,然后利用Numpy函数median取该区域的中数。
def median_blur(image, size):
    # 图像的高宽
    rows, cols = image.shape
    # 窗口的高宽,均为奇数
    win_h, win_w = size
    half_win_h = (win_h-1)/2
    half_win_w = (win_w-1)/2
    # 中值滤波后的输出图像
    median_blur_image = np.zeros(image.shape, image.dtype)
    for r in range(rows):
        for c in range(cols):
            # 判断边界
            r_top = 0 if r-half_win_h < 0 else r-half_win_h  # 镜像扩充
            r_bottom = rows-1 if r+half_win_h > rows-1 else r+half_win_h
            c_left = 0 if r-half_win_w < 0 else c-half_win_w
            c_right = cols-1 if c+half_win_w > cols-1 else c+half_win_w
            # 取领域
            region = image[r_top:r_bottom+1, c_left:c_right+1]
            # 求中值
            median_blur_image[r][c] = np.median(region)
    return median_blur_image
C++
OpenCV没有提供直接计算中数的函数,但是我们可以利用sort函数和reshape函数间接计算中数。sort(InputArray src, OutputArray dst, int flags)的参数解释如下表:
| 参数 | 解释 | 
|---|---|
| src | 输入的单通道矩阵 | 
| dst | 输出矩阵,大小和数据类型与src相同 | 
| flags | SORT_EVERY_ROW:对每一行排序;SORT_EVERY_COLUMN:对每一列排序;SORT_ASCENDING:升序;SORT_DESCENDING:降序 | 
参数flags是可以组合使用的,如SORT_EVERY_ROW+SORT_ASCENDING代表对每行升序排序。
reshape(channels, rows)函数将矩阵变成一行或者一列,然后取中间位置得到中数。
中值平滑要对所有像素点的邻域进行排序,所以一般比卷积慢。
OpenCV提供函数medianBlur(src, dst, int ksize)实现了中值平滑,其中参数ksize代表中值算子尺寸。
双边滤波和导向滤波
双向滤波和导向滤波都有保持边缘的特点。双边滤波是根据每个位置的邻域对该位置构建不同的权重模板,OpenCV提供函数bilateralFilter和adaptiveBilateralFilter实现了双边滤波;而导向滤波是Kaiming He提出的,速度比双边滤波快,效果也更好。作者在文献中给出了较为详细的伪代码。
K. He, J. Sun, and X. Tang. Guided image filtering. In ECCV, pages 1-14. 2010
K. He, J. Sun, and X. Tang. Guided image filtering. TPAMI, 35(6): 1397-1409, 2013
参考
《OpenCV算法精解——基于Python和C++》(张平)第五章