2015-09-26 58 views
2

我想从一张图片中提取红色的球并在图片中获取检测到的椭圆矩阵。如何从OpenCV的fitEllipse函数获取椭圆系数?

这是我的例子: enter image description here

I阈图片,发现的红球通过使用findContour()函数的轮廓,并使用fitEllipse()以适合的椭圆。

但我想要的是得到这个椭圆的系数。因为fitEllipse()返回一个旋转矩形(RotatedRect),所以我需要重写这个函数。

一个椭圆可以表示为Ax^2 + By^2 + Cxy + Dx + Ey + F = 0;所以如果F是1(构造一个椭圆矩阵),我想得到u =(A,B,C,D,E,F)或u =(A,B,C,D,E)。

我读了fitEllipse()的源代码,总共有三个SVD过程,我想我可以从这三个SVD过程的结果中得到上述系数。但是我很困惑每个SVD过程的每个结果(变量cv :: Mat x)代表什么以及为什么这里有三个SVD?

这里是这样的功能:

cv::RotatedRect cv::fitEllipse(InputArray _points) 
{ 
    Mat points = _points.getMat(); 
    int i, n = points.checkVector(2); 
    int depth = points.depth(); 
    CV_Assert(n >= 0 && (depth == CV_32F || depth == CV_32S)); 

    RotatedRect box; 

    if(n < 5) 
     CV_Error(CV_StsBadSize, "There should be at least 5 points to fit the ellipse"); 

    // New fitellipse algorithm, contributed by Dr. Daniel Weiss 
    Point2f c(0,0); 
    double gfp[5], rp[5], t; 
    const double min_eps = 1e-8; 
    bool is_float = depth == CV_32F; 
    const Point* ptsi = points.ptr<Point>(); 
    const Point2f* ptsf = points.ptr<Point2f>(); 

    AutoBuffer<double> _Ad(n*5), _bd(n); 
    double *Ad = _Ad, *bd = _bd; 

    // first fit for parameters A - E 
    Mat A(n, 5, CV_64F, Ad); 
    Mat b(n, 1, CV_64F, bd); 
    Mat x(5, 1, CV_64F, gfp); 

    for(i = 0; i < n; i++) 
    { 
     Point2f p = is_float ? ptsf[i] : Point2f((float)ptsi[i].x, (float)ptsi[i].y); 
     c += p; 
    } 
    c.x /= n; 
    c.y /= n; 

    for(i = 0; i < n; i++) 
    { 
     Point2f p = is_float ? ptsf[i] : Point2f((float)ptsi[i].x, (float)ptsi[i].y); 
     p -= c; 

     bd[i] = 10000.0; // 1.0? 
     Ad[i*5] = -(double)p.x * p.x; // A - C signs inverted as proposed by APP 
     Ad[i*5 + 1] = -(double)p.y * p.y; 
     Ad[i*5 + 2] = -(double)p.x * p.y; 
     Ad[i*5 + 3] = p.x; 
     Ad[i*5 + 4] = p.y; 
    } 

    solve(A, b, x, DECOMP_SVD); 

    // now use general-form parameters A - E to find the ellipse center: 
    // differentiate general form wrt x/y to get two equations for cx and cy 
    A = Mat(2, 2, CV_64F, Ad); 
    b = Mat(2, 1, CV_64F, bd); 
    x = Mat(2, 1, CV_64F, rp); 
    Ad[0] = 2 * gfp[0]; 
    Ad[1] = Ad[2] = gfp[2]; 
    Ad[3] = 2 * gfp[1]; 
    bd[0] = gfp[3]; 
    bd[1] = gfp[4]; 
    solve(A, b, x, DECOMP_SVD); 

    // re-fit for parameters A - C with those center coordinates 
    A = Mat(n, 3, CV_64F, Ad); 
    b = Mat(n, 1, CV_64F, bd); 
    x = Mat(3, 1, CV_64F, gfp); 
    for(i = 0; i < n; i++) 
    { 
     Point2f p = is_float ? ptsf[i] : Point2f((float)ptsi[i].x, (float)ptsi[i].y); 
     p -= c; 
     bd[i] = 1.0; 
     Ad[i * 3] = (p.x - rp[0]) * (p.x - rp[0]); 
     Ad[i * 3 + 1] = (p.y - rp[1]) * (p.y - rp[1]); 
     Ad[i * 3 + 2] = (p.x - rp[0]) * (p.y - rp[1]); 
    } 
    solve(A, b, x, DECOMP_SVD); 

    // store angle and radii 
    rp[4] = -0.5 * atan2(gfp[2], gfp[1] - gfp[0]); // convert from APP angle usage 
    if(fabs(gfp[2]) > min_eps) 
     t = gfp[2]/sin(-2.0 * rp[4]); 
    else // ellipse is rotated by an integer multiple of pi/2 
     t = gfp[1] - gfp[0]; 
    rp[2] = fabs(gfp[0] + gfp[1] - t); 
    if(rp[2] > min_eps) 
     rp[2] = std::sqrt(2.0/rp[2]); 
    rp[3] = fabs(gfp[0] + gfp[1] + t); 
    if(rp[3] > min_eps) 
     rp[3] = std::sqrt(2.0/rp[3]); 

    box.center.x = (float)rp[0] + c.x; 
    box.center.y = (float)rp[1] + c.y; 
    box.size.width = (float)(rp[2]*2); 
    box.size.height = (float)(rp[3]*2); 
    if(box.size.width > box.size.height) 
    { 
     float tmp; 
     CV_SWAP(box.size.width, box.size.height, tmp); 
     box.angle = (float)(90 + rp[4]*180/CV_PI); 
    } 
    if(box.angle < -180) 
     box.angle += 360; 
    if(box.angle > 360) 
     box.angle -= 360; 

    return box; 
} 

的源代码链接:https://github.com/Itseez/opencv/blob/master/modules/imgproc/src/shapedescr.cpp

回答

8

fitEllipse返回RotatedRect包含椭圆的所有参数的函数。

  • XC:x中的中心
  • YC的坐标:中心的y坐标
  • 一个:主要半轴

    的椭圆是由5个参数定义

  • b:小半轴
  • THETA:旋转角度

可以获取这些参数,如:

RotatedRect e = fitEllipse(points); 

float xc = e.center.x; 
float yc = e.center.y; 
float a  = e.size.width/2; // width >= height 
float b  = e.size.height/2; 
float theta = e.angle;    // in degrees 

可以绘制椭圆使用所述RotatedRect功能ellipse

ellipse(image, e, Scalar(0,255,0)); 

,或者等效地使用椭圆参数:

ellipse(res, Point(xc, yc), Size(a, b), theta, 0.0, 360.0, Scalar(0,255,0)); 

如果你需要的隐式方程的系数的值,你可以这样做(从Wikipedia):

enter image description here

所以,你可以让你从RotatedRect需要的参数,您无需更改功能fitEllipsesolve函数用于求解线性系统或最小二乘问题。使用SVD分解方法,系统可以被过度定义和/或矩阵src1可以是单数的。

有关该算法的更多详细信息,可以看到提出了这种拟合椭圆方法的paper of Fitzgibbon