diff --git a/REQUIRE b/REQUIRE index 54190bf..3af5429 100644 --- a/REQUIRE +++ b/REQUIRE @@ -1,2 +1,3 @@ julia 1.0 ImageCore 0.7.4 +Images 0.13.0 \ No newline at end of file diff --git a/src/ImageMorphology.jl b/src/ImageMorphology.jl index d255977..cadaf7c 100644 --- a/src/ImageMorphology.jl +++ b/src/ImageMorphology.jl @@ -3,10 +3,15 @@ __precompile__() module ImageMorphology using ImageCore +using Images + +abstract type AbstractSkeletonizationAlgorithm end +struct MedialAxisTransform <: AbstractSkeletonizationAlgorithm end + include("dilation_and_erosion.jl") include("thinning.jl") - +include("skeletonization.jl") export dilate, erode, @@ -16,7 +21,8 @@ export bothat, morphogradient, morpholaplace, - + MedialAxisTransform, + skeletonize, thinning, GuoAlgo diff --git a/src/skeletonization.jl b/src/skeletonization.jl new file mode 100644 index 0000000..464e8a7 --- /dev/null +++ b/src/skeletonization.jl @@ -0,0 +1,130 @@ +# checks if pixel is an articulation/critical node +function compute_critical_indices(num_configurations::Integer) + critical_index_array = Array{Bool, 1}() + for index = 0:num_configurations - 1 + append!(critical_index_array, maximum(find_components(index)) != maximum(find_components(index & ~2^4))) + end + critical_index_array +end + +function find_components(index::Integer) + label_components(reshape(digits(index, base=2, pad=9), (3, 3)), trues(3, 3)) +end + +# computes degree of cornerness of pixels +function corner_table_lookup(img::AbstractArray{Bool, 2}, corner_table::AbstractArray{T, 1}) where T <: Integer + corner_score = zeros(Integer, size(img) .+ 2) + score_table = [256 128 64; 32 16 8; 4 2 1] + + for i = 2:size(img)[1] + 1 + for j = 2:size(img)[2] + 1 + if img[i - 1, j - 1] + corner_score[i - 1:i + 1, j - 1:j + 1] = corner_score[i - 1:i + 1, j - 1:j + 1] + score_table + end + end + end + corner_table[corner_score[2: size(img)[1] + 1, 2: size(img)[2] + 1] .+ 1] +end + +# traverse foreground and reduce image to skeleton +function inner_skeleton_loop(img::AbstractArray{Bool, 2}, order::AbstractArray{T, 1}, table::AbstractArray{Bool, 1}) where T <: Integer + index = findall(img)[order] + score_table = [256 128 64; 32 16 8; 4 2 1] + result = falses(size(img)) + padded_image = padarray(img, Fill(0, (1, 1), (1, 1))) + for i in index + result[i] = table[sum(padded_image[i[1] - 1:i[1] + 1, i[2] - 1:i[2] + 1] .* score_table) + 1] + padded_image[i[1], i[2]] = result[i] + end + result +end + +""" +``` +skeletonize(img, MedialAxisTransform()) +``` + +Returns Medial axis of binary image that matches the dimensions of the input binary image. + +# Details +Skeletonization is a process for reducing foreground regions in a binary image to a +skeletal remnant that largely preserves the extent and connectivity of the original +region while throwing away most of the original foreground pixels. + +Steps of the algorithm : +- A look up table is built to check whether a pixel is to be removed or not. The pixel +is removed if it has more than one neighbour and if removing it doesn't change the connectedness. +- Distance transform is computed as well as degree of cornerness. +- Foreground pixels are traversed as per distance transform and cornerness. + +# Arguments +The img parameter needs to be a two-dimensional array of bool type where true represents +foreground and false represents background. + +# Example + +Compute the skeleton of rectangular object. +```julia + +using Images, ImageMorphology + +img = Bool.([0 0 0 0 0 + 0 1 1 1 0 + 0 1 1 1 0 + 0 1 1 1 0 + 0 1 1 1 0 + 0 0 0 0 0]) + +skeleton = skeletonize(img, MedialAxisTransform()) + +skeleton = [0 0 0 0 0 + 0 1 0 1 0 + 0 0 1 0 0 + 0 0 1 0 0 + 0 1 0 1 0 + 0 0 0 0 0] +``` + +# References +[1] http://homepages.inf.ed.ac.uk/rbf/HIPR2/skeleton.htm +[2] H. Blum in 1967, and in 1968 by L. Calabi and W. E. Hartnett. +""" +function skeletonize(img::AbstractArray{Bool, 2}, algo::MedialAxisTransform) + + # check if array has only 0s + if all(.~img) + return img + end + + # check if array has only 1s + if sum(img) == length(img) + img[1, 1:end] .= false + img[1:end - 1, end] .= false + img[end, 1:end - 1] .= false + return .~img + end + + # possible number of configurations of 3x3 kernel + num_configurations = 512 + + # build look-up table + patterns = [0:1:num_configurations - 1;] + num_pixels_in_patterns_kernel = sum.(digits.(patterns, base=2)) + center_is_foreground = patterns .& 2^4 .> 0 + table = (center_is_foreground .& compute_critical_indices(num_configurations) .| (num_pixels_in_patterns_kernel .< 3)) + + # store distance transform + distance = distance_transform(feature_transform(.~img)) + + # corners handling + corner_table = 9 .- num_pixels_in_patterns_kernel + corner_score = corner_table_lookup(img, corner_table) + + # generate order for traversing the shape + dist_corner_pair = CartesianIndex.(round.(Integer, distance[img] .* 10), corner_score[img]) + first_sort = sortperm(dist_corner_pair, by=x->x[2]) + second_sort = sortperm(dist_corner_pair[first_sort], by=x->x[1]) + traversal_order = first_sort[second_sort] + + inner_skeleton_loop(img, traversal_order, table) +end diff --git a/test/runtests.jl b/test/runtests.jl index 1fe2d75..7c43258 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -135,4 +135,20 @@ using Test 0 0 0 0 0 0]) @test thinning(img) == img end + @testset "Skeletonization" begin + img = Bool.([0 0 0 0 0 + 0 1 1 1 0 + 0 1 1 1 0 + 0 1 1 1 0 + 0 1 1 1 0 + 0 0 0 0 0]) + ans = Bool.([0 0 0 0 0 + 0 1 0 1 0 + 0 0 1 0 0 + 0 0 1 0 0 + 0 1 0 1 0 + 0 0 0 0 0]) + + @test skeletonize(img, MedialAxisTransform()) == ans + end end