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.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}(3.0000000206172914, 4.0000000621366265, 0.50000000442367)
  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.0000000621366265, 3.0000000206172914), AzimuthalCosine{1, Float64, Float64}((0.50000000442367,), (1.5707963379270948,)))
  Modifiers:
    1. Stretch{Float64, Float64, Float64}
    2. Shift{Float64, Float64, Float64}

fig = triptic(img, t);
fig

This page was generated using Literate.jl.