From: "S. Grudszus" Subject: Q: general ellipse intersection formula Date: Wed, 08 Sep 1999 06:55:06 GMT Newsgroups: sci.image.processing,comp.ai.vison,sci.math.research To: sci-math-research@moderators.isc.org Keywords: computing intersection of two ellipses Hello, perhaps this seems to be a simple question to some readers, but I ran into some trouble solving a formula which should calculate the points of intersection of two ellipses, with different centres and angles. I use following ellipse parameters and equations: ellipse 1: curve parameter t1, angle phi1, centre mx1/my1, axis rx1/ry1 ellipse 2: curve parameter t2, angle phi2, centre mx2/my2, axis rx2/ry2 ellipse 1: x1(t1)= mx1+ rx1* cos(phi1)* cos(t1)+ ry1* sin(phi1)* sin(t1) y1(t1)= my1+ ry1* cos(phi1)* sin(t1)- rx1* sin(phi1)* cos(t1) ellipse 2: x2(t2)= mx2+ rx2* cos(phi2)* cos(t2)+ ry2* sin(phi2)* sin(t2) y2(t2)= my2+ ry2* cos(phi2)* sin(t2)- rx2* sin(phi2)* cos(t2) These formula are not easy to solve, but after repeated usage of some trigonometric formula from BRONSTEIN, especially >> FORMULA 1: >> a1 * sin(t+alpha1) + a2 * sin(t+alpha2) = A * sin(t+alpha) >> >> with: >> A= SQRT(a1*a1 + a2*a2 + 2*a1*a2*cos(alpha2-alpha1)) >> tan(alpha)= (a1*sin(alpha1)+a2*sin(alpha2)) / >> (a1*cos(alpha1)+a2*cos(alpha2)) >> FORMULA 2: >> a * sin(t) + b * cos(t) = A * sin(t+phi) >> >> with: >> A= SQRT(a1*a1 + a2*a2) >> tan(phi)= b / a I could find an equation which solves for the one unknown curve parameter t2 and use this result to solve for unknown t2. Of course, while calculating on and on, there are some places where partly solutions have to be fulfilled. This is what I expected, as there are different cases, which have to be discriminated: a) no points of intersections (poi) b) exactly one poi, where the ellipses touch themselves c) two poi, where they intersect d) three poi, two times intersect, one time touch e) four poi, where they intersect f) trivial solution, when ellipses are equal I would be very lucky, to find out a general solution to this problem. The best case would be two find a unique solution for parameters t1 and t2, describing the cases b, c, d, e, I guess. I surmise, cases a, b, f are easier to calculate directly form ellipse parameters, i.e. there distance, etc. Any help appreciated. Who knows similar solutions to this problem?? Is there any mathematical algorithm to this problem?? Please answer to: S.Grudszus@plettac-electronics.de Regards, Stefan Sent via Deja.com http://www.deja.com/ Share what you know. Learn what you don't. ============================================================================== From: Dave Rusin Subject: Re: Q: general ellipse intersection formula Date: Wed, 8 Sep 1999 03:32:58 -0500 (CDT) Newsgroups: [missing] To: s.grudszus@plettac-electronics.de [deletia --djr] From the data given you can write each ellipse as the solution set to a quadratic polynomial in x and y. You can eliminate one of these two from the two quadratic equations, to get a single polynomial of degree 4 to solve. Solve that quartic and plug in to get the other variable; each real root of the quartic gives one of the points of intersection. Any symbolic algebra program will do elimination. Here is the output from one of them. If the ellipses are described by the euqations f1:=a*x^2+b*x*y+c*y^2+d*x+e*y+f; f2:=a1*x^2+b1*x*y+c1*y^2+d1*x+e1*y+f1; then we may for example eliminate x from these equations and discover y must be the root of the quartic y^4*a^2*c1^2 - y^4*a*b*b1*c1 - 2*y^4*a*c*a1*c1 + y^4*a*c*b1^2 + [much deleted -- see below. --djr] + 2*y*e*f*a1^2 + a*f*d1^2 - d*f*a1*d1 + f^2*a1^2 Solve this however you like (usually people use numerical methods like Newton's method, but in principle you could write down a perfectly awful symbolic solution). For each of the real solutions, compute the x coordinate of the corresponding interesection point with the equation x*a^3*b1*d1*e1 - x*a^3*c1*d1^2 - x*a^2*b*a1*d1*e1 + x*a^2*c*a1*d1^2 - [much deleted -- see below. --djr] b*e*f*a1^3 + c*d*f*a1^3 (Phew!) Of course there are ways to simplify this. You could translate, rotate, and dilate so that one of the ellipses is simply the circle x^2+y^2=1, and you could assume the other ellipse is centered at (h,0) for some h. This would make the formulas much simpler. dave (moderator) ============================================================================== From: "Grudszus, Stefan" Subject: Q: general ellipse intersection formula (2) Date: Thu, 16 Sep 1999 08:31:34 +0000 Newsgroups: [missing] To: Dave Rusin Dear Dave! Firstly, thank you a lot for your email answer, which helped me a lot trying to solve my ellipse intersection point formula problem. Secondly, I do have some questions, which I hopefully adress to you, as I expect you know the answers. > From the data given you can write each ellipse as the solution set to a > quadratic polynomial in x and y. > f1:=a*x^2+b*x*y+c*y^2+d*x+e*y+f > f2:=a1*x^2+b1*x*y+c1*y^2+d1*x+e1*y+f1 I calculated the parameters "a" through "f1" for both rotated and translated ellipses, so that they are expressions of my ellipse parameters mx1, my1, ry1, ry1, phi1, ... phi2.. > You can eliminate one of these two from the two quadratic equations, > to get a single polynomial of degree 4 to solve. I simply used your formula for the quartic q4(y)= 0. > Solve that quartic and plug in to get the other variable; > each real root of the quartic gives one of the points of intersection. Here comes up the first problem. I used your formula for the quartic q4(y) and solvedq4(y)= 0 for y using our language IDL from RSI, which finds the roots using Laguerre's method. But even for very simple ellipses (centered an not rotated) I always get 4 purely real (2 double solutions) or 4 purely imaginary (pairwise conjugated complex) or 4 complex (also pairwise conjugated complex) solutions!!! As you mailed me, I should use real y solutions to q4(y), but in the general case (despite there are 4 points of intersect for test ellipses) there are always pairwise conjugated complex solutions! So, I cannot proceed, as I do not understand this phenomena. > Any symbolic algebra program will do elimination. Here is the output from > one of them. ... Which program did you use. Is the formula correct??? > For each of the real solutions, compute the x coordinate > of the corresponding interesection point with the equation ... I would be lucky, to be able to use this formula again. So can you please help me again?? I am totally confused! But I want to tackle this problem. In expectation of your soon help. Yours sincerefully, Stefan -------------------------------------------------- Dipl. Ing. CCTV Videosensorik Stefan Grudszus plettac electronic security GmbH Würzburger Str. 150 90766 Fürth Germany Phone: +49 (0)911 75884 363 Fax : +49 (0)911 75884 121 eMail: S.Grudszus@plettac-electronics.de -------------------------------------------------- ============================================================================== From: Dave Rusin Subject: Re: Q: general ellipse intersection formula (2) Date: Thu, 16 Sep 1999 13:36:59 -0500 (CDT) Newsgroups: [missing] To: s.grudszus@plettac-electronics.de You're right! I made a very silly mistake. Look what I wrote: > If the ellipses are described by the euqations > f1:=a*x^2+b*x*y+c*y^2+d*x+e*y+f; > f2:=a1*x^2+b1*x*y+c1*y^2+d1*x+e1*y+f1; If you type this into Maple exactly as written you get a very unexpected display. How silly of me to use f1 as the name of the constant term in the second polynomial when f1 is also the name of the first polynomial! If I do this instead, f2:=a1*x^2+b1*x*y+c1*y^2+d1*x+e1*y+fq; then things will work better. I had previously computed with Magma but Maple can do this too, and just as quickly: > readlib(eliminate): > temp:=eliminate({f1,f2},x): The output is very long but the command lprint(temp[2][1]); gives something equivalent to this quartic: f*a*d1^2+a^2*fq^2-d*a*d1*fq+a1^2*f^2-2*a*fq*a1*f-d*d1*a1*f+a1*d^2*fq+(e1*d^2* a1-fq*d1*a*b-2*a*fq*a1*e-f*a1*b1*d+2*d1*b1*a*f+2*e1*fq*a^2+d1^2*a*e-e1*d1*a*d- 2*a*e1*a1*f-f*a1*d1*b+2*f*e*a1^2-fq*b1*a*d-e*a1*d1*d+2*fq*b*a1*d)*y+(e1^2*a^2+ 2*c1*fq*a^2-e*a1*d1*b+fq*a1*b^2-e*a1*b1*d-fq*b1*a*b-2*a*e1*a1*e+2*d1*b1*a*e-c1 *d1*a*d-2*a*c1*a1*f+b1^2*a*f+2*e1*b*a1*d+e^2*a1^2-c*a1*d1*d-e1*b1*a*d+2*f*c*a1 ^2-f*a1*b1*b+c1*d^2*a1+d1^2*a*c-e1*d1*a*b-2*a*fq*a1*c)*y^2+(-2*a*a1*c*e1+e1*a1 *b^2+2*c1*b*a1*d-c*a1*b1*d+b1^2*a*e-e1*b1*a*b-2*a*c1*a1*e-e*a1*b1*b-c1*b1*a*d+ 2*e1*c1*a^2+2*e*c*a1^2-c*a1*d1*b+2*d1*b1*a*c-c1*d1*a*b)*y^3+(a^2*c1^2-2*a*c1* a1*c+a1^2*c^2-b*a*b1*c1-b*b1*a1*c+b^2*a1*c1+c*a*b1^2)*y^4 and the command lprint(temp[1][1]) gives this formula for the x coordinate: x = -(a*fq+a*c1*y^2-a1*c*y^2+a*e1*y-a1*e*y-a1*f)/(a*b1*y+a*d1-a1*b*y-a1*d) Let's test this: here's an ellipse very nearly the circle at the origin with radius 10 (I'm adding some random noise): > su:={a=10,b=1,c=11,d=-1,e=2,f=-1003}: Here's a small pertubation of the ellipse which is twice as wide and half as tall: > su2:={a1=2,b1=-1,c1=39,d1=2,e1=2,fq=-999}: > subs(su union su2, temp[2][1]); 2 4 3 -175252 y - 5898810 y + 136024 y + 9992 y + 63678146 > fsolve("); -4.872380534, -4.472486972, 4.550429218, 4.720980663 giving four y-coordinates roughly where we expect them to be. The x-coordinates are given by > xcoor:=subs(su union su2, temp[1][1]): so that we have an intersection at, for example > subs(y=-4.872380534, xcoor); x = -8.380860050 I hope everything is now in order. Sorry about the silly mistake! dave