C++:给定一个二维点集，找到所有的整体对称轴

这里我列出了两种方法：

 

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);
}

该代码均为笔者手写，如果有不懂可以提问，也能直接拿去使用。