如何使用 Scipy 从 3D 体积中提取任意 2D 切片？

Question

我正在使用 Scipy 从 3D 数据（矢量 200x200x200）渲染平面。 我可以通过 2 个向量或向量和角度来指定所需的平面。 我想从这个 3D 体积中提取这样一个任意切片。 我找到了如何在 Matlab 中做到这一点： http://www.mathworks.com/help/techdoc/ref/slice.html 我如何在 Scipy 中做到这一点？

Answer 1

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

Answer 2

def extract_slice（数据，三元组）：
”“”
算法：
1.找到平面与数据框边缘的交点
2.对于这些点，找到轴导向的b-box
3. 找到从 R2 到 R3 的“后”反式 (A,T)，如 X' = AX + T
使用 (0,0), (0,h), (w,0) 很容易计算
4. 对 2D (w,h) 图像中的每个值使用反式（使用三线性插值）
“”“

我相信作为该项目的一部分，我将在几个月内正确发布代码。

Answer 3

我正在解决与OP类似的任务，所以我想出了基于numpy（不是scipy）的代码，在给定平面上任何点的位置向量和三个正交方向向量的情况下，从体积中提取任何给定的切片。

对于我的回答太长，我深表歉意，但考虑到手头问题的复杂性，我认为最好提供这么多细节。

对于我的特定问题，这些向量以毫米而不是像素为单位定义，因此间距（即每个方向上两个连续体积像素之间的距离）也用作输入。我使用了最近邻方法来插值切片的子像素点。

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 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 体积，因此作为示例，我添加了从体积中生成的切片。

Answer 4

scipy.ndimage.affine_transform 带有旋转矩阵、偏移量和output_shape。

scipy.spatial.transform.Rotation.as_matrix()

可用于创建旋转矩阵。

output_shape

可以是 2D。