Geometry & Recognition

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

$$a{x}^2+b{y}^2+cxy+dx+ey+f=0$$

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

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

$$x^2+y^2+Ax+By+C=0$$

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

주어진 점집합이 원 위의 점이면 우변이 0이 되어야 하나, 실제 데이터를 얻는 과정에서 여러 노이즈에 노출되므로 일반적으로 0이 되지 않는다. 최소자승법은 주어진 점들이 원에서 벗어나는 정도의 제곱 합이 최소가 되도록 하는 계수 $A,B,C$를 결정한다.  원과 점의 편차의 제곱합
$$L=\sum _i{|x_i^2+y_i^2+A{x}_i+B{y}_i+C|}^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,m_y=\frac{1}{N}\sum _i{y}_i$$

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

$$x_i\to x_i-m_x,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_{x{y}^2})$$

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

을 얻을 수 있다. 이를 풀면

$$A=\frac{-S_{y^2}(S_{x^3}+S_{x{y}^2})+S_{xy}(S_{y^3}+S_{x^{2}y})}{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_{x{y}^2})}{S_{x^2}{S}_{y^2}-S_{xy}^2}$$

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

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

$$c_x=-\frac{A}{2},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);
}

https://kipl.tistory.com/357

https://kipl.tistory.com/32

히스토그램처럼 균일한 간격으로 주어진 데이터에서 피크 값의 위치를 찾기 위해서는 먼저 노이즈 제거를 위해 몇 번의 smoothing 과정을 적용해야 한다 (구체적인 방법은 필요에 따라 다양하다). smoothing 과정을 거친 히스토그램에서 피크 위치를 찾는 것은 쉬운 작업이다. 그런데 경우에 따라 그 위치를 sub-order까지 구해야 필요가 생긴다. 예를 들면, 실수 값 데이터 시퀀스 대한 히스토그램을 만들려면 먼저 정수로 양자화시켜야 가능하다. 이 양자화된 정보를 담고 있는 히스토그램에서 피크의 위치를 찾으려 할 때 양자화 과정에서 잃어버린 정밀도를 복원하려면 interpolation을 써야 한다. 간단하게 parabolic 근사로 피크의 위치를 찾는 경우를 생각하자. 이 방법은 수학적으로 단순하지만 영상처리 알고리즘에서 많이 이용이 되고 있다. 데이터 시퀀스가 균일한 간격에서 주어졌으므로 계산은 $-1,0,1$의 세 군데 위치에서 중앙 근방에 피크가 나타나는지를 고려하면 된다. 세 점에서 히스토그램 값이 각각 $h_m,h_0,h_p$일 때 이들을 지나는 이차 곡선의 피크 위치를 찾자. 주어진 이차 곡선을

$$y=a(x-c{)}^2+b$$

꼴로 쓰면 $c$가 0에서 벗어난 정도를 나타낸다.

$$(-1,h_m)\to h_m=a(-1-c{)}^2+b$$

$$(0,h_0)\to h_0=a{c}^2+b$$

$$(+1,h_p)\to h_m=a(+1-c{)}^2+b$$

이므로

$$h_m-h_p=4ac,$$

$$h_m+h_p-2h_0=2a,$$

$$\therefore c=\frac{h_m-h_p}{2(h_m+h_p-2h_0)}$$

아래 코드는 피크의 위치가 중앙점에서 벗어난 정도를 준다.

bool parabolicInterpolate(double hm, double h0, double hp, double *c) {
      double a = hm + hp - 2 * h0;
      if (a >= 0) return false; // not a parabola(a==0), not convex up(a>0);
      *c = 0.5 * (hm - hp) / a;
      if (*c < -0.5 || *c > 0.5) return false; //  too far;
      return true;
}