Guide
Documentation may not be accurate as this is a beta stage package undergoing changes. Raise Github issue if you have questions :)
Installation
We're a registered Julia package, but it's recommended to install the latest revision directly from Github.
using Pkg; Pkg.add(url="https://github.com/aced-differentiate/EquivariantOperators.jl.git")
Scalar & vector fields
Scalar & vector fields are represented as 2d/3d arrays of canonically scalars or vectors (StaticVectors
from StaticArrays.jl
for performance). This vector field representation is consistent with multi-channel images from Julia Images which differs from representations using separate arrays for field components. Most Images
functions are readily applicable. Array values can alternatively be any custom type that Supports addition & multiplication, such as complex numbers and custom structs encoding spherical harmonics.
Customizable grid, interpolation, particle mesh placement
Grid
— TypeGrid(resolutions, origin=ones(length(resolutions)), sz=nothing)
Grid(resolutions, rmax)
Grid(cell, origin=ones(size(cell, 1)), sz=nothing)
Grid(cell, rmax)
Constructs Grid
for interpolation, coordinate retrieval, and particle mesh placement. Grids can be in any dimension and be boundless or bounded. At minimum, a grid stores its discretization cell and origin in pixels.
For an orthogonal grid, supply a list of resolutions in each dimension. For non-orthogonal grid, supply the cell vectors as a column-wise matrix. Origin by default is at index (1, ...).
Bounded grids compute and store additional info eg position vectors and their norms of all grid points. To construct bounded grid, supply the overall integer pixel size sz
. For bounded cubic grid, use convenience constructor with rmax
Params
resolutions
: list of resolutions in each dimension for orthogonal gridsorigin
: indices of the origin, may be decimal valuedsz
: integer size (pixel dimensions) of bounded gridrmax
: halflength (in length units, not pixels) of bounded cubic gridcell
: column-wise matrix of discretization cell vectors
fields
cell
origin
- additional fields of bounded grids
p
:Array
of position vectorsr
:Array
of radial distances (lengths of position vectors)x
,y
,z
:Array
of that coordinate (field aliases callingp
)
Examples
2d interpolation
dx = dy = 0.1
g = Grid((dx, dy))
a = [x^2 + y^2 for x in 0:dx:0.5, y in 0:dy:0.5]
a[g, 0.25, 0.25]
# 0.13
# alternatively the example array can be made using a bounded grid and its coordinate fields
R = 0.5
g = Grid((dx, dy), R)
a = g.x.^2 + g.y.^2 # or g.r.^2
grid construction details
Grid((0.1,)) # 1d grid spaced 0.1 apart with origin at (1,)
Grid((0.1, 0.2), (2, 5)) # 2d grid spaced 0.1 along x and 0.2 along y, with origin at (2 ,5)
Grid((0.1, 0.1), 20.0) # bounded 2d grid spaced 0.1 apart, halflength 20.0, with origin at (201, 201), pixel size (401, 401)
Grid(0.1 * [1 1 0; 1 0 1; 1 1 1]', ones(3), (10, 12, 15)) # bounded 3d grid with cell vectors [.1, .1, 0], [.1, 0, .1], [.1, .1, .1] (note matrix transpose). origin (1, 1, 1), pixel size (10, 12, 15). can construct lattice this way
place!
— Functionplace!(field, grid, rvec, val)
Place a discrete impulse value via particle mesh method, specifically using area weighting according to the CIC (cloud-in-cell) rule. The impulse is thus spread onto the nearest grid points with proximity weighting and discretization scaling.
particle mesh placement and interpolation
dx = dy = dz = 0.1
g = Grid((dx, dy, dz))
a = zeros(5, 5, 5)
v = 1.0
place!(a, g, (0.2, 0.2, 0.2), v)
a[g, 0.2, 0.2, 0.2]
# 1000
v / g.dv
# 1000
Finite difference equivariant operators
Del
— FunctionDel(resolutions; pad = :same, border = :smooth)
Del(cell; pad = :same, border = :smooth)
constructs gradient operator (also divergence, curl) using central difference stencil. By default, boundaries are smooth (C1 or C2 continuous) and output is of same length as input.
Example
1d derivative
dx = 0.1
x = 0:dx:.5
y = x .^ 2
d = Del((dx,))
d(y)
#=
6-element Vector{Float64}:
0.0
0.2
0.4
0.6
0.8
1.0
=#
2d gradient
dy = dx = 0.1
a = [x^2 + y^2 for x in 0:dx:0.5, y in 0:dy:0.5]
▽ = Del((dx, dy))
grad_a = ▽(a)
#=
6×6 Matrix{SVector{2, Float64}}:
[0.2, 0.2] [0.2, 0.2] [0.2, 0.4] [0.2, 0.6] [0.2, 0.8] [0.2, 0.8]
[0.2, 0.2] [0.2, 0.2] [0.2, 0.4] [0.2, 0.6] [0.2, 0.8] [0.2, 0.8]
[0.4, 0.2] [0.4, 0.2] [0.4, 0.4] [0.4, 0.6] [0.4, 0.8] [0.4, 0.8]
[0.6, 0.2] [0.6, 0.2] [0.6, 0.4] [0.6, 0.6] [0.6, 0.8] [0.6, 0.8]
[0.8, 0.2] [0.8, 0.2] [0.8, 0.4] [0.8, 0.6] [0.8, 0.8] [0.8, 0.8]
[0.8, 0.2] [0.8, 0.2] [0.8, 0.4] [0.8, 0.6] [0.8, 0.8] [0.8, 0.8]
=#
2d divergence
dy = dx = 0.1
a = [[2x,3y] for x in 0:dx:0.5, y in 0:dy:0.5]
▽ = Del((dx, dy))
▽ ⋅ (a)
#=
6×6 Matrix{Float64}:
5.0 5.0 5.0 5.0 5.0 5.0
5.0 5.0 5.0 5.0 5.0 5.0
5.0 5.0 5.0 5.0 5.0 5.0
5.0 5.0 5.0 5.0 5.0 5.0
5.0 5.0 5.0 5.0 5.0 5.0
5.0 5.0 5.0 5.0 5.0 5.0
=#
Lap
— FunctionLap(cell; pad = :same, border = :smooth)
constructs Laplacian operator
Examples
# 2d Example
dy = dx = 0.1
a = [x^2 + y^2 for x in 0:dx:0.5, y in 0:dy:0.5]
▽2 = Lap((dx, dy))
▽2(a)
#=
6×6 Matrix{Float64}:
4.0 4.0 4.0 4.0 4.0 4.0
4.0 4.0 4.0 4.0 4.0 4.0
4.0 4.0 4.0 4.0 4.0 4.0
4.0 4.0 4.0 4.0 4.0 4.0
4.0 4.0 4.0 4.0 4.0 4.0
4.0 4.0 4.0 4.0 4.0 4.0
=#
Gauss
— FunctionGauss(resolutions, σ, rmax; kw...)
Gauss(cell, σ, rmax; kw...)
constructs Gaussian diffusion operator with volume Normalized to 1 wrt grid support
Op
— TypeOp(radfunc, rmax, resolutions; l = 0, rmin = 0., pad = :same)
Op(radfunc, rmax, cell; kw...)
Op(radfunc, grid; kw...)
Op
constructs equivariant finite difference operators & custom Green's functions by specifying the radial function of the impulse response. Prebuilt operators like differential operators (▽
) & common Green's functions can be constructed instead using Del
, Lap
.
Args
l
: rotation order,0
for scalar field,1
for vector field
Example
dx = dy = 0.1
resolutions = (dx, dy)
rmin = 1e-9
rmax = 0.2
ϕ = Op(r -> 1 / r, rmax, resolutions; rmin) # 1/r potential
F = Op(r -> 1 / r^2, rmax, resolutions; rmin, l=1) # 1/r^2 field
g = Grid(resolutions,)
a = zeros(5, 5)
a[3, 3] = 1.0 / g.dv # puts discrete value integrating to 1.0 onto array
ϕ(a)
#=
5×5 Matrix{Float64}:
0.0 0.0 5.0 0.0 0.0
0.0 7.07107 10.0 7.07107 0.0
5.0 10.0 0.0 10.0 5.0
0.0 7.07107 10.0 7.07107 0.0
0.0 0.0 5.0 0.0 0.0
=#
F(a)
#=
5×5 Matrix{SVector{2, Float64}}:
[0.0, 0.0] [0.0, 0.0] [-25.0, 0.0] [0.0, 0.0] [0.0, 0.0]
[0.0, 0.0] [-35.3553, -35.3553] [-100.0, 0.0] [-35.3553, 35.3553] [0.0, 0.0]
[0.0, -25.0] [0.0, -100.0] [0.0, 0.0] [0.0, 100.0] [0.0, 25.0]
[0.0, 0.0] [35.3553, -35.3553] [100.0, 0.0] [35.3553, 35.3553] [0.0, 0.0]
[0.0, 0.0] [0.0, 0.0] [25.0, 0.0] [0.0, 0.0] [0.0, 0.0]
=#
function (m::Op)(x::AbstractArray, )
Convolutions
Operators apply to most use cases but you may also use convolution functions directly. We offer feature rich convolution and cross correlation functions with options for padding, stride, boundary conditions, and custom products (tensor field convolutions). We use DSP.conv
as our backend for scalar field convolutions and our own implementations for convolutions involving vector fields or custom products. FFT implementation is automatically invoked when appropriate.
cvconv
— Functioncvconv(x, f; product = *, stride = 1, pad = 0, alg = nothing)
"convolution" in computer vision for any dimension, same as Cross correlation. Automatically uses FFT for big kernels. For convolution in signal processing , use dspconv
instead.
Args
x
input arrayf
filter array
Keywords
product
product in convolution, eg*
,dot
pad
amount of padding or padding option- any integer number of pixels on each boundary
:same
: adds enough padding so ouoverlapdotut is same size as input:outer
: ouoverlapdotut size issize(x) .+ size(f) .- 1
border
type of padding0
value pixels:replicate
repeats edge valuesperiodic
or:circular
: periodic BC:smooth
continuous derivatives at boundaries useful for differential operators:reflect
reflects interior across boundaries which are not repeated:symmetric
same as:reflect
but with boundaries repeated
alg
specifies convolution algorithmnothing
Automatically chooses fastest algorithm:direct
convolution, scales as O(n^2):fft
Fourier convolution, scales as O(n log(n))
Convolutions in other Julia packages, fewer features but perhaps more optimized for speed in their specific use cases
- ImageFiltering.imfilter. Its docs has excellent mathematical explaination of convolutions and correlation as well as padding/border options
DSP.conv
DSP.xcor
Flux.conv
dspconv
— Functiondspconv(a, f; product = *, pad = :outer, border=0)
Convolution in signal processing. For "convolution" in computer vision, use cvconv instead. Automatically uses FFT for big kernels. By default output size is size(x) .+ size(f) .- 1
. See cvconv
for its keyword options which also apply here