Adds medianflow tracker

src/ImageTracking.jl

@@ -10,21 +10,25 @@ using StaticArrays
 include("core.jl")
 include("optical_flow.jl")
 include("haar.jl")
+include("tracker.jl")
 
 export
 
 	# main functions
-  optical_flow,
+	optical_flow,
+	init_tracker,
+	update_tracker,
 
-	# other functions
+  # other functions
 	haar_coordinates,
 	haar_features,
-
-	# other functions
 	polynomial_expansion,
 
 	# optical flow algorithms
 	LK,
-	Farneback
+	Farneback,
+
+	# tracking algorithms
+	TrackerMedianFlow
 
 end

src/medianflow_tracker/medianflow_tracker.jl

@@ -0,0 +1,152 @@
+mutable struct TrackerMedianFlow{I <: Int, F <: Float64} <: Tracker
+    #Initialized
+    points_in_grid::I
+    window_size::I
+    max_level::I
+    termination_criteria::I
+    window_size_for_ncc::I
+    max_median_of_displacement::F
+
+    #Uninitialized
+    model::TrackerModel
+
+    TrackerMedianFlow(points_in_grid::I = 15, window_size::I = 11, max_level::I = 4, termination_criteria::I = 20, window_size_for_ncc::I = 30, max_median_of_displacement::F = 10.0) where {I <: Int, F <: Float64} = new{I, F}(
+                      points_in_grid, window_size, max_level, termination_criteria, window_size_for_ncc, max_median_of_displacement)
+end
+
+function init_impl(tracker::TrackerMedianFlow{}, image::Array{T, 2}, bounding_box::MVector{4, Int}) where T
+    tracker.model = TrackerMedianFlowModel(image, bounding_box)
+end
+
+function update_impl(tracker::TrackerMedianFlow{}, image::Array{T, 2}) where T
+    prev_img = tracker.model.image
+    bounding_box = tracker.model.bounding_box
+    next_img = image
+
+    gray_prev_img = Gray.(prev_img)
+    gray_new_img = Gray.(next_img)
+
+    points_to_track = Vector{SVector{2, Int}}()
+    for i = 1:tracker.points_in_grid
+        for j = 1:tracker.points_in_grid
+            push!(points_to_track, SVector{2}(round(Int, bounding_box[1] + (i+0.5)*(bounding_box[3]/tracker.points_in_grid)), round(Int, bounding_box[2] + (j+0.5)*(bounding_box[4]/tracker.points_in_grid))))
+        end
+    end
+
+    flow, status, err = optical_flow(gray_prev_img, gray_new_img, LK(points_to_track, [SVector{2}(0.0,0.0)], tracker.window_size, tracker.max_level, false, tracker.termination_criteria))
+    new_points = points_to_track .+ flow
+    good_points = findn(status)
+    points_to_track = points_to_track[good_points]
+    new_points = new_points[good_points]
+
+    filter_status = MVector{length(good_points)}(trues(good_points))
+    filter_status = calc_fb_error(tracker, prev_img, next_img, points_to_track, new_points, filter_status)
+    filter_status = calc_ncc_error(tracker, prev_img, next_img, points_to_track, new_points, MVector{length(filter_status)}(filter_status))
+
+    good_points = findn(filter_status)
+    points_to_track = points_to_track[good_points]
+    new_points = new_points[good_points]
+
+    difference = new_points .- points_to_track
+    displacement, bounding_box = get_displacement(tracker, points_to_track, new_points, bounding_box)
+    final_disp = Array{SVector{2, Float64}, 1}()
+    for i = 1:length(difference)
+        push!(final_disp, difference[i] .- displacement)
+    end
+    displacements = sqrt.(dot.(final_disp,final_disp))
+    median_displacements = median(displacements)
+    assert(median_displacements < tracker.max_median_of_displacement)
+
+    tracker.model.image = image
+    tracker.model.bounding_box = bounding_box
+
+    return bounding_box
+end
+
+@inline function get_y(point::SVector{2, T}) where T
+    return point[1]
+end
+
+@inline function get_x(point::SVector{2, T}) where T
+    return point[2]
+end
+
+function get_displacement(tracker::TrackerMedianFlow{}, old_points::Array{SVector{2, Int}, 1}, new_points::Array{SVector{2, Float64}, 1}, bounding_box::MVector{4, Int})
+    new_center = MVector{2}(round(Int, (bounding_box[1] + bounding_box[3] + 1)/2), round(Int, (bounding_box[2] + bounding_box[4] + 1)/2))
+    new_bounding_box = MVector{4}(zeros(Int, 4))
+
+    if length(old_points) == 1
+        new_bounding_box[1] = bounding_box[1] + new_points[1][1] - old_points[1][1]
+        new_bounding_box[2] = bounding_box[2] + new_points[1][2] - old_points[1][2]
+        new_bounding_box[3] = new_bounding_box[1] + bounding_box[3] - bounding_box[1] + 1
+        new_bounding_box[4] = new_bounding_box[2] + bounding_box[4] - bounding_box[2] + 1
+        displacement = MVector{2}(new_points[1][1] - old_points[1][1], new_points[1][2] - old_points[1][2])
+        return displacement, new_bounding_box
+    end
+
+    location_shift = new_points .- old_points
+    y_shift = median(get_y.(location_shift))
+    x_shift = median(get_x.(location_shift))
+    new_center .+= MVector{2}(round(Int, y_shift), round(Int, x_shift))
+    displacement = MVector{2}(round(Int, y_shift), round(Int, x_shift))
+
+    n = length(old_points)
+    scale_buffer = MVector{round(Int, n*(n-1)/2), Float64}()
+    counter = 1
+    for i = 1:length(old_points)
+        for j = 1:i-1
+            new_change = norm(new_points[i] .- new_points[j])
+            old_change = norm(old_points[i] .- old_points[j])
+            scale_buffer[counter] = (old_change == 0.0)? 0.0:(new_change/old_change)
+            counter += 1
+        end
+    end
+    scale = median(scale_buffer)
+    new_bounding_box[1] = round(Int, new_center[1] - scale*(bounding_box[3] - bounding_box[1] + 1)/2)
+    new_bounding_box[2] = round(Int, new_center[2] - scale*(bounding_box[4] - bounding_box[2] + 1)/2)
+    new_bounding_box[3] = round(Int, new_bounding_box[1] + scale*(bounding_box[3] - bounding_box[1] + 1))
+    new_bounding_box[4] = round(Int, new_bounding_box[2] + scale*(bounding_box[4] - bounding_box[2] + 1))
+
+    return displacement, new_bounding_box
+end
+
+function calc_fb_error(tracker::TrackerMedianFlow, prev_img::Array{T, 2}, next_img::Array{T, 2}, old_points::Array{SVector{2, Int}, 1}, new_points::Array{SVector{2, Float64}, 1}, status::MVector{N, Bool}) where N where T
+    #TODO: Add Float64 support for input array in LK optical flow
+    flow, status, err = optical_flow(prev_img, next_img, LK(round(Int, new_points), [SVector{2}(0.0,0.0)], tracker.window_size, tracker.max_level, false, tracker.termination_criteria))
+
+    fb_error = norm.(old_points .- flow)
+    median_fb_error = median(fb_error)
+    status .= status .& (fb_error .<= median_fb_error)
+    return status
+end
+
+function get_patch(image::Array{T, 2}, patch_size::MVector{2, Int}, patch_center::SVector{2, I}) where {T, I <: Real}
+    roi_strat_corner = MVector{2, Int}(round(Int, patch_center[1] - patch_size[1]/2), round(Int, patch_center[2] - patch_size[2]/2))
+    patch_rect = MVector{4, Int}(roi_strat_corner[1], roi_strat_corner[2], roi_strat_corner[1] + patch_size[1] - 1, roi_strat_corner[2] + patch_size[2] - 1)
+
+    if patch_rect[1] >= 1 && patch_rect[2] >= 1 && patch_rect[3] <= size(image)[1] && patch_rect[4] <= size(image)[2]
+        patch = image[patch_rect[1]:patch_rect[3], patch_rect[2]:patch_rect[4]]
+    else
+        itp = interpolate(image, BSpline(Quadratic(Flat())), OnGrid())
+        etp = extrapolate(itp, zero(eltype(image)))
+        patch = etp[patch_rect[1]:patch_rect[3], patch_rect[2]:patch_rect[4]]
+    end
+
+    return patch
+end
+
+function calc_ncc_error(tracker::TrackerMedianFlow, prev_img::Array{T, 2}, next_img::Array{T, 2}, old_points::Array{SVector{2, Int}, 1}, new_points::Array{SVector{2, Float64}, 1}, status::MVector{N, Bool}) where N where T
+    ncc_error = MVector{length(old_points), Float64}()
+    for i = 1:length(old_points)
+        patch_1 = get_patch(prev_img, MVector{2}(tracker.window_size_for_ncc, tracker.window_size_for_ncc), old_points[i])
+        patch_2 = get_patch(next_img, MVector{2}(tracker.window_size_for_ncc, tracker.window_size_for_ncc), new_points[i])
+
+        temp = ncc(Float64.(patch_1), Float64.(patch_2))
+        #TODO: Use DataArrays.jl to handle missing data
+        ncc_error[i] = isnan(temp) ? 0 : temp
+    end
+
+    median_ncc_error = median(ncc_error)
+    status .= status .& (ncc_error .<= median_ncc_error)
+    return status
+end

src/tracker.jl

@@ -0,0 +1,29 @@
+abstract type Tracker end
+
+function init_tracker(tracker::Tracker, image::Array{T, 2}, bounding_box::MVector{4, Int}) where T
+    @assert bounding_box[1] >= 1 && bounding_box[1] <= size(image)[1]
+    @assert bounding_box[2] >= 1 && bounding_box[2] <= size(image)[2]
+    @assert bounding_box[3] >= 1 && bounding_box[3] <= size(image)[1]
+    @assert bounding_box[4] >= 1 && bounding_box[4] <= size(image)[2]
+
+    @assert bounding_box[1] < bounding_box[3]
+    @assert bounding_box[2] < bounding_box[4]
+
+    init_impl(tracker, image, bounding_box)
+end
+
+function update_tracker(tracker::Tracker, image::Array{T, 2}) where T
+    update_impl(tracker, image)
+end
+
+#---------------------
+# TRACKER COMPONENTS
+#---------------------
+
+include("tracker_model.jl")
+
+#------------------
+# IMPLEMENTATIONS
+#------------------
+
+include("medianflow_tracker/medianflow_tracker.jl")

src/tracker_model.jl

@@ -0,0 +1,10 @@
+abstract type TrackerModel end
+
+#-----------------------------
+# MODEL FOR MEDIANFLOW TRACKER
+#-----------------------------
+
+mutable struct TrackerMedianFlowModel{I <: Int} <: TrackerModel
+    image::Array{T, 2} where T
+    bounding_box::MVector{4, I}
+end

test/runtests.jl

@@ -6,3 +6,4 @@ using Base.Test
 
 include("optical_flow.jl")
 include("haar.jl")
+include("tracker.jl")

test/tracker.jl

@@ -0,0 +1,124 @@
+using Images, TestImages, StaticArrays
+
+function find_overlap(tracker_bounding_box::MVector{4, Int}, ground_truth_bounding_box::MVector{4, Int})
+	top = max(tracker_bounding_box[1], ground_truth_bounding_box[1])
+	bottom = min(tracker_bounding_box[3], ground_truth_bounding_box[3])
+	left = max(tracker_bounding_box[2], ground_truth_bounding_box[2])
+	right = min(tracker_bounding_box[4], ground_truth_bounding_box[4])
+
+	overlap_area = (bottom - top + 1)*(right - left + 1)
+	ground_truth_area = (ground_truth_bounding_box[3] - ground_truth_bounding_box[1] + 1)*(ground_truth_bounding_box[4] - ground_truth_bounding_box[2] + 1)
+
+	return (overlap_area/ground_truth_area)*100
+end
+
+@testset "MedianFlow" begin
+
+	#Simple tracking sequence
+	tr = TrackerMedianFlow()
+	a = MVector{4}(248,168,276,196)
+	img = Gray{Float64}.(testimage("lake"))
+	img[250:274, 170:194] .= 0
+	init_tracker(tr, img, a)
+
+	img = Gray{Float64}.(testimage("lake"))
+	img[253:277, 175:199] .= 0
+	bounding_box = update_tracker(tr, img)
+	overlap_percentage = find_overlap(bounding_box, MVector{4, Int}(251, 173, 279, 201))
+	@show overlap_percentage
+	@test overlap_percentage > 50
+
+	img = Gray{Float64}.(testimage("lake"))
+	img[254:278, 177:201] .= 0
+	bounding_box = update_tracker(tr, img)
+	overlap_percentage = find_overlap(bounding_box, MVector{4, Int}(253, 176, 279, 202))
+	@show overlap_percentage
+	@test overlap_percentage > 50
+
+	img = Gray{Float64}.(testimage("lake"))
+	img[257:281, 181:205] .= 0
+	bounding_box = update_tracker(tr, img)
+	overlap_percentage = find_overlap(bounding_box, MVector{4, Int}(256, 180, 282, 206))
+	@show overlap_percentage
+	@test overlap_percentage > 50
+
+	img = Gray{Float64}.(testimage("lake"))
+	img[262:286, 178:202] .= 0
+	bounding_box = update_tracker(tr, img)
+	overlap_percentage = find_overlap(bounding_box, MVector{4, Int}(261, 177, 287, 203))
+	@show overlap_percentage
+	@test overlap_percentage > 50
+
+	img = Gray{Float64}.(testimage("lake"))
+	img[265:289, 176:200] .= 0
+	bounding_box = update_tracker(tr, img)
+	overlap_percentage = find_overlap(bounding_box, MVector{4, Int}(264, 175, 290, 201))
+	@show overlap_percentage
+	@test overlap_percentage > 50
+
+	#Real data tracking sequence
+	tr = TrackerMedianFlow()
+	a = MVector{4}(81,150,120,173)
+	img = Gray{Float64}.(load("test_data/tracking/images/img00000.png"))
+	init_tracker(tr, img, a)
+
+	img = Gray{Float64}.(load("test_data/tracking/images/img00001.png"))
+	bounding_box = update_tracker(tr, img)
+	overlap_percentage = find_overlap(bounding_box, MVector{4, Int}(83, 150, 122, 173))
+	@show overlap_percentage
+	@test overlap_percentage > 50
+
+	img = Gray{Float64}.(load("test_data/tracking/images/img00002.png"))
+	bounding_box = update_tracker(tr, img)
+	overlap_percentage = find_overlap(bounding_box, MVector{4, Int}(84, 151, 123, 174))
+	@show overlap_percentage
+	@test overlap_percentage > 50
+
+	img = Gray{Float64}.(load("test_data/tracking/images/img00003.png"))
+	bounding_box = update_tracker(tr, img)
+	overlap_percentage = find_overlap(bounding_box, MVector{4, Int}(85, 151, 124, 174))
+	@show overlap_percentage
+	@test overlap_percentage > 50
+
+	img = Gray{Float64}.(load("test_data/tracking/images/img00004.png"))
+	bounding_box = update_tracker(tr, img)
+	overlap_percentage = find_overlap(bounding_box, MVector{4, Int}(87, 151, 126, 174))
+	@show overlap_percentage
+	@test overlap_percentage > 50
+
+	img = Gray{Float64}.(load("test_data/tracking/images/img00005.png"))
+	bounding_box = update_tracker(tr, img)
+	overlap_percentage = find_overlap(bounding_box, MVector{4, Int}(89, 152, 128, 175))
+	@show overlap_percentage
+	@test overlap_percentage > 50
+
+	img = Gray{Float64}.(load("test_data/tracking/images/img00006.png"))
+	bounding_box = update_tracker(tr, img)
+	overlap_percentage = find_overlap(bounding_box, MVector{4, Int}(90, 155, 129, 178))
+	@show overlap_percentage
+	@test overlap_percentage > 50
+
+	img = Gray{Float64}.(load("test_data/tracking/images/img00007.png"))
+	bounding_box = update_tracker(tr, img)
+	overlap_percentage = find_overlap(bounding_box, MVector{4, Int}(90, 156, 129, 179))
+	@show overlap_percentage
+	@test overlap_percentage > 50
+
+	img = Gray{Float64}.(load("test_data/tracking/images/img00008.png"))
+	bounding_box = update_tracker(tr, img)
+	overlap_percentage = find_overlap(bounding_box, MVector{4, Int}(91, 158, 130, 181))
+	@show overlap_percentage
+	@test overlap_percentage > 50
+
+	img = Gray{Float64}.(load("test_data/tracking/images/img00009.png"))
+	bounding_box = update_tracker(tr, img)
+	overlap_percentage = find_overlap(bounding_box, MVector{4, Int}(90, 160, 129, 183))
+	@show overlap_percentage
+	@test overlap_percentage > 50
+
+	img = Gray{Float64}.(load("test_data/tracking/images/img00010.png"))
+	bounding_box = update_tracker(tr, img)
+	overlap_percentage = find_overlap(bounding_box, MVector{4, Int}(87, 161, 126, 184))
+	@show overlap_percentage
+	@test overlap_percentage > 50
+end