Rotating an array - and manipulating it afterwards

Question

I am trying to create an array that is (101,101) and the value of each number in the array is given by it's distance from the center.

So far, I have managed to get this working in the 'top quarter' of my array. I found a method for rotating my array, which also works. But what I want to do, is rotate the array and then apply the same piece of code to it, to get the desired effect on the next 'quadrant'.

from pylab import *

close()

z=ones((101,101),dtype=integer)

a=0
b=0
c=100
d=50
while a<=100 and b<=100 and c>=0 and d<=50:
    z[a,b:c]=d
    a=a+1
    b=b+1
    c=c-1
    d=d-1

x=zip(*z[::-1])

when I apply the same code to my new array 'x', I can't get anything to happen.

from pylab import *

close()

z=ones((101,101),dtype=integer)

a=0
b=0
c=100
d=50
while a<=100 and b<=100 and c>=0 and d<=50:
    z[a,b:c]=d
    a=a+1
    b=b+1
    c=c-1
    d=d-1

    x=zip(*z[::-1])


x=zip(*z[::-1])
while a<=100 and b<=100 and c>=0 and d<=50:
    x[a,b:c]=d
    a=a+1
    b=b+1
    c=c-1
    d=d-1
    g=zip(*x[::,-1])

imshow(g)
show()

This just tells me that g is not defined, but it looks like it should work to me... Please can someone advise me on what's going wrong?

Thank you!

Answer 1

I have added to the answer below a beefed up. more general answer and. later, a warning against non vectorized inner loops in numpy

My answer

I have a function for you, it works only when n is odd but it works OK when n is odd...

def radius(n):
    if (n+1)%2 : return None
    a = np.zeros((n,n))
    for r in range(n/2,n):
        x = r - n/2
        for c in range(n/2,n):
            y = c-n/2
            a[r,c] = np.sqrt(x*x+y*y)
    a[n/2-1::-1,n/2:] = a[n/2+1:,n/2:]
    a[:,n/2-1::-1] = a[:,n/2+1:]
    return a

Note that this is suboptimal as you may want to compute n/2 once and reuse the value the following ten times it is used...

Upgrade

This is simply better, especially it works OK also for even values of n.

def radius(n):
    import numpy as np
    if n<1 : return
    if n == 1 : return np.zeros((1,1))
    a = np.zeros((n,n))
    hn = n/2
    C = float(n-1)/2.0
    for r in range(hn,n):
        y2= (r-C)**2
        for c in range(hn,n):
            x2 = (c-C)**2
            a[r,c] = np.sqrt(y2+x2)
    if n%2:
        a[hn-1::-1,hn:] = a[hn+1:,hn:]
        a[:,hn-1::-1] = a[:,hn+1:]
    else:
        a[hn-1::-1,hn:] = a[hn:,hn:]
        a[:,hn-1::-1] = a[:,hn:]
    return a

A quick test

(This is a test of version 1 only, sorry, but I've tested version 2 as well...)

In [1]: import numpy as np

In [2]: def radius(n):
    if (n+1)%2 : return None

    a = np.zeros((n,n))
    for r in range(n/2,n):
        x = r - n/2
        for c in range(n/2,n):
            y = c-n/2
            a[r,c] = np.sqrt(x*x+y*y)
    a[n/2-1::-1,n/2:] = a[n/2+1:,n/2:]
    a[:,n/2-1::-1] = a[:,n/2+1:]
    return a
   ...: 

In [3]: np.set_printoptions(linewidth=120, precision=2)

In [4]: radius(14)

In [5]: radius(15)
Out[5]: 
array([[ 9.9 ,  9.22,  8.6 ,  8.06,  7.62,  7.28,  7.07,  7.  ,  7.07,  7.28,  7.62,  8.06,  8.6 ,  9.22,  9.9 ],
       [ 9.22,  8.49,  7.81,  7.21,  6.71,  6.32,  6.08,  6.  ,  6.08,  6.32,  6.71,  7.21,  7.81,  8.49,  9.22],
       [ 8.6 ,  7.81,  7.07,  6.4 ,  5.83,  5.39,  5.1 ,  5.  ,  5.1 ,  5.39,  5.83,  6.4 ,  7.07,  7.81,  8.6 ],
       [ 8.06,  7.21,  6.4 ,  5.66,  5.  ,  4.47,  4.12,  4.  ,  4.12,  4.47,  5.  ,  5.66,  6.4 ,  7.21,  8.06],
       [ 7.62,  6.71,  5.83,  5.  ,  4.24,  3.61,  3.16,  3.  ,  3.16,  3.61,  4.24,  5.  ,  5.83,  6.71,  7.62],
       [ 7.28,  6.32,  5.39,  4.47,  3.61,  2.83,  2.24,  2.  ,  2.24,  2.83,  3.61,  4.47,  5.39,  6.32,  7.28],
       [ 7.07,  6.08,  5.1 ,  4.12,  3.16,  2.24,  1.41,  1.  ,  1.41,  2.24,  3.16,  4.12,  5.1 ,  6.08,  7.07],
       [ 7.  ,  6.  ,  5.  ,  4.  ,  3.  ,  2.  ,  1.  ,  0.  ,  1.  ,  2.  ,  3.  ,  4.  ,  5.  ,  6.  ,  7.  ],
       [ 7.07,  6.08,  5.1 ,  4.12,  3.16,  2.24,  1.41,  1.  ,  1.41,  2.24,  3.16,  4.12,  5.1 ,  6.08,  7.07],
       [ 7.28,  6.32,  5.39,  4.47,  3.61,  2.83,  2.24,  2.  ,  2.24,  2.83,  3.61,  4.47,  5.39,  6.32,  7.28],
       [ 7.62,  6.71,  5.83,  5.  ,  4.24,  3.61,  3.16,  3.  ,  3.16,  3.61,  4.24,  5.  ,  5.83,  6.71,  7.62],
       [ 8.06,  7.21,  6.4 ,  5.66,  5.  ,  4.47,  4.12,  4.  ,  4.12,  4.47,  5.  ,  5.66,  6.4 ,  7.21,  8.06],
       [ 8.6 ,  7.81,  7.07,  6.4 ,  5.83,  5.39,  5.1 ,  5.  ,  5.1 ,  5.39,  5.83,  6.4 ,  7.07,  7.81,  8.6 ],
       [ 9.22,  8.49,  7.81,  7.21,  6.71,  6.32,  6.08,  6.  ,  6.08,  6.32,  6.71,  7.21,  7.81,  8.49,  9.22],
       [ 9.9 ,  9.22,  8.6 ,  8.06,  7.62,  7.28,  7.07,  7.  ,  7.07,  7.28,  7.62,  8.06,  8.6 ,  9.22,  9.9 ]])
In [6]: 

Addendum et Caveat Emptor

If you need to use the code above repeatedly and/or for large values of n, please have a look at the following timings and change the double loop to something that numpy can vectorize...

In [71]: def do(n):
    a = np.zeros((n,n))
    for r in range(n):
        y2 = r*r
        for c in range(n):
            a[r,c] = np.sqrt(y2+c*c)
    return a
   ....: 

In [72]: do(7)
Out[72]: 
array([[ 0.  ,  1.  ,  2.  ,  3.  ,  4.  ,  5.  ,  6.  ],
       [ 1.  ,  1.41,  2.24,  3.16,  4.12,  5.1 ,  6.08],
       [ 2.  ,  2.24,  2.83,  3.61,  4.47,  5.39,  6.32],
       [ 3.  ,  3.16,  3.61,  4.24,  5.  ,  5.83,  6.71],
       [ 4.  ,  4.12,  4.47,  5.  ,  5.66,  6.4 ,  7.21],
       [ 5.  ,  5.1 ,  5.39,  5.83,  6.4 ,  7.07,  7.81],
       [ 6.  ,  6.08,  6.32,  6.71,  7.21,  7.81,  8.49]])

In [73]: n = 7 ; a = np.arange(n) ; o = np.ones(n) ; np.sqrt(np.outer(o,a*a)+np.outer(a*a,o))
Out[73]: 
array([[ 0.  ,  1.  ,  2.  ,  3.  ,  4.  ,  5.  ,  6.  ],
       [ 1.  ,  1.41,  2.24,  3.16,  4.12,  5.1 ,  6.08],
       [ 2.  ,  2.24,  2.83,  3.61,  4.47,  5.39,  6.32],
       [ 3.  ,  3.16,  3.61,  4.24,  5.  ,  5.83,  6.71],
       [ 4.  ,  4.12,  4.47,  5.  ,  5.66,  6.4 ,  7.21],
       [ 5.  ,  5.1 ,  5.39,  5.83,  6.4 ,  7.07,  7.81],
       [ 6.  ,  6.08,  6.32,  6.71,  7.21,  7.81,  8.49]])

In [74]: %timeit do(1001)
1 loops, best of 3: 2.45 s per loop

In [75]: %timeit n = 1001 ; a = np.arange(n) ; o = np.ones(n) ; np.sqrt(np.outer(o,a*a)+np.outer(a*a,o))
100 loops, best of 3: 13.2 ms per loop

In [76]: 

As you can see, the non-vectorized code (as in my example codes) is almost 200 times slower than the vectorized version!

Answer 2

The problem is that by rotating the array like this x = zip(*z[::-1]) you are turning the numpy.array into a regular Python list, and for those you can not use tuple indexing like z[a,b:c] = d.

Instead of rotating the array and then filling out the same quarter again, why don't you just fill out all the quarters in the same loop? You just have to flip the ranges in the array.

Also, I'd make this a for loop to make things easier. Try this:

z = ones((101, 101), dtype=integer)
for d in range(51):
    a = 50 - d
    b = 50 + d
    z[a, a:b+1] = d
    z[b, a:b+1] = d
    z[a:b+1, a] = d
    z[a:b+1, b] = d
imshow(z)
show()

Or, if you want to keep your "rotating" approach, make sure you turn the rotated array back into an actual numpy.array. And again, I would recommend using two for loops for this.

for i in range(4):
    for d in range(51):
        a = 50 - d
        b = 50 + d
        z[a, a:b] = d
    z = array(zip(*z[::-1]))

Or fill out entire blocks of the array in one go, working your way from the outside inwards:

for d in range(50, 0, -1):
    a = 50 - d
    b = 51 + d
    z[a:b, a:b] = d