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

平面点集的最小包围圆--Cracent整理 2016.5.28目录1、问题背景.... 12、算法及原理.... 13、算法（摘自OPENCV）... 14、基础数学知识.... 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 ¢er, 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 ¢er, 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 ¢er, 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);    radius = (float)norm(Point2f(dx, dy)) / 2.0f + EPS;     for (int k = 0; k < j; ++k)    {        dx = center.x - (float)pts[k].x;        dy = center.y - (float)pts[k].y;        if (norm(Point2f(dx, dy)) < radius)        {            continue;        }        else        {            Point2f ptsf[3];            ptsf[0] = (Point2f)pts[i];            ptsf[1] = (Point2f)pts[j];            ptsf[2] = (Point2f)pts[k];            findEnclosingCircle3pts_orLess_32f(ptsf,3, center, radius);        }    }}  template<typename PT>void findSecondPoint(const PT *pts, int i, Point2f ¢er, float &radius){    center.x = (float)(pts[0].x + pts[i].x) / 2.0f;    center.y = (float)(pts[0].y + pts[i].y) / 2.0f;    float dx = (float)(pts[0].x - pts[i].x);    float dy = (float)(pts[0].y - pts[i].y);    radius = (float)norm(Point2f(dx, dy)) / 2.0f + EPS;     for (int j = 1; j < i; ++j)    {        dx = center.x - (float)pts[j].x;        dy = center.y - (float)pts[j].y;        if (norm(Point2f(dx, dy)) < radius)        {            continue;        }        else        {            findThirdPoint(pts, i, j, center, radius);        }    }}  template<typename PT>static voidfindMinEnclosingCircle(const PT *pts, int count, Point2f ¢er, float &radius){    center.x = (float)(pts[0].x + pts[1].x) / 2.0f;    center.y = (float)(pts[0].y + pts[1].y) / 2.0f;    float dx = (float)(pts[0].x - pts[1].x);    float dy = (float)(pts[0].y - pts[1].y);    radius = (float)norm(Point2f(dx, dy)) / 2.0f + EPS;     for (int i = 2; i < count; ++i)    {        dx = (float)pts[i].x - center.x;        dy = (float)pts[i].y - center.y;        float d = (float)norm(Point2f(dx, dy));        if (d < radius)        {            continue;        }        else        {            findSecondPoint(pts, i, center, radius);        }    }} //see Welzl, Emo. Smallest enclosing disks (balls and ellipsoids). SpringerBerlin Heidelberg, 1991.voidminEnclosingCircle_jbl(InputArray _points, Point2f& _center, float& _radius){    Mat points = _points.getMat();    int count = points.checkVector(2);    int depth = points.depth();    Point2f center;    float radius = 0.f;    CV_Assert(count >= 0 && (depth == CV_32F || depth == CV_32S));     _center.x = _center.y = 0.f;    _radius = 0.f;     if (count == 0)        return;     bool is_float = depth == CV_32F;    const Point* ptsi = points.ptr<Point>();    const Point2f* ptsf = points.ptr<Point2f>();     // point count <= 3    if (count <= 3)    {        Point2f ptsf3[3];        for (int i = 0; i < count; ++i)        {            ptsf3[i] = (is_float) ? ptsf[i] : Point2f((float)ptsi[i].x, (float)ptsi[i].y);        }        findEnclosingCircle3pts_orLess_32f(ptsf3,count, center, radius);        _center = center;        _radius = radius;        return;    }     if (is_float)    {        findMinEnclosingCircle<Point2f>(ptsf, count,center, radius);#if 0        for (size_t m = 0; m < count; ++m)        {            float d = (float)norm(ptsf[m] - center);            if (d > radius)            {                printf("error!\n");            }        }#endif    }    else    {        findMinEnclosingCircle<Point>(ptsi, count,center, radius);#if 0        for (size_t m = 0; m < count; ++m)        {            double dx = ptsi[m].x - center.x;            double dy = ptsi[m].y - center.y;            double d = std::sqrt(dx * dx + dy * dy);            if (d > radius)            {                printf("error!\n");            }        }#endif    }    _center = center;    _radius = radius;}  

4、基础数学知识

三角形的外心三角形的外心是三角形三边的垂直平分线的交点(或三角形外接圆的圆心) 。三角形的三条垂直平分线必交于一点三角形的三条垂直平分线必交于一点已知：△ABC中，AB,AC的垂直平分线DO,EO相交于点O求证：O点在BC的垂直平分线上证明：连结AO,BO,CO，∵DO垂直平分AB，∴AO=BO∵EO垂直平分AC，∴AO=CO∴BO=CO即O点在BC的垂直平分线上三角形的外心的性质1.三角形三条边的垂直平分线的交于一点，该点即为三角形外接圆的圆心.2三角形的外接圆有且只有一个，即对于给定的三角形，其外心是唯一的，但一个圆的内接三角形却有无数个，这些三角形的外心重合。3.锐角三角形的外心在三角形内；钝角三角形的外心在三角形外；直角三角形的外心与斜边的中点重合4.OA=OB=OC=R5.∠BOC=2∠BAC，∠AOB=2∠ACB，∠COA=2∠CBA6.S△ABC=abc/4R三角形的内心三角形的内心是三角形三条内角平分线的交点(或内切圆的圆心)。证明三角形的三条角平分线必交于一点己知：在△ABC中，∠A与∠B的角平分线交于点O，连接OC求证：OC平分∠ACB三角形内心证明：过O点作OD,OE,OF分别垂直于AC,BC,AB,垂足分别为D,E,F∵AO平分∠BAC,∴OD=OF；∵BO平分∠ABC,∴OE=OF ；∴OD=OF∴O在∠ACB角平分线上∴CO平分∠ACB性质1.三角形的三条角平分线交于一点，该点即为三角形的内心2.三角形的内心到三边的距离相等，都等于内切圆半径r3.r=2S/(a+b+c)4.在Rt△ABC中，∠C=90°，r=(a+b-c)/2．5.∠BOC = 90 °+∠A/2 ∠BOA = 90 °+∠C/2 ∠AOC = 90 °+∠B/26.S△=[(a+b+c)r]/2 (r是内切圆半径)三角形的垂心三角形的垂心是三角形三边上的高的交点(通常用H表示)。三角形的三条高必交于一点已知：△ABC中，AD、BE是两条高，AD、BE交于点O，连接CO并延长交AB于点F三角形的三条高必交于一点求证：CF⊥AB证明：连接DE ∵∠ADB=∠AEB=90°,且在AB同旁，∴A、B、D、E四点共圆 ∴∠ADE=∠ABE （同弧上的圆周角相等）∵∠EAO=∠DAC ∠AEO=∠ADC =90°∴△AEO∽△ADC ∴AE/AD=AO/AC 即AE/AO=AD/AC∴ΔEAD∽ΔOAC ∴∠ACF=∠ADE=∠ABE又∵∠ABE+∠BAC=90° ∴∠ACF+∠BAC=90° ∴CF⊥AB三角形的垂心的性质1.锐角三角形的垂心在三角形内；直角三角形的垂心在直角顶点上；钝角三角形的垂心在三角形外2.三角形的垂心是它垂足三角形的内心；或者说，三角形的内心是它旁心三角形的垂心3. 垂心O关于三边的对称点，均在△ABC的外接圆圆上。4.△ABC中，有六组四点共圆，有三组(每组四个)相似的直角三角形，且AO·OD=BO·OE=CO·OF5. H、A、B、C四点中任一点是其余三点为顶点的三角形的垂心(并称这样的四点为一—垂心组)。6.△ABC，△ABO，△BCO，△ACO的外接圆是等圆。7.在非直角三角形中，过O的直线交AB、AC所在直线分别于P、Q，则 AB/AP·tanB+AC/AQ·tanC=tanA+tanB+tanC8.三角形任一顶点到垂心的距离，等于外心到对边的距离的2倍。9.设O，H分别为△ABC的外心和垂心，则∠BAO=∠HAC，∠ABH=∠OBC，∠BCO=∠HCA。10.锐角三角形的垂心到三顶点的距离之和等于其内切圆与外接圆半径之和的2倍。11.锐角三角形的垂心是垂足三角形的内心；锐角三角形的内接三角形(顶点在原三角形的边上)中，以垂足三角形的周长最短。（施瓦尔兹三角形，最早在古希腊时期由海伦发现）12.西姆松(Simson)定理（西姆松线）：从一点向三角形的三边所引垂线的垂足共线的重要条件是该点落在三角形的外接圆上13.设锐角△ABC内有一点P，那么P是垂心的充分必要条件是PB*PC*BC+PB*PA*AB+PA*PC*AC=AB*BC*CA。14.设H为非直角三角形的垂心，且D、E、F分别为H在BC，CA，AB上的射影，H1，H2，H3分别为△AEF，△BDF，△CDE的垂心，则△DEF≌△H1H2H3。15.三角形垂心H的垂足三角形的三边，分别平行于原三角形外接圆在各顶点的切线。三角形的重心三角形的重心是三角形三条中线的交点。三角形的三条中线必交于一点已知：△ABC的两条中线AD、CF相交于点O，连结并延长BO，交AC于点E。三角形的三条中线必交于一点求证：AE=CE证明：延长OE到点G，使OG=OB∵OG=OB,∴点O是BG的中点又∵点D是BC的中点∴OD是△BGC的一条中位线∴AD∥CG∵点O是BG的中点，点F是AB的中点∴OF是△BGA的一条中位线　∴CF∥AG∵AD∥CG，CF∥AG,∴四边形AOCG是平行四边形　∴AC、OG互相平分,∴AE=CE三角形的重心的性质1.重心到顶点的距离与重心到对边中点的距离之比为2：1。2.重心和三角形3个顶点组成的3个三角形面积相等。3.重心到三角形3个顶点距离的平方和最小。4.在平面直角坐标系中，重心的坐标是顶点坐标的算术平均，即其坐标为((X1+X2+X3)/3,(Y1+Y2+Y3)/3)；空间直角坐标系——横坐标：(X1+X2+X3)/3 纵坐标：(Y1+Y2+Y3)/3 竖坐标：（Z1+Z2+Z3）/35.重心和三角形3个顶点的连线的任意一条连线将三角形面积平分。6.重心是三角形内到三边距离之积最大的点。三角形的旁心介绍三角形的一条内角平分线与另两个内角的外角平分线相交于一点，是旁切圆的圆心，称为旁心。旁心常常与内心联系在一起，旁心还与三角形的半周长关系密切，三角形有三个旁心。三角形旁心的性质1、三角形一内角平分线和另外两顶点处的外角平分线交于一点，该点即为三角形的旁心。2、每个三角形都有三个旁心。 3、旁心到三边的距离相等。如图，点A就是△BCD的一个旁心。三角形任意两角的外角平分线和第三个角的内角平分线的交点。一个三角形有三个旁心，而且一定在三角形外。欧拉线非等边三角形的外心、重心、垂心，依次位于同一直线上，这条直线就叫三角形的欧拉线。其中，重心到外心的距离是重心到垂心距离的一半。欧拉线的证法1作△ABC的外接圆，连结并延长BO，交外接圆于点D。连结AD、CD、AH、CH、OH。作中线AM，设AM交OH于点G’∵ BD是直径∴∠BAD、∠BCD是直角∴ AD⊥AB，DC⊥BC∵ CH⊥AB，AH⊥BC∴ DA‖CH，DC‖AH∴四边形ADCH是平行四边形∴ AH=DC∵ M是BC的中点，O是BD的中点∴ OM= 1/2DC∴ OM= 1/2AH∵ OM‖AH∴△OMG’ ∽△HAG’∴AG/GM=2/1∴ G’是△ABC的重心∴ G与G’重合∴ O、G、H三点在同一条直线上如果使用向量,证明过程可以极大的简化,运用向量中的坐标法,分别求出O G H三点的坐标即可.欧拉线的证法2设H,G,O,分别为△ABC的垂心、重心、外心。连接AG并延长交BC于D, 则可知D为BC中点。连接OD ，又因为O为外心，所以OD⊥BC。连接AH并延长交BC于E,因H为垂心,所以 AE⊥BC。所以OD//AE，有∠ODA=∠EAD。由于G为重心，则GA:GD=2:1。连接CG并延长交BA于F,则可知F为AB中点。同理，OF//CM.所以有∠OFC=∠MCF连接FD，有FD平行AC,且有DF:AC=1:2。FD平行AC，所以∠DFC=∠FCA，∠FDA=∠CAD，又∠OFC=∠MCF，∠ODA=∠EAD，相减可得∠OFD=∠HCA,∠ODF=∠EAC,所以有△OFD∽△HCA,所以OD:HA=DF:AC=1:2；又GA:GD=2:1所以OD:HA=GA:GD=2:1又∠ODA=∠EAD，所以△OGD∽△HGA。所以∠OGD=∠AGH,又连接AG并延长，所以∠AGH+∠DGH=180°,所以∠OGD+∠DGH=180°。即O、G、H三点共线。欧拉线的证法3设H,G,O,分别为△ABC的垂心、重心、外心.则向量OH=向量OA+向量+OB+向量OC向量OG=（向量OA+向量OB+向量OC）/3，向量OG*3=向量OH所以O、G、H三点共线转：http://blog.csdn.net/Cracent/article/details/51524354