Geometry & Recognition

주어진 N개의 관측값을 k개의 가우시안 mixture로 모델링한 후 EM 알고리즘을 적용하여 모델을 결정하는 파라미터를 추정하고자 한다. 추정해야 할 모델 파라미터는 k개의 가우시안을 기술하는 평균값(μ)과 공분산 행렬(Covariance Matrix: Σ)과 이들 가우시안 분포의 mixing(α) 정도를 나타내는 파라미터이다.

다음은 이 과정을 수행하는 C++ class와 사용 예제, 및 결과를 보여준다.

//가우시안 kernel Class;
class GaussKernel2D {
    double mx, my;
    double sdx, sdy, sdxy;
    double weight;
public:
    std::vector<double> posterior;
    std::vector<double> wgauss;

    void init(POINT* pts, int np, int sx, int sy, int w) {
        posterior.resize(np);
        wgauss.resize(np);
        weight = w ;
        // initialize model parameters(random하게 선택함-->try another method!!::
        // 주어진 데이터로 전체 분포의 범위와 중심을 알 수 있고, 이를 주어진 클래스 수만큼 임의로 
        // 분배하는 방식으로 초기조건을 설정하면 보다 안정적으로 동작한다)
        mx = sx * double(rand()) / RAND_MAX;
        my = sy * double(rand()) / RAND_MAX;
        sdx = sx / 4 + 1;
        sdy = sy / 4 + 1;
        mx += rand() % 100 ;
        my += rand() % 100 ;
        sdxy = 0;
    }
    double gauss2d(double x, double y) { 
        double varx = sdx * sdx;
        double vary = sdy * sdy;
        double det = varx * vary - sdxy * sdxy;
        double idet = 1.0 / det ;
        double dxx = (x - mx) * (x - mx);
        double dyy = (y - my) * (y - my);
        double dxy = (x - mx) * (y - my);
        return (1.0 / sqrt(det) / 6.28319) * exp(-0.5 * (dxx * vary \
                               + dyy * varx - 2. * dxy * dxy) * idet); 
    }
    void getParams(POINT *pts, int np, double prior = 0) {
        double sx = 0, sy = 0, sxx = 0, syy = 0, sxy = 0;
        for (int j = 0; j < np; j++) {
            double x = pts[j].x, y = pts[j].y;
            sx  += posterior[j] * x;
            sy  += posterior[j] * y;
            sxx += posterior[j] * x * x;
            syy += posterior[j] * y * y;
            sxy += posterior[j] * x * y;
        }
        double denorm = np * weight;
        mx = sx/denorm; my= sy / denorm;
        double devx = sxx / denorm - mx * mx ;
        if (devx <= 0) devx = 0.001;
        sdx = sqrt(devx);
        double devy=syy / denorm - my * my;
        if (devy <= 0) devy = 0.001;
        sdy = sqrt(devy);
        sdxy = sxy / denorm - mx * my;
        // if prior = non-zero -> weight = weight*(1-alpha)+alpha*prior; alpha=0.1?
    };
    // weight; // posterior;
    void estimate(double *px, int np) {
        weight = 0;
        for(int j = 0; j < np; j++) {
            posterior[j] = wgauss[j] / px[j];    
            weight += posterior[j];
        }
        weight /= np; 
    } 
    //P(x|thetal) * prior;
    void setProb(POINT *pts, int np) {
        for(int i = 0; i < np; i++){
            wgauss[i] = weight * gauss2d(pts[i].x, pts[i].y);
        }
    }
    void Draw(CDC* pDC, DWORD color = RGB(0xFF, 0, 0)) {
        CPen pen0(PS_SOLID, 1, color);
        CPen* pOld = pDC->SelectObject(&pen0);
        drawCon(pDC, mx, my, sdx, sdy, sdxy);        // draw ellipses;
        pDC->Ellipse(mx - 2, my - 2, mx + 2, my + 2);// draw center of ellipse;
        pDC->SelectObject(pOld);            
    }
};

void em_main(POINT *pts, int np) {
    GaussKernel2D kernel[NKERNELS];
    int nclass = NKERNELS;
    double weights[20] = {1};
    std::vector<double> px(np);
    double wsum = 0 ;

    for (int i = 0; i < nclass; i++) {
        kernel[i].init(pts, np, 400, 400, weights[i]);
        wsum += weights[i];
    };    
#define MAX_ITER 50
    for (int iter = 0; iter < MAX_ITER; ++iter){   
        for (i = 0; i < np; i++) px[i] = 0;
        for (int k = 0; k < nclass; k++){
            GaussKernel2D & gker = kernel[k];
            gker.setProb(pts, np); 
            for (int i = 0; i < np; i++){
                px[i] += gker.wgauss[i] ;
            }
        }        
        for (k = 0; k < nclass; k++) {
            kernel[k].estimate(&px[0], np);
            kernel[k].getParams(pts, np);
        }
        //또는 log-likelihood를 계산하여서 그 변화가 적으면 loop-끝내면 된다..
    }
}

//참고 : 아래의 데이터는 사전에 라벨링이 된 것이 아니다. 컬러링은 한번 계산한 후에 분포에 맞게 컬러를 조절하여서 다시 계산한 것이다.

 

f(y|θ);

원을 그리기 위해서는 원의 방정식을 이용하면 된다. 반지름이 $r$인 경우에 원점이 중심인 원은 $x^2 + y^2 = r^2$이므로 (중심이 원점이 아닌 경우는 단순 평행 이동만 하면 되기 때문에 따로 취급하지 않는다), $x$를 독립변수로 잡으면 $y= \sqrt {r^2 - y^2}$으로 주어져서 1/2 사분면을 그릴 수 있다. 그리고 부호를 반대로 하면 3/4 분면에 그려진다. 그러나 $x$가 $r$에 근접하는 경우 $y$의 변화율이 급격해져서 $y$값이 연속적으로 그려지지 않는다. 이 문제는 전체 4 사분면을 45도 단위로 8개 영역으로 쪼개서 $x$와 $y$의 역할을 적당히 바꾸면서 그리면 이러한 문제를 피할 수 있다. 그러나, 여전히 실수 연산을 수행해야 하는 문제가 남아 있다. 이를 피하는 방법을 고안하도록 하자. 힌트는 Bresenham의 라인 알고리즘에서 얻을 수 있다. 우선, 영역을 $x = 0$에서 $x = y$인 90-도에서 45-도 사이의 구간만 생각하자. 이 구간에서 $y$의 변화율의 크기는 항상 0에서 1 사이 이므로, 8-방향으로 연결성을 갖도록 그리려면 현재의 점이 $(x, y)$로 주어지는 경우

          (x, y) ---> (x+1, y)     변화율 = 0;
                      (x+1, y-1)   변화율 = 1      (절대값)

의 두 가지 경우만 있다.  그럼 이 두 점 중 어느 것을 선택을 해야 하는가? 이것을 판별하기 위한 식은 이 두 점의 중심인 $(x, y - 1/2)$이 원의 내부에 있는가 또는 외부에 있는가를 보고서 결정하면 된다.

          (x+1, y-1/2) 이  원의 내부점인 경우---->(x+1, y)  선택
                           원의 외부점인 경우---->(x+1, y-1) 선택;

하면 각각의 상황에서 보다 정확한 원에 가까운 점을 선택한 것이 된다. 다음 점에서도 동일한 과정을 거치면 된다.

그러면 $(x+1, y-1/2)$가 원의 내부점인지 외부점인지를 항상 계산해야 하는가? 매 스텝마다 이것을 계산을 할 필요는 없다. 왜냐면 처음 $(0, r-1/2)$에서 계산을 하고

          d = x^2+y^2-r^2 at (0, r-1/2) ;
             = (5-4*r)/4  ; -> 이런 식으로 표현을 하는 것은 정수 연산만으로 
                               하기 위해서 계산 우선 순위를 조절한 것임.

이것과, $F(x, y)=x^2-y^2-r^2$ 라고 할 때,

(1),  $(x+1, y)$를 선택할 때 새로운 $d$값은 $(x+2, y-1/2)$에서 계산하므로

         d_old = F(x+1, y-1/2) = (x+1)^2+(y-1/2)^2-r^2 ;
         d_new = F(x+2, y-1/2) = d_old + 2*(x+1) + 1;

(2),  $(x+1, y-1)$를 선택할 때, 새로운 $d$값은 $(x+2, y-1-1/2)$에서 계산하므로

        d_old = F(x+1, y-1/2) = (x+1)^2+(y-1/2)^2-r^2 ;
        d_new = F(x+2, y-3/2) = d_old + 2 * ((x+1) -(y-1)) + 1;

인 관계를 이용하면 쉽게 정수 연산만으로 $d$값을 업데이트할 수 있다. 따라서 전체적인 코드는

   int x = 0, y = r ;
   int d = (5 - 4 * r) / 4;
   setPixel (x, y) ; 
   while ( x < y ) {
        ++x ;
        if (d < 0)                            // (x+1, y) 선택;
            d +=  2*x + 1;                    // ++x와 연관해서 잘 살펴야 한다.
        else {                                // (x+1, y-1) 선택;
            --y ;
            d += 2 * (x - y) + 1;             // ++x,--y;
        }
        setPixel(x, y);
   };

8 방향 대칭성으로 인해 점 $(x, y)$를 그리면 아래의 점도 그려야 한다 (처음 점 $(0, r)$은 4방향만, 그리고 $x=y$인 경우도 마찬가지나 따로 분리하기보다는 그냥 그리는 것이 더 낫다).

   {(x,-y), (-x,y), (-x,-y), (y,x), (y,-x), (-y,x),(-y,-x)}

단, 원의 중심이 원점이 아닐 경우는 전체적으로 평행이동을 시키면 된다.

이미지 상의 한 점은 픽셀로 표현이 된다. 그러나 실제적으로 색을 표시하기 위해서 픽셀은 유한한 크기를 갖는다. 따라서 픽셀이 한 점을 표시한다는 말을 할 때는 실제상으로는 그 점은 픽셀의 중앙을 의미한다고 생각하면 된다. 예를 들어, 점(1, 3)에 색을 칠한다고 하면 엄밀한 의미로는 [0.5,1.5] x [2.5,3.5]로 주어지는 정사각형 영역을 채우는 것을 의미한다. 픽셀을 가지고 작업을 하는 경우는 항상 정수단위로 움직인다. 이런 사실 때문에 원치 않은 결과가 나오기도 한다. 기울기가 매우 작은 직선을 그리는 작업을 생각해 보자. 예를 들면 기울기가 2/10인 경우에 픽셀상에서의 작업은 0부터 4까지는 같은 y=0인 픽셀에, 5, 9까지는 y=1인 픽셀에 칠해야 한다. 실수 값을 다룰 수 없으므로 4와 5를 지나면서 갑자기 계단을 형성하게 된다. 

이러한 계단 현상은 눈에 보이는 그래픽의 품질을 저하시키는 결과를 초래하므로 이것을 피하는 선 그리기 방법이 있어야 한다. 간단히 생각할 수 있는 방법은 블러링을 주는 방법이다. 우선 기울기가 1보다 작은 경우만 생각하자. 나머지의 경우는 대칭을 생각하면 쉽게 확장이 가능하다. 이 상황에서는 x 좌표값이 독립변수의 역할을 하고, y가 종속변수다. 각각의 픽셀(x=정수)에서 y값을 계산하면 일반적으로 실수 값을 갖는다. 이 y값에 가까운 두 개의 정수 y값에 해당하는 픽셀에 색을 칠하는데, 이때 농도는 거리의 차이에 반하도록 잡는다.

     (x, y)  -> (x, int( y ) ),       명암 가중치 = 1 – (y-int(y));

               -> (x, int( y ) + 1 ),  명암 가중치 = y – int(y) ;

처음점과 끝점 좀 더 특별한 처리가 필요하다. 여기서도 역시 블러링의 효과를 주기 위해서 가중치를 주는데 y에 의한 가중치에 덧붙여 x에 의한 가중치도 준다. x가 정수 값을 취할 경우에 50%의 가중치를 주는 것을 기준으로 하면, x + 0.5를 기준으로 이것의 정수 값이 벗어나는 정도를 취하면 된다.  그리고  x의 경우는 가장 가까운 정수에서부터 시작한다.

     (x, y)     ->  (int(x + 0. 5), int(y)),  명암 가중치 = (1 – (y - int(y))) * (1 - (x + 0.5 - int(x + 0.5)))

                  ->  (int(x + 0. 5), int(y) + 1), 명암 가중치 = (y - int(y)) * (1- ( x + 0.5 - int(x + 0.5)))

오른쪽 끝점은 x에서 오는 가중치는 왼쪽의 complement이므로 x+0.5-int(x+0.5)로 바꾸어야 한다.
아래의 그림은 이 알고리즘으로 그린 선분(위)과 Bresenham 알고리즘(아래)을 이용한 선분을 비교한 것이다( 2배로 확대시킨 그림이다)

이 알고리즘은 1991년도에 Wu에 의해서 제시되었다.

#define FPART(x) ((x) - int(x))     // fractional part of (x)
void DrawLine(double x1, double y1, double x2, double y2) {
    if (x2 < x1) {
    	swap(x1, x2); swap(y1, y2);
    }
    double dx = x2 - x1;
    double dy = y2 - y1;
    double grad = dy / dx ;
    if ( fabs(dx) > fabs(dy)) { // 수평에 가까운 선분;
        //(start point)
        int xs = int(x1 + 0.5);                 //x1 에 제일 가까운 정수 xs ;
        double ys = y1 + grad * (xs - x1);      //픽셀점(xs)에서 y값;
        int ixs = xs ;
        int iys = int(ys) ;
        double xgap = 1 - FPART(x1 + .5);   //x1에 의한 가중치
        setPixel(ixs, iys,     (1 - FPART(ys)) * xgap);
        setPixel(ixs, iys + 1, FPART(ys) * xgap);
        //(end point)
        int xe = int(x2 + 0.5);                 //x2에 제일 가까운 정수 xe;
        double ye = y2 + grad * (xe - x2);      //픽셀점(xe)에서의 y값
        xgap = FPART(x2 + .5);
        int ixe = xe; 
        int iye = int(y2);
        setPixel(ixe, iye,    (1 - FPART(ye)) * xgap);
        setPixel(ixe, iye + 1, FPART(ye) * xgap);
        
        //[ixs+1,..., ixe-1];
        double y = ys + grad ;            // 시작(xs) 다음점에서 y-값;
        for (int ix = ixs + 1; ix < ixe; ix++) {
            setPixel(ix, int(y),      1 - FPART(y));
            setPixel(ix, int(y) + 1,  FPART(y));
            y += grad ;
        }
    } else {
        // 이 경우는 x와 y의 역할을 바꾼다;생략;
    };
};

이 함수는 음수가 아닌 인자에 대해서만 성립한다. 그러나 일반적인 경우로도 쉽게 확장이 가능하다.
참고 :
1. http://en.wikipedia.org/wiki/Xiaolin_Wu's_line_algorithm
2. http://freespace.virgin.net/hugo.elias/graphics/x_wuline.htm
3. http://www.gamedev.net/reference/articles/article382.asp