You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

polygon.cc 11KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362
  1. /*
  2. Inside/Outside polygon filter
  3. Copyright (C) 2002 Robert Lipe, robertlipe+source@gpsbabel.org
  4. This program is free software; you can redistribute it and/or modify
  5. it under the terms of the GNU General Public License as published by
  6. the Free Software Foundation; either version 2 of the License, or
  7. (at your option) any later version.
  8. This program is distributed in the hope that it will be useful,
  9. but WITHOUT ANY WARRANTY; without even the implied warranty of
  10. MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  11. GNU General Public License for more details.
  12. You should have received a copy of the GNU General Public License
  13. along with this program; if not, write to the Free Software
  14. Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111 USA
  15. */
  16. #include "defs.h"
  17. #include "filterdefs.h"
  18. #include <stdio.h>
  19. #if FILTERS_ENABLED
  20. #define MYNAME "Polygon filter"
  21. static char* polyfileopt = NULL;
  22. static char* exclopt = NULL;
  23. /*
  24. * This test for insideness is essentially an odd/even test. The
  25. * traditional (simple) implementation of the odd/even test counts
  26. * intersections along a test ray, and if it should happen to exactly hit
  27. * a vertex of the polygon, it throws away the result and fires a different
  28. * test ray. Since we're potentially testing hundreds of test points against
  29. * a polygon with hundreds of sides, this seemed both wasteful and difficult
  30. * to coordinate, so instead I added extra state to try to figure out what we
  31. * would have seen if we had just missed the vertex. The result is the
  32. * hodgepodge of "limbo" states explained below. On the credit side of the
  33. * ledger, though, the tests for intersection are vastly simplified by always
  34. * having a horizontal test ray.
  35. *
  36. * The general structure of this filter is: we loop over the points in the
  37. * polygon. For each point, we update the state of each of the waypoints in
  38. * the test set. Thus, the state of every waypoint is indeterminate until
  39. * the end of the outer loop, at which point the state of every waypoint
  40. * should theoretically be completely determined.
  41. *
  42. * The bits following this comment encode the current state of the test point
  43. * as we go around the polygon. OUTSIDE clearly isn't a bit; it's just here
  44. * for completeness. If it's not INSIDE, and it's not something else, it's
  45. * clearly OUTSIDE.
  46. *
  47. * INSIDE is self-explanatory. What it means is that the last time we were
  48. * certain of our state, we were inside of the polygon.
  49. *
  50. * LIMBO encodes a bit more state information, to handle the case where our
  51. * test ray (a horizontal line) has intersected one of the vertices of the
  52. * polygon exactly. If the two lines that meet at that vertex are on
  53. * opposite sides of the test ray, it was an intersection. Otherwise, it just
  54. * grazed a local minimum or maximum and so counted as either zero or two
  55. * intersections - not a change in state. Horizontal lines encountered
  56. * while in limbo don't change the limbo state. All other lines do.
  57. * The rest of the bits talk about how we got into limbo, and thus what to do
  58. * when we get out of limbo.
  59. *
  60. * LIMBO_UP means that the last line segment we saw going into limbo was
  61. * headed upward. When we see another non-horizontal line segment, whether
  62. * we flip the INSIDE bit or not depends on whether it also goes upward. If
  63. * it does, we flip the INSIDE bit. Otherwise, we just clear our limbo state
  64. * and go on as if nothing had happened.
  65. *
  66. * LIMBO_BEGIN means that the very first vertex in the polygon put us into a
  67. * limbo state. We won't be able to resolve that limbo state until we get to
  68. * the end of the cycle, and it can actually coexist with another more local
  69. * limbo state. The next two bits talk about the beginning limbo state in
  70. * more detail.
  71. *
  72. * BEGIN_UP means that the first non-horizontal line segment we encountered
  73. * while in a LIMBO_BEGIN state went upward. As with LIMBO_UP, this is used
  74. * to determine the final disposition of the limbo state when we get back
  75. * around to the other end of the cycle.
  76. *
  77. * BEGIN_HOR is fairly temporary. It says that we've encountered one or more
  78. * horizontal line segments at the beginning of the cycle, so we haven't yet
  79. * been able to resolve the state of BEGIN_UP. It's a way of propagating the
  80. * "firstness" forward until we can make a decision, without propagating it
  81. * for every test point.
  82. *
  83. * UP is used to remember which way we were going in case we encounter a
  84. * limbo state.
  85. *
  86. * -- RLP
  87. */
  88. #define OUTSIDE 0
  89. #define INSIDE 1
  90. #define LIMBO 2
  91. #define LIMBO_UP 4
  92. #define LIMBO_BEGIN 8
  93. #define BEGIN_UP 16
  94. #define BEGIN_HOR 32
  95. #define UP 64
  96. typedef struct {
  97. unsigned short state;
  98. unsigned short override;
  99. } extra_data;
  100. static
  101. arglist_t polygon_args[] = {
  102. {
  103. "file", &polyfileopt, "File containing vertices of polygon",
  104. NULL, ARGTYPE_FILE | ARGTYPE_REQUIRED, ARG_NOMINMAX
  105. },
  106. {
  107. "exclude", &exclopt, "Exclude points inside the polygon",
  108. NULL, ARGTYPE_BOOL, ARG_NOMINMAX
  109. },
  110. ARG_TERMINATOR
  111. };
  112. static void polytest(double lat1, double lon1,
  113. double lat2, double lon2,
  114. double wlat, double wlon,
  115. unsigned short* state, int first, int last)
  116. {
  117. if (lat1 == wlat) {
  118. if (lat2 < wlat) {
  119. /* going down */
  120. if (*state & LIMBO) {
  121. if (*state & LIMBO_UP) {
  122. *state = *state ^ INSIDE;
  123. }
  124. *state = *state & ~LIMBO &~LIMBO_UP;
  125. } else if (*state & LIMBO_BEGIN) {
  126. if (*state & BEGIN_HOR) {
  127. *state = *state & ~BEGIN_HOR;
  128. } else if (last) {
  129. if (*state & BEGIN_UP) {
  130. *state = *state ^ INSIDE;
  131. }
  132. *state = *state & ~LIMBO_BEGIN & ~BEGIN_UP;
  133. }
  134. } else if (first && (lon1 > wlon)) {
  135. *state |= LIMBO_BEGIN;
  136. }
  137. } else if (lat2 == wlat) {
  138. if (first & (lon1 > wlon || lon2 > wlon)) {
  139. *state |= LIMBO_BEGIN | BEGIN_HOR;
  140. } else if (last && (*state & LIMBO_BEGIN) && (*state & LIMBO)) {
  141. if ((!!(*state & LIMBO_UP)) != (!!(*state & BEGIN_UP))) {
  142. *state = *state ^ INSIDE;
  143. }
  144. *state = *state & ~LIMBO & ~LIMBO_UP &
  145. ~LIMBO_BEGIN & ~BEGIN_UP;
  146. } else if (*state & LIMBO) {
  147. /* do nothing */
  148. } else {
  149. if (lon1 <= wlon && lon2 > wlon) {
  150. if (*state & UP) {
  151. *state &= ~UP;
  152. *state |= LIMBO_UP;
  153. }
  154. *state = *state | LIMBO;
  155. }
  156. }
  157. } else {
  158. /* going up */
  159. if (*state & LIMBO) {
  160. if (!(*state & LIMBO_UP)) {
  161. *state = *state ^ INSIDE;
  162. }
  163. *state = *state & ~LIMBO & ~LIMBO_UP;
  164. } else if (*state & LIMBO_BEGIN) {
  165. if (*state & BEGIN_HOR) {
  166. *state &= ~BEGIN_HOR;
  167. *state |= BEGIN_UP;
  168. } else if (last) {
  169. if (!(*state & BEGIN_UP)) {
  170. *state = *state ^ INSIDE;
  171. }
  172. *state = *state & ~LIMBO_BEGIN & ~BEGIN_UP;
  173. }
  174. } else if (first && (lon1 > wlon)) {
  175. *state |= LIMBO_BEGIN | BEGIN_UP;
  176. }
  177. }
  178. *state = *state & ~UP;
  179. } else if (lat2 == wlat) {
  180. if (lat1 < wlat) {
  181. if (last) {
  182. if (*state & BEGIN_UP) {
  183. *state = *state ^ INSIDE;
  184. }
  185. *state = *state & ~LIMBO_BEGIN & ~BEGIN_UP;
  186. } else if (lon2 > wlon) {
  187. *state |= LIMBO;
  188. }
  189. }
  190. /* no case for lat1==wlat; that's above */
  191. else {
  192. if (last) {
  193. if (!(*state & BEGIN_UP)) {
  194. *state = *state ^ INSIDE;
  195. }
  196. *state = *state & ~LIMBO_BEGIN & ~BEGIN_UP;
  197. } else if (lon2 > wlon) {
  198. *state |= LIMBO | LIMBO_UP;
  199. } else {
  200. *state |= UP;
  201. }
  202. }
  203. } else {
  204. if ((lat1 > wlat && lat2 < wlat) ||
  205. (lat1 < wlat && lat2 > wlat)) {
  206. /* we only care if the lines might intersect */
  207. if (lon1 > wlon && lon2 > wlon) {
  208. *state = *state ^ INSIDE;
  209. } else if (!(lon1 <= wlon && lon2 <= wlon)) {
  210. /* we're inside the bbox of a diagonal line. math time. */
  211. double loni = lon1+(lon2-lon1)/(lat2-lat1)*(wlat-lat1);
  212. if (loni > wlon) {
  213. *state = *state ^ INSIDE;
  214. }
  215. }
  216. }
  217. }
  218. }
  219. #define BADVAL 999999
  220. void
  221. polygon_process(void)
  222. {
  223. queue* elem, * tmp;
  224. Waypoint* waypointp;
  225. extra_data* ed;
  226. double lat1, lon1, lat2, lon2;
  227. double olat, olon;
  228. int fileline = 0;
  229. int first = 1;
  230. int last = 0;
  231. char* line;
  232. gbfile* file_in;
  233. file_in = gbfopen(polyfileopt, "r", MYNAME);
  234. olat = olon = lat1 = lon1 = lat2 = lon2 = BADVAL;
  235. while ((line = gbfgetstr(file_in))) {
  236. char* pound = NULL;
  237. int argsfound = 0;
  238. fileline++;
  239. pound = strchr(line, '#');
  240. if (pound) {
  241. *pound = '\0';
  242. }
  243. lat2 = lon2 = BADVAL;
  244. argsfound = sscanf(line, "%lf %lf", &lat2, &lon2);
  245. if (argsfound != 2 && strspn(line, " \t\n") < strlen(line)) {
  246. warning(MYNAME
  247. ": Warning: Polygon file contains unusable vertex on line %d.\n",
  248. fileline);
  249. } else if (lat1 != BADVAL && lon1 != BADVAL &&
  250. lat2 != BADVAL && lon2 != BADVAL) {
  251. #if NEWQ
  252. foreach(Waypoint* waypointp, waypt_list) {
  253. #else
  254. QUEUE_FOR_EACH(&waypt_head, elem, tmp) {
  255. waypointp = (Waypoint*)elem;
  256. #endif
  257. if (waypointp->extra_data) {
  258. ed = (extra_data*) waypointp->extra_data;
  259. } else {
  260. ed = (extra_data*) xcalloc(1, sizeof(*ed));
  261. ed->state = OUTSIDE;
  262. ed->override = 0;
  263. waypointp->extra_data = (extra_data*) ed;
  264. }
  265. if (lat2 == waypointp->latitude &&
  266. lon2 == waypointp->longitude) {
  267. ed->override = 1;
  268. }
  269. if (olat != BADVAL && olon != BADVAL &&
  270. olat == lat2 && olon == lon2) {
  271. last = 1;
  272. }
  273. polytest(lat1, lon1, lat2, lon2,
  274. waypointp->latitude,
  275. waypointp->longitude,
  276. &ed->state, first, last);
  277. first = 0;
  278. last = 0;
  279. }
  280. }
  281. if (olat != BADVAL && olon != BADVAL &&
  282. olat == lat2 && olon == lon2) {
  283. olat = BADVAL;
  284. olon = BADVAL;
  285. lat1 = BADVAL;
  286. lon1 = BADVAL;
  287. first = 1;
  288. } else if (lat1 == BADVAL || lon1 == BADVAL) {
  289. olat = lat2;
  290. olon = lon2;
  291. lat1 = lat2;
  292. lon1 = lon2;
  293. } else {
  294. lat1 = lat2;
  295. lon1 = lon2;
  296. }
  297. }
  298. gbfclose(file_in);
  299. #if NEWQ
  300. foreach(Waypoint* wp, waypt_list) {
  301. #else
  302. QUEUE_FOR_EACH(&waypt_head, elem, tmp) {
  303. Waypoint* wp = (Waypoint*) elem;
  304. #endif
  305. ed = (extra_data*) wp->extra_data;
  306. wp->extra_data = NULL;
  307. if (ed) {
  308. if (ed->override) {
  309. ed->state = INSIDE;
  310. }
  311. if (((ed->state & INSIDE) == OUTSIDE) == (exclopt == NULL)) {
  312. waypt_del(wp);
  313. delete wp;
  314. }
  315. xfree(ed);
  316. }
  317. }
  318. }
  319. void
  320. polygon_init(const char* args)
  321. {
  322. /* do nothing */
  323. }
  324. void
  325. polygon_deinit(void)
  326. {
  327. /* do nothing */
  328. }
  329. filter_vecs_t polygon_vecs = {
  330. polygon_init,
  331. polygon_process,
  332. polygon_deinit,
  333. NULL,
  334. polygon_args
  335. };
  336. #endif // FILTERS_ENABLED