원문링크
Gaussian Filter
영상을 부드럽게 하기 위해서는 신호처리분야에서 저역 통과 필터(low pass filter)를 통과시켜야 된다.
저역 통과 필터는 주파수 공간에서 영상의 저주파 성분만을 통과시키는 필터로써 mean filter, weighted mean filter
등 여러가지가 있지만 여기에서는 Gaussian Filter에 대해 알아본다.
가우시안 함수는 자연현상을 가장 잘 표현하는 함수 식 중의 하나이다.
1차원 가우시안 분포를 식으로 나타내면 아래와 같다.
μ는 일반적으로 0으로 놓는다. σ는 값의 분포를 결정지으며, 표준편차를 의미한다.
일반적으로 μ = 0과 σ^2 = 1를 정규분포라고 부른다.
영상처리에서 가우시안 필터는 2차원 가우시안 함수 값을 이용하여 마스크를 생성하고,
입력영상과 마스크 연산을 수행해주는 것을 의미한다.
이 때, σ에 의해 필터의 크기가 결정되는데, 이 때 값의 분포와 그래프는 아래와 같다.
-4 * σ <= x <= 4 * σ
지금까지는 1차원 가우시안 함수를 설명했지만, 실제 영상에서는 2차원 가우시안 함수가 사용되어야 한다.
2차원 가우시안 함수는 아래와 같다.
위의 식에서 μ = 0 이라고 가정했다. (이후부터 μ는 고려하지 않는다.)
아래는 σ에 따른 이미지의 Bluring을 보여주고 있다.
최적화
σ= 1.0으로 놓을 경우 G(x, y)는 총 9 * 9의 filter를 픽셀별로 연산해야 한다.
G(x, y)를 보면 알겠지만 연산량이 상당하다는 것을 알 수 있다.
따라서 가우시안 함수를 최적화 할 필요가 있는데 2차원의 가우시안 함수는 1차원 가우시안 함수로 분리하여
표현할 수 있으며 아래와 같다.
가우시안 필터를 계산하였지만 여전히 픽셀당 9 * 9 를 계산하는것은 부담으로 작용한다.
하지만 2차원 가우시안 함수가 1차원 가우시안 함수로 표현할 수 있다는 것은 마스크 계산량을 줄일 수 있는
여지를 제공한다.
아래 그림은 p4를 기준으로 가우시안 필터를 적용한 상황을 보여주고 있다.
2차원을 1차원으로 풀게되면 3 * 3일 경우 2개의 가우시안 함수값만 계산하면 된다.
그런데 그림을 보면 또다른 계산의 중복이 있음을 알 수 있다.
아래 그림에서 P(-2)를 보도록 하자
그림을 보면 P(-2)와 P4는 각각 G(x-1)P0, G(x)P1, G(x+1)P2를 공통으로 사용하고 있음을
알 수 있다.
따라서 중복된 연산을 제거하여 최적화를 할 수 있다.
여기에서는 일단 모든 이미지의 픽셀에 대해서 자신을 기준으로 x축방향으로 필터링 연산을 한다.
P4일 경우 아래와 같이 연산을 하면 된다.
모든 픽셀에 대해서 x축 기준 가우시안 마스크을 적용하고 나면 그 결과를 가지고
y축 가우시안 마스크를 적용하면된다.
y 축 적용은 이전 픽셀의 x축 적용결과와 다음 픽셀의 y축 적용결과를 합산하기만 하면 되며,
수식을 풀어보면 2차 가우스함수 필터를 적용한 것과 완전히 동일하다.
이런식의 최적화는 열우선 메모리접근에 강점을 보이는 OpenCL에서 더욱더 장점으로 작용한다.
가우시안 필터를 구현한 소스는 아래와 같다.
<![CDATA[
void CVisitorGaussian::inMsg_img(CImg* img)
{
int width, height;
int i;
m_img = img;
for(i=0;i< 3;++i)
{
if(m_tmpData[i])
delete[] m_tmpData[i];
}
width = img->outMsg_img()->outMsg_width();
height = img->outMsg_img()->outMsg_height();
m_tmpData[0] = new float[width * height];
m_tmpData[1] = new float[width * height];
m_tmpData[2] = new float[width * height];
}
void CVisitorGaussian::inMsg_sigma(float sigma, int maskSize)
{
int i, cnt, inx;
float tmpFront, tmpBack, tmp;
m_filterCnt = maskSize;
tmpFront = 1.0f / (sqrt(2.0f * 3.1415926535f) * sigma);
tmp = 2.0f * sigma * sigma;
cnt = maskSize >> 1;
inx = 0;
for(i=-cnt;i< 0;++i)
{
tmpBack = -((i * i) / tmp);
m_filter[inx] = tmpFront * exp(tmpBack);
m_filter[maskSize -1 - inx] = m_filter[inx];
++inx;
}
tmpBack = tmpFront * exp(0.0f);
m_filter[inx] = tmpBack;
}
void CVisitorGaussian::visit(CImg* inst)
{
CScoBmp *dstBmp, *srcBmp;
int width, height, tmpX, tmpY, tmp, i, j;
int maskCnt, a, orientDelta;
int rowInx, colInx;
float sumR, sumG, sumB;
unsigned char uR, uG, uB;
dstBmp = m_img->outMsg_img();
srcBmp = inst->outMsg_img();
width = dstBmp->outMsg_width();
height = dstBmp->outMsg_height();
maskCnt = m_filterCnt;
orientDelta = maskCnt >> 1;
for(j=0;j< height;++j)
{
colInx = width * j;
for(i=0;i< width;++i)
{
sumR = 0.0f;
sumG = 0.0f;
sumB = 0.0f;
tmpX = i - orientDelta;
for(a=0;a< maskCnt;++a)
{
tmp = (tmpX < 0) ? 0 : tmpX;
tmp = (tmp >= width) ? (width - 1) : tmp;
srcBmp->outMsg_pixel(tmp, j, uR, uG, uB);
sumR += (uR * m_filter[a]);
sumG += (uG * m_filter[a]);
sumB += (uB * m_filter[a]);
++tmpX;
}
rowInx = colInx + i;
(m_tmpData[0])[rowInx] = sumR;
(m_tmpData[1])[rowInx] = sumG;
(m_tmpData[2])[rowInx] = sumB;
}
}
for(j=0;j< height;++j)
{
for(i=0;i< width;++i)
{
sumR = 0.0f;
sumG = 0.0f;
sumB = 0.0f;
tmpY = j - orientDelta;
for(a=0;a< maskCnt;++a)
{
tmp = (tmpY < 0) ? 0 : tmpY;
tmp = (tmp >= height) ? (height - 1) : tmp;
rowInx = width * tmp + i;
sumR += ((m_tmpData[0])[rowInx] * m_filter[a]);
sumG += ((m_tmpData[1])[rowInx] * m_filter[a]);
sumB += ((m_tmpData[2])[rowInx] * m_filter[a]);;
++tmpY;
}
sumR += 0.5f;
sumG += 0.5f;
sumB += 0.5f;
Limit(sumR, 0, 255);
Limit(sumG, 0, 255);
Limit(sumB, 0, 255);
dstBmp->inMsg_pixel(i, j, sumR, sumG, sumB);
}
}
}
]]>
