Merge pull request #73 from JuliaImages/teh/polish

Lots of polish

Author: timholy

src/ImageSegmentation.jl

@@ -6,6 +6,8 @@ using LinearAlgebra, Statistics
 using DataStructures, StaticArrays, ImageCore, ImageFiltering, ImageMorphology, LightGraphs, SimpleWeightedGraphs, RegionTrees, Distances, StaticArrays, Clustering, MetaGraphs
 import Clustering: kmeans, fuzzy_cmeans
 
+const PairOrTuple{K,V} = Union{Pair{K,V},Tuple{K,V}}
+
 include("compat.jl")
 include("core.jl")
 include("region_growing.jl")

src/compat.jl

@@ -1,7 +1,7 @@
-if VERSION <= v"1.0.5"
+@static if Base.VERSION <= v"1.0.5"
     # https://github.com/JuliaLang/julia/pull/29442
     _oneunit(::CartesianIndex{N}) where {N} = _oneunit(CartesianIndex{N})
     _oneunit(::Type{CartesianIndex{N}}) where {N} = CartesianIndex(ntuple(x -> 1, Val(N)))
 else
-    _oneunit = Base.oneunit
+    const _oneunit = Base.oneunit
 end

src/core.jl

@@ -95,7 +95,7 @@ function region_adjacency_graph(s::SegmentedImage, weight_fn::Function)
         while !isempty(t)
             temp = pop!(t)
             visited[temp] = true
-            for J in CartesianIndices(_colon(max(Ibegin, temp-I1), min(Iend, temp+I1)))
+            for J in _colon(max(Ibegin, temp-I1), min(Iend, temp+I1))
                 if s.image_indexmap[temp] != s.image_indexmap[J]
                     push!(n,s.image_indexmap[J])
                 elseif !visited[J]
@@ -268,7 +268,7 @@ end
 
 # Once Base has colon defined here we can replace this
 _colon(I::CartesianIndex{N}, J::CartesianIndex{N}) where N =
-    map((i,j) -> i:j, Tuple(I), Tuple(J))
+    CartesianIndices(map((i,j) -> i:j, Tuple(I), Tuple(J)))
 
 
 """

src/fast_scanning.jl

@@ -23,46 +23,10 @@ getscalar(a::Real, i...) = a
 fast_scanning(img::AbstractArray{CT,N}, block::NTuple{N,Int} = ntuple(i->4,Val(N))) where {CT<:ImageCore.NumberLike,N} =
     fast_scanning(img, adaptive_thres(img, block))
 
-"""
-    seg_img = fast_scanning(img, threshold, [diff_fn])
-
-Segments the N-D image using a fast scanning algorithm and returns a
-[`SegmentedImage`](@ref) containing information about the segments.
-
-# Arguments:
-* `img`         : N-D image to be segmented (arbitrary axes are allowed)
-* `threshold`   : Upper bound of the difference measure (δ) for considering
-                  pixel into same segment; an `AbstractArray` can be passed
-                  having same number of dimensions as that of `img` for adaptive
-                  thresholding
-* `diff_fn`     : (Optional) Function that returns a difference measure (δ)
-                  between the mean color of a region and color of a point
-
-# Examples:
-
-```jldoctest; setup = :(using ImageCore, ImageMorphology, ImageSegmentation)
-julia> img = zeros(Float64, (3,3));
-
-julia> img[2,:] .= 0.5;
-
-julia> img[:,2] .= 0.6;
-
-julia> seg = fast_scanning(img, 0.2);
-
-julia> labels_map(seg)
-3×3 $(Matrix{Int}):
- 1  4  5
- 4  4  4
- 3  4  6
-```
-
-# Citation:
-
-Jian-Jiun Ding, Cheng-Jin Kuo, Wen-Chih Hong,
-"An efficient image segmentation technique by fast scanning and adaptive merging"
-"""
-function fast_scanning(img::AbstractArray{CT,N}, threshold::Union{AbstractArray,Real}, diff_fn::Function = default_diff_fn) where {CT<:Union{Colorant,Real},N}
+fast_scanning(img::AbstractArray{CT,N}, threshold::Union{AbstractArray,Real}, diff_fn::Function = default_diff_fn) where {CT<:Union{Colorant,Real},N} =
+    fast_scanning!(fill(-1, axes(img)), img, threshold, diff_fn)
 
+function fast_scanning!(result, img::AbstractArray{CT,N}, threshold::Union{AbstractArray,Real}, diff_fn::DF = default_diff_fn) where {CT<:Union{Colorant,Real},N,DF<:Function}
     if threshold isa AbstractArray
         ndims(img) == ndims(threshold) || error("Dimension count of image and threshold do not match")
     end
@@ -74,7 +38,6 @@ function fast_scanning(img::AbstractArray{CT,N}, threshold::Union{AbstractArray,
 
     # Required data structures
     TM = meantype(CT)
-    result              =   fill(-1, axes(img))                             # Array to store labels
     region_means        =   Dict{Int, TM}()                                 # A map conatining (label, mean) pairs
     region_pix_count    =   Dict{Int, Int}()                                # A map conatining (label, count) pairs
     temp_labels         =   IntDisjointSets(0)                              # Disjoint set to map labels to their equivalence class

src/felzenszwalb.jl

@@ -1,36 +1,17 @@
 """
-```
-segments                = felzenszwalb(img, k, [min_size])
-index_map, num_segments = felzenszwalb(edges, num_vertices, k, [min_size])
-```
+    index_map, num_segments = felzenszwalb(edges, num_vertices, k, min_size=0)
 
-Segments an image using Felzenszwalb's graph-based algorithm. The function can be used in either of two ways -
-
-1. `segments = felzenszwalb(img, k, [min_size])`
-
-Segments an image using Felzenszwalb's segmentation algorithm and returns the result as `SegmentedImage`. The algorithm uses
-euclidean distance in color space as edge weights for the region adjacency graph.
-
-Parameters:
--    img            = input image
--    k              = Threshold for region merging step. Larger threshold will result in bigger segments.
--    min_size       = Minimum segment size
-
-2. `index_map, num_segments = felzenszwalb(edges, num_vertices, k, [min_size])`
-
-Segments an image represented as Region Adjacency Graph(RAG) using Felzenszwalb's segmentation algorithm. Each pixel/region
- corresponds to a node in the graph and weights on each edge measure the dissimilarity between pixels.
+Segment an image represented as Region Adjacency Graph(RAG) using Felzenszwalb's segmentation algorithm. Each pixel/region
+corresponds to a node in the graph and weights on each edge measure the dissimilarity between pixels.
 The function returns the number of segments and index mapping from nodes of the RAG to segments.
 
 Parameters:
--    edges          = Array of edges in RAG. Each edge is represented as `ImageEdge`.
--    num_vertices   = Number of vertices in RAG
--    k              = Threshold for region merging step. Larger threshold will result in bigger segments.
--    min_size       = Minimum segment size
-
-
+- `edges`:        Array of edges in RAG. Each edge is represented as `ImageEdge`.
+- `num_vertices`: Number of vertices in RAG
+- `k`:            Threshold for region merging step. Larger threshold will result in bigger segments.
+- `min_size`:     Minimum segment size (in # pixels)
 """
-function felzenszwalb(edges::Array{ImageEdge}, num_vertices::Int, k::Real, min_size::Int = 0)
+function felzenszwalb(edges::Array{ImageEdge}, num_vertices::Int, k::Float64, min_size::Int = 0)
 
     num_edges = length(edges)
     G = IntDisjointSets(num_vertices)
@@ -61,7 +42,7 @@ function felzenszwalb(edges::Array{ImageEdge}, num_vertices::Int, k::Real, min_s
         end
     end
 
-    segments = OrderedSet()
+    segments = OrderedSet{Int}()
     for i in 1:num_vertices
         push!(segments, find_root!(G, i))
     end
@@ -72,35 +53,59 @@ function felzenszwalb(edges::Array{ImageEdge}, num_vertices::Int, k::Real, min_s
         segments2index[s] = i
     end
 
-    index_map = Array{Int}(undef, num_vertices)
+    index_map = Vector{Int}(undef, num_vertices)
     for i in 1:num_vertices
         index_map[i] = segments2index[find_root!(G, i)]
     end
 
     return index_map, num_sets
 end
+felzenszwalb(edges::Array{ImageEdge}, num_vertices::Integer, k::Real, min_size::Integer = 0) =
+    felzenszwalb(edges, convert(Int, num_vertices)::Int, convert(Float64, k)::Float64, convert(Int, min_size)::Int)
 
 meantype(::Type{T}) where T = typeof(zero(accum_type(T))/2)
 
-function felzenszwalb(img::AbstractArray{T, 2}, k::Real, min_size::Int = 0) where T<:Union{Real,Color}
+"""
+    segments = felzenszwalb(img, k, [min_size])
+
+Segment an image using Felzenszwalb's segmentation algorithm and returns the result as `SegmentedImage`.
+The algorithm uses euclidean distance in color space as edge weights for the region adjacency graph.
+
+Parameters:
+- `img`:      input image
+- `k`:        Threshold for region merging step. Larger threshold will result in bigger segments.
+- `min_size`: Minimum segment size (in # pixels)
+"""
+function felzenszwalb(img::AbstractArray{T}, k::Real, min_size::Int = 0) where T<:Union{Real,Color}
+
+    sz = size(img)
+    num_vertices = prod(sz)
 
-    rows, cols = size(img)
-    num_vertices = rows*cols
-    num_edges = 4*rows*cols - 3*rows - 3*cols + 2
-    edges = Array{ImageEdge}(undef, num_edges)
+    R = CartesianIndices(img)
+    L = LinearIndices(img)
+    Ibegin, Iend = first(R), last(R)
+    I1 = _oneunit(Ibegin)
 
-    R = CartesianIndices(size(img))
-    I1, Iend = first(R), last(R)
-    num = 1
+    # Compute the number of entries per pixel (other than at the image edges)
+    num_edges = 0
+    for I in _colon(-I1, I1)
+        I >= zero(I1) && continue
+        num_edges += 1
+    end
+    num_edges *= num_vertices   # now the number for the whole image
+    edges = Vector{ImageEdge}(undef, num_edges)
+
+    num = 0
     for I in R
-        for J in CartesianIndices(_colon(max(I1, I-I1), min(Iend, I+I1)))
+        imgI = img[I]
+        for J in _colon(max(Ibegin, I-I1), min(Iend, I+I1))
             if I >= J
                 continue
             end
-            edges[num] = ImageEdge((I[2]-1)*rows+I[1], (J[2]-1)*rows+J[1], sqrt(sum(abs2,(img[I])-meantype(T)(img[J]))))
-            num += 1
+            edges[num+=1] = ImageEdge(L[I], L[J], sqrt(abs2(imgI-meantype(T)(img[J]))))
         end
     end
+    deleteat!(edges, num+1:num_edges)   # compensate for the ones we were missing at the image edges
 
     index_map, num_segments = felzenszwalb(edges, num_vertices, k, min_size)
 
@@ -109,12 +114,10 @@ function felzenszwalb(img::AbstractArray{T, 2}, k::Real, min_size::Int = 0) wher
     region_means        = Dict{Int, meantype(T)}()
     region_pix_count    = Dict{Int, Int}()
 
-    for j in axes(img, 2)
-        for i in axes(img, 1)
-            result[i, j] = index_map[(j-1)*rows+i]
-            region_pix_count[result[i,j]] = get(region_pix_count, result[i, j], 0) + 1
-            region_means[result[i,j]] = get(region_means, result[i,j], zero(meantype(T))) + (img[i, j] - get(region_means, result[i,j], zero(meantype(T))))/region_pix_count[result[i,j]]
-        end
+    for I in R
+        result[I] = index_map[L[I]]
+        region_pix_count[result[I]] = get(region_pix_count, result[I], 0) + 1
+        region_means[result[I]] = get(region_means, result[I], zero(meantype(T))) + (img[I] - get(region_means, result[I], zero(meantype(T))))/region_pix_count[result[I]]
     end
 
     return SegmentedImage(result, labels, region_means, region_pix_count)

src/meanshift.jl

@@ -41,7 +41,7 @@ function meanshift(img::Array{CT, 2}, spatial_radius::Real, range_radius::Real;
         I1, Iend = first(R), last(R)
         I = CartesianIndex(rowbin(pt[1]), colbin(pt[2]), colorbin(pt[3]))
 
-        for J in CartesianIndices(_colon(max(I1, I-_oneunit(I)), min(Iend, I+_oneunit(I))))
+        for J in _colon(max(I1, I-_oneunit(I)), min(Iend, I+_oneunit(I)))
             for point in buckets[J]
                 if dist2(pt, SVector(point[1], point[2], img[point])) <= 1
                     num += SVector(point[1], point[2], img[point])