2017-02-09 87 views
3

我想在astropy庫中使用Ellipse2D模型擬合一個橢圓。適合不起作用。建模參數與初始參數相同(可能除幅度參數外)。請參見下面的代碼:使用astropy擬合一個橢圓[Ellipse2d模型]

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() 
+0

您已經嘗試了不同的鉗工?我想知道這是否真的能起作用,因爲Ellipse2D模型看起來是一個堅實的橢圓。 – Iguananaut

+0

是的,我嘗試了所有這些。 SimplexLSQFitter(見下文)給出了某種更好的結果,但它還遠遠不能令人滿意。 –

回答

0

如前所述另一個鉗工可能更適合這項任務,例如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之間的區別模型。

+0

謝謝你的確,SimplexLSQFitter提供了更好的(任何)結果,但是我試圖找到的主要參數是傾斜度(theta),它在那裏失效了。 我認爲Ellipse2D應該是最好的找到橢圓的傾向,但顯然我錯了。你知道有什麼更好的方法來做到這一點嗎? –

2

我在等待這個問題的答案或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) 

雖然它不適合用恆定的幅度給定的橢圓形的高原,但仍提供了正確的theta1.23386185422e-10這實際上爲零。當我將e0theta更改爲某些不同的值時,它會給出正確的值。

希望它有幫助!

+0

它幫助,謝謝。 –

0

由於建議由@ yoonsoo - 對巴赫這裏是工作示例:

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()