由于项目的需要在我还没有搞清楚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;
}