Python numpy 模块,cross() 实例源码

我们从Python开源项目中,提取了以下50个代码示例,用于说明如何使用numpy.cross()

项目:Neural-Networks-for-Inverse-Kinematics    作者:paramrajpura    | 项目源码 | 文件源码
def vector_product(v0, v1, axis=0):
    """Return vector perpendicular to vectors.

    >>> v = vector_product([2, 0, 0], [0, 3, 0])
    >>> numpy.allclose(v, [0, 0, 6])
    True
    >>> v0 = [[2, 0, 0, 2], [0, 2, 0, 2], [0, 0, 2, 2]]
    >>> v1 = [[3], [0], [0]]
    >>> v = vector_product(v0, v1)
    >>> numpy.allclose(v, [[0, 0, 0, 0], [0, 0, 6, 6], [0, -6, 0, -6]])
    True
    >>> v0 = [[2, 0, 0], [2, 0, 0], [0, 2, 0], [2, 0, 0]]
    >>> v1 = [[0, 3, 0], [0, 0, 3], [0, 0, 3], [3, 3, 3]]
    >>> v = vector_product(v0, v1, axis=1)
    >>> numpy.allclose(v, [[0, 0, 6], [0, -6, 0], [6, 0, 0], [0, -6, 6]])
    True

    """
    return numpy.cross(v0, v1, axis=axis)
项目:pybot    作者:spillai    | 项目源码 | 文件源码
def skew(v, return_dv=False):
  """ 
  Returns the skew-symmetric matrix of a vector
  Ref: https://github.com/dreamdragon/Solve3Plus1/blob/master/skew3.m

  Also known as the cross-product matrix [v]_x such that 
  the cross product of (v x w) is equivalent to the 
  matrix multiplication of the cross product matrix of 
  v ([v]_x) and w

  In other words: v x w = [v]_x * w
  """
  sk = np.float32([[0, -v[2], v[1]],
                   [v[2], 0, -v[0]],
                   [-v[1], v[0], 0]])

  if return_dv: 
      dV = np.float32([[0, 0, 0], 
                       [0, 0, -1], 
                       [0, 1, 0], 
                       [0, 0, 1], 
                       [0, 0, 0], 
                       [-1, 0, 0], 
                       [0, -1, 0], 
                       [1, 0, 0], 
                       [0, 0, 0]])
      return sk, dV
  else: 
      return sk
项目:pybot    作者:spillai    | 项目源码 | 文件源码
def points_and_normals(self): 
        """
        Returns the point/normals parametrization for planes, 
        including clipped zmin and zmax frustums

        Note: points need to be in CCW
        """

        nv1, fv1 = self._front_back_vertices
        nv2 = np.roll(nv1, -1, axis=0)
        fv2 = np.roll(fv1, -1, axis=0)

        vx = np.vstack([fv1-nv1, nv2[0]-nv1[0], fv1[2]-fv1[1]])
        vy = np.vstack([fv2-fv1, nv2[1]-nv2[0], fv1[1]-fv1[0]])
        pts = np.vstack([fv1, nv1[0], fv1[1]])

        # vx += 1e-12
        # vy += 1e-12

        vx /= np.linalg.norm(vx, axis=1).reshape(-1,1)
        vy /= np.linalg.norm(vy, axis=1).reshape(-1,1)

        normals = np.cross(vx, vy)
        normals /= np.linalg.norm(normals, axis=1).reshape(-1,1)
        return pts, normals
项目:NeoAnalysis    作者:neoanalysis    | 项目源码 | 文件源码
def faceNormals(self, indexed=None):
        """
        Return an array (Nf, 3) of normal vectors for each face.
        If indexed='faces', then instead return an indexed array
        (Nf, 3, 3)  (this is just the same array with each vector
        copied three times).
        """
        if self._faceNormals is None:
            v = self.vertexes(indexed='faces')
            self._faceNormals = np.cross(v[:,1]-v[:,0], v[:,2]-v[:,0])

        if indexed is None:
            return self._faceNormals
        elif indexed == 'faces':
            if self._faceNormalsIndexedByFaces is None:
                norms = np.empty((self._faceNormals.shape[0], 3, 3))
                norms[:] = self._faceNormals[:,np.newaxis,:]
                self._faceNormalsIndexedByFaces = norms
            return self._faceNormalsIndexedByFaces
        else:
            raise Exception("Invalid indexing mode. Accepts: None, 'faces'")
项目:NeoAnalysis    作者:neoanalysis    | 项目源码 | 文件源码
def faceNormals(self, indexed=None):
        """
        Return an array (Nf, 3) of normal vectors for each face.
        If indexed='faces', then instead return an indexed array
        (Nf, 3, 3)  (this is just the same array with each vector
        copied three times).
        """
        if self._faceNormals is None:
            v = self.vertexes(indexed='faces')
            self._faceNormals = np.cross(v[:,1]-v[:,0], v[:,2]-v[:,0])

        if indexed is None:
            return self._faceNormals
        elif indexed == 'faces':
            if self._faceNormalsIndexedByFaces is None:
                norms = np.empty((self._faceNormals.shape[0], 3, 3))
                norms[:] = self._faceNormals[:,np.newaxis,:]
                self._faceNormalsIndexedByFaces = norms
            return self._faceNormalsIndexedByFaces
        else:
            raise Exception("Invalid indexing mode. Accepts: None, 'faces'")
项目:SLP-Annotator    作者:PhonologicalCorpusTools    | 项目源码 | 文件源码
def vector_product(v0, v1, axis=0):
    """Return vector perpendicular to vectors.

    >>> v = vector_product([2, 0, 0], [0, 3, 0])
    >>> numpy.allclose(v, [0, 0, 6])
    True
    >>> v0 = [[2, 0, 0, 2], [0, 2, 0, 2], [0, 0, 2, 2]]
    >>> v1 = [[3], [0], [0]]
    >>> v = vector_product(v0, v1)
    >>> numpy.allclose(v, [[0, 0, 0, 0], [0, 0, 6, 6], [0, -6, 0, -6]])
    True
    >>> v0 = [[2, 0, 0], [2, 0, 0], [0, 2, 0], [2, 0, 0]]
    >>> v1 = [[0, 3, 0], [0, 0, 3], [0, 0, 3], [3, 3, 3]]
    >>> v = vector_product(v0, v1, axis=1)
    >>> numpy.allclose(v, [[0, 0, 6], [0, -6, 0], [6, 0, 0], [0, -6, 6]])
    True

    """
    return numpy.cross(v0, v1, axis=axis)
项目:code    作者:ActiveState    | 项目源码 | 文件源码
def poly_area(poly):
    if len(poly) < 3: # not a plane - no area
        return 0
    total = [0, 0, 0]
    N = len(poly)
    for i in range(N):
        vi1 = poly[i]
        vi2 = poly[(i+1) % N]
        prod = np.cross(vi1, vi2)
        total[0] += prod[0]
        total[1] += prod[1]
        total[2] += prod[2]
    result = np.dot(total, unit_normal(poly[0], poly[1], poly[2]))
    return abs(result/2)

#unit normal vector of plane defined by points a, b, and c
项目:pypilot    作者:pypilot    | 项目源码 | 文件源码
def __init__(self, plane_fit, gridsize):
        plane = numpy.array(plane_fit)

        origin = -plane / numpy.dot(plane, plane)
        n = numpy.array([plane[1], plane[2], plane[0]])

        u = numpy.cross(plane, n)
        v = numpy.cross(plane, u)

        u /= numpy.linalg.norm(u)
        v /= numpy.linalg.norm(v)

        def project_point(point):
            return origin + point[0]*u + point[1]*v

        vertexes = []

        for x in range(-gridsize+1, gridsize):
            for y in range(-gridsize+1, gridsize):
                vertexes += [project_point((x-1, y-1)),
                             project_point((x, y-1)),
                             project_point((x, y)),
                             project_point((x-1, y))]

        super(self, Plane).__init__(vertexes)
项目:em_examples    作者:geoscixyz    | 项目源码 | 文件源码
def analytic_infinite_wire(obsloc,wireloc,orientation,I=1.):
    """
    Compute the response of an infinite wire with orientation 'orientation'
    and current I at the obsvervation locations obsloc

    Output:
    B: magnetic field [Bx,By,Bz]
    """

    n,d = obsloc.shape
    t,d = wireloc.shape
    d = np.sqrt(np.dot(obsloc**2.,np.ones([d,t]))+np.dot(np.ones([n,d]),(wireloc.T)**2.)
    - 2.*np.dot(obsloc,wireloc.T))
    distr = np.amin(d, axis=1, keepdims = True)
    idxmind = d.argmin(axis=1)
    r = obsloc - wireloc[idxmind]

    orient = np.c_[[orientation for i in range(obsloc.shape[0])]]
    B = (mu_0*I)/(2*np.pi*(distr**2.))*np.cross(orientation,r)

    return B
项目:gym-extensions    作者:Breakend    | 项目源码 | 文件源码
def winding_angle(self, path, point):
        wa = 0
        for i in range(len(path)-1):
            p = np.array([path[i].x, path[i].y])
            pn = np.array([path[i+1].x, path[i+1].y])

            vp = p - point
            vpn = pn - point

            vp_norm = sqrt(vp[0]**2 + vp[1]**2)
            vpn_norm = sqrt(vpn[0]**2 + vpn[1]**2)

            assert (vp_norm > 0)
            assert (vpn_norm > 0)

            z = np.cross(vp, vpn)/(vp_norm * vpn_norm)
            z = min(max(z, -1.0), 1.0)
            wa += asin(z)

        return wa
项目:Tweaker-3    作者:ChristophSchranz    | 项目源码 | 文件源码
def rotate_ascii_stl(self, rotation_matrix, content, filename):
        """Rotate the mesh array and save as ASCII STL."""
        mesh = np.array(content, dtype=np.float64)

        # prefix area vector, if not already done (e.g. in STL format)
        if len(mesh[0]) == 3:
            row_number = int(len(content)/3)
            mesh = mesh.reshape(row_number, 3, 3)

        # upgrade numpy with: "pip install numpy --upgrade"
        rotated_content = np.matmul(mesh, rotation_matrix)

        v0 = rotated_content[:, 0, :]
        v1 = rotated_content[:, 1, :]
        v2 = rotated_content[:, 2, :]
        normals = np.cross(np.subtract(v1, v0), np.subtract(v2, v0)) \
            .reshape(int(len(rotated_content)), 1, 3)
        rotated_content = np.hstack((normals, rotated_content))

        tweaked = list("solid %s" % filename)
        tweaked += list(map(self.write_facett, list(rotated_content)))
        tweaked.append("\nendsolid %s\n" % filename)
        tweaked = "".join(tweaked)

        return tweaked
项目:tensorflow_ocr    作者:BowieHsu    | 项目源码 | 文件源码
def line_cross_point(line1, line2):
    # line1 0= ax+by+c, compute the cross point of line1 and line2
    if line1[0] != 0 and line1[0] == line2[0]:
        print('Cross point does not exist')
        return None
    if line1[0] == 0 and line2[0] == 0:
        print('Cross point does not exist')
        return None
    if line1[1] == 0:
        x = -line1[2]
        y = line2[0] * x + line2[2]
    elif line2[1] == 0:
        x = -line2[2]
        y = line1[0] * x + line1[2]
    else:
        k1, _, b1 = line1
        k2, _, b2 = line2
        x = -(b1-b2)/(k1-k2)
        y = k1*x + b1
    return np.array([x, y], dtype=np.float32)
项目:radar    作者:amoose136    | 项目源码 | 文件源码
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))
项目:Sverchok    作者:Sverchok    | 项目源码 | 文件源码
def vector_product(v0, v1, axis=0):
    """Return vector perpendicular to vectors.

    >>> v = vector_product([2, 0, 0], [0, 3, 0])
    >>> numpy.allclose(v, [0, 0, 6])
    True
    >>> v0 = [[2, 0, 0, 2], [0, 2, 0, 2], [0, 0, 2, 2]]
    >>> v1 = [[3], [0], [0]]
    >>> v = vector_product(v0, v1)
    >>> numpy.allclose(v, [[0, 0, 0, 0], [0, 0, 6, 6], [0, -6, 0, -6]])
    True
    >>> v0 = [[2, 0, 0], [2, 0, 0], [0, 2, 0], [2, 0, 0]]
    >>> v1 = [[0, 3, 0], [0, 0, 3], [0, 0, 3], [3, 3, 3]]
    >>> v = vector_product(v0, v1, axis=1)
    >>> numpy.allclose(v, [[0, 0, 6], [0, -6, 0], [6, 0, 0], [0, -6, 6]])
    True

    """
    return numpy.cross(v0, v1, axis=axis)
项目:bpy_lambda    作者:bcongdon    | 项目源码 | 文件源码
def intersect(self, segment):
        """ point_sur_segment return
            p: point d'intersection
            u: param t de l'intersection sur le segment courant
            v: param t de l'intersection sur le segment segment
        """
        v2d = self.vect_2d
        c2 = np.cross(segment.vect_2d, (0, 0, 1))
        d = np.dot(v2d, c2)
        if d == 0:
            # segments paralleles
            segment._point_sur_segment(self.c0)
            segment._point_sur_segment(self.c1)
            self._point_sur_segment(segment.c0)
            self._point_sur_segment(segment.c1)
            return False, 0, 0, 0
        c1 = np.cross(v2d, (0, 0, 1))
        v3 = self.c0.vect(segment.c0)
        v3[2] = 0.0
        u = np.dot(c2, v3) / d
        v = np.dot(c1, v3) / d
        co = self.lerp(u)
        return True, co, u, v
项目:uncover-ml    作者:GeoscienceAustralia    | 项目源码 | 文件源码
def __init__(self, ROI):
        bounds = [(ROI[0], ROI[2], 0),
                  (ROI[1], ROI[2], 0),
                  (ROI[0], ROI[3], 0),
                  (ROI[1], ROI[3], 0)]
        self.wgs84 = nv.FrameE(name='WGS84')
        lon, lat, hei = np.array(bounds).T
        geo_points = self.wgs84.GeoPoint(longitude=lon, latitude=lat, z=-hei,
                                         degrees=True)
        P = geo_points.to_ecef_vector().pvector.T
        dx = normed(P[1] - P[0])
        dy = P[2] - P[0]
        dy -= dx * dy.dot(dx)
        dy = normed(dy)
        dz = np.cross(dx, dy)
        self.rotation = np.array([dx, dy, dz]).T
        self.mu = np.mean(P.dot(self.rotation), axis=0)[np.newaxis, :]
项目:meshpy    作者:BerkeleyAutomation    | 项目源码 | 文件源码
def _signed_volume_of_tri(self, tri):
        """Return the signed volume of the given triangle.

        Parameters
        ----------
        tri : :obj:`numpy.ndarray` of int
            The triangle for which we wish to compute a signed volume.

        Returns
        -------
        float
            The signed volume associated with the triangle.
        """
        v1 = self.vertices_[tri[0], :]
        v2 = self.vertices_[tri[1], :]
        v3 = self.vertices_[tri[2], :]

        volume = (1.0 / 6.0) * (v1.dot(np.cross(v2, v3)))
        return volume
项目:meshpy    作者:BerkeleyAutomation    | 项目源码 | 文件源码
def _area_of_tri(self, tri):
        """Return the area of the given triangle.

        Parameters
        ----------
        tri : :obj:`numpy.ndarray` of int
            The triangle for which we wish to compute an area.

        Returns
        -------
        float
            The area of the triangle.
        """
        verts = [self.vertices[i] for i in tri]
        ab = verts[1] - verts[0]
        ac = verts[2] - verts[0]
        return 0.5 * np.linalg.norm(np.cross(ab, ac))
项目:BlenderRobotDesigner    作者:HBPNeurorobotics    | 项目源码 | 文件源码
def vector_product(v0, v1, axis=0):
    """Return vector perpendicular to vectors.

    >>> v = vector_product([2, 0, 0], [0, 3, 0])
    >>> numpy.allclose(v, [0, 0, 6])
    True
    >>> v0 = [[2, 0, 0, 2], [0, 2, 0, 2], [0, 0, 2, 2]]
    >>> v1 = [[3], [0], [0]]
    >>> v = vector_product(v0, v1)
    >>> numpy.allclose(v, [[0, 0, 0, 0], [0, 0, 6, 6], [0, -6, 0, -6]])
    True
    >>> v0 = [[2, 0, 0], [2, 0, 0], [0, 2, 0], [2, 0, 0]]
    >>> v1 = [[0, 3, 0], [0, 0, 3], [0, 0, 3], [3, 3, 3]]
    >>> v = vector_product(v0, v1, axis=1)
    >>> numpy.allclose(v, [[0, 0, 6], [0, -6, 0], [6, 0, 0], [0, -6, 6]])
    True

    """
    return numpy.cross(v0, v1, axis=axis)
项目:SPIND    作者:LiuLab-CSRC    | 项目源码 | 文件源码
def calc_rotation_matrix(q1, q2, ref_q1, ref_q2):
  ref_nv = np.cross(ref_q1, ref_q2) 
  q_nv = np.cross(q1, q2)
  if min(norm(ref_nv), norm(q_nv)) == 0.:  # avoid 0 degree including angle
    return np.identity(3)
  axis = np.cross(ref_nv, q_nv)
  angle = rad2deg(acos(ref_nv.dot(q_nv) / (norm(ref_nv) * norm(q_nv))))
  R1 = axis_angle_to_rotation_matrix(axis, angle)
  rot_ref_q1, rot_ref_q2 = R1.dot(ref_q1), R1.dot(ref_q2)  # rotate ref_q1,2 plane to q1,2 plane

  cos1 = max(min(q1.dot(rot_ref_q1) / (norm(rot_ref_q1) * norm(q1)), 1.), -1.)  # avoid math domain error
  cos2 = max(min(q2.dot(rot_ref_q2) / (norm(rot_ref_q2) * norm(q2)), 1.), -1.)
  angle1 = rad2deg(acos(cos1))
  angle2 = rad2deg(acos(cos2))
  angle = (angle1 + angle2) / 2.
  axis = np.cross(rot_ref_q1, q1)
  R2 = axis_angle_to_rotation_matrix(axis, angle)

  R = R2.dot(R1)
  return R
项目:StructEngPy    作者:zhuoju36    | 项目源码 | 文件源码
def __init__(self,origin, pt1, pt2, name=None):
        """
        origin: 3x1 vector
        pt1: 3x1 vector
        pt2: 3x1 vector
        """
        self.__origin=origin    
        vec1 = np.array([pt1[0] - origin[0] , pt1[1] - origin[1] , pt1[2] - origin[2]])
        vec2 = np.array([pt2[0] - origin[0] , pt2[1] - origin[1] , pt2[2] - origin[2]])
        cos = np.dot(vec1, vec2)/np.linalg.norm(vec1)/np.linalg.norm(vec2)
        if  cos == 1 or cos == -1:
            raise Exception("Three points should not in a line!!")        
        self.__x = vec1/np.linalg.norm(vec1)
        z = np.cross(vec1, vec2)
        self.__z = z/np.linalg.norm(z)
        self.__y = np.cross(self.z, self.x)
        self.__name=uuid.uuid1() if name==None else name
项目:StructEngPy    作者:zhuoju36    | 项目源码 | 文件源码
def set_by_3pts(self,origin, pt1, pt2):
        """
        origin: tuple 3
        pt1: tuple 3
        pt2: tuple 3
        """
        self.origin=origin    
        vec1 = np.array([pt1[0] - origin[0] , pt1[1] - origin[1] , pt1[2] - origin[2]])
        vec2 = np.array([pt2[0] - origin[0] , pt2[1] - origin[1] , pt2[2] - origin[2]])
        cos = np.dot(vec1, vec2)/np.linalg.norm(vec1)/np.linalg.norm(vec2)
        if  cos == 1 or cos == -1:
            raise Exception("Three points should not in a line!!")        
        self.x = vec1/np.linalg.norm(vec1)
        z = np.cross(vec1, vec2)
        self.z = z/np.linalg.norm(z)
        self.y = np.cross(self.z, self.x)
项目:PyFRAP    作者:alexblaessle    | 项目源码 | 文件源码
def normalByCross(vec1,vec2):

    r"""Returns normalised normal vectors of plane spanned by two vectors.

    Normal vector is computed by:

    .. math:: \mathbf{n} = \frac{\mathbf{v_1} \times \mathbf{v_2}}{|\mathbf{v_1} \times \mathbf{v_2}|}

    .. note:: Will return zero vector if ``vec1`` and ``vec2`` are colinear.

    Args:
        vec1 (numpy.ndarray): Vector 1.
        vec2 (numpy.ndarray): Vector 2.

    Returns:
        numpy.ndarray: Normal vector.
    """

    if checkColinear(vec1,vec2):

        printWarning("Can't compute normal of vectors, they seem to be colinear. Returning zero.")
        return np.zeros(np.shape(vec1))

    return np.cross(vec1,vec2)/np.linalg.norm(np.cross(vec1,vec2))
项目:geomeTRIC    作者:leeping    | 项目源码 | 文件源码
def calc_e0(self):
        """
        Compute the reference axis for adding dummy atoms. 
        Only used in the case of linear molecules.

        We first find the Cartesian axis that is "most perpendicular" to the molecular axis.
        Next we take the cross product with the molecular axis to create a perpendicular vector.
        Finally, this perpendicular vector is normalized to make a unit vector.
        """
        ysel = self.x0[self.a, :]
        vy = ysel[-1]-ysel[0]
        ev = vy / np.linalg.norm(vy)
        # Cartesian axes.
        ex = np.array([1.0,0.0,0.0])
        ey = np.array([0.0,1.0,0.0])
        ez = np.array([0.0,0.0,1.0])
        self.e0 = np.cross(vy, [ex, ey, ez][np.argmin([np.dot(i, ev)**2 for i in [ex, ey, ez]])])
        self.e0 /= np.linalg.norm(self.e0)
项目:geomeTRIC    作者:leeping    | 项目源码 | 文件源码
def normal_vector(self, xyz):
        xyz = xyz.reshape(-1,3)
        a = np.array(self.a)
        b = self.b
        c = np.array(self.c)
        xyza = np.mean(xyz[a], axis=0)
        xyzc = np.mean(xyz[c], axis=0)
        # vector from first atom to central atom
        vector1 = xyza - xyz[b]
        # vector from last atom to central atom
        vector2 = xyzc - xyz[b]
        # norm of the two vectors
        norm1 = np.sqrt(np.sum(vector1**2))
        norm2 = np.sqrt(np.sum(vector2**2))
        crs = np.cross(vector1, vector2)
        crs /= np.linalg.norm(crs)
        return crs
项目:geomeTRIC    作者:leeping    | 项目源码 | 文件源码
def value(self, xyz):
        xyz = xyz.reshape(-1,3)
        a = np.array(self.a)
        b = self.b
        c = self.c
        d = np.array(self.d)
        xyza = np.mean(xyz[a], axis=0)
        xyzd = np.mean(xyz[d], axis=0)

        vec1 = xyz[b] - xyza
        vec2 = xyz[c] - xyz[b]
        vec3 = xyzd - xyz[c]
        cross1 = np.cross(vec2, vec3)
        cross2 = np.cross(vec1, vec2)
        arg1 = np.sum(np.multiply(vec1, cross1)) * \
               np.sqrt(np.sum(vec2**2))
        arg2 = np.sum(np.multiply(cross1, cross2))
        answer = np.arctan2(arg1, arg2)
        return answer
项目:geomeTRIC    作者:leeping    | 项目源码 | 文件源码
def value(self, xyz):
        xyz = xyz.reshape(-1,3)
        a = self.a
        b = self.b
        c = self.c
        d = self.d
        vec1 = xyz[b] - xyz[a]
        vec2 = xyz[c] - xyz[b]
        vec3 = xyz[d] - xyz[c]
        cross1 = np.cross(vec2, vec3)
        cross2 = np.cross(vec1, vec2)
        arg1 = np.sum(np.multiply(vec1, cross1)) * \
               np.sqrt(np.sum(vec2**2))
        arg2 = np.sum(np.multiply(cross1, cross2))
        answer = np.arctan2(arg1, arg2)
        return answer
项目:measure_lens_alignment    作者:oxford-pcs    | 项目源码 | 文件源码
def getProjectedAngleInXYPlane(self, z=0, ref_axis=[0,1], centre=[0,0], inDeg=True):    
    '''
      Project the OA vector to z=z, calculate the XY position, construct a 
      2D vector from [centre] to this XY and measure the angle subtended by 
      this vector from [ref_axis] (clockwise).
    '''
    ref_axis = np.array(ref_axis)
    centre = np.array(centre)

    point_vector_from_fit_centre = np.array(self.getXY(z=z)) - centre
    dotP = np.dot(ref_axis, point_vector_from_fit_centre)
    crossP = np.cross(ref_axis, point_vector_from_fit_centre)
    angle = np.arccos(dotP/(np.linalg.norm(ref_axis)*np.linalg.norm(point_vector_from_fit_centre)))

    if np.sign(crossP) > 0:
      angle = (np.pi-angle) + np.pi

    if inDeg:
      dir_v = self._eval_direction_vector()
      return np.degrees(angle)
    else:
      return angle
项目:hexmachina    作者:dnkrtz    | 项目源码 | 文件源码
def compute_normals(self):
        """Compute vertex and face normals of the triangular mesh."""

        # Compute face normals, easy as cake.
        for fi, face in enumerate(self.faces):
            self.face_normals[fi] = np.cross(self.vertices[face[2]] - self.vertices[face[0]],
                                             self.vertices[face[1]] - self.vertices[face[0]])

        # Next, compute the vertex normals.
        for fi, face in enumerate(self.faces):
            self.vertex_normals[face[0]] += self.face_normals[fi]
            self.vertex_normals[face[1]] += self.face_normals[fi]
            self.vertex_normals[face[2]] += self.face_normals[fi]

        # Normalize all vectors
        for i, f_norm in enumerate(self.face_normals):
            self.face_normals[i] = normalize(f_norm)
        for i, v_norm in enumerate(self.vertex_normals):
            self.vertex_normals[i] = normalize(v_norm)
项目:hexmachina    作者:dnkrtz    | 项目源码 | 文件源码
def rotate_coord_sys(old_u, old_v, new_norm):
        """Rotate a coordinate system to be perpendicular to the given normal."""
        new_u = old_u
        new_v = old_v
        old_norm = np.cross(old_u, old_v)
        # Project old normal onto new normal
        ndot = np.dot(old_norm, new_norm)
        # If projection is leq to -1, simply reverse
        if ndot <= -1:
            new_u = -new_u
            new_v = -new_v
            return new_u, new_v
        # Otherwise, compute new normal
        perp_old = new_norm - ndot * old_norm
        dperp = (old_norm + new_norm) / (1 + ndot)
        new_u -= dperp * np.dot(new_u, perp_old)
        new_v -= dperp * np.dot(new_v, perp_old)
        return new_u, new_v
项目:MDT    作者:cbclab    | 项目源码 | 文件源码
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
项目:MDT    作者:cbclab    | 项目源码 | 文件源码
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)
项目:MDT    作者:cbclab    | 项目源码 | 文件源码
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
项目:blmath    作者:bodylabs    | 项目源码 | 文件源码
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)
项目:blmath    作者:bodylabs    | 项目源码 | 文件源码
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])
项目:blmath    作者:bodylabs    | 项目源码 | 文件源码
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)
项目:blmath    作者:bodylabs    | 项目源码 | 文件源码
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)
项目:crazyswarm    作者:USC-ACTLab    | 项目源码 | 文件源码
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
项目:Acoustics    作者:fei0324    | 项目源码 | 文件源码
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
项目:Acoustics    作者:fei0324    | 项目源码 | 文件源码
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
项目:albion    作者:Oslandia    | 项目源码 | 文件源码
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)
项目:roboschool    作者:openai    | 项目源码 | 文件源码
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]) ) )
项目:krpcScripts    作者:jwvanderbeck    | 项目源码 | 文件源码
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))
项目:yt    作者:yt-project    | 项目源码 | 文件源码
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)
项目:scipy-lecture-notes-zh-CN    作者:jayleicn    | 项目源码 | 文件源码
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
项目:scipy-lecture-notes-zh-CN    作者:jayleicn    | 项目源码 | 文件源码
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
项目:Splipy    作者:sintefmath    | 项目源码 | 文件源码
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)
项目:Splipy    作者:sintefmath    | 项目源码 | 文件源码
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
项目:kaggle_dsb2017    作者:astoc    | 项目源码 | 文件源码
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
项目:pybot    作者:spillai    | 项目源码 | 文件源码
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)