Lode 的计算机图形学教程

图像滤镜

目录

返回目录

简介

图像滤镜可以让你对照片应用各种特效。本文所介绍的图像滤镜使用一个 2D 滤波矩阵，与 Paint Shop Pro 中的"用户自定义滤镜"以及 Photoshop 中的"自定义滤镜"类似。

卷积

图像滤镜的核心思路是：给定一个 2D 滤波矩阵和一幅 2D 图像，对图像中的每个像素计算乘积之和。每个乘积由当前像素或其邻域像素的颜色值，与滤波矩阵中对应位置的值相乘得到。滤波矩阵的中心元素与当前像素相乘，其余元素与对应的邻域像素相乘。

这种对两个 2D 函数的对应元素求乘积之和，并让其中一个函数在另一个函数的每个位置上滑动的运算，称为卷积（Convolution）或相关（Correlation）。二者的区别在于：卷积需要对滤波矩阵进行镜像翻转，但通常滤波矩阵是对称的，因此二者结果相同。

基于卷积的滤镜相对简单。更复杂的滤镜可以使用更灵活的函数，实现更复杂的效果（例如 Photoshop 中的彩色铅笔滤镜），但本文不作讨论。

2D 卷积运算需要一个四重循环，因此速度并不特别快，除非使用较小的滤波核。本文通常使用 3x3 或 5x5 的滤波核。

关于滤波矩阵，有以下几条规则：
图像的尺寸是有限的，例如在计算左侧边缘的像素时，其左边不存在更多像素，而卷积运算又需要这些像素。此时可以用 0 代替，也可以绕回到图像的另一侧。本教程选择绕回方式，因为这可以用取模运算轻松实现。

应用滤镜后，像素值可能为负数或大于 255。遇到这种情况，可以将其截断：小于 0 的值置为 0，大于 255 的值置为 255。对于负值，也可以取其绝对值。

在傅里叶域（频域）中，卷积运算变为乘法运算，速度更快。在傅里叶域中，可以更快速地应用功能更强大、尺寸更大的滤波器，尤其是配合快速傅里叶变换（FFT）使用时。详细内容请参阅傅里叶变换相关文章。本文将介绍几种典型的小型滤波器，如模糊、边缘检测和浮雕。

图像滤镜目前还不适合实时应用和游戏，但在图像处理领域非常有用。

数字音频和电子滤波器同样使用卷积，只不过是在 1D 场景下。

下面是用于测试各种滤镜的代码。除滤波矩阵外，代码还引入了一个乘法因子（factor）和一个偏置值（bias）。应用滤镜后，结果会先乘以 factor，再加上 bias。因此，如果滤波矩阵中某个元素为 0.25，而 factor 设为 2，则该元素实际上等效于 0.5。bias 可用于提高输出图像的亮度。

单个像素的计算结果先以浮点数 red、green、blue 存储，再转换为整数写入结果缓冲区。

滤波计算本身是一个四重循环，需要遍历图像的每个像素，再遍历滤波矩阵的每个元素。imageX 和 imageY 的计算方式为：对于滤波矩阵的中心元素，坐标为 (x, y)；对于其他元素，坐标为图像中 (x, y) 左、右、上、下方向的某个像素。其坐标通过对图像宽度（w）或高度（h）取模来实现绕回。取模之前还需加上 w 或 h，以确保负数情况下取模结果正确。这样，像素 (-1, -1) 就能正确映射为像素 (w-1, h-1)。

#define filterWidth 3
#define filterHeight 3

double filter[filterHeight][filterWidth] =
{
   0, 0, 0,
   0, 1, 0,
   0, 0, 0
};

double factor = 1.0;
double bias = 0.0;

int main(int argc, char *argv[])
{
  //load the image into the buffer
  unsigned long w = 0, h = 0;
  std::vector<ColorRGB> image;
  loadImage(image, w, h, "pics/photo3.png");
  std::vector<ColorRGB> result(image.size());

  //set up the screen
  screen(w, h, 0, "Filters");

  ColorRGB color; //the color for the pixels

  //apply the filter
  for(int x = 0; x < w; x++)
  for(int y = 0; y < h; y++)
  {
    double red = 0.0, green = 0.0, blue = 0.0;

    //multiply every value of the filter with corresponding image pixel
    for(int filterY = 0; filterY < filterHeight; filterY++)
    for(int filterX = 0; filterX < filterWidth; filterX++)
    {
      int imageX = (x - filterWidth / 2 + filterX + w) % w;
      int imageY = (y - filterHeight / 2 + filterY + h) % h;
      red += image[imageY * w + imageX].r * filter[filterY][filterX];
      green += image[imageY * w + imageX].g * filter[filterY][filterX];
      blue += image[imageY * w + imageX].b * filter[filterY][filterX];
    }

    //truncate values smaller than zero and larger than 255
    result[y * w + x].r = min(max(int(factor * red + bias), 0), 255);
    result[y * w + x].g = min(max(int(factor * green + bias), 0), 255);
    result[y * w + x].b = min(max(int(factor * blue + bias), 0), 255);
  }

  //draw the result buffer to the screen
  for(int y = 0; y < h; y++)
  for(int x = 0; x < w; x++)
  {
    pset(x, y, result[y * w + x]);
  }

  //redraw & sleep
  redraw();
  sleep();
}

如果希望对小于零的值取绝对值而非截断，请改用以下代码：

    //take absolute value and truncate to 255
    result[y * w + x].r = min(abs(int(factor * red + bias)), 255);
    result[y * w + x].g = min(abs(int(factor * green + bias)), 255);
    result[y * w + x].b = min(abs(int(factor * blue + bias)), 255);

当前填入的滤波矩阵为：

[ 0 0 0 ]
[ 0 1 0 ]
[ 0 0 0 ],

该矩阵不做任何处理，只返回原始图像，因为只有中心值为 1，每个像素都乘以 1。

代码会尝试加载图像 "pics/photo3.bmp"。该图像可在此处下载。

原始图像如下所示：



接下来，我们将通过修改滤波矩阵的定义并运行代码，对图像应用多种滤镜。

模糊

模糊效果可以通过取当前像素与其 4 个邻域像素的平均值来实现。将当前像素与其 4 个邻域像素求和后除以 5，即在滤波矩阵中填入 5 个值为 0.2 的元素：

#define filterWidth 3
#define filterHeight 3

double filter[filterHeight][filterWidth] =
{
   0.0, 0.2,  0.0,
   0.2, 0.2,  0.2,
   0.0, 0.2,  0.0
};

double factor = 1.0;
double bias = 0.0;

使用如此小的滤波矩阵，只能产生非常轻微的模糊效果：



使用更大的滤波核可以获得更强的模糊效果（别忘了同步修改 filterWidth 和 filterHeight 的值）：

#define filterWidth 5
#define filterHeight 5

double filter[filterHeight][filterWidth] =
{
  0, 0, 1, 0, 0,
  0, 1, 1, 1, 0,
  1, 1, 1, 1, 1,
  0, 1, 1, 1, 0,
  0, 0, 1, 0, 0,
};

double factor = 1.0 / 13.0;
double bias = 0.0;

滤波矩阵所有元素之和应为 1，但此处并未在矩阵内填入浮点数，而是将 factor 除以所有元素之和（即 13）来等效实现。

这样可以获得更明显的模糊效果：



模糊程度越强，所需的滤波核就越大，或者也可以多次应用同一个小型模糊滤镜。

如果滤波核是一个全部填充相同值的矩形（并配以适当的缩放因子使所有元素之和为 1.0），则该模糊称为均值模糊（box blur）。若需要非常大的均值模糊，本教程中的朴素卷积代码会过慢。但可以用更快的算法来实现：由于每个值的权重相同，可以逐行遍历图像像素，对 N 个值（N 为矩形框的宽度）求和并除以适当的缩放因子。对于每个后续像素，加入矩形框中新出现的像素值，减去从矩形框左侧移出的像素值。对每条扫描线水平处理完毕后，再垂直方向执行相同操作（为优化 CPU 缓存利用率，在垂直方向处理时，应确保实际上仍以扫描线顺序操作，而非按列操作，因此需要为每列维护一个求和值）。这一切都需要仔细处理边界情况（矩形框部分超出图像范围时）以及图像尺寸小于矩形框的情况。本文不提供相关代码，因为这已超出本教程的范围。

高斯模糊

上述模糊所用的滤波核较为生硬。使用高斯核可以获得更平滑的模糊效果。在高斯核中，数值随距中心距离的增大而指数级衰减。公式为：G(x) = exp(-x * x / 2 * sigma * sigma) / sqrt(2 * pi * sigma * sigma)

对于 2D 情形，先在 X 方向应用该公式，再在 Y 方向应用（二者可分离），合并后为： G(x, y) = exp(-(x * x + y * y) / (2 * sigma * sigma)) / (2 * pi * sigma * sigma)

公式中各参数说明：
*) sigma 决定模糊半径（理论上半径无限大，但由于指数衰减，实际上存在一个值小到肉眼不可见的截止点，sigma 越大，截止点越远） *) x 和 y 为坐标值，且以滤波核中心为原点

上述公式可用于构建任意大小的滤波核。以下是一些可直接使用的简单示例：

3x3 滤波核近似：

#define filterWidth 3 #define filterHeight 3 double filter[filterHeight][filterWidth] = { 1, 2, 1, 2, 4, 2, 1, 2, 1, }; double factor = 1.0 / 16.0; double bias = 0.0;

5x5 滤波核近似：

#define filterWidth 5 #define filterHeight 5 double filter[filterHeight][filterWidth] = { 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, }; double factor = 1.0 / 256.0; double bias = 0.0;

精确（非近似）示例：

#define filterWidth 3 #define filterHeight 3 double filter[filterHeight][filterWidth] = { 0.077847, 0.123317, 0.077847, 0.123317, 0.195346, 0.123317, 0.077847, 0.123317, 0.077847, }; double factor = 1.0; double bias = 0.0;

对于较大的模糊半径（例如绘图软件中的高斯模糊），需要更大的滤波核。本教程中的朴素卷积实现对于大半径高斯模糊在实践中会过慢。但有解决方案：使用本系列傅里叶变换教程中介绍的傅里叶变换方法，或更快的近似方法：连续多次执行均值模糊，三次均值模糊已能很好地近似高斯模糊。如何实现快速均值模糊已在上一章节中介绍。该方法有效的原因在于，高斯分布自然地从多个过程的叠加中涌现。

运动模糊

运动模糊通过仅在一个方向上进行模糊来实现。以下是一个 9x9 的运动模糊滤波核：

#define filterWidth 9
#define filterHeight 9

double filter[filterHeight][filterWidth] =
{
  1, 0, 0, 0, 0, 0, 0, 0, 0,
  0, 1, 0, 0, 0, 0, 0, 0, 0,
  0, 0, 1, 0, 0, 0, 0, 0, 0,
  0, 0, 0, 1, 0, 0, 0, 0, 0,
  0, 0, 0, 0, 1, 0, 0, 0, 0,
  0, 0, 0, 0, 0, 1, 0, 0, 0,
  0, 0, 0, 0, 0, 0, 1, 0, 0,
  0, 0, 0, 0, 0, 0, 0, 1, 0,
  0, 0, 0, 0, 0, 0, 0, 0, 1,
};

double factor = 1.0 / 9.0;
double bias = 0.0;



效果犹如相机从左上方向右下方移动，因此得名运动模糊。

边缘检测

用于检测水平边缘的滤波核示例如下：

#define filterWidth 5
#define filterHeight 5

double filter[filterHeight][filterWidth] =
{
   0,  0, -1,  0,  0,
   0,  0, -1,  0,  0,
   0,  0,  2,  0,  0,
   0,  0,  0,  0,  0,
   0,  0,  0,  0,  0,
};

double factor = 1.0;
double bias = 0.0;

此处选用 5x5 而非 3x3 的滤波核，是因为 3x3 滤波核在当前图像上的结果过暗。注意，现在所有元素之和为 0，这将产生一幅非常暗的图像，只有检测到的边缘处才有颜色。



该滤波核能够检测水平边缘的原因在于：使用该滤波核的卷积运算可以看作一种离散微分：取当前像素值减去前一个像素值，得到的差值代表两者之间的变化量，即函数的斜率。

以下滤波核用于检测垂直边缘，同时利用当前像素上方和下方的像素值：

#define filterWidth 5
#define filterHeight 5

double filter[filterHeight][filterWidth] =
{
   0,  0, -1,  0,  0,
   0,  0, -1,  0,  0,
   0,  0,  4,  0,  0,
   0,  0, -1,  0,  0,
   0,  0, -1,  0,  0,
};

double factor = 1.0;
double bias = 0.0;



以下是另一种可能的滤波核，擅长检测 45° 方向的边缘。值 '-2' 的选取没有特别原因，只需确保所有值之和为 0 即可。

#define filterWidth 5
#define filterHeight 5

double filter[filterHeight][filterWidth] =
{
  -1,  0,  0,  0,  0,
   0, -2,  0,  0,  0,
   0,  0,  6,  0,  0,
   0,  0,  0, -2,  0,
   0,  0,  0,  0, -1,
};

double factor = 1.0;
double bias = 0.0;



以下是一个能够检测所有方向边缘的简单边缘检测滤波核：

#define filterWidth 3
#define filterHeight 3

double filter[filterHeight][filterWidth] =
{
  -1, -1, -1,
  -1,  8, -1,
  -1, -1, -1
};

double factor = 1.0;
double bias = 0.0;



锐化

图像锐化与边缘检测非常相似：将原始图像与边缘检测后的图像叠加，结果是一幅边缘得到增强、看起来更清晰的新图像。这两幅图像的叠加，可以通过取上一示例中的边缘检测滤波核并将其中心值加 1 来实现。此时滤波矩阵所有元素之和为 1，输出图像的亮度与原图相同，但更加清晰。

#define filterWidth 3
#define filterHeight 3

double filter[filterHeight][filterWidth] =
{
  -1, -1, -1,
  -1,  9, -1,
  -1, -1, -1
};

double factor = 1.0;
double bias = 0.0;



以下是一个效果更为细腻的锐化滤波核：

#define filterWidth 5
#define filterHeight 5

double filter[filterHeight][filterWidth] =
{
  -1, -1, -1, -1, -1,
  -1,  2,  2,  2, -1,
  -1,  2,  8,  2, -1,
  -1,  2,  2,  2, -1,
  -1, -1, -1, -1, -1,
};

double factor = 1.0 / 8.0;
double bias = 0.0;



以下是一个会过度突出边缘的滤波核：

#define filterWidth 3
#define filterHeight 3

double filter[filterHeight][filterWidth] =
{
   1,  1,  1,
   1, -7,  1,
   1,  1,  1
};

double factor = 1.0;
double bias = 0.0;


浮雕

浮雕滤镜可以为图像添加 3D 阴影效果，非常适合用作图像的凹凸贴图（bumpmap）。其实现方式是：取中心一侧的像素值，减去另一侧的像素值。像素的计算结果可能为正值或负值。为了将负值用作阴影、正值用作高光来生成凹凸贴图，需要为图像添加 128 的偏置值。这样，图像的大部分区域将呈现为灰色，而边缘处则呈现为深灰色/黑色或浅灰色/白色。

以下是一个 45° 角方向的浮雕滤波核示例：

#define filterWidth 3
#define filterHeight 3

double filter[filterHeight][filterWidth] =
{
  -1, -1,  0,
  -1,  0,  1,
   0,  1,  1
};

double factor = 1.0;
double bias = 128.0;



如果确实要将其用作凹凸贴图，可以将其转为灰度图：



以下是一个效果更为夸张的浮雕滤波核：

#define filterWidth 5
#define filterHeight 5

double filter[filterHeight][filterWidth] =
{
  -1, -1, -1, -1,  0,
  -1, -1, -1,  0,  1,
  -1, -1,  0,  1,  1,
  -1,  0,  1,  1,  1,
   0,  1,  1,  1,  1
};

double factor = 1.0;
double bias = 128.0;


均值滤波与中值滤波

均值滤波（Mean Filter）和中值滤波（Median Filter）均可用于去除图像噪声。均值滤波取当前像素与其邻域像素的平均值，例如使用 8 个邻域像素时，其滤波核为：

#define filterWidth 3
#define filterHeight 3

double filter[filterHeight][filterWidth] =
{
  1, 1, 1,
  1, 1, 1,
  1, 1, 1
};

double factor = 1.0 / 9.0;
double bias = 0.0;

这是一个普通的模糊滤波核。我们可以用以下含有"椒盐噪声"（Salt and Pepper Noise）的图像来测试：



应用后，结果图像会变得模糊：



中值滤波的处理方式类似，但不是取均值，而是取中位数。中位数的求法是：将所有值从小到大排序，然后取中间的值。若中间有两个值，则取其平均值。中值滤波在去除椒盐噪声方面效果更好，因为它能完全消除噪声。均值滤波在计算平均值时仍会受到噪声像素颜色值的影响，而取中位数时只保留一个或两个正常像素的颜色值。不过，中值滤波同样会降低图像质量。

中值滤波无法通过卷积实现，需要借助排序算法。本例选用了梳排序（combsort），这是一种相对较快的排序算法。

要对当前像素及其 8 个邻域像素求中位数，可将 filterWidth 和 filterHeight 设为 3，也可以设置更大的值以去除更大的噪声颗粒。

数组 red、green 和 blue 将存储当前像素及其所有邻域像素的值，这些数组将由排序算法排序，以便取得中位数。主函数负责应用滤镜、计算中位数并绘制结果。

#define filterWidth 3
#define filterHeight 3

//color arrays
int red[filterWidth * filterHeight];
int green[filterWidth * filterHeight];
int blue[filterWidth * filterHeight];

int selectKth(int* data, int s, int e, int k);

int main(int argc, char *argv[])
{
  //load the image into the buffer
  unsigned long w = 0, h = 0;
  std::vector<ColorRGB> image;
  loadImage(image, w, h, "pics/noise.png");
  std::vector<ColorRGB> result(image.size());

  //set up the screen
  screen(w, h, 0, "Median Filter");

  ColorRGB color; //the color for the pixels

  //apply the filter
  for(int y = 0; y < h; y++)
  for(int x = 0; x < w; x++)
  {
    int n = 0;
    //set the color values in the arrays
    for(int filterY = 0; filterY < filterHeight; filterY++)
    for(int filterX = 0; filterX < filterWidth; filterX++)
    {
      int imageX = (x - filterWidth / 2 + filterX + w) % w;
      int imageY = (y - filterHeight / 2 + filterY + h) % h;
      red[n] = image[imageY * w + imageX].r;
      green[n] = image[imageY * w + imageX].g;
      blue[n] = image[imageY * w + imageX].b;
      n++;
    }

    int filterSize = filterWidth * filterHeight;
    result[y * w + x].r = red[selectKth(red, 0, filterSize, filterSize / 2)];
    result[y * w + x].g = green[selectKth(green, 0, filterSize, filterSize / 2)];
    result[y * w + x].b = blue[selectKth(blue, 0, filterSize, filterSize / 2)];
  }

  //draw the result buffer to the screen
  for(int y = 0; y < h; y++)
  for(int x = 0; x < w; x++)
  {
    pset(x, y, result[y * w + x]);
  }

  //redraw & sleep
  redraw();
  sleep();
}

数组中存储了所操作矩形区域内每个颜色通道的值，但尚未排序，因此无法直接取中位数。排序后取中间元素是一种方法，但从理论上讲，使用选择算法来选取第 k 大的元素（k = size / 2）更快。以下是一个非常简单的选择算法实现。也可以使用标准 C++ 函数 nth_element，它更简单也更快，但本教程中我们自行实现所有算法。需要注意的是，与统计学上中位数的定义不同，当数组长度为偶数时，本实现不会取两个中间元素的平均值，而是直接取其中一个。

// selects the k-th largest element from the data between start and end (end exclusive)
int selectKth(int* data, int s, int e, int k)  // in practice, use C++'s nth_element, this is for demonstration only
{
  // 5 or less elements: do a small insertion sort
  if(e - s <= 5)
  {
    for(int i = s + 1; i < e; i++)
      for(int j = i; j > 0 && data[j - 1] > data[j]; j--) std::swap(data[j], data[j - 1]);
    return s + k;
  }

  int p = (s + e) / 2; // choose simply center element as pivot

  // partition around pivot into smaller and larger elements
  std::swap(data[p], data[e - 1]); // temporarily move pivot to the end
  int j = s;  // new pivot location to be calculated
  for(int i = s; i + 1 < e; i++)
    if(data[i] < data[e - 1]) std::swap(data[i], data[j++]);
  std::swap(data[j], data[e - 1]);

  // recurse into the applicable partition
  if(k == j - s) return s + k;
  else if(k < j - s) return selectKth(data, s, j, k);
  else return selectKth(data, j + 1, e, k - j + s - 1); // subtract amount of smaller elements from k
}

再次展示含噪声的图像：



3x3 中值滤波可以去除其噪声：



由于代码未经优化，较大尺寸的滤波核运行速度相当慢。针对 2D 中值滤波已有更专门、更快速的算法，但这超出了本教程的范围。较大尺寸的滤波结果颇具艺术感，以下展示不同尺寸的效果：

5x5：


9x9：


15x15：

补充说明：上述中值算法实现非常慢。无论是使用 C++ 的 nth_element 函数，还是这里演示用的 "selectKth"，对于求 9 个或 25 个数的中位数，二者带来的收益都很有限。无论某个算法在大 N 情况下的理论复杂度如何，如果只处理某个固定的小规模输入，就需要选用最适合该输入规模的方案。

如果要实现 3x3 中值滤波，最快的方案是使用一个大小为 9 的硬编码排序网络，取其中间输出即为中位数。然后对每个输出像素，将其对应的 9 个输入像素逐颜色通道地应用该网络。硬编码的优势在于算法无需包含依赖输入规模的条件判断（条件判断，如 if 语句和 for 循环的条件，对 CPU 来说非常慢，因为它们会打断流水线）。本文不提供相关代码，因为高效的实际实现超出了本教程的范围。如有兴趣，可以查阅"排序网络"（sorting network）相关资料——它是一种针对两个数的大小关系进行硬编码交换操作的序列，选用已被证明对所需输入规模最优的方案。由于我们不需要完全排序，只需取中位数，因此可以去掉所有不影响中间输出元素的交换操作，并将那些只有一个输出贡献于中间输出元素的交换操作替换为 min 或 max。这将给出理论上最快的实现，在此基础上的进一步加速只能依靠并行性和/或更优的 CPU 指令。

结语

本文提供了对图像应用卷积滤镜的代码，并展示了几种不同的滤镜及其效果。这些只是图像滤镜的基础入门内容，使用更大的滤波核并进行充分调整，可以获得效果更好的滤镜。

傅里叶变换相关文章介绍了在频域中对图像进行滤波的另一种方式，其中讨论了低通滤波、高通滤波和带通滤波。
版权所有 © 2004-2018 Lode Vandevenne。保留所有权利。