--WIyZ46R2i8wDzkSu
Content-Type: text/plain; charset=us-ascii
Content-Disposition: inline

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


--WIyZ46R2i8wDzkSu
Content-Type: text/plain; charset=us-ascii
Content-Disposition: attachment; filename="smallest_enclosing_circle2.rb"

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
 
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

# ----- 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
  #draw_circle(circle0.center, circle0.radius, :green) if VISUALIZE
  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
  #draw_circle(circle1.center, circle1.radius, :blue) if VISUALIZE
  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
  #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=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

def draw_circle(center, radius, color=false)
  $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=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) }
    draw_circle(circle.center, circle.radius, :orange)
    canvas = Magick::ImageList.new
    canvas.new_image($size, $size)
    $gc.draw(canvas)
    canvas.display
  end
  
end

--WIyZ46R2i8wDzkSu--