dmt@mtunb.ATT.COM (Dave Tutelman) (02/19/89)
Since I received a number of requests, here's a Turbo Pascal program for drawing fast ellipses. (Or drawing ellipses fast :-) A few notes about the program: - I've included a brief description of the math involved, as a lead comment. - While the support stuff (like the pixel function and the clear-screen-set-mode function) are IBM-PC-specific, the ellipse routine itself is perfectly generic Pascal, and easily converted to C. - I made a conscious decision to do it differently from the algorithm posted by Richard Stevens <301@hsi86.hsi.UUCP>, and attributed to J.R. van Aken. I use two multiplications per iteration, but fewer additions and subtractions. Bad trade? I doubt it, since I can do integer arithmetic rather than long arithmetic. That makes it faster on a machine where the penalty for long arithmetic is severe. If the long arithmetic is OK, but multiplication is a dog (not my machine), then van Aken's algorithm is better. +---------------------------------------------------------------+ | Dave Tutelman | | Physical - AT&T Bell Labs - Lincroft, NJ | | Logical - ...att!mtunb!dmt | | Audible - (201) 576 2442 | +---------------------------------------------------------------+ { ELLIPSE - a procedure to draw ellipses, and a demo program. } { Dave Tutelman - Feb 1989 } (******************** THEORY *************************************** We're using Bresenham's algorithm, best known for fast drawing of straight lines. However, it's not hard to adapt for other conic sections. We'll start with the canonical equation for an ellipse, in terms of its half-axes Xh and Yh: x^2 y^2 ---- + ---- = 1 (1) Xh^2 Yh^2 Bresenham's algorithm is a "digital differential analyzer", meaning that it works from the equation for the SLOPE. Let's differentiate (1) for slope: x y dy ---- + ---- -- = 0 (2) Xh^2 Yh^2 dx or, in two "more interesting" forms: (x Yh^2) dx + (y Xh^2) dy = 0 (3) dy x Yh^2 -- = - --- ---- (4) dx y Xh^2 Now, we'll start at <x=0, y=Yh>, and draw the <+,+> quadrant of the ellipse. Each time we find a point, we plot it in all four quadrants, using symmetry. The way the algorithm works is we "step" along, stepping x if the absolute slope is less than 1, and stepping y if the absolute slope is greater than 1. We use (4) to tell us whether to step x or y. If (x * Yh^2) > (y * Xh^2) then step y, otherwise step x. In this quadrant, stepping means we increment x and decrement y. If we didn't know that from visualizing the shape of an ellipse, we could tell by looking at (4); for positive x and y, the slope is negative. At each step, we check an "error" term that tells us if we're above or below the curve. On the basis of the error, we: - Decide whether to go straight (horizontal when stepping x, vertical when stepping y) or diagonal (change BOTH x and y). - Update the error term using (3), where dx is 0 or 1, dy is 0 or -1. If incrementing x, add (x * Yh^2) to error. If decrementing y, sub (y * Xh^2) from error. If moving diagonally, do both. The "ellipse" function below is a straight-forward coding of this algorithm, EXCEPT for our selection of Xh^2 and Yh^2. The equations show that, once we start from the right point, the only thing that matters is the RATIO of these two factors, not their absolute value. However, the magnitude of (x * Yh^2) or (y * Xh^2) can easily exceed the maximum value of an integer in a 16-bit machine. Rather than slow down the operation by doing our iterations in long arithmetic, we SCALE Xh^2 and Yh^2 so that all our products remain integer. We lose a bit of accuracy this way, but not much. (If you need accuracy, it might pay to "draw" it twice; once without putting pixels, and the second with a "corrected" value of Xh. That's probably faster than long arithmetic for machines without an arithmetic accelerator chip.) ********************************************************************) (******************** Pixel Primitives ****************************** * * These elementary routines are necessarily machine dependent. * pixel (x,y, val) put a pixel of color "val" at <x,y> * vidmode (mode) clear & switch the mode of the display. * * In the example below, we've coded for the CGA display of a DOS PC. * Re-do for your computer: *) uses DOS, { need this for BIOS calls } CRT; { just for readKey } { put a bright or dark pixel at <x,y> } procedure pixel (x, y, val : integer); var rr : Registers; begin rr.ah := 12; (* write-dot function *) rr.al := val; rr.cx := x; rr.dx := y; intr ($10, rr); end; { change the video mode } procedure vidmode (mode : byte); var rr : Registers; begin rr.ax := mode; intr ($10, rr); end; (************************ ELLIPSE **************************************) { OK, folks. Here's what you came for. } procedure ellipse (x0, y0, (* center of the ellipse *) xh, yh, (* ellipse full-axis lengths *) val (* 0=dark, 1=bright *) : integer); const BIGINT=22000; (* Choose so that xyh &yxh stay within * the range of an integer. 22,000 * seems to work well on a CGA screen * at medium resolution. *) var xd,yd : integer; (* x & y wrt center of circle *) xh2, yh2 : integer; (* xh squared & yh squared *) maxx,maxy, xfactor,yfactor : longint; maxfactor : longint; (* ;-) *) xyh, yxh : integer; (* x*yh2 and y*xh2 *) err : integer; begin { Actually want half-axis, so div by 2 } xh := abs(xh) div 2; yh := abs(yh) div 2; { It can be tricky to choose xh2 & yh2 so everything stays integer. However, it's worth getting all the non-integer arithmetic OUT OF THE MAIN ITERATION LOOP, for speed. } maxx := xh; maxy := yh; maxx := maxx * maxx; maxy := maxy * maxy; xfactor := maxx * longint(yh) div BIGINT; yfactor := maxy * longint(xh) div BIGINT; if yfactor>xfactor then maxfactor := yfactor else maxfactor := xfactor; if maxfactor<1 then maxfactor := 1; xh2 := maxx div maxfactor; yh2 := maxy div maxfactor; { Set up the iteration variables and go. } xd := 0; yd := yh; err := -yd; xyh:= 0; yxh:= yd*xh2; pixel (x0, y0+yd, val); pixel (x0, y0-yd, val); while yxh > xyh do begin xyh := xd*yh2; yxh := yd*xh2; if err>0 then (* DIAGONAL *) begin err := err + xyh - yxh; xd:=xd+1; yd:=yd-1; end else (* STRAIGHT *) begin err := err + xyh; xd:=xd+1; end; pixel (x0+xd, y0+yd, val); pixel (x0-xd, y0+yd, val); pixel (x0+xd, y0-yd, val); pixel (x0-xd, y0-yd, val); end; while yd>0 do begin xyh := xd*yh2; yxh := yd*xh2; if err<0 then (* DIAGONAL *) begin err := err + xyh - yxh; xd:=xd+1; yd:=yd-1; end else (* STRAIGHT *) begin err := err - yxh; yd:=yd-1; end; pixel (x0+xd, y0+yd, val); pixel (x0-xd, y0+yd, val); pixel (x0+xd, y0-yd, val); pixel (x0-xd, y0-yd, val); end; end; (************************************************************************) { Demo program. Mode setting & screen size are PC- and CGA-specific } const GRAFMODE = 5; ALFAMODE = 3; xxx = 320; yyy = 200; (* x & y size of screen, in pixels *) var xc,yc : integer; (* center of screen *) i : integer; inchar : char; { Make x- and y-axes, with tickmarks. } procedure ticks; var i : integer; xx, yy : integer; x10, y10 : integer; begin xx := xxx-1; x10 := xx div 10; yy := yyy-1; y10 := yy div 10; (* x-axis tickmarks *) for i:= 0 to xx do pixel (i, yc, 1); for i:= 0 to x10 do pixel (10*i+9, yc-1, 1); for i:= 0 to x10 do pixel (10*i+9, yc+1, 1); (* y-axis tickmarks *) for i:= 0 to yy do pixel (xc, i, 1); for i:= 0 to y10 do pixel (xc-1, 10*i+9, 1); for i:= 0 to y10 do pixel (xc+1, 10*i+9, 1); end; { MAIN PROGRAM } begin { Set up the video screen for CGA med-res (320 x 200) } vidmode (GRAFMODE); xc := (xxx-1) div 2; yc := (yyy-1) div 2; { Draw a single ellipse in the upper-left portion of screen } ticks; ellipse (xc div 2, yc div 2, xc-20, yc-20, 1); inchar := readKey; (* wait for keystroke *) vidmode (GRAFMODE); (* clear screen *) { Draw a family of increasingly large ellipses } ticks; for i:=1 to (xc div 10) do ellipse (xc,yc,20*i,10*i, 1); inchar := readKey; vidmode (GRAFMODE); ticks; { Family of constant y-axis } for i:=1 to (xc div 10) do ellipse (xc,yc,20*i,yyy-20, 1); inchar := readKey; { Family of constant x-axis } for i:=1 to (yc div 10) do ellipse (xc,yc,xxx-20,20*i, 1); { wait for keystroke, then back to alpha mode } inchar := ReadKey; vidmode (ALFAMODE); end.