You can keep using the 2D Delaunay algorithm, but do it between two neighboring layers (z1, z2) each time. (z1 < z2)
Assume contour lines are sampled and stored as a set of (x,y,z) in counter-clockwise order on z-plane. You need to construct a set of boundary/hole points for triangulation:
- detect or compute overlapping line segments (x1, y1, z1)(x2, y2, z1), and (x1, y1, z2)(x2, y2, z2)
- move (x1, y1, z1) inwards by a small offset, along the direction of average normal of two neighboring line segments. Similarly, move (x2, y2, z1) inwards; move (x1, y1, z2)(x2,y2,z2) outwards. Now the two (or more) contour lines are separated, and you get two (or more) boundary/hole point set.
- apply Delaunay trangulation within the area defined by boundary and holes (i.e. generate random points, drop points outside boundary or inside holes, construct triangle). I assume Poly2Tri will handle the cases where a triangle partially cover a hole. Once you have all triangles in (x,y) plane, restore the offset and compute z by interpolating original data.
Do triangulation again for (z2, z3) layer, so on and so forth. Note that the hole/boundary points at z2 remains the same, but if overlapping happens, instead of moving inwards, it's now moving outwards. At the end, combine all triangles into one mesh.