【全网最详细】如何实现多边形直角化，判断多边形顺时针or逆时针，凹边形各角度数怎么求（含全部代码），如何为vs工程引入ceres-solve库

概要

看完这篇文章，你将学会：

• 多边形的顺时针逆时针怎么判断
• 在凹边形中如何求得真实的角度数（两边互相垂直的时候，如在多边形中如何知道是90°还是270°）
• 实操实现多边形直角化（如何引入ceres库，实现最小二乘法，为自己的工厂引入非线性优化功能）

所有代码都在最下方，看完之后如果对大家有帮助，大家也可以评论留言，希望不再只是机器人点赞了。

多边形的顺时针逆时针怎么判断

代码：参考AngleOptimizationCostFunction::AngleOptimizationCostFunction(const std::vectorEigen::Vector3d& points) 实现。
原理：封闭图形的顺逆时针方向，有使用有向面积的，OpenGL里也有（我没用这种）。我参考的是极值法，也就是寻找一个极值点，能找到过这个点的一条线，让所有其他点都在这条线的另一边。代码中我直接找某一个分量（x或者y）最大的点。定义相邻边叉乘（不懂就按照我下个问题解答的定义方式），z小于0为逆时针，否则为顺时针。

在凹边形中如何求得真实的角度数

代码：参考AngleOptimizationCostFunction::ComputeAngle(const Eigen::Matrix<T, 3, 1>& v1, const Eigen::Matrix<T, 3, 1>& v2)实现。
原理：如上图，经过观察可以发现，从垂足A为起点，沿着两条垂边向外的两条向量所构成的直角∠BAC，当ABxAC，也就是“前边”（点序号组成更小的边）叉乘“后边”，根据右手螺旋定理，结果朝内（注意是沿着黑色的弧线绕！）。同时，可以观察到整个封闭图形的是逆时针的，也根据右手螺旋定责，朝向向外。多试几组不同的情况后，可以得出结论：当相邻两边的叉乘结果（按照我这种定义，有的算法往往不强调，导致大家BAxAC，什么的都有，这样算出的结果也会出错）和封闭图形的方向相反的时候，这个角度为270°，否则为90°。

解决多边形直角化的问题

代码：参考void BecomeRectangular(std::vectorEigen::Vector3d points, std::vectorEigen::Vector3d & newPnts)实现。
原理：最小二乘法，大家都或多或少听过，大多数印象是在解决线性问题，求最优拟合的时候运用。但实际情况，我们也可以去设置一些预期条件，通过最小二乘找到最优解。解决多边形直角化的关键就是如何设置这些直角化的预期条件，首先我们找到多边形中接近90°或者270°的角，我设置了一个取值范围，75-105和255-285，这些角度我期望都变成90度。然后引入residuals数组，存储每次迭代每个期望垂直的角和实际垂直角的差值，这里还需要额外引入一个多边形内角和，亲试否则求解也会有问题。使用的时候请使用ceres::DynamicAutoDiffCostFunction，这个是动态求解的函数，可支持任意边的多边形。其他的就看代码备注就好。

如何引入ceres库，实现最小二乘法，为自己的工厂引入非线性优化功能

这个放在最后说，我知道对于很多人来说，引入第三方库有时候也是件特别麻烦的事。这里我把我遇到的所有经过和坑都讲一遍，希望有帮助。
首先登入ceres-solve官网，如果你对Cmake很熟练，就想自己尝试各种配置设置，那直接自己编译就好。我找到了官网推荐的专门为Windows生成好的解决方案，如图。https://github.com/tbennun/ceres-windows
找到sln，我用的2017，直接打开ceres-2015.sln，重定向使用140以上版本打开。根据自己工程使用动态库还是静态库的需求，编译对应的库文件即可。
把对应sln同级别 glog\src\windows\下的glog文件夹拷贝到自己工程中；
把对应sln同级别 ceres-solver\include\下的ceres文件夹拷贝到自己工程中；
把对应sln同级别 Eigen\下的Eigen文件夹拷贝到自己工程中；
然后把sln同级别 win\include\ceres\internal下的config.h文件考到刚刚复制到自己工程的ceres\internal文件夹里；
包含libglog_static.lib和ceres.lib两个库文件。
至此，大功告成。

全部代码

// 头文件
#include <Eigen/Dense> 
#include <cmath>
#include <algorithm>
// 定义优化问题的成本函数
class AngleOptimizationCostFunction {

public:
	AngleOptimizationCostFunction(const std::vector<Eigen::Vector3d>& points);

	template <typename T>
	bool operator()(T const* const* parameters, T* residuals) const;

	template <typename T>
	T ComputeAngle(const Eigen::Matrix<T, 3, 1>& v1, const Eigen::Matrix<T, 3, 1>& v2)const;

private:  
	std::vector<Eigen::Vector3d> points_;  	

	bool geoDir_ = false; 

};

// cpp
#include "ceres/ceres.h"
#include "ceres/rotation.h"

AngleOptimizationCostFunction::AngleOptimizationCostFunction(const std::vector<Eigen::Vector3d>& points) 
	: points_(points) {

	//使用极值法
	double maxX = points[0].x();
	int tatNum = 0;
	int num_points = points.size();
	for (int i = 0; i < num_points; i++) {
		if (maxX <= points[i].x()) {
			tatNum = i;
			maxX = points[i].x();
		}
	}
	Eigen::Vector3d v1 = points[(tatNum - 1 + num_points) % num_points] - points[tatNum];
	Eigen::Vector3d v2 = points[(tatNum + 1) % num_points] - points[tatNum];

	if (v1.cross(v2).z() < 0.0) {
		geoDir_ = false;
	}
	else {
		geoDir_ = true; //顺时针朝外
	}

}


template<typename T>
bool AngleOptimizationCostFunction::operator()(T const * const * parameters, T * residuals) const {
	size_t num_points = points_.size();
	std::vector<Eigen::Matrix<T, 3, 1>> points_3d(num_points);

	for (size_t i = 0; i < num_points; ++i) {
		points_3d[i] = Eigen::Matrix<T, 3, 1>(parameters[0][i * 3], parameters[0][i * 3 + 1], parameters[0][i * 3 + 2]);
	}

	std::vector<T> angles(num_points);
	for (size_t i = 0; i < num_points; ++i) {
		Eigen::Matrix<T, 3, 1> v1 = points_3d[(i - 1 + num_points) % num_points] - points_3d[i];
		Eigen::Matrix<T, 3, 1> v2 = points_3d[(i + 1) % num_points] - points_3d[i];
		angles[i] = ComputeAngle(v1, v2);
	}

	const T target_angle_1 = T(90.0);
	const T target_angle_2 = T(270.0);
	for (size_t i = 0; i < num_points; ++i) {
		T angle = angles[i];
		if (angle >= T(75.0) && angle <= T(105.0)) {
			residuals[i] = angle - target_angle_1;
		}
		else if (angle >= T(255.0) && angle <= T(285.0)) {
			residuals[i] = angle - target_angle_2;
		}
		else {
			residuals[i] = T(0.0);
		}
	}

	T angle_sum = T(0.0);
	for (const auto& angle : angles) {
		angle_sum += angle;
	}
	residuals[num_points] = angle_sum - T((num_points - 2) * 180);

	return true;
}

template<typename T>
T AngleOptimizationCostFunction::ComputeAngle(const Eigen::Matrix<T, 3, 1>& v1, const Eigen::Matrix<T, 3, 1>& v2) const
{
	T dot_product = v1.dot(v2);
	T norm_v1 = v1.norm();
	T norm_v2 = v2.norm();
	T cos_theta = dot_product / (norm_v1 * norm_v2);
	//避免值因为double赋值超过边界
	cos_theta = max(-T(1.0), min(T(1.0), cos_theta));

	T angle = acos(cos_theta);
	Eigen::Matrix<T, 3, 1> cross_product = v1.cross(v2);
	bool curDir = cross_product.z() < 0.0 ? false : true;
	if (curDir ^geoDir_) {
		angle = 2 * PI - angle;
	}

	return angle * 180.0 / PI;
}

//多边形直角化代码
void BecomeRectangular(std::vector<Eigen::Vector3d> points, std::vector<Eigen::Vector3d> & newPnts)
{
	
	 // 初始点（展平成一维数组）
	std::vector<double> initial_points;
	for (const auto& p : points) {
		initial_points.push_back(p.x());
		initial_points.push_back(p.y());
		initial_points.push_back(p.z());
	}
	const int num_points = points.size();
	const int num_residuals = num_points + 1; // 每个角度的残差 + 角度和的残差
	const int num_parameters = 3 * num_points;
	// 创建优化问题
	ceres::Problem problem;

	// 创建 DynamicAutoDiffCostFunction 对象,new出来的对象，会在Problem结束时候自动析构
	auto cost_function = new ceres::DynamicAutoDiffCostFunction<AngleOptimizationCostFunction>(
		new AngleOptimizationCostFunction(points));
	cost_function->AddParameterBlock(num_parameters);
	cost_function->SetNumResiduals(num_residuals);

	problem.AddResidualBlock(cost_function, nullptr, initial_points.data());

	// 配置求解器选项
	ceres::Solver::Options options;
	options.minimizer_type = ceres::TRUST_REGION;
	options.linear_solver_type = ceres::DENSE_QR;
	options.max_num_iterations = 500;

	// 执行优化
	ceres::Solver::Summary summary;
	ceres::Solve(options, &problem, &summary);

	for (int i = 0; i < points.size(); i++) {
		newPnts.emplace_back(initial_points[i * 3], initial_points[i * 3 + 1], initial_points[i * 3 + 2]);
	} 

}