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