来自圆柱坐标的Numpy遮罩

nr9pn0ug  于 2022-11-29  发布在  其他
关注(0)|答案(1)|浏览(158)

我生成了一个圆柱体的坐标,它的两个面连接着两个已经给定的任意点。

有没有可能用标准的Python库从坐标中构建一个填充圆柱体的3D numpy遮罩?创建一个2D遮罩看起来很简单,但是我在3D中遇到了一些困难。
下面是用于生成圆柱体的代码,取自herehere

import scipy
import scipy.linalg
import numpy as np
import nibabel as nib
import matplotlib
matplotlib.use('TkAgg')
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# defining mask
shape = (100, 100, 100)
image = np.zeros(shape=shape)

# set radius and centres values
r = 3
start = [30, 45, 60]
end = [40, 58, 70]
p1 = np.array(start)
p2 = np.array(end)

# # calculate p2-p1 distance
# dx = p2[0] - p1[0]
# dy = p2[1] - p1[1]
# dz = p2[2] - p1[2]
# dist = math.sqrt(dx**2 + dy**2 + dz**2)

# vector in direction of axis
v = p2 - p1
# find magnitude of vector
mag = scipy.linalg.norm(v)
# unit vector in direction of axis
v = v / mag
# make some vector not in the same direction as v
not_v = np.array([1, 0, 0])
if (v == not_v).all():
    not_v = np.array([0, 1, 0])
# make vector perpendicular to v
n1 = np.cross(v, not_v)
# normalize n1
n1 /= scipy.linalg.norm(n1)
# make unit vector perpendicular to v and n1
n2 = np.cross(v, n1)
#surface ranges over t from 0 to length of axis and 0 to 2*pi
t = np.linspace(0, mag, 100)
theta = np.linspace(0, 2 * np.pi, 100)
rsample = np.linspace(0, r, 2)

#use meshgrid to make 2d arrays
t, theta2 = np.meshgrid(t, theta)
rsample, theta = np.meshgrid(rsample, theta)

# generate coordinates for surface
# "Tube"
X, Y, Z = [p1[i] + v[i] * t + r * np.sin(theta2) * n1[i] + r * np.cos(theta2) * n2[i] for i in [0, 1, 2]]
# "Bottom"
X2, Y2, Z2 = [p1[i] + rsample[i] * np.sin(theta) * n1[i] + rsample[i] * np.cos(theta) * n2[i] for i in [0, 1, 2]]
# "Top"
X3, Y3, Z3 = [p1[i] + v[i] * mag + rsample[i] * np.sin(theta) * n1[i] + rsample[i] * np.cos(theta) * n2[i] for i in [0, 1, 2]]

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(X, Y, Z)
ax.plot_surface(X2, Y2, Z2)
ax.plot_surface(X3, Y3, Z3)

plt.show()

我需要的3D numpy掩码选择一个3D图像的圆柱体内的所有值。掩码和图像的形状是相同的。

bq3bfh9z

bq3bfh9z1#

最后,我遍历了管和面的坐标,得到了以下链接的坐标:3D points from Numpy meshgrid coordinates

tube = np.stack((X.ravel(), Y.ravel(), Z.ravel()), axis=1)
face1 = np.stack((X2.ravel(), Y2.ravel(), Z2.ravel()), axis=1)
face2 = np.stack((X3.ravel(), Y3.ravel(), Z3.ravel()), axis=1)

# filling numpy mask

for i in range(len(tube)):
    image[int(tube[i][0]), int(tube[i][1]), int(tube[i][2])] = 255

for j in range(len(face1)):
    image[int(face1[j][0]), int(face1[j][1]), int(face1[j][2])] = 255

for k in range(len(face2)):
    image[int(face2[k][0]), int(face2[k][1]), int(face2[k][2])] = 255

mask_new = nib.Nifti1Image(image.astype(np.float32), ctsurg_file.affine)

nib.save(mask_new, os.path.join(currdir, 'mask_cyl.nii.gz'))

相关问题