2011年12月9日星期五

Ellipse fitting and parameter conversion

The third assignment of Computer Vision course involves an interesting problem. Given a large number of points in Cartesian coordinate system, how to find an ellipse that best fit these points. In this paper the illustration and solution of this problem is presented.

The theory and proofs are not trivial and I can`t understand all of them. However, it`s implementation is very simple. Below is a six-line Matlab implementation of the algorithm.


function a = fit_ellipse(x, y)
    D = [x.*x x.*y y.*y x y ones(size(x))];
    S = D' * D;
    C(6, 6) = 0; C(1, 3) = -2; C(2, 2) = 1; C(3, 1) = -2;
    [gevec, geval] = eig(S, C);
    [NegR, NegC] = find(geval < 0 & ~isinf(geval));
    a = gevec(:, NegC);
end

The input of this function is two vectors containing x and y coordinates of the points, respectively. The output is the 6 parameters of implicit equation of the ellipse.

Our TA is very smart and he used a trick to simplify the problem for us to implement. So we can avoid solving the general eigenvector problem, which is not provided by OpenCV library. However, the eigen()function provided by OpenCV can only solve symmetric matrix, which is not our case. Fortunately, there is another linear algebra C++ library called Eigen, and OpenCV provided functions to convert data structures needed by both libraries back and forth. With Eigen, we can get things done, but it leaves me rather messy code. This is why I don`t post my C++ implementation here.

Having only the implicit equation is not sufficient, because we cannot get any useful information from it directly. To visualize the ellipse, we have to get semi-major axis, semi-minor axis, rotation and center from the 6 parameters. Below is the C++ code:

    // Input parameters: a, b, c, d, e, f
    phi = atan(b / (a - c)) * 0.5; // rotation
    double cos_phi = cos(phi);
    double sin_phi = sin(phi);
    double u = (2 * c * d - b * e) / (b * b - 4 * a * c);
    double v = (2 * a * e - b * d) / (b * b - 4 * a * c);
    center = cv::Vec2d(u, v);        // center

    // eliminate rotation and recalculate 6 parameters
    double aa = a * cos_phi * cos_phi - b * cos_phi * sin_phi + c * sin_phi * sin_phi;
    double bb = 0;
    double cc = a * sin_phi * sin_phi + b * cos_phi * sin_phi + c * cos_phi * cos_phi;
    double dd = d * cos_phi - e * sin_phi;
    double ee = d * sin_phi + e * cos_phi;
    double ff = 1 + (d * d) / (4 * a) + (e * e) / (4 * c);

    sx = sqrt(ff / aa);              // semi-major axis
    sy = sqrt(ff / cc);              // semi-minor axis

    if(sx < sy) {
        double temp = sx;
        sx = sy;
        sy = temp;
    }

Now we can use information to visualize the ellipse. Here is an example of my result:

The red dots are cross section points of human body and we want to fit. The ellipse shown is the result. We can see that waist, left and right biceps are quite good, but chest, hips, left and right thigh are not satisfactory. One possible explaination is that there may be many noisy points outside this visual range that affect the result.

没有评论:

发表评论