[QUIZ] The Smallest Circle (#157)

Discussion in 'Ruby' started by Matthew D Moss, Feb 15, 2008.

  1. The three rules of Ruby Quiz 2:

    1. Please do not post any solutions or spoiler discussion for this
    quiz until 48 hours have passed from the time on this message.

    2. Support Ruby Quiz 2 by submitting ideas as often as you can! (A
    permanent, new website is in the works for Ruby Quiz 2. Until then,
    please visit the temporary website at <http://
    matthew.moss.googlepages.com/home>.

    3. Enjoy!

    Suggestion: A [QUIZ] in the subject of emails about the problem
    helps everyone
    on Ruby Talk follow the discussion. Please reply to the original
    quiz message,
    if you can.

    -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
    =-=-=-

    The Smallest Circle
    by Matthew Moss

    Your task this week sounds simple enough, but may be more difficult
    than it first appears. Given a set of points on a plane, your goal is
    to find the smallest circle that encloses those points.

    You are to provide a function, *encircle*, that takes an array of
    points and returns the smallest circle surrounding those points.
    Start with the following base code and extend as needed to solve the
    problem:

    class Point < Struct.new:)x, :y)
    def self.random
    Point.new(rand, rand)
    end

    def to_s
    "(#{x}, #{y})"
    end
    end

    class Circle < Struct.new:)center, :radius)
    def to_s
    "{center:#{center}, radius:#{radius}}"
    end
    end

    def encircle(points) # takes array of Point objects
    # returns a Circle object
    end


    I will be running several tests on the submitted solutions, with
    various point sets, to see how well they perform at this task. I
    recommend you you test your algorithm with a variety of sample sets,
    from small sets consisting of just 1-5 points, up to medium and
    larger sets, containing a few thousand points.

    To generate an array of random points, start with the above code and
    add:

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


    And then you may test your implementation like this:

    # encircle 10 random points
    puts encircle( generate_samples(10) )


    As mentioned, this one may be more difficult than it seems. Feel free
    to estimate the smallest circle, if you get stuck on getting the
    exact solution.
     
    Matthew D Moss, Feb 15, 2008
    #1
    1. Advertising

  2. [Note: parts of this message were removed to make it a legal post.]

    If one of the points is on the enclosing circle, is that a valid solution ?
     
    Horea Raducan, Feb 15, 2008
    #2
    1. Advertising

  3. Horea Raducan wrote:
    > If one of the points is on the enclosing circle, is that a valid solution ?
    >
    >

    I wouldn't care about such details :
    the points are defined by floating point coordinates : I'd say that
    defining a circle with a center and a radius in floating point which
    passes exactly through such points is most often impossible.
    There are obvious exceptions like circles with zero radius and some
    values where there's no approximation in the computations but these
    cases aren't likely to matter here.

    Lionel
     
    Lionel Bouton, Feb 15, 2008
    #3
  4. Matthew D Moss

    Matthew Moss Guest

    Re: The Smallest Circle (#157)

    On Feb 15, 1:06 pm, "Horea Raducan" <> wrote:
    > If one of the points is on the enclosing circle, is that a valid solution ?


    Yes, points on the circle are considered enclosed.

    Which also means that for an input set of one point, your function
    should return a circle with that point as its center and a zero
    radius.
     
    Matthew Moss, Feb 15, 2008
    #4
  5. Matthew D Moss

    Alex Shulgin Guest

    Re: The Smallest Circle (#157)

    On Feb 15, 3:19 pm, Matthew D Moss <> wrote:
    >
    > The Smallest Circle
    > by Matthew Moss
    >
    > Your task this week sounds simple enough, but may be more difficult
    > than it first appears. Given a set of points on a plane, your goal is
    > to find the smallest circle that encloses those points.
    >
    > You are to provide a function, *encircle*, that takes an array of
    > points and returns the smallest circle surrounding those points.
    > Start with the following base code and extend as needed to solve the
    > problem:
    >
    > class Point < Struct.new:)x, :y)
    > def self.random
    > Point.new(rand, rand)
    > end


    May I propose a bit refined random method?

    class Point < Struct.new:)x, :y)
    def self.random
    Point.new(0.5 - rand, 0.5 - rand)
    end

    This should give us negative coordinates as well. :)

    --
    Cheers and welcome new quizmasters!
    Alex
     
    Alex Shulgin, Feb 16, 2008
    #5
  6. Matthew D Moss

    John Joyce Guest

    Interesting, the numbering is continuing from Ruby Quiz ?
    Not starting at 1???

    SUGGESTION: I don't have time to do this quiz, but it looks like fun.
    I for one would encourage anybody to do it with some visual
    representation! (suggested libs: RMagick, Gosu, Rubygame, Ruby/SDL,
    etc...)
    An audio representation could be interesting as well (perhaps a sort
    of sonar?)
     
    John Joyce, Feb 16, 2008
    #6
  7. Matthew D Moss

    Todd Benson Guest

    On Feb 16, 2008 3:21 PM, John Joyce <> wrote:
    > Interesting, the numbering is continuing from Ruby Quiz ?
    > Not starting at 1???
    >
    > SUGGESTION: I don't have time to do this quiz, but it looks like fun.


    I also would love to do this one, but am on vacation :)

    I hope there's someone out there that supplies a multidimensional
    solution (something that gives the user a choice of the dimensional
    space; 3-D, 4-D, etc).

    > I for one would encourage anybody to do it with some visual
    > representation! (suggested libs: RMagick, Gosu, Rubygame, Ruby/SDL,
    > etc...)
    > An audio representation could be interesting as well (perhaps a sort
    > of sonar?)


    Those would be pretty interesting.

    Todd
     
    Todd Benson, Feb 16, 2008
    #7
  8. Matthew D Moss

    Matthew Moss Guest

    Re: The Smallest Circle (#157)

    > May I propose a bit refined random method?
    >
    > =A0 =A0 =A0class Point < Struct.new:)x, :y)
    > =A0 =A0 =A0 =A0 =A0def self.random
    > =A0 =A0 =A0 =A0 =A0 =A0 =A0Point.new(0.5 - rand, 0.5 - rand)
    > =A0 =A0 =A0 =A0 =A0end
    >
    > This should give us negative coordinates as well. :)



    Sure, but your algorithm should work with any set of random points,
    whether they fit in the unit square from 0 to 1, -0.5 to 0.5, or
    coordinates outside the unit square.
     
    Matthew Moss, Feb 17, 2008
    #8
  9. Matthew D Moss

    Bill Kelly Guest

    From: "John Joyce" <>
    >
    > SUGGESTION: I don't have time to do this quiz, but it looks like fun.
    > I for one would encourage anybody to do it with some visual
    > representation! (suggested libs: RMagick, Gosu, Rubygame, Ruby/SDL,
    > etc...)


    Here's an OpenGL/GLUT-based visualizer for the quiz:

    http://tastyspleen.net/~billk/ruby/quiz/157-smallest-circle/157_smallest_circle_visualization.rb

    The visualizer is designed to run with any quiz solution that
    provides the encircle() and generate_samples(n) methods at
    outer scope.


    A screenshot for those curious, but not curious enough to actually
    try installing ruby OpenGL bindings ... <grin>

    http://tastyspleen.net/~billk/ruby/quiz/157-smallest-circle/157_vis.png


    NOTE: I've only tried this program on ruby 1.8.4 so far... I have
    no idea if it might run on 1.9 ...

    Also... I believe I'm using yoshi's opengl-0.32f bindings. I haven't
    tried it with the very latest ruby/OpenGL bindings on rubyforge.


    Regards,

    Bill
     
    Bill Kelly, Feb 17, 2008
    #9
  10. Matthew D Moss

    Alex Shulgin Guest

    Re: The Smallest Circle (#157)

    On Feb 17, 8:59 am, Bill Kelly <> wrote:
    > From: "John Joyce" <>
    >
    >
    >
    > > SUGGESTION: I don't have time to do this quiz, but it looks like fun.
    > > I for one would encourage anybody to do it with some visual
    > > representation! (suggested libs: RMagick, Gosu, Rubygame, Ruby/SDL,
    > > etc...)

    >
    > Here's an OpenGL/GLUT-based visualizer for the quiz:
    >
    > http://tastyspleen.net/~billk/ruby/quiz/157-smallest-circle/157_small...


    I'm trying to utilize gnuplot for visualization. With any luck I can
    post the visual part without spoiling the quiz. ;)

    Anyone interested?

    --
    Alex
     
    Alex Shulgin, Feb 17, 2008
    #10
  11. Re: [QUIZ][SOLUTION] The Smallest Circle (#157)

    Hi,

    here's my solution of this nice quiz. It's heavily inspired by the last
    week's quiz, i.e. solving a non-linear equation. This time the program
    solves the non-smooth convex optimization problem

    min{max{|p_i-x|^2: p_i the given points}: x in R^k},

    where k=1,2,3,... (dimension), using a simple subgradient algorithm.
    Therefore, it works for an arbitrary number of points
    (up to 10000 in about 40s on my 5 year old notebook) and an
    arbitrary dimension (number of coordinates) of the points. The solutions
    are usually a good approximation of the optimal ones.


    ================================================================
    # extends point-class with an arbitrary dimension
    # and simple vector-operations
    class Point
    def initialize( *coords )
    @coords = coords.map{|x| x.to_f}
    end

    def size
    @coords.size
    end

    def []( index )
    @coords[index]
    end

    def []= ( index, x )
    @coords[index] = x
    end

    def self.random( n = 2 )
    Point.new( *(1..n).map{ 0.5 - rand } )
    end

    def +(p)
    return Point.new( *(0...size).map{|i| @coords + p} )
    end

    def -(p)
    return Point.new( *(0...size).map{|i| @coords - p} )
    end

    def *(a)
    return Point.new( *@coords.map{|x| x * a} )
    end

    def /(a)
    return Point.new( *@coords.map{|x| x / a} )
    end

    # calculate the inner product if this point with p
    def ip(p)
    (0...size).inject(0) { |sum, i| sum + @coords * p }
    end

    def to_s
    "(#{@coords.join(', ')})"
    end
    end

    class Circle < Struct.new:)center, :radius)
    def to_s
    "{center:#{center}, radius:#{radius}}"
    end
    end

    def generate_samples(n, dim = 2)
    (1..n).map { Point.random(dim) }
    end


    # calculate value and subgradient of
    # f(x,y) = max{(x-x_i)^2 + (y-y_i)^2: (x_i,y_i) in points}
    def evaluate( m, points )
    y_big = nil
    grad = nil
    points.each do |p|
    d = (m - p)
    y = d.ip( d ) # y = sum (m_i-p_i)^2
    if y_big.nil? or y > y_big
    y_big = y
    grad = d*2
    end
    end

    return [y_big, grad]
    end

    # perform simple subgradient algorithm
    # with given starting-point and number of iterations
    def encircle( points,
    x_start = Point.new( *([0]*points[0].size) ),
    max_iter = 100 )
    x = x_start
    y, g = nil, nil

    for k in 1..max_iter do
    y, g = evaluate( x, points )
    x = x - g/k
    end

    return Circle.new(x, Math.sqrt(y))
    end

    puts encircle( generate_samples(10000, 3) )
    ================================================================

    Bye,
    Frank
     
    Frank Fischer, Feb 17, 2008
    #11
  12. Douglas A. Seifert, Feb 17, 2008
    #12
  13. [Note: parts of this message were removed to make it a legal post.]

    Hello,

    Had a lot of fun with this one :) My solution works by taking the average of
    all points as the center, and the farthest outlying point as the radius. An
    optimal solution is then found by moving the circle center towards the
    outlier as the radius is reduced, while all points still are within the
    radius. The same principles could be applied for higher dimensions as well.

    class Point < Struct.new:)x, :y)
    def self.random
    Point.new(rand, rand)
    end

    def to_s
    "(#{x}, #{y})"
    end
    end

    class Circle < Struct.new:)center, :radius)
    def to_s
    "{center:#{center}, radius:#{radius}}"
    end
    end

    # Calculate distance between points
    def distance(pt_a, pt_b)
    Math.sqrt((pt_a.x - pt_b.x) * (pt_a.x - pt_b.x) +
    (pt_a.y - pt_b.y) * (pt_a.y - pt_b.y))
    end

    # Determine if given points are all within a circle
    def inside_circle?(points, circle)
    for point in points
    dist = distance(point, circle.center)
    return false if dist > circle.radius
    end
    true
    end

    def encircle(points) # takes array of Point objects
    # find the average midpoint of all points
    mid = Point.new(
    points.inject(0){|sum, pt| sum += pt.x} / (points.size * 1.0),
    points.inject(0){|sum, pt| sum += pt.y} / (points.size * 1.0))

    # sort points by longest distance from midpoint
    # longest point (index 0) is the initial radius
    points.sort!{|a,b| distance(mid, a) <=> distance(mid, b) }.reverse!

    # Taking the average midpoint does not necessarily work because the points
    may
    # be weighted more heavily to one side. We correct for this by sliding the
    circle
    # along the line from its original center to the outmost point, decreasing
    the
    # radius as we go. Keep doing this until the circle can be made no
    smaller.
    point = points[0]
    slope = (mid.y - point.y) / (mid.x - point.x)
    new_pt, delta, sign = mid, 0.0, 1.0
    sign = -1.0 if mid.x > point.x
    while inside_circle?(points, Circle.new(new_pt, distance(new_pt, point)))
    mid = new_pt
    delta += 0.000001 * sign
    new_pt = Point.new(mid.x + delta, mid.y + (slope * (delta)))
    end

    return Circle.new(mid, distance(mid, point))
    end

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


    Here are Pasties of the program and a RMagick script to generate images:

    Algorithm: http://pastie.caboo.se/153387
    Graphics code: http://pastie.caboo.se/153388
    Images:
    http://justin.ethier.googlepages.com/157_jae_25_2.jpg
    http://justin.ethier.googlepages.com/157_jae_100_0.jpg
    http://justin.ethier.googlepages.com/157_jae_1000_3.jpg

    Thanks,

    Justin

    On Feb 15, 2008 8:19 AM, Matthew D Moss <> wrote:

    > The three rules of Ruby Quiz 2:
    >
    > 1. Please do not post any solutions or spoiler discussion for this
    > quiz until 48 hours have passed from the time on this message.
    >
    > 2. Support Ruby Quiz 2 by submitting ideas as often as you can! (A
    > permanent, new website is in the works for Ruby Quiz 2. Until then,
    > please visit the temporary website at <http://
    > matthew.moss.googlepages.com/home>.
    >
    > 3. Enjoy!
    >
    > Suggestion: A [QUIZ] in the subject of emails about the problem
    > helps everyone
    > on Ruby Talk follow the discussion. Please reply to the original
    > quiz message,
    > if you can.
    >
    > -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
    > =-=-=-
    >
    > The Smallest Circle
    > by Matthew Moss
    >
    > Your task this week sounds simple enough, but may be more difficult
    > than it first appears. Given a set of points on a plane, your goal is
    > to find the smallest circle that encloses those points.
    >
    > You are to provide a function, *encircle*, that takes an array of
    > points and returns the smallest circle surrounding those points.
    > Start with the following base code and extend as needed to solve the
    > problem:
    >
    > class Point < Struct.new:)x, :y)
    > def self.random
    > Point.new(rand, rand)
    > end
    >
    > def to_s
    > "(#{x}, #{y})"
    > end
    > end
    >
    > class Circle < Struct.new:)center, :radius)
    > def to_s
    > "{center:#{center}, radius:#{radius}}"
    > end
    > end
    >
    > def encircle(points) # takes array of Point objects
    > # returns a Circle object
    > end
    >
    >
    > I will be running several tests on the submitted solutions, with
    > various point sets, to see how well they perform at this task. I
    > recommend you you test your algorithm with a variety of sample sets,
    > from small sets consisting of just 1-5 points, up to medium and
    > larger sets, containing a few thousand points.
    >
    > To generate an array of random points, start with the above code and
    > add:
    >
    > def generate_samples(n)
    > (1..n).map { Point.random }
    > end
    >
    >
    > And then you may test your implementation like this:
    >
    > # encircle 10 random points
    > puts encircle( generate_samples(10) )
    >
    >
    > As mentioned, this one may be more difficult than it seems. Feel free
    > to estimate the smallest circle, if you get stuck on getting the
    > exact solution.
    >
    >
     
    Justin Ethier, Feb 17, 2008
    #13
  14. On 2008-02-17, Justin Ethier <> wrote:
    > Hello,
    >
    > Had a lot of fun with this one :) My solution works by taking the average of
    > all points as the center, and the farthest outlying point as the radius. An
    > optimal solution is then found by moving the circle center towards the
    > outlier as the radius is reduced, while all points still are within the
    > radius. The same principles could be applied for higher dimensions as well.


    Nice solution, this was also one of my first ideas. The problem with
    this approach is if many points are close to each other and only few
    points are far away. In this case the center is near that mass of points
    and the line-search may fail to find an optimal solution. For example:

    puts encircle( ([[0.0,0.0]] * 100 + [[0.0,1.0], [1.0,0.0]]).map{|p|
    Point.new(*p)})

    (100 points at (0,0), one at (1,0) and one at (0,1)) gives a solution of
    {center:(0.00980392156862745, 0.00980392156862745), radius:0.990244611507173}
    but optimal would be (1/2, 1/2) and radius 1/sqrt(2).

    Frank
     
    Frank Fischer, Feb 17, 2008
    #14
  15. Matthew D Moss

    ThoML Guest

    Re: The Smallest Circle (#157)

    > Had a lot of fun with this one :) My solution works by taking the average of
    > all points as the center, and the farthest outlying point as the radius. An
    > optimal solution is then found by moving the circle center towards the
    > outlier as the radius is reduced


    Well, I didn't have the time to check out how to do this the right way
    (as Douglas did it). I nevertheless tried several approaches. This one
    has some imprecision built-in but it turns out (in comparison with
    Douglas' solution) that for random point sets, it performs more or
    less well (for an incorrect semi-solution).

    I also attach some code to run several solution against each other,
    which is fun if you exclude the precise solutions. This requires the
    script lib from http://redshift.sourceforge.net/script/.

    Cool quiz BTW that made me refresh my knowledge on triangles etc. I
    wish a had had more time so that I could also read up on appropriate
    algorithms.

    Regards,
    Thomas.


    My semi-solution that uses a similar approach as Justin:


    def encircle(points) # takes array of Point objects
    points = points.uniq
    return if points.nil? or points.empty?

    x = points.inject(0.0) {|a, p| a += p.x} / points.size
    y = points.inject(0.0) {|a, p| a += p.y} / points.size
    a = Point.new(x, y)
    return Circle.new(a, 0) unless points.size > 1

    f = points.sort_by {|p| distance(a, p)}[-1]
    df = distance(f, a)
    b = med(f, a)
    e = 1e-12
    1000.times do
    db = distance(f, b)
    if points.all? {|p| distance(b, p) <= db + e}
    da = distance(f, a)
    if (da - db).abs < e
    return Circle.new(b, db)
    else
    a, b = b, med(f, b)
    end
    else
    b = med(a, b)
    end
    end
    raise RuntimeError
    end

    def med(a, b)
    Point.new((b.x - a.x) / 2 + a.x, (b.y - a.y) / 2 + a.y)
    end

    def distance(p1, p2)
    Math.sqrt((p1.x - p2.x) ** 2 + (p1.y - p2.y) ** 2)
    end



    The code to run the tests:

    #!/usr/bin/env ruby19

    require 'script'

    files = [
    's_tml.rb',
    's_frank.rb',
    's_justin.rb',
    's_douglas.rb',
    ]


    def generate_samples(n)
    (1..n).map { [rand - 0.5, rand - 0.5] }
    end

    # Numer of samples
    # N = 10
    N = 50
    # N = 100
    # N = 200

    # Tolerance
    E = 1e-10

    # Your favourite solution, you might want to change this
    F = 's_tml.rb'

    sets = Array.new(N) {|i| generate_samples(i + 2)}
    solutions = Hash.new {|h, k| h[k] = {}}

    # Best choice per sample
    winners = Hash.new {|h,k| h[k] = 0}

    # Number of invalid solutions
    losers = Hash.new {|h,k| h[k] = 0}


    puts
    files.each do |f|
    print f
    script = Script.load(f)
    sets.each_with_index do |ps, i|
    if i % 10 == 0
    print '.'; STDOUT.flush
    end
    ps = ps.map {|p| script::point.new(*p)}
    v = script.encircle(ps)
    if ps.any? {|p| distance(v.center, p) > v.radius + E}
    losers[f] += 1
    print 'x'; STDOUT.flush
    v.radius += 1000.0
    end
    solutions[f] = v
    end
    puts
    end
    puts

    f_dist = 0

    puts " %s" % files.map {|f| '%10s' % f}.join(' ')
    solutions.each do |i, s|
    winner = files.sort_by {|f| s[f] ? s[f].radius : 100000000}[0]
    winners[winner] += 1
    f_d = (s[F].radius - s[winner].radius).abs
    f_dist = f_d if f_d > f_dist
    puts '%5d %10s %s (%6.4f)' % [
    i,
    '(%s)' % winner,
    files.map {|f| s[f] && '%10.4f' % s[f].radius}.join(' '),
    f_d,
    ]
    end
    puts


    puts "Winners (best circle):"
    winners.each do |f, n|
    puts '%10s %d' % [f, n]
    end
    puts

    puts "Losers (invalid solutions):"
    losers.each do |f, n|
    puts '%10s %d' % [f, n]
    end
    puts

    puts "Favourite vs Winner:"
    puts 'Max distance: %6.4f' % f_dist
    puts
     
    ThoML, Feb 17, 2008
    #15
  16. Matthew D Moss

    Bill Kelly Guest

    [QUIZ][SOLUTION] The Smallest Circle (#157)

    Greetings,

    Not sure if my solution is worth posting.... I took a naive
    axis-aligned bounding box approach.

    http://tastyspleen.net/~billk/ruby/quiz/157-smallest-circle/157_naive_aabb.rb

    I knew it would be inaccurate, but I was curious how bad it
    would be, so I wrote an OpenGL/GLUT based visualizer:

    http://tastyspleen.net/~billk/ruby/quiz/157-smallest-circle/157_smallest_circle_visualization.rb

    (The visualizer will work with any quiz solution that
    provides the encircle() and generate_samples(n) methods,
    and Point#x and Point#y.)


    Occasionally the naive approach lucks out with a dataset
    that allows it to look good:

    http://tastyspleen.net/~billk/ruby/quiz/157-smallest-circle/157_vis.png

    But more often it looks something like this:

    http://tastyspleen.net/~billk/ruby/quiz/157-smallest-circle/157_vis_naive.png


    I had planned to improve my solution by taking an approach
    similar to Justin Ethier's. . . . By iterating and "moving the
    circle center towards the outlier as the radius is reduced,
    while all points still are within the radius."

    I may still code that if I have time... but then again, it's
    already been done. :)


    Regards,

    Bill
     
    Bill Kelly, Feb 17, 2008
    #16
  17. Matthew D Moss

    ThoML Guest

    Re: The Smallest Circle (#157)

    And again I reply to my own post. The test code requires a distance()
    method:

    def distance(p1, p2)
    Math.sqrt((p1.x - p2.x) ** 2 + (p1.y - p2.y) ** 2)
    end

    In case somebody cares.
     
    ThoML, Feb 17, 2008
    #17
  18. Lionel Bouton, Feb 17, 2008
    #18
  19. 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
     
    Lionel Bouton, Feb 17, 2008
    #19
  20. Matthew D Moss

    Bill Kelly Guest

    Re: [QUIZ][SOLUTION] The Smallest Circle (#157)

    From: "Bill Kelly" <>
    >
    > I had planned to improve my solution by taking an approach
    > similar to Justin Ethier's. . . . By iterating and "moving the
    > circle center towards the outlier as the radius is reduced,
    > while all points still are within the radius."
    >
    > I may still code that if I have time... but then again, it's
    > already been done. :)


    Ah, what the heck... :) Here 'tis..

    http://tastyspleen.net/~billk/ruby/quiz/157-smallest-circle/157_smallest_circle.rb

    # Algorithm used: First enclose the points in a simple
    # axis-aligned bounding box, and choose the center of the
    # box as the initial center of the circle. Then, binary
    # search to determine how far we can move the center point
    # toward the furthest outlying point, while keeping the
    # rest of the points still inside the circle.


    Regards,

    Bill
     
    Bill Kelly, Feb 17, 2008
    #20
    1. Advertising

Want to reply to this thread or ask your own question?

It takes just 2 minutes to sign up (and it's free!). Just click the sign up button to choose a username and then you can ask your own questions on the forum.
Similar Threads
  1. Brent Turner
    Replies:
    0
    Views:
    302
    Brent Turner
    Sep 29, 2003
  2. Reinhard Koenig
    Replies:
    5
    Views:
    531
    Reinhard Koenig
    Nov 17, 2003
  3. Lionel Bouton
    Replies:
    16
    Views:
    197
    M. Edward (Ed) Borasky
    Feb 23, 2008
  4. Matthew Moss

    [QUIZ] Circle Drawing (#166)

    Matthew Moss, Jun 13, 2008, in forum: Ruby
    Replies:
    18
    Views:
    298
    ThoML
    Jul 12, 2008
  5. Daniel Moore
    Replies:
    19
    Views:
    338
Loading...

Share This Page