diff --git a/Project.toml b/Project.toml index a74201a..af3c74e 100644 --- a/Project.toml +++ b/Project.toml @@ -6,6 +6,7 @@ version = "1.6.0" Clustering = "aaaa29a8-35af-508c-8bc3-b662a17a0fe5" DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" Distances = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7" +ImageBase = "c817782e-172a-44cc-b673-b171935fbb9e" ImageCore = "a09fc81d-aa75-5fe9-8630-4744c3626534" ImageFiltering = "6a3955dd-da59-5b1f-98d4-e7296123deb5" ImageMorphology = "787d08f9-d448-5407-9aad-5290dd7ab264" @@ -38,6 +39,8 @@ FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549" ImageMagick = "6218d12a-5da1-5696-b52f-db25d2ecc6d1" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" +TestImages = "5e47fb64-e119-507b-a336-dd2b206d9990" +ImageTransformations = "02fcd773-0e25-5acc-982a-7f6622650795" [targets] -test = ["Documenter", "FileIO", "ImageMagick", "SparseArrays", "Test"] +test = ["Documenter", "FileIO", "ImageMagick", "SparseArrays", "Test", "TestImages", "ImageTransformations"] diff --git a/src/ImageSegmentation.jl b/src/ImageSegmentation.jl index d81cb7f..99e5cfe 100644 --- a/src/ImageSegmentation.jl +++ b/src/ImageSegmentation.jl @@ -5,6 +5,7 @@ import Base: show using LinearAlgebra, Statistics using DataStructures, StaticArrays, ImageCore, ImageFiltering, ImageMorphology, LightGraphs, SimpleWeightedGraphs, RegionTrees, Distances, StaticArrays, Clustering, MetaGraphs using ImageCore.ColorVectorSpace: MathTypes +using ImageBase.ImageCore: GenericGrayImage, GenericImage import Clustering: kmeans, fuzzy_cmeans const PairOrTuple{K,V} = Union{Pair{K,V},Tuple{K,V}} @@ -21,6 +22,7 @@ include("meanshift.jl") include("clustering.jl") include("merge_segments.jl") include("deprecations.jl") +include("chan_vese.jl") export #accessor methods @@ -49,6 +51,7 @@ export kmeans, fuzzy_cmeans, merge_segments, + chan_vese, # types SegmentedImage, diff --git a/src/chan_vese.jl b/src/chan_vese.jl new file mode 100644 index 0000000..a3ec7c0 --- /dev/null +++ b/src/chan_vese.jl @@ -0,0 +1,254 @@ +""" + chan_vese(img; [μ], [λ₁], [λ₂], [tol], [max_iter], [Δt], [reinitial_flag]) + +Segments image `img` by evolving a level set. An active contour model +which can be used to segment objects without clearly defined boundaries. + +# output +Return a `BitMatrix`. + +# Details + +Chan-Vese algorithm deals quite well even with images which are quite +difficult to segment. Since CV algorithm relies on global properties, +rather than just taking local properties under consideration, such as +gradient. Better robustness for noise is one of the main advantages of +this algorithm. See [1], [2], [3] for more details. + +# Options + +The function argument is described in detail below. + +Denote the edge set curve with 𝐶 in the following part. + +## `μ::Real` + +The argument `μ` is a weight controlling the penalty on the total length +of the curve 𝐶; + +For example, if the boundaries of the image are quite smooth, a larger `μ` +can prevent 𝐶 from being a complex curve. + +Default: 0.25 + +## `λ₁::Real`, `λ₂::Real` + +The argument `λ₁` and `λ₂` affect the desired uniformity inside 𝐶 and +outside 𝐶, respectively. + +For example, if set `λ₁` < `λ₂`, we are more possible to get result with +quite uniform background and varying grayscale objects in the foreground. + +Default: λ₁ = 1.0 + λ₂ = 1.0 + +## `tol::Real` + +The argument `tol` controls the level set variation tolerance between +iteration. If the L2 norm difference between two level sets of adjacent +iterations is below `tol`, then the solution will be assumed to be reached. + +Default: 1e-3 + +## `max_iter::Int` + +The argument `max_iter` controls the maximum of iteration number. + +Default: 500 + +## `Δt::Real` + +The argument `Δt` is a multiplication factor applied at calculations +for each step, serves to accelerate the algorithm. Although larger `Δt` +can speed up the algorithm, it might prevent algorithm from converging to +the solution. + +Default: 0.5 + +## normalize::Bool + +The arguement `normalize` controls whether to normalize the input img. + +Default: false + +## init_level_set + +The arguement `init_level_set` allows users to initalize the level set 𝚽⁰, +whose size should be the same as input image. + +If the default init level set is uesd, it will be defined as: +𝚽⁰ₓ = ∏ sin(xᵢ ⋅ π / 5), where xᵢ are pixel coordinates, i = 1, 2, ⋯, ndims(img). +This level set has fast convergence, but may fail to detect implicit edges. + +Default: initial_level_set(size(img)) + +# Examples + +```julia +using TestImages +using ImageSegmentation + +img = testimage("cameraman") + +cv_result = chan_vese(img, max_iter=200) +``` + +# References + +[1] An Active Contour Model without Edges, Tony Chan and Luminita Vese, + Scale-Space Theories in Computer Vision, 1999, :DOI:`10.1007/3-540-48236-9_13` +[2] Chan-Vese Segmentation, Pascal Getreuer Image Processing On Line, 2 (2012), + pp. 214-224, :DOI:`10.5201/ipol.2012.g-cv` +[3] The Chan-Vese Algorithm - Project Report, Rami Cohen, 2011 :arXiv:`1107.2782` +""" +function chan_vese(img::GenericGrayImage; + μ::Real=0.25, + λ₁::Real=1.0, + λ₂::Real=1.0, + tol::Real=1e-3, + max_iter::Int=500, + Δt::Real=0.5, + normalize::Bool=false, + init_level_set=initial_level_set(size(img))) + # Signs used in the codes and comments mainly follow paper[3] in the References. + axes(img) == axes(init_level_set) || throw(ArgumentError("axes of input image and init_level_set should be equal. Instead they are $(axes(img)) and $(axes(init_level_set)).")) + img = float64.(channelview(img)) + N = ndims(img) + iter = 0 + h = 1.0 + del = tol + 1 + if normalize + img .= img .- minimum(img) + if maximum(img) != 0 + img .= img ./ maximum(img) + end + end + + # Precalculation of some constants which helps simplify some integration + area = length(img) # area = ∫H𝚽 + ∫H𝚽ⁱ + ∫u₀ = sum(img) # ∫u₀ = ∫u₀H𝚽 + ∫u₀H𝚽ⁱ + + # Initialize the level set + 𝚽ⁿ = init_level_set + + # Preallocation and initializtion + H𝚽 = trues(size(img)...) + 𝚽ⁿ⁺¹ = similar(𝚽ⁿ) + + Δ = ntuple(d -> CartesianIndex(ntuple(i -> i == d ? 1 : 0, N)), N) + idx_first = first(CartesianIndices(𝚽ⁿ)) + idx_last = last(CartesianIndices(𝚽ⁿ)) + + while (del > tol) & (iter < max_iter) + ϵ = 1e-8 + diff = 0 + + # Calculate the average intensities + @. H𝚽 = 𝚽ⁿ > 0 # Heaviside function + c₁, c₂ = calculate_averages(img, H𝚽, area, ∫u₀) # Compute c₁(𝚽ⁿ), c₂(𝚽ⁿ) + + # Calculate the variation of level set 𝚽ⁿ + @inbounds @simd for idx in CartesianIndices(𝚽ⁿ) + 𝚽₀ = 𝚽ⁿ[idx] # 𝚽ⁿ(x, y) + u₀ = img[idx] # u₀(x, y) + 𝚽₊ = broadcast(i->𝚽ⁿ[i], ntuple(d -> idx[d] != idx_last[d] ? idx + Δ[d] : idx, N)) + 𝚽₋ = broadcast(i->𝚽ⁿ[i], ntuple(d -> idx[d] != idx_first[d] ? idx - Δ[d] : idx, N)) + + # Solve the PDE of equation 9 in paper[3] + center_diff = ntuple(d -> (𝚽₊[d] - 𝚽₋[d])^2 / 4., N) + sum_center_diff = sum(center_diff) + C₊ = ntuple(d -> 1. / sqrt(ϵ + (𝚽₊[d] - 𝚽₀)^2 + sum_center_diff - center_diff[d]), N) + C₋ = ntuple(d -> 1. / sqrt(ϵ + (𝚽₋[d] - 𝚽₀)^2 + sum_center_diff - center_diff[d]), N) + + K = sum(𝚽₊ .* C₊) + sum(𝚽₋ .* C₋) + δₕ = h / (h^2 + 𝚽₀^2) # Regularised Dirac function + difference_from_average = - λ₁ * (u₀ - c₁) ^ 2 + λ₂ * (u₀ - c₂) ^ 2 + + 𝚽ⁿ⁺¹[idx] = 𝚽 = (𝚽₀ + Δt * δₕ * (μ * K + difference_from_average)) / (1. + μ * Δt * δₕ * (sum(C₊) + sum(C₋))) + diff += (𝚽 - 𝚽₀)^2 + end + + del = sqrt(diff / area) + + # Reinitializing the level set is not strictly necessary, so this part of code is commented. + # If you wants to use the reinitialization, just uncommented codes following. + # Function `reinitialize!` is already prepared. + + # reinitialize!(𝚽ⁿ⁺¹, 𝚽ⁿ, Δt, h) # Reinitialize 𝚽 to be the signed distance function to its zero level set + + 𝚽ⁿ .= 𝚽ⁿ⁺¹ + + iter += 1 + end + + return 𝚽ⁿ .> 0 +end + +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 + +function calculate_averages(img::AbstractArray{T, N}, H𝚽::AbstractArray{S, N}, area::Int64, ∫u₀::Float64) where {T<:Real, S<:Bool, N} + ∫u₀H𝚽 = 0 + ∫H𝚽 = 0 + @inbounds for i in eachindex(img) + if H𝚽[i] + ∫u₀H𝚽 += img[i] + ∫H𝚽 += 1 + end + end + ∫H𝚽ⁱ = area - ∫H𝚽 + ∫u₀H𝚽ⁱ = ∫u₀ - ∫u₀H𝚽 + c₁ = ∫u₀H𝚽 / max(1, ∫H𝚽) + c₂ = ∫u₀H𝚽ⁱ / max(1, ∫H𝚽ⁱ) + + return c₁, c₂ +end + +function calculate_reinitial(𝚽::AbstractArray{T, M}, 𝚿::AbstractArray{T, M}, Δt::Float64, h::Float64) where {T<:Real, M} + ϵ = 1e-8 + N = ndims(𝚽) + + Δ = ntuple(d -> CartesianIndex(ntuple(i -> i == d ? 1 : 0, N)), N) + idx_first = first(CartesianIndices(𝚽)) + idx_last = last(CartesianIndices(𝚽)) + + @inbounds @simd for idx in CartesianIndices(𝚽) + 𝚽₀ = 𝚽[idx] # 𝚽ⁿ(x, y) + Δ₊ = ntuple(d -> idx[d] != idx_last[d] ? idx + Δ[d] : idx, N) + Δ₋ = ntuple(d -> idx[d] != idx_first[d] ? idx - Δ[d] : idx, N) + Δ𝚽₊ = broadcast(i -> (𝚽[i] - 𝚽₀) / h, Δ₊) + Δ𝚽₋ = broadcast(i -> (𝚽₀ - 𝚽[i]) / h, Δ₋) + + maxΔ𝚽₊ = max.(Δ𝚽₊, 0) + minΔ𝚽₊ = min.(Δ𝚽₊, 0) + maxΔ𝚽₋ = max.(Δ𝚽₋, 0) + minΔ𝚽₋ = min.(Δ𝚽₋, 0) + + G = 0 + if 𝚽₀ > 0 + G += sqrt(sum(max.(minΔ𝚽₊.^2, maxΔ𝚽₋.^2))) - 1 + elseif 𝚽₀ < 0 + G += sqrt(sum(max.(maxΔ𝚽₊.^2, minΔ𝚽₋.^2))) - 1 + end + sign𝚽 = 𝚽₀ / sqrt(𝚽₀^2 + ϵ) + 𝚿[idx] = 𝚽₀ - Δt * sign𝚽 * G + end + + return 𝚿 +end + +function reinitialize!(𝚽::AbstractArray{T, M}, 𝚿::AbstractArray{T, M}, Δt::Float64, h::Float64, max_reiter::Int=5) where {T<:Real, M} + for i in 1 : max_reiter + 𝚽 .= calculate_reinitial(𝚽, 𝚿, Δt, h) + end +end \ No newline at end of file diff --git a/test/chan_vese.jl b/test/chan_vese.jl new file mode 100644 index 0000000..d9f494e --- /dev/null +++ b/test/chan_vese.jl @@ -0,0 +1,13 @@ +@testset "chan_vese" begin + @info "Test: Chan Vese Segmentation" + + @testset "Gray Image Chan-Vese Segmentation Reference Test" begin + img_gray = imresize(testimage("cameraman"), (64, 64)) + ref = load("references/Chan_Vese_Gray.png") + ref = ref .> 0 + out = chan_vese(img_gray, μ=0.1, tol=1e-2, max_iter=200, normalize=true) + @test eltype(out) == Bool + @test sum(out) == sum(ref) + @test out == ref + end +end \ No newline at end of file diff --git a/test/references/Chan_Vese_Gray.png b/test/references/Chan_Vese_Gray.png new file mode 100644 index 0000000..5e84d54 Binary files /dev/null and b/test/references/Chan_Vese_Gray.png differ diff --git a/test/runtests.jl b/test/runtests.jl index 47c869f..b660757 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,4 +1,5 @@ using ImageSegmentation, ImageCore, ImageFiltering, Test, SimpleWeightedGraphs, LightGraphs, StaticArrays, DataStructures +using FileIO, ImageTransformations, TestImages using RegionTrees: isleaf, Cell, split! using MetaGraphs: MetaGraph, clear_props!, get_prop, has_prop, set_prop!, props, vertices @@ -15,3 +16,4 @@ include("watershed.jl") include("region_merging.jl") include("meanshift.jl") include("merge_segments.jl") +include("chan_vese.jl") \ No newline at end of file