From 5687005579e161546f27833ead40351b358f2f4a Mon Sep 17 00:00:00 2001 From: Johnny Chen Date: Mon, 18 Oct 2021 23:27:09 +0800 Subject: [PATCH] implement FiniteDiffs submodule: gradient, divergence and laplacian operators (#22) New operators: - gradient and adjoint gradient: `fgradient` - divergence: `fdiv` - laplacian: `flaplacian` To better organize the symbols, I deprecate two entrypoints: * `ImageBase.fdiff` => `ImageBase.FiniteDiff.fdiff` * `ImageBase.fdiff!` => `ImageBase.FiniteDiff.fdiff!` * disable meta quality test in old Julia versions Maintaining old compatibility becomes quite troublesome especially when we have 1.6 as the new LTS version now. * CI compatibility fix for Julia < 1.3 Co-authored-by: Tim Holy --- .github/workflows/UnitTest.yml | 8 ++ Project.toml | 3 +- src/ImageBase.jl | 7 -- src/compat.jl | 9 -- src/deprecated.jl | 5 ++ src/diff.jl | 145 ++++++++++++++++++++++++++++++++- test/deprecated.jl | 7 ++ test/diff.jl | 84 ++++++++++++++++++- test/runtests.jl | 24 +++--- 9 files changed, 257 insertions(+), 35 deletions(-) delete mode 100644 src/compat.jl diff --git a/.github/workflows/UnitTest.yml b/.github/workflows/UnitTest.yml index 968caad..8ac1685 100644 --- a/.github/workflows/UnitTest.yml +++ b/.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 diff --git a/Project.toml b/Project.toml index 7357fcc..139d81e 100644 --- a/Project.toml +++ b/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"] diff --git a/src/ImageBase.jl b/src/ImageBase.jl index cd4f0cd..8432a61 100644 --- a/src/ImageBase.jl +++ b/src/ImageBase.jl @@ -5,11 +5,6 @@ export # originally from ImageTransformations.jl restrict, - # finite difference on one-dimension - # originally from Images.jl - fdiff, - fdiff!, - # basic image statistics, from Images.jl minimum_finite, maximum_finite, @@ -26,13 +21,11 @@ using Reexport using Base.Cartesian: @nloops @reexport using ImageCore using ImageCore.OffsetArrays -using ImageCore.MappedArrays: of_eltype include("diff.jl") include("restrict.jl") include("utils.jl") include("statistics.jl") -include("compat.jl") include("deprecated.jl") if VERSION >= v"1.4.2" # work around https://github.com/JuliaLang/julia/issues/34121 diff --git a/src/compat.jl b/src/compat.jl deleted file mode 100644 index 481dddc..0000000 --- a/src/compat.jl +++ /dev/null @@ -1,9 +0,0 @@ -if VERSION < v"1.2" - require_one_based_indexing(A...) = !Base.has_offset_axes(A...) || throw(ArgumentError("offset arrays are not supported but got an array with index other than 1")) -else - const require_one_based_indexing = Base.require_one_based_indexing -end - -if VERSION < v"1.1" - isnothing(x) = x === nothing -end diff --git a/src/deprecated.jl b/src/deprecated.jl index ac74e0a..ba7d43c 100644 --- a/src/deprecated.jl +++ b/src/deprecated.jl @@ -8,4 +8,9 @@ @deprecate maxfinite(A; kwargs...) maximum_finite(A; kwargs...) @deprecate maxabsfinite(A; kwargs...) maximum_finite(abs, A; kwargs...) +# These two symbols are exported by previous ImageBase versions and now organized in the +# FiniteDiff submodule. +@deprecate fdiff(args...; kwargs...) ImageBase.FiniteDiff.fdiff(args...; kwargs...) +@deprecate fdiff!(args...; kwargs...) ImageBase.FiniteDiff.fdiff!(args...; kwargs...) + # END 0.1 deprecation diff --git a/src/diff.jl b/src/diff.jl index c2fa425..1195a75 100644 --- a/src/diff.jl +++ b/src/diff.jl @@ -1,5 +1,36 @@ -# TODO: add keyword `shrink` to give a consistant result on Base -# when this is done, then we can propose this change to upstream Base +# Because there exists the `FiniteDiff.jl` package with quite different purposes, +# this module is not expected to be reexported. +module FiniteDiff + +using ImageCore +using ImageCore.MappedArrays: of_eltype + +""" +Although stored as an array, image can also be viewed as a function from discrete grid space +Zᴺ to continuous space R if it is gray image, to C if it is complex-valued image +(MRI rawdata), to Rᴺ if it is colorant image, etc. +This module provides the discrete +version of gradient-related operators by viewing image arrays as functions. + +This module provides: + +- forward/backward difference [`fdiff`](@ref) are the Images-flavor of `Base.diff` +- gradient operator [`fgradient`](@ref) and its adjoint via keyword `adjoint=true`. +- divergence operator [`fdiv`](@ref) computes the sum of discrete derivatives of vector + fields. +- laplacian operator [`flaplacian`](@ref) is the divergence of the gradient fields. + +Every function in this module has its in-place version. +""" +FiniteDiff + +export + fdiff, fdiff!, + fdiv, fdiv!, + fgradient, fgradient!, + flaplacian, flaplacian! + + """ fdiff(A::AbstractArray; dims::Int, rev=false, boundary=:periodic) @@ -19,7 +50,7 @@ Take vector as an example, it computes `(A[2]-A[1], A[3]-A[2], ..., A[1]-A[end]) # Examples -```jldoctest; setup=:(using ImageBase: fdiff) +```jldoctest; setup=:(using ImageBase.FiniteDiff: fdiff) julia> A = [2 4 8; 3 9 27; 4 16 64] 3×3 $(Matrix{Int}): 2 4 8 @@ -106,3 +137,111 @@ _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...) + +Compute the divergence of vector fields `Vs`. + +See also [`fdiv!`](@ref) for the in-place version. +""" +function fdiv(V₁::AbstractArray{T}, Vs::AbstractArray{T}...) where T + fdiv!(similar(V₁, maybe_floattype(T)), V₁, Vs...) +end +fdiv(Vs::Tuple) = fdiv(Vs...) + +""" + fdiv!(out, Vs...) + +The in-place version of divergence operator [`fdiv`](@ref). +""" +function fdiv!(out::AbstractArray, V₁::AbstractArray, Vs::AbstractArray...) + # This is the optimized version of `sum(v->fgradient(v; ajoint=true), (V₁, Vs...))` + # by removing unnecessary memory allocations. + all(v->axes(v) == axes(out), (V₁, Vs...)) || throw(ArgumentError("All axes of vector fields Vs and X should be the same.")) + + # TODO(johnnychen94): for better performance, we can eliminate this `tmp` allocation by fusing multiple `fdiff` in the inner loop. + out .= fdiff(V₁; dims=1, rev=true, boundary=:periodic) + tmp = similar(out) + for i in 1:length(Vs) + out .+= fdiff!(tmp, Vs[i]; dims=i+1, rev=true, boundary=:periodic) + end + return out +end +fdiv!(out::AbstractArray, Vs::Tuple) = fdiv!(out, Vs...) + +""" + flaplacian(X::AbstractArray) + +The discrete laplacian operator, i.e., the divergence of the gradient fields of `X`. + +See also [`flaplacian!`](@ref) for the in-place version. +""" +flaplacian(X::AbstractArray) = flaplacian!(similar(X, maybe_floattype(eltype(X))), X) + +""" + flaplacian!(out, X) + flaplacian!(out, ∇X::Tuple, X) + +The in-place version of the laplacian operator [`flaplacian`](@ref). + +!!! tip Avoiding allocations + The two-argument method will allocate memory to store the intermediate + gradient fields `∇X`. If you call this repeatedly with images of consistent size and type, + consider using the three-argument form with pre-allocated memory for `∇X`, + which will eliminate allocation by this function. +""" +flaplacian!(out, X::AbstractArray) = fdiv!(out, fgradient(X)) +flaplacian!(out, ∇X::Tuple, X::AbstractArray) = fdiv!(out, fgradient!(∇X, X)) + + +""" + fgradient(X::AbstractArray; adjoint=false) -> (∂₁X, ∂₂X, ..., ∂ₙX) + +Computes the gradient fields of `X`. If `adjoint==true` then it computes the adjoint gradient +fields. + +Each gradient vector is computed as forward difference along specific dimension, e.g., +[`∂ᵢX = fdiff(X, dims=i)`](@ref fdiff). + +Mathematically, the adjoint operator ∂ᵢ' of ∂ᵢ is defined as `<∂ᵢu, v> := `. + +See also the in-place version [`fgradient!(X)`](@ref) to reuse the allocated memory. +""" +function fgradient(X::AbstractArray{T,N}; adjoint::Bool=false) where {T,N} + fgradient!(ntuple(i->similar(X, maybe_floattype(T)), N), X; adjoint=adjoint) +end + +""" + fgradient!(∇X::Tuple, X::AbstractArray; adjoint=false) + +The in-place version of (adjoint) gradient operator [`fgradient`](@ref). + +The input `∇X = (∂₁X, ∂₂X, ..., ∂ₙX)` is a tuple of arrays that are similar to `X`, i.e., +`eltype(∂ᵢX) == eltype(X)` and `axes(∂ᵢX) == axes(X)` for all `i`. +""" +function fgradient!(∇X::NTuple{N, <:AbstractArray}, X; adjoint::Bool=false) where N + all(v->axes(v) == axes(X), ∇X) || throw(ArgumentError("All axes of vector fields ∇X and X should be the same.")) + for i in 1:N + if adjoint + # the negative adjoint of gradient operator for forward difference is the backward difference + # see also + # Getreuer, Pascal. "Rudin-Osher-Fatemi total variation denoising using split Bregman." _Image Processing On Line_ 2 (2012): 74-95. + fdiff!(∇X[i], X, dims=i, rev=true) + # TODO(johnnychen94): ideally we can get avoid flipping the signs for better performance. + @. ∇X[i] = -∇X[i] + else + fdiff!(∇X[i], X, dims=i) + end + end + return ∇X +end + + + +if VERSION < v"1.1" + isnothing(x) = x === nothing +end + +end # module diff --git a/test/deprecated.jl b/test/deprecated.jl index 5aea6ed..b43d249 100644 --- a/test/deprecated.jl +++ b/test/deprecated.jl @@ -14,4 +14,11 @@ @test maxabsfinite(A) == maximum_finite(abs, A) @test maxabsfinite(A, dims=1) == maximum_finite(abs, A, dims=1) end + + @testset "fdiff entrypoints" begin + A = rand(Float32, 5) + @test ImageBase.fdiff(A, rev=true) == ImageBase.FiniteDiff.fdiff(A, rev=true) + out = similar(A) + @test ImageBase.fdiff!(out, A, rev=false) == ImageBase.FiniteDiff.fdiff!(out, A, rev=false) + end end diff --git a/test/diff.jl b/test/diff.jl index 89ed9e0..75b6798 100644 --- a/test/diff.jl +++ b/test/diff.jl @@ -1,8 +1,12 @@ +using ImageBase.FiniteDiff +# TODO(johnnychen94): remove this after we delete `Imagebase.fdiff` and `ImageBase.fdiff!` entrypoints +using ImageBase.FiniteDiff: fdiff, fdiff! + @testset "fdiff" begin # Base.diff doesn't promote integer to float - @test ImageBase.maybe_floattype(Int) == Int - @test ImageBase.maybe_floattype(N0f8) == Float32 - @test ImageBase.maybe_floattype(RGB{N0f8}) == RGB{Float32} + @test ImageBase.FiniteDiff.maybe_floattype(Int) == Int + @test ImageBase.FiniteDiff.maybe_floattype(N0f8) == Float32 + @test ImageBase.FiniteDiff.maybe_floattype(RGB{N0f8}) == RGB{Float32} @testset "API" begin # fdiff! works the same as fdiff @@ -107,3 +111,77 @@ @test fdiff(A, dims=1) == fdiff(float.(A), dims=1) end end + +@testset "fgradient" begin + for T in generate_test_types([N0f8, Float32], [Gray, RGB]) + for sz in [(7, ), (7, 5), (7, 5, 3)] + A = rand(T, sz...) + + ∇A = fgradient(A) + for i in 1:length(sz) + @test ∇A[i] == fdiff(A, dims=i) + end + ∇A = map(similar, ∇A) + @test fgradient!(∇A, A) == fgradient(A) + + ∇tA = fgradient(A, adjoint=true) + for i in 1:length(sz) + @test ∇tA[i] == .-fdiff(A, dims=i, rev=true) + end + ∇tA = map(similar, ∇tA) + @test fgradient!(∇tA, A, adjoint=true) == fgradient(A, adjoint=true) + end + end +end + +@testset "fdiv/flaplacian" begin + ref_laplacian(X) = imfilter(X, Kernel.Laplacian(ntuple(x->true, ndims(X))), "circular") + + X = [ + 5 3 8 1 2 2 3 + 5 5 1 3 3 1 1 + 1 8 2 9 1 2 7 + 7 3 4 5 8 1 5 + 1 4 1 8 8 9 7 + ] + ΔX_ref = [ + -8 10 -26 17 6 7 3 + -8 -3 14 2 -5 4 12 + 23 -21 14 -25 18 2 -19 + -18 11 -5 9 -17 20 2 + 19 -8 20 -17 -5 -18 -10 + ] + ΔX = ref_laplacian(X) + # Base.diff doesn't promote Int to floats so we should probably do the same for laplacian + @test eltype(ΔX) == Int + @test ΔX_ref == ΔX + + for T in generate_test_types([N0f8, Float32], [Gray, RGB]) + for sz in [(7,), (7, 7), (7, 7, 7)] + A = rand(T, sz...) + ∇A = fgradient(A) + out = fdiv(∇A) + @test out ≈ ref_laplacian(A) + + fill!(out, zero(T)) + fdiv!(out, ∇A...) + @test out == fdiv(∇A) + end + end + + for T in generate_test_types([N0f8, Float32], [Gray, RGB]) + for sz in [(7,), (7, 7), (7, 7, 7)] + A = rand(T, sz...) + out = flaplacian(A) + @test out ≈ ref_laplacian(A) + + ∇A = fgradient(A) + foreach((out, ∇A...)) do x + fill!(x, zero(T)) + end + flaplacian!(out, ∇A, A) + @test out == flaplacian(A) + @test ∇A == fgradient(A) + end + end +end diff --git a/test/runtests.jl b/test/runtests.jl index 6617826..4269e49 100644 --- a/test/runtests.jl +++ b/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