Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
import numpy as np
import overlap
from polygon_triangulate import polygon_triangulate, PlotTriangulation
from tetrahedralize_triangles import tetrahedralize
def measure_volume(V_tet, T, origin, r):
sphere = overlap.Sphere(origin, r)
volume = 0.0
for tet in T:
vertices = V_tet[tet]
overlap_tet = overlap.Tetrahedron(vertices)
overlap_volume = overlap.overlap(sphere, overlap_tet)
volume += overlap_volume
return volume
def same_point(x, y):
return np.linalg.norm(x - y) < 1e-10
def remove_duplicate_contour_points(P):
n = len(P)
for i in range(n):
j = (i + 1) % n
if same_point(P[i], P[j]):
return remove_duplicate_contour_points(np.delete(P, i, axis=0))
return P
def measure_cumulative_volume(z, P, thickness, origin, R):
'''
Integration of the intersecting volume of a thick polygon and a sphere of increasing radius.
Parameters
----------
z : float
The z depth-wise coordinate of the thick polygon.
P : (N, 2) array_like
A real array of ordered contour points of the polygon.
thickness : float
The thickness of the polygon.
origin : (3,) array_like
The center of the sphere from which to integrate from.
R : (M,) array_like
The discrete sphere radii from which to measure the polygon in.
Returns
-------
mu: (M,) array_like
The intersecting volumes of the sphere and the thick polygon at radii determined by R.
'''
P = np.array(P)
origin = np.array(origin)
P = remove_duplicate_contour_points(P)
z0 = z - thickness / 2
z1 = z + thickness / 2
_, V_tri, F = polygon_triangulate(P)
V_tet, T = tetrahedralize(V_tri, F, z0, z1)
mu = np.zeros(len(R))
for i, r in enumerate(R):
if r <= 0:
continue
mu[i] = measure_volume(V_tet, T, origin, r)
return mu
def measure_cumulative_points(z, X, thickness, origin, R):
'''
Integration of the points placed in a thick section integrated intersecting with a sphere of increasing size.
Parameters
----------
z : float
The z depth-wise coordinate of the section with the points.
X : (N, 2) array_like
A real array of points inside the section.
thickness : float
The thickness of the section of which the points are in.
origin : (3,) array_like
The center of the sphere from which to integrate from.
R : (M,) array_like
The discrete sphere radii from which to measure the points inside.
Returns
-------
mu: (M,) array_like
The integrated point masses inside the sphere of radii determined by R.
'''
oz = origin[2]
oxy = origin[:2]
dxy = np.array([np.linalg.norm(x - oxy) for x in X])
mu = np.zeros(len(R))
if z == oz:
for i, r in enumerate(R):
s = np.sqrt(np.clip(r**2 - dxy**2, 0, None))
mu[i] = np.sum(np.clip(2*s/thickness, 0, 1))
return mu
# plane is below, mirror it above to simplify calculations
if z < oz:
z = z + 2*(oz - z)
z0 = z - thickness / 2
z1 = z + thickness / 2
d0 = z0 - oz
for i, r in enumerate(R):
h = np.sqrt(np.clip(r**2 - dxy**2, 0, None))
s = h - d0
mu[i] = np.sum(np.clip(s/thickness, 0, 1))
return mu