add the Chan-Vese Segmentation for Gray

[JKay0327]

This is an implementation of Chan-Vese Segmentation without vectorized codes.

[JKay0327]

Different from my previous work: https://github.com/JKay0327/Chan-Vese, I try to modified the vectorized codes by using for cycle. What's more, I find that there appears a lot of repeat calculations when calculate the forward diff and backward diff, so I've done some works to eliminate the repeat calculations.
I couldn't make sure whether CV works well on RGB images, so this version just available for Gray images.

[timholy]

Here's a preliminary review. It took me a little while to understand your intent in several places, so my early comments are a bit confused, particularly about the role of your 3rd dimension. As far as I can tell, there's nothing in the mathematics that limits Chan-Vese to 2d, and in JuliaImages we can transparently support grayscale & color along with arbitrary spatial dimensionality with no ambiguity about what these things mean. That's different from other languages where an m, n, 3 array means "RGB" except when it doesn't (e.g., if it's a stack of 3 slices from a 3d grayscale). ColorVectorSpace defines the arithmetic operations you'll need for RGB.

Another question: this is a very simple stepping solver. I wonder if it would make sense to split out the optimization to allow application of solvers from, e.g., Optim.jl or the DifferentialEquations suite? I haven't thought about what that might look like, so if it's a lot of work we could do that in a separate PR.

[timholy]

Why is this standardization needed?

[JKay0327]

About standardization, I think a standardized image may contribute to a segmentation algorithm since it makes the gray levels of image more dispersive. If the grey levels of image range from 0.2 to 0.8, it seems better to standardize image so that the grey levels will range from 0.0 to 1.0.
The test of chan_vese use a resized 64*64 cameraman as the test image, whose grey levels range from 0.02 to 0.988. The standardization has showed that it will cause some difference :

julia> using TestImages, Images

julia> img_512 = testimage("cameraman");

julia> maximum(img_512)
Gray{N0f8}(1.0)

julia> minimum(img_512)
Gray{N0f8}(0.0)

julia> img_64 = imresize(testimage("cameraman"), (64, 64));

julia> maximum(img)
Gray{N0f8}(0.988)

julia> minimum(img)
Gray{N0f8}(0.02)

Test image shows as following:

julia> res_no_std = chan_vese(img, μ=0.1, tol=1e-2, max_iter=200);

julia> sum(res_no_std)
1075

julia> colorview(Gray, res_no_std)

julia> res_std = chan_vese(img, μ=0.1, tol=1e-2, max_iter=200);

julia> sum(res_std)
1064

julia> colorview(Gray, res_std)

We can find that there are some difference. Without standardization, the part between the man's leg can't be segmented well. So maybe we have to maintain the standardization?

[johnnychen94]

Maybe we can introduce a keyword normalize=false to control this behavior. And we can also refer to ImageContrastAdjustment in the docstring for more complicated normalization workflows.

By hardcoding the standardization in the implementation we lose some possibilities. In many similar cases, we prefer to let the caller rather than the callee do the work. See also https://docs.julialang.org/en/v1/manual/style-guide/#Handle-excess-argument-diversity-in-the-caller

[timholy]

Presumably that is easily compensated by changing the values of λ. The only place in the mathematics where the image magnitude means anything is the λ*sum(abs2, img[componentmask] .- componentmean). So anything you can do by standardizing the image, you can do more quickly by adjusting λ. Alternatively, you can adjust μ since only the ratio of the terms really matters for the final result.

What does standardization even mean for an image with RGB pixels? Do you take the min & max across channels too? Why does the magnitude of the blue channel scale the magnitude of the red one? Aren't these independent pieces of information? Conversely, if my raw image has almost no diversity in the green channel, I'll be really surprised if standardization on the green channel converts something that is perceptually insignificant into something that drives the segmentation.

Basically, there isn't an answer to such questions. So it's better to the let the user be in charge.

[johnnychen94]

I did not review calculate_reinitial part because I think it's quite duplicated to the main body of chan_vese. I believe we need a more generic and programmable way to handle i₊, j₊, i₋, j₋ and a, b, c, d things; the current version is hardcoded to a 2-dimensional case and it seems not necessary and not elegant to me.

[johnnychen94]

For test coverage, also need to test reinitial_flag=true case.

[johnnychen94]

If you want a raw numerical array, then just:

Suggested change

	img = float64.(channelview(img))
	img = Float64.(img)

But generally, we encourage directly handling Array{<:Colorant} without channelview because it permits more generic implementation as Tim said in the previous comment.

[johnnychen94]

Every dependency requires a corresponding compat entry.

[johnnychen94]

We should either remove ImageCore from Project.toml or use

Suggested change

	using ImageBase.ImageCore: GenericGrayImage, GenericImage
	using ImageCore: GenericGrayImage, GenericImage

[JKay0327]

For 3D image, I take image mri in TestImages as an example, here is the demo:

julia> using Images, TestImages
julia> mri = testimage("mri")
julia> mosaicview(mri, ncol=9)

julia> mri_seg = chan_vese(mri, μ=0.005, λ₁=0.8, λ₂=1.0, tol=1e-3, max_iter=2000, Δt=0.4, reinitial_flag=false)

It seems that chan_vese already works well when input image has size 226 * 186 * 27.

I find it reach the convergence when number of iteration comes to 1852. So it works well on 3D image actually!

[johnnychen94]

The extension to 3d case is great and I think it's worth a publication. Here's a very simple illustration on its superiority over frame-wise 2d processing:

using ImageSegmentation, TestImages, ImageShow, ImageCore
using StackViews

mri = Gray{Float64}.(testimage("mri"))
mri .+= 0.1 .* randn(size(mri)) # add small gaussian noise
mosaicview(mri, nrow=3, rowmajor=true)

mri_seg = chan_vese(mri, μ=0.005, λ₁=0.8, λ₂=1.0, tol=1e-3, max_iter=200, Δt=0.4, reinitial_flag=false);
mri_seg_per_frame =  chan_vese.(eachslice(mri, dims=3), μ=0.005, λ₁=0.8, λ₂=1.0, tol=1e-3, max_iter=200, Δt=0.4, reinitial_flag=false);

mosaicview(mri_seg, StackView(mri_seg_per_frame), nrow=6, rowmajor=true) .|> Gray

and if we check a slice:

mosaicview(mri_seg[70, :, :]', StackView(mri_seg_per_frame)[70, :, :]'; npad=10) .|> Gray

Quite obviously that directly processing 3d case produces more consistent and smooth result.

[timholy]

This is starting to look merge-worthy to me! Great work, @JKay0327.

The most important remaining tasks seem to be:

• generalization by eltype, e.g., to RGB images. Mostly, you can do that by deleting code, which is always nice. I've indicated a couple of changes that I think will help, but it is almost certain you'll have to solve a few other problems.
• Adding more tests for different eltype gray images, RGB images, and 2d and 3d images

[timholy]

Suggested change

	function initial_level_set(shape::Tuple{Int64, Int64})
	x₀ = reshape(collect(0:shape[begin]-1), shape[begin], 1)
	y₀ = reshape(collect(0:shape[begin+1]-1), 1, shape[begin+1])
	𝚽₀ = @. sin(pi / 5 * x₀) * sin(pi / 5 * y₀)
	end
	function initial_level_set(shape::Tuple{Int64, Int64, Int64})
	x₀ = reshape(collect(0:shape[begin]-1), shape[begin], 1, 1)
	y₀ = reshape(collect(0:shape[begin+1]-1), 1, shape[begin+1], 1)
	z₀ = reshape(collect(0:shape[begin+2]-1), 1, 1, shape[begin+2])
	𝚽₀ = @. sin(pi / 5 * x₀) * sin(pi / 5 * y₀) * sin(pi / 5 * z₀)
	end
	initial_level_set(sz) = broadcast(*, kernelfactors(ntuple(i -> sin.((pi / 5) * (0:sz[i]-1)), length(sz)))

kernelfactors is from ImageFiltering and will make this truly arbitrary-dimensional.

[JKay0327]

I’ve tried to change the code, but do I need to change the code like following to get a Matrix type output ?

initial_level_set(sz) = broadcast(*, kernelfactors(ntuple(i -> sin.((pi / 5) * (0:sz[i]-1)), length(sz)))...)

Another question is that this changing seems reduce the performance:

# previous code's performance
function initial_level_set(shape::Tuple{Int64, Int64})
    x₀ = reshape(collect(0:shape[begin]-1), shape[begin], 1)
    y₀ = reshape(collect(0:shape[begin+1]-1), 1, shape[begin+1])
    𝚽₀ = @. sin(pi / 5 * x₀) * sin(pi / 5 * y₀)
end
...
julia> @btime chan_vese(img_gray, μ=0.25, λ₁=1.0, λ₂=1.0, tol=1e-3, max_iter=200, Δt=0.5, reinitial_flag=false);
  462.876 ms (119 allocations: 6.46 MiB)

# using kernelfactor
initial_level_set(sz) = broadcast(*, kernelfactors(ntuple(i -> sin.((pi / 5) * (0:sz[i]-1)), length(sz)))...)
...
julia> @btime chan_vese(img_gray, μ=0.25, λ₁=1.0, λ₂=1.0, tol=1e-3, max_iter=200, Δt=0.5, reinitial_flag=false);
  575.990 ms (115 allocations: 6.46 MiB)

It seems that the changing causes more than 100ms's performance reduction.

[timholy]

Shoot, that should have been

initial_level_set(sz) = broadcast(*, kernelfactors(ntuple(i -> sin.((pi / 5) * (0:sz[i]-1)), length(sz)))...)

so that

julia> initial_level_set((3, 5))
3×5 OffsetArray(::Matrix{Float64}, 1:3, 1:5) with eltype Float64 with indices 1:3×1:5:
 0.0  0.0       0.0       0.0       0.0
 0.0  0.345492  0.559017  0.559017  0.345492
 0.0  0.559017  0.904508  0.904508  0.559017

julia> initial_level_set((3, 5, 2))
3×5×2 OffsetArray(::Array{Float64, 3}, 1:3, 1:5, 1:2) with eltype Float64 with indices 1:3×1:5×1:2:
[:, :, 1] =
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0

[:, :, 2] =
 0.0  0.0       0.0       0.0       0.0
 0.0  0.203075  0.328582  0.328582  0.203075
 0.0  0.328582  0.531657  0.531657  0.328582

Sorry about that.

[timholy]

I would recommend deleting this line and the minimum/maximum standardization unless there is strong reason to keep it. Other languages may have done this just to commonize among UInt8 images (ranging from 0:255), UInt16 images (ranging from 0:65535), and Float32/Float64 images (ranging from 0 to 1). But thanks to N0f8 and N0f16, we've already done that. I think it's better to accept the image at face value.

[timholy]

If you like, you can even make this

Suggested change

	Δ = ntuple(d -> CartesianIndex(ntuple(i -> i == d ? 1 : 0, N)), N)
	Δ = ntuple(d -> CartesianIndex(ntuple(i -> i == d, N)), N)

since the Bool will be promoted to Int. (Optional, of course, and you may not like it better). I've probably got places I could use that myself...

[timholy]

Suggested change

	function calculate_averages(img::AbstractArray{T, N}, H𝚽::AbstractArray{S, N}, area::Int64, ∫u₀::Float64) where {T<:Real, S<:Bool, N}
	function calculate_averages(img::AbstractArray{T, N}, H𝚽::AbstractArray{S, N}, area, ∫u₀) where {T<:Real, S<:Bool, N}

Putting these types in doesn't help you, and only prevents you from using, e.g., RGB colors for ∫u₀.

[timholy]

Suggested change

	∫u₀H𝚽 = 0
	∫u₀H𝚽 = zero(∫u₀)

This will make you generic.

[timholy]

Comments and test-coverage would be appreciated here. I'm not quite sure what this is about (I haven't yet seen this in the paper but I haven't read the whole thing carefully, just skimmed it).

[timholy]

Other optional factors to consider:

• do we want to allow the user to pass in an initialization for 𝚽? For grayscale, I could imagine wanting to initialize it using Otsu's method, and a kmeans with 2 clusters for RGB images.
• splitting this into clearly-defined "objective function" and "solver" might make it possible for people to experiment with other solvers. That might facilitate convergence in the 3d case?

[timholy]

@johnnychen94 that's a great comparison. Some might argue that the interpretation of the parameter values might need to differ between 2 and 3 dimensions (for example, the ratio of surface area-to-volume differs), so for a truly fair comparison you might need to tweak the parameters (at least μ) until you get the best match between the two results. Nevertheless, I am confident based on prior experience that the fundamental observation will hold: you'll get better, more consistent results in a natively 3d algorithm.

[JKay0327]

This is a summary for works we've done in this PR:

• Implement function chan_vese. We use for loop to solve the pixel wide problem, so that we got a 7~8 times performance comparing with scikit-image's implementation.
• Make functions in chan_vese available for N dimensions so that we can use chan_vese on 3D images.

The current performance of chan_vese algorithm in Julia is as follow:

julia> using Images, TestImages

julia> using ImageSegmentation

julia> img2d = testimage("cameraman");

julia> @btime chan_vese(img2d, μ=0.25, λ₁=1.0, λ₂=1.0, tol=1e-3, max_iter=200, Δt=0.5, reinitial_flag=false);
462.876 ms (119 allocations: 6.46 MiB)

The performance of sciki-image's implementation is as follow:

image = img_as_float(data.camera())
start =time.clock()
for i in range(100):
    cv = chan_vese(image, mu=0.25, lambda1=1, lambda2=1, tol=1e-3, max_iter=200, dt=0.5, init_level_set="checkerboard", extended_output=True)
end = time.clock()
print('Running time: %s Seconds'%((end-start)/100))

Running time: 3.2634375 Seconds

Using chan_vese on 2D input:

julia> using Images, TestImages, MosaicViews

julia> using ImageSegmentation

julia> img2d = testimage("cameraman");

julia> size(img2d) # 2D image
(512, 512)

julia> img2d_seg = chan_vese(img2d, max_iter=200);

julia> mosaicview(img_gray, segmentation, nrow=1, rowmajor=true)

Using chan_vese on 3D input:

julia> img3d = testimage("mri");

julia> size(img3d) # 3D image
(226, 186, 27)

julia> img3d_seg = chan_vese(img3d, μ=0.005, max_iter=2000);

julia> mosaicview(mri, segmentation, ncol=9, rowmajor=true)

Some problems are still remain, we can solve these problems in next PR:

• Further improve the performance of chan_vese.
• Study whether to split chan_vese into "objective function" and "solver".

[JKay0327]

A list to do before merge:

• Add tests for init_level_set and 3D inputs.
• RGB images should be available. But the algorithm for colorant images is different from the one for 3D image. So we can not just use an RGB type input, we have to change the algorithm.

If there is anything I have omitted, please tell me and I'll extend this todo-list.

Unresolved comments: (Edited by: @johnnychen94)

• Investigate why float64(channelview(img)) is needed compared to float.(img): add the Chan-Vese Segmentation for Gray #84 (comment)
• test both normalize=true and normalize=false and then resolves add the Chan-Vese Segmentation for Gray #84 (comment)
• test reinitial_flag=true case and then resolves add the Chan-Vese Segmentation for Gray #84 (comment) and add the Chan-Vese Segmentation for Gray #84 (comment)
• add missing dependency compat for ImageBase add the Chan-Vese Segmentation for Gray #84 (comment)
• maybe cleanup ImageCore dependency in favor of ImageBase add the Chan-Vese Segmentation for Gray #84 (comment)
• dimension-agnoistic implementation of initial_level_set add the Chan-Vese Segmentation for Gray #84 (comment) and add tests for it
• remove useless type restriction: add the Chan-Vese Segmentation for Gray #84 (comment)
• type-awared zero: add the Chan-Vese Segmentation for Gray #84 (comment)

After merging this:

• add a demo for this algorithm to https://juliaimages.org/stable/examples/ (Edited by: @johnnychen94)

[johnnychen94]

Need to update the keywords for normalize and init_level_set.

[johnnychen94]

Suggested change

	## normalize::Bool
	The arguement `normalize` controls whether to normalize the input img.
	Default: false
	## `normalize::Bool`
	The argument `normalize` controls whether to normalize the input `img` by linear scaling it to `[0.0, 1.0]` range.
	Default: false

[johnnychen94]

@simd won't work if the function is complicated (e.g., has if-branch and uninlined functions)

Suggested change

	@inbounds @simd for idx in CartesianIndices(𝚽ⁿ)
	@inbounds for idx in CartesianIndices(𝚽ⁿ)

[johnnychen94]

Might be worth documenting that this algorithm supports 3d images.

If you don't plan to add RGB support in this PR then might also worth mentioning that RGB image is not supported.

[timholy]

You've brought this to an excellent state, and both the generalization to 3d and the performance advantage over skimage are gratifying. WRT the performance advantage, can you clarify whether both converged after the same number of iterations? It's important because different algorithms might, in principle, simply have different tolerances on convergence (I haven't checked), and it would not be fair to penalize skimage if the only difference is that it chooses to reach a higher standard before termination. Just checking the numeric value of the tolerance is sometimes insufficient unless you're certain that the two are being applied to the same measure (e.g., identical objective function). If anything is unclear about the comparison, you could try forcing early convergence in both after, say, 10 iterations, as in that case it seems almost certain that both will use the full 10 iterations (you could be certain if you can show their timing is linear in the number of iterations you allow) and then you know you'll be using comparable timing data.

we can solve these problems in next PR:
Further improve the performance of chan_vese.
Study whether to split chan_vese into "objective function" and "solver".

I agree that both of these are better left for a future PR.

But the algorithm for colorant images is different from the one for 3D image

Can you elaborate on this? AFAICT this is merely an issue of using ColorVectorSpace.Future.abs2 (#84 (comment)) rather than (x-y)^2, dropping/modifying a couple of type annotations or single lines (#84 (comment), #84 (comment), and deleting img = float64.(channelview(img))). That's just dropping some unnecessarily-specific choices, could be done with a couple minutes of work, and if it suffices it's a big advance. (If it turns out to be more complicated, then I agree it's better left for a future PR.)

If there is anything I have omitted, please tell me and I'll extend this todo-list.

• Any reason not to do https://github.com/JuliaImages/ImageSegmentation.jl/pull/84/files#r716011493? If there is, I'm fine with it, but worth checking whether it's an oversight.
• The whole "reinitialize" framework seems untested and partially commented out. It might be best to fully comment it out and remove the reference in the docstring for now.

Hah, I see @johnnychen94 and I are almost simultaneous again!

[JKay0327]

Any reason not to do https://github.com/JuliaImages/ImageSegmentation.jl/pull/84/files#r716011493? If there is, I'm fine with it, but worth checking whether it's an oversight.

It is shown that this will cause a nearly 100ms's performance reduction, and I have give out a benchmark in https://github.com/JuliaImages/ImageSegmentation.jl/pull/84/files#r716011493.

[timholy]

Sorry if the sequence was unclear. I made a suggestion, you pointed out the problem, I realized I had made an error and amended the suggestion. The performance should be identical because it produces the same output as your solution. for substantially less code while also being fully general across dimensionality.

EDIT: now I see you already added the missing .... Sorry I failed to understand the point you were trying to make.

[JKay0327]

In my test, I've already change the code into:

# using kernelfactor
initial_level_set(sz) = broadcast(*, kernelfactors(ntuple(i -> sin.((pi / 5) * (0:sz[i]-1)), length(sz)))...)
...
julia> @btime chan_vese(img_gray, μ=0.25, λ₁=1.0, λ₂=1.0, tol=1e-3, max_iter=200, Δt=0.5, reinitial_flag=false);
  575.990 ms (115 allocations: 6.46 MiB)

It produces the same output, but gets a worse performance.
Do you means that we have to first find a way to make initial_level_set fully general across dimensionality, then consider the performance of this method?

[timholy]

EDIT: now I see that you had already fixed my first concern in your initial post on this topic, so I failed to realize that the performance problem you were worried about wasn't an artifact. Sorry about the confusion. The text below should still be useful, however.

So this is what makes it surprising:

julia> using ImageFiltering

julia> function initial_level_set_yours(shape::Tuple{Int64, Int64})
           x₀ = reshape(collect(0:shape[begin]-1), shape[begin], 1)
           y₀ = reshape(collect(0:shape[begin+1]-1), 1, shape[begin+1])
           𝚽₀ = @. sin(pi / 5 * x₀) * sin(pi / 5 * y₀)
       end
initial_level_set_yours (generic function with 1 method)

julia> initial_level_set_mine(sz) = broadcast(*, kernelfactors(ntuple(i -> sin.((pi / 5) * (0:sz[i]-1)), length(sz)))...)
initial_level_set_mine (generic function with 1 method)

julia> ϕyours = initial_level_set_yours((5, 8));

julia> ϕmine = initial_level_set_mine((5, 8));

julia> ϕyours == ϕmine
true

However

julia> println(summary(ϕyours))
5×8 Matrix{Float64}

julia> println(summary(ϕmine))
5×8 OffsetArray(::Matrix{Float64}, 1:5, 1:8) with eltype Float64 with indices 1:5×1:8

so the types are different. This suggests you have an inferrability problem (slower and more allocations are typical symptoms). While the best fix is to solve the inferrability problem, in the short term perhaps you can use parent:

julia> println(summary(parent(ϕmine)))
5×8 Matrix{Float64}

which gives the exact same type, too.

Do you means that we have to first find a way to make initial_level_set fully general across dimensionality, then consider the performance of this method?

There's no absolute obligation, but there seems to be no reason to avoid being general (shorter and more general code is typically preferred).