From 96ceac37c67c56b1fcf3909e4bf1d178a5622608 Mon Sep 17 00:00:00 2001
From: JKay <2439951158@qq.com>
Date: Thu, 16 Sep 2021 14:04:32 +0800
Subject: [PATCH 1/6] add the Chan-Vese Segmentation for Gray

---
 src/chan_vese.jl | 252 +++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 252 insertions(+)
 create mode 100644 src/chan_vese.jl

diff --git a/src/chan_vese.jl b/src/chan_vese.jl
new file mode 100644
index 0000000..108e98f
--- /dev/null
+++ b/src/chan_vese.jl
@@ -0,0 +1,252 @@
+using MosaicViews
+using LazyArrays
+using BenchmarkTools
+using Images, TestImages
+using ImageBase.ImageCore: GenericGrayImage, GenericImage
+
+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𝚽ⁱ)
+    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
+    if ∫H𝚽 != 0
+        c₁ = ∫u₀H𝚽 / ∫H𝚽
+    end
+    if ∫H𝚽ⁱ != 0
+        c₂ = ∫u₀H𝚽ⁱ / ∫H𝚽ⁱ
+    end
+
+    return c₁, c₂
+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
+    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₂
+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)
+    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))
+    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
+
+    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)
+        𝚽ⱼ₊ = 0
+
+        for y in 1:n-1
+            𝚽ⱼ₊ = 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.)
+                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")
+
+μ=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)
\ No newline at end of file

From 982cfc4a25d8416eb692b003b254d03153d40905 Mon Sep 17 00:00:00 2001
From: JKay <2439951158@qq.com>
Date: Tue, 21 Sep 2021 02:24:55 +0800
Subject: [PATCH 2/6] modify Chan-Vese Segmentation, add docs and simple test

---
 Project.toml                       |   5 +-
 src/ImageSegmentation.jl           |   3 +
 src/chan_vese.jl                   | 436 +++++++++++++++--------------
 test/chan_vese.jl                  |  13 +
 test/references/Chan_Vese_Gray.png | Bin 0 -> 466 bytes
 test/runtests.jl                   |  20 +-
 6 files changed, 252 insertions(+), 225 deletions(-)
 create mode 100644 test/chan_vese.jl
 create mode 100644 test/references/Chan_Vese_Gray.png

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
index 108e98f..e7a3175 100644
--- a/src/chan_vese.jl
+++ b/src/chan_vese.jl
@@ -1,85 +1,108 @@
-using MosaicViews
-using LazyArrays
-using BenchmarkTools
-using Images, TestImages
-using ImageBase.ImageCore: GenericGrayImage, GenericImage
-
-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𝚽ⁱ)
-    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
-    if ∫H𝚽 != 0
-        c₁ = ∫u₀H𝚽 / ∫H𝚽
-    end
-    if ∫H𝚽ⁱ != 0
-        c₂ = ∫u₀H𝚽ⁱ / ∫H𝚽ⁱ
-    end
+"""
+    chan_vese(img; μ, λ₁, λ₂, tol, max_iter, Δt, reinitial_flag)
 
-    return c₁, c₂
-end
+Segments image `img` by evolving a level set. An active contour model 
+which can be used to segment objects without clearly defined boundaries.
 
-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𝚽)
+# output
+Return a `BitMatrix`.
 
-    if ndims(img) == 2
-        return @. -λ₁ * (img - c₁)^2 + λ₂ * (img - c₂)^2
-    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
+# Details
 
-    return c₁, c₂
-end
+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.
 
-function δₕ(x::AbstractArray{T,N}, h::Float64=1.0) where {T<:Real, N}
-    return @~ @. h / (h^2 + x^2)
-end
+# Options
 
-function initial_level_set(shape::Tuple)
-    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
+The function argument is described in detail below. 
+
+Denote the edge set curve with 𝐶 in the following part.
+
+## `μ::Float64`
+
+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
+
+## `λ₁::Float64`, `λ₂::Float64`
+
+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::Float64`
+
+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::Int64`
+
+The argument `max_iter` controls the maximum of iteration number.
+
+Default: 500
 
+## `Δt::Float64`
+
+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
+
+## reinitial_flag::Bool
+
+The arguement `reinitial_flag` controls whether to reinitialize the
+level set in each step.
+
+Default: false
+
+# Examples
+
+```julia
+using TestImages
+using ImageSegmentation
+
+img = testimage("cameraman")
+
+cv_result = chan_vese(img, μ=0.25, λ₁=1.0, λ₂=1.0, tol=1e-3, max_iter=200, Δt=0.5, reinitial_flag=false)
+```
+
+# 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;
-                    μ::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}
+                   μ::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)
+    # Signs used in the codes and comments mainly follow paper[3] in the References.
     img = float64.(channelview(img))
     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)
 
@@ -87,166 +110,149 @@ function chan_vese(img::GenericGrayImage;
         img .= img ./ maximum(img)
     end
 
-    diff = 0
-    H𝚽 = similar(𝚽ⁿ)
-    u₀H𝚽 = similar(img)
-    ∫u₀ = sum(img)
-    𝚽ᵢ₊ᶜ = zeros(m, 1)
+    # 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
+    𝚽ⁿ = initial_level_set(size(img))
 
+    # Preallocation and initializtion
+    H𝚽 = trues(size(img)...)
+    𝚽ⁿ⁺¹ = similar(𝚽ⁿ)
 
+    # The upper bounds of 𝚽ⁿ's coordinates is `m` and `n`
+    s, t = first(CartesianIndices(𝚽ⁿ))[1], first(CartesianIndices(𝚽ⁿ))[2]
+    m, n = last(CartesianIndices(𝚽ⁿ))[1], last(CartesianIndices(𝚽ⁿ))[2]
+    
     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𝚽
+        ϵ = 1e-8
+        diff = 0
 
-        if ∫H𝚽 != 0
-            c₁ = ∫u₀H𝚽 ./ ∫H𝚽
-        end
-        if ∫H𝚽ⁱ != 0
-            c₂ = ∫u₀H𝚽ⁱ ./ ∫H𝚽ⁱ
+        # 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 𝚽ⁿ
+        for idx in CartesianIndices(𝚽ⁿ) # Denote idx = (x, y)
+            # i₊ ≔ i₊(x, y), denotes 𝚽ⁿ(x, y + 1)'s CartesianIndex
+            # j₊ ≔ j₊(x, y), denotes 𝚽ⁿ(x + 1, y)'s CartesianIndex
+            # i₋ ≔ i₋(x, y), denotes 𝚽ⁿ(x, y - 1)'s CartesianIndex
+            # j₋ ≔ j₋(x, y), denotes 𝚽ⁿ(x - 1, y)'s CartesianIndex
+            # Taking notice that if 𝚽ⁿ(x, y) is the boundary of 𝚽ⁿ, than 𝚽ⁿ(x ± 1, y), 𝚽ⁿ(x, y ± 1) might be out of bound.
+            # So the pixel values of these outbounded terms are equal to 𝚽ⁿ(x, y)
+            i₊ = idx[2] != n ? idx + CartesianIndex(0, 1) : idx
+            j₊ = idx[1] != m ? idx + CartesianIndex(1, 0) : idx
+            i₋ = idx[2] != t ? idx - CartesianIndex(0, 1) : idx
+            j₋ = idx[1] != s ? idx - CartesianIndex(1, 0) : idx
+
+            𝚽₀  = 𝚽ⁿ[idx] # 𝚽ⁿ(x, y)
+            u₀ = img[idx] # u₀(x, y)
+            𝚽ᵢ₊ = 𝚽ⁿ[i₊] # 𝚽ⁿ(x, y + 1)
+            𝚽ⱼ₊ = 𝚽ⁿ[j₊] # 𝚽ⁿ(x + 1, y)
+            𝚽ᵢ₋ = 𝚽ⁿ[i₋] # 𝚽ⁿ(x, y - 1)
+            𝚽ⱼ₋ = 𝚽ⁿ[j₋] # 𝚽ⁿ(x - 1, y)
+
+            # Solve the PDE of equation 9 in paper[3]
+            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 = 𝚽ᵢ₊ * C₁ + 𝚽ᵢ₋ * C₂ + 𝚽ⱼ₊ * C₃ + 𝚽ⱼ₋ * 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 * δₕ * (C₁ + C₂ + C₃ + C₄))
+            diff += (𝚽 - 𝚽₀)^2
         end
 
-        ind = CartesianIndices(reshape(collect(1 : 9), 3, 3)) .- CartesianIndex(2, 2)
-        𝚽ⱼ₊ = 0
-
-        for y in 1:n-1
-            𝚽ⱼ₊ = 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.)
-                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
+        del = sqrt(diff / area)
 
-        𝚽ᵢ₊ = 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  
+        if reinitial_flag
+            # Reinitialize 𝚽 to be the signed distance function to its zero level set
+            reinitialize(𝚽ⁿ⁺¹, 𝚽ⁿ, Δt, h)
+        else
+            𝚽ⁿ .= 𝚽ⁿ⁺¹
         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
+    return 𝚽ⁿ .> 0
 end
 
-img_gray = testimage("cameraman")
+function initial_level_set(shape::Tuple)
+    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
 
-μ=0.25
-λ₁=1.0
-λ₂=1.0
-tol=1e-3
-max_iter=200
-Δt=0.5
+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
+    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𝚽ⁱ)
 
-𝚽, iter_num = chan_vese(img_gray, μ=0.25, λ₁=1.0, λ₂=1.0, tol=1e-3, max_iter=200, Δt=0.5, reinitial_flag=false)
+    return c₁, c₂
+end
 
-@btime chan_vese(img_gray, μ=0.25, λ₁=1.0, λ₂=1.0, tol=1e-3, max_iter=200, Δt=0.5, reinitial_flag=false);
+function calculate_reinitial(𝚽::AbstractArray{T, M}, 𝚿::AbstractArray{T, M}, Δt::Float64, h::Float64) where {T<:Real, M}
+    ϵ = 1e-8
+
+    s, t = first(CartesianIndices(𝚽))[1], first(CartesianIndices(𝚽))[2]
+    m, n = last(CartesianIndices(𝚽))[1], last(CartesianIndices(𝚽))[2]
+
+    for idx in CartesianIndices(𝚽)
+        i₊ = idx[2] != n ? idx + CartesianIndex(0, 1) : idx
+        j₊ = idx[1] != m ? idx + CartesianIndex(1, 0) : idx
+        i₋ = idx[2] != t ? idx - CartesianIndex(0, 1) : idx
+        j₋ = idx[1] != s ? idx - CartesianIndex(1, 0) : idx
+        𝚽₀  = 𝚽[idx]               # 𝚽(i, j)
+        𝚽ᵢ₊ = 𝚽[i₊]                # 𝚽(i + 1, j)
+        𝚽ⱼ₊ = 𝚽[j₊]                # 𝚽(i, j + 1)
+        𝚽ᵢ₋ = 𝚽[i₋]                # 𝚽(i - 1, j)
+        𝚽ⱼ₋ = 𝚽[j₋]                # 𝚽(i, j - 1)
+
+        a = (𝚽₀ - 𝚽ᵢ₋) / h
+        b = (𝚽ᵢ₊ - 𝚽₀) / h
+        c = (𝚽₀ - 𝚽ⱼ₋) / h
+        d = (𝚽ⱼ₊ - 𝚽₀) / h
+
+        a⁺ = max(a, 0)
+        a⁻ = min(a, 0)
+        b⁺ = max(b, 0)
+        b⁻ = min(b, 0)
+        c⁺ = max(c, 0)
+        c⁻ = min(c, 0)
+        d⁺ = max(d, 0)
+        d⁻ = min(d, 0)
+
+        G = 0
+        if 𝚽₀ > 0
+            G += sqrt(max(a⁺^2, b⁻^2) + max(c⁺^2, d⁻^2)) - 1
+        elseif 𝚽₀ < 0
+            G += sqrt(max(a⁻^2, b⁺^2) + max(c⁻^2, d⁺^2)) - 1
+        end
+        sign𝚽 = 𝚽₀ / sqrt(𝚽₀^2 + ϵ)
+        𝚿[idx] = 𝚽₀ - Δt * sign𝚽 * G
+    end
 
-segmentation = 𝚽 .> 0
-print(iter_num)
-𝚽 .= 𝚽 .- minimum(𝚽)
+    return 𝚿
+end
 
-colorview(Gray, segmentation)
\ No newline at end of file
+function reinitialize(𝚽::AbstractArray{T, M}, 𝚿::AbstractArray{T, M}, Δt::Float64, h::Float64, max_reiter::Int64=5) where {T<:Real, M}
+    iter = 0
+    while iter < max_reiter
+        𝚽 .= calculate_reinitial(𝚽, 𝚿, Δt, h)
+        iter += 1
+    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..eb0ab20
--- /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, λ₁=1.0, λ₂=1.0, tol=1e-2, max_iter=200, Δt=0.5, reinitial_flag=false)
+        @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 0000000000000000000000000000000000000000..5e84d54b3beb4cdc6e3d00c0c8b98eb68413046d
GIT binary patch
literal 466
zcmeAS@N?(olHy`uVBq!ia0vp^4j{|{BpCXc^q3eJ7-xFAIEF-UE<MOOX|@9ci^13b
z@^>3q&sg@Yp73aCg=SpWbT}ok`p0^SmG&QFd#)($>uS7p#N04o_763-EHT%fqYl1K
zXOmRquBxPo_D<ZfA~qn%P)NK~#FV>Ny)aCDVc16%zgId}bl4s($Xk?VD6o51cjQrr
znk&sSLV1oXlRLV^P-u456t|-uYhKNk=nSnDR^GT$KS{+dOQS6*T-i`$chw{%x1?}?
zr=uQkUMaa=km!t!b>Y-vd$d^3snauA!S0pPlo=s5-M+D{#}<9mF)J06THWQoamDO_
zjVmfdx_5Dku6AoZ<}SFpO04_V(M27uv8`SGw~nq++qh!&g&U!BuBwQY3Yzdq@jvp9
zWfk4swfjeKxZ}}vkM#asaSeF1;E|eGshCM`bJy&MtDYAuvhN&O*4NZExoh=}D~=a_
zIOe=k&?pQm6r63vAr&qvzHw#qg2LcJvDvp+Ro7{qOA2Lo?3{fpa?LA+j6X_Yrru|p
zj_D_bZoZ;YC%$^t(zK78Bd)4siP?9D`YwWUg!cbtUl-eE+k0{2R8VAjy85}Sb4q9e
E0I1W%p8x;=

literal 0
HcmV?d00001

diff --git a/test/runtests.jl b/test/runtests.jl
index 47c869f..4d64534 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
 
@@ -6,12 +7,13 @@ using Documenter
 Base.VERSION >= v"1.6" && doctest(ImageSegmentation, manual = false)
 @test isempty(detect_ambiguities(ImageSegmentation))
 
-include("core.jl")
-include("region_growing.jl")
-include("felzenszwalb.jl")
-include("fast_scanning.jl")
-include("flood_fill.jl")
-include("watershed.jl")
-include("region_merging.jl")
-include("meanshift.jl")
-include("merge_segments.jl")
+# include("core.jl")
+# include("region_growing.jl")
+# include("felzenszwalb.jl")
+# include("fast_scanning.jl")
+# include("flood_fill.jl")
+# include("watershed.jl")
+# include("region_merging.jl")
+# include("meanshift.jl")
+# include("merge_segments.jl")
+include("chan_vese.jl")
\ No newline at end of file

From 9996e09cf77e1b7ebd46165f00f940e4a7790def Mon Sep 17 00:00:00 2001
From: JKay <2439951158@qq.com>
Date: Tue, 21 Sep 2021 02:29:39 +0800
Subject: [PATCH 3/6] modify runtest.jl

---
 test/runtests.jl | 18 +++++++++---------
 1 file changed, 9 insertions(+), 9 deletions(-)

diff --git a/test/runtests.jl b/test/runtests.jl
index 4d64534..b660757 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -7,13 +7,13 @@ using Documenter
 Base.VERSION >= v"1.6" && doctest(ImageSegmentation, manual = false)
 @test isempty(detect_ambiguities(ImageSegmentation))
 
-# include("core.jl")
-# include("region_growing.jl")
-# include("felzenszwalb.jl")
-# include("fast_scanning.jl")
-# include("flood_fill.jl")
-# include("watershed.jl")
-# include("region_merging.jl")
-# include("meanshift.jl")
-# include("merge_segments.jl")
+include("core.jl")
+include("region_growing.jl")
+include("felzenszwalb.jl")
+include("fast_scanning.jl")
+include("flood_fill.jl")
+include("watershed.jl")
+include("region_merging.jl")
+include("meanshift.jl")
+include("merge_segments.jl")
 include("chan_vese.jl")
\ No newline at end of file

From 7e3828ddbfea4cb880c0ea4d89eafd63a136a185 Mon Sep 17 00:00:00 2001
From: JKay <2439951158@qq.com>
Date: Thu, 23 Sep 2021 23:26:53 +0800
Subject: [PATCH 4/6] generic implementation for N dims

---
 src/chan_vese.jl  | 106 ++++++++++++++++++----------------------------
 test/chan_vese.jl |   2 +-
 2 files changed, 42 insertions(+), 66 deletions(-)

diff --git a/src/chan_vese.jl b/src/chan_vese.jl
index e7a3175..9b115b7 100644
--- a/src/chan_vese.jl
+++ b/src/chan_vese.jl
@@ -1,5 +1,5 @@
 """
-    chan_vese(img; μ, λ₁, λ₂, tol, max_iter, Δt, reinitial_flag)
+    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.
@@ -50,7 +50,7 @@ iterations is below `tol`, then the solution will be assumed to be reached.
 
 Default: 1e-3
 
-## `max_iter::Int64`
+## `max_iter::Int`
 
 The argument `max_iter` controls the maximum of iteration number.
 
@@ -80,7 +80,7 @@ using ImageSegmentation
 
 img = testimage("cameraman")
 
-cv_result = chan_vese(img, μ=0.25, λ₁=1.0, λ₂=1.0, tol=1e-3, max_iter=200, Δt=0.5, reinitial_flag=false)
+cv_result = chan_vese(img, max_iter=200)
 ```
 
 # References
@@ -96,11 +96,12 @@ function chan_vese(img::GenericGrayImage;
                    λ₁::Float64=1.0,
                    λ₂::Float64=1.0,
                    tol::Float64=1e-3,
-                   max_iter::Int64=500,
+                   max_iter::Int=500,
                    Δt::Float64=0.5,
                    reinitial_flag::Bool=false)
     # Signs used in the codes and comments mainly follow paper[3] in the References.
     img = float64.(channelview(img))
+    N = ndims(img)
     iter = 0
     h = 1.0
     del = tol + 1
@@ -111,6 +112,7 @@ function chan_vese(img::GenericGrayImage;
     end
 
     # Precalculation of some constants which helps simplify some integration   
+    # area = length(img) # area = ∫H𝚽 + ∫H𝚽ⁱ
     area = length(img) # area = ∫H𝚽 + ∫H𝚽ⁱ
     ∫u₀ = sum(img)     # ∫u₀ = ∫u₀H𝚽 + ∫u₀H𝚽ⁱ
 
@@ -121,9 +123,9 @@ function chan_vese(img::GenericGrayImage;
     H𝚽 = trues(size(img)...)
     𝚽ⁿ⁺¹ = similar(𝚽ⁿ)
 
-    # The upper bounds of 𝚽ⁿ's coordinates is `m` and `n`
-    s, t = first(CartesianIndices(𝚽ⁿ))[1], first(CartesianIndices(𝚽ⁿ))[2]
-    m, n = last(CartesianIndices(𝚽ⁿ))[1], last(CartesianIndices(𝚽ⁿ))[2]
+    Δ = 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
@@ -134,36 +136,23 @@ function chan_vese(img::GenericGrayImage;
         c₁, c₂ = calculate_averages(img, H𝚽, area, ∫u₀) # Compute c₁(𝚽ⁿ), c₂(𝚽ⁿ)
 
         # Calculate the variation of level set 𝚽ⁿ
-        for idx in CartesianIndices(𝚽ⁿ) # Denote idx = (x, y)
-            # i₊ ≔ i₊(x, y), denotes 𝚽ⁿ(x, y + 1)'s CartesianIndex
-            # j₊ ≔ j₊(x, y), denotes 𝚽ⁿ(x + 1, y)'s CartesianIndex
-            # i₋ ≔ i₋(x, y), denotes 𝚽ⁿ(x, y - 1)'s CartesianIndex
-            # j₋ ≔ j₋(x, y), denotes 𝚽ⁿ(x - 1, y)'s CartesianIndex
-            # Taking notice that if 𝚽ⁿ(x, y) is the boundary of 𝚽ⁿ, than 𝚽ⁿ(x ± 1, y), 𝚽ⁿ(x, y ± 1) might be out of bound.
-            # So the pixel values of these outbounded terms are equal to 𝚽ⁿ(x, y)
-            i₊ = idx[2] != n ? idx + CartesianIndex(0, 1) : idx
-            j₊ = idx[1] != m ? idx + CartesianIndex(1, 0) : idx
-            i₋ = idx[2] != t ? idx - CartesianIndex(0, 1) : idx
-            j₋ = idx[1] != s ? idx - CartesianIndex(1, 0) : idx
-
+        @inbounds @simd for idx in CartesianIndices(𝚽ⁿ)
             𝚽₀  = 𝚽ⁿ[idx] # 𝚽ⁿ(x, y)
-            u₀ = img[idx] # u₀(x, y)
-            𝚽ᵢ₊ = 𝚽ⁿ[i₊] # 𝚽ⁿ(x, y + 1)
-            𝚽ⱼ₊ = 𝚽ⁿ[j₊] # 𝚽ⁿ(x + 1, y)
-            𝚽ᵢ₋ = 𝚽ⁿ[i₋] # 𝚽ⁿ(x, y - 1)
-            𝚽ⱼ₋ = 𝚽ⁿ[j₋] # 𝚽ⁿ(x - 1, y)
+            u₀ = img[idx]  # u₀(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], Δ₊)
+            𝚽₋ = broadcast(i -> 𝚽ⁿ[i], Δ₋)
 
             # Solve the PDE of equation 9 in paper[3]
-            C₁ = 1. / sqrt(ϵ + (𝚽ᵢ₊ - 𝚽₀)^2 + (𝚽ⱼ₊ - 𝚽ⱼ₋)^2 / 4.)
-            C₂ = 1. / sqrt(ϵ + (𝚽₀ - 𝚽ᵢ₋)^2 + (𝚽ⱼ₊ - 𝚽ⱼ₋)^2 / 4.)
-            C₃ = 1. / sqrt(ϵ + (𝚽ᵢ₊ - 𝚽ᵢ₋)^2 / 4. + (𝚽ⱼ₊ - 𝚽₀)^2)
-            C₄ = 1. / sqrt(ϵ + (𝚽ᵢ₊ - 𝚽ᵢ₋)^2 / 4. + (𝚽₀ - 𝚽ⱼ₋)^2)
+            C₊ = ntuple(d -> 1. / sqrt(ϵ + (𝚽₊[d] - 𝚽₀)^2 + (𝚽₊[d % N + 1] - 𝚽₋[d % N + 1])^2 / 4.), N)
+            C₋ = ntuple(d -> 1. / sqrt(ϵ + (𝚽₋[d] - 𝚽₀)^2 + (𝚽₊[d % N + 1] - 𝚽₋[d % N + 1])^2 / 4.), N)
 
-            K = 𝚽ᵢ₊ * C₁ + 𝚽ᵢ₋ * C₂ + 𝚽ⱼ₊ * C₃ + 𝚽ⱼ₋ * C₄
+            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 * δₕ * (C₁ + C₂ + C₃ + C₄))
+            𝚽ⁿ⁺¹[idx] = 𝚽 = (𝚽₀ + Δt * δₕ * (μ * K + difference_from_average)) / (1. + μ * Δt * δₕ * (sum(C₊) + sum(C₋)))
             diff += (𝚽 - 𝚽₀)^2
         end
 
@@ -171,7 +160,7 @@ function chan_vese(img::GenericGrayImage;
 
         if reinitial_flag
             # Reinitialize 𝚽 to be the signed distance function to its zero level set
-            reinitialize(𝚽ⁿ⁺¹, 𝚽ⁿ, Δt, h)
+            reinitialize!(𝚽ⁿ⁺¹, 𝚽ⁿ, Δt, h)
         else
             𝚽ⁿ .= 𝚽ⁿ⁺¹
         end
@@ -191,7 +180,7 @@ 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
-    for i in eachindex(img)
+    @inbounds for i in eachindex(img)
         if H𝚽[i]
             ∫u₀H𝚽 += img[i]
             ∫H𝚽 += 1
@@ -207,40 +196,29 @@ 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, Δ₋)
 
-    s, t = first(CartesianIndices(𝚽))[1], first(CartesianIndices(𝚽))[2]
-    m, n = last(CartesianIndices(𝚽))[1], last(CartesianIndices(𝚽))[2]
-
-    for idx in CartesianIndices(𝚽)
-        i₊ = idx[2] != n ? idx + CartesianIndex(0, 1) : idx
-        j₊ = idx[1] != m ? idx + CartesianIndex(1, 0) : idx
-        i₋ = idx[2] != t ? idx - CartesianIndex(0, 1) : idx
-        j₋ = idx[1] != s ? idx - CartesianIndex(1, 0) : idx
-        𝚽₀  = 𝚽[idx]               # 𝚽(i, j)
-        𝚽ᵢ₊ = 𝚽[i₊]                # 𝚽(i + 1, j)
-        𝚽ⱼ₊ = 𝚽[j₊]                # 𝚽(i, j + 1)
-        𝚽ᵢ₋ = 𝚽[i₋]                # 𝚽(i - 1, j)
-        𝚽ⱼ₋ = 𝚽[j₋]                # 𝚽(i, j - 1)
-
-        a = (𝚽₀ - 𝚽ᵢ₋) / h
-        b = (𝚽ᵢ₊ - 𝚽₀) / h
-        c = (𝚽₀ - 𝚽ⱼ₋) / h
-        d = (𝚽ⱼ₊ - 𝚽₀) / h
-
-        a⁺ = max(a, 0)
-        a⁻ = min(a, 0)
-        b⁺ = max(b, 0)
-        b⁻ = min(b, 0)
-        c⁺ = max(c, 0)
-        c⁻ = min(c, 0)
-        d⁺ = max(d, 0)
-        d⁻ = min(d, 0)
+        maxΔ𝚽₊ = max.(Δ𝚽₊, 0)
+        minΔ𝚽₊ = min.(Δ𝚽₊, 0)
+        maxΔ𝚽₋ = max.(Δ𝚽₋, 0)
+        minΔ𝚽₋ = min.(Δ𝚽₋, 0)
 
         G = 0
         if 𝚽₀ > 0
-            G += sqrt(max(a⁺^2, b⁻^2) + max(c⁺^2, d⁻^2)) - 1
+            G += sqrt(sum(max.(minΔ𝚽₊.^2, maxΔ𝚽₋.^2))) - 1
         elseif 𝚽₀ < 0
-            G += sqrt(max(a⁻^2, b⁺^2) + max(c⁻^2, d⁺^2)) - 1
+            G += sqrt(sum(max.(maxΔ𝚽₊.^2, minΔ𝚽₋.^2))) - 1
         end
         sign𝚽 = 𝚽₀ / sqrt(𝚽₀^2 + ϵ)
         𝚿[idx] = 𝚽₀ - Δt * sign𝚽 * G
@@ -249,10 +227,8 @@ function calculate_reinitial(𝚽::AbstractArray{T, M}, 𝚿::AbstractArray{T, M
     return 𝚿
 end
 
-function reinitialize(𝚽::AbstractArray{T, M}, 𝚿::AbstractArray{T, M}, Δt::Float64, h::Float64, max_reiter::Int64=5) where {T<:Real, M}
-    iter = 0
-    while iter < max_reiter
+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)
-        iter += 1
     end
 end
\ No newline at end of file
diff --git a/test/chan_vese.jl b/test/chan_vese.jl
index eb0ab20..23d93e1 100644
--- a/test/chan_vese.jl
+++ b/test/chan_vese.jl
@@ -5,7 +5,7 @@
         img_gray = imresize(testimage("cameraman"), (64, 64))
         ref = load("references/Chan_Vese_Gray.png")
         ref = ref .> 0
-        out = chan_vese(img_gray, μ=0.1, λ₁=1.0, λ₂=1.0, tol=1e-2, max_iter=200, Δt=0.5, reinitial_flag=false)
+        out = chan_vese(img_gray, μ=0.1, tol=1e-2, max_iter=200)
         @test eltype(out) == Bool
         @test sum(out) == sum(ref)
         @test out == ref

From 9e07ebcb9e4195e9548ffab9824e759f57aec4ee Mon Sep 17 00:00:00 2001
From: JKay <2439951158@qq.com>
Date: Fri, 24 Sep 2021 13:07:52 +0800
Subject: [PATCH 5/6] make the chan_vese function available for 3D image

---
 src/chan_vese.jl | 21 ++++++++++++++-------
 1 file changed, 14 insertions(+), 7 deletions(-)

diff --git a/src/chan_vese.jl b/src/chan_vese.jl
index 9b115b7..ab17770 100644
--- a/src/chan_vese.jl
+++ b/src/chan_vese.jl
@@ -139,14 +139,14 @@ function chan_vese(img::GenericGrayImage;
         @inbounds @simd for idx in CartesianIndices(𝚽ⁿ)
             𝚽₀  = 𝚽ⁿ[idx] # 𝚽ⁿ(x, y)
             u₀ = img[idx]  # u₀(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], Δ₊)
-            𝚽₋ = broadcast(i -> 𝚽ⁿ[i], Δ₋)
+            𝚽₊ = 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]
-            C₊ = ntuple(d -> 1. / sqrt(ϵ + (𝚽₊[d] - 𝚽₀)^2 + (𝚽₊[d % N + 1] - 𝚽₋[d % N + 1])^2 / 4.), N)
-            C₋ = ntuple(d -> 1. / sqrt(ϵ + (𝚽₋[d] - 𝚽₀)^2 + (𝚽₊[d % N + 1] - 𝚽₋[d % N + 1])^2 / 4.), N)
+            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
@@ -171,12 +171,19 @@ function chan_vese(img::GenericGrayImage;
     return 𝚽ⁿ .> 0
 end
 
-function initial_level_set(shape::Tuple)
+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

From 6395478ef69033202f2223835d3cfb3fd7ea0748 Mon Sep 17 00:00:00 2001
From: JKay <2439951158@qq.com>
Date: Tue, 28 Sep 2021 21:22:40 +0800
Subject: [PATCH 6/6] make normalize and init level set be optimal, and some
 small modifications

---
 src/chan_vese.jl  | 69 ++++++++++++++++++++++++++++-------------------
 test/chan_vese.jl |  2 +-
 2 files changed, 42 insertions(+), 29 deletions(-)

diff --git a/src/chan_vese.jl b/src/chan_vese.jl
index ab17770..a3ec7c0 100644
--- a/src/chan_vese.jl
+++ b/src/chan_vese.jl
@@ -21,7 +21,7 @@ The function argument is described in detail below.
 
 Denote the edge set curve with 𝐶 in the following part.
 
-## `μ::Float64`
+## `μ::Real`
 
 The argument `μ` is a weight controlling the penalty on the total length
 of the curve 𝐶;
@@ -31,7 +31,7 @@ can prevent 𝐶 from being a complex curve.
 
 Default: 0.25
 
-## `λ₁::Float64`, `λ₂::Float64`
+## `λ₁::Real`, `λ₂::Real`
 
 The argument `λ₁` and `λ₂` affect the desired uniformity inside 𝐶 and 
 outside 𝐶, respectively. 
@@ -42,7 +42,7 @@ quite uniform background and varying grayscale objects in the foreground.
 Default: λ₁ = 1.0
          λ₂ = 1.0
 
-## `tol::Float64`
+## `tol::Real`
 
 The argument `tol` controls the level set variation tolerance between 
 iteration. If the L2 norm difference between two level sets of adjacent
@@ -56,7 +56,7 @@ The argument `max_iter` controls the maximum of iteration number.
 
 Default: 500
 
-## `Δt::Float64`
+## `Δt::Real`
 
 The argument `Δt` is a multiplication factor applied at calculations 
 for each step, serves to accelerate the algorithm. Although larger `Δt`
@@ -65,13 +65,23 @@ the solution.
 
 Default: 0.5
 
-## reinitial_flag::Bool
+## normalize::Bool
 
-The arguement `reinitial_flag` controls whether to reinitialize the
-level set in each step.
+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
@@ -92,32 +102,34 @@ cv_result = chan_vese(img, max_iter=200)
 [3] The Chan-Vese Algorithm - Project Report, Rami Cohen, 2011 :arXiv:`1107.2782`
 """
 function chan_vese(img::GenericGrayImage;
-                   μ::Float64=0.25,
-                   λ₁::Float64=1.0,
-                   λ₂::Float64=1.0,
-                   tol::Float64=1e-3,
+                   μ::Real=0.25,
+                   λ₁::Real=1.0,
+                   λ₂::Real=1.0,
+                   tol::Real=1e-3,
                    max_iter::Int=500,
-                   Δt::Float64=0.5,
-                   reinitial_flag::Bool=false)
+                   Δ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
-    img .= img .- minimum(img)
-
-    if maximum(img) != 0
-        img .= img ./ maximum(img)
+    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𝚽ⁱ
-    area = length(img) # area = ∫H𝚽 + ∫H𝚽ⁱ
-    ∫u₀ = sum(img)     # ∫u₀ = ∫u₀H𝚽 + ∫u₀H𝚽ⁱ
+    # 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
-    𝚽ⁿ = initial_level_set(size(img))
+    𝚽ⁿ = init_level_set
 
     # Preallocation and initializtion
     H𝚽 = trues(size(img)...)
@@ -158,12 +170,13 @@ function chan_vese(img::GenericGrayImage;
 
         del = sqrt(diff / area)
 
-        if reinitial_flag
-            # Reinitialize 𝚽 to be the signed distance function to its zero level set
-            reinitialize!(𝚽ⁿ⁺¹, 𝚽ⁿ, Δt, h)
-        else
-            𝚽ⁿ .= 𝚽ⁿ⁺¹
-        end
+        # 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
diff --git a/test/chan_vese.jl b/test/chan_vese.jl
index 23d93e1..d9f494e 100644
--- a/test/chan_vese.jl
+++ b/test/chan_vese.jl
@@ -5,7 +5,7 @@
         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)
+        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