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
|
# -*- coding: utf-8 -*-
# Copyright (c) Vispy Development Team. All Rights Reserved.
# Distributed under the (new) BSD License. See LICENSE.txt for more info.
"""Miscellaneous functions"""
import numpy as np
def _cross_2d(x, y):
"""Compute the z-component of the cross product of (arrays of) 2D vectors.
This is meant to replicate the 2D functionality of np.cross(), which is
deprecated in numpy 2.0.
x and y must have broadcastable shapes, with the last dimension being 2.
Parameters
----------
x : array
Input array 1, shape (..., 2).
y : array
Input array 2, shape (..., 2).
Returns
-------
z : array
z-component of cross products of x and y.
See: https://github.com/numpy/numpy/issues/26620
"""
if x.shape[-1] != 2 or y.shape[-1] != 2:
raise ValueError("Input arrays must have shape (..., 2)")
return x[..., 0] * y[..., 1] - x[..., 1] * y[..., 0]
def _fast_cross_3d(x, y):
"""Compute cross product between list of 3D vectors
Much faster than np.cross() when the number of cross products
becomes large (>500). This is because np.cross() methods become
less memory efficient at this stage.
Parameters
----------
x : array
Input array 1.
y : array
Input array 2.
Returns
-------
z : array
Cross product of x and y.
Notes
-----
x and y must both be 2D row vectors. One must have length 1, or both
lengths must match.
"""
assert x.ndim == 2
assert y.ndim == 2
assert x.shape[1] == 3
assert y.shape[1] == 3
assert (x.shape[0] == 1 or y.shape[0] == 1) or x.shape[0] == y.shape[0]
if max([x.shape[0], y.shape[0]]) >= 500:
return np.c_[x[:, 1] * y[:, 2] - x[:, 2] * y[:, 1],
x[:, 2] * y[:, 0] - x[:, 0] * y[:, 2],
x[:, 0] * y[:, 1] - x[:, 1] * y[:, 0]]
else:
return np.cross(x, y)
###############################################################################
# These fast normal calculation routines are adapted from mne-python
def _calculate_normals(rr, tris):
"""Efficiently compute vertex normals for triangulated surface"""
# ensure highest precision for our summation/vectorization "trick"
rr = rr.astype(np.float64)
# first, compute triangle normals
r1 = rr[tris[:, 0], :]
r2 = rr[tris[:, 1], :]
r3 = rr[tris[:, 2], :]
tri_nn = _fast_cross_3d((r2 - r1), (r3 - r1))
# Triangle normals and areas
size = np.sqrt(np.sum(tri_nn * tri_nn, axis=1))
size[size == 0] = 1.0 # prevent ugly divide-by-zero
tri_nn /= size[:, np.newaxis]
npts = len(rr)
# the following code replaces this, but is faster (vectorized):
#
# for p, verts in enumerate(tris):
# nn[verts, :] += tri_nn[p, :]
#
nn = np.zeros((npts, 3))
for verts in tris.T: # note this only loops 3x (number of verts per tri)
for idx in range(3): # x, y, z
nn[:, idx] += np.bincount(verts.astype(np.int32),
tri_nn[:, idx], minlength=npts)
size = np.sqrt(np.sum(nn * nn, axis=1))
size[size == 0] = 1.0 # prevent ugly divide-by-zero
nn /= size[:, np.newaxis]
return nn
def resize(image, shape, kind='linear'):
"""Resize an image
Parameters
----------
image : ndarray
Array of shape (N, M, ...).
shape : tuple
2-element shape.
kind : str
Interpolation, either "linear" or "nearest".
Returns
-------
scaled_image : ndarray
New image, will have dtype np.float64.
"""
image = np.array(image, float)
shape = np.array(shape, int)
if shape.ndim != 1 or shape.size != 2:
raise ValueError('shape must have two elements')
if image.ndim < 2:
raise ValueError('image must have two dimensions')
if not isinstance(kind, str) or kind not in ('nearest', 'linear'):
raise ValueError('mode must be "nearest" or "linear"')
r = np.linspace(0, image.shape[0] - 1, shape[0])
c = np.linspace(0, image.shape[1] - 1, shape[1])
if kind == 'linear':
r_0 = np.floor(r).astype(int)
c_0 = np.floor(c).astype(int)
r_1 = r_0 + 1
c_1 = c_0 + 1
top = (r_1 - r)[:, np.newaxis]
bot = (r - r_0)[:, np.newaxis]
lef = (c - c_0)[np.newaxis, :]
rig = (c_1 - c)[np.newaxis, :]
c_1 = np.minimum(c_1, image.shape[1] - 1)
r_1 = np.minimum(r_1, image.shape[0] - 1)
for arr in (top, bot, lef, rig):
arr.shape = arr.shape + (1,) * (image.ndim - 2)
out = top * rig * image[r_0][:, c_0, ...]
out += bot * rig * image[r_1][:, c_0, ...]
out += top * lef * image[r_0][:, c_1, ...]
out += bot * lef * image[r_1][:, c_1, ...]
else: # kind == 'nearest'
r = np.round(r).astype(int)
c = np.round(c).astype(int)
out = image[r][:, c, ...]
return out
|
