평면 위에 점집합이 주어지고 이들을 잘 기술하는 직선의 방정식을 구해야 할 경우가 많이 발생한다. 이미지의 에지 정보를 이용해 선분을 찾는 경우에 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$의 고윳값과 고유 벡터를 구하면, 큰 고윳값을 갖는 고유 벡터 방향이 선분의 방향이 될 것이다.
잘 피팅이 이루어지려면 두 고윳값의 차이가 커야 한다. 또한 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 결과다.
서로 다른 조건에서 생성이 된 두 개의 영상을 비교하기 위해서는 먼저 영상을 정규화시키어야 한다. 쉽게 생각할 수 있는 정규화의 방법은 두 개의 영상이 동일한 히스토그램을 같도록 만들면 된다. 어떻게 동일한 히스토그램을 만들 수 있는가? 간단히 생각할 수 있는 방법은 주어진 영상의 히스토그램을 전 픽셀에 걸쳐서 균일한 빈도를 보이는 히스토그램이 되도록 변환을 하는 것이다. 즉, 히스토그램을 평탄화시키는 것이다. 이 변환을 $F(x)$라고 하면, 입력 영상의 픽셀 값 $x$는 새로운 픽셀 값 $y = F(x)$로 바뀐다. 영상의 히스토그램은 픽셀 값의 빈도를 나타내므로 확률 밀도 함수로 해석할 수 있고 (픽셀 값을 연속 변수로 볼 때) 주어진 영상의 확률 밀도 함수를 $P(x)$라고 하자. 그러면 $P(x) dx$는 픽셀 값이 $[x, x+dx]$인 픽셀이 영상에서 나타날 확률을 의미한다. 변환 후에서는 임의의 픽셀 값이 균일해야 하므로 확률 밀도 함수는 상수 $a$이고, 픽셀 값의 구간이 $[0,255]$ 이므로 $a = 1/ 255$로 주어진다. 따라서, 변환이 확률을 보존하기 위해서는 $$P(x) dx = a dy$$을 만족해야 하므로 변환 후 픽셀 값 y와 변환 전 픽셀 값 x사이의 관계를 알 수 있다. $$\frac {dy}{ dx} = \frac {1}{a} P(x) $$ 이 식을 적분하면, $$y = F(x) = \frac {1}{a} \int_0^x P(x') dx' = 255\times \text {cumulative sum of } P(x)$$ 즉, 히스토그램을 평탄화시키는 변환은 원래 영상 히스토그램의 누적합을 전체 픽셀 수로 나눈 값, $$k \rightarrow F(k)=255 \times \frac {\sum_{i=0}^{k} H [i]}{ \sum_{i=0}^{255} H [i]}$$으로 결정된다. 영상의 히스토그램은 이산적이기 때문에 실제로 모든 픽셀 값에서 균일하게 나타나지 않지만, 구간 변환은 거의 일정하게 나타난다. 이 변환을 거친 영상의 누적 히스토그램(cumulative histogram)을 $y$-축 픽셀 값을 $x$-축으로 해서 그리면 거의 일직선으로 증가하는 모양을 나타낸다.
이 히스토그램의 평탄화는 픽셀 값이 좁은 영역에 몰려있는 영상을 전 영역에 분포하도록 한다. 인간은 영상의 밝기에 의해서 보다는 밝고 어두움의 변화의 크기에 의해서 인지도가 증가하게 되는데, 평탄화된 영상은 보다 쉽게 인식이 될 수 있는 구조를 가진다.
void histogrm_eq(BYTE *src, int width, int height, BYTE *dst) {
int hist[256] = {0};
int cum[256];
for (int k = width * height; k-- > 0; )
hist[src[k]]++ ;
cum[0] = hist[0];
for (int k = 1; k < 256; k++)
cum[k] = hist[k] + cum[k - 1];
double a = 255. / cum[255];
for (int k = width * height; k-- > 0; ) {
int x = int(a * cum[src[k]] + 0.5);
dst[k] = x > 255 ? 255: x;
}
}