How to use the trimesh.util.diagonal_dot function in trimesh

distance : (n,) float
      [OPTIONAL] Distance from point
    denom : (n,) float
      [OPTIONAL] Denominator
    """

    # check input types
    plane_origins = np.asanyarray(plane_origins, dtype=np.float64)
    plane_normals = np.asanyarray(plane_normals, dtype=np.float64)
    line_origins = np.asanyarray(line_origins, dtype=np.float64)
    line_directions = np.asanyarray(line_directions, dtype=np.float64)

    # vector from line to plane
    origin_vectors = plane_origins - line_origins

    projection_ori = util.diagonal_dot(origin_vectors, plane_normals)
    projection_dir = util.diagonal_dot(line_directions, plane_normals)

    valid = np.abs(projection_dir) > 1e-5

    distance = np.divide(projection_ori[valid],
                         projection_dir[valid])

    on_plane = line_directions[valid] * distance.reshape((-1, 1))
    on_plane += line_origins[valid]

    result = [on_plane, valid]

    if return_distance:
        result.append(distance)
    if return_denom:
        result.append(projection_dir)

edge1 = vert2 - vert0

    # P is a vector perpendicular to the ray direction and one
    # triangle edge.
    P = np.cross(ray_direction, edge1)

    # if determinant is near zero, ray lies in plane of triangle
    det = diagonal_dot(edge0, P)
    candidates[np.abs(det) < tol.zero] = False

    if not candidates.any():
        return candidates
    # remove previously calculated terms which are no longer candidates
    inv_det = 1.0 / det[candidates]
    T = ray_origin - vert0[candidates]
    u = diagonal_dot(T, P[candidates]) * inv_det

    new_candidates = np.logical_not(np.logical_or(u < -tol.zero,
                                                  u > (1 + tol.zero)))
    candidates[candidates] = new_candidates
    if not candidates.any():
        return candidates
    inv_det = inv_det[new_candidates]
    T = T[new_candidates]
    u = u[new_candidates]

    Q = np.cross(T, edge0[candidates])
    v = np.dot(ray_direction, Q.T) * inv_det

    new_candidates = np.logical_not(np.logical_or((v < -tol.zero),
                                                  (u + v > (1 + tol.zero))))

u = u[new_candidates]

    Q = np.cross(T, edge0[candidates])
    v = np.dot(ray_direction, Q.T) * inv_det

    new_candidates = np.logical_not(np.logical_or((v < -tol.zero),
                                                  (u + v > (1 + tol.zero))))

    candidates[candidates] = new_candidates
    if not candidates.any():
        return candidates

    Q = Q[new_candidates]
    inv_det = inv_det[new_candidates]

    t = diagonal_dot(edge1[candidates], Q) * inv_det
    candidates[candidates] = t > tol.zero

    return candidates

points = np.vstack((points[0], points[1:][direction_ok]))
    direction = direction[direction_ok]
    direction_norm = direction_norm[direction_ok]

    # create a vector between every other point, then turn it perpendicular
    # if we have points A B C D
    # and direction vectors A-B, B-C, etc
    # these will be perpendicular to the vectors A-C, B-D, etc
    perp = (points[2:] - points[:-2]).T[::-1].T
    perp_norm = util.row_norm(perp)
    perp_nonzero = perp_norm > tol.merge
    perp[perp_nonzero] /= perp_norm[perp_nonzero].reshape((-1, 1))

    # find the projection of each direction vector
    # onto the perpendicular vector
    projection = np.abs(util.diagonal_dot(perp,
                                          direction[:-1]))

    projection_ratio = np.max((projection / direction_norm[1:],
                               projection / direction_norm[:-1]), axis=0)

    mask = np.ones(len(points), dtype=np.bool)
    # since we took diff, we need to offset by one
    mask[1:-1][projection_ratio < 1e-4 * scale] = False

    merged = points[mask]
    return merged

'''
    candidates = np.ones(len(triangles), dtype=np.bool)

    # edge vectors and vertex locations in (n,3) format
    vert0 = triangles[:, 0, :]
    vert1 = triangles[:, 1, :]
    vert2 = triangles[:, 2, :]
    edge0 = vert1 - vert0
    edge1 = vert2 - vert0

    # P is a vector perpendicular to the ray direction and one
    # triangle edge.
    P = np.cross(ray_direction, edge1)

    # if determinant is near zero, ray lies in plane of triangle
    det = diagonal_dot(edge0, P)
    candidates[np.abs(det) < tol.zero] = False

    if not candidates.any():
        return candidates
    # remove previously calculated terms which are no longer candidates
    inv_det = 1.0 / det[candidates]
    T = ray_origin - vert0[candidates]
    u = diagonal_dot(T, P[candidates]) * inv_det

    new_candidates = np.logical_not(np.logical_or(u < -tol.zero,
                                                  u > (1 + tol.zero)))
    candidates[candidates] = new_candidates
    if not candidates.any():
        return candidates
    inv_det = inv_det[new_candidates]
    T = T[new_candidates]

# each triangle area and mean center
    triangles_area = triangles.area(crosses=crosses, sum=False)
    triangles_center = vertices[faces].mean(axis=1)

    # since the convex hull is (hopefully) convex, the vector from
    # the centroid to the center of each face
    # should have a positive dot product with the normal of that face
    # if it doesn't it is probably backwards
    # note that this sometimes gets screwed up by precision issues
    centroid = np.average(triangles_center,
                          weights=triangles_area,
                          axis=0)
    # a vector from the centroid to a point on each face
    test_vector = triangles_center - centroid
    # check the projection against face normals
    backwards = util.diagonal_dot(normals,
                                  test_vector) < 0.0

    # flip the winding outward facing
    faces[backwards] = np.fliplr(faces[backwards])
    # flip the normal
    normals[backwards] *= -1.0

    # save the work we did to the cache so it doesn't have to be recomputed
    initial_cache = {'triangles_cross': crosses,
                     'triangles_center': triangles_center,
                     'area_faces': triangles_area,
                     'centroid': centroid}

    # create the Trimesh object for the convex hull
    convex = Trimesh(vertices=vertices,
                     faces=faces,

def method_cramer():
        dot00 = util.diagonal_dot(edge_vectors[:, 0], edge_vectors[:, 0])
        dot01 = util.diagonal_dot(edge_vectors[:, 0], edge_vectors[:, 1])
        dot02 = util.diagonal_dot(edge_vectors[:, 0], w)
        dot11 = util.diagonal_dot(edge_vectors[:, 1], edge_vectors[:, 1])
        dot12 = util.diagonal_dot(edge_vectors[:, 1], w)

        inverse_denominator = 1.0 / (dot00 * dot11 - dot01 * dot01)

        barycentric = np.zeros((len(triangles), 3), dtype=np.float64)
        barycentric[:, 2] = (dot00 * dot12 - dot01 *
                             dot02) * inverse_denominator
        barycentric[:, 1] = (dot11 * dot02 - dot01 *
                             dot12) * inverse_denominator
        barycentric[:, 0] = 1 - barycentric[:, 1] - barycentric[:, 2]
        return barycentric

def method_cramer():
        dot00 = util.diagonal_dot(edge_vectors[:, 0], edge_vectors[:, 0])
        dot01 = util.diagonal_dot(edge_vectors[:, 0], edge_vectors[:, 1])
        dot02 = util.diagonal_dot(edge_vectors[:, 0], w)
        dot11 = util.diagonal_dot(edge_vectors[:, 1], edge_vectors[:, 1])
        dot12 = util.diagonal_dot(edge_vectors[:, 1], w)

        inverse_denominator = 1.0 / (dot00 * dot11 - dot01 * dot01)

        barycentric = np.zeros((len(triangles), 3), dtype=np.float64)
        barycentric[:, 2] = (dot00 * dot12 - dot01 *
                             dot02) * inverse_denominator
        barycentric[:, 1] = (dot11 * dot02 - dot01 *
                             dot12) * inverse_denominator
        barycentric[:, 0] = 1 - barycentric[:, 1] - barycentric[:, 2]
        return barycentric

# (2, ) int, list of 2 closest candidate indices
            idxs = distance.argsort()[:2]
            # make sure the two distances are identical
            check_distance = distance[idxs].ptp() < tol.merge
            # make sure the magnitude of both distances are nonzero
            check_magnitude = (np.abs(distance[idxs]) > tol.merge).all()

            # check if query-points are actually off-surface
            if check_distance and check_magnitude:
                # get face normals for two points
                normals = mesh.face_normals[np.array(candidate)[idxs]]
                # compute normalized surface-point to query-point vectors
                vectors = ((points[i] - close_points[idxs]) /
                           distance[idxs, np.newaxis] ** 0.5)
                # compare enclosed angle for both face normals
                dots = util.diagonal_dot(normals, vectors)
                # take the idx with the most positive angle
                idx = idxs[dots.argmax()]

        # take the single closest value from the group of values
        result_close[i] = close_points[idx]
        result_tid[i] = candidate[idx]
        result_distance[i] = distance[idx]

    # we were comparing the distance squared so
    # now take the square root in one vectorized operation
    result_distance **= .5

    return result_close, result_distance, result_tid

[OPTIONAL] Distance from point
    denom : (n,) float
      [OPTIONAL] Denominator
    """

    # check input types
    plane_origins = np.asanyarray(plane_origins, dtype=np.float64)
    plane_normals = np.asanyarray(plane_normals, dtype=np.float64)
    line_origins = np.asanyarray(line_origins, dtype=np.float64)
    line_directions = np.asanyarray(line_directions, dtype=np.float64)

    # vector from line to plane
    origin_vectors = plane_origins - line_origins

    projection_ori = util.diagonal_dot(origin_vectors, plane_normals)
    projection_dir = util.diagonal_dot(line_directions, plane_normals)

    valid = np.abs(projection_dir) > 1e-5

    distance = np.divide(projection_ori[valid],
                         projection_dir[valid])

    on_plane = line_directions[valid] * distance.reshape((-1, 1))
    on_plane += line_origins[valid]

    result = [on_plane, valid]

    if return_distance:
        result.append(distance)
    if return_denom:
        result.append(projection_dir)