更改 3D numpy 数组的基础（分数到笛卡尔坐标）
Change basis of 3D numpy array (fractional to cartesian coordinates)
这是我想要在 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))
我认为这可以通过仿射变换来完成，但我不太确定如何实现。这是我能找到的最接近的实现（仅限 2D），使用 scikit-image；
提前致谢！
完成！貌似还不错，欢迎批评指正：
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()
这是我想要在 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))
我认为这可以通过仿射变换来完成，但我不太确定如何实现。这是我能找到的最接近的实现（仅限 2D），使用 scikit-image；
提前致谢！
完成！貌似还不错，欢迎批评指正：
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()