你可以使用scipy.ndimage.interpolation.rotate来旋转你的3d数组到你想要的任何Angular （它使用样条插值），然后你可以从它中取出一个切片。

def extract_slice（数据，三元组）：
““”
算法：
1.查找平面与数据标注栏边得交点
2.对于这些pts，找到面向轴的b-box
3.求出从R2到R3的“反向”反式（A，T），如X“= AX + T
使用（0，0），（0，h），（w，0）这三个易于计算的参数
4.对2D（w，h）图像中的每个值使用变换（使用三线性插值）
““”
我将在几个月内适当地发布代码，我相信作为这个项目的一部分。

我正在处理与OP类似的任务，因此我提出了基于numpy（而不是scipy）的代码，以从体积中提取任何给定切片，给定平面上任何点的位置矢量和三个正交方向矢量。
我很抱歉我的答复太长，但考虑到问题的复杂性，我认为最好是提供这么多的细节。
对于我的特殊问题，这些矢量以mm而不是像素为单位定义，因此间距（即每个方向上两个连续体积像素之间的距离）也用作输入。我使用了最近邻方法来插值切片的亚像素点。
reslice_volume (volume, spacing, o1, o2, n, pos)
这个算法的主要步骤如下。注意，我可以互换使用平面和切片：

1.获得所需平面和体积边界之间的交线。

def PlaneBoundsIntersectionsLines (n, pos):
  """Outputs points and vectors defining the lines that the view creates by intersecting the volume's bounds.
      Input:
        Normal vector of the given plane and the coords of a point belonging to the plane.
      Output:
        normals_line, points_line
        """
    def intersectionPlanePlane(n1,p1,n2,p2):
        # Get direction of line
        nout = np.cross(n1.reshape((1,3)),n2.reshape((1,3))).reshape(3,1)
        nout = normalizeLength(nout)
        M =  np.concatenate((n1.reshape(1,3),n2.reshape(1,3)), axis=0)
        b = np.zeros((2,1))
        # print(n1.shape, p1.shape)
        b[0,0]=np.dot(n1,p1)
        b[1,0]=np.dot(n2,p2)
        pout,resid,rank,s = np.linalg.lstsq(M,b, rcond=None)
        return pout, nout
      # ... For each face
      normalFaces = np.concatenate((np.eye(3,3),np.eye(3,3)), axis = 1)
      pointsFaces = np.array([[0,0,0],[0,0,0],[0,0,0], [379.9872, 379.9872, 169.5], [379.9872, 379.9872, 169.5], [379.9872, 379.9872, 169.5]]).transpose()
      points_line = np.zeros((3,6))
      normals_line = np.zeros((3,6))
      for face in range(6):
        n1 = normalFaces[:,face].reshape(3,)
        p1 = pointsFaces[:,face].reshape(3,)
        pout, nout = intersectionPlanePlane(n1,p1,n,pos)
        points_line[:,face] = pout.reshape((3,))
        normals_line[:,face] = nout.reshape((3,))
      return normals_line, points_line

2.获得这些线之间的交点，这些交点足够接近体积的边界，以被视为平面和体积之间的交点的角。

def FindPlaneCorners(normals_line, points_line):
  """Outputs the points defined by the intersection of the input lines that 
  are close enough to the borders of the volume to be considered corners of the view plane.
      Input:
        Points and vectors defining lines
      Output:
        p_intersection, intersecting_lines
        """
  def intersectionLineLine(Up,P0,Uq,Q0):
    # Computes the closest point between two lines
    # Must be column points
    b = np.zeros((2,1))
    b[0,0] = -np.dot((P0-Q0),Up)
    b[1,0] = -np.dot((P0-Q0),Uq)
    A = np.zeros((2,2))
    A[0,0] = np.dot(Up,Up)
    A[0,1] = np.dot(-Uq,Up)
    A[1,0] = np.dot(Up,Uq)
    A[1,1] = np.dot(-Uq,Uq)
    if ( np.abs(np.linalg.det(A)) < 10^(-10) ):
      point = np.array([np.nan, np.nan, np.nan]).reshape(3,1)
    else:
      lbd ,resid,rank,s = np.linalg.lstsq(A,b, rcond=None)
      # print('\n')
      # print(lbd)
      P1 = P0 + lbd[0]*Up;
      Q1 = Q0 + lbd[1]*Uq;
      point = (P1+Q1)/2;
    return point
  # ... ... Get closest point for every possible pair of lines and select only the ones inside the box
  npts = 0
  p_intersection = []
  intersecting_lines = []
  # ... Get all possible pairs of lines
  possible_pairs = np.array(list(itertools.combinations(np.linspace(0,5,6), 2)))
  for pair in possible_pairs:
    k = int(pair[0])
    j = int(pair[1])
    Up = normals_line[:,k]
    P0 = points_line[:,k]
    Uq = normals_line[:,j]
    Q0 = points_line[:,j]
    closest_point = intersectionLineLine(Up,P0,Uq,Q0)
    epsilon = 2.2204e-10
    # ... ... Is point inside volume? Is it close to the border?
    if closest_point[0] <= 379.9872 + epsilon and closest_point[0] >= 0 - epsilon and \
                closest_point[1] <= 379.9872 + epsilon and closest_point[1] >= 0 - epsilon and \
                closest_point[2] <= 169.5 + epsilon and closest_point[2] >= 0 - epsilon:
      # ... ...  Is it close to the border? 25 mm?
      th = 25
      if 379.9872 - closest_point[0] <= th or closest_point[0] - 0 <= th or \
          379.9872 - closest_point[1] <= th or closest_point[1] - 0 <= th or \
          169.5 - closest_point[2] <= th or closest_point[2] - 0 <= th:
        # print('It is close to teh border')
        npts += 1
        p_intersection.append(closest_point)
        intersecting_lines.append([k,j])

  p_intersection = np.array(p_intersection).transpose()

  return p_intersection, intersecting_lines

**3.**将找到的点转换到切片的参考系（sRF）（我们可以在切片平面内任意地将RF居中）。

dim = volume.shape
  # ... Get intersection lines between plane and volume bounds
  normals_line, points_line = PlaneBoundsIntersectionsLines (n, pos)

  # ... Get intersections between generated lines to get corners of view plane
  p_intersection, intersecting_lines = FindPlaneCorners (normals_line, points_line)

  # ... Calculate parameters of the 2D slice
  # ... ... Get corners of slice from volume RF (vrf) to slice RF (srf) - in this case centered in the middle of teh slice
  # ... ... ... Define T_vrf2srf
  Pose_slice_vrf = M_creater(o1,o2,n,pos)
  # ... ... ... Apply transform
  p_intersection_slicerf = np.zeros(p_intersection.shape)
  for corner in range(p_intersection.shape[1]):
    pt_arr = np.concatenate((p_intersection[:,corner],np.ones((1,))) ,axis = 0).reshape((4,1))
    p_intersection_slicerf[:,corner] = np.matmul(np.linalg.inv(Pose_slice_vrf), pt_arr)[:-1].reshape((3,))

4.****获得穿过这些点的最小x和y坐标并定义将被用作平面/切片的原点的角点。将该原点变换回体积的RF（vRF）并定义新的变换矩阵，其将RF从vRF变换到sRF，但现在以所述原点为中心。
5.根据这些切片内点坐标，我们可以确定切片的大小，然后使用它来生成目标切片的所有可能的切片内索引。

# ... ... Get slice size based on corners and spacing

  spacing_slice = [1,1,8]
  min_bounds_slice_xy = np.min(p_intersection_slicerf,axis=1)
  max_bounds_slice_xy = np.max(p_intersection_slicerf,axis=1)
  size_slice_x = int(np.ceil((max_bounds_slice_xy[0] - min_bounds_slice_xy[0] - 1e-6) / spacing_slice[0]))
  size_slice_y = int(np.ceil((max_bounds_slice_xy[1] - min_bounds_slice_xy[1] - 1e-6) / spacing_slice[1]))
  slice_size = [size_slice_x, size_slice_y, 1]
  print('slice_size')
  print(slice_size)

  # ... ... Get corner in slice coords and redefine transform mat - make corner origin of the slice
  origin_corner_slice = np.array([min_bounds_slice_xy[0],min_bounds_slice_xy[1],0])
  pt_arr = np.concatenate((origin_corner_slice,np.ones((1,))) ,axis = 0).reshape((4,1))
  origin_corner_slice_vrf = np.matmul(Pose_slice_vrf, pt_arr)[:-1].reshape((3,))
  Pose_slice_origin_corner_vrf = M_creater(o1,o2,n,origin_corner_slice_vrf)

  # ... ... Get every possible inslice coordinates
  xvalues = np.linspace(0,size_slice_x-1,size_slice_x)
  yvalues = np.linspace(0,size_slice_y-1,size_slice_y)
  zvalues = np.linspace(0,0,1)
  xx, yy = np.meshgrid(xvalues, yvalues)
  xx = xx.transpose()
  yy = yy.transpose()
  zz = np.zeros(xx.shape)
  inslice_coords = np.concatenate((xx.reshape(-1,1), yy.reshape(-1,1), zz.reshape(-1,1)), axis = 1)

6.下一步是使用新定义的变换矩阵（步骤4）来将每个可能的切片内索引Map到卷的参考帧。

# ... ... Map every point of slice into volume's RF

  inslice_coords_vrf = np.zeros(inslice_coords.shape)
  for coord_set in range(inslice_coords.shape[0]):
    pt_arr = np.concatenate((inslice_coords[coord_set,:],np.ones((1,))) ,axis = 0).reshape((4,1))
    inslice_coords_vrf[coord_set,:] = np.matmul(Pose_slice_origin_corner_vrf, pt_arr)[:-1].reshape((3,))

7.现在，我们有了切片包含的所有vRF坐标，这些坐标应立即转换为像素值，方法是将它们除以各自的间距。在这一步，我们发现，当切片经过体积的子像素位置时，我们最终得到非整数像素值。我们将像素值四舍五入到其最接近的整数- * 最近邻插值 *。

# ... ... ... Convert to pixel coord - here we used teh resampled spacing

  inslice_coords_vrf_px = inslice_coords_vrf.copy()
  inslice_coords_vrf_px[:,0] = inslice_coords_vrf[:,0] / spacing[0]
  inslice_coords_vrf_px[:,1] = inslice_coords_vrf[:,1] / spacing[1]
  inslice_coords_vrf_px[:,2] = inslice_coords_vrf[:,2] / spacing[2]

  # ... ... Interpolate pixel value at each mapped point - nearest neighbour int
  # ... ... ... Convert pixel value to its closest existing value in the volume
  inslice_coords_vrf_px = np.round(inslice_coords_vrf_px, 0).astype(int)

**8.**接下来，我们确定切片的哪些像素实际上在体积的边界内，并得到它们的值。体积外的像素被填充为0。

# ... ... Find slice voxels within volume bounds

  in_mask = np.zeros((inslice_coords_vrf_px.shape[0], 1))
  idx_in = []
  for vox in range(in_mask.shape[0]):
    if not np.any(inslice_coords_vrf_px[vox,:]<0) and \
                  inslice_coords_vrf_px[vox,0]<dim[0] and \
                  inslice_coords_vrf_px[vox,1]<dim[1] and \
                  inslice_coords_vrf_px[vox,2]<dim[2]:
      in_mask[vox] = 1
      idx_in.append(vox)
  idx_in = np.array(idx_in)
  # ... ... Get pixel value from volume based on interpolated pixel indexes
  extracted_slice = np.zeros((inslice_coords_vrf_px.shape[0], 1))
  for point in range(inslice_coords_vrf_px.shape[0]):
    if point in idx_in:
      vol_idx = inslice_coords_vrf_px[point,:]
      extracted_slice[point] = volume[vol_idx[0], vol_idx[1], vol_idx[2]]
  # ... ... Reshape to slice shape
  extracted_slice = extracted_slice.reshape((slice_size[0], slice_size[1]))

为了更加清晰，我添加了一个图。在这里，体积由黑色的边界框定义。切片与由框/体积的面定义的平面的线交叉点用橙子虚线表示。蓝色表示前面的线之间的交叉点。粉红色的点属于切片，橙色的点属于切片，并且在体积内。
在我的例子中，我处理的是MRI体积，所以作为例子，我添加了体积的结果切片。