[Numpy-discussion] floating point arithmetic issue

[Pauli Virtanen]

Fri, 30 Jul 2010 14:21:23 +0200, Guillaume Chérel wrote:
[clip]
> As for the details about my problem, I'm trying to compute the total
> surface of overlapping disks. I approximate the surface with a grid and
> count how many points of the grid fall into at least one disk.

HTH,

import numpy as np

# some random circles
x = np.random.randn(200)
y = np.random.randn(200)
radius = np.random.rand(200)*0.1

# grid, [-1, 1] x [-1, 1], 200 x 200 points
xmin, xmax = -1, 1
ymin, ymax = -1, 1
grid_x, grid_y = np.mgrid[xmin:xmax:200j, ymin:ymax:200j]

# count
mask = np.zeros((200, 200), bool)
for xx, yy, rr in zip(x, y, radius):
    mask |= (grid_x - xx)**2 + (grid_y - yy)**2 < rr**2

area = mask.sum() * (xmax-xmin) * (ymax-ymin) / float(mask.size)