使用 astropy [Ellipse2d 模型] 拟合椭圆

Question

我正在尝试使用Ellipse2D 模型在astropy 库中拟合椭圆。合身不起作用。建模参数与初始参数相同（可能除了幅度参数）。请看下面的代码：

import numpy as np
from astropy.modeling import models, fitting
import matplotlib.pyplot as pl

# fake data
num = 100
x, y = np.meshgrid(np.linspace(-5., 5., num), np.linspace(-5, 5, num))
e0 = models.Ellipse2D(amplitude=1., x_0=0., y_0=0., a=2, b=1, theta=0.)
z0 = e0(x, y)
print 'DATA:\n', e0, '\n\n'

# initial model
ei = models.Ellipse2D(amplitude=1., x_0=0.0, y_0=0.0, a=2, b=2, theta=0.2)

fi = fitting.LevMarLSQFitter()

#fitted model?
e1 = fi(ei, x, y, z0)
z1 = e1(x, y)
print 'MODEL:\n', e1, '\n\n'

pl.imshow(z0, extent=[-5, 5, -5, 5], alpha=0.5)
pl.imshow(z1, extent=[-5, 5, -5, 5], alpha=0.2)
pl.show()

Answer 1

我一直在在这里或 Astropy 邮件列表中等待这个问题的答案，因为当时我遇到了完全相同的问题。

由于找不到答案，我决定在找出您/我的代码问题之前不使用 Ellipse2D，而是使用 Gaussian2D 获取theta 参数。

您可以尝试以下代码。我只修改了你的代码。

import numpy as np
from astropy.modeling import models, fitting
import matplotlib.pyplot as pl

#%%
# data
num = 100
x, y = np.meshgrid(np.linspace(-5., 5., num), np.linspace(-5, 5, num))
e0 = models.Ellipse2D(amplitude=1., x_0=0., y_0=0., a=2, b=1, theta=0.)
z0 = e0(x, y)
print ('DATA:\n', e0, '\n\n')

#%%
# initial model
ei = models.Ellipse2D(amplitude=1., x_0=0.1, y_0=0.1, a=3, b=2, theta=0.2)
gi = models.Gaussian2D(amplitude=1., x_mean=0.1, y_mean=0.1,
                       x_stddev=3, y_stddev=2, theta=0.2)
fi = fitting.LevMarLSQFitter()

#%% 
# fitted model?
e1 = fi(ei, x, y, z0)
g1 = fi(gi, x, y, z0)
z1 = e1(x, y)
z2 = g1(x, y)
print('MODEL:\n', e1, '\n\n')
print('MODEL:\n', g1, '\n\n')

pl.imshow(z0, extent=[-5, 5, -5, 5], alpha=0.5)
pl.imshow(z1, extent=[-5, 5, -5, 5], alpha=0.2)
pl.imshow(z2, extent=[-5, 5, -5, 5], alpha=0.5)
pl.colorbar()
pl.show()
print(g1.theta.value)

虽然它不适合给定的具有恒定幅度的椭圆形高原，但它仍然给出了正确的theta 值1.23386185422e-10，它实际上为零。当我将e0 的theta 更改为一些不同的值时，它确实给出了正确的值。

希望对您有所帮助！

Answer 2

如前所述，另一个装配工可能更适合此任务，例如SimplexLSQFitter。

这并不完全适合椭圆，但至少 b 参数几乎可以很好地匹配：

...
fi = fitting.SimplexLSQFitter()
e1 = fi(ei, x, y, z0)
z1 = e1(x, y)
print(repr(e1))
# <Ellipse2D(amplitude=0.8765330382805181, 
#            x_0=0.00027076793418705464, 
#            y_0=0.0008061856852329963, 
#            a=2.0019138872185174, 
#            b=1.0985760645823452, 
#            theta=0.22591442574477916)>

但我认为Ellipse2D 不是一个适合拟合的好模型，特别是如果theta 在模型之间存在差异。

Answer 3

感谢@yoonsoo-p-bach 的建议，这里是工作示例：

import numpy as np
from astropy.modeling import models, fitting
import matplotlib.pyplot as pl

# data
num = 100
x, y = np.meshgrid(np.linspace(-5., 5., num), np.linspace(-5, 5, num))
e0 = models.Ellipse2D(amplitude=1., x_0=0.2, y_0=0.3, a=2, b=1, theta=0.4)
z0 = e0(x, y)
print 'DATA:\n', e0, '\n\n'

# fitting procedure
fi = fitting.SimplexLSQFitter() 
#fi = fitting.LevMarLSQFitter()

# gaussian fit (to estimate x_0, y_0 and theta)
gi = models.Gaussian2D(amplitude=1., x_mean=0.1, y_mean=0.2, x_stddev=1, y_stddev=1, theta=0.0)
g1 = fi(gi, x, y, z0, maxiter=1000)
print 'Gaussian:\n', g1, '\n\n'

# initial model
ei = models.Ellipse2D(amplitude=1., x_0=g1.x_mean, y_0=g1.y_mean, a=g1.x_stddev, b=g1.y_stddev, theta=g1.theta, fixed={'x_0': True, 'y_0':True, 'theta':True})

#fitted model
e1 = fi(ei, x, y, z0, maxiter=1000)
z1 = e1(x, y)
print 'MODEL:\n', e1, '\n\n'

pl.imshow(z0-z1, extent=[-5, 5, -5, 5])
pl.show()