JuliaImages / ImageSegmentation.jl
1 12
sharpness(img::AbstractArray{CT,N}) where {CT<:Images.NumberLike,N} = var(imfilter(img, Kernel.Laplacian(ntuple(i->true, Val(N)))))
2

3 6
function adaptive_thres(img::AbstractArray{CT,N}, block::NTuple{N,Int}) where {CT<:Images.NumberLike,N}
4 12
    threshold = zeros(Float64, block)
5 12
    block_length = CartesianIndex(ntuple(i->ceil(Int,length(axes(img,i))/block[i]),Val(N)))
6 12
    net_s = sharpness(img)
7 6
    net_var = var(img)
8 6
    net_end = last(CartesianIndices(axes(img)))
9 12
    for i in CartesianIndices(block)
10 6
        si = CartesianIndex(ntuple(j->(i[j]-1)*block_length[j]+1,Val(N)))
11 6
        ei = min(si + block_length - _oneunit(si), net_end)
12 6
        wi = view(img, map((i,j)->i:j, si.I, ei.I)...)
13 12
        threshold[i] = 0.02 + min(sharpness(wi)/net_s*0.04,0.1) + min(var(wi)/net_var*0.04,0.1)
14
    end
15 12
    threshold
16
end
17

18 6
getscalar(A::AbstractArray{T,N}, i::CartesianIndex{N}, block_length::CartesianIndex{N}) where {T<:Real,N} =
19 6
    A[CartesianIndex(ntuple(j->(i[j]-1)÷block_length[j]+1, Val(N)))]
20

21 6
getscalar(a::Real, i...) = a
22

23 12
fast_scanning(img::AbstractArray{CT,N}, block::NTuple{N,Int} = ntuple(i->4,Val(N))) where {CT<:Images.NumberLike,N} =
24
    fast_scanning(img, adaptive_thres(img, block))
25

26
"""
27
    seg_img = fast_scanning(img, threshold, [diff_fn])
28

29
Segments the N-D image using a fast scanning algorithm and returns a
30
[`SegmentedImage`](@ref) containing information about the segments.
31

32
# Arguments:
33
* `img`         : N-D image to be segmented (arbitrary axes are allowed)
34
* `threshold`   : Upper bound of the difference measure (δ) for considering
35
                  pixel into same segment; an `AbstractArray` can be passed
36
                  having same number of dimensions as that of `img` for adaptive
37
                  thresholding
38
* `diff_fn`     : (Optional) Function that returns a difference measure (δ)
39
                  between the mean color of a region and color of a point
40

41
# Examples:
42

43
```jldoctest; setup = :(using Images, ImageSegmentation)
44
julia> img = zeros(Float64, (3,3));
45

46
julia> img[2,:] .= 0.5;
47

48
julia> img[:,2] .= 0.6;
49

50
julia> seg = fast_scanning(img, 0.2);
51

52
julia> labels_map(seg)
53
3×3 $(Matrix{Int}):
54
 1  4  5
55
 4  4  4
56
 3  4  6
57
```
58

59
# Citation:
60

61
Jian-Jiun Ding, Cheng-Jin Kuo, Wen-Chih Hong,
62
"An efficient image segmentation technique by fast scanning and adaptive merging"
63
"""
64 6
function fast_scanning(img::AbstractArray{CT,N}, threshold::Union{AbstractArray,Real}, diff_fn::Function = default_diff_fn) where {CT<:Union{Colorant,Real},N}
65

66 12
    if threshold isa AbstractArray
67 6
        ndims(img) == ndims(threshold) || error("Dimension count of image and threshold do not match")
68
    end
69

70
    # Neighbourhood function
71 12
    _diagmN = Diagonal([1 for i in 1:N])
72 12
    half_region::NTuple{N,CartesianIndex{N}} = ntuple(i-> CartesianIndex{N}(ntuple(j->_diagmN[j,i], Val(N))), Val(N))
73 6
    neighbourhood(x) = ntuple(i-> x-half_region[i], Val(N))
74

75
    # Required data structures
76 6
    TM = meantype(CT)
77 12
    result              =   fill(-1, axes(img))                             # Array to store labels
78 12
    region_means        =   Dict{Int, TM}()                                 # A map conatining (label, mean) pairs
79 12
    region_pix_count    =   Dict{Int, Int}()                                # A map conatining (label, count) pairs
80 6
    temp_labels         =   IntDisjointSets(0)                              # Disjoint set to map labels to their equivalence class
81 6
    v_neigh             =   MVector{N,Int}(undef)                           # MVector to store valid neighbours
82

83 12
    block_length = CartesianIndex(ntuple(i->ceil(Int,length(axes(img,i))/size(threshold,i)), Val(N)))
84

85 12
    for point in CartesianIndices(axes(img))
86 6
        sz = 0
87 6
        same_label = true
88 6
        prev_label = 0
89 12
        for p in neighbourhood(point)
90 12
            if checkbounds(Bool, img, p)
91 6
                root_p = find_root!(temp_labels, result[p])
92 12
                if diff_fn(region_means[root_p], img[point]) < getscalar(threshold, point, block_length)
93 12
                    if prev_label == 0
94 12
                        prev_label = root_p
95 12
                    elseif prev_label != root_p
96 6
                        same_label = false
97
                    end
98 12
                    sz += 1
99 12
                    v_neigh[sz] = find_root!(temp_labels, root_p)
100
                end
101
            end
102
        end
103

104
        # If no valid label found
105 12
        if sz == 0
106
            # Assign a new label
107 12
            new_label = push!(temp_labels)
108 6
            result[point] = new_label
109 12
            region_means[new_label] = img[point]
110 12
            region_pix_count[new_label] = 1
111

112
        # If all labels are same
113 12
        elseif same_label
114 6
            result[point] = prev_label
115 12
            region_pix_count[prev_label] += 1
116 12
            region_means[prev_label] += (img[point] - region_means[prev_label])/(region_pix_count[prev_label])
117
        else
118
            # Merge segments and assign to this new label
119 6
            union_label = v_neigh[1]
120 12
            for i in 1:sz
121 12
                union_label = union!(temp_labels, union_label, v_neigh[i])
122
            end
123 12
            result[point] = union_label
124 12
            region_pix_count[union_label] += 1
125 12
            region_means[union_label] += (img[point] - region_means[union_label])/(region_pix_count[union_label])
126

127 12
            for i in 1:sz
128 12
                if v_neigh[i] != union_label && haskey(region_pix_count, v_neigh[i])
129 12
                    region_pix_count[union_label] += region_pix_count[v_neigh[i]]
130 12
                    region_means[union_label] += (region_means[v_neigh[i]] - region_means[union_label])*region_pix_count[v_neigh[i]]/region_pix_count[union_label]
131

132
                    # Remove label v_neigh[i] from region_means, region_pix_count
133 6
                    delete!(region_pix_count,v_neigh[i])
134 12
                    delete!(region_means,v_neigh[i])
135
                end
136
            end
137
        end
138
    end
139

140 12
    for point in CartesianIndices(axes(img))
141 12
        result[point] = find_root!(temp_labels, result[point])
142
    end
143

144 12
    SegmentedImage(result, unique(temp_labels.parents), region_means, region_pix_count)
145
end

Read our documentation on viewing source code .

Loading