Changing the values in a numpy array depending on the area of a circle

Question

Context

I need to measure the area of combined circles in python. I came up with way using numpy arrays. First, I populate a grid (numpy array) with zeros, each position in the grid corresponds to a length of 0.5cm. I then drop the centre of the circle onto the grid and change this value to 1 in the grid. I know the radius of the circles so I can calculate the area of the circle, because I know the area of the circle I change the the zeros in the grid to which fall in the area of the circle ones. I then count the frequency of ones in the gird and use this to calculate the area of the combined circles, because I know the length of each position in the grid I can calculate the area. This is currently a very coarse grain method, I plan to change it to a finer grain once I have the algorithm nailed down.

Example

If you look at the image I posted below which describes my idea better. There are two circles on my grid (red lines), the centre of the circles are marked with a blue square and the area the circles occupy are in a light orange. I want to change the the area with marked with orange to ones. I can currently change the the orange squares horizontal and vertical to the circle centre, but the diagonal boxes from the centre are causing me trouble.

enter image description here

Current code

class area():

    def make_grid(self):
        '''
        Each square in the grid represents 0.5 cm
        '''
        import numpy as np
        grid = np.zeros((10,10))
        square_length = 0.5
        circles = {'c1':[[4,2],1.5],'c2':[[5,6],2.0]}

        print grid
        for key,val in circles.iteritems():
            grid[val[0][0]][val[0][1]] = 1
            area = int((val[1] - square_length)/0.5)            
            for i in xrange(1,area+1):
                grid[val[0][0]][val[0][1]+i] = 1 # Change column vals in +ve direction
                grid[val[0][0]][val[0][1]-i] = 1 # Chnage column vals in -ve direction
                grid[val[0][0]+i][val[0][1]] = 1 # Chnage row vals in +ve direction 
                grid[val[0][0]-i][val[0][1]] = 1 # Chnage row vals in -ve direction

        print ''
        print grid

In the dictionary above, the key is the circle name, the first element in the value is the circle centre co-ordinates and the second element is the radius of the circle.

The code outputs:

[[ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.]
 [ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.]
 [ 0.  0.  1.  0.  0.  0.  1.  0.  0.  0.]
 [ 0.  0.  1.  0.  0.  0.  1.  0.  0.  0.]
 [ 1.  1.  1.  1.  1.  0.  1.  0.  0.  0.]
 [ 0.  0.  1.  1.  1.  1.  1.  1.  1.  1.]
 [ 0.  0.  1.  0.  0.  0.  1.  0.  0.  0.]
 [ 0.  0.  0.  0.  0.  0.  1.  0.  0.  0.]
 [ 0.  0.  0.  0.  0.  0.  1.  0.  0.  0.]
 [ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.]]

Shouldn't you instead be stepping through the full grid, calculate the distance to each circle centre for each grid point, then see if it's within the circle and upp the grid point counter by one? That way, some points will have 2 instead of 1, but at least all the points within the circle (rounding problems excepted) will have been counted. It is slower, 'cause you have to step through all the points, but it's probably more straightforward. — , Mar 15 '13 at 12:27
Do you consider geometric solution "to measure the area"? http://mathworld.wolfram.com/Circle-CircleIntersection.html — MBo, Mar 15 '13 at 12:35
I have considered a geometric solution, but a problem arises if I have multiple circles with 'missing areas' as described in this question: http://stackoverflow.com/questions/1667310/combined-area-of-overlapping-circles — Harpal, Mar 15 '13 at 12:38

YXD · Accepted Answer · 2013-03-15T15:35:16.920

Updated according to comments from @Jaime.

I'd probably start with something like this. The main point is the correct calculation of the pixels inside the circle.

import numpy as np
import matplotlib.pyplot as plt

grid = np.zeros((10,10), dtype=np.bool)
square_length = 0.5
circles = {'c1':[[4,2],1.5],'c2':[[5,6],2.0]}

# Generate arrays of indices/coordiates so we can do the
# calculations the Numpy way, without resorting to loops
# I always get the index order wrong so double check...
xx = np.arange(grid.shape[0])
yy = np.arange(grid.shape[1])

for val in circles.itervalues():
    radius = val[1]
    # same index caveat here
    # Calling Mr Pythagoras: Find the pixels that lie inside this circle
    inside = (xx[:,None] - val[0][0]) ** 2 + (yy[None, :] - val[0][1]) ** 2 <= (radius ** 2)
    # do grid & inside and initialize grid with ones for intersection instead of union
    grid = grid | inside 

plt.imshow(grid)
plt.show()

If you are after the intersection, note that your circle centers are sqrt(17) ~= 4.123.. units apart and the sum of the two radii comes to 3.5 so there is actually no overlap.

+1 I think this is the way to go. I would get rid of the `np.sqrt` and compare to `radius**2`. Also, by using `&`, I think you are determining points overlapped by **all** circles. Use a `grid` of `zeros` and `|` instead to get points inside **any** circle. You could also not use `meshgrid`, but have `xx = np.arange(grid.shape[0])` and `yy = np.arange(grid.shape[1])`, and then use broadcasting as `inside = (xx[:, None] - val[0][0])**2 + (yy - val[0][1])**2 <= radius**2` to have a less memory hungry implementation. — Jaime, Mar 15 '13 at 14:33
Ah, yeah I understood "combined circles" to mean the intersection but it makes sense to have the union. Thanks for the comments - will update my answer. — YXD, Mar 15 '13 at 15:27

ghdalum · Answer 2 · 2013-03-15T12:40:50.970

You are only searching directly right and left from the center. For getting the squares diagonal from the center you probably have to use the Pythagorean theorem. With the theorem you can find the hypotenuse of a square by using the offset both horizontally and vertically in relative to the center of the circle as sidelengths. You can then compare the hypotenuse to the radius, and increase the square-value by one if the hypotenuse is the shorter of the two.

Also the usage of radius is a little weird, since it is not using the middle of the square as the center. This makes a circle with a radius of 2 have a diameter of 3.5.

score 1 · Answer 3 · answered Mar 15 '13 at 13:02

You could use Monte carlo method to find the area of the overlapped area. That way it would be more accurate. http://www.chem.unl.edu/zeng/joy/mclab/mcintro.html

Find the bounds of the square, and then fill with random numbers, and for each random number check if the random number is inside the circle.

if (np.sqrt((xx - val[0][0]) ** 2 + (yy - val[0][1]) ** 2) <= radius) :inside=inside+1

Area=total no of pixels inside circle/ Total number of pixels generated*Area of square

score 0 · Answer 4 · answered Mar 15 '13 at 12:58

0

You can use Floodfill algorithm with slight modification of stop condition: take into account both color and distance to circle center

P.S. Hack method for low-size regions: Draw filled circles and count color pixels...

answered Mar 15 '13 at 12:58

MBo

77,366
5
53
86

Changing the values in a numpy array depending on the area of a circle

4 Answers4