Geometry & Recognition

주어진 데이터를 잘 피팅하는 직선을 찾기 위해서는 데이터를 이루는 각 점의 $y$ 값과 같은 위치에서 구하려는 직선과의 차이 (residual)의 제곱을 최소화시키는 조건을 사용한다. 그러나 직선의 기울기가 매우 커지는 경우에는  데이터에 들어있는 outlier에 대해서는 그 차이가 매우 커지는 경우가 발생할 수도 있어 올바른 피팅을 방해할 수 있다. 이를 수정하는 간단한 방법 중 하나는 $y$값의 차이를 이용하지 않고 데이터와 직선 간의 최단거리를 이용하는 것이다. 

수평에서 기울어진 각도가 $\theta$이고 원점에서 거리가 $s$인 직선의 방정식은 

$$ x \sin \theta - y \cos \theta + s=0$$

이고, 한 점 $(x_i, y_i)$에서 이 직선까지 거리는

$$ d_i = | \sin \theta y_i - \cos \theta x_i + s|$$

이다. 따라서 주어진 데이터에서 떨어진 거리 제곱의 합이 최소인 직선을 구하기 위해서는 다음을 최소화시키는 $\theta, s$을 구해야 한다. $$ L = \sum_i \big( \sin \theta x_i - \cos \theta y_i +s \big)^2 \\ (\theta, s)=\text{argmin}(L)$$

$\theta$와 $s$에 대한 극값 조건에서 

$$\frac{1}{2} \frac{\partial L}{\partial \theta} = \frac{1}{2} \sin 2 \theta \sum_i (x_i^2 - y_i^2) - \cos 2 \theta \sum x_i y_i + s \sin \theta \sum_i x_i + s \cos \theta \sum_i y_i = 0$$

$$ \frac{1}{2}\frac{\partial L}{\partial s}=\cos \theta \sum y_i  - \sin \theta \sum_i x_i - N s=0$$

주어진 데이터의 질량중심계에서 계산을 수행하면 $\sum_i x_i = \sum_i y_i =0$ 이므로 데이터의 2차 모멘트를 $$ A= \sum_i (x_i^2 - y_i^2), \qquad B = \sum_i x_i y_i $$로 놓으면 직선의 파라미터를 결정하는 식은

$$ \frac{1}{2} A \sin 2 \theta   - B \cos 2 \theta = 0  \qquad \to \quad  \tan 2\theta = \frac{2B}{A} \\ s = 0$$

두 번째 식은 직선이 질량중심(질량중심계에서 원점)을 통과함을 의미한다. 첫번째 식을 풀면

$$ \tan \theta = \frac{- A \pm \sqrt{A^2 + (2B)^2 }}{2B}$$

두 해 중에서 극소값 조건을 만족시키는 해가 직선을 결정한다. 그런데

$$ \frac{1}{2}\frac{\partial^2 L}{\partial \theta^2}=  A \cos 2 \theta + 2B \sin 2 \theta = \pm \sqrt{A^2 + (2B)^2}  >0$$

이므로 위쪽 부호로 직선($x\sin \theta = y\cos \theta$)이 정해진다. 질량중심계에서는 원점을 지나지만 원좌표계로 돌아오면 데이터의 질량중심을 통과하도록 평행이동시키면 된다.

$$  \left(-A+ \sqrt{A^2+ (2B)^2} \right)  (x-\bar{x}) = 2B (y - \bar{y})  $$

여기서 주어진 데이터의 질량중심은 원좌표계에서

$$ \bar{x} = \frac{1}{N} \sum_i x_i, \quad \bar{y} = \frac{1}{N} \sum_i y_i$$

이다. 또한 원좌표계에서 $A$와 $B$의 계산은 

$$ A = \sum_i [ (x_i - \bar{x})^2 - (y_i - \bar{y})^2], \qquad B = \sum (x_i  - \bar{x})(y_i - \bar{y})$$ 

이 결과는 데이터 분포에 PCA를 적용해서 얻은 결과와 동일하다. PCA에서는 공분산 행렬의 고유값이 큰 쪽에 해당하는 고유벡터가 직선의 방향을 결정했다. https://kipl.tistory.com/211.  또한 통계에서는 Deming regression이라고 불린다.

평면 위에 점집합이 주어지고 이들을 잘 기술하는 직선의 방정식을 구해야 할 경우가 많이 발생한다. 이미지의 에지 정보를 이용해 선분을 찾는 경우에 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 결과다.

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

이미지에서 관찰된 점집합이 $\{(x_i, y_i)| i = 1, 2,\dots, N\}$이 있다. 이 점집합을 직선 $y = a + bx$ 로 피팅을 하고 싶을 때, 보통 최소자승법을 이용하는데, 원리는 직선의 방정식이 예측한 $y$값과 실제 관찰한 $y$값의 차이의 제곱(=square deviation)을 최소화시키는 직선의 기울기 $a$와 절편 $b$를 찾는 것이다:

$$\chi^2(a, b) = \sum_i |y_i - (b x_i +a) |^2 $$

데이터를 얻는 측정 과정에서 측정값 $y_i$는 랜덤 노이즈를 포함하게 되고, 이 노이즈는 참값 $y(x)$ 근방에서 정규분포를 한다고 가정을 할 수 있다. 만약 모든 측정의 노이즈가 동일한 표준편차 $\sigma$를 갖게 된다면, $N$개의 관측 데이터가 나타날 확률(likelihood)은 (개별 측정은 IID 조건을 만족한다고 전제)

$$P = \prod_{i} e^{ -\frac{ | y_i -  (bx_i + a)|^2 }{2\sigma^2 }  }$$

의 형태가 된다. 따라서 최소자승법은 이 likelihood를 최대화시키는 파라미터를 찾는 방법이 된다. 최소자승법은 피팅 파라미터를 주어진 데이터를 이용해서 표현할 수 있는 장점은 있지만, outliers에 대해서 매우 취약하다 (아래의 결과 그림을 참조). 이는 적은 수의 outlier도 χ2에 큰 기여를 할 수 있기 때문이다. 따라서 피팅을 좀 더 robust 하게 만들기 위해서는 outliers가 likelihood에 기여하는 정도를 줄여야 한다. 이를 위해서는 likelihood의 지수 인자를 큰 에러에서 덜 민감하게 반응하는 꼴로 바뀌어야 한다. 이를 만족하는 가장 간단한 것 방법 중 하나가 square-deviation 대신에 absolute-deviation을 이용하는 것이다:

$$\text{absolute deviation} = \sum _i | y_i - (bx_i + a)|   .$$

그러나 이 식을 사용하는 경우에는 최소자승법과 다르게 기울기 $a$와 절편 $b$를 주어진 데이터 $\{(x_i, y_i)\}$로 표현할 수 없고, 반복적인 방법을 이용해서 찾아야 한다. 

수열 $\{c_i\}$에 대해 합 $\sum_{i} |c_i - a|$은 $a$가 수열의 median 값일 때 최솟값을 갖는다는 사실을 이용하면 (증명: 극값을 구하기 위해서 $a$에 대해서 미분하면, $0=(\sum_{c_i > a} 1)-(\sum_{c_i < a} 1)$: 합은 $a$가 $c_i$ 보다 큰 경우와 작은 경우로 분리. 따라서 0이 되기 위해서는 작은 경우와 큰 경우의 수가 같아야 한다. 고로, $a = \text{median}\{c_i\}$ q.e.d.). 고정된 절편 $b$에 대해서 absolute deviation을 최소로 만드는 기울기 $a$는

$$a= \text{median} \{ y_i - b x_i\}$$

임을 알 수 있다. 그리고  absolute deviation 식을 절편 $b$에 대해서 미분해서

$$0 = \sum_i \text{sign} \left( y_i - (bx_i +a) \right)$$

을 얻는데, 위에서 구한 기울기 $a$를 대입한 후 bracketing and bisection 방법을 이용하면 절편 $b$를 얻을 수 있다(불연속 함수이므로 일반적으로 근을 구하는 방법을 사용하는 것은 위험하다). 아래의 코드는 이를 구현한 것이다.

double FitLine_LS(std::vector<double>& x, std::vector<double>& y, double *a, double *b);

// 최소자승법을 이용한 직선 추정:
// return (sigma[dy] / sigma[x]);
double FitLine_LS(std::vector<double>& x, std::vector<double>& y, double& a, double& b) {
    double sx = 0, sy = 0, sxx = 0, sxy = 0;
    for (int i = x.size(); i-->0;) {
        sx  += x[i];        sy  += y[i];
        sxx += x[i] * x[i]; sxy += x[i] * y[i];
    };
    const int n = x.size();
    double det = n * sxx - sx * sx;
    if (det == 0.) return -1;                   // vertical line;
    a = (sxx * sy - sx * sxy) / det;
    b = (n * sxy - sx * sy) / det;
    double chi2 = 0;
    for (int i = x.size(); i-->0;) {
        double t = y[i] - (*a + *b * x[i]);
        chi2 += t * t;
    }
    det /= n;         //det -> var(x) * n;
    // chi2 = var(dy) * n;
    // (dy vs x의 편차비)
    return  sqrt(chi2 / det);
}

// qsort-comparator;
static int comp(const void *a, const void * b) {
    double d = *((double *)a) - *((double *)b);
    return d > 0 ? 1 : d < 0 ? -1 : 0;
}
// 기울기(bb)가 주어진 경우에 y-절편(median = aa)값을 추정하고, 이 때 AD값을 얻는다.
double RhoFunc(std::vector<double>&x, std::vector<double>& y,
               double bb, double& aa, double& abdev) {
    std::vector<double> h(x.size());
    for (int i = x.size(); i-->0;) h[i] = y[i] - bb * x[i];
    qsort(&h[0], h.size(), sizeof(double), comp);
    // median;
    const int med = h.size()/2;
    aa = (h.size() & 1) ? h[med] : (h[med] + h[med-1])/2;

    double sum = 0;
    abdev = 0;
    for (int i = x.size(); i-->0;) {
        double d = y[i] - (bb * x[i] + aa);
        abdev += fabs(d);
        // sign-함수의 원점에서 모호함을 없애기 위해서 증폭을 시킴;
        if (y[i] != 0.) d /= fabs(y[i]);
        if (fabs(d) > DBL_EPSILON) // sign 함수의 모호함에서 벗어나면,
            sum += (d >= 0 ? x[i] : -x[i]);
    }
    return sum; // return sum{xi * sign(yi - (b * xi + a))}
};
// y = a + b * x ;
// Least Absolute Deviation:
double FitLine_MAD (std::vector<double>& x, std::vector<double>& y,
                    double& a, double& b) {
    // least square esimates for (aa, bb);
    double aa, bb, abdev;
    double sigb = FitLine_LS(x, y, aa, bb);   // estimate: y=aa + bb*x;
    double b1 = bb;
    double f1 = RhoFunc(x, y, b1, aa, abdev);
    /* bracket 3-sigma away in the downhill direction;*/
    double b2 = fabs(3 * sigb);
    b2 = bb + (f1 < 0 ? -b2 : b2);
    double f2 = RhoFunc(x, y, b2, aa, abdev);

    // if conditional added to take care of the case of a
    // line input into this function. It is needed to avoid an
    // infinite loop when (b1 == b2) (within floating point error)
    if (fabs(b2 - b1) > (sigb + 0.005)) {
        // bracketing;
        while ((f1 * f2) > 0) {
            bb = 2 * b2 - b1;
            b1 = b2; b2 = bb; 
            f1 = f2; f2 = RhoFunc(x, y, b2, aa, abdev) ;
        }
    }
    // refine until the error is a negligible number of std;
    sigb *= 0.01;
    while (fabs(b2 - b1)> sigb) {
        // bisection;
        bb = (b1 + b2) / 2.;
        if ((bb == b1) || (bb == b2)) break ;
        double f = RhoFunc(x, y, bb, aa, abdev) ;
        if ((f * f1) >= 0) {
            f1 = f; b1 = bb;
        } else {
            f2 = f; b2 = bb;
        }
    }
    a = aa; b = bb; 
    return (abdev/x.size());
}

// 붉은 선--> 최소자승법에 의한 피팅.: outlier에 매우 취약한 구조.
// 파란 선--> least absolute deviation을 이용한 피팅: outlier에 매우 robust 하다.

아래 그림을 보면 우리는 데이터를 연결하는 두 개의 직선을 생각할 수 있을 것이다.  그럼 두 개의 직선을 어떻게 얻을 것인가? 물론, ICA(independent component analysis)를 이용하는 것이 한 가지 방법이 될 것이다. 여기서는 EM 알고리즘을 이용하여서 두 개의 직선을 기술하는 기울기와 $y$-절편의 값을 구하는 방법을 알아보자.

두 직선을 각각 $y=a_1 x + b_1$, $y = a_2  x + b_2$라고 하면, $(a_1, b_1)$, $(a_2, b_2)$를 구하는 문제이다. 만약 각각의 data가 에러가 수반되는 측정에 의해서 얻어졌다고 하자. 에러 분포가 정규분포를 따른다면(여기서는 두 개의 모델 모두 갖은 표준편차를 갖는다고 가정했다) $i$-번째의 데이터가 각각의 직선 모델 1과 2에 속할 확률은 (posterior with equal prior) Bayes 공식에 의해서 

$$ w_1[i] = \frac{ e^{ - \frac{ r_1^2[i]}{2\sigma^2} } } {e^{ - \frac{ r_1^2[i]}{2\sigma^2} } + e^{ - \frac{ r_2^2[i]}{2\sigma^2} } }, \quad w_1[i] = \frac{ e^{ - \frac{ r_2^2[i]}{2\sigma^2} } } {e^{ - \frac{ r_1^2[i]}{2\sigma^2} } + e^{ - \frac{ r_2^2[i]}{2\sigma^2} } }, \quad i = 0,1,2,... $$

로 주어진다. 여기서 $r_1(i)$와 $r_2(i)$는 residual error이다:

$$r_1[i] = a_1 x[i] + b_1 - y[i],\quad r_2[i] = a_2 x[i] + b_2 - y[i], \quad i=0,1,2,...$$

(*이 값 대신에 직선까지 거리=$\frac{|a_k x + b_k - y|}{\sqrt{1+ a_k^2}},~ k=1,2)$로 대체해도 된다)

이제 각각의 데이터에 대해서 posterior를 구하였으므로(E-step) 이 값을 가중치로 하여서 직선의 방정식을 다시 갱신한다. 즉, 각각의 data 점들에 대한 $w_1(i)$를 가중치로 하여서 다시 직선 모델 1의 파라미터를 재계산하고, 동일한 과정을 $w_2(i)$를 가지고 직선 모델 2를 계산하여서 $(a_1, b_1)$, $(a_2, b_2)$를 재계산한다(M-step). 이 갱신된 파라미터를 이용하여서 다시 가중치를 구하는 작업과, 직선의 파라미터를 구하는 작업을 일정하게 수렴할 때까지 반복을 하는 과정을 수행한다.

아래의 그림은 3번 만에 원하는 결과를 얻는 것을 보여준다. 직선의 파라미터에 대한 초기값과 residual error의 표준편차 파라미터에 대한 적절한 값의 선택이 중요하다.

// 코드의 일부...
std::vector<CPoint> data ;                             // data,
std::vector<double> w1(data.size()), w2(data.size());  // weights.
double a1, b1, a2, b2 ;                                // line params;
double sigma ;
// E-step;
void calcWeights() {
    for (int i = data.size(); i-- > 0;) {
        double  x = data[i].x, y = data[i].y ;
        double r1 = a1 * x + b1 - y ;
        double r2 = a2 * x + b2 - y ;
        double n1 = SQR(r1) / SQR(sigma) / 2;
        double n2 = SQR(r2) / SQR(sigma) / 2;
        double p1 = exp( - n1);
        double p2 = exp( - n2);
        w1[i] = p1 / (p1 + p2);
        w2[i] = p2 / (p1 + p2);
    }
};
//  M-step
void estimModels() {
    double s1xx = 0, s1x = 0, s1 = 0, s1y = 0, s1xy = 0;
    double s2xx = 0, s2x = 0, s2 = 0, s2y = 0, s2xy = 0;
    for (int i = data.size(); i-- > 0;) {
        double  x = data[i].x,
                y = data[i].y;
            s1xx += w1[i] * SQR(x);
            s1xy += w1[i] * x * y;
            s1x  += w1[i] * x;
            s1y  += w1[i] * y;
            s1   += w1[i];
            //
            s2xx += w2[i] * SQR(x);
            s2xy += w2[i] * x * y;
            s2x  += w2[i] * x;
            s2y  += w2[i] * y;
            s2   += w2[i];
    };
    double det1 = s1xx * s1 - SQR(s1x);
    double det2 = s2xx * s2 - SQR(s2x);
    a1 = (s1 * s1xy  - s1x * s1y ) / det1;
    b1 = (s1xx * s1y - s1x * s1xy) / det1;
    a2 = (s2 * s2xy  - s2x * s2y ) / det2;
    b2 = (s2xx * s2y - s2x * s2xy) / det2;
}