一 算法思想
网上找内切圆算法没找到理想的，在国外看到一篇文章，作者本人也有源码，编译没通过，自己重新写了下，实现了功能
上图中在三角形形成的二维平面区域中，高度为H，宽度为W。将空间进行等分，高度为n份，宽度为m份，每份为H/n,W/m份
在空间形成网格，将三角区域外的点剔除掉。再剩下的点中，对每个点而言，找出到边最近的距离。这样在很多点中找出最大的距离点。这个点作为种子点。以种子点为中心的区域再剖分，循环迭代
二.代码部分
#define N_CELLS 20
#define M_CELLS 20
int FindInscribedCircleCenter(vector<cv::Point3f> polygon, double* output)
{
if (polygon.size() <= 1)
return 0;
float bounds[4] = {10000,-10000,100000,-10000};
for (int i = 0; i < polygon.size(); i++)
{
cv::Point3f 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;
}
cv::Point3f point_pia;
float flt_tmp = FLT_MAX;
int count = 1;
cv::Point3f 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;
}
找到内切圆函数
cv::Point3f GeometryFindPIA(vector<cv::Point3f> polygon, float bounds[])
{
cv::Point3f pia;
pia.x = (bounds[0] + bounds[1]) / 2;
pia.y = (bounds[2] + bounds[3]) / 2;
cv::Point3f tmp;
float increment_x = (bounds[1] - bounds[0]) / N_CELLS;
float increment_y = (bounds[3] - bounds[2]) / M_CELLS;
float max_distance = 0;
int i, j;
float 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<cv::Point3f>& polyline, cv::Point3f& pt)
{
int count = 0;
cv::Point3f mark1(1000, 1000, 0);
for (int i = 0; i < polyline.size(); i++)
{
int id = (i + 1) % polyline.size();
if (get_line_intersection(polyline[i], polyline[id], pt, mark1) == 1)
{
count++;
}
}
if (count % 2 == 1)
return true;
else
return false;
}
判断两个线段有没有交点 引用了“Unity 数学知识——数学公式汇总”博客中的部分内容
int get_line_intersection(cv::Point3f& p0, cv::Point3f& p1, cv::Point3f& p2, cv::Point3f& p3)
{
float s_numer, t_numer, denom, t;
cv::Point3f s10, s32;
s10 = p1 - p0;
s32 = p3 - p2;
denom = s10.x * s32.y - s32.x * s10.y;
if (denom == 0)//平行或共线
return 0; // Collinear
bool denomPositive = denom > 0;
cv::Point3f s02;
s02 = p0 - p2;
s_numer = s10.x * s02.y - s10.y * s02.x;
if ((s_numer < 0) == denomPositive)//参数是大于等于0且小于等于1的，分子分母必须同号且分子小于等于分母
return 0; // No collision
t_numer = s32.x * s02.y - s32.y * s02.x;
if ((t_numer < 0) == denomPositive)
return 0; // No collision
if (fabs(s_numer) > fabs(denom) || fabs(t_numer) > fabs(denom))
return 0; // No collision
return 1;
}
点到多边形距离
double DistancePointAPolygon(vector<cv::Point3f>& polyline, cv::Point3f& pt)
{
///1 先判断垂足位置
cv::Point3f mark;
double val,val1,val2,val3;
double max_dist = 10000;
for (int i = 0; i < polyline.size(); i++)
{
int id = (i + 1) % polyline.size();
cv::Point3f 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 (get_line_intersection(polyline[i], polyline[id], pt, mark) == 1) // 有交点
{
cv::Point3f 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
{
cv::Point3f dist1 = pt - polyline[i];
cv::Point3f 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;
}
三 效果说明
要求输入的点为有序点（顺时针），代码是c++版，的借助于opencv的Point3f 来表示点。效果如下算法论文
1.“An Efficient Algorithm to Calculate the Center of the Biggest Inscribed Circle in an Irregular Polygon”