@@ -25,11 +25,12 @@ f(x | x̂, Q) = exp[-4ln(2) * (x - x̂)ᵀ Q (x - x̂)]
25
25
where `Q` is the inverse covariance matrix (or precision matrix). This is equivalent to the inverse of the FWHM matrix after squaring each element.
26
26
"""
27
27
struct Gaussian{T,FT,VT<: AbstractVector ,IT<: Tuple } <: PSFModel{T}
28
+ amp:: T
28
29
pos:: VT
29
30
fwhm:: FT
30
31
indices:: IT
31
32
32
- Gaussian {T} (pos:: VT , fwhm:: FT , indices:: IT ) where {T,VT<: AbstractVector ,FT,IT<: Tuple } = new {T,FT,VT,IT} (pos, fwhm, indices)
33
+ Gaussian {T} (amp :: T , pos:: VT , fwhm:: FT , indices:: IT ) where {T,VT<: AbstractVector ,FT,IT<: Tuple } = new {T,FT,VT,IT} (amp, pos, fwhm, indices)
33
34
end
34
35
35
36
# Alias Normal -> Gaussian
@@ -45,36 +46,38 @@ const Normal = Gaussian
45
46
# default type is Float64
46
47
Gaussian (args... ; kwargs... ) = Gaussian {Float64} (args... ; kwargs... )
47
48
# parse indices from maxsize
48
- Gaussian {T} (pos:: AbstractVector , fwhm; maxsize= 3 ) where {T} = Gaussian {T} (pos, fwhm, indices_from_extent (pos, fwhm, maxsize))
49
+ Gaussian {T} (pos:: AbstractVector ; fwhm, amp = one (T), maxsize= 3 ) where {T} = Gaussian {T} (amp, pos, fwhm, indices_from_extent (pos, fwhm, maxsize))
49
50
# default position is [0, 0]
50
- Gaussian {T} (fwhm ; kwargs... ) where {T} = Gaussian {T} (SA[0 , 0 ], fwhm ; kwargs... )
51
+ Gaussian {T} (; kwargs... ) where {T} = Gaussian {T} (SA[0 , 0 ]; kwargs... )
51
52
# # parse position to vector
52
- Gaussian {T} (x:: Number , y:: Number , fwhm ; kwargs... ) where {T} = Gaussian {T} (SA[x, y], fwhm ; kwargs... )
53
- Gaussian {T} (xy:: Tuple , fwhm ; kwargs... ) where {T} = Gaussian {T} (SVector (xy), fwhm ; kwargs... )
53
+ Gaussian {T} (x:: Number , y:: Number ; kwargs... ) where {T} = Gaussian {T} (SA[x, y]; kwargs... )
54
+ Gaussian {T} (xy:: Tuple ; kwargs... ) where {T} = Gaussian {T} (SVector (xy); kwargs... )
54
55
# # translate polar coordinates to cartesian, optionally recentering
55
- Gaussian {T} (p:: Polar , fwhm ; origin= SA[0 , 0 ], kwargs... ) where {T} = Gaussian {T} (CartesianFromPolar ()(p) .+ origin, fwhm ; kwargs... )
56
+ Gaussian {T} (p:: Polar ; origin= SA[0 , 0 ], kwargs... ) where {T} = Gaussian {T} (CartesianFromPolar ()(p) .+ origin; kwargs... )
56
57
57
58
Base. size (g:: Gaussian ) = map (length, g. indices)
58
59
Base. axes (g:: Gaussian ) = g. indices
59
60
61
+ const GAUSS_PRE = - 4 * log (2 )
62
+
60
63
# covers scalar case
61
64
function (g:: Gaussian{T} )(point:: AbstractVector ) where T
62
65
Δ = sqeuclidean (point, g. pos)
63
- val = exp ( - 4 * log ( 2 ) * Δ / g. fwhm^ 2 )
66
+ val = g . amp * exp (GAUSS_PRE * Δ / g. fwhm^ 2 )
64
67
return convert (T, val)
65
68
end
66
69
# covers vector case
67
70
function (g:: Gaussian{T,<:Union{Tuple,AbstractVector}} )(point:: AbstractVector ) where T
68
71
weights = SA[1 / first (g. fwhm)^ 2 , 1 / last (g. fwhm)^ 2 ] # manually invert
69
72
Δ = wsqeuclidean (point, g. pos, weights)
70
- val = exp ( - 4 * log ( 2 ) * Δ)
73
+ val = g . amp * exp (GAUSS_PRE * Δ)
71
74
return convert (T, val)
72
75
end
73
76
74
77
# matrix case
75
78
function (g:: Gaussian{T,<:AbstractMatrix} )(point:: AbstractVector ) where T
76
79
R = point - g. pos
77
80
Δ = R' * ((g. fwhm .^ 2 ) \ R)
78
- val = exp ( - 4 * log ( 2 ) * Δ)
81
+ val = g . amp * exp (GAUSS_PRE * Δ)
79
82
return convert (T, val)
80
83
end
0 commit comments