用通用椭圆方程(包括倾斜和平移)最小二乘拟合椭圆。

huangapple go评论65阅读模式
英文:

Least-Squares Fitting of an Ellipse with General Ellipse Equation (including tilt and translation)

问题

我想使用以下代码将一些二维数据拟合到椭圆上,该代码来自此处 的帖子。

import numpy as np
import matplotlib.pyplot as plt
alpha = 5
beta = 3
N = 500
DIM = 2

np.random.seed(2)

# 通过采样均匀角度在单位圆上生成随机点
theta = np.random.uniform(0, 2*np.pi, (N,1))
eps_noise = 0.2 * np.random.normal(size=[N,1])
circle = np.hstack([np.cos(theta), np.sin(theta)])

# 将圆进行拉伸和旋转,形成具有随机线性变换的椭圆
B = np.random.randint(-3, 3, (DIM, DIM))
noisy_ellipse = circle.dot(B) + eps_noise

# 提取椭圆的 x 坐标和 y 坐标作为列向量
X = noisy_ellipse[:,0:1]
Y = noisy_ellipse[:,1:]

# 构造并解决最小二乘问题 ||Mx - b ||^2
M = np.hstack([X**2, X * Y, Y**2, X, Y])
b = np.ones_like(X)
x = np.linalg.lstsq(M, b)[0].squeeze()

# 打印标准形式的椭圆方程
print('椭圆方程为 {0:.3}x^2 + {1:.3}xy+{2:.3}y^2+{3:.3}x+{4:.3}y = 1'.format(x[0], x[1],x[2],x[3],x[4]))

# 绘制噪声数据
plt.scatter(X, Y, label='数据点')

# 绘制生成数据的原始椭圆
phi = np.linspace(0, 2*np.pi, 1000).reshape((1000,1))
c = np.hstack([np.cos(phi), np.sin(phi)])
ground_truth_ellipse = c.dot(B)
plt.plot(ground_truth_ellipse[:,0], ground_truth_ellipse[:,1], 'k--', label='生成椭圆')

# 绘制最小二乘椭圆
x_coord = np.linspace(-5,5,300)
y_coord = np.linspace(-5,5,300)
X_coord, Y_coord = np.meshgrid(x_coord, y_coord)
Z_coord = x[0] * X_coord ** 2 + x[1] * X_coord * Y_coord + x[2] * Y_coord**2 + x[3] * X_coord + x[4] * Y_coord
plt.contour(X_coord, Y_coord, Z_coord, levels=[1], colors=('r'), linewidths=2)

plt.legend()
plt.xlabel('X')
plt.ylabel('Y')
plt.show()

然而,我需要使用完整的一般方程:Ax^2 + Bxy + Cy^2 + Dx + Ey + F = 1,该方程来自此处 并包括椭圆的倾斜和平移。

在上面的代码中,矩阵 M:

M = np.hstack([X**2, X * Y, Y**2, X, Y])

假定我正在拟合 Ax^2 + Bxy + Cy^2 + Dx + Ey = 1。

问题: 如何“添加” F 系数,或者使用这种方法是否不可能?

英文:

I would like to utilize the following code to fit some 2D data to an ellipse, which I got from this post.

import numpy as np
import matplotlib.pyplot as plt
alpha = 5
beta = 3
N = 500
DIM = 2
np.random.seed(2)
# Generate random points on the unit circle by sampling uniform angles
theta = np.random.uniform(0, 2*np.pi, (N,1))
eps_noise = 0.2 * np.random.normal(size=[N,1])
circle = np.hstack([np.cos(theta), np.sin(theta)])
# Stretch and rotate circle to an ellipse with random linear tranformation
B = np.random.randint(-3, 3, (DIM, DIM))
noisy_ellipse = circle.dot(B) + eps_noise
# Extract x coords and y coords of the ellipse as column vectors
X = noisy_ellipse[:,0:1]
Y = noisy_ellipse[:,1:]
# Formulate and solve the least squares problem ||Mx - b ||^2
M = np.hstack([X**2, X * Y, Y**2, X, Y])
b = np.ones_like(X)
x = np.linalg.lstsq(M, b)[0].squeeze()
# Print the equation of the ellipse in standard form
print('The ellipse is given by {0:.3}x^2 + {1:.3}xy+{2:.3}y^2+{3:.3}x+{4:.3}y = 1'.format(x[0], x[1],x[2],x[3],x[4]))
# Plot the noisy data
plt.scatter(X, Y, label='Data Points')
# Plot the original ellipse from which the data was generated
phi = np.linspace(0, 2*np.pi, 1000).reshape((1000,1))
c = np.hstack([np.cos(phi), np.sin(phi)])
ground_truth_ellipse = c.dot(B)
plt.plot(ground_truth_ellipse[:,0], ground_truth_ellipse[:,1], 'k--', label='Generating Ellipse')
# Plot the least squares ellipse
x_coord = np.linspace(-5,5,300)
y_coord = np.linspace(-5,5,300)
X_coord, Y_coord = np.meshgrid(x_coord, y_coord)
Z_coord = x[0] * X_coord ** 2 + x[1] * X_coord * Y_coord + x[2] * Y_coord**2 + x[3] * X_coord + x[4] * Y_coord
plt.contour(X_coord, Y_coord, Z_coord, levels=[1], colors=('r'), linewidths=2)
plt.legend()
plt.xlabel('X')
plt.ylabel('Y')
plt.show()

However, I need to use the full general equation: Ax^2 + Bxy + Cy^2 + Dx + Ey + F = 1, which I got from here and includes tilt and translation of the ellipse.

In the above code, the M matrix:

M = np.hstack([X**2, X * Y, Y**2, X, Y])

assumes I'm fitting Ax^2 + Bxy + Cy^2 + Dx + Ey = 1.

Question: How do I "add" in the F coefficient, or is that not possible using this method?

答案1

得分: 1

以下是翻译好的部分:

"Indeed it is possible with a little knowledge in optimization. We know that the trivial solution for F is certainly F=1 and other parameters equal to 0, but we don't want that and there is probably another solution other than the trivial solution. So let's get into some theory."

"的确,通过一点优化知识,这是可能的。我们知道F的平凡解肯定是 F=1,而其他参数都为0,但我们不想要那个,很可能除了平凡解之外还有其他解。所以让我们深入一些理论。"

"In the standard least squares method we have"

"在标准最小二乘法中,我们有"

"But there are several ways to extend this method so that it can meet the demands of each problem. In your case, what we know is that the trivial solution resets the parameters A, B, C, D and E, that is, the simplest solution is to think that these parameters cannot be Zero"

"但有几种方法可以扩展这个方法,使其满足每个问题的需求。在你的情况下,我们知道平凡解会重置参数 A, B, C, D 和 E,也就是说,最简单的解决方案是认为这些参数不能Zero"

"So, a solution is to add a regularizer that adds to our objective function, a term that adds an error smaller and smaller if our parameters distance from zero. One of the ways is: to obtain the inverse of the sum of the parameters we have (the higher the sum of the parameters, the smaller the penalty from the regularizer will be.), in the way that Regularizer = alpha/sum(|x|). And sum(|x|) is L1 norm in Linear Algebra. Increase alpha to increase the strength of the regularizer"

"所以,一个解决方案是添加一个正则化项到我们的目标函数中,这个项会在我们的参数远离零时添加一个越来越小的误差。其中一种方法是:获得我们拥有的参数之和的倒数(参数之和越大,正则化项的惩罚就越小),方式是 Regularizer = alpha/sum(|x|)。而 sum(|x|) 是线性代数中的L1范数。增加 alpha 可以增加正则化的强度"

"So our LSTSQ is"

"所以我们的LSTSQ是"

"but be careful, because the norm of x cannot be zero, we will deal with this in code."

"但要小心,因为 x 的范数不能为零,我们将在代码中处理这个问题。"

"So, let's code"

"所以,让我们编写代码"

"The ellipse is given by -0.337x^2 + -0.143xy + -0.543y^2 + -0.0325x + -0.0331y + 5.33 = 1"

"椭圆方程为 -0.337x^2 + -0.143xy + -0.543y^2 + -0.0325x + -0.0331y + 5.33 = 1"

"Minimize of Scipy"

"Scipy的最小化函数"

"Basically the 'minimize' function works by iterations, where it calculates where the cost function has the lowest value, so it adjusts the parameters at each iteration until it finds the best ones given a tolerance and maximum number of iterations allowed."

"基本上,'minimize' 函数通过迭代工作,它计算成本函数的最低值出现在哪里,然后在每次迭代中调整参数,直到找到在给定容差和允许的最大迭代次数下的最佳参数。"

"See more on docs: https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html"

"查看更多文档:https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html"

英文:

Indeed it is possible with a little knowledge in optimization. We know that the trivial solution for F is certainly F=1 and other parameters equal to 0, but we don't want that and there is probably another solution other than the trivial solution. So let's get into some theory.

In the standard least squares method we have

Standard least squares
用通用椭圆方程(包括倾斜和平移)最小二乘拟合椭圆。

But there are several ways to extend this method so that it can meet the demands of each problem. In your case, what we know is that the trivial solution resets the parameters A, B, C, D and E, that is, the simplest solution is to think that these for meters cannot be Zero

So, a solution is to add a regularizer that adds to our objective function, a term that adds an error smaller and smaller if our parameters distance from zero.
One of the ways is: to obtain the inverse of the sum of the parameters we have (the higher the sum of the parameters, the smaller the penalty from the regularizer will be.), in the way that Regularizer = alpha/sum(|x|). And sum(|x|) is L1 norm in Linear Algebra. Increase alpha to increase the strength of the regularizer

So our LSTSQ is

Our modified least squares
用通用椭圆方程(包括倾斜和平移)最小二乘拟合椭圆。

but be careful, because the norm of x cannot be zero, we will deal with this in code.

So, let's code

Step 1

Import minimize from scipy

from scipy.optimize import minimize

Step 2

Create your own minimize function

(the adjusted variable always comes before the other parameters, that's why x comes before)

def target_func(x, M, b):
    norm = np.linalg.norm(x, ord=1)     # get L1 norm of x
    norm = norm if norm != 0 else 1e-6  # Verify Zero, in case Zero, place a tolerance (1e-6 in this example)
    reg = 100/(norm)                    # Evaluate regularizer (using alpha= 100 for this exemple)
    return np.linalg.norm(M@x - b, ord=2) + reg # Write de minimize function based on lstsq

Step 3

Create parameters and apply on minimize function

M = np.hstack([X**2, X * Y, Y**2, X, Y, np.ones_like(X)]) # Add ones for param F
b = np.ones_like(X).flatten()                             # Create b (as flatten for prevent problems)
x0 = np.random.rand(6)                                    # Create Initial solution for our problem

res_log = minimize(target_func,                           # Our minimize target
                   x0,                                    # Initial solution
                   args=(M, b),                           # Args for problem (M and b)
)
x = res_log.x # Get the x 

Step 4

Execute the same code for plot, now with the x[5] (F param).

# Print the equation of the ellipse in standard form
print('The ellipse is given by {0:.3}x^2 + {1:.3}xy + {2:.3}y^2 + {3:.3}x + {4:.3}y + {5:.3} = 1'.format(*x))

# Plot the noisy data
plt.scatter(X, Y, label='Data Points')

# Plot the original ellipse from which the data was generated
phi = np.linspace(0, 2*np.pi, 1000).reshape((1000,1))
c = np.hstack([np.cos(phi), np.sin(phi)])
ground_truth_ellipse = c.dot(B)
plt.plot(ground_truth_ellipse[:,0], ground_truth_ellipse[:,1], 'k--', label='Generating Ellipse')

# Plot the least squares ellipse
x_coord = np.linspace(-5,5,300)
y_coord = np.linspace(-5,5,300)
X_coord, Y_coord = np.meshgrid(x_coord, y_coord)
Z_coord = x[0] * X_coord ** 2 + x[1] * X_coord * Y_coord + x[2] * Y_coord**2 + x[3] * X_coord + x[4] * Y_coord + x[5] # Add x[5]
plt.contour(X_coord, Y_coord, Z_coord, levels=[1], colors=('r'), linewidths=2)

plt.legend()
plt.xlabel('X')
plt.ylabel('Y')
plt.show()

Result

The ellipse is given by -0.337x^2 + -0.143xy + -0.543y^2 + -0.0325x + -0.0331y + 5.33 = 1

用通用椭圆方程(包括倾斜和平移)最小二乘拟合椭圆。

Minimize of Scipy

Basically the "minimize" function works by iterations, where it calculates where the cost function has the lowest value, so it adjusts the parameters at each iteration until it finds the best ones given a tolerance and maximum number of iterations allowed.

See more on docs: https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html

huangapple
  • 本文由 发表于 2023年6月15日 00:54:27
  • 转载请务必保留本文链接:https://go.coder-hub.com/76475905.html
匿名

发表评论

匿名网友

:?: :razz: :sad: :evil: :!: :smile: :oops: :grin: :eek: :shock: :???: :cool: :lol: :mad: :twisted: :roll: :wink: :idea: :arrow: :neutral: :cry: :mrgreen:

确定