--4Ckj6UjgE2iN1+kY
Content-Type: text/plain; charset=us-ascii
Content-Disposition: inline


Here is my solution. Although it uses a brute force approach if the
first attempts (driven by obvious constraints) of determining the
circle fail, it scales on larger sets, because the given equal
distibution of random points form an approximate square on these
sets. In that case the convex hull which is used for the brute force
has very few points.

I also did some visualization with RMagick.

g phil
--4Ckj6UjgE2iN1+kY
Content-Type: text/plain; charset=us-ascii
Content-Disposition: attachment; filename="smallest_enclosing_circle.rb"

class Point < Struct.new(:x, :y) 
  def self.random
    Point.new(rand, rand)
  end
  
  def to_s 
    "(#{x}, #{y})" 
  end

  def distance(p)
    Math.sqrt((x-p.x)**2+(y-p.y)**2)
  end

  def middle(p)
    Point.new((self.x+p.x)/2, (self.y+p.y)/2)
  end
end
 
class Circle < Struct.new(:center, :radius)
  def to_s 
    "{center:#{center}, radius:#{radius}}" 
  end
end
 
def generate_samples(n)
  (1..n).map { Point.random }
end

def encircle(points)

  # only one point? then we're done
  return Circle.new(points.first, 0) if points.size

  # firstly, we calculate the convex hull
  # this will reduce points to calculate
  # especially for larger sets
  hull  onvexHull::convex_hull(points)
  draw_hull(hull.clone, :lightgray) if VISUALIZE

  # find the two spots on the hull with the larges distance
  # these will make a good starting point
  distances0  }
  list_of_points  ull.clone
  until list_of_points.empty?
    p1  ist_of_points.shift
    list_of_points.each do |p2|
      distances0[[p1, p2]]  1.distance(p2)
    end
  end
  max0  istances0.values.max
  p0, p1  istances0.invert[max0]
  circle0  ircle.new(p0.middle(p1), max0/2)
  hull.delete(p0)
  hull.delete(p1)

  # we're done if the hull only had two points
  return circle0 if hull.empty?
  # find the point that has the largest
  # distance to the center of circle0
  distances1  } 
  hull.each do |point|
    distances1[point]  ircle0.center.distance(point)
  end
  max1  istances1.values.max
  # if the farest point from center is
  # closer then the radius we're done
  draw_circle(circle0.center, circle0.radius, :green) if VISUALIZE
  return circle0 if max1 < ircle0.radius

  p2  istances1.invert[max1]
  circle1  OLE::circle_of(p0, p1, p2)
  hull.delete(p2)
  # find the point that has the largest
  # distance to the center of circle1
  distances2  }
  hull.each do |point|
    distances2[point]  ircle1.center.distance(point)
  end
  max2  istances2.values.max
  # if the farest point from center is in the circle we're done
  draw_circle(circle1.center, circle1.radius, :blue) if VISUALIZE
  return circle1 if max2 < ircle1.radius
  
  # ok, now we fall back to brute force on the convex hull
  # we can do this because the random distribution given
  # leads to an approximate square and therefore to convex hulls
  # with very few points, if this weren't the case we should
  # uses clarksons probabilistic approach here, especcially
  # for larger sets of points.
  hull.concat([p0, p1, p2])
  sets_of_three  ubSet::subsets(hull, 3)
  solutions  ]
  sets_of_three.each do |subset|
    circle2  OLE::circle_of(*subset)
    distances3  }
    hull.each do |p3|
      distances3[p3]  ircle2.center.distance(p3) unless subset.include?(p3)
    end
    max3  istances3.values.max
    solutions << circle2 if max3 < ircle2.radius
  end

  circle3  olutions.sort_by { |circle| circle.radius }.first
  draw_circle(circle3.center, circle3.radius, :red) if VISUALIZE
  return circle3

end

# --- calculate subsets

module SubSet
  
  # i admit this could definitly be improved
  def self.subsets(set, ranksl)
    ranks || ..set.size
    ranks  anks..ranks unless ranks.is_a?(Range)
    sets  ll_subsets(set)
    sets.select { |s| ranks.include?(s.size) }
  end
  
  def self.all_subsets(set)
    return [[]] if set.empty?
    set  et.clone
    first  et.shift
    sets  ll_subsets(set)
    sets.concat(sets.collect { |s| [first] + s })
    return sets
  end
  
end

# --- calculate the convex hull

module ConvexHull

  # after graham & andrew
  def self.convex_hull(points)
    lop  oints.sort_by { |p| p.x }
    left  op.shift
    right  op.pop
    lower, upper  left], [left]
    lower_hull, upper_hull  ], []
    det_func  eterminant_function(left, right)
    until lop.empty?
      p  op.shift
      ( det_func.call(p) < 0 ? lower : upper ) << p
    end
    lower << right
    until lower.empty?
      lower_hull << lower.shift
      while (lower_hull.size > ) &&
          !convex?(lower_hull.last(3), true)
        last  ower_hull.pop
        lower_hull.pop
        lower_hull << last
      end
    end
    upper << right
    until upper.empty?
      upper_hull << upper.shift
      while (upper_hull.size > ) &&
          !convex?(upper_hull.last(3), false)
        last  pper_hull.pop
        upper_hull.pop
        upper_hull << last
      end
    end
    upper_hull.shift
    upper_hull.pop
    lower_hull + upper_hull.reverse
  end
  
  private

  def self.determinant_function(p0, p1)
    proc { |p| ((p0.x-p1.x)*(p.y-p1.y))-((p.x-p1.x)*(p0.y-p1.y)) }
  end
  
  def self.convex?(list_of_three, lower)
    p0, p1, p2  ist_of_three
    (determinant_function(p0, p2).call(p1) > 0) ^ lower
  end
  
end

# --- calculate a circle out of three given points

module SOLE
  
  def self.circle_of(*three_points)
    m  ole_matrix_for_circle(three_points)
    solve_sole!(m)
    x  [1].last
    y  [2].last
    r  ath.sqrt(x*x + y*y - m[0].last)
    Circle.new(Point.new(x, y), r)
  end
  
  private

  def self.sole_matrix_for_circle(points)
    matrix  ]
    c  
    points.each do |p|
      matrix[c]  1]
      matrix[c] << -2 * p.x
      matrix[c] << -2 * p.y
      matrix[c] << -p.x**2 - p.y**2
      c + 
    end
    matrix
  end
  
  # solve system of linear equations
  def self.solve_sole!(matrix) 
    num_equations  atrix.size
    num_variables  atrix.first.size
    (num_variables-1).times do |j|
      q  atrix[j][j]
      if q 0
        for i in j+1..num_equations-1 do
          if matrix[i][j] ! 
            for k in 0..num_variables-1 do
              matrix[j][k] + atrix[i][k]
            end
            q  atrix[j][j]
            break
          end
        end
      end
      if q ! 
        for k in 0..num_variables-1
          matrix[j][k]  atrix[j][k] / q
        end
      end
      for i in 0..num_equations-1
        if i ! 
          q  atrix[i][j]
          for k in 0..num_variables-1
            matrix[i][k] - atrix[j][k] * q
          end
        end
      end
    end
  end
  
end

# ----- visualize

def scale(value)
  $size / 6 + value * $size / 1.5
end

def draw_circle(center, radius, colorů±se)
  $gc.stroke_color(color.to_s) if color
  $gc.circle(scale(center.x), scale(center.y), scale(center.x), scale(center.y+radius))
end

def draw_point(center, colorů±se)
  $gc.stroke_color(color.to_s) if color
  $gc.circle(scale(center.x), scale(center.y), scale(center.x), scale(center.y)+2)
end

def draw_line(p0, p1, colorů±se)
  $gc.stroke_color(color.to_s) if color
  $gc.line(scale(p0.x), scale(p0.y), scale(p1.x), scale(p1.y))
end

def draw_hull(points, colorů±se)
  first  oints.shift
  prev  irst
  until points.empty?
    p  oints.shift
    draw_line(prev, p, color)
    prev  
  end
  draw_line(prev, first, color)
end

# ----- if executed directly

if __FILE__ $0
  
  VISUALIZE  rue

  if VISUALIZE
    require 'RMagick'
    $size  00
    $gc  agick::Draw.new
    $gc.fill_opacity(0)
    $gc.stroke_width(2)
    $gc.stroke_opacity(1)
  end
  
  samples  enerate_samples(ARGV[0].to_i)
  circle  ncircle(samples)
  puts circle
  
  if VISUALIZE
    samples.each { |p| draw_point(p, :black) }
    canvas  agick::ImageList.new
    canvas.new_image($size, $size)
    $gc.draw(canvas)
    canvas.display
  end
  
end

--4Ckj6UjgE2iN1+kY--