Skip to content

Commit

Permalink
Adds ability to determine gradient (edge) orientation
Browse files Browse the repository at this point in the history
  • Loading branch information
zygmuntszpak committed Dec 1, 2020
1 parent 073e75d commit 770f34f
Show file tree
Hide file tree
Showing 7 changed files with 402 additions and 8 deletions.
48 changes: 48 additions & 0 deletions docs/src/index.md
Original file line number Diff line number Diff line change
Expand Up @@ -60,6 +60,7 @@ mkpath("images")
```@example Canny
using TestImages, ImageEdgeDetection, MosaicViews
using FileIO # hide
using ImageCore # hide
img = testimage("mandril_gray")
# Detect edges at different scales by adjusting the `spatial_scale` parameter.
img_edges₁ = detect_edges(img, Canny(spatial_scale = 1.4))
Expand Down Expand Up @@ -118,4 +119,51 @@ save("images/demo3.jpg", demo₃); # hide
<p>
```

One can also determine the gradient orientation in an adjustable manner by
defining an [`OrientationConvention`](@ref
ImageEdgeDetection.OrientationConvention). An `OrientationConvention` allows you
to specify the compass direction against which you intend to measure the angle,
and whether you are measuring in a clockwise or counter-clockwise manner.

In the example below, we map the angles `[0...360]` to the unit interval to
visualise the orientation of the circle edge test image using different
orientation conventions. Note that the angle `360` is used as a sentinel value
to demarcate pixels for which the gradient orientation is undefined. The
gradient orientation is undefined when the gradient magnitude is effectively
zero. This corresponds to regions of constant intensity in the image. In the
images that depict the gradient orientation, the undefined orientations are
represent as pure white pixels.

```@example GradientOrientation
# Create a test image (black circle against a white background).
a = 250
b = 250
r = 150
img = Gray.(ones(500, 500))
for i in CartesianIndices(img)
y, x = i.I
img[i] = (x-a)^2 + (y - b)^2 - r^2 < 0 ? 0.0 : 1.0
end
# Determine the image gradients
g₁, g₂ = imgradients(img, KernelFactors.sobel)
orientation_convention₁ = OrientationConvention(in_radians = false, compass_direction = 'S')
orientation_convention₂ = OrientationConvention(in_radians = false, compass_direction = 'N')
orientation_convention₂ = OrientationConvention(in_radians = false, compass_direction = 'E', is_clockwise = true)
angles₁ = detect_gradient_orientation(g₁, g₂, orientation_convention₁) / 360
angles₂ = detect_gradient_orientation(g₁, g₂, orientation_convention₂) / 360
angles₃ = detect_gradient_orientation(g₁, g₂, orientation_convention₃) / 360
demo₃ = mosaicview(img, Gray.(angles₁), Gray.(angles₂), Gray.(angles₃); nrow = 2)
save("images/demo4.jpg", demo₄); # hide
```
```@raw html
<img src="images/demo4.jpg" width="512px" alt="gradient orientation demo image" />
<p>
```


For more advanced usage, please check [function reference](@ref function_reference) page.
15 changes: 12 additions & 3 deletions docs/src/reference.md
Original file line number Diff line number Diff line change
Expand Up @@ -10,12 +10,17 @@ Depth = 3
```@docs
detect_edges
detect_edges!
detect_subpixel_edges
detect_subpixel_edges!
detect_gradient_orientation
detect_gradient_orientation!
thin_edges
thin_edges!
thin_subpixel_edges
thin_subpixel_edges!
```

## Edge Detection Algorithms

```@docs
ImageEdgeDetection.EdgeDetectionAPI.AbstractEdgeDetectionAlgorithm
```
Expand All @@ -26,7 +31,6 @@ ImageEdgeDetection.Canny
```

## Edge Thinning Algorithms

```@docs
ImageEdgeDetection.EdgeDetectionAPI.AbstractEdgeThinningAlgorithm
```
Expand All @@ -38,7 +42,12 @@ ImageEdgeDetection.NonmaximaSuppression

### Non-maxima Suppression (Subpixel)
```@docs
ImageEdgeDetection.NonmaximaSubpixelSuppression
ImageEdgeDetection.SubpixelNonmaximaSuppression
```

### OrientationConvention
```@docs
ImageEdgeDetection.OrientationConvention
```

## Supplementary Types
Expand Down
12 changes: 8 additions & 4 deletions src/ImageEdgeDetection.jl
Original file line number Diff line number Diff line change
Expand Up @@ -18,9 +18,9 @@ import .EdgeDetectionAPI: AbstractEdgeDetectionAlgorithm,
detect_edges, detect_edges!,
detect_subpixel_edges, detect_subpixel_edges!

import .EdgeDetectionAPI: AbstractEdgeThinningAlgorithm,
thin_edges, thin_edges!,
thin_subpixel_edges, thin_subpixel_edges!
import .EdgeDetectionAPI: AbstractEdgeThinningAlgorithm,
thin_edges, thin_edges!,
thin_subpixel_edges, thin_subpixel_edges!

# TODO Relax this to all image color types
const GenericGrayImage = AbstractArray{<:Union{Number, AbstractGray}}
Expand All @@ -39,6 +39,7 @@ end
include("algorithms/nonmaxima_suppression.jl")
include("algorithms/subpixel_nonmaxima_suppression.jl")
include("algorithms/canny.jl")
include("algorithms/gradient_orientation.jl")

# Set the Canny algorithm as the default edge detection algorithm.
detect_edges(img::AbstractArray,
Expand All @@ -53,6 +54,7 @@ export
Canny,
NonmaximaSuppression,
SubpixelNonmaximaSuppression,
OrientationConvention,
Percentile,
thin_edges,
thin_edges!,
Expand All @@ -61,6 +63,8 @@ export
detect_edges,
detect_edges!,
detect_subpixel_edges,
detect_subpixel_edges!
detect_subpixel_edges!,
detect_gradient_orientation,
detect_gradient_orientation!

end # module
204 changes: 204 additions & 0 deletions src/algorithms/gradient_orientation.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,204 @@
"""
```
OrientationConvention(; compass_direction::AbstractChar = 'S', is_clockwise::Bool = false, in_radians = true, tol = sqrt(eps(Float64)))
```
Specifies the coordinate system context for `detect_gradient_orientation` which
determines the meaning of the angles (the gradient orientations).
# Details
You can specify how you want the gradient orientation to be reported. By
default, the orientation is measured counter-clockwise from the south direction.
This is because in a Raster coordinate system, the first spatial dimension
increases as one goes down the image (i.e. it points south), and the second
spatial dimension increases as one moves to the right of the image (i.e. it
points east).
If you wish to interpret the orientation in a canonical Cartesian coordinate
convention you would specify east as the reference compass direction
(`compass_direction = 'E'`) and a counter-clockwise direction (`clockwise =
false`).
If `in_radians = true` the valid angles are reported in the range of `[0...2π)`,
otherwise they are reported in the range `[0...360)`. The values `2π` and `360`
are used as sentinels to designate undefined angles (because the gradient
magnitude was too close to zero). By default, an angle is undefined if `(abs(g₁)
< tol && abs(g₂) < tol)` where `g₁` and `g₂` denote the gradient in the first
and second spatial dimensions, and `tol = sqrt(eps(Float64))`.
# Example
```julia
using TestImages, FileIO, ImageView, ImageEdgeDetection, ImageFiltering
img = testimage("mandril_gray")
# Gradient in the first and second spatial dimension
g₁, g₂ = imgradients(img, KernelFactors.scharr)
# Interpret the angles with respect to a canonical Cartesian coordinate system
# where the angles are measured counter-clockwise from the positive x-axis.
orientation_convention = OrientationConvention(in_radians = true,
is_clockwise = false,
compass_direction = 'E')
angles = detect_gradient_orientation(g₁, g₂, orientation_convention)
```
"""
@with_kw struct OrientationConvention{ T <: AbstractChar}
# 'N', 'S', 'E ', 'W' for (north, south, east, west)
compass_direction::T = 'S'
is_clockwise::Bool = false
in_radians::Bool = true
tol = sqrt(eps(eltype(Float64)))
end


"""
detect_gradient_orientation(g₁::AbstractArray, g₂::AbstractArray, orientation_convention::OrientationConvention, args...; kwargs...)
Given the gradient in the first (`g₁`) and second (`g₂`) spatial dimensions,
returns the gradient orientation, where the orientation is interpreted according
to a supplied [`OrientationConvention`](@ref
ImageEdgeDetection.OrientationConvention).
# Details
You can specify how you want the gradient orientation to be reported by
supplying an [`OrientationConvention`](@ref
ImageEdgeDetection.OrientationConvention). If left unspecified the orientation
is measured counter-clockwise from the south direction. This is because in a
Raster coordinate system, the first spatial dimension increases as one goes down
the image (i.e. it points south), and the second spatial dimension increases as
one moves to the right of the image (i.e. it points east).
If you wish to interpret the orientation in a canonical Cartesian coordinate
convention you would specify east as the reference compass direction
(`compass_direction = 'E'`) and a counter-clockwise direction (`clockwise =
false`).
# Output
Returns a two-dimensional array of angles. If `in_radians = true` the valid
angles are reported in the range of `[0...2π)`, otherwise they are reported in
the range `[0...360)`. The values `2π` and `360` are used as sentinels to
designate undefined angles (because the gradient magnitude was too close to
zero). By default, an angle is undefined if `(abs(g₁) < tol && abs(g₂) < tol)`
where `g₁` and `g₂` denote the gradient in the first and second spatial
dimensions, and `tol = sqrt(eps(Float64))` (as defined in
[`OrientationConvention`](@ref ImageEdgeDetection.OrientationConvention)).
See also: [`detect_gradient_orientation!`](@ref)
"""
function detect_gradient_orientation(g₁::AbstractArray, g₂::AbstractArray, orientation_convention::OrientationConvention)
out = zeros(axes(g₁))
detect_gradient_orientation!(out, g₁, g₂, orientation_convention)
return out
end

function detect_gradient_orientation(g₁::AbstractArray, g₂::AbstractArray)
out = zeros(axes(g₁))
detect_gradient_orientation!(out, g₁, g₂, OrientationConvention())
return out
end

"""
detect_gradient_orientation(out::AbstractArray, g₁::AbstractArray, g₂::AbstractArray, orientation_convention::OrientationConvention, args...; kwargs...)
Given the gradient in the first (`g₁`) and second (`g₂`) spatial dimensions,
returns the gradient orientation in `out`, where the orientation is interpreted
according to a supplied [`OrientationConvention`](@ref
ImageEdgeDetection.OrientationConvention).
# Details
You can specify how you want the gradient orientation to be reported by
supplying an [`OrientationConvention`](@ref
ImageEdgeDetection.OrientationConvention). If left unspecified the orientation
is measured counter-clockwise in radians from the south direction. This is
because in a Raster coordinate system, the first spatial dimension increases as
one goes down the image (i.e. it points south), and the second spatial dimension
increases as one moves to the right of the image (i.e. it points east).
If you wish to interpret the orientation in a canonical Cartesian coordinate
convention you would specify east as the reference compass direction
(`compass_direction = 'E'`) and a counter-clockwise direction (`clockwise =
false`).
# Output
Returns a two-dimensional array of angles. If `in_radians = true` genuine angles are
reported in the range of `[0...2π)`, otherwise they are reported in the range
`[0...360)`. The values `2π` and `360` are used as sentinels to designate
undefined angles (because the gradient magnitude was too close to zero).
By default, an angle is undefined if `(abs(g₁) < tol && abs(g₂) < tol)` where
`g₁` and `g₂` denote the gradient in the first and second spatial dimensions,
and `tol = sqrt(eps(eltype(out)))` (as defined in [`OrientationConvention`](@ref
ImageEdgeDetection.OrientationConvention)).
See also: [`detect_gradient_orientation`](@ref)
"""
function detect_gradient_orientation!(out::AbstractArray, g₁::AbstractArray, g₂::AbstractArray, orientation_convention::OrientationConvention)
@unpack compass_direction, is_clockwise, in_radians = orientation_convention
tol = sqrt(eps(eltype(out)))

# Determine from which compass direction we intend to measure the angle.
if compass_direction == 'S' || compass_direction == 's'
offset = π/2
elseif compass_direction == 'E' || compass_direction == 'e'
offset = 0.0
elseif compass_direction == 'N' || compass_direction == 'n'
offset = -π/2
elseif compass_direction == 'W' || compass_direction == 'w'
offset = π
else
@warn("Unrecognised compass_direction... using a default direction of south (S).")
offset = π/2
end


if !is_clockwise && !in_radians
@inbounds for i in CartesianIndices(g₁)
out[i] = rad2deg(zero_to_2PI(valid_angle(g₁[i], g₂[i], offset, tol)))
end
elseif !is_clockwise && in_radians
@inbounds for i in CartesianIndices(g₁)
out[i] = zero_to_2PI(valid_angle(g₁[i], g₂[i], offset, tol))
end
elseif is_clockwise && !in_radians
@inbounds for i in CartesianIndices(g₁)
out[i] = rad2deg(zero_to_2PI(clockwise(valid_angle(g₁[i], g₂[i], offset, tol))))
end
elseif is_clockwise && in_radians
@inbounds for i in CartesianIndices(g₁)
out[i] = zero_to_2PI(clockwise(valid_angle(g₁[i], g₂[i], offset, tol)))
end
end
return nothing
end

function detect_gradient_orientation!(out::AbstractArray, g₁::AbstractArray, g₂::AbstractArray)
detect_gradient_orientation!(out, g₁, g₂, OrientationConvention())
return out
end

# When the angle is undefined because g₁ and g₂ are almost zero, we return
# a "dummy" value of 2π.
function valid_angle(g₁, g₂, offset, tol)
# The expression for the angle when changing coordinate system from Raster
# to Cartesian coordinates is atan(-g₁, g₂). The offset is used to indicate
# against which compass direction we will measure angles.
is_angle_undefined = (abs(g₁) < tol && abs(g₂) < tol)
return is_angle_undefined ? 2π : atan(-g₁, g₂) + offset
end

function zero_to_2PI(θ)
return>= 0 ? θ : (2*π + θ))
end

function clockwise(x)
return mod(-x, 2π)
end
2 changes: 1 addition & 1 deletion src/algorithms/subpixel_nonmaxima_suppression.jl
Original file line number Diff line number Diff line change
Expand Up @@ -54,7 +54,7 @@ imshow(nms)
1. J. Canny, "A Computational Approach to Edge Detection," in IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. PAMI-8, no. 6, pp. 679-698, Nov. 1986, doi: 10.1109/TPAMI.1986.4767851.
2. F. Devernay, "A non-maxima suppression method for edge detection with sub-pixel accuracy", Tech. Report RR-2724, INRIA, 1995.
"""
@with_kw struct SubpixelNonmaximaSuppression{ T <: Union{Real,AbstractGray, Percentile}} <: AbstractEdgeThinningAlgorithm
@with_kw struct SubpixelNonmaximaSuppression{T <: Union{Real,AbstractGray, Percentile}} <: AbstractEdgeThinningAlgorithm
threshold::T = Percentile(20)
end

Expand Down
Loading

0 comments on commit 770f34f

Please sign in to comment.