8.4 移动凸多边形
(Translating a Convex Polygon)

当机器人形状为凸多边形时，我们会遇到一个棘手的问题。从一个位置移动到另一个位置可能需要旋转机器人。关于旋转的讨论我们放到 8.5 节。在这里，我们仅考虑平移运动。这项任务仍然比移动圆盘更复杂，但幸运的是，通过闵可夫斯基和(Minkowski sums)来膨胀障碍物的想法仍然适用。在讨论算法细节之前，我们先通过一个例子来解释其基本思想。

8.4.1. 闵可夫斯基和示例 (Minkowski Sum Example)

设机器人 \(R\) 为一个正方形，并选择其左下角作为参考点 \(r\)。考虑一个多边形 \(P\)，如图 8.5 所示的五边形。当 \(R\) 沿着 \(P\) 的边界 \(\partial P\) 移动时，参考点 \(r\) 描绘出 \(P^+\) 的边界，\(P^+\) 是膨胀后的障碍区域，定义了 \(r\) 无法进入的范围。请注意，这种情况与圆盘稍有不同，因为我们将参考点选在了 \(R\) 的一个角上。\(P\) 在不同边上膨胀量各不同。例如，沿边 \(e_0\) 移动时， \(P\) 和 \(P^+\) 的边界重合，因为参考点 \(r\) 有可能接触到 \(e_0\)。而且 \(P^+\) 相对于 \(P\) 的边 \(e_1\) 和 \(e_4\) 的偏移量是不同的。相对于边 \(e_1\) 的膨胀量由 \(R\) 的水平宽度决定，相对于 \(e_4\) 取决于 \(R\) 的垂直高度。还要注意在凹顶点(reflex vertex)附近的情况。如果我们沿着与之相连的两条边的全长移动，\(r\) 会描绘出一条自相交的路径，尽管如此，该路径的外侧仍然精确地表征了 \(R\) 靠近 \(P\) 的几何极限。

8.4.2. 闵可夫斯基加法 (Minkowski Addition)

当 \(R\) 是圆盘时，我们曾论证 \(P^+ = P \oplus R\)。但这显然不适用于图 8.5 的情况。例如，\(P \oplus R\) 会凸出 \(P\) 的边 \(e_0\) 之外，但 \(P^+\) 并不会。合适的计算方法应该是取 \(P\) 与 \(R\) 关于参考点 \(r\) 的反射的闵可夫斯基和。take the Minkowski sum of \(P\) with a reflection of \(R\) through the refence point \(r\). 由于 \(r\) 在闵可夫斯基和的公式中作为原点，所以 \(R\) 的反射版本就是 \(-R\)，即 \(R\) 中的每一点都取反。需要反射的原因是，\(R\) 中的任意点 \(p\)，比如说 \(R\) 的右上角，都有将 \(r\) 推离 \(\partial P\) 距离 \(-p\) 的效果。由此可见，图 8.5 中的 \(P^+\) 就是 \(P \oplus -R\)。注意，因为圆盘关于其圆心是中心对称的，所以 \(R = -R\)，因此这个新公式与我们在上一节中的表述是一致的。鉴于"闵可夫斯基剪发(Minkowski subtraction)"这个术语有时用于指代另一个概念(Guggenheimer 1977)，我们将继续称之为闵可夫斯基加法。

我们用一个稍显通用的命题来形式化上述讨论。

定理 8.4.1. 设 \(R\) 是一个区域(机器人)，\(r \in R\) 是一个参考点。\(P\) 是一个障碍物。那么区域 \(P^+ = P \oplus -R\) 是参考点 \(r\) 的禁区，具体表现为:

若平移 \(R\) 使得 \(r\) 严格位于 \(P^+\) 内部，那么 \(R\) 与 \(P\) 就会穿模，或者说是有重叠区域。
若平移 \(R\) 使得 \(r\) 位于 \(\partial P^+\) 上，那么 \(\partial R\) 与 \(\partial P\) 会接触上，或者说是相切。
若平移 \(R\) 使得 \(r\) 严格位于 \(P^+\) 外部，那么 \(R \cap P = \emptyset\)，即不相交。

这是一种更一般的表述，因为 \(R\) 和 \(P\) 都不需要是凸的，甚至不需要是多边形。但在本节中，我们仍然假设两者都是多边形，且 \(R\) 是凸的。为了方便起见，我们将取 \(r\) 设为原点。

8.4.3. 构建闵可夫斯基和 (Constructing the Minkowski Sum)

我们在此概述一种构建两个多边形闵可夫斯基和的方法，其中部分内容将在下文 8.4.4 节中实现。为了将重点保持在运动规划本身，而不迷失于由此引发的各种有趣问题，这里不提供完整细节，而是将其留作习题。

我们继续使用图 8.5 中的例子。图 8.6(b) 展示了 \(P^+ = P \oplus -R\)，其中 \(P\) 的边标记为 \(0, \dots, 4\)，\(-R\) 的边标记为 \(a, \dots, d\)，两者均按逆时针遍历顺序标记，而 \(P^+\) 的边则根据是由 \(P\) 或 \(-R\) 的哪条边"生成"，来标记为 \(\{0, 1, 2, 3, 4, a, b, c, d\}\) 中的元素。因此，当 \(R\) 沿 \(P\) 的边 2 滑动时，参考点 \(r\) 描绘出 \(P^+\) 的一条平行边，我们将其标记为 2。而当 \(R\) 在 \(P\) 的边 1 和 2 交点处沿 \(c\) 垂直滑动时，我们将 \(P^+\) 的那条边标记为 \(c\)。

注意，我们已经标记了包围 \(P^+\) 的整条自相交多边形路径 \(\tau\)，包括在凹顶点(reflex vertex)附近位于 \(P^+\) 内部的边。构造 \(P^+\) 最简单的方法是先构建 \(\tau\)。它有时被称为 \(P\) 和 \(-R\) 的卷积(convolution)(Guibas, Ramshaw & Stolfi 1983)。

\(\tau\) 各边的标记规律，可以通过图 8.6(a) 所示的 \(P\) 和 \(-R\) 边向量的星形图(star diagram)来直观理解。如果我们像在当前示例中那样，将 \(P\) 视为较大的多边形，\(R\) 视为较小的多边形。那么粗略地说，\(\tau\) 包含对应于 \(P\) 各边的边，其间穿插着一些 \(-R\) 的边。实际上，人们可以看到 \(\tau\) 的标签序列 \((0, b, 1, c, d, 2, a, b, 3, c, d, 4, a)\) 包含了 \((0, 1, 2, 3, 4)\) 作为一个子序列。星形图提供了一种预测 \(-R\) 穿插标签的机制。

将每条边视为按逆时针遍历方向的向量，并将它们都移动到一个公共原点，如图 8.6(a) 所示，并归一化为单位长度。这种边向量的排列方式被称为星形图(star diagram)。现在从 0 开始，绕星形图逆时针旋转。在 \(P\) 边的索引 \(i\) 和 \(i+1\) 之间，写下遇到的所有 \(-R\) 的索引。例如，在 0 和 1 之间遇到了 \(b\)，产生了子序列 \((0, b, 1)\)。在 1 和 2 之间遇到了 \(c\) 和 \(d\)，产生了子序列 \((1, c, d, 2)\)。以此类推，当再次回到 0 时，我们就生成了 \(\tau\) 的完整序列。稍后我们将讨论具体的实现方法。

为了从 \(\tau\) 中获得闵可夫斯基和 \(P^+\)，还需要处理卷积的自交问题。这是另一个有趣但我们不打算深入探讨的问题(Guibas et al. (1983); Ramkumar (1996))。虽然细节尚不清楚，但至少可以肯定，存在一个确定的过程，根据 \(P\) 和 \(R\) 构建 \(P^+\)。现在我们直接给出在各种凸性假设下计算 \(P \oplus -R\) 的复杂度结论:

定理 8.4.2. 如果 \(P\) 有 \(n\) 个顶点，且 \(R\) 的顶点数量固定(为常数)，那么闵可夫斯基和 \(P \oplus -R\) 可以在以下时间复杂度内完成构建:

\(R\)	\(P\)	Size of Sum	时间复杂度
凸(convex)	凸(convex)	\(O(n)\)	\(O(n)\)
凸(convex)	非凸(nonconvex)	\(O(n)\)	\(O(n^2 \log n)\)
非凸(nonconvex)	非凸(nonconvex)	\(O(n^2)\)	\(O(n^2 \log n)\)

这些结果由 Guibas et al.(1983), Toussaint (1983b), Sharir (1987) 以及 Kaul, O'Connor & Srinivasan(1991) 得出。请注意，我们将机器人的大小视为具有固定数量顶点的形态，并且仅报告关于障碍物多边形的顶点数 \(n\) 的复杂度。习题 8.4.6[4]–[6] 探讨了作为机器人顶点数函数的复杂度。

8.4.4. 闵可夫斯基卷积的实现 (Implementation of Minkowski Convolution)

尽管构建卷积轨迹仅仅是构建其外边界(即闵可夫斯基和)的一半工作，但星形图(star diagram)方法足够优雅，值得我们实现一遍。所需的大部分精力都在于构建星形图。一旦有了星形图，只需通过重复累加对应的边向量，就可以轻松完成卷积轨迹的追踪。构建星形图本质上是对向量进行角度排序，我们在 3.5.5 节中已经使用 Graham 凸包算法处理过此类任务。我们将尽量遵循为该算法使用的惯例。具体而言，我们将把边向量存储在一个结构体数组中，然后使用 qsort 进行排序。这些结构体以及算法顶层逻辑的 main 函数，如 Code8.1 所示。

展开 Code 8.1. Data structures and main program for Minkowski convolution. Cf. Code 3.6.>>>

        typedef struct tPointStructure tsPoint;
        typedef tsPoint *tPoint;
        
        struct tPointStructure {
            int     vnum;
            tPointi v;
            bool    primary;
        };
        
        #define PMAX 1000                       /* Max # of points */
        typedef tsPoint tPointArray[PMAX];
        static tPointArray P;
        
        int m;      /* Total number of points in both polygons */
        int n;      /* Number of points in primary polygon */
        int s;      /* Number of points in secondary polygon */
        
        main()
        {
            tPointi p0 = {0,0};
            int j0;                               /* index of start point */
        
            j0 = ReadPoints( p0 );
            Vectorize();
            qsort(
                &P[0],                            /* pointer to 1st elem */
                m,                                /* number of elems */
                sizeof( tsPoint ),                /* size of each elem */
                Compare                           /* -1,0,+1 compare function */
            );
            Convolve( j0, p0 );
        }

其中点结构体 tPoint 与 3.5.5 节中使用的结构体仅有一处不同，它增加了一个布尔字段 primary，用于标识该点属于主(非凸)多边形 \(P\) (TRUE) 还是次(凸)多边形 \(R\)(FALSE)。该字段由函数 ReadPoints设置(此处未列出)，该函数还会利用字段 vnum 分别为两组点分配独立的编号序列。这样，当 qsort 将所有点混合在数组 P[] 中，我们仍有足够的信息来识别点的身份。主多边形的顶点数为 \(n\)，次多边形的顶点数为 \(s\)，总顶点数为 \(m = n + s\)。次多边形 \(R\) 在读取时会进行反射变换，因此后续仅存储和使用 \(-R\)。

虽然我们读入的是多边形顶点，但所有计算都是在边向量上进行的。因此，主过程的第一个实质性操作是通过调用 Vectorize(Code8.2)来计算这些边向量。这无需像图 8.6(a) 中那样归一化为单位长度。事实上，保持原样的向量更有用。我们选择将这些向量存储在原本存放顶点的同一个数组 P[] 中。我们需要小心避免误用新计算的向量覆盖了原始点，这也是引入临时变量 last 进行暂存的原因。除此之外，该计算过程还是很直接的。

展开 Code 8.2. Vectorize. See Code 7.6 for SubVec.>>>

        void Vectorize( void )
        {
            int i;
            tPointi last;   /* Holds the last vector difference. */
        
            // --- 处理第一个多边形 (索引 0 到 n-1) ---
            SubVec( P[0].v, P[n-1].v, last );
            for( i = 0; i < n-1; i++ )
                SubVec( P[i+1].v, P[i].v, P[i].v );
            P[n-1].v[X] = last[X];
            P[n-1].v[Y] = last[Y];
        
            // --- 处理第二个多边形 (索引 n 到 n+s-1) ---
            SubVec( P[n].v, P[n+s-1].v, last );
            for( i = 0; i < s-1; i++ )
                SubVec( P[n+i+1].v, P[n+i].v, P[n+i].v );
            P[n+s-1].v[X] = last[X];
            P[n+s-1].v[Y] = last[Y];
        }

对 qsort 的调用与 Code3.6 中的用法类似，不再赘述。与 Graham 算法一样，编写比较函数 Compare 是一项细致的工作。我们无法完全照搬之前的模型(Code3.5)，因为在那里的所有向量相对于原点都位于上半平面。而此处的边向量则向各个方向发散。我们采用如下的角度排序约定:沿 \(x\) 轴向左指向的向量(例如 \((-1, 0)\))具有最小角度，表示从正 \(x\) 轴逆时针旋转 \(-\pi\)。角度值随逆时针旋转而增大，一次经过 \(x\) 轴下方的半平面，经过正 \(x\) 轴，最后穿过上方的半平面。这与我们之前的约定一致，但将其扩展到了完整的 \(2\pi\) 范围。"位于左侧(left-of)"的判定仍然是决定向量先后顺序的关键。在出现平局(共线)的情况下，我们优先选择较短的向量。

尽管 left-of 是 Compare 的核心(Code8.3)，但我们大部分精力都在处理特殊情况。下半平面上的向量优先于上半平面上的向量。这可以通过比较每个向量的 v[Y] 来确定。\(x\) 轴上的向量需要特殊处理。否则，我们将进入 left-of 判定和平局裁决，其处理方式与第 3 章中使用的方法相同。

展开 Code 8.3. Compare. Cf. Code3.5.>>>

        int Compare( const void *tpi, const void *tpj )
        {
            int a;              /* AreaSign result */
            int x, y;           /* projections in 1st quadrant */
            tPoint pi, pj;      /* Recasted points */
            tPointi Origin = {0,0};
        
            pi = (tPoint)tpi;
            pj = (tPoint)tpj;
        
            /* A vector in the open upper halfplane is after
               a vector in the closed lower halfplane. */
            if      ( ( pi->v[Y] > 0 ) && (pj->v[Y] <= 0 ) )
                return 1;
            else if ( ( pi->v[Y] <= 0 ) && (pj->v[Y] > 0 ) )
                return -1;
        
            /* A vector on the x axis and one in the lower
               halfplane are handled by the Left computation below. */
        
            /* Both vectors on the x axis require special handling. */
            else if ( ( pi->v[Y] == 0 ) && (pj->v[Y] == 0 ) ) {
                if      ( ( pi->v[X] < 0 ) && ( pj->v[X] > 0 ) )
                    return -1;
                if      ( ( pi->v[X] > 0 ) && ( pj->v[X] < 0 ) )
                    return 1;
                else if ( abs(pi->v[X]) < abs(pj->v[X]) )
                    return -1;
                else if ( abs(pi->v[X]) > abs(pj->v[X]) )
                    return 1;
                else
                    return 0;
            }
        
            /* Otherwise, both in open upper halfplane, or
               both in closed lower halfplane, but not both on x axis. */
            else {
                a = AreaSign( Origin, pi->v, pj->v );
                if      (a > 0)
                    return -1;
                else if (a < 0)
                    return 1;
                else { /* Begin collinear */
                    x = abs( pi->v[X] ) - abs( pj->v[X] );
                    y = abs( pi->v[Y] ) - abs( pj->v[Y] );
                    if      ( (x < 0) || (y < 0) )
                        return -1;
                    else if ( (x > 0) || (y > 0) )
                        return 1;
                    else /* points are coincident */
                        return 0;
                } /* End collinear */
            }
        }

调用 qsort 之后，带有识别标签的边向量会按角度排序，并存储在 P[] 中。接下来的任务是按照之前描述的方式遍历星形图(star diagram)累加边向量。一个棘手的问题是起始点的选择 —— 初始点 \(p_0\) 应该是什么? 这又取决于我们从哪里开始处理。我们选择从最小角度开始，它对应于向量 \((-1, 0)\)，即使多边形边向量中可能不存在这样的向量。另一种选择是从主多边形的某个特定边向量开始。虽然这并非显而易见的，但这钟选择意味着我们应该从主多边形 \(P\) 的右上角开始，该点还参考次多边形 \(-R\) 的最右上方点进行了偏移。我们不对此进行论证。起点 \(p_0\) 有 ReadPoints 计算得出，并传递给卷积过程 Convolve，该过程首先要"移动至(moveto)" 起点，如 Code8.5 所示。

展开 Code 8.5. Convolve. Subsidiary Postscript output commands are not shown.>>>

        void Convolve( int j0, tPointi p )
        {
            int i;  /* Index into sorted edge vectors P */
            int j;  /* Primary polygon index */
        
            MoveTo_i( p );
        
            i = 0;      /* Start at angle -pi, rightward vector. */
            j = j0;     /* Start searching for j0. */
        
            do {
                /* Advance around secondary edges until next j reached. */
                while ( !(P[i].primary) && P[i].vnum == j ) {
                    if ( !P[i].primary ) {
                        AddVec( p, P[i].v, p );
                        LineTo_i( p );
                    }
                    i = (i + 1) % m;
                }
        
                /* Advance one primary edge. */
                AddVec( p, P[i].v, p );
                LineTo_i( p );
                j = (j + 1) % n;
        
            } while ( j != j0 );
        }

Convolve 随后通过递增索引 i 来遍历星形图，并通过语句 i = (i+1)%m 在必要时进行回绕。这是代码中外层的 do-while 循环。在每次迭代中，它都在寻找主多边形的下一条边向量，由索引 j 标识。该索引 j 被初始化为 \(j_0\)，即基于 \(p_0\) 的主边向量的索引。在搜索主向量 \(j\) 的过程中，它会输出遇到的每一个次向量，这对应于图 8.6(a) 中穿插的 \(-R\) 标签。这是在内层 while 循环中完成的。当内层循环结束时，算法已找到下一个 \(j\)，随即输出一条主边向量并将 \(j\) 递增。不断重复该过程，直到 \(j\) 再次回绕到 \(j_0\)。

图 8.7 展示了代码输出的一个示例。相应的星形图显示在图 8.8 中。角度最小的向量是 \((-20, 0)\)，对应于 \(R\) 的水平底边，经 \(-R\) 变换之后指向左侧。从 \(p_0\)，即 \(P\) 的最上方点，图中画圈处，出发的第一步，就是这个向左的水平"连线(lineto)"。随后，索引 \(i\) 会遍历 P[] 中所有 \(m = n + s = 22 + 8 = 30\) 个点十次之后，回到 \(j_0\)。这十次循环在最终图形的十个环路中清晰可见。从该图可以清晰的看到，从卷积到闵可夫斯基(即图8.7的外边界)和还有很多额外工作。

展开 Figure8.8. 对应图8.7的星形图. 实在太糊了，影响排版美观，默认折叠了，点击展开.>>>

8.4.5. 概念性运动规划算法 (Conceptual Motion Planning Algorithm)

现在我们回到为凸多边形机器人 \(R\) 的规划运动的问题上。同样，这里仅给出算法的粗略概念。设障碍物为 \(P_1, P_2, \ldots, P_m\)，总共有 \(n\) 个顶点。该算法包含四个步骤:

膨胀所有障碍物: \(P_i^+ = P_i \oplus -R\)。
构造它们的并集 \(P^+ = \bigcup_i P_i^+\)。
找到包含 \(s\) 和 \(t\) 的连通分量。
在该连通分量中找到一条连接 \(s\) 和 \(t\) 的路径。

图 8.9 用一个四边形机器人 \(R\) 和八个障碍物(深色阴影)说明了这一过程。

注意，在这个例子中，自由空间有三个连通分量，分别包含点 \(a, b, c\)。从图示的初始位置出发，\(a\) 是可达的，但 \(b, c\) 都不可达。因此，判断目标位置 \(t\) 是否从 \(s\) 可达，被转化为确定 \(s, t\) 是否位于自由空间的同一个连通分量中。而为机器人规划路径，则转化为在该分量内为参考点寻找一条路径。

上述四个步骤都涉及有趣的算法问题，但我们在此不深入探讨。鉴于整个方案的复杂性，Kedem & Sharir(1990) 能够实现卓越的 \(O(n \log n)\) 复杂度，充分彰显了众多致力于该问题各方面研究的学者的贡献。

定理 8.4.3. 寻找一系列平移运动，使凸多边形机器人在给定的起始和终止位置之间移动，同时避开总共有 \(n\) 个顶点的多边形障碍物，这一过程可以在 \(O(n \log n)\) 时间内完成。更准确地说，如果机器人有 \(k\) 个顶点，则复杂度为 \(O(kn \log(kn))\)。

更令人惊喜的是，对于在 3D 空间中移动的凸多面体机器人，避开多面体障碍物的问题，也已确立了大致相同的复杂度(Aronov & Sharir 1997)。

8.4.6. 练习(Exercises)

[easy] 正方形与三角形之和 Sum of square and triangle。正方形和三角形的闵可夫斯基和(Minkowski sum)最多可以有多少条边?
卷积循环 Convolution cycles。根据 \(P\) 的结构，开发一种方法来预测 Convolve(代码 8.5)在星形图上循环的次数。在图 8.6(b) 和 8.7 上测试。
星形图:非凸 Star diagram: nonconvex。探索非凸多边形 \(R\) 的卷积。星形图方法还能找到 \(\tau\) 吗？
[easy] 凸-凸 Convex-Convex。证明分别具有 \(k\) 和 \(n\) 个顶点的两个凸多边形 \(R, P\) 的闵可夫斯基和可以有 \(\Omega(k + n)\) 个顶点。
凸-非凸 Convex-Nonconvex。证明分别具有 \(k\) 和 \(n\) 个顶点的凸多边形 \(R\) 和可能是非凸的多边形 \(P\) 的闵可夫斯基和可以有 \(\Omega(kn)\) 个顶点。
非凸-非凸 Nonconvex-Nonconvex。证明分别具有 \(k\) 和 \(n\) 个顶点的两个可能是非凸的多边形 \(R, P\) 的闵可夫斯基和可以有 \(\Omega(k^2n^2)\) 个顶点。
[difficult] 凸区域的并集 Union of convex regions。\(m\) 个凸多边形共有 \(n\) 个顶点，这些凸多边形的并集，最多有多少顶点? 请用 \(n\) 和 \(m\) 的函数表示答案。
1. 首先猜测一个上界，并用用示例支持。
2. 证明该上界(Kedem, Livne, Pach & Sharir 1986)。