Nonmaximal suppression updates

* Added tests
* Fixed angle discretization for row-major images
* Fixed subpixel offset calc
* Changed angle discretization from (0,pi) to (0,2pi).  By testing
  in the direction of the actual gradient first, we reduce the number
  of interpolations in the other direction by ~1/4
* Added thin_edges, thin_edges_subpix.  Right now, these
  just point to thin_edges_nonmaxsup and thin_edges_nonmaxsup_subpix,
  but could be expanded to use other edge thinning methods in the
  future.

src/Images.jl

@@ -234,6 +234,8 @@ export # types
     sobel,
     ssd,
     ssdn,
+    thin_edges,
+    thin_edges_subpix,
     thin_edges_nonmaxsup,
     thin_edges_nonmaxsup_subpix,
 

src/algorithms.jl

@@ -1571,6 +1571,12 @@ function imedge(img::AbstractArray, method::String="ando3", border::String="repl
     return (grad_x, grad_y, mag, orient)
 end
 
+# Thin edges
+thin_edges{T}(img::AbstractArray{T,2}, gradientangles::AbstractArray=imgradients(img), border::String="replicate") =
+    thin_edges_nonmaxsup(img, gradientangles, border)
+thin_edges_subpix{T}(img::AbstractArray{T,2}, gradientangles::AbstractArray=imgradients(img), border::String="replicate") =
+    thin_edges_nonmaxsup_subpix(img, gradientangles, border)
+
 function imROF{T}(img::Array{T,2}, lambda::Number, iterations::Integer)
     # Total Variation regularized image denoising using the primal dual algorithm
     # Also called Rudin Osher Fatemi (ROF) model

src/nonmaxsup.jl

@@ -1,7 +1,7 @@
 # NONMAXSUP - Non-maxima suppression
 #
 # Usage:
-#          (im,location) = nonmaxsup(img, gradientangle, radius);
+#          (im,location) = nonmaxsup(img, gradientangles, radius);
 #
 # Function for performing non-maxima suppression on an image using an
 # orientation image.  It is assumed that the orientation image gives 
@@ -10,7 +10,7 @@
 # Input:
 #   img - image to be non-maxima suppressed.
 # 
-#   gradientangle  - image containing gradient angles around each pixel in radians
+#   gradientangles  - image containing gradient angles around each pixel in radians
 #             (-pi,pi)
 # 
 #   radius  - Distance in pixel units to be looked at on each side of each
@@ -69,39 +69,43 @@ end
 
 CoordOffset(x::Float64) = ((frac,i) = modf(x); CoordOffset(sign(frac), int(i), abs(frac)))
 (-)(off::CoordOffset) = CoordOffset(-off.s,-off.i, off.f)
+(*)(x::Number, off::CoordOffset) = x*(off.i + off.s*off.f)
+(*)(off::CoordOffset, x::Number) = x*(off.i + off.s*off.f)
+(+)(x::Number, off::CoordOffset) = x + off.i + off.s*off.f
+(+)(off::CoordOffset, x::Number) = x + off.i + off.s*off.f
 
 # Precalculate x and y offsets relative to centre pixel for each orientation angle 
 function _calc_discrete_offsets(θ, radius, transposed)
 
-    θ_count = iround(pi/θ)
-    θ = pi/θ_count
+    θ_count = iround(2π/θ)
+    θ = 2π/θ_count
     angles = (0:θ_count)*θ
 
     # x and y offset of points at specified radius and angles
     # from each reference position.
 
     if transposed
         # θ′ = -π/2 - θ
-        xoffs = [CoordOffset( x) for x in -radius*sin(angles)]
-        yoffs = [CoordOffset(-y) for y in  radius*cos(angles)]
+        xoffs = [CoordOffset(-x) for x in  radius*sin(angles)]
+        yoffs = [CoordOffset( y) for y in  radius*cos(angles)]
     else
-        xoffs = [CoordOffset( x) for x in radius*cos(angles)]
-        yoffs = [CoordOffset(-y) for y in radius*sin(angles)]
+        xoffs = [CoordOffset( x) for x in  radius*cos(angles)]
+        yoffs = [CoordOffset(-y) for y in  radius*sin(angles)]
     end
 
     return xoffs, yoffs
 end
 
 # Discretize an array of orientation angles
-function _discretize_angles(gradientangle::AbstractArray, θ=pi/180)
+function _discretize_angles(gradientangles::AbstractArray, θ=pi/180)
     iθ = 1/θ
-    d_gradientangle = similar(gradientangle, Int)
-    for i = 1:length(d_gradientangle)
-        g = gradientangle[i]
-        angle = g < 0 ? (g+π)*iθ : g*iθ
-        d_gradientangle[i] = iround(angle)+1
+    d_gradientangles = similar(gradientangles, Int)
+    for i = 1:length(d_gradientangles)
+        g = gradientangles[i]
+        angle = g < 0 ? (g+2π)*iθ : g*iθ
+        d_gradientangles[i] = iround(angle)+1
     end
-    d_gradientangle
+    d_gradientangles
 end
 
 # Interpolate the value of an offset from a particular pixel
@@ -131,20 +135,22 @@ function _interp_offset(img::AbstractArray, x::Integer, y::Integer, xoff::CoordO
     return (upperavg + yoff.f * (loweravg - upperavg), min_adjacent)
 end
 
+cv2_count = 0
+
 # Core thinning algorithm using nonmaximal suppression
-function thin_edges_nonmaxsup_core!{T}(location, img::AbstractArray{T,2}, gradientangle, radius, border, theta)
+function thin_edges_nonmaxsup_core!{T}(location, img::AbstractArray{T,2}, gradientangles, radius, border, theta)
     calc_subpixel = !isempty(location)
 
     # Error checking
-    size(img) == size(gradientangle) || error("image and gradient orientation image are of different sizes")
+    size(img) == size(gradientangles) || error("image and gradient orientation image are of different sizes")
     radius < 1.0 && error("radius must be >= 1")
 
     # Precalculate x and y offsets relative to centre pixel for each orientation angle 
     transposed = spatialorder(img)[1] == "x"
     xoffs, yoffs = _calc_discrete_offsets(theta, radius, transposed)
 
     # Discretize orientation image
-    d_gradientangle = _discretize_angles(gradientangle, theta)
+    d_gradientangles = _discretize_angles(gradientangles, theta)
 
     # Allocate output
     (height,width) = size(img)
@@ -159,9 +165,9 @@ function thin_edges_nonmaxsup_core!{T}(location, img::AbstractArray{T,2}, gradie
     # of the centre pixel to be used for the non-maximal suppression.
 
     for x = 1:width, y = 1:height
-        (c = img[y,x]) == 0.0 && continue  # For thresholded images
+        (c = img[y,x]) == 0 && continue  # For thresholded images
 
-        or = d_gradientangle[y,x]   # Disretized orient
+        or = d_gradientangles[y,x]   # Disretized orient
         v1, n1 = _interp_offset(img, x, y, xoffs[or], yoffs[or], Ix, Iy, pad)
 
         if (c > v1) & (c >= n1) # We need to check the value on the other side...
@@ -180,8 +186,8 @@ function thin_edges_nonmaxsup_core!{T}(location, img::AbstractArray{T,2}, gradie
 
                     # location where maxima of fitted parabola occurs
                     r = -b/2a
-                    location[y,x] = transposed ? Point(y - r*yoffs[or].f, x + r*xoffs[or].f) :
-                                                 Point(x + r*xoffs[or].f, y - r*yoffs[or].f)
+                    location[y,x] = transposed ? Point(y + r*yoffs[or], x + r*xoffs[or]) :
+                                                 Point(x + r*xoffs[or], y + r*yoffs[or])
 
                     if T<:FloatingPoint
                         # Store the interpolated value
@@ -200,17 +206,17 @@ function thin_edges_nonmaxsup_core!{T}(location, img::AbstractArray{T,2}, gradie
 end
 
 # Main function call when subpixel location of edges is not needed
-thin_edges_nonmaxsup(img::AbstractArray, gradientangle::AbstractArray;
-                     radius::Float64=1.35, border="replicate", theta=pi/180) =
-    copy(img, thin_edges_nonmaxsup_core!(Array(Point,(0,0)), img, gradientangle, radius, border, theta))
+thin_edges_nonmaxsup(img::AbstractArray, gradientangles::AbstractArray, border::String="replicate";
+                     radius::Float64=1.35, theta=pi/180) =
+    copy(img, thin_edges_nonmaxsup_core!(Array(Point,(0,0)), img, gradientangles, radius, border, theta))
 
 
 # Main function call when subpixel location of edges is desired
-function thin_edges_nonmaxsup_subpix(img::AbstractArray, gradientangle::AbstractArray;
-                                     radius::Float64=1.35, border="replicate", theta=pi/180)
+function thin_edges_nonmaxsup_subpix(img::AbstractArray, gradientangles::AbstractArray,
+                                     border::String="replicate"; radius::Float64=1.35, theta=pi/180)
     (height,width) = size(img)
     location = zeros(Point, height, width)
-    out = thin_edges_nonmaxsup_core!(location, img, gradientangle, radius, border, theta)
+    out = thin_edges_nonmaxsup_core!(location, img, gradientangles, radius, border, theta)
 
     copy(img, out), copy(img, location)
 end

test/algorithms.jl

@@ -453,3 +453,156 @@ for method in ["sobel", "prewitt", "ando3", "ando4", "ando5", "ando4_sep", "ando
                 ((gphase .== 0.0) & (orient .== 0.0)))  # this part is where both are 
                                                         # zero because there is no gradient
 end
+
+# Nonmaximal suppression
+
+function thin_edges(img)
+    # Get orientation
+    gx,gy = imgradients(img)
+    orient = phase(gx,gy)
+
+    # Do NMS thinning
+    thin_edges_nonmaxsup_subpix(img, orient, radius=1.35)
+end
+
+
+function nms_test_horiz_vert(img, which)
+    ## which = :horizontal or :vertical
+
+    # Do NMS thinning
+    t,s = thin_edges(img)
+
+    # Calc peak location by hand
+
+    # Interpolate values 1.35 pixels left and right
+    # Orientation is zero radians -> to the right
+    v1 = 6 - 0.35   # slope on right is 1
+    v2 = 5 - 0.35*2 # slope on left is 2
+    c = 7.0         # peak value
+
+    # solve v = a*r^2 + b*r + c
+    a = (v1 + v2)/2 - c
+    b = a + c - v2
+    r = -b/2a
+
+    @assert abs(r - 1/6) < EPS
+
+    # Location and value at peak
+    peakloc = r*1.35 + 3
+    peakval = a*r^2 + b*r + c
+
+    transposed = spatialorder(img)[1] == "x"
+    horizontal = which == :horizontal
+
+    test_axis1 = transposed $ !horizontal
+
+    @assert test_axis1 ? all(t[:,[1,2,4,5]] .== 0) : all(t[[1,2,4,5],:] .== 0)
+    @assert test_axis1 ? all(t[:,3]   .== peakval) : all(t[3,:]   .== peakval)
+    @assert test_axis1 ? all(s[:,[1,2,4,5]] .== zero(Base.Graphics.Point)) :
+                         all(s[[1,2,4,5],:] .== zero(Base.Graphics.Point))
+
+    if transposed
+        if which == :horizontal
+            @assert     [pt.x for pt in s[:,3]]  == [1:5]
+            @assert all([pt.y for pt in s[:,3]] .== peakloc)
+        else
+            @assert all([pt.x for pt in s[3,:]] .== peakloc)
+            @assert     [pt.y for pt in s[3,:]]  == [1:5]
+        end
+    else
+        if which == :horizontal
+            @assert     [pt.x for pt in s[3,:]]  == [1:5]
+            @assert all([pt.y for pt in s[3,:]] .== peakloc)
+        else
+            @assert all([pt.x for pt in s[:,3]] .== peakloc)
+            @assert     [pt.y for pt in s[:,3]]  == [1:5]
+        end
+    end
+end
+
+# Test image: vertical edge
+m = [3.0  5.0  7.0  6.0  5.0
+     3.0  5.0  7.0  6.0  5.0
+     3.0  5.0  7.0  6.0  5.0
+     3.0  5.0  7.0  6.0  5.0
+     3.0  5.0  7.0  6.0  5.0]
+
+m_xy = grayim(m')
+m_yx = grayim(m)
+m_yx["spatialorder"] = ["y","x"]
+
+nms_test_horiz_vert(m, :vertical)
+nms_test_horiz_vert(m_xy, :vertical)
+nms_test_horiz_vert(m_yx, :vertical)
+
+# Test image: horizontal edge
+m = m'
+m_xy = grayim(m')
+m_yx = grayim(m)
+m_yx["spatialorder"] = ["y","x"]
+
+nms_test_horiz_vert(m, :horizontal)
+nms_test_horiz_vert(m_xy, :horizontal)
+nms_test_horiz_vert(m_yx, :horizontal)
+
+
+function nms_test_diagonal(img)
+    # Do NMS thinning
+    t,s = thin_edges(img)
+
+    # Calc peak location by hand
+
+    # Interpolate values 1.35 pixels up and left, down and right
+    # using bilinear interpolation
+    # Orientation is π/4 radians -> 45 degrees up
+    fr = 1.35*cos(π/4)
+    lower = (7 + fr*(6-7))
+    upper = (6 + fr*(5-6))
+    v1 = lower + fr*(upper-lower)
+
+    lower = (7 + fr*(5-7))
+    upper = (5 + fr*(3-5))
+    v2 = lower + fr*(upper-lower)
+
+    c = 7.0         # peak value
+
+    # solve v = a*r^2 + b*r + c
+    a = (v1 + v2)/2 - c
+    b = a + c - v2
+    r = -b/2a
+
+    @assert (r - 1/6) < EPS
+
+    transposed = spatialorder(img)[1] == "x"
+
+    # Location and value at peak
+
+    x_peak_offset, y_peak_offset = r*fr, -r*fr
+    peakval = a*r^2 + b*r + c
+
+    @assert all(diag(data(t))[2:4] .== peakval)  # Edge pixels aren't interpolated here
+    @assert all(t - diagm(diag(data(t))) .== 0)
+
+    diag_s = copy(s, diagm(diag(data(s))))
+    @assert s == diag_s
+
+    @assert all([pt.x for pt in diag(data(s))[2:4]] - ([2:4] + x_peak_offset) .< EPS)
+    @assert all([pt.y for pt in diag(data(s))[2:4]] - ([2:4] + y_peak_offset) .< EPS)
+
+end
+
+
+# Test image: diagonal edge
+m = [7.0  6.0  5.0  0.0  0.0
+     5.0  7.0  6.0  5.0  0.0
+     3.0  5.0  7.0  6.0  5.0
+     0.0  3.0  5.0  7.0  6.0
+     0.0  0.0  3.0  5.0  7.0]
+
+m_xy = grayim(m')
+m_yx = grayim(m)
+m_yx["spatialorder"] = ["y","x"]
+
+nms_test_diagonal(m)
+nms_test_diagonal(m_xy)
+nms_test_diagonal(m_yx)