-
Notifications
You must be signed in to change notification settings - Fork 0
/
utils.py
61 lines (50 loc) · 2.03 KB
/
utils.py
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
def circumsphere(vertices):
'''
@param vertices: an array of simplices
@returns: radius of circumsphere, coordinates of circumsphere
@rtype: tuple
'''
# based on https://westy31.home.xs4all.nl/Circumsphere/ncircumsphere.htm and
# https://codereview.stackexchange.com/questions/77593/calculating-the-volume-of-a-tetrahedron
from scipy.spatial import distance
# distances
squared_dists = distance.pdist(vertices, metric='sqeuclidean')
n = distance.num_obs_y(squared_dists)
# add ones and a zero to make a Cayley-Menger matrix.
squared_dists_mat = distance.squareform(squared_dists)
with_border = np.insert(np.insert(squared_dists_mat, 0, values=1, axis=1), 0, values=1, axis=0)
np.fill_diagonal(with_border, 0)
#the first diagonal element of its inverse holds -2*r*r
#and first row / first column hold barycentric coordinates of the sphere
inv = np.linalg.inv(with_border)
r = math.sqrt(inv[0][0] / -2)
barycentric_coodinates = inv[1:, 0]
return r, bary2cart(vertices, barycentric_coodinates)
def bary2cart(simplex, point):
simplex, point = np.asarray(simplex), np.asarray(point)
return np.dot(simplex.T, point)
def check_point_in_sphere(sphere, r, point):
'''
@param sphere: coodinates of a shpere
@type sphere: list
@param r: raduius of a sphere
@type r: float
@param point: coordinates of a point to check
@type point: list
@returns: boolean value
@rtype: bool
'''
sphere, point = np.asarray(sphere), np.asarray(point)
if np.sum((sphere - point)**2) <= r ** 2:
return True
else:
return False
'''
Example:
pentachoron = [[1, 1, 1, -1.0 / math.sqrt(5)], [1, -1, -1, -1.0 / math.sqrt(5)], [-1, 1, -1, -1.0 / math.sqrt(
5)], [-1, -1, 1, -1.0 / math.sqrt(5)], [0, 0, 0, math.sqrt(5) - 1.0 / math.sqrt(5)]]
tetrahedron = [[1, 1, 1], [1, -1, -1], [-1, 1, -1], [-1, -1, 1]]
triangle = [(3, 3), (1, 3), (1, 1)]
print circumsphere(pentachoron)
print circumsphere(tetrahedron)
print circumsphere(triangle)'''