Triangulate a set of points with a concave domain

Here is some Python code that does what you want.

First, building the alpha shape (see my previous answer):

def alpha_shape(points, alpha, only_outer=True):
    """
    Compute the alpha shape (concave hull) of a set of points.
    :param points: np.array of shape (n,2) points.
    :param alpha: alpha value.
    :param only_outer: boolean value to specify if we keep only the outer border or also inner edges.
    :return: set of (i,j) pairs representing edges of the alpha-shape. (i,j) are the indices in the points array.
    """
    assert points.shape[0] > 3, "Need at least four points"

    def add_edge(edges, i, j):
        """
        Add a line between the i-th and j-th points,
        if not in the list already
        """
        if (i, j) in edges or (j, i) in edges:
            # already added
            assert (j, i) in edges, "Can't go twice over same directed edge right?"
            if only_outer:
                # if both neighboring triangles are in shape, it's not a boundary edge
                edges.remove((j, i))
            return
        edges.add((i, j))

    tri = Delaunay(points)
    edges = set()
    # Loop over triangles:
    # ia, ib, ic = indices of corner points of the triangle
    for ia, ib, ic in tri.vertices:
        pa = points[ia]
        pb = points[ib]
        pc = points[ic]
        # Computing radius of triangle circumcircle
        # www.mathalino.com/reviewer/derivation-of-formulas/derivation-of-formula-for-radius-of-circumcircle
        a = np.sqrt((pa[0] - pb[0]) ** 2 + (pa[1] - pb[1]) ** 2)
        b = np.sqrt((pb[0] - pc[0]) ** 2 + (pb[1] - pc[1]) ** 2)
        c = np.sqrt((pc[0] - pa[0]) ** 2 + (pc[1] - pa[1]) ** 2)
        s = (a + b + c) / 2.0
        area = np.sqrt(s * (s - a) * (s - b) * (s - c))
        circum_r = a * b * c / (4.0 * area)
        if circum_r < alpha:
            add_edge(edges, ia, ib)
            add_edge(edges, ib, ic)
            add_edge(edges, ic, ia)
    return edges

To compute the edges of the outer boundary of the alpha shape use the following example call:

edges = alpha_shape(points, alpha=alpha_value, only_outer=True)

Now, after the edges of the outer boundary of the alpha-shape of points have been computed, the following function will determine whether a point (x,y) is inside the outer boundary.

def is_inside(x, y, points, edges, eps=1.0e-10):
    intersection_counter = 0
    for i, j in edges:
        assert abs((points[i,1]-y)*(points[j,1]-y)) > eps, 'Need to handle these end cases separately'
        y_in_edge_domain = ((points[i,1]-y)*(points[j,1]-y) < 0)
        if y_in_edge_domain:
            upper_ind, lower_ind = (i,j) if (points[i,1]-y) > 0 else (j,i)
            upper_x = points[upper_ind, 0] 
            upper_y = points[upper_ind, 1]
            lower_x = points[lower_ind, 0] 
            lower_y = points[lower_ind, 1]

            # is_left_turn predicate is evaluated with: sign(cross_product(upper-lower, p-lower))
            cross_prod = (upper_x - lower_x)*(y-lower_y) - (upper_y - lower_y)*(x-lower_x)
            assert abs(cross_prod) > eps, 'Need to handle these end cases separately'
            point_is_left_of_segment = (cross_prod > 0.0)
            if point_is_left_of_segment:
                intersection_counter = intersection_counter + 1
    return (intersection_counter % 2) != 0

On the input shown in the above figure (taken from my previous answer) the call is_inside(1.5, 0.0, points, edges) will return True, whereas is_inside(1.5, 3.0, points, edges) will return False.

Note that the is_inside function above does not handle degenerate cases. I added two assertions to detect such cases (you can define any epsilon value that fits your application). In many applications this is sufficient, but if not and you encounter these end cases, they need to be handled separately. See, for example, here on robustness and precision issues when implementing geometric algorithms.

One of Classic DT algorithms generates first a bounding triangle, then adds all new triangles sorted by x, then prunes out all triangles having a vertex in the supertriangle.

At least from the provided image one can derive the heuristics of pruning out also some triangles having all vertices on the concave hull. Without a proof, the triangles to be pruned out have a negative area when their vertices are sorted in the same order as the concave hull is defined.

This may need the concave hull to be inserted as well, and to be pruned out.