CODEWITHSUNDEEP
pythonnumpygeometrylinear-algebra
Pxxxx96
Pxxxx96

Reputation: 1

how to sample points in 3D in python with origin and normal vector

I have two points p1(x1, y1, z1) and p2(x2, y2, z2) in 3D. And I want to sample points in a circle with radius r that is centered at p1, and the plane which is perpendicular to the vector p2-p1 (so p2-p1 would be the normal vector of that plane). I have the code for sampling in XOY plane using polar system, but suffering on how to generalize to a different normal than (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

Upvotes: 0

Views: 761

Answers (2)

noble_kazoo
noble_kazoo

Reputation: 51

Lets say we have a vector n and we want to find a circle of points around a center p1 with radius r which are orthogonal to n. Here is a working example with code

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)

Upvotes: 0

MBo
MBo

Reputation: 80187

At first you need to define two base vectors in the circle's plane.

The first one is arbitrary vector orthogonal to normal n = p2-p1

Choose component of normal with the largest magnitude and component with the second magnitude.

Exchange their values, negate the largest, and make the third component zero (note that dot product of result with normal is zero, so they are othogonal)

For example, if n.y is the largest and n.z is the second, make

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

Then calculate the second base vector using vector product

u = n x v

Normalize vectors v and u. Circle points using center point p1 on vector form:

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

or in components:

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

Upvotes: 2

Related Questions