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

[Daπid]

A small improvement:

# Dimensions
N = 100
M = 200

# Coordinates of the centre
x0 = 40
y0 = 70

x, y = np.meshgrid(np.arange(N) - x0, np.arange(M) - y0, sparse=True)

# Center at 40, 70, radius 20
out_circle = (x * x + y * y < 20**2)

out_ring = out_circle - (x * x + y * y < 10**2)

plt.imshow(out_circle)
plt.figure()
plt.imshow(out_ring)
plt.show()

Using sparse you can avoid the repetition of the arrays, getting the same
functionality. Also, if the image is big, you can delegate the computations
of out_circle and out_ring to numexpr for speed.

/David.


On 11 February 2014 22:16, Daπid <davidmenhur at gmail.com> wrote:

> Here it is an example:
>
> import numpy as np
> import pylab as plt
>
> N = 100
> M = 200
> x, y = np.meshgrid(np.arange(N), np.arange(M))
>
> # Center at 40, 70, radius 20
> x -= 40
> y -= 70
> out_circle = (x * x + y * y < 20**2)
>
> out_ring = out_circle - (x * x + y * y < 10**2)
>
> plt.imshow(out_circle)
> plt.figure()
> plt.imshow(out_ring)
> plt.show()
>
> Note that I have avoided taking the costly square root of each element by
> just taking the square of the radius. It can also be generalised to
> ellipses or rectangles, if you need it.
>
> Also, don't use 0 as a False value, and don't force it to be a 0. Instead,
> use "if not ring:"
>
>
> /David
>
>
>
>
>
> On 11 February 2014 21:56, Wolfgang Draxinger <
> Wolfgang.Draxinger at physik.uni-muenchen.de> wrote:
>
>> Hi,
>>
>> 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
>>         else:
>>                 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".
>>
>>
>> Cheers,
>>
>> Wolfgang
>> _______________________________________________
>> NumPy-Discussion mailing list
>> NumPy-Discussion at scipy.org
>> http://mail.scipy.org/mailman/listinfo/numpy-discussion
>>
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/numpy-discussion/attachments/20140211/8a680387/attachment.html>