```On Tue, Feb 19, 2008 at 09:32:16AM +0900, Lionel Bouton wrote:
> Philipp Hofmann wrote:
>> 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.
>
> Wow, that was awfully fast for the random distribution, it takes more
> time to display the result than computing it even on 10000 points !
>
> On badly shaped point clouds it suffers greatly, I couldn't get it to
> compute the circle for 1000 points if I feed it with :

What you call a 'badly shaped point clouds' is in fact an approximate
circle. It is true that this is the worst case for convex hull based
approaches because in this case the convex hull has lots of points.

It would be nice to have an algorithm that produces random
distibutions that form arbitrary polygons. But I guess that's a task
for another quiz.

As I mentioned in one of my comments in the code, right before it goes
into brute force, I prefer using Clarkson's probabilistic aproach here.

So here is my improved version, still based on convex hulls but also
able to deal with circular shaped distributions.

I confident it doesn't do any harm to your poor laptop this time. ;)

> ...

g phil

class Point < Struct.new(:x, :y)

# this distribution is the worst case for convex hull
# based algorithms, because it leads to an approximate
# circle, therefore the convex hull has lots of points
def self.random
#Point.new(rand, rand)
ro = rand / 2
theta = Math::PI * 2 * rand
Point.new(ro * Math::cos(theta) + 0.5, ro * Math::sin(theta) + 0.5)
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

def to_s
end
end

def generate_samples(n)
(1..n).map { Point.random }
end

# ----- clarkson's probabilistic approach

def encircle(points)
# this in arbitrary set threshold
return previous_encircle(points) if points.size < 100
circle = nil
pot = {}
points.each { |p| pot[p] = 1 }
finished = false
until finished
# select 13 random points
selection = []
sum = pot.values.inject { |sum, s| sum += s }
13.times { selection << rand(sum) }
selection.uniq!
while selection.size < 13
selection << rand(sum)
selection.uniq!
end
selection.sort!
thirteen = []
sum = 0
threshold = selection.shift
pot.each do |p, s|
sum += s
if threshold <= sum
thirteen << p
break if selection.empty?
threshold = selection.shift
end
end
# use previous method on selection of thirteen
circle = previous_encircle(thirteen)
# double chance of every outsider
finished = true
pot.each do |p, chance|
if !thirteen.include?(p) && p.distance(circle.center) > circle.radius
pot[p] = chance * 2
finished = false
end
end
end
circle
end

# --- convex hull based approach

def previous_encircle(points)

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

# firstly, we calculate the convex hull
# this will reduce points to calculate
# especially for larger sets
hull = ConvexHull::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 = hull.clone
until list_of_points.empty?
p1 = list_of_points.shift
list_of_points.each do |p2|
distances0[[p1, p2]] = p1.distance(p2)
end
end
max0 = distances0.values.max
p0, p1 = distances0.invert[max0]
circle0 = Circle.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] = circle0.center.distance(point)
end
max1 = distances1.values.max
# if the farest point from center is
# closer then the radius we're done
return circle0 if max1 <= circle0.radius

p2 = distances1.invert[max1]
circle1 = SOLE::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] = circle1.center.distance(point)
end
max2 = distances2.values.max
# if the farest point from center is in the circle we're done
return circle1 if max2 <= circle1.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 = SubSet::subsets(hull, 3)
solutions = []
sets_of_three.each do |subset|
circle2 = SOLE::circle_of(*subset)
distances3 = {}
hull.each do |p3|
distances3[p3] = circle2.center.distance(p3) unless subset.include?(p3)
end
max3 = distances3.values.max
solutions << circle2 if max3 <= circle2.radius
end

circle3 = solutions.sort_by { |circle| circle.radius }.first
return circle3

end

# --- calculate subsets

module SubSet

# i admit this could definitly be improved
def self.subsets(set, ranks=nil)
ranks ||= 0..set.size
ranks = ranks..ranks unless ranks.is_a?(Range)
sets = all_subsets(set)
sets.select { |s| ranks.include?(s.size) }
end

def self.all_subsets(set)
return [[]] if set.empty?
set = set.clone
first = set.shift
sets = all_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 = points.sort_by { |p| p.x }
left = lop.shift
right = lop.pop
lower, upper = [left], [left]
lower_hull, upper_hull = [], []
det_func = determinant_function(left, right)
until lop.empty?
p = lop.shift
( det_func.call(p) < 0 ? lower : upper ) << p
end
lower << right
until lower.empty?
lower_hull << lower.shift
while (lower_hull.size >= 3) &&
!convex?(lower_hull.last(3), true)
last = lower_hull.pop
lower_hull.pop
lower_hull << last
end
end
upper << right
until upper.empty?
upper_hull << upper.shift
while (upper_hull.size >= 3) &&
!convex?(upper_hull.last(3), false)
last = upper_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 = list_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 = sole_matrix_for_circle(three_points)
solve_sole!(m)
x = m[1].last
y = m[2].last
r = Math.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 = 0
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 += 1
end
matrix
end

# solve system of linear equations
def self.solve_sole!(matrix)
num_equations = matrix.size
num_variables = matrix.first.size
(num_variables-1).times do |j|
q = matrix[j][j]
if q == 0
for i in j+1..num_equations-1 do
if matrix[i][j] != 0
for k in 0..num_variables-1 do
matrix[j][k] += matrix[i][k]
end
q = matrix[j][j]
break
end
end
end
if q != 0
for k in 0..num_variables-1
matrix[j][k] = matrix[j][k] / q
end
end
for i in 0..num_equations-1
if i != j
q = matrix[i][j]
for k in 0..num_variables-1
matrix[i][k] -= matrix[j][k] * q
end
end
end
end
end

end

# ----- visualize

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

\$gc.stroke_color(color.to_s) if color
end

def draw_point(center, color=false)
\$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=false)
\$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=false)
first = points.shift
prev = first
until points.empty?
p = points.shift
draw_line(prev, p, color)
prev = p
end
draw_line(prev, first, color)
end

# ----- if executed directly

if __FILE__ == \$0

VISUALIZE = true

if VISUALIZE
require 'RMagick'
\$size = 600
\$gc = Magick::Draw.new
\$gc.fill_opacity(0)
\$gc.stroke_width(2)
\$gc.stroke_opacity(1)
end

samples = generate_samples(ARGV[0].to_i)
circle = encircle(samples)
puts circle

if VISUALIZE
samples.each { |p| draw_point(p, :black) }