Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

add the Chan-Vese Segmentation for Gray #84

Open
wants to merge 6 commits into
base: master
Choose a base branch
from
Open

add the Chan-Vese Segmentation for Gray #84

Changes from 1 commit
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
96ceac3
add the Chan-Vese Segmentation for Gray
JKay0327 Sep 16, 2021
982cfc4
modify Chan-Vese Segmentation, add docs and simple test
JKay0327 Sep 20, 2021
9996e09
modify runtest.jl
JKay0327 Sep 20, 2021
7e3828d
generic implementation for N dims
JKay0327 Sep 23, 2021
9e07ebc
make the chan_vese function available for 3D image
JKay0327 Sep 24, 2021
6395478
make normalize and init level set be optimal, and some small modifica…
JKay0327 Sep 28, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
252 changes: 252 additions & 0 deletions src/chan_vese.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,252 @@
using MosaicViews
using LazyArrays
using BenchmarkTools
using Images, TestImages
using ImageBase.ImageCore: GenericGrayImage, GenericImage
JKay0327 marked this conversation as resolved.
Show resolved Hide resolved

function calculate_averages(img::AbstractArray{T, N}, H𝚽::AbstractArray{T, M}) where {T<:Real, N, M}
H𝚽ⁱ = @. 1. - H𝚽
∫H𝚽 = sum(H𝚽)
∫H𝚽ⁱ = sum(H𝚽ⁱ)
JKay0327 marked this conversation as resolved.
Show resolved Hide resolved
if ndims(img) == 2
∫u₀H𝚽 = sum(img .* H𝚽)
∫u₀H𝚽ⁱ = sum(img .* H𝚽ⁱ)
elseif ndims(img) == 3
∫u₀H𝚽 = sum(img .* H𝚽, dims=(1, 2))
∫u₀H𝚽ⁱ = sum(img .* H𝚽ⁱ, dims=(1, 2))
end
johnnychen94 marked this conversation as resolved.
Show resolved Hide resolved
if ∫H𝚽 != 0
c₁ = ∫u₀H𝚽 / ∫H𝚽
end
if ∫H𝚽ⁱ != 0
c₂ = ∫u₀H𝚽ⁱ / ∫H𝚽ⁱ
end

return c₁, c₂
JKay0327 marked this conversation as resolved.
Show resolved Hide resolved
end

function difference_from_average_term(img::AbstractArray{T, N}, H𝚽::AbstractArray{T, M}, λ₁::Float64, λ₂::Float64) where {T<:Real, N, M}
c₁, c₂ = calculate_averages(img, H𝚽)

if ndims(img) == 2
return @. -λ₁ * (img - c₁)^2 + λ₂ * (img - c₂)^2
johnnychen94 marked this conversation as resolved.
Show resolved Hide resolved
elseif ndims(img) == 3
return -λ₁ .* sum((img .- c₁).^2, dims=3) .+ λ₂ .* sum((img .- c₂).^2, dims=3)
end
end
# H𝚽 = LazyArray(@~ @. 1. * (𝚽ⁿ > 0))
function _calculate_averages(img::AbstractArray{T, N}, 𝚽ⁿ::AbstractArray{T, M}) where {T<:Real, N, M}
∫H𝚽 = ∫H𝚽ⁱ = ∫u₀H𝚽 = ∫u₀H𝚽ⁱ = 0

for i in CartesianIndices(img)
H𝚽 = 1. * (𝚽ⁿ[i] > 0)
H𝚽ⁱ = 1. - H𝚽
∫H𝚽 += H𝚽
∫H𝚽ⁱ += H𝚽ⁱ
∫u₀H𝚽 += img[i] * H𝚽
∫u₀H𝚽ⁱ += img[i] * H𝚽ⁱ
end
if ∫H𝚽 != 0
c₁ = ∫u₀H𝚽 ./ ∫H𝚽
end
if ∫H𝚽ⁱ != 0
c₂ = ∫u₀H𝚽ⁱ ./ ∫H𝚽ⁱ
end

return c₁, c₂
JKay0327 marked this conversation as resolved.
Show resolved Hide resolved
end

function δₕ(x::AbstractArray{T,N}, h::Float64=1.0) where {T<:Real, N}
return @~ @. h / (h^2 + x^2)
end

function initial_level_set(shape::Tuple)
johnnychen94 marked this conversation as resolved.
Show resolved Hide resolved
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 chan_vese(img::GenericGrayImage;
μ::Float64=0.25,
λ₁::Float64=1.0,
λ₂::Float64=1.0,
tol::Float64=1e-3,
max_iter::Int64=500,
Δt::Float64=0.5,
reinitial_flag::Bool=false) #where {T<:Real, N}
img = float64.(channelview(img))
Copy link
Member

@timholy timholy Sep 16, 2021 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Now I see where you are getting the 3d handling from: it's not actually 3d, it's your color channel. This is almost never needed in JuliaImages; just use the image as provided, and if you've written everything generically it should "just work." That is, for a grayscale image c₁ will be a number or Gray, but for RGB c₁ will be an RGB. And then you can truly support both 2d and 3d images by writing your code in a manner that generalizes across dimensions. So much easier than how other languages do these things, no?

Keep in mind that Python and Julia have reversed "fast axes" for their arrays (Python has its fastest axis last, Julia first), but that in fact the memory layout is the same. Consequently, in Julia if you do call channelview on an RGB image, the color channel is first, not last.

Copy link
Member

@johnnychen94 johnnychen94 Sep 22, 2021 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

iter = 0
h = 1.0
m, n = size(img)
s = m * n
𝚽ⁿ = initial_level_set((m, n)) # size: m * n
del = tol + 1
img .= img .- minimum(img)

if maximum(img) != 0
img .= img ./ maximum(img)
end
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why is this standardization needed?

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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:
cameraman_64

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)

res_64_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)

res_64_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?

Copy link
Member

@johnnychen94 johnnychen94 Sep 25, 2021 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

Copy link
Member

@timholy timholy Sep 25, 2021 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.


diff = 0
H𝚽 = similar(𝚽ⁿ)
u₀H𝚽 = similar(img)
∫u₀ = sum(img)
𝚽ᵢ₊ᶜ = zeros(m, 1)



while (del > tol) & (iter < max_iter)
ϵ = 1e-16

@. H𝚽 = 1. * (𝚽ⁿ > 0) # size = (m, n)
@. u₀H𝚽 = img * H𝚽 # size = (m, n) or (m, n, 3)

∫H𝚽 = sum(H𝚽)
∫u₀H𝚽 = sum(u₀H𝚽) # (1,)
∫H𝚽ⁱ = s - ∫H𝚽
∫u₀H𝚽ⁱ = ∫u₀ - ∫u₀H𝚽

if ∫H𝚽 != 0
c₁ = ∫u₀H𝚽 ./ ∫H𝚽
end
if ∫H𝚽ⁱ != 0
c₂ = ∫u₀H𝚽ⁱ ./ ∫H𝚽ⁱ
end

ind = CartesianIndices(reshape(collect(1 : 9), 3, 3)) .- CartesianIndex(2, 2)
JKay0327 marked this conversation as resolved.
Show resolved Hide resolved
𝚽ⱼ₊ = 0

for y in 1:n-1
JKay0327 marked this conversation as resolved.
Show resolved Hide resolved
𝚽ⱼ₊ = 0
for x in 1:m-1
i = CartesianIndex(x, y)
𝚽₀ = 𝚽ⁿ[i]
u₀ = img[i]
𝚽ᵢ₋ = 𝚽ᵢ₊ᶜ[i[1]]
𝚽ᵢ₊ᶜ[i[1]] = 𝚽ᵢ₊ = 𝚽ⁿ[i + ind[2, 3]] - 𝚽₀ # except i[2] = n
𝚽ⱼ₋ = 𝚽ⱼ₊
𝚽ⱼ₊ = 𝚽ⁿ[i + ind[3, 2]] - 𝚽₀ # except i[2] = m
𝚽ᵢ = 𝚽ᵢ₊ + 𝚽ᵢ₋
𝚽ⱼ = 𝚽ⱼ₊ + 𝚽ⱼ₋
t1 = 𝚽₀ + 𝚽ᵢ₊
t2 = 𝚽₀ - 𝚽ᵢ₋
t3 = 𝚽₀ + 𝚽ⱼ₊
t4 = 𝚽₀ - 𝚽ⱼ₋

C₁ = 1. / sqrt(ϵ + 𝚽ᵢ₊^2 + 𝚽ⱼ^2 / 4.)
JKay0327 marked this conversation as resolved.
Show resolved Hide resolved
C₂ = 1. / sqrt(ϵ + 𝚽ᵢ₋^2 + 𝚽ⱼ^2 / 4.)
C₃ = 1. / sqrt(ϵ + 𝚽ᵢ^2 / 4. + 𝚽ⱼ₊^2)
C₄ = 1. / sqrt(ϵ + 𝚽ᵢ^2 / 4. + 𝚽ⱼ₋^2)

K = t1 * C₁ + t2 * C₂ + t3 * C₃ + t4 * C₄
δₕ = h / (h^2 + 𝚽₀^2)

𝚽ⁿ[i] = 𝚽 = (𝚽₀ + Δt * δₕ * (μ * K - λ₁ * (u₀ - c₁) ^ 2 + λ₂ * (u₀ - c₂) ^ 2)) / (1. + μ * Δt * δₕ * (C₁ + C₂ + C₃ + C₄))
diff += (𝚽 - 𝚽₀)^2
end
i = CartesianIndex(m, y)
𝚽₀ = 𝚽ⁿ[i]
u₀ = img[i]
𝚽ᵢ₋ = 𝚽ᵢ₊ᶜ[i[1]]
𝚽ᵢ₊ᶜ[i[1]] = 𝚽ᵢ₊ = 𝚽ⁿ[i + ind[2, 3]] - 𝚽₀ # except i[2] = n
𝚽ⱼ₋ = 𝚽ⱼ₊
𝚽ⱼ₊ = 0 # except i[2] = m
𝚽ᵢ = 𝚽ᵢ₊ + 𝚽ᵢ₋
𝚽ⱼ = 𝚽ⱼ₊ + 𝚽ⱼ₋
t1 = 𝚽₀ + 𝚽ᵢ₊
t2 = 𝚽₀ - 𝚽ᵢ₋
t3 = 𝚽₀ + 𝚽ⱼ₊
t4 = 𝚽₀ - 𝚽ⱼ₋

C₁ = 1. / sqrt(ϵ + 𝚽ᵢ₊^2 + 𝚽ⱼ^2 / 4.)
C₂ = 1. / sqrt(ϵ + 𝚽ᵢ₋^2 + 𝚽ⱼ^2 / 4.)
C₃ = 1. / sqrt(ϵ + 𝚽ᵢ^2 / 4. + 𝚽ⱼ₊^2)
C₄ = 1. / sqrt(ϵ + 𝚽ᵢ^2 / 4. + 𝚽ⱼ₋^2)

K = t1 * C₁ + t2 * C₂ + t3 * C₃ + t4 * C₄
δₕ = h / (h^2 + 𝚽₀^2)

𝚽ⁿ[i] = 𝚽 = (𝚽₀ + Δt * δₕ * (μ * K - λ₁ * (u₀ - c₁) ^ 2 + λ₂ * (u₀ - c₂) ^ 2)) / (1. + μ * Δt * δₕ * (C₁ + C₂ + C₃ + C₄))
diff += (𝚽 - 𝚽₀)^2
end

𝚽ᵢ₊ = 0
𝚽ⱼ₊ = 0
for x in 1:m-1
i = CartesianIndex(x, n)
𝚽₀ = 𝚽ⁿ[i]
u₀ = img[i]
𝚽ᵢ₋ = 𝚽ᵢ₊ᶜ[i[1]]
𝚽ᵢ₊ᶜ[i[1]] = 0
𝚽ⱼ₋ = 𝚽ⱼ₊
𝚽ⱼ₊ = 𝚽ⁿ[i + ind[3, 2]] - 𝚽₀ # except i[2] = m
𝚽ᵢ = 𝚽ᵢ₊ + 𝚽ᵢ₋
𝚽ⱼ = 𝚽ⱼ₊ + 𝚽ⱼ₋
t1 = 𝚽₀ + 𝚽ᵢ₊
t2 = 𝚽₀ - 𝚽ᵢ₋
t3 = 𝚽₀ + 𝚽ⱼ₊
t4 = 𝚽₀ - 𝚽ⱼ₋

C₁ = 1. / sqrt(ϵ + 𝚽ᵢ₊^2 + 𝚽ⱼ^2 / 4.)
C₂ = 1. / sqrt(ϵ + 𝚽ᵢ₋^2 + 𝚽ⱼ^2 / 4.)
C₃ = 1. / sqrt(ϵ + 𝚽ᵢ^2 / 4. + 𝚽ⱼ₊^2)
C₄ = 1. / sqrt(ϵ + 𝚽ᵢ^2 / 4. + 𝚽ⱼ₋^2)

K = t1 * C₁ + t2 * C₂ + t3 * C₃ + t4 * C₄
δₕ = h / (h^2 + 𝚽₀^2)

𝚽ⁿ[i] = 𝚽 = (𝚽₀ + Δt * δₕ * (μ * K - λ₁ * (u₀ - c₁) ^ 2 + λ₂ * (u₀ - c₂) ^ 2)) / (1. + μ * Δt * δₕ * (C₁ + C₂ + C₃ + C₄))
diff += (𝚽 - 𝚽₀)^2
end
i = CartesianIndex(m, n)
𝚽₀ = 𝚽ⁿ[i]
u₀ = img[i]
𝚽ᵢ₋ = 𝚽ᵢ₊ᶜ[i[1]]
𝚽ᵢ₊ᶜ[i[1]] = 0
𝚽ⱼ₋ = 𝚽ⱼ₊
𝚽ⱼ₊ = 0
𝚽ᵢ = 𝚽ᵢ₊ + 𝚽ᵢ₋
𝚽ⱼ = 𝚽ⱼ₊ + 𝚽ⱼ₋
t1 = 𝚽₀ + 𝚽ᵢ₊
t2 = 𝚽₀ - 𝚽ᵢ₋
t3 = 𝚽₀ + 𝚽ⱼ₊
t4 = 𝚽₀ - 𝚽ⱼ₋

C₁ = 1. / sqrt(ϵ + 𝚽ᵢ₊^2 + 𝚽ⱼ^2 / 4.)
C₂ = 1. / sqrt(ϵ + 𝚽ᵢ₋^2 + 𝚽ⱼ^2 / 4.)
C₃ = 1. / sqrt(ϵ + 𝚽ᵢ^2 / 4. + 𝚽ⱼ₊^2)
C₄ = 1. / sqrt(ϵ + 𝚽ᵢ^2 / 4. + 𝚽ⱼ₋^2)

K = t1 * C₁ + t2 * C₂ + t3 * C₃ + t4 * C₄
δₕ = h / (h^2 + 𝚽₀^2)

𝚽ⁿ[i] = 𝚽 = (𝚽₀ + Δt * δₕ * (μ * K - λ₁ * (u₀ - c₁) ^ 2 + λ₂ * (u₀ - c₂) ^ 2)) / (1. + μ * Δt * δₕ * (C₁ + C₂ + C₃ + C₄))
diff += (𝚽 - 𝚽₀)^2

del = sqrt(diff / s)
diff = 0

iter += 1
end

return 𝚽ⁿ, iter
end

img_gray = testimage("cameraman")
JKay0327 marked this conversation as resolved.
Show resolved Hide resolved

μ=0.25
λ₁=1.0
λ₂=1.0
tol=1e-3
max_iter=200
Δt=0.5

𝚽, iter_num = chan_vese(img_gray, μ=0.25, λ₁=1.0, λ₂=1.0, tol=1e-3, max_iter=200, Δt=0.5, reinitial_flag=false)

@btime chan_vese(img_gray, μ=0.25, λ₁=1.0, λ₂=1.0, tol=1e-3, max_iter=200, Δt=0.5, reinitial_flag=false);

segmentation = 𝚽 .> 0
print(iter_num)
𝚽 .= 𝚽 .- minimum(𝚽)

colorview(Gray, segmentation)