Well, you could triangulate the boundary, do an intersection test of all the grid points against the triangles, remove those which intersect then retriangulate with all remaining points?
Here I only know the grid points And the boundary points to be subtracted. Boundary points may not fall within the grid points. Then with which coordintes I need to construct the delaunay triangulation? Could you explain?
Wouldn't be particularly efficient but should work I guess...
If using python, Shapely is a nice library for stuff like this.