hi

椭圆的最小二乘拟合

只需要将X换成你的边界坐标即可

x = [59 136

58 137

57 137

56 137

55 138

54 139

53 140

52 141

51 142

51 143

51 144

50 145

50 146

50 147

50 148

49 149

49 150

49 151

49 152

49 153

50 154

50 155

50 156

50 157

51 158

51 159

51 160

52 161

52 162

53 163

54 164

54 165

55 166

56 167

57 168

58 169

59 170

60 171

61 171

62 172

63 172

64 172

65 173

66 173

67 173

68 173

69 173

70 173

71 173

72 173

73 172

74 172

75 172

76 171

77 171

78 170

79 169

79 168

80 167

80 166

80 165

81 164

81 163

81 162

81 161

81 160

81 159

81 158

81 157

81 156

81 155

81 154

81 153

80 152

80 151

80 150

79 149

79 148

79 147

78 146

78 145

77 144

76 143

75 142

74 141

73 140

72 139

71 138

70 138

69 137

68 137

67 137

66 136

65 136

64 136

63 136

62 136

61 136

60 136];

%  p0=[1 1 1 1 1 1];

p0=[0.005 0.005 0.005 0.005 0.005 0.005];

warning off

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)*x(:,2)+p(6);

% 拟合系数，最小二乘方法

p=nlinfit(x,zeros(size(x,1),1),F,p0);

p(1)

p(2)

p(3)

p(4)

p(5)

p(6)

A=p(1)/p(6);

B=p(2)/p(6);

C=p(3)/p(6);

D=p(4)/p(6);

E=p(5)/p(6);

%%椭圆中心

X_center = (B*E-2*C*D)/(4*A*C - B^2);

Y_center = (B*D-2*A*E)/(4*A*C - B^2);

fprintf(' X_center=%g, Y_center=%g\n',X_center,Y_center);

%%长短轴

a= 2*sqrt((2*A*(X_center^2)+2*C*(Y_center^2)+2*B*X_center*Y_center-2)/(A+C+sqrt(((A-C)^2+B^2))));

b= 2*sqrt((2*A*(X_center^2)+2*C*(Y_center^2)+2*B*X_center*Y_center-2)/(A+C-sqrt(((A-C)^2+B^2))));

%%长轴倾角

q=0.5 * atan(B/(A-C));

fprintf(' q=%g\n',q);

fprintf(' a=%g, b=%g\n',a,b);

plot(x(:,1),x(:,2),'ro');

hold on;

xmin=min(x(:,1));

xmax=max(x(:,1));

ymin=min(x(:,2));

ymax=max(x(:,2));

% 作图

ezplot(@(x,y)F(p,[x,y]),[xmin,xmax,ymin,ymax]);

title('曲线拟合');

%legend('样本点','拟合曲线')