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 .