3

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()
Improve this question
2
  • Have you tried a different fitter? I wonder if this would really work since the Ellipse2D model seems to make a solid ellipse. Commented Feb 11, 2017 at 0:36
  • Yes I have tried all of them. The SimplexLSQFitter (see below) gives somehow better results, but it is still far from satisfying. Commented Feb 13, 2017 at 9:53

3 Answers 3

Reset to default
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!

Improve this answer
Sign up to request clarification or add additional context in comments.

Comments

0

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.

Improve this answer

1 Comment

Thank you, indeed SimplexLSQFitter gives better (any) results, but the main parameter I am trying to find is the inclination (theta) and it fails there. I thought that Ellipse2D should be the best to find the inclination of an ellipse, but apparently I was wrong. Do you know any better method to do that?
0

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()
Improve this answer

Comments

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.