1#cython: cdivision=True 2#cython: boundscheck=False 3#cython: nonecheck=False 4#cython: wraparound=False 5import math 6import numpy as np 7 8cimport numpy as cnp 9from libc.math cimport sqrt, sin, cos, floor, ceil, fabs 10from .._shared.geometry cimport point_in_polygon 11 12cnp.import_array() 13 14 15def _coords_inside_image(rr, cc, shape, val=None): 16 """ 17 Return the coordinates inside an image of a given shape. 18 19 Parameters 20 ---------- 21 rr, cc : (N,) ndarray of int 22 Indices of pixels. 23 shape : tuple 24 Image shape which is used to determine the maximum extent of output 25 pixel coordinates. Must be at least length 2. Only the first two values 26 are used to determine the extent of the input image. 27 val : (N, D) ndarray of float, optional 28 Values of pixels at coordinates ``[rr, cc]``. 29 30 Returns 31 ------- 32 rr, cc : (M,) array of int 33 Row and column indices of valid pixels (i.e. those inside `shape`). 34 val : (M, D) array of float, optional 35 Values at `rr, cc`. Returned only if `val` is given as input. 36 """ 37 mask = (rr >= 0) & (rr < shape[0]) & (cc >= 0) & (cc < shape[1]) 38 if val is None: 39 return rr[mask], cc[mask] 40 else: 41 return rr[mask], cc[mask], val[mask] 42 43 44def _line(Py_ssize_t r0, Py_ssize_t c0, Py_ssize_t r1, Py_ssize_t c1): 45 """Generate line pixel coordinates. 46 47 Parameters 48 ---------- 49 r0, c0 : int 50 Starting position (row, column). 51 r1, c1 : int 52 End position (row, column). 53 54 Returns 55 ------- 56 rr, cc : (N,) ndarray of int 57 Indices of pixels that belong to the line. 58 May be used to directly index into an array, e.g. 59 ``img[rr, cc] = 1``. 60 61 See Also 62 -------- 63 line_aa : Anti-aliased line generator 64 """ 65 66 cdef char steep = 0 67 cdef Py_ssize_t r = r0 68 cdef Py_ssize_t c = c0 69 cdef Py_ssize_t dr = abs(r1 - r0) 70 cdef Py_ssize_t dc = abs(c1 - c0) 71 cdef Py_ssize_t sr, sc, d, i 72 73 cdef Py_ssize_t[::1] rr = np.zeros(max(dc, dr) + 1, dtype=np.intp) 74 cdef Py_ssize_t[::1] cc = np.zeros(max(dc, dr) + 1, dtype=np.intp) 75 76 with nogil: 77 if (c1 - c) > 0: 78 sc = 1 79 else: 80 sc = -1 81 if (r1 - r) > 0: 82 sr = 1 83 else: 84 sr = -1 85 if dr > dc: 86 steep = 1 87 c, r = r, c 88 dc, dr = dr, dc 89 sc, sr = sr, sc 90 d = (2 * dr) - dc 91 92 for i in range(dc): 93 if steep: 94 rr[i] = c 95 cc[i] = r 96 else: 97 rr[i] = r 98 cc[i] = c 99 while d >= 0: 100 r = r + sr 101 d = d - (2 * dc) 102 c = c + sc 103 d = d + (2 * dr) 104 105 rr[dc] = r1 106 cc[dc] = c1 107 108 return np.asarray(rr), np.asarray(cc) 109 110 111def _line_aa(Py_ssize_t r0, Py_ssize_t c0, Py_ssize_t r1, Py_ssize_t c1): 112 """Generate anti-aliased line pixel coordinates. 113 114 Parameters 115 ---------- 116 r0, c0 : int 117 Starting position (row, column). 118 r1, c1 : int 119 End position (row, column). 120 121 Returns 122 ------- 123 rr, cc, val : (N,) ndarray (int, int, float) 124 Indices of pixels (`rr`, `cc`) and intensity values (`val`). 125 ``img[rr, cc] = val``. 126 127 References 128 ---------- 129 .. [1] A Rasterizing Algorithm for Drawing Curves, A. Zingl, 2012 130 http://members.chello.at/easyfilter/Bresenham.pdf 131 """ 132 cdef list rr = list() 133 cdef list cc = list() 134 cdef list val = list() 135 136 cdef int dc = abs(c0 - c1) 137 cdef int dc_prime 138 139 cdef int dr = abs(r0 - r1) 140 cdef float err = dc - dr 141 cdef float err_prime 142 143 cdef int c, r, sign_c, sign_r 144 cdef float ed 145 146 if c0 < c1: 147 sign_c = 1 148 else: 149 sign_c = -1 150 151 if r0 < r1: 152 sign_r = 1 153 else: 154 sign_r = -1 155 156 if dc + dr == 0: 157 ed = 1 158 else: 159 ed = sqrt(dc*dc + dr*dr) 160 161 c, r = c0, r0 162 while True: 163 cc.append(c) 164 rr.append(r) 165 val.append(fabs(err - dc + dr) / ed) 166 167 err_prime = err 168 c_prime = c 169 170 if (2 * err_prime) >= -dc: 171 if c == c1: 172 break 173 if (err_prime + dr) < ed: 174 cc.append(c) 175 rr.append(r + sign_r) 176 val.append(fabs(err_prime + dr) / ed) 177 err -= dr 178 c += sign_c 179 180 if 2 * err_prime <= dr: 181 if r == r1: 182 break 183 if (dc - err_prime) < ed: 184 cc.append(c_prime + sign_c) 185 rr.append(r) 186 val.append(fabs(dc - err_prime) / ed) 187 err += dc 188 r += sign_r 189 190 return (np.array(rr, dtype=np.intp), 191 np.array(cc, dtype=np.intp), 192 1. - np.array(val, dtype=float)) 193 194 195def _polygon(r, c, shape): 196 """Generate coordinates of pixels within polygon. 197 198 Parameters 199 ---------- 200 r : (N,) ndarray 201 Row coordinates of vertices of polygon. 202 c : (N,) ndarray 203 Column coordinates of vertices of polygon. 204 shape : tuple 205 Image shape which is used to determine the maximum extent of output 206 pixel coordinates. This is useful for polygons that exceed the image 207 size. If None, the full extent of the polygon is used. 208 209 Returns 210 ------- 211 rr, cc : ndarray of int 212 Pixel coordinates of polygon. 213 May be used to directly index into an array, e.g. 214 ``img[rr, cc] = 1``. 215 """ 216 r = np.atleast_1d(r) 217 c = np.atleast_1d(c) 218 219 cdef Py_ssize_t nr_verts = c.shape[0] 220 cdef Py_ssize_t minr = int(max(0, r.min())) 221 cdef Py_ssize_t maxr = int(ceil(r.max())) 222 cdef Py_ssize_t minc = int(max(0, c.min())) 223 cdef Py_ssize_t maxc = int(ceil(c.max())) 224 225 # make sure output coordinates do not exceed image size 226 if shape is not None: 227 maxr = min(shape[0] - 1, maxr) 228 maxc = min(shape[1] - 1, maxc) 229 230 # make contiguous arrays for r, c coordinates 231 cdef cnp.float64_t[::1] rptr = np.ascontiguousarray(r, 'float64') 232 cdef cnp.float64_t[::1] cptr = np.ascontiguousarray(c, 'float64') 233 cdef cnp.float64_t r_i, c_i 234 235 # output coordinate arrays 236 rr = list() 237 cc = list() 238 239 for r_i in range(minr, maxr+1): 240 for c_i in range(minc, maxc+1): 241 if point_in_polygon(cptr, rptr, c_i, r_i): 242 rr.append(r_i) 243 cc.append(c_i) 244 245 return np.array(rr, dtype=np.intp), np.array(cc, dtype=np.intp) 246 247 248def _circle_perimeter(Py_ssize_t r_o, Py_ssize_t c_o, Py_ssize_t radius, 249 method, shape): 250 """Generate circle perimeter coordinates. 251 252 Parameters 253 ---------- 254 r_o, c_o : int 255 Centre coordinate of circle. 256 radius : int 257 Radius of circle. 258 method : {'bresenham', 'andres'} 259 bresenham : Bresenham method (default) 260 andres : Andres method 261 shape : tuple 262 Image shape which is used to determine the maximum extent of output pixel 263 coordinates. This is useful for circles that exceed the image size. 264 If None, the full extent of the circle is used. 265 266 Returns 267 ------- 268 rr, cc : (N,) ndarray of int 269 Bresenham and Andres' method: 270 Indices of pixels that belong to the circle perimeter. 271 May be used to directly index into an array, e.g. 272 ``img[rr, cc] = 1``. 273 274 Notes 275 ----- 276 Andres method presents the advantage that concentric 277 circles create a disc whereas Bresenham can make holes. There 278 is also less distortions when Andres circles are rotated. 279 Bresenham method is also known as midpoint circle algorithm. 280 Anti-aliased circle generator is available with `circle_perimeter_aa`. 281 282 References 283 ---------- 284 .. [1] J.E. Bresenham, "Algorithm for computer control of a digital 285 plotter", IBM Systems journal, 4 (1965) 25-30. 286 .. [2] E. Andres, "Discrete circles, rings and spheres", Computers & 287 Graphics, 18 (1994) 695-706. 288 """ 289 290 cdef list rr = list() 291 cdef list cc = list() 292 293 cdef Py_ssize_t c = 0 294 cdef Py_ssize_t r = radius 295 cdef Py_ssize_t d = 0 296 297 cdef double dceil = 0 298 cdef double dceil_prev = 0 299 300 cdef char cmethod 301 if method == 'bresenham': 302 d = 3 - 2 * radius 303 cmethod = b'b' 304 elif method == 'andres': 305 d = radius - 1 306 cmethod = b'a' 307 else: 308 raise ValueError('Wrong method') 309 310 while r >= c: 311 rr.extend([r, -r, r, -r, c, -c, c, -c]) 312 cc.extend([c, c, -c, -c, r, r, -r, -r]) 313 314 if cmethod == b'b': 315 if d < 0: 316 d += 4 * c + 6 317 else: 318 d += 4 * (c - r) + 10 319 r -= 1 320 c += 1 321 elif cmethod == b'a': 322 if d >= 2 * (c - 1): 323 d = d - 2 * c 324 c = c + 1 325 elif d <= 2 * (radius - r): 326 d = d + 2 * r - 1 327 r = r - 1 328 else: 329 d = d + 2 * (r - c - 1) 330 r = r - 1 331 c = c + 1 332 333 if shape is not None: 334 return _coords_inside_image(np.array(rr, dtype=np.intp) + r_o, 335 np.array(cc, dtype=np.intp) + c_o, 336 shape) 337 return (np.array(rr, dtype=np.intp) + r_o, 338 np.array(cc, dtype=np.intp) + c_o) 339 340 341def _circle_perimeter_aa(Py_ssize_t r_o, Py_ssize_t c_o, 342 Py_ssize_t radius, shape): 343 """Generate anti-aliased circle perimeter coordinates. 344 345 Parameters 346 ---------- 347 r_o, c_o : int 348 Centre coordinate of circle. 349 radius : int 350 Radius of circle. 351 shape : tuple 352 Image shape which is used to determine the maximum extent of output 353 pixel coordinates. This is useful for circles that exceed the image 354 size. If None, the full extent of the circle is used. 355 356 Returns 357 ------- 358 rr, cc, val : (N,) ndarray (int, int, float) 359 Indices of pixels (`rr`, `cc`) and intensity values (`val`). 360 ``img[rr, cc] = val``. 361 362 Notes 363 ----- 364 Wu's method draws anti-aliased circle. This implementation doesn't use 365 lookup table optimization. 366 367 References 368 ---------- 369 .. [1] X. Wu, "An efficient antialiasing technique", In ACM SIGGRAPH 370 Computer Graphics, 25 (1991) 143-152. 371 """ 372 373 cdef Py_ssize_t c = 0 374 cdef Py_ssize_t r = radius 375 cdef Py_ssize_t d = 0 376 377 cdef double dceil = 0 378 cdef double dceil_prev = 0 379 380 cdef list rr = [r, c, r, c, -r, -c, -r, -c] 381 cdef list cc = [c, r, -c, -r, c, r, -c, -r] 382 cdef list val = [1] * 8 383 384 while r > c + 1: 385 c += 1 386 dceil = sqrt(radius * radius - c * c) 387 dceil = ceil(dceil) - dceil 388 if dceil < dceil_prev: 389 r -= 1 390 rr.extend([r, r - 1, c, c, r, r - 1, c, c]) 391 cc.extend([c, c, r, r - 1, -c, -c, -r, 1 - r]) 392 393 rr.extend([-r, 1 - r, -c, -c, -r, 1 - r, -c, -c]) 394 cc.extend([c, c, r, r - 1, -c, -c, -r, 1 - r]) 395 396 val.extend([1 - dceil, dceil] * 8) 397 dceil_prev = dceil 398 399 if shape is not None: 400 return _coords_inside_image(np.array(rr, dtype=np.intp) + r_o, 401 np.array(cc, dtype=np.intp) + c_o, 402 shape, 403 val=np.array(val, dtype=float)) 404 return (np.array(rr, dtype=np.intp) + r_o, 405 np.array(cc, dtype=np.intp) + c_o, 406 np.array(val, dtype=float)) 407 408 409def _ellipse_perimeter(Py_ssize_t r_o, Py_ssize_t c_o, Py_ssize_t r_radius, 410 Py_ssize_t c_radius, double orientation, shape): 411 """Generate ellipse perimeter coordinates. 412 413 Parameters 414 ---------- 415 r_o, c_o : int 416 Centre coordinate of ellipse. 417 r_radius, c_radius : int 418 Minor and major semi-axes. ``(r/r_radius)**2 + (c/c_radius)**2 = 1``. 419 orientation : double 420 Major axis orientation in clockwise direction as radians. 421 shape : tuple 422 Image shape which is used to determine the maximum extent of output pixel 423 coordinates. This is useful for ellipses that exceed the image size. 424 If None, the full extent of the ellipse is used. 425 426 Returns 427 ------- 428 rr, cc : (N,) ndarray of int 429 Indices of pixels that belong to the ellipse perimeter. 430 May be used to directly index into an array, e.g. 431 ``img[rr, cc] = 1``. 432 433 References 434 ---------- 435 .. [1] A Rasterizing Algorithm for Drawing Curves, A. Zingl, 2012 436 http://members.chello.at/easyfilter/Bresenham.pdf 437 """ 438 439 # If both radii == 0, return the center to avoid infinite loop in 2nd set 440 if r_radius == 0 and c_radius == 0: 441 return np.array(r_o), np.array(c_o) 442 443 # Pixels 444 cdef list rr = list() 445 cdef list cc = list() 446 447 # Compute useful values 448 cdef Py_ssize_t rd = r_radius * r_radius 449 cdef Py_ssize_t cd = c_radius * c_radius 450 451 cdef Py_ssize_t r, c, e2, err 452 453 cdef int ir0, ir1, ic0, ic1, ird, icd 454 cdef double sin_angle, ra, ca, za, a, b 455 456 if orientation == 0: 457 c = -c_radius 458 r = 0 459 e2 = rd 460 err = c * (2 * e2 + c) + e2 461 while c <= 0: 462 # Quadrant 1 463 rr.append(r_o + r) 464 cc.append(c_o - c) 465 # Quadrant 2 466 rr.append(r_o + r) 467 cc.append(c_o + c) 468 # Quadrant 3 469 rr.append(r_o - r) 470 cc.append(c_o + c) 471 # Quadrant 4 472 rr.append(r_o - r) 473 cc.append(c_o - c) 474 # Adjust `r` and `c` 475 e2 = 2 * err 476 if e2 >= (2 * c + 1) * rd: 477 c += 1 478 err += (2 * c + 1) * rd 479 if e2 <= (2 * r + 1) * cd: 480 r += 1 481 err += (2 * r + 1) * cd 482 while r < r_radius: 483 r += 1 484 rr.append(r_o + r) 485 cc.append(c_o) 486 rr.append(r_o - r) 487 cc.append(c_o) 488 489 else: 490 sin_angle = sin(orientation) 491 za = (cd - rd) * sin_angle 492 ca = sqrt(cd - za * sin_angle) 493 ra = sqrt(rd + za * sin_angle) 494 495 a = ca + 0.5 496 b = ra + 0.5 497 za = za * a * b / (ca * ra) 498 499 ir0 = int(r_o - b) 500 ic0 = int(c_o - a) 501 ir1 = int(r_o + b) 502 ic1 = int(c_o + a) 503 504 ca = ic1 - ic0 505 ra = ir1 - ir0 506 za = 4 * za * cos(orientation) 507 w = ca * ra 508 if w != 0: 509 w = (w - za) / (w + w) 510 icd = int(floor(ca * w + 0.5)) 511 ird = int(floor(ra * w + 0.5)) 512 513 # Draw the 4 quadrants 514 rr_t, cc_t = _bezier_segment(ir0 + ird, ic0, ir0, ic0, ir0, ic0 + icd, 1-w) 515 rr.extend(rr_t) 516 cc.extend(cc_t) 517 rr_t, cc_t = _bezier_segment(ir0 + ird, ic0, ir1, ic0, ir1, ic1 - icd, w) 518 rr.extend(rr_t) 519 cc.extend(cc_t) 520 rr_t, cc_t = _bezier_segment(ir1 - ird, ic1, ir1, ic1, ir1, ic1 - icd, 1-w) 521 rr.extend(rr_t) 522 cc.extend(cc_t) 523 rr_t, cc_t = _bezier_segment(ir1 - ird, ic1, ir0, ic1, ir0, ic0 + icd, w) 524 rr.extend(rr_t) 525 cc.extend(cc_t) 526 527 if shape is not None: 528 return _coords_inside_image(np.array(rr, dtype=np.intp), 529 np.array(cc, dtype=np.intp), shape) 530 return np.array(rr, dtype=np.intp), np.array(cc, dtype=np.intp) 531 532 533def _bezier_segment(Py_ssize_t r0, Py_ssize_t c0, 534 Py_ssize_t r1, Py_ssize_t c1, 535 Py_ssize_t r2, Py_ssize_t c2, 536 double weight): 537 """Generate Bezier segment coordinates. 538 539 Parameters 540 ---------- 541 r0, c0 : int 542 Coordinates of the first control point. 543 r1, c1 : int 544 Coordinates of the middle control point. 545 r2, c2 : int 546 Coordinates of the last control point. 547 weight : double 548 Middle control point weight, it describes the line tension. 549 550 Returns 551 ------- 552 rr, cc : (N,) ndarray of int 553 Indices of pixels that belong to the Bezier curve. 554 May be used to directly index into an array, e.g. 555 ``img[rr, cc] = 1``. 556 557 Notes 558 ----- 559 The algorithm is the rational quadratic algorithm presented in 560 reference [1]_. 561 562 References 563 ---------- 564 .. [1] A Rasterizing Algorithm for Drawing Curves, A. Zingl, 2012 565 http://members.chello.at/easyfilter/Bresenham.pdf 566 """ 567 # Pixels 568 cdef list cc = list() 569 cdef list rr = list() 570 571 # Steps 572 cdef double sc = c2 - c1 573 cdef double sr = r2 - r1 574 575 cdef double d2c = c0 - c2 576 cdef double d2r = r0 - r2 577 cdef double d1c = c0 - c1 578 cdef double d1r = r0 - r1 579 cdef double rc = d1c * sr + d1r * sc 580 cdef double cur = d1c * sr - d1r * sc 581 cdef double err 582 583 cdef bint test1, test2 584 585 # If not a straight line 586 if cur != 0 and weight > 0: 587 if (sc * sc + sr * sr > d1c * d1c + d1r * d1r): 588 # Swap point 0 and point 2 589 # to start from the longer part 590 c2 = c0 591 c0 -= <Py_ssize_t>(d2c) 592 r2 = r0 593 r0 -= <Py_ssize_t>(d2r) 594 cur = -cur 595 d1c = 2 * (4 * weight * sc * d1c + d2c * d2c) 596 d1r = 2 * (4 * weight * sr * d1r + d2r * d2r) 597 # Set steps 598 if c0 < c2: 599 sc = 1 600 else: 601 sc = -1 602 if r0 < r2: 603 sr = 1 604 else: 605 sr = -1 606 rc = -2 * sc * sr * (2 * weight * rc + d2c * d2r) 607 608 if cur * sc * sr < 0: 609 d1c = -d1c 610 d1r = -d1r 611 rc = -rc 612 cur = -cur 613 614 d2c = 4 * weight * (c1 - c0) * sr * cur + d1c / 2 + rc 615 d2r = 4 * weight * (r0 - r1) * sc * cur + d1r / 2 + rc 616 617 # Flat ellipse, algo fails 618 if weight < 0.5 and (d2r > rc or d2c < rc): 619 cur = (weight + 1) / 2 620 weight = sqrt(weight) 621 rc = 1. / (weight + 1) 622 # Subdivide curve in half 623 sc = floor((c0 + 2 * weight * c1 + c2) * rc * 0.5 + 0.5) 624 sr = floor((r0 + 2 * weight * r1 + r2) * rc * 0.5 + 0.5) 625 d2c = floor((weight * c1 + c0) * rc + 0.5) 626 d2r = floor((r1 * weight + r0) * rc + 0.5) 627 return _bezier_segment(r0, c0, <Py_ssize_t>(d2r), <Py_ssize_t>(d2c), 628 <Py_ssize_t>(sr), <Py_ssize_t>(sc), cur) 629 630 err = d2c + d2r - rc 631 while d2r <= rc and d2c >= rc: 632 cc.append(c0) 633 rr.append(r0) 634 if c0 == c2 and r0 == r2: 635 # The job is done! 636 return np.array(rr, dtype=np.intp), np.array(cc, dtype=np.intp) 637 638 # Save boolean values 639 test1 = 2 * err > d2r 640 test2 = 2 * (err + d1r) < -d2r 641 # Move (c0, r0) to the next position 642 if 2 * err < d2c or test2: 643 r0 += <Py_ssize_t>(sr) 644 d2r += rc 645 d2c += d1c 646 err += d2c 647 if 2 * err > d2c or test1: 648 c0 += <Py_ssize_t>(sc) 649 d2c += rc 650 d2r += d1r 651 err += d2r 652 653 # Plot line 654 cc_t, rr_t = _line(c0, r0, c2, r2) 655 cc.extend(cc_t) 656 rr.extend(rr_t) 657 658 return np.array(rr, dtype=np.intp), np.array(cc, dtype=np.intp) 659 660 661def _bezier_curve(Py_ssize_t r0, Py_ssize_t c0, 662 Py_ssize_t r1, Py_ssize_t c1, 663 Py_ssize_t r2, Py_ssize_t c2, 664 double weight, shape): 665 """Generate Bezier curve coordinates. 666 667 Parameters 668 ---------- 669 r0, c0 : int 670 Coordinates of the first control point. 671 r1, c1 : int 672 Coordinates of the middle control point. 673 r2, c2 : int 674 Coordinates of the last control point. 675 weight : double 676 Middle control point weight, it describes the line tension. 677 shape : tuple 678 Image shape which is used to determine the maximum extent of output 679 pixel coordinates. This is useful for curves that exceed the image 680 size. If None, the full extent of the curve is used. 681 682 Returns 683 ------- 684 rr, cc : (N,) ndarray of int 685 Indices of pixels that belong to the Bezier curve. 686 May be used to directly index into an array, e.g. 687 ``img[rr, cc] = 1``. 688 689 Notes 690 ----- 691 The algorithm is the rational quadratic algorithm presented in 692 reference [1]_. 693 694 References 695 ---------- 696 .. [1] A Rasterizing Algorithm for Drawing Curves, A. Zingl, 2012 697 http://members.chello.at/easyfilter/Bresenham.pdf 698 """ 699 # Pixels 700 cdef list cc = list() 701 cdef list rr = list() 702 703 cdef int vc, vr 704 cdef double dc, dr, ww, t, q 705 vc = c0 - 2 * c1 + c2 706 vr = r0 - 2 * r1 + r2 707 708 dc = c0 - c1 709 dr = r0 - r1 710 711 if dc * (c2 - c1) > 0: 712 if dr * (r2 - r1): 713 if abs(dc * vr) > abs(dr * vc): 714 c0 = c2 715 c2 = <Py_ssize_t>(dc + c1) 716 r0 = r2 717 r2 = <Py_ssize_t>(dr + r1) 718 if (c0 == c2) or (weight == 1.): 719 t = <double>(c0 - c1) / vc 720 else: 721 q = sqrt(4. * weight * weight * (c0 - c1) * (c2 - c1) + (c2 - c0) * floor(c2 - c0)) 722 if (c1 < c0): 723 q = -q 724 t = (2. * weight * (c0 - c1) - c0 + c2 + q) / (2. * (1. - weight) * (c2 - c0)) 725 726 q = 1. / (2. * t * (1. - t) * (weight - 1.) + 1.0) 727 dc = (t * t * (c0 - 2. * weight * c1 + c2) + 2. * t * (weight * c1 - c0) + c0) * q 728 dr = (t * t * (r0 - 2. * weight * r1 + r2) + 2. * t * (weight * r1 - r0) + r0) * q 729 ww = t * (weight - 1.) + 1. 730 ww *= ww * q 731 weight = ((1. - t) * (weight - 1.) + 1.) * sqrt(q) 732 vc = <int>(dc + 0.5) 733 vr = <int>(dr + 0.5) 734 dr = (dc - c0) * (r1 - r0) / (c1 - c0) + r0 735 736 rr_t, cc_t = _bezier_segment(r0, c0, <int>(dr + 0.5), vc, vr, vc, ww) 737 cc.extend(cc_t) 738 rr.extend(rr_t) 739 740 dr = (dc - c2) * (r1 - r2) / (c1 - c2) + r2 741 r1 = <int>(dr + 0.5) 742 c0 = c1 = vc 743 r0 = vr 744 745 if (r0 - r1) * floor(r2 - r1) > 0: 746 if (r0 == r2) or (weight == 1): 747 t = (r0 - r1) / (r0 - 2. * r1 + r2) 748 else: 749 q = sqrt(4. * weight * weight * (r0 - r1) * (r2 - r1) + (r2 - r0) * floor(r2 - r0)) 750 if r1 < r0: 751 q = -q 752 t = (2. * weight * (r0 - r1) - r0 + r2 + q) / (2. * (1. - weight) * (r2 - r0)) 753 q = 1. / (2. * t * (1. - t) * (weight - 1.) + 1.) 754 dc = (t * t * (c0 - 2. * weight * c1 + c2) + 2. * t * (weight * c1 - c0) + c0) * q 755 dr = (t * t * (r0 - 2. * weight * r1 + r2) + 2. * t * (weight * r1 - r0) + r0) * q 756 ww = t * (weight - 1.) + 1. 757 ww *= ww * q 758 weight = ((1. - t) * (weight - 1.) + 1.) * sqrt(q) 759 vc = <int>(dc + 0.5) 760 vr = <int>(dr + 0.5) 761 dc = (c1 - c0) * (dr - r0) / (r1 - r0) + c0 762 763 rr_t, cc_t = _bezier_segment(r0, c0, vr, <int>(dc + 0.5), vr, vc, ww) 764 cc.extend(cc_t) 765 rr.extend(rr_t) 766 767 dc = (c1 - c2) * (dr - r2) / (r1 - r2) + c2 768 c1 = <int>(dc + 0.5) 769 c0 = vc 770 r0 = r1 = vr 771 772 rr_t, cc_t = _bezier_segment(r0, c0, r1, c1, r2, c2, weight * weight) 773 cc.extend(cc_t) 774 rr.extend(rr_t) 775 776 if shape is not None: 777 return _coords_inside_image(np.array(rr, dtype=np.intp), 778 np.array(cc, dtype=np.intp), shape) 779 return np.array(rr, dtype=np.intp), np.array(cc, dtype=np.intp) 780