diff --git a/Project.toml b/Project.toml index 214974ee..e43fbba4 100644 --- a/Project.toml +++ b/Project.toml @@ -10,6 +10,7 @@ ImageAxes = "2803e5a7-5153-5ecf-9a86-9b4c37f5f5ac" ImageBase = "c817782e-172a-44cc-b673-b171935fbb9e" ImageContrastAdjustment = "f332f351-ec65-5f6a-b3d1-319c6670881a" ImageCore = "a09fc81d-aa75-5fe9-8630-4744c3626534" +ImageCorners = "89d5987c-236e-4e32-acd0-25bd6bd87b70" ImageDistances = "51556ac3-7006-55f5-8cb3-34580c88182d" ImageFiltering = "6a3955dd-da59-5b1f-98d4-e7296123deb5" ImageIO = "82e4d734-157c-48bb-816b-45c225c6df19" @@ -37,6 +38,7 @@ ImageAxes = "0.6.9" ImageBase = "0.1.5" ImageContrastAdjustment = "0.3.3" ImageCore = "0.9.3" +ImageCorners = "0.1.2" ImageDistances = "0.2.5" ImageFiltering = "0.7.1" ImageIO = "0.3, 0.4, 0.5, 0.6" diff --git a/src/Images.jl b/src/Images.jl index 80ba944a..139ecfbc 100644 --- a/src/Images.jl +++ b/src/Images.jl @@ -19,6 +19,7 @@ const is_little_endian = ENDIAN_BOM == 0x04030201 # CHECKME(johnnychen94): is th @reexport using ImageTransformations @reexport using ImageAxes +@reexport using ImageCorners @reexport using ImageMetadata @reexport using ImageFiltering @reexport using ImageMorphology @@ -46,31 +47,19 @@ import FileIO: metadata import Graphics: Point include("compat.jl") -include("misc.jl") include("labeledarrays.jl") include("algorithms.jl") include("deprecations.jl") -include("corner.jl") include("edge.jl") export # types ColorizedArray, - Percentile, # macros @test_approx_eq_sigma_eps, # algorithms - imcorner, - imcorner_subpixel, - corner2subpixel, - harris, - shi_tomasi, - kitchen_rosenfeld, - fastcorners, - meancovs, - gammacovs, imedge, # TODO: deprecate? imgaussiannoise, otsu_threshold, diff --git a/src/corner.jl b/src/corner.jl deleted file mode 100644 index 45aa9de7..00000000 --- a/src/corner.jl +++ /dev/null @@ -1,325 +0,0 @@ -""" -``` -corners = imcorner(img; [method]) -corners = imcorner(img, threshold; [method]) -``` - -Performs corner detection using one of the following methods - - - 1. harris - 2. shi_tomasi - 3. kitchen_rosenfeld - -The parameters of the individual methods are described in their documentation. The -maxima values of the resultant responses are taken as corners. If a threshold is -specified, the values of the responses are thresholded to give the corner pixels. -If `threshold` is a `Percentile` then its type will be preserved. -""" -function imcorner(img::AbstractArray; method::Function = harris, args...) - responses = method(img; args...) - corners = similar(img, Bool) - fill!(corners, false) - maxima = map(CartesianIndex, findlocalmaxima(responses)) - for m in maxima corners[m] = true end - corners -end - -function imcorner(img::AbstractArray, threshold; method::Function = harris, args...) - responses = method(img; args...) - map(i -> i > threshold, responses) -end - -function imcorner(img::AbstractArray, thresholdp::Percentile; method::Function = harris, args...) - responses = method(img; args...) - threshold = StatsBase.percentile(vec(responses), thresholdp.p) - map(i -> i > threshold, responses) -end - -""" -``` -corners = imcorner_subpixel(img; [method]) - -> Vector{HomogeneousPoint{Float64,3}} -corners = imcorner_subpixel(img, threshold, percentile; [method]) - -> Vector{HomogeneousPoint{Float64,3}} -``` - -Same as [`imcorner`](@ref), but estimates corners to sub-pixel precision. - -Sub-pixel precision is achieved by interpolating the corner response values using -the 4-connected neighbourhood of a maximum response value. -See [`corner2subpixel`](@ref) for more details of the interpolation scheme. - -""" -function imcorner_subpixel(img::AbstractArray; method::Function = harris, args...) - responses = method(img; args...) - corner_indicator = similar(img, Bool) - fill!(corner_indicator, false) - maxima = map(CartesianIndex, findlocalmaxima(responses)) - for m in maxima corner_indicator[m] = true end - corners = corner2subpixel(responses,corner_indicator) -end - -function imcorner_subpixel(img::AbstractArray, threshold; method::Function = harris, args...) - responses = method(img; args...) - corner_indicator = map(i -> i > threshold, responses) - corners = corner2subpixel(responses,corner_indicator) -end - -function imcorner_subpixel(img::AbstractArray, thresholdp::Percentile; method::Function = harris, args...) - responses = method(img; args...) - threshold = StatsBase.percentile(vec(responses), thresholdp.p) - corner_indicator = map(i -> i > threshold, responses) - corners = corner2subpixel(responses,corner_indicator) -end - - -""" -``` -corners = corner2subpixel(responses::AbstractMatrix,corner_indicator::AbstractMatrix{Bool}) - -> Vector{HomogeneousPoint{Float64,3}} -``` -Refines integer corner coordinates to sub-pixel precision. - -The function takes as input a matrix representing corner responses -and a boolean indicator matrix denoting the integer coordinates of a corner -in the image. The output is a vector of type [`HomogeneousPoint`](@ref) -storing the sub-pixel coordinates of the corners. - -The algorithm computes a correction factor which is added to the original -integer coordinates. In particular, a univariate quadratic polynomial is fit -separately to the ``x``-coordinates and ``y``-coordinates of a corner and its immediate -east/west, and north/south neighbours. The fit is achieved using a local -coordinate system for each corner, where the origin of the coordinate system is -a given corner, and its immediate neighbours are assigned coordinates of minus -one and plus one. - -The corner and its two neighbours form a system of three equations. For example, -let ``x_1 = -1``, ``x_2 = 0`` and ``x_3 = 1`` denote the local ``x`` coordinates -of the west, center and east pixels and let the vector ``\\mathbf{b} = [r_1, r_2, r_3]`` -denote the corresponding corner response values. With - -```math - \\mathbf{A} = - \\begin{bmatrix} - x_1^2 & x_1 & 1 \\\\ - x_2^2 & x_2 & 1 \\\\ - x_3^2 & x_3 & 1 \\\\ - \\end{bmatrix}, -``` -the coefficients of the quadratic polynomial can be found by solving the -system of equations ``\\mathbf{b} = \\mathbf{A}\\mathbf{x}``. -The result is given by ``x = \\mathbf{A}^{-1}\\mathbf{b}``. - -The vertex of the quadratic polynomial yields a sub-pixel estimate of the -true corner position. For example, for a univariate quadratic polynomial -``px^2 + qx + r``, the ``x``-coordinate of the vertex is ``\\frac{-q}{2p}``. -Hence, the refined sub-pixel coordinate is equal to: - ``c + \\frac{-q}{2p}``, where ``c`` is the integer coordinate. - -!!! note - Corners on the boundary of the image are not refined to sub-pixel precision. - -""" -function corner2subpixel(responses::AbstractMatrix, corner_indicator::AbstractMatrix{Bool}) - row_range, col_range = axes(corner_indicator) - idxs = findall(!iszero, corner_indicator) # findnz - row, col = (getindex.(idxs,1),getindex.(idxs,2)) - ncorners = length(row) - corners = fill(HomogeneousPoint((0.0,0.0,0.0)),ncorners) - invA = @SMatrix [0.5 -1.0 0.5; -0.5 0.0 0.5; 0.0 1.0 -0.0] - for k = 1:ncorners - # Corners on the perimeter of the image will not be interpolated. - if (row[k] == first(row_range) || row[k] == last(row_range) || - col[k] == first(col_range) || col[k] == last(col_range)) - y = convert(Float64,row[k]) - x = convert(Float64,col[k]) - corners[k] = HomogeneousPoint((x,y,1.0)) - else - center, north, south, east, west = - unsafe_neighbourhood_4(responses,row[k],col[k]) - # Solve for the coefficients of the quadratic equation. - a, b, c = invA* @SVector [west, center, east] - p, q, r = invA* @SVector [north, center, south] - # Solve for the first coordinate of the vertex. - u = -b/(2.0a) - v = -q/(2.0p) - corners[k] = HomogeneousPoint((col[k]+u,row[k]+v,1.0)) - end - end - return corners -end - -""" -``` -unsafe_neighbourhood_4(matrix::AbstractMatrix,r::Int,c::Int) -``` - -Returns the value of a matrix at given coordinates together with the values -of the north, south, east and west neighbours. - -This function does not perform bounds checking. It is up to the user to ensure -that the function is not called with indices that are on the boundary of the -matrix. - -""" -function unsafe_neighbourhood_4(matrix::AbstractMatrix,r::Int,c::Int) - center = matrix[r,c] - north = matrix[r-1,c] - south = matrix[r+1,c] - east = matrix[r,c+1] - west = matrix[r,c-1] - return center, north, south, east, west -end - - -""" -``` -harris_response = harris(img; [k], [border], [weights]) -``` - -Performs Harris corner detection. The covariances can be taken using either a mean -weighted filter or a gamma kernel. -""" -function harris(img::AbstractArray; k::Float64 = 0.04, args...) - cov_xx, cov_xy, cov_yy = gradcovs(img, args...) - corner = map((xx, yy, xy) -> xx * yy - xy ^ 2 - k * (xx + yy) ^ 2, cov_xx, cov_yy, cov_xy) - corner -end - -""" -``` -shi_tomasi_response = shi_tomasi(img; [border], [weights]) -``` - -Performs Shi Tomasi corner detection. The covariances can be taken using either a mean -weighted filter or a gamma kernel. -""" -function shi_tomasi(img::AbstractArray; border::AbstractString = "replicate", args...) - cov_xx, cov_xy, cov_yy = gradcovs(img, border; args...) - corner = map((xx, yy, xy) -> ((xx + yy) - (sqrt((xx - yy) ^ 2 + 4 * xy ^ 2))) / 2, cov_xx, cov_yy, cov_xy) - corner -end - -""" -``` -kitchen_rosenfeld_response = kitchen_rosenfeld(img; [border]) -``` - -Performs Kitchen Rosenfeld corner detection. The covariances can be taken using either a mean -weighted filter or a gamma kernel. -""" -function kitchen_rosenfeld(img::AbstractArray; border::AbstractString = "replicate") - meth = KernelFactors.sobel - (grad_x, grad_y) = imgradients(img, meth, border) - (grad_xx, grad_xy) = imgradients(grad_x, meth, border) - (grad_yx, grad_yy) = imgradients(grad_y, meth, border) - map(kr, grad_x, grad_y, grad_xx, grad_xy, grad_yy) -end - -function kr(x::T, y::T, xx::T, xy::T, yy::T) where T<:Real - num = xx*y*y + yy*x*x - 2*xy*x*y - denom = x*x + y*y - ifelse(denom == 0, zero(num)/one(denom), -num/denom) -end - -function kr(x::Real, y::Real, xx::Real, xy::Real, yy::Real) - xp, yp, xxp, xyp, yyp = promote(x, y, xx, xy, yy) - kr(xp, yp, xxp, xyp, yyp) -end - -kr(x::NumberLike, y::NumberLike, xx::NumberLike, xy::NumberLike, yy::NumberLike) = - kr(gray(x), gray(y), gray(xx), gray(xy), gray(yy)) - -function kr(x::AbstractRGB, y::AbstractRGB, xx::AbstractRGB, xy::AbstractRGB, yy::AbstractRGB) - krrgb = RGB(kr(red(x), red(y), red(xx), red(xy), red(yy)), - kr(green(x), green(y), green(xx), green(xy), green(yy)), - kr(blue(x), blue(y), blue(xx), blue(xy), blue(yy))) - gray(convert(Gray, krrgb)) -end - -""" - fastcorners(img, n, threshold) -> corners - -Performs FAST Corner Detection. `n` is the number of contiguous pixels -which need to be greater (lesser) than intensity + threshold (intensity - threshold) -for a pixel to be marked as a corner. The default value for n is 12. -""" -function fastcorners(img::AbstractArray{T}, n::Int = 12, threshold::Float64 = 0.15) where T - img_padded = padarray(img, Fill(0, (3,3))) - corner = falses(size(img)) - R = CartesianIndices(size(img)) - idx = map(CartesianIndex{2}, [(0, 3), (1, 3), (2, 2), (3, 1), (3, 0), (3, -1), (2, -2), (1, -3), - (0, -3), (-1, -3), (-2, -2), (-3, -1), (-3, 0), (-3, 1), (-2, 2), (-1, 3)]) - - idxidx = [1, 5, 9, 13] - for I in R - bright_threshold = img_padded[I] + threshold - dark_threshold = img_padded[I] - threshold - if n >= 12 - sum_bright = 0 - sum_dark = 0 - for k in idxidx - pixel = img_padded[I + idx[k]] - if pixel > bright_threshold - sum_bright += 1 - elseif pixel < dark_threshold - sum_dark += 1 - end - end - if sum_bright < 3 && sum_dark < 3 - continue - end - end - consecutive_bright = 0 - consecutive_dark = 0 - - for i in 1:15 + n - k = mod1(i, 16) - pixel = img_padded[I + idx[k]] - if pixel > bright_threshold - consecutive_dark = 0 - consecutive_bright += 1 - elseif pixel < dark_threshold - consecutive_bright = 0 - consecutive_dark += 1 - end - - if consecutive_dark == n || consecutive_bright == n - corner[I] = true - break - end - end - end - corner -end - -function gradcovs(img::AbstractArray, border::AbstractString = "replicate"; weights::Function = meancovs, args...) - (grad_x, grad_y) = imgradients(img, KernelFactors.sobel, border) - - cov_xx = dotc.(grad_x, grad_x) - cov_xy = dotc.(grad_x, grad_y) - cov_yy = dotc.(grad_y, grad_y) - - weights(cov_xx, cov_xy, cov_yy, args...) -end - -function meancovs(cov_xx, cov_xy, cov_yy, blockSize::Int = 3) - - box_filter_kernel = centered((1 / (blockSize * blockSize)) * ones(blockSize, blockSize)) - - filt_cov_xx = imfilter(cov_xx, box_filter_kernel) - filt_cov_xy = imfilter(cov_xy, box_filter_kernel) - filt_cov_yy = imfilter(cov_yy, box_filter_kernel) - - filt_cov_xx, filt_cov_xy, filt_cov_yy -end - -function gammacovs(cov_xx, cov_xy, cov_yy, gamma::Float64 = 1.4) - kernel = KernelFactors.gaussian((gamma, gamma)) - - filt_cov_xx = imfilter(cov_xx, kernel) - filt_cov_xy = imfilter(cov_xy, kernel) - filt_cov_yy = imfilter(cov_yy, kernel) - - filt_cov_xx, filt_cov_xy, filt_cov_yy -end diff --git a/src/misc.jl b/src/misc.jl deleted file mode 100644 index 31beb19a..00000000 --- a/src/misc.jl +++ /dev/null @@ -1,55 +0,0 @@ -""" - Percentile(x) - -Indicate that `x` should be interpreted as a [percentile](https://en.wikipedia.org/wiki/Percentile) rather than an absolute value. For example, - -- `canny(img, 1.4, (80, 20))` uses absolute thresholds on the edge magnitude image -- `canny(img, 1.4, (Percentile(80), Percentile(20)))` uses percentiles of the edge magnitude image as threshold -""" -struct Percentile{T} <: Real p::T end - - -""" -HomogeneousPoint(x::NTuple{N, T}) - -In projective geometry [homogeneous coordinates](https://en.wikipedia.org/wiki/Homogeneous_coordinates) are the -natural coordinates for describing points and lines. - -For instance, the homogeneous coordinates for a planar point are a triplet of real numbers ``(u, v ,w)``, with ``w \\neq 0``. -This triple can be associated with a point ``P = (x,y)`` in Cartesian coordinates, where ``x = \\frac{u}{w}`` and ``y = \\frac{v}{w}`` -[(more details)](http://www.geom.uiuc.edu/docs/reference/CRC-formulas/node6.html#SECTION01140000000000000000). - -In particular, the `HomogeneousPoint((10.0,5.0,1.0))` is the standardised projective representation of the Cartesian -point `(10.0,5.0)`. -""" -struct HomogeneousPoint{T <: AbstractFloat,N} - coords::NTuple{N, T} -end - -# By overwriting Base.to_indices we can define how to index into an N-dimensional array -# given an (N+1)-dimensional [`HomogeneousPoint`](@ref) type. -# We do this by converting the homogeneous coordinates to Cartesian coordinates -# and rounding to nearest integer. -# -# For homogeneous coordinates of a planar point we return -# a tuple of permuted Cartesian coordinates, (y,x), since matrices -# are indexed according to row and then column. -# For homogeneous coordinates of other dimensions we do not permute -# the corresponding Cartesian coordinates. -Base.to_indices(A::AbstractArray, p::Tuple{<: HomogeneousPoint}) = homogeneous_point_to_indices(p[1]) - -function homogeneous_point_to_indices(p::HomogeneousPoint{T,3}) where T - if p.coords[end] == 1 - return round(Int, p.coords[2]), round(Int, p.coords[1]) - else - return round(Int, p.coords[2] / p.coords[end]), round(Int, p.coords[1] / p.coords[end]) - end -end - -function homogeneous_point_to_indices(p::HomogeneousPoint) - if p.coords[end] == 1 - return round.(Int, p.coords) - else - return round.(Int, p.coords ./ p.coords[end]) - end -end diff --git a/test/corner.jl b/test/corner.jl index 45219a1e..5f1b42ac 100644 --- a/test/corner.jl +++ b/test/corner.jl @@ -148,7 +148,7 @@ using Test, Images @test all(ids[i].coords .≈ corner_pts[i].coords) end - # User specifies percentile. + # User specifies Percentile. corner_pts = imcorner_subpixel(img, Percentile(98), method = harris) @test length(corner_pts) == length(ids) for i = 1:length(ids) @@ -216,7 +216,7 @@ using Test, Images @test all(ids[i].coords .≈ corner_pts_offset[i].coords) end - # User specifies percentile. + # User specifies Percentile. corner_pts_offset = imcorner_subpixel(imgo, Percentile(98), method = harris) @test length(corner_pts_offset) == length(ids) for i = 1:length(ids) @@ -382,6 +382,4 @@ using Test, Images @test all(corners[6:12, 6:12] .== false) end -end - -nothing +end \ No newline at end of file