最小面积圆覆盖

MinimumCircleCover.h

#pragma once
class CMinimumCircleCover
{
public:
CMinimumCircleCover();
~CMinimumCircleCover();

static double getMinimumCircle(AcGePoint3dArray& allpoints, AcGePoint3d& circlecenterpoint);
static AcGePoint3dArray testfilterPoints(AcGePoint3dArray& allpoints, AcGePoint3d& circlecenterpoint,double radiusdis);
};

MinimumCircleCover.cpp

#include "stdafx.h"
#include "MinimumCircleCover.h"

#include<stdio.h>
#include<math.h>
#include<string.h>

#define MAXN 1000
struct pointset {
double x, y;
};
const double precison = 1.0e-8;
pointset maxcic, point[MAXN];
double radius;
int curset[MAXN], posset[3];
int set_cnt, pos_cnt;
inline double dis_2(pointset &from, pointset& to) {
return ((from.x - to.x)*(from.x - to.x) + (from.y - to.y)*(from.y - to.y));
}
int in_cic(int pt) {
if (sqrt(dis_2(maxcic, point[pt])) < radius + precison) return 1;
return 0;
}
int cal_mincic() {
if (pos_cnt == 1 || pos_cnt == 0)
return 0;
else if (pos_cnt == 3) {
double A1, B1, C1, A2, B2, C2;
int t0 = posset[0], t1 = posset[1], t2 = posset[2];
A1 = 2 * (point[t1].x - point[t0].x);
B1 = 2 * (point[t1].y - point[t0].y);
C1 = point[t1].x*point[t1].x - point[t0].x*point[t0].x +
point[t1].y*point[t1].y - point[t0].y*point[t0].y;
A2 = 2 * (point[t2].x - point[t0].x);
B2 = 2 * (point[t2].y - point[t0].y);
C2 = point[t2].x*point[t2].x - point[t0].x*point[t0].x +
point[t2].y*point[t2].y - point[t0].y*point[t0].y;
maxcic.y = (C1*A2 - C2*A1) / (A2*B1 - A1*B2);
maxcic.x = (C1*B2 - C2*B1) / (A1*B2 - A2*B1);
radius = sqrt(dis_2(maxcic, point[t0]));
}
else if (pos_cnt == 2) {
maxcic.x = (point[posset[0]].x + point[posset[1]].x) / 2;
maxcic.y = (point[posset[0]].y + point[posset[1]].y) / 2;
radius = sqrt(dis_2(point[posset[0]], point[posset[1]])) / 2;
}
return 1;
}
int mindisk() {
if (set_cnt == 0 || pos_cnt == 3) {
return cal_mincic();
}
int tt = curset[--set_cnt];
int res = mindisk();
set_cnt++;
if (!res || !in_cic(tt)) {
set_cnt--;
posset[pos_cnt++] = curset[set_cnt];
res = mindisk();
pos_cnt--;
curset[set_cnt++] = curset[0];
curset[0] = tt;
}
return res;
}

CMinimumCircleCover::CMinimumCircleCover()
{
}

CMinimumCircleCover::~CMinimumCircleCover()
{
}

double CMinimumCircleCover::getMinimumCircle(AcGePoint3dArray& allpoints, AcGePoint3d& circlecenterpoint)
{
int n = allpoints.length();

if (n == 0)
return 0;
int i;
for (i = 0; i < n; i++)
{
point[i].x = allpoints[i].x;
point[i].y = allpoints[i].y;
}

if (n == 1) {
maxcic.x = point[0].x;
maxcic.y = point[0].y;
radius = 0;
// acutPrintf(_T("Circle %.2lf,%.2lf,0 %.2lf "), maxcic.x, maxcic.y, radius);
circlecenterpoint.x = maxcic.x;
circlecenterpoint.y = maxcic.y;
circlecenterpoint.z = 0;
return 0;
}
set_cnt = n; pos_cnt = 0;
for (i = 0; i < n; i++)
curset[i] = i;
mindisk();
// acutPrintf(_T("Circle %.2lf,%.2lf,0 %.2lf "), maxcic.x, maxcic.y, radius);
circlecenterpoint.x = maxcic.x;
circlecenterpoint.y = maxcic.y;
circlecenterpoint.z = 0;
return radius;
}

AcGePoint3dArray CMinimumCircleCover::testfilterPoints(AcGePoint3dArray& allpoints, AcGePoint3d& circlecenterpoint, double radiusdis)
{
AcGePoint3dArray pts;

if (radiusdis == 0 && allpoints.length() == 1)
{
pts.append(allpoints[0]);
return pts;
}

AcDbCircle *pCircle = new AcDbCircle(circlecenterpoint, AcGeVector3d(0, 0, 1), radiusdis);
if (pCircle == NULL)
{
return pts;
}


for (int i = 0; i < allpoints.length(); i++)
{
double dis = allpoints[i].distanceTo(circlecenterpoint);
if (abs(dis - radiusdis) <= 1)
{
pts.append(allpoints[i]);
allpoints.remove(allpoints[i]);
i--;
}
}

if (pCircle)
pCircle->close();

return pts;

}