Fitting an ellipse using astropy [Ellipse2d model]

Question

I am trying to fit an ellipse using Ellipse2D model in astropy library. The fit does not work. The modeled parameters are the same as initial parameters (maybe except the amplitude parameter). See the code below:

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

Thanks to suggestion by @yoonsoo-p-bach here is the working example:

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

Answer 2

I was waiting for the answer to this question for here or Astropy mailing list, since I was having exactly the same problem at that time.

As I couldn't find the answer, I decided not to use Ellipse2D until I figure out the problem of your/my code, but to use Gaussian2D for getting the theta parameter.

You may try the following code. I changed only a slight bit of your code.

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)

Although it does not fit the given ellipse-shaped plateau with constant amplitude, but still it gives the correct theta value 1.23386185422e-10 which is effectively zero. It does give correct values when I change theta of e0 to some different values.

Hope it helped!

Answer 3

As already noted another fitter might be better suited for this task, for example the SimplexLSQFitter.

This doesn't perfectly fit the ellipse but at least the b parameter almost gives a good match:

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

But I think Ellipse2D isn't a good model to be fitted, especially if theta differs between the models.