Xref: utzoo comp.graphics:4526 sci.math:5733 comp.sources.wanted:6417
Path: utzoo!utgpu!utstat!jarvis.csri.toronto.edu!mailrus!tut.cis.ohio-state.edu!ucbvax!decwrl!sun!pitstop!sundc!seismo!uunet!brunix!jfh
From: jfh@brunix (John Forbes Hughes)
Newsgroups: comp.graphics,sci.math,comp.sources.wanted
Subject: Re: looking for a fast ellipse algorithm
Keywords: ellipse, graphics, circle
Message-ID: <395@brunix.UUCP>
Date: 17 Feb 89 18:12:30 GMT
References: <399@peritek.UUCP> <1311@ndmath.UUCP>
Sender: news@brunix.UUCP
Reply-To: jfh@hope.UUCP (John Forbes Hughes)
Organization: Brown University Department of Computer Science
Lines: 95

In article <1311@ndmath.UUCP> krueger@ndmath.UUCP (Andreas Krueger) writes:
>In article <399@peritek.UUCP>, dig@peritek.UUCP (David Gotwisner) writes:
>> 
>> 	I am looking for an algorithm (or code) which will allow circle
>> generation onto a graphics device which drives a monitor with non-square
>> pixels.
>
>Drawing a circle on a device with non-square pixels is
>equivalent to drawing an ellipse on a monitor with square pixels.
> [Analytic geometry deleted]
>I've just been teaching this stuff in a third semester calculus course 
> [...]
>Now let x run from 0 to a^2/sqrt(a^2+b^2) in steps of 1 (pixel).
>For each such x, compute y by the above formula, which
>when solved for y yields y = (b/a)*sqrt(a^2 - x^2).
> [More deletion]
>From this, you can write the code yourself, I'd say. [...] If it 
>turns out that the arthmetic is too slow, you'll have to increase
>your steps from 1 to some larger number and draw connecting lines ...

Well, I've just been writing a book (Foley, van Dam, Feiner and Hughes)
about this, and the method above is not too clever. It uses floating point
arithmetic, sqrt, and doesn't handle the case of non-integer center/radius
properly.

Here's a different solution:

Assumptions: Your ellipse is aligned vertically (i.e. not rotated), and
   the center (h, k) is given in rational numbers with denominator d.
   Also, the eccentricity is small (since you say it's for a raster device,
   I assume the pixel aspect ratio (length/width) is near one). I also assume
   that the major and minor radii (a and b) are rational with denominator d.

(1) Write eqn for the ellipse: 

      (b(x - h))^2 + (a(y - k))^2 - (abR)^2 = S(x, y) = 0

(2) Multiply through by 4d^2 and expand out to a form like this

	S(x, y) = Ax^2 + Cy^2 + Dx + Ey + F = 0;

(3) Observe that A...E are all integers divisible by 4. F is not, so round
it off. Important observation: if S(x, y) > 0 at an integer (or half-integer)
point (x, y), it will still be >0 after rounding F, since rounding changes 
the value by less than one. Henceforth F refers to the rounded value of F.

(4) Start at the pixel (x, y) = (ROUND(h), ROUND(k + b)), draw this pixel.
Compute S(x+1, y - 1/2). If this is > 0, then move from (x, y) to (x+1, y-1)
and draw. If it is < 0, then move from (x, y) to (x + 1, y) and draw. Continue
this process until you reach the point where 2Ax + D > 2Cy + E. You have drawn
one eighth of your ellipse. The other 7 octants are a (very tedious) exercise.

IMPORTANT NOTE: recomputing S(x + 1, y - 1/2) for each (x, y) is expensive!
You can save a lot of time as follows: Note that you currently know either
S(x, y - 1/2) or S(x, y+1/2). TO compute S(x + 1, y - 1/2), in the first case
you just add dE = 2Ax + D; in the second case you add dSE = 2Ax + D + 2Cy - E;
A rough approximation of the correct loop now looks like this:

INITIALIZE EVERYTHING IN SIGHT, and then...
while (gradX < gradY)
   plot(x, y)
   dE = dE + 2A               	/* compute new value of dE from previous one */
   gradX = gradX + 2		/* compute new gradient values		*/
   if (last_move == E)		/* Same for dSE				*/
      dSE = dSE + 2A
   else
      dSE = dSE + 2A - 2C
      gradY = gradY - 2		/* moved SE, so gradY changes, too	*/
   
   if (last_move == E)
      S = S + dE
   else
      S = S + dSE
   
   x = x + 1			/* move to correct new point 	*/
   if (S > 0)
      y = y-1
      last_move = SE		/* moved to southeast		*/
   else
      last_move = E		/* moved east -- no change in y	*/
end-while

Note that for each pixel, we do about 3 additions (the numbers 2A, 2C, 
etc should all be precomputed) and a few compares. This should run about a
billion :-) times faster than any routine calling for square roots.

Also note that if your center is at an integer point, then you can actually
draw only 1/4 of the ellipse and use reflections to draw the other points.If
it's not at an integer center, this doesn't work.

All of this is a fairly trivial extension of the Midpoint Algorithm for
scan-converting circles. Foley and van Dam (1st edition) give a description
of another scan converter for circles which can be adapted in the same way.
Now, would you like antialiasing too???
-John Hughes, "the 4th author"