如何匹配不同网格的值？

Question

我试图从圆柱体的半径、大小和方向信息开始创建圆柱体的灰度图像堆栈

这个圆柱体应该包含在一个立方体3D网格中，每个网格点代表一个像素

代码应该是这样的

  #Define meshgrid
  x_ = np.linspace(0,Lx,int(Lx/vox))
  y_ = np.linspace(0,Ly,int(Ly/vox))
  z_ = np.linspace(0,Lz,int(Lz/vox))

  X,Y,Z = np.meshgrid(x_,y_,z_,indexing='ij')

其中vox是我输入的体素大小。这将定义3D网格，我们将从中导出灰度堆栈的矩阵将定义为

  M = np.zeros((x_.size,y_.size,z_.size),dtype=int)

现在我们看到一个圆柱体，它的中心轴从点p1到点p2，半径为r。因此，基于以下链接Numpy mask from cylinder coordinates，我创建了一个额外的网格，其中包含属于该圆柱体的所有点

        #Vector
        v = p2 - p1

        #Normalize vector
        lenght = scipy.linalg.norm(v)
        v = v / lenght

        # make some vector not in the same direction as v
        not_v = np.array([1.0, 0, 0])
        if (v == not_v).all():
              not_v = np.array([0, 1.0, 0])
        # make vector perpendicular to v
        n1 = np.cross(v, not_v)
        # normalize n1
        n1 = n1 / scipy.linalg.norm(n1)
        # make unit vector perpendicular to v and n1
        n2 = np.cross(v, n1)

        #Define gridpoints for the cilinder
        l_ =      np.linspace(0,lenght,100)
        r_ =      np.linspace(0,r,10)
        theeta_ = np.linspace(0,2*np.pi,10)

        #define meshgrid for cilinder
        L,R,Theeta = np.meshgrid(l_,r_,theeta_,indexing='ij')

然后我会从这些柱坐标中得到xyz坐标

#Transform to x, y, z coordinates
Xc, Yc, Zc = [p1[i] + v[i] * L + R * np.sin(Theeta) * n1[i] + r * np.cos(Theeta) * n2[i] for i in [0, 1, 2]]

所以现在的问题如下。我已经在笛卡尔坐标系中定义了构成给定圆柱体的网格点。但是，我现在遇到的问题是尝试将圆柱体的网格网格与定义的网格网格相匹配，以获得图像的灰度堆栈

所以我想做的是取这个Xc，Yc，Zc，坐标，看看它们对应于哪个X，Y，Z坐标，然后点亮矩阵M中相应的像素。但是我看不到一个明显的方法

问候,

Answer 1

如果以错误的方式处理问题，则应直接查找坐标网格中的哪些点位于圆柱体内。这就是你可以做到的：

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# Problem parameters
Lx, Ly, Lz = 50, 20, 30
vox = 2.0
p1 = np.array([12., 8., 4.])
p2 = np.array([37., 14., 22.])
r = 5.
#Define meshgrid
x_ = np.linspace(0, Lx, int(Lx / vox))
y_ = np.linspace(0, Ly, int(Ly / vox))
z_ = np.linspace(0, Lz, int(Lz / vox))
# Grid coordinates
X, Y, Z = np.meshgrid(x_, y_, z_, indexing='ij')
# Stack into an array
coords = np.stack([X, Y, Z], axis=-1)
# Compute distance from each point to cylinder axis
v = p2 - p1
t = np.dot(coords - p1, v) / np.dot(v, v)
p = p1 + np.expand_dims(t, axis=-1) * v
dist = np.linalg.norm(coords - p, axis=-1)
# Select points within cylinder distance and bounds
mask = (dist <= r) & (0 <= t) & (t <= 1)
print(mask.shape)
# (10, 4, 6)
# Select coordinates in cylinder
cyl_coords = coords[mask]
# Plot
ax = plt.figure().add_subplot(111, projection='3d')
ax.scatter3D(cyl_coords[:, 0], cyl_coords[:, 1], cyl_coords[:, 2], s=3)
# Set plot limits for uniform aspect ratio
ax.set_xlim(0, 50)
ax.set_ylim(-15, 35)
ax.set_zlim(-10, 40)
ax.figure.tight_layout()
ax.figure.show()

输出：