Code archives/Graphics/Bresenham like Circle and Ellipse functions

This code has been declared by its author to be Public Domain code.

Download source code

Bresenham like Circle and Ellipse functions by Andy_A2011
After looking at 'furbrain's Circle plotter and MCP asking the question: "How difficult would it be to modify the Circle code above to draw elipses quickly?"
( http://www.blitzmax.com/codearcs/codearcs.php?code=2809 )

During my research into these routines I found that adding scan line fill code was just a couple of lines, so both routines have the option of drawing regular and filled circles/ellipses.

One other thing of note is that the Blitz filled oval is faster than an unfilled oval (approx. 33% faster!)

P.S. There are other Bresenham like circle algos here in the code archives
http://www.blitzbasic.com/codearcs/codearcs.php?code=428
http://www.blitzmax.com/codearcs/codearcs.php?code=1451
But this is the first Bresenham like ellipse plotter

Here are the results:

Circle
;Source:
;	http://homepage.smc.edu/kennedy_john/papers.htm
;	http://homepage.smc.edu/kennedy_john/bcircle.pdf

;	A Fast Bresenham Type Algorithm For Drawing Circles
;	by 
;	John Kennedy

;	Blitzplus/Blitz 3D port by Andy_A


AppTitle "Bresenham circle"
Graphics 800,600,32,2
SetBuffer BackBuffer()

numIterations% = 100

LockBuffer GraphicsBuffer()

st = MilliSecs()
For i = 1 To numIterations
	PlotCircle(400, 300, 290, $FFE020, 1)
Next
et1# = MilliSecs()-st

st = MilliSecs()
For i = 1 To numIterations
	Oval 110, 10, 580, 580, 1
Next
et2# = MilliSecs()-st

UnlockBuffer

Text 5,	5,"Avg/Circle: "+(et1/Float(numIterations))+"ms"
Text 5,20,"  Avg/Oval: "+(et2/Float(numIterations))+"ms"

Flip
WaitMouse()
End

Function PlotCircle(CX%, CY%, R%, colr%, fill% = 0)
	Local X%, Y%
	Local XChange%, YChange%
	Local RadiusError%
	Color (colr And $FF0000) Shr 16, (colr And $FF00) Shr 8, colr And $FF
	X = R
	Y = 0
	XChange = 1 - (R Shl 1)
	YChange = 1
	RadiusError = 0

	While  X >= Y

		If fill <> 0 Then
			Line(CX-X, CY+Y, CX+X, CY+Y)		;used calc'd values to draw
			Line(CX-X, CY-Y, CX+X, CY-Y)		;scan lines from points
			Line(CX-Y, CY+X, CX+Y, CY+X)		;in opposite octants
			Line(CX-Y, CY-X, CX+Y, CY-X)
		Else
			WritePixelFast(CX+X, CY+Y, colr)	; {point in octant 1}
			WritePixelFast(CX-X, CY+Y, colr)	; {point in octant 4}
			WritePixelFast(CX-X, CY-Y, colr)	; {point in octant 5}
			WritePixelFast(CX+X, CY-Y, colr)	; {point in octant 8}
			WritePixelFast(CX+Y, CY+X, colr)	; {point in octant 2}
			WritePixelFast(CX-Y, CY+X, colr)	; {point in octant 3}
			WritePixelFast(CX-Y, CY-X, colr)	; {point in octant 6}
			WritePixelFast(CX+Y, CY-X, colr)	; {point in octant 7}
		End If

		Y = Y + 1
		RadiusError = RadiusError + YChange
		YChange = YChange + 2

		If (RadiusError Shl 1) + XChange > 0 Then
			X = X - 1
			RadiusError = RadiusError + XChange
			XChange = XChange + 2
		End If
	Wend
End Function



Ellipse:
;Source:
;	http://homepage.smc.edu/kennedy_john/papers.htm
;	http://homepage.smc.edu/kennedy_john/belipse.pdf

;	A Fast Bresenham Type Algorithm For Drawing Ellipses
;	by 
;	John Kennedy

;	Blitzplus/Blitz 3D port by Andy_A

AppTitle "Bresenham Ellipse"
Graphics 800,600,32,2
SetBuffer BackBuffer()


numRepeats% = 100

LockBuffer GraphicsBuffer()

st = MilliSecs()
For i = 1 To numRepeats
	Ellipse(400, 300, 390, 290,$FFE020, 0)
Next
et1# = MilliSecs()-st

st = MilliSecs()
For i = 1 To numRepeats
	Oval 10, 10, 780, 580, 1
Next
et2# = MilliSecs()-st

UnlockBuffer


Text 5,	5,"Avg/Ellipse: "+(et1/Float(numRepeats))+"ms"
Text 5,20,"   Avg/Oval: "+(et2/Float(numRepeats))+"ms"


Flip
WaitMouse()
End

Function Ellipse(CX%, CY%, XRadius%, YRadius%, colr%, fill%);
	Local X%, Y%
	Local XChange%, YChange%
	Local EllipseError%
	Local TwoASquare%, TwoBSquare%
	Local StoppingX%, StoppingY%
	Color (colr And $FF0000) Shr 16, (colr And $FF00) Shr 8, colr And $FF

	TwoASquare = (XRadius*XRadius) Shl 1
	TwoBSquare = (YRadius*YRadius) Shl 1
	X = XRadius
	Y = 0
	XChange = YRadius*YRadius*(1-(XRadius Shl 1))
	YChange = XRadius*XRadius
	EllipseError = 0
	StoppingX = TwoBSquare*XRadius
	StoppingY = 0

	While StoppingX >= StoppingY 				; do {1st set of points, y' > -1}

		If fill <> 0 Then
			Line(CX-X, CY+Y, CX+X, CY+Y)		; used calc'd points to draw scan
			Line(CX-X, CY-Y, CX+X, CY-Y)		; lines from opposite quadrants
		Else
			WritePixelFast(CX+X, CY+Y, colr)	; {point in quadrant 1}
			WritePixelFast(CX-X, CY+Y, colr)	; {point in quadrant 2}
			WritePixelFast(CX-X, CY-Y, colr)	; {point in quadrant 3}
			WritePixelFast(CX+X, CY-Y, colr)	; {point in quadrant 4}
		End If

		Y = Y + 1
		StoppingY = StoppingY + TwoASquare
		EllipseError = EllipseError + YChange
		YChange = YChange + TwoASquare

		If (EllipseError Shl 1) + XChange > 0 Then
			X = X - 1
			StoppingX = StoppingX - TwoBSquare
			EllipseError = EllipseError + XChange
			XChange = XChange + TwoBSquare
		End If
	Wend
	
	;{ 1st set of points is done; start the 2nd set of points }
	X = 0
	Y = YRadius
	XChange = YRadius*YRadius
	YChange = XRadius*XRadius*(1 - (YRadius Shl 1))
	EllipseError = 0
	StoppingX = 0
	StoppingY = TwoASquare*YRadius

	While StoppingX <= StoppingY 				;do {2nd set of points, y' < -1}

		If fill <> 0 Then					
			Line(CX-X, CY+Y, CX+X, CY+Y)		; used calc'd points to draw scan
			Line(CX-X, CY-Y, CX+X, CY-Y)		; lines from opposite quadrants
		Else
			WritePixelFast(CX+X, CY+Y, colr)	; {point in quadrant 1}
			WritePixelFast(CX-X, CY+Y, colr)	; {point in quadrant 2}
			WritePixelFast(CX-X, CY-Y, colr)	; {point in quadrant 3}
			WritePixelFast(CX+X, CY-Y, colr)	; {point in quadrant 4}
		End If

		X = X + 1
		StoppingX = StoppingX + TwoBSquare
		EllipseError = EllipseError + XChange
		XChange = XChange + TwoBSquare

		If (EllipseError Shl 1) + YChange > 0 Then
			Y = Y - 1
			StoppingY = StoppingY - TwoASquare
			EllipseError = EllipseError + YChange
			YChange = YChange + TwoASquare
		End If
	Wend
End Function

Comments

MCP2011
Excellent stuff. Thanks for spending your time on this! :)

BTW: With debug off I get 14ms to draw an elipse with your code vs 38ms with Blitz Oval. Nice.


Andy_A2011
No worries,

I was researching some other math's when I stumbled on John Kennedy's site.

I posted because it may help others (until processors are so fast that there isn't any difference in speed).


Jesse2011
not fast but useful I think.
Function Ellipse(CX%, CY%, XRadius%, YRadius%,angle:Float)'
	Local X%, Y%
	Local XChange%, YChange%
	Local EllipseError%
	Local TwoASquare%, TwoBSquare%
	Local StoppingX%, StoppingY%
	Local xdiameter% = XRadius*XRadius
	Local ydiameter% = YRadius*YRadius
	Local xc#,ys#,yc#,xs#
	TwoASquare = (xdiameter) Shl 1
	TwoBSquare = (ydiameter) Shl 1
	X = XRadius
	Y = 0
	XChange = ydiameter*(1-(XRadius Shl 1))
	YChange = xdiameter
	EllipseError = 0
	StoppingX = TwoBSquare*XRadius
	StoppingY = 0
	Local c:Float=Cos(angle)
	Local s:Float=Sin(angle)
	While StoppingX >= StoppingY 				
		xc = c*x
		ys = s*y
		yc = c*y
		xs = s*x
		Plot CX+xc-ys, CY+yc+xs	
		Plot CX-xc-ys, CY+yc-xs	
		Plot CX-xc+ys, CY-yc-xs	
		Plot CX+xc+ys, CY-yc+xs	
		Y = Y + 1
		StoppingY = StoppingY + TwoASquare
		EllipseError = EllipseError + YChange
		YChange = YChange + TwoASquare

		If (EllipseError Shl 1) + XChange > 0 Then
			X = X - 1
			StoppingX = StoppingX - TwoBSquare
			EllipseError = EllipseError + XChange
			XChange = XChange + TwoBSquare
		End If
	Wend
	
	
	X = 0
	Y = YRadius
	XChange = YDiameter
	YChange = XDiameter*(1 - (YRadius Shl 1))
	EllipseError = 0
	StoppingX = 0
	StoppingY = TwoASquare*YRadius
	While StoppingX <= StoppingY 				
		xc = c*x
		ys = s*y
		yc = c*y
		xs = s*x
		Plot CX+xc-ys, CY+yc+xs	
		Plot CX-xc-ys, CY+yc-xs	
		Plot CX-xc+ys, CY-yc-xs	
		Plot CX+xc+ys, CY-yc+xs	
		X = X + 1
		StoppingX = StoppingX + TwoBSquare
		EllipseError = EllipseError + XChange
		XChange = XChange + TwoBSquare

		If (EllipseError Shl 1) + YChange > 0 Then
			Y = Y - 1
			StoppingY = StoppingY - TwoASquare
			EllipseError = EllipseError + YChange
			YChange = YChange + TwoASquare
		End If
	Wend
End Function



Andy_A2011
That's useful if the ovals are small and gaps between pixels is acceptable.

On my machine the WritePixelFast version is about 70x faster!
Function Ellipse(CX%, CY%, XRadius%, YRadius%, colr%, angle#)
	Local X%, Y%
	Local XChange%, YChange%
	Local EllipseError%
	Local TwoASquare%, TwoBSquare%
	Local StoppingX%, StoppingY%
	Local xdiameter% = XRadius*XRadius
	Local ydiameter% = YRadius*YRadius
	Local xc#,ys#,yc#,xs#, ca#, sa#
	
	TwoASquare = (xdiameter) Shl 1
	TwoBSquare = (ydiameter) Shl 1
	X = XRadius
	Y = 0
	XChange = ydiameter*(1-(XRadius Shl 1))
	YChange = xdiameter
	EllipseError = 0
	StoppingX = TwoBSquare*XRadius
	StoppingY = 0
	ca = Cos(angle)
	sa = Sin(angle)
	While StoppingX >= StoppingY 				
		xc = ca*x
		ys = sa*y
		yc = ca*y
		xs = sa*x
		WritePixelFast(CX+xc-ys, CY+yc+xs, colr)	
		WritePixelFast(CX-xc-ys, CY+yc-xs, colr)	
		WritePixelFast(CX-xc+ys, CY-yc-xs, colr)	
		WritePixelFast(CX+xc+ys, CY-yc+xs, colr)	
		Y = Y + 1
		StoppingY = StoppingY + TwoASquare
		EllipseError = EllipseError + YChange
		YChange = YChange + TwoASquare

		If (EllipseError Shl 1) + XChange > 0 Then
			X = X - 1
			StoppingX = StoppingX - TwoBSquare
			EllipseError = EllipseError + XChange
			XChange = XChange + TwoBSquare
		End If
	Wend
	
	
	X = 0
	Y = YRadius
	XChange = YDiameter
	YChange = XDiameter*(1 - (YRadius Shl 1))
	EllipseError = 0
	StoppingX = 0
	StoppingY = TwoASquare*YRadius
	While StoppingX <= StoppingY 				
		xc = ca*x
		ys = sa*y
		yc = ca*y
		xs = sa*x
		WritePixelFast(CX+xc-ys, CY+yc+xs, colr)
		WritePixelFast(CX-xc-ys, CY+yc-xs, colr)	
		WritePixelFast(CX-xc+ys, CY-yc-xs, colr)	
		WritePixelFast(CX+xc+ys, CY-yc+xs, colr)	
		X = X + 1
		StoppingX = StoppingX + TwoBSquare
		EllipseError = EllipseError + XChange
		XChange = XChange + TwoBSquare

		If (EllipseError Shl 1) + YChange > 0 Then
			Y = Y - 1
			StoppingY = StoppingY - TwoASquare
			EllipseError = EllipseError + YChange
			YChange = YChange + TwoASquare
		End If
	Wend
End Function



Jesse2011
I was just trying to illustrate rotation principles. I didn't try to make it exact or fast. if speed would have been my concern I would have given an example with Blitzmax and drawn to an image and simply hardware rotate it. it would have been faster than... and I didn't do it with writepixelfast because blitzbasic don't work on Macs and therefore couldn't test it.


Andy_A2011
I just wrote the B+ version so that those with B+/B3d can use a faster version. But since I didn't mention it in my prior post, thanks for making the modification.


Jimmy2011
Bresenham like code for rotated ellipses

Uses crunches huge numbers more than can be handled which introduces the jumpy bug at some sizes..

But it produces silky smooth ellipses, which is rather scarce to see these days.

Graphics  800,600,16,2:SetBuffer BackBuffer()
Repeat:Cls:xx=MouseX():yy=MouseY()
Color 0,255,0:ellipse(512,384,1+yy/2.0,200,xx/200.0)
Color 255,255,255:Plot 512,384:Plot 512+100,384:Plot 512-100,384:Plot 512,384+200:Plot 512,384-200
Flip
Until MouseDown(2)

Function ellipse (iXC#,iYC#,A#,B#,angle#)
;
; General ellipse
;
; Known bugs:
; jumpy movement sometimes, as it uses huge numbers, too large To fit storage
; It has few lines with trigonometry, these could very possibly be replaced with integer table or 
; tables of degrees Or similiar. It wraps the angle into an acceptable range, and sets xflag accordingly.
;
R2D# = 180.0/Pi:D2R# = Pi/180.0
If angle# < Pi/2.0
xflag#=1
ElseIf angle#>Pi/2 And angle#<Pi
temp#=angle#-Pi/2.0:angle#=Pi/2.0-temp#:xflag#=-1
ElseIf angle#>=Pi And angle#<Pi*1.5
angle#=angle#-Pi:xflag#=1
ElseIf angle#>=Pi*1.5 And angle#<Pi*2.0
angle#=angle#-Pi:temp#=angle#-Pi/2.0:angle#=Pi/2.0-temp#:xflag#=-1
EndIf

; ...To output these 4 integer points. (and the said xflag), These are where major and minor axis of ellipse ends (determines shape as in aligned axis algorithms).
ixa#=Int(Cos(angle#*R2D#)*A#)
iya#=Int(Sin(angle#*R2D#)*A#)
ixb#=Int(Cos((angle#+(Pi/2.0))*R2D#)*B#)
iyb#=Int(Sin((angle#+(Pi/2.0))*R2D#)*B#)

; rest of code uses only variables:
; ixa,iya,ixb,iyb
; ixc, iyc, xflag
; these are as every one else from now on, integers

Plot ixc#-ixa#,iyc#-iya#
Plot ixc#-ixb#,iyc#-iyb#
Plot ixc#-ixb#,iyc#-iya#
Plot ixc#-ixa#,iyc#-iyb#

; From now on, integers only, it uses multiplication much, and needs very large integers (in large resolutions even more than 64bits)
; therefor single precision are used to alloud the space, but theres no fixed point Or floating point involved.

   ixa2#=ixa#*ixa#
   iya2#=iya#*iya#
   ixb2#=ixb#*ixb#
   iyb2#=iyb#*iyb#
   ixaya#=ixa#*iya#
   ixbyb#=ixb#*iyb#
   ila2#=ixa2#+iya2#
   ila4#=ila2#*ila2#
   ilb2#=ixb2#+iyb2#
   ilb4#=ilb2#*ilb2#
   ia#=ixa2#*ilb4#+ixb2#*ila4#
   ib#=ixaya#*ilb4#+ixbyb#*ila4#
   ic#=iya2#*ilb4#+iyb2#*ila4#
   id#=ila4#*ilb4#

   If iYA# <= iXA#
       ; Start AT (-xA,-yA) 
       iX# = -iXA#:iY# = -iYA#:iDx# = -(iB#*iXA#+iC#*iYA#):iDy# = iA#*iXA#+iB#*iYA# 

       ; Arc FROM (-xA,-yA) TO point (x0,y0) where dx/dy = 0 
       While iDx# <= 0
           Plot iXC#+iX#*xflag#,iYC#+iY#
           Plot iXC#-iX#*xflag#,iYC#-iY#
           iY#=iY#+1
           iSigma# = iA#*iX#*iX#+2*iB#*iX#*iY#+iC#*iY#*iY#-iD# 
           If iSigma# < 0 Then iDx#=iDx#-iB#:iDy#=iDy#+iA#:iX#=iX#-1
           iDx#=iDx# + iC# 
           iDy#=iDy# - iB# 
       Wend 

       ; Arc FROM (x0,y0) TO point (x1,y1) where dy/dx = 1 
       While iDx# <= iDy# 
           Plot iXC#+iX#*xflag#,iYC#+iY# 
           Plot iXC#-iX#*xflag#,iYC#-iY# 
           iY#=iY#+1
           iXp1# = iX#+1 
           iSigma# = iA#*iXp1#*iXp1#+2*iB#*iXp1#*iY#+iC#*iY#*iY#-iD# 
           If iSigma# >= 0 Then iDx#=iDx# + iB#:iDy#=iDy# - iA#:iX# = iXp1# 
           iDx#=iDx# + iC# 
           iDy#=iDy# - iB# 
       Wend 

       ; Arc FROM (x1,y1) TO point (x2,y2) where dy/dx = 0 
       While iDy# >= 0
           Plot iXC#+iX#*xflag#,iYC#+iY#
           Plot iXC#-iX#*xflag#,iYC#-iY#
           iX=iX+1 
           iSigma# = iA#*iX#*iX#+2*iB#*iX#*iY#+iC#*iY#*iY#-iD# 
           If iSigma# < 0 Then iDx#=iDx# + iC# : iDy#=iDy# - iB# : iY#=iY#+1 
           iDx#=iDx# + iB# 
           iDy#=iDy# - iA# 
       Wend 

       ; Arc FROM (x2,y2) TO point (x3,y3) where dy/dx = -1 
       While iDy# >= -iDx# 
           Plot iXC#+iX#*xflag#,iYC#+iY#
           Plot iXC#-iX#*xflag#,iYC#-iY#
           iX#=iX#+1 
           iYm1# = iY#-1 
           iSigma# = iA#*iX#*iX#+2*iB#*iX#*iYm1#+iC#*iYm1#*iYm1#-iD# 
           If iSigma# >= 0 Then iDx#=iDx# - iC#:iDy#=iDy# + iB#:iY# = iYm1# 
           iDx#=iDx# + iB# 
           iDy#=iDy# - iA# 
       Wend 

       ; Arc FROM (x3,y3) TO (xa,ya) 
       While iY# >= iYA# 
           Plot iXC#+iX#*xflag#,iYC#+iY#
           Plot iXC#-iX#*xflag#,iYC#-iY#
           iY#=iY#-1
           iSigma# = iA#*iX#*iX#+2*iB#*iX#*iY#+iC#*iY#*iY#-iD# 
           If iSigma# < 0 Then iDx#=iDx# + iB#:iDy#=iDy# - iA#:iX#=iX#+1 
           iDx#=iDx# - iC# 
           iDy#=iDy# + iB# 
       Wend 

   Else 
       ; Start AT (-xa,-ya) 
       iX# = -iXA#:iY# = -iYA#:iDx# = -(iB#*iXA#+iC#*iYA#):iDy# = iA#*iXA#+iB#*iYA#

       ; Arc FROM (-xa,-ya) TO point (x0,y0) where dy/dx = -1 
       While -iDx# >= iDy# 
           Plot iXC#+iX#*xflag#,iYC#+iY#
           Plot iXC#-iX#*xflag#,iYC#-iY# 
           iX#=iX#-1
           iYp1# = iY#+1 
           iSigma# = iA#*iX#*iX#+2*iB#*iX#*iYp1#+iC#*iYp1#*iYp1#-iD# 
           If iSigma# >= 0 Then iDx#=iDx# + iC# : iDy#=iDy# - iB#:iY# = iYp1# 
           iDx#=IDx# - iB# 
           iDy#=Idy# +iA# 
       Wend 

       ; Arc FROM (x0,y0) TO point (x1,y1) where dx/dy = 0 
       While iDx# <= 0
           Plot iXC#+iX#*xflag#,iYC#+iY#
           Plot iXC#-iX#*xflag#,iYC#-iY#
           iY#=iY#+1
           iSigma# = iA#*iX#*iX#+2*iB#*iX#*iY#+iC#*iY#*iY#-iD# 
           If iSigma# < 0 Then iDx#=iDx# - iB#:iDy#=iDy# + iA#:iX#=iX#-1
           iDx#=iDx# + iC# 
           iDy#=iDy# - iB# 
       Wend 

       ; Arc FROM (x1,y1) TO point (x2,y2) where dy/dx = 1 
       While iDx# <= iDy# 
           Plot iXC#+iX#*xflag#,iYC#+iY#
           Plot iXC#-iX#*xflag#,iYC#-iY#
           iY#=iY#+1
           iXp1# = iX#+1 
           iSigma# = iA#*iXp1#*iXp1#+2*iB#*iXp1#*iY#+iC#*iY#*iY#-iD# 
           If iSigma# >= 0 Then iDx#=IDx# + iB# :iDy#=Idy# - iA# : iX# = iXp1# 
           iDx#=iDx# + iC# 
           iDy#=iDy# - iB# 
       Wend 

       ; Arc FROM (x2,y2) TO point (x3,y3) where dy/dx = 0 
       While iDy# >= 0
           Plot iXC#+iX#*xflag#,iYC#+iY#
           Plot iXC#-iX#*xflag#,iYC#-iY# 
           iX#=iX#+1 
           iSigma# = iA#*iX#*iX#+2*iB#*iX#*iY#+iC#*iY#*iY#-iD# 
           If iSigma# < 0 Then iDx#=iDx# + iC# : iDy#=iDy# - iB# : iY#=iY#+1 
           iDx#=iDx# + iB# 
           iDy#=iDy# - iA# 
       Wend 

       ; Arc FROM (x3,y3) TO (xa,ya) 
       While iX# <= iXA#
           Plot iXC#+iX#*xflag#,iYC#+iY# 
           Plot iXC#-iX#*xflag#,iYC#-iY# 
           iX#=iX#+1 
           iYm1# = iY#-1 
           iSigma# = iA#*iX#*iX#+2*iB#*iX#*iYm1#+iC#*iYm1#*iYm1#-iD# 
           If iSigma# >= 0 Then iDx#=iDx# - iC# : iDy#=iDy# + iB# :iY# = iYm1# 
           iDx#=iDx# + iB# 
           iDy#=iDy# - iA# 
       Wend 

   EndIf
End Function



MCP2011
Fantastic stuff Jimmy :)


Andy_A2011
Jimmy that's Brilliant!

I tested to see what a difference it would make to lock the buffer and use WritePixelFast, the results were over 50 times the original speed!

It approaches Bresenham's speed, with rotated ellipses no less.

Here's the WritePixelFast code I used to test:


I used this code at the top of the listing to modify your original post in making the speed comparisons.



Code Archives Forum