Question

matlab拟合离散点曲线

有X=（ 0.48 ， 1.58 ，2.428，3.032，3.36，3.368，2.92，2.172 ，1.368 ，0.512 ，-0.436 ，0.876，-0.48，-1.58，-2.428，-3.032，-3.36，-3.368，-2.92，-2.172，-1.368，-0.512，0.436，0.48 ） Y=（0.752，1.8572，2.39，2.24，1.766，0.8828，0.0572... 有X=（ 0.48 ， 1.58 ，2.428，3.032，3.36，3.368，2.92，2.172 ，1.368 ，0.512 ，-0.436 ，0.876，-0.48，-1.58，-2.428，-3.032，-3.36，-3.368，-2.92，-2.172，-1.368，-0.512，0.436，0.48 ） Y=（0.752，1.8572，2.39，2.24，1.766，0.8828，0.0572，-0.3136，-0.472，-0.6064，-0.7864，-0.2804，-0.752，-1.8572 ，-2.39，-2.24，-1.766，-0.8828，-0.0572 ，0.3136，0.472，0.6064，0.7864，0.752 ）两组数据，画出组成的曲线，怎样用MATLAB将这条曲线拟合成椭圆方程，然后求出这个椭圆的面积，得到这个椭圆的主轴和X轴的夹角。并且得到这个椭圆的X最大值对应的X,Y坐标。求各位大神帮忙啊

Answer 1

问题分析
1、从数据看，与椭圆的形状吻合不是很好，尤其是短轴部分明显凹入。
2、拟合的基本思路是，使用长半轴a、短半轴b和倾角Q作为参数描述椭圆方程（进一步还可以考虑圆心不在原点，再增加两个参数），然后通过优化，目标是使得所有已知点到椭圆的总距离最短。
3、关于点到椭圆的最短距离，有很多讨论，都非常复杂，并不实用。这里采用在椭圆上取足够多的点然后直接计算的距离的方式。
4、拟合得到椭圆参数后，面积S=πab直接计算，主轴和X轴的夹角为拟合得到的数据，而关于椭圆的X最大值对应的X,Y坐标的计算，可以把点取密集一些，直接求最大值，也可以对X求导数，然后解X'=0得到对应的t，再计算X、Y。
 
参考代码
12345678910111213141516171819202122232425262728293031323334353637383940function zdX=[0.48, 1.58,2.428,3.032,3.36,3.368,2.92,2.172,1.368,0.512,-0.436,0.876,-0.48,-1.58,-2.428,-3.032,-3.36,-3.368,-2.92,-2.172,-1.368,-0.512,0.436,0.48];Y=[0.752,1.8572,2.39,2.24,1.766,0.8828,0.0572,-0.3136,-0.472,-0.6064,-0.7864,-0.2804,-0.752,-1.8572,-2.39,-2.24,-1.766,-0.8828,-0.0572,0.3136,0.472,0.6064,0.7864,0.752];clfplot(X,Y,'x')hold onc = lsqcurvefit(@curve,[3.5 1.5 30*pi/180 0 0],[X;Y],X*0)a = c(1);b = c(2);Q = c(3);T=linspace(0,2*pi,100);x1=a*cos(T);y1=b*sin(T);x=x1*cos(Q)-y1*sin(Q);y=x1*sin(Q)+y1*cos(Q);plot(x,y,'r--')% 椭圆的面积S = pi*a*b% 椭圆的X最大值对应的X,Y坐标t=fzero(@(t)-a*cos(Q)*sin(t)-b*sin(Q)*cos(t),0);x=a*cos(t)*cos(Q)-b*sin(t)*sin(Q)y=a*cos(t)*sin(Q)+b*sin(t)*cos(Q)plot(x,y,'mo')legend('原始数据','拟合曲线','X最大值',4)% 椭圆的主轴和X轴的夹角Q = Q * 180/pi function f = curve(c,xy)a = c(1);b = c(2);Q = c(3);t=linspace(0,2*pi,500);x1=a*cos(t);y1=b*sin(t);x=x1*cos(Q)-y1*sin(Q) + c(4);y=x1*sin(Q)+y1*cos(Q) + c(5);f = xy(1,:)*0;for i=1:length(xy)    f(i) = <a
href="http://www.baidu.com/s?wd=min&tn=44039180_cpr&fenlei=mv6quAkxTZn0IZRqIHckPjm4nH00T1YLrj-WnvP-njuBnARsn1m10ZwV5Hcvrjm3rH6sPfKWUMw85HfYnjn4nH6sgvPsT6KdThsqpZwYTjCEQLGCpyw9Uz4Bmy-bIi4WUvYETgN-TLwGUv3EnHTdrjR3njTsPWfYn1n1nHmY"
target="_blank"
class="baidu-highlight">min</a>( sqrt((xy(1,i)-x).^2 + (xy(2,i)-x).^2 ) );end
 
结果
程序输出：
1234567891011121314c =    3.4794    1.4809    0.5155   -0.0066         0 S =   16.1875 x =    3.1141 y =    1.3654 Q =   29.5357
也就是说a=3.48，b=1.48，Q=29.5度，S=16.2，如图所示。