VLBIImagePriors
Documentation for VLBIImagePriors.
VLBIImagePriors.AdditiveLR
VLBIImagePriors.AngleTransform
VLBIImagePriors.CenterImage
VLBIImagePriors.CenteredLR
VLBIImagePriors.CenteredRegularizer
VLBIImagePriors.DiagonalVonMises
VLBIImagePriors.GaussMarkovRandomField
VLBIImagePriors.GaussMarkovRandomField
VLBIImagePriors.GaussMarkovRandomField
VLBIImagePriors.GaussMarkovRandomField
VLBIImagePriors.ImageDirichlet
VLBIImagePriors.ImageSimplex
VLBIImagePriors.ImageSphericalUniform
VLBIImagePriors.ImageUniform
VLBIImagePriors.MarkovRandomFieldCache
VLBIImagePriors.MarkovRandomFieldCache
VLBIImagePriors.MarkovRandomFieldCache
VLBIImagePriors.NamedDist
VLBIImagePriors.SphericalUnitVector
VLBIImagePriors.TDistMarkovRandomField
VLBIImagePriors.TDistMarkovRandomField
VLBIImagePriors.TDistMarkovRandomField
VLBIImagePriors.WrappedUniform
VLBIImagePriors.alr!
VLBIImagePriors.alrinv!
VLBIImagePriors.clr!
VLBIImagePriors.clrinv!
VLBIImagePriors.to_real
VLBIImagePriors.to_simplex
VLBIImagePriors.to_simplex!
VLBIImagePriors.AdditiveLR
— TypeAdditiveLR <: LogRatioTransform
Defines the additive log-ratio transform. The clr
transformation moves from the simplex Sⁿ → $R^{n-1}$ and is given by
alr(x) = [log(x₁/xₙ) ... log(xₙ/xₙ)],
where g(x) = (∏xᵢ)ⁿ⁻¹
is the geometric mean. The inverse transformation is given by
alr⁻¹(x) = exp.(x)./(1 + sum(x[1:n-1])).
VLBIImagePriors.AngleTransform
— TypeAngleTransform
A transformation that moves two vector x
and y
to an angle θ
. Note that is x
and y
are normally distributed then the resulting distribution in θ
is uniform on the circle.
VLBIImagePriors.CenterImage
— MethodCenterImage(X, Y)
CenterImage(img::IntensityMap)
CenterImage(grid::ComradeBase.AbstractDims)
Constructs a projections operator that will take an arbritrary image and return a transformed version whose centroid is at the origin.
Arguments
X
,Y
: the iterators for the image pixel locationsimg
: AIntensityMap
whose grid is the same as the image you want to centergrid
: The grid that the image is defined on.
Example
julia> using ComradeBase: centroid, imagepixels
julia> grid = imagepixels(μas2rad(100.0), μas2rad(100.0), 256, 256)
julia> K = CenterImage(grid)
julia> img = IntensityMap(rand(256, 256), grid)
julia> centroid(img)
(1.34534e-10, 5.23423e-11)
julia> cimg = K(img)
julia> centroid(cimg)
(1.231e-14, -2.43e-14)
Note
This center image works using a linear projection operator. This means that is does not necessarily preserve the image total flux and the positive definiteness of the image. In practise we find that the deviation from the original flux, and the amount of negative flux is small.
Explanation
The centroid is constructed by
X ⋅ I = 0
Y ⋅ I = 0
where I
is the flattened image of length N
, and X
and Y
are the image pixel locations. This can be simplified into the matrix equation
XY ⋅ I = 0
where XY
is the 2×N
matrix whose first for is given by X
and second is Y
. The space of centered images is then given by the null space of XY
. Given this we can form a matrix C
which is the kernel matrix of the nullspace whose columns are the orthogonal basis vectors of the null space of XY
. Using this we can construct a centered image projection opertor by
K = CC'.
Our centered image is then given by I₀ = KI
.
VLBIImagePriors.CenteredLR
— TypeCenteredLR <: LogRatioTransform
Defines the centered log-ratio transform. The clr
transformation moves from the simplex Sⁿ → Rⁿ and is given by
clr(x) = [log(x₁/g(x)) ... log(xₙ/g(x))]
where g(x) = (∏xᵢ)ⁿ⁻¹
is the geometric mean. The inverse transformation is given by the softmax function and is only defined on a subset of the domain otherwise it is not injective
clr⁻¹(x) = exp.(x)./sum(exp, x).
Notes
As mentioned above this transformation is bijective on the entire codomain of the function. However, unlike the additive log-ratio transform it does not treat any pixel as being special.
VLBIImagePriors.CenteredRegularizer
— TypeCenteredRegularizer(x, y, σ, p)
Regularizes a general image prior p
such that the center of light is close the the origin of the imag. After regularization the log density of the prior is modified to
\[ \log p(I) \to \log p(I) - \frac{(x_C^2 + y_C^2)^2}{2\sigma\^2} N_x N_y\]
where N_x
and N_y
are the number of pixels in the x
and y
direction of the image, and $x_C, y_C$ are the center of light of the image I
.
VLBIImagePriors.DiagonalVonMises
— TypeDiagonalVonMises(μ::Real, κ::Real)
DiagonalVonMises(μ::AbstractVector{<:Real}, κ::AbstractVector{<:Real})
Constructs a Von Mises distribution, with mean μ
and concentraion parameter κ
. If μ
and κ
are vectors then this constructs a independent multivariate Von Mises distribution.
Notes
This is a custom implementation since the version in Distributions.jl
has certain properties that do not play well (having an support only between [-π+μ, π+μ]) with usual VLBI problems. Additionally this distribution has a special overloaded product_distribution
method so that concatenating multiple DiagonalVonMises
together preserves the type. This is helpful for Zygote
autodiff.
VLBIImagePriors.GaussMarkovRandomField
— Typestruct GaussMarkovRandomField{T<:Number, C} <: VLBIImagePriors.MarkovRandomField
A image prior based off of the first-order zero mean Gaussian Markov random field. This prior is similar to the combination of total squared variation TSV and L₂ norm, and is given by
(π²Σ)⁻¹ TSV(I) + (ρ^2π²Σ)L₂(I) + lognorm(ρ, Σ)
where ρ and Σ are the correlation length and variance of the random field and lognorm(ρ,Σ)
is the log-normalization of the random field. This normalization is needed to jointly infer I
and the hyperparameters ρ, Σ.
Fields
ρ
The correlation length of the random field.
Σ
The variance of the random field
cache
The Markov Random Field cache used to define the specific Markov random field class used.
Example
julia> ρ, Σ = 2.0, 1.0
julia> d = GaussMarkovRandomField(ρ, Σ, (32, 32))
julia> cache = MarkovRandomFieldCache(Float64, (32, 32)) # now instead construct the cache
julia> d2 = GaussMarkovRandomField(ρ, Σ, cache)
julia> invcov(d) ≈ invcov(d2)
true
VLBIImagePriors.GaussMarkovRandomField
— MethodGaussMarkovRandomField(ρ, Σ, img::AbstractArray)
Constructs a first order zero-mean Gaussian Markov random field with dimensions size(img)
, correlation ρ
and diagonal covariance Σ
.
VLBIImagePriors.GaussMarkovRandomField
— MethodGaussMarkovRandomField(ρ, Σ, cache::MarkovRandomFieldCache)
Constructs a first order Gaussian Markov random field using the precomputed cache cache
. ρ
and diagonal covariance Σ
and the precomputed MarkovRandomFieldCache cache
. If mean
is not included then it is assume the mean is identically zero.
VLBIImagePriors.GaussMarkovRandomField
— MethodGaussMarkovRandomField(ρ, Σ, dims)
Constructs a first order zero-mean Gaussian Markov random field with dimensions dims
, correlation ρ
and diagonal covariance Σ
.
VLBIImagePriors.ImageDirichlet
— TypeImageDirichlet(α::AbstractMatrix)
ImageDirichlet(α::Real, ny, nx)
A Dirichlet distribution defined on a matrix. Samples from this produce matrices whose elements sum to unity. This is a useful image prior when you want to separately constrain the flux. The α parameter defines the usual Dirichlet concentration amount.
Notes
Much of this code was taken from Distributions.jl and it's Dirichlet distribution. However, some changes were made to make it faster. Additionally, we use define a custom rrule
to speed up derivatives.
VLBIImagePriors.ImageSimplex
— TypeImageSimplex(ny,nx)
This defines a transformation from ℝⁿ⁻¹ to the n
probability simplex defined on an matrix with dimension ny×nx
. This is a more natural transformation for rasterized images, which are most naturally represented as a matrix.
Notes
Much of this code was inspired by TransformVariables. However, we have specified custom rrules
using Enzyme as a backend. This allowed the simplex transform to be used with Zygote and we achieved an order of magnitude speedup when computing the pullback of the simplex transform.
VLBIImagePriors.ImageSphericalUniform
— TypeImageSphericalUniform(nx, ny)
Construct a distribution where each image pixel is a 3-sphere uniform variable. This is useful for polarization where the stokes parameters are parameterized on the 3-sphere.
Currently we use a struct of vectors memory layout. That is the image is described by three matrices (X,Y,Z)
grouped together as a tuple, where each matrix is one direction on the sphere, and we require norm((X,Y,Z)) == 1
.
VLBIImagePriors.ImageUniform
— TypeImageUniform(a::Real, b::Real, nx, ny)
A uniform distribution in image pixels where a/b
are the lower/upper bound for the interval. This then concatenates ny×nx uniform distributions together.
VLBIImagePriors.MarkovRandomFieldCache
— Typestruct MarkovRandomFieldCache{A, TD, M}
A cache for a Markov random field.
Fields
Λ
Intrinsic Gaussian Random Field pseudo-precison matrix. This is similar to the TSV regularizer
D
Gaussian Random Field diagonal precision matrix. This is similar to the L2-regularizer
λQ
The eigenvalues of the Λ matrix which is needed to compute the log-normalization constant of the GMRF.
Examples
julia> mean = zeros(2,2) # define a zero mean
julia> cache = GRMFCache(mean)
julia> prior_map(x) = GaussMarkovRandomField(mean, x[1], x[2], cache)
julia> d = HierarchicalPrior(prior_map, product_distribution([Uniform(-5.0,5.0), Uniform(0.0, 10.0)]))
julia> x0 = rand(d)
VLBIImagePriors.MarkovRandomFieldCache
— MethodMarkovRandomFieldCache(grid::AbstractDims)
Create a GMRF cache out of a Comrade
model as the mean image.
Arguments
m
: Comrade model used for the mean image.grid
: The grid of points you want to create the model image on.
Keyword arguments
transform = identity
: A transform to apply to the image when creating the mean image. See the examples
Example
julia> m = MarkovRandomFieldCache(imagepixels(10.0, 10.0, 64, 64))
VLBIImagePriors.MarkovRandomFieldCache
— MethodMarkovRandomFieldCache(mean::AbstractMatrix)
Contructs the MarkovRandomFieldCache
from the mean image mean
. This is useful for hierarchical priors where you change the hyperparameters of the GaussMarkovRandomField
, ρ and Σ
.
VLBIImagePriors.NamedDist
— MethodNamedDist(d::NamedTuple{N})
NamedDist(;dists...)
A Distribution with names N
. This is useful to construct a set of random variables with a set of names attached to them.
julia> d = NamedDist((a=Normal(), b = Uniform(), c = MvNormal(randn(2), rand(2))))
julia> x = rand(d)
(a = 0.13789342, b = 0.2347895, c = [2.023984392, -3.09023840923])
julia> logpdf(d, x)
Note that NamedDist values passed to NamedDist can also be abstract collections of distributions as well
julia> d = NamedDist(a = Normal(),
b = MvNormal(ones(2)),
c = (Uniform(), InverseGamma())
d = (a = Normal(), Beta)
)
How this is done internally is considered an implementation detail and is not part of the public interface.
VLBIImagePriors.SphericalUnitVector
— TypeSphericalUnitVector{N}()
A transformation from a set of N+1
vectors to the N
sphere. The set of N+1
vectors are inherently assumed to be N+1
a distributed according to a unit multivariate Normal distribution.
Notes
For more information about this transformation see the Stan manual. In the future this may be depricated when is merged.
VLBIImagePriors.TDistMarkovRandomField
— Typestruct TDistMarkovRandomField{T<:Number, C} <: VLBIImagePriors.MarkovRandomField
A image prior based off of the first-order Multivariate T distribution Markov random field.
Fields
ρ
The correlation length of the random field.
Σ
The variance of the random field
ν
The student T "degrees of freedom parameter which ≥ 1 for a proper prior
cache
The Markov Random Field cache used to define the specific Markov random field class used.
Examples
julia> ρ, Σ = 2.0, 1.0
julia> d = TDistMarkovRandomField(ρ, Σ, (32, 32))
julia> cache = MarkovRandomFieldCache(Float64, (32, 32)) # now instead construct the cache
julia> d2 = TDistMarkovRandomField(ρ, Σ, cache)
julia> invcov(d) ≈ invcov(d2)
true
VLBIImagePriors.TDistMarkovRandomField
— MethodTDistMarkovRandomField(ρ, Σ, img::AbstractArray)
Constructs a first order TDist Markov random field with mean image mean
and correlation ρ
and diagonal covariance Σ
.
VLBIImagePriors.TDistMarkovRandomField
— MethodTDistMarkovRandomField(ρ, Σ, cache::MarkovRandomFieldCache)
Constructs a first order TDist Markov random field with zero mean ,correlation ρ
, diagonal covariance Σ
, and the precomputed MarkovRandomFieldCache cache
.
VLBIImagePriors.WrappedUniform
— TypeWrappedUniform(period)
Constructs a potentially multivariate uniform distribution that is wrapped a given period
. That is
d = WrappedUniform(period)
logpdf(d, x) ≈ logpdf(d, x+period)
for any x
.
If period
is a vector this creates a multivariate independent wrapped uniform distribution.
VLBIImagePriors.alr!
— Methodalr!(x, y)
Compute the inverse alr transform. That is x
lives in ℜⁿ and y
, lives in Δⁿ
VLBIImagePriors.alrinv!
— Methodalrinv!(x, y)
Computes the additive logit transform inplace. This converts from ℜⁿ → Δⁿ where Δⁿ is the n-simplex
Notes
This function is mainly to transform the GaussMarkovRandomField to live on the simplex. In order to preserve the nice properties of the GRMF namely the log det we only use y[begin:end-1]
elements and the last one is not included in the transform. This shouldn't be a problem since the additional parameter is just a dummy in that case and the Gaussian prior should ensure it is easy to sample from.
VLBIImagePriors.clr!
— Methodclr!(x, y)
Compute the inverse alr transform. That is x
lives in ℜⁿ and y
, lives in Δⁿ
VLBIImagePriors.clrinv!
— Methodclrinv!(x, y)
Computes the additive logit transform inplace. This converts from ℜⁿ → Δⁿ where Δⁿ is the n-simplex
Notes
This function is mainly to transform the GaussMarkovRandomField to live on the simplex.
VLBIImagePriors.to_real
— Methodto_real(t::LogRatioTransform, y)
Transform the value u
that lives on the simplex to a value in the real vector space. See subtypes(LogRatioTransform)
for a list of possible transformations.
The inverse of this transform is given by to_simplex(t, x)
.
Example
```julia julia> y = randn(100) julia> y .= y./sum(y) julia> toreal(CenteredLR(), y) julia> toreal(AdditiveLR(), y)
VLBIImagePriors.to_simplex!
— Methodto_simplex!(t::LogRatioTransform, y, x)
Transform the vector x
assumed to be a real valued array to the simplex using the log-ratio transform t
and stores the value in y
.
The inverse of this transform is given by to_real!(t, x, y)
where y
is a vector that sums to unity, i.e. it lives on the simplex.
Example
julia> x = randn(100)
julia> to_simplex(CenteredLR(), x)
julia> to_simplex(AdditiveLR(), x)
VLBIImagePriors.to_simplex
— Methodto_simplex(t::LogRatioTransform, x)
Transform the vector x
assumed to be a real valued array to the simplex using the log-ratio transform t
. See subtypes(LogRatioTransform)
for a list of possible transformations.
The inverse of this transform is given by to_real(t, y)
where y
is a vector that sums to unity, i.e. it lives on the simplex.
Example
julia> x = randn(100)
julia> to_simplex(CenteredLR(), x)
julia> to_simplex(AdditiveLR(), x)