cursor给的源码

函数声明

#include <geos/geom.h>
#include <geos.h>

// #include <geos_c.h>
#include<iostream>
#include <vector>
#include <memory>

using namespace geos::geom;

// 定义二维点类型（使用 double 坐标）
struct Pointd {
    double x, y;
    Pointd(double x = 0.0, double y = 0.0) : x(x), y(y) {}
};

// 函数声明：计算两个多边形集合的并集
int polygon_union(
    const std::vector<std::vector<Pointd>>& polygons1,
    const std::vector<std::vector<Pointd>>& polygons2,
    std::vector<std::vector<Pointd>>& polygons_out,
    bool& is_union);

函数实现

#include "Union.h"
#include <geos_c.h>
#include <vector>
#include <memory>
#include <iostream> // 用于调试

// 前向声明
static GEOSGeometry* createPolygonWithHoles(
    GEOSContextHandle_t handle,
    const std::vector<Pointd>& shell,
    const std::vector<std::vector<Pointd>>& holes);

static void extractPolygon(
    GEOSContextHandle_t handle, 
    const GEOSGeometry* poly, 
    std::vector<std::vector<Pointd>>& polygons_out);

int polygon_union(
    const std::vector<std::vector<Pointd>>& polygons1,
    const std::vector<std::vector<Pointd>>& polygons2,
    std::vector<std::vector<Pointd>>& polygons_out,
    bool& is_union) {

    // 调试信息
    std::cout << "Input polygons1 count: " << polygons1.size() << std::endl;
    std::cout << "Input polygons2 count: " << polygons2.size() << std::endl;

    // 初始化 GEOS
    GEOSContextHandle_t handle = GEOS_init_r();
    if (!handle) return -1;

    try {
        // 创建第一个多边形（带内环）
        GEOSGeometry* poly1 = nullptr;
        if (!polygons1.empty()) {
            // 第一个vector是外环，其余的都是内环
            const auto& shell = polygons1[0];
            std::vector<std::vector<Pointd>> holes;
            if (polygons1.size() > 1) {
                holes.assign(polygons1.begin() + 1, polygons1.end());
            }
            poly1 = createPolygonWithHoles(handle, shell, holes);
        }

        // 创建第二个多边形（带内环）
        GEOSGeometry* poly2 = nullptr;
        if (!polygons2.empty()) {
            const auto& shell = polygons2[0];
            std::vector<std::vector<Pointd>> holes;
            if (polygons2.size() > 1) {
                holes.assign(polygons2.begin() + 1, polygons2.end());
            }
            poly2 = createPolygonWithHoles(handle, shell, holes);
        }

        // 检查输入是否有效
        if (!poly1 && !poly2) {
            is_union = false;
            GEOS_finish_r(handle);
            return -1;
        }

        // 计算并集
        GEOSGeometry* result = nullptr;
        if (!poly1) {
            result = GEOSGeom_clone_r(handle, poly2);
            is_union = false;
        } 
        else if (!poly2) {
            result = GEOSGeom_clone_r(handle, poly1);
            is_union = false;
        }
        else {
            result = GEOSUnion_r(handle, poly1, poly2);
            // 检查是否有相交
            GEOSGeometry* intersection = GEOSIntersection_r(handle, poly1, poly2);
            is_union = intersection && !GEOSisEmpty_r(handle, intersection);
            if (intersection) GEOSGeom_destroy_r(handle, intersection);
        }

        // 提取结果
        polygons_out.clear();
        if (result) {
            int type = GEOSGeomTypeId_r(handle, result);
            
            if (type == 3) { // GEOS_POLYGON
                // 提取单个多边形
                extractPolygon(handle, result, polygons_out);
            }
            else if (type == 6) { // GEOS_MULTIPOLYGON
                // 处理多个多边形的情况
                for (int i = 0; i < GEOSGetNumGeometries_r(handle, result); i++) {
                    const GEOSGeometry* geom = GEOSGetGeometryN_r(handle, result, i);
                    if (geom) {
                        size_t current_size = polygons_out.size();
                        extractPolygon(handle, geom, polygons_out);
                        
                        // 如果这不是第一个多边形，在结果中添加分隔标记
                        if (i > 0 && polygons_out.size() > current_size) {
                            // 可以在这里添加多边形之间的分隔标记
                            // 比如添加一个空的vector作为分隔符
                            polygons_out.push_back(std::vector<Pointd>());
                        }
                    }
                }
            }
        }

        // 清理资源
        if (poly1) GEOSGeom_destroy_r(handle, poly1);
        if (poly2) GEOSGeom_destroy_r(handle, poly2);
        if (result) GEOSGeom_destroy_r(handle, result);

        GEOS_finish_r(handle);
        return 0;

    } catch (const std::exception& e) {
        std::cout << "Exception: " << e.what() << std::endl;
        if (handle) GEOS_finish_r(handle);
        return -1;
    } catch (...) {
        std::cout << "Unknown exception occurred" << std::endl;
        if (handle) GEOS_finish_r(handle);
        return -1;
    }
}

// 辅助函数：创建带内环的多边形
static GEOSGeometry* createPolygonWithHoles(
    GEOSContextHandle_t handle,
    const std::vector<Pointd>& shell,
    const std::vector<std::vector<Pointd>>& holes) {
    
    if (shell.size() < 3) return nullptr;

    // 创建外环
    GEOSCoordSequence* shell_seq = GEOSCoordSeq_create_r(handle, shell.size() + 1, 2);
    for (size_t i = 0; i < shell.size(); i++) {
        GEOSCoordSeq_setX_r(handle, shell_seq, i, shell[i].x);
        GEOSCoordSeq_setY_r(handle, shell_seq, i, shell[i].y);
    }
    // 闭合多边形
    GEOSCoordSeq_setX_r(handle, shell_seq, shell.size(), shell[0].x);
    GEOSCoordSeq_setY_r(handle, shell_seq, shell.size(), shell[0].y);

    GEOSGeometry* shell_ring = GEOSGeom_createLinearRing_r(handle, shell_seq);
    if (!shell_ring) return nullptr;

    // 创建内环
    std::vector<GEOSGeometry*> hole_rings;
    for (const auto& hole : holes) {
        if (hole.size() < 3) continue;

        GEOSCoordSequence* hole_seq = GEOSCoordSeq_create_r(handle, hole.size() + 1, 2);
        for (size_t i = 0; i < hole.size(); i++) {
            GEOSCoordSeq_setX_r(handle, hole_seq, i, hole[i].x);
            GEOSCoordSeq_setY_r(handle, hole_seq, i, hole[i].y);
        }
        // 闭合内环
        GEOSCoordSeq_setX_r(handle, hole_seq, hole.size(), hole[0].x);
        GEOSCoordSeq_setY_r(handle, hole_seq, hole.size(), hole[0].y);

        GEOSGeometry* hole_ring = GEOSGeom_createLinearRing_r(handle, hole_seq);
        if (hole_ring) {
            hole_rings.push_back(hole_ring);
        }
    }

    // 创建完整的多边形
    GEOSGeometry* polygon = GEOSGeom_createPolygon_r(
        handle,
        shell_ring,
        hole_rings.data(),
        hole_rings.size()
    );

    return polygon;
}

// 辅助函数：提取多边形的坐标
static void extractPolygon(
    GEOSContextHandle_t handle, 
    const GEOSGeometry* poly, 
    std::vector<std::vector<Pointd>>& polygons_out) {
    
    // 提取外环
    const GEOSGeometry* shell = GEOSGetExteriorRing_r(handle, poly);
    if (shell) {
        std::vector<Pointd> shell_points;
        const GEOSCoordSequence* coords = GEOSGeom_getCoordSeq_r(handle, shell);
        unsigned int size;
        GEOSCoordSeq_getSize_r(handle, coords, &size);
        
        for (unsigned int i = 0; i < size - 1; i++) { // size-1 to avoid duplicate point
            double x, y;
            GEOSCoordSeq_getX_r(handle, coords, i, &x);
            GEOSCoordSeq_getY_r(handle, coords, i, &y);
            shell_points.emplace_back(x, y);
        }
        if (!shell_points.empty()) {
            polygons_out.push_back(shell_points);
        }

        // 提取内环
        int num_holes = GEOSGetNumInteriorRings_r(handle, poly);
        for (int i = 0; i < num_holes; i++) {
            const GEOSGeometry* hole = GEOSGetInteriorRingN_r(handle, poly, i);
            if (hole) {
                std::vector<Pointd> hole_points;
                coords = GEOSGeom_getCoordSeq_r(handle, hole);
                GEOSCoordSeq_getSize_r(handle, coords, &size);
                
                for (unsigned int j = 0; j < size - 1; j++) { // size-1 to avoid duplicate point
                    double x, y;
                    GEOSCoordSeq_getX_r(handle, coords, j, &x);
                    GEOSCoordSeq_getY_r(handle, coords, j, &y);
                    hole_points.emplace_back(x, y);
                }
                if (!hole_points.empty()) {
                    polygons_out.push_back(hole_points);
                }
            }
        }
    }
}

使用示例

#include <iostream>
#include <vector>
#include "Union.h"

int main() {
    // 第一个多边形：正方形 (0,0) -> (1,0) -> (1,1) -> (0,1)
    std::vector<std::vector<Pointd>> polygons1 = {
        {   // 外环
            Pointd(0, 0),
            Pointd(1, 0),
            Pointd(1, 1),
            Pointd(0, 1)
        }
        // 可以添加内环，如: {Pointd(...), Pointd(...), ...}
    };

    // 第二个多边形：正方形 (0.5,0.5) -> (1.5,0.5) -> (1.5,1.5) -> (0.5,1.5)
    std::vector<std::vector<Pointd>> polygons2 = {
        {   // 外环
            Pointd(0.5, 0.5),
            Pointd(1.5, 0.5),
            Pointd(1.5, 1.5),
            Pointd(0.5, 1.5)
        }
    };
    
    std::vector<std::vector<Pointd>> result;
    bool is_union;
    
    int status = polygon_union(polygons1, polygons2, result, is_union);
    
    if (status == 0) {
        std::cout << "Union successful!" << std::endl;
        std::cout << "Is union: " << (is_union ? "true" : "false") << std::endl;
        std::cout << "Result has " << result.size() << " rings:" << std::endl;
        for (size_t i = 0; i < result.size(); ++i) {
            std::cout << (i == 0 ? "Exterior ring: " : "Interior ring: ");
            for (const auto& p : result[i]) {
                std::cout << "(" << p.x << "," << p.y << ") ";
            }
            std::cout << std::endl;
        }
    } else {
        std::cerr << "Union failed!" << std::endl;
    }
    
    return 0;
}

cmake文件

cmake_minimum_required(VERSION 3.10)
project(MyProject)

# 使用 pkg-config 查找 GEOS
find_package(PkgConfig REQUIRED)
pkg_check_modules(GEOS REQUIRED geos)

# GEOS 3.13 源码默认安装路径
set(GEOS_INCLUDE_DIR "/usr/local/include/geos")
set(GEOS_LIBRARY_DIR "/usr/local/lib")

# 添加可执行文件
add_executable(my_app main.cpp Union.cpp)
target_include_directories(my_app PRIVATE ${GEOS_INCLUDE_DIR})
target_link_directories(my_app PRIVATE ${GEOS_LIBRARY_DIR})
target_link_libraries(my_app ${GEOS_LIBRARIES})

安装和准备

参考官方Download and Build
先下载源码包geos-3.13.1.tar.bz2
然后执行解压和编译，测试，安装

# Unpack and setup build directory
tar xvfj geos-3.13.1.tar.bz2
cd geos-3.13.1
mkdir _build
cd _build
# Set up the build
cmake \
    -DCMAKE_BUILD_TYPE=Release \
    -DCMAKE_INSTALL_PREFIX=/usr/local \
    ..
# Run the build, test, install
make
ctest
make install