numpy 如何在python中使用原点和法向量在3D中采样点

5us2dqdw  于 2022-11-23  发布在  Python
关注(0)|答案(2)|浏览(202)

我有两个点p1(x1,y1,z1)和p2(x2,y2,z2)。我想在以p1为中心的半径为r的圆上采样点,该圆所在的平面与矢量p2-p1垂直(因此p2-p1将是该平面的法向量)。我有在XOY平面使用极坐标系统采样的代码,但苦于如何推广到不同于(0,0,1)的法线

rho = np.linspace(0, 2*np.pi, 50)
r = 1
x = np.cos(rho) * r
y = np.sin(rho) * r
z = np.zeros(rho.shape)

Sampled points

rqdpfwrv

rqdpfwrv1#

首先,您需要在圆的平面中定义两个基本向量。
第一个是正交于法线n = p2-p1的任意矢量
选择具有最大幅值的法线分量和具有第二幅值的分量。
交换它们的值,取最大值的反,并使第三个分量为零(注意,结果与法线的点积为零,因此它们是正交的)
例如,如果n.y最大,n.z次之,则使

v = (0, n.z, -n.y)

然后使用矢量积计算第二个基矢量

u = n x v

规范化向量vu。使用向量形式上的中心点p1圈出点:

f(rho) = p1 + r * v * cos(rho) + r * u * sin(rho)

或在组件中:

f.x = p1.x + r * v.x * cos(rho) + r * u.x * sin(rho)
and so on
8cdiaqws

8cdiaqws2#

假设我们有一个向量n,我们想找到一个以圆心p1为圆心、半径为r、与n正交的点组成的圆。

p1 = np.array([-21.03181359,   4.54876345,  19.26943601])
n = np.array([-0.06592715,  0.00713031, -0.26809672])
n = n / np.linalg.norm(n) # normalise n
r = 0.5

x = np.array([1,0,0]).astype(np.float64) # take a random vector of magnitude 1
x -= x.dot(n) * n / np.linalg.norm(n)**2  # make it orthogonal to n
x /= np.linalg.norm(x)  # normalize

# find first point on circle (x1). 
# currently it has magnitude of 1, so we multiply it by the r
x1 = p1 + (x*r)

# vector from lumen centre to first circle point
p1x1 = x1 - p1

def rotation_matrix(axis, theta):
    """
    Return the rotation matrix associated with counterclockwise rotation about
    the given axis by theta radians.
    """
    axis = np.asarray(axis)
    axis = axis / math.sqrt(np.dot(axis, axis))
    a = math.cos(theta / 2.0)
    b, c, d = -axis * math.sin(theta / 2.0)
    aa, bb, cc, dd = a * a, b * b, c * c, d * d
    bc, ad, ac, ab, bd, cd = b * c, a * d, a * c, a * b, b * d, c * d
    return np.array([[aa + bb - cc - dd, 2 * (bc + ad), 2 * (bd - ac)],
                     [2 * (bc - ad), aa + cc - bb - dd, 2 * (cd + ab)],
                     [2 * (bd + ac), 2 * (cd - ab), aa + dd - bb - cc]])

# rotate the vector p1x1 around the axis n with angle theta
circle = []
for theta in range(0,360,6):
    circle_i = np.dot(rotation_matrix(n, np.deg2rad(theta)), p1x1)
    circle.append(circle_i+p1)

ax = axes3d.Axes3D(plt.figure(figsize=(10,10)))
ax.scatter3D(*np.array(circle).T, s=10, c='red')
ax.scatter3D(*p1.T, s=10, c='black')
ax.set_xlabel('X', size=40)
ax.set_ylabel('Y', size=40)
ax.set_zlabel('Z', size=40)

ax.set_xlim(-19,-22)
ax.set_ylim(2,5)
ax.set_zlim(18,21)

相关问题