由于項目的需要在我還沒有搞清楚Opencl架構的情況下就開始去用opencl方式去優化一個圖像處理算法--bilateral_filter，可想而知其過程的艱辛，好在昨天已經搞定收工了，感謝這一個月的過程讓我從一個opencl路人變成了親戚，特此寫一篇部落格記錄下一些心得。廢話少說現在開始。
拿到這個任務是我首先感覺相對我之前寫代碼還是比較高大上的，是以心裡也暗暗下決心搞定這個問題，于是先去了解雙邊濾波算法的原理，在網上資料一大堆，看了幾篇部落格也就很快了解了在此總結一下。
雙邊濾波
雙邊濾波是一種非線性濾波器，它可以達到保持邊緣、降噪平滑的效果。和其他的濾波算法一樣，雙邊濾波也是采用了權重平均的方法，隻不過這個權計算的比較複雜，它不僅考慮了像素範圍中周圍像素到中心像素的歐氏距離對中心像素的影響（高斯濾波）還考慮了像素範圍域中的輻射差異（例如卷積核中像素與中心像素之間相似程度、顔色強度，深度距離等，對于C8灰階圖像來說就是灰階值之差，可以了解為頻域，是以在邊緣時值比較大，平坦的地方值比較小）。
數學公式表達：
這裡的權值w就是通過空域和像素範圍域來确定的，公式如下：
像素範圍域：
空域：
于是整合在一起就是：
綜合結果：在圖像的平坦區域，像素值變化很小，對應的像素範圍域權重接近于1，此時空間域權重起主要作用，相當于進行高斯模糊；在圖像的邊緣區域，像素值變化很大，像素範圍域權重變大，進而保持了邊緣的資訊。
盜張圖：來自 https://blog.csdn.net/piaoxuezhong/article/details/78302920
在了解時可以根據公式考慮兩種極端情況，即中心點和周圍點像素值一樣時（平坦地帶）和中心點255周圍是0時（邊緣地帶）。
代碼實作：
for(int i=0;i<gray.rows;i++)
{
for(int j=0;j<gray.cols;j++)
{
fij = *(gray.data + i*gray.cols + j);
for(int K = -ksize / 2; K <= ksize / 2; K++)
{
for(int L = -ksize / 2; L <= ksize / 2; L++)
{
if(((i+K )<0) || ((j + L)<0)|| ((j + L)>(gray.cols-1))|| ((i + K)>(gray.rows-1)))
{
fkl = 0;
}
else
{
fkl = *(gray.data + (i+K)*gray.cols +(j+L));
}
dkl = -(K*K + L*L) / (2 * sigma_d*sigma_d);
rkl = -(fij - fkl)*(fij - fkl) / (2 * sigma_r*sigma_r);
wkl = exp(dkl + rkl);
numerator += fkl * wkl;
denominator += wkl;
}
}
gij = numerator / denominator;
//雙邊濾波後再做一個融合
gij = fij*alpha + gij*(1.0 - alpha);
blurImage.data[i*blurImage.cols + j]= (uchar)gij;
}
}
以上簡單的描述一下原理，下面開始介紹opencl的實作。
Opencl實作
環境：Cyclone V GT 150k le + Atom E3930 +Quatus aoc 編譯
編譯過程比較慢大概2個半小時一次，還是比較痛苦的，下面就總結下opencl的從了解到使用的過程及相關知識脈絡和優化心得。
Opencl（全稱Open Computing Language，開放運算語言），是第一個面向異構系統通用目的并行程式設計的開放式、免費标準，也是一個統一的程式設計環境。這裡注意它是一種程式設計環境，具體還是用C來編碼。
OpenCL裝置有兩部分組成，主控端和OpenCL裝置，主控端負責整體流程控制，一般為CPU，OpenCL裝置主要負責資料運算操作。下面是從實體角度劃分 計算裝置（CD）、計算單元(CU)、處理單元（PE）關系。
主機上的OpenCL 應用程式送出指令(command queue)給裝置中的處理單元以執行計算任務(kernel)。計算單元中的處理單元會作為SIMD 單元（執行指令流的步伐一緻）或SPMD 單元（每個PE 維護自己的程式計數器）執行指令流。
上面是Host端的初始化調用流程和Running階段處理事件的調用過程。
基于Opencl算法優化
優化過程總結幾點優化方法：
1.因為是處理圖像中每一個像素，并且每個像素處理時是基于原圖像互相獨立的，是以我們就可以把每個像素作為一個線程來進行處理，使用NDrange的方式。類似于系統的多線程處理。
2.對于一些乘除法運算比較資源，是以可以考慮用移位的方式。
3.對于那些計算量大但是可以知道範圍的可以提前計算好存入數組中，fpga執行時就可以直接搜尋數組獲得值就可以了。
4.向量計算能夠增加計算效率，如果可以嘗試用向量的方法。
5.對于那些累加計算可以嘗試并行加，這樣可以縮短累加時間。
6.向量讀記憶體會減小讀記憶體時間。
基于opencl的雙邊濾波算法代碼kernel size(11,11)：
kernel void bilateral_filter(const __global uchar16 * restrict sdata,__global unsigned char * restrict dest)
{
//float alpha = 0.4;
//int ksize = 11;
const int Col = get_global_id(0);
const int Row = get_global_id(1);
const int width = get_global_size(0);
const int height = get_global_size(1);
const int index_1 = Row*width+Col;
//uchar fij = sdata[index_1];
uchar fkl[11][11];
int colindex = index_1 >> 4;
int Ltmp = index_1 - (colindex << 4);
c32 vfkl;
for(int K = 0; K < 11; K++)
{
#pragma unroll
for(int i = 0;i<2;i++)
{
vfkl.v[i] = sdata[ofs_space[K]+colindex+i];
}
#pragma unroll
for(int L = 0; L < 11; L++)
{
fkl[K][L] = vfkl.d[L+Ltmp];
}
}
uchar fij = fkl[5][5];
float numerator = (float)(0.0f);
float denominator = (float)(0.0f);
//#pragma unroll
for(int K = 0; K < 11; K++)
{
f16 wkl;
f16 num;
float rkl_a[11];
float dkl_a[11];
#pragma unroll
for(int L = 0; L < 11; L++)
{
uchar temp;
if(fij > fkl[K][L])
temp = fij - fkl[K][L];
else
temp = fkl[K][L] - fij;
rkl_a[L] = value_space[temp];
dkl_a[L] = dis_space[K][L];
}
float16 rkl = (float16)(rkl_a[0],rkl_a[1],rkl_a[2],rkl_a[3],rkl_a[4],rkl_a[5],
rkl_a[6],rkl_a[7],rkl_a[8],rkl_a[9],rkl_a[10],0,0,0,0,0);
float16 dkl = (float16)(dkl_a[0],dkl_a[1],dkl_a[2],dkl_a[3],dkl_a[4],dkl_a[5],
dkl_a[6],dkl_a[7],dkl_a[8],dkl_a[9],dkl_a[10],0,0,0,0,0);
wkl.v = rkl*dkl;
float16 fkl16 = (float16)(fkl[K][0],fkl[K][1],fkl[K][2],fkl[K][3],fkl[K][4],fkl[K][5],
fkl[K][6],fkl[K][7],fkl[K][8],fkl[K][9],fkl[K][10],0,0,0,0,0);
num.v = fkl16 * wkl.v;
float numtmp1[6];
float wkltmp1[6];
#pragma unroll
for(int i = 0;i<5;i++)
{
numtmp1[i] = num.d[i] + num.d[10 - i];
wkltmp1[i] = wkl.d[i] + wkl.d[10 - i];
}
numtmp1[5] = num.d[5];
wkltmp1[5] = wkl.d[5];
float numtmp2[3];
float wkltmp2[3];
#pragma unroll
for(int i = 0;i<3;i++)
{
numtmp2[i] = numtmp1[i] + numtmp1[5 - i];
wkltmp2[i] = wkltmp1[i] + wkltmp1[5 - i];
}
float numtmp3 = numtmp2[0] + numtmp2[1] + numtmp2[2];
float wkltmp3 = wkltmp2[0] + wkltmp2[1] + wkltmp2[2];
numerator += numtmp3;
denominator += wkltmp3;
}
float gij = numerator / denominator;
dest[index_1] = (uchar)gij;
}