Geometry & Recognition

이진 영상에서 전경의 경계를 저장하는 함수다. 전경의 최외각을 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;
    }    
};

MFC를 써서 폴리곤의 내부를 칠하기 위해서는 윈도  GDI region을 다루는  CRgn  클래스의 CRgn::CreatePolygonRgn() 메쏘드를 써서 해당 폴리곤으로 bound 된 영역을 만든 이후에 CDC::FillRgn() 함수를 이용하여서 원하는 색으로 칠할 수 있다. 직접  raster에 칠하려고 할 때는 이러한 방법이 너무 불편하다 (물론 라스터(DIB)를 윈도 비트맵으로 변환한 후에 위의 방법대로 칠하고, 다시 윈도 비트맵을 라스터로 바꾸는 작업을 하면 되는데 번거롭다).

쉽게 생각할 수 있는 방법은  point_in_polygon()  방법을 써서 각각의 픽셀이 주어진 폴리곤 내부점인지를 확인하면서 작업을 할 수 있다. 그러나 이 방법은 한 픽셀을 검사하기 위해서 주어진 폴리곤의 모든 에지를 다 건들어야 하는 낭비를 하게 된다. 보다 효율적인 방법은 point_in_polygon 알고리즘을 구현할 때 사용하는 원리를 이용하면 적어도 한 스캔라인의 모든 픽셀에 대해서는 딱 한 번만 전체 폴리곤의 에지를 검사하게 만들 수 있다. inclusion 테스트는 해당점에서 출발하는 수평 반직선이 폴리곤과 몇 번 교차하는지를 세서 내부점인지 외부점인지를 판별한다. 이것을 생각하면 수평선이 폴리곤과 교차하는 점들을 모두 구하여서(한번 폴리곤의 전체 에지 검사가 필요) 크기의 순서대로 정렬할 때, 짝수번째 구간(0->1, 2->3,.....)에 들어가는 스캔라인의 픽셀들은 모두 폴리곤의 내부점이 된다는 것을 알 수 있다.

알고리즘 순서:

1. 각각의 스캔라인에 대해서 교차점들을 구한다:
   꼭짓점 i와 꼭짓점 j(=i-1)를 연결하는 직선의 방정식이 x = x[i] + (x[j] - x[i])*(y - y[i]) / (y[j] - y[i]) 이므로 현재의 scanline 위치에서 y값이 두 꼭지점의 y값 (y[i], y[j]) 사이에 있으면, 교차점의 x-좌표는 윗식의 좌변 값이 된다.
2. 교차점들을 크기의 순서대로 정렬한다.
3. 짝수번째 구간에 들어가는 스캔라인상의 픽셀들은 해당 색깔로 칠한다.

/*sample code;*/
void FillPolygon(const std::vector<CPoint>& poly,
                 DWORD color,
                 BYTE *image, int width, int height, int bpp) {
    std::vector<int> nodeX(poly.size()); //never exceeds npoly;
    for (int y = 0; y < height; y++) {
        // find intersection nodes;
        int nodes = 0;
        for (int i = 0, j = poly.size()-1; i < poly.size(); j = i++)
            if (poly[i].y < y && poly[j].y >= y || poly[j].y < y && poly[i].y >= y) {
                // 수평인 에지는 그리지 않는다(부등식을 자세히 보라)
                nodeX[nodes++] = (int)(poly[i].x + double(y - poly[i].y) * double(poly[j].x \
                                 - poly[i].x) / double(poly[j].y - poly[i].y) + .5); 
                                 // round to integer!!.
            }
            
        // sort nodes (ascending order);
        std::sort(nodeX.begin(), nodeX.begin() + nodes);
        // fill the pixels between node pairs.
        for (int i = 0; i < nodes; i += 2) {
            if (nodeX[i + 0] >= width) break;
            if (nodeX[i + 1] > 0) {
                if (nodeX[i + 0] < 0)     nodeX[i + 0] = 0 ;
                if (nodeX[i + 1] > width) nodeX[i + 1] = width;
                for (int x = nodeX[i]; x < nodeX[i + 1]; x++) 
                    SetPixel(x, y, color, image, width, height, bpp); 
                    //SetPixel()은 다른 프로토타입을 가질 수 있다.
            }
        }
    } 
};

/**
** http://blog.naver.com/helloktk/80050645334 에서 옮김.
*/

주어진 영상에서 선분을 찾고자 하면 어떻게 해야 할까. 우선 에지 영상을 만들어 선분을 강조하고, 선분이 하나만 있는 경우에는 에지 데이터를 가지고 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$ 인 직선 위의 점들이었다는 것을 추정할 수 있다.

std::vector<Line> HoughTransformLine(BYTE* image, int width, int height, /*background=0*/
                                     double rho/*=2*/, 
                                     double theta/*=Pi/180.*/,
                                     int threshold/*=20*/) 
{
    /* use cosine and sine look-up tables*/
    const double irho = 1 / rho;
    const int numangle = (int) (Pi / theta);
    const int numrho = (int) (((width + height) * 2 + 1) / rho);
    const int rwidth = (numrho + 2);
    std::vector<int> accum((numangle + 2) * (numrho + 2), 0);
    double ang; int n;
    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;
    std::vector<Line> vecLine;
    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 );
            }
        }
    }
    return vecLine;
}

아래 그림에서 첫 번째는 소스 이미지이고, 이 이미지에서 $rθ$-평면으로 Hough Transform 한 것이 두 번째 것이다(rho=2, theta=Pi/360). 원본 이미지에는 8개의 선분이 있고, 4 개씩 같은 방향을 갖고 있다. Hough Transform 된 이미지에서 이러한 특징을 볼 수 있다(가로축이 $θ$이고 세로축이 $r$이다. $r$ 축은 가운데가 $r=0$이다). 누적이 많이 된 부분이 두 군데의 동일한 θ에서 나타나고 있다. 그러나 결과에서 8개의 피크를 분리하기는 쉽지 않다. 위의 코드는 각각의 점에서 4방향으로 체크하여 극대 값을 찾고 있으나, 항상 잘 동작하는 것은 아니다.

local-maxima를 좀 더 잘 잡고 싶으면, 위처럼 주변의 4점만 체크할 것이 아니라, 윈도를 좀 더 크게 잡고, 그 윈도 내에서 최댓값을 찾아보는 것도 한 방법이 된다.

/**
** http://blog.naver.com/helloktk/80051779331
*/