快速高斯模糊(IIR遞迴高斯模糊)
阿新 • • 發佈:2018-12-30
本人是影象處理的初學者,在http://www.cnblogs.com/Imageshop/p/3145216.html博文中看到有關影象柔光的特效處理原理後自己也動手實現了一下,這裡包括一個快速高斯模糊的演算法,具體原理可以參考論文《Recursive Implementation of the gaussian filter》,這個演算法的主要優點是高斯模糊的速度與高斯函式標準差沒有關係,本文主要是給出上述演算法的C++程式實現,這其中參考了utopiaT的http://www.cnblogs.com/utopiaT/p/3305732.html和http://www.cnblogs.com/utopiaT/p/3308991.html
BOOL WINAPI GaussSmoothIIR24(LPSTR lpDIBBits, LONG lWidth, LONG lHeight, float fSigma)
{
int size = 0;
// 影象每行的位元組數
LONG lLineBytes;
// 計算影象每行的位元組數
lLineBytes = WIDTHBYTES(lWidth * 24);
// 指向源影象的指標
unsigned char* lpSrc;
float *w1 = new float[(lWidth + 3 + 3) * 3];//為了消除往右方移動,這裡給w陣列增加三個RGB畫素即9個畫素
float *w2 = new float[(lWidth + 3 + 3) * 3];
float *wB1 = new float[lHeight + 3 + 3];//為了消除往下方移動,這裡給w陣列增加三個畫素
float *wB2 = new float[lHeight + 3 + 3];
float *wG1 = new float[lHeight + 3 + 3];
float *wG2 = new float [lHeight + 3 + 3];
float *wR1 = new float[lHeight + 3 + 3];
float *wR2 = new float[lHeight + 3 + 3];
float *in = new float[lHeight*lWidth * 3];
float *out = new float[lHeight*lWidth * 3];
//----------------計算高斯核---------------------------------------//
float q = 0;
float q2, q3;
double b0, b1, b2, b3;
//double b[4] = { 0 };//這裡聲明瞭一個4個元素的陣列,只賦值一個0表示第一個元素賦值為0,其他三個元素自動為0.
double B = 0;
if (fSigma >= 2.5)
{
q = 0.98711 * fSigma - 0.96330;
}
else if ((fSigma >= 0.5) && (fSigma < 2.5))
{
q = 3.97156 - 4.14554 * (float)sqrt((double)1 - 0.26891 * fSigma);
}
else
{
q = 0.1147705018520355224609375;
}
q2 = q * q;
q3 = q * q2;
b0 = (1.57825 + (2.44413*q) + (1.4281 *q2) + (0.422205*q3));
b1 = ((2.44413*q) + (2.85619*q2) + (1.26661 *q3));
b2 = (-((1.4281*q2) + (1.26661 *q3)));
b3 = ((0.422205*q3));
B = 1.0 - ((b1 + b2 + b3) / b0);
//加速方法 減少迴圈多次/b0
b1 /= b0;
b2 /= b0;
b3 /= b0;
//----------------計算高斯核結束---------------------------------------//
for (int i = 0; i < lHeight; i++)
{
for (int j = 0; j < lWidth * 3; j++)
{
lpSrc = (unsigned char*)lpDIBBits + lLineBytes*(lHeight - 1 - i) + j;
/* 0-255 => 1-256 */
in[i*lWidth * 3 + j] = (float)(*lpSrc + 1.0);
}
}
for (int k = 0; k < lHeight; k++)
{
size = lWidth * 3;
size -= 3;
w1[0 * 3] = in[k * lWidth * 3 + 0];
w1[1 * 3] = in[k * lWidth * 3 + 0];
w1[2 * 3] = in[k * lWidth * 3 + 0];
w1[0 * 3 + 1] = in[k * lWidth * 3 + 1];
w1[1 * 3 + 1] = in[k * lWidth * 3 + 1];
w1[2 * 3 + 1] = in[k * lWidth * 3 + 1];
w1[0 * 3 + 2] = in[k * lWidth * 3 + 2];
w1[1 * 3 + 2] = in[k * lWidth * 3 + 2];
w1[2 * 3 + 2] = in[k * lWidth * 3 + 2];
for (int i = 0, n = 3 * 3; i <= size; i += 3, n += 3)
{
w1[n] = (float)(B*in[k * lWidth * 3 + i] +
((b1 * w1[n - 1 * 3] +
b2 * w1[n - 2 * 3] +
b3 * w1[n - 3 * 3])));
w1[n + 1] = (float)(B*in[k * lWidth * 3 + i + 1] +
((b1 * w1[n - 1 * 3 + 1] +
b2 * w1[n - 2 * 3 + 1] +
b3 * w1[n - 3 * 3 + 1])));
w1[n + 2] = (float)(B*in[k * lWidth * 3 + i + 2] +
((b1 * w1[n - 1 * 3 + 2] +
b2 * w1[n - 2 * 3 + 2] +
b3 * w1[n - 3 * 3 + 2])));
}
//w2[size + 1 * 3] = w1[size + 3 * 3];
//w2[size + 2 * 3] = w1[size + 3 * 3];
//w2[size + 3 * 3] = w1[size + 3 * 3];
//w2[size + 1 * 3 + 1] = w1[size + 3 * 3 + 1];
//w2[size + 2 * 3 + 1] = w1[size + 3 * 3 + 1];
//w2[size + 3 * 3 + 1] = w1[size + 3 * 3 + 1];
//w2[size + 1 * 3 + 2] = w1[size + 3 * 3 + 2];
//w2[size + 2 * 3 + 2] = w1[size + 3 * 3 + 2];
//w2[size + 3 * 3 + 2] = w1[size + 3 * 3 + 2];
w2[size + 9 + 1 * 3] = w1[size + 3 * 3];//為了消除往下方移動,這裡把賦值資料都加了一個9
w2[size + 9 + 2 * 3] = w1[size + 3 * 3];
w2[size + 9 + 3 * 3] = w1[size + 3 * 3];
w2[size + 9 + 1 * 3 + 1] = w1[size + 3 * 3 + 1];
w2[size + 9 + 2 * 3 + 1] = w1[size + 3 * 3 + 1];
w2[size + 9 + 3 * 3 + 1] = w1[size + 3 * 3 + 1];
w2[size + 9 + 1 * 3 + 2] = w1[size + 3 * 3 + 2];
w2[size + 9 + 2 * 3 + 2] = w1[size + 3 * 3 + 2];
w2[size + 9 + 3 * 3 + 2] = w1[size + 3 * 3 + 2];
for (int i = size, n = i + 9; i >= 0; i -= 3, n -= 3)//為了消除往下方移動,這裡n為i+9不是i
{
w2[n] = out[k * lWidth * 3 + i] = (float)(B*w1[n] +
((b1 * w2[n + 1 * 3 + 0] +
b2 * w2[n + 2 * 3 + 0] +
b3 * w2[n + 3 * 3 + 0])));
w2[n + 1] = out[k * lWidth * 3 + i + 1] = (float)(B*w1[n + 1] +
((b1 * w2[n + 1 * 3 + 1] +
b2 * w2[n + 2 * 3 + 1] +
b3 * w2[n + 3 * 3 + 1])));
w2[n + 2] = out[k * lWidth * 3 + i + 2] = (float)(B*w1[n + 2] +
((b1 * w2[n + 1 * 3 + 2] +
b2 * w2[n + 2 * 3 + 2] +
b3 * w2[n + 3 * 3 + 2])));
}
}
//下面的橫向處理直接將資料 out 與 in 對調 省去memcpy
//memcpy(in, out, channelsize * sizeof(float));
//memset(out, 0 , channelsize * sizeof(float));
for (int k = 0; k < lWidth * 3; k += 3)
{
size = lHeight;
size -= 1;
wB1[0] = out[k + 0];
wB1[1] = out[k + 0];
wB1[2] = out[k + 0];
wG1[0] = out[k + 1];
wG1[1] = out[k + 1];
wG1[2] = out[k + 1];
wR1[0] = out[k + 2];
wR1[1] = out[k + 2];
wR1[2] = out[k + 2];
for (int i = 0, n = 3; i <= size; i++, n++)
{
wB1[n] = (float)(B*out[i*lWidth * 3 + k] +
((b1 * wB1[n - 1] +
b2 * wB1[n - 2] +
b3 * wB1[n - 3])));
wG1[n] = (float)(B*out[i*lWidth * 3 + k + 1] +
((b1 * wG1[n - 1] +
b2 * wG1[n - 2] +
b3 * wG1[n - 3])));
wR1[n] = (float)(B*out[i*lWidth * 3 + k + 2] +
((b1 * wR1[n - 1] +
b2 * wR1[n - 2] +
b3 * wR1[n - 3])));
}
//wB2[size + 1] = wB1[size + 3];
//wB2[size + 2] = wB1[size + 3];
//wB2[size + 3] = wB1[size + 3];
//wG2[size + 1] = wG1[size + 3];
//wG2[size + 2] = wG1[size + 3];
//wG2[size + 3] = wG1[size + 3];
//wR2[size + 1] = wR1[size + 3];
//wR2[size + 2] = wR1[size + 3];
//wR2[size + 3] = wR1[size + 3];
wB2[size + 4] = wB1[size + 3];//為了消除往下方移動,這裡本該賦值給size + 1和2和3更改為size + 4和5和6
wB2[size + 5] = wB1[size + 3];
wB2[size + 6] = wB1[size + 3];
wG2[size + 4] = wG1[size + 3];
wG2[size + 5] = wG1[size + 3];
wG2[size + 6] = wG1[size + 3];
wR2[size + 4] = wR1[size + 3];
wR2[size + 5] = wR1[size + 3];
wR2[size + 6] = wR1[size + 3];
for (int i = size, n = i+3; i >= 0; i--, n--)//為了消除往下方移動,這裡n為i+3不是i
{
wB2[n] = in[i * lWidth * 3 + k] = (float)(B*wB1[n] +
((b1 * wB2[n + 1] +
b2 * wB2[n + 2] +
b3 * wB2[n + 3])));
wG2[n] = in[i * lWidth * 3 + k + 1] = (float)(B*wG1[n] +
((b1 * wG2[n + 1] +
b2 * wG2[n + 2] +
b3 * wG2[n + 3])));
wR2[n] = in[i * lWidth * 3 + k + 2] = (float)(B*wR1[n] +
((b1 * wR2[n + 1] +
b2 * wR2[n + 2] +
b3 * wR2[n + 3])));
}
}
//拷貝結果到函式輸出指標
for (int i = 0; i < lHeight; i++)
{
for (int j = 0; j < lWidth * 3; j++)
{
lpSrc = (unsigned char*)lpDIBBits + lLineBytes*(lHeight - 1 - i) + j;
/* 1-256 => 0-255 */
*lpSrc = in[i*lWidth * 3 + j] - 1;
}
}
delete[] w1;
delete[] w2;
delete[] wB1;
delete[] wB2;
delete[] wG1;
delete[] wG2;
delete[] wR1;
delete[] wR2;
delete[] in;
delete[] out;
return TRUE;
}
下面是處理效果:
原圖:
IIR高斯模糊(fSigma=10,用時53ms)的處理效果:
普通二維卷積高斯模糊(fSigma=10,用時13120ms)的處理效果: