znlgis 博客

GIS开发与技术分享

第五章 矩形裁剪与闵可夫斯基操作(C#版)

5.1 引言

除了布尔运算和多边形偏移之外,Clipper2还提供了两个专门的几何操作:矩形裁剪(Rectangle Clipping)和闵可夫斯基运算(Minkowski Operations)。矩形裁剪是一种针对轴对齐矩形优化的高效裁剪算法,在地图瓦片生成、视窗裁剪等场景下非常有用。闵可夫斯基运算则在碰撞检测、机器人路径规划等领域有着重要的应用。本章将详细介绍Clipper2 C#版本中这两种操作的原理和使用方法。

5.2 矩形裁剪

5.2.1 什么是矩形裁剪

矩形裁剪是使用轴对齐矩形(Axis-Aligned Rectangle)作为裁剪区域对多边形进行裁剪的操作。由于矩形具有特殊的几何性质(四条边分别与坐标轴平行),可以设计比通用布尔运算更高效的算法。

     原始多边形              裁剪矩形              裁剪结果
    
        ╱╲                  ┌─────┐               ┌─────┐
       ╱  ╲                 │     │               │     │
      ╱    ╲         ∩      │     │        =      │     │
     ╱      ╲               │     │               │ ╱╲  │
    ╱────────╲              └─────┘               ╱─┘  └─╲

5.2.2 矩形裁剪的优势

相比于使用通用布尔运算进行裁剪,专用的矩形裁剪算法具有以下优势:

更高的性能

矩形裁剪算法的时间复杂度更低,通常比通用布尔运算快2-5倍,对于简单的多边形可能更快。

更简单的实现

由于矩形边与坐标轴对齐,边与边的交点计算更加简单直接。

内存效率

算法不需要构建复杂的事件队列和活动边列表,内存占用更低。

5.2.3 使用RectClip函数

Clipper2提供了RectClip函数用于裁剪闭合多边形:

#include "clipper2/clipper.h"
using namespace Clipper2Lib;

int main() {
    // 创建要裁剪的多边形
    Paths64 subject;
    subject.push_back(MakePath({0, 0, 200, 0, 200, 200, 0, 200}));
    subject.push_back(MakePath({50, 50, 150, 50, 150, 150, 50, 150}));  // 孔洞
    
    // 定义裁剪矩形
    Rect64 clipRect(25, 25, 175, 175);
    
    // 执行矩形裁剪
    Paths64 result = RectClip(clipRect, subject);
    
    // result 包含裁剪后的多边形
    return 0;
}

5.2.4 使用RectClipLines函数

对于开放的折线(而非闭合多边形),使用RectClipLines函数:

// 创建开放折线
Paths64 lines;
lines.push_back(MakePath({0, 100, 200, 100}));  // 水平线
lines.push_back(MakePath({100, 0, 100, 200}));  // 垂直线
lines.push_back(MakePath({0, 0, 200, 200}));    // 对角线

// 定义裁剪矩形
Rect64 clipRect(50, 50, 150, 150);

// 裁剪折线
Paths64 result = RectClipLines(clipRect, lines);
// result 包含被裁剪到矩形内的线段

5.2.5 RectClip64类

对于需要多次使用同一矩形进行裁剪的场景,可以使用RectClip64类以获得更好的性能:

// 创建裁剪对象
Rect64 clipRect(0, 0, 100, 100);
RectClip64 rectClipper(clipRect);

// 裁剪多个多边形
for (const auto& polygon : polygons) {
    Paths64 result = rectClipper.Execute(Paths64{polygon});
    ProcessResult(result);
}

5.2.6 RectClipLines64类

类似地,对于折线裁剪:

Rect64 clipRect(0, 0, 100, 100);
RectClipLines64 lineClipper(clipRect);

for (const auto& line : lines) {
    Paths64 result = lineClipper.Execute(Paths64{line});
    ProcessResult(result);
}

5.2.7 浮点数版本

使用浮点数坐标的版本:

// 浮点数矩形裁剪
PathsD subject;
subject.push_back(MakePathD({0.0, 0.0, 10.0, 0.0, 10.0, 10.0, 0.0, 10.0}));

RectD clipRect(2.5, 2.5, 7.5, 7.5);

PathsD result = RectClip(clipRect, subject);

5.2.8 矩形裁剪的边界情况

边界上的顶点

当多边形顶点恰好在裁剪矩形的边界上时,算法会正确处理:

// 顶点在边界上的多边形
Paths64 subject;
subject.push_back(MakePath({50, 0, 100, 50, 50, 100, 0, 50}));  // 菱形

Rect64 clipRect(0, 0, 100, 100);  // 边界与菱形顶点重合
Paths64 result = RectClip(clipRect, subject);

完全在矩形内的多边形

如果多边形完全在裁剪矩形内部,返回原始多边形:

Paths64 subject;
subject.push_back(MakePath({25, 25, 75, 25, 75, 75, 25, 75}));

Rect64 clipRect(0, 0, 100, 100);
Paths64 result = RectClip(clipRect, subject);
// result 与 subject 相同

完全在矩形外的多边形

如果多边形完全在裁剪矩形外部,返回空结果:

Paths64 subject;
subject.push_back(MakePath({200, 200, 300, 200, 300, 300, 200, 300}));

Rect64 clipRect(0, 0, 100, 100);
Paths64 result = RectClip(clipRect, subject);
// result 为空

5.2.9 矩形裁剪的应用场景

地图瓦片生成

// 生成256x256的地图瓦片
void GenerateTile(int tileX, int tileY, int zoom, const Paths64& mapData) {
    // 计算瓦片的地理范围
    int tileSize = 256;
    int worldSize = tileSize << zoom;
    
    Rect64 tileBounds(
        tileX * tileSize,
        tileY * tileSize,
        (tileX + 1) * tileSize,
        (tileY + 1) * tileSize
    );
    
    // 裁剪地图数据到瓦片范围
    Paths64 tileData = RectClip(tileBounds, mapData);
    
    // 渲染瓦片
    RenderTile(tileData);
}

视窗裁剪

// 裁剪到可见视窗区域
Paths64 ClipToViewport(const Paths64& geometry, 
                        int viewportWidth, int viewportHeight,
                        int scrollX, int scrollY) {
    Rect64 viewport(
        scrollX,
        scrollY,
        scrollX + viewportWidth,
        scrollY + viewportHeight
    );
    
    return RectClip(viewport, geometry);
}

空间分区

// 将几何体分配到四叉树节点
void QuadTreeInsert(QuadTreeNode& node, const Paths64& geometry) {
    Rect64 nodeBounds = node.getBounds();
    Paths64 clipped = RectClip(nodeBounds, geometry);
    
    if (!clipped.empty()) {
        if (node.isLeaf() || Area(clipped) < threshold) {
            node.addGeometry(clipped);
        } else {
            // 分配到子节点
            for (auto& child : node.children()) {
                QuadTreeInsert(child, clipped);
            }
        }
    }
}

5.3 闵可夫斯基运算

5.3.1 什么是闵可夫斯基运算

闵可夫斯基运算是以德国数学家赫尔曼·闵可夫斯基(Hermann Minkowski)命名的几何操作。主要包括两种运算:

闵可夫斯基和(Minkowski Sum)

两个点集A和B的闵可夫斯基和定义为: A ⊕ B = {a + b | a ∈ A, b ∈ B}

直观理解:将形状B的中心放在形状A的每个边界点上,所有这些B形状的并集就是闵可夫斯基和。

     多边形A              多边形B            闵可夫斯基和
    
    ┌───────┐              ○               ╭───────────╮
    │       │              ↓               │   ╭─────╮ │
    │       │      ⊕                  =    │   │     │ │
    │       │                              │   ╰─────╯ │
    └───────┘                              ╰───────────╯
                                          (圆角矩形)

闵可夫斯基差(Minkowski Difference)

A和B的闵可夫斯基差定义为: A ⊖ B = A ⊕ (-B) = {a - b | a ∈ A, b ∈ B}

其中-B是B关于原点的反射。

5.3.2 闵可夫斯基和的几何意义

闵可夫斯基和有几个重要的几何解释:

形态膨胀

当B是一个以原点为中心的圆时,A ⊕ B 等于A的圆形偏移。

可达区域

如果A表示障碍物,B表示移动物体,那么A ⊕ (-B)表示移动物体中心不能到达的区域。

碰撞检测

两个多边形A和B相交,当且仅当(A ⊖ B)包含原点。

5.3.3 使用MinkowskiSum函数

Clipper2提供了闵可夫斯基和函数:

#include "clipper2/clipper.h"
using namespace Clipper2Lib;

int main() {
    // 创建第一个多边形(矩形)
    Path64 pattern = MakePath({-50, -50, 50, -50, 50, 50, -50, 50});
    
    // 创建第二个多边形(三角形,以原点为中心)
    Path64 path = MakePath({0, -30, 26, 15, -26, 15});
    
    // 计算闵可夫斯基和
    Paths64 result = MinkowskiSum(pattern, path, false);
    // 参数:pattern(图案), path(路径), isClosed(路径是否闭合)
    
    return 0;
}

5.3.4 使用MinkowskiDiff函数

// 计算闵可夫斯基差
Path64 polygon1 = MakePath({0, 0, 100, 0, 100, 100, 0, 100});
Path64 polygon2 = MakePath({-20, -20, 20, -20, 20, 20, -20, 20});

Paths64 diff = MinkowskiDiff(polygon1, polygon2, false);

5.3.5 闵可夫斯基和的参数

Paths64 MinkowskiSum(
    const Path64& pattern,    // 图案多边形
    const Path64& path,       // 路径
    bool pathIsClosed,        // 路径是否闭合
    FillRule fillRule = FillRule::NonZero  // 填充规则
);

pattern(图案):被复制到path每个顶点的形状

path(路径):定义图案放置位置的路径

pathIsClosed

5.3.6 开放路径与闭合路径

// 闭合路径的闵可夫斯基和
Path64 circle = MakeCircle(Point64(0, 0), 20);  // 圆形图案
Path64 square = MakePath({0, 0, 100, 0, 100, 100, 0, 100});  // 正方形路径

// 闭合路径:产生圆角矩形
Paths64 closedResult = MinkowskiSum(circle, square, true);

// 开放路径:产生胶囊形
Path64 line = MakePath({0, 0, 100, 0});  // 线段路径
Paths64 openResult = MinkowskiSum(circle, line, false);

5.3.7 实际应用示例

碰撞检测

使用闵可夫斯基差进行碰撞检测:

// 两个多边形
Path64 polygonA = MakePath({0, 0, 50, 0, 50, 50, 0, 50});
Path64 polygonB = MakePath({30, 30, 80, 30, 80, 80, 30, 80});

// 计算闵可夫斯基差
Paths64 minkDiff = MinkowskiDiff(polygonA, polygonB, true);

// 检查原点是否在闵可夫斯基差内
Point64 origin(0, 0);
bool colliding = false;
for (const auto& path : minkDiff) {
    if (PointInPolygon(origin, path) != PointInPolygonResult::IsOutside) {
        colliding = true;
        break;
    }
}

if (colliding) {
    std::cout << "多边形相交!" << std::endl;
}

机器人路径规划

使用闵可夫斯基和计算配置空间障碍物:

// 机器人形状(以中心为原点)
Path64 robotShape = MakePath({-10, -10, 10, -10, 10, 10, -10, 10});

// 障碍物列表
std::vector<Path64> obstacles = {
    MakePath({100, 100, 200, 100, 200, 200, 100, 200}),
    MakePath({300, 50, 400, 50, 400, 150, 300, 150})
};

// 计算配置空间障碍物
Paths64 configSpaceObstacles;
for (const auto& obstacle : obstacles) {
    // 使用机器人形状的反射
    Path64 reflectedRobot;
    for (const auto& pt : robotShape) {
        reflectedRobot.push_back(Point64(-pt.x, -pt.y));
    }
    
    // 闵可夫斯基和
    Paths64 expanded = MinkowskiSum(reflectedRobot, obstacle, true);
    
    // 合并到配置空间障碍物
    configSpaceObstacles.insert(configSpaceObstacles.end(),
                                 expanded.begin(), expanded.end());
}

// 合并重叠的障碍物
Paths64 mergedObstacles = Union(configSpaceObstacles, FillRule::NonZero);

// 现在可以在点空间中进行路径规划,只要路径不穿过这些障碍物

形态学膨胀

使用闵可夫斯基和模拟形态学膨胀:

// 原始形状
Path64 shape = MakePath({...});

// 结构元素(圆形)
Path64 structuringElement = MakeCircle(Point64(0, 0), radius);

// 形态学膨胀
Paths64 dilated = MinkowskiSum(structuringElement, shape, true);

计算可达边界

// 移动物体的形状
Path64 movingObject = MakePath({-5, -5, 5, -5, 5, 5, -5, 5});

// 轨迹路径
Path64 trajectory = MakePath({0, 0, 100, 50, 200, 0, 300, 50});

// 计算移动物体沿轨迹运动时扫过的区域
Paths64 sweptArea = MinkowskiSum(movingObject, trajectory, false);

5.3.8 闵可夫斯基运算的性能考虑

复杂度分析

闵可夫斯基和的时间复杂度大约为O(m * n),其中m和n分别是两个多边形的顶点数。对于复杂的多边形,这可能会很慢。

优化策略

// 1. 简化输入多边形
Path64 simplifiedPattern = SimplifyPath(pattern, tolerance);
Path64 simplifiedPath = SimplifyPath(path, tolerance);
Paths64 result = MinkowskiSum(simplifiedPattern, simplifiedPath, true);

// 2. 对于凸多边形使用专门的算法
if (IsConvex(pattern) && IsConvex(path)) {
    // 凸多边形的闵可夫斯基和可以更高效地计算
    Paths64 result = ConvexMinkowskiSum(pattern, path);
}

// 3. 分解为凸多边形
Paths64 convexPattern = ConvexPartition(pattern);
Paths64 convexPath = ConvexPartition(path);
Paths64 result;
for (const auto& cp : convexPattern) {
    for (const auto& cpp : convexPath) {
        Paths64 partial = MinkowskiSum(cp, cpp, true);
        result.insert(result.end(), partial.begin(), partial.end());
    }
}
result = Union(result, FillRule::NonZero);

5.3.9 闵可夫斯基和的数学性质

交换律

A ⊕ B = B ⊕ A

结合律

(A ⊕ B) ⊕ C = A ⊕ (B ⊕ C)

分配律(对于并集)

A ⊕ (B ∪ C) = (A ⊕ B) ∪ (A ⊕ C)

与原点的关系

如果O是只包含原点的点集,则 A ⊕ O = A

这些性质可以用于优化计算:

// 利用分配律分解计算
Paths64 A = ...;
Paths64 B = ...;  // B由多个分离的多边形组成

// 方法1:直接计算(如果B是一个Paths)
// Paths64 result = MinkowskiSum(A, B);

// 方法2:分解计算(可能更快)
Paths64 result;
for (const auto& bi : B) {
    Paths64 partial = MinkowskiSum(A[0], bi, true);
    result.insert(result.end(), partial.begin(), partial.end());
}
result = Union(result, FillRule::NonZero);

5.4 综合应用案例

5.4.1 地图瓦片系统

class TileSystem {
private:
    int tileSize_;
    int maxZoom_;
    
public:
    TileSystem(int tileSize = 256, int maxZoom = 20)
        : tileSize_(tileSize), maxZoom_(maxZoom) {}
    
    Paths64 ClipToTile(const Paths64& geometry, 
                        int x, int y, int zoom) const {
        // 计算瓦片边界
        int64_t scale = 1LL << zoom;
        int64_t tileX = x * tileSize_;
        int64_t tileY = y * tileSize_;
        
        Rect64 tileBounds(
            tileX, tileY,
            tileX + tileSize_, tileY + tileSize_
        );
        
        // 使用矩形裁剪
        return RectClip(tileBounds, geometry);
    }
    
    std::map<std::pair<int, int>, Paths64> 
    GenerateTiles(const Paths64& geometry, int zoom) const {
        std::map<std::pair<int, int>, Paths64> tiles;
        
        // 获取几何体的边界框
        Rect64 bounds = Bounds(geometry);
        
        // 计算涉及的瓦片范围
        int minTileX = bounds.left / tileSize_;
        int minTileY = bounds.top / tileSize_;
        int maxTileX = bounds.right / tileSize_;
        int maxTileY = bounds.bottom / tileSize_;
        
        // 为每个瓦片裁剪几何体
        for (int ty = minTileY; ty <= maxTileY; ty++) {
            for (int tx = minTileX; tx <= maxTileX; tx++) {
                Paths64 tileGeometry = ClipToTile(geometry, tx, ty, zoom);
                if (!tileGeometry.empty()) {
                    tiles[{tx, ty}] = std::move(tileGeometry);
                }
            }
        }
        
        return tiles;
    }
};

5.4.2 碰撞检测系统

class CollisionDetector {
public:
    // 检测两个多边形是否碰撞
    bool CheckCollision(const Path64& polygonA, 
                        const Path64& polygonB) const {
        // 计算闵可夫斯基差
        Paths64 minkDiff = MinkowskiDiff(polygonA, polygonB, true);
        
        // 检查原点是否在差集内
        Point64 origin(0, 0);
        for (const auto& path : minkDiff) {
            auto result = PointInPolygon(origin, path);
            if (result != PointInPolygonResult::IsOutside) {
                return true;
            }
        }
        return false;
    }
    
    // 计算穿透深度和方向
    std::pair<double, PointD> GetPenetration(
            const Path64& polygonA, 
            const Path64& polygonB) const {
        
        Paths64 minkDiff = MinkowskiDiff(polygonA, polygonB, true);
        
        Point64 origin(0, 0);
        double minDist = std::numeric_limits<double>::max();
        PointD minDirection(0, 0);
        
        for (const auto& path : minkDiff) {
            auto result = PointInPolygon(origin, path);
            if (result == PointInPolygonResult::IsOutside) {
                continue;
            }
            
            // 找到原点到边界的最短距离
            for (size_t i = 0; i < path.size(); i++) {
                Point64 p1 = path[i];
                Point64 p2 = path[(i + 1) % path.size()];
                
                // 计算点到线段的距离
                double dist = DistanceToSegment(origin, p1, p2);
                if (dist < minDist) {
                    minDist = dist;
                    // 计算法线方向
                    PointD edge(p2.x - p1.x, p2.y - p1.y);
                    double len = std::sqrt(edge.x * edge.x + edge.y * edge.y);
                    minDirection = PointD(-edge.y / len, edge.x / len);
                }
            }
        }
        
        return {minDist, minDirection};
    }
    
private:
    double DistanceToSegment(const Point64& pt, 
                             const Point64& a, 
                             const Point64& b) const {
        double dx = b.x - a.x;
        double dy = b.y - a.y;
        double t = std::max(0.0, std::min(1.0,
            ((pt.x - a.x) * dx + (pt.y - a.y) * dy) / (dx * dx + dy * dy)));
        double projX = a.x + t * dx;
        double projY = a.y + t * dy;
        return std::sqrt((pt.x - projX) * (pt.x - projX) + 
                         (pt.y - projY) * (pt.y - projY));
    }
};

5.4.3 扫掠体积计算

class SweptVolume {
public:
    // 计算物体沿路径移动时扫过的区域
    Paths64 ComputeSweptArea(const Path64& objectShape,
                              const Path64& motionPath,
                              bool closedPath = false) const {
        // 使用闵可夫斯基和
        return MinkowskiSum(objectShape, motionPath, closedPath);
    }
    
    // 计算旋转物体的扫掠区域(近似)
    Paths64 ComputeRotationalSweep(const Path64& objectShape,
                                    const Point64& pivot,
                                    double startAngle,
                                    double endAngle,
                                    int numSteps = 36) const {
        Paths64 result;
        
        double angleStep = (endAngle - startAngle) / numSteps;
        
        for (int i = 0; i <= numSteps; i++) {
            double angle = startAngle + i * angleStep;
            Path64 rotated = RotatePath(objectShape, pivot, angle);
            result.push_back(rotated);
        }
        
        // 合并所有位置
        return Union(result, FillRule::NonZero);
    }
    
private:
    Path64 RotatePath(const Path64& path, 
                      const Point64& pivot, 
                      double angle) const {
        Path64 result;
        double cosA = std::cos(angle);
        double sinA = std::sin(angle);
        
        for (const auto& pt : path) {
            double dx = pt.x - pivot.x;
            double dy = pt.y - pivot.y;
            result.push_back(Point64(
                static_cast<int64_t>(pivot.x + dx * cosA - dy * sinA),
                static_cast<int64_t>(pivot.y + dx * sinA + dy * cosA)
            ));
        }
        
        return result;
    }
};

5.4.4 安全区域计算

class SafetyZoneCalculator {
public:
    // 计算机器人的安全导航区域
    Paths64 ComputeSafeZone(const Paths64& obstacles,
                             const Path64& robotShape,
                             double additionalMargin = 0) const {
        Paths64 expandedObstacles;
        
        // 计算机器人形状的反射
        Path64 reflectedRobot;
        for (const auto& pt : robotShape) {
            reflectedRobot.push_back(Point64(-pt.x, -pt.y));
        }
        
        // 如果需要额外边距,先膨胀反射形状
        Path64 marginedRobot = reflectedRobot;
        if (additionalMargin > 0) {
            Paths64 inflated = InflatePaths(Paths64{reflectedRobot}, 
                                            additionalMargin,
                                            JoinType::Round, 
                                            EndType::Polygon);
            if (!inflated.empty()) {
                marginedRobot = inflated[0];
            }
        }
        
        // 扩展每个障碍物
        for (const auto& obstacle : obstacles) {
            Paths64 expanded = MinkowskiSum(marginedRobot, obstacle, true);
            expandedObstacles.insert(expandedObstacles.end(),
                                      expanded.begin(), expanded.end());
        }
        
        // 合并重叠的障碍物
        return Union(expandedObstacles, FillRule::NonZero);
    }
    
    // 计算工作区域内的可用空间
    Paths64 ComputeFreeSpace(const Path64& workspace,
                              const Paths64& obstacles,
                              const Path64& robotShape) const {
        Paths64 expandedObstacles = ComputeSafeZone(obstacles, robotShape);
        return Difference(Paths64{workspace}, expandedObstacles, FillRule::NonZero);
    }
};

5.5 性能对比

5.5.1 矩形裁剪 vs 布尔交集

#include <chrono>

void BenchmarkRectClip() {
    // 创建测试数据
    Paths64 subject;
    for (int i = 0; i < 100; i++) {
        subject.push_back(MakeRandomPolygon());
    }
    
    Rect64 clipRect(0, 0, 500, 500);
    Paths64 clipPolygon;
    clipPolygon.push_back(MakePath({0, 0, 500, 0, 500, 500, 0, 500}));
    
    // 测试矩形裁剪
    auto start1 = std::chrono::high_resolution_clock::now();
    for (int i = 0; i < 1000; i++) {
        Paths64 result = RectClip(clipRect, subject);
    }
    auto end1 = std::chrono::high_resolution_clock::now();
    
    // 测试布尔交集
    auto start2 = std::chrono::high_resolution_clock::now();
    for (int i = 0; i < 1000; i++) {
        Paths64 result = Intersect(subject, clipPolygon, FillRule::NonZero);
    }
    auto end2 = std::chrono::high_resolution_clock::now();
    
    auto rectTime = std::chrono::duration_cast<std::chrono::milliseconds>(end1 - start1);
    auto boolTime = std::chrono::duration_cast<std::chrono::milliseconds>(end2 - start2);
    
    std::cout << "矩形裁剪: " << rectTime.count() << " ms" << std::endl;
    std::cout << "布尔交集: " << boolTime.count() << " ms" << std::endl;
}

典型结果显示,矩形裁剪比布尔交集快2-5倍。

5.6 本章小结

本章我们学习了Clipper2的两个专门几何操作:

  1. 矩形裁剪
    • RectClip和RectClipLines函数
    • RectClip64和RectClipLines64类
    • 性能优势和使用场景
    • 地图瓦片生成、视窗裁剪等应用
  2. 闵可夫斯基运算
    • 闵可夫斯基和与闵可夫斯基差
    • MinkowskiSum和MinkowskiDiff函数
    • 碰撞检测、机器人路径规划等应用
    • 性能考虑和优化策略
  3. 综合应用
    • 地图瓦片系统
    • 碰撞检测系统
    • 扫掠体积计算
    • 安全区域计算

在下一章中,我们将学习Clipper2的高级应用技巧和性能优化方法。