Algorithm for all intersecting geometry parts

Question

I'm trying to understand the ArcGIS "union" method so I can replicate it with an algorithm using Python and Shapely.

I do not want the cumulative intersection of all features where all of them intersect at the same time (ABC):

(Image is from another SO post)

Instead, I want each sub intersection that intersects at least once, repeated how many times it occurs. So in the example image I want ABC returned three times, and AB, AC, and BC returned twice. While the ArcGis method also returns the unique parts A, B, and C that have no intersections, I do not need this.

This is a pretty standard GIS method also available in QGIS, so what is the accepted set theoretic way of achieving this?

Code

In the ArcGis docs for union, they write that they use a "cracks and clusters" method, which is likely a low level modification of the intersection/union clipping algorithm. Since I am not implementing the underlying clipping algorithm, I need a different approach.

Here is pseudocode that I will be trying to implement. It basically means for each feature i find all intersections with others, and then recursively find all intersections between the intersecting parts, adding only the "node" parts that have no intersections with other intersections (i.e. we have cut it up as much as possible):

def isecs(g, geoms):
    for og in geoms:
        if og != g:
            if og.crosses(g):
                yield g.intersection(og)

def process(isecs):
    parts = []
    for g in isecs:
        isecs = getisecs(g, isecs)
        if no isecs:
            # "node" reached
            parts += g
        else:
            # this is the recursive part
            parts += process(isecs)
    return parts

for f in features:
    top_isecs = isecs(f.geom, features.bbox_overlaps(f.bbox))
    parts = process(top_isecs)
    for g in parts:
        addfeat(g)

Correction

@Vince correctly pointed out that the ArcGis operation is called union, instead of intersection as I first wrote. This is not to be confused that we are looking for geometrical unions, ArcGis is simply referring to the fact that they are returning all (hence union) geometrical intersections: "Union calculates the geometric intersection of any number of feature class..." They also allow returning an actual geometrical union through a dissolve option.

Answer 1

I think the key is finding all the overlapping polygons for each polygon, instead of checking for each set of polygons whether they intersect. That way the number of necessary itersect operations can be defined for each polygon (depending on the number of overlapping polygons). I would rewrite (in pseudo-code) your first function as follows (of course you can also use generators instead of lists):

def isecs(featureList):
    overlapList = {}
    for feature in featureList:
        overlappingFeatures = feature.intersects(featureList)
        # depending on code, may be necessary to exclude 'feature' from 'overlappingFeatures'
        overlapList[feature.id()] = overlappingFeatures 
    return overlapList

The overlapList dictionary will contain lists of overlapping polygons for A, B, ... and therefore will yield the correct number of copies of the intersecting parts: i.e. AB and BA for A and B.