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