Geometry & Recognition

이미지 처리에서 각도를 다루게 되면 삼각함수를 이용해야 한다. 일반적인 PC 상에서는 CPU의 성능이 크게 향상이 되어 있고, floating point 연산이 잘 지원이 되므로 큰 문제가 없이 삼각함수를 바로 이용할 수가 있다. 그러나 embeded 환경에서는 연산 속도가 그다지 빠르지 않고 floating point 연산은 에뮬레이션을 하는 경우가 대부분이다. 이런 조건에서 루프 내에서 연속적으로 삼각함수를 호출하는 것은 성능에 많은 영향을 주게 되므로 다른 방법을 고려해야 한다. 일반적으로 각도를 점진적으로 증가시키면 cosine이나 sine 값을 계산해야 할 필요가 있는 경우에는 아래와 같이 처리한다:

for (int i = 0; i < nsteps; ++i) {
    double angle = i * angle_increment;
    double cos_val = cos( angle ) ;
    double sin_val = sin( abgle ) ;
    // do something useful with cos_val, and sin_val;
    // .......................;
}

물론 미리 계산된 cosine/sine의 값을 가지고 처리할 수도 있다. 하지만 매우 작은 각도 단위로 계산하는 경우는 lookup table의 크기가 매우 커지게 되므로 직접 계산을 하는 것보다 이득이 없을 수 있다. 이런 경우에는 다음의  삼각함수의 관계식을 이용하면 매번 cosine이나 sine 함수의 호출 없이도 필요한 값을 얻을 수 있다;

$$\cos(\theta+\delta)-\cos(\theta)=- (2\sin^2(\delta/2) \cos(\theta) + \sin(\delta)\sin(\theta))= -(\alpha \cos(\theta) + \beta \sin (\theta));$$

$$\sin(\theta+\delta)-\sin(\theta)= -(2\sin^2(\delta/2) \sin(\theta) - \sin(\delta)\cos(\theta))= -(\alpha \sin(\theta) - \beta \cos (\theta));$$

여기서 $\alpha = 2 \sin^2(\delta/2)$, $\beta = \sin(\delta)$다. 처음에 한번 계산된 $\alpha$와 $\beta$ 값을 이용하면, 추가적인 계산이 없이 점진적으로 $\cos$과 $\sin$의 값을 계산할 수 있다.

    int i = 0 ;
    double cos_val = 1 ; // = cos(0);
    double sin_val = 0 ; // = sin(0);
    double beta = sin(angle_increment) ;
    double alpha = sin(angle_increment/2) ;
    alpha = 2 * alpha * alpha ; 
    do {
        // do something useful with cos_val and sin_val ;
        // ...........................;
        // calculate next values;

        double cos_inc = (alpha * cos_val + beta * sin_val) ;
        double sin_inc = (alpha * sin_val - beta * cos_val) ;
        cos_val -= cos_inc ;
        sin_val -= sin_inc ;
    } while (++i < nsteps);

이 방법은 이론상 정확한 방법이나 round-off error가 전파가 되어서 반복 회수가 커질수록 에러가 누적이 될 위험이 있다. 이를 방지하기 위해 중간에 카운터를 두어서 일정한 반복 이후에는 cosine/sine값 계산을 다시 하여 에러를 리셋시키면 된다.

 이미지 처리 연산 과정에서 제곱근 값을 계산해야 경우가 많이 생긴다. 이 경우 math-라이브러리의 sqrt() 함수를 사용하는 것보다 원하는 정밀도까지만 계산을 하는 근사적인 제곱근 함수를 만들어서 사용하는 것이 계산 비용 측면에서 유리할 수 있다. 

주어진 양수 $x$의 제곱근 $y = \sqrt{x}$를 구하는 것은 $F(y) = y^2 - x$의 근을 구하는 것과 같다. 프로그램적으로 근은 반복적인 방법에 의해서 구해지는데 많이 사용하는 알고리즘 중 하나가 Newton 방법이다. Newton 방법은 한 단계에서 근의 후보가 정해지면 그 값에서 F에 접하는 직선을 구한 후 그 직선이 가로축과 만나는 점을 다시 새로운 근으로 잡는 반복 과정을 일정한 에러 범위 내에 들어올 때까지 수행한다.

함수 $F(y)$의 근은 아래의 점화식에 의해서 구해지게 된다:

$$\tt y(n + 1) = y(n) - F(y(n)) / F'(y(n)), n = 0,1, 2,...;$$

제곱근의 경우에는 ($F(y) = y^2 - x$) 점화식은 

$$\tt y(n + 1) = (y(n) + x / y(n)) / 2;$$

의 형태로 주어진다. 그러나, 이 점화식에는 계산 부하가 상대적으로 큰 나누기 연산이 들어 있다. 그런데 $\sqrt{x} = x \times (1 / \sqrt{x})$로 표현할 수 있으므로 제곱근의 역수를 쉽게 계산하는 알고리즘을 찾으면 된다. $x$의 제곱근의 역수를 구하는 점화식은 (이 경우 $F(y) = 1 / y^2 - x$)

$$\tt y(n + 1) = y(n) * (3 - y(n) * y(n) * x) / 2;$$

이므로 전체적인 제곱근을 구하는 pseudo 코드는

float fast_sqrt(float x) {
    xhalf = 0.5F * x;
    y = initial_guess_of (1/sqrt(x));
    while (error is not small)
        // find approx inverse sqrt root of (x);
        y = y * 1.5 - y * y * y * xhalf;
    end while

    // finally return sqrt(x) = x * (1 / sqrt(x));
    return (x * y);
}

보통 iteration 횟수는 고정이 되어 있다. 남은 문제는 초기값 추정뿐이다. Newton 방법에서 잘 추정된 초기값은 몇 번의 반복으로도 원하는 정밀도를 얻을 수 있게 하지만 잘못된 추정은 수렴 속도를 더디게 한다. 또한 초기값의 추정 과정도 매우 간단해야 하는데, 좋은 추정을 하기 위해서는 IEEE에서 규정하는 double/float 포맷에 대한 정보와 수치해석 지식이 필요하다.

IEEE의 32비트 float 상수의 비트 구성은 $\tt 0 \sim 22$번째 비트까지는 소수점 아래의 유효숫자(Mantissa: $\tt M$)를 표기하고, $\tt 23$번째부터 $\tt 30$번째까지는 지수 부분(Exponent: $\tt E$)을 표시하는데 $\tt 127$만큼 더해진 값이다(+-지수를 표현하기 위한 것으로 실제 지수는 $\tt e = E - 127$). 마지막 $\tt 31$비트는 부호를 결정한다. 임의의 양의 float 상수 $\tt x$가 

$$\tt x = (1 + M) 2^{e} = (1 + M) 2^{E - 127};$$

로 쓰여지면, 정수로 변환된 비트 데이터는 $\tt ix = *(int~*)\&x = (E + M) 2^{23}$의 값을 갖는다. $\tt x$ 제곱근의 역수는

$$\tt \frac{1}{\sqrt{x}} = \frac{1}{\sqrt{1 + M}} 2^{- e / 2} ;$$

로 표현되므로 비트 데이터에서 지수부는 $\tt (- e / 2) + 127 = 190 - E / 2$로 표현된다. 

이제 필요한 것은 잘 선택된 $\tt 1/\sqrt{x}$의 초기값을 찾는 것으로, 첫 번째로 시도할 근사는 유효숫자 부분 $(\tt 1/\sqrt{1 + M})$을 제거하고, 단순히 지수 부분($\tt 2^{- e / 2}$)만 살리면 된다. 정수 표현으로 바꾸면 $\tt 190$이 $\tt 23$에서 $\tt 30$번째 비트를 채우도록 $\tt 23$을 시프트 하고, $\tt - E / 2$ 부분은 원래 float 상수의 정수 표현을 나누면 유효숫자 부분에 해당하는 정도만 에러가 발생한다:

   정수 표현 = (190 - E / 2) << 23;

                  = (190 << 23) - (x의 정수형 변환 & (0xFF << 23)) >> 1;

                  ~ (190 << 23) - (x의 정수형 변환) >> 1;    // 유효 숫자 / 2 만큼 에러 발생;

                  = 0x5F000000 - (*(int *)&x) >> 1;

유효숫자 부분은 좀 더 높은 초기 정밀도를 갖도록 근사를 할 수 있는데 결과만 쓰면 $\tt 0x5F000000$ 대신에 $\tt 0x5F3759DF$를 이용하면 된다.

//빠른 float 제곱근의 역수;
float fast_invsqrt(float x) {
    float xhalf = 0.5F * x;
    // initial guess of 1 / sqrt(x);
    int i = 0x5F3759DF - (*(int *)&x) >> 1;
    x = *(float *)&i;    
    //iterate
    x = x * (1.5F - x * x * xhalf);
    x = x * (1.5F - x * x * xhalf);
    // end of calculation of 1 / sqrt(x);
    return (x);
};

// 빠른 float 제곱근;
float fast_sqrt(float x) {
    float xhalf = 0.5F * x;
    //initial guess of 1 / sqrt(x);
    int i = 0x5F3759DF - (*(int *)&x) >> 1;
    float y = *(float *)&i;    
    //iterate
    y = y * (1.5F - y * y * xhalf);
    y = y * (1.5F - y * y * xhalf);
    //end of calculation of 1/sqrt(x);
    //return sqrt(x)=x * 1/sqrt(x);
    return (x * y)   
}; 

$ \log_2 (x) = e + \log_2 (1 + M)$이므로 x의 근삿값은 $\log_2 (1+M)$의 근삿값을 잘 선택하면 된다. 함수 $f(M)=\log_2 (1+M) - M$는 $M=1/\ln(2)-1$에서 최댓값 $\tt 0.0860713$을 가지므로

$$0\le \log_2 (1+M)-M \le 0.0860713$$

이다. 따라서 오차를 최댓값의 절반으로 잡아서 근사하면

$$\log_2 (1+M) \approx M+ \delta,\quad \delta = \frac{0.0860713}{2}=0.0430385$$ 

그리고 $\delta$는 $x$에 무관한 값이 된다. $x$의 비트 데이터를 정수로 표현하면

$$ i_x= 2^{23} (E+M) = 2^{23}(e + \delta +M) + 2^{23}(127-\delta) \approx 2^{23}\log_2(x) + 2^{23}(127-\delta).$$ 로 근사할 수 있다. 두 번째 항은 $x$에 무관한 값이다. 따라서, 제곱근의 역수에 해당하는 비트 데이터의 근사적인 정수 표현은

$$i_{\frac{1}{\sqrt{x}}} \approx 2^{23}\log_2 (1/\sqrt{x}) + 2^{23}(127-\delta)\approx \frac{3}{2} 2^{23} (127-\delta ) -\frac{1}{2} i_x$$

로 쓸 수 있고, 위에서 구한 $\delta$을 대입하면

$$\tt  i_{\frac{1}{\sqrt{x}}} \approx  0x5F37BCB6 - i_x>>1$$ 

을 얻을 수 있다. $\delta = 0.0450466$을 대입하면 위에 적힌 코드에서 사용한 $\tt 0x5F3759DF$가 됨을 알 수 있다. 그러나 어떤 근거로 이 값을 선택했는지는...

$$g(x, y; \sigma)= \frac{1}{2\pi \sigma^2}\exp\left( - \frac{x^2 +y^2 }{2\sigma^2 }\right)$$

원점을 중심으로 분산이 1 인 정규분포를 갖는 평면상의 점을 발생시킬 필요가 있는 경우에 아래의 함수를 이용한다: Box-Muller method.

내장 함수 rand()를 이용하면 $[0,1]$ 구간에서 균일한 분포를 만들 수 있다. $[0,1]\times[0,1]$ 영역에서 균일한 분포를 갖는 두 개의 1차원 랜덤 변수 $x_1$, $x_2$을 이용해서 확률 변수 $y_1$, $y_2$의 2차원 Gaussian 분포를 주는 변환을 구하려면, 우선 확률 보존에 의해서 

$$ dx_1 dx_2 = \frac{1}{\sqrt{2\pi}} e^{-y_1^2 /2 } dy_1 \frac{1}{\sqrt{2\pi}} e^{-y_2^2 /2 } dy_2 $$

을 만족시켜야 한다. 극좌표 $\rho^2 = y_1^2 + y_2^2$, $\phi = \tan^{-1}(y_2/y_1)$을 도입하면 

$$ dx_1 dx_2 = e^{-\rho^2/2} \rho d\rho \frac{d\phi}{2\pi}$$로 쓸 수 있으므로 미소 변환은

$$ dx_1 = e^{-\rho^2/2} \frac{d\rho^2}{2},\quad dx_2 = \frac{d\phi}{2\pi}$$ 로 선택할 수 있다. 따라서

$$ x_1 = \exp[-\rho^2/2] = \exp[- (y_1^2 + y_2^2 )/2],$$

$$ x_2 = \frac{1}{2\pi} \tan^{-1} (y_2/y_1).$$

그리고, 역변환은

$$ y_1 = \sqrt{-2 \ln (x_1) } \cos (2\pi x_2),\quad y_2 = \sqrt{-2 \ln(x_1) } \sin (2\pi x_2).$$

삼각함수의 계산을 피하기 위해서 반지름이 1인 원 내부에서 정의되는 두 개의 균일한 랜덤 변수 $v_1, v_2$를 도입하면, 다음 변환에 의해서 반지름 1인 원 영역이 변의 길이가 1인 정사각형 영역으로 보내진다.

$$x_1 = R^2(=S) = v_1^2 + v_2^2, \quad 2\pi x_2 = \tan^{-1} (v_2/v_1).$$

따라서 반지름이 1인 원 내부의 랜덤 변수를 써서 2차원의 Gaussian 분포를 만들어주는 변환은

$$ y_1 = \sqrt{-2\ln S} \frac{v_1}{\sqrt{S}}, \quad y_2 = \sqrt{-2\ln S} \frac{v_2}{\sqrt{S}}$$

// 참고;  Numerical Receipe;
#define drand() ((double)rand() / RAND_MAX)   // [0,1]
// mean = 0; sigma = 1;
void normal_2D(double *x, double *y) {
    while (1) {
        double V1 = -1. + (2. * drand());
        double V2 = -1. + (2. * drand());
        double S = V1 * V1 + V2 * V2;
        if (S < 1.) {
            double rad = sqrt(-2.*log(S) / S);
            *x = V1 * rad;
            *y = V2 * rad;
            /* if (mean, sigma) were given,
            *x = meanx + sigma * (*x);
            *y = meany + sigma * (*y);
            */
            return;
        }
    }
};

평면 위에 점집합이 주어지고 이들을 잘 기술하는 직선의 방정식을 구해야 할 경우가 많이 발생한다. 이미지의 에지 정보를 이용해 선분을 찾는 경우에 hough transform과 같은 알고리즘을 이용하는 할 수도 있지만 수치해석적으로 직접 fitting을 할 수도 있다. 점집합의 데이터를 취합하는 과정은 항상 노이즈에 노출이 되므로 직선 위의 점뿐만 아니라 직선에서 (많이) 벗어난 outlier들이 많이 들어온다. 따라서 line-fitting은 이러한 outlier에 대해서 매우 robust 해야 한다. 데이터 fitting의 경우에 초기에 대략적인 fitting에서 초기 파라미터를 세팅하고, 이것을 이용하여서 점차로 정밀하게 세팅을 해나가는 반복적인 방법을 많이 이용한다. 입력 데이터가 $\{(x_i, y_i)| i=0,..., N-1\}$로 주어지는 경우에 많이 이용하는 최소자승법에서는 각 $x_i$에서 직선상의 $y$ 값과 주어진 $y_i$의 차이(residual)의 제곱을 최소로 하는 직선의 기울기와 $y$ 절편을 찾는다. 그러나 데이터가 $y$축에 평행하게 분포하는 경우를 다루지 못하게 되며, 데이터 점에서 직선까지 거리를 비교하는 것이 아니라 $y$값의 차이만 비교하므로 outlier의 영향을 매우 심하게 받는다.

이러한 문제를 제거 또는 완화하기 위해서는 PCA(principal axis analysis)를 이용할 수 있다. 점들이 선분을 구성하는 경우, 선분 방향으로는 점 위치의 편차가 크지만 수직 방향으로는 편차가 상대적으로 작다. 따라서 평면에서 점 분포에 대한 공분산 행렬 $\tt Cov$의 고윳값과 고유 벡터를 구하면, 큰 고윳값을 갖는 고유 벡터 방향이 선분의 방향이 될 것이다.

$$ {\tt Cov}[\{ (x_i, y_i)\}]=\frac{1}{N} \begin{pmatrix} \sum _i(x_i- \bar{x})^2 & \sum_i (x_i-\bar{x})( y_i-\bar{y}) \\ \sum_ i (x_i-\bar{x})( y_i-\bar{y})  & \sum_i (y_i- \bar{y})^2   \end{pmatrix}$$

잘 피팅이 이루어지려면 두 고윳값의 차이가 커야 한다. 또한 outlier에 robust 한 피팅이 되기 위해서는 각 점에 가중치를 부여해서 공분산 행렬에 기여하는 가중치를 다르게 하는 알고리즘을 구성해야 한다. 처음 방향을 설정할 때는 모든 점에 동일한 가중치를 부여하여 선분의 방향을 구한 후 다음번 계산에서는 직선에서 먼 점이 공분산 행렬에 기여하는 weight를 줄여 주는 식으로 하면 된다. weight는 점과 직선과의 거리에 의존하나 그 형태는 항상 정해진 것이 아니다.

// 점에서 직선까지 거리;
double DistanceToLine(CPoint P, double line[4]) {
    // 중심에서 P까지 변위;
	double dx = P.x - line[2], dy = P.y - line[3]; 
    // 직선의 법선으로 정사영 길이 = 직선까지 거리;
    return fabs(-line[1] * dx + line[0] * dy);
}
// PCA-방법에 의한 line-fitting;
double LineFit_PCA(std::vector<CPoint>& P, std::vector<double>& weight, double line[4]) {
    int res = 1;
    // 초기화 시 weight[i] = 1.;
    double sx = 0, sy = 0, sxx = 0, syy = 0, sxy = 0, sw = 0;
    for (int i = P.size(); i-->0;) {
         int x = P[i].x, y = P[i].y;
         double w = weight[i]; 
         sx += w * x; sy += w * y;
         sxx += w * x * x; syy += w * y * y;
         sxy += w * x * y; 
         sw  += w; 
    }
    // variances;
    double vxx = (sxx - sx * sx / sw) / sw;
    double vxy = (sxy - sx * sy / sw) / sw;
    double vyy = (syy - sy * sy / sw) / sw;
    // principal axis의 기울기;
    double theta = atan2(2 * vxy, vxx - vyy) / 2;
    line[0] = cos(theta); line[1] = sin(theta);
    // center of mass (xc, yc);
    line[2] = sx / sw; line[3] = sy / sw;
    // line-eq:: sin(theta) * (x - xc) = cos(theta) * (y - yc);
    // calculate weights w.r.t the new line;
    std::vector<double> dist(P.size());
    double scale = 0;
    for (int i = P.size(); i-->0;) {
        double d = dist[i] = DistanceToLine(P[i], line);
        if (d > scale) scale = d;
    }
    if (scale == 0) scale = 1;
    for (int i = dist.size(); i-->0; ) {
        double d = dist[i] / scale;
        weight[i] = 1 / (1 + d * d / 2);
    }
    return fitError(P, line);
};
void test_main(std::vector<CPoint>& pts, double line_params[4]) {
    // initial weights = all equal weights;
    std::vector<double> weight(pts.size(), 1); 
    while (1) {
       double err = LineFit_PCA(pts, weight, line_params) ;
       //(1) check goodness of line-fitting; if good enough, break loop;
       //(2) re-calculate weight, normalization not required.
     }
};

아래 그림은 weight를 구하는 함수로 $weight= 1 /\sqrt{1+dist\times dist}$를 이용하고, fitting 과정을 반복하여 얻은 결과다. 상당히 많은 outlier가 있음에도 영향을 덜 받는다. 파란 점이 outlier이고, 빨간 직선은 outlier가 없는 경우 fitting 결과고, 파란 선은 outlier까지 포함한 fitting 결과다.

##: 네이버 블로그에서 이전;