我正在模拟点在倾斜平面上的投影。我想保持与缺少 scipy.interpolate 等效项的 cupy 的兼容性。
一般来说,给定一个点p = [x,y,z]^T
和一个nxn
旋转矩阵m
,输出坐标中的点p'
是p' = m p
。换句话说,m
是从输入到输出空间的前向映射。
在 numpy 中,您可以将其实现为pp = m.dot(p)
or pp = m @ p
。同时对于许多点,使用 tensordot 或 einsum 进行投影要快得多:
# m = 3x3,
# points=16384x3
out = np.tensordot(m, points, axes=(1,1))
m
和的样本points
是
m = np.array([[1.0, -0.0, 0.0, 0.0],
[0.0, 0.5000000000000001, -0.8660254037844386, 0.0],
[-0.0, 0.8660254037844386, 0.5000000000000001, 0.0],
[0.0, 0.0, 0.0, 1.0]])
m = m[:3,:3]
N = 128
x = y = np.arange(N, dtype=np.float64)
x -= N//2
x, y = np.meshgrid(x,y)
z = np.ones_like(x)
x2 = x.ravel()
y2 = y.ravel()
z2 = z.ravel()
points = np.stack((x2,y2,z2), axis=1)
X, Y, Z = points
X = X.reshape(x.shape)
Y = Y.reshape(y.shape)
应用转换后,输入和输出 x,y 点如下所示:
使用 map_coordinates 会导致在输入网格上查找输出域中的点,从而有效地执行逆映射。有没有办法使用 map_coordinates 执行前向映射?不正确的实现看起来像
def regularize(xy, XY, z):
"""Regularize the coordinates XY relative to the frame xy.
This function is used in conjunction with rotate to project
surface figure errors onto tilted planes or other geometries.
Parameters
----------
xy : `numpy.ndarray`
ndarray of dimension (2, m, n)
with [x, y] on the first dimension
represents the input coordinates
implicitly rectilinear
XY : `numpy.ndarray`
ndarray of dimension (2, m, n)
with [x, y] on the first dimension
represents the input coordinates
not necessarily rectilinear
Returns
-------
Z : `numpy.ndarray`
z which exists on the grid XY, looked up at the points xy
"""
xy = np.array(xy)
XY = np.array(XY)
# map coordinates says [0,0] is the upper left corner
# need to adjust XYZ by xyz origin and sample spacing
# d = delta; o = origin
x,y = xy
ox = x[0,0]
oy = y[0,0]
dx = x[0,1] - ox
dy = y[1,0] - oy
XY2 = XY.copy()
X, Y = XY2
X -= ox
Y -= oy
X /= dx
Y /= dy
# ::-1 = reverse X,Y
# ... = leave other axes as-is
XY2 = XY2[::-1,...]
return ndimage.map_coordinates(z, XY2)
不基于 ndimage 的正确实现将使用 interp2 或来自 scipy.interpolate 的 griddata