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

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

 
 
 
 
728x90

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

Fast Float Sqrt  (0) 2020.12.27
2차원 Gaussian 분포 생성  (0) 2020.12.10
Histogram Equalization  (0) 2020.11.12
Least Squares Fitting of Circles  (0) 2020.11.11
Integer Sqrt  (0) 2020.11.11
Posted by helloktk
,

이미지에서 관찰된 점집합이 $\{(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 하다.

728x90

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

RANSAC: Circle Fit  (0) 2008.07.21
KMeans Algorithm  (0) 2008.07.19
EM: Binarization  (0) 2008.07.01
EM Algorithm: Line Fitting  (0) 2008.06.29
Gaussian Mixture Model  (2) 2008.06.07
Posted by helloktk
,

주어진 영상에서 선분을 찾고자 하면 어떻게 해야 할까. 우선 에지 영상을 만들어 선분을 강조하고, 선분이 하나만 있는 경우에는 에지 데이터를 가지고 line-fitting 알고리즘을 적용해서 해당 선분을 기술하는 파라미터를 추출할 수 있다. 그러나 영상이 선분을 여러 개 포함하는 경우에는 이 방법을 이용하기 어렵다(불가능한 것은 아니지만 노이즈 등에 의한 영향을 고려해야 하므로 복잡한 과정을 거쳐야 한다). 평면에서 선분은 원점까지 거리와 그것의 수직(수평)인 방향만 주어지면 결정이 된다. 직선에서 원점까지 거리를 $r$, 법선이 $x$-축과 이루는 각도가 $θ$면

$$r = \cos(θ)  x + \sin(θ)  y;$$

의 형태로 주어진다. 즉, $(r,θ)$ 한쌍이 주어지면 직선 한 개가 정의된다 (주어진 $(r,θ)$ 값이 만드는 직선을 $a  x + b  y = c$  꼴로 쓰면 계수는 $a = \cos θ$, $b = \sin θ$, $c = r$로 표현된다). 

이 직선 방정식은 $rθ$-평면의 stripe=$[0,∞) \times  [0,2π]$을 $xy$-평면으로 보내는 변환으로도 생각할 수 있다. 이 변환은 $rθ$-평면의 한 점을 $xy$-평면의 한 직선으로 보낸다. 역으로는 $xy$-평면의 한 점은 $rθ$-평면의 한 곡선으로 변환이 된다. 직선 상의 점들은 같은 원점 거리=$r$, 같은 기울기=$\theta$를 가지므로, 직선 위의 각 점들이 $rθ$-평면에 만드는 곡선들은 공통의 교점을 가지게 된다. 이 교점의 위치가 $xy$-평면에서 직선을 기술하는 파라미터를 준다.

Hough Transform은 이 변환의 특성을 이용한 것이다. 이미지에서 에지에 해당하는 점들을 위의 변환에 의해서 $rθ$-평면의 곡선으로 보내는 것이다 (프로그램적으로는, 에지점 좌표 $(x, y)$가 주어지면 $θ=0 \rightarrow π$까지 일정 간격으로 증가시키면서 곡선 $r=x \cos(θ) + y\sin(θ)$의 값을 구해서 메모리 상의  $[r,θ]$ 지점의 도수를 1씩 증가시킨다). 만약 에지에 해당되는 점들이 직선의 관계를 가지게 되면, 각 점들에 해당하는 $rθ$-평면에서 곡선들은 한 교점을 형성하게 될 것이다. 긴 선분은 $rθ$-평면의 한 점에서 더 많은 곡선들의 교점을 형성하게 된다. 따라서, $rθ$-평면에서 이러한 교점의 히스토그램을 구하여, 교점 수가 일정 이상 누적이 된 경우를 취하면, $xy$-이미지 평면에서 일정한 길이 이상을 갖는 선분을 골라낼 수 있다.

그러면 왜 $rθ$-평면으로 변환을 생각하는 것일까? 직선이 $y= a x+b$의 형태로 표시되므로 기울기와 $y$축과의 교점을 기술하는 $ab$-공간을 이용할 수도 있지만, 이 경우에 $y$-축에 평행인 직선의 경우에 기울기 $a$ 값이 무한히 커지는 경우가 발생하여 저장공간 할당에 문제가 생긴다. 실제적인 문제에서 되도록이면 compact 한 공간으로 변환을 고려해야 하는 것이다. $rθ$-공간으로의 변환은 유한한 이미지의 경우 항상 유한한 $rθ$-공간의 영역으로 변환된다.

 

아래의 그림은 $x-y=-5$인 직선상의 세 점 $(5,10)$, $(10,15)$, $(15,20)$에 해당하는 $rθ$-평면상의 세 곡선을 보인 것이다: $r = 5 \cos(θ) + 10\sin(θ)$, $r = 10 \cos(θ) + 15  \sin(θ)$, $r = 15 \cos(θ) + 20 \sin(θ)$. 이 세 곡선은 $(r,θ) = (5/\sqrt{2}, 3 \pi /4)$ 지점에서 만남을 알 수 있다. 이 만나는 지점을 구하면, 거꾸로 세 점이  $x-y=-5$ 인 직선 위의 점들이었다는 것을 추정할 수 있다.

 

사용자 삽입 이미지

 

BOOL HoughTransformLine(BYTE* image, int width, int height, /*background=0*/
                                       double rho/*=2*/, 
                                       double theta/*=Pi/180.*/,
                                       int threshold/*=20*/,
                                       std::vector<Line>& vecLine) {
    /* use cosine and sine look-up tables*/
    double irho = 1 / rho;
    int numangle = (int) (Pi / theta);
    int numrho = (int) (((width + height) * 2 + 1) / rho);
    int *accum = (int*)malloc( sizeof(accum[0]) * (numangle + 2) * (numrho + 2) );
    memset( accum, 0, sizeof(accum[0]) * (numangle + 2) * (numrho + 2) );
    double ang;
    int n, rwidth = (numrho + 2);
    for (int  y = 0; y < height; y++ ){
        for (int x = 0; x < width; x++ ){
            if ( image[y * width + x] != 0 ){ // only for foreground pixels;
                for (ang = 0, n = 0; n < numangle; ang += theta, n++ ) {
                    double rho = (x * cos(ang) + y * sin(ang)) * irho ;
                    int r = rho >= 0 ? int(rho + .5) : (-int(-rho + .5));//round to nearest int;
                    // accum을 이차원배열로 생각하고, 1 픽셀의 border를 두는 형태로 한다.
                    // 이는 아래의 local-maxima를 찾을때 경계를 고려할 필요가 없어서 유리하다.
                    r += (numrho - 1) / 2;
                    accum[(n + 1) * rwidth + r + 1]++;
                }
            }
        }
    }
    //find local maxima;
    for (int  r = 0; r < numrho; r++ ) {
        for (int  n = 0; n < numangle; n++ ) {
            int base = (n + 1) * rwidth + r + 1;
            //test whether it is local-maximum(4-way);
            if ( accum[base] > threshold &&
                 accum[base] > accum[base - 1] && accum[base] > accum[base + 1] &&
                 accum[base] > accum[base - rwidth] && accum[base] > accum[base + rwidth] )
            {               
                Line line;
                line.rho = (r - (numrho - 1) * .5) * rho;
                line.angle = n * theta;
                vecLine.push_back( line );
            }
        }
    }
    free(accum);
    return TRUE;
}

아래 그림에서 첫 번째는 소스 이미지이고, 이 이미지에서 $rθ$-평면으로 Hough Transform 한 것이 두 번째 것이다(rho=2, theta=Pi/360). 원본 이미지에는 8개의 선분이 있고, 4 개씩 같은 방향을 갖고 있다. Hough Transform 된 이미지에서 이러한 특징을 볼 수 있다(가로축이 $θ$이고 세로축이 $r$이다. $r$ 축은 가운데가 $r=0$이다). 누적이 많이 된 부분이 두 군데의 동일한 θ에서 나타나고 있다. 그러나 결과에서 8개의 피크를 분리하기는 쉽지 않다. 위의 코드는 각각의 점에서 4방향으로 체크하여 극대 값을 찾고 있으나, 항상 잘 동작하는 것은 아니다.

local-maxima를 좀 더 잘 잡고 싶으면, 위처럼 주변의 4점만 체크할 것이 아니라, 윈도를 좀 더 크게 잡고, 그 윈도 내에서 최댓값을 찾아보는 것도 한 방법이 된다.

사용자 삽입 이미지
사용자 삽입 이미지

 


/**
** http://blog.naver.com/helloktk/80051779331
*/

 

728x90

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

Bright Preserving Histogram Equalization with Maximum Entropy  (0) 2008.07.31
Adaptive Binarization  (2) 2008.07.14
Histogram Equalization  (0) 2008.06.22
FFT2D  (0) 2008.06.10
Otsu Algorithm  (6) 2008.05.30
Posted by helloktk
,