test_polygon.py 27 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676
  1. from sympy.core.numbers import (Float, Rational, oo, pi)
  2. from sympy.core.singleton import S
  3. from sympy.core.symbol import (Symbol, symbols)
  4. from sympy.functions.elementary.complexes import Abs
  5. from sympy.functions.elementary.miscellaneous import sqrt
  6. from sympy.functions.elementary.trigonometric import (acos, cos, sin)
  7. from sympy.functions.elementary.trigonometric import tan
  8. from sympy.geometry import (Circle, Ellipse, GeometryError, Point, Point2D,
  9. Polygon, Ray, RegularPolygon, Segment, Triangle,
  10. are_similar, convex_hull, intersection, Line, Ray2D)
  11. from sympy.testing.pytest import raises, slow, warns
  12. from sympy.core.random import verify_numerically
  13. from sympy.geometry.polygon import rad, deg
  14. from sympy.integrals.integrals import integrate
  15. from sympy.utilities.iterables import rotate_left
  16. def feq(a, b):
  17. """Test if two floating point values are 'equal'."""
  18. t_float = Float("1.0E-10")
  19. return -t_float < a - b < t_float
  20. @slow
  21. def test_polygon():
  22. x = Symbol('x', real=True)
  23. y = Symbol('y', real=True)
  24. q = Symbol('q', real=True)
  25. u = Symbol('u', real=True)
  26. v = Symbol('v', real=True)
  27. w = Symbol('w', real=True)
  28. x1 = Symbol('x1', real=True)
  29. half = S.Half
  30. a, b, c = Point(0, 0), Point(2, 0), Point(3, 3)
  31. t = Triangle(a, b, c)
  32. assert Polygon(Point(0, 0)) == Point(0, 0)
  33. assert Polygon(a, Point(1, 0), b, c) == t
  34. assert Polygon(Point(1, 0), b, c, a) == t
  35. assert Polygon(b, c, a, Point(1, 0)) == t
  36. # 2 "remove folded" tests
  37. assert Polygon(a, Point(3, 0), b, c) == t
  38. assert Polygon(a, b, Point(3, -1), b, c) == t
  39. # remove multiple collinear points
  40. assert Polygon(Point(-4, 15), Point(-11, 15), Point(-15, 15),
  41. Point(-15, 33/5), Point(-15, -87/10), Point(-15, -15),
  42. Point(-42/5, -15), Point(-2, -15), Point(7, -15), Point(15, -15),
  43. Point(15, -3), Point(15, 10), Point(15, 15)) == \
  44. Polygon(Point(-15, -15), Point(15, -15), Point(15, 15), Point(-15, 15))
  45. p1 = Polygon(
  46. Point(0, 0), Point(3, -1),
  47. Point(6, 0), Point(4, 5),
  48. Point(2, 3), Point(0, 3))
  49. p2 = Polygon(
  50. Point(6, 0), Point(3, -1),
  51. Point(0, 0), Point(0, 3),
  52. Point(2, 3), Point(4, 5))
  53. p3 = Polygon(
  54. Point(0, 0), Point(3, 0),
  55. Point(5, 2), Point(4, 4))
  56. p4 = Polygon(
  57. Point(0, 0), Point(4, 4),
  58. Point(5, 2), Point(3, 0))
  59. p5 = Polygon(
  60. Point(0, 0), Point(4, 4),
  61. Point(0, 4))
  62. p6 = Polygon(
  63. Point(-11, 1), Point(-9, 6.6),
  64. Point(-4, -3), Point(-8.4, -8.7))
  65. p7 = Polygon(
  66. Point(x, y), Point(q, u),
  67. Point(v, w))
  68. p8 = Polygon(
  69. Point(x, y), Point(v, w),
  70. Point(q, u))
  71. p9 = Polygon(
  72. Point(0, 0), Point(4, 4),
  73. Point(3, 0), Point(5, 2))
  74. p10 = Polygon(
  75. Point(0, 2), Point(2, 2),
  76. Point(0, 0), Point(2, 0))
  77. p11 = Polygon(Point(0, 0), 1, n=3)
  78. p12 = Polygon(Point(0, 0), 1, 0, n=3)
  79. p13 = Polygon(
  80. Point(0, 0),Point(8, 8),
  81. Point(23, 20),Point(0, 20))
  82. p14 = Polygon(*rotate_left(p13.args, 1))
  83. r = Ray(Point(-9, 6.6), Point(-9, 5.5))
  84. #
  85. # General polygon
  86. #
  87. assert p1 == p2
  88. assert len(p1.args) == 6
  89. assert len(p1.sides) == 6
  90. assert p1.perimeter == 5 + 2*sqrt(10) + sqrt(29) + sqrt(8)
  91. assert p1.area == 22
  92. assert not p1.is_convex()
  93. assert Polygon((-1, 1), (2, -1), (2, 1), (-1, -1), (3, 0)
  94. ).is_convex() is False
  95. # ensure convex for both CW and CCW point specification
  96. assert p3.is_convex()
  97. assert p4.is_convex()
  98. dict5 = p5.angles
  99. assert dict5[Point(0, 0)] == pi / 4
  100. assert dict5[Point(0, 4)] == pi / 2
  101. assert p5.encloses_point(Point(x, y)) is None
  102. assert p5.encloses_point(Point(1, 3))
  103. assert p5.encloses_point(Point(0, 0)) is False
  104. assert p5.encloses_point(Point(4, 0)) is False
  105. assert p1.encloses(Circle(Point(2.5, 2.5), 5)) is False
  106. assert p1.encloses(Ellipse(Point(2.5, 2), 5, 6)) is False
  107. assert p5.plot_interval('x') == [x, 0, 1]
  108. assert p5.distance(
  109. Polygon(Point(10, 10), Point(14, 14), Point(10, 14))) == 6 * sqrt(2)
  110. assert p5.distance(
  111. Polygon(Point(1, 8), Point(5, 8), Point(8, 12), Point(1, 12))) == 4
  112. with warns(UserWarning, \
  113. match="Polygons may intersect producing erroneous output"):
  114. Polygon(Point(0, 0), Point(1, 0), Point(1, 1)).distance(
  115. Polygon(Point(0, 0), Point(0, 1), Point(1, 1)))
  116. assert hash(p5) == hash(Polygon(Point(0, 0), Point(4, 4), Point(0, 4)))
  117. assert hash(p1) == hash(p2)
  118. assert hash(p7) == hash(p8)
  119. assert hash(p3) != hash(p9)
  120. assert p5 == Polygon(Point(4, 4), Point(0, 4), Point(0, 0))
  121. assert Polygon(Point(4, 4), Point(0, 4), Point(0, 0)) in p5
  122. assert p5 != Point(0, 4)
  123. assert Point(0, 1) in p5
  124. assert p5.arbitrary_point('t').subs(Symbol('t', real=True), 0) == \
  125. Point(0, 0)
  126. raises(ValueError, lambda: Polygon(
  127. Point(x, 0), Point(0, y), Point(x, y)).arbitrary_point('x'))
  128. assert p6.intersection(r) == [Point(-9, Rational(-84, 13)), Point(-9, Rational(33, 5))]
  129. assert p10.area == 0
  130. assert p11 == RegularPolygon(Point(0, 0), 1, 3, 0)
  131. assert p11 == p12
  132. assert p11.vertices[0] == Point(1, 0)
  133. assert p11.args[0] == Point(0, 0)
  134. p11.spin(pi/2)
  135. assert p11.vertices[0] == Point(0, 1)
  136. #
  137. # Regular polygon
  138. #
  139. p1 = RegularPolygon(Point(0, 0), 10, 5)
  140. p2 = RegularPolygon(Point(0, 0), 5, 5)
  141. raises(GeometryError, lambda: RegularPolygon(Point(0, 0), Point(0,
  142. 1), Point(1, 1)))
  143. raises(GeometryError, lambda: RegularPolygon(Point(0, 0), 1, 2))
  144. raises(ValueError, lambda: RegularPolygon(Point(0, 0), 1, 2.5))
  145. assert p1 != p2
  146. assert p1.interior_angle == pi*Rational(3, 5)
  147. assert p1.exterior_angle == pi*Rational(2, 5)
  148. assert p2.apothem == 5*cos(pi/5)
  149. assert p2.circumcenter == p1.circumcenter == Point(0, 0)
  150. assert p1.circumradius == p1.radius == 10
  151. assert p2.circumcircle == Circle(Point(0, 0), 5)
  152. assert p2.incircle == Circle(Point(0, 0), p2.apothem)
  153. assert p2.inradius == p2.apothem == (5 * (1 + sqrt(5)) / 4)
  154. p2.spin(pi / 10)
  155. dict1 = p2.angles
  156. assert dict1[Point(0, 5)] == 3 * pi / 5
  157. assert p1.is_convex()
  158. assert p1.rotation == 0
  159. assert p1.encloses_point(Point(0, 0))
  160. assert p1.encloses_point(Point(11, 0)) is False
  161. assert p2.encloses_point(Point(0, 4.9))
  162. p1.spin(pi/3)
  163. assert p1.rotation == pi/3
  164. assert p1.vertices[0] == Point(5, 5*sqrt(3))
  165. for var in p1.args:
  166. if isinstance(var, Point):
  167. assert var == Point(0, 0)
  168. else:
  169. assert var in (5, 10, pi / 3)
  170. assert p1 != Point(0, 0)
  171. assert p1 != p5
  172. # while spin works in place (notice that rotation is 2pi/3 below)
  173. # rotate returns a new object
  174. p1_old = p1
  175. assert p1.rotate(pi/3) == RegularPolygon(Point(0, 0), 10, 5, pi*Rational(2, 3))
  176. assert p1 == p1_old
  177. assert p1.area == (-250*sqrt(5) + 1250)/(4*tan(pi/5))
  178. assert p1.length == 20*sqrt(-sqrt(5)/8 + Rational(5, 8))
  179. assert p1.scale(2, 2) == \
  180. RegularPolygon(p1.center, p1.radius*2, p1._n, p1.rotation)
  181. assert RegularPolygon((0, 0), 1, 4).scale(2, 3) == \
  182. Polygon(Point(2, 0), Point(0, 3), Point(-2, 0), Point(0, -3))
  183. assert repr(p1) == str(p1)
  184. #
  185. # Angles
  186. #
  187. angles = p4.angles
  188. assert feq(angles[Point(0, 0)].evalf(), Float("0.7853981633974483"))
  189. assert feq(angles[Point(4, 4)].evalf(), Float("1.2490457723982544"))
  190. assert feq(angles[Point(5, 2)].evalf(), Float("1.8925468811915388"))
  191. assert feq(angles[Point(3, 0)].evalf(), Float("2.3561944901923449"))
  192. angles = p3.angles
  193. assert feq(angles[Point(0, 0)].evalf(), Float("0.7853981633974483"))
  194. assert feq(angles[Point(4, 4)].evalf(), Float("1.2490457723982544"))
  195. assert feq(angles[Point(5, 2)].evalf(), Float("1.8925468811915388"))
  196. assert feq(angles[Point(3, 0)].evalf(), Float("2.3561944901923449"))
  197. # https://github.com/sympy/sympy/issues/24885
  198. interior_angles_sum = sum(p13.angles.values())
  199. assert feq(interior_angles_sum, (len(p13.angles) - 2)*pi )
  200. interior_angles_sum = sum(p14.angles.values())
  201. assert feq(interior_angles_sum, (len(p14.angles) - 2)*pi )
  202. #
  203. # Triangle
  204. #
  205. p1 = Point(0, 0)
  206. p2 = Point(5, 0)
  207. p3 = Point(0, 5)
  208. t1 = Triangle(p1, p2, p3)
  209. t2 = Triangle(p1, p2, Point(Rational(5, 2), sqrt(Rational(75, 4))))
  210. t3 = Triangle(p1, Point(x1, 0), Point(0, x1))
  211. s1 = t1.sides
  212. assert Triangle(p1, p2, p1) == Polygon(p1, p2, p1) == Segment(p1, p2)
  213. raises(GeometryError, lambda: Triangle(Point(0, 0)))
  214. # Basic stuff
  215. assert Triangle(p1, p1, p1) == p1
  216. assert Triangle(p2, p2*2, p2*3) == Segment(p2, p2*3)
  217. assert t1.area == Rational(25, 2)
  218. assert t1.is_right()
  219. assert t2.is_right() is False
  220. assert t3.is_right()
  221. assert p1 in t1
  222. assert t1.sides[0] in t1
  223. assert Segment((0, 0), (1, 0)) in t1
  224. assert Point(5, 5) not in t2
  225. assert t1.is_convex()
  226. assert feq(t1.angles[p1].evalf(), pi.evalf()/2)
  227. assert t1.is_equilateral() is False
  228. assert t2.is_equilateral()
  229. assert t3.is_equilateral() is False
  230. assert are_similar(t1, t2) is False
  231. assert are_similar(t1, t3)
  232. assert are_similar(t2, t3) is False
  233. assert t1.is_similar(Point(0, 0)) is False
  234. assert t1.is_similar(t2) is False
  235. # Bisectors
  236. bisectors = t1.bisectors()
  237. assert bisectors[p1] == Segment(
  238. p1, Point(Rational(5, 2), Rational(5, 2)))
  239. assert t2.bisectors()[p2] == Segment(
  240. Point(5, 0), Point(Rational(5, 4), 5*sqrt(3)/4))
  241. p4 = Point(0, x1)
  242. assert t3.bisectors()[p4] == Segment(p4, Point(x1*(sqrt(2) - 1), 0))
  243. ic = (250 - 125*sqrt(2))/50
  244. assert t1.incenter == Point(ic, ic)
  245. # Inradius
  246. assert t1.inradius == t1.incircle.radius == 5 - 5*sqrt(2)/2
  247. assert t2.inradius == t2.incircle.radius == 5*sqrt(3)/6
  248. assert t3.inradius == t3.incircle.radius == x1**2/((2 + sqrt(2))*Abs(x1))
  249. # Exradius
  250. assert t1.exradii[t1.sides[2]] == 5*sqrt(2)/2
  251. # Excenters
  252. assert t1.excenters[t1.sides[2]] == Point2D(25*sqrt(2), -5*sqrt(2)/2)
  253. # Circumcircle
  254. assert t1.circumcircle.center == Point(2.5, 2.5)
  255. # Medians + Centroid
  256. m = t1.medians
  257. assert t1.centroid == Point(Rational(5, 3), Rational(5, 3))
  258. assert m[p1] == Segment(p1, Point(Rational(5, 2), Rational(5, 2)))
  259. assert t3.medians[p1] == Segment(p1, Point(x1/2, x1/2))
  260. assert intersection(m[p1], m[p2], m[p3]) == [t1.centroid]
  261. assert t1.medial == Triangle(Point(2.5, 0), Point(0, 2.5), Point(2.5, 2.5))
  262. # Nine-point circle
  263. assert t1.nine_point_circle == Circle(Point(2.5, 0),
  264. Point(0, 2.5), Point(2.5, 2.5))
  265. assert t1.nine_point_circle == Circle(Point(0, 0),
  266. Point(0, 2.5), Point(2.5, 2.5))
  267. # Perpendicular
  268. altitudes = t1.altitudes
  269. assert altitudes[p1] == Segment(p1, Point(Rational(5, 2), Rational(5, 2)))
  270. assert altitudes[p2].equals(s1[0])
  271. assert altitudes[p3] == s1[2]
  272. assert t1.orthocenter == p1
  273. t = S('''Triangle(
  274. Point(100080156402737/5000000000000, 79782624633431/500000000000),
  275. Point(39223884078253/2000000000000, 156345163124289/1000000000000),
  276. Point(31241359188437/1250000000000, 338338270939941/1000000000000000))''')
  277. assert t.orthocenter == S('''Point(-780660869050599840216997'''
  278. '''79471538701955848721853/80368430960602242240789074233100000000000000,'''
  279. '''20151573611150265741278060334545897615974257/16073686192120448448157'''
  280. '''8148466200000000000)''')
  281. # Ensure
  282. assert len(intersection(*bisectors.values())) == 1
  283. assert len(intersection(*altitudes.values())) == 1
  284. assert len(intersection(*m.values())) == 1
  285. # Distance
  286. p1 = Polygon(
  287. Point(0, 0), Point(1, 0),
  288. Point(1, 1), Point(0, 1))
  289. p2 = Polygon(
  290. Point(0, Rational(5)/4), Point(1, Rational(5)/4),
  291. Point(1, Rational(9)/4), Point(0, Rational(9)/4))
  292. p3 = Polygon(
  293. Point(1, 2), Point(2, 2),
  294. Point(2, 1))
  295. p4 = Polygon(
  296. Point(1, 1), Point(Rational(6)/5, 1),
  297. Point(1, Rational(6)/5))
  298. pt1 = Point(half, half)
  299. pt2 = Point(1, 1)
  300. '''Polygon to Point'''
  301. assert p1.distance(pt1) == half
  302. assert p1.distance(pt2) == 0
  303. assert p2.distance(pt1) == Rational(3)/4
  304. assert p3.distance(pt2) == sqrt(2)/2
  305. '''Polygon to Polygon'''
  306. # p1.distance(p2) emits a warning
  307. with warns(UserWarning, \
  308. match="Polygons may intersect producing erroneous output"):
  309. assert p1.distance(p2) == half/2
  310. assert p1.distance(p3) == sqrt(2)/2
  311. # p3.distance(p4) emits a warning
  312. with warns(UserWarning, \
  313. match="Polygons may intersect producing erroneous output"):
  314. assert p3.distance(p4) == (sqrt(2)/2 - sqrt(Rational(2)/25)/2)
  315. def test_convex_hull():
  316. p = [Point(-5, -1), Point(-2, 1), Point(-2, -1), Point(-1, -3), \
  317. Point(0, 0), Point(1, 1), Point(2, 2), Point(2, -1), Point(3, 1), \
  318. Point(4, -1), Point(6, 2)]
  319. ch = Polygon(p[0], p[3], p[9], p[10], p[6], p[1])
  320. #test handling of duplicate points
  321. p.append(p[3])
  322. #more than 3 collinear points
  323. another_p = [Point(-45, -85), Point(-45, 85), Point(-45, 26), \
  324. Point(-45, -24)]
  325. ch2 = Segment(another_p[0], another_p[1])
  326. assert convex_hull(*another_p) == ch2
  327. assert convex_hull(*p) == ch
  328. assert convex_hull(p[0]) == p[0]
  329. assert convex_hull(p[0], p[1]) == Segment(p[0], p[1])
  330. # no unique points
  331. assert convex_hull(*[p[-1]]*3) == p[-1]
  332. # collection of items
  333. assert convex_hull(*[Point(0, 0), \
  334. Segment(Point(1, 0), Point(1, 1)), \
  335. RegularPolygon(Point(2, 0), 2, 4)]) == \
  336. Polygon(Point(0, 0), Point(2, -2), Point(4, 0), Point(2, 2))
  337. def test_encloses():
  338. # square with a dimpled left side
  339. s = Polygon(Point(0, 0), Point(1, 0), Point(1, 1), Point(0, 1), \
  340. Point(S.Half, S.Half))
  341. # the following is True if the polygon isn't treated as closing on itself
  342. assert s.encloses(Point(0, S.Half)) is False
  343. assert s.encloses(Point(S.Half, S.Half)) is False # it's a vertex
  344. assert s.encloses(Point(Rational(3, 4), S.Half)) is True
  345. def test_triangle_kwargs():
  346. assert Triangle(sss=(3, 4, 5)) == \
  347. Triangle(Point(0, 0), Point(3, 0), Point(3, 4))
  348. assert Triangle(asa=(30, 2, 30)) == \
  349. Triangle(Point(0, 0), Point(2, 0), Point(1, sqrt(3)/3))
  350. assert Triangle(sas=(1, 45, 2)) == \
  351. Triangle(Point(0, 0), Point(2, 0), Point(sqrt(2)/2, sqrt(2)/2))
  352. assert Triangle(sss=(1, 2, 5)) is None
  353. assert deg(rad(180)) == 180
  354. def test_transform():
  355. pts = [Point(0, 0), Point(S.Half, Rational(1, 4)), Point(1, 1)]
  356. pts_out = [Point(-4, -10), Point(-3, Rational(-37, 4)), Point(-2, -7)]
  357. assert Triangle(*pts).scale(2, 3, (4, 5)) == Triangle(*pts_out)
  358. assert RegularPolygon((0, 0), 1, 4).scale(2, 3, (4, 5)) == \
  359. Polygon(Point(-2, -10), Point(-4, -7), Point(-6, -10), Point(-4, -13))
  360. # Checks for symmetric scaling
  361. assert RegularPolygon((0, 0), 1, 4).scale(2, 2) == \
  362. RegularPolygon(Point2D(0, 0), 2, 4, 0)
  363. def test_reflect():
  364. x = Symbol('x', real=True)
  365. y = Symbol('y', real=True)
  366. b = Symbol('b')
  367. m = Symbol('m')
  368. l = Line((0, b), slope=m)
  369. p = Point(x, y)
  370. r = p.reflect(l)
  371. dp = l.perpendicular_segment(p).length
  372. dr = l.perpendicular_segment(r).length
  373. assert verify_numerically(dp, dr)
  374. assert Polygon((1, 0), (2, 0), (2, 2)).reflect(Line((3, 0), slope=oo)) \
  375. == Triangle(Point(5, 0), Point(4, 0), Point(4, 2))
  376. assert Polygon((1, 0), (2, 0), (2, 2)).reflect(Line((0, 3), slope=oo)) \
  377. == Triangle(Point(-1, 0), Point(-2, 0), Point(-2, 2))
  378. assert Polygon((1, 0), (2, 0), (2, 2)).reflect(Line((0, 3), slope=0)) \
  379. == Triangle(Point(1, 6), Point(2, 6), Point(2, 4))
  380. assert Polygon((1, 0), (2, 0), (2, 2)).reflect(Line((3, 0), slope=0)) \
  381. == Triangle(Point(1, 0), Point(2, 0), Point(2, -2))
  382. def test_bisectors():
  383. p1, p2, p3 = Point(0, 0), Point(1, 0), Point(0, 1)
  384. p = Polygon(Point(0, 0), Point(2, 0), Point(1, 1), Point(0, 3))
  385. q = Polygon(Point(1, 0), Point(2, 0), Point(3, 3), Point(-1, 5))
  386. poly = Polygon(Point(3, 4), Point(0, 0), Point(8, 7), Point(-1, 1), Point(19, -19))
  387. t = Triangle(p1, p2, p3)
  388. assert t.bisectors()[p2] == Segment(Point(1, 0), Point(0, sqrt(2) - 1))
  389. assert p.bisectors()[Point2D(0, 3)] == Ray2D(Point2D(0, 3), \
  390. Point2D(sin(acos(2*sqrt(5)/5)/2), 3 - cos(acos(2*sqrt(5)/5)/2)))
  391. assert q.bisectors()[Point2D(-1, 5)] == \
  392. Ray2D(Point2D(-1, 5), Point2D(-1 + sqrt(29)*(5*sin(acos(9*sqrt(145)/145)/2) + \
  393. 2*cos(acos(9*sqrt(145)/145)/2))/29, sqrt(29)*(-5*cos(acos(9*sqrt(145)/145)/2) + \
  394. 2*sin(acos(9*sqrt(145)/145)/2))/29 + 5))
  395. assert poly.bisectors()[Point2D(-1, 1)] == Ray2D(Point2D(-1, 1), \
  396. Point2D(-1 + sin(acos(sqrt(26)/26)/2 + pi/4), 1 - sin(-acos(sqrt(26)/26)/2 + pi/4)))
  397. def test_incenter():
  398. assert Triangle(Point(0, 0), Point(1, 0), Point(0, 1)).incenter \
  399. == Point(1 - sqrt(2)/2, 1 - sqrt(2)/2)
  400. def test_inradius():
  401. assert Triangle(Point(0, 0), Point(4, 0), Point(0, 3)).inradius == 1
  402. def test_incircle():
  403. assert Triangle(Point(0, 0), Point(2, 0), Point(0, 2)).incircle \
  404. == Circle(Point(2 - sqrt(2), 2 - sqrt(2)), 2 - sqrt(2))
  405. def test_exradii():
  406. t = Triangle(Point(0, 0), Point(6, 0), Point(0, 2))
  407. assert t.exradii[t.sides[2]] == (-2 + sqrt(10))
  408. def test_medians():
  409. t = Triangle(Point(0, 0), Point(1, 0), Point(0, 1))
  410. assert t.medians[Point(0, 0)] == Segment(Point(0, 0), Point(S.Half, S.Half))
  411. def test_medial():
  412. assert Triangle(Point(0, 0), Point(1, 0), Point(0, 1)).medial \
  413. == Triangle(Point(S.Half, 0), Point(S.Half, S.Half), Point(0, S.Half))
  414. def test_nine_point_circle():
  415. assert Triangle(Point(0, 0), Point(1, 0), Point(0, 1)).nine_point_circle \
  416. == Circle(Point2D(Rational(1, 4), Rational(1, 4)), sqrt(2)/4)
  417. def test_eulerline():
  418. assert Triangle(Point(0, 0), Point(1, 0), Point(0, 1)).eulerline \
  419. == Line(Point2D(0, 0), Point2D(S.Half, S.Half))
  420. assert Triangle(Point(0, 0), Point(10, 0), Point(5, 5*sqrt(3))).eulerline \
  421. == Point2D(5, 5*sqrt(3)/3)
  422. assert Triangle(Point(4, -6), Point(4, -1), Point(-3, 3)).eulerline \
  423. == Line(Point2D(Rational(64, 7), 3), Point2D(Rational(-29, 14), Rational(-7, 2)))
  424. def test_intersection():
  425. poly1 = Triangle(Point(0, 0), Point(1, 0), Point(0, 1))
  426. poly2 = Polygon(Point(0, 1), Point(-5, 0),
  427. Point(0, -4), Point(0, Rational(1, 5)),
  428. Point(S.Half, -0.1), Point(1, 0), Point(0, 1))
  429. assert poly1.intersection(poly2) == [Point2D(Rational(1, 3), 0),
  430. Segment(Point(0, Rational(1, 5)), Point(0, 0)),
  431. Segment(Point(1, 0), Point(0, 1))]
  432. assert poly2.intersection(poly1) == [Point(Rational(1, 3), 0),
  433. Segment(Point(0, 0), Point(0, Rational(1, 5))),
  434. Segment(Point(1, 0), Point(0, 1))]
  435. assert poly1.intersection(Point(0, 0)) == [Point(0, 0)]
  436. assert poly1.intersection(Point(-12, -43)) == []
  437. assert poly2.intersection(Line((-12, 0), (12, 0))) == [Point(-5, 0),
  438. Point(0, 0), Point(Rational(1, 3), 0), Point(1, 0)]
  439. assert poly2.intersection(Line((-12, 12), (12, 12))) == []
  440. assert poly2.intersection(Ray((-3, 4), (1, 0))) == [Segment(Point(1, 0),
  441. Point(0, 1))]
  442. assert poly2.intersection(Circle((0, -1), 1)) == [Point(0, -2),
  443. Point(0, 0)]
  444. assert poly1.intersection(poly1) == [Segment(Point(0, 0), Point(1, 0)),
  445. Segment(Point(0, 1), Point(0, 0)), Segment(Point(1, 0), Point(0, 1))]
  446. assert poly2.intersection(poly2) == [Segment(Point(-5, 0), Point(0, -4)),
  447. Segment(Point(0, -4), Point(0, Rational(1, 5))),
  448. Segment(Point(0, Rational(1, 5)), Point(S.Half, Rational(-1, 10))),
  449. Segment(Point(0, 1), Point(-5, 0)),
  450. Segment(Point(S.Half, Rational(-1, 10)), Point(1, 0)),
  451. Segment(Point(1, 0), Point(0, 1))]
  452. assert poly2.intersection(Triangle(Point(0, 1), Point(1, 0), Point(-1, 1))) \
  453. == [Point(Rational(-5, 7), Rational(6, 7)), Segment(Point2D(0, 1), Point(1, 0))]
  454. assert poly1.intersection(RegularPolygon((-12, -15), 3, 3)) == []
  455. def test_parameter_value():
  456. t = Symbol('t')
  457. sq = Polygon((0, 0), (0, 1), (1, 1), (1, 0))
  458. assert sq.parameter_value((0.5, 1), t) == {t: Rational(3, 8)}
  459. q = Polygon((0, 0), (2, 1), (2, 4), (4, 0))
  460. assert q.parameter_value((4, 0), t) == {t: -6 + 3*sqrt(5)} # ~= 0.708
  461. raises(ValueError, lambda: sq.parameter_value((5, 6), t))
  462. raises(ValueError, lambda: sq.parameter_value(Circle(Point(0, 0), 1), t))
  463. def test_issue_12966():
  464. poly = Polygon(Point(0, 0), Point(0, 10), Point(5, 10), Point(5, 5),
  465. Point(10, 5), Point(10, 0))
  466. t = Symbol('t')
  467. pt = poly.arbitrary_point(t)
  468. DELTA = 5/poly.perimeter
  469. assert [pt.subs(t, DELTA*i) for i in range(int(1/DELTA))] == [
  470. Point(0, 0), Point(0, 5), Point(0, 10), Point(5, 10),
  471. Point(5, 5), Point(10, 5), Point(10, 0), Point(5, 0)]
  472. def test_second_moment_of_area():
  473. x, y = symbols('x, y')
  474. # triangle
  475. p1, p2, p3 = [(0, 0), (4, 0), (0, 2)]
  476. p = (0, 0)
  477. # equation of hypotenuse
  478. eq_y = (1-x/4)*2
  479. I_yy = integrate((x**2) * (integrate(1, (y, 0, eq_y))), (x, 0, 4))
  480. I_xx = integrate(1 * (integrate(y**2, (y, 0, eq_y))), (x, 0, 4))
  481. I_xy = integrate(x * (integrate(y, (y, 0, eq_y))), (x, 0, 4))
  482. triangle = Polygon(p1, p2, p3)
  483. assert (I_xx - triangle.second_moment_of_area(p)[0]) == 0
  484. assert (I_yy - triangle.second_moment_of_area(p)[1]) == 0
  485. assert (I_xy - triangle.second_moment_of_area(p)[2]) == 0
  486. # rectangle
  487. p1, p2, p3, p4=[(0, 0), (4, 0), (4, 2), (0, 2)]
  488. I_yy = integrate((x**2) * integrate(1, (y, 0, 2)), (x, 0, 4))
  489. I_xx = integrate(1 * integrate(y**2, (y, 0, 2)), (x, 0, 4))
  490. I_xy = integrate(x * integrate(y, (y, 0, 2)), (x, 0, 4))
  491. rectangle = Polygon(p1, p2, p3, p4)
  492. assert (I_xx - rectangle.second_moment_of_area(p)[0]) == 0
  493. assert (I_yy - rectangle.second_moment_of_area(p)[1]) == 0
  494. assert (I_xy - rectangle.second_moment_of_area(p)[2]) == 0
  495. r = RegularPolygon(Point(0, 0), 5, 3)
  496. assert r.second_moment_of_area() == (1875*sqrt(3)/S(32), 1875*sqrt(3)/S(32), 0)
  497. def test_first_moment():
  498. a, b = symbols('a, b', positive=True)
  499. # rectangle
  500. p1 = Polygon((0, 0), (a, 0), (a, b), (0, b))
  501. assert p1.first_moment_of_area() == (a*b**2/8, a**2*b/8)
  502. assert p1.first_moment_of_area((a/3, b/4)) == (-3*a*b**2/32, -a**2*b/9)
  503. p1 = Polygon((0, 0), (40, 0), (40, 30), (0, 30))
  504. assert p1.first_moment_of_area() == (4500, 6000)
  505. # triangle
  506. p2 = Polygon((0, 0), (a, 0), (a/2, b))
  507. assert p2.first_moment_of_area() == (4*a*b**2/81, a**2*b/24)
  508. assert p2.first_moment_of_area((a/8, b/6)) == (-25*a*b**2/648, -5*a**2*b/768)
  509. p2 = Polygon((0, 0), (12, 0), (12, 30))
  510. assert p2.first_moment_of_area() == (S(1600)/3, -S(640)/3)
  511. def test_section_modulus_and_polar_second_moment_of_area():
  512. a, b = symbols('a, b', positive=True)
  513. x, y = symbols('x, y')
  514. rectangle = Polygon((0, b), (0, 0), (a, 0), (a, b))
  515. assert rectangle.section_modulus(Point(x, y)) == (a*b**3/12/(-b/2 + y), a**3*b/12/(-a/2 + x))
  516. assert rectangle.polar_second_moment_of_area() == a**3*b/12 + a*b**3/12
  517. convex = RegularPolygon((0, 0), 1, 6)
  518. assert convex.section_modulus() == (Rational(5, 8), sqrt(3)*Rational(5, 16))
  519. assert convex.polar_second_moment_of_area() == 5*sqrt(3)/S(8)
  520. concave = Polygon((0, 0), (1, 8), (3, 4), (4, 6), (7, 1))
  521. assert concave.section_modulus() == (Rational(-6371, 429), Rational(-9778, 519))
  522. assert concave.polar_second_moment_of_area() == Rational(-38669, 252)
  523. def test_cut_section():
  524. # concave polygon
  525. p = Polygon((-1, -1), (1, Rational(5, 2)), (2, 1), (3, Rational(5, 2)), (4, 2), (5, 3), (-1, 3))
  526. l = Line((0, 0), (Rational(9, 2), 3))
  527. p1 = p.cut_section(l)[0]
  528. p2 = p.cut_section(l)[1]
  529. assert p1 == Polygon(
  530. Point2D(Rational(-9, 13), Rational(-6, 13)), Point2D(1, Rational(5, 2)), Point2D(Rational(24, 13), Rational(16, 13)),
  531. Point2D(Rational(12, 5), Rational(8, 5)), Point2D(3, Rational(5, 2)), Point2D(Rational(24, 7), Rational(16, 7)),
  532. Point2D(Rational(9, 2), 3), Point2D(-1, 3), Point2D(-1, Rational(-2, 3)))
  533. assert p2 == Polygon(Point2D(-1, -1), Point2D(Rational(-9, 13), Rational(-6, 13)), Point2D(Rational(24, 13), Rational(16, 13)),
  534. Point2D(2, 1), Point2D(Rational(12, 5), Rational(8, 5)), Point2D(Rational(24, 7), Rational(16, 7)), Point2D(4, 2), Point2D(5, 3),
  535. Point2D(Rational(9, 2), 3), Point2D(-1, Rational(-2, 3)))
  536. # convex polygon
  537. p = RegularPolygon(Point2D(0, 0), 6, 6)
  538. s = p.cut_section(Line((0, 0), slope=1))
  539. assert s[0] == Polygon(Point2D(-3*sqrt(3) + 9, -3*sqrt(3) + 9), Point2D(3, 3*sqrt(3)),
  540. Point2D(-3, 3*sqrt(3)), Point2D(-6, 0), Point2D(-9 + 3*sqrt(3), -9 + 3*sqrt(3)))
  541. assert s[1] == Polygon(Point2D(6, 0), Point2D(-3*sqrt(3) + 9, -3*sqrt(3) + 9),
  542. Point2D(-9 + 3*sqrt(3), -9 + 3*sqrt(3)), Point2D(-3, -3*sqrt(3)), Point2D(3, -3*sqrt(3)))
  543. # case where line does not intersects but coincides with the edge of polygon
  544. a, b = 20, 10
  545. t1, t2, t3, t4 = [(0, b), (0, 0), (a, 0), (a, b)]
  546. p = Polygon(t1, t2, t3, t4)
  547. p1, p2 = p.cut_section(Line((0, b), slope=0))
  548. assert p1 == None
  549. assert p2 == Polygon(Point2D(0, 10), Point2D(0, 0), Point2D(20, 0), Point2D(20, 10))
  550. p3, p4 = p.cut_section(Line((0, 0), slope=0))
  551. assert p3 == Polygon(Point2D(0, 10), Point2D(0, 0), Point2D(20, 0), Point2D(20, 10))
  552. assert p4 == None
  553. # case where the line does not intersect with a polygon at all
  554. raises(ValueError, lambda: p.cut_section(Line((0, a), slope=0)))
  555. def test_type_of_triangle():
  556. # Isoceles triangle
  557. p1 = Polygon(Point(0, 0), Point(5, 0), Point(2, 4))
  558. assert p1.is_isosceles() == True
  559. assert p1.is_scalene() == False
  560. assert p1.is_equilateral() == False
  561. # Scalene triangle
  562. p2 = Polygon (Point(0, 0), Point(0, 2), Point(4, 0))
  563. assert p2.is_isosceles() == False
  564. assert p2.is_scalene() == True
  565. assert p2.is_equilateral() == False
  566. # Equilateral triangle
  567. p3 = Polygon(Point(0, 0), Point(6, 0), Point(3, sqrt(27)))
  568. assert p3.is_isosceles() == True
  569. assert p3.is_scalene() == False
  570. assert p3.is_equilateral() == True
  571. def test_do_poly_distance():
  572. # Non-intersecting polygons
  573. square1 = Polygon (Point(0, 0), Point(0, 1), Point(1, 1), Point(1, 0))
  574. triangle1 = Polygon(Point(1, 2), Point(2, 2), Point(2, 1))
  575. assert square1._do_poly_distance(triangle1) == sqrt(2)/2
  576. # Polygons which sides intersect
  577. square2 = Polygon(Point(1, 0), Point(2, 0), Point(2, 1), Point(1, 1))
  578. with warns(UserWarning, \
  579. match="Polygons may intersect producing erroneous output", test_stacklevel=False):
  580. assert square1._do_poly_distance(square2) == 0
  581. # Polygons which bodies intersect
  582. triangle2 = Polygon(Point(0, -1), Point(2, -1), Point(S.Half, S.Half))
  583. with warns(UserWarning, \
  584. match="Polygons may intersect producing erroneous output", test_stacklevel=False):
  585. assert triangle2._do_poly_distance(square1) == 0