이미지에서 관찰된 점들이 {(xi, yi) | i = 1, 2,..., N}이 있다. 이 점들을 직선 y = a + b*x 로 피팅을 하고 싶을 때, 보통 최소 자승법을 이용하는 데,  원리는 직선의 방정식이 예측한 y값과 실제 관찰한 y값의 차이의 제곱(=square deviation)을 최소화시키는 직선의 파라미터 a, b를 찾는 것과 동일하다:


확률적으로는, 데이터를 얻는 과정인 측정에서 측정값 yi 는 랜덤에러를 가지게 되고, 이 에러는 참값 y(x) 근방에서 정규분포의 모양을 가지게 된다고 가정을 한다. 만약에 모든 측정의 에러가 동일한 표준편차 σ를 갖게 된다면, N개의 관측데이터가 나타날 확률(likelihood)는 (개별측정은 IID조건을 만족한다고 전제)
     

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


그러나 이 식을 이용하는 경우에는 최소자승법과는 달리, 파라미터 a, b를 주어진 데이터 {(xi, yi)}로 표현할 수 없고, 반복적인 방법을 이용해서 찾는 수 밖에 없다. 


주어진 수열 {c
i}와 a 와의 차이의 절대값의 합
i |ci - a| 는 a가 이 수열의 median 값인 경우에 최소된다는 사실을 이용하면 (증명: 극값을 구하기 위해서 a에 대해서 미분하면, 0 = (i 1) - (i 1) : 합은 a가 ci 보다 큰 경우와 작은 경우로 분리=> 따라서 0이 되기 위해서는 작은 경우와 큰 경우의 수가 같아야 한다 => 즉, a = median{ci}), 고정된 b값에 대해서 위 absolute deviation을 최소화 시키는 a값은

 


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


을 얻는데 위의 a 값을 대입해서 위의 식을 만족시키는 b값을 bracketing and bisection 방법을 이용해서 얻을 수 있다(불연속함수여서, 일반적으로 근을 구하는 방법을 사용하는 것은 위험하다). 아래의 코드는 이를 구현한 것이다.


참고 : Numerical Recipes. 

// 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(double *x, double* y, int n, 

               double bb, double *aa, double *abdev) { 

    double *h = new double [n];

    for (int i = 0; i < n; i++)  h[i] = y[i] - bb * x[i];

    qsort(h, n, sizeof(double), comp); 

    // median; 

    *aa = (n & 1) ? h[n / 2] : (h[n / 2] + h[n / 2 - 1]) / 2;

    

    double sum = 0;  

    *abdev = 0; 

    for (i = 0; i < n; i++) { 

        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]); 

    } 

    delete [] h;

    return sum; // return sum{xi * sign(yi - (b * xi + a))} 

}; 

// 최소자승법을 이용한 라인추정: 

// return (sigma[dy] / sigma[x]);

double FitLine_LS(double *x, double *y, int n, double *a, double *b) { 

    double sx = 0, sy = 0, sxx = 0, sxy = 0; 

    for (int i = 0; i < n; i++) { 

        sx += x[i]; 

        sy += y[i]; 

        sxx += x[i] * x[i] ; 

        sxy += x[i] * y[i] ; 

    }; 

    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 (i = 0; i < n; i++) { 

        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);  ; 

}

// y = a + b * x ; 

// Least Absolute Deviation: 

void FitLine_MAD (double *x, double *y, int n, 

                  double *a, double *b, double *mad) 

    // least square esimates for (aa, bb);      

    double aa, bb; 

    double sigb = FitLine_LS(x, y, n, &aa, &bb); 

    double b1 = bb; 

    double abdev; 

    double f1 = RhoFunc(x, y, n, 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, n, 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, n, 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, n, bb, &aa, &abdev) ; 

        if ((f * f1) >= 0) { 

            f1 = f; 

            b1 = bb; 

        } 

        else { 

            f2 = f; 

            b2 = bb; 

        } 

    } 

    *a = aa; 

    *b = bb; 

    *mad = abdev / n; 

}

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

"사용자

'Image Recognition' 카테고리의 다른 글

RANSAC : Circle Fit  (0) 2008.07.21
KMeans Algorithm  (0) 2008.07.19
Robust Line Fitting  (0) 2008.07.08
Bayesian Spam Filtering  (0) 2008.07.03
EM : Binarization  (0) 2008.07.01
EM Algorithm : Line Fitting 예  (0) 2008.06.29
Posted by helloktk