PCA Line Fitting

평면 위에 주어진 점집합을 가장 잘 나타내는 직선을 구해야 하는 경우는 영상 처리와 컴퓨터 비전에서 자주 발생한다. 예를 들어 이미지의 에지 점들로부터 선분을 추출할 때는 Hough Transform을 사용할 수도 있지만, 점집합에 대해 직접 수치적으로 직선을 피팅(fitting)하는 방법도 널리 사용된다.

실제 데이터에는 측정 오차와 잡음이 포함되므로, 직선 위의 점뿐만 아니라 직선에서 크게 벗어난 outlier가 함께 존재하는 경우가 많다. 따라서 직선 피팅 알고리즘은 이러한 outlier에 대해 가능한 한 강인(robust)해야 한다. 일반적으로는 먼저 대략적인 직선을 추정한 뒤, 이를 이용하여 점차 정확도를 높여 가는 반복적인(iterative) 방법을 사용한다.

입력 데이터가 \({(x_i,y_i)\mid i=0,\ldots,N-1}\)로 주어질 때, 가장 널리 사용되는 최소자승법(Least Squares)은 각 \(x_i\)에서 직선이 예측하는 \(y\)값과 실제 \(y_i\)의 차이(residual)의 제곱합을 최소로 하는 직선의 기울기와 절편을 구한다. 그러나 이 방법은 직선을 \(y=ax+b\)의 형태로 가정하므로 수직선이나 이에 가까운 직선을 적절히 표현하지 못한다. 또한 점과 직선 사이의 직교거리가 아니라 \(y\)방향의 오차만을 최소화하며, 잔차를 제곱하기 때문에 outlier의 영향을 크게 받는다.

이러한 문제를 완화하기 위해 PCA(Principal Component Analysis)를 이용할 수 있다. 점들이 하나의 선분을 따라 분포한다면 선분 방향으로는 점들의 분산이 크고, 이에 수직인 방향으로는 분산이 작다. 따라서 점 분포의 공분산 행렬

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

의 고윳값과 고유벡터를 구하면, 가장 큰 고윳값에 대응하는 고유벡터가 직선의 방향을 나타내고, 가장 작은 고윳값에 대응하는 고유벡터는 직선의 법선 방향을 나타낸다. 따라서 PCA는 수직선을 포함한 임의의 방향의 직선을 자연스럽게 피팅할 수 있다.

그러나 PCA 역시 outlier에는 민감하므로 보다 강인한 피팅을 위해서는 각 점에 서로 다른 가중치를 부여하는 반복적인 방법을 사용한다. 먼저 모든 점에 동일한 가중치를 주어 직선의 방향을 추정한 뒤, 이후에는 직선에서 먼 점일수록 작은 가중치를 부여하여 다시 공분산 행렬을 계산한다. 이러한 과정을 반복하면 outlier의 영향을 점차 줄이면서 보다 안정적인 직선을 얻을 수 있다. 가중치는 일반적으로 점과 직선 사이의 직교거리에 따라 결정되며, Huber, Tukey biweight 등 다양한 robust 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를 구하는 함수로 $\text{weight}= 1 /\sqrt{1+(\text{dist/scale})^2} $를 이용하고, fitting 과정을 반복하여 얻은 결과다. 상당히 많은 outlier가 있음에도 영향을 덜 받는다. 파란 점이 outlier이고, 빨간 직선은 outlier가 없는 경우 fitting 결과고, 파란 선은 outlier까지 포함한 fitting 결과다.

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