更改 3D numpy 数组的基础（分数到笛卡尔坐标）

Question

这是我想在 3D 中实现的 2D 示例：

我有一个值数组，A，s.t。 A.shape=(n,m)，例如

>>> A = [[1,    2],
...      [3,    4]]

其索引与沿（任意）基向量的等距步长成比例，例如

>>> v1 = [1,0]
>>> v2 = [cos(pi/4),sin(pi/4)] # [0,1] rotated 45 degrees

对于这个例子，我想要一个应用这个基础的函数来获取

>>> apply_basis2D(A,v1,v2)
[[np.nan,1,    2],
 [3,     4,    np.nan]]

（那么对于 3D 版本，我想要 apply_basis3D(A,v1,v2,v3)），其中 A.shape=(n,m,l))

我认为这可以通过仿射变换来完成，但我不确定如何。这是我能找到的最接近的实现（仅限二维），使用scikit-image；

提前致谢！

Answer 1

完成！似乎工作得很好，但我欢迎批评：

import numpy as np
from scipy.spatial import Delaunay
from scipy.interpolate import interpn
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
%matplotlib notebook

def cartesian_basis2d(A,v1,v2,longest_side=None):
    """convert 2d array in basis v1,v2 to cartesian basis

    Properties
    ----------
    A : array((N,M))
        values in original basis 
    v1 : array((2,))
    v2 : array((2,))
    longest_side : int
        longest side (in terms of indexes) of new array

    Returns
    -------
    B : array((P,Q))
        where P,Q >= longest_side

    """
    longest_side = max(A.shape) if longest_side is None else longest_side

    # assumed
    origin = [0,0]

    # convert to numpy arrays
    origin = np.asarray(origin)
    v1 = np.asarray(v1)
    v2 = np.asarray(v2)

    # pre-compute basis transformation matrix
    M_inv = np.linalg.inv(np.transpose([v1,v2]))

    # only works rigth if transposed before and after?
    A = np.array(A).T
    # add bounding rows/columns for interpolation
    A = np.concatenate((np.array(A[:,0],ndmin=2).T,A,np.array(A[:,-1],ndmin=2).T),axis=1)
    A = np.concatenate((np.array(A[0],ndmin=2),A,np.array(A[-1],ndmin=2)),axis=0)
    # create axes
    axes=[]
    for i,v in enumerate([v1,v2]):
        step = 1./(A.shape[i]-2)
        ax = np.linspace(0,1+step,A.shape[i]) - step/2.
        axes.append(ax)

    # get bounding box and compute it volume and extents
    bbox_pts=np.asarray([origin,v1,v1+v2,v2])
    hull = Delaunay(bbox_pts)

    bbox_x, bbox_y = bbox_pts.T
    new_bounds = bbox_x.min(),bbox_x.max(),bbox_y.min(),bbox_y.max() #l,r,bottom,top
    min_bound, max_bound = min(bbox_x.min(),bbox_y.min()), max(bbox_x.max(),bbox_y.max())

    # compute new array size
    x_length = abs(new_bounds[0]-new_bounds[1])
    y_length = abs(new_bounds[2]-new_bounds[3])
    if x_length>y_length:
        xlen = longest_side
        ylen = int(longest_side*y_length/float(x_length))
    else:
        ylen = longest_side
        xlen = int(longest_side*x_length/float(y_length))

    # compute new array values
    new_array = np.full((xlen,ylen),np.nan)
    xidx, yidx = np.meshgrid(range(new_array.shape[0]),range(new_array.shape[1]))
    xidx=xidx.flatten()
    yidx=yidx.flatten()
    xyidx = np.concatenate((np.array(xidx,ndmin=2).T,np.array(yidx,ndmin=2).T),axis=1)
    xy = min_bound+(xyidx.astype(float)*abs(min_bound-max_bound)/longest_side)
    # find point is inside bounding box 
    inside_mask = hull.find_simplex(xy)>=0
    uv = np.einsum('...jk,...k->...j',M_inv,xy[inside_mask])
    new_array[xyidx[inside_mask][:,0],xyidx[inside_mask][:,1]] = interpn(axes,A,uv,bounds_error=True,method='nearest')
    new_array = new_array.T  
    return new_array

A = np.array(
    [[1,2,3],
     [4,5,6],
     [7,8,9]])
v1 = [2,0]
v2 = [np.cos(np.pi/4),np.sin(np.pi/4)]

new_array = cartesian_basis2d(A,v1,v2,100)

plt.imshow(new_array,origin='lower');

def cartesian_basis3d(A,v1,v2,v3,longest_side=None):
    """convert 3d array in basis v1,v2,v3 to cartesian basis

    Properties
    ----------
    A : array((N,M))
        values in original basis 
    v1 : array((2,))
    v2 : array((2,))
    v3 : array((2,))
    longest_side : int
        longest side (in terms of indexes) of new array

    Returns
    -------
    B : array((P,Q))
        where P,Q >= longest_side

    """
    longest_side = max(A.shape) if longest_side is None else longest_side

    # assumed
    origin = [0,0,0]

    # convert to numpy arrays
    origin = np.asarray(origin)
    v1 = np.asarray(v1)
    v2 = np.asarray(v2)
    v3 = np.asarray(v3)

    # pre-compute basis transformation matrix
    M_inv = np.linalg.inv(np.transpose([v1,v2,v3]))

    # only works rigth if transposed before and after?
    A = np.array(A).T
    # add bounding layers for interpolation
    A = np.concatenate((np.array(A[0],ndmin=3),A,np.array(A[-1],ndmin=3)),axis=0)
    start = np.transpose(np.array(A[:,:,0],ndmin=3),axes=[1,2,0])
    end = np.transpose(np.array(A[:,:,-1],ndmin=3),axes=[1,2,0])
    A = np.concatenate((start,A,end),axis=2)
    start = np.transpose(np.array(A[:,0,:],ndmin=3),axes=[1,0,2])
    end = np.transpose(np.array(A[:,-1,:],ndmin=3),axes=[1,0,2])
    A = np.concatenate((start,A,end),axis=1)

    # create axes
    axes=[]
    for i,v in enumerate([v1,v2,v3]):
        step = 1./(A.shape[i]-2)
        ax = np.linspace(0,1+step,A.shape[i]) - step/2.
        axes.append(ax)

    # get bounding box and compute it volume and extents
    bbox_pts=np.asarray([origin,v1,v2,v3,v1+v2,v1+v3,v1+v2+v3,v2+v3])
    hull = Delaunay(bbox_pts)

    bbox_x, bbox_y, bbox_z = bbox_pts.T
    new_bounds = bbox_x.min(),bbox_x.max(),bbox_y.min(),bbox_y.max(),bbox_z.min(),bbox_z.max() #l,r,bottom,top
    min_bound, max_bound = min(bbox_x.min(),bbox_y.min(),bbox_z.min()), max(bbox_x.max(),bbox_y.max(),bbox_z.min())

    # compute new array size
    x_length = abs(new_bounds[0]-new_bounds[1])
    y_length = abs(new_bounds[2]-new_bounds[3])
    z_length = abs(new_bounds[4]-new_bounds[5])
    if x_length == max([x_length,y_length,z_length]):
        xlen = longest_side
        ylen = int(longest_side*y_length/float(x_length))
        zlen = int(longest_side*z_length/float(x_length))
    elif y_length == max([x_length,y_length,z_length]):
        ylen = longest_side
        xlen = int(longest_side*x_length/float(y_length))
        zlen = int(longest_side*z_length/float(y_length))
    else:
        zlen = longest_side
        xlen = int(longest_side*x_length/float(z_length))
        ylen = int(longest_side*y_length/float(z_length))

    # compute new array values
    new_array = np.full((xlen,ylen,zlen),np.nan)
    xidx, yidx, zidx = np.meshgrid(range(new_array.shape[0]),range(new_array.shape[1]),range(new_array.shape[2]))
    xidx=xidx.flatten()
    yidx=yidx.flatten()
    zidx=zidx.flatten()
    xyzidx = np.concatenate((np.array(xidx,ndmin=2).T,np.array(yidx,ndmin=2).T,np.array(zidx,ndmin=2).T),axis=1)
    xyz = min_bound+(xyzidx.astype(float)*abs(min_bound-max_bound)/longest_side)
    # find point is inside bounding box 
    inside_mask = hull.find_simplex(xyz)>=0
    uvw = np.einsum('...jk,...k->...j',M_inv,xyz[inside_mask])
    new_array[xyzidx[inside_mask][:,0],xyzidx[inside_mask][:,1],xyzidx[inside_mask][:,2]] = interpn(axes,A,uvw,bounds_error=True,method='nearest')
    new_array = new_array.T  
    return new_array

A = np.array(
    [[[1,1],[2,2]],
     [[3,3],[4,4]]])

v1 = [2,0,0]
v2 = [np.cos(np.pi/4),np.sin(np.pi/4),0]
v3 = [0,np.cos(np.pi/4),np.sin(np.pi/4)]

new_array = cartesian_basis3d(A,v1,v2,v3,100)

xs,ys,zs = new_array.nonzero()
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
pcm = ax.scatter(xs, ys, zs, c=new_array[xs,ys,zs],cmap='jet')
plt.show()