边缘崩缺检测——滚球法

应用场景

在进行图像边缘崩缺检测时，我们常使用外接矩形法或凸包法，如果图示：

图中绿色框区域为边缘崩缺位置，红色边框是使用外接矩形法，或凸包法检测出来的物体的边缘轮廓。但当检测物体不是规整的形状，且带有弧度时，如图所示，绿色部分为物体的崩缺，此时外接矩形法和凸包法则不能适用于这种场景。

此时，我们可以使用本文即将介绍的滚球法。如下图所示，黄色点为物体边缘轮廓的下采样点，圆球沿着A点逆时针旋转，第一次经过的是C点，将C点存入最终的轮廓点result。圆球再沿C旋转，第一次经过的是H点，将H点存入最终的轮廓点result。循环以上步骤，最终的 result 中未包含 D,E,F,G 点，此时，连接 result 中的所有点，即可得到没有崩缺时的物体区域。

算法原理

1. 思路

对于任何点集，我们把这些点想象为钉在平面上的钉子。假如拿一个半径大于一定值的球去从边界接近这个钉群，我们可以用这个球在这些钉子群的边界滚一圈，每次滚动，球都能接触到两个钉子而被卡住。这个思路要求一个合法的R，R太小就没法形成一个闭合的图形。由于我们讨论问题的初衷是要形成一个合适的多边形而不是0个或多个，这样对R的选择就应该有一个限制，太小的R必然出不了结果，这里姑且假设给的R值是合适的。此过程若形成一个多边形，则多变形的最长的边一定小于球的直径。如下图显示了这个滚球的过程：

假设刚开始，圆球经过A、Z两点（对应的圆心为E点），圆球沿A点逆时针旋转，先经过C点（对应的圆心为F点），再经过B点（对应的圆心为G点）。我们可以通过比较向量AE到向量AF、AG逆时针旋转的角度，来判断圆球先经过哪个点。

2. 步骤

根据上面的思路，我们可以将滚球法算法归纳为一下步骤：

1. 对所有轮廓点进行逆时针排序，并从小到大分配编号。根据设定的半径r，统计每个轮廓点在半径为 2r 范围内的所有点（因为当两个点距离大于2R时，滚球不可能同时经过这两个点）；
2. 选择一个起始点Point(x,y)，这里，我选择了Y坐标最大时，X坐标最大的点，此时最初的圆心坐标为Point(x, y-r)（也可以选择其他的点作为起始点，对应的起始圆心坐标需要调整）；
3. 根据设定的半径，沿着起始点逆时针旋转圆球，通过比较两个向量之间沿逆时针的夹角，判断圆球先经过哪个点。将第一次经过的点存入最终的轮廓点集result，并设为新的旋转点;
4. 迭代执行步骤3，当圆球第一次经过的点编号小于旋转点编号，或第一次经过的点已经在 result 中时，停止迭代，得到的 result 点集即为最终的边缘轮廓点集。

将 result 点集（所有红色点）用蓝线连接，即可得没有崩缺时的物体区域。其中，黄色点为最初物体轮廓下采样点。用蓝线包围区域与原物体区域相减，即可得到崩缺区域。

重点代码详解

1. 判断两个点是否逆时针排序

通过判断两个点A，B与轮廓中心点O，组成的向量OA，OB的关系，即可判断从A到B是否是逆时针排序

// 逆时针flag
bool inverseFlag = true;

// OA , OB 向量
int oax = points[0].X - center.X;
int oay = points[0].Y - center.Y;
int obx = points[1].X - center.X;
int oby = points[1].Y - center.Y;

auto atan2a = std::atan2(oay, oax);
auto atan2b = std::atan2(oby, obx);
if (atan2a < 0) atan2a = atan2a + CV_PI * 2;
if (atan2b < 0) atan2b = atan2b + CV_PI * 2;   // 小于0的角度需要变换下，使得它可以进行总体的角度排序；
if (atan2a > atan2b)
    inverseFlag = true;

2. 根据两点和半径，计算圆心

假设圆球在沿轮廓中心点 Center 逆时针滚动过程中，圆心从点O1滚动到O2，再到O3，圆心在O1时，同时通过A，C两点，圆心在O2，O3时，同时通过A，B两点。 而O1才是我们想要的圆心坐标，此时可以通过计算向量AO1与AO2，AO3之间的夹角，夹角小于180度时对应的圆心即为所得。

// 根据两点及半径计算圆心
CPoint RollingBall::CalCenterByPtsAndRadius(CPoint a, CPoint b, double r)
{
    CPoint circle(-1, -1);

    double x1 = a.X;
    double y1 = a.Y;
    double x2 = b.X;
    double y2 = b.Y;

    if (x1 == x2) {
        double y_c = (y1 + y2) / 2;
        double x_c1 = x1 + sqrt(pow(r, 2) - pow((y1 - y2) / 2, 2));
        double x_c2 = x1 - sqrt(pow(r, 2) - pow((y1 - y2) / 2, 2));
        circle.Y = y_c;
        circle.X = x_c1;
        double angle = GetVectorRadian(a, circle, b);
        if (angle < 180) { return circle; }
        else {
            circle.X = x_c2;
            return circle;
        }
    }

    if (y1 == y2) {
        double x_c = (x1 + x2) / 2;
        double y_c1 = y1 + sqrt(pow(r, 2) - pow((x1 - x2) / 2, 2));
        double y_c2 = y1 - sqrt(pow(r, 2) - pow((x1 - x2) / 2, 2));
        circle.X = x_c;
        circle.Y = y_c1;
        double angle = GetVectorRadian(a, circle, b);
        if (angle < 180) { return circle; }
        else {
            circle.Y = y_c2;
            return circle;
        }
    }

    double x_mid = (x1 + x2) / 2;
    double y_mid = (y1 + y2) / 2;
    double k = (y1 - y2) / (x1 - x2);
    double m = -1 / k;
    double A = 1 + pow(m, 2);
    double B = 2 * m * (y_mid - m * x_mid - y1) - 2 * x1;
    double C = pow((y_mid - m * x_mid - y1), 2) + pow(x1, 2) - pow(r, 2);
    double x_c1 = (-B + sqrt(B * B - 4 * A * C)) / (2 * A);
    double x_c2 = (-B - sqrt(B * B - 4 * A * C)) / (2 * A);
    double y_c1 = m * (x_c1 - x_mid) + y_mid;
    double y_c2 = m * (x_c2 - x_mid) + y_mid;

    circle.X = x_c1;
    circle.Y = y_c1;
    double angle = GetVectorRadian(a, circle, b);

    if (angle < 180) {
        return circle;
    }
    else {
        circle.X = x_c2;
        circle.Y = y_c2;
        return circle;
    }
}


// 计算两个向量逆时针旋转角度
double RollingBall::GetVectorRadian(CPoint origin, CPoint p1, CPoint p2)
{
    CPoint op1(p1.X - origin.X, origin.Y - p1.Y);
    CPoint op2(p2.X - origin.X, origin.Y - p2.Y);

    double op1_dot_op2 = op1.X * op2.X + op1.Y * op2.Y;
    double op1_cross_op2 = op1.X * op2.Y - op1.Y * op2.X;

    double angle = atan2(op1_cross_op2, op1_dot_op2) * 180 / CV_PI;
    if (angle < 0) { angle = 360 + angle; }

    return angle;
}

3. 获取下一点轮廓

// 获取下一个轮廓点
int RollingBall::GetNextCPoint(int pre, int cur, vector<int> list, int radius)
{
    CPoint prePt = _points[pre];
    CPoint curPt = _points[cur];
    if (pre == -1) {
        prePt.X = curPt.X - 1;
        prePt.Y = curPt.Y;
    }

    CPoint originCenter = CalCenterByPtsAndRadius(prePt, curPt, radius);    //计算起始点对应的圆心

    int index = -1;
    float angle = 360.0;
    for (int j = 0; j < list.size(); j++)
    {
        if (_signs[list[j]])
        {
            continue;
        }
        CPoint tmpPt = _points[list[j]];
        CPoint circleCenter = CalCenterByPtsAndRadius(curPt, tmpPt, radius);
        double tmpAngle = GetVectorRadian(curPt, originCenter, circleCenter);
        if (tmpAngle < angle) {
            index = list[j];
            angle = tmpAngle;
        }
    }
    return index;
}

参考：
https://blog.csdn.net/u013279723/article/details/108085334
https://blog.csdn.net/KiterPAN/article/details/134700455