Merge pull request #685 from zygmuntszpak/imcorner_subpixel

Add sub-pixel corner detection capability #682

REQUIRE

@@ -13,6 +13,7 @@ ImageMetadata
 IndirectArrays
 MappedArrays
 SIUnits
+StaticArrays
 Graphics 0.1
 FileIO
 Compat 0.19

src/Images.jl

@@ -13,11 +13,13 @@ import Base: abs, atan2, clamp, convert, copy, copy!, ctranspose, delete!,
              strides, sum, write, zero
 
 export float32, float64
+export HomogeneousPoint
 
 using Base: depwarn
 using Base.Order: Ordering, ForwardOrdering, ReverseOrdering
 
 using Compat
+using StaticArrays
 
 # "deprecated imports" are below
 
@@ -64,6 +66,52 @@ Indicate that `x` should be interpreted as a [percentile](https://en.wikipedia.o
 """
 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
+
 include("labeledarrays.jl")
 include("algorithms.jl")
 include("exposure.jl")
@@ -154,6 +202,8 @@ export # types
     forwarddiffx,
     forwarddiffy,
     imcorner,
+    imcorner_subpixel,
+    corner2subpixel,
     harris,
     shi_tomasi,
     kitchen_rosenfeld,
@@ -244,7 +294,7 @@ Algorithms:
     - Exposure : `imhist`, `histeq`, `adjust_gamma`, `histmatch`, `imadjustintensity`, `imstretch`, `imcomplement`, `clahe`, `cliphist`
     - Gradients: `backdiffx`, `backdiffy`, `forwarddiffx`, `forwarddiffy`, `imgradients`
     - Edge detection: `imedge`, `imgradients`, `thin_edges`, `magnitude`, `phase`, `magnitudephase`, `orientation`, `canny`
-    - Corner detection: `imcorner`, `harris`, `shi_tomasi`, `kitchen_rosenfeld`, `meancovs`, `gammacovs`, `fastcorners`
+    - Corner detection: `imcorner`,`imcorner_subpixel`, `harris`, `shi_tomasi`, `kitchen_rosenfeld`, `meancovs`, `gammacovs`, `fastcorners`
     - Blob detection: `blob_LoG`, `findlocalmaxima`, `findlocalminima`
     - Morphological operations: `dilate`, `erode`, `closing`, `opening`, `tophat`, `bothat`, `morphogradient`, `morpholaplace`, `feature_transform`, `distance_transform`, `convexhull`
     - Connected components: `label_components`, `component_boxes`, `component_lengths`, `component_indices`, `component_subscripts`, `component_centroids`

src/corner.jl

@@ -35,6 +35,142 @@ function imcorner(img::AbstractArray, thresholdp::Percentile; method::Function =
     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 = indices(corner_indicator)
+    row, col, _ = findnz(corner_indicator)
+    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])