[974] | 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 | |
---|