Article 71337 of comp.sys.cbm:
Path: news.acns.nwu.edu!newsfeed.acns.nwu.edu!news.luc.edu!chi-news.cic.net!ddsw1!news.mcs.net!hammer.uoregon.edu!infeed1.internetmci.com!newsfeed.internetmci.com!europa.clark.net!dispatch.news.demon.net!demon!delos.dra.hmg.gb!server1.netnews.ja.net!lyra.csx.cam.ac.uk!news.ox.ac.uk!news
From: imc@ecs.ox.ac.uk (Ian Collier)
Newsgroups: comp.sys.cbm,comp.sys.sinclair
Subject: Re: Elite line routine (was Re: spectrum faster in 3D)
Date: 11 Jul 1997 17:57:41 GMT
Organization: Oxford University Computing Laboratory, UK
Message-ID: <11449.imc@comlab.ox.ac.uk>
References: <33A9DD07.153E@skynet.be> <33C3830B.446B@wi.leidenuniv.nl> <33C4CD48.41C6@wi.leidenuniv.nl> <5q4jpr$4cu@news.acns.nwu.edu>
NNTP-Posting-Host: boothp2.ecs.ox.ac.uk
X-Local-Date: Friday, 11th July 1997 at 6:57pm BST
Lines: 118
Xref: news.acns.nwu.edu comp.sys.cbm:71337 comp.sys.sinclair:41239
In article <5q4jpr$4cu@news.acns.nwu.edu>, sjudd@nwu.edu (Stephen Judd) wrote:
>Now, I have a question for anyone interested: What is the Bresenham
>algorithm?
>I always assumed it was the strategy I use (along with just about everyone
>who sits down to think up an algorithm), but I keep seeing these 2DX-DY
>kinds of things that really mystify me. What's the thinking behind the
>algorithm?
They are, of course, more or less equivalent (it's down to which way
you go when the answer is exactly zero) and so I guess they can both
be referred to as Bresenham's algorithm.
The derivation is something like this (actually I just went to the library
so at least one text gives this derivation (-: ).
| ___ B
* T ___/
P |___/
___/
___/ * S
___/ |
A x=x0
You are drawing an approximation to the mathematical straight line AB (which
appears here as a wavy line, but that wasn't intentional!), and you have
reached a particular x co-ordinate, x0. The vertical line shown is the
line x=x0 and there are two adjacent pixels on it, called S and T. Your
job is to choose whichever one of S and T is closest to the intersection
between AB and the line x=x0, labelled P. The idea is to maintain a
variable d which is proportional to the value of PS-PT. Then clearly
you should choose point S if d<0 and T if d>0. If d=0 then either point
will do.
Initially, you have plotted point A and you are considering which pixel to
choose on the column immediately to the right of A. Then PT is 1-dy/dx and
PS is dy/dx so d is some constant times 2dy/dx-1. We may as well choose
that constant to be dx, so d is initially 2dy-dx.
If you move east, the new values of PS and PT are PS+dy/dx and PT-dy/dx
respectively, so the new value of d should be d+2dy.
If you move north, the new values of PS and PT are PS-1 and PT+1
respectively, so the new value of d would be d-2dx. So for a northeast
move therefore, the new value of d should be d+2(dy-dx).
This gives the Bresenham algorithm which Alvin quoted. We can of course
now make optimisations such as dividing the whole thing by two (we can show
formally that the lost half that might result from dividing an odd number by
two makes no difference provided it is rounded the correct way) and moving
the additions around slightly, to give our favourite algorithm.
>(I've also heard of Bresenham circle algorithms, but have never seen
>one of those either).
It can be derived in a manner rather similar to the above reasoning, but
instead here is my take on the matter. This may not be exactly Bresenham's
algorithm but it's similar in style and totally amazing in its simplicity if
you have never seen one before.
The first thing to notice is that you only need plot an eighth of the circle
since the rest can be done with mirrors. I will arbitrarily choose the arc
which rises from the east to the northeast (apparently Bresenham did the
whole 360 degrees, but anyway...). I will also simplify by assuming that the
circle is centred on the origin. The radius will be R.
Let e be the value x*x + y*y - R*R. For each value of y in the octant
that we are interested in, the perfect circle has an x such that e=0.
Call this x0. We would like to choose our x to be at most half a pixel
away from x0. So:
x - 1/2 <= x0 < x + 1/2
(square each term, noting that in our octant x-1/2 is positive for R>=1)
<=> x*x - x + 1/4 <= x0*x0 < x*x + x + 1/4
(add y*y-R*R to each term)
<=> e - x + 1/4 <= 0 < e + x + 1/4
<=> e - x + 1/4 <= 0 and 0 < e + x + 1/4
<=> e < x and -x <= e
<=> -x <= e < x
Let E be the value e+x. Then we want to choose a pixel such that
0 <= E < 2x.
Initially we have that x=R, y=0 and E=R. When we move north, the new
value of E would be:
x*x + (y+1)*(y+1) - R*R + x
= x*x + y*y - R*R - 2y - 1 + x
= E + 2y + 1.
If we then have to move west, the new value of E would be:
(x-1)*(x-1) + y*y - R*R + (x-1)
= x*x + y*y - R*R - x
= E - 2x.
Therefore we have to decide whether or not to move west depending on which
of E and E-2x is between 0 and 2x. In other words, we move west if E >= 2x.
So we have the following algorithm.
x=R; y=0; E=R;
while (y<=x) {
plot8(x,y);
E=E+2y+1; y=y+1;
if (E>=2x) {
E=E-2x; x=x-1;
}
}
If we let f = 2x-1-E then the check E>=2x becomes f<0 which is slightly
easier.
Interestingly I get a slightly different program each time I try this,
but I think the above one is correct this time. :-)
--
---- Ian Collier : imc@comlab.ox.ac.uk : WWW page (including Spectrum section):
------ http://www.comlab.ox.ac.uk/oucl/users/ian.collier/imc.html
Article 71378 of comp.sys.cbm:
Path: news.acns.nwu.edu!merle!judd
From: judd@merle.acns.nwu.edu (Stephen Judd)
Newsgroups: comp.sys.cbm,comp.sys.sinclair
Subject: Re: Elite line routine (was Re: spectrum faster in 3D)
Date: 12 Jul 1997 22:05:48 GMT
Organization: Northwestern University, Evanston, IL
Lines: 83
Message-ID: <5q8v3s$m38@news.acns.nwu.edu>
References: <33A9DD07.153E@skynet.be> <33C4CD48.41C6@wi.leidenuniv.nl> <5q4jpr$4cu@news.acns.nwu.edu> <11449.imc@comlab.ox.ac.uk>
Reply-To: sjudd@nwu.edu (Stephen Judd)
NNTP-Posting-Host: merle.acns.nwu.edu
Xref: news.acns.nwu.edu comp.sys.cbm:71378 comp.sys.sinclair:41299
In article <11449.imc@comlab.ox.ac.uk>, Ian Collier wrote:
>In article <5q4jpr$4cu@news.acns.nwu.edu>, sjudd@nwu.edu (Stephen Judd) wrote:
>>Now, I have a question for anyone interested: What is the Bresenham
>>algorithm?
>
>They are, of course, more or less equivalent (it's down to which way
>you go when the answer is exactly zero) and so I guess they can both
>be referred to as Bresenham's algorithm.
>
>The derivation is something like this (actually I just went to the library
>so at least one text gives this derivation (-: ).
Thanks -- I've gone through it, still find it confusing, and will hence go
through it again later :).
For what it's worth, my own derivation goes like:
From the definition of a line: dy=m*dx where m=slope. That is, if I
take a step of size dx in the x direction, I will take a corresponding
step dy in the y-direction.
Since we're on a computer, steps of interest are of size 1. Let's say we
move in the x-direction by one at each step:
x=x+1
y=y+m
where m=Dy/Dx i.e. (y2-y1)/(x2-x1). So, if we just keep a running sum
of m at each step, once that sum is >1 it's time to take a step in y.
The running sum just looks like
m + m + m + ... = Dy/Dx + Dy/Dx + Dy/Dx + ...
= (Dy + Dy + Dy + ...)/Dx
So, the idea is to just keep adding Dy to itself, and once it is larger
than Dx, the fraction is >1 and it is time to take a step in y. By then
subtracting off Dx, we retain the remainder. The only thing to realize
is that by starting this counter off at Dx/2 we are in effect rounding
the numbers up (we've just added 0.5 to the fraction), which in turn
has the effect of splitting one of the horizontal line segments evenly
between the first and the last point.
So, by my reckoning, the 2*dy-dx junk is all just a complicated way
of starting a counter at dx/2.
>>(I've also heard of Bresenham circle algorithms, but have never seen
>>one of those either).
>
>It can be derived in a manner rather similar to the above reasoning, but
>instead here is my take on the matter. This may not be exactly Bresenham's
>algorithm but it's similar in style and totally amazing in its simplicity if
>you have never seen one before.
Again, thanks. The circle routine I came up with was
Y = Radius
X = 0
A = Radius/2
:loop
Plot4(x,y) ;Could just use Plot8
X = X + 1
A = A - X
if A<0 then A=A+Y : Y=Y-1
if X < Y then :loop
As you say, you can use plot8, but in BLARG I actually computed the next
eighth by basically reversing the above steps, and using plot4. That way
I could maintain four pointers and make the plotting very fast. I think I
put a derivation in C=Hacking somewhere; what it does is solve the
differential equation for a circle
dy/dx = -x/y
using a running sum, similar to the line routine above. Actually when
I first wrote it I had a very strange geometric/pictoral reasoning -- the
above derivation came a day later. I actually have a secret belief that
there's something deeper going on, because the routine is just too beautiful,
and works too well across all values of R.
Anyways, it looks pretty similar to your routine, and it continues
to amaze me that it actually works, and works incredibly well.
-S
Article 71468 of comp.sys.cbm:
Path: news.acns.nwu.edu!newsfeed.acns.nwu.edu!news.ece.nwu.edu!news.cse.psu.edu!rutgers!cam-news-feed2.bbnplanet.com!cam-news-hub1.bbnplanet.com!cpk-news-hub1.bbnplanet.com!news.bbnplanet.com!europa.clark.net!dispatch.news.demon.net!demon!delos.dra.hmg.gb!server1.netnews.ja.net!server5.netnews.ja.net!server3.netnews.ja.net!server4.netnews.ja.net!server2.netnews.ja.net!news.ox.ac.uk!news
From: imc@ecs.ox.ac.uk (Ian Collier)
Newsgroups: comp.sys.cbm,comp.sys.sinclair
Subject: Re: Elite line routine (was Re: spectrum faster in 3D)
Date: 14 Jul 1997 17:48:56 GMT
Organization: Oxford University Computing Laboratory, UK
Message-ID: <11471.imc@comlab.ox.ac.uk>
References: <33A9DD07.153E@skynet.be> <5q4jpr$4cu@news.acns.nwu.edu> <11449.imc@comlab.ox.ac.uk> <5q8v3s$m38@news.acns.nwu.edu>
NNTP-Posting-Host: boothp2.ecs.ox.ac.uk
X-Local-Date: Monday, 14th July 1997 at 6:48pm BST
Lines: 77
Xref: news.acns.nwu.edu comp.sys.cbm:71468 comp.sys.sinclair:41435
In article <5q8v3s$m38@news.acns.nwu.edu>, sjudd@nwu.edu (Stephen Judd) wrote:
>For what it's worth, my own derivation goes like:
>From the definition of a line: dy=m*dx where m=slope. That is, if I
>take a step of size dx in the x direction, I will take a corresponding
>step dy in the y-direction.
[and so on]
I once (for the purposes of a computing exercise set by someone else for
which I had to give a class) derived our favourite algorithm from this
program:
VAR yt: REAL;
xd, yd, xi, yi: INTEGER;
FOR xi:=0 TO xd DO BEGIN
yt := (yd/xd)*xi;
yi := trunc(yt + 1/2);
write_pixel(xi, yi);
END;
(and don't ask me what language this is supposed to be in (-: ).
>So, by my reckoning, the 2*dy-dx junk is all just a complicated way
>of starting a counter at dx/2.
Well as I showed it does have a mathematical reason (which perhaps boils
down to starting a counter at dx/2).
>Again, thanks. The circle routine I came up with was
> Y = Radius
> X = 0
> A = Radius/2
>:loop
> Plot4(x,y) ;Could just use Plot8
> X = X + 1
> A = A - X
> if A<0 then A=A+Y : Y=Y-1
> if X < Y then :loop
We should perhaps test these, since I know mine is correct and yours isn't
the same. :-) Incidentally, the one I looked up at the same time as the
line algorithm had some more complicated formulae in it, which is why I
didn't pay any attention to it.
I used to use the old y=sqrt(R*R-x*x) but this is obviously rather quicker.
> I actually have a secret belief that
>there's something deeper going on, because the routine is just too beautiful,
>and works too well across all values of R.
Not really; after all the solution to dy/dx=-x/y does seem to be the right
thing to calculate. I did observe that mine was similar to a line routine
for dx=-y and dy=x.
>Anyways, it looks pretty similar to your routine, and it continues
>to amaze me that it actually works, and works incredibly well.
Mine is equivalent to this:
Y = Radius
X = 0
A = (Radius-1)/2
:loop
Plot8(x,y)
X = X + 1
A = A - X + 1/2
if A<0 then A=A+Y-1 : Y=Y-1
if X <= Y then :loop
so not so different. And I have given a proof for mine so it shouldn't be
surprising that it works. Yours loses 1/2 for each straight step and gains
1/2 for each diagonal step.
--
---- Ian Collier : imc@comlab.ox.ac.uk : WWW page (including Spectrum section):
------ http://www.comlab.ox.ac.uk/oucl/users/ian.collier/imc.html
Article 71466 of comp.sys.cbm:
Path: news.acns.nwu.edu!newsfeed.acns.nwu.edu!news.ece.nwu.edu!news.cse.psu.edu!uwm.edu!vixen.cso.uiuc.edu!howland.erols.net!europa.clark.net!dispatch.news.demon.net!demon!delos.dra.hmg.gb!server1.netnews.ja.net!server5.netnews.ja.net!server6.netnews.ja.net!server4.netnews.ja.net!server2.netnews.ja.net!news.ox.ac.uk!news
From: imc@ecs.ox.ac.uk (Ian Collier)
Newsgroups: comp.sys.cbm,comp.sys.sinclair
Subject: Re: Elite line routine (was Re: spectrum faster in 3D)
Date: 14 Jul 1997 17:54:08 GMT
Organization: Oxford University Computing Laboratory, UK
Message-ID: <11472.imc@comlab.ox.ac.uk>
References: <33A9DD07.153E@skynet.be> <33C4CD48.41C6@wi.leidenuniv.nl> <5q4jpr$4cu@news.acns.nwu.edu> <33C67EF8.3BA9@chebucto.ns.ca>
NNTP-Posting-Host: boothp2.ecs.ox.ac.uk
X-Local-Date: Monday, 14th July 1997 at 6:54pm BST
Lines: 30
Xref: news.acns.nwu.edu comp.sys.cbm:71466 comp.sys.sinclair:41419
In article <33C67EF8.3BA9@chebucto.ns.ca>, aa601@chebucto.ns.ca wrote:
>x=-1
>y=r+1
>loop:
>x=x+1
>r=r+x+1
>if r>=0 then y=y-1:r=r-2*y
>plot8(x,y,c)
>if y>x then goto loop:
This can't be right, for the simple reason that you are adding x and
subtracting 2y, so the approximate gradient of the curve will be -x/2y.
>Judd version:
>a=r/2
>x=r
>y=0
>loop:
>y=y+1
>a=a-y
>if a<=0 then x=x-1:a=a+x
^^^^^^^^^^^
>plot8(x,y)
>if y
References: <33A9DD07.153E@skynet.be> <11449.imc@comlab.ox.ac.uk> <5q8v3s$m38@news.acns.nwu.edu> <11471.imc@comlab.ox.ac.uk>
Reply-To: sjudd@nwu.edu (Stephen Judd)
NNTP-Posting-Host: merle.acns.nwu.edu
Xref: news.acns.nwu.edu comp.sys.cbm:71656 comp.sys.sinclair:41637
In article <11471.imc@comlab.ox.ac.uk>, Ian Collier wrote:
>In article <5q8v3s$m38@news.acns.nwu.edu>, sjudd@nwu.edu (Stephen Judd) wrote:
>
>> Y = Radius
>> X = 0
>> A = Radius/2
>>:loop
>> Plot4(x,y) ;Could just use Plot8
>> X = X + 1
>> A = A - X
>> if A<0 then A=A+Y : Y=Y-1
>> if X < Y then :loop
>
>We should perhaps test these, since I know mine is correct and yours isn't
(How's that for a misrepresented quotation :)
>the same. :-) Incidentally, the one I looked up at the same time as the
I know the above works, because it is what I used in BLARG. It passes
the all-important eye test at all radii R>0. I think I also checked it
against (r*cos(t),r*sin(t)), and it was pretty much dead-on.
Interestingly, it is very robust to mangling. For instance, you could
interchange the A=A+Y:Y=Y-1 and still get a decent looking circle
(at most radii, but a few of the smaller ones don't come out right).
>line algorithm had some more complicated formulae in it, which is why I
>didn't pay any attention to it.
Every circle routine I've seen (up until yours :) has always been quite
complicated.
>> I actually have a secret belief that
>>there's something deeper going on, because the routine is just too beautiful,
>>and works too well across all values of R.
>
>Not really; after all the solution to dy/dx=-x/y does seem to be the right
>thing to calculate. I did observe that mine was similar to a line routine
Well, if you buy that it is just a solution of the diffeq, then you are
accepting that a wholly trivial discrete addition routine solves a nonlinear,
nonautonomous differential equation for all initial conditions!
Besides, like I said, the diffeq story was just a cover :). Originally I
imagined stacking blocks.
Anyways, it's something interesting to lay awake at night thinking about.
-S
Article 71783 of comp.sys.cbm:
Path: news.acns.nwu.edu!newsfeed.acns.nwu.edu!news.luc.edu!chi-news.cic.net!infeed1.internetmci.com!newsfeed.internetmci.com!cpk-news-hub1.bbnplanet.com!news.bbnplanet.com!baron.netcom.net.uk!netcom.net.uk!server3.netnews.ja.net!news.ox.ac.uk!news
From: imc@ecs.ox.ac.uk (Ian Collier)
Newsgroups: comp.sys.cbm,comp.sys.sinclair
Subject: Re: Elite line routine (was Re: spectrum faster in 3D)
Date: 18 Jul 1997 16:19:46 GMT
Organization: Oxford University Computing Laboratory
Lines: 81
Message-ID: <11516.imc@comlab.ox.ac.uk>
References: <33A9DD07.153E@skynet.be> <33C4CD48.41C6@wi.leidenuniv.nl> <5q4jpr$4cu@news.acns.nwu.edu> <33C67EF8.3BA9@chebucto.ns.ca>
NNTP-Posting-Host: gruffle.comlab.ox.ac.uk
X-Local-Date: Friday, 18th July 1997 at 5:19pm BST
Xref: news.acns.nwu.edu comp.sys.cbm:71783 comp.sys.sinclair:41774
In article <33C67EF8.3BA9@chebucto.ns.ca>, aa601@chebucto.ns.ca wrote:
>I can't recall it right now. Here are the bresenham routines:
>(at least some variation)
>void circle(im,x0,y0,r,c)
> Image *im;
> int x0,y0,r,c;
>{
> register int x=-1,y=r+1;
> do
> { r += ((++x - ( (r >= 0) ? --y : 0) )<<1) + 1;
[pixel writes elided]
[ } ]
> while (y > x);
>}
>/* and semi translated to basic:
>x=-1
>y=r+1
>loop:
>x=x+1
>r=r+x+1
>if r>=0 then y=y-1:r=r-2*y
>plot8(x,y,c)
>if y>x then goto loop: */
OK, it seems that you lost a 2 somewhere. The exact translation is...
x=-1: y=r+1
loop:
if r>=0 then y=y-1: r=r-2*y
x=x+1: r=r+2*x+1
plot8(x,y,c)
if y>x then goto loop:
and it turns out to be equivalent to the program I posted earlier.
To find this out you need to move the if backwards past the loop,
x=-1: y=r: r=-r
loop:
x=x+1: r=r+2*x+1
plot8(x,y,c)
if r>=0 then y=y-1: r=r-2*y
if y>=x+1 then goto loop:
let the new x be the old x+1
x=0: y=r: r=-r
loop:
plot8(x,y,c)
x=x+1: r=r+2*x-1
if r>=0 then y=y-1: r=r-2*y
if y>=x then goto loop:
and let e be r+2y.
x=0: y=r: e=r
loop:
plot8(x,y,c)
e=e+2*x+1: x=x+1
if r>=2*y then e=e-2*y: y=y-1
if y>=x then goto loop:
Not that this interests you very much of course. :-)
>Judd version:
>a=r/2
>x=r
>y=0
>loop:
>y=y+1
>a=a-y
>if a<=0 then x=x-1:a=a+x
>plot8(x,y)
>if y