gaussian2d_example.py

gistfile1.py

from astropy.modeling import models, fitting
import numpy as np
from scipy.optimize import leastsq

def gaussian2d(p, x, y):
    amplitude, x_mean, y_mean, x_stddev, y_stddev, theta = p[:]
    a = 0.5 * ((np.cos(theta) / x_stddev) ** 2 + (np.sin(theta) / y_stddev) ** 2)
    b = 0.5 * (np.cos(theta) * np.sin(theta) * (1. / x_stddev ** 2 - 1. / y_stddev ** 2))
    c = 0.5 * ((np.sin(theta) / x_stddev) ** 2 + (np.cos(theta) / y_stddev) ** 2)
    return amplitude * np.exp(-a * (x - x_mean) ** 2 - b * (x - x_mean) * (y - y_mean) - c * (y - y_mean) ** 2)


def error_function(p, x, y, z):
    return gaussian2d(p, x, y) - z

X,Y = np.meshgrid(np.arange(11),np.arange(18))
p0 = [137., 5.1, 5.4, 1.5, 2., np.pi/4]
data = gaussian2d(p0, X.ravel(), Y.ravel()).reshape(X.shape)
data += np.random.normal(0., 1., size=data.shape)

gauss = models.Gaussian2DModel(amplitude=10., x_mean=5., y_mean=5.,
x_stddev=4., y_stddev=4., theta=0.5)
gauss_fit=fitting.NonLinearLSQFitter(gauss)
gauss_fit(X, Y, data, estimate_jacobian=True)

gauss1= models.Gaussian2DModel(amplitude=10., x_mean=5., y_mean=5.,
x_stddev=4., y_stddev=4., theta=0.5)
gauss_fit1=fitting.NonLinearLSQFitter(gauss1)
gauss_fit1(X, Y, data)

p_opt, success = leastsq(error_function,
[10., 5., 5., 4., 4., 0.5],
args=(X.ravel(),Y.ravel(),data.ravel()),maxfev=1000, col_deriv=1)

print(gauss.parameters)
print(gauss1.parameters)
print(p_opt)