Grid(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 grids
• origin: indices of the origin, may be decimal valued
• sz: integer size (pixel dimensions) of bounded grid
• rmax: halflength (in length units, not pixels) of bounded cubic grid
• cell: column-wise matrix of discretization cell vectors

fields

• cell
• origin
• additional fields of bounded grids
  • p: Array of position vectors
  • r: Array of radial distances (lengths of position vectors)
  • x, y, z: Array of that coordinate (field aliases calling p)

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

coordinate retrieval

# alternatively the previous 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

function (m::Op)(x::AbstractArray, )

Op(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

• radfunc: radial function

Keywords

• 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]
=#

Del(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.0, 0.0]  [0.0, 0.2]  [0.2, 0.4]  [0.2, 0.6]  [0.2, 0.8]  [0.2, 0.8]
 [0.2, 0.0]  [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, 1.0]
 [0.8, 0.2]  [0.8, 0.2]  [0.8, 0.4]  [0.8, 0.6]  [1.0, 0.8]  [1.0, 1.0]
=#

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
=#

Gauss(resolutions, σ, rmax; kw...)
Gauss(cell, σ, rmax; kw...)

constructs Gaussian diffusion operator with volume Normalized to 1 wrt grid support

Lap(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
=#

cvconv(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 array
• f 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 is size(x) .+ size(f) .- 1
• border type of padding
  • 0 value pixels
  • :replicate repeats edge values
  • periodic 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 algorithm
  • nothing 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(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

fftconv(x, y; product=*)

Signal processing convolution via FFT in any dimension.

function nae(yhat, y; s = sum(abs.(y)))

Normalized absolute error. eg 0.1 denotes 10% absolute error

place!(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