图形算法库研发（六）：通用的二维曲线求交算法

wangzy1982

于 2024-04-04 22:49:54 发布

阅读量1.9k

点赞数 41

分类专栏：从零开始做一个世界级的三维几何算法库文章标签：算法

本文链接：https://blog.csdn.net/weixin_45387973/article/details/137375251

版权

从零开始做一个世界级的三维几何算法库专栏收录该内容

10 篇文章 1 订阅

订阅专栏

通用的二维参数曲线稳定高效求交算法

曲线求交算法存在的问题
WGP中通用二维参数曲线求交算法实现流程
求交结果的设计
测试代码
和OpenCaseCade的对比
源代码地址

曲线求交算法存在的问题

通用曲线求交算法是很多二维算法的基础，它看似简单，但要做到稳定高效，难度还是非常大的。
作者试验了OpenCaseCade和SISL两个比较著名的库，在某些病态情况下，它们都会出现低效和计算结果错误的情况。
一些病态情况，例如接近重合，在实际项目中还是非常常见的。
作者就曾经遇到过大量的CAD图纸因为绘制的问题导致看似重合实际不重合的情况，这些图纸在后续处理中，因为求交算法的不稳定，往往得到错误的结果，造成严重的损失。
总的来说，求交算法主要要解决下面几个问题。
1、如何得到所有交点，这是根本问题，但大多数数值算法只能得到单一交点。
2、曲线重合如何快速检测。
3、曲线接近重合时如何高效处理，这种情况下需要能得到满足误差要求的重合区域，这样在后续处理多边形和多段线的时候（例如多边形的布尔运算）更稳定。

WGP中通用二维参数曲线求交算法实现流程

本文只介绍通用二维曲线求交算法，特殊求交（例如直线与圆弧、圆弧与圆弧等）虽然是提速的好办法，但它们相对简单，只是工作量的问题，本文就不做针对性描述了。

步骤1：曲线初步分割
分割的原则如下：

如果是分段曲线（例如Nurbs）则在分段位置进行分割。
如果曲线绕过的角度大于Pi，则将曲线进行分割。

步骤2：快速解方程，确定根的情况
利用图形算法库研发（五）：非线性方程组求解中提到的方程组求解算法，借助如下方程，快速求出两条曲线根的情况。

X1(t1) = X2(t2)
Y1(t1) = Y2(t2)

方程组的结果分为两类：

精确解：交点
模糊解：可能相交的区域

步骤3：合并精确解
因为在有多个交点时候，方程组求解器采用的是分割求交。如果分割的位置正好是精确解的位置则可能存在某个解重复出现的情况，这时候需要对结果进行合并。
合并的结果也分两种情况：

精确解单独存在或者精确解只与精确解合并：这时候可以记录该合并后的精确解为交点（能合并成一个点）或者一段重合线（不能合并成一个点）。
精确解和模糊解合并：这时候，精确解成为模糊解的边界，并记录该边界端点是一个根。

步骤4：合并模糊解
合并模糊解的目的是将分割开的可能重合的区间合并成更大区域。
所以，合并的条件是两个区间相邻位置是根。

步骤5：计算模糊解边界
该步骤用来计算模糊解未知边界是否为根，如果一个模糊解两端为根，则很可能重合，后续开始重合检测。
这步的算法的思想如下：
根据第一根线边界点和该点垂线方向做一条直线，与第二条线相交。如果有交点，则判断交点和第一根线边界点是否重合；如果没有交点，怎用第二根线的端点和刚才的直线方向和第一根线相交计算交点并比较交点的重合情况。
交点的计算方法：
因为这种情况，直线和曲线接近垂直，初值也很好确定，所以可以先采用牛顿迭代快速计算，如果牛顿迭代失败（极小概率事件）再采用通用解方程组的方式来计算交点。

步骤6：重合段检测
两端点都为根的曲线段进行重合检测。
该步采用采样的方式，这种方式可以理解成通过高次拟合来逼近原直线，如果采样个数够多，准确率是有保证的。
采样点的计算方法同步骤4中的交点计算方法。
如果采样点不重合，则从采样点分割成两段进入后续步骤。

步骤7：精确求交
不重合的段，按接近重合处理。
这块是难点，一个好的方法是通过非线性变化，将接近重合在变化后界限明显。
本系统采用的算法是根据边界点和中间采样点的中心确定一个圆弧，如下图所示。将曲线上的点带入该圆弧的隐式方程，如果两条接近重合的曲线不相交，理论上会一条曲线上所有点都大于0，而另一条都小于0。这样就清楚的隔离开了两条曲线。上述方程我们添加到步骤1的方程组中一起求解，可以极大的提高收敛速度（这块就是作者在方程组求解器中提到的求解器要支持方程组个数大于变量个数）。
黑线是两条接近重合的曲线，红线是构造的圆弧线
步骤8：重复3到7步，直到所有的模糊解都被处理掉

步骤9：对交点和交线进行合并
需要合并的原因主要是计算过程中的分割，分割边界点可能重复出现。
例如曲线A被分成A1和A2，曲线B被分为B1和B2，如果交点恰巧在分割点处，则会计算出4个交点：A1B1、A1B2、A2B1、A2B2，需要对这些点进项合并。

求交结果的设计

交点的数据类型如下：

class WGP_API Curve2dCurve2dInt {
    public:
        Curve2dCurve2dIntType Type;
        void* Tag1;
        void* Tag2;
        Variable T1;
        Variable T2;
        Vector2d Point1;
        Vector2d Point2;
    };

需要说明的是：

Tag用来标记曲线信息，这在后续运算中（例如布尔运算）很有用。
交点类型分为交叉点、重合段起点、重合段内部采样点、重合段终点。本系统会保存采样点供后续算法做稳定性处理，起点、内部采样点和终点是有序的。
该类型只是当前测试类型，后续可能会加入更多属性，例如曲线交叉关系等。

除了交点外，系统还提供了一个数据结构和算法来对多组交点排序。
数据结构如下：

class WGP_API Curve2dCurve2dIntIndex {
    public:
        Array<Curve2dCurve2dInt>* Array;
        int StartIndex;
        int EndIndex;
    };
WGP_API Array<Curve2dCurve2dIntIndex> SortIntersections(Array<Curve2dCurve2dInt>* int_array_list, int int_array_count, CompareTagFunction compare_tag_function, bool is_sorted_by_first);

测试代码

本次同样采用圆弧来测试，好处是圆弧的形状更好控制。再次强调，这个测试比较的是通用求交算法，不是圆弧与圆弧的专用算法。
同时这个测试也演示了如何使用该求交算法
测试代码如下：

inline void TestCurve2dCurve2dInt1() {
	ArcCurve2d arc1(Vector2d(0, 0), 100, 0, 0, g_pi * 2);
	ArcCurve2d arc2(Vector2d(50, 0), 100, 0, 0, g_pi * 2);
	ArcCurve2d arc3(Vector2d(0, 0), 100, 0, 0, g_pi * 2);
	ArcCurve2d arc4(Vector2d(0.001, 0), 100, 0, 0, g_pi * 2);
	int n;
	Array<Curve2dCurve2dInt> result;
	Stopwatch sw;
	//cross
	//*
	std::cout << "=======cross=======\n";
	sw.Resume();
	for (int i = 0; i < 1000; ++i) {
		Array<Curve2dCurve2dInt> result;
		Intersect(&arc1, &arc2, nullptr, nullptr, 1E-6, result);
	}
	sw.Suspend();
	std::cout << sw.GetTotalTime() << "\n";
	std::cout << "==============\n";
	Intersect(&arc1, &arc2, nullptr, nullptr, 1E-6, result);
	n = 0;
	for (int i = 0; i < result.GetCount(); ++i) {
		if (result.GetPointer(i)->Type == Curve2dCurve2dIntType::Cross || result.GetPointer(i)->Type == Curve2dCurve2dIntType::OverlapBegin) {
			++n;
		}
	}
	std::cout << n << "\n";
	n = 0;
	while (n < result.GetCount()) {
		std::cout << "==============\n";
		if (result.GetPointer(n)->Type == Curve2dCurve2dIntType::Cross) {
			std::cout << result.Get(n).T1.Value << "\n";
			std::cout << result.Get(n).T2.Value << "\n";
			++n;
		}
		else {
			std::cout << result.Get(n).T1.Value << "\n";
			std::cout << result.Get(n).T2.Value << "\n";
			++n;
			while (n < result.GetCount()) {
				if (result.GetPointer(n)->Type == Curve2dCurve2dIntType::OverlapEnd) {
					std::cout << "->\n";
					std::cout << result.Get(n).T1.Value << "\n";
					std::cout << result.Get(n).T2.Value << "\n";
					++n;
					break;
				}
				++n;
			}
		}
	}
	sw.Reset();
	//*/
	//*
	//overlapping
	std::cout << "=======overlapping=======\n";
	sw.Resume();
	for (int i = 0; i < 1000; ++i) {
		Array<Curve2dCurve2dInt> result;
		Intersect(&arc1, &arc3, nullptr, nullptr, 1E-6, result);
	}
	sw.Suspend();
	std::cout << sw.GetTotalTime() << "\n";
	std::cout << "==============\n";
	result.Clear();
	Intersect(&arc1, &arc3, nullptr, nullptr, 1E-6, result);
	n = 0;
	for (int i = 0; i < result.GetCount(); ++i) {
		if (result.GetPointer(i)->Type == Curve2dCurve2dIntType::Cross || result.GetPointer(i)->Type == Curve2dCurve2dIntType::OverlapBegin) {
			++n;
		}
	}
	std::cout << n << "\n";
	n = 0;
	while (n < result.GetCount()) {
		std::cout << "==============\n";
		if (result.GetPointer(n)->Type == Curve2dCurve2dIntType::Cross) {
			std::cout << result.Get(n).T1.Value << "\n";
			std::cout << result.Get(n).T2.Value << "\n";
			++n;
		}
		else {
			std::cout << result.Get(n).T1.Value << "\n";
			std::cout << result.Get(n).T2.Value << "\n";
			++n;
			while (n < result.GetCount()) {
				if (result.GetPointer(n)->Type == Curve2dCurve2dIntType::OverlapEnd) {
					std::cout << "->\n";
					std::cout << result.Get(n).T1.Value << "\n";
					std::cout << result.Get(n).T2.Value << "\n";
					++n;
					break;
				}
				++n;
			}
		}
	}
	sw.Reset();
	//*/
	//*
	//close to overlapping
	std::cout << "=======close to overlapping=======\n";
	sw.Resume();
	for (int i = 0; i < 1000; ++i) {
		Array<Curve2dCurve2dInt> result;
		Intersect(&arc1, &arc4, nullptr, nullptr, 1E-6, result);
	}
	sw.Suspend();
	std::cout << sw.GetTotalTime() << "\n";
	std::cout << "==============\n";
	result.Clear();
	Intersect(&arc1, &arc4, nullptr, nullptr, 1E-6, result);
	n = 0;
	for (int i = 0; i < result.GetCount(); ++i) {
		if (result.GetPointer(i)->Type == Curve2dCurve2dIntType::Cross || result.GetPointer(i)->Type == Curve2dCurve2dIntType::OverlapBegin) {
			++n;
		}
	}
	std::cout << n << "\n";
	n = 0;
	while (n < result.GetCount()) {
		std::cout << "==============\n";
		if (result.GetPointer(n)->Type == Curve2dCurve2dIntType::Cross) {
			std::cout << result.Get(n).T1.Value << "\n";
			std::cout << result.Get(n).T2.Value << "\n";
			++n;
		}
		else {
			std::cout << result.Get(n).T1.Value << "\n";
			std::cout << result.Get(n).T2.Value << "\n";
			++n;
			while (n < result.GetCount()) {
				if (result.GetPointer(n)->Type == Curve2dCurve2dIntType::OverlapEnd) {
					std::cout << "->\n";
					std::cout << result.Get(n).T1.Value << "\n";
					std::cout << result.Get(n).T2.Value << "\n";
					++n;
					break;
				}
				++n;
			}
		}
	}
	sw.Reset();
	//*/
}

运行结果如下：

=======cross=======
0.0195802
==============
2
==============
1.31812
1.82348
==============
4.96507
4.45971
=======overlapping=======
0.0927182
==============
1
==============
0
0
->
6.28319
6.28319
=======close to overlapping=======
0.146492
==============
2
==============
1.57079
1.5708
->
1.5708
1.57081
==============
4.71239
4.71238
->
4.7124
4.71239

每种情况运行1000次，可见在相交、重合、接近重合情况下，效率和准确性都很高。

和OpenCaseCade的对比

下表比较的是在通用求交算法方面和OpenCaseCade进行的对比。
通用求交算法和OpenCaseCade的对比
从对比结果看：
相交情况下，WGP算法效率是OpenCaseCade的3倍。
重合情况下，WGP算法效率是OpenCaseCade的30倍。
接近重合情况，OpenCaseCade会得到错误的结果，效率依然很低。