From 6f785c71f56e8cf93e54b3c96d95d86e0e587a7b Mon Sep 17 00:00:00 2001 From: Paul Date: Tue, 22 Dec 2020 13:28:15 -0500 Subject: [PATCH 1/2] Created log spiral branch --- src/filters.jl | 53 ++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 53 insertions(+) diff --git a/src/filters.jl b/src/filters.jl index 15ee250..cdb62a8 100644 --- a/src/filters.jl +++ b/src/filters.jl @@ -38,6 +38,59 @@ function _update(θ::AbstractFilter, p) return typeof(θ)(p) end + +@doc """ + $(SIGNATURES) +Filter type for a logarithmic spiral that also has a +Gaussian falloff for the inner and outer radius of the arms +""" +@with_kw struct LogSpiral{T<:Real} <: AbstractImageFilter + """ Unit curvature of the logarithmic spiral """ + κ::T + """ Radius of the logarithmic spiral """ + r0::T + """ width of the Gaussian spiral arm """ + σ::T + """ peak brightness location """ + ϕmax::T + """ angular extent of arc """ + σϕ::T + """ azimuthal angle of the Gaussian spiral arm """ + ξ::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) == 8 + LogSpiral{T}(p[1],p[2],p[3],p[4],p[5],p[6],p[7],p[8]) +end +size(::Type{LogSpiral}) = 8 +@inline function (θ::LogSpiral)(x,y) + @unpack κ, r0, σ, ϕmax, σϕ, ξ, x0, y0 = θ + r = hypot(x-x0,y-y0) + α = atan(y-y0,x-x0) + x′,y′ = rotate(x-x0,y-y0,ξ) + k = sqrt(1-κ*κ)/κ + n = (log(r/r0)/k - α)/(2π) + nc = ceil(n) + nf = floor(n) + rc = r0*exp(k*(α + nc*2π)) + rf = r0*exp(k*(α + nf*2π)) + dist = min(abs(rc-r),abs(rf-r)) + return exp(-dist^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 From 3e0cdba47a9878fd6fadfb0caf686f67bcb4e635 Mon Sep 17 00:00:00 2001 From: Paul Tiede Date: Wed, 23 Dec 2020 10:55:25 -0500 Subject: [PATCH 2/2] Logspiral finished --- src/VIDA.jl | 4 ++-- src/filters.jl | 60 ++++++++++++++++++++++++++++++++----------------- test/common.jl | 2 ++ test/filters.jl | 12 ++++++++++ 4 files changed, 55 insertions(+), 23 deletions(-) diff --git a/src/VIDA.jl b/src/VIDA.jl index 8792fcf..d7fc847 100644 --- a/src/VIDA.jl +++ b/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 diff --git a/src/filters.jl b/src/filters.jl index cdb62a8..bd4dcb2 100644 --- a/src/filters.jl +++ b/src/filters.jl @@ -41,21 +41,21 @@ end @doc """ $(SIGNATURES) -Filter type for a logarithmic spiral that also has a -Gaussian falloff for the inner and outer radius of the arms +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 - """ Radius of the logarithmic spiral """ - r0::T - """ width of the Gaussian spiral arm """ + """ thickness of the Gaussian spiral arm """ σ::T + """ Azimuthal extent of the spiral arm """ + δϕ::T """ peak brightness location """ - ϕmax::T - """ angular extent of arc """ - σϕ::T - """ azimuthal angle of the Gaussian spiral arm """ ξ::T """ x location of disk center in μas """ x0::T @@ -63,23 +63,40 @@ Gaussian falloff for the inner and outer radius of the arms y0::T end function LogSpiral(p::Vector{T}) where {T<:Real} - @assert length(p) == 8 - LogSpiral{T}(p[1],p[2],p[3],p[4],p[5],p[6],p[7],p[8]) + @assert length(p) == 7 + LogSpiral{T}(p[1],p[2],p[3],p[4],p[5],p[6],p[7]) end -size(::Type{LogSpiral}) = 8 +size(::Type{LogSpiral{T}}) where {T} = 7 + + @inline function (θ::LogSpiral)(x,y) - @unpack κ, r0, σ, ϕmax, σϕ, ξ, x0, y0 = θ - r = hypot(x-x0,y-y0) - α = atan(y-y0,x-x0) - x′,y′ = rotate(x-x0,y-y0,ξ) + @unpack κ, σ, r0, δϕ, ξ, x0, y0 = θ + x′,y′ = x-x0,y-y0 + #Set up the spiral k = sqrt(1-κ*κ)/κ - n = (log(r/r0)/k - α)/(2π) + 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 = r0*exp(k*(α + nc*2π)) - rf = r0*exp(k*(α + nf*2π)) - dist = min(abs(rc-r),abs(rf-r)) - return exp(-dist^2/(2*σ^2)) + 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 """ @@ -894,3 +911,4 @@ function filter_image(θ::AbstractFilter, end return (xitr,yitr,img) end + \ No newline at end of file diff --git a/test/common.jl b/test/common.jl index b7d59eb..9043ebc 100644 --- a/test/common.jl +++ b/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] diff --git a/test/filters.jl b/test/filters.jl index ea2f3a8..aa2af3f 100644 --- a/test/filters.jl +++ b/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))