implement LBP rotation-invariant version

.gitignore

@@ -3,3 +3,4 @@
 *.jl.mem
 Manifest.toml
 /docs/build/
+/benchmark/tune.json

Project.toml

@@ -5,10 +5,12 @@ version = "0.1.0"
 
 [deps]
 ImageCore = "a09fc81d-aa75-5fe9-8630-4744c3626534"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 TiledIteration = "06e1c1a7-607b-532d-9fad-de7d9aa2abac"
 
 [compat]
 ImageCore = "0.8, 0.9"
+StaticArrays = "0.12, 1"
 TiledIteration = "0.2, 0.3"
 julia = "1"
 

benchmark/benchmarks.jl

@@ -14,7 +14,9 @@ tst_types = (N0f8, Float32, Gray{N0f8}, Gray{Float32})
 
 const SUITE = BenchmarkGroup()
 
-alg_list = (( "Original", lbp_original),)
+alg_list = (( "Original", lbp_original),
+            ( "Rotation-Invariant", lbp_ri),
+            )
 
 function add_algorithm_benchmark!(suite, img, alg_name, alg;
                                   tst_sizes, tst_types)

src/LocalBinaryPatterns.jl

@@ -2,9 +2,16 @@ module LocalBinaryPatterns
 
 using TiledIteration: EdgeIterator
 using ImageCore
+using ImageCore.OffsetArrays
+using StaticArrays
 
-export lbp_original, lbp_original!
+export
+    lbp_original, lbp_original!,
+    lbp_ri, lbp_ri!
 
 include("lbp_original.jl")
+include("lbp_ri.jl")
+include("utils.jl")
+include("compat.jl")
 
 end

src/compat.jl

@@ -0,0 +1,6 @@
+if VERSION < v"1.5"
+    # https://github.com/JuliaLang/julia/blob/70cc57cb36d839afe6ce56ea48ff6ed01bc262c4/base/int.jl#L524-L551
+    function bitrotate(x::T, k::Integer) where {T <: Integer}
+        (x << ((sizeof(T) << 3 - 1) & k)) | (x >>> ((sizeof(T) << 3 - 1) & -k))
+    end
+end

src/lbp_original.jl

@@ -2,7 +2,7 @@
     lbp_original(X)
     lbp_original!(out, X)
 
-Compute the local binary pattern, the original version [1], of gray image using 3x3
+Compute the local binary pattern, the original version, of gray image using 3x3
 neighborhood matrix.
 
 ```text
@@ -12,6 +12,8 @@ neighborhood matrix.
 2  1  7         1  0  1              4  16 128         0   0  128
 ```
 
+# Examples
+
 ```jldoctest; setup=:(using LocalBinaryPatterns)
 julia> X = [6 7 9; 5 6 3; 2 1 7]
 3×3 $(Matrix{Int}):
@@ -32,7 +34,7 @@ julia> lbp_original(X)
 - [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.
 """
 lbp_original(X::AbstractArray) = lbp_original!(similar(X, UInt8), X)
-function lbp_original!(out, X::AbstractMatrix{T}) where T<:Union{Real, Gray}
+function lbp_original!(out, X::AbstractMatrix{T}) where T<:Union{Real,Gray}
     # nearest interpolation, 3x3 neighborhood
 
     # The original version [1] uses clockwise order; here we use anti-clockwise order

src/lbp_ri.jl

@@ -0,0 +1,49 @@
+"""
+    lbp_ri(X)
+    lbp_ri!(out, X)
+
+Compute the local binary pattern, the rotation invariant version, of gray image using 3x3
+neighborhood matrix.
+
+```text
+3x3 block     center-thresholded     weights         multiplied by weights      sum      rotation invariant encoding
+6  7  9         1  1  1              1  8  32          1   8  32
+5  6  3  ==>    0  x  0     ==>      2  x  64  ==>     0   x  0            ==>  169  ==>    53
+2  1  7         1  0  1              4  16 128         0   0  128
+```
+
+The rotation encoding is to map all elements in the bitrotation equivalent class to the
+minimal value of this class. For example, `0b11010000` and `0b01000011` belongs to the same
+class because `bitrotate(0b01000011, -2) == 0b11010000`, thus both values are mapped to
+`0b00001101`, which is `minimum(bitrotate.(0b11010000, 0:7))`.
+
+# 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_ri(X)
+3×3 $(Matrix{UInt8}):
+ 0x03  0x01  0x00
+ 0x0d  0x35  0x1b
+ 0x05  0x5b  0x00
+```
+
+# References
+
+- [1] Pietikäinen, Matti, Timo Ojala, and Zelin Xu. "Rotation-invariant texture classification using feature distributions." _Pattern recognition_ 33.1 (2000): 43-52.
+- [2] 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_ri(X::AbstractArray) = lbp_ri!(similar(X, UInt8), X)
+
+function lbp_ri!(out, X::AbstractMatrix{T}) where T<:Union{Real, Gray}
+    # nearest interpolation, 3x3 neighborhood, rotation invariant
+    encoding_table = build_circular_invariant_encoding_table(UInt8)
+    lbp_original!(out, X)
+    @. out = encoding_table[out + 1]
+    return out
+end

src/utils.jl

@@ -0,0 +1,35 @@
+# This implements the solver for Eq. (8) in the Ojala 2002 as lookup table.
+# Computing the encoding table using naive implementation is time-consuming,
+# since it's read-only, we simply cache the encoding table.
+const _circular_invariant_encoding_tables = Dict()
+const SupportedEncodingTypes = Union{UInt8, UInt16}
+function build_circular_invariant_encoding_table(::Type{T}) where T<:SupportedEncodingTypes
+    d = _circular_invariant_encoding_tables
+    haskey(d, T) && return d[T]
+    d[T] = _freeze(T, _build_inverse_table(_bitrotate_quotation_space(T)))
+    return d[T]
+end
+
+# Mathematically, this is the quotation space under circular bitshift of the N-bits binary
+# pattern space. For instance, 0b00001101 and 0b01000011 belong to the same equivalance
+# class.
+# TODO(johnnychen94): maybe support UInt32 and beyond by providing a more efficient implementation
+function _bitrotate_quotation_space(::Type{T}) where T<:SupportedEncodingTypes
+    s = Vector{T}[]
+    for x in typemin(T):typemax(T)
+        # without the following skipping mechanism actually runs faster
+        # any(c->x in c, s) && continue
+        push!(s, sort!(bitrotate.(x, 0:8sizeof(T)-1)))
+    end
+    Dict(minimum(c) => c for c in unique!(s))
+end
+function _build_inverse_table(d::Dict{T,<:AbstractVector{T}}) where {T}
+    id = Vector{T}(undef, typemax(T)-typemin(T)+1)
+    for k in keys(d)
+        for q in d[k]
+            id[q+1] = k
+        end
+    end
+    return id
+end
+_freeze(::Type{T}, v::Vector) where T = SVector{typemax(T)-typemin(T)+1}(v)

test/lbp_ri.jl

@@ -0,0 +1,33 @@
+@testset "lbp_ri" begin
+    # reference result comes from [1]
+    # - [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.
+    X = [6 7 9; 5 6 3; 2 1 7]
+    ref_out = [3 1 0; 13 53 27; 5 91 0]
+
+    out = lbp_ri(X)
+    @test eltype(out) == UInt8
+    @test size(out) == (3, 3)
+    @test out == ref_out
+
+    fill!(out, 0)
+    lbp_ri!(out, X)
+    @test out == ref_out
+
+    # intensity-invariant
+    @test lbp_ri(X./255) == out
+    # Gray inputs
+    @test lbp_ri(Gray.(X./255)) == out
+
+    # LBP is only defined for scalar values
+    @test_throws MethodError lbp_ri(RGB.(Gray.(X./255)))
+
+    # not yet ready for N-dimensional array (although it's doable)
+    @test_throws MethodError lbp_ri(rand(3, 3, 3))
+
+    @testset "OffsetArrays" begin
+        Xo = OffsetArray(X, -1, -1)
+        out = lbp_ri(Xo)
+        @test axes(out) == axes(Xo)
+        @test OffsetArrays.no_offset_view(out) == ref_out
+    end
+end

test/runtests.jl

@@ -13,5 +13,6 @@ using Aqua, Documenter
 end
 
 include("lbp_original.jl")
+include("lbp_ri.jl")
 
 end