设: \(f(x), g(x)\) 是 \(\mathbb{R}1\) 上的两个可积函数， \(f(x)\) 和 \(g(x)\) 的卷积定义如下:

\( (f * g)(t) \, \stackrel{\mathrm{def}}{=}\ \int_{-\infty}^\infty f(\tau)\, g(t-\tau) \, d\tau = \int_{-\infty}^\infty f(t-\tau) \, g(\tau)\, d\tau \)

按公式来，卷积计算的等效C++代码如下

/*
 * m - a 的长度
 * n - b 的长度
 */
void convolution(double *a, double *b, double *output, int m, int n)
{ 
    double *x = new double[m+n-1]; 
    // do convolution 
    for (int i = 0; i < m+n-1; i++)
    { 
        x[i] = 0.0;
        for (int j = 0; j < m; j++) 
        { 
            if (i-j > 0 && i-j < n)
            x[i] += a[j] * b[i-j]; 
        } 
     }
     // set value to the output array 
     for (int i = 0; i < m; i++) 
     output[i] = x[i + (n-1) / 2];
     delete[] x;
}

R 中有现成的关于向量的卷积计算函数 convolve()。使用 R 中的 convolve() 时，默认 使用的是 type=”circular”，这时要求输入参数 a、b 长度相同。长度不同，则只能使用 type=”open” 或者 “filter”。需要注意的是，在 R 中使用 convolve() 只能做一维向量 的计算，并且计算的结果与其它的库得到的有些差异，关键差异还是在边界处。 下面是在 R 中使用 convolve() 计算向量的卷积的例子:

> a <- 1:9
> convolve(a, c(1,2,1,rep(0,6))) # circular type convolution
[1]  8 12 16 20 24 28 32 27 13
> convolve(a, c(1,2,1), type="open") # open type convolution
[1]  1  4  8 12 16 20 24 28 32 26  9
> convolve(a,c(1,2,1), type="filter") # filter type convolution
[1]  8 12 16 20 24 28 32

为了对比，我们也来看一下例如 python 中的 mahotas 的 convolve 的计算结果。 默认情况下， mahotas 中的 convolve 使用的 mode 为 reflect。mahotas 中的 mode 是处理边界处的方式，有: reflect, nearest, wrap, mirror, constant, ignore。

>>> import numpy as np
>>> import mahotas as mh 
>>> a = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9])
>>> b = np.array([1, 2, 1])
>>> mh.convolve(a, b)
array([ 5,  8, 12, 16, 20, 24, 28, 32, 35])
>>> mh.convolve(a, b, mode='nearest')
array([ 5,  8, 12, 16, 20, 24, 28, 32, 35])
>>> mh.convolve(a, b, mode='wrap')
array([13,  8, 12, 16, 20, 24, 28, 32, 27])
>>> mh.convolve(a, b, mode='mirror')
array([ 6,  8, 12, 16, 20, 24, 28, 32, 34])
>>> mh.convolve(a, b, mode='constant')
array([ 4,  8, 12, 16, 20, 24, 28, 32, 26])
>>> mh.convolve(a, b, mode='ignore')
array([ 4,  8, 12, 16, 20, 24, 28, 32, 26])

一般地，为了加快计算，一般不大会使用上面提到的直接展开卷积定义式进行循环计算的 方法，在工程上我们会使用傅立叶的卷积定理来计算卷积:

\( (f*g)(t) = \mathcal{F}^{-1}( \mathcal{F}(f) \times \mathcal{F}(g) )\)

其中 \(\mathcal{F}(f)\) 为 \(f(t)\) 的傅立叶变换。

R 中的 convolve 的实现使用的正是傅立叶变换，并用卷积定理来求解的，其代码如下:

convolve <- function (x, y, conj = TRUE, type = c("circular", "open", "filter"))
{
    type <- match.arg(type)
    n <- length(x)
    ny <- length(y)
    Real <- is.numeric(x) && is.numeric(y)
    if (type == "circular") {
        if (ny != n)
            stop("length mismatch in convolution")
    }
    else {
        n1 <- ny - 1
        x <- c(rep.int(0, n1), x)
        n <- length(y <- c(y, rep.int(0, n - 1)))
    }
    x <- fft(fft(x) * (if (conj)
        Conj(fft(y))
    else fft(y)), inverse = TRUE)
    if (type == "filter")
        (if (Real)
            Re(x)
        else x)[-c(1L:n1, (n - n1 + 1L):n)]/n
    else (if (Real)
        Re(x)
    else x)/n
}

为了在 R 中提供可用于图像的2维卷积计算，可以对这个实现进行改进。这段代码已经够 短了，而如果我们只考虑 open 边界条件，还可以进一步简化为:

convolve.open <- function (x, y)
{
    m <- length(x)
    n <- length(y)

    x <- c(rep.int(0, n - 1), x)  # x前面补 n-1 个 0
    y <- c(y, rep.int(0, m - 1))  # y后面补 m-1 个 0

    c <- fft(fft(x) * Conj(fft(y)), inverse = TRUE)
    Re(c)/length(y)
}

这时，结构就非常清晰了，可以进一步把它改造成2维的卷积，我的代码如下:

conv2d <- function (x, y)
{
    if (is.vector(x) | is.vector(y)) 
        convolve(x, y)
    else {
        m <- dim(x)
        n <- dim(y)

        # x前面补 n-1 个 0
        a <- matrix(0, m[1]+n[1]-1, m[2]+n[2]-1)
        a[n[1]:(m[1]+n[1]-1), n[2]:(m[2]+n[2]-1)] <- x

        # y后面补 m-1 个 0
        b <- matrix(0, m[1]+n[1]-1, m[2]+n[2]-1)
        b[1:n[1], 1:n[2]] <- y

        c <- fft(fft(a) * Conj(fft(b)), inverse = TRUE)
        c <- Re(c)/((m[1]+n[1]-1)*(m[2]+n[2]-1))

        # 裁成 m 的尺寸并输出
        c[ceiling(n[1]/2):(ceiling(n[1]/2)-1+m[1]),
          ceiling(n[2]/2):(ceiling(n[2]/2)-1+m[2])]
    }
}

下面，可以来试一下这个卷积的图像处理效果了。先做两个差分模板:

ptn.a <- matrix(0, 4, 4)
ptn.a[,1:2] <- -1
ptn.a[,3:4] <- 1
ptn.a
#  [,1] [,2] [,3] [,4]
#  [1,]   -1   -1    1    1
#  [2,]   -1   -1    1    1
#  [3,]   -1   -1    1    1
#  [4,]   -1   -1    1    1
ptn.b <- t(ptn.a)
ptn.b
#  [,1] [,2] [,3] [,4]
#  [1,]   -1   -1   -1   -1
#  [2,]   -1   -1   -1   -1
#  [3,]    1    1    1    1
#  [4,]    1    1    1    1

再打开昨天的那张图像来试试，参见R操作图像数据，并且会用到这里的一些 函数，不再在这里重复。

KonikHorses 原图

使用 ptn.a 做一次卷积，并显示一下效果:

pic <- read.yuv('horses-640x360.jpg')  # 读取图像
a <- conv2d(pic$y, ptn.a)              # 使用模板做卷积
b <- conv2d(pic$y, ptn.b)
# 有必要做一下正规化操作
y.r <- range(pic$y)
a.r <- range(a)
b.r <- range(b)
a.y <- (a - a.r[1]) * (y.r[2] - y.r[1]) / (a.r[2] - a.r[1])
b.y <- (b - b.r[1]) * (y.r[2] - y.r[1]) / (b.r[2] - b.r[1])
a.yuv <- list(y=a.y, u=pic$u, v=pic$v) # 做成一幅 yuv 图像
b.yuv <- list(y=b.y, u=pic$u, v=pic$v)
show.yuv(a.yuv)                        # 显示效果
show.yuv(b.yuv)

r-conv-h 与ptn.a卷积的效果

r-conv-v 与ptn.b卷积的效果