图像卷积和 C 代码实现(Image Convolution)

2020/02/07 Image

图像卷积是计算机视觉中非常基础的一个问题。本文作者在 Tesla 面试时就遇到了要求实现无 bug 的快速版本(即可拆分卷积)的高斯核图像卷积的问题(氧化钙, 面试官你 45 分钟能实现并且无 bug 的话我叫你大爷)。很遗憾面试失败了,痛定思痛,在这里整理一下基础知识,并参考网上牛人的代码自己实现一遍。

参考和扩展:


一维卷积(Convolution in 1D)

一维卷积的定义比较简单。给定一个信号 x 和另一个信号 h,它们的卷积是:

可以称 h 是核函数。

注意

  • 里面的符号 * 是代表卷积的意思;
  • 虽然 k 的范围是整个值域,但是其实我们只需要计算两个信号相乘结果不是 0 的那些项就足够了。
  • 两个信号的遍历顺序是相反的。
  • 根据定义,对于两个信号来说,卷积结果信号 y 的长度是 m + n - 1,其中 m 和 n 分别是 x 和 h 信号的长度。

举例

we have two signals:

x = { 3, 4, 5 } 
h = { 2, 1 }

For instance, the signal x means x[0] = 3, x[1] = 4, x[2] = 5, and x[.] = 0 for all other places.

Their convolution is:

y[0] = x[0] · h[0 - 0] + x[1] · h[0 - 1] + x[2] · h[0-2] = x[0] · h[0 - 0] = 6
y[1] = x[0] · h[1 - 0] + x[1] · h[1 - 1] + x[2] · h[1-2] = x[1] · h[0] + x[0] · h[1]  = 11
y[2] = x[0] · h[2 - 0] + x[1] · h[2 - 1] + x[2] · h[2-2] = x[2] · h[0] + x[1] · h[1]  + x[0] · h[2] = 14
y[3] = 5

Then, y[n] = {6, 11, 14, 5}

卷积时两个数组对齐的样子大概就是这样:

image1

可以看到,核函数 h 是翻转的,然后将翻转后的 h 从结尾开始,对准每一个 x 的值,计算所有对应项的乘积的和。

代码实现

实现时,通常不是去翻转 h,而是翻转 x,这里用了卷积的交换律性质:x * h = h * x。好处是,通常 h 的长度很小而 x 很长,因此 x[i] 通常在外循环,h[j] 是内循环。因此,对于任意一个 h[j], 更容易找到与它对应的 x 的元素,其实就是 x[i-j]。即,最开始的是始终是 x[i] 和 h[0] 相对应,接着 j 增大,x 的索引减小。

for ( i = 0; i < sampleCount; i++ )
{
    y[i] = 0;                       // set to zero before sum
    for ( j = 0; j < kernelCount; j++ )
    {
        if (i >= j)
            y[i] += x[i - j] * h[j];    // convolve: multiply and accumulate
    }
}

二维卷积(Convolution in 2D)

定义

二维卷积的定义类似一维卷积的扩展:

这里 I 是待卷积的对象(例如是图片),H 是核函数(Kernel)。二位情况下,H 的上下左右方向都要翻转。当然,交换律也是成立的:I * H = H * I

这里有非常详细的图文举例

计算某个位置 I(i, j) 的二维卷积时,其实就是把 H 左右上下翻转后,再将 H 的中心点和 I(i, j) 对齐,再计算对应项的乘积的和。例如,我们想要计算下图中的 (1,1) 点的卷积:

image1

结果就是 y(1,1) = 1 * 1 + 2 * 2 + 3 * 1 + 4 * 0 + 5 * 0 + 6 * 0 + 7 * (-1) + 8 * (-2) + 9 * (-1) = -24

不过有一点要注意,如果把上面这个例子直接套用定义的话,那么核函数的坐标范围是 (-1,1) 而不是 (0,2)。这其实牵涉到我们通常使用卷积的方法:通常情况下,我们使用的核函数都是对称的,即它的长和宽都是奇数,正中心被看成是 (0,0) 点(因此代码实现时一定注意)。尺寸是偶数的核函数非常罕见,此时的中心点是中间两个点中任意一个都行,并没有严格定义。

Brute Force 代码实现

到了臭名昭度的二维卷积实现了。它很讨厌的地方有两个:

(1) 因为核函数要翻转,因此在寻找核函数 H 的每个元素对应的图片 I 中的元素时,很容易出错。实现时的细节是这样:

  • 首先,框架依然是两重循环,图片元素 I(i,j) 在外循环,核函数 H(x,y) 是内循环,内循环中 x = 0...m-1, y = 0...n-1,其中m = H.rows, n = H.cols,即核函数的宽和高。很多时候核函数都是一个正方形,即 m = n
  • 考虑到核函数的翻转,对于内循环 (x,y) 位置,我们要考虑的核函数元素其实是 H(m-x-1, n-y-1)
  • 然后,我们要找到 (x,y) 这个位置对应的 I 的元素,它其实是 I(i+x-hx, j+y-hy), 其中 (hx, hy) 是核函数的一半尺寸,即 hx = m/2, hy = n/2。不过注意,如果是偶数尺寸的核函数的话,这里的 hx 和 hy 的定义显然就认为核函数中间偏右边的元素是中心点(就像前面所述,偶数尺寸核函数的中心点并没有严格定义)。
  • 累计每一对对应元素的乘积和到最终的卷积结果:
    C(x,y) += H(m - x - 1, n - y - 1) * I(i + x - hx, j + y - hy);
    

一个 5x5 的核函数的卷积例子:

image2

(2) 输入图片 I 的边界部分的卷积容易出错。当计算边界处的卷积时,核函数 H 有一部分在边界之外。此时,最直观的方法自然是,给输入图片 I 增加一部分 0 padding(即边界周围增加 0),增加的尺寸其实就是 m/2 和 n/2,这样的好处是,I 中每个元素都一样了,无需特殊判断。不过,如果使用 C 或 C++ 实现,这样显然很麻烦,因为这和可能要额外分配一块内存用于存储这个新的图片。因此,只能按照普通方法,对于每个点 I(i,j),判断它是否为边界点,然后特别处理。不过这样似乎会有很多种 if 情况要考虑。其实,根据上面给出的步骤,如果确定了 I(i+x-hx, j+y-hy),其实就很容易了:只需要判断这两个坐标是否越界即可。如果超出图片界限,说明就不考虑它就行了。即,只需要增加一句 if 语句判断一下就行了。

整个部分的代码如下:

int m = H.rows;
int n = H.cols;
int hx = m / 2;
int hy = n / 2;
for (int i = 0; i < I.rows; ++i)  // rows
{
    for (int j = 0; j < I.cols; ++j)  // columns
    {
        for (int x = 0; x < m; ++x)  // kernel rows
        {
            for (int y = 0; y < n; ++y)  // kernel columns
            {
                int ii = i + x - hx;
                int jj = j + y - hy;
                // ignore input samples which are out of bound
                if (ii >= 0 && ii < I.rows && jj >= 0 && jj <= I.cols)
                    C(i, j) += H(m - x - 1, n - y - 1) * I(ii, jj);
            }
        }
    }
}

复杂度上,每个点的卷积计算都需要 mn 次乘法和 mn 次加法。

加速方法:可拆分卷积(Separable Convolution 2D)

参考:

简介

计算一个 m x n 的核函数的卷积时,每个点的卷积计算都需要 mn 次乘法和 mn 次加法。例如对于 3x3 的 kernel,分别需要 9 次乘法和加法。这是复杂度很高的。但是,kernel 有一种性质可以降低时间复杂度。

如果一个 kernel 是可拆分的(separable),那么计算卷积可以降低为每个元素只需要 m+n 次乘法。例如,对于 3x3 的 kernel,就变成了 6 次乘法。

它依照的性质是这样的。一个如下所示的核函数矩阵可以被拆分成列向量和行向量相乘:

于是,和该核函数的卷积就可以被拆分成:

即,两个 2D 矩阵之间的卷积就变成了一个 2D 矩阵和两个 1D 向量的卷积。而矩阵和向量的卷积的话,每个元素只需要 t 次乘法,t 就是向量长度。证明并不难,这个链接给出了详细的证明和例子

常见可拆分的核函数

Gaussian kernel

最常见的高斯滤波。对应 OpenCV 中的 cv::GaussianBlur() 函数。这里是一个常见的 3x3 的全是整数的高斯核函数:

以及 5x5 的:

Smoothing Kernel (Average kernel)

参考:wiki主页介绍

通常用于计算平均值。即 OpenCV 中的 cv::blur() 函数和 box filter 函数 cv::boxFilter().

以及另一个加强中心点的 smoothing kernel/filter:

Sobel Kernel

通常用于 edge detector,对应 OpenCV 中的 cv::Sobel() 函数。

复杂度分析

一个矩阵 I 同一个长度为 m 的列向量 V 做卷积的方法很简单,就是计算每个数 I(i,j) 和 V 的卷积,方法和普通的 1D 卷积相同:首先翻转列向量 V,然后把 V 的中心点对准 I(i,j),再计算对应项乘积和。同理,和行向量的卷积类似。显然,和一个列向量卷积只需要 m 次乘法和 m-1 次加法,和行向量卷积需要 n 次乘法和 n-1 次加法,因此每个数的卷积一共需要 m+n 次乘法和 m+n-2 次加法。这通常比普通的卷积所需的 mn 次乘法和 mn - 1 次加法要好得多。

不过,拆分核函数的一个缺点是,需要用额外的和输入矩阵 I 相同大小的空间来记录第一次和向量卷积的结果(如果不能修改输入的 I 本身的话)。当 I 很大时,这个缺点很明显。不过,考虑到通常的 I 只是类似一张图片之类,这个额外空间可以接受。

代码实现

输入图片是一维数组

如果输入图片 I 是二维矩阵,几乎没什么好说的,直接按照流程计算和列向量卷积,再计算和行向量卷积即可。唯一要注意的是边界情况。可以采用像前面 1D 情况的方法,对于每个 kernel 元素对应的图片元素,先判断它是否越界,不越界才会使用它。

但是,在计算机视觉应用中,通常一个图片矩阵的输入是一个 row-major 1D 行向量数组,而不是 2D 矩阵,例如 OpenCV 中的 Mat。这就给实现带来了一些难度。

首先,两个向量的卷积顺序是可以互换的,例如先计算行向量卷积,再计算列向量的卷积。而由于输入时一维数组,那么先计算行向量卷积,再计算列向量卷积会更方便。因为一维数组和行向量的卷积就比较直观了,唯一要注意的就是边界情况了。

那么,如何计算一维数组和列向量之间的卷积?对于 I(i,j),直观上,计算它的卷积其实就是找到它上下几个元素和对应的 kernel 元素的乘积和。因此,似乎不可避免的要使用乘法计算二维坐标了,例如类似 i*cols+j 这样。而乘法肯定会比较慢的。其实,有一种巧妙的实现方法可以避开计算二维坐标并加速:

  • 首先,依然是三层循环。最外层循环肯定是第 i 行。如果是直观法,肯定是第 j 列在第二层循环,kernel[k] 在最里层。而巧妙方法是,我们可以反过来,将 kernel[k] 放在第二层,把第 j 列放在最里层(最外层依然是 i)。这样的意思就变成了,我们在平移 kernel 列向量并累加计算卷积 C(i, j) += (*ptr) * kernel[k],这里的 ptr 就是图像元素的指针,它始终会 ptr++ 递增。
  • 接着,第三层循环走完,j 到头后,第二层循环设置 k–(因为 kernel 要翻转后才使用),而 j 此时又会变成从 0 开始了。因为 ptr 始终在递增,它自动会指向图片的下一行的第一个元素。然后,由于更新了 k,因此 kernel[k] 也就是 ptr 这一行对应的 kernel 元素。因此,接着重复上面的累加计算卷积的步骤,新的这一行的第 j 列的累计和结果还是被累加到了 C(i, j) 中。这正是我们想要的。
  • 当 k = 0 后,第二层循环停止,C(i, j) 中就是卷积结果了。
  • 接着就回到了最外层循环,i++,此时只需要更新 ptr 为上一个 i 开始时的 ptr+cols 就行了。因此,额外使用一个变量 ptr_init,在第二层循环开始前记录一下 ptr 的初始位置即可。

上述加速方法巧妙使用了 C 指针,避免使用乘法,而指针操作速度是飞快的。

边界情况的处理

另一个繁琐的地方自然又是臭名昭度的边界情况了。当然我们可以采用前面给出的普通卷积中使用的方法,增加对每个位置的判断条件,这样最简单直接。不过,显然 if 语句也是耗时的,并且繁多的判断语句会让代码更加繁琐。为了加速,实现上其实可以不用判断语句,而是将整个 I 分成 3 个区间来分别处理。例如,在和列向量做卷积时,可以将 I 的高度方向(x 方向)分成 3 个区间:

  • 上边界部分 [0, hx-1];
  • 中间部分 [hx, rows-hx-1];
  • 下边界部分[rows-hx, rows-1];

区间全是闭区间。这里的 hx 就是高的一半,依然是之前定义的 hx = m/2。实现时,对上边界部分和下边界部分特殊处理,中间部分就是普通情况。同理,在和行向量卷积时,对 I 的宽度方向(y 方向)的区间划分完全类似。举个例子就很明显了,对于 5x5 的图片 I 以及 3x3 的核函数,三个高度方向的区间就是 [0,0],[1,3] 和 [4,4]。对 [0,0] 和 [4,4] 做特殊处理就行了。这样就完全不用任何 if 语句即可,可以很大程度提升效率。

上边界区间和下边界区间的特殊处理也不难,无非是 kernel 中能够使用的元素少了几个。上边界区间的话,kernel 范围其实就是[hx+offset, 0](区间是倒过来的,考虑到翻转),这里 offset 是一个新增的偏差值,初始是 0,每一行后就 offset++。下边界正好相反,kernel 范围是 [n-1, offset],其中 offset 初值是 1,依然是每一行后 offset++,n 是 kernel 长度。

最终代码

下面给出了笔者参考了教程后自己实现的针对作为一维数组输入的图片的可拆分卷积函数,来自这个链接。里面使用了所有尽可能的加速方法。因为要区分多个区间,代码显得很长,但是运行速度飞快。

//! 使用可拆分卷积实现图像卷积
/*!
    @param in 输入图像的指针,这里输入必须是一个 1-channel 的图片,例如灰度图
    @param rows 图像行数
    @param cols 图像列数
    @param kernelX 核函数拆分后的列向量
    @param kernelY 核函数拆分后的行向量
    @param sizeX 核函数列向量长度
    @param sizeY 核函数行向量长度
    @param out 输出图像的指针
*/
bool convolve2DSeparable(
    unsigned char* in, int rows, int cols, float* kernelX, int sizeX, float* kernelY, int sizeY, unsigned char* out)
{
    if (!in || !out || !rows || !cols || !kernelX || !kernelY)
        return false;
    int N = rows * cols;
    float* tmp = new float[N];  // 额外分配空间用于存储第一次卷积的结果
    for (int i = 0; i < N; ++i)
        tmp[i] = 0;

    // 行向量卷积
    int hy = (sizeY >> 1);
    int jump_gap = ((sizeY - 1) >> 1);  // 每一行遍历后,指向图片元素的指针需要调过的元素个数
    int left_border = hy;
    int right_border = cols - hy;
    float* ptr_tmp = tmp;  // 用指针运算可以加快速度
    unsigned char* ptr_in = in;
    int offset = 0;  // 用于每次向右移动 kernel 一位
    for (int i = 0; i < rows; ++i)
    {
        // 左边界区间 [0, hy - 1],kernel 右边有一部分用不到
        offset = 0;  // 用于每次向右移动 kernel 一位
        for (int j = 0; j < left_border; ++j)
        {
            *ptr_tmp = 0;  // 其实最开始已经初始化全部的 tmp 为 0 了,为了保险还是再初始化一下
            for (int k = hy + offset, t = 0; k >= 0; t++, k--)
                *ptr_tmp += *(ptr_in + t) * kernelY[k];  // 这里第一个 * 是取地址操作符
            ptr_tmp++;
            offset++;
        }
        // 中区间 [hy, cols - hy - 1],一般情况
        for (int j = left_border; j < right_border; ++j)
        {
            *ptr_tmp = 0;
            for (int k = sizeY - 1, t = 0; k >= 0; t++, k--)
                *ptr_tmp += *(ptr_in + t) * kernelY[k];
            ptr_tmp++;
            ptr_in++;  // 注意这里也递增了图片指针,始终让它指向 Kernel 的左边界对应的图片的点
        }
        // 右边界区间 [colos - hy, cols - 1],kernel 左边有一部分用不到
        offset = 1;  // 这次它表示在右边界处理中,kernel 最右边的元素位置
        for (int j = right_border; j < cols; ++j)
        {
            *ptr_tmp = 0;
            for (int k = sizeY - 1, t = 0; k >= offset; t++, k--)
                *ptr_tmp += *(ptr_in + t) * kernelY[k];
            ptr_tmp++;
            ptr_in++;
            offset++;
        }
        // 一行的行向量卷积结束后,输入指针会停在离右边界的距离是 sizeY/2 再加上 1 的位置。
        // 因此需要跳过一部分数据到达下一行。
        ptr_in += jump_gap;
    }

    // 列向量卷积。实现上很巧妙。
    int hx = (sizeX >> 1);  // 类似上面的 hy
    int up_border = hx;     // 区间的两个边界
    int bottom_border = rows - hx;
    ptr_tmp = tmp;  // 从头开始
    unsigned char* ptr_out = out;
    float* sum = new float[cols];  // 这一步额外需要一个行向量记录每一行的结果,因为输出的 out 数组是 unsigned char,
                                   // 但是卷积结果是 float,因此最好还是额外定义一个数组先存储结果,再赋值到 out 中
    for (int i = 0; i < cols; ++i)
        sum[i] = 0;
    // 上边界区间 [0, hx - 1]
    offset = 0;
    for (int i = 0; i < up_border; ++i)
    {
        for (int k = hx + offset; k >= 0; k--)
        {
            // 当 k 减小后,j 又回到了从 0 开始增大,因此,同一列的乘积和结果累积到了 sum[j] 中
            for (int j = 0; j < cols; ++j)
            {
                sum[j] += *ptr_tmp * kernelX[k];
                ptr_tmp++;  // 注意它是始终递增的
            }
        }
        for (int j = 0; j < cols; ++j)
        {
            *ptr_out = (unsigned char)(sum[j]);  // 这里其实默认将 float 转成了 unsigned char
            ptr_out++;
            sum[j] = 0;  // 记得清零供下次使用
        }
        ptr_tmp = tmp;  // 上边界区间中,计算完一行卷积后,ptr_tmp 始终指向开头
        offset++;
    }
    // 中间区间 [hx, rows - hx - 1]。做法和上边界区间很类似。
    float* ptr_init = ptr_tmp; // 拷贝一份初始指针
    for (int i = up_border; i < bottom_border; ++i)
    {
        for (int k = sizeX - 1; k >= 0; k--)
        {
            for (int j = 0; j < cols; ++j)
            {
                sum[j] += *ptr_tmp * kernelX[k];
                ptr_tmp++;
            }
        }
        for (int j = 0; j < cols; ++j)
        {
            *ptr_out = sum[j];
            ptr_out++;
            sum[j] = 0; 
        }
        ptr_init += cols;  // 左上角指针要向下移动一行
        ptr_tmp = ptr_init;
    }
    // 下边界区间 [rows - hx, rows - 1]
    offset = 1;
    for (int i = bottom_border; i < rows; ++i)
    {
        for (int k = sizeX - 1; k >= offset; k--)
        {
            for (int j = 0; j < cols; ++j)
            {
                sum[j] += *ptr_tmp * kernelX[k];
                ptr_tmp++;
            }
        }
        for (int j = 0; j < cols; ++j)
        {
            *ptr_out = sum[j];
            ptr_out++;
            sum[j] = 0;
        }
        offset++;
        ptr_init += cols;
        ptr_tmp = ptr_init;
    }
    // 最后记得释放内存
    delete[] sum;
    delete[] tmp;
    return true;
}

Search

    Table of Contents