basis = np.array([[a1, a2, a3],
                  [b1, b2, b3],
                  [c1, c2, c3]])

for i in range(vox.shape[0]):
    for j in range(vox.shape[1]):
        for k in range(vox.shape[2]):
            p = np.dot(basis, np.array([i, j, k]))

if np.dot(p, q) > 0.0:
                desired_vox_indices.append([i, j, k])