Adding a Custom Template

If you want to add your own template you just need to define a new:

subtype of VIDA.AbstractImageTemplate
implement the interface desribed in VIDA.AbstractImageTemplate

In this tutorial we will implement a custom template from scratch.

using VIDA



struct SlashedExponentialRing{T} <: VLBISkyModels.AbstractImageTemplate
    αouter::T
    αinner::T #standard deviation of the Gaussian
    s::T
end

function VIDA.intensity_point(m::SlashedExponentialRing, p)
    (;X, Y) = p
    (;αinner, αouter, s) = m
    r = hypot(X, Y)
    ϕ = atan(X, -Y)

    n = (1-s*cos(ϕ))

    return n*r^αinner/(1 + r^(αouter + αinner+1))
end

  Activating project at `~/work/VIDA.jl/VIDA.jl/example`

We can also add a convienence constructor

function SlashedExponentialRing(r0, αouter, αinner, s, ξs, x0, y0)
    return modify(SlashedExponentialRing(αouter, αinner, s),
            Stretch(r0, r0), Rotate(ξs), Shift(x0, y0))
end

Main.SlashedExponentialRing

Then you can simply call the same optimizing functions and plotting functions. For example lets create a fake image and fit it

template = SlashedExponentialRing(μas2rad(20.0), 3.0, 4.0, 0.5, π/2, 0.0, 0.0)

ModifiedModel
  base model: Main.SlashedExponentialRing{Float64}(3.0, 4.0, 0.5)
  Modifiers:
    1. Stretch{Float64, Float64, Float64}
    2. Rotate{Float64}
    3. Shift{Float64, Float64, Float64}

VIDA uses ComradeBase and VLBISkyModels interface. Therefore, we can create an image using intensitymap

g = imagepixels(μas2rad(128.0), μas2rad(128.0), 64, 64)
img = intensitymap(template, g)

┌ 64×64 IntensityMap{Float64, 2} ┐
├────────────────────────────────┴─────────────────────────────────────── dims ┐
  ↓ X Sampled{Float64} LinRange{Float64}(-3.0543261909900767e-10, 3.0543261909900767e-10, 64) ForwardOrdered Regular Points,
  → Y Sampled{Float64} LinRange{Float64}(-3.0543261909900767e-10, 3.0543261909900767e-10, 64) ForwardOrdered Regular Points
└──────────────────────────────────────────────────────────────────────────────┘
  ↓ →          -3.05433e-10  -2.95736e-10  …  2.95736e-10  3.05433e-10
 -3.05433e-10   1.64145e-5    1.73373e-5      1.73373e-5   1.64145e-5
 -2.95736e-10   1.76458e-5    1.86754e-5      1.86754e-5   1.76458e-5
 -2.8604e-10    1.89748e-5    2.0124e-5       2.0124e-5    1.89748e-5
 -2.76344e-10   2.04082e-5    2.16915e-5      2.16915e-5   2.04082e-5
  ⋮                                        ⋱               
  2.66648e-10   4.34643e-5    4.69293e-5   …  4.69293e-5   4.34643e-5
  2.76344e-10   4.10121e-5    4.41773e-5      4.41773e-5   4.10121e-5
  2.8604e-10    3.86798e-5    4.15695e-5      4.15695e-5   3.86798e-5
  2.95736e-10   3.64663e-5    3.91032e-5      3.91032e-5   3.64663e-5
  3.05433e-10   3.43693e-5    3.67748e-5      3.67748e-5   3.43693e-5

We can also plot the image

using CairoMakie
imageviz(img)

Now lets see if we can get the correct parameters. First we construct our optimization divergence the Bhattacharyya divergence.

bh = Bhattacharyya(img);

To fit we need to define a fitting function. For this our template function needs to accept a named tuple.

temp(θ) = SlashedExponentialRing(θ.r0, θ.αout, θ.αin, θ.s, θ.ξs, θ.x0, θ.y0)

temp (generic function with 1 method)

Additionally we need to define the search region for our template extraction

upper = (r0=μas2rad(40.0), αout=10.0, αin = 10.0, s=0.999, ξs=1π, x0= μas2rad(60.0), y0= μas2rad(60.0))
lower = (r0=μas2rad(5.0), αout=1.0, αin = 0.0, s=0.001, ξs=-1π, x0= -μas2rad(60.0), y0= -μas2rad(60.0))

(r0 = 2.42406840554768e-11, αout = 1.0, αin = 0.0, s = 0.001, ξs = -3.141592653589793, x0 = -2.908882086657216e-10, y0 = -2.908882086657216e-10)

We can now create our problem

prob = VIDAProblem(bh, temp, lower, upper)

VIDAProblem{VIDA.InplaceDivergence{Bhattacharyya, SpatialIntensityMap{Float64, RectiGrid{Tuple{X{DimensionalData.Dimensions.Lookups.Sampled{Float64, LinRange{Float64, Int64}, DimensionalData.Dimensions.Lookups.ForwardOrdered, DimensionalData.Dimensions.Lookups.Regular{Float64}, DimensionalData.Dimensions.Lookups.Points, DimensionalData.Dimensions.Lookups.NoMetadata}}, Y{DimensionalData.Dimensions.Lookups.Sampled{Float64, LinRange{Float64, Int64}, DimensionalData.Dimensions.Lookups.ForwardOrdered, DimensionalData.Dimensions.Lookups.Regular{Float64}, DimensionalData.Dimensions.Lookups.Points, DimensionalData.Dimensions.Lookups.NoMetadata}}}, Serial, DimensionalData.Dimensions.Lookups.NoMetadata, Float64}, Matrix{Float64}, Tuple{X{DimensionalData.Dimensions.Lookups.Sampled{Float64, LinRange{Float64, Int64}, DimensionalData.Dimensions.Lookups.ForwardOrdered, DimensionalData.Dimensions.Lookups.Regular{Float64}, DimensionalData.Dimensions.Lookups.Points, DimensionalData.Dimensions.Lookups.NoMetadata}}, Y{DimensionalData.Dimensions.Lookups.Sampled{Float64, LinRange{Float64, Int64}, DimensionalData.Dimensions.Lookups.ForwardOrdered, DimensionalData.Dimensions.Lookups.Regular{Float64}, DimensionalData.Dimensions.Lookups.Points, DimensionalData.Dimensions.Lookups.NoMetadata}}}, Tuple{}, Symbol}}, typeof(Main.temp), Nothing, @NamedTuple{r0::Float64, αout::Float64, αin::Float64, s::Float64, ξs::Float64, x0::Float64, y0::Float64}}(VIDA.InplaceDivergence{Bhattacharyya, SpatialIntensityMap{Float64, RectiGrid{Tuple{X{DimensionalData.Dimensions.Lookups.Sampled{Float64, LinRange{Float64, Int64}, DimensionalData.Dimensions.Lookups.ForwardOrdered, DimensionalData.Dimensions.Lookups.Regular{Float64}, DimensionalData.Dimensions.Lookups.Points, DimensionalData.Dimensions.Lookups.NoMetadata}}, Y{DimensionalData.Dimensions.Lookups.Sampled{Float64, LinRange{Float64, Int64}, DimensionalData.Dimensions.Lookups.ForwardOrdered, DimensionalData.Dimensions.Lookups.Regular{Float64}, DimensionalData.Dimensions.Lookups.Points, DimensionalData.Dimensions.Lookups.NoMetadata}}}, Serial, DimensionalData.Dimensions.Lookups.NoMetadata, Float64}, Matrix{Float64}, Tuple{X{DimensionalData.Dimensions.Lookups.Sampled{Float64, LinRange{Float64, Int64}, DimensionalData.Dimensions.Lookups.ForwardOrdered, DimensionalData.Dimensions.Lookups.Regular{Float64}, DimensionalData.Dimensions.Lookups.Points, DimensionalData.Dimensions.Lookups.NoMetadata}}, Y{DimensionalData.Dimensions.Lookups.Sampled{Float64, LinRange{Float64, Int64}, DimensionalData.Dimensions.Lookups.ForwardOrdered, DimensionalData.Dimensions.Lookups.Regular{Float64}, DimensionalData.Dimensions.Lookups.Points, DimensionalData.Dimensions.Lookups.NoMetadata}}}, Tuple{}, Symbol}}(Bhattacharyya(), [5.068628969357313e-6 5.353567410504405e-6 … 5.353567410504405e-6 5.068628969357313e-6; 5.448839085356065e-6 5.76676872803105e-6 … 5.76676872803105e-6 5.448839085356065e-6; … ; 1.1260411251307359e-5 1.2074669819611665e-5 … 1.2074669819611665e-5 1.1260411251307359e-5; 1.0612879435517843e-5 1.135568292615902e-5 … 1.135568292615902e-5 1.0612879435517843e-5], [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0]), Main.temp, nothing, (r0 = 2.42406840554768e-11, αout = 1.0, αin = 0.0, s = 0.001, ξs = -3.141592653589793, x0 = -2.908882086657216e-10, y0 = -2.908882086657216e-10), (r0 = 1.939254724438144e-10, αout = 10.0, αin = 10.0, s = 0.999, ξs = 3.141592653589793, x0 = 2.908882086657216e-10, y0 = 2.908882086657216e-10))

The vida method can use any optimizer that works with Optimization.jl For this work we will use BlackBoxOptim.jl.

using OptimizationBBO
xopt, θopt, divmin = vida(prob, BBO_adaptive_de_rand_1_bin(); maxiters=50_000)

@show θopt

ModifiedModel
  base model: Main.SlashedExponentialRing{Float64}(2.9999999753751374, 4.000000111851719, 0.5000000117680279)
  Modifiers:
    1. Stretch{Float64, Float64, Float64}
    2. Rotate{Float64}
    3. Shift{Float64, Float64, Float64}

Let's also plot the results

fig = triptic(img, θopt)
fig

Now with all of this said this template actually already exists in VIDA using the flexible VIDA.RingTemplate.

rad = RadialDblPower(xopt.αin, xopt.αout)
azi = AzimuthalCosine(xopt.s, xopt.ξs)
t   = modify(RingTemplate(rad, azi), Stretch(xopt.r0), Shift(xopt.x0, xopt.y0))

ModifiedModel
  base model: RingTemplate{RadialDblPower{Float64, Float64}, AzimuthalCosine{1, Float64, Float64}}(RadialDblPower{Float64, Float64}(4.000000111851719, 2.9999999753751374), AzimuthalCosine{1, Float64, Float64}((0.5000000117680279,), (1.570796363695795,)))
  Modifiers:
    1. Stretch{Float64, Float64, Float64}
    2. Shift{Float64, Float64, Float64}

fig = triptic(img, t);
fig

This page was generated using Literate.jl.