卷積與平滑濾波器的圖像處理應用

卷積的介紹

卷積（convolution）是泛函分析里的一個概念，不過泛函分析一般都是數學系才學的，計算機系的學生大多在概率統計課本里了解到。它分為兩種形式，一個是離散形式，一個是連續（積分）形式。在圖像處理中我們更關心離散卷積，不過也先看看積分形式的卷積。現在假設我們有兩個函數 $f(x)$ 和 $g(x)$ ，這里 $g(x)$ 又叫做平滑函數或者卷積核，那么它們在連續空間的卷積是：

$(f*g)(x)=\int_{-\infty}^{\infty}f(t)g(x-t)dt$

一般我們有一個這樣的結論，就是當 $f(x)$ 經過足夠多次相同平滑函數 $g(x)$ 卷積，就會足夠接近高斯函數，也就是正態分布的函數形式。卷積就是一種平滑操作，這說明高斯函數就是“最平滑的函數”。引入熱力學中熵的概念，高斯函數就是擁有最高熵的函數，最穩定的狀態，以至于自然界大多數的統計規律都呈現出正態分布：

$((\cdots((f*g)*g)\cdots)*g)(x) \rightarrow \frac 1{\sigma\sqrt{2\pi}} e^{-x^2/{\sigma^2}}$

下面介紹離散形式的卷積。這卷積，首先是由有限項的多項式體現。神奇的是，而它們的乘積就是卷積。首先我們設有兩個多項式 $p = a_0 + a_1 x + a_2 x^2$ 以及 $q = b_0 + b_1 x + b_2 x^2 + b_3 x^3$ 。計算它們的乘積：

$\begin{align*} r = p\cdot q &= (a_0 b_0) \\ &+ (a_0 b_1 + a_1 b_0) x \\ &+ (a_0 b_2 + a_1 b_1+ a_2 b_0) x^2 \\ &+ (a_0 b_3 + a_1 b_2 + a_2 b_1) x^3 \\ &+ (a_1 b_3 + a_2 b_2) x^4 \\ &+ (a_2 b_3) x^5 \end{align*}$

再引入離散形式卷積（向量卷積）的定義，大家比較一下這個定義和上面多項式的計算。稍微說明一下，中括號的意義是p[n]代表向量第n個元素。將兩個多項式的系數寫成向量形式然后進行向量卷積，也就是例如 $p = [a_0, a_1, a_2]$ ，而沒定義的地方當作0。可以發現，兩者是完全一致的：

$\begin{align*} (p * q)[n] &= \sum_{m=-\infty}^\infty p[m]\cdot q[n-m] \\ r[1] &= \sum_{m=0}^1 p[m]\cdot q[1-m] &&= a_0 b_1 + a_1 b_0 \\ r[2] &= \sum_{m=0}^2 p[m]\cdot q[2-m] &&= a_0 b_2 + a_1 b_1 + a_2 b_0 \\ &\cdots \end{align*}$

知道了多項式的乘積就是其相應的卷積，我們甚至可以直接得出兩個冪級數卷積的結果。因為泰勒級數就是冪級數的一種，所以我們可以將幾乎所有的連續函數轉換成離散形式，避免了繁復的積分運算：比如我們希望得到 $r(x) = p(x) * q(x)$ ，其中 $p(x) = \sum a_i x^i,\ q(x) = \sum b_i x^i$ ，只需要簡單地計算這兩個級數的柯西乘積，所得結果就是 $r(x)$ 的卷積。當然了，這是后話，與本文的主題無關。

卷積與圖像處理

在開始講圖像處理之前，我希望先理解一下卷積的整個過程是怎樣的。從上面的公式看得還是有點懵懵懂懂，從直覺上去理解一下很有必要。觀察卷積的公式以及下面的圖片，這個過程可以看作，當你想求一個r[n]的時候：

你先把卷積核q疊在p上面，盡量使左端靠近（如果左對齊就再好不過了），然后看看在[0, n]內p, q重疊的部分是從哪里到哪里，分別寫成向量，那么r[n]就等于其中一個向量與另一個向量的逆序的內積。

比如當n = 2時，兩個向量是[a_0, a_1, a_2]和[b_2, b_1, b_0]；n = 4時，兩個向量是[a_1, a_2, a_3, a_4]和[b_3, b_2, b_1, b_0]。至于求內積，一定難不倒你。下圖說明了這一點：

$\begin{align*} && a_0 && a_1 && a_2 && a_3 && a_4 \\ && a_0b_0 && a_0b_1 && a_0b_2 && a_0b_3 \\ && && a_1b_0 && a_1b_1 && a_1b_2 && a_1b_3 \\ && && && a_2b_0 && a_2b_1 && a_2b_2 && a_2b_3 \\ && && && && a_3b_0 && a_3b_1 && a_3b_2 && a_3b_3 \\ && && && && && a_4b_0 && a_4b_1 && a_4b_2 && a_4b_3 \\ \hline && c_0 && c_1 && c_2 && c_3 && c_4 && c_5 && c_6 && c_7 \end{align*}$

上面是對某一點上卷積的理解。對整個域的卷積，則可以看成是將卷積核（除了開頭幾個外）不停向右移動，每移動一格就將重疊部分拿出來求內積。

這時我們可以把圖像處理和卷積聯系起來了。圖像處理是，將一副“源圖像”（Source），通過一些算法，變成一副“目標圖像”（Destination）。當我們進行平滑處理的時候，用到一個叫做濾波器（filter）的東西，也叫做濾鏡。想想我們現實生活中放大鏡是怎么用的：拿著放大鏡，從報紙的左上角開始，一直掃啊掃到右下角，掃的過程中一直望著放大鏡和報紙的重疊區域（其實就是望著放大鏡，因為它比報紙小多了），這樣你就瀏覽完了一張放大過的報紙。平滑濾鏡也是同樣的使用方法，從源圖的左上角開始掃到右下角，掃的過程中一直取出重疊部分進行內積計算，然后將結果存放到目標圖像中 —— 顯然這個操作跟卷積是一致的，只不過定義在二維空間內。

為了方便量化表示，我們把圖像抽象成定義在 $R \cap [0, 1]$ 數環內的二維矩陣，其意義是灰度值，顏色信息我們暫且忽略。卷積核，也就是濾波器同樣也是定義在 $R \cap [0, 1]$ 內的二維矩陣。這樣，二維的卷積我們這樣定義它的離散形式：

$\text{Dest}[i, j] = \sum_{y=-\infty}^\infty \sum_{x=-\infty}^\infty \text{Src}[y,x] \cdot \text{Ker}[i - y, j - x]$

我們的卷積核大小并不是無限的，它一個半徑r，這樣它的大小就是2r+1。規定了這個r使得，當|x| > r 或 |y| > r，都有Ker[y, x] = 0。規定過大小之后，由 |i-y| < r; |j-x| < r得到 i-r < y < i+r; j-r < x < j+r。同時我們規定Dest和Src的大小是 $m \times n$ 。于是我們得到了濾波器的算法：

$\text{Dest}[i, j] = \sum_{y=\max\{0,i-r\}}^{\min\{i+r,n\}} \left( \sum_{x=\max\{0,j-r\}}^{\min\{j+r,m\}} \text{Src}[y, x] \cdot \text{Ker}[i - y, j - x] \right)$

高斯濾波器

二維的高斯函數（俗稱避孕套函數）

高斯濾波器是最常用的平滑濾波器之一，在Photoshop里面它被用作高斯模糊濾鏡。高斯濾波器的定義很經典，就是簡單地把正態分布離散開來。二維形式只是單純把x替代成(x² + y²)，然后修改系數令實數域上的積分為1：

$\begin{align*} \text{Ker}_1[x] &= \frac 1{\sigma\sqrt{2\pi}} e^{-x^2/{\sigma^2}} \\ \text{Ker}_2[i, j] &= \frac 1{2\sigma^2\pi} e^{-(i^2+j^2)/{\sigma^2}} \end{align*}$

也許你已經發現了一個這樣的規律，這一規律，這在更高維上仍然是滿足的，也就是在連續空間里同樣滿足。這將成為我們優化算法的關鍵。將這個規律代回到二維離散卷積的公式里，因為y在第二個連加中相當于常數系數可以提出來，我們發現：

$\begin{align*} \text{Ker}_2[i, j] &= \text{Ker}_1[i] \cdot \text{Ker}_1[j] \\ \text{Dest}[i, j] &= \sum_{y=-\infty}^\infty \sum_{x=-\infty}^\infty \text{Src}[y,x] \cdot \text{Ker}_2[i - y, j - x] \\ &= \sum_{y=-\infty}^\infty \sum_{x=-\infty}^\infty \text{Src}[y,x] \cdot \text{Ker}_1[i-y] \cdot \text{Ker}_1[j-x]\\ &= \sum_{y=-\infty}^\infty \left( \sum_{x=-\infty}^\infty \text{Src}[y,x] \cdot \text{Ker}_1[j-x]\right) \text{Ker}_1[i-y] \end{align*}$

如果x連加所表示是卷積是從右上角開始按照文字書寫順序從左到右，然后從上到下的順序進行一維卷積，那么y連加表示的卷積就是先從上到下，再從左到有的順序卷積。在OpenCV提供的數據結構的基礎上，不用imgproc提供的算法，我寫了一個示例：

1 // cflags: -lopencv_highgui -lopencv_core
2
3 #include <iostream>
4 #include <cmath>
5 #include <opencv2/highgui/highgui.hpp>
6
7 using namespace cv;
8 using namespace std;
9
10 const char* title = "gaussian-filter";
11
12 Mat kernelMatrix(int radius, double sigma)
13 {
14         int d = radius * 2 + 1;
15         Mat kernel(2, d, CV_64F);
16
17         double coef = 0;
18         for(int i = 0; i <= radius; i++) {
19                 // f(x) = 1/(sigma * sqrt(2 pi)) * e ^ -x^2/(2 s^2)
20                 int dx = i - radius;
21                 double dx_2 = dx * dx;
22                 double w = pow(M_E, - dx_2 / (2 * sigma * sigma));
23
24                 coef += w;
25                 kernel.at<double>(0, i) = w;
26                 kernel.at<double>(0, d - i - 1) = w;
27                 // when you used values from i to j (j>i), the sum of them is:
28                 // kernel[1, j] - (i ? kernel[1, i-1] : 0)
29                 kernel.at<double>(1, i) = coef;
30         }
31
32         for(int i = radius + 1; i < d; i++) {
33                 coef += kernel.at<double>(0, i);
34                 kernel.at<double>(1, i) = coef;
35         }
36
37         return kernel;
38 }
39
40
41 void convolution(const Mat& img, const Mat& kernel, Mat& output, bool t = true)
42 {
43         for(int y = 0, x = 0; y < img.rows; x = (++x<img.cols)? x : (y++, 0)) {
44                 Vec3d r(0, 0, 0);
45
46                 int ideal = x - int(kernel.cols / 2),
47                     ran_beg = max(ideal, 0) - ideal,
48                     ran_end = min(ideal + kernel.cols, img.cols) - ideal;
49
50                 for(int i = ran_beg; i < ran_end; i++) {
51                         double weight = kernel.at<double>(0, i);
52                         Vec3b pixel = img.at<Vec3b>(y, ideal + i);
53
54                         r[0] += pixel[0] * weight;
55                         r[1] += pixel[1] * weight;
56                         r[2] += pixel[2] * weight;
57                 }
58
59                 double coef = kernel.at<double>(1, ran_end - 1);
60                 if(ran_beg) coef -= kernel.at<double>(1, ran_beg - 1);
61
62                 output.at<Vec3b>(t?x:y, t?y:x) = Vec3b(
63                         saturate_cast<uchar>(r[0]/coef),
64                         saturate_cast<uchar>(r[1]/coef),
65                         saturate_cast<uchar>(r[2]/coef));
66         }
67 }
68
69
70 int main()
71 {
72         namedWindow(title, WINDOW_AUTOSIZE);
73
74         const int r = 10;
75
76         Mat img = imread("ai-sample.jpg"),
77             kernel = kernelMatrix(r, (r - 1) * 0.3 + 0.8);
78
79         Mat product_v = Mat(img.cols, img.rows, img.type());
80         Mat product_h = Mat(img.rows, img.cols, img.type());
81
82         convolution(img, kernel, product_v);
83         convolution(product_v, kernel, product_h);
84
85         imshow(title, product_h);
86         for(; waitKey(0) > 0;);
87
88         destroyWindow(title);
89 }

渲染效果圖

posted on 2015-08-12 00:35 Shihira 閱讀(3229) 評論(0) 編輯收藏引用所屬分類: 圖形編程

只有注冊用戶登錄后才能發表評論。
【推薦】100%開源！大型工業跨平臺軟件C++源碼提供，建模，組態！

相關文章: 卷積與平滑濾波器的圖像處理應用 AGG入門（八） - 渲染 AGG入門（七） - 頂點源 AGG入門（六） - 練習和細節 AGG入門（五） - 基礎渲染器 AGG入門（四） - 渲染緩存和混合器 AGG入門（三） - 渲染器介紹 AGG入門（一） - 配置開發環境 AGG入門（二） - 平臺支持

網站導航: 博客園 IT新聞 BlogJava 博問 Chat2DB 管理

卷積與平滑濾波器的圖像處理應用

卷積的介紹

卷積與圖像處理

高斯濾波器

導航

統計

公告

留言簿(2)

隨筆分類

搜索

最新隨筆

最新評論

閱讀排行榜

評論排行榜