implement DiffView and fdiv

Author: johnnychen94

Project.toml

@@ -14,6 +14,7 @@ julia = "1"
 [extras]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+ImageFiltering = "6a3955dd-da59-5b1f-98d4-e7296123deb5"
 ImageIO = "82e4d734-157c-48bb-816b-45c225c6df19"
 ImageMagick = "6218d12a-5da1-5696-b52f-db25d2ecc6d1"
 OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
@@ -23,4 +24,4 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 TestImages = "5e47fb64-e119-507b-a336-dd2b206d9990"
 
 [targets]
-test = ["Aqua", "Documenter", "Test", "ImageIO", "ImageMagick", "OffsetArrays", "Statistics", "StackViews", "TestImages"]
+test = ["Aqua", "Documenter", "Test", "ImageFiltering", "ImageIO", "ImageMagick", "OffsetArrays", "Statistics", "StackViews", "TestImages"]

src/ImageBase.jl

@@ -9,6 +9,9 @@ export
     # originally from Images.jl
     fdiff,
     fdiff!,
+    fdiv,
+    fdiv!,
+    DiffView,
 
     # basic image statistics, from Images.jl
     minimum_finite,

src/diff.jl

@@ -1,3 +1,94 @@
+abstract type BoundaryCondition end
+struct Periodic <: BoundaryCondition end
+struct ZeroFill <: BoundaryCondition end
+
+"""
+    DiffView(A::AbstractArray, [rev=Val(false)], [bc::BoundaryCondition=Periodic()]; dims)
+
+Lazy version of finite difference [`fdiff`](@ref).
+
+!!! tip
+    For performance, `rev` should be stable type `Val(false)` or `Val(true)`.
+
+# Arguments
+
+- `rev::Bool`
+  If `rev==Val(true)`, then it computes the backward difference
+  `(A[end]-A[1], A[1]-A[2], ..., A[end-1]-A[end])`.
+- `boundary::BoundaryCondition`
+  By default it computes periodically in the boundary, i.e., `Periodic()`.
+  In some cases, one can fill zero values with `ZeroFill()`.
+"""
+struct DiffView{T,N,AT<:AbstractArray,BC,REV} <: AbstractArray{T,N}
+    data::AT
+    dims::Int
+end
+function DiffView(
+        data::AbstractArray{T,N},
+        bc::BoundaryCondition=Periodic(),
+        rev::Union{Val, Bool}=Val(false);
+        dims=_fdiff_default_dims(data)) where {T,N}
+    isnothing(dims) && throw(UndefKeywordError(:dims))
+    rev = to_static_bool(rev)
+    DiffView{maybe_floattype(T),N,typeof(data),typeof(bc),typeof(rev)}(data, dims)
+end
+function DiffView(
+        data::AbstractArray,
+        rev::Union{Val, Bool},
+        bc::BoundaryCondition = Periodic();
+        kwargs...)
+    DiffView(data, bc, rev; kwargs...)
+end
+
+to_static_bool(x::Union{Val{true},Val{false}}) = x
+function to_static_bool(x::Bool)
+    @warn "Please use `Val($x)` for performance"
+    return Val(x)
+end
+
+Base.size(A::DiffView) = size(A.data)
+Base.axes(A::DiffView) = axes(A.data)
+Base.IndexStyle(::DiffView) = IndexCartesian()
+
+Base.@propagate_inbounds function Base.getindex(A::DiffView{T,N,AT,Periodic,Val{true}}, I::Vararg{Int, N}) where {T,N,AT}
+    data = A.data
+    I_prev = map(ntuple(identity, N), I, axes(data)) do i, p, r
+        i == A.dims || return p
+        p == first(r) && return last(r)
+        p - 1
+    end
+    return convert(T, data[I...]) - convert(T, data[I_prev...])
+end
+Base.@propagate_inbounds function Base.getindex(A::DiffView{T,N,AT,Periodic,Val{false}}, I::Vararg{Int, N}) where {T,N,AT}
+    data = A.data
+    I_next = map(ntuple(identity, N), I, axes(data)) do i, p, r
+        i == A.dims || return p
+        p == last(r) && return first(r)
+        p + 1
+    end
+    return convert(T, data[I_next...]) - convert(T, data[I...])
+end
+Base.@propagate_inbounds function Base.getindex(A::DiffView{T,N,AT,ZeroFill,Val{false}}, I::Vararg{Int, N}) where {T,N,AT}
+    data = A.data
+    I_next = I .+ ntuple(i->i==A.dims, N)
+    if checkbounds(Bool, data, I_next...)
+        vi = convert(T, data[I...]) # it requires the caller to pass @inbounds
+        @inbounds convert(T, data[I_next...]) - vi
+    else
+        zero(T)
+    end
+end
+Base.@propagate_inbounds function Base.getindex(A::DiffView{T,N,AT,ZeroFill,Val{true}}, I::Vararg{Int, N}) where {T,N,AT}
+    data = A.data
+    I_prev = I .- ntuple(i->i==A.dims, N)
+    if checkbounds(Bool, data, I_prev...)
+        vi = convert(T, data[I...]) # it requires the caller to pass @inbounds
+        @inbounds vi - convert(T, data[I_prev...])
+    else
+        zero(T)
+    end
+end
+
 # TODO: add keyword `shrink` to give a consistant result on Base
 #       when this is done, then we can propose this change to upstream Base
 """
@@ -52,13 +143,16 @@ julia> fdiff(A, dims=2, boundary=:zero) # fill boundary with zeros
 ```
 
 See also [`fdiff!`](@ref) for the in-place version.
+
+See also [`adjoint_fdiff`](@ref), [`fdiv`](@ref), and [`flaplacian`](@ref) for other related
+discrete functions.
 """
 fdiff(A::AbstractArray; kwargs...) = fdiff!(similar(A, maybe_floattype(eltype(A))), A; kwargs...)
 
 """
     fdiff!(dst::AbstractArray, src::AbstractArray; dims::Int, rev=false, boundary=:periodic)
 
-The in-place version of [`ImageBase.fdiff`](@ref)
+The in-place version of [`fdiff`](@ref).
 """
 function fdiff!(dst::AbstractArray, src::AbstractArray;
         dims=_fdiff_default_dims(src),
@@ -106,3 +200,48 @@ _fdiff_default_dims(A::AbstractVector) = 1
 maybe_floattype(::Type{T}) where T = T
 maybe_floattype(::Type{T}) where T<:FixedPoint = floattype(T)
 maybe_floattype(::Type{CT}) where CT<:Color = base_color_type(CT){maybe_floattype(eltype(CT))}
+
+
+"""
+    fdiv(Vs::AbstractArray...; boundary=:periodic)
+
+Discrete divergence operator for vector field (V₁, V₂, ..., Vₙ).
+
+# Example
+
+Laplacian operator of array `A` is the divergence of its gradient vector field (∂₁A, ∂₂A, ..., ∂ₙA):
+
+```jldoctest
+julia> using ImageFiltering, ImageBase
+
+julia> X = Float32.(rand(1:9, 7, 7));
+
+julia> laplacian(X) = fdiv(ntuple(i->DiffView(X, dims=i), ndims(X))...)
+laplacian (generic function with 1 method)
+
+julia> laplacian(X) == imfilter(X, Kernel.Laplacian(), "circular")
+true
+```
+
+See also [`fdiv!`](@ref) for the in-place version.
+"""
+function fdiv(V₁::AbstractArray, Vs...; kwargs...)
+    fdiv!(similar(V₁, floattype(eltype(V₁))), V₁, Vs...; kwargs...)
+end
+
+"""
+    fdiv!(dst::AbstractArray, Vs::AbstractArray...)
+
+The in-place version of [`fdiv`](@ref).
+"""
+function fdiv!(dst::AbstractArray, Vs::AbstractArray...)
+    ∇ = map(ntuple(identity, length(Vs)), Vs) do n, V
+        DiffView(V, Val(true), dims=n)
+    end
+    @inbounds for i in CartesianIndices(dst)
+        dst[i] = sum(x->_inbound_getindex(x, i), ∇)
+    end
+    return dst
+end
+
+@inline _inbound_getindex(x, i) = @inbounds x[i]

test/diff.jl

@@ -107,3 +107,48 @@
         @test fdiff(A, dims=1) == fdiff(float.(A), dims=1)
     end
 end
+
+@testset "DiffView" begin
+    A = rand(6)
+    @test DiffView(A) ===
+        DiffView(A, Val(false), ImageBase.Periodic()) ===
+        DiffView(A, ImageBase.Periodic(), Val(false)) ===
+        @test_logs((:warn, "Please use `Val(false)` for performance"), DiffView(A, false))
+
+    for T in generate_test_types([N0f8, Float32], [Gray, RGB])
+        A = rand(T, 6)
+        Av = DiffView(A)
+        @test Av == DiffView(A, ImageBase.Periodic(), Val(false))
+        @test eltype(Av) == floattype(T)
+        @test axes(Av) == axes(A)
+        @test Av == fdiff(A)
+        @test DiffView(A, Val(true)) == fdiff(A; rev=true)
+        @test DiffView(A, ImageBase.ZeroFill()) == fdiff(A; boundary=:zero)
+        @test DiffView(A, ImageBase.ZeroFill(), Val(true)) == fdiff(A; boundary=:zero, rev=true)
+
+        A = rand(T, 6, 6)
+        Av = DiffView(A, dims=1)
+        @test eltype(Av) == floattype(T)
+        @test axes(Av) == axes(A)
+        @test Av == fdiff(A, dims=1)
+        @test DiffView(A, Val(true), dims=1) == fdiff(A; dims=1, rev=true)
+        @test DiffView(A, ImageBase.ZeroFill(), dims=1) == fdiff(A; boundary=:zero, dims=1)
+        @test DiffView(A, ImageBase.ZeroFill(), Val(true), dims=1) == fdiff(A; boundary=:zero, rev=true, dims=1)
+    end
+
+    A = OffsetArray(rand(6, 6), -1, -1)
+    Av = DiffView(A, dims=1)
+    @test axes(Av) == axes(A)
+    @test Av == fdiff(A, dims=1)
+end
+
+@testset "fdiv" begin
+    laplacian(X) = fdiv(ntuple(i->DiffView(X, dims=i), ndims(X))...)
+    ref_laplacian(X) = imfilter(X, Kernel.Laplacian(ntuple(x->true, ndims(X))), "circular")
+    for T in generate_test_types([N0f8, Float32], [Gray, RGB])
+        for sz in [(7,), (7, 7), (7, 7, 7)]
+            A = rand(T, sz...)
+            @test laplacian(A) ≈ ref_laplacian(A)
+        end
+    end
+end

test/runtests.jl

@@ -1,4 +1,5 @@
 using ImageBase, OffsetArrays, StackViews
+using ImageFiltering
 using Test, TestImages, Aqua, Documenter
 
 using OffsetArrays: IdentityUnitRange