Smallest circle covering set of points
in [General Programming]
|
Prev: Heathfield: You're wrong: here's proof: (WAS: Re: Welcome To The Word Wide Mirror Sci.Math)
Next: SETP-10 Call for papers
From: Torben =?iso-8859-1?Q?=C6gidius?= Mogensen on 18 Nov 2009 08:04 Chip Eastham <hardmath(a)gmail.com> writes: > Evidently the Welzl algorithm generalizes to > higher dimensions more readily than Meggido's. And also more easily to other kinds of enclosing figures, such as pairs of parallel lines/planes/hyperplanes. It is also a lot simpler: You can describe Welzl's algorithm in just a few lines of pseudocode: minCircle(ps) = mC(ps, {}) mC(ps, qs) = the circle defined by (ps ++ qs) if |ps|+|qs| <= 3 mC({p} ++ ps, qs) = let c = mC(ps, qs) in if p is inside c then c else mC(ps, {p} ++ qs) ++ is disjount union of sets. If you have three points or less, it is easy to find the minimum circle: - One point gives a zero-radius circle centered at this point. - Two points define a diameter of a circle. - Three points define a unique circle that touches all three, but a circle defined by two of these (as above) might be smaller. This is quickly determined. Torben
From: Ronald Bruck on 20 Nov 2009 12:06 In article <181120090046541409%bruck(a)math.usc.edu>, Ronald Bruck <bruck(a)math.usc.edu> wrote: > In article > <5b980a85-0f88-440e-aabc-a0faf1345ca9(a)w19g2000yqk.googlegroups.com>, > Chip Eastham <hardmath(a)gmail.com> wrote: > > > On Nov 14, 11:54�pm, Andrew Tomazos <and...(a)tomazos.com> wrote: > > > On Nov 15, 5:29�am, Andrew Tomazos <and...(a)tomazos.com> wrote: > > > > > > > Given a set of points (a1,b1) (a2,b2)..(an,bn) in the plane - what's a > > > > good algorithm for determining the centerpoint and radius of the > > > > smallest circle that will cover all the points? �Thanks, Andrew. > Hmmm. I think I've found a way to solve this using semi-definite > programming. > > Call the known points (u_i, vi) with known radius ri. Let (x,y) be the > center of the spanning circle and r its radius. The problem is then to > minimize r subject to the conditions > > (*) (ui - x)^2 + (vi - y)^2 <= (r-ri)^2, r >= ri. > > This is clearly a quadratic programming problem loosely of SDP type. > We can make it exactly into an SDP problem with a little manipulation. > > We consider a family of two-dimensional vectors as the unknowns. (There > are obvious generalizations to balls in R^N.) These include the > vectors e1 and e2 --- the standard orthonormal basis --- which we > regard as unknowns by the fiction that they are unknown but satisfy the > relations > > (1) <ei, ej> = Kronecker-delta (i,j). > > The only other unknowns are the (really!) unknown X = (x,y) and R = > (r,0). These satisfy the equations > > (2) <R,e1> = r, <R,e2> = 0. > > Then the constraint (*) becomes > > (**) ui^2 + vi^2 - ri^2 <= r^2 - 2 r ri + 2 ui x + 2 vi y - (x^2 + y^2) > > = <R,R> - 2ri <R,e1> + 2 ui <X,e1> + 2 vi <X,e2> - <X,X> > > The problem is to minimize r = <R,e1> subject to (**) and the > additional constraints that > > <X,e1> >= ri for all i. > > and that's why it's an SDP problem: the Gramian of the vectors e1, e2, > R and X is positive semi-definite, and the problem is one of minimizing > one element of the Gramian (<R,e1>) subject to linear constraints on > the elements of the Gramian. > > I'm currently finishing a paper where I did a similar thing for > intersecting spheres (closed balls) in R^N. It turns out to be rather > efficient. THAT application even allowed us to find proximinal points > on intersections of spheres, COMPLEMENTS of spheres, and closed > half-spaces (a problem which is not generally convex). > > One catch with SDP is that you can't really specify the dimension in > which the problem lives. In particular, you would have to check > whether > > <X,X> = <X,e1>^2 + <X,e2>^2. > > I suspect you'll find solutions where this isn't true, that the SDP > lives in dimension 3 or 4; but that with a little manipulation you can > bring it down to dimensions 2. (You can't just toss that last equation > in as a constraint, because it's not a LINEAR constraint.) > > I'll try programming a few problems and see what happens. I should mention that I'm solving a slightly more general problem: find the smallest circle which contains as subsets a finite number of CIRCLES (not just points). I'm not convinced this is very different from the original problem. -- Ron Bruck
From: Chip Eastham on 21 Nov 2009 11:17 On Nov 18, 3:46 am, Ronald Bruck <br...(a)math.usc.edu> wrote: > In article > <5b980a85-0f88-440e-aabc-a0faf1345...(a)w19g2000yqk.googlegroups.com>, > > > > Chip Eastham <hardm...(a)gmail.com> wrote: > > On Nov 14, 11:54 pm, Andrew Tomazos <and...(a)tomazos.com> wrote: > > > On Nov 15, 5:29 am, Andrew Tomazos <and...(a)tomazos.com> wrote: > > > > > Given a set of points (a1,b1) (a2,b2)..(an,bn) in the plane - what's a > > > > good algorithm for determining the centerpoint and radius of the > > > > smallest circle that will cover all the points? Thanks, Andrew. > > > > (some thoughts) > > > > If the centerpoint is (x,y), > > > > than the distance to each point k is given by sqrt((x-ak)^2 + (y-bk) > > > ^2) > > > > so we want to select (x,y) such as to minimize the largest element of > > > the set: > > > > { (x-a1)^2 + (y-b1)^2 , > > > (x-a2)^2 + (y-b1)^2 , > > > ... , > > > (x-an)^2 + (y-bn)^2 > > > } > > > > But how to do that? > > > -Andrew. > > > The short answer is to check out Dave Eberly's > > implementation linked below in the thread! > > > From what I remember there is a "naive" strategy > > that involves pivoting. First find the pair of > > points farthest apart; the circle must have at > > least this diameter. If that diameter (between > > two points farthest apart) yields a containing > > circle, we're done. If not, the minimum circle > > will have (at least) three points on its > > circumference. Pick a point outside the given > > "two point" circle, and see if the circle one > > that and the first two points contains all the > > points. If so, done. > > > Then you have a sequence of three-point circles > > to consider until you find one that contains > > all the points. Here's where "pivoting" comes > > into the picture. We have three points A,B,C > > that determine my current circle, and a point > > D not contained by the circle. Which point > > A,B,C is to be replaced by D? A simple way > > to answer this is by trying all three and > > comparing which has the smallest radius yet > > still contains the replaced point. IIRC a > > slightly more efficient check is based on the > > angles of the resulting triangle. > > > The inscribed triangle (of three points on > > the circle) must be acute, else the circle's > > radius isn't minimal (having already ruled > > out the two-point circle using two farthest > > apart). The right choice for D to replace > > is that point which results in the triangle > > with the minimal maximum angle (whose largest > > angle is the smallest among 3 possibilities). > > Since the largest angle opposes the longest > > side, this check can be done via the cosine > > formula for triangles. > > Hmmm. I think I've found a way to solve this using semi-definite > programming. > > Call the known points (u_i, vi) with known radius ri. Let (x,y) be the > center of the spanning circle and r its radius. The problem is then to > minimize r subject to the conditions > > (*) (ui - x)^2 + (vi - y)^2 <= (r-ri)^2, r >= ri. > > This is clearly a quadratic programming problem loosely of SDP type. > We can make it exactly into an SDP problem with a little manipulation. > > We consider a family of two-dimensional vectors as the unknowns. (There > are obvious generalizations to balls in R^N.) These include the > vectors e1 and e2 --- the standard orthonormal basis --- which we > regard as unknowns by the fiction that they are unknown but satisfy the > relations > > (1) <ei, ej> = Kronecker-delta (i,j). > > The only other unknowns are the (really!) unknown X = (x,y) and R = > (r,0). These satisfy the equations > > (2) <R,e1> = r, <R,e2> = 0. > > Then the constraint (*) becomes > > (**) ui^2 + vi^2 - ri^2 <= r^2 - 2 r ri + 2 ui x + 2 vi y - (x^2 + y^2) > > = <R,R> - 2ri <R,e1> + 2 ui <X,e1> + 2 vi <X,e2> - <X,X> > > The problem is to minimize r = <R,e1> subject to (**) and the > additional constraints that > > <X,e1> >= ri for all i. > > and that's why it's an SDP problem: the Gramian of the vectors e1, e2, > R and X is positive semi-definite, and the problem is one of minimizing > one element of the Gramian (<R,e1>) subject to linear constraints on > the elements of the Gramian. > > I'm currently finishing a paper where I did a similar thing for > intersecting spheres (closed balls) in R^N. It turns out to be rather > efficient. THAT application even allowed us to find proximinal points > on intersections of spheres, COMPLEMENTS of spheres, and closed > half-spaces (a problem which is not generally convex). > > One catch with SDP is that you can't really specify the dimension in > which the problem lives. In particular, you would have to check > whether > > <X,X> = <X,e1>^2 + <X,e2>^2. > > I suspect you'll find solutions where this isn't true, that the SDP > lives in dimension 3 or 4; but that with a little manipulation you can > bring it down to dimensions 2. (You can't just toss that last equation > in as a constraint, because it's not a LINEAR constraint.) > > I'll try programming a few problems and see what happens. Hi, Ron: This is similar to what I read about Meggido's approach, in that his linear programming for the smallest circle gives a problem in 3D, and IIRC for minimum spheres a problem in 4D. regards, chip
From: Dave Eberly on 22 Nov 2009 13:43 "Ronald Bruck" <bruck(a)math.usc.edu> wrote in message news:201120090906363755%bruck(a)math.usc.edu... > I should mention that I'm solving a slightly more general problem: > find the smallest circle which contains as subsets a finite number of > CIRCLES (not just points). > > I'm not convinced this is very different from the original problem. See publication lists at: http://www.inf.ethz.ch/personal/fischerk/ A couple of implementations at: http://miniball.sourceforge.net/ http://www.compgeom.com/meb/ -- Dave Eberly http://www.geometrictools.com
From: keith on 23 Nov 2009 20:48
Andrew Tomazos pravi: > On Nov 15, 5:29 am, Andrew Tomazos <and...(a)tomazos.com> wrote: >> Given a set of points (a1,b1) (a2,b2)..(an,bn) in the plane - what's a >> good algorithm for determining the centerpoint and radius of the >> smallest circle that will cover all the points? Thanks, Andrew. > > (some thoughts) > > If the centerpoint is (x,y), > > than the distance to each point k is given by sqrt((x-ak)^2 + (y-bk) > ^2) > > so we want to select (x,y) such as to minimize the largest element of > the set: > > { (x-a1)^2 + (y-b1)^2 , > (x-a2)^2 + (y-b1)^2 , > ... , > (x-an)^2 + (y-bn)^2 > } > > But how to do that? It's very easy actually. I just wanted to ask a question related to yours. You have a point set P with elements rp, you're looking for a center rc of a minimum circle of radius r. What I do in my algorithm is, I first derive an AABB of the point set and then put a starting rc in the center of this AABB. That's my initial guess. Then I use the hooke-jeeves algorithm (paper is online, but better check Schaum's Operations Research for an easy explanation) to minimize: max({||rp-rc||^2}) with rc as the argument and all rps from P taken into account. You'll get close to the minimum ball center this way and will also get the radius (returned by the cost function). Even so, the welzl algorithm will usually outdo this approach, due finite floating point precision. What I wanted to ask is: Is the cost function max({||rp-rc||^2}) over all rp continuous? Differentiable? Maybe one could get to the rc even faster using derivatives? This approach works, I have working code. |