Bounded Hough Transform

.gitignore

*.png
*.avi
.*.un~
.*.swp

bht.rb

#!/usr/bin/env ruby
require 'hornetseye_ffmpeg'
require 'hornetseye_rmagick'
require 'hornetseye_xorg'
include Hornetseye
class Node
  GRAD_SIGMA = 1
  COV_SIGMA = 1
  K = 0.05
  def corners
    x = gauss_gradient GRAD_SIGMA, 0
    y = gauss_gradient GRAD_SIGMA, 1
    a = (x * x).gauss_blur COV_SIGMA
    b = (y * y).gauss_blur COV_SIGMA
    c = (x * y).gauss_blur COV_SIGMA
    tr = a + b
    det = a * b - c * c
    det - tr * tr * K
  end
  THRESHOLD = 10000
  def nms(threshold = THRESHOLD)
    self >= dilate.major(threshold)
  end
end
raise "Syntax: #{$0} <input video>" if ARGV.size != 1
name = ARGV.first
input = AVInput.new name
output = AVOutput.new "test#{name}", 100000, input.width, input.height, 5.quo(2)
px = 192
py = 144
pw = 0
w2 = 40
h2 = 40
img = MultiArray.load_ubyte 'obj.png'
rect = [px - w2 .. px + w2 , py - h2 .. py + h2]
obj = img[*rect]
obj_corners = obj.corners.nms
corner_x = lazy(*obj.shape) { |x, y| x }.mask obj_corners
corner_y = lazy(*obj.shape) { |x, y| y }.mask obj_corners
xr2 = 5
yr2 = 5
wr2 = 4
hxr2 = 2 * xr2
hyr2 = 2 * yr2
dx = MultiArray.byte(2 * wr2 + 1, *obj.shape).fill!
dy = MultiArray.byte(2 * wr2 + 1, *obj.shape).fill!
valid = MultiArray.bool(2 * wr2 + 1, *obj.shape).fill!
warp_x = lazy(2 * xr2 + 1, 2 * xr2 + 1) do |i, j|
  i - xr2
end
warp_y = lazy(2 * xr2 + 1, 2 * xr2 + 1) do |i, j|
  j - yr2
end
ramp = lazy(2 * xr2 + 1, 2 * yr2 + 1) do |i, j|
  (i - xr2).abs.major((j - yr2).abs).minor xr2
end.normalise 1 .. 0
for k in -wr2 .. wr2
  w = -k * Math::PI / 180
  cw = Math.cos w
  sw = Math.sin w
  for j in -yr2 .. yr2
    for i in -xr2 .. xr2
      for c in 0 ... corner_x.size
        x = corner_x[c] - w2 + i
        y = corner_y[c] - h2 + j
        xs = cw * x - sw * y - (x - i)
        ys = sw * x + cw * y - (y - j)
        if (-w2 .. w2).member?(x) and (-h2 .. h2).member?(y)
          valid[k + wr2, x + w2, y + h2] = true
          dx[k + wr2, x + w2, y + h2] = xs.round
          dy[k + wr2, x + w2, y + h2] = ys.round
        end
      end
    end
  end
end
o_x = corner_x
o_y = corner_y
X11Display.show 2 * input.width do
  puts "px = #{px}, py = #{py}, pw = #{pw}"
  begin
    img = input.read_ubyte
  rescue RuntimeError
    exit
  end
  w = pw * Math::PI / 180
  cw = Math.cos w
  sw = Math.sin w
  warp_xs = finalise do |i, j, k|
    (o_x[k] - w2) * cw - (o_y[k] - h2) * sw + warp_x[i, j] + px
  end.round.to_int
  warp_ys = finalise do |i, j, k|
    (o_x[k] - w2) * sw + (o_y[k] - h2) * cw + warp_y[i, j] + py
  end.round.to_int
  #corners = img.corners.nms
  #box = img.warp warp_xs, warp_ys
  #cbox = corners.warp warp_xs, warp_ys
  box = (img.warp(warp_xs, warp_ys).unroll * ramp).roll
  cbox = MultiArray.dfloat *box.shape
  (0 ... o_x.size).each do |k|
    patch = box[k].corners
    cbox[k] = patch >= [patch.max, 100000].max
  end
  corner_x = warp_xs.mask cbox
  corner_y = warp_ys.mask cbox
  hist = MultiArray.sint(2 * wr2 + 1, 2 * hxr2 + 1, 2 * hyr2 + 1).fill!
  for c in 0 ... corner_x.size
    x = corner_x[c] - px
    y = corner_y[c] - py
    xs = cw * x + sw * y
    ys = -sw * x + cw * y
    x = xs.round
    y = ys.round
    if (-w2 .. w2).member?(x) and (-h2 .. h2).member?(y)
      for k in -wr2 .. wr2
        if valid[k + wr2, x + w2, y + h2]
          hist[k + wr2,
               dx[k + wr2, x + w2, y + h2] + hxr2,
               dy[k + wr2, x + w2, y + h2] + hyr2] += 1
        end
      end
    end
  end
  hist = hist.gauss_blur 1.0
  if hist.max > 0
    m = argmax { |k, i, j| hist[k, i, j] }
    sx = m[1] - hxr2
    sy = m[2] - hyr2
    w = pw * Math::PI / 180
    cw = Math.cos w
    sw = Math.sin w
    px += (cw * sx - sw * sy).round
    py += (sw * sx + cw * sy).round
    pw += m[0] - wr2
  end
  result = img.to_ubytergb
  for i in 0 ... corner_x.size
    result[corner_x[i], corner_y[i]] = RGB 0, 0, 255
  end
  for i in 0 ... o_x.size
    x = o_x[i] - w2
    y = o_y[i] - h2
    w = pw * Math::PI / 180
    cw = Math.cos w
    sw = Math.sin w
    xs = cw * x - sw * y
    ys = sw * x + cw * y
    x = (xs + px).round
    y = (ys + py).round
    if (1 ... img.width - 1).member?(x) and (1 ... img.height - 1).member?(y)
      result[x - 1 .. x + 1, y - 1 .. y + 1] |= RGB 255, 0, 0
    end
  end
  ox = 2 * xr2 + 1
  oy = 2 * yr2 + 1
  for i in 0 ... o_x.size
    k = 2 * i + 1
    result[k * ox ... (k + 1) * ox, oy ... 2 * oy] =
      box[i] / 2 + 64
    result[k * ox ... (k + 1) * ox, 3 * oy ... 4 * oy] =
      cbox[i].conditional 255, 128
  end
  output.write result
  result
end

Makefile

test: testw.avi testx.avi testy.avi

testx.avi: bht.rb x.avi obj.png
	HORNETSEYE_PRELOAD_CACHE=1 ruby bht.rb x.avi

testy.avi: bht.rb y.avi obj.png
	HORNETSEYE_PRELOAD_CACHE=1 ruby bht.rb y.avi

testw.avi: bht.rb w.avi obj.png
	HORNETSEYE_PRELOAD_CACHE=1 ruby bht.rb w.avi

obj.png: x00.png
	cp $< $@

x.avi: x00.png
	mencoder mf://x*.png -ovc x264 -x264encopts crf=23 -o $@

y.avi: y00.png
	mencoder mf://y*.png -ovc x264 -x264encopts crf=23 -o $@

w.avi: w00.png
	mencoder mf://w*.png -ovc x264 -x264encopts crf=23 -o $@

x00.png: x.ini x.pov
	povray x

y00.png: y.ini y.pov
	povray y

w00.png: w.ini w.pov
	povray w

clean:
	rm -f *.png *.avi

w.ini

;Width=768
;Height=576
Width=384
Height=288

Display=On

Antialias=On
Jitter=On
Sampling_Method=2
Antialias_Threshold=0.5
Antialias_Depth=2
Test_Abort_Count=100
Input_File_Name=w.pov
Output_File_Type=N

Initial_Frame = 0
Final_Frame = 49
Initial_Clock = 0.0
Final_Clock = 1.0
Cyclic_Animation = true

Continue_Trace=Off
;Continue_Trace=On
Test_Abort=On

Display=Off
;Display=On

Pause_when_Done=Off

w.pov

#include "colors.inc"

global_settings {
  assumed_gamma 2.2
  max_trace_level 5
  ambient_light White
}

background { color Black }

#declare dist = 10000;

camera {
   location < 0.0, 0.0, -dist >
   look_at 0
   angle atan(20.0 / dist) * 180.0 / pi
}

polygon {
   4,
   < -1, -1 >, < 1, -1 >, < 1, 1 >, < -1, 1 >
   texture {
      finish { ambient 1 diffuse 0 }
      pigment { color White }
  }
  rotate z * clock * 90
}

x.ini

;Width=768
;Height=576
Width=384
Height=288

Display=On

Antialias=On
Jitter=On
Sampling_Method=2
Antialias_Threshold=0.5
Antialias_Depth=2
Test_Abort_Count=100
Input_File_Name=x.pov
Output_File_Type=N

Initial_Frame = 0
Final_Frame = 49
Initial_Clock = 0.0
Final_Clock = 1.0
Cyclic_Animation = true

Continue_Trace=Off
;Continue_Trace=On
Test_Abort=On

Display=Off
;Display=On

Pause_when_Done=Off

x.pov

#include "colors.inc"

global_settings {
  assumed_gamma 2.2
  max_trace_level 5
  ambient_light White
}

background { color Black }

#declare dist = 10000;

camera {
   location < 0.0, 0.0, -dist >
   look_at 0
   angle atan(20.0 / dist) * 180.0 / pi
}

polygon {
   4,
   < -1, -1 >, < 1, -1 >, < 1, 1 >, < -1, 1 >
   texture {
      finish { ambient 1 diffuse 0 }
      pigment { color White }
  }
  translate 3 * x * clock
}

y.ini

;Width=768
;Height=576
Width=384
Height=288

Display=On

Antialias=On
Jitter=On
Sampling_Method=2
Antialias_Threshold=0.5
Antialias_Depth=2
Test_Abort_Count=100
Input_File_Name=y.pov
Output_File_Type=N

Initial_Frame = 0
Final_Frame = 49
Initial_Clock = 0.0
Final_Clock = 1.0
Cyclic_Animation = true

Continue_Trace=Off
;Continue_Trace=On
Test_Abort=On

Display=Off
;Display=On

Pause_when_Done=Off

y.pov

#include "colors.inc"

global_settings {
  assumed_gamma 2.2
  max_trace_level 5
  ambient_light White
}

background { color Black }

#declare dist = 10000;

camera {
   location < 0.0, 0.0, -dist >
   look_at 0
   angle atan(20.0 / dist) * 180.0 / pi
}

polygon {
   4,
   < -1, -1 >, < 1, -1 >, < 1, 1 >, < -1, 1 >
   texture {
      finish { ambient 1 diffuse 0 }
      pigment { color White }
  }
  translate 3 * y * clock
}