def tensor_spherical_to_cartesian(theta, phi, psi):
"""Calculate the eigenvectors for a Tensor given the three angles.
This will return the eigenvectors unsorted, since this function knows nothing about the eigenvalues. The caller
of this function will have to sort them by eigenvalue if necessary.
Args:
theta (ndarray): matrix of list of theta's
phi (ndarray): matrix of list of phi's
psi (ndarray): matrix of list of psi's
Returns:
tuple: The three eigenvector for every voxel given. The return matrix for every eigenvector is of the given
shape + [3].
"""
v0 = spherical_to_cartesian(theta, phi)
v1 = rotate_orthogonal_vector(v0, spherical_to_cartesian(theta + np.pi / 2.0, phi), psi)
v2 = np.cross(v0, v1)
return v0, v1, v2
python类cross()的实例源码
def rotate_vector(basis, to_rotate, psi):
"""Uses Rodrigues' rotation formula to rotate the given vector v by psi around k.
If a matrix is given the operation will by applied on the last dimension.
Args:
basis: the unit vector defining the rotation axis (k)
to_rotate: the vector to rotate by the angle psi (v)
psi: the rotation angle (psi)
Returns:
vector: the rotated vector
"""
cross_product = np.cross(basis, to_rotate)
dot_product = np.sum(np.multiply(basis, to_rotate), axis=-1)[..., None]
cos_psi = np.cos(psi)[..., None]
sin_psi = np.sin(psi)[..., None]
return to_rotate * cos_psi + cross_product * sin_psi + basis * dot_product * (1 - cos_psi)
def rotate_orthogonal_vector(basis, to_rotate, psi):
"""Uses Rodrigues' rotation formula to rotate the given vector v by psi around k.
If a matrix is given the operation will by applied on the last dimension.
This function assumes that the given two vectors (or matrix of vectors) are orthogonal for every voxel.
This assumption allows for some speedup in the rotation calculation.
Args:
basis: the unit vector defining the rotation axis (k)
to_rotate: the vector to rotate by the angle psi (v)
psi: the rotation angle (psi)
Returns:
vector: the rotated vector
"""
cross_product = np.cross(basis, to_rotate)
cos_psi = np.cos(psi)[..., None]
sin_psi = np.sin(psi)[..., None]
return to_rotate * cos_psi + cross_product * sin_psi
def signed_angle(v1, v2, look):
'''
Compute the signed angle between two vectors.
Returns a number between -180 and 180. A positive number indicates a
clockwise sweep from v1 to v2. A negative number is counterclockwise.
'''
# The sign of (A x B) dot look gives the sign of the angle.
# > 0 means clockwise, < 0 is counterclockwise.
sign = np.sign(np.cross(v1, v2).dot(look))
# 0 means collinear: 0 or 180. Let's call that clockwise.
if sign == 0:
sign = 1
return sign * angle(v1, v2, look)
def rotation_from_up_and_look(up, look):
'''
Rotation matrix to rotate a mesh into a canonical reference frame. The
result is a rotation matrix that will make up along +y and look along +z
(i.e. facing towards a default opengl camera).
Note that if you're reorienting a mesh, you can use its `reorient` method
to accomplish this.
up: The foot-to-head direction.
look: The direction the eyes are facing, or the heel-to-toe direction.
'''
up, look = np.array(up, dtype=np.float64), np.array(look, dtype=np.float64)
if np.linalg.norm(up) == 0:
raise ValueError("Singular up")
if np.linalg.norm(look) == 0:
raise ValueError("Singular look")
y = up / np.linalg.norm(up)
z = look - np.dot(look, y)*y
if np.linalg.norm(z) == 0:
raise ValueError("up and look are colinear")
z = z / np.linalg.norm(z)
x = np.cross(y, z)
return np.array([x, y, z])
def from_points(cls, p1, p2, p3):
'''
If the points are oriented in a counterclockwise direction, the plane's
normal extends towards you.
'''
from blmath.numerics import as_numeric_array
p1 = as_numeric_array(p1, shape=(3,))
p2 = as_numeric_array(p2, shape=(3,))
p3 = as_numeric_array(p3, shape=(3,))
v1 = p2 - p1
v2 = p3 - p1
normal = np.cross(v1, v2)
return cls(point_on_plane=p1, unit_normal=normal)
def from_points_and_vector(cls, p1, p2, vector):
'''
Compute a plane which contains two given points and the given
vector. Its reference point will be p1.
For example, to find the vertical plane that passes through
two landmarks:
from_points_and_normal(p1, p2, vector)
Another way to think about this: identify the plane to which
your result plane should be perpendicular, and specify vector
as its normal vector.
'''
from blmath.numerics import as_numeric_array
p1 = as_numeric_array(p1, shape=(3,))
p2 = as_numeric_array(p2, shape=(3,))
v1 = p2 - p1
v2 = as_numeric_array(vector, shape=(3,))
normal = np.cross(v1, v2)
return cls(point_on_plane=p1, unit_normal=normal)
def rpy(self):
acc = self.acceleration()
yaw = self.yaw()
norm = np.linalg.norm(acc)
# print(acc)
if norm < 1e-6:
return (0.0, 0.0, yaw)
else:
thrust = acc + np.array([0, 0, 9.81])
z_body = thrust / np.linalg.norm(thrust)
x_world = np.array([math.cos(yaw), math.sin(yaw), 0])
y_body = np.cross(z_body, x_world)
x_body = np.cross(y_body, z_body)
pitch = math.asin(-x_body[2])
roll = math.atan2(y_body[2], z_body[2])
return (roll, pitch, yaw)
# "private" methods
def triangleArea(triangleSet):
"""
Calculate areas of subdivided triangles
Input: the set of subdivided triangles
Output: a list of the areas with corresponding idices with the the triangleSet
"""
triangleAreaSet = []
for i in range(len(triangleSet)):
v1 = triangleSet[i][1] - triangleSet[i][0]
v2 = triangleSet[i][2] - triangleSet[i][0]
area = np.linalg.norm(np.cross(v1, v2))/2
triangleAreaSet.append(area)
return triangleAreaSet
def crossArea(forceVecs,triangleAreaSet,triNormVecs):
"""
Preparation for Young's Modulus
Calculate the cross sectional areas perpendicular to the force vectors
Input: forceVecs = a list of force vectors
triangleAreaSet = area of triangles
triNormVecs = a list of normal vectors for each triangle (should be given by the stl file)
Output: A list of cross sectional area, approximated by the area of the triangle perpendicular to the force vector
"""
crossAreaSet = np.zeros(len(triangleAreaSet))
for i in range(len(forceVecs)):
costheta = np.dot(forceVecs[i],triNormVecs[i])/(np.linalg.norm(forceVecs[i])*np.linalg.norm(triNormVecs[i]))
crossAreaSet[i] = abs(costheta*triangleAreaSet[i])
return crossAreaSet
def computeNormals(vtx, idx):
nrml = numpy.zeros(vtx.shape, numpy.float32)
# compute normal per triangle
triN = numpy.cross(vtx[idx[:,1]] - vtx[idx[:,0]], vtx[idx[:,2]] - vtx[idx[:,0]])
# sum normals at vtx
nrml[idx[:,0]] += triN[:]
nrml[idx[:,1]] += triN[:]
nrml[idx[:,2]] += triN[:]
# compute norms
nrmlNorm = numpy.sqrt(nrml[:,0]*nrml[:,0]+nrml[:,1]*nrml[:,1]+nrml[:,2]*nrml[:,2])
return nrml/nrmlNorm.reshape(-1,1)
def convex_hull(points, vind, nind, tind, obj):
"super ineffective"
cnt = len(points)
for a in range(cnt):
for b in range(a+1,cnt):
for c in range(b+1,cnt):
vec1 = points[a] - points[b]
vec2 = points[a] - points[c]
n = np.cross(vec1, vec2)
n /= np.linalg.norm(n)
C = np.dot(n, points[a])
inner = np.inner(n, points)
pos = (inner <= C+0.0001).all()
neg = (inner >= C-0.0001).all()
if not pos and not neg: continue
obj.out.write("f %i//%i %i//%i %i//%i\n" % (
(vind[a], nind[a], vind[b], nind[b], vind[c], nind[c])
if (inner - C).sum() < 0 else
(vind[a], nind[a], vind[c], nind[c], vind[b], nind[b]) ) )
#obj.out.write("f %i/%i/%i %i/%i/%i %i/%i/%i\n" % (
# (vind[a], tind[a], nind[a], vind[b], tind[b], nind[b], vind[c], tind[c], nind[c])
# if (inner - C).sum() < 0 else
# (vind[a], tind[a], nind[a], vind[c], tind[c], nind[c], vind[b], tind[b], nind[b]) ) )
def test_broadcasting_shapes(self):
u = np.ones((2, 1, 3))
v = np.ones((5, 3))
assert_equal(np.cross(u, v).shape, (2, 5, 3))
u = np.ones((10, 3, 5))
v = np.ones((2, 5))
assert_equal(np.cross(u, v, axisa=1, axisb=0).shape, (10, 5, 3))
assert_raises(ValueError, np.cross, u, v, axisa=1, axisb=2)
assert_raises(ValueError, np.cross, u, v, axisa=3, axisb=0)
u = np.ones((10, 3, 5, 7))
v = np.ones((5, 7, 2))
assert_equal(np.cross(u, v, axisa=1, axisc=2).shape, (10, 5, 3, 7))
assert_raises(ValueError, np.cross, u, v, axisa=-5, axisb=2)
assert_raises(ValueError, np.cross, u, v, axisa=1, axisb=-4)
# gh-5885
u = np.ones((3, 4, 2))
for axisc in range(-2, 2):
assert_equal(np.cross(u, u, axisc=axisc).shape, (3, 4))
def _setup_normalized_vectors(self, normal_vector, north_vector):
normal_vector, north_vector = _validate_unit_vectors(normal_vector,
north_vector)
mylog.debug('Setting normalized vectors' + str(normal_vector)
+ str(north_vector))
# Now we set up our various vectors
normal_vector /= np.sqrt(np.dot(normal_vector, normal_vector))
if north_vector is None:
vecs = np.identity(3)
t = np.cross(normal_vector, vecs).sum(axis=1)
ax = t.argmax()
east_vector = np.cross(vecs[ax, :], normal_vector).ravel()
# self.north_vector must remain None otherwise rotations about a fixed axis will break.
# The north_vector calculated here will still be included in self.unit_vectors.
north_vector = np.cross(normal_vector, east_vector).ravel()
else:
if self.steady_north or (np.dot(north_vector, normal_vector) != 0.0):
north_vector = north_vector - np.dot(north_vector,normal_vector)*normal_vector
east_vector = np.cross(north_vector, normal_vector).ravel()
north_vector /= np.sqrt(np.dot(north_vector, north_vector))
east_vector /= np.sqrt(np.dot(east_vector, east_vector))
self.normal_vector = normal_vector
self.north_vector = north_vector
self.unit_vectors = YTArray([east_vector, north_vector, normal_vector], "")
self.inv_mat = np.linalg.pinv(self.unit_vectors)
def base_vectors(self):
""" Returns 3 orthognal base vectors, the first one colinear to
the axis of the loop.
"""
# normalize n
n = self.direction / (self.direction**2).sum(axis=-1)
# choose two vectors perpendicular to n
# choice is arbitrary since the coil is symetric about n
if np.abs(n[0])==1 :
l = np.r_[n[2], 0, -n[0]]
else:
l = np.r_[0, n[2], -n[1]]
l /= (l**2).sum(axis=-1)
m = np.cross(n, l)
return n, l, m
def base_vectors(n):
""" Returns 3 orthognal base vectors, the first one colinear to n.
"""
# normalize n
n = n / np.sqrt(np.square(n).sum(axis=-1))
# choose two vectors perpendicular to n
# choice is arbitrary since the coil is symetric about n
if abs(n[0]) == 1 :
l = np.r_[n[2], 0, -n[0]]
else:
l = np.r_[0, n[2], -n[1]]
l = l / np.sqrt(np.square(l).sum(axis=-1))
m = np.cross(n, l)
return n, l, m
def normal(self, t, above=True):
""" Evaluate the normal of the curve at the given parametric value(s).
This function returns an *n* × 3 array, where *n* is the number of
evaluation points.
The normal is computed as the cross product between the binormal and
the tangent of the curve.
:param t: Parametric coordinates in which to evaluate
:type t: float or [float]
:param bool above: Evaluation in the limit from above
:return: Derivative array
:rtype: numpy.array
"""
# error test input
if self.dimension != 3:
raise RuntimeError('Normals require dimension = 3')
# compute derivative
T = self.tangent(t, above=above)
B = self.binormal(t, above=above)
return np.cross(B,T)
def test_curvature(self):
# linear curves have zero curvature
crv = Curve()
self.assertAlmostEqual(crv.curvature(.3), 0.0)
# test multiple evaluation points
t = np.linspace(0,1, 10)
k = crv.curvature(t)
self.assertTrue(np.allclose(k, 0.0))
# test circle
crv = CurveFactory.circle(r=3) + [1,1]
t = np.linspace(0,2*pi, 10)
k = crv.curvature(t)
self.assertTrue(np.allclose(k, 1.0/3.0)) # circles: k = 1/r
# test 3D (np.cross has different behaviour in 2D/3D)
crv.set_dimension(3)
k = crv.curvature(t)
self.assertTrue(np.allclose(k, 1.0/3.0)) # circles: k = 1/r
def thru_plane_position(dcm):
"""Gets spatial coordinate of image origin whose axis
is perpendicular to image plane.
"""
orientation = tuple((float(o) for o in dcm.ImageOrientationPatient))
position = tuple((float(p) for p in dcm.ImagePositionPatient))
rowvec, colvec = orientation[:3], orientation[3:]
normal_vector = np.cross(rowvec, colvec)
slice_pos = np.dot(position, normal_vector)
return slice_pos
def read_pose(gt):
cam_dir, cam_up = gt.cam_dir, gt.cam_up
z = cam_dir / np.linalg.norm(cam_dir)
x = np.cross(cam_up, z)
y = np.cross(z, x)
R = np.vstack([x, y, z]).T
t = gt.cam_pos / 1000.0
return RigidTransform.from_Rt(R, t)
