타원은 원뿔을 평면으로 잘랐을 때 생기는 곡선 중의 하나로 다음 이차 형식(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개의 점을 선택해야 한다.
// 2024.05.30에 새롭게 작성;
double alg_distance(const CfPt& pts, double conic_param[6]) {
return conic_param[0] * pts.x * pts.x +
conic_param[1] * pts.x * pts.y +
conic_param[2] * pts.y * pts.y +
conic_param[3] * pts.x +
conic_param[4] * pts.y +
conic_param[5];
}
int num_sampling5(double prob_fail, double inlier_ratio) {
return int(log(prob_fail)/log(1-pow(inlier_ratio, 5)));
}
std::vector<int> Ransac_EllipseFit(std::vector<CfPt>& points, double ellipse_param[5]) {
if (points.size() < 5) return std::vector<int> (); //return null_vector;
CfPt center; double inv_scale;
// normalize input points for the sake of numerical stability;
std::vector<CfPt> nor_pts = normalize(points, inv_scale, center);
// RANSAC;
const double prob_fail = 0.01;
int sample_num = 1000; //number of sample
int ransac_count = 0;
const double distance_thresh = sqrt(3.84) * inv_scale;
double best_eparam[5] = {0};
std::vector<int> best_inliers;
while (sample_num > ransac_count) {
// pick random 5 indices:[0,points.size()-1];
int quintet[5];
random_quintet(points.size()-1, quintet);
CfPt selected[5];
for (int i = 0; i < 5; i++)
selected[i] = nor_pts[quintet[i]];
// ellipse parameters with 5 points;
double conic_param[6];
solve_conic(selected, conic_param);
// find inliers;
std::vector<int> inliers;
inliers.reserve(points.size());
for (int i = nor_pts.size(); i-->0;) {
// error measure = algebric distance;
double deviation = alg_distance(nor_pts[i], conic_param);
if (fabs(deviation) < distance_thresh)
inliers.push_back(i);
}
if ((inliers.size() > best_inliers.size()) &&
solve_ellipse(conic_param, ellipse_param)) {
// update sampling_num;
sample_num = num_sampling5(prob_fail, double(inliers.size())/points.size());
// update best_inliers;
best_inliers.swap(inliers);
// update best ellipse param;
for (int i = 0; i < 5; i++)
best_eparam[i] = ellipse_param[i];
}
if (++ransac_count > 1500) {
TRACE("error! ransac_count exceed!\n");
break;
}
}
// recover original scale and coordinate;
denormalize(best_eparam, best_eparam, inv_scale, center);
if (best_eparam[0] > 0 && best_eparam[1] > 0) {
for (int i = 0; i < 5; i++)
ellipse_param[i] = best_eparam[i];
TRACE("ellipse found(%d, %d)\n", sample_num, ransac_count);
} else
best_inliers.clear();
return best_inliers;
}'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 |
