1

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

该圆柱体应包含在立方 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中的相应像素。但我看不到明显的方法做它。

问候，

4

1 回答 1

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()

输出：

输出图像

于 2020-07-13T14:18:07.760 回答