任意多边形的最大内切圆算法

admin1

[code]//任意多边形的最大内切圆算法https://blog.csdn.net/u011533238/article/details/88844092[/code]
//输入的点为有序点（顺时针）
#define N_CELLS 20
#define M_CELLS 20


int  FindInscribedCircleCenter(vector<AcGePoint3d> polygon, double* output)
{

         if (polygon.size() <= 1)
                 return 0;
         double bounds[4] = {10000,-10000,100000,-10000};
         for (int i = 0; i < polygon.size(); i++)
         {
                 AcGePoint3d pt = polygon[i];
                 if (pt.x < bounds[0]) bounds[0] = pt.x;
                 if (pt.x > bounds[1]) bounds[1] = pt.x;
                 if (pt.y < bounds[2]) bounds[2] = pt.y;
                 if (pt.y > bounds[3]) bounds[3] = pt.y;
         }

         AcGePoint3d point_pia;

         double flt_tmp = FLT_MAX;
         int count = 1;
         AcGePoint3d point_tmp;// = new cv::Point3f;
         while (count++)
         {

                 // find new candidate PIA

                 point_tmp = GeometryFindPIA(polygon, bounds);


                 // update current PIA
                 point_pia.x = point_tmp.x;
                 point_pia.y = point_tmp.y;

                 // update the bounds
                 flt_tmp = (bounds[1] - bounds[0]) / (sqrtf(2) * 2);
                 bounds[0] = point_pia.x - flt_tmp;
                 bounds[1] = point_pia.x + flt_tmp;
                 flt_tmp = (bounds[3] - bounds[2]) / (sqrtf(2) * 2);
                 bounds[2] = point_pia.y - flt_tmp;
                 bounds[3] = point_pia.y + flt_tmp;


                 if (bounds[1] - bounds[0] < 0.001 || bounds[3] - bounds[2] < 0.001) break;

                 //        printf("Candidate PIA: (%f,%f)\n", point_pia.x, point_pia.y);
         }


         double tmp_distance = DistancePointAPolygon(polygon, point_pia);
         output[0] = point_pia.x;
         output[1] = point_pia.y;
         output[2] = polygon[0].z;
         output[3] = tmp_distance;
         return 1;

}


//找到内切圆函数
AcGePoint3d GeometryFindPIA(vector<AcGePoint3d> polygon, double bounds[])
{
         AcGePoint3d pia;

         pia.x = (bounds[0] + bounds[1]) / 2;
         pia.y = (bounds[2] + bounds[3]) / 2;
         AcGePoint3d tmp;


         double increment_x = (bounds[1] - bounds[0]) / N_CELLS;
         double increment_y = (bounds[3] - bounds[2]) / M_CELLS;


         double max_distance = 0;


         int i, j;
         double tmp_distance = FLT_MAX;
         for (i = 0; i <= N_CELLS; i++)
         {
                 tmp.x = bounds[0] + i * increment_x;

                 for (j = 0; j <= M_CELLS; j++)
                 {
                         tmp.y = bounds[2] + j * increment_y;

                         // compare with candidate PIA if point is in polygon
                         if (IsPointInPolygon(polygon, tmp))
                         {
                                 tmp_distance = DistancePointAPolygon(polygon, tmp);

                                 if (tmp_distance > max_distance)
                                 {
                                         max_distance = tmp_distance;
                                         pia.x = tmp.x;
                                         pia.y = tmp.y;
                                 }
                         }
                 }
         }
         return pia;
}





//判断点在多边形内部的函数 ，算法思路是这个点和多边形的边去求交点如果有奇数个交点为多边形内部。
bool IsPointInPolygon(vector<AcGePoint3d>& polyline, AcGePoint3d& pt)
{
         int count = 0;
         AcGePoint3d mark1(1000, 1000, 0);
         for (int i = 0; i < polyline.size(); i++)
         {
                 int id = (i + 1) % polyline.size();
                 if (Four_point_intersection(polyline[i], polyline[id], pt, mark1) == 1)
                 {
                         count++;
                 }

         }
         if (count % 2 == 1)
                 return true;
         else
                 return false;

}



//四点交点 判断两个线段有没有交点
static  int Four_point_intersection(AcGePoint3d& p0, AcGePoint3d& p1, AcGePoint3d& p2, AcGePoint3d& p3)
{

         double  s_numer, t_numer, denom;
         AcGeVector3d s10 = p1 - p0;
         AcGeVector3d s32 = p3 - p2;


         denom = s10.x * s32.y - s32.x * s10.y;
         if (denom == 0)//平行或共线
                 return 0; // Collinear
         bool denomPositive = denom > 0;

         AcGeVector3d s02 = p0 - p2;

         s_numer = s10.x * s02.y - s10.y * s02.x;
         if ((s_numer < 0) == denomPositive)//参数是大于等于0且小于等于1的，分子分母必须同号且分子小于等于分母
                 return 0; // 无碰撞


         t_numer = s32.x * s02.y - s32.y * s02.x;
         if ((t_numer < 0) == denomPositive)
                 return 0; // 无碰撞

         if (fabs(s_numer) > fabs(denom) || fabs(t_numer) > fabs(denom))
                 return 0; // 无碰撞
         return 1;
}



//点到多边形距离
static double DistancePointAPolygon(vector<AcGePoint3d>& polyline, AcGePoint3d& pt)
{
         ///1 先判断垂足位置
         AcGePoint3d  mark;
         double val,val1,val2;
         double max_dist = 10000;
         for (int i = 0; i < polyline.size(); i++)
         {
                 int id = (i + 1) % polyline.size();

                 AcGeVector3d diff = polyline[id] - polyline[i];  ///  -b a

                 mark.x = -diff.y * 10000 + pt.x;
                 mark.y = diff.x * 10000 + pt.y;
                 double dist = 0;

                 if (Four_point_intersection(polyline[i], polyline[id], pt, mark) == 1)  // 有交点
                 {
                         AcGeVector3d  pt_edge = pt - polyline[i];
                         val = pt_edge.x*diff.x + pt_edge.y*diff.y;
                         val1 = pt_edge.x*pt_edge.x + pt_edge.y*pt_edge.y;
                         if (val1 < 0.00001)
                                 dist = 0;
                         else
                         {
                                 val1 = sqrt(val1);
                                 val2 = diff.x*diff.x + diff.y*diff.y;
                                 val2 = sqrt(val2);
                                 double temp = val / (val1*val2);
                                 if (temp > 1) temp = 1.0;
                                 if (temp < -1) temp = -1.0;

                                 dist = val1*sin(acos(temp));
                                 // double tmep_ = val / (val2*val1);
                                 // dist = sqrt(1 - temp*temp)*val1;
                         }
                 }
                 else
                 {
                         AcGeVector3d dist1 = pt - polyline[i];
                         AcGeVector3d dist2 = pt - polyline[id];
                         val1 = dist1.x*dist1.x + dist1.y*dist1.y;
                         val2 = dist2.x*dist2.x + dist2.y*dist2.y;
                         if (val1 > val2)
                         {
                                 val2 = sqrt(val2);
                                 dist = val2;
                         }
                         else
                         {
                                 val1 = sqrt(val1);
                                 dist = val1;
                         }

                 }
                 if (dist < max_dist)
                         max_dist = dist;
         }
         return max_dist;
}