Move extremization, blob_LoG from Images by timholy · Pull Request #223 · JuliaImages/ImageFiltering.jl

timholy · 2021-09-01T12:13:58Z

I think (?) the best way to do this is to make a minor bump. So we should make sure we move everything at once. Am I missing things? What about imROF?

xref JuliaImages/ImageSegmentation.jl#72

codecov · 2021-09-01T12:14:12Z

Codecov Report

Merging #223 (db0068a) into master (f85e7b0) will decrease coverage by 1.01%.
The diff coverage is 71.76%.

❗ Current head db0068a differs from pull request most recent head 88a9e89. Consider uploading reports for the commit 88a9e89 to get more accurate results

@@            Coverage Diff             @@
##           master     #223      +/-   ##
==========================================
- Coverage   91.31%   90.29%   -1.02%     
==========================================
  Files           9       12       +3     
  Lines        1450     1535      +85     
==========================================
+ Hits         1324     1386      +62     
- Misses        126      149      +23

Impacted Files	Coverage Δ
src/ImageFiltering.jl	`77.27% <ø> (ø)`
src/findall_window.jl	`0.00% <0.00%> (ø)`
src/compat.jl	`50.00% <50.00%> (ø)`
src/extrema.jl	`100.00% <100.00%> (ø)`
src/kernel.jl	`99.17% <0.00%> (+0.82%)`	⬆️

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 f85e7b0...88a9e89. Read the comment docs.

timholy · 2021-09-01T12:31:05Z

Is the doc failure (segfault 😱) any reason to hold off on merging this?

johnnychen94 · 2021-09-01T12:34:46Z

The documentation segfault might be related to Plots, see also JuliaImages/juliaimages.github.io#210. I haven't had the time to investigate it so the best strategy is to temporarily remove it.

src/extrema.jl

johnnychen94 · 2021-09-01T12:47:33Z

src/extrema.jl

+            # If i is in the interior, we don't have to worry about i+j being out-of-bounds
+            for j in Rregion
+                j == z && continue
+                if !Base.Order.lt(order, img[i+j], img_i)


I was actually thinking if it's better to expose this as argument f and make

findlocalextrema(f, img; [dims,] [edges])

And furthermore, maybe it's better to have it work similarly to mapwindow:

findall_window(f, img, window; [border,] [indices])

That's pretty interesting. From a pragmatic standpoint one would probably want to smooth with a separable filter and then use a 3x3 window rather than an nxn window (the latter is O(n^2), the former is O(n)). But I do see the point.

johnnychen94 · 2021-09-01T13:12:31Z

Am I missing things? What about imROF?

I guess gaussian_pyramid (I've commented on JuliaImages/Images.jl#971 (comment))

This makes more dramatic depatures from our old implementations: - add a `show` method for `BlobLog` - encode the entire multidimensional `σ` (after multiplying by `σshape`) - add a threshold to dispose of spurious maxima by default - separate the filtering of `blob_LoG` into an independent (but non-exported) function - change the API of `findlocalmaxima` etc to be `window`-based rather than `dims`-based. This is strictly more flexible. Co-authored-by: Johnny Chen <johnnychen94@hotmail.com>

timholy · 2021-09-01T20:42:25Z

OK, I've made the suggested changes and more. Much nicer IMO, but see what you think. Thanks for suggesting the more thorough revision, it's really the perfect time. I didn't quite mimic the mapwindow interface but it's closer. (I don't think border and indices are relevant, but the window suggestion is brilliant.)

The commit comment 8d854bb explains the changes.

johnnychen94

I'm thinking of an even generic window version of findall that iterates over all the interior points of img and return all p that satisfies f(window_p) == true. Just like mapwindow is the window version of map, I'm thinking it as the window version of findall.

Since we need this PR for JuliaImages/ImageSegmentation.jl#72, maybe we can leave findall_window as future work if that sounds too much coding work.

johnnychen94 · 2021-09-02T05:43:14Z

src/ImageFiltering.jl

-    freqkernel, spacekernel
+    freqkernel, spacekernel,
+    findlocalminima, findlocalmaxima,
+    blob_LoG, BlobLoG


Do we really need to export the struct BlobLoG?

Together with my kwarg method, it makes the output cut/pasteable ("round-trippable"), but maybe that's not important.

johnnychen94 · 2021-09-02T05:44:26Z

src/extrema.jl

+    Rinterior = clippedinds(R0, halfwindow)
+    Rwindow = _colon(-halfwindow, halfwindow)
+    z = zero(halfwindow)
+    for i in R


It looks valid for me to skip the default bounds checking here:

Suggested change

for i in R

@inbounds for i in R

timholy · 2021-09-02T06:36:33Z

I'm thinking of an even generic window version of findall that iterates over all the interior points of img and return all p that satisfies f(window_p) == true. Just like mapwindow is the window version of map, I'm thinking it as the window version of findall...maybe we can leave findall_window as future work if that sounds too much coding work.

An excellent design concept and a worthy goal. No time like the present...

johnnychen94 · 2021-09-02T06:38:13Z

No time like the present...

Yeah... I'm also running short of time..

Merge when you ready.

timholy · 2021-09-02T06:39:08Z

Sorry, I meant "there is no better time than the present." So I'll try it now. (Yes, I am busy too, but your suggestion is too good to leave for later.)

EDIT: https://idioms.thefreedictionary.com/There%27s+no+time+like+the+present

timholy · 2021-09-02T06:47:36Z

One interesting question: when you're within a window's width of the boundary, the center point of the valid box is not the base point. Should f be f(V, base_index) where V is the view of img? That makes it a bit different from findall.

johnnychen94 · 2021-09-02T06:57:32Z

Should f be f(V, base_index) where V is the view of img?

My understanding on findall is that we are checking the values and not the index. If we don't consider border condition, it might be as simple as:

results = CartesianIndex{N}[]
for (window_ax, i) in zip(TileIterator(axes(X), window)), CartesianIndices(X))
    window = @view X[window_ax...]
    f(window) && push!(results, i)
end

Well. Some efforts are needed to achieve the current performance of findlocalextrema, however. I did try this before I commented #223 (review) but it seems 2x slower.

For findlocalminima I defined f(window) = minimal(window) == window[OffsetArrays.center(window)...]

timholy · 2021-09-02T07:04:37Z

Currently we exclude windows where A[j] == A[i] when i != j, and several of our tests fail if we change that. For example, using > for f, clearly the meaning is "the base point has a higher value for all other points" whereas using >= wouldn't care. findall_window(>=, zeros(m, n), window) would return all points in the array!

timholy · 2021-09-02T07:44:18Z

OK, see if you like this API.

Performance:

julia> A = rand(100, 101);

julia> @btime findlocalmaxima($A);
  132.642 μs (11 allocations: 64.47 KiB)

julia> @btime ImageFiltering.findmax_window($A, $(ImageFiltering.default_window(A)));
  182.121 μs (11 allocations: 64.47 KiB)

So not too shabby, but probably enough of a difference that we still want findlocalmaxima.

johnnychen94 · 2021-09-02T07:56:55Z

Will f(img[centerpoint], img[otherpoint]) too strict assumption here? I think f(Rview) == true should be the most generic case. Recall my slow mapwindow example, JuliaImages/ImageSegmentation.jl#72 (comment) I feel it's more natural to write the do-block here.

Of course, the performance is really an issue.

Edit: let me also take a try :) This function is quite interesting.

timholy · 2021-09-02T08:03:50Z

My docs were lying a bit, I combined two things I shouldn't have combined. This is more accurate (the code itself hasn't changed).

timholy · 2021-09-02T08:05:15Z

Did you see this (edited) comment:

findall_window(>=, zeros(m, n), window) would return all points in the array!

That's the core problem with your implementation in JuliaImages/ImageSegmentation.jl#72 (comment).

johnnychen94 · 2021-09-02T08:11:38Z

OK, see if you like this API.

Right, I think I'm thinking of the same implementation as your current version. I'll try more if there are better ideas to optimize the performance and then will comment.

If we allow paddings like mapwindow then we don't need to pass extra i to f(OffsetArray(view(img, Rview), Rview.indices), i); we can simply pass the local coordinates f(Rview). I'm not sure of it, the current version is also good to me.

timholy · 2021-09-02T08:20:59Z

Not sure I understand this:

we can simply pass the local coordinates f(Rview)

That only passes the indices, and leaves the array itself implicit. That's pretty different from how mapwindow works. Is that really what you mean? If so, then you don't even need to pass img to mapwindow and findall_window, and that's different from map and findall.

timholy · 2021-09-02T08:30:35Z

Ooh, I got it down a little further:

julia> @btime ImageFiltering.findmax_window($A, $(ImageFiltering.default_window(A)));
  158.325 μs (6 allocations: 43.30 KiB)

Do you think that's a small enough gap (~15%) that we could do away with findlocalmaxima and findlocalminima? Or is it still better to have those? If so, I think we could defer the decision about how findall_window works.

johnnychen94 · 2021-09-02T08:54:20Z

src/findall_window.jl

+@inline function all_window(f::F, V::AbstractArray, basepoint::CartesianIndex, excludepoint::Union{Nothing,CartesianIndex}=nothing) where F
+    @inbounds ref = V[basepoint]
+    @inbounds for (i, v) in pairs(V)
+        i == excludepoint || f(ref, v) || return false


I'm expecting f(ref, v) is scalar operation and is type stable, but if I change to

i == excludepoint || ref == v || return false

I get

julia> @btime ImageFiltering.findmax_window($A, $(ImageFiltering.default_window(A))); 74.347 μs (1 allocation: 80 bytes)

Surprisingly less allocations, do you have any ideas?

Edit: hmmm, maybe is due to the push! operation in findall_window.

Presumably the output array is empty? So it's probably just push! reallocating the output buffer.

With

function findall_window(f::F, img::AbstractArray{T,N}, window::Dims{N}; allow_edges::NTuple{N,Bool}=ntuple(_->true, N)) where {F,T<:Union{Gray,Number},N} buffer_size = min(512, length(img)) # a heuristic threshold to work around the overhead of `push!` basepoints = Vector{CartesianIndex{N}}(undef, buffer_size) Iedge = CartesianIndex(map(!, allow_edges)) R0 = CartesianIndices(img) R = clippedinds(R0, Iedge) halfwindow = CartesianIndex(map(x -> x >> 1, window)) i_positive = 0 @inbounds for i in R Rview = _colon(i-halfwindow, i+halfwindow) ∩ R0 if f(OffsetArray(view(img, Rview), Rview.indices), i) i_positive += 1 if length(basepoints) < i_positive resize!(basepoints, min(2length(basepoints), length(img))) end basepoints[i_positive] = i end end return resize!(basepoints, i_positive) end

I'm getting closer:

julia> @btime ImageFiltering.findmax_window($A, $(ImageFiltering.default_window(A))); 146.643 μs (4 allocations: 60.31 KiB) # after 167.950 μs (11 allocations: 64.47 KiB) # before

Will this be too tricky?

sizehint! does the same thing and is much simpler.

And these optimizations could be applied to findlocalextrema, widening the gap once again. (The same comment applies to my own i == excludepoint || ... change.)

I think the key question is any overhead from creating the wrapper, which findlocalextrema doesn't do.

timholy · 2021-09-02T11:42:26Z

Should we split out the findall_window into a separate PR? Or settle on it here?

johnnychen94 · 2021-09-02T11:50:54Z

Since the API is settled we can have both here and then try to improve the performance in future PRs.

Edit: I'll travel back to Shanghai tomorrow for the new semester so won't be available to tackle the performance issue this weekend.

johnnychen94

Since the API is settled we can have both here and then try to improve the performance in future PRs.

Well... Sorry for the back and forth. Let's separate it to another PR so that we have more time to think about the implementation.

johnnychen94 · 2021-09-02T13:13:36Z

src/findall_window.jl

@@ -0,0 +1,51 @@
+"""
+    findall_window(f, img::AbstractArray{T,N}, window::Dims{N}; allow_edges::(true...,))


I just realize that allow_edges doesn't work correctly when window != (3, 3).

It depends on what you mean by "correctly." The sense here is, "will it return a solution on the edge of the image?" In that sense it works fine (I think). What sense did you mean?

Move extremization, blob_LoG from Images

bad0c70

johnnychen94 reviewed Sep 1, 2021

View reviewed changes

johnnychen94 mentioned this pull request Sep 1, 2021

🧹 code organization and dead code elimination JuliaImages/Images.jl#971

Merged

1 task

Fix & test IO; test rthresh

88a9e89

johnnychen94 approved these changes Sep 2, 2021

View reviewed changes

johnnychen94 mentioned this pull request Sep 2, 2021

introduce a findall version of mapwindow #224

Open

timholy force-pushed the teh/move_extrema branch from 5bd10ea to d84449d Compare September 2, 2021 08:03

johnnychen94 reviewed Sep 2, 2021

View reviewed changes

timholy force-pushed the teh/move_extrema branch from db0068a to 88a9e89 Compare September 2, 2021 18:09

timholy merged commit 02407fa into master Sep 2, 2021

timholy deleted the teh/move_extrema branch September 2, 2021 20:45

timholy mentioned this pull request Sep 3, 2021

Add findall_window #225

Draft

		@@ -0,0 +1,51 @@
		"""
		findall_window(f, img::AbstractArray{T,N}, window::Dims{N}; allow_edges::(true...,))

Conversation

timholy commented Sep 1, 2021

Uh oh!

codecov bot commented Sep 1, 2021 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Codecov Report

Uh oh!

timholy commented Sep 1, 2021

Uh oh!

johnnychen94 commented Sep 1, 2021

Uh oh!

Uh oh!

Uh oh!

Uh oh!

johnnychen94 Sep 1, 2021 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

johnnychen94 Sep 1, 2021 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

johnnychen94 commented Sep 1, 2021

Uh oh!

timholy commented Sep 1, 2021 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

johnnychen94 left a comment • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

timholy commented Sep 2, 2021

Uh oh!

johnnychen94 commented Sep 2, 2021

Uh oh!

timholy commented Sep 2, 2021 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

timholy commented Sep 2, 2021 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

johnnychen94 commented Sep 2, 2021 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

timholy commented Sep 2, 2021 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

timholy commented Sep 2, 2021

Uh oh!

johnnychen94 commented Sep 2, 2021 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

timholy commented Sep 2, 2021

Uh oh!

timholy commented Sep 2, 2021

Uh oh!

johnnychen94 commented Sep 2, 2021 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

timholy commented Sep 2, 2021

Uh oh!

timholy commented Sep 2, 2021

Uh oh!

johnnychen94 Sep 2, 2021 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

johnnychen94 Sep 2, 2021 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

codecov bot commented Sep 1, 2021 •

edited

Loading

johnnychen94 Sep 1, 2021 •

edited

Loading

johnnychen94 Sep 1, 2021 •

edited

Loading

timholy commented Sep 1, 2021 •

edited

Loading

johnnychen94 left a comment •

edited

Loading

timholy commented Sep 2, 2021 •

edited

Loading

timholy commented Sep 2, 2021 •

edited

Loading

johnnychen94 commented Sep 2, 2021 •

edited

Loading

timholy commented Sep 2, 2021 •

edited

Loading

johnnychen94 commented Sep 2, 2021 •

edited

Loading

johnnychen94 commented Sep 2, 2021 •

edited

Loading

johnnychen94 Sep 2, 2021 •

edited

Loading

johnnychen94 Sep 2, 2021 •

edited

Loading

johnnychen94 commented Sep 2, 2021 •

edited

Loading