타원은 원뿔을 평면으로 잘랐을 때 생기는 곡선 중의 하나로 다음 이차 형식(quadratic form)으로 표현된다.
$$\text{conic eq: } ax^2 + bxy + cy^2 +d x+ ey + f=0$$.
이차 형식의 계수를 구하기 위해서는 타원 위의 서로 다른 5개의 점이 필요하다 (타원이기 위해서는 $b^2 - 4a c < 0$). 주어진 5개의 평면상의 점 $\{(x_i, y_i)|i=1,2,..,5\}$가 결정하는 타원은
\begin{gather} ax_1^2 + bx_1 y_1 + c y_1^2 + d x_1 +e y_1 + f =0 \\ ax_2^2 + bx_2 y_2 + c y_2^2 + d x_2 +e y_2 + f =0 \\ ax_3^2 + bx_3 y_3 + c y_3^2 + d x_3 +e y_3 + f =0 \\ ax_4^2 + bx_4 y_4 + c y_4^2 + d x_4 +e y_4 + f =0 \\ ax_5^2 + bx_5 y_5 + c y_5^2 + d x_5 +e y_5 + f =0 \end{gather}
이를 만족시키는 해는 다음 행렬식으로 쓰인 방정식의 해와 같다:
$$ \begin{pmatrix}x_1^ 2& x_1 y_1 & y_1 ^2 & x_1 & y_1 & 1\\x_2^ 2& x_2 y_2 & y_2 ^2 & x_2 & y_2 & 1\\ x_3^ 2& x_3 y_3 & y_3 ^2 & x_3 & y_3 & 1\\ x_4^ 2& x_4 y_4 & y_4 ^2 & x_4 & y_4 & 1\\ x_5^ 2& x_5 y_5 & y_5 ^2 & x_5 & y_5 & 1 \end{pmatrix} \begin{pmatrix} a\\b\\c\\d\\e\\f \end{pmatrix}=0 $$
또는 간단히
$$ \mathbf{A}.\mathbf{x} =\mathbf{0}$$
의 형태로 쓰여진다. 6개의 conic-parameters를 결정해야 하므로, $\mathbf{A}$는 마지막 행을 0으로 채워서 $6\times6$ 행렬로 만들어서 사용한다. 이 경우 $\mathbf{A}_{6\times 6}$는 determinant가 0인 singular 행렬이 되고 적어도 null 공간의 차원은 1 이상이다(1인 경우만 의미있는 해를 준다). $\mathbf{A}$의 singular value decomposition(SVD)을 이용하면위 방정식의 nontrivial한 해를 얻을 수 있다. SVD에 의해 $\mathbf{A}$는 다음처럼 쓸 수 있다( singular value 중 1개가 0인 경우).
$$\mathbf{A}= \mathbf{U}. \text{diag}(w_0,w_1, w_2, w_3, w_4, 0).\mathbf{V}^T= \mathbf{U}.\begin{pmatrix} w_0 \mathbf{v}_0^T \\ w_1 \mathbf{v}_1^T\\ w_2 \mathbf{v}_2^T\\ w_3 \mathbf{v}_3^T\\ w_4 \mathbf{v}_4^T\\ 0\end{pmatrix}. $$
행렬 $\mathbf{V}$의 각 열벡터 $\{\mathbf{v}_i\}$는 6차원 벡터 공간의 기저를 형성한다. 그리고 $\mathbf{U}$는 orthogonal 행렬이므로 오른쪽에서 곱해지는 벡터의 크기를 변화시키지 않고 단지 회전만 시킨다. singular value가 $0$에 해당하는 $\mathbf{V}$의 열벡터 $\mathbf{v}_5$를 대입하면
$$ \mathbf{A} \cdot \mathbf{v}_5 = \mathbf{U}. (0,0,0,0,0)^T = 0$$
이어서 $\mathbf{A}.\mathbf{x}=0$의 nontrivial solution이 됨을 알 수 있다. 물론 singular value가 0인 경우가 여럿이 생기면 해당하는 열벡터의 선형결합도 해가 되므로 원하는 답이 아니다. 이 경우는 주어진 5개의 점들의 일부가 일직선상에 있는 경우에 해당한다(또는 겹치는 경우). 이때는 다시 5개의 점을 선택해야 한다.
// 주어진 점에서 타원까지의 대수적 거리: conic-eq의 좌변의 값;
static double fitting_error(double conic_params[6], double x, double y);
static double fitting_error(CfPt q, double conic_params[6]) {
double x = q.x, y = q.y;
return conic_params[0] * x * x +
conic_params[1] * x * y +
conic_params[2] * y * y +
conic_params[3] * x +
conic_params[4] * y +
conic_params[5];
}//
static int ransac_sample_number(int ndata, int ngood, double prob, int ntry);
static int ransac_sample_number(int ndata, int ngood, double prob, int ntry) {
if (ndata == ngood) return 1;
return (int)(log(1. - prob)/log(1. - pow((double)ngood/ndata, ntry)));
}int ransac_ellipse_fitter(std::vector<CPoint>& points, int width, int height,
double ellipse_params[5],/*(a, b, cx, cy, theta)*/
std::vector<int> &best_index) {
if (points.size() < 5) return 0;
// 안정적인 conic fitting을 위해서 데이터를 (-1,1)구간으로 normalize함;
std::vector<CfPt> nor_pts(points.size());
double ellipse_scale; //원점에서 타원점들까지의 평균거리;
CfPt ellipse_center; //타원의 중심 예측;
normalize_points(points, nor_pts, &ellipse_scale, &ellipse_center);
//
best_index.resize(nor_pts.size());
// inverse_scale을 기준으로 만들어졌음;
const double dist_threshold = sqrt(3.84) / ellipse_scale;
double best_params[5];
#define MAX_ITERATION (1500)
int sample_num = 1000; //number of sample
int ransac_iter = 0;
int max_inliers = 0;
while (sample_num > ransac_iter) {
int rand_idx[5];
// 5개 점을 랜덤하게 선택;
random_quintet((nor_pts.size() - 1), rand_idx);
// 5개 점을 이용해서 conic parameter estimate;
double conic_params[6];
get_conic_params(nor_pts, rand_idx, conic_params);
// conic parameter를 이용해서 inlier 후보들을 모음;
std::vector<int> inliers; // inlier index을 저장;
for (int i = nor_pts.size(); i-->0;) {
double error = fitting_error(nor_pts[i], conic_params);
if (fabs(error) < dist_threshold) inliers.push_back(i);
}
if (inliers.size() > max_inliers) {
if (conic_to_ellipse(conic_params, ellipse_params)) {
get_real_params(ellipse_params, ellipse_scale, ellipse_center, ellipse_params);
// 주어진 영역안의 정상적인 타원만 선택함;
if (regular_ellipse(ellipse_params, width, height, 0.5, 2.0)) {
max_inliers = inliers.size();
//이 단계에서 inliers를 이용해서 다시 타원을 예측한다;
//if (inliers.size() >= 6) {
// std::vector<CfPt> pts(inliers.size());
// for (int i = inliers.size(); i-->0;) pts[i] = nor_pts[inliers[i]];
// EllipseFit(pts, conic_params);
// conic_to_ellipse(conic_params, ellipse_params);
// get_real_params(ellipse_params, ellipse_scale, ellipse_center, ellipse_params);
//}
for (int i = 0; i < 5; i++) best_params[i] = ellipse_params[i];
for (int i = inliers.size(); i-->0;) best_index[i] = inliers[i];
sample_num = ransac_sample_number(nor_pts.size(), inliers.size(), 0.99, 5);
}
}
}
ransac_iter++;
if (ransac_iter > MAX_ITERATION) {
TRACE("Error! ransac_iter exceeds!\n");
max_inliers = 0;
break ;
}
}
if (best_params[0] > 0 && best_params[1] > 0)
for (int i = 0; i < 5; i++) ellipse_params[i] = best_params[i];
else
max_inliers = 0;
best_index.resize(max_inliers);
return best_index.size();
}
'Image Recognition' 카테고리의 다른 글
| Expectation Maximization Algorithm for Two-Component Gaussian Mixture (0) | 2017.01.02 |
|---|---|
| Union-Find Connected Component Labeling (0) | 2012.11.01 |
| Autofocus Algorithm (0) | 2012.06.03 |
| Statistical Region Merging (2) | 2012.03.25 |
| Local Histogram Equalization (0) | 2012.03.10 |
