def f(a, b):
    return a + b
a = np.array([1, 3, 9, 15])
b = np.array([6, 18, 32])
f_ab = f(a[:, (np.newaxis)], b)
f_ab = np.empty((len(a), len(b)), dtype=np.array(f(a[0], b[0])).dtype)
for i, a_ in enumerate(a):
    for j, b_ in enumerate(b):
        f_ab[i, j] = f(a_, b_)
Лучшая структура данных для 3D-данных?