Initial (trivial) implementation of flood_fill (#75)

Also adds window-iterators to traverse different neighborhood patterns.

Author: timholy

Project.toml

@@ -36,8 +36,10 @@ julia = "1"
 
 [extras]
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
+ImageMagick = "6218d12a-5da1-5696-b52f-db25d2ecc6d1"
 Images = "916415d5-f1e6-5110-898d-aaa5f9f070e0"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Documenter", "Images", "Test"]
+test = ["Documenter", "FileIO", "ImageMagick", "Images", "Test"]

src/ImageSegmentation.jl

@@ -13,11 +13,13 @@ 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("clustering.jl")
 include("merge_segments.jl")
+include("deprecations.jl")
 
 export
     #accessor methods
@@ -31,6 +33,9 @@ export
     unseeded_region_growing,
     felzenszwalb,
     fast_scanning,
+    fast_scanning!,
+    flood_fill,
+    flood_fill!,
     watershed,
     hmin_transform,
     region_adjacency_graph,
@@ -43,7 +48,7 @@ export
     kmeans,
     fuzzy_cmeans,
     merge_segments,
-    
+
     # types
     SegmentedImage,
     ImageEdge

src/compat.jl

@@ -5,3 +5,7 @@
 else
     const _oneunit = Base.oneunit
 end
+
+# Once Base has colon defined here we can replace this
+_colon(I::CartesianIndex{N}, J::CartesianIndex{N}) where N =
+    CartesianIndices(map((i,j) -> i:j, Tuple(I), Tuple(J)))

src/core.jl

@@ -1,7 +1,8 @@
-accum_type(::Type{T}) where {T<:Integer} = Int
-accum_type(::Type{Float32})    = Float32
-accum_type(::Type{T}) where {T<:Real} = Float64
-accum_type(::Type{C}) where {C<:Colorant} = base_colorant_type(C){accum_type(eltype(C))}
+accum_type(::Type{T}) where {T<:Integer}    = Int
+accum_type(::Type{Float32})                 = Float32
+accum_type(::Type{T}) where {T<:Real}       = Float64
+accum_type(::Type{T}) where {T<:FixedPoint} = floattype(T)
+accum_type(::Type{C}) where {C<:Colorant}   = base_colorant_type(C){accum_type(eltype(C))}
 
 """
 `SegmentedImage` type contains the index-label mapping, assigned labels,
@@ -266,10 +267,68 @@ function prune_segments(s::SegmentedImage, is_rem::Function, diff_fn::Function)
 
 end
 
-# Once Base has colon defined here we can replace this
-_colon(I::CartesianIndex{N}, J::CartesianIndex{N}) where N =
-    CartesianIndices(map((i,j) -> i:j, Tuple(I), Tuple(J)))
 
+"""
+    box_iterator(window)
+
+Return a function that constructs a box-shaped iterable region.
+
+# Examples
+```jldoctest; setup=:(using ImageSegmentation), filter=r"#\\d+"
+julia> fiter = ImageSegmentation.box_iterator((3, 3))
+#17 (generic function with 1 method)
+
+julia> center = CartesianIndex(17, 24)
+CartesianIndex(17, 24)
+
+julia> fiter(center)
+3×3 CartesianIndices{2, Tuple{UnitRange{$Int}, UnitRange{$Int}}}:
+ CartesianIndex(16, 23)  CartesianIndex(16, 24)  CartesianIndex(16, 25)
+ CartesianIndex(17, 23)  CartesianIndex(17, 24)  CartesianIndex(17, 25)
+ CartesianIndex(18, 23)  CartesianIndex(18, 24)  CartesianIndex(18, 25)
+```
+"""
+function box_iterator(window::Dims{N}) where N
+    for dim in window
+        dim > 0 || error("Dimensions of the window must be positive")
+        isodd(dim) || error("Dimensions of the window must be odd")
+    end
+    halfwindow = CartesianIndex(map(x -> x ÷ 2, window))
+    return function(center::CartesianIndex{N})
+        _colon(center-halfwindow, center+halfwindow)
+    end
+end
+
+"""
+    diamond_iterator(window)
+
+Return a function that constructs a diamond-shaped iterable region.
+
+# Examples
+```jldoctest; setup=:(using ImageSegmentation), filter=r"#\\d+"
+julia> fiter = ImageSegmentation.diamond_iterator((3, 3))
+#18 (generic function with 1 method)
+
+julia> center = CartesianIndex(17, 24)
+CartesianIndex(17, 24)
+
+julia> fiter(center)
+(CartesianIndex(18, 24), CartesianIndex(17, 25), CartesianIndex(16, 24), CartesianIndex(17, 23))
+```
+"""
+function diamond_iterator(window::Dims{N}) where N
+    for dim in window
+        dim > 0 || error("Dimensions of the window must be positive")
+        isodd(dim) || error("Dimensions of the window must be odd")
+    end
+    halfwindow = CartesianIndex(map(x -> x ÷ 2, window))
+    return function(center::CartesianIndex{N})
+        (ntuple(i -> center + CartesianIndex(ntuple(j -> i == j, N)), N)...,
+         ntuple(i -> center - CartesianIndex(ntuple(j -> i == j, N)), N)...)
+    end
+end
+
+window_neighbors(img::AbstractArray{T,N}) where {T,N} = ntuple(_ -> 3, N)
 
 """
     G, vertex2cartesian = region_adjacency_graph(img, weight_fn, R)

src/deprecations.jl

@@ -0,0 +1,11 @@
+@deprecate seeded_region_growing(
+    img::AbstractArray,
+    seeds::AbstractVector,
+    kernel_dim::Vector{Int},
+    diff_fn::Function = default_diff_fn) seeded_region_growing(img, seeds, (kernel_dim...,), diff_fn)
+
+@deprecate unseeded_region_growing(
+    img::AbstractArray,
+    threshold::Real,
+    kernel_dim::Vector{Int},
+    diff_fn::Function = default_diff_fn) unseeded_region_growing(img, threshold, (kernel_dim...,), diff_fn)

src/flood_fill.jl

@@ -0,0 +1,37 @@
+flood_fill(f, src::AbstractArray, idx::Union{Integer,CartesianIndex}, nbrhood_function=diamond_iterator(window_neighbors(src))) =
+    flood_fill!(f, falses(axes(src)) #=fill!(similar(src, Bool), false)=#, src, idx, nbrhood_function)
+
+function flood_fill(src::AbstractArray, idx::Union{Integer,CartesianIndex},
+                    nbrhood_function=diamond_iterator(window_neighbors(src)); thresh)
+    validx = src[idx]
+    validx = accum_type(typeof(validx))(validx)
+    return let validx=validx
+        flood_fill(val -> default_diff_fn(val, validx) < thresh, src, idx, nbrhood_function)
+    end
+end
+
+function flood_fill!(f, dest, src::AbstractArray, idx::Union{Int,CartesianIndex}, nbrhood_function=diamond_iterator(window_neighbors(src)))
+    R = CartesianIndices(src)
+    axes(dest) == R.indices || throw(DimensionMismatch("$(axes(dest)) do not match $(Tuple(R))"))
+    idx = R[idx]  # ensure cartesian indexing
+    f(src[idx]) || throw(ArgumentError("starting point fails to meet criterion"))
+    q = [idx]
+    _flood_fill!(f, dest, src, R, q, nbrhood_function)
+    return dest
+end
+flood_fill!(f, dest, src::AbstractArray, idx::Integer, nbrhood_function=diamond_iterator(window_neighbors(src))) =
+    flood_fill!(f, dest, src, Int(idx)::Int, nbrhood_function)
+
+# This is a trivial implementation (just to get something working), better would be a raster implementation
+function _flood_fill!(f::F, dest, src, R::CartesianIndices{N}, q, nbrhood_function::FN) where {F,N,FN}
+    while !isempty(q)
+        idx = pop!(q)
+        dest[idx] = true
+        @inbounds for j in nbrhood_function(idx)
+            j ∈ R || continue
+            if f(src[j]) && !dest[j]
+                push!(q, j)
+            end
+        end
+    end
+end

src/region_growing.jl

@@ -1,5 +1,6 @@
 
-default_diff_fn(c1::CT1,c2::CT2) where {CT1<:Union{Colorant,Real}, CT2<:Union{Colorant,Real}} = sqrt(sum(abs2,(c1)-accum_type(CT2)(c2)))
+default_diff_fn(c1::CT1,c2::CT2) where {CT1<:Union{Colorant,Real}, CT2<:Union{Colorant,Real}} = sqrt(_abs2((c1)-accum_type(CT2)(c2)))
+_abs2(c) = mapreducec(v->float(v)^2, +, 0, c)/length(c)
 
 """
     seg_img = seeded_region_growing(img, seeds, [kernel_dim], [diff_fn])
@@ -53,19 +54,12 @@ Albert Mehnert, Paul Jackaway (1997), "An improved seeded region growing algorit
 Pattern Recognition Letters 18 (1997), 1065-1071
 """
 function seeded_region_growing(img::AbstractArray{CT,N}, seeds::AbstractVector{<:PairOrTuple{CartesianIndex{N},Int}},
-    kernel_dim::Union{Vector{Int}, NTuple{N, Int}} = ntuple(i->3,N), diff_fn::Function = default_diff_fn) where {CT<:Union{Colorant,Real}, N}
-    length(kernel_dim) == N || error("Dimension count of image and kernel_dim do not match")
-    for dim in kernel_dim
-        dim > 0 || error("Dimensions of the kernel must be positive")
-        isodd(dim) || error("Dimensions of the kernel must be odd")
-    end
-    pt = CartesianIndex(ntuple(i->kernel_dim[i]÷2, N))
-    neighbourhood_gen(t) = c->_colon(c-t,c+t)
-    seeded_region_growing(img, seeds, neighbourhood_gen(pt), diff_fn)
+    kernel_dim::Dims{N} = ntuple(i->3,N), diff_fn::Function = default_diff_fn) where {CT<:Union{Colorant,Real}, N}
+    seeded_region_growing(img, seeds, box_iterator(kernel_dim), diff_fn)
 end
 
 function seeded_region_growing(img::AbstractArray{CT,N}, seeds::AbstractVector{<:PairOrTuple{CartesianIndex{N},Int}},
-    neighbourhood::NF, diff_fn::DF = default_diff_fn) where {CT<:Union{Colorant,Real}, N, NF, DF}
+    neighbourhood::NF, diff_fn::DF = default_diff_fn) where {CT<:Union{Colorant,Real}, N, NF<:Function, DF<:Function}
 
     # Check if labels are positive integers
     for seed in seeds
@@ -243,18 +237,11 @@ julia> labels_map(seg)
 
 """
 function unseeded_region_growing(img::AbstractArray{CT,N}, threshold::Real,
-    kernel_dim::Union{Vector{Int}, NTuple{N, Int}} = ntuple(i->3,N), diff_fn::Function = default_diff_fn) where {CT<:Colorant, N}
-    length(kernel_dim) == N || error("Dimension count of image and kernel_dim do not match")
-    for dim in kernel_dim
-        dim > 0 || error("Dimensions of the kernel must be positive")
-        isodd(dim) || error("Dimensions of the kernel must be odd")
-    end
-    pt = CartesianIndex(ntuple(i->kernel_dim[i]÷2, N))
-    neighbourhood_gen(t) = c->_colon(c-t,c+t)
-    unseeded_region_growing(img, threshold, neighbourhood_gen(pt), diff_fn)
+    kernel_dim::Dims{N} = ntuple(i->3,N), diff_fn::Function = default_diff_fn) where {CT<:Colorant, N}
+    unseeded_region_growing(img, threshold, box_iterator(kernel_dim), diff_fn)
 end
 
-function unseeded_region_growing(img::AbstractArray{CT,N}, threshold::Real, neighbourhood::Function, diff_fn = default_diff_fn) where {CT<:Colorant,N}
+function unseeded_region_growing(img::AbstractArray{CT,N}, threshold::Real, neighbourhood::Function, diff_fn::Function = default_diff_fn) where {CT<:Colorant,N}
     TM = meantype(CT)
 
     # Fast linear<->cartesian indexing lookup

test/flood_fill.jl

@@ -0,0 +1,65 @@
+using ImageSegmentation
+using ImageSegmentation.Colors
+using ImageSegmentation.FixedPointNumbers
+using FileIO
+using Statistics
+using Test
+
+@testset "flood_fill" begin
+    # 0d
+    a = reshape([true])
+    @test flood_fill(identity, a, CartesianIndex()) == a
+    @test_throws ArgumentError flood_fill(!, a, CartesianIndex())
+    # 1d
+    a = 1:7
+    @test flood_fill(==(2), a, CartesianIndex(2)) == (a .== 2)
+    @test_throws ArgumentError flood_fill(==(2), a, CartesianIndex(3))
+    @test flood_fill(x -> 1 < x < 4, a, CartesianIndex(2)) == [false, true, true, false, false, false, false]
+    @test flood_fill(isinteger, a, CartesianIndex(2)) == trues(7)
+    # 2d
+    ab = [true false false false;
+         true true false false;
+         true false false true;
+         true true true true]
+    an0f8 = N0f8.(ab)
+    agray = Gray.(an0f8)
+        for (f, a) in ((identity, ab), (==(1), an0f8), (==(1), agray))
+        for idx in CartesianIndices(a)
+            if f(a[idx])
+                @test flood_fill(f, a, idx) == a
+            else
+                @test_throws ArgumentError flood_fill(f, a, idx)
+            end
+        end
+    end
+    @test flood_fill(identity, ab, Int16(1)) == ab
+    # 3d
+    k = 10
+    a = falses(k, k, k)
+    idx = CartesianIndex(1,1,1)
+    incs = [CartesianIndex(1,0,0), CartesianIndex(0,1,0), CartesianIndex(0,0,1)]
+    a[idx] = true
+    while any(<(k), Tuple(idx))
+        d = rand(1:3)
+        idx += incs[d]
+        idx = min(idx, CartesianIndex(k,k,k))
+        a[idx] = true
+    end
+    for idx in eachindex(a)
+        if a[idx]
+            @test flood_fill(identity, a, idx) == a
+        else
+            @test_throws ArgumentError flood_fill(identity, a, idx)
+        end
+    end
+    # Colors
+    path = download("https://github.com/JuliaImages/juliaimages.github.io/raw/source/docs/src/pkgs/segmentation/assets/flower.jpg")
+    img = load(path)
+    seg = flood_fill(img, CartesianIndex(87,280); thresh=0.3)
+    @test 0.2*length(seg) <= sum(seg) <= 0.25*length(seg)
+    c = mean(img[seg])
+    # N0f8 makes for easier approximate testing
+    @test N0f8(red(c)) ≈ N0f8(0.855)
+    @test N0f8(green(c)) ≈ N0f8(0.161)
+    @test N0f8(blue(c)) ≈ N0f8(0.439)
+end

test/region_growing.jl

@@ -156,29 +156,33 @@ of_mean_type(p::Colorant) = ImageSegmentation.meantype(typeof(p))(p)
     expected = ones(Int, (3,3))
     expected[1:3,3] .= 2
     expected_labels = [1,2]
-    expected_means = Dict(1=>RGB{Float64}(0.3,1.0,0.0), 2=>RGB{Float64}(0.0,0.0,0.0))
+    expected_means = Dict(1=>RGB{Float32}(0.3,1.0,0.0), 2=>RGB{Float32}(0.0,0.0,0.0))
     expected_count = Dict(1=>6, 2=>3)
 
     seg = seeded_region_growing(img, seeds)
     @test all(label->(label in expected_labels), seg.segment_labels)
     @test all(label->(label in seg.segment_labels), expected_labels)
     @test expected_count == seg.segment_pixel_count
-    @test expected_means == seg.segment_means
+    @test expected_means[1] ≈ seg.segment_means[1]
+    @test expected_means[2] ≈ seg.segment_means[2]
     @test seg.image_indexmap == expected
 
     expected = ones(Int, (3,3))
     expected[1:3,2] .= 0
     expected[1:3,3] .= 2
     expected_labels = [0,1,2]
-    expected_means = Dict(1=>RGB{Float64}(0.4,1.0,0.0), 2=>RGB{Float64}(0.0,0.0,0.0))
+    expected_means = Dict(1=>RGB{Float32}(0.4,1.0,0.0), 2=>RGB{Float32}(0.0,0.0,0.0))
     expected_count = Dict(0=>3, 1=>3, 2=>3)
 
-    seg = seeded_region_growing(img, seeds, [3,3], (c1,c2)->abs(of_mean_type(c1).r - of_mean_type(c2).r))
+    seg = seeded_region_growing(img, seeds, (3,3), (c1,c2)->abs(of_mean_type(c1).r - of_mean_type(c2).r))
     @test all(label->(label in expected_labels), seg.segment_labels)
     @test all(label->(label in seg.segment_labels), expected_labels)
     @test expected_count == seg.segment_pixel_count
     @test expected_means == seg.segment_means
     @test seg.image_indexmap == expected
+    @info "The deprecation warning below is expected"   # but can be deleted eventually!
+    segd = seeded_region_growing(img, seeds, [3,3], (c1,c2)->abs(of_mean_type(c1).r - of_mean_type(c2).r))
+    @test labels_map(segd) == labels_map(seg)
 end
 
 @testset "Unseeded Region Growing" begin
@@ -296,11 +300,13 @@ end
     expected_means = Dict(1=>of_mean_type(img[1,1]), 3=>of_mean_type(img[1,3]), 2=>of_mean_type(img[1,2]))
     expected_count = Dict(1=>3, 2=>3, 3=>3)
 
-    seg = unseeded_region_growing(img, 0.2, [3,3], (c1,c2)->abs(of_mean_type(c1).r - of_mean_type(c2).r))
+    seg = unseeded_region_growing(img, 0.2, (3,3), (c1,c2)->abs(of_mean_type(c1).r - of_mean_type(c2).r))
     @test all(label->(label in expected_labels), seg.segment_labels)
     @test all(label->(label in seg.segment_labels), expected_labels)
     @test expected_count == seg.segment_pixel_count
     @test expected_means == seg.segment_means
     @test seg.image_indexmap == expected
-
+    @info "The deprecation warning below is expected"   # but can be deleted eventually!
+    segd = unseeded_region_growing(img, 0.2, [3,3], (c1,c2)->abs(of_mean_type(c1).r - of_mean_type(c2).r))
+    @test labels_map(segd) == labels_map(seg)
 end

test/runtests.jl

@@ -1,14 +1,16 @@
 using ImageSegmentation, Images, Test, SimpleWeightedGraphs, LightGraphs, StaticArrays, DataStructures
-using RegionTrees: isleaf, Cell, split! 
+using RegionTrees: isleaf, Cell, split!
 using MetaGraphs: MetaGraph, clear_props!, get_prop, has_prop, set_prop!, props, vertices
 
 using Documenter
-doctest(ImageSegmentation, manual = false)
+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")