poly_point_isect_py2py3.py 42 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911912913914915916917918919920921922923924925926927928929930931932933934935936937938939940941942943944945946947948949950951952953954955956957958959960961962963964965966967968969970971972973974975976977978979980981982983984985986987988989990991992993994995996997998999100010011002100310041005100610071008100910101011101210131014101510161017101810191020102110221023102410251026102710281029103010311032103310341035103610371038103910401041104210431044104510461047104810491050105110521053105410551056105710581059106010611062106310641065106610671068106910701071107210731074107510761077107810791080108110821083108410851086108710881089109010911092109310941095109610971098109911001101110211031104110511061107110811091110111111121113111411151116111711181119112011211122112311241125112611271128112911301131113211331134113511361137113811391140114111421143114411451146114711481149115011511152115311541155115611571158115911601161116211631164116511661167116811691170117111721173117411751176117711781179118011811182118311841185118611871188118911901191119211931194119511961197119811991200120112021203120412051206120712081209121012111212121312141215121612171218121912201221122212231224122512261227122812291230123112321233123412351236123712381239124012411242124312441245124612471248124912501251125212531254125512561257125812591260126112621263126412651266126712681269127012711272127312741275127612771278127912801281128212831284128512861287128812891290129112921293129412951296
  1. # BentleyOttmann sweep-line implementation
  2. # (for finding all intersections in a set of line segments)
  3. __all__ = (
  4. "isect_segments",
  5. "isect_polygon",
  6. # same as above but includes segments with each intersections
  7. "isect_segments_include_segments",
  8. "isect_polygon_include_segments",
  9. # for testing only (correct but slow)
  10. "isect_segments__naive",
  11. "isect_polygon__naive",
  12. )
  13. # ----------------------------------------------------------------------------
  14. # Main Poly Intersection
  15. # Defines to change behavior.
  16. #
  17. # Whether to ignore intersections of line segments when both
  18. # their end points form the intersection point.
  19. USE_IGNORE_SEGMENT_ENDINGS = True
  20. USE_DEBUG = True
  21. USE_VERBOSE = False
  22. # checks we should NOT need,
  23. # but do them in case we find a test-case that fails.
  24. USE_PARANOID = False
  25. # Support vertical segments,
  26. # (the bentley-ottmann method doesn't support this).
  27. # We use the term 'START_VERTICAL' for a vertical segment,
  28. # to differentiate it from START/END/INTERSECTION
  29. USE_VERTICAL = True
  30. # end defines!
  31. # ------------
  32. # ---------
  33. # Constants
  34. X, Y = 0, 1
  35. # -----------------------------------------------------------------------------
  36. # Switchable Number Implementation
  37. NUMBER_TYPE = 'native'
  38. if NUMBER_TYPE == 'native':
  39. Real = float
  40. ########################################################################
  41. # decreased 1e-10 to 1e-4 here, otherwise large float values of 10k+
  42. # caused not found errors in remove()
  43. NUM_EPS = Real("1e-4")
  44. ########################################################################
  45. NUM_INF = Real(float("inf"))
  46. elif NUMBER_TYPE == 'decimal':
  47. # Not passing tests!
  48. import decimal
  49. Real = decimal.Decimal
  50. decimal.getcontext().prec = 80
  51. NUM_EPS = Real("1e-10")
  52. NUM_INF = Real(float("inf"))
  53. elif NUMBER_TYPE == 'numpy':
  54. import numpy
  55. Real = numpy.float64
  56. del numpy
  57. NUM_EPS = Real("1e-10")
  58. NUM_INF = Real(float("inf"))
  59. elif NUMBER_TYPE == 'gmpy2':
  60. # Not passing tests!
  61. import gmpy2
  62. gmpy2.set_context(gmpy2.ieee(128))
  63. Real = gmpy2.mpz
  64. NUM_EPS = Real(float("1e-10"))
  65. NUM_INF = gmpy2.get_emax_max()
  66. del gmpy2
  67. else:
  68. raise Exception("Type not found")
  69. NUM_EPS_SQ = NUM_EPS * NUM_EPS
  70. NUM_ZERO = Real(0.0)
  71. NUM_ONE = Real(1.0)
  72. class Event:
  73. __slots__ = (
  74. "type",
  75. "point",
  76. "segment",
  77. # this is just cache,
  78. # we may remove or calculate slope on the fly
  79. "slope",
  80. "span",
  81. ) + (() if not USE_DEBUG else (
  82. # debugging only
  83. "other",
  84. "in_sweep",
  85. ))
  86. class Type:
  87. END = 0
  88. INTERSECTION = 1
  89. START = 2
  90. if USE_VERTICAL:
  91. START_VERTICAL = 3
  92. def __init__(self, type, point, segment, slope):
  93. assert(isinstance(point, tuple))
  94. self.type = type
  95. self.point = point
  96. self.segment = segment
  97. # will be None for INTERSECTION
  98. self.slope = slope
  99. if segment is not None:
  100. self.span = segment[1][X] - segment[0][X]
  101. if USE_DEBUG:
  102. self.other = None
  103. self.in_sweep = False
  104. # note that this isn't essential,
  105. # it just avoids non-deterministic ordering, see #9.
  106. def __hash__(self):
  107. return hash(self.point)
  108. def is_vertical(self):
  109. # return self.segment[0][X] == self.segment[1][X]
  110. return self.span == NUM_ZERO
  111. def y_intercept_x(self, x):
  112. # vertical events only for comparison (above_all check)
  113. # never added into the binary-tree its self
  114. if USE_VERTICAL:
  115. if self.is_vertical():
  116. return None
  117. if x <= self.segment[0][X]:
  118. return self.segment[0][Y]
  119. elif x >= self.segment[1][X]:
  120. return self.segment[1][Y]
  121. # use the largest to avoid float precision error with nearly vertical lines.
  122. delta_x0 = x - self.segment[0][X]
  123. delta_x1 = self.segment[1][X] - x
  124. if delta_x0 > delta_x1:
  125. ifac = delta_x0 / self.span
  126. fac = NUM_ONE - ifac
  127. else:
  128. fac = delta_x1 / self.span
  129. ifac = NUM_ONE - fac
  130. assert(fac <= NUM_ONE)
  131. return (self.segment[0][Y] * fac) + (self.segment[1][Y] * ifac)
  132. @staticmethod
  133. def Compare(sweep_line, this, that):
  134. if this is that:
  135. return 0
  136. if USE_DEBUG:
  137. if this.other is that:
  138. return 0
  139. current_point_x = sweep_line._current_event_point_x
  140. this_y = this.y_intercept_x(current_point_x)
  141. that_y = that.y_intercept_x(current_point_x)
  142. # print(this_y, that_y)
  143. if USE_VERTICAL:
  144. if this_y is None:
  145. this_y = this.point[Y]
  146. if that_y is None:
  147. that_y = that.point[Y]
  148. delta_y = this_y - that_y
  149. assert((delta_y < NUM_ZERO) == (this_y < that_y))
  150. # NOTE, VERY IMPORTANT TO USE EPSILON HERE!
  151. # otherwise w/ float precision errors we get incorrect comparisons
  152. # can get very strange & hard to debug output without this.
  153. if abs(delta_y) > NUM_EPS:
  154. return -1 if (delta_y < NUM_ZERO) else 1
  155. else:
  156. this_slope = this.slope
  157. that_slope = that.slope
  158. if this_slope != that_slope:
  159. if sweep_line._before:
  160. return -1 if (this_slope > that_slope) else 1
  161. else:
  162. return 1 if (this_slope > that_slope) else -1
  163. delta_x_p1 = this.segment[0][X] - that.segment[0][X]
  164. if delta_x_p1 != NUM_ZERO:
  165. return -1 if (delta_x_p1 < NUM_ZERO) else 1
  166. delta_x_p2 = this.segment[1][X] - that.segment[1][X]
  167. if delta_x_p2 != NUM_ZERO:
  168. return -1 if (delta_x_p2 < NUM_ZERO) else 1
  169. return 0
  170. def __repr__(self):
  171. if self.segment is not None:
  172. return ("Event(0x%x, s0=%r, s1=%r, p=%r, type=%d, slope=%r)" % (
  173. id(self),
  174. self.segment[0], self.segment[1],
  175. self.point,
  176. self.type,
  177. self.slope,
  178. ))
  179. else:
  180. return ("Event(0x%x, s0=%r, s1=%r, p=%r, type=%d, slope=%r)" % (
  181. id(self),
  182. "None", "None",
  183. self.point,
  184. self.type,
  185. self.slope,
  186. ))
  187. class SweepLine:
  188. __slots__ = (
  189. # A map holding all intersection points mapped to the Events
  190. # that form these intersections.
  191. # {Point: set(Event, ...), ...}
  192. "intersections",
  193. "queue",
  194. # Events (sorted set of ordered events, no values)
  195. #
  196. # note: START & END events are considered the same so checking if an event is in the tree
  197. # will return true if its opposite side is found.
  198. # This is essential for the algorithm to work, and why we don't explicitly remove START events.
  199. # Instead, the END events are never added to the current sweep, and removing them also removes the start.
  200. "_events_current_sweep",
  201. # The point of the current Event.
  202. "_current_event_point_x",
  203. # A flag to indicate if we're slightly before or after the line.
  204. "_before",
  205. )
  206. def __init__(self):
  207. self.intersections = {}
  208. self._current_event_point_x = None
  209. self._events_current_sweep = RBTree(cmp=Event.Compare, cmp_data=self)
  210. self._before = True
  211. def get_intersections(self):
  212. """
  213. Return a list of unordered intersection points.
  214. """
  215. if Real is float:
  216. return list(self.intersections.keys())
  217. else:
  218. return [(float(p[0]), float(p[1])) for p in self.intersections.keys()]
  219. # Not essential for implementing this algorithm, but useful.
  220. def get_intersections_with_segments(self):
  221. """
  222. Return a list of unordered intersection '(point, segment)' pairs,
  223. where segments may contain 2 or more values.
  224. """
  225. if Real is float:
  226. return [
  227. (p, [event.segment for event in event_set])
  228. for p, event_set in self.intersections.items()
  229. ]
  230. else:
  231. return [
  232. (
  233. (float(p[0]), float(p[1])),
  234. [((float(event.segment[0][0]), float(event.segment[0][1])),
  235. (float(event.segment[1][0]), float(event.segment[1][1])))
  236. for event in event_set],
  237. )
  238. for p, event_set in self.intersections.items()
  239. ]
  240. # Checks if an intersection exists between two Events 'a' and 'b'.
  241. def _check_intersection(self, a, b):
  242. # Return immediately in case either of the events is null, or
  243. # if one of them is an INTERSECTION event.
  244. if ((a is None or b is None) or
  245. (a.type == Event.Type.INTERSECTION) or
  246. (b.type == Event.Type.INTERSECTION)):
  247. return
  248. if a is b:
  249. return
  250. # Get the intersection point between 'a' and 'b'.
  251. p = isect_seg_seg_v2_point(
  252. a.segment[0], a.segment[1],
  253. b.segment[0], b.segment[1])
  254. # No intersection exists.
  255. if p is None:
  256. return
  257. # If the intersection is formed by both the segment endings, AND
  258. # USE_IGNORE_SEGMENT_ENDINGS is true,
  259. # return from this method.
  260. if USE_IGNORE_SEGMENT_ENDINGS:
  261. if ((len_squared_v2v2(p, a.segment[0]) < NUM_EPS_SQ or
  262. len_squared_v2v2(p, a.segment[1]) < NUM_EPS_SQ) and
  263. (len_squared_v2v2(p, b.segment[0]) < NUM_EPS_SQ or
  264. len_squared_v2v2(p, b.segment[1]) < NUM_EPS_SQ)):
  265. return
  266. # Add the intersection.
  267. events_for_point = self.intersections.pop(p, set())
  268. is_new = len(events_for_point) == 0
  269. events_for_point.add(a)
  270. events_for_point.add(b)
  271. self.intersections[p] = events_for_point
  272. # If the intersection occurs to the right of the sweep line, OR
  273. # if the intersection is on the sweep line and it's above the
  274. # current event-point, add it as a new Event to the queue.
  275. if is_new and p[X] >= self._current_event_point_x:
  276. event_isect = Event(Event.Type.INTERSECTION, p, None, None)
  277. self.queue.offer(p, event_isect)
  278. def _sweep_to(self, p):
  279. if p[X] == self._current_event_point_x:
  280. # happens in rare cases,
  281. # we can safely ignore
  282. return
  283. self._current_event_point_x = p[X]
  284. def insert(self, event):
  285. assert(event not in self._events_current_sweep)
  286. assert(not USE_VERTICAL or event.type != Event.Type.START_VERTICAL)
  287. if USE_DEBUG:
  288. assert(event.in_sweep == False)
  289. assert(event.other.in_sweep == False)
  290. self._events_current_sweep.insert(event, None)
  291. if USE_DEBUG:
  292. event.in_sweep = True
  293. event.other.in_sweep = True
  294. def remove(self, event):
  295. try:
  296. self._events_current_sweep.remove(event)
  297. if USE_DEBUG:
  298. assert(event.in_sweep == True)
  299. assert(event.other.in_sweep == True)
  300. event.in_sweep = False
  301. event.other.in_sweep = False
  302. return True
  303. except KeyError:
  304. if USE_DEBUG:
  305. assert(event.in_sweep == False)
  306. assert(event.other.in_sweep == False)
  307. return False
  308. def above(self, event):
  309. return self._events_current_sweep.succ_key(event, None)
  310. def below(self, event):
  311. return self._events_current_sweep.prev_key(event, None)
  312. '''
  313. def above_all(self, event):
  314. while True:
  315. event = self.above(event)
  316. if event is None:
  317. break
  318. yield event
  319. '''
  320. def above_all(self, event):
  321. # assert(event not in self._events_current_sweep)
  322. return self._events_current_sweep.key_slice(event, None, reverse=False)
  323. def handle(self, p, events_current):
  324. if len(events_current) == 0:
  325. return
  326. # done already
  327. # self._sweep_to(events_current[0])
  328. assert(p[0] == self._current_event_point_x)
  329. if not USE_IGNORE_SEGMENT_ENDINGS:
  330. if len(events_current) > 1:
  331. for i in range(0, len(events_current) - 1):
  332. for j in range(i + 1, len(events_current)):
  333. self._check_intersection(
  334. events_current[i], events_current[j])
  335. for e in events_current:
  336. self.handle_event(e)
  337. def handle_event(self, event):
  338. t = event.type
  339. if t == Event.Type.START:
  340. # print(" START")
  341. self._before = False
  342. self.insert(event)
  343. e_above = self.above(event)
  344. e_below = self.below(event)
  345. self._check_intersection(event, e_above)
  346. self._check_intersection(event, e_below)
  347. if USE_PARANOID:
  348. self._check_intersection(e_above, e_below)
  349. elif t == Event.Type.END:
  350. # print(" END")
  351. self._before = True
  352. e_above = self.above(event)
  353. e_below = self.below(event)
  354. self.remove(event)
  355. self._check_intersection(e_above, e_below)
  356. if USE_PARANOID:
  357. self._check_intersection(event, e_above)
  358. self._check_intersection(event, e_below)
  359. elif t == Event.Type.INTERSECTION:
  360. # print(" INTERSECTION")
  361. self._before = True
  362. event_set = self.intersections[event.point]
  363. # note: events_current aren't sorted.
  364. reinsert_stack = [] # Stack
  365. for e in event_set:
  366. # Since we know the Event wasn't already removed,
  367. # we want to insert it later on.
  368. if self.remove(e):
  369. reinsert_stack.append(e)
  370. self._before = False
  371. # Insert all Events that we were able to remove.
  372. while reinsert_stack:
  373. e = reinsert_stack.pop()
  374. self.insert(e)
  375. e_above = self.above(e)
  376. e_below = self.below(e)
  377. self._check_intersection(e, e_above)
  378. self._check_intersection(e, e_below)
  379. if USE_PARANOID:
  380. self._check_intersection(e_above, e_below)
  381. elif (USE_VERTICAL and
  382. (t == Event.Type.START_VERTICAL)):
  383. # just check sanity
  384. assert(event.segment[0][X] == event.segment[1][X])
  385. assert(event.segment[0][Y] <= event.segment[1][Y])
  386. # In this case we only need to find all segments in this span.
  387. y_above_max = event.segment[1][Y]
  388. # self.insert(event)
  389. for e_above in self.above_all(event):
  390. if e_above.type == Event.Type.START_VERTICAL:
  391. continue
  392. y_above = e_above.y_intercept_x(
  393. self._current_event_point_x)
  394. if USE_IGNORE_SEGMENT_ENDINGS:
  395. if y_above >= y_above_max - NUM_EPS:
  396. break
  397. else:
  398. if y_above > y_above_max:
  399. break
  400. # We know this intersects,
  401. # so we could use a faster function now:
  402. # ix = (self._current_event_point_x, y_above)
  403. # ...however best use existing functions
  404. # since it does all sanity checks on endpoints... etc.
  405. self._check_intersection(event, e_above)
  406. # self.remove(event)
  407. class EventQueue:
  408. __slots__ = (
  409. # note: we only ever pop_min, this could use a 'heap' structure.
  410. # The sorted map holding the points -> event list
  411. # [Point: Event] (tree)
  412. "events_scan",
  413. )
  414. def __init__(self, segments, line):
  415. self.events_scan = RBTree()
  416. # segments = [s for s in segments if s[0][0] != s[1][0] and s[0][1] != s[1][1]]
  417. for s in segments:
  418. assert(s[0][X] <= s[1][X])
  419. slope = slope_v2v2(*s)
  420. if s[0] == s[1]:
  421. pass
  422. elif USE_VERTICAL and (s[0][X] == s[1][X]):
  423. e_start = Event(Event.Type.START_VERTICAL, s[0], s, slope)
  424. if USE_DEBUG:
  425. e_start.other = e_start # FAKE, avoid error checking
  426. self.offer(s[0], e_start)
  427. else:
  428. e_start = Event(Event.Type.START, s[0], s, slope)
  429. e_end = Event(Event.Type.END, s[1], s, slope)
  430. if USE_DEBUG:
  431. e_start.other = e_end
  432. e_end.other = e_start
  433. self.offer(s[0], e_start)
  434. self.offer(s[1], e_end)
  435. line.queue = self
  436. def offer(self, p, e):
  437. """
  438. Offer a new event ``s`` at point ``p`` in this queue.
  439. """
  440. existing = self.events_scan.setdefault(
  441. p, ([], [], [], []) if USE_VERTICAL else
  442. ([], [], []))
  443. # Can use double linked-list for easy insertion at beginning/end
  444. '''
  445. if e.type == Event.Type.END:
  446. existing.insert(0, e)
  447. else:
  448. existing.append(e)
  449. '''
  450. existing[e.type].append(e)
  451. # return a set of events
  452. def poll(self):
  453. """
  454. Get, and remove, the first (lowest) item from this queue.
  455. :return: the first (lowest) item from this queue.
  456. :rtype: Point, Event pair.
  457. """
  458. assert(len(self.events_scan) != 0)
  459. p, events_current = self.events_scan.pop_min()
  460. return p, events_current
  461. def isect_segments_impl(segments, include_segments=False):
  462. # order points left -> right
  463. if Real is float:
  464. segments = [
  465. # in nearly all cases, comparing X is enough,
  466. # but compare Y too for vertical lines
  467. (s[0], s[1]) if (s[0] <= s[1]) else
  468. (s[1], s[0])
  469. for s in segments]
  470. else:
  471. segments = [
  472. # in nearly all cases, comparing X is enough,
  473. # but compare Y too for vertical lines
  474. (
  475. (Real(s[0][0]), Real(s[0][1])),
  476. (Real(s[1][0]), Real(s[1][1])),
  477. ) if (s[0] <= s[1]) else
  478. (
  479. (Real(s[1][0]), Real(s[1][1])),
  480. (Real(s[0][0]), Real(s[0][1])),
  481. )
  482. for s in segments]
  483. sweep_line = SweepLine()
  484. queue = EventQueue(segments, sweep_line)
  485. while len(queue.events_scan) > 0:
  486. if USE_VERBOSE:
  487. print(len(queue.events_scan), sweep_line._current_event_point_x)
  488. p, e_ls = queue.poll()
  489. for events_current in e_ls:
  490. if events_current:
  491. sweep_line._sweep_to(p)
  492. sweep_line.handle(p, events_current)
  493. if include_segments is False:
  494. return sweep_line.get_intersections()
  495. else:
  496. return sweep_line.get_intersections_with_segments()
  497. def isect_polygon_impl(points, include_segments=False):
  498. n = len(points)
  499. segments = [
  500. (tuple(points[i]), tuple(points[(i + 1) % n]))
  501. for i in range(n)]
  502. return isect_segments_impl(segments, include_segments=include_segments)
  503. def isect_segments(segments):
  504. return isect_segments_impl(segments, include_segments=False)
  505. def isect_polygon(segments):
  506. return isect_polygon_impl(segments, include_segments=False)
  507. def isect_segments_include_segments(segments):
  508. return isect_segments_impl(segments, include_segments=True)
  509. def isect_polygon_include_segments(segments):
  510. return isect_polygon_impl(segments, include_segments=True)
  511. # ----------------------------------------------------------------------------
  512. # 2D math utilities
  513. def slope_v2v2(p1, p2):
  514. if p1[X] == p2[X]:
  515. if p1[Y] < p2[Y]:
  516. return NUM_INF
  517. else:
  518. return -NUM_INF
  519. else:
  520. return (p2[Y] - p1[Y]) / (p2[X] - p1[X])
  521. def sub_v2v2(a, b):
  522. return (
  523. a[0] - b[0],
  524. a[1] - b[1])
  525. def dot_v2v2(a, b):
  526. return (
  527. (a[0] * b[0]) +
  528. (a[1] * b[1]))
  529. def len_squared_v2v2(a, b):
  530. c = sub_v2v2(a, b)
  531. return dot_v2v2(c, c)
  532. def line_point_factor_v2(p, l1, l2, default=NUM_ZERO):
  533. u = sub_v2v2(l2, l1)
  534. h = sub_v2v2(p, l1)
  535. dot = dot_v2v2(u, u)
  536. return (dot_v2v2(u, h) / dot) if dot != NUM_ZERO else default
  537. def isect_seg_seg_v2_point(v1, v2, v3, v4, bias=NUM_ZERO):
  538. # Only for predictability and hashable point when same input is given
  539. if v1 > v2:
  540. v1, v2 = v2, v1
  541. if v3 > v4:
  542. v3, v4 = v4, v3
  543. if (v1, v2) > (v3, v4):
  544. v1, v2, v3, v4 = v3, v4, v1, v2
  545. div = (v2[0] - v1[0]) * (v4[1] - v3[1]) - (v2[1] - v1[1]) * (v4[0] - v3[0])
  546. if div == NUM_ZERO:
  547. return None
  548. vi = (((v3[0] - v4[0]) *
  549. (v1[0] * v2[1] - v1[1] * v2[0]) - (v1[0] - v2[0]) *
  550. (v3[0] * v4[1] - v3[1] * v4[0])) / div,
  551. ((v3[1] - v4[1]) *
  552. (v1[0] * v2[1] - v1[1] * v2[0]) - (v1[1] - v2[1]) *
  553. (v3[0] * v4[1] - v3[1] * v4[0])) / div,
  554. )
  555. fac = line_point_factor_v2(vi, v1, v2, default=-NUM_ONE)
  556. if fac < NUM_ZERO - bias or fac > NUM_ONE + bias:
  557. return None
  558. fac = line_point_factor_v2(vi, v3, v4, default=-NUM_ONE)
  559. if fac < NUM_ZERO - bias or fac > NUM_ONE + bias:
  560. return None
  561. # vi = round(vi[X], 8), round(vi[Y], 8)
  562. return vi
  563. # ----------------------------------------------------------------------------
  564. # Simple naive line intersect, (for testing only)
  565. def isect_segments__naive(segments):
  566. """
  567. Brute force O(n2) version of ``isect_segments`` for test validation.
  568. """
  569. isect = []
  570. # order points left -> right
  571. if Real is float:
  572. segments = [
  573. (s[0], s[1]) if s[0][X] <= s[1][X] else
  574. (s[1], s[0])
  575. for s in segments]
  576. else:
  577. segments = [
  578. (
  579. (Real(s[0][0]), Real(s[0][1])),
  580. (Real(s[1][0]), Real(s[1][1])),
  581. ) if (s[0] <= s[1]) else
  582. (
  583. (Real(s[1][0]), Real(s[1][1])),
  584. (Real(s[0][0]), Real(s[0][1])),
  585. )
  586. for s in segments]
  587. n = len(segments)
  588. for i in range(n):
  589. a0, a1 = segments[i]
  590. for j in range(i + 1, n):
  591. b0, b1 = segments[j]
  592. if a0 not in (b0, b1) and a1 not in (b0, b1):
  593. ix = isect_seg_seg_v2_point(a0, a1, b0, b1)
  594. if ix is not None:
  595. # USE_IGNORE_SEGMENT_ENDINGS handled already
  596. isect.append(ix)
  597. return isect
  598. def isect_polygon__naive(points):
  599. """
  600. Brute force O(n2) version of ``isect_polygon`` for test validation.
  601. """
  602. isect = []
  603. n = len(points)
  604. if Real is float:
  605. pass
  606. else:
  607. points = [(Real(p[0]), Real(p[1])) for p in points]
  608. for i in range(n):
  609. a0, a1 = points[i], points[(i + 1) % n]
  610. for j in range(i + 1, n):
  611. b0, b1 = points[j], points[(j + 1) % n]
  612. if a0 not in (b0, b1) and a1 not in (b0, b1):
  613. ix = isect_seg_seg_v2_point(a0, a1, b0, b1)
  614. if ix is not None:
  615. if USE_IGNORE_SEGMENT_ENDINGS:
  616. if ((len_squared_v2v2(ix, a0) < NUM_EPS_SQ or
  617. len_squared_v2v2(ix, a1) < NUM_EPS_SQ) and
  618. (len_squared_v2v2(ix, b0) < NUM_EPS_SQ or
  619. len_squared_v2v2(ix, b1) < NUM_EPS_SQ)):
  620. continue
  621. isect.append(ix)
  622. return isect
  623. # ----------------------------------------------------------------------------
  624. # Inline Libs
  625. #
  626. # bintrees: 2.0.2, extracted from:
  627. # http://pypi.python.org/pypi/bintrees
  628. #
  629. # - Removed unused functions, such as slicing and range iteration.
  630. # - Added 'cmp' and and 'cmp_data' arguments,
  631. # so we can define our own comparison that takes an arg.
  632. # Needed for sweep-line.
  633. # - Added support for 'default' arguments for prev_item/succ_item,
  634. # so we can avoid exception handling.
  635. # -------
  636. # ABCTree
  637. from operator import attrgetter
  638. _sentinel = object()
  639. class _ABCTree(object):
  640. def __init__(self, items=None, cmp=None, cmp_data=None):
  641. """T.__init__(...) initializes T; see T.__class__.__doc__ for signature"""
  642. self._root = None
  643. self._count = 0
  644. if cmp is None:
  645. def cmp(cmp_data, a, b):
  646. if a < b:
  647. return -1
  648. elif a > b:
  649. return 1
  650. else:
  651. return 0
  652. self._cmp = cmp
  653. self._cmp_data = cmp_data
  654. if items is not None:
  655. self.update(items)
  656. def clear(self):
  657. """T.clear() -> None. Remove all items from T."""
  658. def _clear(node):
  659. if node is not None:
  660. _clear(node.left)
  661. _clear(node.right)
  662. node.free()
  663. _clear(self._root)
  664. self._count = 0
  665. self._root = None
  666. @property
  667. def count(self):
  668. """Get items count."""
  669. return self._count
  670. def get_value(self, key):
  671. node = self._root
  672. while node is not None:
  673. cmp = self._cmp(self._cmp_data, key, node.key)
  674. if cmp == 0:
  675. return node.value
  676. elif cmp < 0:
  677. node = node.left
  678. else:
  679. node = node.right
  680. raise KeyError(str(key))
  681. def pop_item(self):
  682. """T.pop_item() -> (k, v), remove and return some (key, value) pair as a
  683. 2-tuple; but raise KeyError if T is empty.
  684. """
  685. if self.is_empty():
  686. raise KeyError("pop_item(): tree is empty")
  687. node = self._root
  688. while True:
  689. if node.left is not None:
  690. node = node.left
  691. elif node.right is not None:
  692. node = node.right
  693. else:
  694. break
  695. key = node.key
  696. value = node.value
  697. self.remove(key)
  698. return key, value
  699. popitem = pop_item # for compatibility to dict()
  700. def min_item(self):
  701. """Get item with min key of tree, raises ValueError if tree is empty."""
  702. if self.is_empty():
  703. raise ValueError("Tree is empty")
  704. node = self._root
  705. while node.left is not None:
  706. node = node.left
  707. return node.key, node.value
  708. def max_item(self):
  709. """Get item with max key of tree, raises ValueError if tree is empty."""
  710. if self.is_empty():
  711. raise ValueError("Tree is empty")
  712. node = self._root
  713. while node.right is not None:
  714. node = node.right
  715. return node.key, node.value
  716. def succ_item(self, key, default=_sentinel):
  717. """Get successor (k,v) pair of key, raises KeyError if key is max key
  718. or key does not exist. optimized for pypy.
  719. """
  720. # removed graingets version, because it was little slower on CPython and much slower on pypy
  721. # this version runs about 4x faster with pypy than the Cython version
  722. # Note: Code sharing of succ_item() and ceiling_item() is possible, but has always a speed penalty.
  723. node = self._root
  724. succ_node = None
  725. while node is not None:
  726. cmp = self._cmp(self._cmp_data, key, node.key)
  727. if cmp == 0:
  728. break
  729. elif cmp < 0:
  730. if (succ_node is None) or self._cmp(self._cmp_data, node.key, succ_node.key) < 0:
  731. succ_node = node
  732. node = node.left
  733. else:
  734. node = node.right
  735. if node is None: # stay at dead end
  736. if default is _sentinel:
  737. raise KeyError(str(key))
  738. return default
  739. # found node of key
  740. if node.right is not None:
  741. # find smallest node of right subtree
  742. node = node.right
  743. while node.left is not None:
  744. node = node.left
  745. if succ_node is None:
  746. succ_node = node
  747. elif self._cmp(self._cmp_data, node.key, succ_node.key) < 0:
  748. succ_node = node
  749. elif succ_node is None: # given key is biggest in tree
  750. if default is _sentinel:
  751. raise KeyError(str(key))
  752. return default
  753. return succ_node.key, succ_node.value
  754. def prev_item(self, key, default=_sentinel):
  755. """Get predecessor (k,v) pair of key, raises KeyError if key is min key
  756. or key does not exist. optimized for pypy.
  757. """
  758. # removed graingets version, because it was little slower on CPython and much slower on pypy
  759. # this version runs about 4x faster with pypy than the Cython version
  760. # Note: Code sharing of prev_item() and floor_item() is possible, but has always a speed penalty.
  761. node = self._root
  762. prev_node = None
  763. while node is not None:
  764. cmp = self._cmp(self._cmp_data, key, node.key)
  765. if cmp == 0:
  766. break
  767. elif cmp < 0:
  768. node = node.left
  769. else:
  770. if (prev_node is None) or self._cmp(self._cmp_data, prev_node.key, node.key) < 0:
  771. prev_node = node
  772. node = node.right
  773. if node is None: # stay at dead end (None)
  774. if default is _sentinel:
  775. raise KeyError(str(key))
  776. return default
  777. # found node of key
  778. if node.left is not None:
  779. # find biggest node of left subtree
  780. node = node.left
  781. while node.right is not None:
  782. node = node.right
  783. if prev_node is None:
  784. prev_node = node
  785. elif self._cmp(self._cmp_data, prev_node.key, node.key) < 0:
  786. prev_node = node
  787. elif prev_node is None: # given key is smallest in tree
  788. if default is _sentinel:
  789. raise KeyError(str(key))
  790. return default
  791. return prev_node.key, prev_node.value
  792. def __repr__(self):
  793. """T.__repr__(...) <==> repr(x)"""
  794. tpl = "%s({%s})" % (self.__class__.__name__, '%s')
  795. return tpl % ", ".join(("%r: %r" % item for item in self.items()))
  796. def __contains__(self, key):
  797. """k in T -> True if T has a key k, else False"""
  798. try:
  799. self.get_value(key)
  800. return True
  801. except KeyError:
  802. return False
  803. def __len__(self):
  804. """T.__len__() <==> len(x)"""
  805. return self.count
  806. def is_empty(self):
  807. """T.is_empty() -> False if T contains any items else True"""
  808. return self.count == 0
  809. def set_default(self, key, default=None):
  810. """T.set_default(k[,d]) -> T.get(k,d), also set T[k]=d if k not in T"""
  811. try:
  812. return self.get_value(key)
  813. except KeyError:
  814. self.insert(key, default)
  815. return default
  816. setdefault = set_default # for compatibility to dict()
  817. def get(self, key, default=None):
  818. """T.get(k[,d]) -> T[k] if k in T, else d. d defaults to None."""
  819. try:
  820. return self.get_value(key)
  821. except KeyError:
  822. return default
  823. def pop(self, key, *args):
  824. """T.pop(k[,d]) -> v, remove specified key and return the corresponding value.
  825. If key is not found, d is returned if given, otherwise KeyError is raised
  826. """
  827. if len(args) > 1:
  828. raise TypeError("pop expected at most 2 arguments, got %d" % (1 + len(args)))
  829. try:
  830. value = self.get_value(key)
  831. self.remove(key)
  832. return value
  833. except KeyError:
  834. if len(args) == 0:
  835. raise
  836. else:
  837. return args[0]
  838. def prev_key(self, key, default=_sentinel):
  839. """Get predecessor to key, raises KeyError if key is min key
  840. or key does not exist.
  841. """
  842. item = self.prev_item(key, default)
  843. return default if item is default else item[0]
  844. def succ_key(self, key, default=_sentinel):
  845. """Get successor to key, raises KeyError if key is max key
  846. or key does not exist.
  847. """
  848. item = self.succ_item(key, default)
  849. return default if item is default else item[0]
  850. def pop_min(self):
  851. """T.pop_min() -> (k, v), remove item with minimum key, raise ValueError
  852. if T is empty.
  853. """
  854. item = self.min_item()
  855. self.remove(item[0])
  856. return item
  857. def pop_max(self):
  858. """T.pop_max() -> (k, v), remove item with maximum key, raise ValueError
  859. if T is empty.
  860. """
  861. item = self.max_item()
  862. self.remove(item[0])
  863. return item
  864. def min_key(self):
  865. """Get min key of tree, raises ValueError if tree is empty. """
  866. return self.min_item()[0]
  867. def max_key(self):
  868. """Get max key of tree, raises ValueError if tree is empty. """
  869. return self.max_item()[0]
  870. def key_slice(self, start_key, end_key, reverse=False):
  871. """T.key_slice(start_key, end_key) -> key iterator:
  872. start_key <= key < end_key.
  873. Yields keys in ascending order if reverse is False else in descending order.
  874. """
  875. return (k for k, v in self.iter_items(start_key, end_key, reverse=reverse))
  876. def iter_items(self, start_key=None, end_key=None, reverse=False):
  877. """Iterates over the (key, value) items of the associated tree,
  878. in ascending order if reverse is True, iterate in descending order,
  879. reverse defaults to False"""
  880. # optimized iterator (reduced method calls) - faster on CPython but slower on pypy
  881. if self.is_empty():
  882. return []
  883. if reverse:
  884. return self._iter_items_backward(start_key, end_key)
  885. else:
  886. return self._iter_items_forward(start_key, end_key)
  887. def _iter_items_forward(self, start_key=None, end_key=None):
  888. for item in self._iter_items(left=attrgetter("left"), right=attrgetter("right"),
  889. start_key=start_key, end_key=end_key):
  890. yield item
  891. def _iter_items_backward(self, start_key=None, end_key=None):
  892. for item in self._iter_items(left=attrgetter("right"), right=attrgetter("left"),
  893. start_key=start_key, end_key=end_key):
  894. yield item
  895. def _iter_items(self, left=attrgetter("left"), right=attrgetter("right"), start_key=None, end_key=None):
  896. node = self._root
  897. stack = []
  898. go_left = True
  899. in_range = self._get_in_range_func(start_key, end_key)
  900. while True:
  901. if left(node) is not None and go_left:
  902. stack.append(node)
  903. node = left(node)
  904. else:
  905. if in_range(node.key):
  906. yield node.key, node.value
  907. if right(node) is not None:
  908. node = right(node)
  909. go_left = True
  910. else:
  911. if not len(stack):
  912. return # all done
  913. node = stack.pop()
  914. go_left = False
  915. def _get_in_range_func(self, start_key, end_key):
  916. if start_key is None and end_key is None:
  917. return lambda x: True
  918. else:
  919. if start_key is None:
  920. start_key = self.min_key()
  921. if end_key is None:
  922. return (lambda x: self._cmp(self._cmp_data, start_key, x) <= 0)
  923. else:
  924. return (lambda x: self._cmp(self._cmp_data, start_key, x) <= 0 and
  925. self._cmp(self._cmp_data, x, end_key) < 0)
  926. # ------
  927. # RBTree
  928. class Node(object):
  929. """Internal object, represents a tree node."""
  930. __slots__ = ['key', 'value', 'red', 'left', 'right']
  931. def __init__(self, key=None, value=None):
  932. self.key = key
  933. self.value = value
  934. self.red = True
  935. self.left = None
  936. self.right = None
  937. def free(self):
  938. self.left = None
  939. self.right = None
  940. self.key = None
  941. self.value = None
  942. def __getitem__(self, key):
  943. """N.__getitem__(key) <==> x[key], where key is 0 (left) or 1 (right)."""
  944. return self.left if key == 0 else self.right
  945. def __setitem__(self, key, value):
  946. """N.__setitem__(key, value) <==> x[key]=value, where key is 0 (left) or 1 (right)."""
  947. if key == 0:
  948. self.left = value
  949. else:
  950. self.right = value
  951. class RBTree(_ABCTree):
  952. """
  953. RBTree implements a balanced binary tree with a dict-like interface.
  954. see: http://en.wikipedia.org/wiki/Red_black_tree
  955. """
  956. @staticmethod
  957. def is_red(node):
  958. if (node is not None) and node.red:
  959. return True
  960. else:
  961. return False
  962. @staticmethod
  963. def jsw_single(root, direction):
  964. other_side = 1 - direction
  965. save = root[other_side]
  966. root[other_side] = save[direction]
  967. save[direction] = root
  968. root.red = True
  969. save.red = False
  970. return save
  971. @staticmethod
  972. def jsw_double(root, direction):
  973. other_side = 1 - direction
  974. root[other_side] = RBTree.jsw_single(root[other_side], other_side)
  975. return RBTree.jsw_single(root, direction)
  976. def _new_node(self, key, value):
  977. """Create a new tree node."""
  978. self._count += 1
  979. return Node(key, value)
  980. def insert(self, key, value):
  981. """T.insert(key, value) <==> T[key] = value, insert key, value into tree."""
  982. if self._root is None: # Empty tree case
  983. self._root = self._new_node(key, value)
  984. self._root.red = False # make root black
  985. return
  986. head = Node() # False tree root
  987. grand_parent = None
  988. grand_grand_parent = head
  989. parent = None # parent
  990. direction = 0
  991. last = 0
  992. # Set up helpers
  993. grand_grand_parent.right = self._root
  994. node = grand_grand_parent.right
  995. # Search down the tree
  996. while True:
  997. if node is None: # Insert new node at the bottom
  998. node = self._new_node(key, value)
  999. parent[direction] = node
  1000. elif RBTree.is_red(node.left) and RBTree.is_red(node.right): # Color flip
  1001. node.red = True
  1002. node.left.red = False
  1003. node.right.red = False
  1004. # Fix red violation
  1005. if RBTree.is_red(node) and RBTree.is_red(parent):
  1006. direction2 = 1 if grand_grand_parent.right is grand_parent else 0
  1007. if node is parent[last]:
  1008. grand_grand_parent[direction2] = RBTree.jsw_single(grand_parent, 1 - last)
  1009. else:
  1010. grand_grand_parent[direction2] = RBTree.jsw_double(grand_parent, 1 - last)
  1011. # Stop if found
  1012. if self._cmp(self._cmp_data, key, node.key) == 0:
  1013. node.value = value # set new value for key
  1014. break
  1015. last = direction
  1016. direction = 0 if (self._cmp(self._cmp_data, key, node.key) < 0) else 1
  1017. # Update helpers
  1018. if grand_parent is not None:
  1019. grand_grand_parent = grand_parent
  1020. grand_parent = parent
  1021. parent = node
  1022. node = node[direction]
  1023. self._root = head.right # Update root
  1024. self._root.red = False # make root black
  1025. def remove(self, key):
  1026. """T.remove(key) <==> del T[key], remove item <key> from tree."""
  1027. if self._root is None:
  1028. raise KeyError(str(key))
  1029. head = Node() # False tree root
  1030. node = head
  1031. node.right = self._root
  1032. parent = None
  1033. grand_parent = None
  1034. found = None # Found item
  1035. direction = 1
  1036. # Search and push a red down
  1037. while node[direction] is not None:
  1038. last = direction
  1039. # Update helpers
  1040. grand_parent = parent
  1041. parent = node
  1042. node = node[direction]
  1043. direction = 1 if (self._cmp(self._cmp_data, node.key, key) < 0) else 0
  1044. # Save found node
  1045. if self._cmp(self._cmp_data, key, node.key) == 0:
  1046. found = node
  1047. # Push the red node down
  1048. if not RBTree.is_red(node) and not RBTree.is_red(node[direction]):
  1049. if RBTree.is_red(node[1 - direction]):
  1050. parent[last] = RBTree.jsw_single(node, direction)
  1051. parent = parent[last]
  1052. elif not RBTree.is_red(node[1 - direction]):
  1053. sibling = parent[1 - last]
  1054. if sibling is not None:
  1055. if (not RBTree.is_red(sibling[1 - last])) and (not RBTree.is_red(sibling[last])):
  1056. # Color flip
  1057. parent.red = False
  1058. sibling.red = True
  1059. node.red = True
  1060. else:
  1061. direction2 = 1 if grand_parent.right is parent else 0
  1062. if RBTree.is_red(sibling[last]):
  1063. grand_parent[direction2] = RBTree.jsw_double(parent, last)
  1064. elif RBTree.is_red(sibling[1-last]):
  1065. grand_parent[direction2] = RBTree.jsw_single(parent, last)
  1066. # Ensure correct coloring
  1067. grand_parent[direction2].red = True
  1068. node.red = True
  1069. grand_parent[direction2].left.red = False
  1070. grand_parent[direction2].right.red = False
  1071. # Replace and remove if found
  1072. if found is not None:
  1073. found.key = node.key
  1074. found.value = node.value
  1075. parent[int(parent.right is node)] = node[int(node.left is None)]
  1076. node.free()
  1077. self._count -= 1
  1078. # Update root and make it black
  1079. self._root = head.right
  1080. if self._root is not None:
  1081. self._root.red = False
  1082. if not found:
  1083. raise KeyError(str(key))