#!/usr/bin/env ruby

class Point < Struct.new(:x, :y)
  def self.random
    Point.new(rand, rand)
  end

  def dist(p2)
    xd = self.x - p2.x
    yd = self.y - p2.y
    Math.sqrt( xd**2 + yd**2 )
  end

  def -(p2)
    Point.new(self.x - p2.x, self.y - p2.y)
  end

  def +(v2)
    Point.new(self.x + v2.x, self.y + v2.y)
  end

  def *(scalar)
    Point.new(self.x * scalar, self.y * scalar)
  end

  def to_s
    "(#{x}, #{y})"
  end
end

class Circle < Struct.new(:center, :radius)
  def to_s
    "{center:#{center}, radius:#{radius}}"
  end
end

class AABB  # axis-aligned bounding box
  def initialize
    @p1 = Point.new
    @p2 = Point.new
  end
  
  def enclose(points)
    points.each {|point| add(point)}
  end
  
  def add(point)
    @p1.x = point.x if @p1.x.nil? or @p1.x > point.x
    @p1.y = point.y if @p1.y.nil? or @p1.y > point.y
    @p2.x = point.x if @p2.x.nil? or @p2.x < point.x
    @p2.y = point.y if @p2.y.nil? or @p2.y < point.y
  end

  def width
    @p2.x - @p1.x
  end

  def height
    @p2.y - @p1.y
  end
  
  def center
    Point.new(@p1.x + width/2.0, @p1.y + height/2.0)
  end
  
  def to_s  
    "#@p1,#@p2"
  end
end

# $SEED = 0

def generate_samples(n)
  # $SEED += 1 ; srand($SEED) ; puts "SEED=#$SEED"
  (1..n).map { Point.random }
end

SEARCH_EPSILON = 0.00001

# Binary search... might be better to use a fixed number of
# iterations, rather than the epsilon...
# All parameters are constant except low, high.
def search_closest_fit(outlier, vec, center, points, low, high)
  pivot = low + ((high - low) / 2.0)
  new_center = center + (vec * pivot)
  r = new_center.dist(outlier)
  too_far = points.any? {|pt| new_center.dist(pt) > r}
  if too_far
    high = pivot
  else
    low = pivot
  end
  done = (high - low) < SEARCH_EPSILON
  if done
    [new_center, r]
  else
    search_closest_fit(outlier, vec, center, points, low, high)
  end
end

# Algorithm used: First enclose the points in a simple
# axis-aligned bounding box, and choose the center of the
# box as the initial center of the circle. Then, binary
# search to determine how far we can move the center point
# toward the furthest outlying point, while keeping the
# rest of the points still inside the circle.
def encircle(points)        # takes array of Point objects
  bbox = AABB.new
  bbox.enclose(points)
  c = bbox.center
  points = points.sort_by {|pt| -pt.dist(c)}
  outlier = points.shift
  if points.empty?
    r = 0.0  # c.dist(outlier)
  else
    vec = outlier - c
    c, r = search_closest_fit(outlier, vec, c, points, 0.0, 1.0)
  end
  Circle.new(c, r)
end


if $0 == __FILE__

require 'test/unit'

class TestEncircle < Test::Unit::TestCase

  EPSILON = 0.001

  def assert_in_delta(expected, actual, delta=EPSILON)
    super(expected, actual, delta)
  end

  def test_single_point
    circ = encircle( [ Point.new(0.8, 0.3) ] )
    assert_in_delta( 0.8, circ.center.x )
    assert_in_delta( 0.3, circ.center.y )
    assert_in_delta( 0.0, circ.radius )
  end

  def test_two_points
    p1, p2 = Point.new(0.25, 0.25), Point.new(0.75, 0.75)
    circ = encircle( [p1, p2] )
    assert_in_delta( 0.5, circ.center.x )
    assert_in_delta( 0.5, circ.center.y )
    assert_in_delta( circ.center.dist(p1), circ.radius )
  end

end

end