WIP: add ajoint_fdiff and div

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

src/diff.jl

@@ -52,13 +52,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), 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),
@@ -101,3 +104,67 @@ end
 
 _fdiff_default_dims(A) = nothing
 _fdiff_default_dims(A::AbstractVector) = 1
+
+"""
+    adjoint_fdiff(A::AbstractArray; dims::Int, kwargs...)
+
+The adjoint operator of finite difference [`fdiff`](@ref).
+
+The keywords are the same as [`fdiff`](@ref). See also
+[`ImageBaseadjoint_fdiff!`](@ref) for the in-place version.
+"""
+adjoint_fdiff(A::AbstractArray; kwargs...) = adjoint_fdiff!(similar(A), A; kwargs...)
+
+"""
+    adjoint_fdiff!(dst::AbstractArray, src::AbstractArray; kwargs...)
+
+The in-place version of [`adjoint_fdiff`](@ref).
+"""
+function adjoint_fdiff!(dst::AbstractArray, src::AbstractArray; rev=false, kwargs...)
+    # Mathematically, the negative adjoint is the backward difference. See section 4 in the
+    # following paper:
+    # Getreuer, Pascal. "Rudin-Osher-Fatemi total variation denoising using split Bregman." Image Processing On Line 2 (2012): 74-95.
+    fdiff!(dst, src; rev=!rev, kwargs...)
+    @. dst = -dst
+end
+
+"""
+    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->fdiff(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, V₁, Vs::AbstractArray...; boundary=:periodic)
+    tmp = fdiff(V₁; dims=1, rev=true, boundary=boundary)
+    dst .= tmp
+    for i in 1:length(Vs)
+        dst .+= fdiff!(tmp, Vs[i]; dims=i+1, rev=true, boundary=boundary)
+    end
+    return dst
+end

test/diff.jl

@@ -71,7 +71,7 @@
                 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)
@@ -97,3 +97,21 @@
         @test A_out.parent == fdiff(parent(A), dims=1, rev=true)
     end
 end
+
+@testset "adjoint_fdiff" begin
+    for T in generate_test_types([N0f8, Float32], [Gray, RGB])
+        A = rand(T, 6)
+        @test adjoint_fdiff(A) == .-fdiff(A, rev=true)
+        @test adjoint_fdiff(A, rev=true) == .-fdiff(A, rev=false)
+        @test eltype(adjoint_fdiff(A)) == T
+
+        A = rand(T, 6, 6)
+        for n in 1:ndims(A)
+            @test adjoint_fdiff(A, dims=n) == .-fdiff(A, rev=true, dims=n)
+            @test adjoint_fdiff(A, dims=n, rev=true) == .-fdiff(A, rev=false, dims=n)
+        end
+    end
+
+    A = rand(4, 4)
+    @test_throws UndefKeywordError adjoint_fdiff(A)
+end