这里我列出了两种方法:
1.基于质心的解法。(这种方法要求除非你拥有“分数”的数据结构,不然质心的求解会出现误差)
2.基于向量的解法。
该算法的核心思想是:先求出凸包,凸包的对称轴才可能是整体的对称轴,以此减小搜索范围。在这过程中,可以使用各种方法(例如向量和为0等等)来不断减少不可能的情况,最终求出结果。
按理说,质心求解更快,但是苦于有精度误差,所以还是方案2更好。
#include"iostream"
#include"math.h"
#include"vector"
#include"algorithm"
#include "iomanip"
#include"mex.h"
using namespace std;
#define PI 3.1415926
class Point//点类
{
public:
Point()
{
}
Point(double x, double y)//构造函数的重载
{
this->x = x;
this->y = y;
}
Point& operator=(const Point &A)//重载赋值运算符
{
this->x = A.x;
this->y = A.y;
return *this;
}
Point& operator+(const Point &A)//重载加号
{
static Point B;
B.x = this->x + A.x;
B.y = this->y + A.y;
return B;
}
Point& operator/(const double &A)//重载除号
{
static Point B;
B.x = this->x / A;
B.y = this->y / A;
return B;
}
double operator*(const Point &A)//重载向量的点乘
{
return x * A.y - y * A.x;
}
bool operator==(const Point&A)//重载相等比较符号
{
return (this->x == A.x && this->y == A.y);
}
bool operator!=(const Point&A)//重载不相等比较符号
{
return !(this->x == A.x && this->y == A.y);
}
void set(double x, double y)//赋值函数
{
this->x = x;
this->y = y;
}
Point To(Point& A)//返回向量
{
Point OA;
OA.x = A.x - this->x;
OA.y = A.y - this->y;
return OA;
}
double DistanceTo(Point B)//求距离函数
{
return sqrt((this->x - B.x)*(this->x - B.x) + (this->y - B.y)*(this->y - B.y));
}
double thea()//求到原点的极角
{
if (x == 0) return PI / 2;
else if (x > 0)
{
return atan(y / x);
}
else if (y > 0)
return PI + atan(y / x);
else
return atan(y / x) - PI;
}
double theaFrom(Point N)//求到点N的极角
{
double x = this->x - N.x;
double y = this->y - N.y;
if (x == 0) return PI / 2;
else if (x > 0)
{
return atan(y / x);
}
else if (y > 0)
return PI + atan(y / x);
else
return atan(y / x) - PI;
}
double x;
double y;
};
class Line
{
public:
Line() {}
Line(double a, double b, double c)//三值式
{
this->a = a;
this->b = b;
this->c = c;
if (a == -0) a = 0;
if (b == -0) b = 0;
if (c == -0) c = 0;
}
Line(Point A, double k)//点斜式
{
this->a = k;
this->b = -1;
this->c = -k * A.x + A.y;
if (a == -0) a = 0;
if (b == -0) b = 0;
if (c == -0) c = 0;
}
Line(Point A, Point B)//两点式
{
if (A.x == B.x)
{
a = 1;
b = 0;
c = -A.x;
}
else
{
double k = (A.y - B.y) / (A.x - B.x);
a = k;
b = -1;
c = -k * A.x + A.y;
}
if (a == -0) a = 0;
if (b == -0) b = 0;
if (c == -0) c = 0;
OA = A.To(B);//得到特征向量OA,用于辅助判断是否垂直
O = A;
}
Line(Point A, Line L)//点垂式
{
if (L.b == 0)
{
a = 0;
b = 1;
c = -A.y;
}
else
if (L.a == 0)
{
a = 1;
b = 0;
c = -A.x;
}
else
{
double k = L.b / L.a;
a = k;
b = -1;
c = -k * A.x + A.y;
}
if (a == -0) a = 0;
if (b == -0) b = 0;
if (c == -0) c = 0;
}
Line& operator=(const Line &A)
{
this->a = A.a;
this->b = A.b;
this->c = A.c;
return *this;
}
bool operator==(Line &A)
{
if (this->a == A.a&&this->b == A.b &&this->c == A.c)
return true;
else
return false;
}
double DistanceFrom(Point A)//求点到直线的距离
{
double x, y;
if (a*A.x + b * A.y + c == 0)
return 0;
else
if (a == 0)
{
x = A.x;
y = -c / b;
}
else
{
y = (a*a*A.y - a * b*A.x - b * c) / (a*a + b * b);
x = (-c - b * y) / a;
}
return A.DistanceTo(Point(x, y));
}
int PointOnLine(Point A)//判断点与直线的位置关系
{
double d = a * A.x + b * A.y + c;
if (d > 0) return 1;
else if (d < 0) return -1;
else return 0;
}
Point& SymmetryPoint(const Point &A)//对称点
{
static Point A1;
if (a*A.x + b * A.y + c == 0)
A1 = A;
else
if (b == 0)
{
A1.x = 2 * (-c / a) - A.x;//a和b不可能同时为0
A1.y = A.y;
}
else if (a == 0)
{
A1.x = A.x;
A1.y = 2 * (-c / b) - A.y;
}
else
{
A1.y = (-2 * a*b*A.x + (a*a - b * b)*A.y - 2 * b*c) / (a*a + b * b);
A1.x = a / b * (A1.y - A.y) + A.x;
}
return A1;
}
double a = 0, b = 0, c = 0;
Point OA = Point(0, 0);
Point O = Point(0, 0);
};
Point startPoint;//起点
Point centroid;//定义质心
vector<Point>::iterator it;//定义一个点循环指针,专门用于循环
vector<Point> Up;//点集的上半部分(由对称轴确定的)
vector<Point> Down;//点集的下半部分(由对称轴确定的)
vector<Line> Diad;//对称轴的容器
bool on_segment(Point p, Point p1, Point p2)
{
double min_x = p1.x < p2.x ? p1.x : p2.x;
double max_x = p1.x > p2.x ? p1.x : p2.x;
double min_y = p1.y < p2.y ? p1.y : p2.y;
double max_y = p1.y > p2.y ? p1.y : p2.y;
if (p.x >= min_x && p.x <= max_x
&& p.y >= min_y && p.y <= max_y) return true;
else return false;
}
bool sortByPolorAngle(Point & p1, Point & p2)
{
Point op1 = startPoint.To(p1);
Point op2 = startPoint.To(p2);
double d = op1 * op2;//向量点乘,也就是叉积
if (d < 0) return true;//叉积小于0,说明是顺时针旋转了,我们需要
if (d > 0) return false;//叉积大于0,说明是逆时针旋转了,我们不需要
if (d == 0 && on_segment(startPoint, p1, p2))return true;//共线的两种情况,我们需要
if (d == 0 && on_segment(p2, startPoint, p1)) return true;
return false;
}
vector<Point>& GetConvexHull(vector<Point> & point)//求凸壳
{
Point p0 = point[0];//起点
int k = 0;
for (int i = 1; i < point.size(); i++)
{
if (point[i].y < p0.y ||
point[i].y == p0.y && point[i].x < p0.x)
{
p0 = point[i];
k = i;
}
}
point.erase(point.begin() + k);
point.insert(point.begin(), p0);//把起点调到开头
static vector<Point> convex_hull;
do {
convex_hull.push_back(point[0]);
startPoint = point[0];//确定起点
point.erase(point.begin());//消去点集中的第一个点
sort(point.begin(), point.end(), sortByPolorAngle);//按照规则排序
point.push_back(convex_hull[convex_hull.size() - 1]);//把原来的开头压入到现在的结尾
if (point[0] == convex_hull[0]) break;//冒泡到了开头了
} while (1);
for (vector<Point>::iterator it = convex_hull.begin(); it < convex_hull.end(); it++)
//将点集中的凸壳去除,便于之后的“对凸壳内部的点是否符合对称轴”的验证
{
vector<Point>::iterator it1 = point.begin();
while (it1 < point.end())
{
if ((*it) == (*it1))
{
point.erase(it1);
break;
}
else
it1++;
}
}
return convex_hull;
}
bool PointsCanSplit(vector<Point> v, Line L)//判断直线是否可以划分两边的点
{
Up.clear(); Down.clear();
int one = 0, two = 0;
for (vector<Point>::iterator it = v.begin(); it < v.end(); it++)
{
Point OB = L.O.To(*it);
if (L.OA*(OB) > 0)//利用叉积将点集其分为上下两部分
{
Up.push_back(*it);
one++;
}
if (L.OA*(OB) < 0)
{
Down.push_back(*it);
two++;
}
}
if (one == two) //如果上下两部分点数相同,则返回真
return true;
else
return false;
}
bool find(vector<Point> v, Point A, vector<Point>::iterator &it1)
//查找函数,用于在点集中查找需要的点,并通过指针it1返回位置
{
vector<Point>::iterator it;
for (it = v.begin(); it < v.end(); it++)
{
if (*it == A)
{
it1 = it;
return true;
}
}
return false;
}
bool find(vector<Line> L, Line P, vector<Line>::iterator &it1)
{
//查找函数,用于在线集中查找需要的线,并通过指针it1返回位置
vector<Line>::iterator it;
for (it = L.begin(); it < L.end(); it++)
{
if (*it == P)
{
it1 = it;
return true;
}
}
return false;
}
bool IsSymmetric(vector<Point> A, Line L)//直接法判断是否对称(对称点的求解有误差)
{
vector<Point>::iterator it1;
for (it = Up.begin(); it < Up.end(); it++)//对凸壳的上半部分每个点
{
if (!find(Down, L.SymmetryPoint(*it), it1))
{
Point P = L.SymmetryPoint(*it);
return false;
}
}
return true;
}
bool IsSymmetric2(vector<Point> A, Line L)//用“垂线+距离”判断法来判断是否对称。(垂线的求解有误差)
{
vector<Point>::iterator it1;
for (it = Up.begin(); it < Up.end(); it++)//对凸壳的上半部分每个点
{
Line L1((*it), L);//点垂式求出直线
bool flag = 0;
for (it1 = Down.begin(); it1 < Down.end(); it1++)
{
if (L1.PointOnLine(*it1) == 0)
{
if (L.DistanceFrom(*it) == L.DistanceFrom(*it1))
flag = 1;
}
}
if (flag != 1) return false;
}
return true;
}
bool IsSymmetric3(vector<Point> A, Line L)//基于向量的垂直判断+距离
{
vector<Point>::iterator it1;
for (it = Up.begin(); it < Up.end(); it++)//对凸壳的上半部分每个点
{
bool flag = 0;
for (it1 = Down.begin(); it1 < Down.end(); it1++)
{
Line L1(*it, *it1);
if (L1.OA.x*L.OA.x + L1.OA.y*L.OA.y == 0)//向量垂直的判断
{
Point OB = L.O.To(((*it) + (*it1)) / 2);
if (L.OA*OB == 0)//向量平行的判断
flag = 1;
}
}
if (flag != 1)
return false;
}
return true;
}
bool IsBalance(vector<Point> A, Line LX)//判断点集内的点是否均匀
{
vector<Point>::iterator it;
double sum_x = 0, sum_y = 0, x, y;
Line LY(centroid, LX);
for (it = A.begin(); it < A.end(); it++)
{
double x_, y_;
x_ = LX.a*(*it).x + LX.b*(*it).y + LX.c;
y_ = LY.a*(*it).x + LY.b*(*it).y + LY.c;
if (x_ >= 0)
x = LY.DistanceFrom(*it);
else
x = -LY.DistanceFrom(*it);
if (y_ >= 0)
y = LX.DistanceFrom(*it);
else
y = -LX.DistanceFrom(*it);
sum_x = sum_x + x;
sum_y = sum_y + y;
}
if (sum_x == 0 && sum_y == 0)
return true;
else
return false;
}
void DealWithDiad(vector<Point> A, Line L)//对于对称轴的处理
{
vector<Line>::iterator it1;
if (PointsCanSplit(A, L))
if (IsSymmetric3(A, L))
if (!find(Diad, L, it1))
Diad.push_back(L);
}
void Design1(vector<Point> A)//基于质心的求凸壳对称轴方法,若质心的求解有误差,则转入枚举法
{
double sum_x = 0, sum_y = 0;
int N = A.size();
for (int i = 0; i < N; i++)
{
sum_x = sum_x + A[i].x; sum_y = sum_y + A[i].y;
}
centroid.x = sum_x / N; centroid.y = sum_y / N;//求出质心
vector<Point> convex_hull = GetConvexHull(A);//求出凸壳
int Number = convex_hull.size() / 2;//只需要前一半的点数
double *AdjacentEdge = new double[Number + 1];//定义相邻边数组(只需要前半部分)
Point *EdgeMidpoint = new Point[Number + 1];//定义边的中点数组(只需要前半部分)
vector<Line>::iterator it1;
if (convex_hull.size() == 2)//凸壳只有两个点,对称轴肯定是两点连线和它的垂线了
{
Line L1(convex_hull[0], convex_hull[1]);
Line L2((convex_hull[0] + convex_hull[1]) / 2, L1);
Diad.push_back(L1);
Diad.push_back(L2);
return;
}
//有两个点以上的情况:
for (int i = 0; i < Number + 1; i++)//求出前一半的中点,相邻边长
{
EdgeMidpoint[i] = (convex_hull[i] + convex_hull[i + 1]) / 2.0;
AdjacentEdge[i] = convex_hull[i].DistanceTo(convex_hull[i + 1]);
}
//判断相邻边的中点与质心的连线是否为凸壳对称轴
for (int i = 0; i < Number + 1; i++)
{
Line L(EdgeMidpoint[i], centroid);
DealWithDiad(convex_hull, L);
}
//判断相邻边的角平分线是否为凸壳对称轴
for (int i = 0; i < Number + 1; i++)
{
if (AdjacentEdge[i] == AdjacentEdge[i + 1])//当前边与下一条边相等
{
Line L(convex_hull[i + 1], centroid);
DealWithDiad(convex_hull, L);
}
}
}
void Design2(vector<Point> A)//用枚举法求凸壳的对称轴
{
vector<Point> convex_hull = GetConvexHull(A);//求出凸壳
int N = convex_hull.size();
if (N == 2)//凸壳只有两个点,肯定是两点连线和它的垂线了
{
Line L1(convex_hull[0], convex_hull[1]);
Line L2((convex_hull[0] + convex_hull[1]) / 2, L1);
Diad.push_back(L1);
Diad.push_back(L2);
return;
}
Point *EdgeMidpoint = new Point[N];//定义边的中点数组
for (int i = 0; i < N; i++)
{
if (i == convex_hull.size() - 1)
EdgeMidpoint[i] = (convex_hull[i] + convex_hull[0]) / 2.0;
else
EdgeMidpoint[i] = (convex_hull[i] + convex_hull[i + 1]) / 2.0;
}
vector<Point> convex_hull_new;//新的包含中点的凸壳
for (int i = 0; i < N; i++)
{
convex_hull_new.push_back(convex_hull[i]);
}
for (int i = 0; i < N; i++)
{
convex_hull_new.push_back(EdgeMidpoint[i]);
}
N = convex_hull_new.size();
for (int i = 0; i < N; i++)
{
for (int j = i; j < N; j++)
{
if (i == j)
continue;
else
{
Line L(convex_hull_new[i], convex_hull_new[j]);//用两点式求得直线
DealWithDiad(convex_hull, L);
}
}
}
}
/*
int main()
{
vector<Point> A;//定义点集
Point a;//定义一个点
cout << "请输入点的个数:" << endl;
int N;
cin >> N;
cout << "请输入点的横纵坐标:" << endl;
double x, y;
for (int i = 0; i < N; i++)
{
cin >> x;
cin >> y;
a.set(x, y);
A.push_back(a);
}
// Design1(A);//方案一,优化之后的方案
Design2(A);//方案二
if (A.size() != 0)//内部点集不为空则进行验证
{
for (vector<Line>::iterator it = Diad.begin(); it < Diad.end(); it++)//验证内部点
{
if (IsBalance(A, (*it)))
{
if (IsSymmetric3(A, (*it)))
continue;
else
Diad.erase(it);
}
}
}
cout << "对称轴的条数为:" << Diad.size() << endl;
cout << "对称轴为:" << endl;
if (Diad.size() == 0)
cout << "无" << endl;
else
{
for (vector<Line>::iterator it = Diad.begin(); it < Diad.end(); it++)
{
cout << setw(20) << (*it).a << "*x + " << setw(20) << (*it).b << "*y + " << setw(20) << (*it).c << "=0" << endl;
}
}
system("pause");
return 0;
}*/
void Main(int N, double *x, double *y)
{
vector<Point> A;//定义点集
Point a;//定义一个点
for (int i = 0; i < N; i++)
{
a.set(x[i], y[i]);
A.push_back(a);
}
// Design1(A);//方案一,优化之后的方案
Design2(A);//方案二
if (A.size() != 0)//内部点集不为空则进行验证
{
for (vector<Line>::iterator it = Diad.begin(); it < Diad.end(); it++)//验证内部点
{
if (IsBalance(A, (*it)))
{
if (IsSymmetric3(A, (*it)))
continue;
else
Diad.erase(it);
}
}
}
cout << "对称轴的条数为:" << Diad.size() << endl;
cout << "对称轴为:" << endl;
if (Diad.empty())
cout << "无" << endl;
else
{
for (auto it = Diad.begin(); it < Diad.end(); ++it)
{
cout << setw(20) << (*it).a << "*x + " << setw(20) << (*it).b << "*y + " << setw(20) << (*it).c << "=0" << endl;
}
}
}
两种方法均列在这里面。去掉注释打开main即可使用。这其中Main函数就是执行函数,目前选中的执行方案是方案一,方案二大家可以自己改良。
这里给出了Matlab的兼容性代码,将这个和上面的代码放在一起,在matlab中编译,即可运行。
void mexFunction(int nlhs, mxArray*plhs[], int nrhs, const mxArray*prhs[])
{
int N;
double *inMatrix1;
double *inMatrix2;
if (nrhs != 4)
{
mexErrMsgIdAndTxt("MyToolbox:arrayProduct:nrhs", "There inputs required.");
}
if (nlhs != 0)
{
mexErrMsgIdAndTxt("MyToolbox:arrayProduct:nlhs", "Zero output requird.");
}
if (!mxIsDouble(prhs[0]) || mxIsComplex(prhs[10]) || mxGetNumberOfElements(prhs[0]) != 1)
{
mexErrMsgIdAndTxt("MyToolbox:arrayProduct:notScalar", "Input multiplier must be a scalar.");
}
if (mxGetM(prhs[1]) != 1)
{
mexErrMsgIdAndTxt("MyToolbox:arrayProduct:notRowVector", "Input must be a row vector.");
}
if (mxGetM(prhs[2]) != 1)
{
mexErrMsgIdAndTxt("MyToolbox:arrayProduct:notRowVector", "Input must be a row vector.");
}
N = mxGetScalar(prhs[0]);
inMatrix1 = mxGetPr(prhs[1]);
inMatrix2 = mxGetPr(prhs[2]);
Main(N, inMatrix1, inMatrix2);
}
该代码均为笔者手写,如果有不懂可以提问,也能直接拿去使用。