image matching using SIFT

sift_match.rb

#!/usr/bin/env ruby

require 'rubygems'
require 'rmagick'

include Math

class KDTreeNode  
  attr_accessor :data, :left_child, :right_child

  def initialize
    @data = nil
    @left_child = nil
    @right_child = nil
  end

  def is_leaf
    return (self.left_child.nil? and self.right_child.nil?)
  end  
end

# Public:
#
# Examples
# 
#   kdtree = KDTree.construct_from_keypt(keypt)
# 
class KDTree
  attr_accessor :root_node

  def initialize(kpdata)
    def kdtree_build(kpdata, depth=0)
      if kpdata.size == 0 then
        return nil
      end

      # use only feature elements

      header = 4 # (x, y, scale, orientation)
      dim = kpdata.size - header
      axis = depth % dim
      sortedKpdata = kpdata.sort{|p,q| p[header + axis] <=> q[header + axis]}
      median = sortedKpdata.size / 2
      
      node = KDTreeNode.new
      node.data = sortedKpdata[median]
      node.left_child = kdtree_build(sortedKpdata[0, median], depth+1)
      node.right_child = kdtree_build(sortedKpdata[median+1..-1], depth+1)  

      return node
    end    
    
    @root_node = kdtree_build(kpdata, 0)
  end
  
  # class method
  def self.construct_from_keypoint(keypt)
    kdtree = self.new(keypt.data)
    return kdtree
  end

  def search_all(keypt)
    def search(qt, node, depth=0)
      if node.nil? then
        return nil
      end

      # if we have reached a leaf
      if node.is_leaf then
        return node.data
      end

      # this node is no leaf

      header = 4 
      dim = node.data.size - header
      axis = depth % dim

      # compare query point and point of current node in selected dimension
      # and figure out which subtree is farther than the other
      if qt[header + axis].to_i < node.data[header + axis].to_i then
        near_subtree = node.left_child
        far_subtree = node.right_child
      else
        near_subtree = node.right_child
        far_subtree = node.left_child
      end
      
      # recursively search through the tree until a leaf is found
      if near_subtree then
        if near_subtree.is_leaf then
          leaf = near_subtree.data
        else
          leaf = search(qt, near_subtree, depth+1)
        end
      else
        leaf = search(qt, far_subtree, depth+1)
      end

      # while unwinding the recursion, check whether there could be 
      # any points on the other side of splitting plane that are closer to the query point
      # than the current best
      if far_subtree &&
          square_distance(qt, leaf) > (qt[header + axis].to_i - node.data[header + axis].to_i)**2 then
        leaf = closer(qt, leaf, search(qt, far_subtree, depth+1))
      end

      return closer(qt, leaf, node.data)
    end

    # Finds an image feature's nearest neighbors in a kd tree
    pairs = Array.new
    for kpdata in keypt.data do
      pairs.push(search(kpdata, @root_node))
      p pairs.size
    end
    return pairs
  end

end

# Public: SIFT keypoint
# 
# Examples
# 
#   keypt = Keypoint.read_from_file(filename)
#
class Keypoint
  attr_accessor :num, :data

  def initialize(filename)
    @num = 0
    @data = Array.new
    IO.foreach(filename) do |s|
      @data.push(s.split(" ")[0..-1])
      @num += 1
    end
  end

  # class method
  def self.read_from_file(filename)
    keypt = self.new(filename)
    return keypt
  end

  def draw_circle(im, coord, scale)
    gc = Magick::Draw.new
    gc.fill_opacity(0)
    gc.stroke('red')

    origin_x = coord[0].to_i
    origin_y = coord[1].to_i
    scl = scale.to_f
    perim_x = (coord[0].to_f - scl * cos(0.0)).to_i
    perim_y = origin_y
    gc.circle(origin_x, origin_y, perim_x, perim_y)
    gc.draw(im)
  end
  private :draw_circle

  def draw_orientation(im, coord, scale, orientation)
    gc = Magick::Draw.new
    gc.stroke('blue')

    here_x = coord[0].to_i
    here_y = coord[1].to_i
    scl = scale.to_f
    ori = orientation.to_f
    len = scl * cos(0.0)  
    there_x = here_x + (len * cos(ori)).to_i
    there_y = here_y + (len * sin(ori)).to_i
    gc.line(here_x, here_y, there_x, there_y)
    gc.draw(im)
  end
  private :draw_orientation

  def plot_features(im)
    range = Range.new(0, @num, true)
    for i in range
      coord = @data[i][0..1]
      scale = @data[i][2]
      orientation = @data[i][3]
      # circle
      draw_circle(im, coord, scale)

      # orientation
      draw_orientation(im, coord, scale, orientation)
    end    
  end
end

def square_distance(ptA, ptB)
  header = 4
  dist = 0
  dims = Range.new(header, ptA.size, true)
  for dim in dims do
    dist += (ptA[dim].to_i - ptB[dim].to_i)**2
  end
  return dist
end

def closer(ref, ptA, ptB)
  distA = square_distance(ref, ptA)
  distB = square_distance(ref, ptB)
  return distA < distB ? ptA : ptB
end

def process_image(imagename, resultname, params="--edge-thresh 10 --peak-thresh 5")
  img = Magick::Image.from_blob( File.read(imagename) ).shift
  if img.format() != "PGM" 
    img.format = "PGM"
    File.open("tmp.pgm", "wb") {|f| f << img.to_blob}
    imagename = "tmp.pgm"
  end

  cmmd = "sift " + imagename + " --output=" + resultname + " " + params
  system(cmmd)
  puts "processed " + imagename + " to " + resultname
end

def plot_matches(im, im2, keypt, pairs)
  im3 = Magick::ImageList.new
  im3.push(im)
  im3.push(im2)
  im4 = im3.append(false)

  ofs_x = im.columns

  range = Range.new(0, pairs.size, true)
  scores = Array.new(pairs.size, 0)
  for i in range do 
    scores[i] = (square_distance(keypt.data[i], pairs[i]))
  end
  # plot X % of matched points
  x = pairs.size * 0.1
  print "plot ", x.to_i, " points\n"
  th_score = scores.sort[x.to_i]
  
  for i in range do
    next if scores[i] > th_score
    
    from_x = keypt.data[i][0].to_i
    from_y = keypt.data[i][1].to_i
    to_x = ofs_x + pairs[i][0].to_i
    to_y = pairs[i][1].to_i

    gc = Magick::Draw.new
    gc.stroke('green')
    gc.stroke_width(2)
    gc.line(from_x, from_y, to_x, to_y)
    gc.draw(im4)
  end
  im4.write("matched.jpg")
end

if __FILE__ == $0
  process_image('box.pgm', 'box.sift')
  keypt = Keypoint.read_from_file('box.sift')
  im = Magick::Image.read('box.pgm').first
  keypt.plot_features(im)
  im.write("box_feat.jpg")

  process_image('scene.pgm', 'scene.sift')
  keypt2 = Keypoint.read_from_file('scene.sift')
  im2 = Magick::Image.read('scene.pgm').first
  keypt2.plot_features(im2)
  im2.write("scene_feat.jpg")

  kdtree = KDTree.construct_from_keypoint(keypt2)
  pairs = kdtree.search_all(keypt)

  plot_matches(im, im2, keypt, pairs)
end

sadly this it does not say which program is used to generate the SIFT (line 216)

http://www.vlfeat.org/man/sift.html