In [5]:
%matplotlib
import numpy as np
import matplotlib.pyplot as plt
from pyhull.convex_hull import ConvexHull
X, Y = np.random.uniform(1, 9, size=(2, 30))
ch = ConvexHull(np.transpose([X, Y]))
plt.plot(X, Y, "ro")
for v in ch.vertices:
plt.plot(X[v], Y[v], "g")
plt.axis([0, 10, 0, 10])
plt.show()
In [6]:
%matplotlib
import numpy as np
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
from pyhull.convex_hull import ConvexHull
X, Y, Z = np.random.uniform(1, 9, size=(3, 30))
ch = ConvexHull(np.transpose([X + [10], Y + [10] , Z + [10]]))
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(X, Y, Z, c='r')
ax.plot_trisurf(X, Y, Z, alpha=0.1, triangles=ch.vertices)
plt.show()
In [10]:
%matplotlib
import numpy as np
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
with open("resources/bun_zipper_res4.off") as inp:
inp.readline() # OFF
v, f, _ = map(int, inp.readline().split()) # vertices count, faces count, edges count
X, Y, Z = np.transpose([list(map(float, inp.readline().split())) for _ in range(v)]) # vertices
surface = [list(map(int, inp.readline().split()[1:])) for _ in range(f)] # edges
fig = plt.figure()
ax = fig.add_subplot(121, projection='3d', aspect="equal")
ax.scatter(Z, X, Y, c='g')
ax = fig.add_subplot(122, projection='3d', aspect="equal")
ax.plot_trisurf(Z, X, Y, triangles=surface, color='r', linewidth=1)
plt.show()
Один из самых простых для реализации алгоритмов: "A simple algorithm for homeomorphic surface reconstruction" за авторством N. Amenta, S. Choi, T.K. Dey и N. Leekha (2000).
Ячейка диаграммы Вороного для точки $p_i$ получается пересечением полуплоскостей, образованных серединными перпендикулярами пар точек $(p_i, p_j),\space j \ne i$.
In [8]:
%matplotlib
import numpy as np
import matplotlib.pyplot as plt
from scipy.spatial import Voronoi, voronoi_plot_2d
import matplotlib
np.random.seed(239)
points = np.random.rand(200,2)
points[72] = (0, 0)
vor = Voronoi(points)
ax = plt.subplot()
voronoi_plot_2d(vor, ax, line_width=10, show_vertices=False, show_points=False)
for region in vor.regions:
if not -1 in region:
polygon = [vor.vertices[i] for i in region]
plt.fill(*zip(*polygon), alpha=0.3)
def plot_segment(ps, color='black', **kwargs):
ax.plot([p[0] for p in ps],
[p[1] for p in ps], color=color, **kwargs)
size = 0.003
min_x, min_y, max_x, max_y = 0.3, 0.3, 0.6, 0.6
for (i, j), (va, vb) in vor.ridge_dict.items():
if -1 in vor.regions[vor.point_region[i]] + vor.regions[vor.point_region[j]]:
continue # ignore infinite regions
pi, pj = points[i], points[j]
# if np.any(pi < [min_x, min_y]) or np.any(pi > [max_x, max_y])\
# or np.any(pj < [min_x, min_y]) or np.any(pj > [max_x, max_y]):
# continue # working only if both vertices are in ROI
if 74 not in (i, j):
continue
if 74 not in (i, j):
continue
#print(i, j)
plot_segment([pi, pj], linestyle='-', linewidth=0.3) # plotting line between vertices
m = (pi + pj) / 2
a = (pj - pi) / np.linalg.norm(pj - pi)
va, vb = vor.vertices[va], vor.vertices[vb]
b = (vb - va) / np.linalg.norm(vb - va)
if (np.linalg.norm(m - b*size - va) > np.linalg.norm(m + b*size - va)):
plot_segment([va, m - b*size], linestyle='-', linewidth=0.5)
else:
plot_segment([va, m + b*size], linestyle='-', linewidth=0.5)
plot_segment([m + b * size, m + b * size + a * size, m + a * size], linewidth=0.3) # 90 degree angle
plot_segment([m - b * size, m - b * size - a * size, m - a * size], linewidth=0.3) # 90 degree angle
c = np.random.rand(3, 1)
plot_segment([(pi + m) / 2 - b * size/2, (pi + m) / 2 + b * size/2], linewidth=3.0, color=c) # equal segments stroke
plot_segment([(pj + m) / 2 - b * size/2, (pj + m) / 2 + b * size/2], linewidth=3.0, color=c) # equal segments stroke
ax.axis([min_x, max_x, min_y, max_y])
plt.axis('off')
ax.set_aspect('equal')
plt.show()
Рассмотрим множество полуплоскостей, заданных нормалью $n$ и точкой $p$ на прямой, ограничивающей соответствующую полуплоскость (напомню неравенство для $i$-й полуплоскости: $n_i \cdot (p - p_i) \le 0$). Пусть пересечение полуплоскостей $CP$ непусто и нам известна точка $q \in int(CP)$. Передвинем систему координат в точку $q$. Уравнения полуплоскостей изменятся на величину $n_i \cdot p_i$, которая станет положительной для всех полуплоскостей, что позволит переписать все уравнения в виде $a_i x + b_i y + c_i \ge 0, \space c_i \gt 0$.
Полуплоскостям $(a_i, b_i, c_i)$ в исходном пространстве будут соответствовать точки $d_i = (\frac{a_i}{c_i}, \frac{b_i}{c_i})$ в двойственном пространстве. Построим в двойственном пространстве выпуклую оболочку $CH$ множества точек $\{d_i\} \cup \{0\}$.
Рассмотрим множество полуплоскостей, заданных нормалью $n$ и точкой $p$ на прямой, ограничивающей соответствующую полуплоскость (напомню неравенство для $i$-й полуплоскости: $n_i \cdot (p - p_i) \le 0$). Потребуем, чтобы все нормали были ориентированы вертикально, то есть, чтобы $n_y \gt 0$. Точка $(x,y)$ пересечения прямых, ограничивающих полуплоскости, является решением системы уравнений:
$$ \left\{ \begin{array}{ll} n_{1x} x + n_{1y} y = n_1 \cdot p_1\\ n_{2x} x + n_{2y} y = n_2 \cdot p_2 \end{array} \right. $$$$ x = \frac{ \begin{vmatrix} n_1 \cdot p_1 && n_{1y} \\ n_2 \cdot p_2 && n_{2y} \end{vmatrix}}{ \begin{vmatrix} n_{1x} && n_{1y} \\ n_{2x} && n_{2y} \end{vmatrix}}\space\space\space y = \frac{ \begin{vmatrix} n_{1y} && n_1 \cdot p_1 \\ n_{2y} && n_2 \cdot p_2 \end{vmatrix}}{ \begin{vmatrix} n_{1x} && n_{1y} \\ n_{2x} && n_{2y} \end{vmatrix}} $$Заметим, что если нормали $n_1, n_3$ и $n_2$ упорядочены по повороту и точка пересечения прямых $l_1$ и $l_2$ лежит в полуплоскости $h_3$, то эту полуплоскость можно выкинуть из рассмотрения: вклад в пересечение полуплоскостей она не даёт.
Домножим на него, получим
$$ n_{3x} \begin{vmatrix} n_1 \cdot p_1 && n_{1y} \\ n_2 \cdot p_2 && n_{2y} \end{vmatrix} + n_{3y} \begin{vmatrix} n_{1y} && n_1 \cdot p_1 \\ n_{2y} && n_2 \cdot p_2 \end{vmatrix} - n_3 \cdot p_3 \begin{vmatrix} n_{1x} && n_{1y} \\ n_{2x} && n_{2y} \end{vmatrix} \gt 0 $$Что эквивалентно
$$ \begin{vmatrix} n_{1x} && n_1 \cdot p_1 && n_{1y} \\ n_{2x} && n_2 \cdot p_2 && n_{2y} \\ n_{3x} && n_3 \cdot p_3 && n_{3y} \end{vmatrix} \gt 0 $$Что эквивалентно проверке поворота в двойственном пространстве точек с однородными координатами $(n_{ix}, n_i \cdot p_i, n_{iy})$ или, после деления на $n_{iy}$:
$$ \left(\frac{n_{ix}}{n_{iy}}, \frac{n_i \cdot p_i}{n_{iy}}, 1\right) $$Отметим, что точкам, упорядоченным по координате $x$ соответствуют полуплоскости, упорядоченные по повороту нормалей.
Таким образом, верхняя оболочка множества точек в двойственном пространстве даст последовательность полуплоскостей, составляющих границу пересечения.
In [7]:
%matplotlib
import numpy as np
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
from pyhull.convex_hull import ConvexHull
from math import sqrt
def process_data(Nx, r):
Ny = np.ones(len(Nx)) - Nx**2
X, Y = r[:, 0], r[:, 1]
Xd, Yd = Nx / Ny, -(Nx * X + Ny * Y) / Ny
ch = ConvexHull(np.transpose([Xd, Yd]))
ax = plt.subplot(1, 2, 1, aspect="equal")
for i, (x, y, nx, ny) in enumerate(zip(X, Y, Nx, Ny), 1):
ax.arrow(x, y, nx, ny, head_width=0.05, head_length=0.1, fc="k", ec="k")
ax.plot([x + ny * 100, x - ny * 100],
[y - nx * 100, y + nx * 100], "g", alpha=0.3)
ax.text(x, y, i)
ax.axis([0, 10, 0, 10])
ax = plt.subplot(1, 2, 2)
ax.plot(Xd, Yd, "ro")
for i, (x, y) in enumerate(zip(Xd, Yd), 1):
ax.text(x, y, i)
for v in ch.vertices:
if Xd[v[0]] > Xd[v[1]]:
ax.plot(Xd[v], Yd[v], "g")
plt.show()
process_data(np.array([-0.9, -0.5, -0.25, 0, 0.45, 0.8]),
np.array([[1, 2], [2, 4], [3, 6], [4, 5], [6, 4], [7, 3]]))
#npoints = 9
#np.random.seed(42)
#process_data(np.random.uniform(-1, 1, size=npoints),
# np.random.uniform(1, 9, size=(npoints, 2)))
Рассмотрим триангуляцию $T$ множества точек $P$. $T$ называется триангуляцией Делоне, если внутренность проведённого около любого треугольника $t \in T$ круга $c$ не содержит точек из $P$. Сопоставим каждой точке $p \in P$ точку $p' = (p_x, p_y, p^2)$ (спроецируем точки на параболоид). Из уравнения полупространства $n \cdot (p - p_0) \le 0$ и параболоида $z=x^2+y^2$ следует, что множество точек пересечения параболоида полупространством будет удовлетворять уравнению $ax+by+c(x^2+y^2)+d \le 0$, что является уравнением круга $(p - p_0)^2 \le r^2$ в случае, если нормаль к плоскости $n=(a, b, c)$ будет направлена вверх ($c \gt 0$).
Чему будет эквивалентна триангуляция Делоне?
In [5]:
%matplotlib
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter
import matplotlib.pyplot as plt
import numpy as np
fig = plt.figure()
ax = fig.gca(projection='3d', aspect="equal")
X = np.arange(1, 9, 0.25)
Y = np.arange(1, 9, 0.25)
X, Y = np.meshgrid(X, Y)
Z = X**2 + Y**2
surf = ax.plot_surface(X, Y, Z, rstride=1, cstride=1, alpha=0.2, color="g",
linewidth=0, antialiased=False)
ax.set_zlim(0, 200)
ax.axis([0, 10, 0, 10])
r = np.random.uniform(2, 8, size=(20, 2))
X, Y = r[:, 0], r[:, 1]
Z = X ** 2 + Y ** 2
ax.scatter(X, Y, np.zeros(len(r)), c="r", s=3)
ax.scatter(X, Y, Z, c="b", s=3)
ch = ConvexHull(np.transpose([X, Y, Z]))
triangles = [v for v in ch.vertices if np.linalg.det([X[v], Y[v], np.ones(3)]) > 0]
ax.plot_trisurf(X, Y, Z, alpha=0.4, triangles=triangles, color="b")
ax.plot_trisurf(X, Y, np.zeros(len(r)), alpha=0.4, triangles=triangles, color="b")
for x, y in r:
ax.plot([x, x], [y, y], "gray", alpha=0.3, zs=[0, x**2 + y**2])
plt.show()
In [4]:
%matplotlib
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter
import matplotlib.pyplot as plt
import numpy as np
from math import pi
from pyhull.convex_hull import ConvexHull
np.random.seed(321)
X, Y, Z = np.random.uniform(1, 9, size=(3, 15))
W = X**2 + Y**2 + Z**2
ch = ConvexHull(np.transpose([X, Y, Z, W]))
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d', aspect="equal")
for v in ch.vertices:
if np.linalg.det([X[v], Y[v], Z[v], np.ones(4)]) < 0:
t = v + [v[1], v[0], v[2], v[0], v[3]]
ax.plot(X[t], Y[t], "g", zs=Z[t])
ax.scatter(X, Y, Z, c="b")
ax.axis([0, 10, 0, 10])
ax.set_zlim(0, 10)
plt.show()
Аналогично введём триангуляцию Делоне на сфере. Рассмотрим множество точек $P$ на сфере $S$. Скажем, что триангуляция множества точек $P$ является триангуляцией Делоне, если внутренность описанного около каждого треугольника $t$ круга $c$ не содержит точек из $P$. Существование триангуляции Делоне следует из того, что окружность — это пересечение плоскости, а круг — это пересечение полупространства со сферой $S$.
In [3]:
%matplotlib
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter
import matplotlib.pyplot as plt
import numpy as np
from math import pi
from pyhull.convex_hull import ConvexHull
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d', aspect="equal")
a, b, c, d = 0.5, 0.5, 0.5, 5
x, y = np.transpose([[-10, -10], [-10, 10], [10, 10], [10, -10]])
z = (d - a * x - b * y) / c
ax.plot_trisurf(x, y, z, triangles=[[0, 1, 2], [0, 2, 3]], color='r', alpha=0.5, linewidth=0)
u = np.linspace(0, 2 * np.pi, 100)
v = np.linspace(0, np.pi, 100)
x = 10 * np.outer(np.cos(u), np.sin(v))
y = 10 * np.outer(np.sin(u), np.sin(v))
z = 10 * np.outer(np.ones(np.size(u)), np.cos(v))
ax.plot_surface(x, y, z, alpha=0.2, rstride=4, cstride=4, color='b', linewidth=0.5)
ax.axis([-11, 11, -11, 11])
ax.set_zlim(-11, 11)
plt.show()
Поэтому условие "точка не лежит внутри круга $c$" переформулируется как "точка не лежит в полупространстве, пересекающем сферу $S$". Чем является триангуляция Делоне в этом случае?
In [1]:
%matplotlib
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter
import matplotlib.pyplot as plt
import numpy as np
from math import pi
from pyhull.convex_hull import ConvexHull
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d', aspect="equal")
u = np.linspace(0, 2 * np.pi, 100)
v = np.linspace(0, np.pi, 100)
x = 10 * np.outer(np.cos(u), np.sin(v))
y = 10 * np.outer(np.sin(u), np.sin(v))
z = 10 * np.outer(np.ones(np.size(u)), np.cos(v))
ax.plot_surface(x, y, z, alpha=0.2, rstride=4, cstride=4, color='b', linewidth=0)
u, t = np.random.uniform(-1, 1, 120), np.random.uniform(0, 2*pi, 120)
ones = np.ones(len(u))
X = 10 * np.sqrt(ones - u ** 2) * np.cos(t)
Y = 10 * np.sqrt(ones - u ** 2) * np.sin(t)
Z = 10 * u
ax.scatter(X, Y, Z, c="b", s=3)
ch = ConvexHull(np.transpose([X, Y, Z]))
ax.plot_trisurf(X, Y, Z, alpha=0.9, triangles=ch.vertices)
ax.axis([-11, 11, -11, 11])
ax.set_zlim(-11, 11)
plt.show()