728x90

원뿔을 평면으로 잘랐을 때 나타나는 곡선인 conic section은 직교 좌표계에서 $(x, y)$에 대한 2차 형식으로 쓰인다.

이 conic section이 타원을 기술할 때 parameter {$a, b, c, d, e, f$}를 이용해서 타원의 중심, 장축과 단축의 길이, 그리고 회전각에 대한 공식을 구해보자. 2차 항을 행렬을 써서 표현하면

따라서, 적절한 회전 변환을 하면 두 좌표의 곱으로 주어지는 $xy$-항을 없앨 수 있다. 회전 변환을 시행하면 행렬은 eigenvalue를 성분으로 하는 대각 행렬이 된다. 회전 변환이 determinant를 보존하므로 determinant는 행렬의 두 eigenvalue의 곱으로 주어짐을 알 수 있다.

회전 후의 방정식이 타원의 방정식(원점이 이동된)을 기술하기 위해서는 $det>0$ 이어야 한다 (회전시킨 후 식에서 $x^2$과 $y^2$의 계수는 두 eigenvalue로 주어지므로 같은 부호를 가져야 한다.)

conic section $F$가 타원을 기술한다면, 다음과 같이 평행 이동 및 회전 변환을

$F$에 적용해서 $xy$-항이 없어지는 조건과 1차 항이 사라지도록 하는 조건을 찾으면 타원의 중심 $(x_0, y_0)$와 회전각 $\theta$가 아래처럼 주어짐을 알 수 있다:

  

그리고, $x^2$와 $y^2$의 계수는 두 eigenvalue를 주어진다.

이 eigenvalue는 위에서 구한 회전각 $\theta$를 이용해서 표현하면,

나머지 상수항은

로 주어짐을 알 수 있다. 따라서 타원의 두 축의 길이는

로 주어진다.

// conic_params(a, b, c, d, e, f);
// ellipse_params(major_axis, minor_axis, center_x, center_y, tilt_angle);
bool conic_to_ellipse(double conic_param[6], double ellipse_param[5]) {
    const double a = conic_param[0];
    const double b = conic_param[1];
    const double c = conic_param[2];
    const double d = conic_param[3];
    const double e = conic_param[4];
    const double f = conic_param[5];
    //get ellipse orientation
    const double theta = 0.5 * atan2(b, a - c);
    //get scaled major/minor axes
    const double ct = cos(theta);
    const double st = sin(theta);
    const double ap = a * ct * ct + b * ct * st + c * st * st;
    const double cp = a * st * st - b * ct * st + c * ct * ct;
    //get translations
    const double cx = (2 * c * d - b * e) / (b * b - 4 * a * c);
    const double cy = (2 * a * e - b * d) / (b * b - 4 * a * c);
    //get scale factor
    const double val = a * cx * cx + b * cx * cy + c * cy * cy;
    const double scale_inv = val - f;
    if (scale_inv / ap <= 0 || scale_inv / cp <= 0) {
        //TRACE("Error! ellipse parameters are imaginary\n");
        return 0;
    }
    ellipse_param[0] = sqrt(scale_inv / ap);
    ellipse_param[1] = sqrt(scale_inv / cp);
    ellipse_param[2] = cx;
    ellipse_param[3] = cy;
    ellipse_param[4] = theta;
    return 1;
};

Posted by helloktk

댓글을 달아 주세요

728x90

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

.

이 이차 형식의 계수를 구하기 위해서는 타원 위의 5개의 서로 다른 점이 필요하다 $(a c > 0)$. 그러나 주어진 점 정보를 얻는 과정에서 노이즈가 포함되어 conic eq에 대입했을 때 우변이 0이 안될 수 있다. 이런 경우를 고려해서, 5개 점 $\{(x_i,y_i)|i=1,...,5\}$에 대해서 conic eq의 좌변 (대수적인 거리^2) 최소가 되도록 하는 conic parameter 구하도록 하자:

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

또는 간단히

의 형태로 쓰여진다. 6개의 conic-parameters를 결정해야 하므로, A는 마지막 행을 0으로 채워서 6x6 행렬로 만들어서 사용한다.  이 경우 $\bf A$(->6x6)는 determinatn가 0인 singular 행렬이 된다. 근사해를 구하기 위해 $\bf A$를 singular value decomposition(SVD)해서 아래처럼 쓰면

가장 작은 eigenvalue$(w)$에 해당하는 $\bf V$의 열벡터가 에러를 최소로 하는 conic parameter를 준다 (행렬 $\bf V$의 각 열들이 6차원 벡터 공간의 기저를 형성한다. $\bf U$는 orthogonal 행렬이어서 오른쪽에서 곱해지는 벡터의 크기를 변화시키지 않고 단지 회전만 시킨다. 따라서 가장 작은 고윳값에 해당하는 $\bf V$의 열벡터가 $|\bf A \cdot x|$를  최소화시킨다)

 

// 주어진 점에서 타원까지의 대수적 거리: conic-eq의 좌변의 값;

static double fitting_error(double conic_params[6], double x, double y);

더보기
static double fitting_error(double conic_params[6], double x, double 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) {
    const int npts = points.size(); 
    if (npts < 5) {
        return 0;
    }
    // 안정적인 conic fitting을 위해서 데이터를 (-1,1)구간으로 normalize함;
    std::vector<CfPt> nor_pts(npts);
    double ellipse_scale; //원점에서  타원점들까지의 평균거리;
    CfPt ellipse_center;  //타원의 중심 예측;
    normalize_points(points, nor_pts, &ellipse_scale, &ellipse_center);
    //
    std::vector<int> inliers_index(npts);
    best_index.resize(npts);
    // 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((npts - 1), rand_idx);
  //5 점을 이용해서 conic parameter estimate;
 double conic_params[6];
 get_conic_params(nor_pts, rand_idx, conic_params);
        //conic parameter를 이용해서 inlier후보들을 모음;
        int ninliers = 0;
        for (int i = 0; i < npts; i++) {
            double error = fitting_error(conic_params, nor_pts[i].x, nor_pts[i].y);
            if (fabs(error) < dist_threshold) {
                inliers_index[ninliers] = i;
                ninliers++;
            }
        }
        if (ninliers > 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 = ninliers;
                    //이 단계에서 inliers를 이용해서 다시 타원을 예측한다;
                    //if (ninliers >= 6) { 
                    // std::vector<double> pts(2 * ninliers);
                    // for (int i = 0; i < ninliers; i++) {
                    // pts[2 * i + 0] = points[inliers_index[i]].x ;
                    // pts[2 * i + 1] = points[inliers_index[i]].y ;
                    // }
                    // EllipseFit(&pts[0], ninliers, conic_params);
                    // for (int i = 0; i < 5; i++) best_params[i] = conic_params[i];
                    //}
                    // if (ninliers == 5) 
                    for (int i = 0; i < 5; i++) best_params[i] = ellipse_params[i];
                    for (int i = 0; i < ninliers; i++) best_index[i] = inliers_index[i];
                    sample_num = ransac_sample_number(npts, ninliers, 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();
}

Posted by helloktk

댓글을 달아 주세요

  1. 오현경 2019.05.05 22:59  댓글주소  수정/삭제  댓글쓰기

    죄송하지만 전체 소스를 볼 수 있을까요??