Add finite difference: fdiff and fdiff! (#11)

Co-authored-by: JKay <[email protected]>

Author: johnnychen94

Project.toml

@@ -12,10 +12,14 @@ Reexport = "0.2, 1"
 julia = "1"
 
 [extras]
+Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
+Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+ImageIO = "82e4d734-157c-48bb-816b-45c225c6df19"
 ImageMagick = "6218d12a-5da1-5696-b52f-db25d2ecc6d1"
 OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
+StackViews = "cae243ae-269e-4f55-b966-ac2d0dc13c15"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 TestImages = "5e47fb64-e119-507b-a336-dd2b206d9990"
 
 [targets]
-test = ["Test", "ImageMagick", "OffsetArrays", "TestImages"]
+test = ["Aqua", "Documenter", "Test", "ImageIO", "ImageMagick", "OffsetArrays", "StackViews", "TestImages"]

README.md

@@ -14,4 +14,10 @@ This package can be seen as an experimental package inside JuliaImages:
 2. is very conservative to external dependencies outside JuliaImages unless there is a real need, in which case,
    it may just fit the first case.
 
+Functions provided by this package:
+
+- `restrict` for two-fold image downsample. (Originally from ImageTransformations.jl)
+- finite difference operator `fdiff`/`fdiff!` (Originally from Images.jl)
+
+
 [ImageCore]: https://github.com/JuliaImages/ImageCore.jl

src/ImageBase.jl

@@ -1,13 +1,23 @@
 module ImageBase
 
-export restrict
+export
+    # two-fold downsampling
+    # originally from ImageTransformations.jl
+    restrict,
+
+    # finite difference on one-dimension
+    # originally from Images.jl
+    fdiff,
+    fdiff!
+
 
 using Reexport
 
 using Base.Cartesian: @nloops
 @reexport using ImageCore
 using ImageCore.OffsetArrays
 
+include("diff.jl")
 include("restrict.jl")
 include("compat.jl")
 include("deprecated.jl")

src/compat.jl

@@ -3,3 +3,7 @@ if VERSION < v"1.2"
 else
     const require_one_based_indexing = Base.require_one_based_indexing
 end
+
+if VERSION < v"1.1"
+    isnothing(x) = x === nothing
+end

src/diff.jl

@@ -0,0 +1,103 @@
+# TODO: add keyword `shrink` to give a consistant result on Base
+#       when this is done, then we can propose this change to upstream Base
+"""
+    fdiff(A::AbstractArray; dims::Int, rev=false, boundary=:periodic)
+
+A one-dimension finite difference operator on array `A`. Unlike `Base.diff`, this function doesn't
+shrink the array size.
+
+Take vector as an example, it computes `(A[2]-A[1], A[3]-A[2], ..., A[1]-A[end])`.
+
+# Keywords
+
+- `rev::Bool`
+  If `rev==true`, then it computes the backward difference
+  `(A[end]-A[1], A[1]-A[2], ..., A[end-1]-A[end])`.
+- `boundary`
+  By default it computes periodically in the boundary, i.e., `:periodic`.
+  In some cases, one can fill zero values with `boundary=:zero`.
+
+# Examples
+
+```jldoctest; setup=:(using ImageBase: fdiff)
+julia> A = [2 4 8; 3 9 27; 4 16 64]
+3×3 $(Matrix{Int}):
+ 2   4   8
+ 3   9  27
+ 4  16  64
+
+julia> diff(A, dims=2) # this function exists in Base
+3×2 $(Matrix{Int}):
+  2   4
+  6  18
+ 12  48
+
+julia> fdiff(A, dims=2)
+3×3 $(Matrix{Int}):
+  2   4   -6
+  6  18  -24
+ 12  48  -60
+
+julia> fdiff(A, dims=2, rev=true) # reverse diff
+3×3 $(Matrix{Int}):
+  -6   2   4
+ -24   6  18
+ -60  12  48
+
+julia> fdiff(A, dims=2, boundary=:zero) # fill boundary with zeros
+3×3 $(Matrix{Int}):
+  2   4  0
+  6  18  0
+ 12  48  0
+```
+
+See also [`fdiff!`](@ref) for the in-place version.
+"""
+fdiff(A::AbstractArray; kwargs...) = fdiff!(similar(A), A; kwargs...)
+
+"""
+    fdiff!(dst::AbstractArray, src::AbstractArray; dims::Int, rev=false, boundary=:periodic)
+
+The in-place version of [`ImageBase.fdiff`](@ref)
+"""
+function fdiff!(dst::AbstractArray, src::AbstractArray;
+        dims=_fdiff_default_dims(src),
+        rev=false,
+        boundary::Symbol=:periodic)
+    isnothing(dims) && throw(UndefKeywordError(:dims))
+    axes(dst) == axes(src) || throw(ArgumentError("axes of all input arrays should be equal. Instead they are $(axes(dst)) and $(axes(src))."))
+    N = ndims(src)
+    1 <= dims <= N || throw(ArgumentError("dimension $dims out of range (1:$N)"))
+
+    r = axes(src)
+    r0 = ntuple(i -> i == dims ? UnitRange(first(r[i]), last(r[i]) - 1) : UnitRange(r[i]), N)
+    r1 = ntuple(i -> i == dims ? UnitRange(first(r[i])+1, last(r[i])) : UnitRange(r[i]), N)
+
+    d0 = ntuple(i -> i == dims ? UnitRange(last(r[i]), last(r[i])) : UnitRange(r[i]), N)
+    d1 = ntuple(i -> i == dims ? UnitRange(first(r[i]), first(r[i])) : UnitRange(r[i]), N)
+
+    if rev
+        dst[r1...] .= view(src, r1...) .- view(src, r0...)
+        if boundary == :periodic
+            dst[d1...] .= view(src, d1...) .- view(src, d0...)
+        elseif boundary == :zero
+            dst[d1...] .= zero(eltype(dst))
+        else
+            throw(ArgumentError("Wrong boundary condition $boundary"))
+        end
+    else
+        dst[r0...] .= view(src, r1...) .- view(src, r0...)
+        if boundary == :periodic
+            dst[d0...] .= view(src, d1...) .- view(src, d0...)
+        elseif boundary == :zero
+            dst[d0...] .= zero(eltype(dst))
+        else
+            throw(ArgumentError("Wrong boundary condition $boundary"))
+        end
+    end
+
+    return dst
+end
+
+_fdiff_default_dims(A) = nothing
+_fdiff_default_dims(A::AbstractVector) = 1

test/diff.jl

@@ -0,0 +1,99 @@
+@testset "fdiff" begin
+    @testset "API" begin
+        # fdiff! works the same as fdiff
+        mat_in = rand(3, 3, 3)
+        mat_out = similar(mat_in)
+        fdiff!(mat_out, mat_in, dims = 2)
+        @test mat_out == fdiff(mat_in, dims = 2)
+
+        mat_in = rand(3, 3, 3)
+        mat_out = similar(mat_in)
+        fdiff!(mat_out, mat_in, dims = 3, rev=true)
+        @test mat_out == fdiff(mat_in, dims = 3, rev=true)
+
+        # default `dims` are only available when input is a vector
+        v = rand(9)
+        @test fdiff(v) == fdiff(v; dims=1)
+        m = rand(3, 3)
+        @test_throws UndefKeywordError fdiff(m)
+    end
+
+    @testset "NumericalTests" begin
+        a = reshape(collect(1:9), 3, 3)
+        b_fd_1 = [1 1 1; 1 1 1; -2 -2 -2]
+        b_fd_2 = [3 3 -6; 3 3 -6; 3 3 -6]
+        b_bd_1 = [-2 -2 -2; 1 1 1; 1 1 1]
+        b_bd_2 = [-6 3 3; -6 3 3; -6 3 3]
+        out = similar(a)
+
+        @test fdiff(a, dims = 1) == b_fd_1
+        @test fdiff(a, dims = 2) == b_fd_2
+        @test fdiff(a, dims = 1, rev=true) == b_bd_1
+        @test fdiff(a, dims = 2, rev=true) == b_bd_2
+        fdiff!(out, a, dims = 1)
+        @test out == b_fd_1
+        fdiff!(out, a, dims = 2)
+        @test out == b_fd_2
+        fdiff!(out, a, dims = 1, rev=true)
+        @test out == b_bd_1
+        fdiff!(out, a, dims = 2, rev=true)
+        @test out == b_bd_2
+
+        @testset "Boundary" begin
+            # These originally lives in Images and we use it to check the compatibility
+            forwarddiffy(u::AbstractMatrix) = [u[2:end,:]; u[end:end,:]] - u
+            forwarddiffx(u::AbstractMatrix) = [u[:,2:end] u[:,end:end]] - u
+            backdiffy(u::AbstractMatrix) = u - [u[1:1,:]; u[1:end-1,:]]
+            backdiffx(u::AbstractMatrix) = u - [u[:,1:1] u[:,1:end-1]]
+
+            X = rand(3, 3)
+            @test forwarddiffx(X) == fdiff(X, dims=2, boundary=:zero)
+            @test forwarddiffy(X) == fdiff(X, dims=1, boundary=:zero)
+            @test backdiffx(X) == fdiff(X, dims=2, rev=true, boundary=:zero)
+            @test backdiffy(X) == fdiff(X, dims=1, rev=true, boundary=:zero)
+        end
+
+        if VERSION >= v"1.1"
+            # dims keyword for diff only exists after Julia 1.1
+            # check numerical results with base implementation
+            drop_last_slice(X, dims) = collect(StackView(collect(eachslice(X, dims=dims))[1:end-1]..., dims=dims))
+            function drop_last_slice(X::AbstractVector, dims)
+                @assert dims==1
+                X[1:end-1]
+            end
+            drop_first_slice(X, dims) = collect(StackView(collect(eachslice(X, dims=dims))[2:end]..., dims=dims))
+            function drop_first_slice(X::AbstractVector, dims)
+                @assert dims==1
+                X[2:end]
+            end
+
+            for N in 1:3
+                sz = ntuple(_->5, N)
+                A = rand(sz...)
+                A_out = similar(A)
+                
+                for dims = 1:N
+                    out_base = diff(A; dims=dims)
+                    out = fdiff(A; dims=dims)
+                    @test out_base == drop_last_slice(out, dims)
+
+                    out_base = .-reverse(diff(reverse(A; dims=dims); dims=dims); dims=dims)
+                    out = fdiff(A; dims=dims, rev=true)
+                    @test out_base == drop_first_slice(out, dims)
+                end
+            end
+        end
+    end
+
+    @testset "OffsetArrays" begin
+        A = OffsetArray(rand(3, 3), -1, -1)
+        A_out = fdiff(A, dims=1)
+        @test axes(A_out) == (0:2, 0:2)
+        @test A_out.parent == fdiff(parent(A), dims=1)
+
+        A = OffsetArray(rand(3, 3), -1, -1)
+        A_out = fdiff(A, dims=1, rev=true)
+        @test axes(A_out) == (0:2, 0:2)
+        @test A_out.parent == fdiff(parent(A), dims=1, rev=true)
+    end
+end

test/runtests.jl

@@ -1,9 +1,30 @@
-using ImageBase, ImageCore, OffsetArrays
-using Test, TestImages
+using ImageBase, OffsetArrays, StackViews
+using Test, TestImages, Aqua, Documenter
 
 using OffsetArrays: IdentityUnitRange
 
 @testset "ImageBase.jl" begin
+
+    @testset "Project meta quality checks" begin
+        # Not checking compat section for test-only dependencies
+        ambiguity_exclude_list = [
+            # https://github.com/JuliaDiff/ChainRulesCore.jl/pull/367#issuecomment-869071000
+            Base.:(==),
+        ]
+        Aqua.test_ambiguities([ImageCore, Base, Core], exclude=ambiguity_exclude_list)
+        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)
+        end
+    end
+    
+    include("diff.jl")
     include("restrict.jl")
 
     @info "deprecations are expected"