Merge branch 'master' of https://github.com/ptiede/VIDA.jl

src/VIDA.jl

@@ -1,4 +1,4 @@
-__precompile__()
+#__precompile__()
 """
     VIDA
 Is a image feature extraction tool for use with EHT images of black holes.
@@ -32,7 +32,7 @@ export
     #Filters
     GaussianRing,SlashedGaussianRing,EllipticalGaussianRing,
     TIDAGaussianRing,GeneralGaussianRing, Constant, AsymGaussian,
-    CosineRing,Disk,ImageFilter,
+    CosineRing,Disk,ImageFilter,LogSpiral,
     #Filter helper functions
     stack,split,unpack,
     #Image functions

src/filters.jl

@@ -38,6 +38,76 @@ function _update(θ::AbstractFilter, p)
     return typeof(θ)(p)
 end
 
+
+@doc """
+    $(SIGNATURES)
+Filter type for a logarithmic spiral segment
+
+## Fields
+$(FIELDS)
+"""
+@with_kw struct LogSpiral{T<:Real} <: AbstractImageFilter
+    """ Radius of the spiral peak brightness """
+    r0::T
+    """ Unit curvature of the logarithmic spiral """
+    κ::T
+    """ thickness of the Gaussian spiral arm """
+    σ::T
+    """ Azimuthal extent of the spiral arm """
+    δϕ::T
+    """ peak brightness location """
+    ξ::T
+    """ x location of disk center in μas """
+    x0::T
+    """ y location of disk center in μas """
+    y0::T
+end
+function LogSpiral(p::Vector{T}) where {T<:Real}
+    @assert length(p) == 7
+    LogSpiral{T}(p[1],p[2],p[3],p[4],p[5],p[6],p[7])
+end
+size(::Type{LogSpiral{T}}) where {T} = 7
+
+
+@inline function (θ::LogSpiral)(x,y)
+    @unpack κ, σ, r0, δϕ, ξ, x0, y0 = θ
+    x′,y′ = x-x0,y-y0
+    #Set up the spiral
+    k = sqrt(1-κ*κ)/κ
+    rc = exp(k*10π) #This finds where we should start our spiral arm from
+    a = r0/rc #Get on the correct logspiral
+
+    r = hypot(x′,y′)
+    α = (atan(y′,x′)) - ξ
+
+    #Now I need to find the distance from the closest spiral arm
+    n = (log(r/a)/k - α)/(2π)
+    nc = ceil(n)
+    nf = floor(n)
+    rc = a*exp(k*(α + nc*2π))
+    rf = a*exp(k*(α + nf*2π))
+    r1,r2 = abs(rc-r),abs(rf-r)
+    if r1 < r2
+        nn = nc
+        dist = r1
+    else
+        nn = nf
+        dist = r2
+    end
+    #Get the angular extent
+    dtheta = (10π - (α + nn*2π))
+    return exp(-dist^2/(2*σ^2) -dtheta^2/(2*(δϕ/2)^2))
+end
+
+@doc """
+    $(SIGNATURES)
+Type for an image filter. This takes an EHTImage or any image object
+and creates a filter out of it. The parameters of the image are
+the center of the image `x0`, `y0`.
+
+# Fields
+$(FIELDS)
+"""
 @with_kw struct ImageFilter{T} <: AbstractImageFilter
     x0::Float64
     y0::Float64
@@ -841,3 +911,4 @@ function filter_image(θ::AbstractFilter,
     end
     return (xitr,yitr,img)
 end
+

test/common.jl

@@ -15,6 +15,8 @@ const ξs_1 = 0.25
 const ξs_2 = -0.567
 const x0 = -5.0
 const y0 = 2.5
+const δϕ = 0.1π
+
 
 const npix = 512
 const xlim = [-60.0,60.0]

test/filters.jl

@@ -107,6 +107,18 @@ end
     @test length(fieldnames(GeneralGaussianRing)) == VIDA.size(GeneralGaussianRing)
 end
 
+@testset "FilterLogSpiral" begin
+    θ = LogSpiral(r0, τ, σ, δϕ, ξs, x0, y0)
+    θ1 = LogSpiral(r0=r0,x0=x0,σ=σ,y0=y0,κ=τ,ξ=ξs, δϕ=δϕ)
+    θ2 = LogSpiral([r0, τ, σ, δϕ, ξs, x0, y0])
+
+    @test unpack(θ) == unpack(θ1)
+    @test unpack(θ1) == unpack(θ2)
+    fimg = VIDA.filter_image(θ,npix,xlim,ylim)
+    @test length(fieldnames(LogSpiral)) == VIDA.size(typeof(θ))
+end
+
+
 @testset "FilterImageFilter" begin
     tmp = GaussianRing(r0, σ, x0, y0)
     img = VIDA.make_image(tmp, 60, (-60.0, 60.0), (-60.0, 60.0))