support up to 32 circular points (#12)

Author: johnnychen94

docs/src/reference.md

@@ -1,5 +1,13 @@
 # API References
 
-```@autodocs
-Modules = [LocalBinaryPatterns]
+```@docs
+lbp_original
+```
+
+## Internal utilities
+
+```@docs
+LocalBinaryPatterns.rotation_encoding_table
+LocalBinaryPatterns.RotationEncodingTable
+LocalBinaryPatterns.uniform_encoding_table
 ```

src/LocalBinaryPatterns.jl

@@ -6,11 +6,12 @@ using ImageCore
 using ImageCore.OffsetArrays
 using StaticArrays
 
-export
-    lbp_original, lbp_original!
+export lbp_original
 
 include("lbp_original.jl")
 
+include("bitrotate_encoding.jl")
+include("uniform_encoding.jl")
 include("utils.jl")
 include("compat.jl")
 

src/bitrotate_encoding.jl

@@ -0,0 +1,104 @@
+const _RotationEncodingTables = Dict()
+
+"""
+    rotation_encoding_table(T, nbits)
+    rotation_encoding_table(T, X::AbstractUnitRange)
+
+Build the rotation-invariant encoding table. If global runtime cache exists, then reuse it.
+"""
+rotation_encoding_table(::Type{T}, nbits::Int) where T = rotation_encoding_table(T, 0:2^nbits-1)
+function rotation_encoding_table(::Type{T}, X::AbstractUnitRange) where T
+    d = _RotationEncodingTables
+    k = (T, minimum(X), maximum(X))
+    haskey(d, k) && return d[k]
+    lookup = d[k] = RotationEncodingTable{T}(X)
+    return lookup
+end
+
+"""
+    RotationEncodingTable{T<:Unsigned}(X)
+
+The lookup table for the quotient space constructed using `bitrotate` equivalance relation
+on bit representation of datatype `T`.
+
+By definition, ``a`` is equivalant to ``b``(``a ~ b``) if there exists ``n`` such that
+`bitrotate(a, n) == b`.
+
+This equivalance class gives a unique encoding map from ``x∈X`` to ``a``, where ``a`` is the
+minimal value of the equivalance class ``[x]``. For instance, the equivalance class of
+`0b10101001` (UInt8) is
+
+```julia
+0b10101001 # 169
+0b01010011 # 83
+0b10100110 # 166
+0b01001101 # 77
+0b10011010 # 154
+0b00110101 # 53
+0b01101010 # 106
+0b11010100 # 212
+```
+
+thus `f(169) == 0x35 == 53`.
+
+```jldoctest
+julia> using LocalBinaryPatterns: RotationEncodingTable
+
+julia> X = RotationEncodingTable{UInt8}(0:255);
+
+julia> X[169] # 53
+0x35
+
+julia> X = RotationEncodingTable{UInt8}(128:255);
+
+julia> X[169] # 154
+0x9a
+```
+
+The values in the lookup table is implemented in lazy cache manner so the real computation
+only happens once at the first retrieval.
+
+This lookup table is used to build an efficient implementation of rotation-invariant local
+binary pattern [1].
+
+# References
+
+- [1] T. Ojala, M. Pietikäinen, and T. Mäenpää, “A Generalized Local Binary Pattern Operator for Multiresolution Gray Scale and Rotation Invariant Texture Classification,” in _Advances in Pattern Recognition — ICAPR 2001, vol. 2013, S. Singh, N. Murshed, and W. Kropatsch, Eds. Berlin, Heidelberg: Springer Berlin Heidelberg_, 2001, pp. 399–408. doi: 10.1007/3-540-44732-6_41.
+"""
+struct RotationEncodingTable{T<:Unsigned,AT<:AbstractUnitRange,LAT} <: AbstractVector{T}
+    X::AT
+    Y::LAT
+end
+function RotationEncodingTable{T}(X::AbstractUnitRange) where T
+    if !((maximum(X) <= typemax(T)) && (minimum(X) >= typemin(T)))
+        throw(ArgumentError("range $(first(X)):$(last(X)) can't be uniquely represented by type $T."))
+    end
+
+    if log2(length(X)) >= 26
+        datasize = round(length(X)/1024/1024/1024, digits=3)
+        @warn "$(datasize)GB data will be created."
+    end
+
+    Y = Vector{Union{T,Missing}}(missing, length(X))
+    RotationEncodingTable{T,typeof(X),typeof(Y)}(X, Y)
+end
+
+@inline Base.axes(X::RotationEncodingTable) = (first(X.X):last(X.X), )
+@inline Base.size(X::RotationEncodingTable) = (length(X.X)-1)
+@inline function Base.getindex(X::RotationEncodingTable{T}, x::Int) where T
+    i = x-first(axes(X, 1))+1
+    v = X.Y[i]
+    if ismissing(v)
+        v = minimum_bitrotate(T(x), X.X)
+        X.Y[i] = v
+    end
+    return T(v)
+end
+
+function minimum_bitrotate(x::T, X::AbstractUnitRange) where T<:Unsigned
+    function bit_rotate_project_to(x, n, X)
+        v = bitrotate(x, n)
+        return v in X ? v : x
+    end
+    minimum(bit_rotate_project_to(x, n, X) for n in 0:8sizeof(T)-1)
+end

src/lbp_original.jl

@@ -1,38 +1,13 @@
 """
-    lbp_original(X; kwargs...)
-    lbp_original(X, npoints, radius, interpolation=Linear(); kwargs...)
-    lbp_original!(out, X, offsets; kwargs...)
+    lbp_original(X; [rotation], [uniform_degree])
 
 Compute the local binary pattern of gray image `X` using the original method.
 
 # Arguments
 
 - `X::AbstractMatrix`: the input image matrix. For colorful images, one can manually convert
-  it to some monochrome space, e.g., `Gray` or the L-channel of `Lab`.
-
-[2,4] proposes a generalized interpolation-based version using circular neighborhood matrix,
-it produces better result for `rotation=true` case but is usually much slower than the plain
-3x3 matrix version. The arguments for this version are:
-
-- `npoints::Int`(4 ≤ npoints ≤ 8): the number of (uniform-spaced) neighborhood points.
-- `radius::Real`(radius ≥ 1.0): the radius of the circular.
-- `interpolation::Union{Degree, InterpolationType}=Linear()`: the interpolation method used
-  to generate non-grid pixel value. In most cases, `Linear()` are good enough. One can also
-  try other costly interpolation methods, e.g., `Cubic(Line(OnGrid()))`(also known as
-  "bicubic"), `Lanczos()`. See also Interpolations.jl for more choices.
-
-!!! note "neighborhood order differences"
-    Different implementation might use different neighborhood orders; this will change the
-    encoding result but will not change the overall distribution. For instance,
-    `lbp_original(X)` differs from `lbp_original(X, 8, 1, Constant())` only by how `offsets`
-    (see below) are ordered; the former uses column-major top-left to bottom-right 3x3 matrix
-    order and the latter uses circular order.
-
-Arguments for in-place version:
-
-- `offsets::Tuple`: tuple of neighborhood matrix, e.g., `((0, 1), (0, -1), (1, 0), (-1, 0))`
-  specifies the 4-neighborhood matrix. If `X isa Interpolations.AbstractInterpolation`
-  holds, then the position values can be float numbers, e.g, `(0.7, 0.7)`.
+  it to some monochrome space, e.g., `Gray`, the L-channel of `Lab`. One could also do
+  channelwise LBP and then concatenate together.
 
 # Parameters
 
@@ -61,12 +36,6 @@ julia> lbp_original(X)
  0x68  0xa9  0x1b
  0x28  0x6b  0x00
 
-julia> lbp_original(X, 4, 1) # 4-neighbor with circular radius 1
-3×3 $(Matrix{UInt8}):
- 0x01  0x01  0x00
- 0x03  0x02  0x0e
- 0x02  0x07  0x00
-
 julia> lbp_original(X; rotation=true)
 3×3 $(Matrix{UInt8}):
  0x03  0x01  0x00
@@ -128,21 +97,90 @@ unchanged because it only has ``2`` bit transitions.
 - [3] Pietikäinen, Matti, Timo Ojala, and Zelin Xu. "Rotation-invariant texture classification using feature distributions." _Pattern recognition_ 33.1 (2000): 43-52.
 - [4] T. Ojala, M. Pietikainen, and T. Maenpaa, “Multiresolution gray-scale and rotation invariant texture classification with local binary patterns,” _IEEE Trans. Pattern Anal. Machine Intell._, vol. 24, no. 7, pp. 971–987, Jul. 2002, doi: 10.1109/TPAMI.2002.1017623.
 """
-function lbp_original(X::AbstractArray; kwargs...)
+function lbp_original(X::AbstractArray; rotation::Bool=false, uniform_degree::Union{Nothing,Int}=nothing)
     # The original version [1] uses clockwise order; here we use anti-clockwise order
     # because Julia is column-major order. If we consider memory order differences then they
     # are equivalent.
     offsets = ((-1, -1), (0, -1), (1, -1), (-1, 0), (1, 0), (-1, 1), (0, 1), (1, 1))
-    lbp_original!(similar(X, UInt8), X, offsets; kwargs...)
+    lookups = _build_lbp_original_lookup(UInt8, 8; rotation=rotation, uniform_degree=uniform_degree)
+    lbp_original!(similar(X, UInt8), X, offsets, lookups)
+end
+
+"""
+    lbp_original(X, npoints, radius, interpolation=Linear(); [rotation], [uniform_degree])
+
+Compute the local binary pattern of gray image `X` using the interpolation-based original
+method with circular neighborhood matrix.
+
+This produces better result for `rotation=true` case but is usually slower than the plain
+3x3 matrix version `lbp_original(X)`.
+
+# Arguments
+
+- `npoints::Int`(4 ≤ npoints ≤ 32): the number of (uniform-spaced) neighborhood points. It
+    is recommended to use one of {4, 8, 12, 16, 24}.
+- `radius::Real`(radius ≥ 1.0): the radius of the circular. Larger radius computes the
+    pattern of a larger local window/block.
+- `interpolation::Union{Degree, InterpolationType}=Linear()`: the interpolation method used
+    to generate non-grid pixel value. In most cases, `Linear()` are good enough. One can
+    also try other costly interpolation methods, e.g., `Cubic(Line(OnGrid()))`(also known as
+    "bicubic"), `Lanczos()`. See also Interpolations.jl for more choices.
+
+!!! note "neighborhood order differences"
+    Different implementation might use different neighborhood orders; this will change the
+    encoding result but will not change the overall distribution. For instance,
+    `lbp_original(X)` differs from `lbp_original(X, 8, 1, Constant())` only by how
+    `offsets` (see below) are ordered; the former uses column-major top-left to
+    bottom-right 3x3 matrix order and the latter uses circular order.
+
+# Examples
+
+```jldoctest; setup=:(using LocalBinaryPatterns)
+julia> X = [6 7 9; 5 6 3; 2 1 7]
+3×3 $(Matrix{Int}):
+ 6  7  9
+ 5  6  3
+ 2  1  7
+
+julia> lbp_original(X, 4, 1) # 4-neighbor with circular radius 1
+3×3 $(Matrix{UInt32}):
+ 0x00000001  0x00000001  0x00000000
+ 0x00000003  0x00000002  0x0000000e
+ 0x00000002  0x00000007  0x00000000
+
+julia> lbp_original(X, 4, 1; rotation=true)
+3×3 $(Matrix{UInt32}):
+ 0x00000001  0x00000001  0x00000000
+ 0x00000003  0x00000001  0x00000007
+ 0x00000001  0x00000007  0x00000000
+```
+
+"""
+function lbp_original(
+        X::AbstractArray, npoints, radius, interpolation=Linear();
+        rotation::Bool=false, uniform_degree::Union{Nothing,Int}=nothing)
+    interpolation = wrap_BSpline(interpolation)
+    offsets = _circular_neighbor_offsets(npoints, radius)
+    if interpolation == BSpline(Constant())
+        offsets = map(offsets) do o
+            round.(Int, o, RoundToZero)
+        end
+        itp = X
+    else
+        itp = interpolate(X, interpolation)
+    end
+    # For the sake of type-stability, hardcode to UInt32 here.
+    lookups = _build_lbp_original_lookup(UInt32, npoints, rotation=rotation, uniform_degree=uniform_degree)
+    lbp_original!(similar(X, UInt32), itp, offsets, lookups)
 end
+
 function lbp_original!(
         out,
         X::AbstractMatrix{T},
-        offsets::Tuple;
-        rotation::Bool=false,
-        uniform_degree::Union{Nothing,Int}=nothing,
+        offsets::Tuple,
+        lookups::Vector
         ) where T<:Union{Real,Gray}
-    length(offsets) > 8 && throw(ArgumentError("length(offsets) >= 8 is not supported."))
+    length(offsets) > 32 && throw(ArgumentError("length(offsets)=$(length(offsets)) is expected to be no larger than 32."))
     outerR = CartesianIndices(X)
     r = CartesianIndex(
         ceil(Int, maximum(abs.(extrema(first.(offsets))))),
@@ -175,28 +213,26 @@ function lbp_original!(
         out[I] = rst
     end
 
-    # The building is cached and chained(if there are multiple encoding passes) thus the
-    # cost is decreased to one brocasting to `getindex` and `setindex!`.
-    encoding_table = build_LBP_encoding_table(UInt8; rotation=rotation, uniform_degree=uniform_degree)
-    if !isnothing(encoding_table)
-        @. out = encoding_table[out + 1]
+    for F in lookups
+        @. out = F[out]
     end
 
     return out
 end
 
-function lbp_original(X::AbstractArray, npoints, radius, interpolation=Linear(); kwargs...)
-    # TODO(johnnychen94): support npoints=24 as [2, 4] indicate
-    4<= npoints <= 8 || throw(ArgumentError("currently only support 4<= npoints <=8, instead it is $(npoints)."))
-    interpolation = wrap_BSpline(interpolation)
-    offsets = _circular_neighbor_offsets(npoints, radius)
-    if interpolation == BSpline(Constant())
-        offsets = map(offsets) do o
-            round.(Int, o, RoundToZero)
-        end
-        itp = X
-    else
-        itp = interpolate(X, interpolation)
+function _build_lbp_original_lookup(
+        ::Type{T}, nbits;
+        rotation::Bool,
+        uniform_degree::Union{Nothing,Int}
+    ) where T <: Unsigned
+    # TODO(johnnychen94): fuse these operation into one broadcasting could potentially give
+    # better performance.
+    lookups = []
+    if rotation
+        push!(lookups, rotation_encoding_table(T, nbits))
+    end
+    if !isnothing(uniform_degree)
+        push!(lookups, uniform_encoding_table(T, nbits, uniform_degree))
     end
-    lbp_original!(similar(X, UInt8), itp, offsets; kwargs...)
+    return lookups
 end

src/uniform_encoding.jl

@@ -0,0 +1,62 @@
+const _UniformEncodingTables = Dict()
+"""
+    uniform_encoding_table(T, X::AbstractUnitRange, degree::Int)
+
+Get the encoding table of `X` using uniform degree filter.
+
+A uniform degree filter is defined as `f(x) = U(x) > degree ? P+1 : x`, where `P` is number
+of bits of `x`, and `U(x)` is the circular bit transition of `x` with bit representation of
+datatype `T`.
+
+For example, `0b0000001` has ``2`` bit transitions thus is unchanged, while `0b00000101` has
+``4`` bit transitions. The transition is count in circular sense, i.e., the highest and
+lowest bits are also compared.
+
+```jldoctest
+julia> using LocalBinaryPatterns: uniform_encoding_table
+
+julia> X = uniform_encoding_table(UInt8, 0:255, 2);
+
+julia> X[0b0000001]
+0x01
+
+julia> X[0b0000101] # miscellaneous
+0x09
+```
+
+This function is used to distinguish local binary patterns from texture patterns and
+miscellaneous patterns. See [1] for more information.
+
+!!! info "Runtime caches"
+    For better performance, this function uses a runtime cache to store the lookup and
+    shared among all callers. The result is expected to be used in read-only mode.
+
+# References
+
+- [1] T. Ojala, M. Pietikäinen, and T. Mäenpää, “A Generalized Local Binary Pattern Operator for Multiresolution Gray Scale and Rotation Invariant Texture Classification,” in _Advances in Pattern Recognition — ICAPR 2001, vol. 2013, S. Singh, N. Murshed, and W. Kropatsch, Eds. Berlin, Heidelberg: Springer Berlin Heidelberg_, 2001, pp. 399–408. doi: 10.1007/3-540-44732-6_41.
+"""
+function uniform_encoding_table(::Type{T}, X::AbstractUnitRange, degree::Int) where T
+    d = _UniformEncodingTables
+    k = (T, degree, minimum(X), maximum(X))
+    haskey(d, k) && return d[k]
+    lookup = d[k] = _uniform_encoding_table(T, X, degree)
+    return lookup
+end
+uniform_encoding_table(::Type{T}, nbits::Int, degree::Int) where T = uniform_encoding_table(T, 0:2^nbits-1, degree)
+
+function _uniform_encoding_table(::Type{T}, X::AbstractUnitRange, degree::Int) where T<:Unsigned
+    if !((maximum(X) <= typemax(T)) && (minimum(X) >= typemin(T)))
+        throw(ArgumentError("range $(first(X)):$(last(X)) can't be uniquely represented by type $T."))
+    end
+
+    function _count_bit_transitions(x)
+        # Eq. (10) in Ojala 2002
+        count_ones(x << 7 & typemax(typeof(x))) + count_ones(x ⊻ (x >> 1))
+    end
+    # Eq. (9) for Ojala 2002
+    lookup = map(X) do x
+        n = _count_bit_transitions(T(x))
+        T(ifelse(n > degree, 8sizeof(T)+1, x))
+    end
+    OffsetVector(lookup, X)
+end