matlab axis函数_OpenCV与MATLAB中的椭圆拟合（四十七）

一、OpenCV椭圆拟合

//创建一个用于绘制图像的空白图    cv::Mat image = cv::Mat::ones(480, 640, CV_8UC3);  //设置蓝色背景  image.setTo(cv::Scalar(100, 0, 0));    //输入拟合点    std::vector<:point> points;    points.push_back(cv::Point(200, 240));  points.push_back(cv::Point(300, 400));  points.push_back(cv::Point(400, 360));  points.push_back(cv::Point(500, 300));  points.push_back(cv::Point(500, 200));  points.push_back(cv::Point(300, 150));    //将拟合点绘制到空白图上    for (int i = 0; i < points.size(); i++)  {      cv::circle(image, points[i], 5, cv::Scalar(0, 0, 255), 2, 8, 0);  }    //获取拟合椭圆的外包围矩形  cv::RotatedRect rotate_rect = cv::fitEllipse(points);  //绘制拟合椭圆  cv::ellipse(image, rotate_rect, cv::Scalar(0, 255, 255), 2, 8);    cv::imshow("image", image);  cv::waitKey(0);

opencv中的椭圆拟合函数是：

C++: fitEllipse

C  : cvFitEllipse2

函数原型：

C++: RotatedRect fitEllipse(InputArray points)

C  : CvBox2D cvFitEllipse2(const CvArr* points)

对应C++接口，fitEllipse的输入2维点集可以以std::vector<> or Mat形式存储。函数的返回是RotatedRect 类型,它具有3个成员变量：

center(块中心(x,y)), size(宽和高), angle(旋转角)，刚好与椭圆具备的5个参数[a,b,x,y,w,h]吻合，即分别对应椭圆的中心，短轴和长轴，旋转角度。 

类RotatedRect的原型声明：

class CV_EXPORTS RotatedRect{   public:  //构造函数    RotatedRect();    RotatedRect(const Point2f& center, const Size2f& size, float angle);    RotatedRect(const CvBox2D& box);    void points(Point2f pts[]) const;    //!返回矩形的4个顶点    Rect boundingRect() const;  //返回包含旋转矩形的最小矩形    operator CvBox2D() const;  //!转换到旧式的cvbox2d结构    Point2f center; //矩形的质心    Size2f size;    //矩形的边长    float angle;    //旋转角度，当角度为0、90、180、270等时，矩形就成了一个直立的矩形};

     需要注意的是angle对应角度是角度,如果使用三角函数sin,cos,tan之类的需要先转换成弧度,boundingRect对应水平方向外接矩形而不是椭圆的外接矩形,如果需要获取椭圆的两个半轴长度可以通过size.width和size.height得到,其余的参数看一下上面代码对应注释.

     下面就是怎么通过一个RotatedRect来获取椭圆参数,具体原理可以参考百度文库的一个文档:点击打开链接,直接截图看一下:

因为这里默认的是顺时针旋转为正,而图像里面默认是逆时针旋转为正,只需要把B乘以一个负号就行,感兴趣的同学自己推导一下.

好了，到了上代码的时间了，代码如下：

#include #include #include using namespace cv;using namespace std; struct EllipsePara{  Point2f c;  float A;  float B;  float C;  float F;}; void getEllipsePara(RotatedRect & ellipsemege, EllipsePara& EP_t){  float theta = ellipsemege.angle * CV_PI / 180.0 ;  float a = ellipsemege.size.width / 2.0;  float b = ellipsemege.size.height / 2.0;   EP_t.c.x = ellipsemege.center.x;  EP_t.c.y = ellipsemege.center.y;   EP_t.A = a * a * sin(theta) * sin(theta) + b * b * cos(theta) * cos(theta);  EP_t.B = (-2.0) * (a * a - b * b) * sin(theta) * cos(theta);  EP_t.C = a * a * cos(theta) * cos(theta) + b * b * sin(theta) * sin(theta);  EP_t.F = (-1.0) * a * a * b * b;  //cout << "theta: " << theta << " x: " << EP.c.x << " y: " << EP.c.y << " C: " << EP.C << " F: " << EP.F << endl;}void main(){  //RotatedRect(const Point2f& center, const Size2f& size, float angle);    RotatedRect Rec_t = RotatedRect(Point2f(0,0), Size2f(4,1), 90);  cout << Rec_t.size.width << " " << Rec_t.size.height << endl;   cout << Rec_t.boundingRect().width <<  " " << Rec_t.boundingRect().height <<  endl;  EllipsePara EP_t;  getEllipsePara(Rec_t, EP_t);   cout << EP_t.A << " " << EP_t.B << " " << EP_t.C << " " << EP_t.F << endl;}

得到的椭圆方程为：

二、MATLAB椭圆拟合

设平面任意位置椭圆方程为： 

x^2+Axy+By^2+Cx+Dy+E=0 

设P_i (x_i,y_i )(i=1,2,…,N)为椭圆轮廓上的N(N≥5) 个测量点，依据最小二乘原理，所拟合的目标函数为： 

 

欲使F为最小，需使 

 

由此可以得方程： 

 

解方程可以得到A，B，C，D，E的值。 

根据椭圆的几何知识，可以计算出椭圆的五个参数：位置参数(θ,x_0,y_0 )以及形状参数(a,b)。 

matlab程序

function [J,x0, y0, L_axis, S_axis, phi] = ellipse_fit_x(x,y)[N,N1] = size(x);x = x(:);y = y(:);% % x^2+a*x*y+b*y^2+c*x+d*y+e = 0;M = [sum(x.^2.*y.^2) sum(x.*y.^3) sum(x.^2.*y) sum(x.*y.^2) sum(x.*y)     sum(x.*y.^3) sum(y.^4) sum(x.*y.^2) sum(y.^3) sum(y.^2)     sum(x.^2.*y) sum(x.*y.^2) sum(x.^2) sum(x.*y) sum(x)     sum(x.*y.^2) sum(y.^3) sum(x.*y) sum(y.^2) sum(y)     sum(x.*y) sum(y.^2) sum(x) sum(y) N] ; B = -[sum(x.^3.*y) sum(x.^2.*y.^2) sum(x.^3) sum(x.^2.*y) sum(x.^2)]';g = M\B;[a,b,c,d,e] = deal(g(1),g(2),g(3),g(4),g(5));J = [a,b,c,d,e] delta = a^2-4*b; x0 = (2*b*c-a*d)/deltay0 = (2*d-a*c)/deltanom = 2 * (a*c*d - b*c^2 - d^2 + 4*b*e - a^2*e);s = sqrt(a^2+(1-b)^2);% s1 = nom/(delta*(b-s+1));% s2 = nom/(delta*(b+s+1));% if s1<0 ||s2<0%    return% enda_prime = sqrt(nom/(delta*(b-s+1))); b_prime = sqrt(nom/(delta*(b+s+1))); % nom = b*c^2-a*c*d +d^2 - 4*b*e+ a^2*e;% s = sqrt(a^2+(1-b)^2);% a_prime = sqrt((nom/delta)*(2*(b+s+1)/delta)); % b_prime = sqrt((nom/delta)*(2*(b-s+1)/delta)); L_axis = max(a_prime, b_prime) S_axis = min(a_prime, b_prime)% phi = atan(sqrt((a_prime^2-b_prime^2*b)/(a_prime^2*b-b_prime^2)));phi = 0.5*atan(a/(1-b))%  if (a_prime%       phi = pi/2 - phi;% end

X= [200,300,400,500,500,300]Y = [240,400,360,300,200,150]ellipse_fit_x(X,Y)

输出如下：

J =   -1.1271   -0.0522 -429.8358  470.6007 -651.5050x0 =  388.9551y0 =  308.8295L_axis =   0.0000e+00 + 1.9743e+02iS_axis =   96.4370phi =   -0.4099ans =   -1.1271   -0.0522 -429.8358  470.6007 -651.5050

三、绘制椭圆

1、

  ellipse(img,    Point( WINDOW_WIDTH/2, WINDOW_WIDTH/2 ),    Size( WINDOW_WIDTH/4, WINDOW_WIDTH/16 ),    angle,    0,    360,    Scalar( 255, 129, 0 ),    1,    lineType );

2、

  Mat img；  RotatedRect ellipsebox;  ellipse(img,   ellipsebox,    Scalar(255, 129, 0),     1,     8);

四、MATLAB最小二乘拟合

clc;%% 鼠标输入点，enter键结束axis([-10 10 -10 10]);[x,y]=ginput;   %读取坐标直到按下回车键,返回坐标点的x,y坐标num=length(x);  %计算点的个数%% 直接用最小二乘进行拟合%通过最小化误差的平方和寻找数据的最佳函数匹配[p1,s1]=polyfit(x,y,1);     %n=1为直线拟合 x,y为数据点，n为多项式阶数，返回p为幂次从高到低的多项式系数向量[p2,s2]=polyfit(x,y,num-2); %n>1为曲线拟合，找到次数为n的多项式，对于数据点集(x,y),满足差的平方和最小[p3,s3]=polyfit(x,y,num-1); %x必须是单调的。矩阵s用于在polyval中来估计误差。xcurve=-10:2:10;            %在这些点处求多项式的值p1curve=polyval(p1,xcurve); %多项式曲线求值,返回对应自变量xcurve在给定系数P的多项式的值p2curve=polyval(p2,xcurve);p3curve=polyval(p3,xcurve);%% 画图plot(xcurve,p1curve,'g-',xcurve,p2curve,'b-',...    xcurve,p3curve,'r-',x,y,'*');title('不同次数的最小二乘拟合');legend('degree 1','degree num-2','degree num-1','points');

五、基于RANSAC的直线拟合

clc;clear;%% 随机点生成g_NumOfPoints = 500;   % 点数g_ErrPointPart = 0.4;  % 噪声g_NormDistrVar = 1;    % 标准偏差% 生成随机点theta = (rand(1) + 1) * pi/6;R = ( rand([1 g_NumOfPoints]) - 0.5) * 100;DIST = randn([1 g_NumOfPoints]) * g_NormDistrVar;Data = [cos(theta); sin(theta)] * R + [-sin(theta); cos(theta)] * DIST;Data(:, 1:floor(g_ErrPointPart * g_NumOfPoints)) = 2 * [max(abs(Data(1,:))), 0; 0, max(abs(Data(2,:)))] *...                                                        (rand([2 floor(g_ErrPointPart * g_NumOfPoints)]) - 0.5);plot(Data(1, :), Data(2, :), '.', 'Tag', 'DATA');hold on;%% RANSAC拟合% 拟合模型初始化nSampLen = 2;                       %设定模型所依据的点数nIter = 50;                         %最大循环次数dThreshold = 2;                     %残差阈值nDataLen = size(Data, 2);           %数据长度RANSAC_model = NaN;                 %跳过缺失模型RANSAC_mask = zeros([1 nDataLen]);  %全0矩阵，1表示符合模型，0表示不符合nMaxInlyerCount = -1;%% 主循环for i = 1:nIter     %  抽样，选取两个不同的点    SampleMask = zeros([1 nDataLen]);      while sum( SampleMask ) ~= nSampLen        ind = ceil(nDataLen .* rand(1, nSampLen - sum(SampleMask))); %rand产生随机数，ceil向离它最近的大整数圆整        SampleMask(ind) = 1;    end        Sample = find( SampleMask );        %找出非零元素的索引值，即建立模型的点    %% 建立模型,并查找符合模型的点    ModelSet = feval(@TLS, Data(:, Sample));    %计算所有符合模型的点的最小二乘    for iModel = 1:size(ModelSet, 3)       CurModel = ModelSet(:, :, iModel);        %当前模型对应的直线参数         CurMask =( abs( CurModel * [Data; ones([1 size(Data, 2)])])< dThreshold);%到直线距离小于阈值的点符合模型,标记为1      nCurInlyerCount = sum(CurMask);           %计算符合直线模型的点的个数        %% 选取最佳模型        if nCurInlyerCount > nMaxInlyerCount    %符合模型的点数最多的模型即为最佳模型            nMaxInlyerCount = nCurInlyerCount;            RANSAC_mask = CurMask;            RANSAC_model = CurModel;        end    endend%% 画最佳模型的拟合结果MinX=min(Data(1, :));MaxX=max(Data(1, :));MinX_Y=-(RANSAC_model(1).*MinX+RANSAC_model(3))./RANSAC_model(2);MaxX_Y=-(RANSAC_model(1).*MaxX+RANSAC_model(3))./RANSAC_model(2);plot([MinX MaxX],[MinX_Y MaxX_Y],'r-');title('ransac在噪声情况下的直线拟合');

六、基于RANSAC的椭圆拟合

%%clc;clear;%%  生成 带噪声的椭圆% 参数初始化g_NumOfPoints = 500;   % 点数g_ErrPointPart = 0.4;  % 噪声g_NormDistrVar = 3;    % 标准偏差a=10;b=20;             %长轴短轴angle=60;              %倾斜角%% 椭圆生成beta = angle * (pi / 180);alpha = linspace(0, 360, g_NumOfPoints) .* (pi / 180); X = (a * cos(alpha) * cos(beta)- b * sin(alpha) * sin(beta) )+wgn(1,length(alpha),g_NormDistrVar^2,'linear');    Y = (a * cos(alpha) * sin(beta)+ b * sin(alpha) * cos(beta) )+wgn(1,length(alpha),g_NormDistrVar^2,'linear');Data=[X;Y];plot(Data(1, :), Data(2, :), '.', 'Tag', 'DATA');hold on;%% RANSAC椭圆拟合%椭圆一般方程：Ax2+Bxy+Cy2+Dx+Ey+F=0%F=@(p,x)p(1)*x(:,1).^2+p(2)*x(:,1).*x(:,2)+p(3)*x(:,2).^2+p(4)*x(:,1)+p(5)%%  参数初始化nSampLen = 3;               %设定模型所依据的点数nDataLen = size(Data, 2);   %数据长度nIter = 50;                 %最大循环次数dThreshold = 2;             %残差阈值nMaxInlyerCount=-1;         %点数下限A=zeros([2 1]);B=zeros([2 1]);P=zeros([2 1]);%%  主循环for i = 1:nIter    SampleMask = zeros([1 nDataLen]);       while sum( SampleMask ) ~= nSampLen        ind = ceil(nDataLen .* rand(1, nSampLen - sum(SampleMask))); %抽样，选取nSampLen个不同的点        SampleMask(ind) = 1;     end        Sample = find( SampleMask );                                    %找出非零元素的索引值，即建立模型的点    %%  建立模型，存储建模需要的坐标点,焦点和过椭圆的一个点      %椭圆定义方程：到两定点之间距离和为常数      A(:,1)=Data(:,ind(1));    %焦点      B(:,1)=Data(:,ind(2));    %焦点      P(:,1)=Data(:,ind(3));    %椭圆上一点      DIST=sqrt((P(1,1)-A(1,1)).^2+(P(2,1)-A(2,1)).^2)+sqrt((P(1,1)-B(1,1)).^2+(P(2,1)-B(2,1)).^2);      xx=[];       nCurInlyerCount=0;        %初始化点数为0个     %%  是否符合模型？     for k=1:g_NumOfPoints        CurModel=[A(1,1)   A(2,1)  B(1,1)  B(2,1)  DIST ];        pdist=sqrt((Data(1,k)-A(1,1)).^2+(Data(2,k)-A(2,1)).^2)+sqrt((Data(1,k)-B(1,1)).^2+(Data(2,k)-B(2,1)).^2);        CurMask =(abs(DIST-pdist)< dThreshold);     %到直线距离小于阈值的点符合模型,标记为1        nCurInlyerCount =nCurInlyerCount+CurMask;             %计算符合椭圆模型的点的个数        if(CurMask==1)             xx =[xx,Data(:,k)];        end     end       %% 选取最佳模型        if nCurInlyerCount > nMaxInlyerCount   %符合模型的点数最多的模型即为最佳模型            nMaxInlyerCount = nCurInlyerCount;            Ellipse_mask = CurMask;             Ellipse_model = CurModel;             Ellipse_points = [A B P];             Ellipse_x =xx;        end     end%% 由符合点拟合椭圆%椭圆一般方程：Ax2+Bxy+Cy2+Dx+Ey+F=0F=@(p,x)p(1)*x(:,1).^2+p(2)*x(:,1).*x(:,2)+p(3)*x(:,2).^2+p(4)*x(:,1)+p(5)*x(:,2)+p(6);p0=[1 1 1 1 1 1];x=Ellipse_x';pr=nlinfit(x,zeros(size(x,1),1),F,p0);  % 拟合系数，最小二乘方法xmin=min(x(:,1));xmax=max(x(:,1));ymin=min(x(:,2));ymax=max(x(:,2));%% 画点作图plot(Ellipse_points(1,:),Ellipse_points(2,:),'r*');hold on;plot(Ellipse_x(1,:),Ellipse_x(2,:),'yo');hold on;ezplot(@(x,y)F(pr,[x,y]),[-1+xmin,1+xmax,-1+ymin,1+ymax]);title('RANSAC椭圆拟合');legend('样本点','抽取点','符合点','拟合曲线')

最小二乘拟合&基于RANSAC的直线拟合&椭圆拟合

https://blog.csdn.net/weixin_34061555/article/details/85521018?utm_medium=distribute.pc_relevant.none-task-blog-BlogCommendFromMachineLearnPai2-2.nonecase&depth_1-utm_source=distribute.pc_relevant.none-task-blog-BlogCommendFromMachineLearnPai2-2.nonecase

matlab椭圆拟合

​https://blog.csdn.net/qq_30282427/java/article/details/82117365

获取椭圆参数

https://blog.csdn.net/sinat_31425585/article/details/75103239

拟合椭圆

​https://blog.csdn.net/i_chaoren/java/article/details/78358991

机器人实时纠偏系统(一)机器人实时纠偏(二)OpenCV+VS开发环境配置(三)基于结构光视觉的焊接机器人纠偏(四)基于结构光视觉的机器人焊接(五)基于结构光视觉的机器人焊接(六)机器人初始点导引(七)MATLAB摄像机工具箱标定相机参数(八)机器人的手眼标定(九)机器人坐标获取(十)机器人调试(十一)TCP/IP客户端API编程(十二)结构光传感器上位机界面多线程编程(十三)TCP&UDP(十四)C/C++ Programing(十五)机器人扫描与跟踪调试(十六)结构光传感器库函数(十七)结构光传感器编程(十八)C/C++ Programing(十九)C/C++ Programing(二十)结构光传感器编程(二十一)DX200操作要领(二十二)DX200操作要领(二十三)DX200独立协调—工装轴协调(二十四)DX200独立协调—无夹具协调(二十五)图像处理调试(二十六)STM32MODBUS_CRC编程(二十七)

在C++中调用Matlab函数(二十八)

机器人手眼标定MATLAB及C++实现

机器人位姿运算及Eigen的使用(三十)

OpenCV与Eigen矩阵运算(三十一)

VS中数据读写及OpenCV拟合(三十二)

VS2013配置OpenGL库(三十三)

曲线拟合/插值C++MATLAB实现(三十四)

曲线拟合绘制滤波及机器人平移(三十五)

DX200操作要领—示教1(三十六)

直接打开与平移变换(三十七)PAM与镜像平移变换(三十八)

修改与编辑程序(三十九)

YRC1000 宏程序命令(四十)

程序编辑与试运行(四十一)

程序编辑与再现(四十二)

再现(四十三)

程序管理(四十四)

便捷功能(四十五)

便捷功能(四十六)