Randomized Convex Hull

説明

乱択逐次構成法で点列から凸包を求めます。

• 三角形から始めて順次凸包を大きくしていく。
• 実行時間はO(nlogn)。

実装

→Ideone

class ConvexHull2D; end
class << ConvexHull2D
  private
  def norm(v)
    Math.sqrt(v[0]*v[0]+v[1]*v[1])
  end
  def dot(v0,v1)
    v0[0]*v1[0]+v0[1]*v1[1]
  end
  def det(v0,v1)
    v0[0]*v1[1]-v0[1]*v1[0]
  end
  def judge(p0,p1,p2)
    v0 = [p1[0]-p0[0],p1[1]-p0[1]]
    v1 = [p2[0]-p0[0],p2[1]-p0[1]]
    d = det(v0,v1)
    if d > 0
      2
    elsif d < 0
      -2
    else
      (f,n = dot(v0,v1),norm(v0);
      (0 < f && f < n*n)) ? 1 : -1
    end
  end
  def is_inner(xs, x)
    (0...xs.length).each do |i|
      si,di = i,i+1==xs.length ? 0 : i+1
      return false if judge(xs[si],xs[di],x) == 2
    end
    true
  end
  def prepare(xs)
    (2...xs.length).each do |i|
      d = judge(xs[0],xs[1],xs[i])
      if d == 2
        xs[2],xs[i] = xs[i],xs[2]
        return
      end
    end
    raise "invalid data"
  end
  def reconstruct(xs, x)
    acc = []
    t = judge(xs[-1],xs[0],x)
    prev_state = [2,1].include?(t)
    (0...xs.length).each do |i|
      si,di = i,i+1==xs.length ? 0 : i+1
      if [2,1].include?(judge(xs[si],xs[di],x))
        prev_state = true
        acc << xs[si]
      else
        if prev_state == true
          acc << xs[si] if t != 1
          acc << x
        end
        prev_state = false
      end
    end
    acc
  end
  public
  def calc(input)
    return nil if input.length < 3
    xs = input.shuffle
    prepare(xs)

    acc = xs[0..2]
    (3...xs.length).each do |i|
      unless is_inner(acc, xs[i])
        acc = reconstruct(acc, xs[i])
      end
    end
    acc
  end
end

if __FILE__ == $0
  def time(mes, &blk)
    start_time = Time.now
    blk.call
    puts "time [%s] : %.2f (sec)" % [mes, Time.now-start_time]
  end

  a = [[0,0],[2,0],[1,1],[0,2]]
  n = 1_000
  n.times do
    a << [Random.rand,Random.rand]
  end

  time "ConvexHull" do
    b = a.shuffle
    p ConvexHull2D.calc(b)
  end
end

実行結果

d3.jsライブラリを用いて生成。

所感

点と凸多角形の内外判定には「重心」と「点と直線の距離」を用いれば高速化できるかもしれない。