多边形裁剪：Sutherland Hodgman算法

前言

最近一段时间一直在专研物理引擎，这真是一个集代数，平面几何，物理和计算机于一身的交叉领域，要一层层拨开里面的迷雾真的不容易。

前段时间我阅读了某个仿box2d的c++物理引擎clib的源码，虽说里面的代码写得乱七八糟的，但是也有一点收获。后来就转去看matter.js的源码了，感叹没有早点发现matter.js，里面的代码真的堪称工程典范，条理清晰，井井有条。

跑题了，说回今天要介绍的内容。在看完clib的碰撞检测部分后，我意识到我之前有一个地方说错了：

最近在看box2d的源码，看得好累。发现box2d的碰撞检测不止用到了SAT（分离轴）算法，还有GJK算法和V-clip算法（连google都找不到的冷门算法，不知道具体原理）。box2d貌似把3种算法糅合起来了，根本看不懂。

额，其实并没有V-clip这个算法（怪不得google搜不到），box2d里面的V-clip函数其实是指Vertex Clip，即多边形裁剪，而该多边形裁剪算法的真正名字叫Sutherland Hodgman算法。

Sutherland Hodgman算法

Sutherland Hodgman是一个用于在指定区域裁剪多边形的算法，注意，裁剪不是分割，这是两回事。sh的基本思想其实很简单，核心就是根据上个顶点和当前顶点于裁剪区域的关系进行一步步迭代裁剪。Sutherland Hodgman把上个顶点和当前顶点于裁剪区域的关系分成四种情况：

上一个顶点在裁剪区域外，当前顶点在裁剪区域内
上一个顶点在裁剪区域内，当前顶点在裁剪区域外
上一个和当前顶点都在顶点在裁剪区域外
上一个和当前顶点都在顶点在裁剪区域内

而根据上面的四种情况，Sutherland Hodgman对多边形所有顶点进行遍历，选择哪些顶点需要保留，哪些顶点需要遗弃：

情况1：保留current和last连线与裁剪线的交点和current
情况2：保留current和last连线与裁剪线的交点和last
情况3：遗弃顶点current
情况3：保留顶点current

对剪辑区域的每条剪辑边进行以上的操作，便可以完成多边形的剪辑。

求两条线段的交点

可以看到，Sutherland Hodgman在情况1和情况2下需要求两条线段的交点，然而其实求线段交点并不是一件简单的事。

假如有两条线段a1a2，b1b2，a1，a2，b1，b2为已知顶点，要求其交点，通常情况下都是联立两条直线方程解方程组，联立之前还要把线段化成一般式，这是最常规最稳固的方法，但是这不是对计算机友好的方法，因为计算机不会联立方程，需要先人为地将x，y化为系数的表达式，然后在代码中塞入长长的一串算式。

当然这没有问题，但是不够优雅。我使用的是一种利用向量和相似三角形的求解法。首先将线段看成两个向量：向量a1a2和向量b1b2，然后建立辅助线，如下图：

首先，利用向量叉积的性质：两二维向量叉积的几何意义是两个向量围城的平行四边形的面积，求出b1点到直线a1a2的距离d1, 和b2点到直线a1a2的距离d2：

由于三角形b1b2e和三角形b1oc相似，所以有：

求出了向量b2o，又已知b2的值，便可以求出交点o的值

具体代码很简单，没有什么逻辑性的东西，只要理解了就能写出来，详细注释就不写了，意义不大：

/**
 * 求交点
 * @param line1 第一条线段
 * @param line2 第二条线段
 */
function intersection(line1: polygonVex, line2: polygonVex): vector {
    let v1 = Vector.sub(line1[1], line1[0]),
        v2 = Vector.sub(line2[1], line2[0]),
        tv1 = Vector.sub(line1[0], line2[0]),
        tv2 = Vector.sub(line1[1], line2[1]),
        d1 = Math.abs(Vector.cor(tv1, v1)/Vector.len(v1)),
        d2 = Math.abs(Vector.cor(tv2, v1)/Vector.len(v1)),
        tv3 = Vector.scl(d1/(d1 + d2), v2);

    return Vector.add(line2[0], tv3);
}

算法实现

/**
 * Sutherland Hodgman算法
 * @param polygon 要裁剪的多边形
 * @param clipPlane 裁剪区域
 */
export function SutherlandHodgman(polygon: polygonVex, clipPlane: polygonVex): polygonVex {
        // 裁剪完成的多边形
    let resultPolygon: polygonVex,
        // 保存当前在裁剪的多边形
        tmpPolygon: polygonVex,
        // 上一个顶点
        lastVertex: vector,
        // 当前顶点
        curVertex: vector,
        // 上一个顶点是否在裁剪区域内
        isLastVertexInside: boolean,
        i: number,
        j: number;

    // 初始化裁剪完成的多边形为原始多边形
    resultPolygon = polygon.slice(0);
    // 当前正在裁剪的多边形为空
    tmpPolygon = [];
    i = 0;

    // 遍历裁剪区域的边
    while(i < clipPlane.length - 1) {
        // 当前裁剪边
        let plane: polygonVex = [clipPlane[i], clipPlane[i + 1]];

        // 初始化上一个顶点为多边形的最后一个顶点
        lastVertex = resultPolygon[resultPolygon.length - 1];
        // 记录上一个顶点与裁剪区域的关系
        isLastVertexInside = isInside(plane, lastVertex);

        // 遍历多边形顶点
        for(j = 0; j < resultPolygon.length; j++) {
            // 当前顶点
            curVertex = resultPolygon[j];

            // 若当前顶点在裁剪区域内
            if(isInside(plane, curVertex)) {
                // 而上一个顶点不在裁剪区域内 => 情况1
                if(!isLastVertexInside) {
                    // 添加相交点
                    tmpPolygon.push(intersection(plane, [lastVertex, curVertex]));
                }

                // 添加当前顶点
                tmpPolygon.push(curVertex);

                // 修改标志
                isLastVertexInside = true;
            }
            // 若当前顶点不在剪辑区域内
            else {
                // 而上一顶点在剪辑内 => 情况2
                if(isLastVertexInside) {
                    // 添加相交点
                    tmpPolygon.push(intersection(plane, [lastVertex, curVertex]));
                }

                // 修改标志位
                isLastVertexInside = false;
            }

            // 当前顶点设为上一顶点
            lastVertex = curVertex;
        }

        // 完成了一条裁剪边的裁剪，将当前完成裁剪的多边形保存到裁剪完成的多边形数组，使用该数组进行下一次裁剪
        resultPolygon = tmpPolygon.slice(0);
        // 清空当前在裁剪的多边形数组
        tmpPolygon = [];
        // 下一个裁剪边
        i++;
    }

    return resultPolygon;
}   

/**
 * 判断点在线段的哪一侧
 * @param line 线段
 * @param point 被检测点
 */
function isInside(line: polygonVex, point: vector): boolean {
    let v1 = Vector.sub(line[1], line[0]),
        v2 = Vector.sub(point, line[0]);

    // < 0：左侧；> 0：右侧；= 0：点在线上
    return Vector.cor(v2, v1) <= 0;
}

/**
 * 求交点
 * @param line1 第一条线段
 * @param line2 第二条线段
 */
function intersection(line1: polygonVex, line2: polygonVex): vector {
    let v1 = Vector.sub(line1[1], line1[0]),
        v2 = Vector.sub(line2[1], line2[0]),
        tv1 = Vector.sub(line1[0], line2[0]),
        tv2 = Vector.sub(line1[1], line2[1]),
        d1 = Math.abs(Vector.cor(tv1, v1)/Vector.len(v1)),
        d2 = Math.abs(Vector.cor(tv2, v1)/Vector.len(v1)),
        tv3 = Vector.scl(d1/(d1 + d2), v2);

    return Vector.add(line2[0], tv3);
}

判断点在剪辑线的哪一侧也用到了向量叉积，看见叉积的应用在图形学还是很广泛的，但是不知道为什么中学数学没有叉积的内容。

效果展示

Sutherland Hodgman算法在碰撞检测中的应用

看到这里的小伙伴可能会问：说了这么多，这个算法跟物理引擎究竟有什么关系呢？在阅读了几个项目的源码，在结合自己踩过的一些坑后，我得出了结论：

Sutherland Hodgman主要用于求解碰撞点

如图，两个图形因为碰撞发生穿透，红点为两个图形的交点，这时不能直接判定相交点即为两个图形的碰撞点，因为这不符合现实规则，所以必须进行一次穿透修正后，相交点才是碰撞点。Sutherland Hodgman就是用作求这两个相交点的。这里的碰撞点可能与真实物理世界情况下有些误差，可以看到两个图形在修正后依然相交，这是因为在相交修正时必须保持一个允许微小穿透的slop（通常 < 0.1），用作消除图形堆叠产生的抖动，但是这种精度一般已经足够了。

---EOF---

phenomLi / Blog