（轉）卷積在C++中的實現
@author：http://www.cnblogs.com/yabin/p/6664011.html
白袍小夫轉，分類在Graphy,後面再放一個C++And HLSL or GLSL的。
? ?
1 信號處理中的卷積
無論是信號處理、圖像處理還是其他一些領域，我們經常會在一些相互關聯的數據處理中使用卷積。卷積可以說是算法中一個非常重要的概念。這個概念最早起
源於信號處理之中。
假設對於一個線性系統其在單位脈沖δ(t)的響應下，輸出為h(t).那麽在Aδ(t)的響應下輸出為Ah(t).而所有的信號都可以用δ(t)乘以一個系數的和來表示。即
.於是對於線性系統而言，我們可以將當前及過去所有時刻信號（看出無數個脈沖乘以系數）產生的輸出進行疊加來得到當前時刻的輸出。即對於任意時刻
的輸出為
.關於卷積的來源和其物理意義教課書和網上都有很多的介紹，我這裏也不再啰嗦了，卷積的物理意義主要體現在如下四個方面：
- 信號響應的計算
- 和多項式的關系
- 濾波器的應用
- 圖像處理中的應用
2 matlab中卷積函數
matlab中提供了兩個常用的卷積函數conv和conv2.分別用來計算一維卷積和二維卷積。這兩個函數都是對離散的數據進行卷積，接下來我們就分別討論這兩種卷積。
2.1一維卷積conv函數
通過對之前信號處理中卷積的介紹，我們將其進行離散化處理可以得到如下離散的一維卷積表達式：
.其中u
卷積和多項式
其實兩個多項式之間的乘法我們可以看作是它們系數矩陣之間的卷積。例如
之間的乘積就可以看成u=[1 0 1]和v = [2 7]之間的卷積。這時卷積後的結果序列長度為size(u)+size(v)-1.對應matlab代碼如下：
clear
clc
u = [1 0 1];
v = [2 7];
w = conv(u,v)
執行結果為w=[2 7 2 7],表示u,v得到的對應多項式為：
.由於卷積滿足交換律，所以
卷積和濾波器
從另一方面我們也可以將卷積看作是一個濾波的過程（FIR），將u看成被濾波的信號，將v看成一個濾波器。此時濾波後得到的結果序列長度和u相同。
在v裏面有一個濾波中心點，在卷積乘加的過程中，中心點對應"時刻"的濾波值始終和對應卷積結果時刻的值相乘，通過平移濾波器序列改變中心點的位
置，後面圖示會直觀的解釋這一概念。對於信號處理由於只能知道過去時刻的信號，所以中心點應該選擇v的第一個元素上（即v[0]為中心點）。對於已
知整個序列的信號（如圖像處理）則中心點的位置一般選擇在size(v)/2+1(如果不能整除則選在這個值的左邊或右邊)上。特殊的如果選擇的濾波器各個值
都相等且和為1則這是一個均值濾波器，下面是一個中心點在v[0]，size(v)=3的一個均值濾波器濾波示意圖：
上圖中為了方便觀測卷積濾波器的物理意義，這裏直接求w[2]和w[3].圖中的連線表示被連接的兩個數需要相乘。可以看出此時v就相當於一個濾波因子對
u的對應值進行濾波,我們將v看出一個單位脈沖響應，因此也稱為FIR(有限脈沖響應)濾波器。
對於圖像處理或者已知整個信號序列值的場合，濾波器的中心點往往選在在v的中心（或左右），matlab中conv函數在調用的過程中如果加入了same參
數則表示把該卷積看成一個濾波器，且濾波器的中心點在size(v)/2+1或者round(size(v)/2)(v長度為奇數時)。下面是matlab中將卷積看成濾波器且濾
波器中心點在中間的計算代碼和計算過程的示意圖：
clear
clc
u = [1,2,3,4,3,1];
v = [1/3 1/3 1/3];
w = conv(u,v,‘same‘)
執行結果：
w =
1.0000 2.0000 3.0000 3.3333 2.6667 1.3333
我們發現將卷積看成是一個濾波器（濾波算子）後，在乘加的過程中我們總需要將v進行交叉。這也是為什麽在圖像處理中為什麽各種算子和對應的卷積核總是
倒轉了180°的原因。針對這個問題有兩種解決方法，一個是改變在濾波過程中的卷積公式，將其中的減改成加即定義濾波卷積公式為：
。第二個解決方法是，就接受這個相反的事實，在畫相關示意圖時手動將它調轉180°。我們對上面兩個圖進行反轉也行看的能更順眼一些：
? ?
在信號濾波中濾波器中心點一般選擇在最左邊，而在圖像處理中濾波器中心一般選擇在中心位置。其他位置的使用較少。事實上，可以通過對上節通用的
卷積結果進行移位和截斷能得到任意中心位置的卷積結果。
除了same參數以外，matlab還提供了valid參數用來計算有效的濾波點數。所謂的有效點數就是中心點在濾波器中心的條件下，能全部使用濾波器數
據進行計算的數據點。比如上面中的w[0]就不是一個有效的數據點。用圖像處理的角度來說就是，算子在計算邊界時沒有足夠的數據進行計算然後舍棄
這些點。計算後的w長度=max(length(u)-length(v)+1,0).下面是一個matlab示例，可以體會一下：
clear
clc
u = [1,2,3,4,3,1];
v = [1/3 1/3 1/3];
w = conv(u,v,‘valid‘)
執行結果如下：
w =
2.0000 3.0000 3.3333 2.6667
2.2二維卷積conv2函數
和前面一維卷積一樣，我們可以將卷積的概念推廣到二維卷積上。二維卷積最常用的地方就是圖像的濾波處理。這裏我們將卷積所用的濾波器稱為卷積核，
一些常用的卷積核和我們又稱之為某某算子。比如常用於邊沿檢測的Sobel算子。事實上，由於卷積的概念來自於數字信號處理，所以在圖像處理中算子
和卷積核剛好關於中心點成中心對稱（前面我們已經解釋了）。但是不引起混淆的情況下，我默認它們是一樣的不再進行區分了。
首先我們我們先看一下二維卷積的定義，可以使用doc conv2參考matlab中的說明文檔，但需要註意matlab中角標是從1開始而不是0開始，公式會
有些不一致。
兩個二維矩陣a,b的卷積結果c可以用如下的公式表示：
對於C=conv2(A,B),如果有[ma,na]=size(A),[mb,nb]=size(C),[mc,nc]=size(C).那麽mc=max([ma+mb-1],ma,mb),nc=max(na+nb-1,na,nb).
對於二維卷積結果中的一點c(m,n),它其實就是將a(m,n)周圍的點和b中所有的點對應相乘然後求和（假設b矩陣的維數小於a）。在圖像處理中a(m,n)
周圍的點是指以a(m,n)為中心周圍的點。具體的意義可以百度"二維卷積"一般都有配圖，我這裏就不畫圖了，我在網上隨便截了個圖，大概就是這個意思：
在matlab在調用二維卷積的時候同樣也提供了三個可選參數:
full: 計算兩個二維矩陣的常規卷積（默認參數），和一維卷積一樣
same: 和一維卷積一樣
valid: 和一維卷積一樣
上面三個參數中當參數為full時表示的意義和信號處理中卷積的意義相同。當參數為same時，對應的卷積公式為：
上面這個公式實際上是在原來的卷積公式中，將卷積核向右，向下分別平移了卷積核行列一半的向下取整個單位。也就是說讓卷積核的中心和計算卷積
的點對應。下面是full參數下和same參數下卷積的物理意義對比，以計算c(0,0)為例。
full參數對應公式：
same參數對應公式:
可以看出full參數計算卷積求c(0,0)時，從圖像上來看首先是將b矩陣旋轉180°(註意b矩陣左上角為(2,2),右下角為(0,0)),然後讓a(0,0)和
b(0,0)對齊，最後對應乘加計算c(0,0)的值，這裏有可能會涉及到a矩陣邊界以外的點。
對應same參數來說大體計算流程並沒有改變多少，只是這裏不是讓a(0,0)和b(0,0)對齊，而是讓a(0,0)和b矩陣的中心點b(mb/2,nb/2)對
齊(行列為奇數時向下取整)，同樣這裏有可能會涉及到a矩陣邊界以外的點。
由於邊界以外的點的情況我們往往並不知道，所以matlab中提供了一個valid參數，該參數下計算卷積時，如果計算某個卷積點c(i,j)用到了a
矩陣邊界外的點，那麽就會去掉這個卷積結果。其對應的公式為：
該公式的思想和same參數下公式的思想相同。只是這裏計算c(0,0)時，將a(0,0)和b(2,2)對齊。由於a(0,0)和b(2,2)都是矩陣在計算時的右
上角元素，所以計算c(0,0)時並不會用到a矩陣邊界外的點。接下來計算的時候只要保證卷積核b不移動到a矩陣外面即可。
除了上面和一維卷積相同的調用形式以外，matlab中還提供了另外一種調用形式，即C=conv2(h1,h2,A).它首先讓A的每一列和h1進行卷積，
然後再讓結果的每一行和h2進行卷積，如果n1=length(h1),n2=length(h2),那麽mc=max([ma+n1-1,ma,n1]),nc=max([na+n2-1],n
a,n2).一般情況下這種調用形式我們用的比較少。
2.3 邊界的擴充
matlab的卷積計算過程中使用的是信號處理中求系統響應的標準計算公式（圖像處理等實際上就是在標準卷積結果中進行截取），這時候針對圖像
處理中的"邊界"以外的點都是認為為0的。然而在計算機圖形學中有時候我們並不認為這是一個合理的假設。比如這幅圖是在一個大的圖像中截取的
我們認為邊界點的以外的點約等於邊界上的點更合適。或者這幅圖是一個周期循環的圖，我們在這裏截取的是一個周期上的圖像，這時候認為邊界以
外的點等於另一邊界上的點更合適。
此外，當卷積核一部分跑到邊界以外的時候，我們想舍去這些點。這時候我們就可以用valid參數進行舍棄。
3 卷積的C++實現
為了使我們的代碼在計算系統響應、多項式乘法、濾波器設計、圖像處理中有更廣泛的通用性，除了實現上面matlab中提供的卷積實現以外，我們還
要增加matlabconv2函數中沒有的邊界處理條件選擇。
可選的邊界條件有：邊界以外的點為0（默認），邊界以外的點和邊界值相等，邊界以外的點和另一邊界上的點組成周期信號。
這裏實現C++卷積算法的時候使用的是C++矩陣計算庫Eigen.
需要註意的是這裏卷積實現的時候將conv和conv2合在一起了。事實上，在matlab中的大部分場合使用conv的地方都可以使用conv2代替，因此
在C++實現裏不再區分conv和conv2.但是需要註意的是，使用conv計算一個行向量和一個列向量的卷積得到的是一個向量，而使用conv2得到的
卻是一個矩陣，在這一點上conv和conv2是有區別的。因此在使用時如果你確保你求的是一維卷積，請保證兩個向量要麽都是行向量要麽都是列向量。
下面是matlab的一個示例：
clear
clc
a = rand(1,3);
b = rand(3,1);
c = conv(a,b)
d = conv2(a,b)
e = conv2(a,b‘)
執行結果：
c =
0.1409
0.4958
0.2941
0.1608
0.0380
d =
0.1409 0.0357 0.0384
0.4601 0.1165 0.1255
0.1391 0.0352 0.0380
e =
0.1409 0.4958 0.2941 0.1608 0.0380
3.1 C++中卷積的簡單實現
有了上面的分析之後，首先我們先編寫一個較為簡單的卷積函數來驗證一下功能。之後再將這個函數封裝成一個類。簡單的C++卷積代碼如下：
#include <iostream>
using namespace std;
enum BoundaryCondition
{
zero
,bound
,period
};
enum Method
{
full
,same
,valid
};
//參數設置
const int A_Row = 4;
const int A_Col = 5;
const int B_Row = 2;
const int B_Col = 3;
BoundaryCondition Bc = zero;
Method method = full;
float A[A_Row][A_Col] = {{1,2,3,4,5},
{2,3,4,5,6},
{3,4,5,6,7},
{4,5,6,7,8}};
float B[B_Row][B_Col] = {{1,2,1},
{2,3,2}};
float** C = 0;
int cR;
int cC;
template <typename T>
T GetA_Ele(int row,int col);
template <typename T>
void conv(const T a[][A_Col],const T b[][B_Col],T**& c);
int main()
{
conv(A,B,C);
for(int i = 0;i < cR;i++)
{
for(int j = 0;j < cC;j++)
cout<<C[i][j]<<"\t";
cout<<endl;
}
delete[] C;
}
//計算A和B的卷積
template<typename T>
void conv(const T a[][A_Col],const T b[][B_Col],T**& c)
{
int offsetR = 0;
int offsetC = 0;
switch(method)
{
case full:
cR = A_Row + B_Row - 1;
cC = A_Col + B_Col - 1;
break;
case same:
cR = A_Row;
cC = A_Col;
offsetR = B_Row/2;
offsetC = B_Col/2;
break;
case valid:
cR = A_Row - B_Row + 1;
cC = A_Col - B_Col + 1;
if((cR < 1)|(cC < 1))
return;
offsetR = B_Row - 1;
offsetC = B_Col - 1;
break;
default:
return;
}
c = new T*[cR]; //給二維數組分配空間
for(int i = 0;i < cR;i++)
c[i] = new T[cC];
for(int i = 0;i < cR;i++)
{
for(int j = 0;j < cC;j++)
{
c[i][j] = 0;
for(int k1 = 0;k1 < B_Row;k1++)
{
for(int k2 = 0;k2 < B_Col;k2++)
c[i][j] += b[k1][k2]*GetA_Ele<float>(i - k1 + offsetR,j - k2 + offsetC);
}
}
}
}
//根據邊界條件獲取A矩陣的元素
template <typename T>
T GetA_Ele(int row,int col)
{
switch(Bc)
{
case zero: //索引超出界限認為0
if((row < 0)|(row > A_Row)|(col < 0)|(col > A_Col))
return 0;
case bound: //超出索引部分和邊界值相等
if(row < 0)
row = 0;
else if(row >= A_Row)
row = A_Row - 1;
if(col < 0)
col = 0;
else if(col >= A_Col)
col = A_Col - 1;
return A[row][col];
case period:
while((row < 0)|(row >= A_Row))
{
if(row < 0)
row += A_Row;
else
row -= A_Row;
}
while((col < 0)|(col >= A_Col))
{
if(col < 0)
col += A_Col;
else
col -= A_Col;
}
return A[row][col];
default:
return T(0);
}
}
? ?
上面函數執行結果為:
1 4 8 12 16 19 15
4 14 26 37 48 56 43
7 22 37 48 59 67 51
10 30 48 59 70 78 59
12 35 55 66 77 85 64
由於沒有使用矩陣運行庫且算法較為簡單，所以參數設置起來比較麻煩。需要手動給出矩陣的行和列。這裏方法中的三個參數full,same,valid和
MATLAB中相同。此外為了更方便的進行圖像處理，有提供了3種可選的卷積邊界條件。由於在卷積計算的過程中，對A矩陣的取值有可能超出
索引，因此必須使用一個Get方法進行封裝。
為了使用matlab對我們的算法進行驗證，這裏取Bc = zero;，然後判斷三個不同的方法計算結果是否和matlab中相同。matlab對應代碼如下：
clear
clc
a = [1,2,3,4,5;...
2,3,4,5,6;...
3,4,5,6,7;...
4 5 6 7 8];
b = [1 2 1;2 3 2];
c = conv2(a,b,‘full‘)
有了上面的簡單移植代碼做鋪墊，我們就可以編寫一個較為正式的conv模板來供我們以後使用了。對應的C++類模板conv使用方法如下
（源碼見文章的最後面）：
#include "Eigen\core"
#include <iostream>
#include "Conv.h"
using namespace std;
using namespace Eigen;
using namespace ConvSpace;
int main()
{
MatrixXd a(4,5);
MatrixXd b(2,3);
a << 1, 2, 3, 4, 5,
2, 3, 4, 5, 6,
3, 4, 5, 6, 7,
4, 5, 6, 7, 8;
b << 1, 2, 1,
2, 3, 2;
try
{
Conv<MatrixXd> myConv(a, b); //第一種調用形式
//Conv<MatrixXd> myConv(a, b,ConvType::Full,BoundaryConditon::Period); //不使用默認參數
cout << myConv.Eval() << endl; //直接輸出結果
//const MatrixXd &r = myConv.Eval(); //第二種獲取結果的形式
//cout << r << endl;
//Conv<MatrixXd> myConv; //第二種調用形式
//myConv.SetA(a);
//myConv.SetB(b);
//cout << myConv.Eval() << endl;
}
catch(ConvError ce)
{
cout<<ce.GetMessage()<<endl;
}
int c;
cin >> c; //程序暫停
}
執行結果如下:
該算法需要使用到Eigen矩陣運行庫，因此需要包含Eigen相關文件和命名空間。此外這裏將該算法在ConvSpace命名空間中實現，因此在使用
函數之前必須要包含該命名空間和相關文件。
該函數大體上有兩種調用形式，一種是在構造函數中給出被卷積矩陣然後進行計算，另外一種是使用空構造函數，然後通過Set方法設置相關矩陣，
最後計算。第一種方式調用比較簡單如下：
Conv<MatrixXd> myConv(a, b); //第一種調用形式
cout << myConv.Eval() << endl; //直接輸出結果
在使用上面代碼時應該將所有的代碼放到一個try語句塊中，並且捕獲ConvError異常。當出現參數配置錯誤或其他錯誤時，會拋出該異常通過查看
該異常的說明可以判斷異常出現的原因。
上面使用了默認的邊界條件zero和默認的卷積方法full。當然也可以自己選擇邊界條件和卷積方法（如上面註釋）。
前面在介紹原理的時候提到過，matlab中conv2還有另外一種調用形式，是計算兩個向量和一個矩陣的卷積。這裏並沒有實現這個算法。實際上向量
也是一種特殊的矩陣，有興趣的同學可以在代碼中補全這一部分內容，最好將補全後的代碼再分享出來。
源代碼
4 關於圖像卷積的思考
由前面卷積的定義，我們可以發現把卷積的概念引入圖像處理中並不能直接的看出其物理意義，主要由以下三個方面的移植問題：
- 它將其中的一個輸入序列"向前"進行了平移，這就要求計算某個點k的卷積結果時，需要知道k+i"時刻"的輸入序列值。而這在時間序列作為輸入
- 時往往是不可能做到的。
- 有圖像卷積的圖示可以知道，卷積實際上是將卷積核先進行中心對稱旋轉後，再和另外一個矩陣對應乘加。由於並不是直接乘加因此沒有濾波算子
- 的所表達的物理意義明確。
- 在圖像處理中卷積往往有一個"邊界的概念"。而在數字信號處理中並不會有這種概念，因為在數字信號零狀態下，線性系統沒有輸入就沒有輸出。
由上面幾點可以看出，卷積的概念推廣到圖像處理中，已經有何很大的改變，也賦予了卷積更多的意義。由於卷積的概念在圖像處理中並實現完全 - 的"本地化"給其賦予一個直觀的物理意義。因此我們可以對卷積的定義稍作修改，定義一個在圖像處理中意義更加明確的圖像卷積如下：
其對應的示意圖如下：
上面定義的圖像卷積和數組信號處理中的卷積略有不同，這樣將卷積的概念移植到圖像處理中並不會將卷積核進行顛倒，註意b矩陣左上角為b(0,0).
而且可以通過兩個可調的偏移量l和m調整卷積核在乘加過程中的中心位置。上面的圖像卷積公式有著更加直觀的物理意義。
當然有利有弊，在數字信號處理中我們關於卷積的討論已經有了很多非常成熟的結論。因此如果我們使用的是普通的卷積概念，這些討論和結論在圖像
處理中是可以直接使用的。如卷積等於頻域中的乘積、卷積的交換律和分配率等。雖然上面卷積通過討論也能得到相似的性質，但是畢竟增加了許多的工作量。
總之，有利有弊，如果我們需要直觀的在圖像處理中使用卷積的話，上面的公式在某種方面也有助於我們對基本卷積公式的理解。
? ?
來自 <http://www.cnblogs.com/yabin/p/6664011.html>
（轉）卷積在C++中的實現