--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, ranks l)
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--