def align_vectors(a, b, return_angle=False):
    """
    Find a transform between two 3D vectors.

    Implements the method described here:
    http://ethaneade.com/rot_between_vectors.pdf

    Parameters
    --------------
    a : (3,) float
      Source vector
    b : (3,) float
      Target vector
    return_angle : bool
      If True return the angle between the two vectors

    Returns
    -------------
    transform : (4, 4) float
      Homogenous transform from a to b
    angle : float
      Angle between vectors in radians
      Only returned if return_angle
    """
    # copy of input vectors
    a = np.array(a, dtype=np.float64, copy=True)
    b = np.array(b, dtype=np.float64, copy=True)

    # make sure vectors are 3D
    if a.shape != (3,) or b.shape != (3,):
        raise ValueError('only works for (3,) vectors')

    # unitize input vectors
    a /= np.linalg.norm(a)
    b /= np.linalg.norm(b)

    # projection of a onto b
    dot = np.dot(a, b)

    # are vectors just reversed
    if dot < (tol.zero - 1):
        # a reversed vector is 180 degrees
        angle = np.pi

        # get an arbitrary perpendicular vector to a
        perp = util.generate_basis(a)[0] * np.eye(3)

        # (3, 3) rotation from a to b
        rotation = (2 * np.dot(perp, perp.T)) - np.eye(3)

    # are vectors already the same
    elif dot > (1 - tol.zero):
        angle = 0.0
        # no rotation
        rotation = np.eye(3)

    # vectors are at some angle to each other
    else:
        # we already handled values out of the range [-1.0, 1.0]
        angle = np.arccos(dot)

        # (3,) vector perpendicular to both a and b
        w = np.cross(a, b)

        # a float between 0.5 and 1.0
        c = 1.0 / (1.0 + dot)

        # (3, 3) skew- symmetric matrix from the (3,) vector w
        # the matrix has the property: wx == -wx.T
        wx = np.array([[0, -w[2], w[1]],
                       [w[2], 0, -w[0]],
                       [-w[1], w[0], 0]])

        # (3, 3) rotation from a to b
        rotation = np.eye(3) + wx + (np.dot(wx, wx) * c)

    # put rotation into homogenous transformation matrix
    transform = np.eye(4)
    transform[:3, :3] = rotation

    if return_angle:
        return transform, angle

    return transform
Find a transform between two 3D vectors.

    Implements the method described here:
    http://ethaneade.com/rot_between_vectors.pdf

    Parameters
    --------------
    a : (3,) float
      Source vector
    b : (3,) float
      Target vector
    return_angle : bool
      If True return the angle between the two vectors

    Returns
    -------------
    transform : (4, 4) float
      Homogenous transform from a to b
    angle : float
      Angle between vectors in radians
      Only returned if return_angle