地理信息系统算法基础

.Vector

已于 2023-10-08 19:08:09 修改

阅读量755

点赞数 3

分类专栏： GIS相关文章标签：算法 c++

于 2022-07-03 17:17:23 首次发布

本文链接：https://blog.csdn.net/QomoL_/article/details/125374487

版权

GIS相关专栏收录该内容

3 篇文章 0 订阅

订阅专栏

储备知识

已知两点坐标A(a1,b1),B(a2,b2,),则向量AB=(a2-a1,b2-b1)。

(写向量时前后两个字母的顺序,前面的字母为向量的起点,后面的字母为向量的终点)

已知向量OA(x1,y1)、OB(x2,y2)

向量加法(满足平行四边形法则和三角形法则)：OA+OB=OC;OA+OB=(x1+x2,y1+y2);

向量减法(共同起点，指向被减)：OA-OB=BA;OA-OB=(x1-x2,y1-y2)

向量点乘(向量的内积):

①数学定义：OA·OB=x1·x2+y1·y2；②几何定义：OA·OB=|OA||OB|cosθ；

(向量点乘的结果是一个标量，几何意义为一个向量在另一个向量上的投影的长度)

向量叉乘(向量的外积)：

①数学定义：OA x OB=x1y2-x2y1;②几何定义：OA x OB=|OA||OB|sinθ；

(向量叉乘的结果是一个向量，数学上为垂直于OA、OB向量所在平面的向量，符合右手螺旋法则。几何意义为OA、OB向量构成的平行四边形的面积)

一、类&全局变量

class point //定义二维点
{
public:
	double x;
	double y;
	point() {};
	point(double x, double y)
	{
		this->x = x;
		this->y = y;
	}
};

class segment //定义线段
{
public:
	point startpoint;
	point endpoint;
	segment(point p1,point p2)
	{
		startpoint = p1;    //线段起点
		endpoint = p2;		//线段终点
	}	
};

class my_vector  //定义矢量
{
public:
	double x; 
	double y;
	my_vector() {}
	my_vector(double mx,double my)   //坐标构成矢量
	{
		this->x = mx;
		this->y = my;
	}
	my_vector(point pA, point pB)    //两点构成矢量
	{
		this->x = pB.x - pA.x;
		this->y = pB.y - pA.y;
	}
	my_vector add(my_vector v)//向量加法
	{ 
		my_vector temp;
		temp.x = this->x + v.x;
		temp.y = this->y + v.y;
		return temp;
	}
	my_vector sub(my_vector v)//向量减法
	{
		my_vector temp;
		temp.x = this->x - v.x;
		temp.y = this->y - v.y;
		return temp;
	}
	double dotmulit(my_vector v)//向量点乘
	{
		double temp;
		temp = this->x*v.x+this->y*v.y;		
		return temp;
	}
	double foxmulit(my_vector v)//向量叉乘
	{
		double temp;
		temp = this->x*v.y - v.x*this->y;
		return temp;//向量叉乘是还是向量，但是这里主要是利用计算结果，所以返回数值
	}
};

class grid//定义栅格像元
{
public:
	int row;
	int col;
	grid() {};
	grid(int r,int c)
	{
		this->row = r;
		this->col = c;
	}
};

double Ymax = 10000;//矢量数据栅格化所需参数
double Xmin = 0;
double deltaY = 2;
double deltaX = 2;

二、辅助方法

vector<segment> point2segment(vector<point> pointset) //点构成折线段
{
	vector<segment> mulitseg;
	for (int i = 0; i < pointset.size() - 1; i++)
	{
		segment s(pointset[i], pointset[i + 1]);
		mulitseg.push_back(s);
	}
	return mulitseg;
}

double pointsdistance(point p1, point p2)//计算两点之间的距离
{
	return sqrt((p1.x - p2.x)*(p1.x - p2.x) + (p1.y - p2.y)*(p1.y - p2.y));
}

double min(double a,double b)//求最小值
{
	return a < b ? a : b;
}

double max(double a, double b)//求最大值
{
	return a > b ? a : b;
}

grid VexsRowCol(point p)//矢量转栅格——求点所在行列号
{	
	grid g;
	g.row = (Ymax - p.y) / deltaY;
	g.col = (p.x - Xmin) / deltaX;
	return g;//返回栅格单元
}

三、实现方法

1.判断线段和线段是否相交

int issegmentcross(segment s1,segment s2)//判断线段和线段是否相交
{	
        //快速排斥
	if (max(s1.startpoint.x,s1.endpoint.x)>=min(s2.startpoint.x,s2.endpoint.x)&&
		max(s2.startpoint.x, s2.endpoint.x)>=min(s1.startpoint.x, s1.endpoint.x)&&
		max(s1.startpoint.y, s1.endpoint.y) >= min(s2.startpoint.y, s2.endpoint.y) &&
		max(s2.startpoint.y, s2.endpoint.y) >= min(s1.startpoint.y, s1.endpoint.y) )
	{	
        //跨立实验
		my_vector v1(s1.startpoint,s1.endpoint);
		my_vector v2(s1.startpoint,s2.startpoint);
		my_vector v3(s1.startpoint, s2.endpoint);
		double res1 = v1.foxmulit(v2)*v1.foxmulit(v3);

		my_vector v4(s2.startpoint, s2.endpoint);
		my_vector v5(s2.startpoint, s1.startpoint);
		my_vector v6(s2.startpoint, s1.endpoint);
		double res2 = v4.foxmulit(v5)*v5.foxmulit(v6);
		if (res1<=0&&res2<=0)//s2跨立s1，s1跨立s2
		{
			return 1;//相交
		}
		else
		{			
			return 0;//不相交
		}
	}
	else
	{
		return 0;//不相交
	}
}

2.判断点是否在线段上

int ispointinsegment(point p,segment s)//判断点是否在线段上
{	
	//1.判断两向量是否共线
	my_vector seg(s.startpoint,s.endpoint);
	my_vector ptemp(p,s.endpoint);
	double a = seg.foxmulit(ptemp);
	if (a==0)//共线
	{
		//2.判断点是否在矩形框内	
		if (min(s.startpoint.x,s.endpoint.x)<=p.x&&p.x<= max(s.startpoint.x, s.endpoint.x)
			&& min(s.startpoint.y, s.endpoint.y) <=p.y&&p.y<= max(s.startpoint.x, s.endpoint.x))
		{
			return 1;//在线段上
		}		
	}
	else
	{
		return 0;//不在线段上
	}
}

3.判断点是否在折线上

int ispointinpolyline(point p,vector<point> pointset)//判断点是否在折线上
{
	//1.点构成折线
	vector<segment> polyline=point2segment(pointset);
	//2.依次判断点是否在折线的组成线段上
	int count=0;
	for (int i = 0; i < polyline.size(); i++)
	{
		int a = ispointinsegment(p, polyline[i]);
		if (a==1)
		{
			count++; //如果点在组成线段上，计数器+1
		}
	}
	if (count!=0) //只要计数器不等于0
	{
		return 1;//点在折线上
	}
	else
	{
		return 0;//其他情况，点不在折线上
	}
}

4.判断点是否在多边形内

int ispointinpolygon(point p, vector<point> pointset)//判断点是否在多边形内
{
	if (pointset.size()<3) return 3;//如果小于3个点，不是多边形就返回
		
		int count=0;//计数器
		point rayp(999999999, p.y); segment ray(p, rayp); //做p的射线L(点rapy的x值要设置一个很大很大的数值，这里暂定999999999)

		//点构成多边形外围线段	
		pointset.push_back(pointset[0]);//多边形的起始点加入到尾部，连成环
		vector<segment> polygon = point2segment(pointset);

		for (int i = 0; i < polygon.size(); i++)
		{
			if (ispointinsegment(p, polygon[i]) == 1)//点在构成多边形的线段上
			{
				return 1;
			}
			else if (polygon[i].startpoint.y != polygon[i].endpoint.y)//多边形的边不是水平的
			{
				//多边形的边的一个端点在射线L上
				double avg = (polygon[i].startpoint.y*polygon[i].endpoint.y) / 2;

				//边s的一个端点在L上，且该端点是s两端点中纵坐标较大的端点
				if (ispointinsegment(polygon[i].startpoint, ray) == 1&& polygon[i].startpoint.y> avg 
					|| ispointinsegment(polygon[i].endpoint, ray) == 1 && polygon[i].endpoint.y > avg)
				{
					count++;				
				}
				//多边形的边s和射线L相交
				else if (issegmentcross(polygon[i], ray) == 1)
				{
					count++;
				}
			}
		}
		//射线法规则：交点个数为奇数则在多边形内，为偶数则在多边形外
		if (count % 2 == 0)
		{
			return 0;
		}
		else
		{
			return 1;
		}
}

5.求任意多边形面积和几何中心

double polygonAandC(vector<point> pointset)//求多边形面积和几何中心
{
	if (pointset.size() < 3) return 3;//如果不是多边形，返回

	//已知有序线段，则多边形面积为：任意一个点与顺序相邻两点组成的三角形面积之和
	point anypoint(pointset[0].x, pointset[0].y); //这里默认把任意点设为起始点
	double AreaSum=0, cenx=0, ceny=0;
	for (int i = 1; i < pointset.size()-1; i++)
	{	
		my_vector v1(anypoint, pointset[i]);
		my_vector v2(anypoint, pointset[i + 1]);
		double AreaPart = v1.foxmulit(v2);
		AreaSum += AreaPart; 

		//求几何中心		
		my_vector v = v1.add(v2);
		cenx += v.x*AreaPart;
		ceny += v.y*AreaPart;
	}
	AreaSum = AreaSum / 2;//向量叉乘是平行四边形的面积，除以2才是三角形的面积	
	cenx = cenx / (AreaSum*3);
	ceny = ceny / (AreaSum*3);
	return AreaSum;
}

6.过点，求线段AB垂足

point verticalpoint(point c,segment s)//过已知点C，做已知线段AB的垂足
{
	my_vector AC(s.startpoint,c);
	my_vector AB(s.startpoint,s.endpoint);
	double ap = AC.dotmulit(AB);//得到AC在AB上的投影长度
	double L = pointsdistance(s.startpoint,s.endpoint);//线段AB长度

	point Vpoint;
	Vpoint.x= s.startpoint.x + (s.endpoint.x - s.startpoint.x)*ap / L;
	Vpoint.y= s.startpoint.y + (s.endpoint.y - s.startpoint.y)*ap / L;
	return Vpoint;
}

7.过点，求到线段AB的垂距

double verticaldis(segment s,point p)//已知线段AB，已知有点P，求垂距
{
	my_vector AB(s.startpoint,s.endpoint);
	my_vector AP(s.startpoint,p);
	double area = AB.foxmulit(AP); //向量AP、AB叉乘，得到平行四边形面积
	double ab = pointsdistance(s.startpoint,s.endpoint);//求线段长度，向量AB模长
	double Vdis = area / ab;//得到垂距
	return Vdis;
}

8.求线段AB延长距离d点坐标

point EXtendsegment(segment s,point p,double d)//已知一线段 AB，在某个方向选择一点 P，输入距离 d，求线段 AB延长距离 d后的点坐标；
{
	//求取线段AB的长度L
	double L = pointsdistance(s.startpoint,s.endpoint);

	//判断点A和点B谁距P最近
	double AP = pointsdistance(s.startpoint, p);
	double BP = pointsdistance(s.endpoint, p);

	if (AP>BP)//距离B点近
	{
		point newB;
		newB.x = s.endpoint.x + (s.endpoint.x - s.startpoint.x)*d / L;
		newB.y = s.endpoint.y + (s.endpoint.y - s.startpoint.y)*d / L;
		return newB;
	}
	else  //距离A点近
	{
		point newA;
		newA.x = s.startpoint.x + (s.startpoint.x - s.endpoint.x)*d / L;
		newA.y = s.startpoint.y + (s.startpoint.y - s.endpoint.y)*d / L;
		return newA;
	}
}

9.求点到线段L的最短距离

double minDisPoint2segment(segment s,point p)//求点 P(x,y)到线段 L的最短距离；
{
	if (ispointinsegment(p,s)==1)//点在线段上
	{
		return 0;//最短距离为0
	}

	point czu = verticalpoint(p,s);//求垂足
	if (ispointinsegment(czu, s) == 1)//垂足在线段上
	{
		return verticaldis(s,p);//最短距离为垂距
	}
	else
	{
		//p点离哪个端点近就返回该最短距离
		return min(pointsdistance(s.startpoint,p),pointsdistance(s.endpoint,p));
	}
}

10.过点，做线段AB平行线

segment drawparallelLine(segment s,point c)//已知线段AB，过线段外的点C做平行线
{
	point p = verticalpoint(c,s);//求垂足
	double dx = s.endpoint.x - s.startpoint.x;
	double dy = s.endpoint.y - s.startpoint.y;
	double ap = pointsdistance(s.startpoint,p);
	double bp = pointsdistance(s.endpoint, p);
	if (ap<=bp)//距离A点近
	{
		point d;
		d.x = c.x + dx;
		d.y = c.y + dy;
		segment s(c, d);
		return s;
	}
	else//距离B点近
	{
		point d;
		d.x = c.x - dx;
		d.y = c.y - dy;
		segment s(c, d);
		return s;
	}
}

11.间隔取点实现矢量数据压缩

vector<segment> getIntervalpoints(vector<point> pointset)//编程实现间隔取点法达到矢量数据的压缩；
{
	vector<point> newpointset;
	for (int i = 0; i < pointset.size(); i+=2)
	{
		newpointset.push_back(pointset[i]);
	}
	if (pointset.size()%2==0) //输入点集为偶数个
	{
		newpointset.push_back(pointset[pointset.size()-1]);//首尾点都要保留
	}
	vector<segment> compressLine = point2segment(newpointset);//点转线
	return compressLine;
}

12.链式编码实现栅格数据压缩

vector<int> ChainCodes(vector<grid> gridset)//编写链式编码程序实现栅格数据的压缩；(默认输入的栅格是有顺序的)
{
	vector<int> result;//存储编码
	result.push_back(gridset[0].row);
	result.push_back(gridset[0].col);

	for (int i = 0; i < gridset.size() - 1; i++)
	{
		if (gridset[i + 1].row == gridset[i].row&&gridset[i + 1].col == gridset[i].col + 1)
		{
			result.push_back(0);
		}
		else if (gridset[i + 1].row == gridset[i].row + 1 && gridset[i + 1].col == gridset[i].col + 1)
		{
			result.push_back(1);
		}
		else if (gridset[i + 1].row == gridset[i].row + 1 && gridset[i + 1].col == gridset[i].col)
		{
			result.push_back(2);
		}
		else if (gridset[i + 1].row == gridset[i].row + 1 && gridset[i + 1].col == gridset[i].col - 1)
		{
			result.push_back(3);
		}
		else if (gridset[i + 1].row == gridset[i].row && gridset[i + 1].col == gridset[i].col - 1)
		{
			result.push_back(4);
		}
		else if (gridset[i + 1].row == gridset[i].row - 1 && gridset[i + 1].col == gridset[i].col - 1)
		{
			result.push_back(5);
		}
		else if (gridset[i + 1].row == gridset[i].row - 1 && gridset[i + 1].col == gridset[i].col)
		{
			result.push_back(6);
		}
		else if (gridset[i + 1].row == gridset[i].row - 1 && gridset[i + 1].col == gridset[i].col + 1)
		{
			result.push_back(7);
		}
	}
	return result;
}

13.实现八方向矢量线栅格化

（还有点问题）

void EightDircVector2Grid(vector<point> pointset)//八方向矢量线栅格化
{
	vector<segment> mulitline = point2segment(pointset); //点转线
	vector<grid> resultGrid;
	for (int i = 0; i < mulitline.size(); i++)//逐线段
	{
		grid g1 = VexsRowCol(mulitline[i].startpoint);
		grid g2 = VexsRowCol(mulitline[i].endpoint);
		int derow = abs(g1.row - g2.row);//行差
		int decol = abs(g1.col - g2.col);//列差

		if (derow>decol)//行差大于列差
		{
			resultGrid.push_back(g1);//把首端点加入栅格容器
			for (int k = g2.row+1; k < g1.row; k++)//(暂时假设前一个点的行号大于后一个点)
			{
				point p;//求交点
				p.y = Ymax - (2 * k + 1) / 2 * deltaY;
				p.x = (p.y - mulitline[i].startpoint.y)*
					((mulitline[i].endpoint.x - mulitline[i].startpoint.x) /
					(mulitline[i].endpoint.y - mulitline[i].startpoint.y) +
						mulitline[i].startpoint.x);
				grid temp = VexsRowCol(p);
				resultGrid.push_back(temp);			
			}
			resultGrid.push_back(g2);//把尾端点加入栅格容器
		}
		else if (derow < decol)//列差大于行差
		{
			resultGrid.push_back(g1);//把首端点加入栅格容器
			for (int j = g1.col+1; j < g2.col; j++)//(暂时假设前一个点的列号小于后一个点)
			{
				point p;//求交点
				p.x = (2 * j + 1) / 2 * deltaX;
				p.y = (p.x - mulitline[i].startpoint.x)*
					((mulitline[i].endpoint.y - mulitline[i].startpoint.y) / 
					(mulitline[i].endpoint.x - mulitline[i].startpoint.x) +
						mulitline[i].startpoint.y);
				grid temp = VexsRowCol(p);//把交点转为栅格
				resultGrid.push_back(temp);//把该栅格加入栅格容器
			}
			resultGrid.push_back(g2);//把尾端点加入栅格容器
		}	
	}
}

14.递归求N的阶乘及阶层之和

int Nfactorial(int n)//求n的阶乘
{
	if (n==1)
	{
		return 1;
	}
	return n * Nfactorial(n-1);
}

int sumNfactorial(int n)//求1！+2！+3！+...+n!
{
	int sum = 0;
	for (int i = 1; i <= n; i++)
	{		
		sum += Nfactorial(i);		
	}
	return sum;
}

15.求斐波那契数列前n项之和

double Fibonacci(int n) //求1+2+3/2+5/3+8/5+···的前n项之和
{
	double a = 1, b = 1, c = 0, res = 1;
	for (int i = 2; i <= n; i++)//从第二项开始
	{
		c = a + b;
		a = b;
		b = c;
		res += c / a;
	}
	return res;
}

16.判断线段是否在多边形内

（还有点问题，实现代码下有附上算法描述）

int issegmentinpolygon(segment s, vector<point> pointset)//判断线段是否在多边形内
{
	if (pointset.size() < 3) return 3;//如果小于3个点，不是多边形就返回

	//点构成多边形外围线段	
	pointset.push_back(pointset[0]);//多边形的起始点加入到尾部，连成环
	vector<segment> polygon = point2segment(pointset);

	if (ispointinpolygon(s.startpoint, pointset)!=1 ||  //如果线段两端有任何一点不在多边形内
		ispointinpolygon(s.endpoint, pointset) != 1)
	{
		return 0;//线段不在多边形内
	}
	else
	{
		vector<point> temp;
		for (int i = 0; i < polygon.size(); i++)//遍历多边形边a
		{
			if (ispointinsegment(s.startpoint, polygon[i]) == 1) //如果线段的首端点在a上
			{
				temp.push_back(s.startpoint);
			}
			else if (ispointinsegment(s.endpoint, polygon[i]) == 1)//如果线段的尾端点在a上
			{
				temp.push_back(s.endpoint);
			}
			else if (ispointinsegment(polygon[i].startpoint, s))//如果边a的首点在线段s上
			{
				temp.push_back(polygon[i].startpoint);
			}
			else if (ispointinsegment(polygon[i].endpoint, s))//如果边a的尾点在线段s上
			{
				temp.push_back(polygon[i].endpoint);
			}
			else if (issegmentcross(polygon[i], s) == 1)//如果边a和线段s相交（内交）
			{
				return 0;//线段不在多边形内
			}
		}

		//将temp中的点按照X-Y坐标排序
		std::sort(temp.begin(), temp.end(), //按x排序
			[](point pt1, point pt2) {return pt1.x < pt2.x; });
		for (int i = 0; i < temp.size() - 1; i++)
		{
			point halfp;//求相邻两点的中点
			halfp.x = (temp[i].x + temp[i + 1].x) / 2;
			halfp.y = (temp[i].y + temp[i + 1].y) / 2;
			if (ispointinpolygon(halfp, pointset) != 1)//如果中点不在多边形内
			{
				return 0;//线段不在多边形内
			}
		}
		return 1;//线段在多边形内
	}	
}

if 线端PQ的端点不都在多边形内
      then return false;
    点集pointSet初始化为空;
    for 多边形的每条边s
      do if 线段的某个端点在s上
           then 将该端点加入pointSet;
         else if s的某个端点在线段PQ上
           then 将该端点加入pointSet;
         else if s和线段PQ相交 // 这时候已经可以肯定是内交了
           then return false;
    将pointSet中的点按照X-Y坐标排序;
    for pointSet中每两个相邻点 pointSet[i] , pointSet[ i+1]
      do if pointSet[i] , pointSet[ i+1] 的中点不在多边形中
           then return false;
    return true;

17.判断圆是否在多边形内

int iscircleinpolygon(point centre, double r, vector<point> pointset)//判断圆是否在多边形内
{
	if (pointset.size() < 3) return 3;//如果小于3个点，不是多边形就返回

	//点构成多边形外围线段	
	pointset.push_back(pointset[0]);//多边形的起始点加入到尾部，连成环
	vector<segment> polygon = point2segment(pointset);

	for (int i = 0; i < polygon.size(); i++)//遍历多边形的边
	{
		if (minDisPoint2segment(polygon[i],centre)<r)//圆心到边的最短距离小于半径
		{
			return 0;//圆不在多边形内
		}
	}
	return 1;
}

待更新，待完善。

.Vector

关注

3
点赞
踩
16

收藏

觉得还不错? 一键收藏
0
评论
地理信息系统算法基础

二、辅助方法三、实现方法1.判断线段和线段是否相交2.判断点是否在线段上3.判断点是否在折线上4.判断点是否在多边形内5.求任意多边形面积和几何中心6.过点，求线段AB垂足7.过点，求到线段AB的垂距8.求线段AB延长距离d点坐标9.求点到线段L的最短距离10.过点，做线段AB平行线11.间隔取点实现矢量数据压缩12.链式编码实现栅格数据压缩13.实现
复制链接

扫一扫