From 90cede0ddd729645c34d988ab3c48c2f240473c7 Mon Sep 17 00:00:00 2001 From: Deepank308 Date: Mon, 15 Apr 2019 03:12:08 +0530 Subject: [PATCH 1/2] Add skeletonization function --- REQUIRE | 3 +- src/ImageMorphology.jl | 10 ++++-- src/skeletonization.jl | 81 ++++++++++++++++++++++++++++++++++++++++++ test/runtests.jl | 16 +++++++++ 4 files changed, 107 insertions(+), 3 deletions(-) create mode 100644 src/skeletonization.jl diff --git a/REQUIRE b/REQUIRE index c8e8cf8..9812907 100644 --- a/REQUIRE +++ b/REQUIRE @@ -1,2 +1,3 @@ julia 0.7 -ImageCore 0.7.0 +Images 0.13.0 +ImageCore 0.7.0 diff --git a/src/ImageMorphology.jl b/src/ImageMorphology.jl index d255977..322c103 100644 --- a/src/ImageMorphology.jl +++ b/src/ImageMorphology.jl @@ -3,10 +3,15 @@ __precompile__() module ImageMorphology using ImageCore +using Images + +abstract type SkeletonizationAlgo end +struct MedialAxis <: SkeletonizationAlgo end + include("dilation_and_erosion.jl") include("thinning.jl") - +include("skeletonization.jl") export dilate, erode, @@ -16,7 +21,8 @@ export bothat, morphogradient, morpholaplace, - + MedialAxis, + skeletonize, thinning, GuoAlgo diff --git a/src/skeletonization.jl b/src/skeletonization.jl new file mode 100644 index 0000000..c027c0e --- /dev/null +++ b/src/skeletonization.jl @@ -0,0 +1,81 @@ +# checks if element is essential for connectivity +function compute_essential_indices() + essential_index_array = Array{Bool, 1}() + for index = 0:511 + append!(essential_index_array, maximum(find_components(index)) != maximum(find_components(index & ~2^4))) + end + essential_index_array +end + +function find_components(index::Int64) + label_components(reshape(digits(index, base=2, pad=9), (3, 3)), trues(3, 3)) +end + +function corner_table_lookup(img::AbstractArray{Bool, 2}, corner_table::AbstractArray{Int, 1}) + corner_score = zeros(Int, 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 + +function inner_skeleton_loop(img::AbstractArray{Bool, 2}, order::Array{Int, 1}, table::AbstractArray{Bool, 1}) + index = findall(img)[order] + score_table = [256 128 64; 32 16 8; 4 2 1] + result = zeros(Int, size(img)) + padded_image = zeros(Int, size(img) .+ 2) + padded_image[2: size(img)[1] + 1, 2: size(img)[2] + 1] = img + for i in index + result[i] = table[sum(padded_image[i[1]:i[1] + 2, i[2]:i[2] + 2] .* score_table) + 1] + padded_image[i[1] + 1, i[2] + 1] = result[i] + end + result +end + +function skeletonize(img::AbstractArray{Bool}, algo::SkeletonizationAlgo) + skeletonize_impl(img, algo) +end + +function skeletonize_impl(img::AbstractArray{Bool, 2}, algo::MedialAxis) + + # check if array has only 0s + if sum(img) == 0 + 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 convert(Array{Bool}, .~img) + end + + # build look-up table + patterns = [0:1:511;] + num_pixels_in_patterns_kernel = sum.(digits.(patterns, base=2)) + center_is_foreground = convert(Array{Bool, 1}, patterns .& 2^4 .> 0) + table = (center_is_foreground .& compute_essential_indices() .| (num_pixels_in_patterns_kernel .< 3)) + table = convert(Array{Bool, 1}, table) + + # 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.(Int, 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]) + order = first_sort[second_sort] + + convert(Array{Bool, 2}, inner_skeleton_loop(img, order, table)) +end diff --git a/test/runtests.jl b/test/runtests.jl index 1fe2d75..5cea871 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, MedialAxis()) == ans + end end From d4f8638bed04f073d3ae3bf53b3b65b2a42f47b1 Mon Sep 17 00:00:00 2001 From: Deepank308 Date: Tue, 16 Apr 2019 00:12:56 +0530 Subject: [PATCH 2/2] Cleans up and add doc-string --- src/ImageMorphology.jl | 6 +- src/skeletonization.jl | 149 +++++++++++++++++++++++++++-------------- test/runtests.jl | 2 +- 3 files changed, 103 insertions(+), 54 deletions(-) diff --git a/src/ImageMorphology.jl b/src/ImageMorphology.jl index 322c103..cadaf7c 100644 --- a/src/ImageMorphology.jl +++ b/src/ImageMorphology.jl @@ -5,8 +5,8 @@ module ImageMorphology using ImageCore using Images -abstract type SkeletonizationAlgo end -struct MedialAxis <: SkeletonizationAlgo end +abstract type AbstractSkeletonizationAlgorithm end +struct MedialAxisTransform <: AbstractSkeletonizationAlgorithm end include("dilation_and_erosion.jl") @@ -21,7 +21,7 @@ export bothat, morphogradient, morpholaplace, - MedialAxis, + MedialAxisTransform, skeletonize, thinning, GuoAlgo diff --git a/src/skeletonization.jl b/src/skeletonization.jl index c027c0e..464e8a7 100644 --- a/src/skeletonization.jl +++ b/src/skeletonization.jl @@ -1,81 +1,130 @@ -# checks if element is essential for connectivity -function compute_essential_indices() - essential_index_array = Array{Bool, 1}() - for index = 0:511 - append!(essential_index_array, maximum(find_components(index)) != maximum(find_components(index & ~2^4))) +# 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 - essential_index_array + critical_index_array end -function find_components(index::Int64) +function find_components(index::Integer) label_components(reshape(digits(index, base=2, pad=9), (3, 3)), trues(3, 3)) end -function corner_table_lookup(img::AbstractArray{Bool, 2}, corner_table::AbstractArray{Int, 1}) - corner_score = zeros(Int, size(img) .+ 2) +# 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 end corner_table[corner_score[2: size(img)[1] + 1, 2: size(img)[2] + 1] .+ 1] end -function inner_skeleton_loop(img::AbstractArray{Bool, 2}, order::Array{Int, 1}, table::AbstractArray{Bool, 1}) - index = findall(img)[order] - score_table = [256 128 64; 32 16 8; 4 2 1] - result = zeros(Int, size(img)) - padded_image = zeros(Int, size(img) .+ 2) - padded_image[2: size(img)[1] + 1, 2: size(img)[2] + 1] = img - for i in index - result[i] = table[sum(padded_image[i[1]:i[1] + 2, i[2]:i[2] + 2] .* score_table) + 1] - padded_image[i[1] + 1, i[2] + 1] = result[i] - end - result -end - -function skeletonize(img::AbstractArray{Bool}, algo::SkeletonizationAlgo) - skeletonize_impl(img, algo) +# 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 -function skeletonize_impl(img::AbstractArray{Bool, 2}, algo::MedialAxis) +""" +``` +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 sum(img) == 0 + 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 convert(Array{Bool}, .~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:511;] + patterns = [0:1:num_configurations - 1;] num_pixels_in_patterns_kernel = sum.(digits.(patterns, base=2)) - center_is_foreground = convert(Array{Bool, 1}, patterns .& 2^4 .> 0) - table = (center_is_foreground .& compute_essential_indices() .| (num_pixels_in_patterns_kernel .< 3)) - table = convert(Array{Bool, 1}, table) - - # 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.(Int, 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]) - order = first_sort[second_sort] - - convert(Array{Bool, 2}, inner_skeleton_loop(img, order, table)) + 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 5cea871..7c43258 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -149,6 +149,6 @@ using Test 0 1 0 1 0 0 0 0 0 0]) - @test skeletonize(img, MedialAxis()) == ans + @test skeletonize(img, MedialAxisTransform()) == ans end end