根据已知坐标点画椭圆曲线

最新推荐文章于 2024-05-19 23:31:32 发布

Arvin457

最新推荐文章于 2024-05-19 23:31:32 发布

阅读量105

点赞数

文章标签： python python3.11

本文链接：https://blog.csdn.net/weixin_66361419/article/details/132911765

版权

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize

# 输入坐标点数据
x_data = np.array([3.941480293, 3.475033914, 3.073941549, 2.805914826, 2.734021537,
                   2.942555288, 3.354419877, 3.937493254, 4.441122968, 4.998814452,
                   5.490826688, 5.911054241, 6.336517315, 6.774639331, 7.146489382,
                   7.42027032, 7.434313787, 7.155838222, 6.594675284, 6.009223856,
                   5.487949163, 4.960735454, 4.425830988])
y_data = np.array([4.074905924, 3.591903876, 3.042621541, 2.485151431, 1.78094537,
                   1.409438513, 1.21628921, 1.218700552, 1.347703702, 1.587670099,
                   1.875068964, 2.178580205, 2.549639346, 3.023493594, 3.555617652,
                   4.194067893, 4.810987605, 5.191017166, 5.342526898, 5.267524226,
                   5.089825813, 4.825749131, 4.474423764])

# 定义椭圆参数化函数
def custom_ellipse(params, x):
    h, k, a, b, theta = params
    angle = np.radians(theta)
    x_rotated = (x - h) * np.cos(angle) + (y_data - k) * np.sin(angle)
    y_rotated = (y_data - k) * np.cos(angle) - (x - h) * np.sin(angle)
    return ((x_rotated / a) ** 2 + (y_rotated / b) ** 2 - 1)

# 初始参数估计值
params_initial = (5, 3, 4, 2, -45)

# 优化椭圆参数以最小化残差
result = minimize(lambda params: np.sum(custom_ellipse(params, x_data) ** 2), params_initial, method='Nelder-Mead')

# 获取最优参数
params_optimal = result.x

# 解析最优参数
h, k, a, b, theta = params_optimal

# 计算椭圆的一般多项式方程的系数
A = (np.cos(np.radians(theta)) ** 2) / (a ** 2) + (np.sin(np.radians(theta)) ** 2) / (b ** 2)
B = 2 * ((np.sin(np.radians(theta)) * np.cos(np.radians(theta))) / (a ** 2) - (np.sin(np.radians(theta)) * np.cos(np.radians(theta))) / (b ** 2))
C = (np.sin(np.radians(theta)) ** 2) / (a ** 2) + (np.cos(np.radians(theta)) ** 2) / (b ** 2)
D = -2 * A * h - B * k
E = -B * h - 2 * C * k
F = A * (h ** 2) + B * h * k + C * (k ** 2) - 1

# 输出椭圆的一般多项式方程
print("椭圆的一般多项式方程：")
print(f"{A:.6f}x^2 + {B:.6f}xy + {C:.6f}y^2 + {D:.6f}x + {E:.6f}y + {F:.6f} = 0")

# 绘制描点连线法得到的椭圆
t = np.linspace(0, 2 * np.pi, 100)
x_fit = h + a * np.cos(t) * np.cos(np.radians(theta)) - b * np.sin(t) * np.sin(np.radians(theta))
y_fit = k + a * np.cos(t) * np.sin(np.radians(theta)) + b * np.sin(t) * np.cos(np.radians(theta))

# 绘制原始数据和拟合结果
plt.scatter(x_data, y_data, label='Data Points', color='blue')
plt.plot(x_fit, y_fit, label='Fitted Ellipse', color='red')
plt.xlabel('X')
plt.ylabel('Y')
plt.legend()
plt.gca().set_aspect('equal', adjustable='box')  # 设置坐标轴比例为等比例

# 输出拟合结果
print("Center (h, k):", h, k)
print("Semi-major axis (a):", a)
print("Semi-minor axis (b):", b)
print("Rotation angle (theta):", theta)

# 显示拟合结果
plt.show()