def METHOD_NAME(tknots, ts, points, quats):
    base_indexes = tknots.searchsorted(ts, side='left')-1
    t1s, t2s = tknots[base_indexes], tknots[base_indexes+1]
    dts = (ts - t1s) / (t2s - t1s)
    #dts = np.clip(dts, 0.0, 1.0) # Just in case...
    matrix_out = []
    # TODO: ideally this should be vectorized with numpy;
    # but that would require implementation of quaternion
    # interpolation in numpy.
    for dt, base_index, point in zip(dts, base_indexes, points):
        q1, q2 = quats[base_index], quats[base_index+1]
        # spherical linear interpolation.
        # TODO: implement `squad`.
        if dt < 0:
            q = q1
        elif dt > 1.0:
            q = q2
        else:
            q = q1.slerp(q2, dt)
        matrix = q.to_matrix().to_4x4()
        matrix.translation = Vector(point)
        matrix_out.append(matrix)
    return matrix_out