Geometry & Recognition

이미지 처리 과정에서 미분에 해당하는 그래디언트 필드(gradient field: $g_x$, $g_y$ )를 이용하면 이미지 상의 특징인 corner, edge, ridge 등의 정보를 쉽게 얻을 수 있다. 이미지의 한 지점이 이러한 특징을 가지는 특징점이 되기 위해서는 그래디언트 필드의 값이 그 점 주변에서 (3x3나 5x5정도 크기의 window) 일정한 패턴을 유지해야 한다. 이 패턴을 찾기 위해서 그래디언트 필드에 PCA를 적용해보자. 수평과 수직방향의 그래디언트 field인 $g_x$와 $g_y$ 사이의 covariance 행렬은 다음 식으로 정의된다:

$$ \Sigma = \left [ \begin {array}{cc} < g_x^2 > & < g_x g_y > \\    <g_x g_y> & <g_y^2 > \end {array}\right] =\left [ \begin {array}{cc} s_{xx} & s_{xy} \\ s_{xy} & s_{yy}\end {array}\right];$$

$<...> = \int_{W}(...) dxdy$는 픽셀 윈도에 대한 적분을 의미한다. $\Sigma$의 eigenvalue는 항상 음이 아닌 값을 갖게 되는데 (matrix 자체가 positive semi-definitive), 두 eigenvalue이 $λ_1$, $λ_2$면

$$λ_1 + λ_2 = s_{xx} + s_{yy} \ge 0, \quad \quad    λ_1  λ_2 = s_{xx} s_{yy} - s_{xy}^2 \ge0 $$

을 만족한다 (완전히 상수 이미지를 배제하면 0인 경우는 없다). eigenvalue $λ_1$, $λ_2$는 principal axis 방향으로 그래디언트 필드의 변동(분산)의 크기를 의미한다. edge나 ridge의 경우는 그 점 주변에서 잘 정의된 방향성을 가져야 하고, corner의 경우는 방향성이 없어야 한다. edge나 ridge처럼 일방향성의 그래디언트 특성을 갖거나 corner처럼 방향성이 없는 특성을 서로 구별할 수 있는 measure가 필요한데, $λ_1$과 $λ_2$를 이용하면 차원이 없는 measure을 만들 수 있다. 가장 간단한 차원이 없는 측도(dimensionless measure)는  eigenvalue의 기하평균과 산술평균의 비를 비교하는 것이다.

$$ Q = \frac { {λ_{1} λ_{2}} }{ \left( \frac {λ_{1}+λ_{2}}{2} \right)^2} = 4\frac { s_{xx} s_{yy} - s_{xy}^2}{(s_{xx} + s_{yy})^2};$$

기하평균은 산술평균보다도 항상 작으므로

$$ 0 \le Q \le 1 $$

의 범위를 갖는다. 그리고 $Q$의 complement로

$$P = 1-Q = \frac{(s_{xx}-s_{yy})^2 + 4 s_{xy}^2}{(s_{xx}+s_{yy})^2};$$를 정의할 수 있는 데 $0 \le P \le 1$이다. $Q$와 $P$의 의미는 무엇인가? 자세히 증명을 할 수 있지만 간단히 살펴보면 한 지점에서 $Q \rightarrow 1$이려면 $λ_{1} \approx λ_{2}$이어야 하고, 이는 두 주축이 동등하다는 의미이므로 그 점에서는 방향성이 없는 코너의 특성을 갖게 된다. 반대로 $Q \rightarrow 0$이면 강한 방향성을 갖게 되어서 edge(ridge) 특성을 갖게 된다.

실제적인 응용으로는 지문 인식에서 지문 영역을 알아내거나 (이 경우는 상당이 큰 윈도를 사용해야 한다) 또는 이미지 텍스쳐 특성을 파악하기 위해서는 이미지를 작은 블록으로 나누고 그 블록 내의 미분 연산자의 균일성을 파악할 필요가 있는데 이 차원이 없는 측도는 이미지의 상태에 상관없이 좋은 기준을 주게 된다.

참고 논문:

Image field categorization and edge/corner detection from gradient covariance
Ando, S.
Pattern Analysis and Machine Intelligence, IEEE Transactions on
Volume 22, Issue 2, Feb 2000 Page(s):179 - 190

** 네이버 블로그 이전;

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

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

점집합을 일반적인 2차 곡선으로 피팅하는 경우에 방정식은

$$ a x^2 + by^2 + cxy +d x + ey +f = 0$$

의 계수를 주어진 데이터를 이용하여서 구해야 한다. 실제 문제에서는 타원, 포물선 쌍곡 선등의 타입에 따라 몇 가지 제약 조건을 넣어 피팅을 한다. 원은 타원의 특별한 경우로 일반적으로 $a = b$, $c = 0$의 제약 조건이 필요하다. 그러나 보다 엄밀하게 제약을 하게 되면 $a = b = 1$의 추가 조건을 줄 수 있다. 이 경우는 점들이 모두 일직선에 있는 경우를 ($a = b = 0$) 취급할 수 없게 된다. 이 예외적인 경우를 제외하고는 최소자승법을 사용하면 계수를 매우 쉽게 구할 수 있기 때문에 많이 이용된다.

문제: 주어진 데이터를 fitting 하는 이차곡선(원)

$$x^2  + y^2 + A x + B  y + C = 0$$

의 계수 $A, B, C$를 최소자승법을 사용해서 구하라. 

주어진 점집합이 원 위의 점이면 우변이 0이 되어야 하나, 실제 데이터를 얻는 과정에서 여러 노이즈에 노출되므로 일반적으로 0이 되지 않는다. 최소자승법은 주어진 점들이 원에서 벗어나는 정도의 제곱 합이 최소가 되도록 하는 계수 $A, B, C$를 결정한다.  원과 점의 편차의 제곱합
$$ L=\sum_ i   \left |x_i^2 + y_i^2 + A x_i + B y_i + C \right|^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, \quad m_y = \frac{1}{N} \sum_i y_i$$

계에서 계산을 하자. 이를 위해 입력점의 $x$, $y$ 성분에서 각각 $m_x$, $m_y$만큼을 빼준 값을 좌표값으로 대입하면 된다: 

$$ x_i \to x_i - m_x,\quad 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_{xy^2})  $$

$$ S_{xy}  A  + S_{y^2} B = - (S_{y^3} + S_{x^2 y} )$$

을 얻을 수 있다. 이를 풀면

$$ A = \frac{- S_{y^2} ( S_{x^3} + S_{xy^2}) + S_{xy} (S_{y^3} + S_{x^2y})  }{ 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_{xy^2}) }{S_{x^2} S_{y^2}- S_{xy}^2} $$

여기서 주어진 데이터의 각 차수에 해당하는 moment는 아래처럼 계산된다:

추정된 원의 중심 $(c_x, c_y)$는 

$$ c_x = - \frac{A}{2},   \qquad 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);
}

히스토그램처럼 균일한 간격으로 주어진 데이터에서 피크 값의 위치를 찾기 위해서는 먼저 노이즈 제거를 위해 몇 번의 smoothing 과정을 적용해야 한다 (구체적인 방법은 필요에 따라 다양하다). smoothing 과정을 거친 히스토그램에서 피크 위치를 찾는 것은 쉬운 작업이다. 그런데 경우에 따라 그 위치를 sub-order까지 구해야 필요가 생긴다. 예를 들면, 실수 값 데이터 시퀀스 대한 히스토그램을 만들려면 먼저 정수로 양자화시켜야 가능하다. 이 양자화된 정보를 담고 있는 히스토그램에서 피크의 위치를 찾으려 할 때 양자화 과정에서 잃어버린 정밀도를 복원하려면 interpolation을 써야 한다. 간단하게 parabolic 근사로 피크의 위치를 찾는 경우를 생각하자. 이 방법은 수학적으로 단순하지만 영상처리 알고리즘에서 많이 이용이 되고 있다. 데이터 시퀀스가 균일한 간격에서 주어졌으므로 계산은 $-1, 0, 1$의 세 군데 위치에서 중앙 근방에 피크가 나타나는지를 고려하면 된다. 세 점에서 히스토그램 값이 각각 $h_m, h_0, h_p$일 때 이들을 지나는 이차 곡선의 피크 위치를 찾자. 주어진 이차 곡선을

$$y = a  (x - c)^2 + b$$

꼴로 쓰면 $c$가 0에서 벗어난 정도를 나타낸다.

$$(-1, h_m) \quad \rightarrow \quad h_m = a(-1-c)^2 +b; \\(0, h_0) \quad \rightarrow \quad h_0 = a c^2 +b; \\ (+1, h_p) \quad \rightarrow \quad h_m = a(+1-c)^2 +b; $$ 이므로 $$h_m - h_p = 4ac, \\ h_m + h_p -2 h_0=2 a,\\ \therefore ~c = \frac { h_m - h_p}{2(h_m + h_p -2 h_0 )}$$

아래 코드는 피크의 위치가 중앙점에서 벗어난 정도를 준다.

bool parabolicInterpolate(double hm, double h0, double hp, double *c) {
      double a = hm + hp - 2 * h0;
      if (a >= 0) return false; // not a parabola(a==0), not convex up(a>0);
      *c = 0.5 * (hm - hp) / a;
      if (*c < -0.5 || *c > 0.5) return false; //  too far;
      return true;
}