影像處理會依照[1]的章節進行，但不直接使用OpenCV而是實作原理，主要了解影像處理原理的特性，雖然速度不比OpenCV好，但在看OpenCV原始碼你會更好理解運作原理並且可以加以修改自己所需要模式，這樣做可讓自己更有機會發展出更多不同的處理方式，這次會介紹4種平滑濾波器平滑濾波、高斯濾波、中值濾波和雙邊濾波，首先來介紹一些基本概念。
相信大家對常態分佈不陌生如下圖。[3]提到在自然界或數學當中，一個不明的隨機變量我們會使用常態分佈來去表示或計算，而在人工智慧當中高斯分布也是其中的種要角色之一。主要將它當做計算一個機率的分布，接著介紹它的數學公式。
來源：[2]
一維公式：
來源[3]
這裡sigma在控制高低如下圖一，mi在控制偏移量如下圖二。但通常影像處理的維度是二維，因此要帶入二維常態分佈。
圖一
圖二
二維公式：
來源[4]
二維公式推導為兩個機率分布相乘，一維公式 * 一維公式，結果為下圖。
以矩陣的方式求矩陣的總和平均值，例如下圖一範例3x3的矩陣，將全部數值相加除以陣列的大小等於5，再將計算後的結果放入矩陣的中心位置(1, 1)，這樣就完成一個像素的處理，下圖二為結果圖。
圖一
圖二
標頭檔：Library.h加入
public：均值平滑
Kevin
/*
Blur8bit Parameter:
src = source of image
pur = purpose of image
width = Image's width
height = Image's height
size = 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 image
pur = purpose of image
width = Image's width
height = Image's height
size = gauss temp size
sigma = 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 image
pur = purpose of image
width = Image's width
height = Image's height
size = 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 image
pur = purpose of image
width = Image's width
height = Image's height
spaceSigma = space sigma
colorSigma = space sigma
size = 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若是偶數則加一。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;
}
}
}
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).
iThome鐵人賽
看影片追技術