-
Notifications
You must be signed in to change notification settings - Fork 3
/
hermite.py
190 lines (166 loc) · 7.84 KB
/
hermite.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
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
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
"""
Functions that calculate Hermite curves and surfaces.
"""
import numpy as np
# Function to visualize a Hermite Curve
def HermiteCurve(p, num):
points = np.linspace(0, 1, num)
m = [[2, -2, 1, 1], [-3, 3, -2, -1], [0, 0, 1, 0], [1, 0, 0, 0]]
curve = np.array((points**3, points**2, points, points**0)).transpose() @ m @ p
return curve
# Function to visualize a Hermite Surface
def HermiteSurface(p, num_u, num_v):
points_u = np.linspace(0, 1, num_u)
points_v = np.linspace(0, 1, num_v)
m = np.array([[2, -2, 1, 1], [-3, 3, -2, -1], [0, 0, 1, 0], [1, 0, 0, 0]])
surface = (
np.array((points_u**3, points_u**2, points_u, points_u**0)).transpose()
@ m @ p @ m.transpose()
@ np.array((points_v**3, points_v**2, points_v, points_v**0))
)
return surface
# Function to calculate continuity between Hermite Curves
def HermiteCurveContinuity(p1, p2):
"""Return the continuity of two Hermite curves as a tuple (str, int), or return None if no continuity exists."""
M = np.array([[0, 0, 0, 0], [0, 0, 0, 0], [12, -12, 6, 6], [-6, 6, -4, -2]])
p1u0 = np.dot(np.dot(np.array([0, 0, 0, 1]), M), p1[:, :, 0].transpose())
p1u1 = np.dot(np.dot(np.array([1, 1, 1, 1]), M), p1[:, :, 0].transpose())
p2u0 = np.dot(np.dot(np.array([0, 0, 0, 1]), M), p2[:, :, 0].transpose())
p2u1 = np.dot(np.dot(np.array([1, 1, 1, 1]), M), p2[:, :, 0].transpose())
if (p1[:, 1, :] == p2[:, 0, :]).all() or (p2[:, 1, :] == p1[:, 0, :]).all():
if (p1[:, 3, :] == p2[:, 2, :]).all() or (p2[:, 3, :] == p1[:, 2, :]).all():
if (p1u1 == p2u0).all() or (p1u0 == p2u1).all():
return ('C', 2)
return ('C', 1)
if 0 in p1[:, 2, :] or 0 in p1[:, 3, :]:
return ('CG', 0)
if (p1[:, 3, :] / p2[:, 2, :] == (p1[:, 3, :] / p2[:, 2, :])[0]).all() or (p2[:, 3, :] / p1[:, 2, :] == (p2[:, 3, :] / p1[:, 2, :])[0]).all():
return ('G', 1)
return ('CG', 0)
return None
# Function to calculate continuity between Hermite Curves
def HermiteSurfaceContinuity(p1, p2):
"""Return the continuity of two Hermite surfaces as a tuple (str, int), or return None if no continuity exists."""
p1sides = np.zeros((4, len(p1), len(p1[0])))
for i in range(len(p1)):
p1sides[0][i] = np.array([p1[i][0][0], p1[i][0][1], p1[i][2][0], p1[i][2][1]])
p1sides[1][i] = np.array([p1[i][0][0], p1[i][1][0], p1[i][0][2], p1[i][1][2]])
p1sides[2][i] = np.array([p1[i][1][0], p1[i][1][1], p1[i][3][0], p1[i][3][1]])
p1sides[3][i] = np.array([p1[i][0][1], p1[i][1][1], p1[i][0][3], p1[i][1][3]])
p2sides = np.zeros((4, len(p2), len(p2[0])))
for i in range(len(p1)):
p2sides[0][i] = np.array([p2[i][0][0], p2[i][0][1], p2[i][2][0], p2[i][2][1]])
p2sides[1][i] = np.array([p2[i][0][0], p2[i][1][0], p2[i][0][2], p2[i][1][2]])
p2sides[2][i] = np.array([p2[i][1][0], p2[i][1][1], p2[i][3][0], p2[i][3][1]])
p2sides[3][i] = np.array([p2[i][0][1], p2[i][1][1], p2[i][0][3], p2[i][1][3]])
for i in range(len(p1sides)):
for j in range(len(p2sides)):
if (p1sides[i].T[0] == p2sides[j].T[0]).all() and (p1sides[i].T[1] == p2sides[j].T[1]).all():
p1 = p1sides[i]
p2 = p2sides[j]
if (p1sides[i].T[2] == p2sides[j].T[2]).all() and (p1sides[i].T[3] == p2sides[j].T[3]).all():
M = np.array([[0, 0, 0, 0], [0, 0, 0, 0], [12, -12, 6, 6], [-6, 6, -4, -2]])
p1u0 = np.dot(np.dot(np.array([0, 0, 0, 1]), M), p1.transpose())
p1u1 = np.dot(np.dot(np.array([1, 1, 1, 1]), M), p1.transpose())
p2u0 = np.dot(np.dot(np.array([0, 0, 0, 1]), M), p2.transpose())
p2u1 = np.dot(np.dot(np.array([1, 1, 1, 1]), M), p2.transpose())
if (p1u1 == p2u0).all() or (p1u0 == p2u1).all():
return ('C', 2)
return ('C', 1)
if (p1[:, 3, :] / p2[:, 2, :] == (p1[:, 3, :] / p2[:, 2, :])[0]).all() or (p2[:, 3, :] / p1[:, 2, :] == (p2[:, 3, :] / p1[:, 2, :])[0]).all():
return ('G', 1)
return ('CG', 0)
return None
# # Initialize control points
# cp = np.array([[-2, 3, 5, 11],
# [-2, 9, 4, 10]])
# cpp = np.array([[[1, 3, 6, 8],
# [1, 3, 6, 8],
# [1, 3, 6, 8],
# [1, 3, 6, 8]],
# [[20, 21, 22, 23],
# [17, 17, 17, 17],
# [14, 14, 14, 14],
# [11, 11, 11, 11]],
# [[2, 5, 4, 3],
# [2, 6, 5, 5],
# [2, 6, 5, 4],
# [2, 3, 4, 3]]])
# Create Curve and Surface with 100 discrete points
# curve = BezierCurve(cp, 100)
# surface = BezierSurface(cpp, 100, 100)
# Plot curve and control points
# figure = plt.figure()
# ax = figure.add_subplot(projection='3d')
# ax.set_xlabel('X-axis')
# ax.set_ylabel('Y-axis')
# ax.set_zlabel('Z-axis')
# ax.set_title('Bézier Curve')
# for i in range(len(cp[0])):
# ax.scatter(cp[0][i], cp[1][i])
# ax.plot(curve[0], curve[1])
# Plot surface and control points
# figure = plt.figure()
# ax = figure.add_subplot(projection='3d')
# ax.set_xlabel('X-axis')
# ax.set_ylabel('Y-axis')
# ax.set_zlabel('Z-axis')
# ax.set_title('Bézier Surface')
# for i in range(len(cpp[0])):
# for j in range(len(cpp[0])):
# ax.scatter(cpp[0][i][j], cpp[1][i][j], cpp[2][i][j])
# ax.plot_surface(surface[0], surface[1], surface[2])
# Starting and Ending points and tangents
# p0, p1, p0u, p1u for x, y and z
# p1 = np.array([[1, 10, 1, 2], [20, 22, 2, 2], [3, 2, 0, -1]])
# p2 = np.array([[10, 20, 4, 2], [22, 24, 4, 2], [2, 1, -2, -1]])
#
# cp_1 = np.array([[[1, 5, 0], [3, 8, 0], [3, 3, 0], [1.9286, -1.2321, 0]]]).transpose()
# cp_2 = np.array([[[3, 8, 0], [6, 4, 0], [1.9286, -1.2321, 0], [4.2857, -1.0714, 0]]]).transpose()
# Create Hermite Curve with 100 discrete points
# curve = HermiteCurve(p1, 100)
# print(HermiteCurveContinuity(cp_1, cp_2))
# Plot curve and starting and ending points
# figure = plt.figure()
# ax = figure.add_subplot(projection='3d')
# ax.set_xlabel('X-axis')
# ax.set_ylabel('Y-axis')
# ax.set_zlabel('Z-axis')
# ax.set_title('Hermite Curve')
# for i in range(2):
# ax.scatter(p[0][i], p[1][i], p[2][i])
# ax.plot(curve[0], curve[1], curve[2])
# plt.show()
# import matplotlib.pyplot as plt
# from mpl_toolkits import mplot3d
# # Starting and Ending points and tangents
# # p00 p01 p00w p01w
# # p10 p11 p10w p11w
# # p00u p01u p00uw p01uw
# # p10u p11u p10uw p11uw
# p = [[[-50, -50, 0, 0], [50, 50, 0, 0], [50, 5, 0, .5], [5, 5, .5, 0]],
# [[0, -50, 20, -5], [-50, 0, 5, -5], [50, 5, 0, .5], [-5, -5, -.5, 0]],
# [[50, -50, -20, -5], [50, -50, -5, -5], [0, 0, 0, .5], [0, 0, -.5, 0]]]
# # Create Hermite Surface with 100 discrete points for u and v
# surface = HermiteSurface(p, 100, 100)
# # Plot surface and corner points
# figure = plt.figure()
# ax = figure.add_subplot(projection='3d')
# ax.set_xlabel('X-axis')
# ax.set_ylabel('Y-axis')
# ax.set_zlabel('Z-axis')
# ax.set_title('Hermite Surface')
# for i in range(2):
# for j in range(2):
# ax.scatter(p[0][i][j], p[1][i][j], p[2][i][j])
# ax.plot_surface(surface[0], surface[1], surface[2])
# plt.show()
# cp_1 = np.array([[[0, 0, 0], [0, 10, 0], [0, 1, 0], [0, 1, 0]],
# [[10, 0, 0], [10, 10, 0], [0, 1, 0], [0, 1, 0]],
# [[1, 0, 0], [1, 0, 0], [0, 0, 1], [0, 0, 1]],
# [[1, 0, 0], [1, 0, 0], [0, 0, 1], [0, 0, 1]]]).transpose((2, 0, 1))
# cp_2 = np.array([[[10, 0, 0], [10, 10, 0], [0, 1, 0], [0, 1, 0]],
# [[20, 0, 0], [20, 10, 0], [0, 1, 0], [0, 1, 0]],
# [[1, 0, 0], [1, 0, 0], [0, 0, 1], [0, 0, 1]],
# [[1, 0, 0], [1, 0, 0], [0, 0, 1], [0, 0, 1]]]).transpose((2, 0, 1))
# print(HermiteSurfaceContinuity(cp_1, cp_2))