How to generate a matrix with circle of ones in numpy/scipy

Question

There are some signal generation helper functions in python's scipy, but these are only for 1 dimensional signal.

I want to generate a 2-D ideal bandpass filter, which is a matrix of all zeros, with a circle of ones to remove some periodic noise from my image.

I am now doing with:

def unit_circle(r):
    def distance(x1, y1, x2, y2):
        return math.sqrt((x1 - x2) ** 2 + (y1 - y2) ** 2)
    d = 2*r + 1
    mat = np.zeros((d, d))
    rx , ry = d/2, d/2
    for row in range(d):
        for col in range(d):
            dist = distance(rx, ry, row, col)
            if abs(dist - r) < 0.5:
                mat[row, col] = 1
    return mat

result:

In [18]: unit_circle(6)
Out[18]:
array([[ 0.,  0.,  0.,  0.,  1.,  1.,  1.,  1.,  1.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  1.,  1.,  0.,  0.,  0.,  0.,  0.,  1.,  1.,  0.,  0.],
       [ 0.,  1.,  1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  1.,  1.,  0.],
       [ 0.,  1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  1.,  0.],
       [ 1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  1.],
       [ 1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  1.],
       [ 1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  1.],
       [ 1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  1.],
       [ 1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  1.],
       [ 0.,  1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  1.,  0.],
       [ 0.,  1.,  1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  1.,  1.,  0.],
       [ 0.,  0.,  1.,  1.,  0.,  0.,  0.,  0.,  0.,  1.,  1.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  1.,  1.,  1.,  1.,  1.,  0.,  0.,  0.,  0.]])

Is there a more direct way to generate a matrix of circle of ones, all else zeros?

Edit: Python 2.7.12

Answer 1

Here is a pure NumPy alternative that should run significantly faster and looks cleaner, imho. Basically, we vectorise your code by replacing built-in sqrt and abs with their NumPy alternatives and working on matrices of indices.

Updated to replace distance with np.hypot(courtesy of James K)

In [5]: import numpy as np

In [6]: def my_unit_circle(r):
   ...:     d = 2*r + 1
   ...:     rx, ry = d/2, d/2
   ...:     x, y = np.indices((d, d))
   ...:     return (np.abs(np.hypot(rx - x, ry - y)-r) < 0.5).astype(int)
   ...: 

In [7]: my_unit_circle(6)
Out[7]: 
array([[ 0.,  0.,  0.,  0.,  1.,  1.,  1.,  1.,  1.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  1.,  1.,  0.,  0.,  0.,  0.,  0.,  1.,  1.,  0.,  0.],
       [ 0.,  1.,  1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  1.,  1.,  0.],
       [ 0.,  1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  1.,  0.],
       [ 1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  1.],
       [ 1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  1.],
       [ 1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  1.],
       [ 1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  1.],
       [ 1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  1.],
       [ 0.,  1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  1.,  0.],
       [ 0.,  1.,  1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  1.,  1.,  0.],
       [ 0.,  0.,  1.,  1.,  0.,  0.,  0.,  0.,  0.,  1.,  1.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  1.,  1.,  1.,  1.,  1.,  0.,  0.,  0.,  0.]])

Benchmarks

In [12]: %timeit unit_circle(100)
100 loops, best of 3: 17.7 ms per loop

In [13]: %timeit my_unit_circle(100)
1000 loops, best of 3: 480 µs per loop

Answer 2

Here's a vectorized approach -

def unit_circle_vectorized(r):
    A = np.arange(-r,r+1)**2
    dists = np.sqrt(A[:,None] + A)
    return (np.abs(dists-r)<0.5).astype(int)

Runtime test -

In [165]: %timeit unit_circle(100) # Original soln
10 loops, best of 3: 31.1 ms per loop

In [166]: %timeit my_unit_circle(100) #@Eli Korvigo's soln
100 loops, best of 3: 2.68 ms per loop

In [167]: %timeit unit_circle_vectorized(100)
1000 loops, best of 3: 582 µs per loop

Answer 3

result of code execution

def gen_circle(img: np.ndarray, center: tuple, diameter: int) -> np.ndarray:


  
    """
        Creates a matrix of ones filling a circle.
    """

    # gets the radious of the image
    radious  = diameter//2

    # gets the row and column center of the image
    row, col = center 

    # generates theta vector to variate the angle
    theta = np.arange(0, 360)*(np.pi/180)

    # generates the indexes of the column
    y = (radious*np.sin(theta)).astype("int32") 

    # generates the indexes of the rows
    x = (radious*np.cos(theta)).astype("int32") 

    # with:
    # img[x, y] = 1
    # you can draw the border of the circle 
    # instead of the inner part and the border. 

    # centers the circle at the input center
    rows = x + (row)
    cols  = y + (col)

    # gets the number of rows and columns to make 
    # to cut by half the execution
    nrows = rows.shape[0] 
    ncols = cols.shape[0]

    # makes a copy of the image
    img_copy = copy.deepcopy(img)

    # We use the simetry in our favour
    # does reflection on the horizontal axes 
    # and in the vertical axes

    for row_down, row_up, col1, col2 in zip(rows[:nrows//4],
                            np.flip(rows[nrows//4:nrows//2]),
                            cols[:ncols//4],
                            cols[nrows//2:3*ncols//4]):
    
        img_copy[row_up:row_down, col2:col1] = 1

 
    return img_copy

center = (30,40)
ones = np.zeros((center[0]*2, center[1]*2))
diameter = 30

circle = gen_circle(ones, center, diameter)
plt.imshow(circle)