이미지에서 Voronoi diagram으로 영역을 분할할 때 각 픽셀이 어느 Voronoi cell에 포함되는가를 알아야 하는 경우가 있다. 보통은 Voronoi 다이어그램으로 구한 cell을 폴리곤으로 표현하고, 해당 픽셀이 어느 폴리곤에 들어가는 가는 체크 해야 한다. 그러나, 이 과정은 복잡하고 계산이 많이 발생한다. 이미지에 만들어진 Voronoi diagram의 경우 cell mask를 이용하면 해당 픽셀이 어느 cell에 들어있는지를 바로 판단할 수 있다. 특히, cell의 개수가 적은 경우 mask를 gray 이미지로 처리할 수 있어서 메모리 사용도 줄일 수 있다. Voronoi diagram의 이미지화 과정은 Voronoi 알고리즘을 이용할 필요는 없고 단지 각 cell을 형성하는 픽셀들은 그 cell의 중심까지 거리가 다른 cell보다 가깝다는 사실만 이용한다.
void rasterize_voronoi(std::vector<CPoint>& vorocenter,
BYTE *image, int width, int height) {
std::vector<BYTE> red(vorocenter.size()),
green(vorocenter.size()),
blue(vorocenter.size());
for (int i = vorocenter.size(); i-->0;) {
red[i] = rand() % 256;
green[i] = rand() % 256;
blue[i] = rand() % 256;
}
for (int y = 0; y < height; y++) {
for (int x = 0; x < width; x++) {
int min_id = 0;
int dist2_min = INT_MAX;
for (int k = vorocenter.size(); k-->0;) {
int dx = x - vorocenter[k].x;
int dy = y - vorocenter[k].y;
int dist2 = dx * dx + dy * dy;
if (dist2 < dist2_min) {
dist2_min = dist2;
min_id = k;
}
}
*image++ = blue[min_id];
*image++ = green[min_id];
*image++ = red[min_id];
}
}
// draw cell center;
}
어떤 현상을 설명하는 이론 모델을 만들려면 현상에서 관측 데이터를 모아야 한다. 그런데 관측 데이터에는 모델에 대한 잘못된 가정이나 측정의 오차 등에서 생기는 여러 가지 형태의 노이즈가 들어 있게 마련이다. 이러한 노이즈는 모델을 정확히 예측하는 것을 방해한다. 모델 파라미터의 예측을 방해하는 데이터(outliners)가 들어 있는 관측 데이터로부터 어떻게 하면 모델을 구축할 수 있는가라는 질문에 대한 한 해결책을 RANSAC("RANdomSAmple Consensus") 알고리즘이 제공한다. 이 알고리즘은 1981년 Fischler과 Bolles에 의해서 제안되었다.
RANSAC 알고리즘에 들어가는 입력은 관측된 데이터, 관측 결과를 피팅하거나 설명할 수 있는 모델 파라미터, 그리고 신뢰성을 나타내는 파라미터를 필요로 한다.
RANSAC 알고리즘은 주어진 원본 데이터에서 일부를 임의로 선택하는 과정을 반복하여서 최적의 파라미터를 예측하는 프러시저의 형태를 가진다. 전체의 관측 데이터가 M개 있고, 모델 파라미터를 예측하는데 N개의 데이터가 필요한 경우, 알고리즘의 동작은 1. 임의로 관측 데이터에서 N개의 부분 데이터를 선택한다. 2. 선택 데이터를 (가상의) inlier로 생각하고 모델을 예측한다. 3. 원본 데이터(M) 중에서 예측된 모델에 잘 맞는가를 체크한다. 잘 맞는 데이터 수를 K라고 한다. 4. 만약에 K가 충분하면, 이 모델을 확정하고 알고리즘을 종료한다. 5. 모델이 좋지 않으면 1→4 과정을 L번 반복한다. 6. 충분한 반복 후에도 좋은 모델을 찾지 못하면 예측에 실패한 것으로 간주한다. 1. 임계값 K는 주어진 전체 데이터 중에서 어느 정도의 비율로 찾는 모델에 잘 맞는가로 판단하는 기준인데, 사용자가 결정해야 한다. 대략적으로 주어진 샘플에서 inlier의 비율을 Pg정도라고 생각되면 다음 정도로 잡으면 된다:
K = M * (1- Pg) 2. 얼마나 많은 반복을 해야하는가? 주어진 관측 데이터에서 inlier일 확률이 Pg인 경우에 L번의 모델 예측 시도가 실패할 확률을 계산하여서 이것이 주어진 설정값, pfail 보다도 작은 경우에 모델 예측의 실패로 처리하므로, pfail = L번의 모델예측 실패 확률 = (한 번 시도가 실패할 확률)L = (1 - 한 번 시도가 성공할 확률))L = (1 - (랜덤 데이터가 모델을 예측할 확률)N))L = (1- (Pg) N)L 이 사실로부터 최대 반복 횟수는 L = log(pfail)/log(1-(Pg)N) 로 주어진다. inlier가 반이 섞인 경우 Pg (=주어진 데이터중에서 inlier일 확률) = 0.5, pfail = 0.01 인 경우, N = 3 (세 점만 있으면 모델 구성이 가능한 원의 피팅이 한 예)이면 최대 반복 횟수는 윗 식에 적용하면, L = 35 회 RANSAC 알고리즘은 주어진 관측 데이터에 많은 outlier가 존재하더라도 매우 정확히 모델 파라미터를 예측할 수 있는 능력이 있는 robust 한 알고리즘이지만, 얼마나 많은 반복을 해야 하는가에 대한 상한 값이 제공되지 않는다(사용자가 pfail 를 설정한다). 따라서 사용자 설정에 의한 상한 값은 최적의 모델이 아닐 수 있다.
"The RANSAC procedure is opposite to that of conventional smoothing techniques: Rather than using as much of the data as possible to obtain an initial solution and then attempting to eliminate the invalid data points, RANSAC uses as small an initial data set as feasible and enlarges this set with consistent data when possible".
Martin A. Fischler and Robert C. Bolles (June 1981). "Random Sample Consensus: A Paradigm for Model Fitting with Applications to Image Analysis and Automated Cartography". Comm. of the ACM24: 381–395.
이진 영상에서 전경의 경계를 저장하는 함수다. 전경의 최외각을 8-방향 연결성을 체크하면서 추적하도록 설계되었다(chain code를 참조하면 된다). 전경의 두께가 1 픽셀이더라도 항상 닫힌 윤곽선을 형성하게 만들었다(다르게 행동하도록 바꿀 수 있다). 전경에 서로 연결되지 않은 blob가 여러 개 있는 경우도 쉽게 처리할 수 있게 출력은 벡터 컨테이너를 사용했다.
struct Cntr {
int x, y ; //position;
int idirn ;//chain-code;
Cntr() {}
Cntr(int x_, int y_, int idirn_) : x(x_), y(y_), idirn(idirn_){}
};
typedef std::vector<Cntr> CntrVector;
#define FGVAL 255
#define BGVAL 0
#define VISITED 33 /* pixel value on accepted contour */
int GetNextCntr (BYTE** image, int *x, int *y, int *idirn);
int ContourTrace (BYTE** image/*binary image(FGVAL,BGVAL)*/, int width, int height,
std::vector<CntrVector* > &CntList) {
/* CAUTION:: one-pixel border should have BGVAL!!!!*/
for (int x = width; x-->0;) image[0][x] = image[height-1][x] = BGVAL;
for (int y = height; y-->0;) image[y][0] = image[y][width-1] = BGVAL;
for (int y = height-1; y-->1;) {
for (int x = width-1; x-->1;) {
if (image[y][x] == FGVAL && image[y][x - 1] == BGVAL) {
CntrVector *pCntrVec = new CntrVector ;
int idirn = 2;
int xStart = x, yStart = y;
pCntrVec->push_back(Cntr(x, y, idirn));
do {
image[y][x] = VISITED; /* set the default value to VISITED */
GetNextCntr (image, &x, &y, &idirn);
pCntrVec->push_back(Cntr(x, y, idirn));
} while (!(x == xStart && y == yStart));
CntList.push_back(pCntrVec);
}
}
}
return CntList.size();
};
/* 다음 픽셀을 검사. 체인 코드; */ int GetNextCntr (BYTE** image, int *x, int *y, int *idirn);
void GetNextCntr (BYTE** image, int *x, int *y, int *idirn) {
/*idirn::direction from last pixel on contour */
int i, j;
BYTE ring[16];
/* find neighborhood ring of pixels */
ring[0] = ring[ 8] = image[*y - 1][*x ];
ring[1] = ring[ 9] = image[*y - 1][*x + 1];
ring[2] = ring[10] = image[*y ][*x + 1];
ring[3] = ring[11] = image[*y + 1][*x + 1];
ring[4] = ring[12] = image[*y + 1][*x ];
ring[5] = ring[13] = image[*y + 1][*x - 1];
ring[6] = ring[14] = image[*y ][*x - 1];
ring[7] = ring[15] = image[*y - 1][*x - 1];
i = (*idirn + 4 + 1) % 8;
if (i == 0) i = 8; /* this never happens */
for (j = 0; j < 8; j++, i++)
if (ring[i] != BGVAL && ring[i - 1] == BGVAL) break;
if (j == 8) return; /* isolated pixel */
*idirn = i % 8;
switch (i) {
case 1:
case 9:
*x = *x + 1; *y = *y - 1;
break;
case 2:
case 10:
*x = *x + 1;
break;
case 3:
case 11:
*x = *x + 1; *y = *y + 1;
break;
case 4:
case 12:
*y = *y + 1;
break;
case 5:
case 13:
*x = *x - 1; *y = *y + 1;
break;
case 6:
case 14:
*x = *x - 1;
break;
case 7:
case 15:
*x = *x - 1; *y = *y - 1;
break;
case 8:
case 0:
*y = *y - 1;
break;
default:
break;
}
};
주어진 영상에서 선분을 찾고자 하면 어떻게 해야 할까. 우선 에지 영상을 만들어 선분을 강조하고, 선분이 하나만 있는 경우에는 에지 데이터를 가지고 line-fitting 알고리즘을 적용해서 해당 선분을 기술하는 파라미터를 추출할 수 있다. 그러나 영상이 선분을 여러 개 포함하는 경우에는 이 방법을 이용하기 어렵다(불가능한 것은 아니지만 노이즈 등에 의한 영향을 고려해야 하므로 복잡한 과정을 거쳐야 한다). 평면에서 선분은 원점까지 거리와 그것의 수직(수평)인 방향만 주어지면 결정이 된다. 직선에서 원점까지 거리를 $r$, 법선이 $x$-축과 이루는 각도가 $θ$면
$$r = \cos(θ) x + \sin(θ) y;$$
의 형태로 주어진다. 즉, $(r,θ)$ 한쌍이 주어지면 직선 한 개가 정의된다 (주어진 $(r,θ)$ 값이 만드는 직선을 $a x + b y = c$ 꼴로 쓰면 계수는 $a = \cos θ$, $b = \sin θ$, $c = r$로 표현된다).
이 직선 방정식은 $rθ$-평면의 stripe=$[0,∞) \times [0,2π]$을 $xy$-평면으로 보내는 변환으로도 생각할 수 있다. 이 변환은 $rθ$-평면의 한 점을 $xy$-평면의 한 직선으로 보낸다. 역으로는 $xy$-평면의 한 점은 $rθ$-평면의 한 곡선으로 변환이 된다. 직선 상의 점들은 같은 원점 거리=$r$, 같은 기울기=$\theta$를 가지므로, 직선 위의 각 점들이 $rθ$-평면에 만드는 곡선들은 공통의 교점을 가지게 된다. 이 교점의 위치가 $xy$-평면에서 직선을 기술하는 파라미터를 준다.
Hough Transform은 이 변환의 특성을 이용한 것이다. 이미지에서 에지에 해당하는 점들을 위의 변환에 의해서 $rθ$-평면의 곡선으로 보내는 것이다 (프로그램적으로는, 에지점 좌표 $(x, y)$가 주어지면 $θ=0 \rightarrow π$까지 일정 간격으로 증가시키면서 곡선 $r=x \cos(θ) + y\sin(θ)$의 값을 구해서 메모리 상의 $[r,θ]$ 지점의 도수를 1씩 증가시킨다). 만약 에지에 해당되는 점들이 직선의 관계를 가지게 되면, 각 점들에 해당하는 $rθ$-평면에서 곡선들은 한 교점을 형성하게 될 것이다. 긴 선분은 $rθ$-평면의 한 점에서 더 많은 곡선들의 교점을 형성하게 된다. 따라서, $rθ$-평면에서 이러한 교점의 히스토그램을 구하여, 교점 수가 일정 이상 누적이 된 경우를 취하면, $xy$-이미지 평면에서 일정한 길이 이상을 갖는 선분을 골라낼 수 있다.
그러면 왜 $rθ$-평면으로 변환을 생각하는 것일까? 직선이 $y= a x+b$의 형태로 표시되므로 기울기와 $y$축과의 교점을 기술하는 $ab$-공간을 이용할 수도 있지만, 이 경우에 $y$-축에 평행인 직선의 경우에 기울기 $a$ 값이 무한히 커지는 경우가 발생하여 저장공간 할당에 문제가 생긴다. 실제적인 문제에서 되도록이면 compact 한 공간으로 변환을 고려해야 하는 것이다. $rθ$-공간으로의 변환은 유한한 이미지의 경우 항상 유한한 $rθ$-공간의 영역으로 변환된다.
아래의 그림은 $x-y=-5$인 직선상의 세 점 $(5,10)$, $(10,15)$, $(15,20)$에 해당하는 $rθ$-평면상의 세 곡선을 보인 것이다: $r = 5 \cos(θ) + 10\sin(θ)$, $r = 10 \cos(θ) + 15 \sin(θ)$, $r = 15 \cos(θ) + 20 \sin(θ)$. 이 세 곡선은 $(r,θ) = (5/\sqrt{2}, 3 \pi /4)$ 지점에서 만남을 알 수 있다. 이 만나는 지점을 구하면, 거꾸로 세 점이 $x-y=-5$ 인 직선 위의 점들이었다는 것을 추정할 수 있다.
BOOL HoughTransformLine(BYTE* image, int width, int height, /*background=0*/
double rho/*=2*/,
double theta/*=Pi/180.*/,
int threshold/*=20*/,
std::vector<Line>& vecLine) {
/* use cosine and sine look-up tables*/
double irho = 1 / rho;
int numangle = (int) (Pi / theta);
int numrho = (int) (((width + height) * 2 + 1) / rho);
int *accum = (int*)malloc( sizeof(accum[0]) * (numangle + 2) * (numrho + 2) );
memset( accum, 0, sizeof(accum[0]) * (numangle + 2) * (numrho + 2) );
double ang;
int n, rwidth = (numrho + 2);
for (int y = 0; y < height; y++ ){
for (int x = 0; x < width; x++ ){
if ( image[y * width + x] != 0 ){ // only for foreground pixels;
for (ang = 0, n = 0; n < numangle; ang += theta, n++ ) {
double rho = (x * cos(ang) + y * sin(ang)) * irho ;
int r = rho >= 0 ? int(rho + .5) : (-int(-rho + .5));//round to nearest int;
// accum을 이차원배열로 생각하고, 1 픽셀의 border를 두는 형태로 한다.
// 이는 아래의 local-maxima를 찾을때 경계를 고려할 필요가 없어서 유리하다.
r += (numrho - 1) / 2;
accum[(n + 1) * rwidth + r + 1]++;
}
}
}
}
//find local maxima;
for (int r = 0; r < numrho; r++ ) {
for (int n = 0; n < numangle; n++ ) {
int base = (n + 1) * rwidth + r + 1;
//test whether it is local-maximum(4-way);
if ( accum[base] > threshold &&
accum[base] > accum[base - 1] && accum[base] > accum[base + 1] &&
accum[base] > accum[base - rwidth] && accum[base] > accum[base + rwidth] )
{
Line line;
line.rho = (r - (numrho - 1) * .5) * rho;
line.angle = n * theta;
vecLine.push_back( line );
}
}
}
free(accum);
return TRUE;
}
아래 그림에서 첫 번째는 소스 이미지이고, 이 이미지에서 $rθ$-평면으로 Hough Transform 한 것이 두 번째 것이다(rho=2, theta=Pi/360). 원본 이미지에는 8개의 선분이 있고, 4 개씩 같은 방향을 갖고 있다. Hough Transform 된 이미지에서 이러한 특징을 볼 수 있다(가로축이 $θ$이고 세로축이 $r$이다. $r$ 축은 가운데가 $r=0$이다). 누적이 많이 된 부분이 두 군데의 동일한 θ에서 나타나고 있다. 그러나 결과에서 8개의 피크를 분리하기는 쉽지 않다. 위의 코드는 각각의 점에서 4방향으로 체크하여 극대 값을 찾고 있으나, 항상 잘 동작하는 것은 아니다.
local-maxima를 좀 더 잘 잡고 싶으면, 위처럼 주변의 4점만 체크할 것이 아니라, 윈도를 좀 더 크게 잡고, 그 윈도 내에서 최댓값을 찾아보는 것도 한 방법이 된다.
주어진 이미지의 gaussian scale space의 표현은 이미지와 gaussian kernel의 convolution 결과로 주어진다. $$L(x, y; σ)= G(x, y; σ) * I(x, y);$$
이 표현은 임의의 $σ$에 대해서 정의되나 실제적으로는 유한한 구간의 값들만 취급한다. 그리고 스케일 값이 한 옥타브를 넘어서는 경우에는 원본 이미지의 크기를 유지할 필요가 없다. 주어진 스케일보다 작은 특징들은 고려 대상이 되지 않는다는 점을 인식하면, 스케일이 두 배로 커지는 경우에는 원본 이미지를 반으로 줄이고 스케일은 동일하게 유지하는 것이 훨씬 더 이득이다. 따라서 scale space를 이미지 pyramid로 표현하는 것이 자연스럽다. 아래 예는 3 옥타브까지 scale space를 표현한 것이다. 각 옥타브에서 3 개의 구간(수직 방향)을 갖는다(각 옥타브에서 마지막 이미지는 비교를 위해 넣은 것이다). 한 단계 위 옥타브의 이미지는 전 단계 옥타브의 동일한 레벨의 이미지를 1/2 다운 샘플링하면 바로 얻을 수 있다. *** $L(x; σ) = G(x; σ) * I(x)$ 의 1/2 down-sampling 이미지 $L'(x;\sigma)$는 (1-dim만 고려해도 충분하다) $$\begin{align} L'(x; σ)\equiv L(2x; σ) &= \int \frac{1}{σ} \exp\Big[-\frac{(2x-x')^2}{2σ^2}\Big] I(x') dx', \quad (x'\rightarrow 2x') \\&=\int \frac{1}{σ} \exp\Big[- \frac{(2x-2x')^2}{2σ^2} \Big] I(2x') d(2x')\\&= \int \frac{1}{\sigma/2} \exp \Big[ -\frac{(x-x')^2 }{2(σ/2)^2}\Big] I(2x') dx' \\ &=G(x; σ/2) * I(2 x);\end{align}$$ 임을 알 수 있다. 즉, 스케일 $σ$로 gaussian convolution 이미지를 1/2 다운샘플링한 것은 1/2 다운샘플링한 이미지를 $σ/2$ 스케일로 gaussian convolution 결과와 같다. 따라서 스케일이 2배인 이미지가 필요하면 현재의 스케일에서 이미지를 반으로 줄이면 된다.
