In [2]:
import numpy as np
def circumscribed_circle(X, Y, Z, W=None):
if W is None:
W = X ** 2 + Y ** 2 + Z ** 2
O = np.ones(4)
Dx = + np.linalg.det([W, Y, Z, O])
Dy = - np.linalg.det([W, X, Z, O])
Dz = + np.linalg.det([W, X, Y, O])
a = np.linalg.det([X, Y, Z, O])
c = np.linalg.det([W, X, Y, Z])
C = np.array([Dx, Dy, Dz]) / 2 / a
r = sqrt(Dx ** 2 + Dy ** 2 + Dz ** 2 - 4 * a * c) / 2 / fabs(a)
return C, r
def draw_sphere(C, r, ax):
u = np.linspace(0, 2 * np.pi, 100)
v = np.linspace(0, np.pi, 100)
x = r * np.outer(np.cos(u), np.sin(v)) + C[0]
y = r * np.outer(np.sin(u), np.sin(v)) + C[1]
z = r * np.outer(np.ones(np.size(u)), np.cos(v)) + C[2]
ax.plot_surface(x, y, z, alpha=0.2, rstride=4, cstride=4, color='b', linewidth=0.5)
In [3]:
%matplotlib
import numpy as np
from math import sqrt, fabs
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
from pyhull.convex_hull import ConvexHull
class Pole:
__slots__ = ("r", "c")
def __init__(self, r, c):
self.r = r
self.c = c
class Triangle:
__slots__ = ("v0", "v1", "check_vertices")
def __init__(self, v0, check_vertex):
self.v1 = None
self.v0 = v0
self.check_vertices = (check_vertex, )
def add_simplex(self, v1, check_vertex):
assert self.v1 is None
self.check_vertices = (self.check_vertices[0], check_vertex)
self.v1 = v1
X, Y, Z = np.transpose(np.reshape(np.fromfile("resources/bun_zipper_res4.xyz", sep=" "), (-1, 3)))
# X, Y, Z = np.transpose(np.reshape(np.fromfile("resources/bun_zipper.xyz", sep=" "), (-1, 3)))
W = X ** 2 + Y ** 2 + Z ** 2
O = np.ones(4)
ch = ConvexHull(np.transpose([X, Y, Z, W]))
simplices = [s for s in ch.vertices if np.linalg.det([X[s], Y[s], Z[s], O]) < 0]
# for s in simplices:
# C, r = circumscribed_circle(X[s], Y[s], Z[s], W[s])
# for i in range(len(X)):
# if not i in s:
# v = [X[i], Y[i], Z[i]]
# if np.linalg.norm(v - C) < r:
# print(s, i, np.linalg.norm(v - C), r)
poles = {}
triangles = {}
for s in simplices:
C, r = circumscribed_circle(X[s], Y[s], Z[s], W[s])
for v in s:
if not v in poles or poles[v].r < r:
poles[v] = Pole(r, C)
for a, b, c, d in (0, 1, 2, 3), (0, 2, 3, 1), (0, 1, 3, 2), (1, 2, 3, 0):
t_idx = tuple(sorted((s[a], s[b], s[c])))
if t_idx in triangles:
triangles[t_idx].add_simplex(C, s[d])
else:
triangles[t_idx] = Triangle(C, s[d])
Nx, Ny, Nz = np.transpose([(np.array(poles[i].c) - np.array([X[i], Y[i], Z[i]])) / poles[i].r for i in range(len(X))])
def intersection_check(triangle, vertex, echo=False):
if triangle.v1 is None:
return False
v = np.array([X[vertex], Y[vertex], Z[vertex]])
vn = np.array([Nx[vertex], Ny[vertex], Nz[vertex]])
v0, v1 = map(np.array, (triangle.v0, triangle.v1))
d0, d1 = sorted(np.dot(vn, r / np.linalg.norm(r)) for r in (v0 - v, v1 - v))
if echo:
print(d0, d1)
if d1 <= -0.38 or d0 >= 0.38:
return False
return True
# candidate_triangles = {idx: triangle for idx, triangle in triangles.items()
# # if any(intersection_check(triangle, v) for v in idx + triangle.check_vertices)}
# if all(intersection_check(triangle, v) for v in idx)}
# assert len(poles) == len(X)
# with open("test.off", "w") as out:
# out.write("OFF\n")
# out.write("{} {} 0\n".format(len(X), len(candidate_triangles)))
# for v in zip(X, Y, Z):
# out.write("{} {} {}\n".format(*v))
# for f in candidate_triangles.keys():
# out.write("3 {} {} {}\n".format(*f))
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d', aspect="equal")
# def perimeter(a, b, c):
# v = [np.array([X[i], Y[i], Z[i]]) for i in (a, b, c)]
# return sum(np.linalg.norm(v[i] - v[i - 1]) for i in range(3))
# candidate_triangles_keys = sorted(candidate_triangles.keys())
# max_trg = max(candidate_triangles.keys(), key=lambda a: perimeter(*a))
# ax.plot_trisurf(X, Y, Z, triangles=[max_trg], color='r', alpha=1, linewidth=1)
# ax.plot_trisurf(X, Y, Z, triangles=list(candidate_triangles.keys()), color='r', alpha=1, linewidth=1)
# trg = candidate_triangles[max_trg]
# t = list(max_trg) + list(trg.check_vertices)
# c1 = list(max_trg) + [trg.check_vertices[0]]
# c2 = list(max_trg) + [trg.check_vertices[1]]
# draw_sphere(*circumscribed_circle(X[c1], Y[c1], Z[c1], W[c1]), ax)
# draw_sphere(*circumscribed_circle(X[c2], Y[c2], Z[c2], W[c2]), ax)
ax.scatter(Z, X, Y, c='g')
ax.quiver(Z, X, Y, Nz, Nx, Ny, length=0.01, pivot="tail")
# ax.plot(*np.transpose([trg.v0, trg.v1]))
# ax.scatter(*trg.v0)
# ax.scatter(*trg.v1)
ax.axis([-0.2, -0.02, 0.02, 0.2])
ax.set_zlim(-0.08, 0.08)
# for i in max_trg + trg.check_vertices:
# print(i, intersection_check(trg, i, True))
plt.show()
In [5]:
%matplotlib
import numpy as np
from math import sqrt, fabs
X, Y, Z = np.random.uniform(1, 9, size=(3, 4))
C, r = circumscribed_circle(X, Y, Z)
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d', aspect="equal")
ax.scatter(X, Y, Z, c="b")
ax.scatter(*C, c="r")
draw_sphere(C, r, ax)
ax.axis([0, 10, 0, 10])
ax.set_zlim(0, 10)
plt.show()