如何匹配不同网格的值?
How to match values of different meshgrids?
我正在尝试根据有关圆柱体的半径、大小和方向的信息创建圆柱体的灰度图像堆栈。
此圆柱体应包含在立方体 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。因此,基于以下 link ,我创建了一个额外的网格,其中包含属于此类圆柱体的所有点。
#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中对应的像素点。但是我看不到明显的方法。
此致,
您以错误的方式解决问题,您应该做的是直接找出坐标网格中的哪些点落在圆柱体内。您可以这样做:
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()
输出:
我正在尝试根据有关圆柱体的半径、大小和方向的信息创建圆柱体的灰度图像堆栈。
此圆柱体应包含在立方体 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。因此,基于以下 link
#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中对应的像素点。但是我看不到明显的方法。
此致,
您以错误的方式解决问题,您应该做的是直接找出坐标网格中的哪些点落在圆柱体内。您可以这样做:
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()
输出: