From: Lionel Bouton <lionel-subscription@...>
Date: 2008-02-18T03:40:27+09:00
Subject: Re: [QUIZ] The Smallest Circle (#157)

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