foundamental support for structuring element concept

[johnnychen94]

This PR tries to unify the structuring element concept via strel function and a few internal trait types. For the proof of concept, will work on #11 and #63 on top of this.

This function converts from one structuring element representation to
the other. Currently, three structuring element representations are
supported:

• SEOffsets: Vector{CartesianIndex{N}}
• SEMask: BitArray{Bool}
  - SEAlias

MorphologySE and the above subtypes are internal types. For user
interface, they are:

• Offsets: strel(CartesianIndex, se)
• Mask: strel(Bool, se)
  - Alias: strel([T], symbol, N)

cc: @ThomasRetornaz

[johnnychen94]

A second thought: the SEAlias seems a bad design because it seems quite inconvenient to use in practice.

A better idea is perhaps introducing two new types:

• SEDiamond: e.g., Bool[0 1 0; 1 1 1; 0 1 0]
• SEWindow: e.g., trues(3, 3)

These two are special cases of the above SEOffsets and SEMask, but they can be used for dispatch purposes to provide optimized implementations.

[johnnychen94]

This seems to be quite convenient to use so I think it's stable.

In a summary, this PR introduces a few helper functions and types to properly support the structuring element concept:

• strel_type(x) to infer the SE type of x, where x can be either
  • AbstractArray{Bool} => SEMask
  • AbstractArray{<:CartesianIndex} => SEOffsets
• strel(T, x) function to convert between different SE types
  • strel(CartesianIndex, x): convert x to SEOffsets representation
  • strel(Bool, x): convert x to mask representation
• strel_size helper to infer the size of the minimal enclosing connectivity mask: this can be very useful to infer the radius of given SE

To support multiple dispatches on special SE cases, two generators and corresponding types are provided:

• SEDiamond/SEDiamondArray: can be generated from strel_diamond. For instance, strel_diamond((3, 3)) is C4 connectivity in 2D and strel_diamond((3, 3, 3)) is C6 connectivity in 3D.
• SEWindow/SEWindowArray: can be generated from strel_window. For instance, strel_window((3, 3)) is C8 connectivity in 2D and strel_window((3, 3, 3)) is C26 connectivity in 3D.

---

To develop generic morphology operation without lossing the opportunities to provide optimized version, we can use the following dispatching pattern:

f(img, Ω=strel_diamond(img)) = _f(strel_type(Ω), img, Ω)
f(img, dims::Dims) = f(img, strel_diamond(img, dims))

function _f(::MorphologySE, img, Ω)
... # generic implementation
end

function _f(::SEDiamond, img, Ω)
... # optimized implementation for diamond shape SE
end

function _f(::SEDiamond, img::AbstractArray{Bool}, Ω)
... # optimized implementation for diamond shape SE and Bool array
end

---

Before merging this PR, I'll first verify this design by solving #11 and upgrading #60 to see if this works well.

[ThomasRetornaz]

Why? the list of offset could contains the center, sometimes we use it sometimes not

[ThomasRetornaz]

Awesome :) my level of julia is not at the level to properly review this PR-> I will reads some docs
But i know how i could use it :)
I just have an objection on offset array representation: why remove the center (ok in all my PR i skeep it) but its not true in all cases
Many thanks for your help
Hope the situation will become less stressful on your side

[ThomasRetornaz]

the center is just a displacement of (0,....)

[ThomasRetornaz]

great, we could speed up some process by distinguish border point than non border point along iteration

[ThomasRetornaz]

this surprise me but ok

[johnnychen94]

I just have an objection on offset array representation: why remove the center (ok in all my PR i skeep it) but its not true in all cases

I feel the center pixel is usually handled separately, for instance, to get the local maxima using the offset representation, it is:

function _maxima(A, p::CartesianIndex, offsets)
    out = A[p] # center pixel is retrived separately.
    for o in offsets
        out = max(out, A[p+o])
    end
    return out
end

If we include the center pixel by default, then you'll have to 1) either do one more max operation (but this might not apply for generic function), or 2) filter out the center pixel before the loop.

Maybe I should add a keyword strel(T, x; include_center=false) to allow the center point included in the SE.

[ThomasRetornaz]

Ok i just thinking of algorithm which need the center INSIDE the Neighbordlist maybee some skeleton one
But in this case i probably made a custom algo because skeleton have special cases
So i m with you on this :)
Thanks for your time

[johnnychen94]

Changes summary:

• rebase on master
• a few interface changes for convenience
• rename SEOffsets to SEOffset
• add a detailed explanation doc page for the SE concept and how to use it in ImageMorphology

Going to merge this in a few days if without further comments. This is a long PR, but most of the contents are documentation and test.

[johnnychen94]

@andrewjradcliffe Thank you for improving the docs!

[johnnychen94]

Rebased. Plan to do a merge commit without squash later today.

[johnnychen94]

slack message from @timholy

my only concern is that strel is kind of an ugly name. But I can go with it if you prefer it!

I disliked this name at first glance because I wasn't aware of the concept of "structuring element" but when I started to play around with it to generalize erode function #11,
I feel quite comfortable typing strel, strel_type, and others. For new users, my best bet is that they will feel the same.

To quickly introduce new users to this concept and make them familiar with strel, I decided to set up the documentation with a page docs/src/structuring_element.md explaining this.
Hopefully, that'll be good reading to close the cognition gap.

[johnnychen94]

As a proof of concept, #72 introduces extreme_filter, the new extremefilt! implementation based on strel work in this PR.

[timholy]

I like it but I think there are a few conceptual issues we should consider:

• how to handle anisotropic imaging where individual dimensions have different scaling? This is common in biomedical 3d imaging, where the 3rd dimension often has reduced resolution. The natural metric on such spaces is something likeΔx^2/sx^2 + Δy^2/sy^2 +Δz^2/sz^2 but this does not fit well with notions of "radius" which are used heavily here unless you also specify (sx, sy, sz).
• from a conceptual standpoint, AbstractArray{Bool} elements should probably have their axes centered on zero. This is also consistent with how imfilter works, and of course many morphological operations are equivalent to filtering (e.g., dilate is just mapwindow(maximum, kernel).
• we should think about color images. https://www.hindawi.com/journals/mpe/2012/678326/ has a nice overview. I favor the "RO" approach, where the user can pass in an f mapping colors to real values and then ordering operations are performed on this real value.

[timholy]

Structuring element is definitely a name we want to introduce (it's very widespread), but in many operations they are simply a binary kernel (e.g., as used by imfilter). I think that's worth pointing out.

[johnnychen94]

I'll add another example to illustrate the difference between extreme_filter and mapwindow in #72, if that's what you're expecting.

[timholy]

What is dims for? Is it a list of dimensions over which this operates? How does that differ from having size[d] = 0? Why is r=1 the default rather than r[d] = size[d] ÷ 2? For anisotropic imaging (commonplace in 3d biomedical imaging), presumably you might want to look several pixels in the high-resolution axes but only at nearest-neighbors in the low-resolution axes.

[johnnychen94]

For the r::Int decision, this is mainly because I'm lazy... because if we extend r to be a tuple, implementing an efficient algorithm to check if a given point is within the N-dimensional diamond(polytope) isn't something trivial.

[johnnychen94]

Separating out dims allows me to generalize the definition a bit to skip certain dimensions as extremefilt! allows, yet without the need to play with the is-point-in-polytope algorithm.

[timholy]

Personally I'd rather think of this as

1
1
1

and specify it as strel_window((3,1)). Or in terms of half-sizes, strel_window((1,0)).

[johnnychen94]

I would think both are equivalent representations so I won't be too worried about it. It totally depends on how the user creates it.

Actually, if we play with the strel conversion, we'll see they're eventually unified to the same:

julia> strel(Bool, strel(CartesianIndex, strel_box((3,1))))
3×1 BitMatrix:
 1
 1
 1

julia> strel(Bool, strel(CartesianIndex, strel_box((3,3), (1,))))
3×1 BitMatrix:
 1
 1
 1

But I'll make the default strel_box(img, dims) to use the smaller version:

  function strel_box(img::AbstractArray{T,N}, dims=coords_spatial(img); kw...) where {T,N}
-     return strel_box(ntuple(i->3, N), dims; kw...)
+     return strel_box(ntuple(i->in(i, dims) ? 3 : 1, N), dims; kw...)
  end

[johnnychen94]

The new commits address some of the comments from @timholy

• rename strel_window to strel_box and simplify its usage
• assume the connectivity mask is zero-centered. For one-based indexing masks, a permanent deprecation is given, and for other cases errors will be throw.
• more documentation is added to the trait types

To relieve me from the maintenance burden, I decided to drop support for Julia < 1.6.

---

how to handle anisotropic imaging where individual dimensions have different scaling? This is common in biomedical 3d imaging, where the 3rd dimension often has reduced resolution. The natural metric on such spaces is something likeΔx^2/sx^2 + Δy^2/sy^2 +Δz^2/sz^2 but this does not fit well with notions of "radius" which are used heavily here unless you also specify (sx, sy, sz).

The current strel_diamond and strel_window only meant to provide some room for optimization using the trait trick. If you check the #72 benchmark, you'll notice that it's perhaps not that bad to use a generic mask SE or displacement offsets SE.
Thus the answer for this is to build some simple function to generate the needed SE.

We might want to introduce a way to cartesian product two low-dimensional SEs into a high dimensional one just like ImageFiltering.KernelFactors does. But I'll leave it to future work.

we should think about color images.

Yes, but it's orthogonal to the structuring element as far as I understand. The color reduced-order provides a "meaningful" way to define the max/min operation for colors, but that's not relevant to SE.

---

This becomes a long PR. Since this design works quite well on #72, I'm going to merge this. If there are more comments, I'll try to address them in PRs.

[timholy]

Thanks for adding what is sure to become well used functionality!

[ThomasRetornaz]

@johnnychen94 pretty good job :)
Do you think its mature enough? If true, i will rebase my PR and update my functions ?(i will add func! variant also)
Do you prefer i try to optimize extreme_filter first for some weel known cases (eg 2D or 3D, with diamond/box strel)?
FYI i port code from http://fulguro.sourceforge.net/

[johnnychen94]

Do you think its mature enough? If true, i will rebase my PR and update my functions ?(i will add func! variant also)

@ThomasRetornaz Yes! I think this becomes quite convenient to use based on my extreme_filter development in #72. Since this is a new feature without breaking old codes, I can tag a release if you want.

Do you prefer i try to optimize extreme_filter first for some weel known cases (eg 2D or 3D, with diamond/box strel)?

I've already made a few optimizations for the diamond shape and for the Boolean array case. Not sure if it's the best version, though. If you want to take it a try, it would be fantastic!
I'm currently working on the docs and the extremfilt! deprecation.

FYI i port code from http://fulguro.sourceforge.net/

I noticed that this library uses GPL licenses. Thus if you port the code here, you'll have to clearly keep the license. And there will need to be some extra work if we want to make ImageMorphology a pure MIT license.
It's okay if you want to port the codes, just that I personally prefer a clean implementation.

[ThomasRetornaz]

Yes maybee a tag could be useful
Don't be afraid of the licence, i just take the ideas of fulguro (i have already port this code to C++ in over organisations) and the algorithm were published in hardware aware conferences.
Yes it will pure new code in pure julia (not a wrapper around fulguro)
Have a good day o/

[johnnychen94]

@ThomasRetornaz I just asked this in slack, and unfortunately, this part of code, if ported, should still be GPL code unless you have a clean implementation according to the papers.

As far as I've concerned, it's okay to have a portion of the codebase being GPL licensed, but here we just need to be clear about the license to avoid unnecessary issues in the future.

From @vchuravy slack handle: Valentin (NOT vchuravy)

looking at code & translating it without the cleanroom dance would require GPL, yes

ironically, a wrapper can be whatever license, since it links dynamically

obligatory not a lawyer, this isn't legal advice, yada yada

if they're implementing off of the papers, it should be fine, if they're referencing the GPL C++ for how to write the code, the copyleft nature applies

(as long as the papers don't have the exact same code verbatim...)

well, I also don't think legal action is a concern. It can just lead to really weird & bad situations when others assume something is MIT and then it turns out it wasn't

at least that's my reasoning for why I try to avoid GPL and put everything I don't plan to make money off of out there with MIT

[ThomasRetornaz]

@johnnychen94 thanks i will implement algorithms with this in mind