JuliaImages / ImageSegmentation.jl
1

2 12
default_diff_fn(c1::CT1,c2::CT2) where {CT1<:Union{Colorant,Real}, CT2<:Union{Colorant,Real}} = sqrt(sum(abs2,(c1)-Images.accum(CT2)(c2)))
3

4
"""
5
    seg_img = seeded_region_growing(img, seeds, [kernel_dim], [diff_fn])
6
    seg_img = seeded_region_growing(img, seeds, [neighbourhood], [diff_fn])
7

8
Segments the N-D image `img` using the seeded region growing algorithm
9
and returns a [`SegmentedImage`](@ref) containing information about the segments.
10

11
# Arguments:
12
* `img`             :  N-D image to be segmented (arbitrary axes are allowed)
13
* `seeds`           :  `Vector` containing seeds. Each seed is a Tuple of a
14
                       CartesianIndex{N} and a label. See below note for more
15
                       information on labels.
16
* `kernel_dim`      :  (Optional) `Vector{Int}` having length N or a `NTuple{N,Int}`
17
                       whose ith element is an odd positive integer representing
18
                       the length of the ith edge of the N-orthotopic neighbourhood
19
* `neighbourhood`   :  (Optional) Function taking CartesianIndex{N} as input and
20
                       returning the neighbourhood of that point.
21
* `diff_fn`         :  (Optional) Function that returns a difference measure(δ)
22
                       between the mean color of a region and color of a point
23

24
!!! note
25
    The labels attached to points must be positive integers, although multiple
26
    points can be assigned the same label. The output includes a labelled array
27
    that has same indexing as that of input image. Every index is assigned to
28
    either one of labels or a special label '0' indicating that the algorithm
29
    was unable to assign that index to a unique label.
30

31
# Examples
32

33
```jldoctest; setup = :(using Images, ImageSegmentation)
34
julia> img = zeros(Gray{N0f8},4,4);
35

36
julia> img[2:4,2:4] .= 1;
37

38
julia> seeds = [(CartesianIndex(3,1),1),(CartesianIndex(2,2),2)];
39

40
julia> seg = seeded_region_growing(img, seeds);
41

42
julia> labels_map(seg)
43
4×4 $(Matrix{Int}):
44
 1  1  1  1
45
 1  2  2  2
46
 1  2  2  2
47
 1  2  2  2
48
```
49

50
# Citation:
51

52
Albert Mehnert, Paul Jackaway (1997), "An improved seeded region growing algorithm",
53
Pattern Recognition Letters 18 (1997), 1065-1071
54
"""
55 6
function seeded_region_growing(img::AbstractArray{CT,N}, seeds::AbstractVector{Tuple{CartesianIndex{N},Int}},
56 6
    kernel_dim::Union{Vector{Int}, NTuple{N, Int}} = ntuple(i->3,N), diff_fn::Function = default_diff_fn) where {CT<:Union{Colorant,Real}, N}
57 12
    length(kernel_dim) == N || error("Dimension count of image and kernel_dim do not match")
58 12
    for dim in kernel_dim
59 12
        dim > 0 || error("Dimensions of the kernel must be positive")
60 12
        isodd(dim) || error("Dimensions of the kernel must be odd")
61
    end
62 12
    pt = CartesianIndex(ntuple(i->kernel_dim[i]÷2, N))
63 12
    neighbourhood_gen(t) = c->CartesianIndices(_colon(c-t,c+t))
64 12
    seeded_region_growing(img, seeds, neighbourhood_gen(pt), diff_fn)
65
end
66

67 6
function seeded_region_growing(img::AbstractArray{CT,N}, seeds::AbstractVector{Tuple{CartesianIndex{N},Int}},
68
    neighbourhood::Function, diff_fn::Function = default_diff_fn) where {CT<:Union{Colorant,Real}, N}
69

70 12
    _QUEUE_SZ = 10
71

72
    # Check if labels are positive integers
73 12
    for seed in seeds
74 12
        seed[2] > 0 || error("Seed labels need to be positive integers!")
75
    end
76

77 6
    TM = meantype(CT)
78

79
    # Required data structures
80 12
    result              =   fill(-1, axes(img))                                     # Array to store labels
81 12
    nhq                 =   Queue{CartesianIndex{N}}(_QUEUE_SZ)                     # Neighbours holding queue
82 12
    pq                  =   PriorityQueue{CartesianIndex{N}, Float64}()             # Priority Queue to hold the pixels of same δ value
83 12
    labelsq             =   Queue{Int}(_QUEUE_SZ)                                   # Queue to hold labels
84 12
    holdingq            =   Queue{CartesianIndex{N}}(_QUEUE_SZ)                     # Queue to hold points corresponding to the labels in `labelsq`
85 12
    region_means        =   Dict{Int, TM}()                                         # A map containing (label, mean) pairs
86 12
    region_pix_count    =   Dict{Int, Int}()                                        # A map containing (label, pixel_count) pairs
87 6
    labels              =   Vector{Int}()                                           # A vector containing list of labels
88

89
    # Labelling initial seeds and initialising `region_means` and `region_pix_count`
90 12
    for seed in seeds
91 6
        result[seed[1]] = seed[2]
92 12
        region_pix_count[seed[2]] = get(region_pix_count, seed[2], 0) + 1
93 12
        region_means[seed[2]] = get(region_means, seed[2], zero(TM)) + (img[seed[1]] - get(region_means, seed[2], zero(TM)))/(region_pix_count[seed[2]])
94 12
        if ! (seed[2] in labels)
95 12
            push!(labels, seed[2])
96
        end
97
    end
98

99
    #=  Labeling scheme for the Array "result"-
100
            Unlabelled => -1
101
            In nhq => -2
102
            In pq => -3
103
            Tied => 0
104
            Labelled => Seed value
105
    =#
106

107
    # Enqueue all the neighbours of seeds into `nhq`
108 12
    for seed in seeds
109 12
        for point in neighbourhood(seed[1])
110 12
            if point != seed[1] && checkbounds(Bool, img, point) && result[point] == -1
111 6
                enqueue!(nhq, point)
112 12
                @inbounds result[point] = -2
113
            end
114
        end
115
    end
116

117 12
    while !isempty(pq) || !isempty(nhq)
118

119 12
        while !isempty(nhq)
120
            # For every neighbouring point, get the minimum δ
121 6
            p = dequeue!(nhq)
122 12
            δ = Inf
123 12
            for point in neighbourhood(p)
124 12
                if point != p && checkbounds(Bool, img, point) && result[point] > 0
125 12
                    curr_diff = diff_fn(region_means[result[point]], img[p])
126 12
                    if δ > curr_diff
127 12
                        δ = curr_diff
128
                    end
129
                end
130
            end
131 12
            pq[p] = δ
132 12
            @inbounds result[p] = -3
133
        end
134

135
        # Get the pixels with minimum δ from `pq` and add them to `holdingq` and their labels to `labelsq`
136 12
        if !isempty(pq)
137 6
            δ_min = peek(pq)[2]
138
        end
139 12
        while (!isempty(pq) && isapprox(peek(pq)[2], δ_min)) #, atol=1e-8))
140 12
            p = dequeue!(pq)
141 6
            mindifflabel = -1
142 12
            mindiff = Inf
143 12
            for point in neighbourhood(p)
144 12
                if point!=p && checkbounds(Bool, img, point) && result[point] > 0
145 12
                    if mindifflabel < 0
146 6
                        mindifflabel = result[point]
147 12
                        mindiff = diff_fn(region_means[result[point]], img[p])
148 12
                    elseif mindifflabel != result[point]
149 12
                        curr_diff = diff_fn(region_means[result[point]], img[p])
150 12
                        if curr_diff < mindiff
151 6
                            mindiff = curr_diff
152 12
                            mindifflabel = result[point]
153 12
                        elseif isapprox(curr_diff, mindiff) #, atol=1e-8)
154 12
                            mindifflabel = 0
155
                        end
156
                    end
157
                end
158
            end
159 12
            if mindifflabel != 0                # new labels, including resolved ties in last step
160 6
                enqueue!(labelsq, mindifflabel)
161 6
                enqueue!(holdingq, p)
162 12
                @inbounds result[p] = -2
163 12
            elseif result[p] == -3              # new ties are enqueued in pq
164 12
                pq[p] = Inf
165 12
                result[p] = 0
166
            end                                 # old ties are not enqueued again
167
        end
168

169
        # Add label to each point in `holdingq` and add their neighbours to `nhq`
170 12
        while !isempty(holdingq)
171 6
            label = dequeue!(labelsq)
172 6
            p = dequeue!(holdingq)
173 6
            result[p] = label
174 12
            region_pix_count[label] += 1
175 12
            region_means[label] += (img[p] - region_means[label])/(region_pix_count[label])
176 12
            for point in neighbourhood(p)
177 12
                if point!=p && checkbounds(Bool, img, point) && (result[point] == -1 || result[point] == -3)
178 6
                    enqueue!(nhq, point)
179 12
                    @inbounds result[point] = -2
180
                end
181
            end
182
        end
183

184
    end
185

186 6
    c0 = count(i->(i==0),result)
187 12
    if c0 != 0
188 6
        push!(labels, 0)
189 12
        region_pix_count[0] = c0
190
    end
191 12
    SegmentedImage(result, labels, region_means, region_pix_count)
192
end
193

194

195
"""
196
    seg_img = unseeded_region_growing(img, threshold, [kernel_dim], [diff_fn])
197
    seg_img = unseeded_region_growing(img, threshold, [neighbourhood], [diff_fn])
198

199
Segments the N-D image using automatic (unseeded) region growing algorithm
200
and returns a [`SegmentedImage`](@ref) containing information about the segments.
201

202
# Arguments:
203
* `img`             :  N-D image to be segmented (arbitrary axes are allowed)
204
* `threshold`       :  Upper bound of the difference measure (δ) for considering
205
                       pixel into same segment
206
* `kernel_dim`      :  (Optional) `Vector{Int}` having length N or a `NTuple{N,Int}`
207
                       whose ith element is an odd positive integer representing
208
                       the length of the ith edge of the N-orthotopic neighbourhood
209
* `neighbourhood`   :  (Optional) Function taking CartesianIndex{N} as input and
210
                       returning the neighbourhood of that point.
211
* `diff_fn`         :  (Optional) Function that returns a difference measure (δ)
212
                       between the mean color of a region and color of a point
213

214
# Examples
215

216
```jldoctest; setup = :(using Images, ImageSegmentation)
217
julia> img = zeros(Gray{N0f8},4,4);
218

219
julia> img[2:4,2:4] .= 1;
220

221
julia> seg = unseeded_region_growing(img, 0.2);
222

223
julia> labels_map(seg)
224
4×4 $(Matrix{Int}):
225
 1  1  1  1
226
 1  2  2  2
227
 1  2  2  2
228
 1  2  2  2
229
```
230

231
"""
232 6
function unseeded_region_growing(img::AbstractArray{CT,N}, threshold::Real,
233 6
    kernel_dim::Union{Vector{Int}, NTuple{N, Int}} = ntuple(i->3,N), diff_fn::Function = default_diff_fn) where {CT<:Colorant, N}
234 12
    length(kernel_dim) == N || error("Dimension count of image and kernel_dim do not match")
235 12
    for dim in kernel_dim
236 12
        dim > 0 || error("Dimensions of the kernel must be positive")
237 12
        isodd(dim) || error("Dimensions of the kernel must be odd")
238
    end
239 12
    pt = CartesianIndex(ntuple(i->kernel_dim[i]÷2, N))
240 12
    neighbourhood_gen(t) = c->CartesianIndices(_colon(c-t,c+t))
241 12
    unseeded_region_growing(img, threshold, neighbourhood_gen(pt), diff_fn)
242
end
243

244 6
function unseeded_region_growing(img::AbstractArray{CT,N}, threshold::Real, neighbourhood::Function, diff_fn = default_diff_fn) where {CT<:Colorant,N}
245 12
    TM = meantype(CT)
246

247
    # Required data structures
248 12
    result                  =   fill(-1, axes(img))                             # Array to store labels
249 12
    neighbours              =   PriorityQueue{CartesianIndex{N},Float64}()      # Priority Queue containing boundary pixels with δ as the priority
250 12
    region_means            =   Dict{Int, TM}()                                 # A map containing (label, mean) pairs
251 12
    region_pix_count        =   Dict{Int, Int}()                                # A map containing (label, pixel_count) pairs
252 6
    labels                  =   Vector{Int}()                                   # Vector containing assigned labels
253

254
    # Initialize data structures
255 6
    start_point = first(CartesianIndices(axes(img)))
256 6
    result[start_point] = 1
257 6
    push!(labels, 1)
258 12
    region_means[1] = img[start_point]
259 12
    region_pix_count[1] = 1
260

261
    # Enqueue neighouring points of `start_point`
262 12
    for p in neighbourhood(start_point)
263 12
        if p != start_point && checkbounds(Bool, img, p) && result[p] == -1
264 12
            enqueue!(neighbours, p, diff_fn(region_means[result[start_point]], img[p]))
265
        end
266
    end
267

268 12
    while !isempty(neighbours)
269 12
        point = dequeue!(neighbours)
270 12
        δ = Inf
271 6
        minlabel = -1
272 6
        pixelval = img[point]
273 12
        for p in neighbourhood(point)
274 12
            if p != point && checkbounds(Bool, img, p) && result[p] > 0
275 6
                curδ = diff_fn(region_means[result[p]], pixelval)
276 12
                if curδ < δ
277 6
                    δ = curδ
278 12
                    minlabel = result[p]
279
                end
280
            end
281
        end
282

283 12
        if δ < threshold
284
            # Assign point to minlabel
285 12
            result[point] = minlabel
286
        else
287
            # Select most appropriate label
288 12
            δ = Inf
289 6
            minlabel = -1
290 12
            for label in labels
291 6
                curδ = diff_fn(region_means[label], pixelval)
292 12
                if curδ < δ
293 6
                    δ = curδ
294 12
                    minlabel = label
295
                end
296
            end
297

298 12
            if δ < threshold
299 12
                result[point] = minlabel
300
            else
301
                # Assign point to a new label
302 6
                minlabel = length(labels) + 1
303 6
                push!(labels, minlabel)
304 6
                result[point] = minlabel
305
            end
306
        end
307

308
        #Update region_means
309 12
        region_pix_count[minlabel] = get(region_pix_count, minlabel, 0) + 1
310 12
        region_means[minlabel] = get(region_means, minlabel, zero(TM)) + (pixelval - get(region_means, minlabel, zero(TM)))/(region_pix_count[minlabel])
311

312
        # Enqueue neighbours of `point`
313 12
        for p in neighbourhood(point)
314 12
            if checkbounds(Bool, img, p) && result[p] == -1
315 12
                if haskey(neighbours, p)
316 12
                    continue
317
                end
318 12
                δ = Inf
319 12
                for tp in neighbourhood(p)
320 12
                    if checkbounds(Bool, img, tp) && tp != p && result[tp] > 0
321 12
                        δ = min(δ, diff_fn(region_means[result[tp]], img[p]))
322
                    end
323
                end
324 12
                enqueue!(neighbours, p, δ)
325
            end
326
        end
327

328
    end
329 12
    SegmentedImage(result, labels, region_means, region_pix_count)
330
end

Read our documentation on viewing source code .

Loading