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

Daπid davidmenhur at gmail.com
Tue Feb 11 16:16:46 EST 2014

```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:
>
>         in_segment = None
>         if ring == 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 )
>                 in_segment = in_center_circle
>         else:
>                 angle = ( self.a_sector(sector, ring),
>                           self.a_sector( (sector+1) % self.n_sectors(),
> ring) )
>                 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 )
>                            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):
>
> 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/cb50a7bb/attachment.html>
```