[计算几何]半平面交 POJ_2451 1474 1279 3335
半平面交,一开始以为很神奇的东西,看zzy的论文,看了好久,发现思路还好懂,不会写代码= =,只是觉得和melkman有点像.后来网上膜拜了各位牛人的代码,自己终于写出来了,A了2451,不过1474&1279的时候发现了自己代码的错误,如果有多条直线平行那么我的代码会错,后来又在3335发现如果有反向平行且不共线的两条直线那么我的会错.于是各种改过之后,觉得现在的代码应该没啥问题了= =||
主要思路就是先对边进行排序[按极角序],然后去重,把平行同向的直线只保留规约最紧的一条.然后对留下来的进行半平面交.现在觉得单独的半平面交不难,难的是把题目转化为半平面交
2451代码如下[如有错误还请指出]
1. /
2. 注意事项:
3. 1.单独的点是否是符合条件的
4. 2.对于输入边有没有两条反向的
5.
6. 对于1需要改变的是check里面判断时是>还是>=
7. 如果是>那么会忽略掉点 否则会留下
8.
9. 对于2的话需要特判一下是不是有两条边方向且不共线
10. 这里需要在排序然后去重之后做
11. /
12. #include
13. #include
14. #include
15. #include
16. #include
17.
18. using namespace std;
19. /点的结构体/
20. typedef struct
21. {
22. double x,y;
23. }TPoint;
24. /直线的结构体 angle是和X轴正方向的夹角/
25. typedef struct
26. {
27. TPoint p1,p2;
28. double angle;
29. }TLine;
30. / ax+by+c=0直线中的a b c/
31. typedef struct
32. {
33. double a,b,c;
34. }abcLine;
35.
36. const int maxnum=20016
37. const double eps=1e-8
38.
39. int n;
40. /pg存半平面最后的点/
41. TPoint pg[maxnum];
42. /line存输入的边 dq是双端队列/
43. TLine line[maxnum],dq[maxnum+2];
44. int idx;/the index of pg;/
45.
46. /加边 所有的边都表示的是
47. void Add_line(int idx,double x1,double y1,double x2,double y2)
48. {
49. line[idx].p1.x = x1;
50. line[idx].p1.y = y1;
51. line[idx].p2.x = x2;
52. line[idx].p2.y = y2;
53. line[idx].angle = atan2(y2-y1,x2-x1);
54. }
55. /输入/
56. void input()
57. {
58. int i;
59. double x1,y1,x2,y2;
60. for(i=0i
61. {
62. scanf(“%lf%lf%lf%lf”,x1,&y1,&x2,&y2);
63. Add_line(i,x1,y1,x2,y2);
64. }
65. /加外面的限制区域/
66. Add_line(n,0,0,10000,0);
67. Add_line(n+1,10000,0,10000,10000);
68. Add_line(n+2,10000,10000,0,10000);
69. Add_line(n+3,0,10000,0,0);
70. n+=4
71. }
72. /
73. 把一个double转化到一个int上
74. <-esp 的转化为-1
75. >eps 的转化为1
76. 其他的为0
77. /
78. int dou2int(double x)
79. {
80. if(x>eps)
81. return 1
82. if(x<-eps)
83. return -1
84. return 0
85. }
86.
87. /
88. 计算
89. /
90. double cross(TPoint a,TPoint b,TPoint c)
91. {
92. return (b.x-a.x)(c.y-a.y)-(b.y-a.y)(c.x-a.x);
93. }
94.
95. /
96. 多边形排序的模板函数
97. 首先按照极角序排列,如果极角序一样的话
98. 那么限制宽的排在前面
99. /
100. bool cmp(TLine l1,TLine l2)
101. {
102. int ret = dou2int(l1.angle-l2.angle);
103. if(ret > 0)
104. return true
105. if(ret < 0)
106. return false
107. if(dou2int(cross(l2.p1,l2.p2,l1.p2))>0)
108. return false
109. return true
110. }
111.
112. /
113. 由直线上的两点得到ax+by+c=0中的a b c
114. /
115. abcLine Get_line(TPoint p1,TPoint p2)
116. {
117. abcLine ret;
118. ret.a = p1.y-p2.y;
119. ret.b = p2.x-p1.x;
120. ret.c = p1.xp2.y-p2.xp1.y;
121. return ret;
122. }
123. /
124. 得到TLine结构体定义的两条直线的交点
125. /
126. TPoint Get_jiao(TLine l1,TLine l2)
127. {
128. TPoint ret;
129. abcLine A,B;
130. A = Get_line(l1.p1,l1.p2);
131. B = Get_line(l2.p1,l2.p2);
132. ret.x = (B.bA.c-A.bB.c)/(B.aA.b-A.aB.b);
133. ret.y = (A.aB.c-B.aA.c)/(B.aA.b-A.aB.b);
134. return ret;
135. }
136. /
137. 半平面交的时候判断当前平面是否合法
138. 是否需要退栈等
139. 如果单独的点不符合条件的话那么if语句
140. 的判断直接>0就行了,否则写成>=0
141. /
142. bool check(TLine l1,TLine l2,TLine l3)
143. {
144. TPoint cro = Get_jiao(l1,l2);
145. if(dou2int(cross(l3.p1,l3.p2,cro))>0)//加上等于0那么后来的直线过交点的时候不会退出
146. return false
147. return true
148. }
149. /
150. 判断两条直线是否平行反向且不共线
151. 下面的代码中需要加上不共线 需要对出去之后的直线进行判断
152. /
153. bool overline(TLine l1,TLine l2)
154. {
155. int ret;
156. int a,b,c,d;
157. a = (l1.p1.y-l1.p2.y)(l2.p1.x-l2.p2.x);
158. b = (l1.p1.x-l1.p2.x)(l2.p1.y-l2.p2.y);
159. c = l1.p1.x-l1.p2.x;
160. d = l2.p1.x-l2.p2.x;
161. ret = dou2int(cross(l1.p1,l1.p2,l2.p1));
162. if((a==b) & (cd<0) && (0 != ret))
163. {
164. //printf(“%f %f %f %f\n”,l1.p1.x,l1.p1.y,l1.p2.x,l1.p2.y);
165. //printf(“%f %f %f %f\n”,l2.p1.x,l2.p1.y,l2.p2.x,l2.p2.y);
166. return true
167. }
168. return false
169. }
170.
171. /
172. work 求半平面交 最后的顶点存在pg数组里面
173. 用到上面的TPoint TLine abcLine结构体
174. 会用到上面的
175. check :判断当前测试的直线是否在一直区域的范围内
176. dou2int :double转int <-eps对应-1 >eps 对应1 否则为0
177. cross :算叉积
178. Get_jiao :由两条直线得到交点
179. Get_line :由直线上的两个点得到直线的一般式
180. cmp等函数:sort用的比较函数,首先是旋转角,如果旋转角
181. 一样的话,那么要求松的放到前面,这样后面好写点
182. /
183. void work()
184. {
185. sort(line,line+n,cmp);
186. /首先对所有直线进行排序/
187. /
188. 把平行且限制宽的直线去掉
189. /
190. for(int i=1i
191. {
192. if(dou2int(line[i].angle-line[i-1].angle)==0)
193. {
194. memmove(line+i-1,line+i,(n-i)sizeof(line[0]));
195. n–;
196. i–;/防止多条直线平行/
197. }
198. }
199. /上面的for循环对直线进行”去重”/
200. int top,bot;
201. bot = 1
202. top=2
203. dq[bot]=line[0];
204. dq[top]=line[1];
205.
206. /
207. 半平面交核心算法
208. 和melkman算法有点像,每次考虑当前直线
209. 然后把这条直线加入到双端队列中
210. /
211. for(int i=2i
212. {
213. while(top>bot & check(dq[top-1],dq[top],line[i]))
214. top–;
215. while(top>bot & check(dq[bot+1],dq[bot],line[i]))
216. bot++;
217.
218. dq[++top]=line[i];
219.
220. }
221.
222. while(top>bot & check(dq[top-1],dq[top],dq[bot]))
223. top–;
224. while(top>bot & check(dq[bot+1],dq[bot],dq[top]))
225. bot++;
226. dq[–bot]=dq[top];/把最后一条直线加到最前面 这样方便计算/
227.
228. /
229. 输出由题目决定 这里表示无闭合区域
230. /
231. if(top-bot<3)
232. {/因为上面加了一条直线 所以这里是3/
233. printf(“0.0\n”);
234. return ;
235. }
236. idx=0
237. for(int i=bot;i
238. pg[idx++] = Get_jiao(dq[i],dq[i+1]);
239. idx–;
240. /下面用叉积算面积/
241. double ans=pg[idx].xpg[0].y-pg[0].xpg[idx].y;
242. for(int i=0i
243. ans += pg[i].xpg[i+1].y-pg[i+1].xpg[i].y;
244. ans /= 2.0
245. if(ans < 0)/如果是负的话转为正的/
246. ans = -ans;
247. /
248. 这里ans = -ans要注意别把0转化为-0了
249. 比如写成ans=ans>0?ans:-ans;的话会把0变成-0;
250. /
251. printf(“%.1f\n”,ans);
252. }
253.
254. int main(void)
255. {
256. #ifndef ONLINE_JUDGE
257. freopen(“2451.in”,“r”,stdin);
258. freopen(“2451.out”,“w”,stdout);
259. #endif
260. scanf(“%d”,n);
261. input();
262. work();
263. return 0
264. }
