任意离散点曲线求交点c++实现

本文介绍了一种使用C++实现的离散曲线交点计算算法。通过遍历多条曲线上的线段并求交点,适用于地理信息系统或计算机图形学中曲线分析任务。文章提供了完整的代码示例和matlab绘图验证。

摘要生成于 C知道 ,由 DeepSeek-R1 满血版支持, 前往体验 >

已知两条曲线上的点坐标(xi,yi),求二者交点。只需对曲线上线段进行遍历求线段交点即可,效果如下

                               

下面是c++代码实现

头文件 CalLineCrossPt.h

#include "stdafx.h"
#include <vector>
using namespace std;

typedef struct tagPosition
{
	double  x;
	double  y;
	tagPosition(double _x,double _y) { x=_x; y=_y;}
	tagPosition() {};
	bool operator==(const tagPosition & pt) { return (x==pt.x && y==pt.y);} 
}CPosition;

CPosition IsLineCross(CPosition pt1,CPosition pt2,CPosition pt3,CPosition pt4);
void CalTwoLineCrossPoint(CPosition *_coord1,int _num1, CPosition *_coord2, int _num2,CPosition mincal,CPosition maxcal, vector<CPosition> &_crosspts);
void CalMulLinesCrossPoint(CPosition **_coord, int _numline, int *_numpt, vector<CPosition> &_crosspts);

实现文件 CalLineCrossPt.cpp

//================================================================
// 功能:    离散曲线求交点
//
// 作者:    jiangjp   1034378054@qq.com
// 单位:    China University of Geosciences (Wuhan)
// 日期:    2013/8/13
//================================================================

#include "stdafx.h"
#include "CalLineCrossPt.h"

void CalMulLinesCrossPoint(CPosition **_coord, int _numline, int *_numpt, vector<CPosition> &_crosspts)
{
	float xmin1,xmax1,xmin2,xmax2;    // 存储两条线所在矩形区域的两个端点(最大最小坐标)
	float ymin1,ymax1,ymin2,ymax2;
	CPosition mincal,maxcal;          // 存储两条线所在区域重叠形成的矩形区域两端点(最大最小坐标)
	CPosition *linemin=new CPosition[_numline];   // 开辟空间存储每条测线的x,y值最小点
	CPosition *linemax=new CPosition[_numline];   // 开辟空间存储每条测线的x,y值最大点

	vector<CPosition> tmpcrosspt;
	for(int i=0;i<_numline;i++){            // 查找每条线最大最小值
		linemin[i].x=_coord[i][0].x;
		linemin[i].y=_coord[i][0].y;
		linemax[i].x=_coord[i][0].x;
		linemax[i].y=_coord[i][0].y;

		for(int j=1;j<_numpt[i];j++){        
			if(_coord[i][j].x>linemax[i].x)  linemax[i].x=_coord[i][j].x;
			if(_coord[i][j].x<linemin[i].x)  linemin[i].x=_coord[i][j].x;
			if(_coord[i][j].y>linemax[i].y)  linemax[i].y=_coord[i][j].y;
			if(_coord[i][j].y<linemin[i].y)  linemin[i].y=_coord[i][j].y;	
		}
	}

	for(int i=0;i<_numline-1;i++){        // 对所有线两两循环	
		xmin1=linemin[i].x;               // 获取第一条线最大最小值点
		xmax1=linemax[i].x;
		ymin1=linemin[i].y;
		ymax1=linemax[i].y;

		for(int j=i+1;j<_numline;j++){ 
			xmin2=linemin[j].x;           // 获取第二条线最大最小值点
			xmax2=linemax[j].x;
			ymin2=linemin[j].y;
			ymax2=linemax[j].y;

			if(xmin1>xmax2 || xmin2>xmax1 || ymin1>ymax2 || ymin2>ymax2)   // 如果两线所在区域重叠则不需要计算交点
				continue;
			else
			{
				mincal.x=xmin1>xmin2 ? xmin2 : xmin1;
				mincal.y=ymin1>ymin2 ? ymin2 : ymin1;

				maxcal.x=xmax1>xmax2 ? xmax1 : xmax2;
				maxcal.y=ymax1>ymax2 ? ymax1 : ymax2;
			}

			CalTwoLineCrossPoint(_coord[i],_numpt[i], _coord[j],_numpt[j],mincal, maxcal,tmpcrosspt);

			int numcrosspt=tmpcrosspt.size();
			for (int k=0;k<numcrosspt;k++){
				_crosspts.push_back(tmpcrosspt.at(k));
			}

			tmpcrosspt.clear();
		}
	}

	delete []linemin;
	delete []linemax;
}


void CalTwoLineCrossPoint(CPosition *_coord1,int _num1, CPosition *_coord2, int _num2,CPosition mincal,CPosition maxcal, vector<CPosition> &_crosspts)         // 计算两条线交点
{     //  _coord1, _coord2分别为2条线坐标,_num1,_num2分别为2条线上点数目,mincal ,maxcal为两条线相交区域的矩形区域两个角点,_crosspts为所求的交点

	CPosition pos1,pos2,pos3,pos4;
	bool Found=false;

	for(int i=0;i<_num1-1;i++){        // 对条线相邻两点组成的线段依次求交点
		for(int j=0;j<_num2-1;j++){
			pos1=_coord1[i];           //  获取相邻点线段
			pos2=_coord1[i+1];
			pos3=_coord2[j];
			pos4=_coord2[j+1];
			// 判断所求线段是否在两测线相交区域内,如果在则进行求交点,否则不求
			if(pos1.x>=mincal.x && pos1.x<=maxcal.x && pos1.y>=mincal.y && pos1.y<=maxcal.y      
				&& pos2.x>=mincal.x && pos2.x<=maxcal.x && pos2.y>=mincal.y && pos2.y<=maxcal.y
				&& pos3.x>=mincal.x && pos3.x<=maxcal.x && pos3.y>=mincal.y && pos3.y<=maxcal.y
				&& pos4.x>=mincal.x && pos4.x<=maxcal.x && pos4.y>=mincal.y && pos4.y<=maxcal.y )
			{
				CPosition cpt=IsLineCross(pos1,pos2,pos3,pos4);      //  对两线段求交点

				if(cpt.x!=-1 && cpt.y!=-1){
					_crosspts.push_back(cpt);
				}
			}
		}
	}
}


CPosition IsLineCross(CPosition pt1,CPosition pt2,CPosition pt3,CPosition pt4)
{                         // 线段pt1pt2用1标记,线段pt3pt4用2标记,如k1,k2
	float xmin1=pt1.x>pt2.x ? pt2.x: pt1.x;
	float xmax1=pt1.x>pt2.x ? pt1.x: pt2.x;
	float ymin1=pt1.y>pt2.y ? pt2.y: pt1.y;
	float ymax1=pt1.y>pt2.y ? pt1.y: pt2.y;

	float xmin2=pt3.x>pt4.x ? pt4.x: pt3.x;
	float xmax2=pt3.x>pt4.x ? pt3.x: pt4.x;
	float ymin2=pt3.y>pt4.y ? pt4.y: pt3.y;
	float ymax2=pt3.y>pt4.y ? pt3.y: pt4.y;

	if(xmin1>xmax2 || xmin2>xmax1 || ymin1>ymax2 || ymin2>ymax2)
		return CPosition(-1,-1);

	if(pt1.x==pt2.x)      // 当线段1斜率不存在
	{
		if(pt3.x==pt4.x){  // 当线段2斜率不存在	
			return CPosition(-1,-1);   // 当交点不存在返回(-1,-1)点
		}
		else              // 当线段2斜率存在  
		{
			float k2=(pt4.y-pt3.y)/(pt4.x-pt3.x);   // 计算线段2斜率
			float x=pt1.x;
			float y=k2*(x-pt3.x)+pt3.y;         // 将线段1的点x坐标带入线段2
			float xmin2=pt3.x>pt4.x ? pt4.x : pt3.x;
			float xmax2=pt3.x>pt4.x ? pt3.x : pt4.x;
			float ymin1=pt1.y>pt2.y ? pt2.y :pt1.y;
			float ymax1=pt1.y>pt2.y ? pt1.y :pt2.y;

			if( x>=xmin2 && x<=xmax2 ){     // 不考虑交点在两线段端点
				if( y>=ymin1 && y<=ymax1 ){ // 如果所求y在线段1的范围之内则有交点
					return CPosition(x,y);    // 返回交点
				}
				else{
					return CPosition(-1,-1);
				}
			}
			else{
				return CPosition(-1,-1);
			}
		}
	}
	else   // 如果线段1斜率存在
	{ 
		if(pt3.x==pt4.x){    // 如果线段2斜率不存在
		
			float k1=(pt2.y-pt1.y)/(pt2.x-pt1.x);
			float x=pt3.x;
			float y=k1*(x-pt1.x)+pt1.y;
			float xmin1=pt1.x>pt2.x ? pt2.x :pt1.x;
			float xmax1=pt1.x>pt2.x ? pt1.x :pt2.x;
			float ymin2=pt3.y>pt4.y ? pt4.y :pt3.y;
			float ymax2=pt3.y>pt4.y ? pt3.y :pt4.y;

			if( x>=xmin1 && x<=xmax1 ){  //如果交点位于线段端点视为没有交点
				if( y>=ymin2 && y<=ymax2 ){   
					return CPosition(x,y);
				}
				else{
					return CPosition(-1,-1);
				}
			}
			else{
				return CPosition(-1,-1);
			}
		}
		else   // 如果线段1和线段2斜率都存在
		{
			float k1=(pt2.y-pt1.y)/(pt2.x-pt1.x);   // 计算线段1和线段2的斜率
			float k2=(pt4.y-pt3.y)/(pt4.x-pt3.x);

			if(k1==k2)
				return CPosition(-1,-1);

			float x=((pt3.y-pt1.y+k1*pt1.x-k2*pt3.x))/(k1-k2);   // 计算两线段的交点
			float y=k1*(x-pt1.x)+pt1.y;
			float xmin1=pt1.x>pt2.x ? pt2.x : pt1.x;
			float xmax1=pt1.x>pt2.x ? pt1.x : pt2.x;
			float xmin2=pt3.x>pt4.x ? pt4.x : pt3.x;
			float xmax2=pt3.x>pt4.x ? pt3.x : pt4.x;

			if((x>=xmin1 && x<=xmax1) && (x>=xmin2 && x<=xmax2)){  // 判断线段交点是否在两线段x范围之内
				return CPosition(x,y);
			}
			else{
				return CPosition(-1,-1);
			}
		}
	}
}

主程序 main

// CurveCrossTest.cpp : 定义控制台应用程序的入口点。
//

#include "stdafx.h"

#include <vector>
using namespace std;

#include "CalLineCrossPt.h"

int _tmain(int argc, _TCHAR* argv[])
{
	vector<CPosition> a;

	float line1_x[9]={1, 3, 4.5, 5.5, 6, 6,   8, 9,   10};
	float line1_y[9]={1, 3, 0,   0.5, 2, 2.7, 0, 1.5, 1.3};

	float line2_x[10]={0, 1, 1.3, 3,   4,   5.5, 7, 7, 8, 9};
	float line2_y[10]={1, 3, 3,   1.5, 1.9, 8,   0, 8, 3, 0};

	float line3_x[6]={1, 2, 3, 4, 5, 6};
	float line3_y[6]={6, 6, 6, 6, 6, 5.5};

	int line_num=3;
	int *numpt_line=new int[line_num];
	numpt_line[0]=9;
	numpt_line[1]=10;
	numpt_line[2]=6;
	CPosition **coord=new CPosition *[line_num];
	for (int i=0;i<line_num;i++){
		coord[i]=new CPosition[numpt_line[i]];
	}

	for (int j=0;j<numpt_line[0];j++){
		coord[0][j]=CPosition(line1_x[j],line1_y[j]);
		printf("line1: %f %f\n",coord[0][j].x,coord[0][j].y);
	}

	for (int j=0;j<numpt_line[1];j++){
		coord[1][j]=CPosition(line2_x[j],line2_y[j]);
		printf("line2: %f %f\n",coord[1][j].x,coord[1][j].y);
	}

	for (int j=0;j<numpt_line[2];j++){
		coord[2][j]=CPosition(line3_x[j],line3_y[j]);
		printf("line3: %f %f\n",coord[2][j].x,coord[2][j].y);
	}



	vector<CPosition> crosspts;
    CalMulLinesCrossPoint(coord, line_num, numpt_line, crosspts);

	int N=crosspts.size();

	FILE *fp_m = fopen("crosspt.txt", "wt");

	for (int i = 0; i < N; i++){
		fprintf(fp_m, "%lf %lf\n", crosspts.at(i).x,crosspts.at(i).y);
	}
	fclose(fp_m);

	return 0;
}

程序运行生成交点坐标文件crosspt.txt,可用matlab绘制图像如下

                                   

 

附matlab绘图代码

clear all;
clc;
close all;

line1_x=[1, 3, 4.5, 5.5, 6, 6,   8, 9,   10];
line1_y=[1, 3, 0,   0.5, 2, 2.7, 0, 1.5, 1.3];

line2_x=[0, 1, 1.3, 3,   4,   5.5, 7, 7, 8, 9];
line2_y=[1, 3, 3,   1.5, 1.9, 8,   0, 8, 3, 0];

line3_x=[1, 2, 3, 4, 5, 6];
line3_y=[6, 6, 6, 6, 6, 5.5];

crosspt=textread('crosspt.txt');
[numcrosspt,n]=size(crosspt);

figure
h1=plot(line1_x,line1_y,'*-k'); hold on
h2=plot(line2_x,line2_y,'*-b');
h3=plot(line3_x,line3_y,'*-g');
h4=plot(crosspt(:,1),crosspt(:,2),'or','LineWidth',3);
legend([h1 h2 h3 h4],'Line1','Line2','Line3','CrossPoints');

 

matlab离散点连成的两曲线交点-intersections.m 本帖最后由 kastin 于 2012-12-29 11:47 编辑 引言     曾经思考过曲面交,结果发现是学术界的一个难题,并且也想出了一个当前广泛使用方法原理一样的近似解法(追踪法)。当然网上也有很多方法,只不过那些方法非常粗糙,无非就是meshgrid出离散网格,比较两曲面在某位置的坐标是否在某一精度范围内,然后标记显示之。这个方法仅仅当离散网格非常细的时候才比较精确。除此之外,还有个非常严重的问题:上面的“精度范围”不是你随心所欲给的,而且也没规律寻找,当给得不恰当的时候,在格点处两曲面点作比较,会出很多个符合要的点,或者一个也没有。这样就会使得交线非常曲折,甚至断裂等,严重影响精确度。 ———————————————————分割线————————————————————————     当然,既然有曲面交,那么也有曲线交,其基本结构就是两曲线交。只是曲线交问题,事先得澄清一些注意点:     1. 数学分析层面曲线交点,其实就是方程组解;     2. “曲线”概念包括“直线”(处处曲率半径为无穷大);     3. Matlab的重点是离散点 矩阵运算,因此所有运算都是基于离散的,因而这里的曲线并不是绝对光滑的。     4. 近似试探与未知函数表达式。 对于1,我想说的是,如果你想要得两曲线的精确交点,并且一个不漏,那就直接解方程组,不用看本帖下文; 对于2,直线在Matlab里面是两个点确定,因此交点如果是一段线(无穷个点)的情况,可能只是显示两端点为交点; 对于3,很简单的例子,参数方程 x=cos,y=sin 在数学分析(即连续空间)层面上是个圆,但是如果你在离散t的时候,间距比较大,那么最后Matlab绘制的图像不是圆,而是正多边形了。因此,此时我们讨论曲线交点是这个离散点连线的图形与其他图形的交点,而非圆与其他交点。这也是我在标题中加了“离散点连成”的修饰词,防止被误会。 对于4,既然是曲线交点,那么本方法可以作为方程组的近似解。当然,如果离散点够多,解的精确度可以保证,不过不能保证一个不漏。另外就是,对于一组离散点构成的曲线,很难知道它们的解析表达式,因此想通过非线性方程组解的方法交点,就不大可能了(不过你可以用曲线拟合出函数解析式),因此,本帖的方法将会是一个较为有效交点方法。     废话了那么多,下面就说说曲线交点方法吧。除了解方程组,很多人想到的方法就是“离散点 判断距离是否足够接近”,这个方法原理跟引言中曲面交的方法是一样的。因此缺点也是一样的——太粗糙了。网上这种方法的代码也很多,这里就不上了。 下面将阐述我的方法以及给出例子代码。     我有两种思路,一种是高级绘图层面的(不涉及到底层操作),一种是底层的。我只给出了第一种的代码,因为我不会底层操作。     思路一:既然matlab曲线绘图是通过有序离散点依次连线形成,也就是说,通过“以直代曲”的过程,那么曲线交点无非就是离散点(结点)或者两线段交点。这比上面直接用交点附近的结点替代交点方法要精确得多了。而两直线交点很容易,只要知道四个点坐标,那么交点精确坐标自然可以表示出来。这就是交点的原理。只是还有一些细节处理和要注意的地方,我会留到后面再详细说。     思路二:仔细观察两曲线交点的特性,很容易发现,其实交点就是操作系统底层绘图重叠的那些像素点。因此,只要给要绘制的像素点做个标记,将那些重合的点突出显示(比如换个颜色),那么就相当于显示出交点了。这种方法由于是本质性的,因此不会遗漏任何交点,而且精确度极高,适用范围广。Matlab提供的plot plot3 surf等绘图函数都属于高级绘图,底层绘图(或称低级绘图)只有line surface以及patch等少数函数。但是,这里的“底层”并非真正的底层,因为它还是经过封装了的,而C 的MFC里面直接用刷子绘图,那才是依靠操作系统完成的真正的“底层”绘图操作(包括所有窗口都是操作系统绘制的)。这里扯远了,想要说明的就是底层绘图的概念而已。只是我不会用matlab实现这些底层绘图。     上面说了思路,下面就详细说说一些注意点和需要处理的细节。     为了算法的健壮性,就必须考虑各种奇异的情况,防止bug。我们要考虑曲线有分支(很多代数曲线是这样的,代数几何里面研究的东西)、间断跳跃(有绝对值函数或者存在渐近线情况)、首尾是交点、在切点相交,等等这些情况。而且对于定位交点处附近的四个最近端点也是个问题(因为这里存在一个情况,如果曲线1上的一条线段与曲线2上的两条或者以上的线段相交,我的程序因为这个问题没能有效解决,出现在一些非常特殊的情况下会遗漏部分交点)。上面的情况如果不考虑,那么你的程序就会出现各种各样的问题。     对于通常情况,我考虑使用变号法则来判断交点(也就是高数里面“连续函数变号端点内存在零点”),对于上面说的特殊情况,那么预先处理,比如先看是否存在eps内的,或者为零的结点,有则直接记录,没有的话,通过两线段交来确定交点。至于遍历顺序的问题,为了简便,我指考虑两曲线离散点个数相同的情况(因为不同的话,会出现一些无法处理的情况),而且优先考虑离散点的坐标值中x或者y都相同的情况(比如x=0:0.1:pi; y1=sin, y2=x.^2这两条曲线的x值相同分布)。 下面是曲线y=cos.*exp)与y2=sin.^2 cos在[0:pi/18:2*pi]区间内的交点的代码: 注意:我没有写成接口的形式,虽然对于比那些较懒的人来说不太方便,但是这样做是为了让你能更好弄懂原理,并能自己改造代码。因此,下面的代码可以稍作修改,就能解决别的曲线交点。这样,不愿思考的懒人就没法达到自己的目的了~% 绘制两离散曲线交点 % 注意: %   1. 这里的“交点”指的是离散点连线绘出的图形的交点,而非函数或者方程理论分析上的交点, %      因此,这个程序不能作为根来用。 %   2. 要曲线离散点的个数一样。 %   3. 两个曲线出现参数方程的话,大多数情况正常。但是经测试发现,对于某些非常特殊的情况会出现bug, %      除非调用ezplot的数据(xdata,ydata)。 % %   by kastin @Mar 21, 2012 clear; debug=false; %关闭显示交点过程 % 曲线1 x=0:pi/18:2*pi; y=cos.*exp); % 曲线2 [x1 N]=sort;  %此处对于C1参数方程,C2为显式函数;或者均为参数方程时候有用 % 下面几句代码在本个案下没有什么特殊作用,但是当出现参数方程的时候,下面的方法改动一下就会有用。 y1=sin.^2 cos; %用于作图 x2=x; y2=sin.^2 cos; %用于寻点 h=plot; y<=eps)=0; y20; neg=cy<=0; %确定变号位置 fro=diff~=0; %变号的前导位置 rel=diff~=0; %变号的尾巴位置 zpf=find; %记录索引 zpr=find 1; %记录索引 zpfr=[zpf; zpr]; hold on % 观看交点过程 if debug, hp=plot,y,'r.-',x2,y2,'g.-'); end %线性交 x0=.*-y)-x.*-y))./ y2-y-y2); y0=y ).*-y)./-x); if any), y0=y2; end %加入已经判断为零的位置 x0=[x<=eps) x0].'; y0=[y<=eps) y0].'; hc=plot; %绘制交点 if debug, legend;hp],'C1','C2','交点','微线段1','微线段2',0); end legend xlabel, ylabel, zlabel; title axis equal hold off disp disp) %排除重复的点复制代码经测试十几种奇怪的曲线相交(包括参数方程形式的曲线),目前发现上述代码的方法有四种情况会出现遗漏一两个交点。(其实上面代码本意是显式函数的曲线交点,或者未知表达式的离散点曲线交点,并未针对参数方程,隐函数方程做优化,但是可以凑合着用用。)
评论 3
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值