Geometry & Recognition

점집합을 일반적인 2차 곡선으로 피팅하는 경우에 방정식은

$$ a x^2 + by^2 + cxy +d x + ey +f = 0$$

의 계수를 주어진 데이터를 이용하여서 구해야 한다. 실제 문제에서는 타원, 포물선 쌍곡 선등의 타입에 따라 몇 가지 제약 조건을 넣어 피팅을 한다. 원은 타원의 특별한 경우로 일반적으로 $a = b$, $c = 0$의 제약 조건이 필요하다. 그러나 보다 엄밀하게 제약을 하게 되면 $a = b = 1$의 추가 조건을 줄 수 있다. 이 경우는 점들이 모두 일직선에 있는 경우를 ($a = b = 0$) 취급할 수 없게 된다. 이 예외적인 경우를 제외하고는 최소자승법을 사용하면 계수를 매우 쉽게 구할 수 있기 때문에 많이 이용된다.

문제: 주어진 데이터를 fitting 하는 이차곡선(원)

$$x^2  + y^2 + A x + B  y + C = 0$$

의 계수 $A, B, C$를 최소자승법을 사용해서 구하라. 

주어진 점집합이 원 위의 점이면 우변이 0이 되어야 하나, 실제 데이터를 얻는 과정에서 여러 노이즈에 노출되므로 일반적으로 0이 되지 않는다. 최소자승법은 주어진 점들이 원에서 벗어나는 정도의 제곱 합이 최소가 되도록 하는 계수 $A, B, C$를 결정한다.  원과 점의 편차의 제곱합
$$ L=\sum_ i   \left |x_i^2 + y_i^2 + A x_i + B y_i + C \right|^2 , $$

의 극값을 찾기 위해서 $A, B,$ 그리고 $C$에 대해 미분을 하면

$$\frac{\partial L}{\partial A} = 2 \sum_i (x_i^2 + y_i^2 + A x_i + B y_i + C) x_i = 0, $$

$$\frac{\partial L}{\partial B} = 2 \sum_i (x_i^2 + y_i^2 + A x_i + B y_i + C) y_i = 0, $$

$$\frac{\partial L}{\partial C} = 2 \sum_i (x_i^2 + y_i^2 + A x_i + B y_i + C) = 0. $$

이 연립방정식을 풀면  $A, B, C$를 구할 수 있다. 계산을 단순하게 만들고 수치적 안정성을 높이기 위해 입력점들의 질량중심 

$$ m_x = \frac{1}{N} \sum_i x_i, \quad m_y = \frac{1}{N} \sum_i y_i$$

계에서 계산을 하자. 이를 위해 입력점의 $x$, $y$ 성분에서 각각 $m_x$, $m_y$만큼을 빼준 값을 좌표값으로 대입하면 된다: 

$$ x_i \to x_i - m_x,\quad y_i \to y_i - m_y$$

그러면 질량중심 좌표계에서는 $S_x = \sum_i x_i =0$, $S_y= \sum_i y_i =0$이 된다.

우선 세 번째 식에서 

$$ C = -\frac{S_{x^2} + S_{y^2}}{N} $$

을 얻을 수 있고, 첫번째와 두 번째 식에서는 각각

$$ S_{x^2} A +  S_{xy} B = -  (S_{x^3} + S_{xy^2})  $$

$$ S_{xy}  A  + S_{y^2} B = - (S_{y^3} + S_{x^2 y} )$$

을 얻을 수 있다. 이를 풀면

$$ A = \frac{- S_{y^2} ( S_{x^3} + S_{xy^2}) + S_{xy} (S_{y^3} + S_{x^2y})  }{ S_{x^2} S_{y^2} - S_{xy}^2 } \\ B= \frac{-S_{x^2}(S_{y^3} + S_{x^2 y}) +S_{xy}  (S_{x^3} + S_{xy^2}) }{S_{x^2} S_{y^2}- S_{xy}^2} $$

여기서 주어진 데이터의 각 차수에 해당하는 moment는 아래처럼 계산된다:

추정된 원의 중심 $(c_x, c_y)$는 

$$ c_x = - \frac{A}{2},   \qquad c_y = - \frac{B}{2} $$

로 주어지고, 반지름은 

$$r^2 =  c_x^2 +c_y^2 - C = c_x^2 + c_y^2 + \frac{1}{N}( S_{x^2}+S_{y^2})$$

로 주어진다.

Ref: I. Kasa, A curve fitting procedure and its error analysis. IEEE Trans. Inst. Meas., 25:8-14, 1976

/* 구현 코드: 2024.04.01, typing error 수정 & 질량중심계 계산으로 수정;*/
double circleFit_LS(std::vector<CPoint>& Q, double& cx, double& cy, double& radius) {
    if (Q.size() < 3) return -1;
    double sx2 = 0.0, sy2 = 0.0, sxy  = 0.0;
    double sx3 = 0.0, sy3 = 0.0, sx2y = 0.0, sxy2 = 0.0;
    double mx = 0, my = 0;            /* center of mass;*/
    for (int k = Q.size(); k-->0;)
        mx += Q[k].x, my += Q[k].y;
    mx /= Q.size(); my /= Q.size();
    /* compute moments; */
    for (int k = Q.size(); k-->0;) { /* offset (mx, my)*/
        double x = Q[k].x - mx, xx = x * x;
        double y = Q[k].y - my, yy = y * y;
        sx2  += xx;      sy2  += yy;      sxy  += x * y;
        sx3  += x * xx;  sy3  += y * yy;
        sx2y += xx * y;  sxy2 += yy * x;
    }
    double det = sx2 * sy2 - sxy * sxy;
    if (fabs(det) < 1.e-10) return -1;    /*collinear한 경우임;*/
    /* center in cm frame; */
    double a = sx3 + sxy2;
    double b = sy3 + sx2y;
    cx = (sy2 * a - sxy * b) / det / 2;
    cy = (sx2 * b - sxy * a) / det / 2;
    /* radius squared */
    double radsq = cx * cx + cy * cy + (sx2 + sy2) / Q.size();
    radius = sqrt(radsq);
    cx += mx; cy += my; /* recover offset; */
    return fitError(Q, cx, cy, radius);
}

정수의 제곱근을 정수 연산만 가지고 계산해야 하는 경우가 많이 있는데 (특히, 임베디드 환경에서 프로그래밍을 하는 경우) 이를 구현하는 방법에 대해서 알아보자. 32비트 정수의 제곱근은 16비트 이하이므로 가장 최상위의 비트 값이 세팅이 된 정수를 가지고 제곱근 조건을 만족시키는지 확인하고, 만족시키면 다음의 하위 비트 값이 세팅이 된 정수 값을 기존의 값에 더해서 다시 제곱근 조건을 만족시키는지 확인하는 과정을  마지막 비트까지 반복하면 된다. 그런데 가능한 가장 최상위의 비트 값을 갖는 정수는 0x8000 (= 1<<15) 이므로 양의 정수 x의 제곱근을 구하는 과정은 아래와 같다.

           unsigned s = 0x8000;  //= (1 << 15);
           unsigned t = 0;
           do {                          //새로 세팅된 값이 제곱근 조건을 만족하는가?
                if ((t + s) * (t + s) <= x ) t = t + s;  //만족하면 하위 비트값을 더함;
                s = s >> 1;                                //다음 하위 비트로 내려감;
           } while (s != 0) ;
           return (t);

이 알고리즘은 16번의 반복으로 항상 끝이 나지만, 중간에 곱셈 연산이 들어가 있다. 이 곱셈 연산을 없애는 방법을 알아보자. 먼저           

 (t + s) * (t + s) = t * t + 2 * t * s + s * s;

이므로 if 문 ( ) 내부 식 양변에서 t^2를 빼도록 하자. 이 경우 t의 값이 갱신이 되면

 x - t^2 => x - (told + sold)^2 => x - told * told - (2 * told * sold + sold * sold) ;

이므로, x 값에서도 역시 2 * t * s + s * s의 값을 제거해야 한다. 따라서

                  unsigned nr = 2 * t * s + s * s;
                  if (nr <= x) {
                         x = x - nr;
                         t = t + s;
                  }

로 고쳐야 한다. 이제, s * s 와 2 * t * s을 별도의 변수를 도입하여 컨트롤하면,

          unsigned m = s * s;    // m = 0x40000000;
          unsigned r = 2 * t * s ;// => 2 * (told + sold) * sold
                                        // = 2 * told * sold + 2sold * sold 
                                        // = r + m + m = nr + m;
          do {
                  nr = m + r ;
                  if (nr <= x) {
                         x = x - nr ;
                         r = nr + m;
                         t = t + s;
                  } ;
                  s = s >> 1;
                  m = m >> 2;
                  r = r >> 1;
          } while (s != 0);

로 사용이 가능하다. 그런데 여기서 r = 2 * t * s인데, 마지막으로 루프를 빠져나올 때 s가 1인 상태에서 s / 2를 하므로 (실수 연산을 고려하면) 이때 r 값은 t 값과 같다. 따라서 t 변수를 사용할 필요가 없고, 마찬가지로 s 변수도 사용할 필요가 없어진다 (m을 컨트롤 변수로 이용). 

//16-step integer_sqrt;
unsigned isqrt( unsigned x ) {
       unsigned m = 0x40000000 ;   //=(1 << 30)
       unsigned r = 0, nr;
       do {
            nr = m + r ;
            if (nr <= x) {
                   x -= nr ;
                   r = nr + m;
            }
            m >>= 2;
            r >>= 1;
       } while (m != 0);
       return (r) ;
};

결과는 항상 <= sqrt(x) 이다. 실제로 빠를까?

한번 더 다듬도록 하자. 나중에 r을 시프트 하지 말고 먼저 하면, r = 2 * t * s => r / 2 = t * s 이므로, if 문의 r은 r =  r + m이 된다. 정리하면

unsigned isqrt( unsigned x ) {
       unsigned m = 0 x 40000000 ;   //= (1 << 30)
       unsigned r = 0, nr;
       do {
            nr = m + r ;
            r >>= 1;
            if (nr <= x) {
                   x -= nr ;
                   r += m;
            }
            m >>= 2;
       } while (m != 0);
       return (r) ;
};

이 된다. 이는 잘 알려진 integer sqrt의 값을 구하는 함수이다.

다시 한번 더 형태를 바꾸면, r = r + m 항을 보면, r은 m을 누적하는 변수이고, m은 하나의 비트만 세팅이 된 변수이므로 덧셈을 비트 or 연산으로 바꿀수 있다.

             r += m   ==> r |= m;

그리고 매크로를 도입하여 루프를 풀어쓰기 하면, 아래의 코드를 얻을 수 있다.

unsigned isqrt( unsigned x ) {
       unsigned m = 0x40000000 ;   //= (1 << 30)
       unsigned r = 0, nr;
#define STEP(k)             \
           nr = m + r;      \
           r <<= 1;         \
           if (nr <= x) {   \
                x -= nr;    \
                r |= m;     \
           }                 
       STEP(15) ;  m>>=2;
       STEP(14) ;  m>>=2;
       STEP(13) ;  m>>=2;
       STEP(12) ;  m>>=2;
       STEP(11) ;  m>>=2;
       STEP(10) ;  m>>=2;
       STEP(9)  ;  m>>=2;
       STEP(8)  ;  m>>=2;
       STEP(7)  ;  m>>=2;
       STEP(6)  ;  m>>=2;      
       STEP(5)  ;  m>>=2;
       STEP(4)  ;  m>>=2;
       STEP(3)  ;  m>>=2;
       STEP(2)  ;  m>>=2;
       STEP(1)  ;  m>>=2;
       STEP(0)  ;
       return (r) ;
};

이 코드 또한 인터넷 상에 ARM 개발자 사이트나 게임 개발자 사이트에 널리 알려진 코드이다.