implement FiniteDiffs submodule: gradient, divergence and laplacian operators #22

johnnychen94 · 2021-10-14T11:23:35Z

This is a rework of #19 on top of the fdiff function. The advantages of this version are that it's simpler and that it is friendly to GPU.

To better organize the symbols, I introduce a submodule FiniteDiffs (definitely not the FiniteDiffs.jl). I also deprecate the entrypoints ImageBase.fdiff and ImageBase.fdiff!.

If we don't consider the memory limit, benchmark on the Laplacian operator shows that this version is more performant than ImageFiltering's version on gray images, and less performant on RGB images.

using ImageFiltering
using ImageBase
using ImageBase.FiniteDiffs

ref_laplacian(X) = imfilter(X, Kernel.Laplacian(ntuple(x->true, ndims(X))), "circular")

X = rand(Float32, 256, 256);
flaplacian(X) ≈ ref_laplacian(X) # true
@btime flaplacian($X); # 115.390 μs (10 allocations: 1.25 MiB)
@btime ref_laplacian($X); # 319.879 μs (16 allocations: 521.02 KiB)

X = rand(RGB{Float32}, 256, 256);
flaplacian(X) ≈ ref_laplacian(X) # true
@btime flaplacian($X); # 528.354 μs (10 allocations: 3.75 MiB)
@btime ref_laplacian($X); # 520.816 μs (16 allocations: 1.52 MiB)

X = rand(Float32, 1024, 1024);
flaplacian(X) ≈ ref_laplacian(X) # true
@btime flaplacian($X); # 3.134 ms (10 allocations: 20.00 MiB)
@btime ref_laplacian($X); # 5.430 ms (16 allocations: 8.03 MiB)

X = rand(RGB{Float32}, 1024, 1024);
flaplacian(X) ≈ ref_laplacian(X) # true
@btime flaplacian($X); # 12.943 ms (10 allocations: 60.00 MiB)
@btime ref_laplacian($X); # 9.759 ms (16 allocations: 24.06 MiB)

closes #19
closes #20

src/diff.jl

codecov · 2021-10-14T11:28:52Z

Codecov Report

Merging #22 (412952e) into master (bd3db5b) will increase coverage by 1.47%.
The diff coverage is 96.29%.

@@            Coverage Diff             @@
##           master      #22      +/-   ##
==========================================
+ Coverage   90.59%   92.07%   +1.47%     
==========================================
  Files           6        5       -1     
  Lines         202      227      +25     
==========================================
+ Hits          183      209      +26     
+ Misses         19       18       -1

Impacted Files	Coverage Δ
src/diff.jl	`96.55% <96.29%> (-0.23%)`	⬇️

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update bd3db5b...412952e. Read the comment docs.

johnnychen94 · 2021-10-14T13:01:48Z

Benchmark to the optimized Images.div shows that there is still more room for performance if we rewrite using for loops.

X = rand(256, 256);
∇X = fgradient(X);
p = cat(∇X..., dims=3);

@btime Images.div($p); # 48.573 μs (2 allocations: 512.05 KiB)
@btime fdiv($∇X); # 146.490 μs (6 allocations: 1.50 MiB)

X = rand(1024, 1024);
∇X = fgradient(X);
p = cat(∇X..., dims=3);

@btime Images.div($p); # 1.314 ms (2 allocations: 8.00 MiB)
@btime fdiv($∇X); # 4.385 ms (6 allocations: 24.00 MiB)

I'll work out imROF first and then come back revisit this PR.

johnnychen94 · 2021-10-14T19:56:26Z

I think this and #24 are ready. There is still some room for better performance but I'd like to keep it a future work #25 so that we can move forward to a new version #23 and use it in Images.jl

timholy · 2021-10-15T18:03:05Z

I hope to get to this tomorrow morning, sorry for the delay.

johnnychen94 · 2021-10-15T18:47:04Z

No worries, I start a lecture series to introduce Julia to students in our school every Sunday, I'll task switch to it in the meantime. BTW, your new "Why Julia" introduction https://www.youtube.com/watch?v=x4oi0IKf52w is super insightful 😆

timholy · 2021-10-15T20:42:31Z

Next lecture in the series is Monday and then they start coming quickly. Steal anything from https://github.com/timholy/AdvancedScientificComputing that's useful to you (see the schedule link to see what's ahead).

timholy

It's quite telling that almost all my comments are about documentation and tests. Nice work!

src/diff.jl

timholy · 2021-10-16T14:07:15Z

src/diff.jl

+
+- 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) is the negative sum of the adjoint gradient operator of


I wonder if people will be surprised by "negative sum", given that divergence is typically written without appeal to the adjoint. That said, from a technical standpoint you are correct, with one quibble: it's not really the gradient because that would produce n components from a scalar-valued input. Technically, I guess fdiv is defined as <∇u, ∇v> = -<u, fdiv(∇v)>, right?

I wonder if it would be better to be a bit vague here, possibly linking to something? Perhaps something like "the divergence operator fdiv is a sum of discrete derivatives of vector fields" and specify the more precise meaning below?

This is definitely an oversight, I should make this a comment and implementation details instead of the docstring. Thanks for catching it.

timholy · 2021-10-16T14:14:52Z

src/diff.jl

+!!! tips Non-allocating
+    This function will allocate a new set of memories to store the intermediate
+    gradient fields `∇X`, if you pre-allcoate the memory for `∇X`, then this function
+    will use it and is thus non-allcating.


Suggested change

!!! tips Non-allocating

This function will allocate a new set of memories to store the intermediate

gradient fields `∇X`, if you pre-allcoate the memory for `∇X`, then this function

will use it and is thus non-allcating.

!!! 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.

Note it's !!! tip and not !!! tips: https://juliadocs.github.io/Documenter.jl/stable/showcase/#Tip-admonition

Every now and then I find myself evil to take your limited open-source time on language checks. 😢

timholy · 2021-10-16T14:15:36Z

src/diff.jl

+
+The in-place version of divergence operator [`fdiv`](@ref).
+"""
+function fdiv!(out::AbstractArray, V₁::AbstractArray, Vs::AbstractArray...)


I'm guessing this should be @inline, otherwise I think you'll get a bit of splat penalty from callers like fdiv(::Tuple) and flaplacian!`.

I didn't observe a clear performance difference by trying the flaplacian benchmark, so I choose to keep it unchanged for now.

However, I'm quite curious about what makes you trying to think this way. Are there any references or information outside that I can take a look at?

Not sure about references, other than maybe checking the MethodInstances that get created via MethodAnalysis.jl. If you see abstract instances then it's possible that things will be better with forced-inlining. However, varargs typically work out well when the types are homogenous, and that seems likely to be true in this case, which may explain why you didn't see this.

src/diff.jl

timholy · 2021-10-16T19:07:12Z

test/deprecated.jl

+
+    @testset "fdiff entrypoints" begin
+        A = rand(Float32, 5)
+        @test ImageBase.fdiff(A, rev=true) == ImageBase.FiniteDiff.fdiff(A, rev=true)


You could even just check ImageBase.fdiff === ImageBase.FiniteDiff.tdiff and so on.

Unfortunately, this doesn't hold true if I do

@deprecate fdiff ImageBase.FiniteDiff.fdiff

without inspecting the deprecation details, it seems that @deprecate only passes positional arguments and not keyword arguments.

julia> A = rand(Float32, 5); julia> ImageBase.fdiff(A, rev=true) ERROR: MethodError: no method matching fdiff(::Vector{Float32}; rev=true) Closest candidates are: fdiff(::Any...) at deprecated.jl:45 got unsupported keyword argument "rev"

timholy · 2021-10-16T19:08:38Z

test/diff.jl

@@ -107,3 +111,76 @@
        @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])


Suggested change

for T in generate_test_types([N0f8, Float32], [Gray, RGB])

for T in generate_test_types(Any[N0f8, Float32], Any[Gray, RGB])

(tiny bit easier on compiler and speeds up the tests)

I'll fix this and elsewhere in a separate PR.

Edit:

Actually, this doesn't change anything by comparing the first run time; in both cases, they take about 0.07s to compile the function. I'll keep them unchanged for now.

That's fine. If you're fixing inference triggers in SnoopCompile these may show up without the Anys but that's just "analysis noise" and not anything concerning.

timholy · 2021-10-16T19:09:58Z

test/diff.jl

+        19   -8   20  -17   -5  -18  -10
+    ]
+    ΔX = ref_laplacian(X)
+    @test eltype(ΔX) == Int


Is this important to test here?

Yep, to ensure that we don't promote the Int to float types. This gives a consistent result to Base's diff.

To clarify, my point is that if it doesn't pass it's actually a bug in ImageFiltering, but that isn't something that can be fixed here.

You are right, I just keep it here as a detector for potential bug tracing when the test fails. If it turns out to be annoying we can always remove it.

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!`

Maintaining old compatibility becomes quite troublesome especially when we have 1.6 as the new LTS version now.

timholy

Yay!

Just responding to a couple of your excellent points.

timholy · 2021-10-19T07:44:42Z

src/diff.jl

+
+The in-place version of divergence operator [`fdiv`](@ref).
+"""
+function fdiv!(out::AbstractArray, V₁::AbstractArray, Vs::AbstractArray...)


Not sure about references, other than maybe checking the MethodInstances that get created via MethodAnalysis.jl. If you see abstract instances then it's possible that things will be better with forced-inlining. However, varargs typically work out well when the types are homogenous, and that seems likely to be true in this case, which may explain why you didn't see this.

timholy · 2021-10-19T07:46:10Z

test/diff.jl

@@ -107,3 +111,76 @@
        @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])


That's fine. If you're fixing inference triggers in SnoopCompile these may show up without the Anys but that's just "analysis noise" and not anything concerning.

timholy · 2021-10-19T07:47:21Z

test/diff.jl

+        19   -8   20  -17   -5  -18  -10
+    ]
+    ΔX = ref_laplacian(X)
+    @test eltype(ΔX) == Int


To clarify, my point is that if it doesn't pass it's actually a bug in ImageFiltering, but that isn't something that can be fixed here.

Four legacy methods are deprecated in favor of `ImageBase.fdiff`: - `forwarddiffx` - `forwarddiffy` - `backdiffx` - `backdiffy` See the following two PRs for more information: - JuliaImages/ImageBase.jl#11 - JuliaImages/ImageBase.jl#22 This commit bumps ImageBase compatibility to v0.1.5

johnnychen94 commented Oct 14, 2021

View reviewed changes

src/diff.jl Outdated Show resolved Hide resolved

johnnychen94 force-pushed the jc/fdiv_iter branch 2 times, most recently from baae567 to 2dadaee Compare October 14, 2021 12:33

johnnychen94 force-pushed the jc/fdiv_iter branch from 2dadaee to d2a42b8 Compare October 14, 2021 13:07

johnnychen94 mentioned this pull request Oct 14, 2021

ImageBase v0.1.5 #23

Merged

3 tasks

johnnychen94 force-pushed the jc/fdiv_iter branch from d2a42b8 to a86e651 Compare October 14, 2021 15:45

johnnychen94 mentioned this pull request Oct 14, 2021

implement generic ROF model using Chambolle04 primal-dual method #24

Closed

1 task

johnnychen94 force-pushed the jc/fdiv_iter branch from a86e651 to 2eea7af Compare October 14, 2021 19:16

johnnychen94 mentioned this pull request Oct 14, 2021

tweak fdiff, fdiv performance #25

Open

johnnychen94 requested a review from timholy October 14, 2021 19:57

timholy approved these changes Oct 16, 2021

View reviewed changes

johnnychen94 added 3 commits October 18, 2021 23:26

disable meta quality test in old Julia versions

ff56525

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

0c8ffb7

johnnychen94 force-pushed the jc/fdiv_iter branch from 412952e to 0c8ffb7 Compare October 18, 2021 15:26

johnnychen94 merged commit 5687005 into master Oct 18, 2021

johnnychen94 deleted the jc/fdiv_iter branch October 18, 2021 15:27

timholy reviewed Oct 19, 2021

View reviewed changes

johnnychen94 mentioned this pull request Nov 15, 2021

Use IndirectArrays JuliaImages/DitherPunk.jl#47

Merged

adrhill mentioned this pull request Nov 21, 2021

Separate inner loop of error diffusion JuliaImages/DitherPunk.jl#49

Closed

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

implement FiniteDiffs submodule: gradient, divergence and laplacian operators #22

implement FiniteDiffs submodule: gradient, divergence and laplacian operators #22

johnnychen94 commented Oct 14, 2021 •

edited

codecov bot commented Oct 14, 2021 •

edited

johnnychen94 commented Oct 14, 2021 •

edited

johnnychen94 commented Oct 14, 2021

timholy commented Oct 15, 2021

johnnychen94 commented Oct 15, 2021

timholy commented Oct 15, 2021

timholy left a comment

timholy Oct 16, 2021

johnnychen94 Oct 18, 2021

timholy Oct 16, 2021

johnnychen94 Oct 18, 2021

timholy Oct 16, 2021

johnnychen94 Oct 18, 2021 •

edited

timholy Oct 19, 2021

timholy Oct 16, 2021

johnnychen94 Oct 18, 2021

timholy Oct 16, 2021

johnnychen94 Oct 18, 2021 •

edited

timholy Oct 19, 2021

timholy Oct 16, 2021

johnnychen94 Oct 18, 2021

timholy Oct 19, 2021

johnnychen94 Oct 19, 2021

timholy left a comment

timholy Oct 19, 2021

timholy Oct 19, 2021

timholy Oct 19, 2021

-!!! tips Non-allocating
-    This function will allocate a new set of memories to store the intermediate
-    gradient fields `∇X`, if you pre-allcoate the memory for `∇X`, then this function
-    will use it and is thus non-allcating.
+!!! 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.

	for T in generate_test_types([N0f8, Float32], [Gray, RGB])
	for T in generate_test_types(Any[N0f8, Float32], Any[Gray, RGB])

implement FiniteDiffs submodule: gradient, divergence and laplacian operators #22

implement FiniteDiffs submodule: gradient, divergence and laplacian operators #22

Conversation

johnnychen94 commented Oct 14, 2021 • edited

codecov bot commented Oct 14, 2021 • edited

Codecov Report

johnnychen94 commented Oct 14, 2021 • edited

johnnychen94 commented Oct 14, 2021

timholy commented Oct 15, 2021

johnnychen94 commented Oct 15, 2021

timholy commented Oct 15, 2021

timholy left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

johnnychen94 Oct 18, 2021 • edited

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

johnnychen94 Oct 18, 2021 • edited

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

timholy left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

johnnychen94 commented Oct 14, 2021 •

edited

codecov bot commented Oct 14, 2021 •

edited

johnnychen94 commented Oct 14, 2021 •

edited

johnnychen94 Oct 18, 2021 •

edited

johnnychen94 Oct 18, 2021 •

edited