平面几何

向量(Vector)
在几乎所有的几何问题中,向量(有时也称矢量)是一个基本点。向量的定义包含方向和一个数(长度)。在二维空间中,一个向量可以用一对x和y来表示。例如由点(1,3)到(5,1的向量可以用(4,-2)来表示。这里大家要特别注意,我这样说并不代表向量定义了起点和终点。向量仅仅定义方向和长度。

向量加法
向量也支持各种数学运算。最简单的就是加法。我们可以对两个向量相加,得到的仍然是一个向量。我们有:
    V1(x1, y1)+V2(x2, y2)=V3(x1+x2, y1+y2)
下图表示了四个向量相加。注意就像普通的加法一样,相加的次序对结果没有影响(满足交换律),减法也是一样的。
点看全图

点乘(Dot Product)
如果说加法是凭直觉就可以知道的,另外还有一些运算就不是那么明显的,比如点乘和叉乘。
点乘比较简单,是相应元素的乘积的和
    V1( x1, y1)   V2(x2, y2) = x1*x2 + y1*y2
注意结果不是一个向量,而是一个标量(Scalar)。点乘有什么用呢,我们有:
    A   B = |A||B|Cos(θ)
θ是向量A和向量B见的夹角。这里|A|我们称为向量A的模(norm),也就是A的长度, 在二维空间中就是|A| = sqrt(x2+y2)。这样我们就和容易计算两条线的夹角
    Cos(θ) = A   B /(|A||B|)

当然你知道要用一下反余弦函数acos()啦。(回忆一下cos(90)=0 和cos(0) = 1还是有好处的,希望你没有忘记。)这可以告诉我们如果点乘的结果,简称点积,为0的话就表示这两个向量垂直。当两向量平行时,点积有最大值
另外,点乘运算不仅限于2维空间,他可以推广到任意维空间。(译注:不少人对量子力学中的高维空间无法理解,其实如果你不要试图在视觉上想象高维空间,而仅仅把它看成三维空间在数学上的推广,那么就好理解了)
点看全图

叉乘(cross product)
相对于点乘,叉乘可能更有用吧。2维空间中的叉乘是:
    V1(x1, y1) X V2(x2, y2) = x1y2 – y1x2
看起来像个标量,事实上叉乘的结果是个向量,方向在z轴上。上述结果是它的模。在二维空间里,让我们暂时忽略它的方向,将结果看成一个向量,那么这个结果类似于上述的点积,我们有:
    A x B = |A||B|Sin(θ)
然而角度 θ和上面点乘的角度有一点点不同,他是有正负的,是指从A到B的角度。下图中 θ为负。
另外还有一个有用的特征那就是叉积的绝对值就是A和B为两边说形成的平行四边形的面积。也就是AB所包围三角形面积的两倍。在计算面积时,我们要经常用到叉积。
(译注:三维及以上的叉乘参看维基:http://en.wikipedia.org/wiki/Cross_product)
点看全图

点-线距离
找出一个点和一条线间的距离是经常遇见的几何问题之一。假设给出三个点,A,B和C,你想找出点C到点A、B定出的直线间距离。第一步是找出A到B的向量AB和A到C的向量AC,现在我们用该两向量的叉积除以|AB|,这就是我们要找的的距离了(下图中的红线)。
    d = (AB x AC)/|AB|
点看全图

如果你有基础的高中几何知识,你就知道原因了。上一节我们知道(AB X AC)/2是三角形ABC的面积,这个三角形的底是|AB|,高就是C到AB的距离。有时叉积得到的是一个负值,这种情况下距离就是上述结果的绝对值。
当我们要找点到线段的距离时,情况变得稍稍复杂一些。这时线段与点的最短距离可能是点到线段的某一端点,而不是点到直线的垂线。例如上图中点C到线段AB的最短距离应该是线段BC。我们有集中不同的方法来判断这种特殊情况。第一种情况是计算点积AB BC来判定两线段间夹角。如果点积大于等于零,那么表示AB到BC是在-90到90度间,也就是说C到AB的垂线在AB外,那么AB上到C距离最近的点就是B。同样,如果BAAC大于等于零,那么点A就是距离C最近的点。如果两者均小于零,那么距离最近的点就在线段AB中的莫一点。

 

代码
1 3. 浮点函数
2  //浮点几何函数库
3 #include <math.h>
4 #define eps 1e-8
5 #define zero(x) (((x)>0?(x):-(x))<eps)
6 struct point{double x,y;};
7 struct line{point a,b;};
8
9 //计算cross product (P1-P0)x(P2-P0)
10 double xmult(point p1,point p2,point p0){
11 return (p1.x-p0.x)*(p2.y-p0.y)-(p2.x-p0.x)*(p1.y-p0.y);
12 }
13 double xmult(double x1,double y1,double x2,double y2,double x0,double y0){
14 return (x1-x0)*(y2-y0)-(x2-x0)*(y1-y0);
15 }
16
17 //计算dot product (P1-P0).(P2-P0)
18 double dmult(point p1,point p2,point p0){
19 return (p1.x-p0.x)*(p2.x-p0.x)+(p1.y-p0.y)*(p2.y-p0.y);
20 }
21 double dmult(double x1,double y1,double x2,double y2,double x0,double y0){
22 return (x1-x0)*(x2-x0)+(y1-y0)*(y2-y0);
23 }
24
25 //两点距离
26 double distance(point p1,point p2){
27 return sqrt((p1.x-p2.x)*(p1.x-p2.x)+(p1.y-p2.y)*(p1.y-p2.y));
28 }
29 double distance(double x1,double y1,double x2,double y2){
30 return sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
31 }
32
33 //判三点共线
34 int dots_inline(point p1,point p2,point p3){
35 return zero(xmult(p1,p2,p3));
36 }
37 int dots_inline(double x1,double y1,double x2,double y2,double x3,double y3){
38 return zero(xmult(x1,y1,x2,y2,x3,y3));
39 }
40
41 //判点是否在线段上,包括端点
42 int dot_online_in(point p,line l){
43 return zero(xmult(p,l.a,l.b))&&(l.a.x-p.x)*(l.b.x-p.x)<eps&&(l.a.y-p.y)*(l.b.y-p.y)<eps;
44 }
45 int dot_online_in(point p,point l1,point l2){
46 return zero(xmult(p,l1,l2))&&(l1.x-p.x)*(l2.x-p.x)<eps&&(l1.y-p.y)*(l2.y-p.y)<eps;
47 }
48 int dot_online_in(double x,double y,double x1,double y1,double x2,double y2){
49 return zero(xmult(x,y,x1,y1,x2,y2))&&(x1-x)*(x2-x)<eps&&(y1-y)*(y2-y)<eps;
50 }
51
52 //判点是否在线段上,不包括端点
53 int dot_online_ex(point p,line l){
54 return dot_online_in(p,l)&&(!zero(p.x-l.a.x)||!zero(p.y-l.a.y))&&(!zero(p.x-l.b.x)||!zero(p.y-l.b.y));
55 }
56 int dot_online_ex(point p,point l1,point l2){
57 return dot_online_in(p,l1,l2)&&(!zero(p.x-l1.x)||!zero(p.y-l1.y))&&(!zero(p.x-l2.x)||!zero(p.y-l2.y));
58 }
59 int dot_online_ex(double x,double y,double x1,double y1,double x2,double y2){
60 return dot_online_in(x,y,x1,y1,x2,y2)&&(!zero(x-x1)||!zero(y-y1))&&(!zero(x-x2)||!zero(y-y2));
61 }
62
63 //判两点在线段同侧,点在线段上返回0
64 int same_side(point p1,point p2,line l){
65 return xmult(l.a,p1,l.b)*xmult(l.a,p2,l.b)>eps;
66 }
67 int same_side(point p1,point p2,point l1,point l2){
68 return xmult(l1,p1,l2)*xmult(l1,p2,l2)>eps;
69 }
70
71 //判两点在线段异侧,点在线段上返回0
72 int opposite_side(point p1,point p2,line l){
73 return xmult(l.a,p1,l.b)*xmult(l.a,p2,l.b)<-eps;
74 }
75 int opposite_side(point p1,point p2,point l1,point l2){
76 return xmult(l1,p1,l2)*xmult(l1,p2,l2)<-eps;
77 }
78
79 // 点关于直线的对称点 // by lyt
80 // 缺点:用了斜率
81 // 也可以利用"点到直线上的最近点"来做,避免使用斜率。
82 point symmetric_point(point p1, point l1, point l2) {
83 point ret;
84 if (l1.x > l2.x - eps && l1.x < l2.x + eps) {
85 ret.x = (2 * l1.x - p1.x);
86 ret.y = p1.y;
87 } else {
88 double k = (l1.y - l2.y ) / (l1.x - l2.x);
89 ret.x = (2*k*k*l1.x + 2*k*p1.y - 2*k*l1.y - k*k*p1.x + p1.x) / (1 + k*k);
90 ret.y = p1.y - (ret.x - p1.x ) / k;
91 }
92 return ret;
93 }
94
95 //判两直线平行
96 int parallel(line u,line v){
97 return zero((u.a.x-u.b.x)*(v.a.y-v.b.y)-(v.a.x-v.b.x)*(u.a.y-u.b.y));
98 }
99 int parallel(point u1,point u2,point v1,point v2){
100 return zero((u1.x-u2.x)*(v1.y-v2.y)-(v1.x-v2.x)*(u1.y-u2.y));
101 }
102
103 //判两直线垂直
104 int perpendicular(line u,line v){
105 return zero((u.a.x-u.b.x)*(v.a.x-v.b.x)+(u.a.y-u.b.y)*(v.a.y-v.b.y));
106 }
107 int perpendicular(point u1,point u2,point v1,point v2){
108 return zero((u1.x-u2.x)*(v1.x-v2.x)+(u1.y-u2.y)*(v1.y-v2.y));
109 }
110
111 //判两线段相交,包括端点和部分重合
112 int intersect_in(line u,line v){
113 if (!dots_inline(u.a,u.b,v.a)||!dots_inline(u.a,u.b,v.b))
114 return !same_side(u.a,u.b,v)&&!same_side(v.a,v.b,u);
115 return dot_online_in(u.a,v)||dot_online_in(u.b,v)||dot_online_in(v.a,u)||dot_online_in(v.b,u);
116 }
117 int intersect_in(point u1,point u2,point v1,point v2){
118 if (!dots_inline(u1,u2,v1)||!dots_inline(u1,u2,v2))
119 return !same_side(u1,u2,v1,v2)&&!same_side(v1,v2,u1,u2);
120 return dot_online_in(u1,v1,v2)||dot_online_in(u2,v1,v2)||dot_online_in(v1,u1,u2)||dot_online_in(v2,u1,u2);
121 }
122
123 //判两线段相交,不包括端点和部分重合
124 int intersect_ex(line u,line v){
125 return opposite_side(u.a,u.b,v)&&opposite_side(v.a,v.b,u);
126 }
127 int intersect_ex(point u1,point u2,point v1,point v2){
128 return opposite_side(u1,u2,v1,v2)&&opposite_side(v1,v2,u1,u2);
129 }
130
131 //计算两直线交点,注意事先判断直线是否平行!
132 //线段交点请另外判线段相交(同时还是要判断是否平行!)
133 point intersection(line u,line v){
134 point ret=u.a;
135 double t=((u.a.x-v.a.x)*(v.a.y-v.b.y)-(u.a.y-v.a.y)*(v.a.x-v.b.x))
136 /((u.a.x-u.b.x)*(v.a.y-v.b.y)-(u.a.y-u.b.y)*(v.a.x-v.b.x));
137 ret.x+=(u.b.x-u.a.x)*t;
138 ret.y+=(u.b.y-u.a.y)*t;
139 return ret;
140 }
141 point intersection(point u1,point u2,point v1,point v2){
142 point ret=u1;
143 double t=((u1.x-v1.x)*(v1.y-v2.y)-(u1.y-v1.y)*(v1.x-v2.x))
144 /((u1.x-u2.x)*(v1.y-v2.y)-(u1.y-u2.y)*(v1.x-v2.x));
145 ret.x+=(u2.x-u1.x)*t;
146 ret.y+=(u2.y-u1.y)*t;
147 return ret;
148 }
149
150 //点到直线上的最近点
151 point ptoline(point p,line l){
152 point t=p;
153 t.x+=l.a.y-l.b.y,t.y+=l.b.x-l.a.x;
154 return intersection(p,t,l.a,l.b);
155 }
156 point ptoline(point p,point l1,point l2){
157 point t=p;
158 t.x+=l1.y-l2.y,t.y+=l2.x-l1.x;
159 return intersection(p,t,l1,l2);
160 }
161
162 //点到直线距离
163 double disptoline(point p,line l){
164 return fabs(xmult(p,l.a,l.b))/distance(l.a,l.b);
165 }
166 double disptoline(point p,point l1,point l2){
167 return fabs(xmult(p,l1,l2))/distance(l1,l2);
168 }
169 double disptoline(double x,double y,double x1,double y1,double x2,double y2){
170 return fabs(xmult(x,y,x1,y1,x2,y2))/distance(x1,y1,x2,y2);
171 }
172
173 //点到线段上的最近点
174 point ptoseg(point p,line l){
175 point t=p;
176 t.x+=l.a.y-l.b.y,t.y+=l.b.x-l.a.x;
177 if (xmult(l.a,t,p)*xmult(l.b,t,p)>eps)
178 return distance(p,l.a)<distance(p,l.b)?l.a:l.b;
179 return intersection(p,t,l.a,l.b);
180 }
181 point ptoseg(point p,point l1,point l2){
182 point t=p;
183 t.x+=l1.y-l2.y,t.y+=l2.x-l1.x;
184 if (xmult(l1,t,p)*xmult(l2,t,p)>eps)
185 return distance(p,l1)<distance(p,l2)?l1:l2;
186 return intersection(p,t,l1,l2);
187 }
188
189 //点到线段距离
190 double disptoseg(point p,line l){
191 point t=p;
192 t.x+=l.a.y-l.b.y,t.y+=l.b.x-l.a.x;
193 if (xmult(l.a,t,p)*xmult(l.b,t,p)>eps)
194 return distance(p,l.a)<distance(p,l.b)?distance(p,l.a):distance(p,l.b);
195 return fabs(xmult(p,l.a,l.b))/distance(l.a,l.b);
196 }
197 double disptoseg(point p,point l1,point l2){
198 point t=p;
199 t.x+=l1.y-l2.y,t.y+=l2.x-l1.x;
200 if (xmult(l1,t,p)*xmult(l2,t,p)>eps)
201 return distance(p,l1)<distance(p,l2)?distance(p,l1):distance(p,l2);
202 return fabs(xmult(p,l1,l2))/distance(l1,l2);
203 }
204
205 //矢量V以P为顶点逆时针旋转angle并放大scale倍
206 point rotate(point v,point p,double angle,double scale){
207 point ret=p;
208 v.x-=p.x,v.y-=p.y;
209 p.x=scale*cos(angle);
210 p.y=scale*sin(angle);
211 ret.x+=v.x*p.x-v.y*p.y;
212 ret.y+=v.x*p.y+v.y*p.x;
213 return ret;
214 }
215
216

 

posted @ 2010-11-08 09:39  Cranny  阅读(800)  评论(0编辑  收藏  举报