Документация Engee

Reference

Страница в процессе перевода.

Structuring element

strel([T], X::AbstractArray)

Convert structuring element (SE) X to appropriate presentation format with element type T. This is a useful tool to generate SE that most ImageMorphology functions understand.

ImageMorphology currently supports two commonly used representations:

  • T=CartesianIndex: offsets to its center point. The output type is Vector{CartesianIndex{N}}.

  • T=Bool: connectivity mask where true indicates connected to its center point. The output type is BitArray{N}.

julia> se_mask = centered(Bool[1 1 0; 1 1 0; 0 0 0]) # connectivity mask
3×3 OffsetArray(::Matrix{Bool}, -1:1, -1:1) with eltype Bool with indices -1:1×-1:1:
 1  1  0
 1  1  0
 0  0  0

julia> se_offsets = strel(CartesianIndex, se_mask) # displacement offsets to its center point
3-element Vector{CartesianIndex{2}}:
 CartesianIndex(-1, -1)
 CartesianIndex(0, -1)
 CartesianIndex(-1, 0)

julia> se = strel(Bool, se_offsets)
3×3 OffsetArray(::BitMatrix, -1:1, -1:1) with eltype Bool with indices -1:1×-1:1:
 1  1  0
 1  1  0
 0  0  0

See also strel_diamond and strel_box for SE constructors for two special cases.

strel_chain(A, B, ...)
strel_chain(As)

For structuring elements of the same dimensions, chain them together to build a bigger one.

The output dimension is the same as the inputs dimensions. See also strel_product that cartesian producting each SE.

For some morphological operations f such as dilation and erosion, if se can be decomposed into smaller ones se1, se2, ..., seN, then f(img, se) is equivalent to f(...f(f(img, se1), se2), ..., seN). Because applying f to smaller SEs is more efficient than to the original big one, this trick is used widely in image morphology.

 
julia> img = rand(512, 512);

julia> se1, se2 = [centered(rand(Bool, 3, 3)) for _ in 1:2];

julia> se = strel_chain(se1, se2);

julia> out_se = dilate(img, se);

julia> out_pipe = dilate(dilate(img, se1), se2);

julia> out_se[2:end-1, 2:end-1] == out_pipe[2:end-1, 2:end-1] # except for the boundary
true
strel_product(A, B, ...)
strel_product(se_list)

Cartesian product of multiple structuring elements; the output dimension ndims(out) == sum(ndims, se_list).

See also strel_chain that chains SEs in the same dimension.

julia> strel_product(strel_diamond((5, 5)), centered(Bool[1, 1, 1]))
5×5×3 SEChainArray{3, OffsetArrays.OffsetArray{Bool, 3, BitArray{3}}} with indices -2:2×-2:2×-1:1:
[:, :, -1] =
 0  0  1  0  0
 0  1  1  1  0
 1  1  1  1  1
 0  1  1  1  0
 0  0  1  0  0

[:, :, 0] =
 0  0  1  0  0
 0  1  1  1  0
 1  1  1  1  1
 0  1  1  1  0
 0  0  1  0  0

[:, :, 1] =
 0  0  1  0  0
 0  1  1  1  0
 1  1  1  1  1
 0  1  1  1  0
 0  0  1  0  0
strel_box(A; r=1)
strel_box(size; r=size .÷ 2)

Construct the N-dimensional structuring element (SE) with all elements in the local window connected.

If image A is provided, then the SE size will be (2r+1, 2r+1, ...) with default half-size r=1. If size is provided, the default r will be size .÷ 2. The default dims will be all dimensions, that is, (1, 2, ..., length(size)).

julia> img = rand(64, 64);

julia> strel_box(img)
3×3 SEBoxArray{2, UnitRange{Int64}} with indices -1:1×-1:1:
 1  1  1
 1  1  1
 1  1  1

julia> strel_box(img; r=2)
5×5 SEBoxArray{2, UnitRange{Int64}} with indices -2:2×-2:2:
 1  1  1  1  1
 1  1  1  1  1
 1  1  1  1  1
 1  1  1  1  1
 1  1  1  1  1

julia> strel_box((5,5); r=(1,2))
5×5 SEBoxArray{2, UnitRange{Int64}} with indices -2:2×-2:2:
 0  0  0  0  0
 1  1  1  1  1
 1  1  1  1  1
 1  1  1  1  1
 0  0  0  0  0

The box shape SEBox is a special type for which many morphology algorithms may provide efficient implementations. For this reason, if one tries to collect an SEBoxArray into other array types (e.g. Array{Bool} via collect), then a significant performance drop is very likely to occur.

See also strel and strel_box.

strel_diamond(A::AbstractArray, [dims]; r=1)
strel_diamond(size, [dims]; [r])

Construct the N-dimensional structuring element (SE) for a diamond shape.

If image A is provided, then the SE size will be (2r+1, 2r+1, ...) with default half-size r=1. If size is provided, the default r will be maximum(size)÷2. The default dims will be all dimensions, that is, (1, 2, ..., length(size)).

julia> img = rand(64, 64);

julia> strel_diamond(img) # default size for image input is (3, 3)
3×3 SEDiamondArray{2, 2, UnitRange{Int64}, 0} with indices -1:1×-1:1:
 0  1  0
 1  1  1
 0  1  0

julia> strel_diamond(img; r=2) # equivalent to `strel_diamond((5,5))`
5×5 SEDiamondArray{2, 2, UnitRange{Int64}, 0} with indices -2:2×-2:2:
 0  0  1  0  0
 0  1  1  1  0
 1  1  1  1  1
 0  1  1  1  0
 0  0  1  0  0

julia> strel_diamond(img, (1,)) # mask along dimension 1
3×1 SEDiamondArray{2, 1, UnitRange{Int64}, 1} with indices -1:1×0:0:
 1
 1
 1

julia> strel_diamond((3,3), (1,)) # 3×3 mask along dimension 1
3×3 SEDiamondArray{2, 1, UnitRange{Int64}, 1} with indices -1:1×-1:1:
 0  1  0
 0  1  0
 0  1  0

The diamond shape SEDiamond is a special type for which many morphology algorithms may provide much more efficient implementations. For this reason, if one tries to collect an SEDiamondArray into other array types (e.g. Array{Bool} via collect), then a significant performance drop is very likely to occur.

See also strel and strel_box.

strel_type(x)

Infer the structuring element type for x.

This function is used to dispatch special SE types, e.g., SEBoxArray, to optimized implementation of particular morphology filter. In this sense it is required for custom SE array types to define this method.

 
strel_size(x)

Calculate the minimal block size that contains the structuring element. The result will be a tuple of odd integers.

julia> se = strel_diamond((5, 5); r=1)
5×5 SEDiamondArray{2, 2, UnitRange{Int64}, 0} with indices -2:2×-2:2:
 0  0  0  0  0
 0  0  1  0  0
 0  1  1  1  0
 0  0  1  0  0
 0  0  0  0  0

julia> strel_size(se) # is not (5, 5)
(3, 3)

julia> strel(Bool, strel(CartesianIndex, se)) # because it only checks the minimal enclosing block
3×3 OffsetArray(::BitMatrix, -1:1, -1:1) with eltype Bool with indices -1:1×-1:1:
 0  1  0
 1  1  1
 0  1  0

julia> se = [CartesianIndex(1, 1), CartesianIndex(-2, -2)];

julia> strel_size(se) # is not (4, 4)
(5, 5)

julia> strel(Bool, se) # because the connectivity mask has to be odd size
5×5 OffsetArray(::BitMatrix, -2:2, -2:2) with eltype Bool with indices -2:2×-2:2:
 1  0  0  0  0
 0  0  0  0  0
 0  0  1  0  0
 0  0  0  1  0
 0  0  0  0  0

julia> se = strel_diamond((5, 5), (1, ); r=1)
5×5 SEDiamondArray{2, 1, UnitRange{Int64}, 1} with indices -2:2×-2:2:
 0  0  0  0  0
 0  0  1  0  0
 0  0  1  0  0
 0  0  1  0  0
 0  0  0  0  0

julia> strel_size(se)
(3, 1)
strel_ndims(x)::Int

Infer the dimension of the structuring element x

upper, lower = strel_split([T], se)

Split a symmetric structuring element into its upper and lower half parts based on its center point.

For each element o in strel(CartesianIndex, upper), its negative -o is an element of strel(CartesianIndex, lower). This function is not the inverse of strel_chain.

The splited non-symmetric SE parts will be represented as array of T, where T is either a Bool or CartesianIndex. By default, T = eltype(se).

julia> se = strel_diamond((3, 3))
3×3 SEDiamondArray{2, 2, UnitRange{Int64}, 0} with indices -1:1×-1:1:
 0  1  0
 1  1  1
 0  1  0

julia> upper, lower = strel_split(se);

julia> upper
3×3 OffsetArray(::Matrix{Bool}, -1:1, -1:1) with eltype Bool with indices -1:1×-1:1:
 0  1  0
 1  1  0
 0  0  0

julia> lower
3×3 OffsetArray(::Matrix{Bool}, -1:1, -1:1) with eltype Bool with indices -1:1×-1:1:
 0  0  0
 0  1  1
 0  1  0

If the se is represented as displacement offset array, then the splited result will also be displacement offset array:

julia> se = strel(CartesianIndex, se)
4-element Vector{CartesianIndex{2}}:
 CartesianIndex(0, -1)
 CartesianIndex(-1, 0)
 CartesianIndex(1, 0)
 CartesianIndex(0, 1)

julia> upper, lower = strel_split(se);

julia> upper
2-element Vector{CartesianIndex{2}}:
 CartesianIndex(0, -1)
 CartesianIndex(-1, 0)

julia> lower
2-element Vector{CartesianIndex{2}}:
 CartesianIndex(1, 0)
 CartesianIndex(0, 1)
centered(A, cp=center(A)) -> Ao

Shift the center coordinate/point cp of array A to (0, 0, ..., 0). Internally, this is equivalent to OffsetArray(A, .-cp).

Совместимость: OffsetArrays 1.9

This method requires at least OffsetArrays 1.9.

Examples

julia> A = reshape(collect(1:9), 3, 3)
3×3 Matrix{Int64}:
 1  4  7
 2  5  8
 3  6  9

julia> Ao = OffsetArrays.centered(A); # axes (-1:1, -1:1)

julia> Ao[0, 0]
5

julia> Ao = OffsetArray(A, OffsetArrays.Origin(0)); # axes (0:2, 0:2)

julia> Aoo = OffsetArrays.centered(Ao); # axes (-1:1, -1:1)

julia> Aoo[0, 0]
5

Users are allowed to pass cp to change how "center point" is interpreted, but the meaning of the output array should be reinterpreted as well. For instance, if cp = map(last, axes(A)) then this function no longer shifts the center point but instead the bottom-right point to (0, 0, ..., 0). A commonly usage of cp is to change the rounding behavior when the array is of even size at some dimension:

julia> A = reshape(collect(1:4), 2, 2) # Ideally the center should be (1.5, 1.5) but OffsetArrays only support integer offsets
2×2 Matrix{Int64}:
 1  3
 2  4

julia> OffsetArrays.centered(A, OffsetArrays.center(A, RoundUp)) # set (2, 2) as the center point
2×2 OffsetArray(::Matrix{Int64}, -1:0, -1:0) with eltype Int64 with indices -1:0×-1:0:
 1  3
 2  4

julia> OffsetArrays.centered(A, OffsetArrays.center(A, RoundDown)) # set (1, 1) as the center point
2×2 OffsetArray(::Matrix{Int64}, 0:1, 0:1) with eltype Int64 with indices 0:1×0:1:
 1  3
 2  4

See also center.

center(A, [r::RoundingMode=RoundDown])::Dims

Return the center coordinate of given array A. If size(A, k) is even, a rounding procedure will be applied with mode r.

Совместимость: OffsetArrays 1.9

This method requires at least OffsetArrays 1.9.

Examples

julia> A = reshape(collect(1:9), 3, 3)
3×3 Matrix{Int64}:
 1  4  7
 2  5  8
 3  6  9

julia> c = OffsetArrays.center(A)
(2, 2)

julia> A[c...]
5

julia> Ao = OffsetArray(A, -2, -2); # axes (-1:1, -1:1)

julia> c = OffsetArrays.center(Ao)
(0, 0)

julia> Ao[c...]
5

To shift the center coordinate of the given array to (0, 0, ...), you can use centered.

is_symmetric(se)

Check if a given structuring element array se is symmetric with respect to its center pixel.

More formally, this checks if mask[I] == mask[-I] for any valid I ∈ CartesianIndices(mask) in the connectivity mask represetation mask = strel(Bool, se).

SEMask{N}()

A (holy) trait type for representing structuring element as connectivity mask. This connectivity mask SE is a bool array where true indicates that pixel position is connected to the center point.

julia> se = centered(Bool[0 1 0; 1 1 1; 0 1 0]) # commonly known as C4 connectivity
3×3 OffsetArray(::Matrix{Bool}, -1:1, -1:1) with eltype Bool with indices -1:1×-1:1:
 0  1  0
 1  1  1
 0  1  0

julia> strel_type(se)
ImageMorphology.StructuringElements.SEMask{2}()

See also SEOffset for the displacement offset representation. More details can be found on he documentation page Structuring Element.

SEOffset{N}()

A (holy) trait type for representing structuring element as displacement offsets. This displacement offsets SE is an array of CartesianIndex where each element stores the displacement offset from the center point.

julia> se = [CartesianIndex(-1, 0), CartesianIndex(0, -1), CartesianIndex(1, 0), CartesianIndex(0, 1)]
4-element Vector{CartesianIndex{2}}:
 CartesianIndex(-1, 0)
 CartesianIndex(0, -1)
 CartesianIndex(1, 0)
 CartesianIndex(0, 1)

julia> strel_type(se)
ImageMorphology.StructuringElements.SEOffset{2}()

See also SEMask for the connectivity mask representation. More details can be found on he documentation page Structuring Element.

SEDiamond{N}(axes, [dims]; [r])

A (holy) trait type for the N-dimensional diamond shape structuring element. This is a special case of SEMask that ImageMorphology algorithms might provide optimized implementation.

It is recommended to use strel_diamond and strel_type:

julia> se = strel_diamond((3, 3)) # C4 connectivity
3×3 SEDiamondArray{2, 2, UnitRange{Int64}, 0} with indices -1:1×-1:1:
 0  1  0
 1  1  1
 0  1  0

julia> strel_type(se)
SEDiamond{2, 2, UnitRange{Int64}}((-1:1, -1:1), (1, 2), 1)

julia> se = centered(collect(se)) # converted to normal centered array
3×3 OffsetArray(::Matrix{Bool}, -1:1, -1:1) with eltype Bool with indices -1:1×-1:1:
 0  1  0
 1  1  1
 0  1  0

julia> strel_type(se)
SEMask{2}()
SEBox{N}(axes; [r])

The N-dimensional structuring element with all elements connected. This is a special case of SEMask that ImageMorphology algorithms might provide optimized implementation.

It is recommended to use strel_box and strel_type:

julia> se = strel_box((3, 3)) # C8 connectivity
3×3 SEBoxArray{2, UnitRange{Int64}} with indices -1:1×-1:1:
 1  1  1
 1  1  1
 1  1  1

julia> strel_type(se)
SEBox{2, UnitRange{Int64}}((-1:1, -1:1), (1, 1))

julia> se = centered(collect(se)) # converted to normal centered array
3×3 OffsetArray(::Matrix{Bool}, -1:1, -1:1) with eltype Bool with indices -1:1×-1:1:
 1  1  1
 1  1  1
 1  1  1

julia> strel_type(se)
SEMask{2}()
SEDiamondArray(se::SEDiamond)

The instantiated array object of SEDiamond.

SEBoxArray(se::SEBox)

The instantiated array object of SEBox.

Morphological operations

extreme_filter(f, A; r=1, [dims]) -> out
extreme_filter(f, A, Ω) -> out

Filter the array A using select function f(x, y) for each Ω-neighborhood. The name "extreme" comes from the fact that typical select function f choice is min and max.

For each pixel p in A, the select function f is applied to its Ω-neighborhood iteratively in a f(...(f(f(A[p], A[p+Ω[1]]), A[p+Ω[2]]), ...) manner. For instance, in the 1-dimensional case, out[p] = f(f(A[p], A[p-1]), A[p+1]) for each p is the default behavior.

The Ω-neighborhood is defined by the dims or Ω argument. The r and dims keywords specifies the box shape neighborhood Ω using strel_box. The Ω is also known as structuring element (SE), it can be either displacement offsets or bool array mask, please refer to strel for more details.

Examples

julia> M = [4 6 5 3 4; 8 6 9 4 8; 7 8 4 9 6; 6 2 2 1 7; 1 6 5 2 6]
5×5 Matrix{Int64}:
 4  6  5  3  4
 8  6  9  4  8
 7  8  4  9  6
 6  2  2  1  7
 1  6  5  2  6

julia> extreme_filter(max, M) # max-filter using 4 direct neighbors along both dimensions
5×5 Matrix{Int64}:
 8  9  9  9  8
 8  9  9  9  9
 8  9  9  9  9
 8  8  9  9  9
 6  6  6  7  7

julia> extreme_filter(max, M; dims=1) # max-filter along the first dimension (column)
5×5 Matrix{Int64}:
 8  6  9  4  8
 8  8  9  9  8
 8  8  9  9  8
 7  8  5  9  7
 6  6  5  2  7

Ω can be either an AbstractArray{Bool} mask array with true element indicating connectivity, or a AbstractArray{<:CartesianIndex} array with each element indicating the displacement offset to its center element.

julia> Ω_mask = centered(Bool[1 1 0; 1 1 0; 1 0 0]) # custom neighborhood in mask format
3×3 OffsetArray(::Matrix{Bool}, -1:1, -1:1) with eltype Bool with indices -1:1×-1:1:
 1  1  0
 1  1  0
 1  0  0

julia> out = extreme_filter(max, M, Ω_mask)
5×5 Matrix{Int64}:
 4  8  6  9  4
 8  8  9  9  9
 8  8  9  9  9
 7  8  8  9  9
 6  6  6  5  7

julia> Ω_offsets = strel(CartesianIndex, Ω_mask) # custom neighborhood as displacement offset
4-element Vector{CartesianIndex{2}}:
 CartesianIndex(-1, -1)
 CartesianIndex(0, -1)
 CartesianIndex(1, -1)
 CartesianIndex(-1, 0)

julia> out == extreme_filter(max, M, Ω_offsets) # both versions work equivalently
true

See also the in-place version extreme_filter!. Another function in ImageFiltering package ImageFiltering.mapwindow provides similar functionality.

extreme_filter!(f, out, A; [r], [dims])
extreme_filter!(f, out, A, Ω)

The in-place version of extreme_filter where out is the output array that gets modified.

dilate(img; dims=coords_spatial(img), r=1)
dilate(img, se)

Perform a max-filter over the neighborhood of img, specified by structuring element se.

se is the structuring element that defines the neighborhood of the image. See strel for more details. If se is not specified, then it will use the strel_box with an extra keyword dims to control the dimensions to filter, and half-size r to control the diamond size.

Examples

julia> img = falses(5, 5); img[3, [2, 4]] .= true; img
5×5 BitMatrix:
 0  0  0  0  0
 0  0  0  0  0
 0  1  0  1  0
 0  0  0  0  0
 0  0  0  0  0

julia> dilate(img)
5×5 BitMatrix:
 0  0  0  0  0
 1  1  1  1  1
 1  1  1  1  1
 1  1  1  1  1
 0  0  0  0  0

julia> dilate(img; dims=1)
5×5 BitMatrix:
 0  0  0  0  0
 0  1  0  1  0
 0  1  0  1  0
 0  1  0  1  0
 0  0  0  0  0

julia> dilate(img, strel_diamond(img)) # use diamond shape SE
5×5 BitMatrix:
 0  0  0  0  0
 0  1  0  1  0
 1  1  1  1  1
 0  1  0  1  0
 0  0  0  0  0

See also

  • dilate! is the in-place version of this function

  • erode is the dual operator of dilate in the sense that complement.(dilate(img)) == erode(complement.(img)).

If se is symmetric with repsect to origin, i.e., se[b] == se[-b] for any b, then dilation becomes the Minkowski sum: A⊕B={a+b|a∈A, b∈B}.

 
dilate!(out, img; [dims], [r])
dilate!(out, img, se)

The in-place version of dilate with input image img and output image out.

out = erode(img; dims=coords_spatial(img), r=1)
out = erode(img, se)

Perform a min-filter over the neighborhood of img, specified by structuring element se.

se is the structuring element that defines the neighborhood of the image. See strel for more details. If se is not specified, then it will use the strel_box with an extra keyword dims to control the dimensions to filter, and half-size r to control the diamond size.

Examples

julia> img = trues(5, 5); img[3, [2, 4]] .= false; img
5×5 BitMatrix:
 1  1  1  1  1
 1  1  1  1  1
 1  0  1  0  1
 1  1  1  1  1
 1  1  1  1  1

julia> erode(img)
5×5 BitMatrix:
 1  1  1  1  1
 0  0  0  0  0
 0  0  0  0  0
 0  0  0  0  0
 1  1  1  1  1

julia> erode(img; dims=1)
5×5 BitMatrix:
 1  1  1  1  1
 1  0  1  0  1
 1  0  1  0  1
 1  0  1  0  1
 1  1  1  1  1

julia> erode(img, strel_diamond(img)) # use diamond shape SE
5×5 BitMatrix:
 1  1  1  1  1
 1  0  1  0  1
 0  0  0  0  0
 1  0  1  0  1
 1  1  1  1  1

See also

  • erode! is the in-place version of this function

  • dilate is the dual operator of erode in the sense that complement.(dilate(img)) == erode(complement.(img)).

If se is symmetric with repsect to origin, i.e., se[b] == se[-b] for any b, then erosion becomes the Minkowski difference: A⊖B={a-b|a∈A, b∈B}.

 
erode!(out, img; [dims], [r])
erode!(out, img, se)

The in-place version of erode with input image img and output image out.

opening(img; dims=coords_spatial(img), r=1)
opening(img, se)

Perform the morphological opening on img. The opening operation is defined as erosion followed by a dilation: dilate(erode(img, se), se).

se is the structuring element that defines the neighborhood of the image. See strel for more details. If se is not specified, then it will use the strel_box with an extra keyword dims to control the dimensions to filter, and half-size r to control the diamond size.

Examples

julia> img = trues(7,7); img[2, 2] = false; img[3:5, 3:5] .= false; img[4, 4] = true; img
7×7 BitMatrix:
 1  1  1  1  1  1  1
 1  0  1  1  1  1  1
 1  1  0  0  0  1  1
 1  1  0  1  0  1  1
 1  1  0  0  0  1  1
 1  1  1  1  1  1  1
 1  1  1  1  1  1  1

julia> opening(img)
7×7 BitMatrix:
 0  0  1  1  1  1  1
 0  0  1  1  1  1  1
 1  1  0  0  0  1  1
 1  1  0  0  0  1  1
 1  1  0  0  0  1  1
 1  1  1  1  1  1  1
 1  1  1  1  1  1  1

julia> opening(img, strel_diamond(img)) # use diamond shape SE
7×7 BitMatrix:
 1  1  1  1  1  1  1
 1  0  1  1  1  1  1
 1  1  0  0  0  1  1
 1  1  0  0  0  1  1
 1  1  0  0  0  1  1
 1  1  1  1  1  1  1
 1  1  1  1  1  1  1

See also

  • opening! is the in-place version of this function.

  • closing is the dual operator of opening in the sense that complement.(opening(img)) == closing(complement.(img)).

opening!(out, img, buffer; [dims], [r])
opening!(out, img, se, buffer)

The in-place version of opening with input image img and output image out. The intermediate erosion result is stored in buffer.

closing(img; dims=coords_spatial(img), r=1)
closing(img, se)

Perform the morphological closing on img. The closing operation is defined as dilation followed by an erosion: erode(dilate(img, se), se).

se is the structuring element that defines the neighborhood of the image. See strel for more details. If se is not specified, then it will use the strel_box with an extra keyword dims to control the dimensions to filter, and half-size r to control the diamond size.

Examples

julia> img = falses(7,7); img[2, 2] = true; img[3:5, 3:5] .= true; img[4, 4] = false; img
7×7 BitMatrix:
 0  0  0  0  0  0  0
 0  1  0  0  0  0  0
 0  0  1  1  1  0  0
 0  0  1  0  1  0  0
 0  0  1  1  1  0  0
 0  0  0  0  0  0  0
 0  0  0  0  0  0  0

julia> closing(img)
7×7 BitMatrix:
 1  1  0  0  0  0  0
 1  1  0  0  0  0  0
 0  0  1  1  1  0  0
 0  0  1  1  1  0  0
 0  0  1  1  1  0  0
 0  0  0  0  0  0  0
 0  0  0  0  0  0  0

julia> closing(img, strel_diamond(img)) # # use diamond shape SE
7×7 BitMatrix:
 0  0  0  0  0  0  0
 0  1  0  0  0  0  0
 0  0  1  1  1  0  0
 0  0  1  1  1  0  0
 0  0  1  1  1  0  0
 0  0  0  0  0  0  0
 0  0  0  0  0  0  0

See also

  • opening! is the in-place version of this function.

  • closing is the dual operator of opening in the sense that complement.(opening(img)) == closing(complement.(img)).

closing!(out, img, buffer; [dims], [r])
closing!(out, img, se, buffer)

The in-place version of closing with input image img and output image out. The intermediate dilation result is stored in buffer.

tophat(img; dims=coords_spatial(img), r=1)
tophat(img, se)

Performs morphological top-hat transform for given image, i.e., img - opening(img, se).

se is the structuring element that defines the neighborhood of the image. See strel for more details. If se is not specified, then it will use the strel_box with an extra keyword dims to control the dimensions to filter, and half-size r to control the diamond size.

This white top-hat transform can be used to extract small white elements and details from an image. To extract black details, the black top-hat transform, also known as bottom-hat transform, bothat can be used.

Examples

julia> img = falses(5, 5); img[1, 1] = true; img[3:5, 3:5] .= true; img
5×5 BitMatrix:
 1  0  0  0  0
 0  0  0  0  0
 0  0  1  1  1
 0  0  1  1  1
 0  0  1  1  1

julia> Int.(tophat(img))
5×5 Matrix{Int64}:
 1  0  0  0  0
 0  0  0  0  0
 0  0  0  0  0
 0  0  0  0  0
 0  0  0  0  0

julia> Int.(tophat(img, strel_diamond(img))) # use diamond shape SE
5×5 Matrix{Int64}:
 1  0  0  0  0
 0  0  0  0  0
 0  0  1  0  0
 0  0  0  0  0
 0  0  0  0  0
tophat!(out, img, buffer; [dims], [r])
tophat!(out, img, se, buffer)

The in-place version of tophat with input image img and output image out. The intermediate erosion result is stored in buffer.

bothat(img; dims=coords_spatial(img), r=1)
bothat(img, se)

Performs morphological bottom-hat transform for given image, i.e., closing(img, se) - img.

se is the structuring element that defines the neighborhood of the image. See strel for more details. If se is not specified, then it will use the strel_box with an extra keyword dims to control the dimensions to filter, and half-size r to control the diamond size.

This bottom-hat transform, also known as black top-hat transform, can be used to extract small black elements and details from an image. To extract white details, the white top-hat transform tophat can be used.

Examples

julia> img = falses(7, 7); img[3:5, 3:5] .= true; img[4, 6] = true; img[4, 4] = false; img
7×7 BitMatrix:
 0  0  0  0  0  0  0
 0  0  0  0  0  0  0
 0  0  1  1  1  0  0
 0  0  1  0  1  1  0
 0  0  1  1  1  0  0
 0  0  0  0  0  0  0
 0  0  0  0  0  0  0

julia> Int.(bothat(img))
7×7 Matrix{Int64}:
 0  0  0  0  0  0  0
 0  0  0  0  0  0  0
 0  0  0  0  0  0  0
 0  0  0  1  0  0  1
 0  0  0  0  0  0  0
 0  0  0  0  0  0  0
 0  0  0  0  0  0  0

julia> Int.(bothat(img, strel_diamond(img))) # use diamond shape SE
7×7 Matrix{Int64}:
 0  0  0  0  0  0  0
 0  0  0  0  0  0  0
 0  0  0  0  0  0  0
 0  0  0  1  0  0  0
 0  0  0  0  0  0  0
 0  0  0  0  0  0  0
 0  0  0  0  0  0  0

See also bothat! for the in-place version.

bothat!(out, img, buffer; [dims], [r])
bothat!(out, img, se, buffer)

The in-place version of bothat with input image img and output image out. The intermediate dilation result is stored in buffer.

mgradient(img; mode=:beucher, dims=coords_spatial(img), r=1)
mgradient(img, se; mode=:beucher)

Calculate morphological gradient of the image using given mode.

There are three widely used modes[1]:

  • :beucher: the default mode. It calculates the arithmetic difference between the dilation and the erosion — dilate(img, se) - erode(img, se).

  • :internal: also known as half-gradient by erosion. It calculates the arithmetic difference between the original image and its erosion — img - erode(img, se).

  • :external: also known as half-gradient by dilation. It calculates the arithmetic difference between dilation and the original image — dilate(img, se) - se.

se is the structuring element that defines the neighborhood of the image. See strel for more details. If se is not specified, then it will use the strel_box with an extra keyword dims to control the dimensions to filter, and half-size r to control the diamond size.

Examples

julia> img = falses(7, 7); img[3:5, 3:5] .= true; img
7×7 BitMatrix:
 0  0  0  0  0  0  0
 0  0  0  0  0  0  0
 0  0  1  1  1  0  0
 0  0  1  1  1  0  0
 0  0  1  1  1  0  0
 0  0  0  0  0  0  0
 0  0  0  0  0  0  0

julia> Int.(mgradient(img)) # default mode :beucher always creates a two-pixel wide boundary
7×7 Matrix{Int64}:
 0  0  0  0  0  0  0
 0  1  1  1  1  1  0
 0  1  1  1  1  1  0
 0  1  1  0  1  1  0
 0  1  1  1  1  1  0
 0  1  1  1  1  1  0
 0  0  0  0  0  0  0

julia> Int.(mgradient(img; mode=:internal)) # half-gradient -- the boundary is internal to original image
7×7 Matrix{Int64}:
 0  0  0  0  0  0  0
 0  0  0  0  0  0  0
 0  0  1  1  1  0  0
 0  0  1  0  1  0  0
 0  0  1  1  1  0  0
 0  0  0  0  0  0  0
 0  0  0  0  0  0  0

julia> Int.(mgradient(img; mode=:external)) # half-gradient -- the boundary is external to original image
7×7 Matrix{Int64}:
 0  0  0  0  0  0  0
 0  1  1  1  1  1  0
 0  1  0  0  0  1  0
 0  1  0  0  0  1  0
 0  1  0  0  0  1  0
 0  1  1  1  1  1  0
 0  0  0  0  0  0  0

julia> Int.(mgradient(img, strel_diamond(img))) # use diamond shape SE
7×7 Matrix{Int64}:
 0  0  0  0  0  0  0
 0  0  1  1  1  0  0
 0  1  1  1  1  1  0
 0  1  1  0  1  1  0
 0  1  1  1  1  1  0
 0  0  1  1  1  0  0
 0  0  0  0  0  0  0

The beucher operator is a self-complementary operator in the sense that mgradient(img, se; mode=:beucher) == mgradient(complement.(img), se; mode=:beucher). When r>1, it is usually called thick gradient. If a line segment is used as se, then the gradient becomes the directional gradient.

See also

  • mgradient! is the in-place version of this function.

  • mlaplacian for the laplacian operator.

  • ImageBase.FiniteDiff also provides a few finite difference operators, including fdiff, fgradient, etc.

References

  • [1] Rivest, Jean-Francois, Pierre Soille, and Serge Beucher. "Morphological gradients." Journal of Electronic Imaging 2.4 (1993): 326-336.

mgradient!(out, img, buffer; [dims], [r], [mode])
mgradient!(out, img, se, buffer; [mode])

The in-place version of mgradient with input image img and output image out.

The buffer array is required for :beucher mode. For :internal and :external modes, buffer is not needed and can be nothing.

mlaplacian(img; dims=coords_spatial(img), r=1)
mlaplacian(img, se)

Calculate morphological laplacian of the image.

The morphological lapalacian operator is defined as ∇⁺A - ∇⁻A where ∇⁺A is the external gradient A - erode(A, se) and ∇⁻A is the internal gradient dilate(A, se) - A. Thus is dilate(A, se) + erode(A, se) - 2A.

se is the structuring element that defines the neighborhood of the image. See strel for more details. If se is not specified, then it will use the strel_box with an extra keyword dims to control the dimensions to filter.

Examples

julia> img = falses(7, 7); img[3:5, 3:5] .= true; img[4, 4] = false; img
7×7 BitMatrix:
 0  0  0  0  0  0  0
 0  0  0  0  0  0  0
 0  0  1  1  1  0  0
 0  0  1  0  1  0  0
 0  0  1  1  1  0  0
 0  0  0  0  0  0  0
 0  0  0  0  0  0  0

julia> Int.(mlaplacian(img))
7×7 Matrix{Int64}:
 0  0   0   0   0  0  0
 0  1   1   1   1  1  0
 0  1  -1  -1  -1  1  0
 0  1  -1   1  -1  1  0
 0  1  -1  -1  -1  1  0
 0  1   1   1   1  1  0
 0  0   0   0   0  0  0

julia> Int.(mlaplacian(img, strel_diamond(img))) # use diamond shape SE
7×7 Matrix{Int64}:
 0  0   0   0   0  0  0
 0  0   1   1   1  0  0
 0  1  -1  -1  -1  1  0
 0  1  -1   1  -1  1  0
 0  1  -1  -1  -1  1  0
 0  0   1   1   1  0  0
 0  0   0   0   0  0  0

See also

  • mlaplacian! is the in-place version of this function.

  • mgradient for the gradient operator.

  • ImageBase.FiniteDiff also provides a few finite difference operators, including fdiff, fgradient, etc.

mlaplacian!(out, img, buffer; [dims], [r])
mlaplacian!(out, img, se, buffer)

The in-place version of mlaplacian with input image img and output image out. The intermediate erosion result is stored in buffer.

Geodesic operations

mreconstruct(op, marker, mask; [dims])
mreconstruct(op, marker, mask, se)

Morphological reconstruction of marker image by operation op.

The op argument is either erode or dilate, indicating reconstruction by erosion or by dilation. The mask argument has the same shape as marker and is used to restrict the output value range.

The dims keyword is used to specify the dimension to process by constructing the box shape structuring element strel_box(marker; dims). For generic structuring element, the half-size is expected to be either 0 or 1 along each dimension.

By definition, the reconstruction is done by applying marker = select.(op(marker; dims), mask) repeatly until reaching stability. For dilation op, select = dilate, min and for erosion op, select = erode, max.

Examples

julia> marker = [0 0 0 0 0; 0 9 0 0 0; 0 0 0 0 0; 0 0 0 5 0; 0 0 0 0 0; 0 9 0 0 0]
6×5 Matrix{Int64}:
 0  0  0  0  0
 0  9  0  0  0
 0  0  0  0  0
 0  0  0  5  0
 0  0  0  0  0
 0  9  0  0  0

julia> mask = [9 0 0 0 0; 0 8 7 1 0; 0 9 0 4 0; 0 0 0 4 0; 0 0 6 5 6; 0 0 9 8 9]
6×5 Matrix{Int64}:
 9  0  0  0  0
 0  8  7  1  0
 0  9  0  4  0
 0  0  0  4  0
 0  0  6  5  6
 0  0  9  8  9

julia> mreconstruct(dilate, marker, mask) # equivalent to underbuild(marker, mask)
6×5 Matrix{Int64}:
 8  0  0  0  0
 0  8  7  1  0
 0  8  0  4  0
 0  0  0  4  0
 0  0  4  4  4
 0  0  4  4  4

See also

The inplace version of this function is mreconstruct!. There are also aliases underbuild for reconstruction by dilation and overbuild for reconstruction by erosion.

References

  • [1] L. Vincent, “Morphological grayscale reconstruction in image analysis: applications and efficient algorithms,” IEEE Trans. on Image Process., vol. 2, no. 2, pp. 176—​201, Apr. 1993, doi: 10.1109/83.217222.

  • [2] P. Soille, Morphological Image Analysis. Berlin, Heidelberg: Springer Berlin Heidelberg, 2004. doi: 10.1007/978-3-662-05088-0.

mreconstruct!(op, out, marker, mask; [dims])

The in-place version of morphological reconstruction mreconstruct.

underbuild(marker, mask; [dims])
underbuild(marker, mask, se)

Reconstruction by dilation. This is an alias for mreconstruct with op=dilate.

See also the in-place version underbuild!, and the dual operator overbuild.

underbuild!(out, marker, mask; [dims])
underbuild!(out, marker, mask, se)

The in-place version of underbuild with output image out being modified in place.

overbuild(marker, mask; [dims])
overbuild(marker, mask, se)

Reconstruction by erosion. This is an alias for mreconstruct with op=erode.

See also the in-place version overbuild!, and the dual operator underbuild.

overbuild!(out, marker, mask; [dims])
overbuild!(out, marker, mask, se)

The in-place version of overbuild with output image out being modified in place.

Components and segmentation

label = label_components(A; [dims=coords_spatial(A)], [r=1], [bkg])
label = label_components(A, se; [bkg])

Find and label the connected components of array A where the connectivity is defined by structuring element se. Each component is assigned a unique integer value as its label with 0 representing the background defined by bkg.

se is the structuring element that defines the neighborhood of the image. See strel for more details. If se is not specified, then it will use the strel_box with an extra keyword dims to control the dimensions to filter, and half-size r to control the diamond size.

Examples

julia> A = [false true false true  false;
            true false false  true  true]
2×5 Matrix{Bool}:
 0  1  0  1  0
 1  0  0  1  1

julia> label_components(A) # default diamond shape C4 connectivity
2×5 Matrix{Int64}:
 0  2  0  3  0
 1  0  0  3  3

julia> label_components(A; dims=2) # only the rows are considered
2×5 Matrix{Int64}:
 0  2  0  3  0
 1  0  0  4  4

julia> label_components(A, strel_box((3, 3))) # box shape C8 connectivity
2×5 Matrix{Int64}:
 0  1  0  2  0
 1  0  0  2  2

The in-place version is label_components!. See also component_boxes, component_lengths, component_indices, component_centroids for basic properties of the labeled components.

label_components!(out, A; [dims], [r] [bkg])
label_components!(out, A, se; [bkg])

The in-place version of label_components.

boxes = component_boxes(labeled_array)

Calculates the minimal bounding boxes for each label including the background label. The labels can be computed by label_components.

Each bounding box is represented as a CartesianIndices. boxes is shifted to 0-based indexing vector so that background region is boxes[0].

julia> A = [2 2 2 2 2; 1 1 1 0 1; 1 0 2 1 1; 1 1 2 2 2; 1 0 2 2 2]
5×5 Matrix{Int64}:
 2  2  2  2  2
 1  1  1  0  1
 1  0  2  1  1
 1  1  2  2  2
 1  0  2  2  2

julia> label = label_components(A) # four disjoint components
5×5 Matrix{Int64}:
 1  1  1  1  1
 2  2  2  0  4
 2  0  3  4  4
 2  2  3  3  3
 2  0  3  3  3

julia> boxes = component_boxes(label) # get bounding boxes of all regions
5-element OffsetArray(::Vector{CartesianIndices{2, Tuple{UnitRange{Int64}, UnitRange{Int64}}}}, 0:4) with eltype CartesianIndices{2, Tuple{UnitRange{Int64}, UnitRange{Int64}}} with indices 0:4:
 CartesianIndices((2:5, 2:4))
 CartesianIndices((1:1, 1:5))
 CartesianIndices((2:5, 1:3))
 CartesianIndices((3:5, 3:5))
 CartesianIndices((2:3, 4:5))

julia> A[boxes[1]] # crop the image region with label 1
1×5 Matrix{Int64}:
 2  2  2  2  2

julia> A[boxes[4]] # crop the image region with label 4
2×2 Matrix{Int64}:
 0  1
 1  1
counts = component_lengths(labeled_array)

Count the number of each labels in the input labeled array. counts is shifted to 0-based indexing vector so that the number of background pixels is counts[0].

julia> A = [2 2 2 2 2; 1 1 1 0 1; 1 0 2 1 1; 1 1 2 2 2; 1 0 2 2 2]
5×5 Matrix{Int64}:
 2  2  2  2  2
 1  1  1  0  1
 1  0  2  1  1
 1  1  2  2  2
 1  0  2  2  2

julia> label = label_components(A) # four disjoint components
5×5 Matrix{Int64}:
 1  1  1  1  1
 2  2  2  0  4
 2  0  3  4  4
 2  2  3  3  3
 2  0  3  3  3

julia> component_lengths(label)
5-element OffsetArray(::Vector{Int64}, 0:4) with eltype Int64 with indices 0:4:
 3
 5
 7
 7
 3

For gray images, labels can be computed by label_components.

indices = component_indices([T], labeled_array)

Get the indices of each label in the input labeled array. indices is shifted to 0-based indexing vector so that the indices of background pixels is indices[0].

The optional type T can be either Int/IndexLinear() or CartesianIndex/IndexCartesian() that is used to specify the type of the indices. The default choice is IndexStyle(labeled_array).

julia> A = [2 2 2 2 2; 1 1 1 0 1; 1 0 2 1 1; 1 1 2 2 2; 1 0 2 2 2]
5×5 Matrix{Int64}:
 2  2  2  2  2
 1  1  1  0  1
 1  0  2  1  1
 1  1  2  2  2
 1  0  2  2  2

julia> label = label_components(A) # four disjoint components
5×5 Matrix{Int64}:
 1  1  1  1  1
 2  2  2  0  4
 2  0  3  4  4
 2  2  3  3  3
 2  0  3  3  3

julia> indices = component_indices(label)
5-element OffsetArray(::Vector{Vector{Int64}}, 0:4) with eltype Vector{Int64} with indices 0:4:
 [8, 10, 17]
 [1, 6, 11, 16, 21]
 [2, 3, 4, 5, 7, 9, 12]
 [13, 14, 15, 19, 20, 24, 25]
 [18, 22, 23]

julia> indices = component_indices(CartesianIndex, label)
5-element OffsetArray(::Vector{Vector{CartesianIndex{2}}}, 0:4) with eltype Vector{CartesianIndex{2}} with indices 0:4:
 [CartesianIndex(3, 2), CartesianIndex(5, 2), CartesianIndex(2, 4)]
 [CartesianIndex(1, 1), CartesianIndex(1, 2), CartesianIndex(1, 3), CartesianIndex(1, 4), CartesianIndex(1, 5)]
 [CartesianIndex(2, 1), CartesianIndex(3, 1), CartesianIndex(4, 1), CartesianIndex(5, 1), CartesianIndex(2, 2), CartesianIndex(4, 2), CartesianIndex(2, 3)]
 [CartesianIndex(3, 3), CartesianIndex(4, 3), CartesianIndex(5, 3), CartesianIndex(4, 4), CartesianIndex(5, 4), CartesianIndex(4, 5), CartesianIndex(5, 5)]
 [CartesianIndex(3, 4), CartesianIndex(2, 5), CartesianIndex(3, 5)]

For gray images, labels can be computed by label_components.

centroids = component_centroids(labeled_array)

Compute the centroid of each label in the input labeled array. centroids is shifted to 0-based indexing vector so that the centroid of background pixels is centroids[0].

The centroid of a finite set X, also known as geometric center, is calculated using sum(X)/length(X). For label i, all (Cartesian) indices of pixels with label i are used to build the set X

julia> A = [2 2 2 2 2; 1 1 1 0 1; 1 0 2 1 1; 1 1 2 2 2; 1 0 2 2 2]
5×5 Matrix{Int64}:
 2  2  2  2  2
 1  1  1  0  1
 1  0  2  1  1
 1  1  2  2  2
 1  0  2  2  2

julia> label = label_components(A) # four disjoint components
5×5 Matrix{Int64}:
 1  1  1  1  1
 2  2  2  0  4
 2  0  3  4  4
 2  2  3  3  3
 2  0  3  3  3

julia> component_centroids(label)
5-element OffsetArray(::Vector{Tuple{Float64, Float64}}, 0:4) with eltype Tuple{Float64, Float64} with indices 0:4:
 (3.3333333333333335, 2.6666666666666665)
 (1.0, 3.0)
 (3.142857142857143, 1.5714285714285714)
 (4.285714285714286, 3.857142857142857)
 (2.6666666666666665, 4.666666666666667)

For gray images, labels can be computed by label_components.

Max tree

Max-tree morphological representation of an image.

Details

Let’s consider a thresholding operation,

    mask = [val ≥ threshold for val in image]

One can identify the connected components (the sets of neighboring true values) in mask. When image thresholding is sequentially applied for all possible thresholds, it generates a collection of connected components that could be organized into a hierarchical structure called component tree. Consider 1D "image" with values 1, 2 and 3:

       2233233312223322

The connected components would be

    1: AAAAAAAAAAAAAAAA
    2: BBBBBBBB.CCCCCCC
    3: ..DD.EEE....FF..

Here, the letters are the labels of the resulting connected components, and . specifies that the pixel value is below the threshold. In this example, the corresponding component tree is:

      A
     ⭩ ⭨
    B   C
   ⭩ ⭨   ⭨
  D   E   F

A max-tree is an efficient representation of the component tree. A connected component at threshold level is represented by the single reference pixel from this level (image[r] == t), which is the parent to all other pixels of and also to the reference pixels of the connected components at higher thresholds, which are the children of . In our example, the reference pixels (denoted by the letter of the corresponding component) would be:

    1: ........A.......
    2: B........C......
    3: ..D..E......F...

I.e.

Comp Ref.Pixel

A

9

B

1

C

10

D

3

E

6

F

13

So the whole max-tree could be encoded as a vector of indices of parent pixels:

9  1  1  3  1  1  6  6  9  9 10 10 10 13 10 10

The max-tree is the basis for many morphological operators, namely connected operators. Unlike morphological openings and closings, these operators do not require a fixed structuring element, but rather act with a flexible structuring element that meets a certain criterion.

See also

area_opening, area_closing, diameter_opening, diameter_closing.

References

  1. Salembier, P., Oliveras, A., & Garrido, L. (1998). Antiextensive Connected Operators for Image and Sequence Processing. IEEE Transactions on Image Processing, 7(4), 555-570.

    https://doi.org/10.1109/83.663500

  2. Berger, C., Geraud, T., Levillain, R., Widynski, N., Baillard, A., Bertin, E. (2007). Effective Component Tree Computation with Application to Pattern Recognition in Astronomical Imaging. In International Conference on Image Processing (ICIP), 41-44.

    https://doi.org/10.1109/ICIP.2007.4379949

  3. Najman, L., & Couprie, M. (2006). Building the component tree in quasi-linear time. IEEE Transactions on Image Processing, 15(11), 3531-3539.

    https://doi.org/10.1109/TIP.2006.877518

  4. Carlinet, E., & Geraud, T. (2014). A Comparative Review of Component Tree Computation Algorithms. IEEE Transactions on Image Processing, 23(9), 3885-3895.

    https://doi.org/10.1109/TIP.2014.2336551

areas(maxtree::MaxTree) -> Array{Int}

Computes the areas of all maxtree components.

Returns

The array of the same shape as the original image. The i-th element is the area (in pixels) of the component that is represented by the reference pixel with index i.

See also

diameters, area_opening, area_closing.

boundingboxes(maxtree::MaxTree) -> Array{NTuple{2, CartesianIndex}}

Computes the minimal bounding boxes of all maxtree components.

Returns

The array of the same shape as the original image. The i-th element is the tuple of the minimal and maximal cartesian indices for the bounding box of the component that is represented by the reference pixel with index i.

See also

diameters.

diameters(maxtree::MaxTree) -> Array{Int}

Computes the "diameters" of all maxtree components.

"Diameter" of the max-tree connected component is the length of the widest side of the component’s bounding box.

Returns

The array of the same shape as the original image. The i-th element is the "diameter" of the component that is represented by the reference pixel with index i.

See also

boundingboxes, areas, diameter_opening, diameter_closing.

area_opening(image, [maxtree]; min_area=64, connectivity=1) -> Array

Performs an area opening of the image.

Area opening replaces all bright components of an image that have a surface smaller than min_area with the darker value taken from their first ancestral component (in max-tree representation of image) that has the area no smaller than min_area.

Details

Area opening is similar to morphological opening (see opening), but instead of using a fixed structuring element (e.g. disk) it employs small (less than min_area) components of the max-tree. Consequently, the area_opening with min_area = 1 is the identity transformation.

In the binary case, area opening is equivalent to remove_small_objects; this operator is thus extended to gray-level images.

Arguments

  • image::GenericGrayImage: the -dimensional input image

  • min_area::Number=64: the smallest size (in pixels) of the image component to keep intact

  • connectivity::Integer=1: the neighborhood connectivity. The maximum number of orthogonal steps to reach a neighbor of the pixel. In 2D, it is 1 for a 4-neighborhood and 2 for a 8-neighborhood.

  • maxtree::MaxTree: optional pre-built max-tree. Note that maxtree and connectivity optional parameters are mutually exclusive.

Returns

An array of the same type and shape as the image.

See also

area_opening!, area_closing, diameter_opening, MaxTree, opening

References

  1. Vincent, L. (1993). Grayscale area openings and closings, their efficient implementation and applications, Proc. of EURASIP Workshop on Mathematical Morphology and its Applications to Signal Processing, Barcelona, Spain, 22-27

  2. Soille, P. (2003). Chapter 6 Geodesic Metrics of Morphological Image Analysis: Principles and Applications, 2nd edition, Springer.

    https://doi.org/10.1007/978-3-662-05088-0

  3. Salembier, P., Oliveras, A., & Garrido, L. (1998). Antiextensive Connected Operators for Image and Sequence Processing. IEEE Transactions on Image Processing, 7(4), 555-570.

    https://doi.org/10.1109/83.663500

  4. Najman, L., & Couprie, M. (2006). Building the component tree in quasi-linear time. IEEE Transactions on Image Processing, 15(11), 3531-3539.

    https://doi.org/10.1109/TIP.2006.877518

  5. Carlinet, E., & Geraud, T. (2014). A Comparative Review of Component Tree Computation Algorithms. IEEE Transactions on Image Processing, 23(9), 3885-3895.

    https://doi.org/10.1109/TIP.2014.2336551

Examples

Creating a test image f (quadratic function with a maximum in the center and 4 additional local maxima):

julia> w = 12;

julia> f = [20 - 0.2*((x - w/2)^2 + (y-w/2)^2) for x in 0:w, y in 0:w];

julia> f[3:4, 2:6] .= 40; f[3:5, 10:12] .= 60; f[10:12, 3:5] .= 80;

julia> f[10:11, 10:12] .= 100; f[11, 11] = 100;

julia> f_aopen = area_opening(f, min_area=8, connectivity=1);

The peaks with a surface smaller than 8 are removed.

area_opening!(output, image, [maxtree];
              min_area=64, connectivity=1) -> output

Performs in-place area opening of the image and stores the result in output. See area_opening for the detailed description of the method.

area_closing(image, [maxtree]; min_area=64, connectivity=1) -> Array

Performs an area closing of the image.

Area closing replaces all dark components of an image that have a surface smaller than min_area with the brighter value taken from their first ancestral component (in max-tree representation of image) that has the area no smaller than min_area.

Details

Area closing is the dual operation to area opening (see area_opening). It is similar to morphological closings (see closing), but instead of using a fixed structuring element (e.g. disk) it employs small (less than min_area) components of the max-tree. Consequently, the area_closing with min_area = 1 is the identity transformation.

In the binary case, area closing is equivalent to remove_small_holes; this operator is thus extended to gray-level images.

Arguments

  • image::GenericGrayImage: the -dimensional input image

  • min_area::Number=64: the smallest size (in pixels) of the image component to keep intact

  • connectivity::Integer=1: the neighborhood connectivity. The maximum number of orthogonal steps to reach a neighbor of the pixel. In 2D, it is 1 for a 4-neighborhood and 2 for a 8-neighborhood.

  • maxtree::MaxTree: optional pre-built max-tree. Note that maxtree and connectivity optional parameters are mutually exclusive.

Returns

An array of the same type and shape as the image.

See also

area_closing!, area_opening, diameter_closing, MaxTree, closing

References

  1. Vincent, L. (1993). Grayscale area openings and closings, their efficient implementation and applications, Proc. of EURASIP Workshop on Mathematical Morphology and its Applications to Signal Processing, Barcelona, Spain, 22-27

  2. Soille, P. (2003). Chapter 6 Geodesic Metrics of Morphological Image Analysis: Principles and Applications, 2nd edition, Springer.

    https://doi.org/10.1007/978-3-662-05088-0

  3. Salembier, P., Oliveras, A., & Garrido, L. (1998). Antiextensive Connected Operators for Image and Sequence Processing. IEEE Transactions on Image Processing, 7(4), 555-570.

    https://doi.org/10.1109/83.663500

  4. Najman, L., & Couprie, M. (2006). Building the component tree in quasi-linear time. IEEE Transactions on Image Processing, 15(11), 3531-3539.

    https://doi.org/10.1109/TIP.2006.877518

  5. Carlinet, E., & Geraud, T. (2014). A Comparative Review of Component Tree Computation Algorithms. IEEE Transactions on Image Processing, 23(9), 3885-3895.

    https://doi.org/10.1109/TIP.2014.2336551

Examples

Creating a test image f (quadratic function with a minimum in the center and 4 additional local minima):

julia> w = 12;

julia> f = [180 + 0.2*((x - w/2)^2 + (y-w/2)^2) for x in 0:w, y in 0:w];

julia> f[3:4, 2:6] .= 40; f[3:5, 10:12] .= 60; f[10:12, 3:5] .= 80;

julia> f[10:11, 10:12] .= 100; f[11, 11] = 100;

julia> f_aclose = area_closing(f, min_area=8, connectivity=1);

All small minima are removed, and the remaining minima have at least a size of 8.

area_closing!(output, image, [maxtree];
              min_area=64, connectivity=1) -> output

Performs in-place area closing of the image and stores the result in output. See area_closing for the detailed description of the method.

diameter_opening(image, [maxtree]; min_diameter=8, connectivity=1) -> Array

Performs a diameter opening of the image.

Diameter opening replaces all bright structures of an image that have the diameter (the widest dimension of their bounding box) smaller than min_diameter with the darker value taken from their first ancestral component (in max-tree representation of image) that has the diameter no smaller than min_diameter.

The operator is also called Bounding Box Opening. In practice, the result is similar to a morphological opening, but long and thin structures are not removed.

Arguments

  • image::GenericGrayImage: the -dimensional input image

  • min_diameter::Number=8: the minimal length (in pixels) of the widest dimension of the bounding box of the image component to keep intact

  • connectivity::Integer=1: the neighborhood connectivity. The maximum number of orthogonal steps to reach a neighbor of the pixel. In 2D, it is 1 for a 4-neighborhood and 2 for a 8-neighborhood.

  • maxtree::MaxTree: optional pre-built max-tree. Note that maxtree and connectivity optional parameters are mutually exclusive.

Returns

An array of the same type and shape as the image.

See also

diameter_opening!, diameter_closing, area_opening, MaxTree, opening

References

  1. Walter, T., & Klein, J.-C. (2002). Automatic Detection of Microaneurysms in Color Fundus Images of the Human Retina by Means of the Bounding Box Closing. In A. Colosimo, P. Sirabella, A. Giuliani (Eds.), Medical Data Analysis. Lecture Notes in Computer Science, vol 2526, 210-220. Springer Berlin Heidelberg.

    https://doi.org/10.1007/3-540-36104-9_23

  2. Carlinet, E., & Geraud, T. (2014). A Comparative Review of Component Tree Computation Algorithms. IEEE Transactions on Image Processing, 23(9), 3885-3895.

    https://doi.org/10.1109/TIP.2014.2336551

Examples

Creating a test image f (quadratic function with a maximum in the center and 4 additional local maxima):

julia> w = 12;

julia> f = [20 - 0.2*((x - w/2)^2 + (y-w/2)^2) for x in 0:w, y in 0:w];

julia> f[3:4, 2:6] .= 40; f[3:5, 10:12] .= 60; f[10:12, 3:5] .= 80;

julia> f[10:11, 10:12] .= 100; f[11, 11] = 100;

julia> f_dopen = diameter_opening(f, min_diameter=3, connectivity=1);

The peaks with a maximal diameter of 2 or less are removed. For the remaining peaks the widest side of the bounding box is at least 3.

diameter_opening!(output, image, [maxtree];
                  min_diameter=8, connectivity=1) -> output

Performs in-place diameter opening of the image and stores the result in output. See diameter_opening for the detailed description of the method.

diameter_closing(image, [maxtree]; min_diameter=8, connectivity=1) -> Array

Performs a diameter closing of the image.

Diameter closing replaces all dark structures of an image that have the diameter (the widest dimension of their bounding box) smaller than min_diameter with the brighter value taken from their first ancestral component (in max-tree representation of image) that has the diameter no smaller than min_diameter.

Arguments

  • image::GenericGrayImage: the -dimensional input image

  • min_diameter::Number=8: the minimal length (in pixels) of the widest dimension of the bounding box of the image component to keep intact

  • connectivity::Integer=1: the neighborhood connectivity. The maximum number of orthogonal steps to reach a neighbor of the pixel. In 2D, it is 1 for a 4-neighborhood and 2 for a 8-neighborhood.

  • maxtree::MaxTree: optional pre-built max-tree. Note that maxtree and connectivity optional parameters are mutually exclusive.

Returns

An array of the same type and shape as the image.

See also

diameter_closing!, diameter_opening, area_closing, MaxTree, closing

References

  1. Walter, T., & Klein, J.-C. (2002). Automatic Detection of Microaneurysms in Color Fundus Images of the Human Retina by Means of the Bounding Box Closing. In A. Colosimo, P. Sirabella, A. Giuliani (Eds.), Medical Data Analysis. Lecture Notes in Computer Science, vol 2526, 210-220. Springer Berlin Heidelberg.

    https://doi.org/10.1007/3-540-36104-9_23

  2. Carlinet, E., & Geraud, T. (2014). A Comparative Review of Component Tree Computation Algorithms. IEEE Transactions on Image Processing, 23(9), 3885-3895.

    https://doi.org/10.1109/TIP.2014.2336551

Examples

Creating a test image f (quadratic function with a minimum in the center and 4 additional local minima):

julia> w = 12;

julia> f = [180 + 0.2*((x - w/2)^2 + (y-w/2)^2) for x in 0:w, y in 0:w];

julia> f[3:4, 2:6] .= 40; f[3:5, 10:12] .= 60; f[10:12, 3:5] .= 80;

julia> f[10:11, 10:12] .= 100; f[11, 11] = 100;

julia> f_dclose = diameter_closing(f, min_diameter=3, connectivity=1);

All small minima with a diameter of 2 or less are removed. For the remaining minima the widest bounding box side is at least 3.

diameter_closing!(output, image, [maxtree];
                  min_diameter=8, connectivity=1) -> output

Performs in-place diameter closing of the image and stores the result in output. See diameter_closing for the detailed description of the method.

local_maxima!(output, image, [maxtree]; connectivity=1) -> output

Detects the local maxima of image and stores the result in output. See local_maxima for the detailed description of the method.

local_maxima(image, [maxtree]; connectivity=1) -> Array

Determines and labels all local maxima of the image.

Details

The local maximum is defined as the connected set of pixels that have the same value, which is greater than the values of all pixels in direct neighborhood of the set.

Technically, the implementation is based on the max-tree representation of an image. It’s beneficial if the max-tree is already computed, otherwise ImageFiltering.findlocalmaxima would be more efficient.

Arguments

  • image::GenericGrayImage: the -dimensional input image

  • connectivity::Integer=1: the neighborhood connectivity. The maximum number of orthogonal steps to reach a neighbor of the pixel. In 2D, it is 1 for a 4-neighborhood and 2 for a 8-neighborhood.

  • maxtree::MaxTree: optional pre-built max-tree. Note that maxtree and connectivity optional parameters are mutually exclusive.

Returns

An integer array of the same shape as the image. Pixels that are not local maxima have 0 value. Pixels of the same local maximum share the same positive value (the local maximum id).

See also

MaxTree, local_maxima!, local_minima, ImageFiltering.findlocalmaxima

Examples

Create f (quadratic function with a maximum in the center and 4 additional constant maxima):

julia> w = 10;

julia> f = [20 - 0.2*((x - w/2)^2 + (y-w/2)^2) for x in 0:w, y in 0:w];

julia> f[3:5, 3:5] .= 40; f[3:5, 8:10] .= 60; f[8:10, 3:5] .= 80; f[8:10, 8:10] .= 100;

julia> f_maxima = local_maxima(f); # Get all local maxima of `f`

The resulting image contains the 4 labeled local maxima.

local_minima!(output, image, [maxtree]; connectivity=1) -> output

Detects the local minima of image and stores the result in output. See local_minima for the detailed description of the method.

local_minima(image, [maxtree]; connectivity=1) -> Array

Determines and labels all local minima of the image.

Details

The local minimum is defined as the connected set of pixels that have the same value, which is less than the values of all pixels in direct neighborhood of the set.

Technically, the implementation is based on the max-tree representation of an image. It’s beneficial if the max-tree is already computed, otherwise ImageFiltering.findlocalminima would be more efficient.

Arguments

  • image::GenericGrayImage: the -dimensional input image

  • connectivity::Integer=1: the neighborhood connectivity. The maximum number of orthogonal steps to reach a neighbor of the pixel. In 2D, it is 1 for a 4-neighborhood and 2 for a 8-neighborhood.

  • maxtree::MaxTree: optional pre-built max-tree. Note that maxtree and connectivity optional parameters are mutually exclusive.

Returns

An integer array of the same shape as the image. Pixels that are not local minima have 0 value. Pixels of the same local minimum share the same positive value (the local minimum id).

See also

MaxTree, local_minima!, local_maxima, ImageFiltering.findlocalminima

Examples

Create f (quadratic function with a minimum in the center and 4 additional constant minimum):

julia> w = 10;

julia> f = [180 + 0.2*((x - w/2)^2 + (y-w/2)^2) for x in 0:w, y in 0:w];

julia> f[3:5, 3:5] .= 40; f[3:5, 8:10] .= 60; f[8:10, 3:5] .= 80; f[8:10, 8:10] .= 100;

julia> f_minima = local_minima(f); # Calculate all local minima of `f`

The resulting image contains the labeled local minima.

rebuild!(maxtree::MaxTree, image::GenericGrayImage,
         neighbors::AbstractVector{CartesianIndex}) -> maxtree

Rebuilds the maxtree for the image using neighbors as the pixel connectivity specification.

Details

The pixels in the connected components generated by the method should be connected to each other by a path through neighboring pixels. The pixels and are neighbors, if neighbors array contains , such that .

See also

MaxTree

filter_components!(output::GenericGrayImage, image::GenericGrayImage,
                   maxtree::MaxTree, attrs::AbstractVector,
                   min_attr, all_below_min) -> output

Filters the connected components of the image and stores the result in output.

The is the copy of the exluding the connected components, whose attribute value is below min_attr. That is, the pixels of the exluded component are reset to the value of the reference pixel of its first valid ancestor (the connected component with the attribute value greater or equal to min_attr).

Arguments

  • maxtree::MaxTree: pre-built max-tree representation of the image

  • attrs::AbstractVector: attrs[i] is the attribute value for the -th component of the tree ( being the linear index of its reference pixel)

  • all_below_min: the value to fill the output if all attributes of all components (including the root one) are below min_attr

Details

This function is the basis for area_opening, diameter_opening and similar transformations. E.g. for area_opening the attribute is the area of the components. In this case, the max-tree components of the output have area no smaller than min_attr pixels.

The method assumes that the attribute values are monotone with respect to the components hieararchy, i.e. <= attrs[maxtree.parentindices[i]]] for each i.

Feature transform

feature_transform(img::AbstractArray{Bool, N};
                  weights=nothing, nthreads=Threads.nthreads()) -> F

Compute the feature transform of a binary image I, finding the closest "feature" (positions where I is true) for each location in I. Specifically, F[i] is a CartesianIndex encoding the position closest to i for which I[F[i]] is true. In cases where two or more features in I have the same distance from i, an arbitrary feature is chosen. If I has no true values, then all locations are mapped to an index where each coordinate is typemin(Int).

Optionally specify the weight w assigned to each coordinate. For example, if I corresponds to an image where voxels are anisotropic, w could be the voxel spacing along each coordinate axis. The default value of nothing is equivalent to w=(1,1,...).

See also: distance_transform.

Citation

  • [1] Maurer, Calvin R., Rensheng Qi, and Vijay Raghavan. "A linear time algorithm for computing exact Euclidean distance transforms of binary images in arbitrary dimensions." IEEE Transactions on Pattern Analysis and Machine Intelligence 25.2 (2003): 265-270.

distance_transform(F::AbstractArray{CartesianIndex}, [w=nothing]) -> D

Compute the distance transform of F, where each element F[i] represents a "target" or "feature" location assigned to i. Specifically, D[i] is the distance between i and F[i]. Optionally specify the weight w assigned to each coordinate; the default value of nothing is equivalent to w=(1,1,...).

See also: feature_transform.

cleared_img = clearborder(img)
cleared_img = clearborder(img, width)
cleared_img = clearborder(img, width, background)

Returns a copy of the original image after clearing objects connected to the border of the image. Parameters:

  • img = Binary/Grayscale input image

  • width = Width of the border examined (Default value is 1)

  • background = Value to be given to pixels/elements that are cleared (Default value is 0)

Misc

chull = convexhull(img)

Computes the convex hull of a binary image and returns the vertices of convex hull as a CartesianIndex array.

isboundary(img::AbstractArray; background = 0, dims = coords_spatial(A), kwargs...)

Finds the boundaries that are just within each object. background is the scalar value of the background pixels which will not be marked as boundaries. Keyword arguments are passed to extremefilt! which include dims indicating the dimension(s) over which to discover boundaries.

See also its in-place version isboundary! and the alternative version that finds thick boundaries, isboundary_thick.

Examples

DocTestSetup = quote
    import ImageMorphology: isboundary
end
julia> A = zeros(Int64, 16, 16); A[4:8, 4:8] .= 5; A[4:8, 9:12] .= 6; A[10:12,13:15] .= 3; A[10:12,3:6] .= 9; A
16×16 Matrix{Int64}:
 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
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  5  5  5  5  5  6  6  6  6  0  0  0  0
 0  0  0  5  5  5  5  5  6  6  6  6  0  0  0  0
 0  0  0  5  5  5  5  5  6  6  6  6  0  0  0  0
 0  0  0  5  5  5  5  5  6  6  6  6  0  0  0  0
 0  0  0  5  5  5  5  5  6  6  6  6  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  9  9  9  9  0  0  0  0  0  0  3  3  3  0
 0  0  9  9  9  9  0  0  0  0  0  0  3  3  3  0
 0  0  9  9  9  9  0  0  0  0  0  0  3  3  3  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  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  0

julia> isboundary(A)
16×16 Matrix{Int64}:
 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
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  1  1  1  1  1  1  1  1  1  0  0  0  0
 0  0  0  1  0  0  0  1  1  0  0  1  0  0  0  0
 0  0  0  1  0  0  0  1  1  0  0  1  0  0  0  0
 0  0  0  1  0  0  0  1  1  0  0  1  0  0  0  0
 0  0  0  1  1  1  1  1  1  1  1  1  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  1  1  1  1  0  0  0  0  0  0  1  1  1  0
 0  0  1  0  0  1  0  0  0  0  0  0  1  0  1  0
 0  0  1  1  1  1  0  0  0  0  0  0  1  1  1  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  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  0

julia> isboundary(A .!= 0)
16×16 BitMatrix:
 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
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  1  1  1  1  1  1  1  1  1  0  0  0  0
 0  0  0  1  0  0  0  0  0  0  0  1  0  0  0  0
 0  0  0  1  0  0  0  0  0  0  0  1  0  0  0  0
 0  0  0  1  0  0  0  0  0  0  0  1  0  0  0  0
 0  0  0  1  1  1  1  1  1  1  1  1  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  1  1  1  1  0  0  0  0  0  0  1  1  1  0
 0  0  1  0  0  1  0  0  0  0  0  0  1  0  1  0
 0  0  1  1  1  1  0  0  0  0  0  0  1  1  1  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  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  0

julia> isboundary(A .!= 0; dims = 1)
16×16 BitMatrix:
 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
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  1  1  1  1  1  1  1  1  1  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  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  1  1  1  1  1  1  1  1  1  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  1  1  1  1  0  0  0  0  0  0  1  1  1  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  1  1  1  1  0  0  0  0  0  0  1  1  1  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  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  0

julia> isboundary(A .!= 0; dims = 2)
16×16 BitMatrix:
 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
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  1  0  0  0  0  0  0  0  1  0  0  0  0
 0  0  0  1  0  0  0  0  0  0  0  1  0  0  0  0
 0  0  0  1  0  0  0  0  0  0  0  1  0  0  0  0
 0  0  0  1  0  0  0  0  0  0  0  1  0  0  0  0
 0  0  0  1  0  0  0  0  0  0  0  1  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  1  0  0  1  0  0  0  0  0  0  1  0  1  0
 0  0  1  0  0  1  0  0  0  0  0  0  1  0  1  0
 0  0  1  0  0  1  0  0  0  0  0  0  1  0  1  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  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  0
isboundary!(img::AbstractArray; background = 0, dims = coords_spatial(A), kwargs...)

Finds the boundaries that are just within each object, replacing the original image. background is the scalar value of the background pixels which will not be marked as boundaries. Keyword arguments are passed to extreme_filter which include dims indicating the dimension(s) over which to discover boundaries.

See out-of-place version, isboundary, for examples.

isboundary_thick(img::AbstractArray; dims = coords_spatial(img), kwargs...)

Find thick boundaries that are just outside and just inside the objects. This is a union of the inner and outer boundaries. Keyword dims indicates over which dimensions to look for boundaries. This dims and additional keywords kwargs are passed to extreme_filter.

See also isboundary which just yields the inner boundaries.

Examples

DocTestSetup = quote
    import ImageMorphology: isboundary_thick
end
julia> A = zeros(Int64, 16, 16); A[4:8, 4:8] .= 5; A[4:8, 9:12] .= 6; A[10:12,13:15] .= 3; A[10:12,3:6] .= 9; A 16×16 Matrix{Int64}:  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  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  5  5  5  5  5  6  6  6  6  0  0  0  0  0  0  0  5  5  5  5  5  6  6  6  6  0  0  0  0  0  0  0  5  5  5  5  5  6  6  6  6  0  0  0  0  0  0  0  5  5  5  5  5  6  6  6  6  0  0  0  0  0  0  0  5  5  5  5  5  6  6  6  6  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  9  9  9  9  0  0  0  0  0  0  3  3  3  0  0  0  9  9  9  9  0  0  0  0  0  0  3  3  3  0  0  0  9  9  9  9  0  0  0  0  0  0  3  3  3  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  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  0

julia> isboundary_thick(A) 16×16 BitMatrix:  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  0  0  1  1  1  1  1  1  1  1  1  1  1  0  0  0  0  0  1  1  1  1  1  1  1  1  1  1  1  0  0  0  0  0  1  1  0  0  0  1  1  0  0  1  1  0  0  0  0  0  1  1  0  0  0  1  1  0  0  1  1  0  0  0  0  0  1  1  0  0  0  1  1  0  0  1  1  0  0  0  0  0  1  1  1  1  1  1  1  1  1  1  1  0  0  0  0  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  0  1  1  1  1  1  1  0  0  0  0  1  1  1  1  1  0  1  1  0  0  1  1  0  0  0  0  1  1  0  1  1  0  1  1  1  1  1  1  0  0  0  0  1  1  1  1  1  0  1  1  1  1  1  1  0  0  0  0  1  1  1  1  1  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  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0

julia> isboundary_thick(A) .& (A .!= 0) 16×16 BitMatrix:  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  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  1  1  1  1  1  1  1  1  1  0  0  0  0  0  0  0  1  0  0  0  1  1  0  0  1  0  0  0  0  0  0  0  1  0  0  0  1  1  0  0  1  0  0  0  0  0  0  0  1  0  0  0  1  1  0  0  1  0  0  0  0  0  0  0  1  1  1  1  1  1  1  1  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  1  1  1  1  0  0  0  0  0  0  1  1  1  0  0  0  1  0  0  1  0  0  0  0  0  0  1  0  1  0  0  0  1  1  1  1  0  0  0  0  0  0  1  1  1  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  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  0

julia> isboundary_thick(A) == isboundary(A; background = -1) true

julia> isboundary_thick(A) .& (A .!= 0) == isboundary(A) # inner boundaries true

julia> isboundary_thick(A .!= 0) .& (A .== 0)  == isboundary(A .== 0) # outer boundaries true