Calculate vertex distances of a mesh

I am using Numpy arrays to express a triangular mesh.

I have two matrices: coordinates is a 3 x n matrix, and connectivity is an n x n matrix that uses 0s and 1s to store vertex connectivity.

Now I want to calculate a n x n matrix named distances that stores vertex distances. Only the positions where connectivity[i,j] == 1 are calculated. Anywhere else are not.

What is the most elegant way to calculate this in Python?

For example, I have a mesh of 4 vertices like this:

Then

import numpy as np

coordinates = np.array(
    [
        [-1, -1, 0],    # A
        [1, -1, 0],     # B
        [1, 1, 0],      # C
        [-1, 1, 0]      # D
    ], 
    dtype=np.float32
)

connectivity = np.array(
    [
        [0, 1, 1, 1],    # A-B, A-C, A-D
        [1, 0, 1, 0],    # B-A, B-C
        [1, 1, 0, 1],    # C-A, C-B, C-D
        [1, 0, 1, 0],    # D-A, D-C
    ],
    dtype=np.int32
)

# For this example, expected `distances` is like this
distances = np.array(
    [
        [0, 2, 2.828, 2],   # A-B, A-C, A-D
        [2, 0, 2, 0],       # B-A, B-C
        [2.828, 2, 0, 2],   # C-A, C-B, C-D
        [2, 0, 2, 0],       # D-A, D-C
    ],
    dtype=np.float32
)

Here is a solution using np.tile to efficiently generate all pairs of distances, which are then masked out by the connectivity matrix.

The idea is to generate two matrices containing the tiled vertex coordinates, in the following pattern:

ABCD 
ABCD
ABCD
ABCD

and

AAAA
BBBB
CCCC
DDDD

The difference of those two matrices then gives you all combinations. If you only want the entries that actually are connected, then you can simply multiply these matrices by the connectivity matrix.

mat = np.tile(coordinates, (4,1,1)) * connectivity[..., None]
mat_t = np.transpose(mat, axes=(1,0,2))

distances = np.sqrt(((mat - mat_t)**2).sum(axis=2))

Which gives you the expected result:

array([[0.      , 2.      , 2.828427, 2.      ],
       [2.      , 0.      , 2.      , 0.      ],
       [2.828427, 2.      , 0.      , 2.      ],
       [2.      , 0.      , 2.      , 0.      ]], dtype=float32)

Alternatively, assuming you generated the connectivity matrix from a triangle index array, you could use it to generate the distances, and then generate the matrix in a similar fashion as you generate the connectivity matrix.

Compute distance between coordinates with scipy.spatial.distance.cdist
and ensure it with connectivity array at the end:

from scipy.spatial.distance import cdist

dist = cdist(coordinates, coordinates, 'euclidean') * connectivity
print(dist) 

[[0.         2.         2.82842712 2.        ]
 [2.         0.         2.         0.        ]
 [2.82842712 2.         0.         2.        ]
 [2.         0.         2.         0.        ]]