xright = max(ax, bx, cx)
ytop = min(ay, by, cy)
ybottom = max(ay, by, cy)
barytransform = np.linalg.inv([[ax, bx, cx], [ay, by, cy], [1, 1, 1]])
grid = np.mgrid[xleft:xright, ytop:ybottom].reshape(2, -1)
grid = np.vstack((grid, np.ones((1, grid.shape[1]))))
barycoords = np.dot(barytransform, grid)
barycoords = barycoords[:, (np.all(barycoords >= 0, axis=0))]
Increasing efficiency of barycentric coordinate calculation in python