Stencils

mapstencil(f, A::StencilArray, args::AbstractArray...)
mapstencil(f, stencil::Stencil, A::AbstractArray, args::AbstractArray...; kw...)

Stencil mapping where f is passed a Stencil centered at each index in A, followed by the values from args at each stencil center.

Keywords

• boundary: a BoundaryCondition like Wrap. Remove() is the default.
• padding: Padding like Conditional or Halo{:in}. Conditional() is the default.

The result is returned as a new array.

mapstencil!(f, dest::AbstractArray, source::StencilArray, args::AbstractArray...)
mapstencil!(f, A::SwitchingStencilArray, args::AbstractArray...)

Stencil mapping where f is passed a stencil centered at each index in src, followed by the values from args at each stencil center. The result of f is written to dest.

For SwitchingStencilArray the internal source and dest arrays are used, returning a switched version of the array.

dest must either be smaller than src by the stencil radius on all sides, or be the same size, in which case it is assumed to also be padded.

Stencil <: StaticVector

Stencils define a pattern of neighboring cells around the current cell. They reduce the structure and dimensions of the neighborhood into a StaticVector of values.

Stencil objects are updated to contain the neighbors for an array index.

This design is so that user functions can be passed a single object from which they can retrieve center, neighbors, offsets, distances to neighbors and other information.

Stencils also provide a range of compile-time utility funcitons like distances and offsets.

NamedStencil <: AbstractStencil

NamedStencil(; kw...)
NamedStencil(values::NamedTuple)
NamedStencil{Keys}(values)

A named stencil that can take arbitrary shapes where each offset position is named. This can make stencil code easier to read by removing magic numbers.

Example

julia> using Stencils

julia> ns = NamedStencil(; west=(0, -1), north=(1, 0), south=(-1, 0), east=(0, 1)) 
NamedStencil{(:west, :north, :south, :east), ((0, -1), (1, 0), (-1, 0), (0, 1)), 1, 2, 4, Nothing}
▄▀▄
 ▀ 

julia> A = StencilArray((1:10) * (1:10)', ns);

julia> stencil(A, (5, 5)).east # we can access values by name
30

julia> mapstencil(s -> s.east + s.west, A); # and use them in `mapstencil` functions

We can also take some shortcuts, and just name an existing stencil:

julia> ns = NamedStencil{(:w,:n,:s,:e)}(VonNeumann(1)) 

The stencil radius is calculated from the most distant coordinate, and the dimensionality N of the stencil is taken from the length of the first coordinate, e.g. 1, 2 or 3.

AngledCross <: Stencil

AngledCross(; radius=1, ndims=2)
AngledCross(radius, ndims)
AngledCross{R,N}()

A neighboorhood where all diagonals are included, for 2:N dimensions

BackSlash <: Stencil

BackSlash(; radius=1, ndims=2)
BackSlash(radius, ndims)
BackSlash{R,N}()

A neighboorhood along the 'backwards' diagonal. Contains 2R+1 neighbors for for 2:N dimensions

Circle <: Stencil

Circle(; radius=1, ndims=2)
Circle(radius, ndims)
Circle{R,N}()

A circular or spherical stencil

Cross <: Stencil

Cross(; radius=1, ndims=2)
Cross(radius, ndims)
Cross{R,N}()

A cross-shaped neighboorhood where positions with offsets of 0 on least N-1 axes are included

Diamond <: Stencil

Diamond(; radius=1, ndims=2)
Diamond(radius, ndims)
Diamond{R,N}()

A diamond or regular octahedron

ForwardSlash <: Stencil

ForwardSlash(; radius=1, ndims=2)
ForwardSlash(radius, ndims)
ForwardSlash{R,N}()

A neighboorhood where only 'forward' diagonals are included. Contains 2R+1 neighbors for 2:N dimensions

Kernel <: AbstractKernelStencil

Kernel(stencil::Stencil, kernel::AbstractArray)
Kernel(f::Function, stencil::Stencil)

Wrap any other stencil object, and includes a kernel array of the same length and positions as the stencil. A function of the stencil and kernel, like kernelproduct can be used in mapstencil.

A function f may be passed as the first argument, and a kernel array will be calculated with map(f, distances(stencil)).

As an example, Kernels can be convolved with an Array

using Stencils

# Define a random array that the kernel will be convolved with
r = rand(1000, 1000)

# Define kernel array
sharpen = [0 -1 0;
           -1 5 -1;
           0 -1 0]

# Define a stencil that is the same size as the kernel array
stencil = Window(1)

# Create a stencil Kernel from the stencil and the kernel array
k = Kernel(stencil, sharpen)

# Wrap the random array and the Kernel in a StencilArray
A = StencilArray(r, k)

# use `mapstencil` with the `kernelproduct` function to convolve the Kernel with array. 
# Note: `kernelproduce is similar to `Base.dot` but `kernelproduct` 
# lets you use an array of StaticArray and it will still work (dot is recursive).
mapstencil(kernelproduct, A) 

Horizontal <: Stencil

Horizontal(; radius=1, ndims=2)
Horizontal(radius, ndims)
Horizontal{R,N}()

A horizontal bar or plane

Layered <: Abstract

Layered(layers::Union{Stencil,Layered}...)
Layered(; layer_keywords...)
Layered(layers::Union{Tuple,NamedTuple})

Tuple or NamedTuple of stencils that can be used together.

neighbors for Layered returns a tuple of iterators for each stencil layer.

Moore <: Stencil

Moore(; radius=1, ndims=2)
Moore(radius, ndims)
Moore{R,N}()

Moore stencils define the stencil as all cells within a horizontal or vertical distance of the central cell. The central cell is omitted.

Positional <: AbstractPositionalStencil

Positional(offsets::Tuple{Vararg{Int}}...)
Positional(offsets::Tuple{Tuple{Vararg{Int}}})
Positional{O}()

Stencils that can take arbitrary shapes by specifying each coordinate, as Tuple{Int,Int} of the row/column distance (positive and negative) from the central point.

The stencil radius is calculated from the most distant coordinate, and the dimensionality N of the stencil is taken from the length of the first coordinate, e.g. 1, 2 or 3.

See NamedStencil for a similar stencil with named offsets.

Example

julia> p = Positional((0, -1), (2, 1), (-1, 1), (0, 1)) 
Positional{((0, -1), (2, 1), (-1, 1), (0, 1)), 2, 2, 4, Nothing}
   ▄ 
 ▀ ▀ 
   ▀ 

Rectangle <: Stencil

Rectangle(offsets::Tuple{Tuple}...)
Rectangle{O}()

Rectanglar stencils of arbitrary shapes. These are specified with pulles of offsets around the center point, one for each dimension.

Vertical <: Stencil

Vertical(; radius=1, ndims=2)
Vertical(radius, ndims)
Vertical{R,N}()

A vertical bar or plane

VonNeumann <: Stencil

VonNeumann(; radius=1, ndims=2)
VonNeumann(radius, ndims)
VonNeumann{R,N}()

Diamond-shaped neighborhood (in 2 dimensions), without the central cell In 1 dimension it is identical to Moore.

Window <: Stencil

Window(; radius=1, ndims=2)
Window(radius, ndims)
Window{R,N}()

A neighboorhood of radius R that includes the central cell.

Annulus <: Stencil

Annulus(; outerradius=2, innerradius=1, ndims=2) Annulus(outerradius, innerradius, ndims) Annulus{RO,RI,N}()

A donut or hollowed spherical stencil

Cardinal <: Stencil

Cardinal(; radius=1, ndims=2)
Cardinal(radius, ndims)
Cardinal{R,N}()

Cardinal (as in N,S,W,E compass directions) stencil

Ordinal <: Stencil

Ordinal(; radius=1, ndims=2)
Ordinal(radius, ndims)
Ordinal{R,N}()

Ordinal (as in NE,SE,SW,NW wind directions) stencil

stencil(x) -> Stencil

Returns a stencil object.

neighbors(x::Stencil) -> iterable

Returns a basic SVector of all cells in the stencil.

center(x::Stencil)

Return the value of the central cell a stencil is offset around. It may or may not be part of the stencil itself.

offsets(x)

Return an SVector of NTuple{N,Int}, containing all positions in the stencil as offsets from the central cell.

Custom Stencils must define this method.

indices(x::Stencil, I::Union{Tuple,CartesianIndex})
indices(x::AbstractStencilArray, I::Union{Tuple,CartesianIndex})

Returns an SVector of CartesianIndices for each neighbor around I.

indices for Stencil do not know about array boundaries and wil not wrap or reflect. On AbstractStencilArray they will wrap and reflect depending on the boundary condition of the array.

distances(hood::Stencil)

Returns an SVector of center-to-center distance of each stencil position from the central cell, so that horizontally or vertically adjacent cells have a distance of 1.0, and a diagonally adjacent cell has a distance of sqrt(2.0).

Values are calculated at compile time, so distances can be used with little overhead.

radius(stencil) -> Int

Return the radius of a stencil.

diameter(rule) -> Int

The diameter of a stencil is 2r + 1 where r is the radius.

merge(stencils::Stencil...)

Merge multiple stencils into a single Positional stencil.

Dimensionality must be the same for all stencils but radii may differ.

merge(stencils::NamedStencil...)

Merge multiple named stencils into a single NamedStencil.

Dimensionality must be the same for all stencils, and all named stencils must have unique names and offsets.

distance_zones(hood::Stencil)

Returns an SVector of Int distance zones for each offset in the Stencil.

kernel(hood::AbstractKernelStencil) => iterable

Returns the kernel object, an array or iterable matching the length of the stencil.

kernelproduct(hood::AbstractKernelStencil)
kernelproduct(hood::Stencil, kernel)

Returns the vector dot product of the stencil and the kernel, although differing from dot in that it is not taken iteratively for members of the stencil - they are treated as scalars.

unsafe_neighbors([stencil::Stencil,] A::AbstractStencilArray, I::CartesianIndex) => SArray

Get stencil neighbors from A around center I as an SVector, without checking bounds of I.

unsafe_stencil(x, A::AbstractArray, I) => Stencil

Update the neighbors of a stencil to values from the array A around index I, without checking bounds of I. Bounds checking of neighbors still occurs, but with the assumption that I is inbounds.

rebuild(x::Stencil, neighbors::StaticArray)

Rebuild a Stencil, returning an stencil of the same size and shape, with new neighbor values.

getneighbor(A::AbstractStencilArray, I::CartesianIndex)

Get an array value from a stencil neighborhood.

This method handles boundary conditions.

BoundaryCondition

Abstract supertype for flags that specify the boundary conditions. These determine what happens when a stencil extends outside of the grid.

Remove <: BoundaryCondition

Remove()

BoundaryCondition flag that specifies to assign padval to cells that overflow grid boundaries.

Use <: BoundaryCondition

Use()

BoundaryCondition flag that specifies to use the existing padding, which is only possible when Halo{:in} is used for padding.

Wrap <: BoundaryCondition

Wrap()

BoundaryCondition flag to wrap cordinates that boundary boundaries back to the opposite side of the grid.

Reflect <: BoundaryCondition

Reflect()

BoundaryCondition flag to mirror cordinates that boundary boundaries back to the source cell of the grid.

Padding

Abstract supertype for padding modes, e.g. Conditional and Halo.

Conditional <: Padding

Padding that doesn't change the array size, but checks getindex for out-of-bounds indexing, and inserts padval with Remove or values from the other side of the array with Wrap.

Halo{X} <: Padding

Padding that uses an in-memory halo around the array so that parts of a stencil that go off the edge of the array can index directly into it without a bounds check or any conditional. This has the benefit of possibly better performance during window broadcasts, but some downsides.

In :out mode, a whole new array is alocated, larger than the original. This may not be worth doing unless you are using it multiple times. with :in mode, the outside edge of the array is used as padding. This may be more accurate as there are no boundary effects from using a padding value.:w

Example

halo_in = Halo(:in)
halo_out = Halo(:out)

AbstractStencilArray <: StaticArray

Supertype for arrays with a Stencil, a BoundaryCondition, and Padding.

StencilArray <: AbstractStencilArray

StencilArray(A::AbstractArray, stencil::Stencil; kw...)

An array with a Stencil and a BoundaryCondition, and Padding.

For most uses a StencilArray works exactly the same as a regular array.

Except it can be indexed at any point with stencil to return a filled Stencil object, or neighbors to return an SVector of neighbors.

Arguments

• A: an AbstractArray
• stencil: a Stencil.

Keywords

• boundary: a BoundaryCondition like Wrap. Remove() is the default.
• padding: Padding like Conditional or Halo{:in}. Conditional() is the default.

Example

using Stencils, Statistics
sa = StencilArray((1:10) * (10:20)', Moore(1); boundary=Wrap())
sa .*= 2 # Broadcast works as usual
means = mapstencil(mean, sa) # mapstencil works
stencil(sa, 5, 6) # manually reading a stencil works too

# output

Moore{1, 2, 8, Int64}
█▀█
▀▀▀

with neighbors:
8-element StaticArraysCore.SVector{8, Int64} with indices SOneTo(8):
 112
 140
 168
 120
 180
 128
 160
 192

AbstractSwitchingStencilArray

Abstract supertype for AbstractStencilArray that wrap two arrays that switch places with each broadcast.

SwitchingStencilArray <: AbstractSwitchingStencilArray

SwitchingStencilArray(A::AbstractArray, stencil; 
    boundary=Remove(zero(eltype(parent))),
    padding=Conditional(),
)

An AbstractArray with a Stencil, a BoundaryCondition, Padding, and two array layers that are switched with each broadcast_stencil operation.

The use case for this operation is in simulations where stencil operations are repeatedly run over the same data, or where a filter (such as a blur) needs to be applied many times.

For most uses a SwitchingStencilArray works exactly the same as a regular array - the dest array can be safely ignored.

However, when using mapstencil! you need to use the output, not the original array. Switching does not happen in-place, but as a new returned array.

Arguments

• A: an AbstractArray
• stencil: a Stencil.

Keywords

• boundary: a BoundaryCondition like Wrap. Remove() is the default.
• padding: Padding like Conditional or Halo{:in}. Conditional() is the default.

Example

using Stencils, Statistics

sa = SwitchingStencilArray(rand(10, 10), Moore(1); boundary=Wrap())
sa .*= 2 # Broadcast works as usual
mapstencil(mean, sa) # As does runing `mapstencils
s = stencil(sa, 5, 10) # And retreiving a single stencil
# But we can also run it in-place, here doing 10 iterations of mean blur:
# Note: if you dont assign new variable with `A =`, the array will
# not switch and will not be blurred.
let sa = sa
    for i in 1:10
        sa = mapstencil!(mean, sa)
    end
end
# output

boundary(A::AbstractStencilArray)

Get the BoundaryCondition object from an AbstractStencilArray.

padding(A::AbstractStencilArray)

Get the Padding object from an AbstractStencilArray.

switch(A::SwitchingStencilArray)

Swap the source and dest of a SwitchingStencilArray.