Add Moffat Kernel

docs/src/function_reference.md

@@ -23,6 +23,7 @@ Kernel.gaussian
 Kernel.DoG
 Kernel.LoG
 Kernel.Laplacian
+Kernel.moffat
 ```
 
 # KernelFactors

src/ImageFiltering.jl

@@ -64,7 +64,7 @@ ArrayLike{T} = Union{ArrayType{T}, AnyIIR{T}}
 
 include("kernel.jl")
 using .Kernel
-using .Kernel: Laplacian, reflect, ando3, ando4, ando5, scharr, bickley, prewitt, sobel, gabor
+using .Kernel: Laplacian, reflect, ando3, ando4, ando5, scharr, bickley, prewitt, sobel, gabor, moffat
 
 NDimKernel{N,K} = Union{AbstractArray{K,N},ReshapedOneD{K,N},Laplacian{N}}
 

src/kernel.jl

@@ -390,6 +390,31 @@ function validate_gabor(σ::Real,λ::Real,γ::Real)
     end
 end
 
+"""
+    moffat(α, β, ls) -> k
+
+Constructs a 2D, symmetric Moffat kernel `k` with core width, `α`, and power, `β`.
+Size of kernel defaults to 4 * full-width-half-max or as specified in `ls`.
+See [this notebook](https://nbviewer.jupyter.org/github/ysbach/AO_2017/blob/master/04_Ground_Based_Concept.ipynb#1.2.-Moffat) for details.
+
+# Citation
+Moffat, A. F. J. "A theoretical investigation of focal stellar images in the photographic emulsion and application to photographic photometry." Astronomy and Astrophysics 3 (1969): 455.
+"""
+function moffat(α::Real, β::Real, ls::Tuple{Integer, Integer})
+    ws = map(n->(ceil(Int,n)>>1), ls)
+    R = CartesianIndices(map(w->IdentityUnitRange(-w:w), ws))
+    α2 = α^2
+    amp = (β - 1)/(π * α2)
+    @. amp*((1+df(R)/α2)^-β)
+end
+moffat(α::Real, β::Real, ls::Integer)    = moffat(α, β, (ls,ls))
+moffat(α::Real, β::Real)                 = moffat(α, β, ceil(Int, (α*2*sqrt(2^(1/β) - 1))*4))
+
+@inline function df(I::CartesianIndex)
+    x = SVector(Tuple(I))
+    sum(x.^2)
+end
+
 """
     reflect(kernel) --> reflectedkernel
 

test/specialty.jl

@@ -216,6 +216,15 @@ using ImageFiltering: IdentityUnitRange
         ImageFiltering.fillbuf_nan[] = false
         @test Kernel.LoG(2.5) == Kernel.LoG((2.5,2.5))
     end
+
+    @testset "moffat" begin
+        α = rand()
+        β = rand()
+        @test Kernel.moffat(α,β,2) == Kernel.moffat(α,β,(2,2))
+
+        fwhm = ceil(Int, (α*2*sqrt(2^(1/β) - 1)) * 4)
+        @test Kernel.moffat(α, β) == Kernel.moffat(α, β, (fwhm, fwhm))
+    end
 end
 
 nothing