add LBP circular implementation

Project.toml

@@ -5,11 +5,13 @@ version = "0.1.0"
 
 [deps]
 ImageCore = "a09fc81d-aa75-5fe9-8630-4744c3626534"
+Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 TiledIteration = "06e1c1a7-607b-532d-9fad-de7d9aa2abac"
 
 [compat]
 ImageCore = "0.8, 0.9"
+Interpolations = "0.12, 0.13"
 StaticArrays = "0.12, 1"
 TiledIteration = "0.2, 0.3"
 julia = "1"

src/LocalBinaryPatterns.jl

@@ -1,6 +1,7 @@
 module LocalBinaryPatterns
 
 using TiledIteration: EdgeIterator
+using Interpolations
 using ImageCore
 using ImageCore.OffsetArrays
 using StaticArrays

src/lbp_original.jl

@@ -1,16 +1,46 @@
 """
     lbp_original(X; kwargs...)
-    lbp_original!(out, X; kwargs...)
+    lbp_original(X, npoints, radius, interpolation=Linear(); kwargs...)
+    lbp_original!(out, X, offsets; kwargs...)
 
 Compute the local binary pattern, the original version, of gray image `X` using 3x3
 neighborhood matrix.
 
+# 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)`.
+
 # Parameters
 
-The parameters control whether and the degree additional encoding passses are used to get
-patterns that are more robust to certain changes, e.g., rotation. The following lists are
-ordered as encoding order. For example, if `rotation=true` and `uniform_degree=2`, then
-rotation encoding will be applied first.
+The parameters control whether and what degree additional encoding passses are used to
+compute more patterns that are more robust/invariant to certain changes, e.g., rotation. The
+following lists are ordered as encoding order. For example, if `rotation=true` and
+`uniform_degree=2`, then rotation encoding will be applied first.
 
 - `rotation=false`: set `true` to generate patterns that are invariant to rotation [3]. For
   example, pattern `0b00001101` is equivalent to `0b01000011` when `rotation=true`.
@@ -32,6 +62,12 @@ 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
@@ -71,6 +107,9 @@ mapped to `0b00001101`. See also Eq.(8) in [2].
 For 3x3 neighborhood matrix, applying rotation-invariant encoding decreases the possible
 number of binary patterns from ``256`` to ``36``.
 
+The interpolation-based version provides more robust result for rotation-invariant pattern,
+see [2,4] for more details.
+
 ## Uniform encoding
 
 Authors of [2] states that certain local binary patterns are fundamental properties of
@@ -88,23 +127,29 @@ because it only has ``2`` bit transitions.
 - [1] T. Ojala, M. Pietikäinen, and D. Harwood, “A comparative study of texture measures with classification based on featured distributions,” _Pattern Recognition_, vol. 29, no. 1, pp. 51–59, Jan. 1996, doi: 10.1016/0031-3203(95)00067-4.
 - [2] 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.
 - [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.
 """
-lbp_original(X::AbstractArray; kwargs...) = lbp_original!(similar(X, UInt8), X; kwargs...)
+function lbp_original(X::AbstractArray; kwargs...)
+    # 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...)
+end
 function lbp_original!(
         out,
-        X::AbstractMatrix{T};
+        X::AbstractMatrix{T},
+        offsets::Tuple;
         rotation::Bool=false,
         uniform_degree::Union{Nothing,Int}=nothing,
         ) where T<:Union{Real,Gray}
-    # nearest interpolation, 3x3 neighborhood
-
-    # 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))
-    offsets = CartesianIndex.(offsets)
+    length(offsets) > 8 && throw(ArgumentError("length(offsets) >= 8 is not supported."))
     outerR = CartesianIndices(X)
-    innerR = first(outerR)+oneunit(first(outerR)):last(outerR)-oneunit(last(outerR))
+    r = CartesianIndex(
+        ceil(Int, maximum(abs.(extrema(first.(offsets))))),
+        ceil(Int, maximum(abs.(extrema(last.(offsets)))))
+    )
+    innerR = first(outerR)+r:last(outerR)-r
 
     # TODO(johnnychen94): use LoopVectorization
     @inbounds for I in innerR
@@ -113,8 +158,8 @@ function lbp_original!(
         # better performance.
         rst = 0
         for i in 1:length(offsets)
-            o = offsets[i]
-            rst += ifelse(gc <= X[I+o], 1, 0) << (i-1)
+            p = I.I .+ offsets[i]
+            rst += ifelse(gc <= _inbounds_getindex(X, p), 1, 0) << (i-1)
         end
         out[I] = rst
     end
@@ -124,9 +169,9 @@ function lbp_original!(
         gc = X[I]
         rst = 0
         for i in 1:length(offsets)
-            o = offsets[i]
-            checkbounds(Bool, X, I+o) || continue
-            rst += ifelse(gc <= X[I+o], 1, 0) << (i-1)
+            p = I.I .+ offsets[i]
+            checkbounds(Bool, X, p...) || continue
+            rst += ifelse(gc <= _inbounds_getindex(X, p), 1, 0) << (i-1)
         end
         out[I] = rst
     end
@@ -140,3 +185,19 @@ function lbp_original!(
 
     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)
+    end
+    lbp_original!(similar(X, UInt8), itp, offsets; kwargs...)
+end

src/utils.jl

@@ -1,4 +1,4 @@
-const SupportedEncodingTypes = Union{UInt8, UInt16}
+const SupportedEncodingTypes = Union{UInt8}
 const _LBP_encoding_table = Dict()
 function build_LBP_encoding_table(::Type{T};
         rotation::Bool,
@@ -58,3 +58,21 @@ function _build_inverse_table(d::Dict{T,<:AbstractVector{T}}) where {T}
     return id
 end
 _freeze(::Type{T}, v::Vector) where T = SVector{typemax(T)-typemin(T)+1}(v)
+
+function _circular_neighbor_offsets(npoints::Integer, radius::Real)
+    radius >= 1 || throw(ArgumentError("radius >= 1.0 is expected."))
+
+    ntuple(npoints) do p
+        θ=2π*(p-1)/npoints
+        pos = (-1, 1) .* radius .* sincos(θ)
+        # rounding the positions for more stable result when applying interpolation
+        round.(pos; digits=8)
+    end
+end
+
+# A helper function to let `interpolation` keyword correctly understands `Degree` inputs
+@inline wrap_BSpline(itp::Interpolations.InterpolationType) = itp
+@inline wrap_BSpline(degree::Interpolations.Degree) = BSpline(degree)
+
+@inline _inbounds_getindex(A, I) = @inbounds A[I...]
+@inline _inbounds_getindex(A::Interpolations.AbstractInterpolation, I) = @inbounds A(I...)

test/lbp_original.jl

@@ -10,10 +10,11 @@
     @test out == ref_out
 
     # ensure default parameter values are not changed accidently
-    @test out == lbp_original(X; rotation=false)
+    @test out == lbp_original(X; rotation=false, uniform_degree=nothing)
 
     fill!(out, 0)
-    lbp_original!(out, X)
+    offsets = ((-1, -1), (0, -1), (1, -1), (-1, 0), (1, 0), (-1, 1), (0, 1), (1, 1))
+    lbp_original!(out, X, offsets)
     @test out == ref_out
 
     # intensity-invariant
@@ -64,3 +65,83 @@
         @test out == ref_out
     end
 end
+
+@testset "lbp_original, Interpolation-based" begin
+    X = [6 7 9; 5 6 3; 2 1 7]
+    ref_out = [1 1 0; 3 2 14; 2 7 0]
+    out = lbp_original(X, 4, 1)
+    @test eltype(out) == UInt8
+    @test size(out) == (3, 3)
+    @test out == ref_out
+
+    X = [6 7 9; 5 6 3; 2 1 7]
+    ref_out = [129 1 0; 7 14 124; 6 31 0]
+    out = lbp_original(X, 8, 1)
+    @test eltype(out) == UInt8
+    @test size(out) == (3, 3)
+    @test out == ref_out
+
+    @testset "OffsetArrays" begin
+        X = [6 7 9; 5 6 3; 2 1 7]
+        ref_out = [129 1 0; 7 14 124; 6 31 0]
+        Xo = OffsetArray(X, -1, -1)
+        out = lbp_original(Xo, 8, 1)
+        @test axes(out) == axes(Xo)
+        @test OffsetArrays.no_offset_view(out) == ref_out
+    end
+
+    # ensure default parameter values are not changed accidently
+    out = lbp_original(X, 8, 1)
+    @test lbp_original(X, 8, 1, Linear(); rotation=false, uniform_degree=nothing) == out
+    @test lbp_original(X, 8, 1, Linear()) == lbp_original(X, 8, 1, BSpline(Linear()))
+
+    # The pattern encoding is different because the offset orders are not the same
+    @test lbp_original(X, 8, 1, Constant()) != lbp_original(X)
+
+    @testset "Rotation Invariant" begin
+        X = [6 7 9; 5 6 3; 2 1 7]
+        ref_out = [3 1 0; 7 7 31; 3 31 0]
+        out = lbp_original(X, 8, 1; rotation=true)
+        @test eltype(out) == UInt8
+        @test size(out) == (3, 3)
+        @test out == ref_out
+    end
+
+    @testset "Uniform encoding" begin
+        X = [6 7 9; 5 6 3; 2 1 7]
+        ref_out = [9 1 0; 7 14 124; 6 31 0]
+        out = lbp_original(X, 8, 1; uniform_degree=2)
+        @test eltype(out) == UInt8
+        @test size(out) == (3, 3)
+        @test out == ref_out
+    end
+
+    @testset "Rotation Invariant, Uniform encoding" begin
+        X = [6 7 9; 5 6 3; 2 1 7]
+        ref_out = [3 1 0; 7 7 31; 3 31 0]
+        out = lbp_original(X, 8, 1; rotation=true, uniform_degree=2)
+        @test eltype(out) == UInt8
+        @test size(out) == (3, 3)
+        @test out == ref_out
+    end
+
+    @testset "radius" begin
+        X = [3  1  1  2  2
+             1  2  1  5  5
+             1  1  2  4  5
+             3  3  5  1  2
+             3  3  5  2  2]
+        ref_out = [0   9   13  8   8
+                   9   9   13  0   0
+                   11  11  9   0   0
+                   1   0   0   6   6
+                   1   0   0   6   6]
+        out = lbp_original(X, 4, 2)
+        @test eltype(out) == UInt8
+        @test size(out) == (5, 5)
+        @test out == ref_out
+
+        @test_nowarn lbp_original(X, 4, 1.5)
+        @test_throws ArgumentError lbp_original(X, 4, 0.5)
+    end
+end

test/runtests.jl

@@ -1,6 +1,7 @@
 using LocalBinaryPatterns
 using ImageCore
 using ImageCore.OffsetArrays
+using Interpolations
 using Test
 using Aqua, Documenter