[SciPy-User] Integration over Voronoi cells

Eraldo Pomponi eraldo.pomponi@gmail....
Thu Nov 24 08:46:30 CST 2011


Dear Friedrich,

First of all let me thank you for your suggestions. It took a little bit to
get into them
but I still have some BIG doubts.


> 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:
>

Not so clear why the Gaussian is approximately one (it is true just when
r is small enough)


> 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.


Let me rewrite the function (it is just the kernel not the full function):

K(r,t) = 1/(2piDr)exp[-r^2/Dt]                       [1]

This function is rotational invariant [a] so a clever way to integrate it
is to consider
a disk of radius R0. The perimeter of this disk is equal to:

l=2*pi*R0             [2]

The term  1/(2*pi*D*r) gives a constant contribution along the path eq.2,
equal to 1/D so
what is left to do is to integrate the term:

exp[-r^2/Dt]

in the interval [R0/D,R0/D+delta] and multiply by the length of the chose
path, i.e. eq.2. and sum
on a 1D grid in the interval [0,R0].

------>>>>>   Is it correct?  <<<<<-----------

The term exp[-r^2/Dt] can be analytically integrated (removing t for
convenience) :

g(r) = sqrt(pi)*erf(r*sqrt(1/D))/(2*sqrt(1/D))   [3]

Its defined integral in [R0/D,R0/D+delta],  multiplied by the length of the
choose path,
give us:

g'(r) = 2*pi*r * [ g(r/D + delta) - g(r/D)]    [4]

summing in the grid k=1..R0, we obtain the integral of  eq.1 in [0,R0]:

integral([1]) =  sum([g'(k) for k in np.arange(0,R0,0.001)] )

What do you think ?


>

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.


I didn't understood what you mean here. Sorry


> 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.
>
>
Yes , I hope so. When the path in eq.2 is no more circular, the integration
along it is more complicated. A bit difficult to figure out how to do that
programmatically but this is the next problem.


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

I guess so but I think that I will fully understand this sentence just when
my
doubts will be cleared. Thanks a lot for your help.

Cheers,
Eraldo
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mail.scipy.org/pipermail/scipy-user/attachments/20111124/1fd374c3/attachment.html 


More information about the SciPy-User mailing list