地理信息系统软件开发中经常需要求取点、线、面之间的交点、交线、封闭区域面积和闭合集等结果，采用以矢量点乘和叉乘为基础的求取算法符合实际工作中已给出点位置和法向量等条件的情况，效率较高。
首先给出基本公式的推导。
矢量的结合率和交换率：
U+(V+W) = (U+V) + W;
U+V=V+U;
而设线段的起点和终点为P,Q;则中间任意一点R可以表示为:
R = P + r(Q - P); (0<r<1)
由P0,P1和P2三点构成的共面参数方程可表示为如下一种形式：
P(s,t) = (1-s-t)P0+sP1+tP2，其中s,t为参数
矢量的长度定义为：
设 V = (a,b); 则 |V| = sqrt(a*a + b*b); （即各轴数据平方和的开方）
而向量的规范化是 U = V/|V|，使得向量的长度为单位1。
矢量的点乘定义：
VW = v1w1+v2w2，V= (v1,v2) , W = (w1,w2);
此定义可从二维空间类推到三维和其他高次空间。
我们知道，二维平面上面一点P（x,y）在坐标轴转动某角度α后，新坐标为P(xnew,ynew)
设P点距离原点为r，此点与原X轴夹角为δ，则
x = rcosδ, y = rsinδ；
可求得新坐标为：xnew = rcos(δ-α)=rcosδcosα + rsinδsinα = xcosα + ysinα；
ynew =rsin(δ-α)=rsinδcosα –rcosδsinα = ycosα –xsinα；
下面推导广泛运用的一个公式：
VW = |V||W|cosϕ，ϕ为两向量的交角；
若V = (v1,0),既此点落在X轴上，则VW = v1w1 + 0w2 = v1w1；
因为V点落在轴上，则 |V| = v1,|W| cosϕ = w1；所以 VW = |V||W|cosϕ；
若为一般情况，即V= (v1,v2),则我们可旋转原始坐标轴角度δ后使得V点落在新X轴上，
由于： (v1w1 + v2w2)/|V| = w1cosδ+ w2sinδ = W点在新坐标系下的X轴数值 = |W| cosϕ
因此 v1w1 + v2w2 = |V||W|cosϕ；
在平时的运用中，计算cosϕ = VW / |V||W| ，可对两向量的交角大小进行判断，若VW为0，则说明两向量互相垂直；若VW大于0，则|ϕ|<90度，说明交角为锐角。
设向量V(a,b),对其旋转90度后，得到U= (acos90 + bsin90, ycos90-xsin90) = (b,-a)；可由此方法得到向量的垂直向量。
由VW = |V||W|cosϕ公式的推导中可以类推|V||W|sinϕ等于什么呢？推导如下：
|W|sinϕ = W点在新坐标系下的Y轴数值 = w2cosδ-w1sinδ= w2v1/|V|-w1v2/|V|=
(v1w2-v2w1)/|V|,即
v1w2-v2w1 = |V||W|sinϕ，ϕ为两向量的交角；
在平时的运用中，可采用此公式计算三角形和四边形的面积。也能判断两向量的位置关系，如果计算v1w2-v2w1为0，则说明两向量交角为0度或180度，即共线。如果计算结果大于0，则说明两向量交角0度到180度之间。
求取向量之间位置关系和三角形面积的参考代码：
// Input: three points P0, P1, and P2
// Return: >0 for P2 left of the line through P0 and P1
// =0 for P2 on the line
// <0 for P2 right of the line
double CDEMAlgorithm::isLeft( Point P0, Point P1, Point P2 )
{
return ( (P1.x - P0.x) * (P2.y - P0.y)
- (P2.x - P0.x) * (P1.y - P0.y) ); //使用右手法则确定叉乘的方向
}
double CDEMAlgorithm::orientation2D_Triangle( Point V0, Point V1, Point V2 )
{
return isLeft(V0, V1, V2); //判断位置关系
}
double CDEMAlgorithm::area2D_Triangle( Point V0, Point V1, Point V2 )
{
return isLeft(V0, V1, V2) / 2.0; //求三角形的面积
}
double CDEMAlgorithm::dist3D_Line_to_Line( Line L1, Line L2)
{
Vector u = L1.P1 - L1.P0;
Vector v = L2.P1 - L2.P0;
Vector w = L1.P0 - L2.P0;
double a = Dot(u,u); // always >= 0
double b = Dot(u,v);
double c = Dot(v,v); // always >= 0
double d = Dot(u,w);
double e = Dot(v,w);
double D = a*c - b*b; // always >= 0
double sc, tc;
// compute the line parameters of the two closest points
if (D < EPS) { // 几乎平行
sc = 0.0;
tc = (b>c ? d/b : e/c); //
}
else {
sc = (b*e - c*d) / D;
tc = (a*e - b*d) / D;
}
//sc,tc分别指示了在两条 直线上的 比例参数
Vector dP = w + (sc * u) - (tc * v); // = L1(sc) - L2(tc)
return sqrt(Dot(dP,dP)); // return the closest distance
}
附图:
