python类cross()的实例源码

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)

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

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

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'")

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'")

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)

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

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)

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

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

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

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)

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 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)

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

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, :]

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

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))

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)

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

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

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)

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))

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)

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

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

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

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

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)

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