Geometry & Recognition

평면 위에 점집합이 주어지고 이들을 잘 기술하는 직선의 방정식을 구해야 할 경우가 많이 발생한다. 이미지의 에지 정보를 이용해 선분을 찾는 경우에 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)를 이용할 수 있다. 점들이 선분을 구성하는 경우, 선분 방향으로는 점 위치의 편차가 크지만 수직 방향으로는 편차가 상대적으로 작다. 따라서 평면에서 점 분포에 대한 공분산 행렬 $\mathtt{C}\mathtt{o}\mathtt{v}$의 고윳값과 고유 벡터를 구하면, 큰 고윳값을 갖는 고유 벡터 방향이 선분의 방향이 될 것이다.

$$\mathtt{C}\mathtt{o}\mathtt{v}[\{(x_i,y_i)\}]=\frac{1}{N}\left(\begin{array}{cc}\sum (x_i-\overline{x}{)}^2& \sum (x_i-\overline{x})(y_i-\overline{y})\\ \sum (x_i-\overline{x})(y_i-\overline{y})& \sum (y_i-\overline{y}{)}^2\end{array}\right)$$

잘 피팅이 이루어지려면 두 고윳값의 차이가 커야 한다. 또한 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]) {
    // initial setting: weight[i] = 1.;
    // compute weighted moments;
    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 결과다.

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

서로 다른 조건에서 생성이 된 두 개의 영상을 비교하기 위해서는 먼저 영상을 정규화시키어야 한다. 쉽게 생각할 수 있는 정규화의 방법은 두 개의 영상이 동일한 히스토그램을 같도록 만들면 된다. 어떻게 동일한 히스토그램을 만들 수 있는가? 간단히 생각할 수 있는 방법은 주어진 영상의 히스토그램을 전 픽셀에 걸쳐서 균일한 빈도를 보이는 히스토그램이 되도록 변환을 하는 것이다. 즉, 히스토그램을 평탄화시키는 것이다. 이 변환을 $F(x)$라고 하면, 입력 영상의 픽셀 값 $x$는 새로운 픽셀 값 $y=F(x)$로 바뀐다. 영상의 히스토그램은 픽셀 값의 빈도를 나타내므로 확률 밀도 함수로 해석할 수 있고 (픽셀 값을 연속 변수로 볼 때) 주어진 영상의 확률 밀도 함수를 $P(x)$라고 하자. 그러면 $P(x)dx$는 픽셀 값이 $[x,x+dx]$인 픽셀이 영상에서 나타날 확률을 의미한다. 변환 후에서는 임의의 픽셀 값이 균일해야 하므로 확률 밀도 함수는 상수 $a$이고, 픽셀 값의 구간이 $[0,255]$ 이므로 $a=1/255$로 주어진다. 따라서, 변환이 확률을 보존하기 위해서는 $$P(x)dx=ady$$을 만족해야 하므로 변환 후 픽셀 값 y와 변환 전 픽셀 값 x사이의 관계를 알 수 있다. $$\frac{dy}{dx}=\frac{1}{a}P(x)$$ 이 식을 적분하면, $$y=F(x)=\frac{1}{a}{\int }_0^{x}P(x^{\prime })d{x}^{\prime }=255\times \text{cumulative sum of }P(x)$$ 즉, 히스토그램을 평탄화시키는 변환은 원래 영상 히스토그램의 누적합을 전체 픽셀 수로 나눈 값,

$$k\to F(k)=255\times \frac{\sum _{i=0}^{k}H[i]-H[0]}{\sum _{i=0}^{255}H[i]-H[0]}$$으로 결정된다($0\to 0$이도록 $H[0]$을 뺀다). 영상의 히스토그램은 이산적이기 때문에 실제로 모든 픽셀 값에서 균일하게 나타나지 않지만, 구간 변환은 거의 일정하게 나타난다. 이 변환을 거친 영상의 누적 히스토그램(cumulative histogram)을 $y$-축 픽셀 값을 $x$-축으로 해서 그리면 거의 일직선으로 증가하는 모양을 나타낸다.

이 히스토그램의 평탄화는 픽셀 값이 좁은 영역에 몰려있는 영상을 전 영역에 분포하도록 한다. 인간은 영상의 밝기에 의해서 보다는 밝고 어두움의 변화의 크기에 의해서 인지도가 증가하게 되는데, 평탄화된 영상은 보다 쉽게 인식이 될 수 있는 구조를 가진다.

std::vector<BYTE> histogram_eq(const std::vector<BYTE> &src) {
    int hist[256] = {0}, cum[256], map[256];
    for (int k = src.size(); k-->0;) ++hist[src[k]];
    for (int k = 0, s = 0; k < 256; k++) {
    	s += hist[k]; cum[k] = s;
    }
    double factor = double(255) / (cum[255] - cum[0]) ;
    for (int k = 0; k < 256; k++) {
        int x = int(factor * (cum[k] - cum[0]));
        map[k] = x > 255 ? 255: x;
    }
    std::vector<BYTE> dst(src.size());
    for (int k = src.size(); k-->0;) dst[k] = map[src[k]];
    return dst;
}

점집합을 일반적인 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