Skip to content

Commit

Permalink
Merge pull request #2038 from larrybradley/gauss2d
Browse files Browse the repository at this point in the history 
Fixed Gaussian2D model
  • Loading branch information
@nden
nden committed Feb 21, 2014
2 parents b3239f1 + c5011c9 commit 5d2c722
Show file tree
Hide file tree
Showing 3 changed files with 127 additions and 74 deletions.
2 changes: 2 additions & 0 deletions CHANGES.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -313,6 +313,8 @@ Bug Fixes

- Raise a `NotImplementedError` when fitting composite models. [#1915]

- Fixed bug in computation of ``Gaussian2D`` model. [#2038]

- ``astropy.nddata``

- ``astropy.sphinx``
Expand Down
163 changes: 89 additions & 74 deletions astropy/modeling/functional_models.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -34,11 +34,11 @@ class Gaussian1D(Parametric1DModel):
Parameters
----------
amplitude : float
Amplitude of the gaussian
Amplitude of the Gaussian.
mean : float
Mean of the gaussian
Mean of the Gaussian.
stddev : float
Standard deviation of the gaussian
Standard deviation of the Gaussian.
Notes
-----
Expand Down Expand Up @@ -102,17 +102,18 @@ def __init__(self, amplitude, mean, stddev, **constraints):
super(Gaussian1D, self).__init__(param_dim=param_dim,
amplitude=amplitude, mean=mean,
stddev=stddev, **constraints)

@staticmethod
def eval(x, amplitude, mean, stddev):
"""
Model function Gauss1D
Gaussian1D model function.
"""
return amplitude * np.exp(- 0.5 * (x - mean) ** 2 / stddev ** 2)

@staticmethod
def fit_deriv(x, amplitude, mean, stddev):
"""
Model function derivatives Gauss1D
Gaussian1D model function derivatives.
"""

d_amplitude = np.exp(-0.5 / stddev ** 2 * (x - mean) ** 2)
Expand All @@ -128,22 +129,24 @@ class Gaussian2D(Parametric2DModel):
Parameters
----------
amplitude : float
Amplitude of the gaussian
Amplitude of the Gaussian.
x_mean : float
Mean of the gaussian in x
Mean of the Gaussian in x.
y_mean : float
Mean of the gaussian in y
Mean of the Gaussian in y.
x_stddev : float
Standard deviation of the gaussian in x
Either x_fwhm or x_stddev must be specified
Standard deviation of the Gaussian in x.
x_stddev and y_stddev must be specified unless a covariance
matrix (cov_matrix) is input.
y_stddev : float
Standard deviation of the gaussian in y
Either y_fwhm or y_stddev must be specified
theta : float
Rotation angle in radians. Note: increases clockwise.
cov_matrix : ndarray
A 2x2 covariance matrix. If specified, overrides stddev, fwhm, and
theta specification.
Standard deviation of the Gaussian in y.
x_stddev and y_stddev must be specified unless a covariance
matrix (cov_matrix) is input.
theta : float, optional
Rotation angle in radians. The rotation angle increases clockwise.
cov_matrix : ndarray, optional
A 2x2 covariance matrix. If specified, overrides the x_stddev,
y_stddev, and theta specification.
Notes
-----
Expand All @@ -157,14 +160,14 @@ class Gaussian2D(Parametric2DModel):
Using the following definitions:
.. math::
a = \\left(- \\frac{\\sin^{2}{\\left (\\theta \\right )}}{2 \\sigma_{y}^{2}} -
\\frac{\\cos^{2}{\\left (\\theta \\right )}}{2 \\sigma_{x}^{2}}\\right)
a = \\left(\\frac{\\cos^{2}{\\left (\\theta \\right )}}{2 \\sigma_{x}^{2}} +
\\frac{\\sin^{2}{\\left (\\theta \\right )}}{2 \\sigma_{y}^{2}}\\right)
b = \\left(\\frac{\\sin{\\left (2 \\theta \\right )}}{2 \\sigma_{y}^{2}} -
\\frac{\\sin{\\left (2 \\theta \\right )}}{2 \\sigma_{x}^{2}}\\right)
b = \\left(\\frac{-\\sin{\\left (2 \\theta \\right )}}{2 \\sigma_{x}^{2}} +
\\frac{\\sin{\\left (2 \\theta \\right )}}{2 \\sigma_{y}^{2}}\\right)
c = \\left(\\frac{\\cos^{2}{\\left (\\theta \\right )}}{2 \\sigma_{y}^{2}} +
\\frac{\\sin^{2}{\\left (\\theta \\right )}}{2 \\sigma_{x}^{2}}\\right)
c = \\left(\\frac{\\sin^{2}{\\left (\\theta \\right )}}{2 \\sigma_{x}^{2}} +
\\frac{\\cos^{2}{\\left (\\theta \\right )}}{2 \\sigma_{y}^{2}}\\right)
See Also
Expand All @@ -183,11 +186,11 @@ def __init__(self, amplitude, x_mean, y_mean, x_stddev=None, y_stddev=None,
theta=0.0, cov_matrix=None, **constraints):
if y_stddev is None and cov_matrix is None:
raise InputParameterError(
"Either y_stddev must be specified, or a "
"Either x/y_stddev must be specified, or a "
"covariance matrix.")
elif x_stddev is None and cov_matrix is None:
raise InputParameterError(
"Either x_stddev must be specified, or a "
"Either x/y_stddev must be specified, or a "
"covariance matrix.")
elif cov_matrix is not None and (x_stddev is not None or
y_stddev is not None):
Expand All @@ -211,54 +214,65 @@ def __init__(self, amplitude, x_mean, y_mean, x_stddev=None, y_stddev=None,
def eval(x, y, amplitude, x_mean, y_mean, x_stddev, y_stddev, theta):
"""Two dimensional Gaussian function"""

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))
cost2 = np.cos(theta) ** 2
sint2 = np.sin(theta) ** 2
sin2t = np.sin(2. * theta)
xstd2 = x_stddev ** 2
ystd2 = y_stddev ** 2
xdiff = x - x_mean
ydiff = y - y_mean
a = 0.5 * ((cost2 / xstd2) + (sint2 / ystd2))
b = 0.5 * (-(sin2t / xstd2) + (sin2t / ystd2))
c = 0.5 * ((sint2 / xstd2) + (cost2 / ystd2))
return amplitude * np.exp(-((a * xdiff ** 2) + (b * xdiff * ydiff) +
(c * ydiff ** 2)))

@staticmethod
def fit_deriv(x, y, amplitude, x_mean, y_mean, x_stddev, y_stddev, theta):
"""Two dimensional Gaussian function derivative with respect to parameters"""

# Helper quantities
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)

da_dtheta = np.sin(theta) * np.cos(theta) * (1/y_stddev**2 - 1/x_stddev**2)
da_dx_stddev = -np.cos(theta)**2/x_stddev**3
da_dy_stddev = -np.sin(theta)**2/y_stddev**3
db_dtheta = 0.5*np.cos(2*theta) * (1/x_stddev**2 - 1/y_stddev**2)
db_dx_stddev = -np.cos(theta)*np.sin(theta)/x_stddev**3
db_dy_stddev = np.cos(theta)*np.sin(theta)/y_stddev**3
dc_dtheta = np.cos(theta)*np.sin(theta)*(1/x_stddev**2 - 1/y_stddev**2)
dc_dx_stddev = -np.sin(theta)**2/x_stddev**3
dc_dy_stddev = -np.cos(theta)**2/y_stddev**3

d_A = np.exp(- a * (x - x_mean) ** 2
- b * (x - x_mean) * (y - y_mean)
- c * (y - y_mean) ** 2)
d_theta = amplitude * (-(x - x_mean) ** 2 * da_dtheta - (x - x_mean) *
(y - y_mean) * db_dtheta - (y - y_mean) ** 2 * dc_dtheta) * d_A

d_x_stddev = amplitude * (-(x - x_mean) ** 2 *da_dx_stddev -
(x - x_mean) * (y - y_mean) * db_dx_stddev
- (y - y_mean) ** 2 * dc_dx_stddev) * d_A
d_y_stddev = amplitude * (-(x - x_mean) ** 2 *da_dy_stddev -
(x - x_mean) * (y - y_mean) * db_dy_stddev
- (y - y_mean) ** 2 * dc_dy_stddev) * d_A
d_y_mean = amplitude * (+(x - x_mean) * b + 2 * (y - y_mean) * c) * d_A
d_x_mean = amplitude * ((y - y_mean) * b + 2 * (x - x_mean) * a) * d_A
return [d_A, d_x_mean, d_y_mean, d_x_stddev, d_y_stddev, d_theta]
cost = np.cos(theta)
sint = np.sin(theta)
cost2 = np.cos(theta) ** 2
sint2 = np.sin(theta) ** 2
cos2t = np.cos(2. * theta)
sin2t = np.sin(2. * theta)
xstd2 = x_stddev ** 2
ystd2 = y_stddev ** 2
xstd3 = x_stddev ** 3
ystd3 = y_stddev ** 3
xdiff = x - x_mean
ydiff = y - y_mean
xdiff2 = xdiff ** 2
ydiff2 = ydiff ** 2
a = 0.5 * ((cost2 / xstd2) + (sint2 / ystd2))
b = 0.5 * (-(sin2t / xstd2) + (sin2t / ystd2))
c = 0.5 * ((sint2 / xstd2) + (cost2 / ystd2))
g = amplitude * np.exp(-((a * xdiff2) + (b * xdiff * ydiff) +
(c * ydiff2)))
da_dtheta = (sint * cost * ((1. / ystd2) - (1. / xstd2)))
da_dx_stddev = -cost2 / xstd3
da_dy_stddev = -sint2 / ystd3
db_dtheta = (-cos2t / xstd2) + (cos2t / ystd2)
db_dx_stddev = sin2t / xstd3
db_dy_stddev = -sin2t / ystd3
dc_dtheta = -da_dtheta
dc_dx_stddev = -sint2 / xstd3
dc_dy_stddev = -cost2 / ystd3
dg_dA = g / amplitude
dg_dx_mean = g * ((2. * a * xdiff) + (b * ydiff))
dg_dy_mean = g * ((b * xdiff) + (2. * c * ydiff))
dg_dx_stddev = g * (-(da_dx_stddev * xdiff2 +
db_dx_stddev * xdiff * ydiff +
dc_dx_stddev * ydiff2))
dg_dy_stddev = g * (-(da_dy_stddev * xdiff2 +
db_dy_stddev * xdiff * ydiff +
dc_dy_stddev * ydiff2))
dg_dtheta = g * (-(da_dtheta * xdiff2 +
db_dtheta * xdiff * ydiff +
dc_dtheta * ydiff2))
return [dg_dA, dg_dx_mean, dg_dy_mean, dg_dx_stddev, dg_dy_stddev,
dg_dtheta]


class Shift(Model):
Expand Down Expand Up @@ -423,7 +437,7 @@ class Linear1D(Parametric1DModel):

def __init__(self, slope, intercept, **constraints):
super(Linear1D, self).__init__(slope=slope, intercept=intercept,
**constraints)
**constraints)

@staticmethod
def eval(x, slope, intercept):
Expand Down Expand Up @@ -472,7 +486,7 @@ class Lorentz1D(Parametric1DModel):

def __init__(self, amplitude, x_0, fwhm, **constraints):
super(Lorentz1D, self).__init__(amplitude=amplitude, x_0=x_0,
fwhm=fwhm, **constraints)
fwhm=fwhm, **constraints)

@staticmethod
def eval(x, amplitude, x_0, fwhm):
Expand Down Expand Up @@ -738,7 +752,7 @@ def eval(x, amplitude, x_0, width):

return np.select([np.logical_and(x >= x_0 - width / 2.,
x <= x_0 + width / 2.)],
[amplitude], 0)
[amplitude], 0)

@classmethod
def fit_deriv(cls, x, amplitude, x_0, width):
Expand Down Expand Up @@ -801,8 +815,10 @@ def __init__(self, amplitude, x_0, y_0, x_width, y_width, **constraints):
@staticmethod
def eval(x, y, amplitude, x_0, y_0, x_width, y_width):
"""Two dimensional Box model function"""
x_range = np.logical_and(x >= x_0 - x_width / 2., x <= x_0 + x_width / 2.)
y_range = np.logical_and(y >= y_0 - y_width / 2., y <= y_0 + y_width / 2.)
x_range = np.logical_and(x >= x_0 - x_width / 2.,
x <= x_0 + x_width / 2.)
y_range = np.logical_and(y >= y_0 - y_width / 2.,
y <= y_0 + y_width / 2.)
return np.select([np.logical_and(x_range, y_range)], [amplitude], 0)


Expand Down Expand Up @@ -889,7 +905,6 @@ def __init__(self, amplitude, x_0, y_0, R_0, slope, **constraints):
x_0=x_0, y_0=y_0, R_0=R_0,
slope=slope, **constraints)


@staticmethod
def eval(x, y, amplitude, x_0, y_0, R_0, slope):
"""Two dimensional Trapezoid Disk model function"""
Expand Down
36 changes: 36 additions & 0 deletions astropy/modeling/tests/test_functional_models.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,13 +6,15 @@
import numpy as np
from numpy.testing import assert_allclose
from .. import models
from astropy.coordinates import Angle

try:
from scipy import optimize
HAS_SCIPY = True
except ImportError:
HAS_SCIPY = False


def test_Trapezoid1D():
"""Regression test for
https://github.com/astropy/astropy/issues/1721
Expand All @@ -22,3 +24,37 @@ def test_Trapezoid1D():
yy = model(xx)
yy_ref = [0., 1.41428571, 3.12857143, 4.2, 4.2, 3.12857143, 1.41428571, 0.]
assert_allclose(yy, yy_ref, rtol=0, atol=1e-6)


def test_Gaussian2D():
"""
Test rotated elliptical Gaussian2D model.
https://github.com/astropy/astropy/pull/2038
"""
model = models.Gaussian2D(100, 1.7, 3.1, x_stddev=3.3, y_stddev=5.0,
theta=np.pi/6.)
y, x = np.mgrid[0:5, 0:5]
g = model(x, y)
g_ref = [[77.86787722, 79.8450774, 75.66323816, 66.26262987, 53.6289606],
[86.01743889, 90.20335997, 87.419011, 78.29536082, 64.80568981],
[90.11888627, 96.6492347, 95.79172412, 87.74139348, 74.27249818],
[89.54601432, 98.21442091, 99.55228375, 93.25543692, 80.73169335],
[84.38744554, 94.65710923, 98.12408063, 94.00369635, 83.22642276]]
assert_allclose(g, g_ref, rtol=0, atol=1e-6)


def test_Gaussian2DRotation():
amplitude = 42
x_mean, y_mean = 0, 0
x_stddev, y_stddev = 2, 3
theta = Angle(10, 'deg')
pars = dict(amplitude=amplitude, x_mean=x_mean, y_mean=y_mean,
x_stddev=x_stddev, y_stddev=y_stddev)
rotation_matrix = models.MatrixRotation2D(angle=theta.degree)
point1 = (x_mean + 2 * x_stddev, y_mean + 2 * y_stddev)
point2 = rotation_matrix(*point1)
g1 = models.Gaussian2D(theta=0, **pars)
g2 = models.Gaussian2D(theta=theta.radian, **pars)
value1 = g1(*point1)
value2 = g2(*point2)
assert_allclose(value1, value2)

0 comments on commit 5d2c722

Please sign in to comment.