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)