Opencv源码之平面点集的最小包围圆

平面点集的最小包围圆

--Cracent整理 2016.5.28

目录

1、问题背景.... 1

2、算法及原理.... 1

3、算法（摘自OPENCV）... 1

4、基础数学知识.... 7

三角形的外心.... 7

三角形的三条垂直平分线必交于一点.... 7

三角形的外心的性质.... 7

三角形的内心.... 8

证明.... 8

性质.... 8

三角形的垂心.... 8

三角形的三条高必交于一点.... 9

三角形的垂心的性质.... 9

三角形的重心.... 10

三角形的三条中线必交于一点.... 10

三角形的重心的性质.... 11

三角形的旁心.... 11

介绍.... 11

三角形旁心的性质.... 11

欧拉线.... 12

欧拉线的证法1.. 12

欧拉线的证法2.. 13

欧拉线的证法3.. 13

 

 

1、问题背景

考察固定在工作平台上的一直机械手，要捡起散落在不同位置的多个零件，并送到别的地方。那么，这只机械手的底座应该选在哪里呢？根据直觉，应该选在机械手需够着的那些位置的“中心”。准确地讲，也就是包围这些点的那个最小圆的圆心----该位置的好处是，可使机械手的底座到它需要够着的那些点的最大距离最小化。于是可得如下问题：给定由平面上n个点所组成的一个集合P（对应于机械手需要够着的工作平台的那些位置），试找出P的最小包围圆（smallest enclosing disc）----亦即，包含P中所有点、半径最小的那个圆。这个最小包围圆必然是唯一的。

 

 2、算法及原理

算法介绍：我们本次算法的设计是基于这样一个简单直观的性质：在既定的给定点条件下，如果引入一张新的半平面，只要此前的最优解顶点（即唯一确定最小包围圆的几个关键顶点）能够包含于其中，则不必对此最优解进行修改，亦即此亦为新点集的最优解；否则，新的最优解顶点必然位于这个新的半空间的边界上。

定理可以通过反证法证明。

于是，基于此性质，我们便可得到一个类似于线性规划算法的随机增量式算法。定义Di为相对于pi的最小包围圆。此算法实现的关键在于对于pi∉Di-1时的处理。显然，如果pi∈Di-1，则Di= Di-1；否则，需要对Di另外更新。而且，Di的组成必然包含了pi；因此，此种情况下的最小包围圆是过pi点且覆盖点集{ p1，p2，p3……pi-1}的最小包围圆。则仿照上述处理的思路，Di={ p1，pi}，逐个判断点集{ p2，p3……pi-1}，如果存在pj∉Di，则Di={pj，pi}。同时，再依次对点集{ p1，p2，p3……pj-1}判断是否满足pk∈Di，若有不满足，则Di={pk，pj，pi}。由于，三点唯一地确定一个圆，故而，只需在此基础上判断其他的点是否位于此包围圆内，不停地更新pk。当最内层循环完成时，退出循环，转而更新pj；当次内层循环结束时，退出循环，更新pi。当i=n时，表明对所有的顶点均已处理过，此时的Dn即表示覆盖了给定n个点的最小包围圆。

 

3、算法（摘自OPENCV）

主函数：minEnclosingCircle_jbl(InputArray _points, Point2f& _center, float& _radius)

 

static float innerProduct(Point2f &v1, Point2f &v2)

{

    return v1.x * v2.y - v1.y * v2.x;

}

 

//三点定圆法

static void findCircle3pts(Point2f *pts, Point2f &center, float &radius)

{

    // two edges of the triangle v1, v2

    Point2f v1 = pts[1] - pts[0];

    Point2f v2 = pts[2] - pts[0];

 

    if (innerProduct(v1, v2) == 0.0f)

    {

        // v1, v2 colineation, can not determine a uniquecircle

        // find the longtest distance as diameter line

        float d1 = (float)norm(pts[0] - pts[1]);

        float d2 = (float)norm(pts[0] - pts[2]);

        float d3 = (float)norm(pts[1] - pts[2]);

        if (d1 >= d2 && d1 >= d3)

        {

            center = (pts[0] + pts[1]) / 2.0f;

            radius = (d1 / 2.0f);

        }

        else if (d2 >= d1 && d2 >= d3)

        {

            center = (pts[0] + pts[2]) / 2.0f;

            radius = (d2 / 2.0f);

        }

        else if (d3 >= d1 && d3 >= d2)

        {

            center = (pts[1] + pts[2]) / 2.0f;

            radius = (d3 / 2.0f);

        }

    }

    else

    {

        // center is intersection of midperpendicular linesof the two edges v1, v2

        // a1*x + b1*y = c1 where a1 = v1.x, b1 = v1.y

        // a2*x + b2*y = c2 where a2 = v2.x, b2 = v2.y

        Point2f midPoint1 = (pts[0] + pts[1]) / 2.0f;

        float c1 = midPoint1.x * v1.x + midPoint1.y * v1.y;

        Point2f midPoint2 = (pts[0] + pts[2]) / 2.0f;

        float c2 = midPoint2.x * v2.x + midPoint2.y * v2.y;

        float det = v1.x * v2.y - v1.y * v2.x;

        float cx = (c1 * v2.y - c2 * v1.y) / det;

        float cy = (v1.x * c2 - v2.x * c1) / det;

        center.x = (float)cx;

        center.y = (float)cy;

        cx -= pts[0].x;

        cy -= pts[0].y;

        radius = (float)(std::sqrt(cx *cx + cy * cy));

    }

}

 

const float EPS = 1.0e-4f;

 

static voidfindEnclosingCircle3pts_orLess_32f(Point2f *pts, int count, Point2f &center, float &radius)

{

    switch (count)

    {

    case 1:

        center = pts[0];

        radius = 0.0f;

        break;

    case 2:

        center.x = (pts[0].x + pts[1].x) / 2.0f;

        center.y = (pts[0].y + pts[1].y) / 2.0f;

        radius = (float)(norm(pts[0] - pts[1]) / 2.0);

        break;

    case 3:

        findCircle3pts(pts, center, radius);

        break;

    default:

        break;

    }

 

    radius += EPS;

}

 

template<typename PT>

static void findThirdPoint(const PT *pts, int i, int j, Point2f &center, float &radius)

{

    center.x = (float)(pts[j].x + pts[i].x) / 2.0f;

    center.y = (float)(pts[j].y + pts[i].y) / 2.0f;

    float dx = (float)(pts[j].x - pts[i].x);

    float dy = (float)(pts[j].y - pts[i].y);