U235U235

mat_A,_=np.mgrid[0:100:5, 0:100:5]

Еще вариант без for:

b=np.array([256*256,256,1], dtype=np.uint32) # or [1, 256, 256*256] BGR<->RGB
result=image@b

Можно так: arr[bool_mask]=42
Старайтесь не использовать циклы for, если работаете с массивами.

Можно использовать meshgrid и reshape. Что-то подобное:

import numpy as np
x = [1, 2, 3]
y = [10, 20, 30]
z = [100, 200, 300]
res1=np.array(np.meshgrid(x,y)).T.reshape(-1,3,2)
print(res1)
res2=np.array(np.meshgrid(x,y,z)).T.reshape(-1,3,3)
print(res2)

Если у вас уже скелет есть, то, самое простое, это взять среднее по окрестности 3x3, а дальше порогом выделять разветления или конечные точки.

В scipy есть scipy.spatial.distance.pdist¶
По-моему это то, что вам нужно.
Поправлю: cdist, конечно же.

Простите, но если умножить вектор на транспонированый вектор обратных значений разве результатом не будет N - число элементов в векторе?

Зачем Вы присваиваете grid случайные числа, а затем меняете их на default_angle?

grid=np.ndarray((num_columns, num_rows))
grid[:,:]=math.pi * 0.25

и цикл не нужен.