WIP: add lazy finite difference DiffView and fdiv

.github/workflows/UnitTest.yml

@@ -48,6 +48,14 @@ jobs:
             ${{ runner.os }}-test-
             ${{ runner.os }}-
       - uses: julia-actions/julia-buildpkg@v1
+      - name: "Compat fix for Julia < v1.3.0"
+        if: ${{ matrix.version == '1.0' }}
+        run: |
+          using Pkg
+          Pkg.add([
+            PackageSpec(name="AbstractFFTs", version="0.5"),
+          ])
+        shell: julia --project=. --startup=no --color=yes {0}
       - uses: julia-actions/julia-runtest@v1
       - uses: julia-actions/julia-processcoverage@v1
       - uses: codecov/codecov-action@v1

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,11 @@ export
     # originally from Images.jl
     fdiff,
     fdiff!,
+    fdiv,
+    fdiv!,
+    flaplacian,
+    flaplacian!,
+    DiffView,
 
     # basic image statistics, from Images.jl
     minimum_finite,

src/diff.jl

@@ -1,3 +1,90 @@
+abstract type BoundaryCondition end
+struct Periodic <: BoundaryCondition end
+struct ZeroFill <: BoundaryCondition end
+
+"""
+    DiffView(A::AbstractArray, dims::Val{D}, [bc::BoundaryCondition=Periodic()], [rev=Val(false)])
+
+Lazy version of finite difference [`fdiff`](@ref).
+
+!!! tip
+    For performance, both `dims` and `rev` require `Val` types.
+
+# Arguments
+
+- `dims::Val{D}`
+  Specify the dimension D that dinite difference is applied to.
+- `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,D,BC,REV,AT<:AbstractArray} <: AbstractArray{T,N}
+    data::AT
+end
+function DiffView(
+        data::AbstractArray{T,N},
+        ::Val{D},
+        bc::BoundaryCondition=Periodic(),
+        rev::Val = Val(false)
+    ) where {T,N,D}
+    DiffView{maybe_floattype(T),N,D,typeof(bc),typeof(rev),typeof(data)}(data)
+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,D,Periodic,Val{true}}, I::Vararg{Int, N}) where {T,N,D}
+    data = A.data
+    r = axes(data, D)
+    x = I[D]
+    x_prev = first(r) == x ? last(r) : x - 1
+    I_prev = update_tuple(I, x_prev, Val(D))
+    return convert(T, data[I...]) - convert(T, data[I_prev...])
+end
+Base.@propagate_inbounds function Base.getindex(A::DiffView{T,N,D,Periodic,Val{false}}, I::Vararg{Int, N}) where {T,N,D}
+    data = A.data
+    r = axes(data, D)
+    x = I[D]
+    x_next = last(r) == x ? first(r) : x + 1
+    I_next = update_tuple(I, x_next, Val(D))
+    return convert(T, data[I_next...]) - convert(T, data[I...])
+end
+Base.@propagate_inbounds function Base.getindex(A::DiffView{T,N,D,ZeroFill,Val{false}}, I::Vararg{Int, N}) where {T,N,D}
+    data = A.data
+    x = I[D]
+    if last(axes(data, D)) == x
+        zero(T)
+    else
+        I_next = update_tuple(I, x+1, Val(D))
+        convert(T, data[I_next...]) - convert(T, data[I...])
+    end
+end
+Base.@propagate_inbounds function Base.getindex(A::DiffView{T,N,D,ZeroFill,Val{true}}, I::Vararg{Int, N}) where {T,N,D}
+    data = A.data
+    x = I[D]
+    if first(axes(data, D)) == x
+        zero(T)
+    else
+        I_prev = update_tuple(I, x-1, Val(D))
+        convert(T, data[I...]) - convert(T, data[I_prev...])
+    end
+end
+
+@generated function update_tuple(A::NTuple{N, T}, x::T, ::Val{i}) where {T, N, i}
+    # This is equivalent to `ntuple(j->j==i ? x : A[j], N)` but is optimized by moving
+    # the if branches to compilation time.
+    ex = :()
+    for j in Base.OneTo(N)
+        new_x = i == j ? :(x) : :(A[$j])
+        ex = :($ex..., $new_x)
+    end
+    return ex
+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
 """
@@ -51,14 +138,15 @@ julia> fdiff(A, dims=2, boundary=:zero) # fill boundary with zeros
  12  48  0
 ```
 
-See also [`fdiff!`](@ref) for the in-place version.
+See also [`fdiff!`](@ref) for the in-place version, and [`DiffView`](@ref) for the
+non-allocating version.
 """
 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 +194,69 @@ _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...)
+
+Discrete divergence operator for vector field (V₁, V₂, ..., Vₙ).
+
+See also [`fdiv!`](@ref) for the in-place version.
+"""
+fdiv(V₁::AbstractArray, Vs...) = fdiv!(similar(V₁, floattype(eltype(V₁))), V₁, Vs...)
+
+"""
+    fdiv!(dst::AbstractArray, Vs::AbstractArray...)
+
+The in-place version of [`fdiv`](@ref).
+"""
+function fdiv!(dst::AbstractArray, Vs::AbstractArray...)
+    # negative adjoint of gradient is equivalent to the reversed finite difference
+    ∇ = fnegative_adjoint_gradient(Vs...)
+    @inbounds for i in CartesianIndices(dst)
+        dst[i] = heterogeneous_getindex_sum(i, ∇...)
+    end
+    return dst
+end
+
+@generated function heterogeneous_getindex_sum(i, Vs::Vararg{<:AbstractArray, N}) where N
+    # This method is equivalent to `sum(V->V[i], Vs)` but is optimized for heterogeneous arrays
+    ex = :(zero(eltype(Vs[1])))
+    for j in Base.OneTo(N)
+        ex = :($ex + Vs[$j][i])
+    end
+    return ex
+end
+
+"""
+    flaplacian(X::AbstractArray)
+
+The Laplacian operator ∇² is the divergence of the gradient operator.
+"""
+flaplacian(X::AbstractArray) = flaplacian!(similar(X, maybe_floattype(eltype(X))), X)
+
+"""
+    flaplacian!(dst::AbstractArray, X::AbstractArray)
+
+The in-place version of the Laplacian operator [`laplacian`](@ref).
+"""
+flaplacian!(dst::AbstractArray, X::AbstractArray) = fdiv!(dst, fgradient(X)...)
+
+# These two functions pass dimension information `Val(i)` to DiffView so that
+# we can move computations to compilation time.
+@generated function fgradient(X::AbstractArray{T, N}) where {T, N}
+    ex = :()
+    for i in Base.OneTo(N)
+        new_x = :(DiffView(X, Val($i), Periodic(), Val(false)))
+        ex = :($ex..., $new_x)
+    end
+    return ex
+end
+@generated function fnegative_adjoint_gradient(Vs::Vararg{<:AbstractArray, N}) where N
+    ex = :()
+    for i in Base.OneTo(N)
+        new_x = :(DiffView(Vs[$i], Val($i), Periodic(), Val(true)))
+        ex = :($ex..., $new_x)
+    end
+    return ex
+end

test/diff.jl

@@ -107,3 +107,41 @@
         @test fdiff(A, dims=1) == fdiff(float.(A), dims=1)
     end
 end
+
+@testset "DiffView" begin
+    for T in generate_test_types([N0f8, Float32], [Gray, RGB])
+        A = rand(T, 6)
+        Av = DiffView(A, Val(1))
+        @test Av == DiffView(A, Val(1), ImageBase.Periodic(), Val(false))
+        @test eltype(Av) == floattype(T)
+        @test axes(Av) == axes(A)
+        @test Av == fdiff(A)
+        @test DiffView(A, Val(1), ImageBase.Periodic(), Val(true)) == fdiff(A; rev=true)
+        @test DiffView(A, Val(1), ImageBase.ZeroFill()) == fdiff(A; boundary=:zero)
+        @test DiffView(A, Val(1), ImageBase.ZeroFill(), Val(true)) == fdiff(A; boundary=:zero, rev=true)
+
+        A = rand(T, 6, 6)
+        Av = DiffView(A, Val(1))
+        @test eltype(Av) == floattype(T)
+        @test axes(Av) == axes(A)
+        @test Av == fdiff(A, dims=1)
+        @test DiffView(A, Val(1), ImageBase.Periodic(), Val(true)) == fdiff(A; dims=1, rev=true)
+        @test DiffView(A, Val(1), ImageBase.ZeroFill()) == fdiff(A; boundary=:zero, dims=1)
+        @test DiffView(A, Val(1), ImageBase.ZeroFill(), Val(true)) == fdiff(A; boundary=:zero, rev=true, dims=1)
+    end
+
+    A = OffsetArray(rand(6, 6), -1, -1)
+    Av = DiffView(A, Val(1))
+    @test axes(Av) == axes(A)
+    @test Av == fdiff(A, dims=1)
+end
+
+@testset "fdiv/flaplacian" begin
+    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 flaplacian(A) ≈ ref_laplacian(A)
+        end
+    end
+end

test/runtests.jl

@@ -1,23 +1,23 @@
 using ImageBase, OffsetArrays, StackViews
+using ImageFiltering
 using Test, TestImages, Aqua, Documenter
 
 using OffsetArrays: IdentityUnitRange
 include("testutils.jl")
 
 @testset "ImageBase.jl" begin
-
     @testset "Project meta quality checks" begin
-        # Not checking compat section for test-only dependencies
-        Aqua.test_ambiguities(ImageBase)
-        Aqua.test_all(ImageBase;
-            ambiguities=false,
-            project_extras=true,
-            deps_compat=true,
-            stale_deps=true,
-            project_toml_formatting=true
-        )
-        if VERSION >= v"1.2"
-            doctest(ImageBase,manual = false)
+        if VERSION >= v"1.3"
+            # Not checking compat section for test-only dependencies
+            Aqua.test_ambiguities(ImageBase)
+            Aqua.test_all(ImageBase;
+                ambiguities=false,
+                project_extras=true,
+                deps_compat=true,
+                stale_deps=true,
+                project_toml_formatting=true
+            )
+            doctest(ImageBase, manual = false)
         end
     end