diff --git a/REQUIRE b/REQUIRE index 80a6fdd..1abc461 100644 --- a/REQUIRE +++ b/REQUIRE @@ -3,6 +3,3 @@ Images 0.5 Colors 0.6 ColorVectorSpace 0.1 FixedPointNumbers 0.1 -FileIO -Compat 0.7.15 -StatsBase diff --git a/src/ImageFeatures.jl b/src/ImageFeatures.jl index 3f50310..2b99262 100644 --- a/src/ImageFeatures.jl +++ b/src/ImageFeatures.jl @@ -1,5 +1,4 @@ module ImageFeatures - # package code goes here using Images, ColorTypes, FixedPointNumbers, Distributions @@ -12,8 +11,9 @@ include("brief.jl") include("orb.jl") include("freak.jl") include("brisk.jl") +include("censure.jl") -export Keypoint, Keypoints, Feature, Features, Params, BRIEF, ORB, FREAK, BRISK +export Keypoint, Keypoints, BRIEF, DescriptorParams, ORB, CENSURE, OctagonFilter, BoxFilter, BiFilter, Detector export #Core @@ -56,5 +56,8 @@ export random_coarse, gaussian, gaussian_local, - centered + centered, + + #CENSURE + censure end diff --git a/src/censure.jl b/src/censure.jl new file mode 100644 index 0000000..6e0fd7f --- /dev/null +++ b/src/censure.jl @@ -0,0 +1,246 @@ +abstract BiFilter + +type BoxFilter <: BiFilter + + scale :: Int + in_length :: Int + out_length :: Int + in_area :: Float64 + out_area :: Float64 + in_weight :: Float64 + out_weight :: Float64 + +end + +type OctagonFilter <: BiFilter + + m_out :: Int + m_in :: Int + n_out :: Int + n_in :: Int + in_area :: Float64 + out_area :: Float64 + in_weight :: Float64 + out_weight :: Float64 + +end + +type CENSURE{F} <: Detector + + smallest :: Int + largest :: Int + filter_type :: Type{F} + filter_stack :: Array{F} + response_threshold :: Float64 + line_threshold :: Float64 + +end + +function OctagonFilter(mo, mi, no, ni) + OF = OctagonFilter(mo, mi, no, ni, 0.0, 0.0, 0.0, 0.0) + OF.out_area = OF.m_out ^ 2 + 2 * OF.n_out ^ 2 + 4 * OF.m_out * OF.n_out + OF.in_area = OF.m_in ^ 2 + 2 * OF.n_in ^ 2 + 4 * OF.m_in * OF.n_in + OF.out_weight = 1.0 / (OF.out_area - OF.in_area) + OF.in_weight = 1.0 / OF.in_area + OF +end + +function BoxFilter(s) + BF = BoxFilter(s, 2 * s + 1, 4 * s + 1, (2 * s + 1) ^ 2, (4 * s + 1) ^ 2, 0.0, 0.0) + BF.in_weight = 1.0 / BF.in_area + BF.out_weight = 1.0 / (BF.out_area - BF.in_area) + BF +end + +const octagon_filter_kernels = [[5, 3, 2, 0], + [5, 3, 3, 1], + [7, 3, 3, 2], + [9, 5, 4, 2], + [9, 5, 7, 3], + [13, 5, 7, 4], + [15, 5, 10, 5]] + +const box_filter_kernels = [1, 2, 3, 4, 5, 6, 7] + +_getkernel(::Type{BoxFilter}) = box_filter_kernels +_getkernel(::Type{OctagonFilter}) = octagon_filter_kernels + +function _get_filter_stack(filter_type::Type, smallest::Integer, largest::Integer) + k = _getkernel(filter_type) + filter_stack = map(f -> filter_type(f...), k[smallest : largest]) +end + +_get_integral_image(img, filter_type::BoxFilter) = integral_image(img) + +function _get_integral_image(img, filter_type::OctagonFilter) + img_shape = size(img) + int_img = zeros(img_shape) + right_slant_img = zeros(img_shape) + left_slant_img = zeros(img_shape) + + int_img[1, :] = cumsum(img[1, :]) + right_slant_img[1, :] = int_img[1, :] + left_slant_img[1, :] = int_img[1, :] + + for i in 2:img_shape[1] + sum = 0.0 + for j in 1:img_shape[2] + sum += img[i, j] + int_img[i, j] = sum + int_img[i - 1, j] + left_slant_img[i, j] = sum + right_slant_img[i, j] = sum + + if j > 1 left_slant_img[i, j] += left_slant_img[i - 1, j - 1] end + right_slant_img[i, j] += j < img_shape[2] ? right_slant_img[i - 1, j + 1] : right_slant_img[i - 1, j] + end + end + int_img, right_slant_img, left_slant_img +end + +function _filter_response{T}(int_img::AbstractArray{T, 2}, BF::BoxFilter) + margin = BF.scale * 2 + n = BF.scale + img_shape = size(int_img) + response = zeros(T, img_shape) + R = CartesianRange(CartesianIndex((margin + 2, margin + 2)), CartesianIndex((img_shape[1] - margin, img_shape[2] - margin))) + + for I in R + topleft = I + CartesianIndex(- n - 1, - n - 1) + topright = I + CartesianIndex(- n - 1, n) + bottomleft = I + CartesianIndex(n, - n - 1) + bottomright = I + CartesianIndex(n, n) + A = checkbounds(Bool, int_img, topleft) ? int_img[topleft] : zero(T) + B = checkbounds(Bool, int_img, topright) ? int_img[topright] : zero(T) + C = checkbounds(Bool, int_img, bottomleft) ? int_img[bottomleft] : zero(T) + D = checkbounds(Bool, int_img, bottomright) ? int_img[bottomright] : zero(T) + in_sum = A + D - B - C + + topleft = I + CartesianIndex(- 2 * n - 1, - 2 * n - 1) + topright = I + CartesianIndex(- 2 * n - 1, 2 * n) + bottomleft = I + CartesianIndex(2 * n, - 2 * n - 1) + bottomright = I + CartesianIndex(2 * n, 2 * n) + A = checkbounds(Bool, int_img, topleft) ? int_img[topleft] : zero(T) + B = checkbounds(Bool, int_img, topright) ? int_img[topright] : zero(T) + C = checkbounds(Bool, int_img, bottomleft) ? int_img[bottomleft] : zero(T) + D = checkbounds(Bool, int_img, bottomright) ? int_img[bottomright] : zero(T) + out_sum = A + D - B - C - in_sum + response[I] = out_sum * BF.out_weight - BF.in_weight * in_sum + end + + response +end + +function _filter_response(int_imgs::Tuple, OF::OctagonFilter) + int_img = int_imgs[1] + rs_img = int_imgs[2] + ls_img = int_imgs[3] + + T = eltype(int_img) + + margin = Int(floor(OF.m_out / 2 + OF.n_out)) + m_in2 = Int(floor(OF.m_in / 2)) + m_out2 = Int(floor(OF.m_out / 2)) + + img_shape = size(int_img) + response = zeros(T, img_shape) + R = CartesianRange(CartesianIndex((margin + 2, margin + 2)), CartesianIndex((img_shape[1] - margin, img_shape[2] - margin))) + + for I in R + topleft = I + CartesianIndex(- m_in2 - 1, - m_in2 - OF.n_in - 1) + topright = I + CartesianIndex(m_in2, - m_in2 - OF.n_in - 1) + bottomleft = I + CartesianIndex(- m_in2 - 1, m_in2 + OF.n_in) + bottomright = I + CartesianIndex(m_in2, m_in2 + OF.n_in) + A = checkbounds(Bool, int_img, topleft) ? int_img[topleft] : zero(T) + B = checkbounds(Bool, int_img, topright) ? int_img[topright] : zero(T) + C = checkbounds(Bool, int_img, bottomleft) ? int_img[bottomleft] : zero(T) + D = checkbounds(Bool, int_img, bottomright) ? int_img[bottomright] : zero(T) + in_sum = A + D - B - C + + trap_top_right = bottomright + trap_bot_right = I + CartesianIndex(m_in2 + OF.n_in, m_in2) + trap_top_left = I + CartesianIndex(m_in2, - m_in2 - OF.n_in) + trap_bot_left = I + CartesianIndex(m_in2 + OF.n_in, - m_in2 - 1) + A = checkbounds(Bool, ls_img, trap_top_left) ? ls_img[trap_top_left] : zero(T) + B = checkbounds(Bool, rs_img, trap_top_right) ? rs_img[trap_top_right] : zero(T) + C = checkbounds(Bool, ls_img, trap_bot_left) ? ls_img[trap_bot_left] : zero(T) + D = checkbounds(Bool, rs_img, trap_bot_right) ? rs_img[trap_bot_right] : zero(T) + in_sum += A + D - B - C + + trap_top_right = I + CartesianIndex(- m_in2 - OF.n_in - 1, m_in2 - 1) + trap_top_left = I + CartesianIndex(- m_in2 - OF.n_in - 1, - m_in2) + trap_bot_right = I + CartesianIndex(- m_in2 - 1, m_in2 + OF.n_in - 1) + trap_bot_left = I + CartesianIndex(- m_in2 - 1, - m_in2 - OF.n_in) + A = checkbounds(Bool, rs_img, trap_top_left) ? rs_img[trap_top_left] : zero(T) + B = checkbounds(Bool, ls_img, trap_top_right) ? ls_img[trap_top_right] : zero(T) + C = checkbounds(Bool, rs_img, trap_bot_left) ? rs_img[trap_bot_left] : zero(T) + D = checkbounds(Bool, ls_img, trap_bot_right) ? ls_img[trap_bot_right] : zero(T) + in_sum += A + D - B - C + + topleft = I + CartesianIndex(- m_out2 - 1, - m_out2 - OF.n_out - 1) + topright = I + CartesianIndex(m_out2, - m_out2 - OF.n_out - 1) + bottomleft = I + CartesianIndex(- m_out2 - 1, m_out2 + OF.n_out) + bottomright = I + CartesianIndex(m_out2, m_out2 + OF.n_out) + A = checkbounds(Bool, int_img, topleft) ? int_img[topleft] : zero(T) + B = checkbounds(Bool, int_img, topright) ? int_img[topright] : zero(T) + C = checkbounds(Bool, int_img, bottomleft) ? int_img[bottomleft] : zero(T) + D = checkbounds(Bool, int_img, bottomright) ? int_img[bottomright] : zero(T) + out_sum = A + D - B - C + + trap_top_right = bottomright + trap_bot_right = I + CartesianIndex(m_out2 + OF.n_out, m_out2) + trap_top_left = I + CartesianIndex(m_out2, - m_out2 - OF.n_out) + trap_bot_left = I + CartesianIndex(m_out2 + OF.n_out, - m_out2 - 1) + A = checkbounds(Bool, ls_img, trap_top_left) ? ls_img[trap_top_left] : zero(T) + B = checkbounds(Bool, rs_img, trap_top_right) ? rs_img[trap_top_right] : zero(T) + C = checkbounds(Bool, ls_img, trap_bot_left) ? ls_img[trap_bot_left] : zero(T) + D = checkbounds(Bool, rs_img, trap_bot_right) ? rs_img[trap_bot_right] : zero(T) + out_sum += A + D - B - C + + trap_top_right = I + CartesianIndex(- m_out2 - OF.n_out - 1, m_out2 - 1) + trap_top_left = I + CartesianIndex(- m_out2 - OF.n_out - 1, - m_out2) + trap_bot_right = I + CartesianIndex(- m_out2 - 1, m_out2 + OF.n_out - 1) + trap_bot_left = I + CartesianIndex(- m_out2 - 1, - m_out2 - OF.n_out) + A = checkbounds(Bool, rs_img, trap_top_left) ? rs_img[trap_top_left] : zero(T) + B = checkbounds(Bool, ls_img, trap_top_right) ? ls_img[trap_top_right] : zero(T) + C = checkbounds(Bool, rs_img, trap_bot_left) ? rs_img[trap_bot_left] : zero(T) + D = checkbounds(Bool, ls_img, trap_bot_right) ? ls_img[trap_bot_right] : zero(T) + out_sum += A + D - B - C + out_sum = out_sum - in_sum + + response[I] = in_sum * OF.in_weight - out_sum * OF.out_weight + end + response + +end + +function CENSURE(; smallest::Integer = 1, largest::Integer = 7, filter::Type = BoxFilter, response_threshold::Number = 0.15, line_threshold::Number = 10) + CENSURE{filter}(smallest, largest, filter, _get_filter_stack(filter, smallest, largest), response_threshold, line_threshold) +end + +function censure{T, F}(img::AbstractArray{T, 2}, params::CENSURE{F}) + int_img = _get_integral_image(img, params.filter_stack[1]) + responses = map(f -> _filter_response(int_img, f), params.filter_stack) + response_matrix = reshape(hcat(responses...), size(img)..., size(responses)...) + minima, maxima = extrema_filter(convert(Array{Float64}, padarray(response_matrix, [1, 1, 1], [1, 1, 1], "replicate")), [3, 3, 3]) + features = map(i -> (minima[i] == response_matrix[i] || maxima[i] == response_matrix[i]) && ( response_matrix[i] > params.response_threshold ), CartesianRange(size(response_matrix))) + (grad_x, grad_y) = imgradients(img, "sobel", "replicate") + cov_xx = grad_x .* grad_x + cov_xy = grad_x .* grad_y + cov_yy = grad_y .* grad_y + for i in 1:params.largest - params.smallest + 1 + gamma = (1 + (params.smallest + i - 1) / 3.0) + filt_cov_xx = imfilter_gaussian(cov_xx, [gamma, gamma]) + filt_cov_xy = imfilter_gaussian(cov_xy, [gamma, gamma]) + filt_cov_yy = imfilter_gaussian(cov_yy, [gamma, gamma]) + + features[:, :, i] = map((xx, yy, xy, f) -> (xx + yy) ^ 2 > params.line_threshold * (xx * yy - xy ^ 2) ? false : f, filt_cov_xx, filt_cov_yy, filt_cov_xy, features[:, :, i]) + end + keypoints = Array{Keypoint}([]) + scales = Array{Integer}([]) + for scale in 1:(params.largest - params.smallest + 1) + rows, cols, _ = findnz(features[:, :, scale]) + append!(keypoints, map((r, c) -> Keypoint(r, c), rows, cols)) + append!(scales, ones(length(rows)) * scale) + end + keypoints, scales +end \ No newline at end of file diff --git a/test/censure.jl b/test/censure.jl new file mode 100644 index 0000000..09c0e94 --- /dev/null +++ b/test/censure.jl @@ -0,0 +1,110 @@ +using FactCheck, Base.Test, Images, Colors, FixedPointNumbers, ImageFeatures + +facts("CENSURE") do + + context("Filters") do + + bf = BoxFilter(1) + @fact bf.scale --> 1 + @fact bf.in_length --> 3 + @fact bf.out_length --> 5 + @fact bf.in_area --> 9.0 + @fact bf.out_area --> 25.0 + @fact isapprox(bf.in_weight, 1 / 9) --> true + @fact isapprox(bf.out_weight, 1 / 16) --> true + bf = BoxFilter(5) + @fact bf.scale --> 5 + @fact bf.in_length --> 11 + @fact bf.out_length --> 21 + @fact bf.in_area --> 121.0 + @fact bf.out_area --> 441.0 + @fact isapprox(bf.in_weight, 1 / 121) --> true + @fact isapprox(bf.out_weight, 1 / 320) --> true + + of = OctagonFilter(5, 3, 2, 0) + @fact of.m_out --> 5 + @fact of.m_in --> 3 + @fact of.n_out --> 2 + @fact of.n_in --> 0 + @fact of.in_area --> 9.0 + @fact of.out_area --> 73.0 + @fact isapprox(of.in_weight, 1 / 9) --> true + @fact isapprox(of.out_weight, 1 / 64) --> true + + of = OctagonFilter(13, 5, 7, 4) + @fact of.m_out --> 13 + @fact of.m_in --> 5 + @fact of.n_out --> 7 + @fact of.n_in --> 4 + @fact of.in_area --> 137.0 + @fact of.out_area --> 631.0 + @fact isapprox(of.in_weight, 1 / 137) --> true + @fact isapprox(of.out_weight, 1 / 494) --> true + + end + + context("Integral Image") do + + img = ones(5, 5) + bf = BoxFilter(1) + @fact all(integral_image(img) .== ImageFeatures._get_integral_image(img, bf)) --> true + + of = OctagonFilter(5, 3, 2, 0) + i, rs, ls = ImageFeatures._get_integral_image(img, of) + @fact all(i .== integral_image(img)) --> true + r_check = [ 1.0 2.0 3.0 4.0 5.0 + 3.0 5.0 7.0 9.0 10.0 + 6.0 9.0 12.0 14.0 15.0 + 10.0 14.0 17.0 19.0 20.0 + 15.0 19.0 22.0 24.0 25.0 ] + @fact all(rs .== r_check) --> true + l_check = [ 1.0 2.0 3.0 4.0 5.0 + 1.0 3.0 5.0 7.0 9.0 + 1.0 3.0 6.0 9.0 12.0 + 1.0 3.0 6.0 10.0 14.0 + 1.0 3.0 6.0 10.0 15.0 ] + @fact all(ls .== l_check) --> true + + end + + context("Filter Response") do + + img = [4.0 2.0 6.0 1.0 1.0 8.0 2.0; + 9.0 3.0 8.0 3.0 5.0 8.0 4.0; + 4.0 10.0 6.0 2.0 5.0 1.0 5.0; + 10.0 5.0 8.0 6.0 6.0 3.0 3.0; + 5.0 6.0 4.0 8.0 3.0 3.0 9.0; + 3.0 1.0 5.0 6.0 2.0 2.0 2.0; + 9.0 3.0 4.0 1.0 10.0 8.0 6.0] + bf = BoxFilter(1) + response = ImageFeatures._filter_response(ImageFeatures._get_integral_image(img, bf), bf) + @fact isapprox(response[4, 4], -0.895833, rtol = 0.001) --> true + @fact isapprox(response[4, 5], 0.888889, rtol = 0.001) --> true + @fact isapprox(response[5, 4], -0.958333, rtol = 0.001) --> true + @fact isapprox(response[5, 5], 0.604167, rtol = 0.001) --> true + @fact all(response[:, 1:3] .== 0) --> true + @fact all(response[1:3, :] .== 0) --> true + @fact all(response[:, 6:6] .== 0) --> true + @fact all(response[6:7, :] .== 0) --> true + + img = [ 3.0 5.0 8.0 8.0 5.0 4.0 1.0 + 3.0 8.0 6.0 10.0 1.0 5.0 4.0 + 10.0 1.0 6.0 4.0 10.0 4.0 5.0 + 2.0 10.0 9.0 4.0 5.0 3.0 7.0 + 1.0 7.0 5.0 9.0 7.0 6.0 3.0 + 5.0 6.0 9.0 9.0 1.0 4.0 8.0 + 2.0 4.0 6.0 9.0 8.0 4.0 10.0] + response = ImageFeatures._filter_response(ImageFeatures._get_integral_image(img, bf), bf) + @fact isapprox(response[4, 4], -0.93056, rtol = 0.001) --> true + @fact isapprox(response[4, 5], -0.02777, rtol = 0.001) --> true + @fact isapprox(response[5, 4], -0.69444, rtol = 0.001) --> true + @fact isapprox(response[5, 5], 1.35417, rtol = 0.001) --> true + @fact all(response[:, 1:3] .== 0) --> true + @fact all(response[1:3, :] .== 0) --> true + @fact all(response[:, 6:6] .== 0) --> true + @fact all(response[6:7, :] .== 0) --> true + + + end + +end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 643fc42..b7be5f6 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -55,6 +55,7 @@ include("corner.jl") include("orb.jl") include("freak.jl") include("brisk.jl") +include("censure.jl") isinteractive() || FactCheck.exitstatus()