Sample a 2D matrix over a line (non necessarily horizontal or vertical)

Question

TL;DR How to find the matrix coefficients sampled "over a line" with Numpy? (nearest neighbour)

Is there a ready-to-use Numpy/Scipy function for this (maybe Bresenham)?
I'm not 100% sure the following naive version works in all cases.

import numpy as np

a = np.arange(36).reshape((6, 6))
# [[ 0  1  2  3  4  5]
#  [ 6  7  8  9 10 11]
#  [12 13 14 15 16 17]
#  [18 19 20 21 22 23]
#  [24 25 26 27 28 29]
#  [30 31 32 33 34 35]]

def sample_over_line(a, i0, j0, i1, j1):
    di = i1 - i0
    dj = j1 - j0
    if abs(di) > abs(dj):
        if i0 > i1:
            i0, j0, i1, j1 = i1, j1, i0, j0
        b = [a[i, round((j1 - j0) / (i1 - i0) * (i - i0) + j0)] for i in range(i0, i1 + 1)]
    else:
        if j0 > j1:
            i0, j0, i1, j1 = i1, j1, i0, j0
        b = [a[round((i1 - i0) / (j1 - j0) * (j - j0) + i0), j] for j in range(j0, j1 + 1)]
    return b

print(sample_over_line(a=a, i0=1, j0=0, i1=3, j1=5))  # line from [1,0] to [3,5]
# [6, 7, 14, 15, 22, 23]

print(sample_over_line(a=a, i0=1, j0=2, i1=5, j1=3))  # line from [1,2] to [5,3]
# [8, 14, 20, 27, 33]

print(sample_over_line(a=a, i0=5, j0=3, i1=1, j1=2)) # line from [5,3] to [1,2]
# [8, 14, 20, 27, 33]