[SciPy-User] Integration over Voronoi cells

[Friedrich Romstedt]

2011/11/23 Eraldo Pomponi <eraldo.pomponi at gmail.com>:
> Hi folks,
> I'm working on the integration of a function like:
> K(r,t) = 1/(2piDr)exp[-r^2/Dt]                       [1]
> over Voronoi cells (r is the distance from the point at which is associated
> the cell).
> I googled a lot and I found this two useful hints:
> http://stackoverflow.com/questions/5941113/looking-for-python-package-for-numerical-integration-over-a-tessellated-domain
> http://mathforum.org/kb/message.jspa?messageID=4963570&tstart=0
> but I'm still not able to understand how I should do this integration.
> I have the function that returns the segments ( list of
> [(x_star,y_start),(x_sop,y,stop)] ) necessary to construct the Voronoi cells
> associated to a set of points.
> Could someone suggest how to proceed ?
> There's also a numerical problem connected with the integration due to the
> singularity in r==0.
> Could you suggest which is a reasonable stable integration method available
> in scipy
> that could handle the function [1]?

Well, AISI, your function is radially symmetric, so the integration
over 1/r can be done analytically within a small disk where the
Gaussian is approximately one.  Since:

2 \pi r * [1 / (2 \pi D r)]  is just 1/D constantly.  So you're left
with [the integration up to] R_0/D when R_0 is the radius of that disk
(where the integration is done analytically).  [No more divergence,
since you cut out the centre region.]  You could also just integrate
with the full Gaussian, since the kernel just integrates up the
lengthes of the disk perimeters, if you understand.

All you need to do, AISI atm, is to calculate a function that gives
you the perimeter length of the circle which is inside your Voronoi
cell.  This is zero at r = 0 and will stay bounded for all finite
radii.  So integration should be fairly straightforward.  You could
even just do a sum over a 1D grid for the radius, since the function
varies slowly, this should be fast and easy (both in runtime as in
coding time).  You might not even need scipy for this.

Additionally, the Gaussian just introduces a suppression of radii
which are farther outside.  It makes the complete function bounded
even for an infinitely large Voronoi cell :-)

The kernel 1/r makes things easy and convenient, instead of making
things troublesome:

1)  because it exhibits rotational invariance;
2)  because it converges nicely in 2D for an integral over a circle line

I might overlook something obvious.  Is the function you gave really
the full function or only the "kernel", the weighting function?

Friedrich

P.S.: You divide your tringles resulting from the centre point and the
boundary lines into two parts, separated by the closest point on the
boundary, which is always in the bounds of that boundary line.  Then
you just have a perimeter length which is linear in r until you touch
the line, and then it'll be a littel more complicated.  I think you'll
figure it out.

P.P.S.: You're really lucky that your kernel is 1/r and not 1/r^2, iirc ;-)