Re: Convex Hull
- To: mathgroup at smc.vnet.net
- Subject: [mg20947] Re: [mg20839] Convex Hull
- From: Daniel Lichtblau <danl at wolfram.com>
- Date: Wed, 1 Dec 1999 01:50:14 -0500 (EST)
- References: <199911170840.DAA02595@smc.vnet.net>
- Sender: owner-wri-mathgroup at wolfram.com
Virgil Stokes wrote:
>
> My problem is this:
>
> I have two ellipses that in general are different (lengths of minor/major
> axes may differ) and are positioned at random in the XY-plane
> with the condition that they do not overlap or touch at any point along their
> boundaries.
>
> The problem is to find the convex hull for these two ellipses -- fast.
> Any ideas/suggestions would be appreciated
>
> Please reply direct to my email since I am still unable to get messages
> generated
> from the MathGroup list.
>
> -- Virgil
Here is a reasonable method. It works when the ellipses are presented in
implicit form. I checked via e-mail that specifically they will have the
form
(x-a)^2/e^2 + (y-b)^2/f^2 == 1
The code given will actually work more generally with pairs of implicit
curves. I'll set this up using indeterminate coefficients, then fill
them in with random values. The first thing to realize is that the
convex hull will be comprised of two arcs (one from each ellipse) and
two line segments that each intersect both ellipses tangentially. To
find these we'll need the pairs of points. These we get just by solving
a system of polynomial equations.
ellipse1 = (x-a1)^2/p1^2 + (y-b1)^2/q1^2 - 1;
ellipse2 = (x-a2)^2/p2^2 + (y-b2)^2/q2^2 - 1;
slope1 = First[yprime /. Solve[
(D[ellipse1/.y->y[x],x] /. {y[x]->y, y'[x]->yprime}) == 0, yprime]];
slope2 = First[yprime /. Solve[
(D[ellipse2/.y->y[x],x] /. {y[x]->y, y'[x]->yprime}) == 0, yprime]];
rules1 = {x->x1, y->y1};
rules2 = {x->x2, y->y2};
polys = Numerator[Together[{ellipse1/.rules1, ellipse2/.rules2,
(slope1/.rules1)-(slope2/.rules2),
(y2-y1)-(slope1/.rules1)*(x2-x1)}]];
These polynomials are exactly what must be satisfied in order for the
pair of points {{x1,y1}, {x2,y2}} to satisfy:
(i) {xj,yj} lies on ellipsej.
(ii) The tangents to ellipsej through {xj,yj} have the same slopes.
(iii) {x2,y2} lies on the tangent to ellipse1 through {x1,y1}.
Now we'll make random substitutions to get two actual ellipses.
SeedRandom[1111]
replacements = Join[
Thread[{a1,b1}->Table[Random[Real,{2,5}],{2}]],
Thread[{a2,b2}->Table[Random[Real,{-5,-2}],{2}]],
Thread[{p1,q2}->Table[Random[Real,{.3,.9}],{2}]],
Thread[{p2,q1}->Table[Random[Real,{1.5,3}],{2}]]];
{ell1, ell2} = {ellipse1,ellipse2} /. replacements;
Next we find all such pairs of points. For two quadratic plane curves we
get four solutions. In some cases two might be complex valued. In
particular this will happen if the ellipses intersect. While it is
stipulated that this situation does not occur, we'll explicitly take
only real-valued solutions to make this code fairly general.
In[145]:= pairs = Cases[
{{x1,y1},{x2,y2}} /. NSolve[polys/.replacements, {x1,y1,x2,y2}],
{{_Real,_Real}, {_Real,_Real}}]
Out[145]= {{{4.49236, 6.30875}, {-5.76366, -4.49148}},
{{4.24555, 5.8998}, {-0.0761669, -4.60452}},
{{5.18463, 2.67326}, {-0.0907271, -4.63022}},
{{4.99884, 2.48789}, {-5.70571, -4.44488}}}
Now we need to determine which two pairs give the convex hull segments
and which do not. One way to find this is as follows: the midpoints of
the other three segments must all lie on the same side of the line
through a given hull segment. The code below will use this fact to find
the correct two segments.
midpoint[{p1_,p2_}] := (p1+p2)/2
slope[{p1_,p2_}] := (p2[[2]]-p1[[2]]) / (p2[[1]]-p1[[1]])
sameside[pair:{p1_,p2_}, midpts:{mid1_,mid2_,mid3_}] := Module[
{x, y, m=slope[pair], line, signs},
line = (y-p1[[2]]) - m*(x-p1[[1]]);
Apply[Equal, Map[Sign,
line /. Map[Thread[{x,y}->#]&, midpts]]]
]
getHullPoints[e1_, e2_, pairs:{pair1_, pair2_, pair3_, pair4_}] :=
Module[
{midpts, pointdata},
pointdata = MapIndexed[{#1, midpoint[#1],#2}&, pairs];
midpts = Map[#[[2]]&, pointdata];
Map[First, Select[pointdata, sameside[#[[1]],
Drop[midpts,#[[3]]]]&]]
]
hullpts = getHullPoints[ell1, ell2, pairs]
One can readily put this together in a picture.
plot1 = ContourPlot[ell1, {x,-8,8}, {y,-8,8}, Contours->{0},
PlotPoints->200, ContourShading->False,
DisplayFunction->Identity];
plot2 = ContourPlot[ell2, {x,-8,8}, {y,-8,8}, Contours->{0},
PlotPoints->200, ContourShading->False,
DisplayFunction->Identity];
lines = Map[Graphics[Line[#]]&, hullpts];
Show[plot1, plot2, Apply[Sequence,lines],
DisplayFunction->$DisplayFunction];
A final step might be to convert the ellipses to parametric form and
only use the two arcs that lie on the hull. These in turn can be deduced
either by checking angles with the tangent segments or by a slight
variation of the midpoints-on-one-side technique.
Daniel Lichtblau
Wolfram Research