How is GEOS UnaryUnion function implemented?

08 Jan 2013

At work I encountered the use of gaiaUnaryUnion function to merge (union) polygons, and I wondered how it is actually implemented, my guess was it must have built some kind of tree structure out of all the input polygons.

Let’s dive in to some code in GEOS:

/* In CascadedPolygonUnion.cpp */

geom::Geometry* CascadedPolygonUnion::Union()
    if (inputPolys->empty())
        return NULL;

    geomFactory = inputPolys->front()->getFactory();

     * A spatial index to organize the collection
     * into groups of close geometries.
     * This makes unioning more efficient, since vertices 
     * are more likely to be eliminated on each round.
    index::strtree::STRtree index(STRTREE_NODE_CAPACITY);

    typedef std::vector<geom::Polygon*>::iterator iterator_type;
    iterator_type end = inputPolys->end();
    for (iterator_type i = inputPolys->begin(); i != end; ++i) {
        geom::Geometry* g = dynamic_cast<geom::Geometry*>(*i);
        index.insert(g->getEnvelopeInternal(), g);

        itemTree (index.itemsTree());

    return unionTree(itemTree.get());

So indeed, it uses a R-Tree to collect spatially neighboring polygons together and forms a tree! After the tree is built, it is flattened out into a list (itemTree; needs to be verified), which is then thrown to unionTree() to do the real union/merge work.

Now let’s have a look of how unionTree is realized:

/* Still in CascadedPolygonUnion.cpp */

geom::Geometry* CascadedPolygonUnion::unionTree(
    index::strtree::ItemsList* geomTree)
     * Recursively unions all subtrees in the list into 
     * single geometries. The result is a list of Geometry's only
    std::auto_ptr<GeometryListHolder> geoms(reduceToGeometries(geomTree));
    return binaryUnion(geoms.get());

/* reduceToGeometries() calls unionTree() from within; 
 * this is actually inorder-traversal over the R-Tree or 
 * walking the tree in a depth-first fashion. It has to 
 * be this way, since we want to begin merging from the 
 * smallest geometries.
CascadedPolygonUnion::reduceToGeometries (index::strtree::ItemsList* geomTree)
        geoms (new GeometryListHolder());

    typedef index::strtree::ItemsList::iterator iterator_type;
    iterator_type end = geomTree->end();
    for (iterator_type i = geomTree->begin(); i != end; ++i) {
        if ((*i).get_type() == index::strtree::ItemsListItem::item_is_list) {
            std::auto_ptr<geom::Geometry> geom (unionTree((*i).get_itemslist()));
        else if ((*i).get_type() == index::strtree::ItemsListItem::item_is_geometry) {
        else {
            assert(!"should never be reached");

    return geoms.release();

/* Divide-and-conquer union of all the geometries */
geom::Geometry* CascadedPolygonUnion::binaryUnion(GeometryListHolder* geoms, 
    std::size_t start, std::size_t end)
    if (end - start <= 1) {
        return unionSafe(geoms->getGeometry(start), NULL);
    else if (end - start == 2) {
        return unionSafe(geoms->getGeometry(start), geoms->getGeometry(start + 1));
    else {
        // recurse on both halves of the list
        std::size_t mid = (end + start) / 2;
        std::auto_ptr<geom::Geometry> g0 (binaryUnion(geoms, start, mid));
        std::auto_ptr<geom::Geometry> g1 (binaryUnion(geoms, mid, end));
        return unionSafe(g0.get(), g1.get());

In conclusion, the implementation walks the tree in a depth-first fashion and merges geometries along the way with a divide-and-conquer operation - binaryUnion(). So this is it! Conceptually it’s quite simple, isn’t it?!

to flatten a tree structure efficiently. In a follow-up to this post, you might want to figure out how actually a R-Tree is implemented and how