0

我正在模拟点在倾斜平面上的投影。我想保持与缺少 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

4

0 回答 0