Is there an easy way to evaluate an outer product with an "arbitrary function" in python or its standard libraries?

Question

I have a 20x20 grid of pixels. Each of these pixels has coordinates (i,j), where i and j both go from 0 to 19.

I would like to make a 400x400 similarity matrix C of these pixels, where if pixel p has coordinates (i,j) and pixel q has coordinates (k,l), then C[p,q] has the similarity score of p and q.

For my similarity score, I'm using the rbf kernel.

It feels like there should be some clever way using broadcasting to create this matrix C without using for loops.

For example, if pixels was a 400x2 array consisting of the coordinates of the 400 pixels, then the "outer product" pixels @ pixels.T would have dimension 400x400 (thought it would consist of the wrong values).

The only different between pixels @ pixels.T and what I would like to do is the similarity function I'm applying to the different rows of pixels. In the matrix multiplication case, it's an inner product. I would just like to apply the rbf kernel (which takes in two vectors, and outputs a scalar).

Is there a particularly elegant or efficient way to do this in python or its core libraries (numpy, scipy)?

Here's what I'm currently doing, but the for loops annoy me:

def sqeuclidean_distance(x, y):
    return np.sum((x-y)**2)

# RBF Kernel
def rbf_kernel(x, y, length_scale=1e1):
    return np.exp( -sqeuclidean_distance(x, y) / (length_scale ** 2))

def make_cov_mat(d):
    x = np.arange(d)
    pixels = np.array([e for e in itertools.product(*[x,x])])
    cov_mat = np.zeros((d**2, d**2))
    for i in range(d**2):
        for j in range(d**2):
            cov_mat[i,j] = rbf_kernel(pixels[i, :], pixels[j, :])
    return cov_mat