CODEWITHSUNDEEP
pythongaussianringastropymodel-fitting
Isaac Wong
Isaac Wong

Reputation: 23

How to make a fittable Ring Gaussian model in python?

The purpose of my experiment is to fit a ring gaussian model to image data and find out the parameters of elliptical or ring Gaussian object in an image.

I've tried Astropy to make a Ring Gaussian model for simplicity and trial. Unfortunately, it's completely different from my artificial data.

from astropy.modeling import Fittable2DModel, Parameter, fitting
import matplotlib.pyplot as plt
import numpy as np

class ringGaussian(Fittable2DModel):

    background = Parameter(default=5.)
    amplitude = Parameter(default=1500.)
    x0 = Parameter(default=15.)
    y0 = Parameter(default=15.)
    radius = Parameter(default=2.)
    width = Parameter(default=1.)

    @staticmethod
    def evaluate(x, y, background, amplitude, x0, y0, radius, width):
        z = background + amplitude * np.exp( -((np.sqrt((x-x0)**2. + (y-y0)**2.) - radius) / width)**2 )
        return z

Then I made some artificial data (initial parameters) to test the fitting function of ring gaussian class.

back = 10 # background
amp = 2000 # highest value of the Gaussian
x0 = 10 # x coordinate of center
y0 = 10 # y coordinate of center
radius = 3
width = 1
xx, yy = np.mgrid[:30, :30]
z = back + amp * np.exp( -((np.sqrt((xx-x0)**2. + (yy-y0)**2.) - radius) / width)**2 )

I plotted xx, yy and z using contourf:

fig = plt.subplots(figsize=(7,7))
plt.contourf(xx,yy,z)
plt.show()

It is what I got: enter image description here

Then I tried to fit the z using my fittable class:

p_init = ringGaussian() #bounds={"x0":[0., 20.], "y0":[0., 20.]}
fit_p = fitting.LevMarLSQFitter()
p = fit_p(p_init, xx, yy, z)
# It is the parameter I got:
<ringGaussian(background=133.0085329497139, amplitude=-155.53652181827655, x0=25.573499373946227, y0=25.25813520725603, radius=8.184302497405568, width=-7.273935403490675)>

I plotted the model:

fig = plt.subplots(figsize=(7,7))
plt.contourf(xx,yy,p(xx,yy))
plt.show()

It is what I got: enter image description here

Originally, I also tried to include the derivative in my class:

@staticmethod
def fit_deriv(x, y, background, amplitude, x0, y0, radius, width):
    g = (np.sqrt((x-x0)**2. + (y-y0)**2.) - radius) / width
    z = amplitude * np.exp( -g**2 )

    dg_dx0 = - (x-x0)/np.sqrt((x-x0)**2. + (y-y0)**2.)
    dg_dy0 = - (y-y0)/np.sqrt((x-x0)**2. + (y-y0)**2.)
    dg_dr0 = - 1/width
    dg_dw0 = g * -1/width

    dz_dB = 1.
    dz_dA = z / amplitude
    dz_dx0 = -2 * z * g**3 * dg_dx0
    dz_dy0 = -2 * z * g**3 * dg_dy0
    dz_dr0 = -2 * z * g**3 * dg_dr0
    dz_dw0 = -2 * z * g**3 * dg_dw0

    return [dz_dB, dz_dA, dz_dx0, dz_dy0, dz_dr0, dz_dw0]

But it returned "ValueError: setting an array element with a sequence."

I'm quite desperate now. Can anyone suggest some possible solutions? or alternative ways to implement the ring gaussian fit in python?

Many many thanks~~~

Upvotes: 0

Views: 909

Answers (1)

Christoph
Christoph

Reputation: 2970

Your implementation is OK (I did not try or check your fit_deriv implementation).

The issue is just that your fit doesn't converge because the initial model parameters are too far from the true values, so the optimiser fails. When I run your code, I get this warning:

WARNING: The fit may be unsuccessful; check fit_info['message'] for more information. [astropy.modeling.fitting]

If you change to model parameters so that your model roughly matches the data, the fit succeeds:

p_init = ringGaussian(x0=11, y0=11)

To check how your model compares with the data, you can use imshow to display data and model images (or also e.g. residual images):

plt.imshow(z)
plt.imshow(p_init(xx,yy))
plt.imshow(p(xx,yy))

Upvotes: 2

Related Questions