RANSAC: Ellipse Fitting

타원은 원뿔을 평면으로 잘랐을 때 생기는 곡선 중의 하나로 다음 이차 형식(quadratic form)으로 표현된다.

$$\text{conic eq: }a{x}^2+bxy+c{y}^2+dx+ey+f=0$$

이차 형식의 계수를 구하기 위해서는 타원 위의 서로 다른 5개의 점이 필요하다 (타원이기 위해서는 $b^2-4ac<0$). 주어진 5개의 평면상의 점 $\{(x_i,y_i)|i=1,2,..,5\}$가 결정하는 타원은

$$\begin{array}{c}a{x}_1^2+b{x}_1{y}_1+c{y}_1^2+d{x}_1+e{y}_1+f=0\\ a{x}_2^2+b{x}_2{y}_2+c{y}_2^2+d{x}_2+e{y}_2+f=0\\ a{x}_3^2+b{x}_3{y}_3+c{y}_3^2+d{x}_3+e{y}_3+f=0\\ a{x}_4^2+b{x}_4{y}_4+c{y}_4^2+d{x}_4+e{y}_4+f=0\\ a{x}_5^2+b{x}_5{y}_5+c{y}_5^2+d{x}_5+e{y}_5+f=0\end{array}$$

이를 만족시키는 해는 다음 행렬식으로 쓰인 방정식의 해와 같다:

$$\left(\begin{array}{cccccc}{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{array}\right)\left(\begin{array}{c}a\\ b\\ c\\ d\\ e\\ f\end{array}\right)=0$$

또는 간단히

$$\mathbf{A}.\mathbf{x}=\mathbf{0}$$

의 형태로 쓰여진다. 6개의 conic-parameters를 결정해야 하므로, $\mathbf{A}$는 마지막 행을 0으로 채워서 $6\times 6$ 행렬로 만들어서 사용한다.  이 경우 ${\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}.\left(\begin{array}{c}{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{array}\right)$$

행렬 $\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개의 점을 선택해야 한다.

red = inliers

// 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;
}