[Numpy-discussion] Geometrically defined masking arrays; how to optimize?

Wolfgang Draxinger Wolfgang.Draxinger at physik.uni-muenchen.de
Tue Feb 11 15:56:01 EST 2014


I implemented the following helper function to create masking arrays:

def gen_mask(self, ring, sector):
	in_segment = None
	if ring == 0:
		radius = self.scales()[0]
		def in_center_circle(xy):
			dx = xy[0] - self.xy[0]
			dy = xy[1] - self.xy[1]
			r = math.sqrt( dx**2 + dy**2 )
			return r < radius
		in_segment = in_center_circle
		angle = ( self.a_sector(sector, ring),
			  self.a_sector( (sector+1) % self.n_sectors(), ring) )
		radii = self.scales()[ring:ring+1]
		def in_segment_ring(xy):
			dx = xy[0] - self.xy[0]
			dy = xy[1] - self.xy[1]
			r = math.sqrt( dx**2 + dy**2 )
			a = math.atan2( dy, dx )
			return r >= radii[0] and r < radii[1] \
			   and a >= angle[0] and a < angle[1]
		in_segment = in_segment_ring

	width,height = self.arr.shape

	mask = numpy.zeros(shape=(width, height), dtype=numpy.bool)
	for y in range(height):
		for x in range(width):
			mask[x][y] = in_segment((x,y))
	return mask

self.scales() generates a list of radius scaling factors.
self.a_sector() gives the dividing angle between sectors "sector" and 
"sector+1" on the given ring.

The function works, it generates the masks as I need it. The problem is 
- of course - that it's quite slow due to the inner loops the perform 
the geometrical test if the given element of the array self.arr is 
withing the bounds of the ring-sector for which the mask is generated.

I wonder if you guys have some ideas, on how I could accelerate this. 
This can be perfectly well constructed by boolean combination of 
elementary geometrical objects. For example a ring would be

ring(p, r0, r1): disk(p, r1) xor disk(p, r0) # where r0 < r1

The sector would then be

ring_sector(p, r, s): ring(p, r, r + ...) and sector(p, s, s+1)

I'd like to avoid doing this using some C helper routine. I'm looking 
for something that is the most efficient when it comes to 
"speedup"/"time invested to develop this".



More information about the NumPy-Discussion mailing list