JuliaImages / ImageFiltering.jl

Julia implementations of multidimensional array convolution and nonlinear stencil operations
Other
99 stars 51 forks source link

Move extremization, blob_LoG from Images #223

Closed timholy closed 2 years ago

timholy commented 2 years ago

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 https://github.com/JuliaImages/ImageSegmentation.jl/pull/72

codecov[bot] commented 2 years ago

Codecov Report

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

:exclamation: Current head db0068a differs from pull request most recent head 88a9e89. Consider uploading reports for the commit 88a9e89 to get more accurate results Impacted file tree graph

@@            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%) :arrow_up:

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 commented 2 years ago

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

johnnychen94 commented 2 years ago

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

johnnychen94 commented 2 years ago

Am I missing things? What about imROF?

I guess gaussian_pyramid (I've commented on https://github.com/JuliaImages/Images.jl/pull/971#issuecomment-910271594)

timholy commented 2 years ago

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.

timholy commented 2 years ago

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 commented 2 years ago

No time like the present...

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

Merge when you ready.

timholy commented 2 years ago

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 commented 2 years ago

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 commented 2 years ago

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 https://github.com/JuliaImages/ImageFiltering.jl/pull/223#pullrequestreview-744624172 but it seems 2x slower.


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

timholy commented 2 years ago

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 commented 2 years ago

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 commented 2 years ago

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, https://github.com/JuliaImages/ImageSegmentation.jl/pull/72#issuecomment-908251756 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 commented 2 years ago

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 commented 2 years ago

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 https://github.com/JuliaImages/ImageSegmentation.jl/pull/72#issuecomment-908251756.

johnnychen94 commented 2 years ago

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 commented 2 years ago

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 commented 2 years ago

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.

timholy commented 2 years ago

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

johnnychen94 commented 2 years ago

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.