I've just pushed some new changes to libigl and gptoolbox to expose the variadic implementation of our robust mesh boolean operations. By variadic, I mean that the boolean function can take one mesh as input or two meshes as input (usual binary case) or three or four and so on. This means you can easily take the union or intersection of n objects, but it also means that complex operations can be reduced to single call. For example, identifying regions inside at least k out of n objects. I first became aware of this variadic concept when reading "QuickCSG: Arbitrary and Faster Boolean Combinations of N Solids".
In libigl, you can call igl::copyleft::cgal::mesh_boolean with an n-long list of mesh vertex arrays and an n-long list of mesh face arrays. Rather than the usual union, intersection, minus etc., you can also pass a C++ function handle. This function handle will filter the winding number vector at any point in space. For example, if you'd like to extract the region inside object 0 and not inside object 1 and either inside object 2 or not inside object 3. The filter would return w(0) > 0 && w(1)<=0 && (w(2)>0 || w(3)<= 0).
After a bit of pointer tweaking, I've also exposed this interface to the matlab wrapper in gptoolbox. You can pass a matlab function handle after the optional argument 'WindingNumberFilter'. For example, assuming V and F are cell arrays containing vertex position arrays and face indices respectively, the following command will extract the region inside at least 3 of the inputs:
[VV,FF,J] = mesh_boolean(V,F,'','WindingNumberFilter',@(w) sum(w>0)>=3);
You could also call it with
[VV,FF,J] = mesh_boolean(V{1},F{1},V{2},F{2},V{3},F{3},V{4},F{4},'','WindingNumberFilter',@(w) sum(w>0)>=3);
Here's an example of extracting min-k results from 4 spheres:
Since extraction from our mesh arrangement is cheap, each operation takes the same amount of time. And fortunately it seems that the mex function handle overhead is not so bad.
Update: Variadic operations are also useful for condensing entire binary CSG trees: (A union B) minus (C union D) could be a single extraction function on (A,B,C,D) return (w(0)>0 || w(1)>0) && !(w(2)>0 || w(3)>0)