Hi,

here’s my solution. I initially thought that the best circle would be

either defined by a diameter (2 points) or one combinations of 3 points

through which it passes. The idea was that given three points defining a

circle a fourth point is either in the circle or defines a circle with 2

from the 3 originals that includes the remaining point.

This would be neat as I wasn’t so happy about a gradient walking

solution. In this case it can be proven that you can deterministically

get the best solution by gradient walking (I don’t think there are local

maximums) but it’s not the general case. Unfortunately it can be proven

that such a circle can encircle all points, but not that it’s the

smallest (and my experiments proved that it isn’t). Finally this is even

slower than my proposed solution below, …

So here’s a more brute force method : minimizing the maximum distance

with some optimizations (heuristic for initial values and a not so naive

2D-space walker). There’s probably some more optimizations for walking

the gradient but I’m happy with the current speed (~10s for 10000 points

on my 1.06GHz Core Duo).

This seems similar in spirit to Frank’s solution, but less clean from

the Math PoV and more optimized for speed.

It outputs the solution and creates a gnuplot plotting file for you to

verify (you have the path used by the algorithm to find the center). I’m

no Gnuplot expert so unfortunately fancy colors are out and crappy

legends in…

Good reading,

Lionel

=============================================

include Math

# Uses 10 points by default

POINT_COUNT = (ARGV[0].to_i == 0) ? 10 : ARGV[0].to_i

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

def self.random

Point.new(rand, rand)

end

def self.generate_samples(n)

(1…n).map { self.random }

end

def to_s

“(#{x}, #{y})”

end

end

class PotentialCenter < Struct.new(:x, :y, :points)

# Square distance with a point

def square_distance(point)

((point.x - x) ** 2) + ((point.y - y) ** 2)

end

# Maximum square distance with my points

# It’s cached because it’s called several times and can be very costly

def max_square_distance

@sq_dist ||= (points.map { |point| square_distance(point) }.max)

end

# Compute the best move with given distance (ro)

# and angles (thetas)

# returns nil if none is better than the current position

def best_move(ro, thetas)

potential_moves = thetas.map { |theta|

new_x = x + (ro * cos(theta))

new_y = y + (ro * sin(theta))

PotentialCenter.new(new_x,new_y,points)

}

choose_move_1(potential_moves)

end

# Dumber, faster

def choose_move_1(potential_moves)

potential_moves.detect { |move|

max_square_distance > move.max_square_distance

}

end

# Look for the shortest path, slower (~1.4x on average with many

points)

def choose_move_2(potential_moves)

potential_moves.reject! { |move|

max_square_distance <= move.max_square_distance

}

potential_moves.min { |c1,c2|

c1.max_square_distance <=> c2.max_square_distance

}

end

# Check that the given offset is big enough to move

# the current position

# (floating point rounding errors prevent infinite precision…)

def can_be_moved_by(offset)

((x + offset) != x) || ((y + offset) != y)

end

def to_s

“(#{x}, #{y})”

end

end

class Circle < Struct.new(:center, :radius)

def to_s

“{center:#{center}, radius:#{radius}}”

end

end

def chose_original_center_and_move_amount(points)

# Start with a good potential (around the middle of the point cloud)

max_x = points.map{|p| p.x}.max

min_x = points.map{|p| p.x}.min

x = (max_x + min_x) / 2

max_y = points.map{|p| p.y}.max

min_y = points.map{|p| p.y}.min

y = (max_y + min_y) / 2

center = PotentialCenter.new(x, y, points)

# The original displacement value is based on the fact that a truly

# random distribution gets a mean value with a precision of this

order.

# The actual center is probably reachable with a few of these steps

move_amount = ((max_y - min_y) + (max_x - min_x)) / points.size

return [center, move_amount]

end

# Look for the best center by minimizing

# its maximum distance with the given points

# This is the same as minimizing its square

# (the maximum of the square distances)

# This can be seen as a dichotomy on a 2-dimensional space

def encircle(points)

center, ro = chose_original_center_and_move_amount(points)

# We’ll move with polar coordinates (distance + angle : ro/theta)

# How many directions do we follow?

# each new angle creates a new (costly) max_square_distance call

# so we take the minimum needed to move around and climb the

# square_distance gradient: 3

angle_count = 3

thetas = (1…angle_count).map { |c| (2 * c * Math::PI) / angle_count}

iter = 0

# Trace our moves

center_list = []

loop do

center_list << center

iter += 1

puts “#{iter}: #{center.max_square_distance}”

# Is there a best move available ?

if new_center = center.best_move(ro, thetas)

center = new_center

else

ro /= 2

# We stop when ro becomes too small to move the center

break unless center.can_be_moved_by(ro)

end

end

puts “Circle found in #{iter} iterations”

return [ Circle.new(center, sqrt(center.max_square_distance)),

center_list ]

end

points = Point.generate_samples(POINT_COUNT)

t = Time.now

circle, center_list = encircle(points)

puts circle

puts “Found in: #{Time.now - t}s”

# Generate a GnuPlot command file to see the result

f = File.open(‘plotcommand.cmd’, ‘w’)

f.print <<EOF

set parametric

set size square

set xrange [-0.2:1.2]

set yrange [-0.2:1.2]

set multiplot

plot “-”

#{points.map { |p| “#{p.x} #{p.y}” }.join("\n")}

end

plot “-” with lines

#{center_list.map { |p| “#{p.x} #{p.y}” }.join("\n")}

end

plot [0:2*pi]*

(sin(t)#{circle.radius})+#{circle.center.x},(cos(t)*#{circle.radius})+#{circle.center.y}

set nomultiplot

EOF

f.close

puts <<EOF

use ‘gnuplot -persist plotcommand.cmd’ to see

the circle, points and path to the center

EOF