前言

影像处理会依照[1]的章节进行，但不直接使用OpenCV而是实作原理，主要了解影像处理原理的特性，虽然速度不比OpenCV好，但在看OpenCV原始码你会更好理解运作原理并且可以加以修改自己所需要模式，这样做可让自己更有机会发展出更多不同的处理方式，这次会介绍4种平滑滤波器平滑滤波、高斯滤波、中值滤波和双边滤波，首先来介绍一些基本概念。

---

高斯分布(常态分布)

相信大家对常态分布不陌生如下图。[3]提到在自然界或数学当中，一个不明的随机变量我们会使用常态分布来去表示或计算，而在人工智慧当中高斯分布也是其中的种要角色之一。主要将它当做计算一个机率的分布，接着介绍它的数学公式。

来源：[2]

一维公式： 来源[3]
这里sigma在控制高低如下图一，mi在控制偏移量如下图二。但通常影像处理的维度是二维，因此要带入二维常态分布。

图一

图二

二维公式： 来源[4]
二维公式推导为两个机率分布相乘，一维公式 * 一维公式，结果为下图。

均值平滑滤波

以矩阵的方式求矩阵的总和平均值，例如下图一範例3x3的矩阵，将全部数值相加除以阵列的大小等于5，再将计算后的结果放入矩阵的中心位置(1, 1)，这样就完成一个像素的处理，下图二为结果图。

图一

图二

程式码

标头档：Library.h加入
public：均值平滑

/*Blur8bit Parameter:src= source of imagepur= purpose of imagewidth= Image's widthheight= Image's heightsize= blur block*/void Blur8bit(C_UCHAE* src, UCHAE* pur, C_UINT32 width, C_UINT32 height, C_UINT32 size);

实作档：Library.cpp加入

size若是偶数则加一。将图片依公式填补到，输出与输入大小相同。走访每一块block后计算平均赋予目的图片值。

void Library::Blur8bit(C_UCHAE* src, UCHAE* pur, C_UINT32 width, C_UINT32 height, C_UINT32 size){assert(src != nullptr && pur != nullptr);assert(width > 0 && height > 0);if (!(size & 1)){const_cast<UINT32&>(size) = size + 1;}C_UINT32 blockSize = size * size;C_INT32 pad = (size >> 1);C_UINT32 padWidth = width + pad * 2;C_UINT32 padHeight = height + pad * 2;UCHAE* padSrc = new UCHAE[padWidth * padWidth];ImagePadding8bit(src, padSrc, width, height, pad);Image srcImage(padSrc, padWidth, padHeight, MNDT::ImageType::GRAY_8BIT);Image purImage(pur, width, height, MNDT::ImageType::GRAY_8BIT);for (UINT32 row = pad; row < padHeight - pad; row++){for (UINT32 col = pad; col < padWidth - pad; col++){UINT32 sum = 0;for (int32_t blockRow = -pad; blockRow <= pad; blockRow++){for (int32_t blockCol = -pad; blockCol <= pad; blockCol++){sum += srcImage.image[row + blockRow][col + blockCol];}}purImage.image[row - pad][col - pad] = sum / blockSize;}}delete[] padSrc;padSrc = nullptr;}

高斯平滑滤波

而高斯平滑和均值平滑差别在于高斯平滑是依照离中心点位置不同来去分配高斯分布计算平滑，以3x3模板为例带入二为模板如下图一，再来看图二，黄色为依照图一带入的高斯分布机算结果，相对位置乘上以红色为中心的绿色像素模板，计算出来结果再塞入红色的位置当中即可，结果如图三，可以清楚看到使用高斯可以模糊又较不失去焦点。

图一　来源：[4]

图二

图三

程式码

标头档：Library.h加入
public：高斯平滑

/*BlurGauss8bit Parameter:src= source of imagepur= purpose of imagewidth= Image's widthheight= Image's heightsize= gauss temp sizesigma= gauss temp sigma*/void BlurGauss8bit(C_UCHAE* src, UCHAE* pur, C_UINT32 width, C_UINT32 height, C_UINT32 size, C_FLOAT sigma);

private：高斯模板

void Gaussian2DTemp(float** const temp, C_INT32 size, C_FLOAT sigma);

实作档：Library.cpp加入

size若是偶数则加一。将图片依公式填补到，输出与输入大小相同。设定高斯模板。走访每一块block后计算平均赋予目的图片值。

void Library::BlurGauss8bit(C_UCHAE* src, UCHAE* pur, C_UINT32 width, C_UINT32 height, C_UINT32 size, C_FLOAT sigma){assert(src != nullptr && pur != nullptr);assert(width > 0 && height > 0);if (!(size & 1)){const_cast<UINT32&>(size) = size + 1;}C_UINT32 blockSize = size * size;C_INT32 pad = (size >> 1);C_UINT32 padWidth = width + pad * 2;C_UINT32 padHeight = height + pad * 2;UCHAE* padSrc = new UCHAE[padWidth * padWidth];ImagePadding8bit(src, padSrc, width, height, pad);Image srcImage(padSrc, padWidth, padHeight, MNDT::ImageType::GRAY_8BIT);Image purImage(pur, width, height, MNDT::ImageType::GRAY_8BIT);float** temp = new float*[size];for (UINT32 index = 0; index < size; index++){temp[index] = new float[size];}Gaussian2DTemp(temp, size, sigma);for (UINT32 row = pad; row < padHeight - pad; row++){for (UINT32 col = pad; col < padWidth - pad; col++){float sum = 0;for (int32_t blockRow = -pad; blockRow <= pad; blockRow++){for (int32_t blockCol = -pad; blockCol <= pad; blockCol++){sum += (srcImage.image[row + blockRow][col + blockCol] * temp[pad + blockRow][pad + blockCol]);}}purImage.image[row - pad][col - pad] = static_cast<UCHAE>(sum);}}for (UINT32 index = 0; index < size; index++){delete[] temp[index];temp[index] = nullptr;}delete[] temp;temp = nullptr;delete[] padSrc;padSrc = nullptr;}

将二维高斯分布带入指定size大小，sigma起伏设定为可调，mi这边设定为0可忽略。

void Library::Gaussian2DTemp(float** const temp, C_INT32 size, C_FLOAT sigma){float sum = 0;C_INT32 center = size >> 1;C_FLOAT expBase = -(2.0f * sigma * sigma);C_FLOAT scaleBase = (2.0f * sigma * sigma) *  static_cast<float>(MNDT::PI);for (int32_t row = 0; row < size; row++){C_INT32 y = center - row;for (int32_t col = 0; col < size; col++){C_INT32 x = col - center;float gaussNum = 0;gaussNum = exp((x * x + y * y) / expBase);gaussNum = (gaussNum / scaleBase);temp[row][col] = gaussNum;sum += gaussNum;}}for (int32_t row = 0; row < size; row++){for (int32_t col = 0; col < size; col++){temp[row][col] /= sum;}}}

中值滤波

中值滤波也是依照block处理，将block内排序后取中间值后即是中值滤波，下为结果图，使用用5x5的滤镜，可以很明显看到影像中的白点被过滤掉了，但也变得较模糊。

程式码

标头档：Library.h加入
public：中值滤波

/*MedianBlur8bit Parameter:src= source of imagepur= purpose of imagewidth= Image's widthheight= Image's heightsize= blur block*/void MedianBlur8bit(C_UCHAE* src, UCHAE* pur, C_UINT32 width, C_UINT32 height, C_UINT32 size);

实作档：Library.cpp加入

size若是偶数则加一。将图片依公式填补到，输出与输入大小相同。走访每一块block把值塞入指标内后，使用sort排序取出中间值。

void Library::MedianBlur8bit(C_UCHAE* src, UCHAE* pur, C_UINT32 width, C_UINT32 height, C_UINT32 size){assert(src != nullptr && pur != nullptr);assert(width > 0 && height > 0);if (!(size & 1)){const_cast<UINT32&>(size) = size + 1;}C_UINT32 blockSize = size * size;C_UINT32 medianBlockSize = blockSize >> 1;C_INT32 pad = (size >> 1);C_UINT32 padWidth = width + pad * 2;C_UINT32 padHeight = height + pad * 2;UCHAE* padSrc = new UCHAE[padWidth * padWidth];ImagePadding8bit(src, padSrc, width, height, pad);Image srcImage(padSrc, padWidth, padHeight, MNDT::ImageType::GRAY_8BIT);Image purImage(pur, width, height, MNDT::ImageType::GRAY_8BIT);UCHAE* allPix = new UCHAE[blockSize];for (UINT32 row = pad; row < padHeight - pad; row++){for (UINT32 col = pad; col < padWidth - pad; col++){UINT32 index = 0;for (int32_t blockRow = -pad; blockRow <= pad; blockRow++){for (int32_t blockCol = -pad; blockCol <= pad; blockCol++){allPix[index++] = srcImage.image[row + blockRow][col + blockCol];}}std::sort(allPix, allPix + blockSize);purImage.image[row - pad][col - pad] = allPix[medianBlockSize];}}delete[] allPix;allPix = nullptr;delete[] padSrc;padSrc = nullptr;}

双边滤波

[5]提到双边滤波是由一个几何空间滤波器决定系数再由一个像素差来决定滤波器系数，而它和中值滤波差别在于较能保留边缘(中值滤波易模糊)，主要提取高斯分布的距离和Alpha截尾的像素差%来做运算，若想更加了解公式由来请点图解原理，接着来看公式。

高斯距离几何函数

来源：[5]

Alpha截尾灰度差函数

来源：[5]

相乘结果函数

来源：[5]

Block计算平均函数

来源：[5]

可以看到距离和像素的公式是很值观的，现在就让我们来看一下结果图滤波大小为30x30sigma为21。

程式码

标头档：Library.h加入
public：双边滤波

/*MedianBlur8bit Parameter:src= source of imagepur= purpose of imagewidth= Image's widthheight= Image's heightspaceSigma= space sigmacolorSigma= space sigmasize= blur block*/void BilateralBlur8bit(C_UCHAE* src, UCHAE* pur, C_UINT32 width, C_UINT32 height, C_FLOAT spaceSigma, C_FLOAT colorSigma, C_UINT32 size);

private：高斯距离样板、Alpha截尾像素差样板

void BilateralSpaceTemp(float** const temp, C_INT32 size, C_FLOAT sigma);void BilateralColorTemp(float* const temp, C_FLOAT sigma);

实作档：Library.cpp加入

size若是偶数则加一。将图片依公式填补到，输出与输入大小相同。设定高斯距离样板和Alpha截尾像素差样板走访每一块block计算，并用两个变数纪录分子和分母，在带入最后一个公式。

void Library::BilateralBlur8bit(C_UCHAE* src, UCHAE* pur, C_UINT32 width, C_UINT32 height, C_FLOAT spaceSigma, C_FLOAT colorSigma, C_UINT32 size){assert(src != nullptr && pur != nullptr);assert(width > 0 && height > 0);if (!(size & 1)){const_cast<UINT32&>(size) = size + 1;}C_UINT32 blockSize = size * size;C_INT32 pad = (size >> 1);C_UINT32 padWidth = width + pad * 2;C_UINT32 padHeight = height + pad * 2;UCHAE* padSrc = new UCHAE[padWidth * padWidth];ImagePadding8bit(src, padSrc, width, height, pad);Image srcImage(padSrc, padWidth, padHeight, MNDT::ImageType::GRAY_8BIT);Image purImage(pur, width, height, MNDT::ImageType::GRAY_8BIT);C_FLOAT colorBase = -2.0f * colorSigma * colorSigma;C_FLOAT spaceBase = -2.0f * spaceSigma * spaceSigma;float* colorTemp = new float[256];float** spaceTemp = new float *[size];for (UINT32 index = 0; index < size; index++){spaceTemp[index] = new float[size];}BilateralSpaceTemp(spaceTemp, size, spaceSigma);BilateralColorTemp(colorTemp, colorSigma);for (UINT32 row = pad; row < padHeight - pad; row++){for (UINT32 col = pad; col < padWidth - pad; col++){float sum = 0;float base = 0;C_CHAE& nowPix = srcImage.image[row][col];for (int32_t blockRow = -pad; blockRow <= pad; blockRow++){for (int32_t blockCol = -pad; blockCol <= pad; blockCol++){C_CHAE& blockPix = srcImage.image[row + blockRow][col + blockCol];C_UINT32 diff = abs(nowPix - blockPix);C_FLOAT num = colorTemp[diff] * spaceTemp[pad + blockRow][pad + blockCol];base += num;sum += (num * blockPix);}}purImage.image[row - pad][col - pad] = static_cast<UCHAE>(sum / base);}}for (UINT32 index = 0; index < size; index++){delete[] spaceTemp[index];spaceTemp[index] = nullptr;}delete[] spaceTemp;spaceTemp = nullptr;delete[] colorTemp;colorTemp = nullptr;delete[] padSrc;padSrc = nullptr;}

带入高斯距离公式计算size * size大小样板，sigma为可调整，mi为0所以忽略。

void Library::BilateralSpaceTemp(float** const temp, C_INT32 size, C_FLOAT sigma){C_INT32 pad = (size >> 1);C_FLOAT base = -2.0f * sigma * sigma;for (int32_t row = 0; row < size; row++){C_INT32 y = pad - row;for (int32_t col = 0; col < size; col++){C_INT32 x = col - pad;temp[row][col] = ((x * x) + (y * y)) / base;}}}

带入Alpha截尾公式计算256种可能的样板，sigma为可调整。

void Library::BilateralColorTemp(float* const temp, C_FLOAT sigma){C_FLOAT base = -2.0f * sigma * sigma;for (int32_t index = 0; index < 256; index++){temp[index] = exp((index * index) / base);}}

---

结论

虽然实作出来的函数比OpenCV慢，但我想这样的程式码算是很值观的，也可以让自己在未来回来看时不会想那么久，这次最慢的大概是在双边滤波的处理，处理速度与OpenCV差了大概4倍，而去翻OpenCV的过滤器和处理都是用一维去处理，但我想它快的原因主要与parallel_for_这函数有关。
而这次C#部分因为没有特殊处理跟之前文章呼叫方式相同，所以就没有多介绍了，若有问题欢迎提问，有错误的地方麻烦大大纠正谢谢。

参考文献

[1]阿洲(2015). OpenCV教学 | 阿洲的程式教学 from: http://monkeycoding.com/?page_id=12 (2018.10.09).
[2]科男品管圈(2012). 常态分配 (Normal distribution) from:http://lobogaw.pixnet.net/blog/post/90548816-%E5%B8%B8%E6%85%8B%E5%88%86%E9%85%8D-%28normal-distribution%29 (2018.10.10).
[3]维基百科(2018). 常态分布 from:https://zh.wikipedia.org/wiki/%E6%AD%A3%E6%80%81%E5%88%86%E5%B8%83 (2018.10.10).
[4]kancloud(2015). 高斯模糊的算法 · 看云 from https://www.kancloud.cn/kancloud/gaussian_blur/49498 (2018.10.09).
[5]taotao1233(2013). 双边滤波器的原理及实现 - Where there is life, there is hope - CSDN博客 from:https://blog.csdn.net/jinshengtao/article/details/17529617 (2018.10.09).