Python Binning using triangular bins

I'm trying to find a simple python module/package that has implemented 2D triangular bins so that it can be use in a similar fashion to scipy binned_statistic_dd. Is anyone aware of such a tool? I've searched but not found anything: the closest I've found is matplotlib's hexbin.

If I have to create a home-made solution, generating the vertex points for the triangular grid is easy, but how would you efficiently (need to avoid slow loops if possible as datasets are about 100K points) search which triangle a point lies in?

import matplotlib.pyplot as plt
import matplotlib.tri as tri
import numpy as np
def plot_triangular_bin_freq(x,y,Vx,Vy):
    X, Y = np.meshgrid(x, y)
    Ny, Nx = X.shape
    
    iy,ix = np.indices((Ny-1, Nx-1))
    # max vertice is supposed to be 
    # max(iy)*Nx + max(ix) + (Nx+1)
    # = (Ny-2)*Nx + (Nx-2) + (Nx+1)
    # = Ny * Nx - 1
    assert iy.max() == Ny-2
    assert ix.max() == Nx-2
    
    # build square grid and split it in a lower-left, upper-right triangles
    # and construct the triangulation
    vertices = (((iy * Nx) + ix)[:,:,None] + np.array([0,1,Nx,Nx,Nx+1,1])[None,None,:]).reshape(-1, 3)
    triangles = tri.Triangulation(X.flatten(), Y.flatten(), vertices)
    
    # Normalized point coordinates
    Vx = (np.asarray(Vx).flatten() - x[0]) * ((Nx-1) / (x[-1] - x[0]))
    Vy = (np.asarray(Vy).flatten() - y[0]) * ((Ny-1) / (y[-1] - y[0]))
    
    m = (0 <= Vx) & (Vx < Nx-1) & (0 <= Vy) & (Vy < Ny-1)
    
    # get indices on the x,y boxes
    Ix, Rx = divmod(Vx[m], 1)
    Iy, Ry = divmod(Vy[m], 1)
    
    # (Rx+Ry)=1 is the boundary between the two triangles
    # w indicates the index of the triangle where the point lies on
    w = ((Rx+Ry)>=1) +  2*(Ix + (Nx-1)*Iy)
    assert max(Ix) < Nx-1
    assert max(Iy) < Ny-1
    assert max(Ix + Iy*(Nx-1)) < (Nx-1)*(Ny-1)
    
    # z[i] is the number of points that lies inside z[i]
    z = np.bincount(w.astype(np.int64), minlength=2*(Nx-1)*(Ny-1))
    plt.tripcolor(triangles, z, shading='flat')

x = np.arange(15)/2.
y = np.arange(10)/2.
Vx = np.random.randn(1000) + 3
Vy = np.random.randn(1000) + 1

plot_triangular_bin_freq(x,y,Vx,Vy)