Printable View
I need to draw an ellipse through point P2, which is the apex of either the minor or major axis, and two other points P1 and P3. Am I correct in believing there can only be one such ellipse? And if so, how can I find its centre, eccentricity and major axis length?
Attachment 108763
BB
EDIT: It would be useful to know if the apex is the end of the major or minor axis, although I could calculate that if necessary from the above data.
You can draw a circle through those three points. This circle is an ellipse with zero eccentricity. There are many other ellipses with varying eccentricity that can be draw through these three points.
Thank you Logophobic. I forgot that there must be a circle through any 3 points, including my seasonally coloured P1, P2, P3. I am still wondering though, given the constraint that P2 must lie on either the minor or the major axis, can there be more than one non-circular ellipse through the three points?
Yes. Choose any arbitrary point Q and draw line QP2. For any point A on line QP2 there exists points B and C such that |AP1| + |BP1| = |AP2| + |BP2| and |AP3| + |CP3| = |AP2| + |CP2|. For some A on QP2, points B and C coincide. Then points A and B/C are the focal points of an ellipse through points P1, P2 and P3 with point P2 being the apex of the major axis.
I cannot prove this, but I think you'll agree with my logic here.
Actually, I'm certain that this is not universally true. Still, there is definitely more than one ellipse here.
Just to be clear, |AP1| is the distance from point A to point P1.
Thanks for that, but without proof I'm not quite convinced yet. On playing around with my home-made ellipse tool, it seems to me that when P1 and P3 are symmetrically disposed either side of P2, there will be many (if not infinite) solutions. But I'm not sure when they are not symmetrical.Quote:
Originally Posted by Logophobic
I'd like to try something with the well-known ellipse formula
where a and b are the horizontal and vertical semi-axes, and x and y are the coordinates of any point on the ellipse relative to the centre. The ellipse can start in the canonical position (with the centre as the origin and the major axis horizontal), because the solution, if found, can be translated/rotated to any position.Code:(x/a)^2 + (y/b)^2 = 1
Attachment 108833
P2 now lies at (a, 0) relative to the origin, but a is itself initially unknown.
For P1 and P3 we can write:
I'd appreciate any advice on how to solve these simultaneous equations for a and b. If we assume a = b, we ought to get the circle solution. I suspect there will be many or infinite solutions when x3 = x1 and y3 = - y1. But in other cases, will there be a unique non-circular solution?Code:(x1 / a)2 + (y1 / b) ^ 2 = 1
(x3 / a) ^ 2 + (y3 / b) ^ 2 = 1
EDIT: after a night's sleep I realize that the above equations are incorrect. x1 an x3 are known in the same coordinate space as x2, and it is the centre of the eillipse which is unknown, at least in the x direction. So the origin is now ((x2 - a), 0) and the equations for P1 and P2 become:
Any ideas about how to proceed from here?Code:((x1 - x2 + a) / a) ^ 2 + (y1 / b) ^ 2 = 1
((x3 - x2 + a) / a) ^ 2 + (y3 / b) ^ 2 = 1
EDIT2
I'm about 50 years out of practice with this kind of thing, but here's my own attempt at the above equations:
Attachment 108877
I hope I haven't made any mistakes in there. The outcome looks to me like an ordinary quadratic equation, but my algebra is too rusty to solve it or draw conclusions from it. All opinions welcome.
BB
Translate the points so P2 is at the origin. Suppose you wish to make an ellipse with center (0, A). The resulting ellipse is
((x-A)/A)^2 + (y/B)^2 = 1,
for some B^2>0. Let P1=(u, v), P3=(r, s), and for convenience set u'=u-A, r'=r-A. P1 and P3 are on this ellipse if and only if
(u'/A)^2 + (v/B)^2 = 1
(r'/A)^2 + (s/B)^2 = 1
Here's a mess of algebra to solve for A in terms of u, v, r, s. Multiplying the first equation by s^2, the second equation by v^2, and subtracting gives
(u's/A)^2 + (vs/B)^2 - (r'v/A)^2 - (vs/B)^2 = s^2 - v^2
=> ((u's)^2-(r'v)^2)/A^2 = s^2 - v^2
=> (s^2 - v^2) A^2 = ((u's)^2-(r'v)^2)
=> (s-v)(s+v) A^2 = (u's-r'v)(u's+r'v)
= (us-As-rv+Av)(us-As+rv-Av)
= (us-rv-A(s-v))(us+rv-A(s+v))
=> A^2 = ((us-rv)/(s-v) - A)((us+rv)/(s+v) - A)
= (us-rv)(us+rv)/((s-v)(s+v)) - A[(us-rv)/(s-v) + (us+rv)/(s+v)] + A^2
=> A = [((us)^2-(rv)^2)/(s^2-v^2)]/[2(us^2 - rv^2)/(s^2-v^2)]
=> A = ((us)^2-(rv)^2)/(2(us^2 - rv^2)).
Given A, it's easy to solve for B^2, since there's no cross term. For instance, B^2 = (Av)^2 / (2Au - u^2) works. You can check (I told Mathematica to do it) that these values for A and B^2 indeed give equality in general. Hence we can compute the unique solution A, B^2 to the above system of equations. There's some degeneracy when we divide by zero in the above steps, but I'll ignore that for now.
Using this, we can produce an example of Logophobic's behavior, and indeed solve the problem in general. For the example, start with P2=(0, 0), P1=(1, 0), P3=(2cos(1), 2sin(1)). Rotate counterclockwise about the origin an angle t to get P2=(0, 0), P1=(cos(t), sin(t)), P3=(2cos(t+1), 2sin(t+1)). For this t, the above A and B^2 are
A = sin(1) sin(1+2t)/(2cos(t) sin(1+t)^2 - cos(1+t) sin(t)^2)
B^2 = 2sec(t)^2 sin(1)^2 sin(1 + 2 t)^2 / ((-2+cos(1)-2cos(2+2t)+cos(1+2t)) (2sin(1+t)^2-cos(1+t) sin(t) tan(t)))
These are both positive for 0.1 <= t <= 0.5, say, giving an ellipse (x/A-1)^2 + y^2/B^2 = 1 passing through (0, 0), (cos(t), sin(t)), and (2cos(t+1), 2sin(t+1)). Hence rotating (x, y) by t, i.e. replacing x with x cos(t) - y sin(t) and y with x sin(t) + y cos(t), we'll get an ellipse passing through (0, 0), (1, 0), (2cos(1), 2sin(1)). Here's a .gif animating the solutions for a full 0 <= t <= pi using these methods.
So, it seems you'll have to add another constraint if you want uniqueness.
Thanks very much for your insights Jemediah. Unfortunately there are several things I do not understand about your explanation. When you say that we can compute "the unique solution A, B^2 to the above system of equations", I take it this does not mean there is a unique solution for given values of u, v, r and s. Never mind, the gif animation is convincing enough: there is no unique solution to the conditions I stated, so an extra constraint is needed.
What I am trying to do is to design a curve drawing tool which provides different kinds of arc: parabola, ellipse, circle, catenary and possibly others. For the moment an elliptical arc is my priority, since I suspect the others will give less problems (and you already helped me with the catenary some time back). I would like the tool to work if possible with only three control points per arc.
The cubic spline provided by GDI+ plots a curve through three points, but the result tends to look like a parabola. Besides, the apex of the parabola is not necessarily at at a user-dragged point. The curve has a mind of its own, so to speak, making it less intuitive to use. Here is an example of how it looks. The three points P1, P2 and P3 are draggable but the apex of the curve moves away from the dragged point P2.
Attachment 108991
I would like the curve to be elliptical with P2 as a major or minor vertex. The point m in the above image is the midpoint between P1 and P3. I wonder if a requirement for the major or minor axis of the ellipse to pass through point m as well as P2 would be adequate as an extra constraint for a unique (non-circular) solution. If so, I would be grateful for a suggestion for adapting the equations posted so far.
Here is the VB.Net code which I used to produce the above interactive image:
Wishing everyone a happy new year,Code:Public Class Form1
Private P1, P2, P3, midpoint As DraggablePoint
Private P1dragging, P2dragging, P3dragging As Boolean
Private Sub Form3_Load(sender As System.Object, e As System.EventArgs) Handles MyBase.Load
Me.DoubleBuffered = True
Dim centre As New PointF(Me.ClientSize.Width / 2.0F, Me.ClientSize.Height / 2.0F)
P1 = New DraggablePoint(New PointF(centre.X - 100, centre.Y), 10, Color.Red, "P1", Me.Font)
P2 = New DraggablePoint(centre, 10, Color.Red, "P2", Me.Font, -20, -25)
P3 = New DraggablePoint(New PointF(centre.X + 100, centre.Y), 10, Color.Red, "P3", Me.Font)
midpoint = New DraggablePoint(centre, 5, Color.Gray, "m", Me.Font, -0, 0)
End Sub
Private Sub Form3_MouseDown(sender As Object, e As System.Windows.Forms.MouseEventArgs) Handles Me.MouseDown
For Each dp As DraggablePoint In {P1, P2, P3}
Dim r As New RectangleF(dp.X - 5, dp.Y - 5, 10, 10)
dp.Dragging = r.Contains(e.Location)
Next
End Sub
Private Sub Form3_MouseMove(sender As Object, e As System.Windows.Forms.MouseEventArgs) Handles Me.MouseMove
For Each dp As DraggablePoint In {P1, P2, P3}
If dp.Dragging Then
dp.X = e.X : dp.Y = e.Y
End If
Next
midpoint.X = (P1.X + P3.X) / 2
midpoint.Y = (P1.Y + P3.Y) / 2
Me.Invalidate()
End Sub
Private Sub Form3_MouseUp(sender As Object, e As System.Windows.Forms.MouseEventArgs) Handles Me.MouseUp
For Each dp As DraggablePoint In {P1, P2, P3}
dp.dragging = False
Next
End Sub
Private Sub Form3_Paint(sender As Object, e As System.Windows.Forms.PaintEventArgs) Handles Me.Paint
e.Graphics.SmoothingMode = Drawing2D.SmoothingMode.HighQuality
Using pn As New Pen(Brushes.Gray) With {.DashStyle = Drawing2D.DashStyle.Dot}
e.Graphics.DrawLine(pn, P2.Location, midpoint.Location)
End Using
Using pn As New Pen(Brushes.Black) With {.DashStyle = Drawing2D.DashStyle.Dash}
e.Graphics.DrawLine(pn, P1.Location, P3.Location)
End Using
e.Graphics.DrawCurve(Pens.Blue, {P1.Location, P2.Location, P3.Location})
For Each dp As DraggablePoint In {midpoint, P1, P2, P3}
dp.Draw(e.Graphics)
Next
End Sub
Private Class DraggablePoint
Public ReadOnly Property Location As PointF
Get
Return New PointF(X, Y)
End Get
End Property
Public X, Y As Single
Public Property Diameter As Single
Public Property FillColor As Color
Public Property Lable As String
Public Property Font As Font
Public Property OffsetX As Single
Public Property OffsetY As Single
Public Property Dragging As Boolean
Public Sub New(location As PointF, diameter As Single, fillcolor As Color, lable As String, font As Font)
X = location.X : Y = location.Y : _Diameter = diameter : _FillColor = fillcolor
_Lable = lable : _Font = font : _OffsetX = 5 : _OffsetY = 5
Dragging = False
End Sub
Public Sub New(location As PointF, diameter As Single, fillcolor As Color, lable As String, font As Font, offsetX As Single, offsetY As Single)
Me.New(location, diameter, fillcolor, lable, font)
_OffsetX = offsetX + 5 : _OffsetY = offsetY + 5
End Sub
Public Sub Draw(g As Graphics)
Dim r As New RectangleF(X - Diameter / 2.0F, Y - Diameter / 2.0F, Diameter, Diameter)
Using sbr As New SolidBrush(FillColor)
g.FillEllipse(sbr, r)
End Using
g.DrawEllipse(Pens.Black, r)
g.DrawString(Lable, Font, Brushes.Black, X + OffsetX, Y + OffsetY)
End Sub
End Class
End Class
BB
Actually it means exactly the opposite: given u, v, r, s, there is a unique solution A, B^2 to the system of two equations I listed, which can be computed using the formulas I gave. Physically, that means there is a unique ellipse (sometimes hyperbola, which is "essentially" the same thing; it's just an ellipse passing through infinity) passing through the origin and the two points (u, v), (r, s), where the origin is at the tip of a major or minor axis aligned with the x-axis. As mentioned, this glosses over some "degenerate cases", eg. when u=v=r=s=0. If you care enough you can trace through my derivation and check when I divide by zero and handle those cases separately, though they're certainly "rare" in practice. My guess is that non-colinearity is all that's needed.Quote:
When you say that we can compute "the unique solution A, B^2 to the above system of equations", I take it this does not mean there is a unique solution for given values of u, v, r and s.
Solving that system is the hard part. In general, there is a unique ellipse (or hyperbola) passing through three fixed points where the second one is at the tip of a major or minor axis with specified (though arbitrary) orientation. To compute it, first subtract P2.x and P2.y from all three points' x and y coordinates; then rotate the resulting points by your arbitrary angle t; then compute A and B^2 from these u, v, r, s; then plug this A and B^2 into the ellipse equation at the top of my first post; then replace x and y with x cos(t) - y sin(t) and x sin(t) + y cos(t), respectively; finally replace x with x-P2.x and y with y-P2.y. That is the equation of your ellipse for this angle t, passing through your original three points. This is essentially what I did to create the animation.
That would almost give uniqueness. It would specify the arbitrary angle above up to a multiple of pi. It would also have some degenerate cases.Quote:
I wonder if a requirement for the major or minor axis of the ellipse to pass through point m as well as P2 would be adequate as an extra constraint for a unique (non-circular) solution. If so, I would be grateful for a suggestion for adapting the equations posted so far.
Sure. Here's pseudocode; I can write some real code if necessary, it's just more time-consuming. I'll assume "m" is the midpoint of P1 and P2. There will undoubtedly be edge cases not covered here, like when m=P2. Hopefully it's not "buggy".Quote:
If so, I would be grateful for a suggestion for adapting the equations posted so far.
INPUT: P1, P2, P3
OUTPUT: Coordinates of a rectangle bounding an ellipse passing through P1, P2, P3, with P2 being a major or minor axis vertex, and with the midpoint of P1 and P3 being on the axis touching P2.
ALGORITHM:
// Recenters problem with P2 at the origin
1. Set (u, v) = (P1.x - P2.x, P2.y - P2.y); (r, s) = (P3.x - P2.x, P3.y - P2.y).
2. Let m be the midpoint of (u, v) and (r, s).
3. Let t=atan2(m.y, m.x) = the angle m makes with the +x axis measured counterclockwise.
// Fixes major/minor axis touching P2 as the x-axis. Note m is rotated to the +x-axis.
4. Rotate (u, v) and (r, s) counterclockwise about the origin by -t.
// Computes coefficients of (unique) ellipse/hyperbola solving this special case
// The center is (0, A), the x-axis-aligned axis has length A, the y-axis-aligned axis has length sqrt(B2)
5. Set A = ((us)^2-(rv)^2)/(2(us^2 - rv^2)), B2 = (Av)^2 / (2Au - u^2).
// Translate to the original situation
6. Let R be a rectangle centered at (0, A) with width 2A and height 2*sqrt(B2).
7. Rotate each vertex of R counterclockwise about the origin by t. Add P2 to each vertex of R.
8. Return R.
Thanks again Jemediah. I've made an initial try at implementing the algorithm. Since a rectangle structure can't be rotated in dotnet, it seemed more convenient to make a function to return the values of A and B plus the location of the ellipse centre. Since the angle of inclination (t) is known, it should be possible to plot the ellipse arc from this data in GDI+. However, I haven't been able to get it working convincingly yet. Here's what I have been trying:
To test this code, I added this to the Paint event in the code of post #7 to draw the calculated centre.Code:Private Structure EllipseArguments
Public A As Single
Public B As Single
Public centre As PointF
Public angle As Single
End Structure
Private Function RotatePoint(X As Single, Y As Single, radians As Double) As PointF
Dim rX As Single = CSng(X * Math.Cos(radians) - Y * Math.Sin(radians))
Dim rY As Single = CSng(X * Math.Sin(radians) + Y * Math.Cos(radians))
Return New PointF(rX, rY)
End Function
Private Function FindEllipse(p1 As PointF, p2 As PointF, p3 As PointF, mid As PointF) As EllipseArguments
'Recenter problem with P2 at the origin
Dim u As Single = p1.X - p2.X
Dim v As Single = p1.Y - p2.Y
Dim r As Single = p3.X - p2.X
Dim s As Single = p3.Y - p2.Y
Dim m As New PointF(mid.X - p2.X, mid.Y - p2.Y)
Dim t As Double = Math.Atan2(m.X, m.Y)
'Fix major/minor axis touching P2 as the x-axis.
u = RotatePoint(u, v, -t).X
v = RotatePoint(u, v, -t).Y
r = RotatePoint(r, s, -t).X
s = RotatePoint(r, s, -t).Y
'Compute coefficients of (unique) ellipse/hyperbola solving this special case
'calculate A = ((us)^2-(rv)^2)/(2(us^2 - rv^2)), B2 = (Av)^2 / (2Au - u^2)
Dim A As Single = CSng((((u * s) ^ 2 - (r * v) ^ 2)) / (2 * (u * (s ^ 2) - r * (v ^ 2))))
Dim B As Single = CSng(Math.Sqrt(((A * v) ^ 2) / (2 * A * u - u ^ 2)))
'Translate to the original situation
Dim centre As New PointF(A, 0)
Dim eArgs As EllipseArguments
eArgs.centre = RotatePoint(centre.X, centre.Y, t)
eArgs.centre.X += p2.X
eArgs.centre.Y += p2.Y
'etc.
eArgs.A = A
eArgs.B = B
eArgs.angle = CSng(t * 180 / Math.Pi)
Return eArgs
End Function
But the dot painted in the Paint sub clearly isn't the centre of a plausible ellipse through the points. So I wonder if my code is a correct interpretation of the algorithm. I've tried to follow the algorithm as closely as possible apart from an apparent typo in Step 1 (I assume P2.y-P2.y should be P1.y - P2.y).Code:Dim ea As EllipseArguments = FindEllipse(P1.Location, P2.Location, P3.Location, midpoint.Location)
e.Graphics.FillEllipse(Brushes.Cyan, ea.centre.X - 6, ea.centre.Y - 6, 12, 12)
BB
For what its worth, three points defining an ellipse would be F1, F2 and any P on the curve, as in the image in post #5. I think a tool with those parameters would be most useful.
After looking through the code, it's a faithful implementation (thanks for correcting the typo) except you've reversed the order of the arguments to atan2. Instead of "Dim t As Double = Math.Atan2(m.X, m.Y)" it should be "Dim t As Double = Math.Atan2(m.Y, m.X)". [I had actually typed "Let t=atan2(m.y, m.x)" since this is the more common, but not only, convention.] If that doesn't fix it, I'll install VB.NET and debug my method.
Reversing the atan2 arguments didn't change anything obvious. I won't be able to debug this without help, so I'd appreciate assistance. I've attached all the code: 1 file VB.Net + 2 files C# in a zip. You only have to start a Forms project and drop your preferred version into the form code. For C#, build the 2 files. Then run. This version draws the ellipse and its bounding rectangle (when visible) as well as the ellipse centre. Just drag the red points around, an ellipse will flit by sooner or later!
@Logophobic: Once this is sorted out it shouldn't be too hard to make a version with 2F+P. I suspect it will in fact be easier.
Attachment 109027
BB
Hold it! I just discovered I forgot to convert the output of Atan2 to degrees for use in Graphics.RotateTransform. Let's see if it fixes the problem. BB (+forhead-slapping)
[edit] It still doesn't work correctly, but at least the angle is now correct which I hope will make it easier to spot the problem.
The reason I noticed the angle conversion error was that I decided to have a quick go at making a 2F+P ellipse tool as suggested by Logophobic. I used a VB.Net version of the pins-and-string method. It worked on the first attempt, except (for now obvious reasons) the angle. The completed code for the 2F+P ellipse is as follows:
The Distance function calculates Math.Sqrt(dX^2+dY^2). See the attached zip for full VB.Net code with a test form. The ellipse tools and the LabelledPoint are now in separate files from the test form.Code:Public Shared Function FindEllipse2(F1 As PointF, F2 As PointF, P As PointF) As EllipseTools.EllipseArguments
Dim ea As EllipseTools.EllipseArguments
ea.centre = New PointF((F1.X + F2.X) / 2.0F, (F1.Y + F2.Y) / 2.0F)
Dim f As Double = Distance(F1, F2) / 2
ea.A = CSng((Distance(F1, P) + Distance(P, F2)) / 2.0F) 'half the length of the "string"
ea.B = CSng(Math.Sqrt(ea.A ^ 2 - f ^ 2))
ea.angle = CSng(Math.Atan2(F2.Y - F1.Y, F2.X - F1.X) * 180 / Math.PI)
Return ea
End Function
edit several weeks later: on reviewing the thread, I just noticed that I accidentally deleted the Dim f line from the code posted above at some point. Now I have put it back. My apologies to anyone who tried running it. The zip should still be OK.
I wrote (0, A) when I meant (A, 0), so "Dim centre As New PointF(0, A)" should be "Dim centre As New PointF(A, 0)". With that change, when the line connecting P2 and m is horizontal, it generally works fine (though when the solution is a hyperbola or degenerate, the ellipse drawing routine gets confused, as you'd expect). There's definitely something wrong with rotation. Continuing to debug; I'll update this when I have more.
Edit: Oh of course,
and the r, s version makes no sense. In the second line, u will have been changed, screwing up v. It should be something likeCode:u = RotatePoint(u, v, -t).X
v = RotatePoint(u, v, -t).Y
The results are much more reasonable, though not perfect. Rotation is still screwed up.Code:Dim up As Single = RotatePoint(u, v, -t).X
Dim vp As Single = RotatePoint(u, v, -t).Y
Dim rp As Single = RotatePoint(r, s, -t).X
Dim sp As Single = RotatePoint(r, s, -t).Y
...
Dim A As Single = CSng((((up * sp) ^ 2 - (rp * vp) ^ 2)) / (2 * (up * (sp ^ 2) - rp * (vp ^ 2))))
Dim B As Single = CSng(Math.Sqrt(((A * vp) ^ 2) / (2 * (A * up) - (up ^ 2))))
Edit: Oh, it was just the issue you mentioned earlier, the RotateTransform line should be "e.Graphics.RotateTransform(.angle * 180 / Math.PI)" rather than "e.Graphics.RotateTransform(.angle)". With all of these changes, I now have a beautiful ellipse with the desired properties, except when the solution is a hyperbola or degenerate.
A few other comments: sorry about the (0, A)/(A, 0) bug; brain fart. I should have caught the uv/rs bug earlier when I looked through your code but it just didn't occur to me. It seems m is always the center of the resulting ellipse, which I'm sure could be proved quickly, though I won't take the time. I mention it since this observation might yield a cleaner derivation of A and B^2. The "That would almost give uniqueness" comment in post #8 is untrue; both choices give the same ellipse in the end.
No apologies are called for from you jemediah. I thought of declaring new variables instead of recycling u, v, r and s, but then forgot to try it. I wondered about 0,A but the effect of reversing the arguments wasn't clear so I left them as they were. If I hadn't made the other errors the effect would have been more self-evident.
Now that the algorithm is working more predictably I have to look at it critically as a drawing tool. With m-P2 as the major axis it is perfect. When m-P2 is the minor axis, it does not work as I hoped. The major vertices of the ellipse shoot out sideways towards infinity as P2 approaches m, making the curve uncontrollable from a user viewpoint. The problem is that the algorithm is fixing the centre of the ellipse at m, which is too strict a constraint. My idea was that the centre only needs to be collinear with P2 and m. In other words, as P2 approaches m, the ellipse centre should slide past m along the line P2-m.
Here is an illustration of what it might look like:
Attachment 109089
What I have done here is to draw an arbitrary ellipse, and dragged points P1, P2 and P3 onto appropriate points of the ellipse. It appears empirically (by dragging the points around) that solutions of this kind are only possible when P1 and P3 are symmetrically disposed on either side of P2, and that the line P1-P3 is parallel to the major axis. I'm still not sure this means there is a unique solution for the ellipse in this symmetrical configuration of p1, p2, p3.
In any case, it's clear I have to decide what kind of behaviour I want the curve to have when P2 is the vertex of the minor axis. Maybe the constraint of the centre of the ellipse lying on line P2-m was not a good choice.
BB
With P3 the reflection of P1 across the P2 axis, there are infinite ellipses:
Attachment 109101
Blue - P2 is apex of major axis
Yellow - Circle
Green - P2 is apex of minor axis
Edit - Added pink ellipses
In the "normalized" form of post #6, a configuration with the line from P1 to P3 perpendicular to the axis from P2 corresponds to P2=(0, 0), P1=(u, v), P3=(u, -v), i.e. u=r, v=-s. Indeed, then us^2 - rv^2 = us^2 - u(-s)^2 = 0, so we divide by 0 to get A, which indicates this configuration is degenerate. The formula above is actually 0/0, so it can still give reasonable results "in the limit", but in any case we lose uniqueness. The two equations from #6 are then
(u/A-1)^2 + (v/B)^2 = 1
(u/A-1)^2 + (-v/B)^2 = 1,
so the second is redundant. We have simply v^2/(1 - (u/A-1)^2) = B^2 for any value of A whatsoever, so there's a whole "ray" (possibly line segment) of solutions (more if you allow hyperbolas), as Logophobic's drawing illustrates. It's somewhat unfortunate that "half the time" there is a unique solution to your problem in post #7, which is a hyperbola, and the "other half" of the time, there is a unique solution, which is an ellipse centered on m (as usual ignoring degenerate cases).
If you drew the hyperbolas, it would at least be more predictable from a user standpoint. It's clearer "why" the major axis gets arbitrarily large sometimes with hyperbolas added in; otherwise it just feels arbitrary and broken. Oh, you could also give the user control of "m" rather than forcing it to be the midpoint of P1 and P3. The algorithm works fine as-is with m not necessarily this midpoint.
Thanks for the suggestions. Plotting the hyperbolas would be interesting a geometrical demonstration but I think it would be a distraction for my present purpose. I only want to bend the line P1-P2-P3 into an arc in various ways for graphic purposes. It ought be possible to create asymmetrical elliptical arcs like the one shown in red below, by dragging on P2. The image below was again made by first drawing the ellipse and then dragging the points on to it -- the reverse of the required procedure. But after several tryouts this way, I am convinced that solutions must be possible.
Attachment 109117
For this purpose, I would like to consider abandoning the requirement for P2 to be at the tip of the minor axis (it's fine for the major axis). Visually, the minor vertex turns out not to be a very interesting point anyway on flattish arcs like this. But P2 must remain on the ellipse: the user drags the curve itself, not some remote control point. The ellipse centre c is some point on the line P2-m which is effectively "controlled" by P2. If an extra constraint is needed to get a unique solution, I would like to consider letting P2 be the point on the arc furthest perpendicularly from the line P1P2. It wouldn't surprise me if that raises some more mathematical difficulties. So unless someone is keen to rush in with an answer, I think I ought to take time to think about it myself.
I don't think the constraint you've said makes complete sense. You want P2 to be as far as possible perpendicularly from line P1P3 (typo), but on "the arc": you're using this constraint to define which arc you mean, so referencing "the arc" here is circular. Perhaps your hope is that among all possible ellipses passing through the three points, there is a unique one where P2 is a point on the ellipse furthest perpendicularly from the line P1P3? I can try to compute such a thing if you're sure it's what you want.
You could also let the user specify the eccentricity of the ellipse (the ratio A/B, say). Eccentricity 1 would give a circle, etc.
Yes, that's what I meant.Or you could say choose the elliptical arc whose furthest point from line P1P3 is closest to P2. But I suspect this won't be necessary.Quote:
Originally Posted by jemidiah
In the symmetrical situation which Logophobic illustrated in post #16, the arc can be circular.
I can't find any ambiguity in the asymmetrical situation I illustrated in post #18. Once P1 and P3 are positioned on the ellipse, there seems to be only one position where I can place P2 on the arc so that P2-m-c forms a straight line. I realize this isn't a proof that there is only one solution, but it seems reasonable evidence. What is more, the resulting position of P2 always appears to be (roughly) the furthest point of the arc from line P1P3; so I suspect that this will turn out to be a consequence rather than a constraint!
I'll try to sum up the desired behaviour of the tool (square brackets mean line length):
1. While [P2m] is less than [P1m]
- a. if [P1P2] = [P3P2] the arc is part of the circle through P1, P2 and P3.
- b. otherwise, P2, P2 and P3 are on the arc and the centre of the ellipse c is collinear with P2m and on the side of m opposite to P2.
2. When [P2m] = [P1m], the arc is a semicircle of radius [p1m] through P1, P2 and P3. The centre c coincides with m.
3. While [P2m] > [P1m] the arc is an ellipse through P1, P3 and P3 with P2 at the vertex of the major axis. The centre c coincides with m.
Situation 1b is the only one which needs resolving.
That's another one I would save for a more mathematical version of the tool.Quote:
You could also let the user specify the eccentricity of the ellipse (the ratio A/B, say). Eccentricity 1 would give a circle, etc.
regards, BB
Hmm. I think your problem is still underdetermined. I'll give a system of equations that describe your conditions and then have Mathematica solve them. This shows there is a single free parameter, so another constraint is needed. This'll be the last significant time I spend on this problem, I think, since the quarter is starting and I should focus on my research :).
I'll switch to the matrix representation of conic sections; it seems more natural here since I can't see an obvious "normal" form in the spirit of the "special case" I brute forced in post #6. An equivalent condition to P2 being at the furthest point in the m-P2 direction from P1-P3 on the ellipse is for the tangent line to P2 to be parallel to the line P1-P3; this is essentially the idea behind Lagrange multipliers; the key is that the normal vector, i.e. the gradient of the defining equation determining level sets, is then perpendicular to the line P1-P3. You've suggested this condition is equivalent to the P2-m-c colinear condition, which seems at least plausible, so I'll go with that.
Let A_Q be the matrix of the conic, which will be determined; this means the conic is given by the set of solutions (x, y) to X^T A_Q X = 0 where X=(x; y; 1) (a column vector) (see Wiki page). For the three points to be on the ellipse, we need P_i^T A_Q P_i = 0 for i=1, 2, 3 (where eg. P_1 = (P1.x; P1.y; 1)), three conditions. The matrix has 6 variables, but only 5 degrees of freedom, since we can multiply it by an arbitrary non-zero scalar without affecting anything. One way to kill this extra degree of freedom is, in the notation of the Wiki article, to set A+C=1, i.e. replace C with 1-A (rigorously this can be justified by dividing by A+C and substituting; A+C will roughly correspond to the sum of the major and minor axes, which is in any case non-zero for ellipses). Since the three conditions above each take care of 1 degree of freedom in the "general" case, we have 2 more degrees of freedom to kill.
In post #8, you had given two additional constraints, that P2 be at a vertex (one degree) and that m be on a major or minor axis (another degree), so heuristically there's a unique solution in the "general" case; my algebra showed this to be true.
Now, you've just given one additional constraint. I'll cast it in algebraic terms. The gradient of X^T A_Q X as a function of x, y at some point (x0, y0) is normal to the level set through (x0, y0), in the sense that it is perpendicular to the tangent line to the level set at P. Hence the condition is just that the gradient is perpendicular to P1-P2. The gradient at P2, again in the notation of the Wiki article, is just
(2A P2.x + B P2.y + D, 2C P2.y + B P2.x + E)
The algebra simplifies quite a bit if we just take P2=(0, 0) by translating it beforehand. Then P_2^T A_Q P_2 = 0 if and only if F=0, and the gradient is just (D, E). Letting P1=(u, v), P3=(r, s), the conditions then simplify to:
Au^2 + Buv + (1-A)v^2 + Du + Ev = 0,
Ar^2 + Brs + (1-A)s^2 + Dr + Es = 0,
Du + Ev = Dr + Es
Mathematica says there is one free parameter, which I chose to be E, with the remaining variables satisfying
A = (s(r-u)v + E(rs-uv))/((r-u) (ru + sv))
B = (-(r - u)(su + rv) + E (-r^2 + s^2 + u^2 - v^2))/((r - u) (ru + sv))
D = E(-s + v)/(r - u)
I won't take the time to plot solutions, though it could be done from the bolded bits above. Again I'd expect "half" of them to be hyperbolas.
Jemediah, I greatly appreciate all the time you have putting into this. I've tried following your reasoning in conjunction with the Wiki page but I have to admit there is a lot I don't understand. I guess the best thing I can do now is to try to plot the solution in code. I hope you still have time to help me with a few questions about your latest post.
Can I take it that A and B still represent the lengths of the major and minor semi-axes respectively? If so, are the axes oriented vertically/horizontally as illustrated in post #18?
Is the gradient referred to as "just (D, E)" the normal to the tangent at P2, which is not something I need to calculate since it is perpendicular to P1P3. In that case, why do I need the value of D?
If E is a free parameter, I suppose I could just try substituting literal values to see what happens. E clearly does not represent eccentricity, because otherwise E=0 would give equal values for A and B (the circle). But maybe a value range 0-1 would be a good place to start.
BB
Ah, sorry, I should have explained my notation more explicitly. A, B, C, D, E, and F are all parameters to the general quadratic curve equation,
Ax^2 + Bxy + Cy^2 + Dx + Ey + F = 0,
listed on the linked page. The matrix A_Q (I don't know why they chose "A" there; it has nothing to do with the A above) organizes these in a nice general way, though I didn't end up needing that after I refactored part of my post.
A and C are no longer major/minor axis lengths, unfortunately. Here's how to use the info in that post:
1. Subtract P2 from P1 and P3; set P1=(u, v); P3=(r, s). Pick any value for E (only some will give reasonable answers). Compute A, B, D using the bolded formulas. Set C=1-A, F=0. The equation of your conic is
Ax^2 + Bxy + Cy^2 + Dx + Ey + F = 0.
2. Replace x with x+P2.x and y with P2.y.
3. Plot the resulting conic. There are a number of ways to do this. My first thought was diagonalizing the final A_Q matrix, but that's probably too much overhead. A quick google search yielded this set of notes, in particular section 4 on polar plots. That should be pretty straightforward to implement, though you'll have to complete the square in x and y to get it in a form without Dx or Ey.
To answer the rest of your questions: the axes are no longer oriented vertically/horizontally as illustrated in post #18; the gradient (D, E) is indeed normal to the tangent at P2, at least that's what the third constraint equation I wrote says. You need to know D because it's used in the general quadratic curve equation above. E is definitely not eccentricity. You should probably try both positive and negative values; I have no idea what a reasonable range would be. It's just an algebraic parameter at this point without a lot of physical meaning.
Thanks for the clarification. I think I should be able to implement polar formula given in section 4 of the document you linked to.
In your step 3, however, I'm not sure how to get from the "big letter" formula:
Ax^2 + Bxy + Cy^2 + Dx + Ey + F =0 (C=1-A, F=0)
to the "small letter" one used in that section:
ax^2 + bxy + cy^2 = +/-1
Is it correct to assume that in this case a = A, b = B and so on? It's been a long time since I completed a square, but it would probably do me good to revive my knowledge in that respect.
BB
Whoops, the cross term is more complicated than just a simple complete-the-square problem. Noting that
A(x-D/(2A))^2 = Ax^2 - Dx + D^2/(4A)
A(y-E/(2C))^2 = Ay^2 - Ey + E^2/(4C)
will get rid of the Dx and Ey terms (at the cost of changing F), but then the Bxy term is
B(x-D/(2A) + D/(2A))(y-E/(2C) + E/(2C)),
which will reintroduce them when expanded. To handle the cross term, one should probably consider replacing x with x' = p(x-x0), y with y'=q(y-y0) for some constants x0, y0, p, and q; tuning them right should get rid of the Dx and Ey terms entirely. But that's overcomplicated; instead, just follow the derivation given in that section but use it on the more general form (the "big letter" formula); you'll end up with a quadratic equation in r, which can be solved with the quadratic formula anyway, no extra linear algebra needed.
If the midpoint (m) of P1-P2 lies on either axis, then it must be the center of the ellipse. This should be easy to prove.
Attachment 109223
@Logophobic: Yup. I think bb and I both observed that somewhere above; the working version of the code discussed in post #14 makes it clear empirically. Here's a proof using my formula from post #8. Here P2=(0, 0), P1=(u, v), P3=(r, s), the midpoint of P1-P3 lies on the x-axis, and the center of the ellipse is at (A, 0) for (excepting degenerate cases) A = ((us)^2-(rv)^2)/(2(us^2 - rv^2)). Since v=-s, this simplifies to
A = s^2(u^2-r^2)/(2s^2(u-r))
= (u-r)(u+r)/(2(u-r))
= (u+r)/2,
so indeed the center is ((u+r)/2, 0), which is precisely the midpoint of P1-P3. My hope was that a simpler derivation might exist for this fact which might give the formulas in my post #8 more cleanly, but I didn't explore it.
This post suggests that bb was still having trouble with P2 as apex of minor axis. Case 1 above is incorrect. It should be:Quote:
Originally Posted by boops boops
1. While [P2m] < [P1m] the arc is an ellipse through P1, P2 and P3 with P2 at the vertex of the minor axis. The centre c coincides with m.
Note that P1P3 is the major axis when perpendicular to P2m (1a in quoted post). Similarly, P1P3 is the minor axis when perpendicular to P2m in case 3.
Edit: I haven't been following the math you've shown, but it occurred to me that you haven't been moving the center to the origin.
Attachment 109271
From (x/a)2 + (y/b)2 = 1 we get b2 = y2 / (1 - (x/a)2), which is valid for all a > x.
The behaviour described in your edit is intentional.Quote:
Originally Posted by Logophobic
The above diagram illustrates the result of the (corrected) code of post #9. That method works satisfactorily when P2 is on the major axis (case 3) but not on the minor axis (case 1). The problem is that the tool becomes unusable as the user pulls point P2 towards the line P1-P3. The ellipse becomes excessively eccentric and eventually transforms into a hyperbola, so from the user's point of view it is unstable.
So I proposed abandoning the requirement for P2 to be on the minor ellipse axis. Sorry about deviating from the thread title, but that happens:o. As a result, the centre of the ellipse is no longer at the midpoint of P1-P3. Instead, it must slide along line P2-m in response to movements of P2, as illustrated in post #18. That is the intended behaviour of case 1b. Jemediah has analysed this situation mathematically. Now it's up to me to try implementing his solution, which may take a while yet!
BB
OK, now I understand your problem. I also have a solution. When |P2m| < |P1m|, move c along P2m so that |P2c| = |P1m|. This should provide a unique solution, and will also allow a smooth transition between case 1 and case 3. You don't need case 1a, it has a noticeably different behavior than the other cases.
1. While [P2m] < [P1m], the arc is an ellipse through P1, P2 and P3. The centre c is collinear with P2m on the side of m opposite to P2; [P2c] = [P1m]
2. When [P2m] = [P1m], the arc is a semicircle of radius [p1m] through P1, P2 and P3. The centre c coincides with m.
3. While [P2m] > [P1m], the arc is an ellipse through P1, P2 and P3 with P2 at the vertex of the major axis. The centre c coincides with m.
Keeping [P2c] = [P1m] strikes me as a very good idea. It supplies the missing constraint, and above all it makes it easy to calculate where the centre of the ellipse will be.Quote:
Originally Posted by Logophobic
This leaves the question of how to plot an ellipse with a given centre through three points (given that c is on the other side of m compared to P2). I presume it can be done analytically similarly to Jemediah's recommended method. But I would certainly like to hear if there are other approaches in this situation, particularly ones which directly yield the major and minor axes and the rotation of the ellipse.
EDIT: with a bit of Googling, it appears D=E=0 in the general equation Ax^2 + Bxy + Cy^2 + Dx + Ey + F = 0 when the centre is at the origin. This answers how to apply the known centre to Jemediah's method, since the centre can be translated to the origin.
BB
EDIT: In the light of day, it was clear that the absurdly large absolute values of A, B and C must mean an error in my coding of the formula. Now testing with missing brackets added as in red below.
Hi Jemediah,
At last I got round to trying to put your advice of posts #21 and 23 (with the polar coordinate function) into practice. The following function returns an array of PointFs which can be plotted in GDI+ with DrawPolygon or DrawCurve:
The value of E ranges from -1000 to 1000. ScaleFactor is a fudge needed to make the output visible; it typically takes a value of around 1 million. In other words, the values generated are on the small side (speaking logarithmically ;)) The typical output forms were a hyperbola or a "ninja star" like those below:Code:Public Shared Function FindEllipse3(P1 As PointF, P2 As PointF, P3 As PointF, E As Double, scaleFactor As Integer) As PointF()
'Translate the origin to P2:
Dim u As Single = P1.X - P2.X
Dim v As Single = P1.Y - P2.Y
Dim r As Single = P3.X - P2.X
Dim s As Single = P3.Y - P2.Y
'Calculate coefficients for the general ellipse equation:
'A = (s(r-u)v + E(rs-uv))/((r-u) (ru + sv))
Dim A As Double = (s * (r - u) * v + E * (r * s - u * v)) / ((r - u) * (r * u + s * v))
'B = (-(r - u)(su + rv) + E (-r^2 + s^2 + u^2 - v^2))/((r - u) (ru + sv))
Dim B As Double = (-(r - u) * (s * u + r * v) + E * (-r ^ 2 + s ^ 2 + u ^ 2 - v ^ 2)) / ((r - u) * (r * u + s * v))
Dim C As Double = 1 - A
'D = E(-s + v)/(r - u)
Dim D As Double = E * (-s + v) / (r - u)
'F = 0
'Plot points with the polar formula: r = sqrt( 1 / | a*cos2(th) + b*sin(th)*cos(th) + c*sin2(th) | )
Dim pfs As New List(Of PointF)
For i As Integer = 0 To 359 Step 5
Dim theta As Double = i * Math.PI / 180
Dim sin As Double = Math.Sin(theta)
Dim cos As Double = Math.Cos(theta)
Dim radius As Double = Math.Sqrt(1 / Math.Abs(A * (cos ^ 2) + B * sin * cos + C * (sin ^ 2))) * scaleFactor
Dim radX As Single = CSng(radius * cos)
Dim radY As Single = CSng(radius * sin)
pfs.Add(New PointF(radX + P2.X, radY + P2.Y))
Next
Return pfs.ToArray
End Function
Attachment 109499Attachment 109501
No sign of an ellipse through P1 , P2 and P3, alas, at any value of E or relative positions of P1, P2 and P3. The value of D is not in fact being used, so maybe that indicates something wrong in the procedure or of my interpretation of it.
Now I am unsure what quadratic equation you refer to. Is it the one used to derive the polar coordinate formula I have used in the code above, with r representing the polar radius? I assume the quadratic formula is the well known (or well forgotten in my case) formula for solving ax^2 + bx + c = 0.Quote:
you'll end up with a quadratic equation in r, which can be solved with the quadratic formula
I hope you can find time to look at this.
cheers, BB
I'll try to take a look in the next couple of days.
[duplicate post deleted]
[I seem to be getting a rash of duplicate posts. Deleted text.]
Thanks, there's no hurry but I hope to get this working some day:rolleyes:.Quote:
Originally Posted by jemidiah
Meanwhile, with brackets added as in the edit to post #32, I now get an ellipse for some value ranges of E. The range depends on the position of P2 and is positive or negative according to which side it lies of P1P3. As a rough example, the plot is elliptical when E=- 150 to -50 for some positions of P2 north of P1P3. However:
1. The ellipse is always centred on P2.
2. It has no obvious other relation to P1 and P3.
3. A scaling factor of 100 is still needed to get a plot at what looks like a suitable size.
I'll report further progress if any as edits to this post.
EDIT1: I found advice here for solving the general ellipse equation by completing the square. It looks like just what I need, so I'll see if I can code it up.
EDIT1a: but the equation there doen't have Bxy term; now I see the relevance of post #25.
EDIT2: found a possibly better approach for the complete general conic equation here.
EDIT2a: curiously it "inverts" B and C: Ax^2 + By^2 + Cxy + Dx + Ey + F = 0. But the page gives ways of deriving the ellipse axes, rotation and centre from the formula: everything I need for my EllipseArguments type.
EDIT3: an ellipse that fits all the requirements except one in Post #18 is Steiner's circumellipse. It is a unique ellipse which passes through 3 points and has its centre at the centroid of the triangle formed. The tangent at P2 will be parallel to P1P3, and the centre will lie on the median P2m -- but on the "wrong" side of it. Still, I should like to figure out how to plot the circumellipse and see how it behaves in practice. This paper looks useful. Theorem 2.1 derives the foci of a triangle's inellipse although I need to learn to interpret the (vector?) notation. I can construct a triangle with the midpoints of its sides at P1, P2 and P3; and the inellipse of that is the circumellipse of triangle P1P2P3.
Sorry for not getting to this. I came down with the flu for the last week and am now playing catch-up. Hopefully soon.
Don't worry. I haven't had much time either. I think Steiner's circumellipse is going to be an interesting approach and I should try it out before going further. I anticipate it will show similar "undesirable" behaviour to your solution at #14, when P2 comes close to P1 or P3 and the triangle flattens out. But it could throw up some useful ideas.
BB
The formula given here for the Steiner inellipse looks simple enough, but I am having trouble with understanding the notation:
Attachment 109831
The formula for g seems clear enough because I can treat it as two formulas:
and its counterpart for y.Code:g.X = (z1.X + z2.X + z3.X) / 3
That won't work for the main square root formula however (I tried it and it gave incorrect results). I considered interpreting g^2, z1z2 etc. as dot products of the vectors g, z1 etc. but the square root would then only produce a single number which can't be added or subtracted from g. As for the possibility that they are cross products, I don't know how to calculate them and I wouldn't know what meaning to give to the square root.
Can anyone enlighten me on this?
EDIT: the article mentions complex analysis, so presumably g and the z's are treated as complex numbers. I still can't resolve the formula -- what do I do with all the imaginary terms? But I'll do a bit of reading up.
EDIT 2: OK, I've sorted out the complex numbers with the help of WikiPedia. Now I have a working inellipse + circumellipse! My code needs a lot of tidying up but I'll post a demo soon.
BB
There's a demo of the CircumEllipse tool in the attached zip. Here's how it looks:
Attachment 109959
You can drag the red points P1, P2 and P3 to any position. The point m is the midpoint of [P1P3], and c is the centroid of triangle [P1P2P3]. The yellow ellipse is the Steiner circumellipse of the inner triangle. Its axes and focal points F1 and F2 marked in green. The ellipse is at the same time the inellipse of the grey "circumtriangle" (maybe another name exists for this) of [P1P2P3], with midpoints at P1, P2 and P3.
The point P2R is the reflection of point P2 in the centroid i.e. the intersection of the median [P2m] with the opposite side of the ellipse. The length ratio of [P2m] to [mP2R] is 3 to 1. This is clear empirically, but it shouldn't be hard to prove geometrically. This ratio is important because it makes it possible find P2R from P2 and m without calculating the ellipse first; and conversely P2 can be found from P2R. In turn this makes it possible for P2 and P2R to swap roles as control points.
The red arc shown through P1, P2 and P3 above is what I call the "major arc" (>180 degrees). If you deselect the lower checkbox, the red arc becomes the "minor arc" through P2R, which then becomes the middle control point instead of P2.
Uncheck the upper checkbox to see the arc and control points only. The minor arc looks to me more useful as the basis for a drawing tool, but the major arc may turn out to have its applications too.
I don't think this solution perfect as a tool. For example, it would be more intuitive to drag the apex of the major axis as in the solution of post #12. There may be some way to combine the two approaches, but I'll put that off for some other time. In any case it's better than dragging the apex of the minor axis, so I'll take satisfaction with this solution and mark the thread resolved. But I'd be glad to hear any comments or suggested improvements.
Huge thanks to Jemediah and Logophobic, without whose help I would have no idea where to start.
BB
That's beautiful and it seems reasonably intuitive. Glad it worked out. Sorry I never got to debugging your earlier post; life got too busy for me to put in the time and effort.
Welcome to the forum and thanks for taking an interest. My question was about non-circular ellipses. I thought that would be obvious because a circle doesn't have an apex, but the potential confusion wasn't cleared up until post #3 of this thread. BBQuote:
Originally Posted by jimmyjjohn