znlgis 博客

GIS开发与技术分享

第18章：多边形偏移算法实现

18.1 概述

本章将深入分析 ClipperOffset 的核心算法实现，包括法向量计算、各种连接类型的处理以及圆弧近似等技术细节。

18.2 GetUnitNormal 方法

计算边的单位法向量：

internal static DoublePoint GetUnitNormal(IntPoint pt1, IntPoint pt2)
{
    double dx = (pt2.X - pt1.X);
    double dy = (pt2.Y - pt1.Y);
    
    if ((dx == 0) && (dy == 0)) 
        return new DoublePoint();

    double f = 1 * 1.0 / Math.Sqrt(dx * dx + dy * dy);
    dx *= f;
    dy *= f;

    return new DoublePoint(dy, -dx);
}

18.2.1 数学原理

边方向向量: (dx, dy)
单位化: (dx/|d|, dy/|d|)  其中 |d| = √(dx² + dy²)
旋转90°得法向量: (dy/|d|, -dx/|d|)

      法向量
         ↑
         │
    pt1 ─────► pt2
         边方向

法向量指向边的左侧，用于确定偏移方向。

18.3 DoOffset 主方法

private void DoOffset(double delta)
{
    m_destPolys = new Paths();
    m_delta = delta;

    // 计算圆弧参数
    if (ArcTolerance <= 0.0)
        m_miterLim = def_arc_tolerance;
    else if (ArcTolerance > Math.Abs(delta) * def_arc_tolerance)
        m_miterLim = Math.Abs(delta) * def_arc_tolerance;
    else
        m_miterLim = ArcTolerance;

    // 计算每弧度的步数
    double steps = Math.PI / Math.Acos(1 - m_miterLim / Math.Abs(delta));
    m_sin = Math.Sin(two_pi / steps);
    m_cos = Math.Cos(two_pi / steps);
    m_StepsPerRad = steps / two_pi;
    
    if (delta < 0.0) m_sin = -m_sin;

    // 计算斜接限制
    m_miterLim = MiterLimit;
    if (m_miterLim < 2.0) m_miterLim = 2.0;
    m_miterLim = 2.0 / (m_miterLim * m_miterLim);

    // 处理每个输入多边形
    for (int i = 0; i < m_polyNodes.ChildCount; i++)
    {
        PolyNode node = m_polyNodes.Childs[i];
        m_srcPoly = node.m_polygon;

        int len = m_srcPoly.Count;
        if (len == 0 || (delta <= 0 && (len < 3 || 
            node.m_endtype != EndType.etClosedPolygon)))
            continue;

        m_destPoly = new Path();

        if (len == 1)
        {
            // 单点特殊处理
            if (node.m_jointype == JoinType.jtRound)
            {
                // 创建圆形
                double X = 1.0, Y = 0.0;
                for (int j = 1; j <= steps; j++)
                {
                    m_destPoly.Add(new IntPoint(
                        Round(m_srcPoly[0].X + X * delta),
                        Round(m_srcPoly[0].Y + Y * delta)));
                    double X2 = X;
                    X = X * m_cos - m_sin * Y;
                    Y = X2 * m_sin + Y * m_cos;
                }
            }
            else
            {
                // 创建方形
                double X = -1.0, Y = -1.0;
                for (int j = 0; j < 4; ++j)
                {
                    m_destPoly.Add(new IntPoint(
                        Round(m_srcPoly[0].X + X * delta),
                        Round(m_srcPoly[0].Y + Y * delta)));
                    if (X < 0) X = 1;
                    else if (Y < 0) Y = 1;
                    else X = -1;
                }
            }
            m_destPolys.Add(m_destPoly);
            continue;
        }

        // 构建法向量
        m_normals.Clear();
        m_normals.Capacity = len;
        for (int j = 0; j < len - 1; j++)
            m_normals.Add(GetUnitNormal(m_srcPoly[j], m_srcPoly[j + 1]));
        
        if (node.m_endtype == EndType.etClosedLine || 
            node.m_endtype == EndType.etClosedPolygon)
            m_normals.Add(GetUnitNormal(m_srcPoly[len - 1], m_srcPoly[0]));
        else
            m_normals.Add(new DoublePoint(m_normals[len - 2]));

        // 根据端点类型处理
        if (node.m_endtype == EndType.etClosedPolygon)
        {
            int k = len - 1;
            for (int j = 0; j < len; j++)
                OffsetPoint(j, ref k, node.m_jointype);
            m_destPolys.Add(m_destPoly);
        }
        else if (node.m_endtype == EndType.etClosedLine)
        {
            int k = len - 1;
            for (int j = 0; j < len; j++)
                OffsetPoint(j, ref k, node.m_jointype);
            m_destPolys.Add(m_destPoly);
            
            // 第二遍（反向）
            m_destPoly = new Path();
            DoublePoint n = m_normals[len - 1];
            for (int j = len - 1; j > 0; j--)
                m_normals[j] = new DoublePoint(-m_normals[j - 1].X, -m_normals[j - 1].Y);
            m_normals[0] = new DoublePoint(-n.X, -n.Y);
            
            k = 0;
            for (int j = len - 1; j >= 0; j--)
                OffsetPoint(j, ref k, node.m_jointype);
            m_destPolys.Add(m_destPoly);
        }
        else
        {
            // 开放路径
            // ... 处理开放端点 ...
        }
    }
}

18.4 OffsetPoint 方法

处理单个顶点的偏移：

internal void OffsetPoint(int j, ref int k, JoinType jointype)
{
    // 计算角度的正弦值
    m_sinA = (m_normals[k].X * m_normals[j].Y - m_normals[j].X * m_normals[k].Y);

    if (Math.Abs(m_sinA * m_delta) < 1.0)
    {
        // 角度太小，检查是否共线或 180° 转向
        double cosA = (m_normals[k].X * m_normals[j].X + m_normals[j].Y * m_normals[k].Y);
        
        if (cosA > 0)
        {
            // 角度 < 90°
            m_destPoly.Add(new IntPoint(Round(m_srcPoly[j].X + m_normals[k].X * m_delta),
                                        Round(m_srcPoly[j].Y + m_normals[k].Y * m_delta)));
            return;
        }
    }
    else if (m_sinA > 1.0) m_sinA = 1.0;
    else if (m_sinA < -1.0) m_sinA = -1.0;

    if (m_sinA * m_delta < 0)
    {
        // 凸角（外侧）
        m_destPoly.Add(new IntPoint(Round(m_srcPoly[j].X + m_normals[k].X * m_delta),
                                    Round(m_srcPoly[j].Y + m_normals[k].Y * m_delta)));
        m_destPoly.Add(m_srcPoly[j]);
        m_destPoly.Add(new IntPoint(Round(m_srcPoly[j].X + m_normals[j].X * m_delta),
                                    Round(m_srcPoly[j].Y + m_normals[j].Y * m_delta)));
    }
    else
    {
        // 凹角（内侧）
        switch (jointype)
        {
            case JoinType.jtMiter:
                DoMiter(j, k);
                break;
            case JoinType.jtSquare:
                DoSquare(j, k);
                break;
            case JoinType.jtRound:
                DoRound(j, k);
                break;
        }
    }
    k = j;
}

18.4.1 m_sinA 计算

m_sinA = (m_normals[k].X * m_normals[j].Y - m_normals[j].X * m_normals[k].Y);

这是两个法向量的叉积，等于 sin(θ)，其中 θ 是两条边的夹角。

        n[k]
         ↑
         │
    ─────●───── 
         │
         ↓
        n[j]

sinA > 0: 凹角（内侧需要连接）
sinA < 0: 凸角（外侧需要连接）
sinA ≈ 0: 接近共线

18.5 DoMiter 方法

斜接连接：

internal void DoMiter(int j, int k)
{
    double q = m_delta / (1 + (m_normals[j].X * m_normals[k].X + 
                               m_normals[j].Y * m_normals[k].Y));
    
    // 检查斜接限制
    double miterDistance = q / m_delta;
    if (miterDistance > m_miterLim)
    {
        // 超出限制，改用方形连接
        DoSquare(j, k);
        return;
    }
    
    m_destPoly.Add(new IntPoint(
        Round(m_srcPoly[j].X + (m_normals[k].X + m_normals[j].X) * q),
        Round(m_srcPoly[j].Y + (m_normals[k].Y + m_normals[j].Y) * q)));
}

18.5.1 斜接点计算

两条边的偏移线相交：

        偏移线1
           ╲
            ╲ ← 交点 (斜接点)
             ╲
    原点 ●────────
              ╱
             ╱
            ╱
         偏移线2

交点 = 原点 + (n1 + n2) * q
其中 q = delta / (1 + n1·n2) = delta / (1 + cosθ)

18.6 DoSquare 方法

方形连接：

internal void DoSquare(int j, int k)
{
    double dx = Math.Tan(Math.Atan2(m_sinA,
        m_normals[k].X * m_normals[j].X + m_normals[k].Y * m_normals[j].Y) / 4);
    
    m_destPoly.Add(new IntPoint(
        Round(m_srcPoly[j].X + m_delta * (m_normals[k].X - m_normals[k].Y * dx)),
        Round(m_srcPoly[j].Y + m_delta * (m_normals[k].Y + m_normals[k].X * dx))));
    m_destPoly.Add(new IntPoint(
        Round(m_srcPoly[j].X + m_delta * (m_normals[j].X + m_normals[j].Y * dx)),
        Round(m_srcPoly[j].Y + m_delta * (m_normals[j].Y - m_normals[j].X * dx))));
}

18.6.1 方形连接图解

            pt1
           ╱
         ╱   ← 角平分线方向
       ╱
    ● ───────────────
      ╲           pt2
       ╲
        ╲

在角平分线方向创建两个点，形成 45° 斜切

18.7 DoRound 方法

圆形连接：

internal void DoRound(int j, int k)
{
    double a = Math.Atan2(m_sinA,
        m_normals[k].X * m_normals[j].X + m_normals[k].Y * m_normals[j].Y);
    int steps = Math.Max((int)Round(m_StepsPerRad * Math.Abs(a)), 1);

    double X = m_normals[k].X, Y = m_normals[k].Y, X2;
    for (int i = 0; i < steps; ++i)
    {
        m_destPoly.Add(new IntPoint(
            Round(m_srcPoly[j].X + X * m_delta),
            Round(m_srcPoly[j].Y + Y * m_delta)));
        X2 = X;
        X = X * m_cos - m_sin * Y;
        Y = X2 * m_sin + Y * m_cos;
    }
    m_destPoly.Add(new IntPoint(
        Round(m_srcPoly[j].X + m_normals[j].X * m_delta),
        Round(m_srcPoly[j].Y + m_normals[j].Y * m_delta)));
}

18.7.1 圆弧近似

使用旋转矩阵生成圆弧上的点：

[X']   [cos  -sin] [X]
[Y'] = [sin   cos] [Y]

每一步旋转 2π/steps 弧度

    ●────●────●
   ╱            ╲
  ●              ●
  │              │
  ●              ●
   ╲            ╱
    ●────●────●

18.8 开放路径端点处理

// 开放路径的端点处理
int k = 0;
for (int j = 1; j < len - 1; ++j)
    OffsetPoint(j, ref k, node.m_jointype);

IntPoint pt1;
int m;
// 处理末端
switch (node.m_endtype)
{
    case EndType.etOpenButt:
        {
            m = len - 1;
            pt1 = new IntPoint((cInt)Round(m_srcPoly[m].X + m_normals[m].X * m_delta),
                              (cInt)Round(m_srcPoly[m].Y + m_normals[m].Y * m_delta));
            m_destPoly.Add(pt1);
            pt1 = new IntPoint((cInt)Round(m_srcPoly[m].X - m_normals[m].X * m_delta),
                              (cInt)Round(m_srcPoly[m].Y - m_normals[m].Y * m_delta));
            m_destPoly.Add(pt1);
        }
        break;
        
    case EndType.etOpenSquare:
        DoSquare(len - 1, len - 1);
        break;
        
    case EndType.etOpenRound:
        DoRound(len - 1, len - 1);
        break;
}

18.8.1 平头端点

// 简单地在末端添加两个点
pt1 = srcPoly[m] + normal[m] * delta;  // 一侧
pt2 = srcPoly[m] - normal[m] * delta;  // 另一侧

18.8.2 方头端点

// 复用 DoSquare，传入相同的索引
DoSquare(len - 1, len - 1);

18.8.3 圆头端点

// 复用 DoRound，创建半圆
DoRound(len - 1, len - 1);

18.9 圆弧步数计算

double steps = Math.PI / Math.Acos(1 - m_miterLim / Math.Abs(delta));

18.9.1 数学推导

为了使圆弧近似误差不超过 ArcTolerance：

误差 = r * (1 - cos(θ/2))
其中 r = delta, θ = 步进角度

要求: r * (1 - cos(θ/2)) ≤ tolerance
解得: θ ≤ 2 * arccos(1 - tolerance/r)

步数 = π / (θ/2) = π / arccos(1 - tolerance/r)

18.10 旋转矩阵优化

// 预计算旋转参数
m_sin = Math.Sin(two_pi / steps);
m_cos = Math.Cos(two_pi / steps);

// 在圆弧生成中使用
for (int i = 0; i < steps; ++i)
{
    // 添加当前点
    m_destPoly.Add(new IntPoint(
        Round(m_srcPoly[j].X + X * m_delta),
        Round(m_srcPoly[j].Y + Y * m_delta)));
    
    // 旋转到下一个位置
    X2 = X;
    X = X * m_cos - m_sin * Y;
    Y = X2 * m_sin + Y * m_cos;
}

这比每次调用三角函数更高效。

18.11 性能优化

18.11.1 法向量预计算

m_normals.Clear();
m_normals.Capacity = len;
for (int j = 0; j < len - 1; j++)
    m_normals.Add(GetUnitNormal(m_srcPoly[j], m_srcPoly[j + 1]));

18.11.2 避免重复计算

// 在 OffsetPoint 中使用预计算的法向量
// 而不是每次重新计算

18.11.3 使用 Clipper 清理

// 偏移可能产生自相交，使用 Clipper 清理
Clipper clpr = new Clipper();
clpr.AddPaths(m_destPolys, PolyType.ptSubject, true);
clpr.Execute(ClipType.ctUnion, solution, ...);

18.12 本章小结

本章详细分析了 ClipperOffset 的核心算法：

  1. 法向量计算：
    • GetUnitNormal 计算边的单位法向量
    • 法向量用于确定偏移方向
  2. 连接类型实现：
    • DoMiter：斜接连接（带限制）
    • DoSquare：方形连接
    • DoRound：圆弧连接
  3. 端点处理：
    • 平头、方头、圆头
    • 闭合路径的循环处理
  4. 圆弧近似：
    • 基于 ArcTolerance 计算步数
    • 使用旋转矩阵优化性能
  5. 后处理：
    • 使用 Clipper 清理自相交
上一章：ClipperOffset详解 返回目录 下一章：辅助函数与工具
↑