## Fitting an axis-aligned ellipse to four points

Discussions related to mathematics, numerical methods, graph plotting etc.
RichardRussell
Posts: 533
Joined: Tue 15 Oct 2019, 09:10

### Fitting an axis-aligned ellipse to four points

As announced elsewhere, the ellipsefit example program supplied with BBC BASIC for SDL 2.0 (version 1.08a or later) is also fully compatible with BBC BASIC for Windows. Although there's nothing to stop a BB4W user from downloading BBCSDL (free) just to get this example, I've listed it here to save him/her the trouble:

Code: Select all

``````      REM Program to demonstrate deriving an axis-aligned ellipse from four points
REM Version 1.0 by Richard Russell, http://www.rtrussell.co.uk/, 17-Nov-2019
REM Compatible with both 'BBC BASIC for Windows' and 'BBC BASIC for SDL 2.0'

BB4W% = INKEY\$(-256) == "W"
VDU 23,22,640;500;8,20,16,128+8
ON ERROR PROCcleanup : IF ERR=17 CHAIN @lib\$+"../examples/tools/touchide" ELSE MODE 3 : PRINT REPORT\$ : END
ON CLOSE PROCcleanup : QUIT

IF BB4W% THEN
INSTALL @lib\$ + "GDIPLIB"
PROC_gdipinit
*FONT Arial, 18
ELSE
INSTALL @lib\$ + "aagfxlib"
OSCLI "FONT """ + @lib\$ + "DejaVuSans"", 18"
ENDIF
INSTALL @lib\$ + "arraylib"
OFF : VDU 5

DIM X(3), Y(3)

Cx = 500 + RND(300)
Cy = 350 + RND(300)
Rx = 200 + RND(100)
Ry = 200 + RND(100)

FOR I% = 0 TO 3
X(I%) = Rx * COS(Alpha) + Cx
Y(I%) = Ry * SIN(Alpha) + Cy
NEXT

title\$ = "Drag the points to resize the axis-aligned ellipse"
title% = (1280 - WIDTH(title\$)) / 2
fault\$ = "No ellipse passes through all four points"
fault% = (1280 - WIDTH(fault\$)) / 2
range\$ = "Ellipse is too big to draw"
range% = (1280 - WIDTH(range\$)) / 2

*REFRESH OFF
selected% = -1
oldX% = 0
oldY% = 0
oldB% = 0
REPEAT
CLS
GCOL 0 : VDU 30 : PLOT 0, title%, 0 : PRINT title\$;
code% = FNellipse4p(X(), Y())

FOR I% = 0 TO 3
IF BB4W% THEN
IF I% = selected% C% = &FF008000 ELSE C% = &FF0000FF
brush% = FN_gdipcreatebrush(C%)
PROC_gdipsector(brush%, X(I%), Y(I%), 16, 16, 0, 360)
PROC_gdipdeletebrush(brush%)
ELSE
IF I% = selected% C% = &FF008000 ELSE C% = &FFFF0000
PROC_aasector(X(I%), Y(I%), 16, 16, 0, 360, C%)
ENDIF
NEXT

CASE code% OF
WHEN 1: GCOL 1 : VDU 30,11 : PLOT 0, fault%, 0 : PRINT fault\$;
WHEN 2: GCOL 1 : VDU 30,11 : PLOT 0, range%, 0 : PRINT range\$;
ENDCASE

change% = FALSE
REPEAT
MOUSE X%, Y%, B%

IF B% IF B% <> oldB% THEN
near% = -1
dist% = &7FFFFFFF
oldsel% = selected%
FOR I% = 0 TO 3
D% = (X% - X(I%))^2 + (Y% - Y(I%))^2
IF D% < dist% dist% = D% : near% = I%
NEXT
IF dist% < 50*50 selected% = near% ELSE selected% = -1
IF selected% <> oldsel% change% = TRUE
oldX% = X% : oldY% = Y%
ENDIF

IF B% IF selected% >= 0 THEN
IF X% <> oldX% X(selected%) += X% - oldX% : change% = TRUE
IF Y% <> oldY% Y(selected%) += Y% - oldY% : change% = TRUE
ENDIF

oldX% = X% : oldY% = Y% : oldB% = B%

CASE INKEY(1) OF
WHEN 9:   selected% = (selected% + 3) MOD 4   : change% = TRUE
WHEN 155: selected% = (selected% + 1) MOD 4   : change% = TRUE
WHEN 136: IF selected% >= 0 X(selected%) -= 4 : change% = TRUE
WHEN 137: IF selected% >= 0 X(selected%) += 4 : change% = TRUE
WHEN 138: IF selected% >= 0 Y(selected%) -= 4 : change% = TRUE
WHEN 139: IF selected% >= 0 Y(selected%) += 4 : change% = TRUE
OTHERWISE:
ENDCASE

IF BB4W% SYS "InvalidateRect", @hwnd%, FALSE
*REFRESH

UNTIL change%

UNTIL FALSE
END

REM Calculate and draw an axis-aligned ellipse from four points:
DEF FNellipse4p(x(), y())
LOCAL a, b, c, x, y, pen%, lhs(), rhs(), pqu()
DIM lhs(2,2), rhs(2), pqu(2)
IF DIM(x(),1) <> 3 OR DIM(y(),1) <> 3 THEN = 3

lhs() = y(0)^2-y(1)^2, x(0)^2-x(1)^2, -2*(x(0)-x(1)), \
\       y(0)^2-y(2)^2, x(0)^2-x(2)^2, -2*(x(0)-x(2)), \
\       y(0)^2-y(3)^2, x(0)^2-x(3)^2, -2*(x(0)-x(3))

rhs() = 2*(y(0)-y(1)), 2*(y(0)-y(2)), 2*(y(0)-y(3))

PROC_invert(lhs())

pqu() = lhs() . rhs()

REM Find centre of ellipse:
x = pqu(2) / pqu(1)
y = 1 / pqu(0)

x() -= x
y() -= y

ON ERROR LOCAL x() += x : y() += y : = 1
c = (x(0)^2 * y(3)^2 - x(3)^2 * y(0)^2)
a = SQR(c / (y(3)^2 - y(0)^2))
b = SQR(c / (x(0)^2 - x(3)^2))

x() += x
y() += y

IF a > 32767 OR b > 1000 THEN = 2
IF BB4W% THEN
pen% = FN_gdipcreatepen(&FFFF0000, 0, 6)
PROC_gdiparc(pen%, x, y, a, b, 0, 360)
PROC_gdipdeletepen(pen%)
ELSE
PROC_aaarc(x, y, a, b, 0, 360, 6, &FF0000FF, 0)
ENDIF
= 0

DEF PROCcleanup
ON ERROR OFF
*REFRESH ON
IF BB4W% PROC_gdipexit
ENDPROC
``````
You can drag the points with the mouse (or your finger/stylus on a touchscreen) or you can use the keyboard: Tab and Shift+Tab select the required point and the four arrow keys move it. Four points are sufficient to define an axis-aligned ellipse uniquely, but there may be no ellipse that passes through them all (they may then define a hyperbola or parabola).

If you have a comment about the style or tone of this message please report it to the moderators by clicking the exclamation mark icon, rather than complaining on the public forum.