import torch


def laplacian_beltrami_matrix(vertices, faces):
    """
    Return the Laplacian-Beltrami matrix
    See https://en.wikipedia.org/wiki/Discrete_Laplace_operator#Mesh_Laplacians
    How we compute cotangent :
        -> Cross product ABxAC = AB*AC*sin(BAC)
        -> Dot product AB*AC = AB*AC*cos(BAC)
        -> cot(BAC) = dot_product(AB, AC)/cross_product(AB, AC)
    We compute the 3 cotangents of each faces. Then to sum opposing angles of each edge, we use the
    properties of sparse tensors to sum overlapping values with coalesce()
    :param vertices: torch.float, [V, 3]
    :param faces: torch.long, [F, 3]
    :return: torch.float, [V, V] sparse COO tensor
    """
    # First we get the cotangents
    pts = vertices[faces]
    vectors = pts[:, [1, 2, 0]] - pts  # [V, 3] vectors 01, 21, 20
    dot = torch.bmm(vectors.view(-1, 3).unsqueeze(1),
                    -vectors[:, [2, 0, 1]].view(-1, 3).unsqueeze(2)
                    ).view(faces.shape[0], 3)
    len_cross = torch.cross(vectors[:, 0], -vectors[:, 1], dim=1).square().sum(1).sqrt().expand(3, -1).T
    half_cot_alpha = dot / len_cross / 2.  # cotangent / 2

    # Then we get their position in the matrix
    indices = faces[:, [[1, 2], [2, 0], [0, 1]]].view(-1, 2).T  # Vertices i and j for i != j
    ind_diag = torch.cat([faces[:, [1, 2, 0]].view(-1).expand(2, -1),
                          faces[:, [2, 0, 1]].view(-1).expand(2, -1)], dim=1)  # Vertices i and j for i == j
    indices = torch.cat([indices, indices[[1, 0], :], ind_diag], dim=1)
    values = half_cot_alpha.view(-1).repeat(2)  # Values for i != j
    values_diag = - half_cot_alpha.view(-1).repeat(2)  # Values for i == j
    values = torch.cat([values, values_diag], dim=0)

    # We fill the sparse matrix, and add overlapping values
    laplacian = torch.sparse_coo_tensor(indices=indices, values=values, size=[vertices.shape[0], vertices.shape[0]])
    laplacian = laplacian.coalesce()

    return laplacian