simetri.geometry.ellipse
defs for working with ellipses.
1"""defs for working with ellipses.""" 2 3import cmath 4from copy import deepcopy 5from math import cos, sin, pi, atan2, sqrt, degrees, ceil 6 7import numpy as np 8 9from ..graphics.shape import Shape, custom_attributes 10from ..graphics.batch import Batch 11from ..graphics.points import Points 12from ..graphics.affine import rotation_matrix 13from ..graphics.common import Point 14from ..graphics.all_enums import Types 15from ..geometry.geometry import ( 16 line_angle, 17 distance, 18 positive_angle, 19 rotate_point, 20 homogenize, 21) 22from ..canvas.style_map import shape_style_map 23from ..settings.settings import defaults 24from ..helpers.utilities import solve_quadratic_eq 25 26isclose = np.isclose 27 28 29class Arc(Shape): 30 """A circular or elliptic arc defined by a center, radius_x, radius_y, start angle, 31 and span angle. If radius_y is not provided, the arc is a circular arc.""" 32 33 def __init__( 34 self, 35 center: Point, 36 radius_x: float, 37 radius_y: float = None, 38 start_angle: float = 0, 39 span_angle: float = pi / 2, 40 rot_angle: float = 0, 41 n_points: int = None, 42 xform_matrix: "ndarray" = None, 43 **kwargs, 44 ): 45 """ 46 Args: 47 center (Point): The center of the arc. 48 radius_x (float): The x radius of the arc. 49 radius_y (float): The y radius for elliptical arcs. 50 start_angle (float): The starting angle of the arc. 51 span_angle (float): The span angle of the arc. 52 rot_angle (float, optional): Rotation angle. Defaults to 0. If negative, the arc is drawn clockwise. 53 xform_matrix (ndarray, optional): Transformation matrix. Defaults to None. 54 **kwargs: Additional keyword arguments. 55 """ 56 if radius_y is None: 57 radius_y = radius_x 58 if n_points is None: 59 n = defaults["n_arc_points"] 60 n_points = ceil(n * abs(span_angle) / (2 * pi)) 61 62 vertices = elliptic_arc_points( 63 center, radius_x, radius_y, start_angle, span_angle, n_points=n_points 64 ) 65 if rot_angle: 66 rot_matrix = rotation_matrix(rot_angle, center) 67 if xform_matrix is not None: 68 xform_matrix = np.dot(rot_matrix, xform_matrix) 69 else: 70 xform_matrix = rot_matrix 71 72 super().__init__(vertices, xform_matrix=xform_matrix, **kwargs) 73 self.subtype = Types.ARC 74 self.n_points = n_points 75 self.__dict__["start_angle"] = start_angle 76 self.__dict__["span_angle"] = span_angle 77 cx, cy = center[:2] 78 self._c = [cx, cy, 1] 79 _a = [radius_x, 0, 1] 80 _b = [0, radius_y, 1] 81 self._orig_triangle = [self._c[:], _a, _b] 82 83 def __setattr__(self, name, value): 84 """Set an attribute of the arc. 85 86 Args: 87 name (str): The name of the attribute. 88 value (Any): The value of the attribute. 89 """ 90 if name == "center": 91 diff = np.array(value[:2]) - np.array(self.center[:2]) 92 self.translate(diff[0], diff[1], reps=0) 93 elif name == "radius_x": 94 c, a, _ = self._orig_triangle @ self.xform_matrix 95 cur_radius = distance(c, a) 96 ratio = value / cur_radius 97 self.scale(ratio, 1, about=self.center) 98 elif name == "radius_y": 99 c, _, b = self._orig_triangle @ self.xform_matrix 100 cur_radius = distance(c, b) 101 ratio = value / cur_radius 102 self.scale(1, ratio, about=self.center) 103 elif name == "start_angle": 104 center, a, b = self._orig_triangle @ self.xform_matrix 105 a = distance(center, a) 106 b = distance(center, b) 107 span = self.span_angle 108 n_points = self.n_points 109 points = elliptic_arc_points(center, a, b, value, span, n_points) 110 self.primary_points = Points(points) 111 self.__dict__["start_angle"] = value 112 elif name == "span_angle": 113 center, a, b = self._orig_triangle @ self.xform_matrix 114 a = distance(center, a) 115 b = distance(center, b) 116 start = self.start_angle 117 n_points = self.n_points 118 points = elliptic_arc_points(center, a, b, start, value, n_points) 119 self.primary_points = Points(points) 120 self.__dict__["span_angle"] = value 121 else: 122 super().__setattr__(name, value) 123 124 @property 125 def center(self): 126 """Return the center of the arc. 127 128 Returns: 129 Point: The center of the arc. 130 """ 131 return (self._c @ self.xform_matrix).tolist()[:2] 132 133 @property 134 def radius_x(self): 135 """Return the x radius of the arc. 136 137 Returns: 138 float: The x radius of the arc. 139 """ 140 c, a, _ = self._orig_triangle @ self.xform_matrix 141 return distance(a, c) 142 143 @property 144 def radius_y(self): 145 """Return the y radius of the arc. 146 147 Returns: 148 float: The y radius of the arc. 149 """ 150 c, _, b = self._orig_triangle @ self.xform_matrix 151 return distance(b, c) 152 153 def copy(self): 154 """Return a copy of the arc.""" 155 center = self.center 156 start_angle = self.start_angle 157 span_angle = self.span_angle 158 radius_x = self.radius_x 159 radius_y = self.radius_y 160 161 arc = Arc(center, radius_x, radius_y, start_angle, span_angle, rot_angle=0) 162 arc.primary_points = self.primary_points.copy() 163 arc.xform_matrix = self.xform_matrix.copy() 164 arc._orig_triangle = deepcopy(self._orig_triangle) 165 arc._c = self._c[:] 166 arc.n_points = self.n_points 167 168 for attrib in shape_style_map: 169 setattr(arc, attrib, getattr(self, attrib)) 170 arc.subtype = self.subtype 171 custom_attribs = custom_attributes(self) 172 arc_attribs = ["center", "start_angle", "span_angle", "radius_x", "radius_y"] 173 for attrib in custom_attribs: 174 if attrib not in arc_attribs: 175 setattr(arc, attrib, getattr(self, attrib)) 176 177 return arc 178 179 180class Ellipse(Shape): 181 """An ellipse defined by center, width, and height.""" 182 183 def __init__( 184 self, 185 center: Point, 186 width: float, 187 height: float, 188 angle: float = 0, 189 xform_matrix: "ndarray" = None, 190 **kwargs, 191 ) -> None: 192 """ 193 Args: 194 center (Point): The center of the ellipse. 195 width (float): The width of the ellipse. 196 height (float): The height of the ellipse. 197 angle (float, optional): Rotation angle. Defaults to 0. 198 xform_matrix (ndarray, optional): Transformation matrix. Defaults to None. 199 **kwargs: Additional keyword arguments. 200 """ 201 n_points = defaults["n_ellipse_points"] 202 vertices = [ 203 tuple(p) 204 for p in ellipse_points(center, width / 2, height / 2, angle, n_points) 205 ] 206 super().__init__(vertices, closed=True, xform_matrix=xform_matrix, **kwargs) 207 a = width / 2 208 b = height / 2 209 self.a = a 210 self.b = b 211 self.center = center 212 self.width = width 213 self.height = height 214 self.angle = angle 215 self.smooth = True 216 self.closed = True 217 self.subtype = Types.ELLIPSE 218 219 @property 220 def closed(self): 221 """Return True ellipse is always closed. 222 223 Returns: 224 bool: Always returns True. 225 """ 226 return True 227 228 @closed.setter 229 def closed(self, value: bool): 230 pass 231 232 def _update(self, xform_matrix: np.array, reps: int = 0) -> Batch: 233 """Used internally. Update the shape with a transformation matrix. 234 235 Args: 236 xform_matrix (array): The transformation matrix. 237 reps (int, optional): The number of repetitions, defaults to 0. 238 239 Returns: 240 Batch: The updated shape or a batch of shapes. 241 """ 242 if reps == 0: 243 center = list(self.center[:2]) + [1] 244 start = list(self.vertices[0][:2]) + [1] 245 end = list(self.vertices[-1][:2]) + [1] 246 points = [center, start, end] 247 center2, start2, end2 = np.dot(points, xform_matrix).tolist() 248 self.center = center2[:2] 249 self.start_point = start2[:2] 250 self.end_point = end2[:2] 251 self.start_angle = line_angle(center2, start2) 252 253 return super()._update(xform_matrix, reps) 254 255 def copy(self): 256 """Return a copy of the ellipse. 257 258 Returns: 259 Ellipse: A copy of the ellipse. 260 """ 261 center = self.center 262 width = self.width 263 height = self.height 264 ellipse = Ellipse(center, width, height) 265 custom_attribs = custom_attributes(self) 266 for attrib in custom_attribs: 267 setattr(ellipse, attrib, getattr(self, attrib)) 268 269 return ellipse 270 271 272def ellipse_tangent(a, b, x, y, tol=0.001): 273 """Calculates the angle of the tangent line to an ellipse at the point (x, y). 274 275 Args: 276 a (float): Semi-major axis of the ellipse. 277 b (float): Semi-minor axis of the ellipse. 278 x (float): x-coordinate of the point. 279 y (float): y-coordinate of the point. 280 tol (float, optional): Tolerance for point on ellipse check. Defaults to .001. 281 282 Returns: 283 float: Angle of the tangent line in radians. 284 """ 285 if abs((x**2 / a**2) + (y**2 / b**2) - 1) >= tol: 286 res = False 287 else: 288 # res = atan2(-(b**2 * x), (a**2 * y)) 289 res = atan2((b**2 * x), -(a**2 * y)) 290 291 return res 292 293 294def r_central(a, b, theta): 295 """Return the radius (distance between the center and the intersection point) 296 of the ellipse at the given angle. 297 298 Args: 299 a (float): Semi-major axis of the ellipse. 300 b (float): Semi-minor axis of the ellipse. 301 theta (float): Angle in radians. 302 303 Returns: 304 float: Radius at the given angle. 305 """ 306 return (a * b) / sqrt((b * cos(theta)) ** 2 + (a * sin(theta)) ** 2) 307 308 309def ellipse_line_intersection(a, b, point): 310 """Return the intersection points of an ellipse and a line segment 311 connecting the given point to the ellipse center at (0, 0). 312 313 Args: 314 a (float): Semi-major axis of the ellipse. 315 b (float): Semi-minor axis of the ellipse. 316 point (tuple): Point coordinates (x, y). 317 318 Returns: 319 list: Intersection points. 320 """ 321 # adapted from http:# mathworld.wolfram.com/Ellipse-LineIntersection.html 322 # a, b is the ellipse width/2 and height/2 and (x_0, y_0) is the point 323 324 x_0, y_0 = point[:2] 325 x = ((a * b) / (sqrt(a**2 * y_0**2 + b**2 * x_0**2))) * x_0 326 y = ((a * b) / (sqrt(a**2 * y_0**2 + b**2 * x_0**2))) * y_0 327 328 return [(x, y), (-x, -y)] 329 330 331def elliptic_arc_points( 332 center, radius_x, radius_y, start_angle, span_angle, n_points=None 333): 334 """Generate points on an elliptic arc. 335 These are generated from the parametric equations of the ellipse. 336 They are not evenly spaced. 337 338 Args: 339 center (tuple): (x, y) coordinates of the ellipse center. 340 radius_x (float): Length of the semi-major axis. 341 radius_y (float): Length of the semi-minor axis. 342 start_angle (float): Starting angle of the arc. 343 span_angle (float): Span angle of the arc. 344 n_points (int): Number of points to generate. 345 346 Returns: 347 numpy.ndarray: Array of (x, y) coordinates of the ellipse points. 348 """ 349 rx = radius_x 350 if radius_y is None: 351 radius_y = radius_x 352 ry = radius_y 353 if n_points is None: 354 n = defaults["n_arc_points"] 355 n_points = ceil(n * abs(span_angle) / (2 * pi)) 356 start_angle = positive_angle(start_angle) 357 clockwise = span_angle < 0 358 if clockwise: 359 if start_angle + span_angle < 0: 360 end_angle = positive_angle(start_angle + span_angle) 361 t0 = get_ellipse_t_for_angle(end_angle, rx, ry) 362 t1 = get_ellipse_t_for_angle(2 * pi, rx, ry) 363 t = np.linspace(t0, t1, n_points) 364 x = center[0] + rx * np.cos(t) 365 y = center[1] + ry * np.sin(t) 366 slice1 = np.column_stack((x, y)) 367 t0 = get_ellipse_t_for_angle(0, rx, ry) 368 t1 = get_ellipse_t_for_angle(start_angle, rx, ry) 369 t = np.linspace(t0, t1, n_points) 370 x = center[0] + rx * np.cos(t) 371 y = center[1] + ry * np.sin(t) 372 slice2 = np.column_stack((x, y)) 373 res = np.flip(np.concatenate((slice1, slice2)), axis=0) 374 else: 375 end_angle = start_angle + span_angle 376 t0 = get_ellipse_t_for_angle(end_angle, rx, ry) 377 t1 = get_ellipse_t_for_angle(start_angle, rx, ry) 378 t = np.linspace(t0, t1, n_points) 379 x = center[0] + rx * np.cos(t) 380 y = center[1] + ry * np.sin(t) 381 res = np.flip(np.column_stack((x, y)), axis=0) 382 383 else: 384 if start_angle + span_angle > 2 * pi: 385 t0 = get_ellipse_t_for_angle(start_angle, rx, ry) 386 t1 = get_ellipse_t_for_angle(2 * pi, rx, ry) 387 t = np.linspace(t0, t1, n_points) 388 x = center[0] + rx * np.cos(t) 389 y = center[1] + ry * np.sin(t) 390 slice1 = np.column_stack((x, y)) 391 392 t0 = get_ellipse_t_for_angle(0, rx, ry) 393 t1 = get_ellipse_t_for_angle(start_angle + span_angle - 2 * pi, rx, ry) 394 t = np.linspace(t0, t1, n_points) 395 x = center[0] + rx * np.cos(t) 396 y = center[1] + ry * np.sin(t) 397 slice2 = np.column_stack((x, y)) 398 399 res = np.concatenate((slice1, slice2)) 400 else: 401 end_angle = start_angle + span_angle 402 t0 = get_ellipse_t_for_angle(start_angle, rx, ry) 403 t1 = get_ellipse_t_for_angle(end_angle, rx, ry) 404 t = np.linspace(t0, t1, n_points) 405 x = center[0] + rx * np.cos(t) 406 y = center[1] + ry * np.sin(t) 407 res = np.column_stack((x, y)) 408 409 return res 410 411 412def ellipse_points( 413 center: Point, 414 a: float, 415 b: float, 416 angle: float, 417 n_points: int = None, 418) -> np.ndarray: 419 """Generate points on an ellipse. 420 These are generated from the parametric equations of the ellipse. 421 They are not evenly spaced. 422 423 Args: 424 center (tuple): (x, y) coordinates of the ellipse center. 425 a (float): Length of the semi-major axis. 426 b (float): Length of the semi-minor axis. 427 angle (float): Rotation angle of the ellipse. 428 n_points (int): Number of points to generate. 429 430 Returns: 431 numpy.ndarray: Array of (x, y) coordinates of the ellipse points. 432 """ 433 if n_points is None: 434 n_points = defaults["n_ellipse_points"] 435 436 t = np.linspace(0, 2 * pi, n_points) 437 x = center[0] + a * np.cos(t) 438 y = center[1] + b * np.sin(t) 439 440 points = homogenize(np.column_stack((x, y))) @ rotation_matrix(angle, center) 441 442 return points[:, :2].tolist() 443 444def elliptic_arclength(t_0, t_1, a, b): 445 """Return the arclength of an ellipse between the given parametric angles. 446 The ellipse has semi-major axis a and semi-minor axis b. 447 448 Args: 449 t_0 (float): Starting parametric angle. 450 t_1 (float): Ending parametric angle. 451 a (float): Semi-major axis of the ellipse. 452 b (float): Semi-minor axis of the ellipse. 453 454 Returns: 455 float: Arclength of the ellipse. 456 """ 457 from scipy.special import ellipeinc # this takes too long to import 458 459 m = 1 - (b / a) ** 2 460 t1 = ellipeinc(t_1 - 0.5 * pi, m) 461 t0 = ellipeinc(t_0 - 0.5 * pi, m) 462 return a * (t1 - t0) 463 464 465def central_to_parametric_angle(a, b, phi): 466 """ 467 Converts a central angle to a parametric angle on an ellipse. 468 469 Args: 470 a (float): Semi-major axis of the ellipse. 471 b (float): Semi-minor axis of the ellipse. 472 phi (float): Angle of the line intersecting the center and the point. 473 474 Returns: 475 float: Parametric angle (in radians). 476 """ 477 t = atan2((a / b) * sin(phi), cos(phi)) 478 if t < 0: 479 t += 2 * pi 480 481 return t 482 483 484def parametric_to_central_angle(a, b, t): 485 """ 486 Converts a parametric angle on an ellipse to a central angle. 487 488 Args: 489 a (float): Semi-major axis of the ellipse. 490 b (float): Semi-minor axis of the ellipse. 491 t (float): Parametric angle (in radians). 492 493 Returns: 494 float: Angle of the line intersecting the center and the point. 495 """ 496 phi = atan2((b / a) * sin(t), cos(t)) 497 if phi < 0: 498 phi += 2 * pi 499 500 return phi 501 502 503def ellipse_point(a, b, angle): 504 """Return a point on an ellipse with the given a=width/2, b=height/2, and angle. 505 angle is the central-angle and in radians. 506 507 Args: 508 a (float): Semi-major axis of the ellipse. 509 b (float): Semi-minor axis of the ellipse. 510 angle (float): Central angle in radians. 511 512 Returns: 513 tuple: Coordinates of the point on the ellipse. 514 """ 515 r = r_central(a, b, angle) 516 517 return (r * cos(angle), r * sin(angle)) 518 519 520def ellipse_param_point(a, b, t): 521 """Return a point on an ellipse with the given a=width/2, b=height/2, and parametric angle. 522 t is the parametric angle and in radians. 523 524 Args: 525 a (float): Semi-major axis of the ellipse. 526 b (float): Semi-minor axis of the ellipse. 527 t (float): Parametric angle in radians. 528 529 Returns: 530 tuple: Coordinates of the point on the ellipse. 531 """ 532 return (a * cos(t), b * sin(t)) 533 534 535def get_ellipse_t_for_angle(angle, a, b): 536 """ 537 Calculates the parameter t for a given angle on an ellipse. 538 539 Args: 540 angle (float): The angle in radians. 541 a (float): Semi-major axis of the ellipse. 542 b (float): Semi-minor axis of the ellipse. 543 544 Returns: 545 float: The parameter t. 546 """ 547 t = atan2(a * sin(angle), b * cos(angle)) 548 if t < 0: 549 t += 2 * pi 550 return t 551 552 553def ellipse_central_angle(t, a, b): 554 """ 555 Calculates the central angle of an ellipse for a given parameter t. 556 557 Args: 558 t (float): The parameter value. 559 a (float): The semi-major axis of the ellipse. 560 b (float): The semi-minor axis of the ellipse. 561 562 Returns: 563 float: The central angle in radians. 564 """ 565 theta = atan2(a * sin(t), b * cos(t)) 566 567 return theta 568 569 570def ellipse_intersection(x1, y1, a, b, phi, x2, y2, c, d, phi2): 571 """Calculate the intersection points of two ellipses. 572 The ellipses are defined by their center, radii, and rotation angle.""" 573 # Taken from https:# github.com/VoyakaGOD/intersection-of-two-ellipses/blob/master/geometry.js 574 phi0 = phi2 - phi 575 576 p0 = rotate_point(((x2 - x1), (y2 - y1)), -phi) 577 578 # cos(t1) = A + Bcos(t2) + Csin(t2) 579 # sin(t1) = D + Ecos(t2) + Fsin(t2) 580 A = p0[0] / a 581 B = c * cos(phi0) / a 582 C = -d * sin(phi0) / a 583 D = p0[1] / b 584 E = c * sin(phi0) / b 585 F = d * cos(phi0) / b 586 587 # Gx^2 + Hx + Ixy + Jy + K = 0, x = cos(t2), y = sin(t2) 588 G = B * B + E * E - C * C - F * F 589 H = 2 * (A * B + D * E) 590 I = 2 * (B * C + E * F) 591 J = 2 * (A * C + D * F) 592 K = A * A + D * D + C * C + F * F - 1 593 594 roots = [] 595 L = G * G + I * I 596 if isclose(L, 0, 0, 1e-7): 597 # Gx + Hy + K = 0 598 roots = solve_quadratic_eq(G, H, K) 599 600 elif isclose(I, 0, 1e-7): 601 # Gx^2 + Hx + K = 0 602 roots = solve_quadratic_eq(G, H, K) 603 604 elif isclose(G, 0, 1e-7): 605 # Hx + Jy + K = 0 606 roots = solve_quadratic_eq(H * H + J * J, 2 * K * H, K * K - J * J) 607 608 else: 609 # Lx^4 + Mx^3 + Nx^2 + Ox + P = 0 610 M = 2 * (G * H + I * J) 611 N = H * H + 2 * G * K + J * J - I * I 612 O = 2 * (K * H - I * J) 613 P = K * K - J * J 614 615 iL = 1 / L 616 roots = solve_quartic_equation(M * iL, N * iL, O * iL, P * iL) 617 618 points = [] 619 # for(i = 0 i < roots.length i++) 620 # for i in range(len(roots)): 621 for i, x in enumerate(roots): 622 # x = roots[i] 623 if isclose(I * x + J, 0, 1e-7): 624 y = sqrt(1 - x * x) 625 points.append((x, y)) 626 627 if not isclose(y, 0, 1e-7): 628 points.append((x, -y)) 629 else: 630 y = -(G * x * x + H * x + K) / (I * x + J) 631 if abs(y) < 1e-2: 632 y = sqrt(1 - x * x) 633 points.append((x, y)) 634 points.append((x, -y)) 635 else: 636 points.append((x, y)) 637 638 for i, pnt in enumerate(points): 639 x, y = pnt 640 points[i] = (A + B * x + C * y) * a, (D + E * x + F * y) * b 641 642 for i, pnt in enumerate(points): 643 x, y = pnt 644 x, y = rotate_point((x, y), phi) 645 points[i] = (x + x1, y + y1) 646 647 return points 648 649 650# def IsCloseToZero(num): 651# return abs(num) < 1e-7 652 653 654def inverse_complex_number(z: complex) -> complex: 655 """ 656 Calculate the inverse of a complex number. 657 658 Args: 659 z (complex): The complex number. 660 661 Returns: 662 complex: The inverse of the complex number. 663 """ 664 a = z.real 665 b = z.imag 666 667 if a == 0 and b == 0: 668 raise ZeroDivisionError("Cannot calculate the inverse of 0") 669 670 return complex(a / (a**2 + b**2), -b / (a**2 + b**2)) 671 672 673def Re(num): 674 return complex(num, 0) 675 # return Complex(num, 0) 676 677 678def Im(num): 679 return complex(0, num) 680 # return Complex(0, num) 681 682 683def Sqrt(complex_): 684 cmath.sqrt(complex_) 685 # return complex.sqrt 686 687 688def Qbrt(complex_): 689 # Taken from https:# github.com/VoyakaGOD/intersection-of-two-ellipses/blob/master/quartic.js 690 691 angle = cmath.phase(complex_) * 0.33333333333 692 # angle = atan2(complex.y, complex.x) * 0.33333333333 693 l = pow( 694 complex_.real * complex_.real + complex_.imag * complex_.imag, 0.16666666666 695 ) 696 # return Complex(l * cos(angle), l * sin(angle)) 697 return complex(l * cos(angle), l * sin(angle)) 698 699 700# x^2 + bx + c = 0 701def solve_complex_quadratic_equation(b, c): 702 # Taken from https:# github.com/VoyakaGOD/intersection-of-two-ellipses/blob/master/quartic.js 703 704 # sqrtD = Sqrt(b.MulComplex(b).Sub(c.Mul(4))) 705 sqrtD = cmath.sqrt(b * b - (c * 4)) 706 return [(b + sqrtD) * (-0.5), (b - sqrtD) * (-0.5)] 707 708 709# x^3 + ax^2 + bx + c = 0 710def get_one_cubic_equation_root(a, b, c): 711 # Taken from https:# github.com/VoyakaGOD/intersection-of-two-ellipses/blob/master/quartic.js 712 713 p = b - a * a * 0.33333333333 714 q = (2 / 27) * a * a * a - a * b * 0.33333333333 + c 715 sqrt_Q = cmath.sqrt(Re(0.03703703703 * p * p * p + 0.25 * q * q)) 716 A = Qbrt(Re(-0.5 * q) - sqrt_Q) 717 if A.real == 0 and A.imag == 0: 718 A = Qbrt(Re(-0.5 * q) + sqrt_Q) 719 B = inverse_complex_number(A) * (-p * 0.33333333333) 720 721 return (A + B) - (Re(a * 0.33333333333)) 722 723# x^4 + ax^3 + bx^2 + cx + d = 0 724def solve_quartic_equation(a, b, c, d): 725 # Taken from https:# github.com/VoyakaGOD/intersection-of-two-ellipses/blob/master/geometry.js 726 727 a2 = a * a 728 a3 = a2 * a 729 a4 = a3 * a 730 p = b - (3/8) * a2 731 q = (1 / 8) * a3 - 0.5 * a * b + c 732 r = d - 0.25 * a * c + (1 / 16) * a2 * b - (3 / 256) * a4 733 734 result = [] 735 736 if(isclose(q, 0, 1e-7)): 737 D = p * p - 4 * r 738 if(abs(D) < 1e-5): 739 m = -0.5 * p 740 if (m >= 0): 741 result.append(sqrt(m)) 742 if (m > 0): 743 result.append(-sqrt(m)) 744 elif (D > 0): 745 sqrt_D = sqrt(D) 746 m1 = (-p - sqrt_D) * 0.5 747 m2 = (-p + sqrt_D) * 0.5 748 sqrt_m1 = sqrt(m1) 749 sqrt_m2 = sqrt(m2) 750 if (m1 >= 0): 751 result.append(sqrt_m1) 752 if (m1 > 0): 753 result.append(-sqrt_m1) 754 if (m2 >= 0): 755 result.append(sqrt_m2) 756 if (m2 > 0): 757 result.append(-sqrt_m2) 758 else: 759 t = get_one_cubic_equation_root(2 * p, p * p-4 * r, -q * q) 760 z = cmath.sqrt(t) 761 u = (Re(p)+(t)) * 0.5 762 v = inverse_complex_number(z) * (q * 0.5) 763 x12 = solve_complex_quadratic_equation(z, u - v) 764 x34 = solve_complex_quadratic_equation(-z, u + v) 765 766 if (isclose(x12[0].imag, 0, 1e-7)): 767 result.append(x12[0].real) 768 if (isclose(x12[1].imag, 0, 1e-7)): 769 result.append(x12[1].real) 770 if (isclose(x34[0].imag, 0, 1e-7)): 771 result.append(x34[0].real) 772 if (isclose(x34[1].imag, 0, 1e-7)): 773 result.append(x34[1].real) 774 775 # for(i = 0 i < result.length i++) 776 for i in range(len(result)): 777 result[i] -= 0.25 * a 778 779 return result
@array_function_dispatch(_isclose_dispatcher)
def
isclose(a, b, rtol=1e-05, atol=1e-08, equal_nan=False):
2316@array_function_dispatch(_isclose_dispatcher) 2317def isclose(a, b, rtol=1.e-5, atol=1.e-8, equal_nan=False): 2318 """ 2319 Returns a boolean array where two arrays are element-wise equal within a 2320 tolerance. 2321 2322 The tolerance values are positive, typically very small numbers. The 2323 relative difference (`rtol` * abs(`b`)) and the absolute difference 2324 `atol` are added together to compare against the absolute difference 2325 between `a` and `b`. 2326 2327 .. warning:: The default `atol` is not appropriate for comparing numbers 2328 with magnitudes much smaller than one (see Notes). 2329 2330 Parameters 2331 ---------- 2332 a, b : array_like 2333 Input arrays to compare. 2334 rtol : array_like 2335 The relative tolerance parameter (see Notes). 2336 atol : array_like 2337 The absolute tolerance parameter (see Notes). 2338 equal_nan : bool 2339 Whether to compare NaN's as equal. If True, NaN's in `a` will be 2340 considered equal to NaN's in `b` in the output array. 2341 2342 Returns 2343 ------- 2344 y : array_like 2345 Returns a boolean array of where `a` and `b` are equal within the 2346 given tolerance. If both `a` and `b` are scalars, returns a single 2347 boolean value. 2348 2349 See Also 2350 -------- 2351 allclose 2352 math.isclose 2353 2354 Notes 2355 ----- 2356 .. versionadded:: 1.7.0 2357 2358 For finite values, isclose uses the following equation to test whether 2359 two floating point values are equivalent.:: 2360 2361 absolute(a - b) <= (atol + rtol * absolute(b)) 2362 2363 Unlike the built-in `math.isclose`, the above equation is not symmetric 2364 in `a` and `b` -- it assumes `b` is the reference value -- so that 2365 `isclose(a, b)` might be different from `isclose(b, a)`. 2366 2367 The default value of `atol` is not appropriate when the reference value 2368 `b` has magnitude smaller than one. For example, it is unlikely that 2369 ``a = 1e-9`` and ``b = 2e-9`` should be considered "close", yet 2370 ``isclose(1e-9, 2e-9)`` is ``True`` with default settings. Be sure 2371 to select `atol` for the use case at hand, especially for defining the 2372 threshold below which a non-zero value in `a` will be considered "close" 2373 to a very small or zero value in `b`. 2374 2375 `isclose` is not defined for non-numeric data types. 2376 :class:`bool` is considered a numeric data-type for this purpose. 2377 2378 Examples 2379 -------- 2380 >>> np.isclose([1e10,1e-7], [1.00001e10,1e-8]) 2381 array([ True, False]) 2382 >>> np.isclose([1e10,1e-8], [1.00001e10,1e-9]) 2383 array([ True, True]) 2384 >>> np.isclose([1e10,1e-8], [1.0001e10,1e-9]) 2385 array([False, True]) 2386 >>> np.isclose([1.0, np.nan], [1.0, np.nan]) 2387 array([ True, False]) 2388 >>> np.isclose([1.0, np.nan], [1.0, np.nan], equal_nan=True) 2389 array([ True, True]) 2390 >>> np.isclose([1e-8, 1e-7], [0.0, 0.0]) 2391 array([ True, False]) 2392 >>> np.isclose([1e-100, 1e-7], [0.0, 0.0], atol=0.0) 2393 array([False, False]) 2394 >>> np.isclose([1e-10, 1e-10], [1e-20, 0.0]) 2395 array([ True, True]) 2396 >>> np.isclose([1e-10, 1e-10], [1e-20, 0.999999e-10], atol=0.0) 2397 array([False, True]) 2398 """ 2399 # Turn all but python scalars into arrays. 2400 x, y, atol, rtol = ( 2401 a if isinstance(a, (int, float, complex)) else asanyarray(a) 2402 for a in (a, b, atol, rtol)) 2403 2404 # Make sure y is an inexact type to avoid bad behavior on abs(MIN_INT). 2405 # This will cause casting of x later. Also, make sure to allow subclasses 2406 # (e.g., for numpy.ma). 2407 # NOTE: We explicitly allow timedelta, which used to work. This could 2408 # possibly be deprecated. See also gh-18286. 2409 # timedelta works if `atol` is an integer or also a timedelta. 2410 # Although, the default tolerances are unlikely to be useful 2411 if (dtype := getattr(y, "dtype", None)) is not None and dtype.kind != "m": 2412 dt = multiarray.result_type(y, 1.) 2413 y = asanyarray(y, dtype=dt) 2414 elif isinstance(y, int): 2415 y = float(y) 2416 2417 with errstate(invalid='ignore'), _no_nep50_warning(): 2418 result = (less_equal(abs(x-y), atol + rtol * abs(y)) 2419 & isfinite(y) 2420 | (x == y)) 2421 if equal_nan: 2422 result |= isnan(x) & isnan(y) 2423 2424 return result[()] # Flatten 0d arrays to scalars
Returns a boolean array where two arrays are element-wise equal within a tolerance.
The tolerance values are positive, typically very small numbers. The
relative difference (rtol * abs(b)) and the absolute difference
atol are added together to compare against the absolute difference
between a and b.
The default atol is not appropriate for comparing numbers
with magnitudes much smaller than one (see Notes).
Parameters
a, b : array_like
Input arrays to compare.
rtol : array_like
The relative tolerance parameter (see Notes).
atol : array_like
The absolute tolerance parameter (see Notes).
equal_nan : bool
Whether to compare NaN's as equal. If True, NaN's in a will be
considered equal to NaN's in b in the output array.
Returns
y : array_like
Returns a boolean array of where a and b are equal within the
given tolerance. If both a and b are scalars, returns a single
boolean value.
See Also
allclose math.isclose
Notes
New in version 1.7.0.
For finite values, isclose uses the following equation to test whether two floating point values are equivalent.::
absolute(a - b) <= (atol + rtol * absolute(b))
Unlike the built-in math.isclose, the above equation is not symmetric
in a and b -- it assumes b is the reference value -- so that
isclose(a, b) might be different from isclose(b, a).
The default value of atol is not appropriate when the reference value
b has magnitude smaller than one. For example, it is unlikely that
a = 1e-9 and b = 2e-9 should be considered "close", yet
isclose(1e-9, 2e-9) is True with default settings. Be sure
to select atol for the use case at hand, especially for defining the
threshold below which a non-zero value in a will be considered "close"
to a very small or zero value in b.
isclose is not defined for non-numeric data types.
bool is considered a numeric data-type for this purpose.
Examples
>>> np.isclose([1e10,1e-7], [1.00001e10,1e-8])
array([ True, False])
>>> np.isclose([1e10,1e-8], [1.00001e10,1e-9])
array([ True, True])
>>> np.isclose([1e10,1e-8], [1.0001e10,1e-9])
array([False, True])
>>> np.isclose([1.0, np.nan], [1.0, np.nan])
array([ True, False])
>>> np.isclose([1.0, np.nan], [1.0, np.nan], equal_nan=True)
array([ True, True])
>>> np.isclose([1e-8, 1e-7], [0.0, 0.0])
array([ True, False])
>>> np.isclose([1e-100, 1e-7], [0.0, 0.0], atol=0.0)
array([False, False])
>>> np.isclose([1e-10, 1e-10], [1e-20, 0.0])
array([ True, True])
>>> np.isclose([1e-10, 1e-10], [1e-20, 0.999999e-10], atol=0.0)
array([False, True])
30class Arc(Shape): 31 """A circular or elliptic arc defined by a center, radius_x, radius_y, start angle, 32 and span angle. If radius_y is not provided, the arc is a circular arc.""" 33 34 def __init__( 35 self, 36 center: Point, 37 radius_x: float, 38 radius_y: float = None, 39 start_angle: float = 0, 40 span_angle: float = pi / 2, 41 rot_angle: float = 0, 42 n_points: int = None, 43 xform_matrix: "ndarray" = None, 44 **kwargs, 45 ): 46 """ 47 Args: 48 center (Point): The center of the arc. 49 radius_x (float): The x radius of the arc. 50 radius_y (float): The y radius for elliptical arcs. 51 start_angle (float): The starting angle of the arc. 52 span_angle (float): The span angle of the arc. 53 rot_angle (float, optional): Rotation angle. Defaults to 0. If negative, the arc is drawn clockwise. 54 xform_matrix (ndarray, optional): Transformation matrix. Defaults to None. 55 **kwargs: Additional keyword arguments. 56 """ 57 if radius_y is None: 58 radius_y = radius_x 59 if n_points is None: 60 n = defaults["n_arc_points"] 61 n_points = ceil(n * abs(span_angle) / (2 * pi)) 62 63 vertices = elliptic_arc_points( 64 center, radius_x, radius_y, start_angle, span_angle, n_points=n_points 65 ) 66 if rot_angle: 67 rot_matrix = rotation_matrix(rot_angle, center) 68 if xform_matrix is not None: 69 xform_matrix = np.dot(rot_matrix, xform_matrix) 70 else: 71 xform_matrix = rot_matrix 72 73 super().__init__(vertices, xform_matrix=xform_matrix, **kwargs) 74 self.subtype = Types.ARC 75 self.n_points = n_points 76 self.__dict__["start_angle"] = start_angle 77 self.__dict__["span_angle"] = span_angle 78 cx, cy = center[:2] 79 self._c = [cx, cy, 1] 80 _a = [radius_x, 0, 1] 81 _b = [0, radius_y, 1] 82 self._orig_triangle = [self._c[:], _a, _b] 83 84 def __setattr__(self, name, value): 85 """Set an attribute of the arc. 86 87 Args: 88 name (str): The name of the attribute. 89 value (Any): The value of the attribute. 90 """ 91 if name == "center": 92 diff = np.array(value[:2]) - np.array(self.center[:2]) 93 self.translate(diff[0], diff[1], reps=0) 94 elif name == "radius_x": 95 c, a, _ = self._orig_triangle @ self.xform_matrix 96 cur_radius = distance(c, a) 97 ratio = value / cur_radius 98 self.scale(ratio, 1, about=self.center) 99 elif name == "radius_y": 100 c, _, b = self._orig_triangle @ self.xform_matrix 101 cur_radius = distance(c, b) 102 ratio = value / cur_radius 103 self.scale(1, ratio, about=self.center) 104 elif name == "start_angle": 105 center, a, b = self._orig_triangle @ self.xform_matrix 106 a = distance(center, a) 107 b = distance(center, b) 108 span = self.span_angle 109 n_points = self.n_points 110 points = elliptic_arc_points(center, a, b, value, span, n_points) 111 self.primary_points = Points(points) 112 self.__dict__["start_angle"] = value 113 elif name == "span_angle": 114 center, a, b = self._orig_triangle @ self.xform_matrix 115 a = distance(center, a) 116 b = distance(center, b) 117 start = self.start_angle 118 n_points = self.n_points 119 points = elliptic_arc_points(center, a, b, start, value, n_points) 120 self.primary_points = Points(points) 121 self.__dict__["span_angle"] = value 122 else: 123 super().__setattr__(name, value) 124 125 @property 126 def center(self): 127 """Return the center of the arc. 128 129 Returns: 130 Point: The center of the arc. 131 """ 132 return (self._c @ self.xform_matrix).tolist()[:2] 133 134 @property 135 def radius_x(self): 136 """Return the x radius of the arc. 137 138 Returns: 139 float: The x radius of the arc. 140 """ 141 c, a, _ = self._orig_triangle @ self.xform_matrix 142 return distance(a, c) 143 144 @property 145 def radius_y(self): 146 """Return the y radius of the arc. 147 148 Returns: 149 float: The y radius of the arc. 150 """ 151 c, _, b = self._orig_triangle @ self.xform_matrix 152 return distance(b, c) 153 154 def copy(self): 155 """Return a copy of the arc.""" 156 center = self.center 157 start_angle = self.start_angle 158 span_angle = self.span_angle 159 radius_x = self.radius_x 160 radius_y = self.radius_y 161 162 arc = Arc(center, radius_x, radius_y, start_angle, span_angle, rot_angle=0) 163 arc.primary_points = self.primary_points.copy() 164 arc.xform_matrix = self.xform_matrix.copy() 165 arc._orig_triangle = deepcopy(self._orig_triangle) 166 arc._c = self._c[:] 167 arc.n_points = self.n_points 168 169 for attrib in shape_style_map: 170 setattr(arc, attrib, getattr(self, attrib)) 171 arc.subtype = self.subtype 172 custom_attribs = custom_attributes(self) 173 arc_attribs = ["center", "start_angle", "span_angle", "radius_x", "radius_y"] 174 for attrib in custom_attribs: 175 if attrib not in arc_attribs: 176 setattr(arc, attrib, getattr(self, attrib)) 177 178 return arc
A circular or elliptic arc defined by a center, radius_x, radius_y, start angle, and span angle. If radius_y is not provided, the arc is a circular arc.
Arc( center: Sequence[float], radius_x: float, radius_y: float = None, start_angle: float = 0, span_angle: float = 1.5707963267948966, rot_angle: float = 0, n_points: int = None, xform_matrix: 'ndarray' = None, **kwargs)
34 def __init__( 35 self, 36 center: Point, 37 radius_x: float, 38 radius_y: float = None, 39 start_angle: float = 0, 40 span_angle: float = pi / 2, 41 rot_angle: float = 0, 42 n_points: int = None, 43 xform_matrix: "ndarray" = None, 44 **kwargs, 45 ): 46 """ 47 Args: 48 center (Point): The center of the arc. 49 radius_x (float): The x radius of the arc. 50 radius_y (float): The y radius for elliptical arcs. 51 start_angle (float): The starting angle of the arc. 52 span_angle (float): The span angle of the arc. 53 rot_angle (float, optional): Rotation angle. Defaults to 0. If negative, the arc is drawn clockwise. 54 xform_matrix (ndarray, optional): Transformation matrix. Defaults to None. 55 **kwargs: Additional keyword arguments. 56 """ 57 if radius_y is None: 58 radius_y = radius_x 59 if n_points is None: 60 n = defaults["n_arc_points"] 61 n_points = ceil(n * abs(span_angle) / (2 * pi)) 62 63 vertices = elliptic_arc_points( 64 center, radius_x, radius_y, start_angle, span_angle, n_points=n_points 65 ) 66 if rot_angle: 67 rot_matrix = rotation_matrix(rot_angle, center) 68 if xform_matrix is not None: 69 xform_matrix = np.dot(rot_matrix, xform_matrix) 70 else: 71 xform_matrix = rot_matrix 72 73 super().__init__(vertices, xform_matrix=xform_matrix, **kwargs) 74 self.subtype = Types.ARC 75 self.n_points = n_points 76 self.__dict__["start_angle"] = start_angle 77 self.__dict__["span_angle"] = span_angle 78 cx, cy = center[:2] 79 self._c = [cx, cy, 1] 80 _a = [radius_x, 0, 1] 81 _b = [0, radius_y, 1] 82 self._orig_triangle = [self._c[:], _a, _b]
Arguments:
- center (Point): The center of the arc.
- radius_x (float): The x radius of the arc.
- radius_y (float): The y radius for elliptical arcs.
- start_angle (float): The starting angle of the arc.
- span_angle (float): The span angle of the arc.
- rot_angle (float, optional): Rotation angle. Defaults to 0. If negative, the arc is drawn clockwise.
- xform_matrix (ndarray, optional): Transformation matrix. Defaults to None.
- **kwargs: Additional keyword arguments.
center
125 @property 126 def center(self): 127 """Return the center of the arc. 128 129 Returns: 130 Point: The center of the arc. 131 """ 132 return (self._c @ self.xform_matrix).tolist()[:2]
Return the center of the arc.
Returns:
Point: The center of the arc.
radius_x
134 @property 135 def radius_x(self): 136 """Return the x radius of the arc. 137 138 Returns: 139 float: The x radius of the arc. 140 """ 141 c, a, _ = self._orig_triangle @ self.xform_matrix 142 return distance(a, c)
Return the x radius of the arc.
Returns:
float: The x radius of the arc.
radius_y
144 @property 145 def radius_y(self): 146 """Return the y radius of the arc. 147 148 Returns: 149 float: The y radius of the arc. 150 """ 151 c, _, b = self._orig_triangle @ self.xform_matrix 152 return distance(b, c)
Return the y radius of the arc.
Returns:
float: The y radius of the arc.
def
copy(self):
154 def copy(self): 155 """Return a copy of the arc.""" 156 center = self.center 157 start_angle = self.start_angle 158 span_angle = self.span_angle 159 radius_x = self.radius_x 160 radius_y = self.radius_y 161 162 arc = Arc(center, radius_x, radius_y, start_angle, span_angle, rot_angle=0) 163 arc.primary_points = self.primary_points.copy() 164 arc.xform_matrix = self.xform_matrix.copy() 165 arc._orig_triangle = deepcopy(self._orig_triangle) 166 arc._c = self._c[:] 167 arc.n_points = self.n_points 168 169 for attrib in shape_style_map: 170 setattr(arc, attrib, getattr(self, attrib)) 171 arc.subtype = self.subtype 172 custom_attribs = custom_attributes(self) 173 arc_attribs = ["center", "start_angle", "span_angle", "radius_x", "radius_y"] 174 for attrib in custom_attribs: 175 if attrib not in arc_attribs: 176 setattr(arc, attrib, getattr(self, attrib)) 177 178 return arc
Return a copy of the arc.
Inherited Members
181class Ellipse(Shape): 182 """An ellipse defined by center, width, and height.""" 183 184 def __init__( 185 self, 186 center: Point, 187 width: float, 188 height: float, 189 angle: float = 0, 190 xform_matrix: "ndarray" = None, 191 **kwargs, 192 ) -> None: 193 """ 194 Args: 195 center (Point): The center of the ellipse. 196 width (float): The width of the ellipse. 197 height (float): The height of the ellipse. 198 angle (float, optional): Rotation angle. Defaults to 0. 199 xform_matrix (ndarray, optional): Transformation matrix. Defaults to None. 200 **kwargs: Additional keyword arguments. 201 """ 202 n_points = defaults["n_ellipse_points"] 203 vertices = [ 204 tuple(p) 205 for p in ellipse_points(center, width / 2, height / 2, angle, n_points) 206 ] 207 super().__init__(vertices, closed=True, xform_matrix=xform_matrix, **kwargs) 208 a = width / 2 209 b = height / 2 210 self.a = a 211 self.b = b 212 self.center = center 213 self.width = width 214 self.height = height 215 self.angle = angle 216 self.smooth = True 217 self.closed = True 218 self.subtype = Types.ELLIPSE 219 220 @property 221 def closed(self): 222 """Return True ellipse is always closed. 223 224 Returns: 225 bool: Always returns True. 226 """ 227 return True 228 229 @closed.setter 230 def closed(self, value: bool): 231 pass 232 233 def _update(self, xform_matrix: np.array, reps: int = 0) -> Batch: 234 """Used internally. Update the shape with a transformation matrix. 235 236 Args: 237 xform_matrix (array): The transformation matrix. 238 reps (int, optional): The number of repetitions, defaults to 0. 239 240 Returns: 241 Batch: The updated shape or a batch of shapes. 242 """ 243 if reps == 0: 244 center = list(self.center[:2]) + [1] 245 start = list(self.vertices[0][:2]) + [1] 246 end = list(self.vertices[-1][:2]) + [1] 247 points = [center, start, end] 248 center2, start2, end2 = np.dot(points, xform_matrix).tolist() 249 self.center = center2[:2] 250 self.start_point = start2[:2] 251 self.end_point = end2[:2] 252 self.start_angle = line_angle(center2, start2) 253 254 return super()._update(xform_matrix, reps) 255 256 def copy(self): 257 """Return a copy of the ellipse. 258 259 Returns: 260 Ellipse: A copy of the ellipse. 261 """ 262 center = self.center 263 width = self.width 264 height = self.height 265 ellipse = Ellipse(center, width, height) 266 custom_attribs = custom_attributes(self) 267 for attrib in custom_attribs: 268 setattr(ellipse, attrib, getattr(self, attrib)) 269 270 return ellipse
An ellipse defined by center, width, and height.
Ellipse( center: Sequence[float], width: float, height: float, angle: float = 0, xform_matrix: 'ndarray' = None, **kwargs)
184 def __init__( 185 self, 186 center: Point, 187 width: float, 188 height: float, 189 angle: float = 0, 190 xform_matrix: "ndarray" = None, 191 **kwargs, 192 ) -> None: 193 """ 194 Args: 195 center (Point): The center of the ellipse. 196 width (float): The width of the ellipse. 197 height (float): The height of the ellipse. 198 angle (float, optional): Rotation angle. Defaults to 0. 199 xform_matrix (ndarray, optional): Transformation matrix. Defaults to None. 200 **kwargs: Additional keyword arguments. 201 """ 202 n_points = defaults["n_ellipse_points"] 203 vertices = [ 204 tuple(p) 205 for p in ellipse_points(center, width / 2, height / 2, angle, n_points) 206 ] 207 super().__init__(vertices, closed=True, xform_matrix=xform_matrix, **kwargs) 208 a = width / 2 209 b = height / 2 210 self.a = a 211 self.b = b 212 self.center = center 213 self.width = width 214 self.height = height 215 self.angle = angle 216 self.smooth = True 217 self.closed = True 218 self.subtype = Types.ELLIPSE
Arguments:
- center (Point): The center of the ellipse.
- width (float): The width of the ellipse.
- height (float): The height of the ellipse.
- angle (float, optional): Rotation angle. Defaults to 0.
- xform_matrix (ndarray, optional): Transformation matrix. Defaults to None.
- **kwargs: Additional keyword arguments.
closed
220 @property 221 def closed(self): 222 """Return True ellipse is always closed. 223 224 Returns: 225 bool: Always returns True. 226 """ 227 return True
Return True ellipse is always closed.
Returns:
bool: Always returns True.
def
copy(self):
256 def copy(self): 257 """Return a copy of the ellipse. 258 259 Returns: 260 Ellipse: A copy of the ellipse. 261 """ 262 center = self.center 263 width = self.width 264 height = self.height 265 ellipse = Ellipse(center, width, height) 266 custom_attribs = custom_attributes(self) 267 for attrib in custom_attribs: 268 setattr(ellipse, attrib, getattr(self, attrib)) 269 270 return ellipse
Return a copy of the ellipse.
Returns:
Ellipse: A copy of the ellipse.
Inherited Members
def
ellipse_tangent(a, b, x, y, tol=0.001):
273def ellipse_tangent(a, b, x, y, tol=0.001): 274 """Calculates the angle of the tangent line to an ellipse at the point (x, y). 275 276 Args: 277 a (float): Semi-major axis of the ellipse. 278 b (float): Semi-minor axis of the ellipse. 279 x (float): x-coordinate of the point. 280 y (float): y-coordinate of the point. 281 tol (float, optional): Tolerance for point on ellipse check. Defaults to .001. 282 283 Returns: 284 float: Angle of the tangent line in radians. 285 """ 286 if abs((x**2 / a**2) + (y**2 / b**2) - 1) >= tol: 287 res = False 288 else: 289 # res = atan2(-(b**2 * x), (a**2 * y)) 290 res = atan2((b**2 * x), -(a**2 * y)) 291 292 return res
Calculates the angle of the tangent line to an ellipse at the point (x, y).
Arguments:
- a (float): Semi-major axis of the ellipse.
- b (float): Semi-minor axis of the ellipse.
- x (float): x-coordinate of the point.
- y (float): y-coordinate of the point.
- tol (float, optional): Tolerance for point on ellipse check. Defaults to .001.
Returns:
float: Angle of the tangent line in radians.
def
r_central(a, b, theta):
295def r_central(a, b, theta): 296 """Return the radius (distance between the center and the intersection point) 297 of the ellipse at the given angle. 298 299 Args: 300 a (float): Semi-major axis of the ellipse. 301 b (float): Semi-minor axis of the ellipse. 302 theta (float): Angle in radians. 303 304 Returns: 305 float: Radius at the given angle. 306 """ 307 return (a * b) / sqrt((b * cos(theta)) ** 2 + (a * sin(theta)) ** 2)
Return the radius (distance between the center and the intersection point) of the ellipse at the given angle.
Arguments:
- a (float): Semi-major axis of the ellipse.
- b (float): Semi-minor axis of the ellipse.
- theta (float): Angle in radians.
Returns:
float: Radius at the given angle.
def
ellipse_line_intersection(a, b, point):
310def ellipse_line_intersection(a, b, point): 311 """Return the intersection points of an ellipse and a line segment 312 connecting the given point to the ellipse center at (0, 0). 313 314 Args: 315 a (float): Semi-major axis of the ellipse. 316 b (float): Semi-minor axis of the ellipse. 317 point (tuple): Point coordinates (x, y). 318 319 Returns: 320 list: Intersection points. 321 """ 322 # adapted from http:# mathworld.wolfram.com/Ellipse-LineIntersection.html 323 # a, b is the ellipse width/2 and height/2 and (x_0, y_0) is the point 324 325 x_0, y_0 = point[:2] 326 x = ((a * b) / (sqrt(a**2 * y_0**2 + b**2 * x_0**2))) * x_0 327 y = ((a * b) / (sqrt(a**2 * y_0**2 + b**2 * x_0**2))) * y_0 328 329 return [(x, y), (-x, -y)]
Return the intersection points of an ellipse and a line segment connecting the given point to the ellipse center at (0, 0).
Arguments:
- a (float): Semi-major axis of the ellipse.
- b (float): Semi-minor axis of the ellipse.
- point (tuple): Point coordinates (x, y).
Returns:
list: Intersection points.
def
elliptic_arc_points(center, radius_x, radius_y, start_angle, span_angle, n_points=None):
332def elliptic_arc_points( 333 center, radius_x, radius_y, start_angle, span_angle, n_points=None 334): 335 """Generate points on an elliptic arc. 336 These are generated from the parametric equations of the ellipse. 337 They are not evenly spaced. 338 339 Args: 340 center (tuple): (x, y) coordinates of the ellipse center. 341 radius_x (float): Length of the semi-major axis. 342 radius_y (float): Length of the semi-minor axis. 343 start_angle (float): Starting angle of the arc. 344 span_angle (float): Span angle of the arc. 345 n_points (int): Number of points to generate. 346 347 Returns: 348 numpy.ndarray: Array of (x, y) coordinates of the ellipse points. 349 """ 350 rx = radius_x 351 if radius_y is None: 352 radius_y = radius_x 353 ry = radius_y 354 if n_points is None: 355 n = defaults["n_arc_points"] 356 n_points = ceil(n * abs(span_angle) / (2 * pi)) 357 start_angle = positive_angle(start_angle) 358 clockwise = span_angle < 0 359 if clockwise: 360 if start_angle + span_angle < 0: 361 end_angle = positive_angle(start_angle + span_angle) 362 t0 = get_ellipse_t_for_angle(end_angle, rx, ry) 363 t1 = get_ellipse_t_for_angle(2 * pi, rx, ry) 364 t = np.linspace(t0, t1, n_points) 365 x = center[0] + rx * np.cos(t) 366 y = center[1] + ry * np.sin(t) 367 slice1 = np.column_stack((x, y)) 368 t0 = get_ellipse_t_for_angle(0, rx, ry) 369 t1 = get_ellipse_t_for_angle(start_angle, rx, ry) 370 t = np.linspace(t0, t1, n_points) 371 x = center[0] + rx * np.cos(t) 372 y = center[1] + ry * np.sin(t) 373 slice2 = np.column_stack((x, y)) 374 res = np.flip(np.concatenate((slice1, slice2)), axis=0) 375 else: 376 end_angle = start_angle + span_angle 377 t0 = get_ellipse_t_for_angle(end_angle, rx, ry) 378 t1 = get_ellipse_t_for_angle(start_angle, rx, ry) 379 t = np.linspace(t0, t1, n_points) 380 x = center[0] + rx * np.cos(t) 381 y = center[1] + ry * np.sin(t) 382 res = np.flip(np.column_stack((x, y)), axis=0) 383 384 else: 385 if start_angle + span_angle > 2 * pi: 386 t0 = get_ellipse_t_for_angle(start_angle, rx, ry) 387 t1 = get_ellipse_t_for_angle(2 * pi, rx, ry) 388 t = np.linspace(t0, t1, n_points) 389 x = center[0] + rx * np.cos(t) 390 y = center[1] + ry * np.sin(t) 391 slice1 = np.column_stack((x, y)) 392 393 t0 = get_ellipse_t_for_angle(0, rx, ry) 394 t1 = get_ellipse_t_for_angle(start_angle + span_angle - 2 * pi, rx, ry) 395 t = np.linspace(t0, t1, n_points) 396 x = center[0] + rx * np.cos(t) 397 y = center[1] + ry * np.sin(t) 398 slice2 = np.column_stack((x, y)) 399 400 res = np.concatenate((slice1, slice2)) 401 else: 402 end_angle = start_angle + span_angle 403 t0 = get_ellipse_t_for_angle(start_angle, rx, ry) 404 t1 = get_ellipse_t_for_angle(end_angle, rx, ry) 405 t = np.linspace(t0, t1, n_points) 406 x = center[0] + rx * np.cos(t) 407 y = center[1] + ry * np.sin(t) 408 res = np.column_stack((x, y)) 409 410 return res
Generate points on an elliptic arc. These are generated from the parametric equations of the ellipse. They are not evenly spaced.
Arguments:
- center (tuple): (x, y) coordinates of the ellipse center.
- radius_x (float): Length of the semi-major axis.
- radius_y (float): Length of the semi-minor axis.
- start_angle (float): Starting angle of the arc.
- span_angle (float): Span angle of the arc.
- n_points (int): Number of points to generate.
Returns:
numpy.ndarray: Array of (x, y) coordinates of the ellipse points.
def
ellipse_points( center: Sequence[float], a: float, b: float, angle: float, n_points: int = None) -> numpy.ndarray:
413def ellipse_points( 414 center: Point, 415 a: float, 416 b: float, 417 angle: float, 418 n_points: int = None, 419) -> np.ndarray: 420 """Generate points on an ellipse. 421 These are generated from the parametric equations of the ellipse. 422 They are not evenly spaced. 423 424 Args: 425 center (tuple): (x, y) coordinates of the ellipse center. 426 a (float): Length of the semi-major axis. 427 b (float): Length of the semi-minor axis. 428 angle (float): Rotation angle of the ellipse. 429 n_points (int): Number of points to generate. 430 431 Returns: 432 numpy.ndarray: Array of (x, y) coordinates of the ellipse points. 433 """ 434 if n_points is None: 435 n_points = defaults["n_ellipse_points"] 436 437 t = np.linspace(0, 2 * pi, n_points) 438 x = center[0] + a * np.cos(t) 439 y = center[1] + b * np.sin(t) 440 441 points = homogenize(np.column_stack((x, y))) @ rotation_matrix(angle, center) 442 443 return points[:, :2].tolist()
Generate points on an ellipse. These are generated from the parametric equations of the ellipse. They are not evenly spaced.
Arguments:
- center (tuple): (x, y) coordinates of the ellipse center.
- a (float): Length of the semi-major axis.
- b (float): Length of the semi-minor axis.
- angle (float): Rotation angle of the ellipse.
- n_points (int): Number of points to generate.
Returns:
numpy.ndarray: Array of (x, y) coordinates of the ellipse points.
def
elliptic_arclength(t_0, t_1, a, b):
445def elliptic_arclength(t_0, t_1, a, b): 446 """Return the arclength of an ellipse between the given parametric angles. 447 The ellipse has semi-major axis a and semi-minor axis b. 448 449 Args: 450 t_0 (float): Starting parametric angle. 451 t_1 (float): Ending parametric angle. 452 a (float): Semi-major axis of the ellipse. 453 b (float): Semi-minor axis of the ellipse. 454 455 Returns: 456 float: Arclength of the ellipse. 457 """ 458 from scipy.special import ellipeinc # this takes too long to import 459 460 m = 1 - (b / a) ** 2 461 t1 = ellipeinc(t_1 - 0.5 * pi, m) 462 t0 = ellipeinc(t_0 - 0.5 * pi, m) 463 return a * (t1 - t0)
Return the arclength of an ellipse between the given parametric angles. The ellipse has semi-major axis a and semi-minor axis b.
Arguments:
- t_0 (float): Starting parametric angle.
- t_1 (float): Ending parametric angle.
- a (float): Semi-major axis of the ellipse.
- b (float): Semi-minor axis of the ellipse.
Returns:
float: Arclength of the ellipse.
def
central_to_parametric_angle(a, b, phi):
466def central_to_parametric_angle(a, b, phi): 467 """ 468 Converts a central angle to a parametric angle on an ellipse. 469 470 Args: 471 a (float): Semi-major axis of the ellipse. 472 b (float): Semi-minor axis of the ellipse. 473 phi (float): Angle of the line intersecting the center and the point. 474 475 Returns: 476 float: Parametric angle (in radians). 477 """ 478 t = atan2((a / b) * sin(phi), cos(phi)) 479 if t < 0: 480 t += 2 * pi 481 482 return t
Converts a central angle to a parametric angle on an ellipse.
Arguments:
- a (float): Semi-major axis of the ellipse.
- b (float): Semi-minor axis of the ellipse.
- phi (float): Angle of the line intersecting the center and the point.
Returns:
float: Parametric angle (in radians).
def
parametric_to_central_angle(a, b, t):
485def parametric_to_central_angle(a, b, t): 486 """ 487 Converts a parametric angle on an ellipse to a central angle. 488 489 Args: 490 a (float): Semi-major axis of the ellipse. 491 b (float): Semi-minor axis of the ellipse. 492 t (float): Parametric angle (in radians). 493 494 Returns: 495 float: Angle of the line intersecting the center and the point. 496 """ 497 phi = atan2((b / a) * sin(t), cos(t)) 498 if phi < 0: 499 phi += 2 * pi 500 501 return phi
Converts a parametric angle on an ellipse to a central angle.
Arguments:
- a (float): Semi-major axis of the ellipse.
- b (float): Semi-minor axis of the ellipse.
- t (float): Parametric angle (in radians).
Returns:
float: Angle of the line intersecting the center and the point.
def
ellipse_point(a, b, angle):
504def ellipse_point(a, b, angle): 505 """Return a point on an ellipse with the given a=width/2, b=height/2, and angle. 506 angle is the central-angle and in radians. 507 508 Args: 509 a (float): Semi-major axis of the ellipse. 510 b (float): Semi-minor axis of the ellipse. 511 angle (float): Central angle in radians. 512 513 Returns: 514 tuple: Coordinates of the point on the ellipse. 515 """ 516 r = r_central(a, b, angle) 517 518 return (r * cos(angle), r * sin(angle))
Return a point on an ellipse with the given a=width/2, b=height/2, and angle. angle is the central-angle and in radians.
Arguments:
- a (float): Semi-major axis of the ellipse.
- b (float): Semi-minor axis of the ellipse.
- angle (float): Central angle in radians.
Returns:
tuple: Coordinates of the point on the ellipse.
def
ellipse_param_point(a, b, t):
521def ellipse_param_point(a, b, t): 522 """Return a point on an ellipse with the given a=width/2, b=height/2, and parametric angle. 523 t is the parametric angle and in radians. 524 525 Args: 526 a (float): Semi-major axis of the ellipse. 527 b (float): Semi-minor axis of the ellipse. 528 t (float): Parametric angle in radians. 529 530 Returns: 531 tuple: Coordinates of the point on the ellipse. 532 """ 533 return (a * cos(t), b * sin(t))
Return a point on an ellipse with the given a=width/2, b=height/2, and parametric angle. t is the parametric angle and in radians.
Arguments:
- a (float): Semi-major axis of the ellipse.
- b (float): Semi-minor axis of the ellipse.
- t (float): Parametric angle in radians.
Returns:
tuple: Coordinates of the point on the ellipse.
def
get_ellipse_t_for_angle(angle, a, b):
536def get_ellipse_t_for_angle(angle, a, b): 537 """ 538 Calculates the parameter t for a given angle on an ellipse. 539 540 Args: 541 angle (float): The angle in radians. 542 a (float): Semi-major axis of the ellipse. 543 b (float): Semi-minor axis of the ellipse. 544 545 Returns: 546 float: The parameter t. 547 """ 548 t = atan2(a * sin(angle), b * cos(angle)) 549 if t < 0: 550 t += 2 * pi 551 return t
Calculates the parameter t for a given angle on an ellipse.
Arguments:
- angle (float): The angle in radians.
- a (float): Semi-major axis of the ellipse.
- b (float): Semi-minor axis of the ellipse.
Returns:
float: The parameter t.
def
ellipse_central_angle(t, a, b):
554def ellipse_central_angle(t, a, b): 555 """ 556 Calculates the central angle of an ellipse for a given parameter t. 557 558 Args: 559 t (float): The parameter value. 560 a (float): The semi-major axis of the ellipse. 561 b (float): The semi-minor axis of the ellipse. 562 563 Returns: 564 float: The central angle in radians. 565 """ 566 theta = atan2(a * sin(t), b * cos(t)) 567 568 return theta
Calculates the central angle of an ellipse for a given parameter t.
Arguments:
- t (float): The parameter value.
- a (float): The semi-major axis of the ellipse.
- b (float): The semi-minor axis of the ellipse.
Returns:
float: The central angle in radians.
def
ellipse_intersection(x1, y1, a, b, phi, x2, y2, c, d, phi2):
571def ellipse_intersection(x1, y1, a, b, phi, x2, y2, c, d, phi2): 572 """Calculate the intersection points of two ellipses. 573 The ellipses are defined by their center, radii, and rotation angle.""" 574 # Taken from https:# github.com/VoyakaGOD/intersection-of-two-ellipses/blob/master/geometry.js 575 phi0 = phi2 - phi 576 577 p0 = rotate_point(((x2 - x1), (y2 - y1)), -phi) 578 579 # cos(t1) = A + Bcos(t2) + Csin(t2) 580 # sin(t1) = D + Ecos(t2) + Fsin(t2) 581 A = p0[0] / a 582 B = c * cos(phi0) / a 583 C = -d * sin(phi0) / a 584 D = p0[1] / b 585 E = c * sin(phi0) / b 586 F = d * cos(phi0) / b 587 588 # Gx^2 + Hx + Ixy + Jy + K = 0, x = cos(t2), y = sin(t2) 589 G = B * B + E * E - C * C - F * F 590 H = 2 * (A * B + D * E) 591 I = 2 * (B * C + E * F) 592 J = 2 * (A * C + D * F) 593 K = A * A + D * D + C * C + F * F - 1 594 595 roots = [] 596 L = G * G + I * I 597 if isclose(L, 0, 0, 1e-7): 598 # Gx + Hy + K = 0 599 roots = solve_quadratic_eq(G, H, K) 600 601 elif isclose(I, 0, 1e-7): 602 # Gx^2 + Hx + K = 0 603 roots = solve_quadratic_eq(G, H, K) 604 605 elif isclose(G, 0, 1e-7): 606 # Hx + Jy + K = 0 607 roots = solve_quadratic_eq(H * H + J * J, 2 * K * H, K * K - J * J) 608 609 else: 610 # Lx^4 + Mx^3 + Nx^2 + Ox + P = 0 611 M = 2 * (G * H + I * J) 612 N = H * H + 2 * G * K + J * J - I * I 613 O = 2 * (K * H - I * J) 614 P = K * K - J * J 615 616 iL = 1 / L 617 roots = solve_quartic_equation(M * iL, N * iL, O * iL, P * iL) 618 619 points = [] 620 # for(i = 0 i < roots.length i++) 621 # for i in range(len(roots)): 622 for i, x in enumerate(roots): 623 # x = roots[i] 624 if isclose(I * x + J, 0, 1e-7): 625 y = sqrt(1 - x * x) 626 points.append((x, y)) 627 628 if not isclose(y, 0, 1e-7): 629 points.append((x, -y)) 630 else: 631 y = -(G * x * x + H * x + K) / (I * x + J) 632 if abs(y) < 1e-2: 633 y = sqrt(1 - x * x) 634 points.append((x, y)) 635 points.append((x, -y)) 636 else: 637 points.append((x, y)) 638 639 for i, pnt in enumerate(points): 640 x, y = pnt 641 points[i] = (A + B * x + C * y) * a, (D + E * x + F * y) * b 642 643 for i, pnt in enumerate(points): 644 x, y = pnt 645 x, y = rotate_point((x, y), phi) 646 points[i] = (x + x1, y + y1) 647 648 return points
Calculate the intersection points of two ellipses. The ellipses are defined by their center, radii, and rotation angle.
def
inverse_complex_number(z: complex) -> complex:
655def inverse_complex_number(z: complex) -> complex: 656 """ 657 Calculate the inverse of a complex number. 658 659 Args: 660 z (complex): The complex number. 661 662 Returns: 663 complex: The inverse of the complex number. 664 """ 665 a = z.real 666 b = z.imag 667 668 if a == 0 and b == 0: 669 raise ZeroDivisionError("Cannot calculate the inverse of 0") 670 671 return complex(a / (a**2 + b**2), -b / (a**2 + b**2))
Calculate the inverse of a complex number.
Arguments:
- z (complex): The complex number.
Returns:
complex: The inverse of the complex number.
def
Re(num):
def
Im(num):
def
Sqrt(complex_):
def
Qbrt(complex_):
689def Qbrt(complex_): 690 # Taken from https:# github.com/VoyakaGOD/intersection-of-two-ellipses/blob/master/quartic.js 691 692 angle = cmath.phase(complex_) * 0.33333333333 693 # angle = atan2(complex.y, complex.x) * 0.33333333333 694 l = pow( 695 complex_.real * complex_.real + complex_.imag * complex_.imag, 0.16666666666 696 ) 697 # return Complex(l * cos(angle), l * sin(angle)) 698 return complex(l * cos(angle), l * sin(angle))
def
solve_complex_quadratic_equation(b, c):
def
get_one_cubic_equation_root(a, b, c):
711def get_one_cubic_equation_root(a, b, c): 712 # Taken from https:# github.com/VoyakaGOD/intersection-of-two-ellipses/blob/master/quartic.js 713 714 p = b - a * a * 0.33333333333 715 q = (2 / 27) * a * a * a - a * b * 0.33333333333 + c 716 sqrt_Q = cmath.sqrt(Re(0.03703703703 * p * p * p + 0.25 * q * q)) 717 A = Qbrt(Re(-0.5 * q) - sqrt_Q) 718 if A.real == 0 and A.imag == 0: 719 A = Qbrt(Re(-0.5 * q) + sqrt_Q) 720 B = inverse_complex_number(A) * (-p * 0.33333333333) 721 722 return (A + B) - (Re(a * 0.33333333333))
def
solve_quartic_equation(a, b, c, d):
725def solve_quartic_equation(a, b, c, d): 726 # Taken from https:# github.com/VoyakaGOD/intersection-of-two-ellipses/blob/master/geometry.js 727 728 a2 = a * a 729 a3 = a2 * a 730 a4 = a3 * a 731 p = b - (3/8) * a2 732 q = (1 / 8) * a3 - 0.5 * a * b + c 733 r = d - 0.25 * a * c + (1 / 16) * a2 * b - (3 / 256) * a4 734 735 result = [] 736 737 if(isclose(q, 0, 1e-7)): 738 D = p * p - 4 * r 739 if(abs(D) < 1e-5): 740 m = -0.5 * p 741 if (m >= 0): 742 result.append(sqrt(m)) 743 if (m > 0): 744 result.append(-sqrt(m)) 745 elif (D > 0): 746 sqrt_D = sqrt(D) 747 m1 = (-p - sqrt_D) * 0.5 748 m2 = (-p + sqrt_D) * 0.5 749 sqrt_m1 = sqrt(m1) 750 sqrt_m2 = sqrt(m2) 751 if (m1 >= 0): 752 result.append(sqrt_m1) 753 if (m1 > 0): 754 result.append(-sqrt_m1) 755 if (m2 >= 0): 756 result.append(sqrt_m2) 757 if (m2 > 0): 758 result.append(-sqrt_m2) 759 else: 760 t = get_one_cubic_equation_root(2 * p, p * p-4 * r, -q * q) 761 z = cmath.sqrt(t) 762 u = (Re(p)+(t)) * 0.5 763 v = inverse_complex_number(z) * (q * 0.5) 764 x12 = solve_complex_quadratic_equation(z, u - v) 765 x34 = solve_complex_quadratic_equation(-z, u + v) 766 767 if (isclose(x12[0].imag, 0, 1e-7)): 768 result.append(x12[0].real) 769 if (isclose(x12[1].imag, 0, 1e-7)): 770 result.append(x12[1].real) 771 if (isclose(x34[0].imag, 0, 1e-7)): 772 result.append(x34[0].real) 773 if (isclose(x34[1].imag, 0, 1e-7)): 774 result.append(x34[1].real) 775 776 # for(i = 0 i < result.length i++) 777 for i in range(len(result)): 778 result[i] -= 0.25 * a 779 780 return result