Merge pull request #11 from JuliaImages/brisk

BRISK Descriptors

Author: mronian

src/ImageFeatures.jl

@@ -11,8 +11,9 @@ include("glcm.jl")
 include("brief.jl")
 include("orb.jl")
 include("freak.jl")
+include("brisk.jl")
 
-export Keypoint, Keypoints, BRIEF, DescriptorParams, ORB, FREAK
+export Keypoint, Keypoints, Feature, Features, Params, BRIEF, ORB, FREAK, BRISK
 
 export 
 	#Core

src/brief.jl

@@ -12,7 +12,7 @@ brief_params = BRIEF([size = 128], [window = 9], [sigma = 2 ^ 0.5], [sampling_ty
 | **seed** | Int | Random seed used for generating the sampling pairs. For matching two descriptors, the seed used to build both should be same. | 
 
 """
-type BRIEF{F} <: DescriptorParams
+type BRIEF{F} <: Params
     size::Int
     window::Int
     sigma::Float64

src/brisk.jl

@@ -0,0 +1,106 @@
+"""
+```
+brisk_params = BRISK([pattern_scale = 1.0])
+```
+
+| Argument | Type | Description |
+|----------|------|-------------|
+| `pattern_scale` | `Float64` | Scaling factor for the sampling window | 
+"""
+type BRISK <: Params
+    threshold::Float64
+    octaves::Int
+    pattern_scale::Float64
+    pattern_table::Vector{Vector{SamplePair}}
+    smoothing_table::Vector{Vector{Float16}}
+    orientation_weights::Vector{OrientationWeights}
+    short_pairs::Vector{OrientationPair}
+    long_pairs::Vector{OrientationPair}
+end
+
+function BRISK(; threshold::Float64 = 0.25, octaves::Int = 4, pattern_scale = 1.0)
+    pattern_table, smoothing_table = _brisk_tables(pattern_scale)
+    orientation_weights = OrientationWeights[]
+    short_pairs = OrientationPair[]
+    long_pairs = OrientationPair[]
+    dminsq = (brisk_dmin * pattern_scale) ^ 2
+    dmaxsq = (brisk_dmax * pattern_scale) ^ 2
+    for i in 2:brisk_points
+        for j in 1:i - 1
+            dy = pattern_table[1][j][1] - pattern_table[1][i][1]
+            dx = pattern_table[1][j][2] - pattern_table[1][i][2]
+            norm = dy ^ 2 + dx ^ 2
+            if norm > dminsq
+                push!(long_pairs, OrientationPair((j, i)))
+                push!(orientation_weights, OrientationWeights((dy / norm, dx / norm)))
+            elseif norm < dmaxsq
+                push!(short_pairs, OrientationPair((j, i)))
+            end
+        end
+    end
+    BRISK(threshold, octaves, pattern_scale, pattern_table, smoothing_table, orientation_weights, short_pairs, long_pairs)
+end
+
+function _brisk_orientation{T<:Gray}(int_img::AbstractArray{T, 2}, keypoint::Keypoint, pattern::Array{SamplePair}, 
+                                        orientation_weights::Array{OrientationWeights}, sigmas::Array{Float16}, long_pairs::Array{OrientationPair})
+    direction_sum_y = 0.0
+    direction_sum_x = 0.0
+    for (i, o) in enumerate(long_pairs)
+        offset_1 = pattern[o[1]]
+        offset_2 = pattern[o[2]]
+        intensity_diff = ImageFeatures._freak_mean_intensity(int_img, keypoint, offset_1, sigmas[o[1]]) - ImageFeatures._freak_mean_intensity(int_img, keypoint, offset_2, sigmas[o[2]])
+        direction_sum_y += orientation_weights[i][1] * intensity_diff
+        direction_sum_x += orientation_weights[i][2] * intensity_diff
+    end
+    angle = atan2(direction_sum_y, direction_sum_x)
+    scaled_angle = ceil(Int, (angle + pi) * brisk_orientation_steps / (2 * pi))
+    scaled_angle
+end
+
+function _brisk_tables(pattern_scale::Float64)
+    pattern_table = Vector{SamplePair}[]
+    smoothing_table = Vector{Float16}[]
+    for ori in 0:brisk_orientation_steps - 1
+        theta = ori * 2 * pi / brisk_orientation_steps 
+        pattern = SamplePair[]
+        sigmas = Float16[]
+        for (i, n) in enumerate(brisk_num_circular_pattern)
+            for circle_number in 0:n - 1
+                angle = (circle_number * 2 * pi / n) + theta
+
+                push!(pattern, SamplePair((brisk_radii[i] * sin(angle) * pattern_scale * 0.85, 
+                                            brisk_radii[i] * cos(angle) * pattern_scale * 0.85)))
+                push!(sigmas, brisk_sigma[i] * pattern_scale * 0.85)
+            end
+        end
+        push!(pattern_table, pattern)
+        push!(smoothing_table, sigmas)
+    end
+    pattern_table, smoothing_table
+end
+
+function create_descriptor{T<:Gray}(img::AbstractArray{T, 2}, features::Array{Feature}, params::BRISK)
+    int_img = integral_image(img)
+    descriptors = BitArray{1}[]
+    ret_features = Feature[]
+    window_size = ceil(Int, (brisk_radii[end] + brisk_sigma[end]) * params.pattern_scale * 0.85) + 1
+    lim = CartesianIndex(window_size, window_size)
+    for feature in features
+        keypoint = Keypoint(feature)
+        checkbounds(Bool, img, keypoint - lim) && checkbounds(Bool, img, keypoint + lim) || continue
+        orientation = _brisk_orientation(int_img, keypoint, params.pattern_table[1], params.orientation_weights, params.smoothing_table[1], params.long_pairs)
+        sampled_intensities = T[]
+        for (i, p) in enumerate(params.pattern_table[orientation])
+            push!(sampled_intensities, ImageFeatures._freak_mean_intensity(int_img, keypoint, p, params.smoothing_table[orientation][i]))
+        end
+        descriptor = falses(size(params.short_pairs, 1))
+        for (i, f) in enumerate(params.short_pairs)
+            point_1 = sampled_intensities[f[1]]
+            point_2 = sampled_intensities[f[2]]
+            descriptor[i] = point_1 < point_2
+        end
+        push!(descriptors, descriptor)
+        push!(ret_features, Feature(keypoint, orientation))
+    end
+    descriptors, ret_features   
+end

src/const.jl

@@ -96,3 +96,17 @@ const freak_radii = [freak_largest_radius, freak_largest_radius - 6 * freak_circ
 const freak_sigma = vcat(freak_radii[1:end - 1] / 2, freak_radii[end - 1] / 2)
 
 const freak_orientation_steps = 256;
+
+const brisk_radii = [0, 2.9, 4.9, 7.4, 10.8]
+
+const brisk_sigma = [0.65, 1.16, 1.41, 2.00, 2.19]
+
+const brisk_num_circular_pattern = [1, 10, 14, 15, 20]
+
+const brisk_orientation_steps = 1024
+
+const brisk_dmin = 8.2
+
+const brisk_dmax = 5.85
+
+const brisk_points = 60

src/core.jl

@@ -1,29 +1,70 @@
-abstract Detector
-abstract DescriptorParams
+abstract Params
 
 """
 ```
 keypoint = Keypoint(y, x)
+keypoint = Keypoint(feature)
 ```
 
-The `Keypoint` type.
+A `Keypoint` may be created by passing the coordinates of the point or from a feature.
 """
 typealias Keypoint CartesianIndex{2}
 
 """
 ```
 keypoints = Keypoints(boolean_img)
+keypoints = Keypoints(features)
 ```
 
-Creates a `Vector{Keypoint}` of the true values in the `boolean_img`.
+Creates a `Vector{Keypoint}` of the `true` values in a boolean image or from a list of features.
 """
 typealias Keypoints Vector{CartesianIndex{2}}
 
+"""
+```
+feature = Feature(keypoint, orientation = 0.0, scale = 0.0)
+```
+
+The `Feature` type has the keypoint, its orientation and its scale.
+"""
+immutable Feature
+    keypoint::Keypoint
+    orientation::Float64
+    scale::Float64
+end
+
+"""
+```
+features = Features(boolean_img)
+features = Features(keypoints)
+```
+
+Returns a `Vector{Feature}` of features generated from the `true` values in a boolean image or from a 
+list of keypoints.
+"""
+typealias Features Vector{Feature}
+
+Feature(k::Keypoint) = Feature(k, 0.0, 0.0)
+
+Feature(k::Keypoint, ori::Number) = Feature(k, ori, 0.0)
+
+Features(keypoints::Keypoints) = map(k -> Feature(k), keypoints)
+
+Features(img::AbstractArray) = Features(Keypoints(img))
+
+Keypoint(feature::Feature) = feature.keypoint
+
 function Keypoints(img::AbstractArray)
     r, c, _ = findnz(img)
     map((ri, ci) -> Keypoint(ri, ci), r, c)
 end
 
+Keypoints(features::Features) = map(f -> f.keypoint, features)        
+
+typealias OrientationPair Tuple{Int16, Int16}
+typealias OrientationWeights Tuple{Float16, Float16}
+typealias SamplePair Tuple{Float16, Float16}
+
 """
 ```
 distance = hamming_distance(desc_1, desc_2)

src/freak.jl

@@ -7,30 +7,28 @@ freak_params = FREAK([pattern_scale = 22.0])
 |----------|------|-------------|
 | **pattern_scale** | Float64 | Scaling factor for the sampling window | 
 """
-type FREAK{S, T, O} <: DescriptorParams
+type FREAK <: Params
     pattern_scale::Float64
-    pattern_table::Array{Array{S, 1}, 1}
-    smoothing_table::Array{T, 1}
-    orientation_weights::Array{O, 1}
+    pattern_table::Vector{Vector{SamplePair}}
+    smoothing_table::Vector{Vector{Float16}}
+    orientation_weights::Vector{OrientationWeights}
 end
 
-typealias OrientationPair Vector{Int}
-typealias SamplePair Vector{Float64}
-
 function FREAK(; pattern_scale::Float64 = 22.0)
     pattern_table, smoothing_table = _freak_tables(pattern_scale)
-    orientation_weights = OrientationPair[]
+    orientation_weights = OrientationWeights[]
     for o in freak_orientation_sampling_pattern
         offset_1 = pattern_table[1][o[1]]
         offset_2 = pattern_table[1][o[2]]
-        dy, dx = offset_1 - offset_2
+        dy = offset_1[1] - offset_2[1]
+        dx = offset_1[2] - offset_2[2]
         norm = (dx ^ 2 + dy ^ 2)
-        push!(orientation_weights, OrientationPair([round(Int, dy * 4096 / norm), round(Int, dx * 4096 / norm)]))
+        push!(orientation_weights, OrientationWeights((dy / norm, dx / norm)))
     end
     FREAK(pattern_scale, pattern_table, smoothing_table, orientation_weights)
 end
 
-function _freak_mean_intensity{T<:Gray}(int_img::AbstractArray{T, 2}, keypoint::Keypoint, offset::SamplePair, sigma::Float64)
+function _freak_mean_intensity{T<:Gray}(int_img::AbstractArray{T, 2}, keypoint::Keypoint, offset::SamplePair, sigma::Float16)
     y = keypoint[1] + offset[1]
     x = keypoint[2] + offset[2]
     if sigma < 0.5
@@ -45,15 +43,15 @@ function _freak_mean_intensity{T<:Gray}(int_img::AbstractArray{T, 2}, keypoint::
 end
 
 function _freak_orientation{T<:Gray}(int_img::AbstractArray{T, 2}, keypoint::Keypoint, pattern::Array{SamplePair}, 
-                                        orientation_weights::Array{OrientationPair}, sigmas::Array{Float64})
+                                        orientation_weights::Array{OrientationWeights}, sigmas::Array{Float16})
     direction_sum_y = 0.0
     direction_sum_x = 0.0
     for (i, o) in enumerate(freak_orientation_sampling_pattern)
         offset_1 = pattern[o[1]]
         offset_2 = pattern[o[2]]
         intensity_diff = _freak_mean_intensity(int_img, keypoint, offset_1, sigmas[o[1]]) - _freak_mean_intensity(int_img, keypoint, offset_2, sigmas[o[2]])
-        direction_sum_y += orientation_weights[i][1] * intensity_diff / 4096
-        direction_sum_x += orientation_weights[i][2] * intensity_diff / 4096
+        direction_sum_y += orientation_weights[i][1] * intensity_diff
+        direction_sum_x += orientation_weights[i][2] * intensity_diff
     end
     angle = atan2(direction_sum_y, direction_sum_x)
     scaled_angle = ceil(Int, (angle + pi) * freak_orientation_steps / (2 * pi))
@@ -62,18 +60,18 @@ end
 
 function _freak_tables(pattern_scale::Float64)
     pattern_table = Vector{SamplePair}[]
-    smoothing_table = Vector{Float64}[]
+    smoothing_table = Vector{Float16}[]
     for ori in 0:freak_orientation_steps - 1
         theta = ori * 2 * pi / freak_orientation_steps 
         pattern = SamplePair[]
-        sigmas = Float64[]
+        sigmas = Float16[]
         for (i, n) in enumerate(freak_num_circular_pattern)
             for circle_number in 0:n - 1
                 alt_offset = (pi / n) * ((i - 1) % 2)
                 angle = (circle_number * 2 * pi / n) + alt_offset + theta
 
-                push!(pattern, SamplePair([freak_radii[i] * sin(angle) * pattern_scale, 
-                                            freak_radii[i] * cos(angle) * pattern_scale]))
+                push!(pattern, SamplePair((freak_radii[i] * sin(angle) * pattern_scale, 
+                                            freak_radii[i] * cos(angle) * pattern_scale)))
                 push!(sigmas, freak_sigma[i] * pattern_scale)
             end
         end
@@ -88,10 +86,9 @@ function create_descriptor{T<:Gray}(img::AbstractArray{T, 2}, keypoints::Keypoin
     descriptors = BitArray{1}[]
     ret_keypoints = Keypoint[]
     window_size = ceil(Int, (freak_radii[1] + freak_sigma[1]) * params.pattern_scale) + 1
-    tl_lim = CartesianIndex(-window_size, -window_size)
-    br_lim = CartesianIndex(window_size, window_size)
+    lim = CartesianIndex(window_size, window_size)
     for k in keypoints
-        checkbounds(Bool, img, k + tl_lim) && checkbounds(Bool, img, k + br_lim) || continue
+        checkbounds(Bool, img, k - lim) && checkbounds(Bool, img, k + lim) || continue
         orientation = _freak_orientation(int_img, k, params.pattern_table[1], params.orientation_weights, params.smoothing_table[1])
         sampled_intensities = T[]
         for (i, p) in enumerate(params.pattern_table[orientation])

src/orb.jl

@@ -13,7 +13,7 @@ orb_params = ORB([num_keypoints = 500], [n_fast = 12], [threshold = 0.25], [harr
 | **levels** | Int | Number of levels in the gaussian pyramid.  See [`gaussian_pyramid`] in Images.jl | 
 | **sigma** | Float64 | Used for gaussian smoothing in each level of the gaussian pyramid.  See [`gaussian_pyramid`] in Images.jl | 
 """
-type ORB <: DescriptorParams
+type ORB <: Params
     num_keypoints::Int
     n_fast::Int
     threshold::Float64