Assigning surfaces to zones based on the 3D region

Given a set of surfaces in three-dimensional space, I am attempting to assign each surface to a zone referring to the smallest 3D region the set encloses, or no zone if this is not applicable. I also want to determine if a surface is an interface between two zones. So, for example, if we had 11 surfaces representing two cubes stacked on top of each other, the surfaces in the top cube would be in the same zone and the surfaces in the bottom would be in a different zone (with the interface surface being in both zones).

As an example, I want to take in a set of surfaces such as this and turn it in to this. Each color here represents a zone, with gray being no zone associated (as in the flap at the bottom).

I have done some searching around trying to find if someone has already come up with an algorithm to do this, but I have not found anything (most seem to identify regions rather than link surfaces to the region they enclose). As such I am trying to come up with my own algorithm and am wondering if there are any other alternatives or if my method would work.

I am assuming all surfaces are connected.

My idea is the following:

Select a random surface whose sides each touch exactly one other surface, and add this to zone 1.
Add each connected surface to zone 1 provided each of its sides touch exactly one other surface.
For those connected surfaces that touch more than one surface on at least one of its sides, add it to the "maybe" list.
For each new surface in zone 1, repeat steps 2-3.
Once a surface has been added to the "maybe" list twice, add it to zone 1 and remove from the "maybe" list. Mark this surface as a zone interface.
Add the zone interface to zone 2.
Select one random surface from the "maybe" list and assign it to zone 2 and clear the "maybe" list.
Repeat steps 2-7 (updating the zone number of course) until there are no surfaces that are unassigned.

This seems to work for simple scenarios (e.g., two cubes stacked on top of one another), but I am not sure if there are any tricky conditions I need to watch out for, or if it falls apart once there are more than two zones that share a side.

Any improvement on my rough algorithm/alternate ideas for implementation would be appreciated. Thanks!

EDIT: Here are some more details in response to some comments. A zone by my definition is simply a group of surfaces that completely bound a 3D region with no gaps. So if I had two cubes, A and B, that do not touch, I would have two zones: one consisting of all the surfaces of cube A and the other of all the surfaces for cube B. If I had a cube that was missing one side, there would be no zone associated with those surfaces.

My end goal is to make an automated process for grouping surfaces in a modeling tool I am creating. The specifics are classified, but essentially I am dealing with models where certain properties are common only between surfaces in the same "zone" as described above. I want to make an automated process that creates these zones so that the user can apply these properties to all surfaces in the zone at once instead of doing it manually.

Essentially the problem boils down to finding the smallest 3D regions that are completely enclosed by an arbitrary set of surfaces, and keeping track of which surfaces belong to which regions. I hope this makes my question more clear.

What you are interested in, then, would be discovering closed surface (volume) mesh topology from a set of input polygons; in other words - polytopes. This is common to pretty much every 3d modeling package. I would guess that Blender has code that does that. There are different ways of doing this, commonly however, some version of half-edge graph is used. See wiki link here: Doubly Linked half edge graph. The idea is to walk your input polies, and build these graphs. Once done, you can easily query each graph to see if there are holes (edges missing, etc).

I attached a picture explaining how to use a half-edge structure to get what you want: Say you are given a soup of five rectangles (they make up a cube with out a top. U process your first rectangle say ABCD, this creates your first graph, say G1. Now you process second polygon, say FEHG, none of these vertices you have seen yet, so you create second graph, G2. Now say you process polygon CDGH. You have seen these vertices before, so instead of creating a new graph, you merge(connect) existing graphs that share these nodes. Proceed until you process all polygons. You get graph in picture.

Now, to query the graph to get your information. Once you walk the graph, you will see that there are exactly four vertices (nodes) that are missing edges. Those verts correspond to the missing top of the box (the edges are red in the illustration). Hence you know that this graph is not a closed manifold. If you had another box, that did not share nodes with this one, you would have another graph. So each graph, once you done processing your polygons, is a "zone" for you.

Note, if you have two say intersecting shapes, you can track those too using these graphs, but its much more complicated. Basically when processing a new polygon, you would not only have to see if any of its verts belong to already processed graphs, but also see if this polygon intersects any of the previously processed polygons, if so, split this polygon and add all this to the intersected graph.