Hi all,
I did not see any voronoi diagram in matplotlib examples so I created a simple
one from the available tri.Triangulation function (I hope I did not miss
something evident).
Nicolas
#!/usr/bin/env python
# -----------------------------------------------------------------------------
# Voronoi diagram from a list of points
# Copyright (C) 2011 Nicolas P. Rougier
#
# Distributed under the terms of the BSD License.
# -----------------------------------------------------------------------------
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
def circumcircle(P1, P2, P3):
'''
Return center of the circle containing P1, P2 and P3
If P1, P2 and P3 are colinear, return None
Adapted from:
http://local.wasp.uwa.edu.au/~pbourke/geometry/circlefrom3/Circle.cpp
'''
delta_a = P2 - P1
delta_b = P3 - P2
if np.abs(delta_a[0]) <= 0.000000001 and np.abs(delta_b[1]) <= 0.000000001:
center_x = 0.5*(P2[0] + P3[0])
center_y = 0.5*(P1[1] + P2[1])
else:
aSlope = delta_a[1]/delta_a[0]
bSlope = delta_b[1]/delta_b[0]
if np.abs(aSlope-bSlope) <= 0.000000001:
return None
center_x= (aSlope*bSlope*(P1[1] - P3[1]) + bSlope*(P1[0] + P2 [0]) \
- aSlope*(P2[0]+P3[0]))/(2.*(bSlope-aSlope))
center_y = -(center_x - (P1[0]+P2[0])/2.)/aSlope + (P1[1]+P2[1])/2.
return center_x, center_y
def voronoi(X,Y):
''' Return line segments describing the voronoi diagram of X and Y '''
P = np.zeros((X.size+4,2))
P[:X.size,0], P[:Y.size,1] = X, Y
# We add four points at (pseudo) "infinity"
m = max(np.abs(X).max(), np.abs(Y).max())*1e5
P[X.size:,0] = -m, -m, +m, +m
P[Y.size:,1] = -m, +m, -m, +m
D = matplotlib.tri.Triangulation(P[:,0],P[:,1])
T = D.triangles
n = T.shape[0]
C = np.zeros((n,2))
for i in range(n):
C[i] = circumcircle(P[T[i,0]],P[T[i,1]],P[T[i,2]])
X,Y = C[:,0], C[:,1]
segments = []
for i in range(n):
for k in D.neighbors[i]:
if k != -1:
segments.append([(X[i],Y[i]), (X[k],Y[k])])
return segments
# -----------------------------------------------------------------------------
if __name__ == '__main__':
P = np.random.random((2,256))
X,Y = P[0],P[1]
fig = plt.figure(figsize=(10,10))
axes = plt.subplot(1,1,1)
plt.scatter(X,Y, s=5)
segments = voronoi(X,Y)
lines = matplotlib.collections.LineCollection(segments, color='0.75')
axes.add_collection(lines)
plt.axis([0,1,0,1])
plt.show()
------------------------------------------------------------------------------
The ultimate all-in-one performance toolkit: Intel(R) Parallel Studio XE:
Pinpoint memory and threading errors before they happen.
Find and fix more than 250 security defects in the development cycle.
Locate bottlenecks in serial and parallel code that limit performance.
http://p.sf.net/sfu/intel-dev2devfeb
_______________________________________________
Matplotlib-devel mailing list
Matplotlib-devel@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/matplotlib-devel
