Rewritten answer, didn't understand the question correctly the first time around.
Your goal should be to divide the plane into polygons in such a way that above each polygon, there is a constant number of surfaces, i.e. at every point within a polygon, a vertical ray will intersect the same number of triangles. Then you can simply alternate between “hole” and “solid”, and in that way sum up the holes. In order to obtain this division of the plane, look at “vertical” parts of the mesh. On the one hand, this includes triangles which are actually vertical. On the other hand, this also includes edges between triangles if the normal of one triangle is below horizontal and the other above.
So I'd suggest the following approach:
- Identify all “vertical” components of the mesh, and project them into the ground plane.
- Combine these line segments to obtain a partition of the plane into polygons.
- Iterate over all triangles and assign each to one of these polygons. Polygons at the boundary between polygons will be divided into several triangles.
- For the triangles of each polygon of the partition, identify the different layers. This can be done by simply following adjacency between triangles, with the exception of vertical components.
- Order the layers for each polygon, e.g. by intersecting a single ray with them. Then compute volumnes between layers and the ground layer using the shoelace formula you already referenced.
- Compute the correct sums and differences: the first layer is positive, the second negative, and so on, up to the last layer which does not count at all. Do so for all polygons, and sum the results. Edit: Thinking about it, triangle orientations will already take care of the alternating signs, so disregarding the top layer is the crucial point here.
I don't have any ready ideas for dealing with the slightly slanted scenarios you mention. Looking at the holes found by the above might help. You can combine them accross boundaries of the partition, then investigate them one at a time to see whether they can be made any smaller.
Edit: After some more thought, I realized that since you only have to disregard the topmost layer, you could also tackle this aspect from a different angle: Iterate over all triangles in order of decreasing height (of any vertex, or sume of three vertices, or whatever; details don't matter if the mesh doesn't intersect itself). While iterating, maintain a 2d idea of what areas above ground are already below some triangle. You might represent this as a set of polygons, perhaps even as a set of convex polygons.
For every triangle, if it is completely inside one region, you have to count it. So add the voume it spans with the origin (o.e. its determinant) to the accumulated volume. If it is completely outside all regions, don't count its volume but add its projection as a new polygon. If it is partially inside a region and partially outside, only compute and add the volume for the part which is outside, and add the projection to the region it did intersect. If it intersects more than one region, only count the part which is outside all of them, and consider joining the regions into a bigger one to simplify your data structures.
