More pythonic circle?

[Pythor]

I wrote the following code for a personal project. I need a function
that will plot a filled circle in a two dimensional array. I found
Bresenham's algorithm, and produced this code. Please tell me there's
a better way to do this.

import numpy

def circle(field=None,radius,center=(0,0),value=255,):
'''Returns a list of points within 'radius' distance from point
'center'.'''
if field==None:
field=numpy.zeros((radius,radius),'u')
cx,cy=center
filllist=[]
dy,dx=0,radius
r2=radius**2
while dy<=(radius*.71): # sin of 45 degrees
if dx**2+dy**2<=r2:
for x in range(dx):
filllist.append((x,dy))
dy+=1
else:
dx-=1
if dx<dy : break
resultlist=[]
for (i,j) in filllist:
field[cx+i][cy+j]=value
field[cx+j][cy+i]=value
field[cx-j][cy+i]=value
field[cx-i][cy+j]=value
field[cx-i][cy-j]=value
field[cx-j][cy-i]=value
field[cx+j][cy-i]=value
field[cx+i][cy-j]=value
return field

 

[John Machin]

OTTOMH, in a rush to go out: never mind Pythonic, following apply to
any language:
(1) accuracy: (a) sue me if I'm wrong, but I think you need range(dx+1)
so that the dx pixel is filled in (b) a few more digits after 0.71
might be useful
(2) efficiency: seems that range(dy, dx+1) would save some wear & tear
on the circuitry
(3) legibility: there's no prize for the script with the absolutely
minimum number of space characters
I think I've got an article on better Bresenham somewhere in the
archives; will dig it out later.
Cheers,
John

 

[Pythor]

OTTOMH, in a rush to go out: never mind Pythonic, following apply to
any language:
(1) accuracy: (a) sue me if I'm wrong, but I think you need range(dx+1)
so that the dx pixel is filled in

Hmm. I think you're right. Thanks.

(b) a few more digits after 0.71
might be useful

Sine of 45 degrees is actually .707... I rounded up, since I was using
<=. Figured that would make it clear.

(2) efficiency: seems that range(dy, dx+1) would save some wear & tear
on the circuitry

It took me a few minutes to figure out waht you meant here. This will
certainly help reduce the repeated coordinates. Thanks, again.

(3) legibility: there's no prize for the script with the absolutely
minimum number of space characters

True. I assume your saying I should make cx,cy,dx, and dy better
names. I probably will. Up to now I was just playing around with
this, and not really expecting anyone else to read it.

I think I've got an article on better Bresenham somewhere in the
archives; will dig it out later.

I'd definitely appreciate it. In fact, I'm trying to find a decent
sphere version, and getting nowhere with google. I tried to figure it
out on my own, and ended up with 48 coordinates for each valid test.
I'm not sure if that's right.

Lee

 

[Felipe Almeida Lessa]

Em Sáb, 2006-04-08 às 21:17 -0700, Pythor escreveu:

True. I assume your saying I should make cx,cy,dx, and dy better
names. I probably will. Up to now I was just playing around with
this, and not really expecting anyone else to read it.

This is kinda funny because you just posted the code to a list with
hundreds of developers. =)

 

[Michael Tobis]

Proving yet again that it's possible to write Fortran in any language.

You aren't getting any benefit from numpy or python here. Are you
aiming for speed or legibility?

Also, with this code, you are using radius for the dimensions of the
enclosing box, as well as the radius of the circle, so it's guaranteed
to not to actually produce a whole circle. Recall what python does with
negative indices!

I'll bet this does the trick for you and runs faster than what you've
got

def circle(rad = 5,max_x = 20, max_y = 20,cx = 10, cy= 10, value=255):
radsq = rad * rad
return numpy.array([[((x - cx) ** 2 + (y - cy) ** 2 < radsq) and
value or 0 for x in range(max_x)] for y in range(max_y)],'u')

I think the pure numpy solution should be something like (untested)

def circle(rad = 5,max_x = 20, max_y = 20,cx = 10, cy= 10, value=255):
def sqdist(x,y):
return (x - cx) * (x - cx) + (y - cy) * (y - cy)
distarray = numpy.fromfunction(sqdist,(max_y,max_x))
return
numpy.asarray(numpy.choose(greater(distarray,rad*rad),(0,value),'u')

Neither approach will get you the eightfold speedup that the messy code
was aimed at, but in practice they will spend less time at the
interpreter level and will likely run faster.

mt

 

[Steven D'Aprano]

True. I assume your saying I should make cx,cy,dx, and dy better
names. I probably will. Up to now I was just playing around with
this, and not really expecting anyone else to read it.

No, "minimum number of space characters" means "you don't use enough
spaces", not "your variable names are too short" *wink*

Generally speaking, a little bit of whitespace helps makes text more
readable, at least for native speakers of languages derived from Latin and
other Indo-European languages. Graphic designers tend to manually adjust
the spacing between paragraphs, lines, words and even characters, but it
isn't necessary to go to that extreme to increase readability.

People's tastes vary, and these are not rules, only guidelines, but it is
usual to leave one (and sometimes more) blank lines between functions and
methods, and even between logical sections of a single function.

Within a single line, a good guideline is to leave a single space on
either side of pluses and minuses (e.g. x**2 + 5*x - 3). Likewise, a
single space on both sides of an equal sign and a single space after
commas tend to be usual.

As I said, none of these are hard and fast rules, but breaking up the flow
of tokens will increase readability immensely.

See Guido van Rossum's style guide and the official Python style guide. As
usual, Guido is always right, and if his style guide contradicts me, he's
wrong but you should follow his guide anyway *smiles*

http://www.python.org/doc/essays/styleguide.html
http://www.python.org/dev/peps/pep-0008/

As for variables cx, cy, dx and dy, I don't believe that they are unclear.
If your function is highly mathematical in nature, I believe it is
acceptable if not expected to follow standard mathematical conventions
such as r for radius, x and y for real numbers, z for complex, dx for
delta (change of) x, etc. If in doubt, when initialising the variable add
a comment spelling it out in full.

On the other hand, you do have an argument "value" with default 255, with
not even hint for what it is.

 

[Pythor]

Proving yet again that it's possible to write Fortran in any language.
Ouch...

You aren't getting any benefit from numpy or python here. Are you
aiming for speed or legibility?

Speed will be a necessity, eventually. I was just really aiming for
something that works, and that I am capable of writing.

Also, with this code, you are using radius for the dimensions of the
enclosing box, as well as the radius of the circle, so it's guaranteed
to not to actually produce a whole circle. Recall what python does with
negative indices!

I'm not sure what you mean here. It produces an eighth-circle, and
then plots each point in the 8 symmetrical positions on the circle.
Except for the (dx+1) point made above, what piece of the circle is
missing?

I'll bet this does the trick for you and runs faster than what you've
got

def circle(rad = 5,max_x = 20, max_y = 20,cx = 10, cy= 10, value=255):
radsq = rad * rad
return numpy.array([[((x - cx) ** 2 + (y - cy) ** 2 < radsq) and
value or 0 for x in range(max_x)] for y in range(max_y)],'u')

I think the pure numpy solution should be something like (untested)

def circle(rad = 5,max_x = 20, max_y = 20,cx = 10, cy= 10, value=255):
def sqdist(x,y):
return (x - cx) * (x - cx) + (y - cy) * (y - cy)
distarray = numpy.fromfunction(sqdist,(max_y,max_x))
return
numpy.asarray(numpy.choose(greater(distarray,rad*rad),(0,value),'u')

I'll take a look at both of these. At this point, I can't quite wrap
my head around what you're doing for either one.

 

[Pythor]

No, "minimum number of space characters" means "you don't use enough
spaces", not "your variable names are too short" *wink*

Hmm. Guess I can't read too well.

Within a single line, a good guideline is to leave a single space on
either side of pluses and minuses (e.g. x**2 + 5*x - 3). Likewise, a
single space on both sides of an equal sign and a single space after
commas tend to be usual.

What I produced was natural for my fingers, but I can see that it's
difficult on the eyes. I'll try to remember that.

As for variables cx, cy, dx and dy, I don't believe that they are unclear.
If your function is highly mathematical in nature, I believe it is
acceptable if not expected to follow standard mathematical conventions
such as r for radius, x and y for real numbers, z for complex, dx for
delta (change of) x, etc. If in doubt, when initialising the variable add
a comment spelling it out in full.

On the other hand, you do have an argument "value" with default 255, with
not even hint for what it is.

Well, value is just a value. It's the number that get's punched into
the array for any points within the circle. I didn't have any better
name I could think of.

 

[Fredrik Lundh]

Speed will be a necessity, eventually. I was just really aiming for
something that works, and that I am capable of writing.

any special reason you cannot use an existing graphics library ?

</F>

 

[Pythor]

any special reason you cannot use an existing graphics library ?

</F>

Well, I'm not really interested in pretty pictures, but in the
resulting array. It might be worth using a graphics library and then
converting from an image to an array when I'm done. I've been playing
with PIL for that. But I wanted to see what I could do on my own, too.
In addition, I'll eventually need some functions that are definitely
not in a standard graphics library, and I wanted to get started with
something I thought was relatively simple.

 

[John Machin]

It's also possible to write microprocessor assembly language in any
other language.
The following code generates the OP's list of points with nothing more
complicated than integer addition/subtraction inside the loop. It also
does the right thing if the radius is not an integer, and avoids the
OP's "0.71 approx == sin(45 deg) aka 1/sqrt(2)" caper.

Wrt your Pythonic suggestion: "I'll bet this does the trick for you and
runs faster than what you've got": You lose on "does the trick" (should
be <= radsq). As for the second clause, a prerequisite to testing that
is to get the OP to say what his typical radius and typical enclosing
box size are (and get your point that they should not be the same).

Cheers,
John

# def octant(radius):
# assert radius >= 0
# filllist = []
# dx = int(radius)
# dy = 0
# trigger = dx * dx - int(radius * radius)
# dx_squared_delta = dx + dx - 1
# dy_squared_delta = 1
# while dy <= dx:
# if trigger <= 0:
# for x in range(dy, dx+1):
# filllist.append((x, dy))
# dy += 1
# trigger += dy_squared_delta
# dy_squared_delta += 2
# else:
# dx -= 1
# trigger -= dx_squared_delta
# dx_squared_delta -= 2
# filllist.sort()
# print "%.2f %r" % (radius, filllist)

# if __name__ == "__main__":
# octant(3.99)
# octant(4)
# octant(4.01)
# octant(3.60)
# octant(3.61)
# octant(0.01)
# octant(0)

 

[John Machin]

[Michael Tobis]
Also, with this code, you are using radius for the dimensions of the
enclosing box, as well as the radius of the circle, so it's guaranteed
to not to actually produce a whole circle. Recall what python does with
negative indices!

[Pythor]
I'm not sure what you mean here. It produces an eighth-circle, and
then plots each point in the 8 symmetrical positions on the circle.
Except for the (dx+1) point made above, what piece of the circle is
missing?
======
What Michael means is that for example a circle of radius 5 is 11
pixels wide and 11 pixels high. You are trying to cram it into a box of
5 x 5 pixels [or maybe 6x6 (I'm numpy-challenged)]. The result will
resemble a train smash. Have you considered *TESTING* your code? It's
not very difficult at all to draw the expected results for radii of
about 4 or 5 pixels on the back of an envelope ...

By the way, there are signs of a benchmark war breaking out. What are
typical sizes you would be using in practice for the radius and the
enclosing box?

Cheers,
John

 

[Scott David Daniels]

I wrote the following code for a personal project. I need a function
that will plot a filled circle in a two dimensional array. I found
Bresenham's algorithm, and produced this code. Please tell me there's
a better way to do this.

import numpy

def circle(field=None,radius,center=(0,0),value=255,):

...

Break this code into two functions. Roughly:

def wedge_nz(radius):
'''points in octant 2 of 0-origin radius-sized circle, no zeroes'''
x, y = int(radius), 1
dr2 = x ** 2 + y ** 2
r2 = radius ** 2
dy_limit = int(radius * .71) # sin of 45 degrees + a smidge
while y <= dy_limit:
if r2 >= dr2: # dx**2 + dy**2
for tx in range(1, x + 1):
yield tx, y
dr2 += y * 2 + 1 # dr2 = x ** 2 + (y + 1) ** 2
y += 1
else:
dr2 -= x * 2 - 1 # dr2 = (x - 1) ** 2 + y ** 2
x -= 1
if x < y:
break

and:

def circle2(field=None, radius=3.2, center=(0, 0), entry=255):
'''Fill field at all points within 'radius' of center with entry.

center is assumed integral.
'''
x, y = center
if field is None:
field = numpy.zeros((int(x + radius), int(y + radius)), 'u')

for px, py in wedge_nz(radius):
for dx in (px, -px):
for dy in (py, -py):
# Done once per quadrant. Do both octants.
field[max(0, y + dy)][max(0, x + dx)] = entry
field[max(0, y + dx)][max(0, x + dy)] = entry
# do the diameters
for dx in range(-radius, radius + 1):
field[y][max(0, x + dx)] = entry
field[max(0, y + dx)][x] = entry
return field

There is still overlap done at x = y and x = -y; you could squeeze that
out as well by changing wedge_nz not to put it out, and making circle2
do diagonal diameters as well (leaving the center the sole overwrite).

--Scott David Daniels
(e-mail address removed)

 

[Pythor]

[Michael Tobis]
Also, with this code, you are using radius for the dimensions of the
enclosing box, as well as the radius of the circle, so it's guaranteed
to not to actually produce a whole circle. Recall what python does with
negative indices!

[Pythor]
I'm not sure what you mean here. It produces an eighth-circle, and
then plots each point in the 8 symmetrical positions on the circle.
Except for the (dx+1) point made above, what piece of the circle is
missing?
======
What Michael means is that for example a circle of radius 5 is 11
pixels wide and 11 pixels high. You are trying to cram it into a box of
5 x 5 pixels [or maybe 6x6 (I'm numpy-challenged)]. The result will
resemble a train smash. Have you considered *TESTING* your code? It's
not very difficult at all to draw the expected results for radii of
about 4 or 5 pixels on the back of an envelope ...

Sure, I tested it. And it works, for the most part. I did miss those
the dx+1 points, which John pointed out. On the other hand, I'm not
having any trouble producing a whole circle, while you seem to think
I'm only producing half a circle. The code that limits itself to a 5x5
box is only expected to produce an eighth of the circle. The actual
assignment portion uses reflection to plot points in the whole area.
If there's some other problem with it, I haven't noticed.

By the way, there are signs of a benchmark war breaking out. What are
typical sizes you would be using in practice for the radius and the
enclosing box?

Well, I'm not really asking someone else to write my code for me, but
here goes.

The desire is to have an array of about 1000 X 1000, bigger would be
better. I want to fill this array with regions of values, with each
region being between .5% and 5% of the total area. I'm using random
circular regions placed around the array, allowing new regions to
overlap old ones. So, I'm using circles of roughly 5 to 50 radius, and
throwing a few thousand into the array.
Actually, this is a prelude to trying the same thing with spheres in a
3-dimensional array, but a 1000X1000x1000 array is larger than my
memory can handle.

 

[John Machin]

[Pythor]
Sure, I tested it.
===
I don't think that word means what you think it means

[Pythor]
On the other hand, I'm not
having any trouble producing a whole circle, while you seem to think
I'm only producing half a circle. The code that limits itself to a 5x5
box is only expected to produce an eighth of the circle. The actual
assignment portion uses reflection to plot points in the whole area.
If there's some other problem with it, I haven't noticed.
===
Here's the problem:
if field==None:
field=numpy.zeros((radius,radius),'u')

If radius is 5, field is 5x5 => train smash.