1 | import numpy as np |
---|
2 | from numpy import arcsin, arctan, tan, sqrt, cross |
---|
3 | from matplotlib.collections import LineCollection |
---|
4 | import matplotlib.pyplot as plt |
---|
5 | |
---|
6 | import numba |
---|
7 | jit=numba.jit(nopython=True, parallel=True) |
---|
8 | |
---|
9 | ############### spherical primitives ################# |
---|
10 | |
---|
11 | def dotprod(a,b): return a[0]*b[0]+a[1]*b[1]+a[2]*b[2] |
---|
12 | @jit |
---|
13 | def norm2(a): return a[0]*a[0]+a[1]*a[1]+a[2]*a[2] |
---|
14 | @jit |
---|
15 | def norm(a): return np.sqrt(a[0]*a[0]+a[1]*a[1]+a[2]*a[2]) |
---|
16 | def unit_vector(a): return a/np.sqrt((a*a).sum()) |
---|
17 | def arc(a,b): return 2.*np.arcsin(.5*norm(a-b)) |
---|
18 | |
---|
19 | @jit |
---|
20 | def cross_jit(a,b, c): |
---|
21 | c[0]=a[1]*b[2]-b[1]*a[2] |
---|
22 | c[1]=a[2]*b[0]-b[2]*a[0] |
---|
23 | c[2]=a[0]*b[1]-b[0]*a[1] |
---|
24 | |
---|
25 | def area(a,b,c): |
---|
26 | # http://mathworld.wolfram.com/LHuiliersTheorem.html |
---|
27 | a,b,c = .25*arc(b,c), .25*arc(c,a), .25*arc(a,b) |
---|
28 | area = sqrt(tan(a+b+c)*tan(a+b-c)*tan(b+c-a)*tan(a+c-b)) |
---|
29 | # print("area(a,b,c):",a,b,c, a+b-c, b+c-a, a+c-b, area) |
---|
30 | return 4.*arctan(area) |
---|
31 | |
---|
32 | def areas(points, simplices): |
---|
33 | for simplex in simplices: |
---|
34 | a,b,c = points[simplex,:] |
---|
35 | yield area(a,b,c) |
---|
36 | |
---|
37 | def circum(A,B,C): |
---|
38 | # see Dubos et al., GMD 2015, appendix |
---|
39 | # the computation is stable, and valid whether a,b,c is clockwise or not |
---|
40 | A,B,C=map(unit_vector, (A,B,C)) |
---|
41 | b,c = B-A, C-A |
---|
42 | p1,p2 = cross(b,c), norm2(c)*b-norm2(b)*c |
---|
43 | p = A+.5*cross(p1,p2)/norm2(p1) |
---|
44 | p = p/norm(p) |
---|
45 | a,b,c=norm(p-A),norm(p-B),norm(p-C) |
---|
46 | if abs(a-b)>1e-15: |
---|
47 | print("Check circum", a-b, b-c, a-c) |
---|
48 | raise(ValueError) |
---|
49 | return p |
---|
50 | |
---|
51 | @jit |
---|
52 | def circum_fast(simp, points): |
---|
53 | """ simp : int[M,3] describing M simplices""" |
---|
54 | M=simp.shape[0] |
---|
55 | circ = np.zeros((M,3)) |
---|
56 | p1, p2, p = np.zeros(3), np.zeros(3), np.zeros(3) |
---|
57 | for n in range(M): |
---|
58 | i,j,k = simp[n,:] |
---|
59 | A,B,C = points[i,:], points[j,:], points[k,:] |
---|
60 | b,c = B-A, C-A |
---|
61 | cross_jit(b,c, p1) |
---|
62 | p2 = norm2(c)*b-norm2(b)*c |
---|
63 | cross_jit(p1,p2, p) |
---|
64 | p = A+p*(.5/norm2(p1)) |
---|
65 | circ[n,:]=p*(1./norm(p)) |
---|
66 | return circ |
---|
67 | |
---|
68 | def barycenters_fast(regions_size, vert, points, vertices, weight): |
---|
69 | # print("== 4 x number of edges =%d"%len(vert)) |
---|
70 | N = points.shape[0] |
---|
71 | bary = np.zeros((N,3)) |
---|
72 | i=0 # running index in vert1, vert2 |
---|
73 | A=0. |
---|
74 | for n in range(N): |
---|
75 | point = points[n,:] |
---|
76 | b = 0. |
---|
77 | for j in range(regions_size[n]): |
---|
78 | v1 = vertices[vert[i],:] |
---|
79 | i = i+1 |
---|
80 | v2 = vertices[vert[i],:] |
---|
81 | i = i+1 |
---|
82 | dA=area(point, v1, v2) |
---|
83 | A = A+dA |
---|
84 | bb = unit_vector(point+v1+v2) |
---|
85 | b = b + dA*bb*weight(bb) |
---|
86 | bary[n,:] = unit_vector(b) |
---|
87 | print("== Total area = 4*pi+%f"%(A-4*np.pi)) |
---|
88 | return bary |
---|
89 | |
---|
90 | def barycenters(regions, points, edges, vertices): |
---|
91 | N = points.shape[0] |
---|
92 | bary = np.zeros((N,3)) |
---|
93 | A=0. |
---|
94 | for n in range(N): |
---|
95 | point, region = points[n,:], regions[n] # generator and Voronoi region |
---|
96 | b = 0. |
---|
97 | for edge in region: |
---|
98 | (n1,p1), (n2,p2) = edges[edge] |
---|
99 | v1, v2 = vertices[n1,:], vertices[n2,:] |
---|
100 | dA=area(point, v1, v2) |
---|
101 | if N<10 : |
---|
102 | print("n, n1, n2", n,n1,n2) |
---|
103 | print("point,v1,v2,dA=",point,v1,v2,dA) |
---|
104 | A = A+dA |
---|
105 | bb = unit_vector(point+v1+v2) |
---|
106 | b = b + dA*bb |
---|
107 | bary[n,:] = unit_vector(b) |
---|
108 | if N<100: print("==== Total area = 4*pi+%f"%(A-4*np.pi)) |
---|
109 | return bary |
---|
110 | |
---|
111 | def plot_mesh(points,edges,orient, color): |
---|
112 | xy=points[:,0:2] |
---|
113 | segments=[ xy[[i, j],:] for i,j in edges if dotprod(orient,points[i,:])>0. ] |
---|
114 | ax = plt.axes(xlim=(-1.1, 1.1), ylim=(-1.1, 1.1)) |
---|
115 | ax.add_collection(LineCollection(segments, color=color)) |
---|
116 | |
---|