-
Notifications
You must be signed in to change notification settings - Fork 184
/
functions.py
254 lines (189 loc) · 6.46 KB
/
functions.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
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
from __future__ import annotations
import numpy as np
from numpy import cos, float64, ndarray, sin
RAD2DEG = 180.0 / np.pi
DEG2RAD = 1 / RAD2DEG
def sign_shape(pts: ndarray) -> float64:
pts2 = np.roll(pts, 1, axis=0)
dx = pts2[:, 0] - pts[:, 0]
y = pts2[:, 1] + pts[:, 1]
return np.sign((dx * y).sum())
def area(pts: ndarray) -> float64:
"""Returns the area."""
pts2 = np.roll(pts, 1, axis=0)
dx = pts2[:, 0] - pts[:, 0]
y = pts2[:, 1] + pts[:, 1]
return (dx * y).sum() / 2
def manhattan_direction(p0, p1, tol=1e-5):
"""Returns manhattan direction between 2 points."""
dp = p1 - p0
dx, dy = dp[0], dp[1]
if abs(dx) < tol:
sx = 0
elif dx > 0:
sx = 1
else:
sx = -1
if abs(dy) < tol:
sy = 0
elif dy > 0:
sy = 1
else:
sy = -1
return np.array((sx, sy))
def remove_flat_angles(points: ndarray) -> ndarray:
a = angles_deg(np.vstack(points))
da = a - np.roll(a, 1)
da = np.mod(np.round(da, 3), 180)
# To make sure we do not remove points at the edges
da[0] = 1
da[-1] = 1
to_rm = list(np.where(np.abs(da[:-1]) < 1e-9)[0])
if isinstance(points, list):
while to_rm:
i = to_rm.pop()
points.pop(i)
else:
points = points[da != 0]
return points
def remove_identicals(
pts: ndarray, grids_per_unit: int = 1000, closed: bool = True
) -> ndarray:
if len(pts) > 1:
identicals = np.prod(abs(pts - np.roll(pts, -1, 0)) < 0.5 / grids_per_unit, 1)
if not closed:
identicals[-1] = False
pts = np.delete(pts, identicals.nonzero()[0], 0)
return pts
def centered_diff(a: ndarray) -> ndarray:
d = (np.roll(a, -1, axis=0) - np.roll(a, 1, axis=0)) / 2
return d[1:-1]
def centered_diff2(a: ndarray) -> ndarray:
d = (np.roll(a, -1, axis=0) - a) - (a - np.roll(a, 1, axis=0))
return d[1:-1]
def curvature(points: ndarray, t: ndarray) -> ndarray:
"""Args are the points and the tangents at each point.
points : numpy.array shape (n, 2)
t: numpy.array of size n
Return:
The curvature at each point.
Computes the curvature at every point excluding the first and last point.
For a planar curve parametrized as P(t) = (x(t), y(t)), the curvature is given
by (x' y'' - x'' y' ) / (x' **2 + y' **2)**(3/2)
"""
# Use centered difference for derivative
dt = centered_diff(t)
dp = centered_diff(points)
dp2 = centered_diff2(points)
dx = dp[:, 0] / dt
dy = dp[:, 1] / dt
dx2 = dp2[:, 0] / dt**2
dy2 = dp2[:, 1] / dt**2
return (dx * dy2 - dx2 * dy) / (dx**2 + dy**2) ** (3 / 2)
def radius_of_curvature(points, t):
return 1 / curvature(points, t)
def path_length(points: ndarray) -> float64:
"""Returns: The path length.
Args:
points: With shape (N, 2) representing N points with coordinates x, y.
"""
dpts = points[1:, :] - points[:-1, :]
_d = dpts**2
return np.sum(np.sqrt(_d[:, 0] + _d[:, 1]))
def snap_angle(a: float64) -> int:
"""Returns angle snapped along manhattan angle (0, 90, 180, 270).
a: angle in deg
Return angle snapped along manhattan angle
"""
a = a % 360
if -45 < a < 45:
return 0
elif 45 < a < 135:
return 90
elif 135 < a < 225:
return 180
elif 225 < a < 315:
return 270
else:
return 0
def angles_rad(pts: ndarray) -> ndarray:
"""Returns the angles (radians) of the connection between each point and the next."""
_pts = np.roll(pts, -1, 0)
return np.arctan2(_pts[:, 1] - pts[:, 1], _pts[:, 0] - pts[:, 0])
def angles_deg(pts: ndarray) -> ndarray:
"""Returns the angles (degrees) of the connection between each point and the next."""
return angles_rad(pts) * RAD2DEG
def extrude_path(
points: ndarray,
width: float,
with_manhattan_facing_angles: bool = True,
spike_length: float64 | int | float = 0,
start_angle: int | None = None,
end_angle: int | None = None,
grid: float | None = None,
) -> ndarray:
"""Deprecated. Use gdsfactory.path.Path.extrude() instead.
Extrude a path of `width` along a curve defined by `points`.
Args:
points: numpy 2D array of shape (N, 2).
width: of the path to extrude.
with_manhattan_facing_angles: snaps to manhattan angles.
spike_length: in um.
start_angle: in degrees.
end_angle: in degrees.
grid: in um.
Returns:
numpy 2D array of shape (2*N, 2).
"""
from gdsfactory.pdk import get_grid_size
grid = grid or get_grid_size()
if isinstance(points, list):
points = np.stack([(p[0], p[1]) for p in points], axis=0)
a = angles_deg(points)
if with_manhattan_facing_angles:
_start_angle = snap_angle(a[0] + 180)
_end_angle = snap_angle(a[-2])
else:
_start_angle = a[0] + 180
_end_angle = a[-2]
start_angle = start_angle if start_angle is not None else _start_angle
end_angle = end_angle if end_angle is not None else _end_angle
a2 = angles_rad(points) * 0.5
a1 = np.roll(a2, 1)
a2[-1] = end_angle * DEG2RAD - a2[-2]
a1[0] = start_angle * DEG2RAD - a1[1]
a_plus = a2 + a1
cos_a_min = np.cos(a2 - a1)
offsets = np.column_stack((-sin(a_plus) / cos_a_min, cos(a_plus) / cos_a_min)) * (
0.5 * width
)
points_back = np.flipud(points - offsets)
if spike_length != 0:
d = spike_length
a_start = start_angle * DEG2RAD
a_end = end_angle * DEG2RAD
p_start_spike = points[0] + d * np.array([[cos(a_start), sin(a_start)]])
p_end_spike = points[-1] + d * np.array([[cos(a_end), sin(a_end)]])
pts = np.vstack((p_start_spike, points + offsets, p_end_spike, points_back))
else:
pts = np.vstack((points + offsets, points_back))
pts = np.round(pts / grid) * grid
return pts
def polygon_grow(polygon: ndarray, offset: float) -> ndarray:
"""Returns a grown closed shaped polygon by an offset."""
s = remove_identicals(polygon)
s = remove_flat_angles(s)
s = np.vstack([s, s[0]])
if len(s) <= 1:
return s
# Make sure the shape is oriented in the correct direction for scaling
ss = sign_shape(s)
offset *= -ss
a2 = angles_rad(s) * 0.5
a1 = np.roll(a2, 1)
a2[-1] = a2[0]
a1[0] = a1[-1]
a = a2 + a1
c_minus = cos(a2 - a1)
offsets = np.column_stack((-sin(a) / c_minus, cos(a) / c_minus)) * offset
return s + offsets