1""" 2.. module:: helpers 3 :platform: Unix, Windows 4 :synopsis: Provides spline evaluation helper functions 5 6.. moduleauthor:: Onur Rauf Bingol <orbingol@gmail.com> 7 8""" 9 10import os 11from copy import deepcopy 12from . import linalg 13from .exceptions import GeomdlException 14try: 15 from functools import lru_cache 16except ImportError: 17 from .functools_lru_cache import lru_cache 18 19 20def find_span_binsearch(degree, knot_vector, num_ctrlpts, knot, **kwargs): 21 """ Finds the span of the knot over the input knot vector using binary search. 22 23 Implementation of Algorithm A2.1 from The NURBS Book by Piegl & Tiller. 24 25 The NURBS Book states that the knot span index always starts from zero, i.e. for a knot vector [0, 0, 1, 1]; 26 if FindSpan returns 1, then the knot is between the half-open interval [0, 1). 27 28 :param degree: degree, :math:`p` 29 :type degree: int 30 :param knot_vector: knot vector, :math:`U` 31 :type knot_vector: list, tuple 32 :param num_ctrlpts: number of control points, :math:`n + 1` 33 :type num_ctrlpts: int 34 :param knot: knot or parameter, :math:`u` 35 :type knot: float 36 :return: knot span 37 :rtype: int 38 """ 39 # Get tolerance value 40 tol = kwargs.get('tol', 10e-6) 41 42 # In The NURBS Book; number of knots = m + 1, number of control points = n + 1, p = degree 43 # All knot vectors should follow the rule: m = p + n + 1 44 n = num_ctrlpts - 1 45 if abs(knot_vector[n + 1] - knot) <= tol: 46 return n 47 48 # Set max and min positions of the array to be searched 49 low = degree 50 high = num_ctrlpts 51 52 # The division could return a float value which makes it impossible to use as an array index 53 mid = (low + high) / 2 54 # Direct int casting would cause numerical errors due to discarding the significand figures (digits after the dot) 55 # The round function could return unexpected results, so we add the floating point with some small number 56 # This addition would solve the issues caused by the division operation and how Python stores float numbers. 57 # E.g. round(13/2) = 6 (expected to see 7) 58 mid = int(round(mid + tol)) 59 60 # Search for the span 61 while (knot < knot_vector[mid]) or (knot >= knot_vector[mid + 1]): 62 if knot < knot_vector[mid]: 63 high = mid 64 else: 65 low = mid 66 mid = int((low + high) / 2) 67 68 return mid 69 70 71def find_span_linear(degree, knot_vector, num_ctrlpts, knot, **kwargs): 72 """ Finds the span of a single knot over the knot vector using linear search. 73 74 Alternative implementation for the Algorithm A2.1 from The NURBS Book by Piegl & Tiller. 75 76 :param degree: degree, :math:`p` 77 :type degree: int 78 :param knot_vector: knot vector, :math:`U` 79 :type knot_vector: list, tuple 80 :param num_ctrlpts: number of control points, :math:`n + 1` 81 :type num_ctrlpts: int 82 :param knot: knot or parameter, :math:`u` 83 :type knot: float 84 :return: knot span 85 :rtype: int 86 """ 87 span = 0 # Knot span index starts from zero 88 while span < num_ctrlpts and knot_vector[span] <= knot: 89 span += 1 90 91 return span - 1 92 93 94def find_spans(degree, knot_vector, num_ctrlpts, knots, func=find_span_linear): 95 """ Finds spans of a list of knots over the knot vector. 96 97 :param degree: degree, :math:`p` 98 :type degree: int 99 :param knot_vector: knot vector, :math:`U` 100 :type knot_vector: list, tuple 101 :param num_ctrlpts: number of control points, :math:`n + 1` 102 :type num_ctrlpts: int 103 :param knots: list of knots or parameters 104 :type knots: list, tuple 105 :param func: function for span finding, e.g. linear or binary search 106 :return: list of spans 107 :rtype: list 108 """ 109 spans = [] 110 for knot in knots: 111 spans.append(func(degree, knot_vector, num_ctrlpts, knot)) 112 return spans 113 114 115def find_multiplicity(knot, knot_vector, **kwargs): 116 """ Finds knot multiplicity over the knot vector. 117 118 Keyword Arguments: 119 * ``tol``: tolerance (delta) value for equality checking 120 121 :param knot: knot or parameter, :math:`u` 122 :type knot: float 123 :param knot_vector: knot vector, :math:`U` 124 :type knot_vector: list, tuple 125 :return: knot multiplicity, :math:`s` 126 :rtype: int 127 """ 128 # Get tolerance value 129 tol = kwargs.get('tol', 10e-8) 130 131 mult = 0 # initial multiplicity 132 133 for kv in knot_vector: 134 if abs(knot - kv) <= tol: 135 mult += 1 136 137 return mult 138 139 140def basis_function(degree, knot_vector, span, knot): 141 """ Computes the non-vanishing basis functions for a single parameter. 142 143 Implementation of Algorithm A2.2 from The NURBS Book by Piegl & Tiller. 144 Uses recurrence to compute the basis functions, also known as Cox - de 145 Boor recursion formula. 146 147 :param degree: degree, :math:`p` 148 :type degree: int 149 :param knot_vector: knot vector, :math:`U` 150 :type knot_vector: list, tuple 151 :param span: knot span, :math:`i` 152 :type span: int 153 :param knot: knot or parameter, :math:`u` 154 :type knot: float 155 :return: basis functions 156 :rtype: list 157 """ 158 left = [0.0 for _ in range(degree + 1)] 159 right = [0.0 for _ in range(degree + 1)] 160 N = [1.0 for _ in range(degree + 1)] # N[0] = 1.0 by definition 161 162 for j in range(1, degree + 1): 163 left[j] = knot - knot_vector[span + 1 - j] 164 right[j] = knot_vector[span + j] - knot 165 saved = 0.0 166 for r in range(0, j): 167 temp = N[r] / (right[r + 1] + left[j - r]) 168 N[r] = saved + right[r + 1] * temp 169 saved = left[j - r] * temp 170 N[j] = saved 171 172 return N 173 174 175def basis_function_one(degree, knot_vector, span, knot): 176 """ Computes the value of a basis function for a single parameter. 177 178 Implementation of Algorithm 2.4 from The NURBS Book by Piegl & Tiller. 179 180 :param degree: degree, :math:`p` 181 :type degree: int 182 :param knot_vector: knot vector 183 :type knot_vector: list, tuple 184 :param span: knot span, :math:`i` 185 :type span: int 186 :param knot: knot or parameter, :math:`u` 187 :type knot: float 188 :return: basis function, :math:`N_{i,p}` 189 :rtype: float 190 """ 191 # Special case at boundaries 192 if (span == 0 and knot == knot_vector[0]) or \ 193 (span == len(knot_vector) - degree - 2) and knot == knot_vector[len(knot_vector) - 1]: 194 return 1.0 195 196 # Knot is outside of span range 197 if knot < knot_vector[span] or knot >= knot_vector[span + degree + 1]: 198 return 0.0 199 200 N = [0.0 for _ in range(degree + span + 1)] 201 202 # Initialize the zeroth degree basis functions 203 for j in range(0, degree + 1): 204 if knot_vector[span + j] <= knot < knot_vector[span + j + 1]: 205 N[j] = 1.0 206 207 # Computing triangular table of basis functions 208 for k in range(1, degree + 1): 209 # Detecting zeros saves computations 210 saved = 0.0 211 if N[0] != 0.0: 212 saved = ((knot - knot_vector[span]) * N[0]) / (knot_vector[span + k] - knot_vector[span]) 213 214 for j in range(0, degree - k + 1): 215 Uleft = knot_vector[span + j + 1] 216 Uright = knot_vector[span + j + k + 1] 217 218 # Zero detection 219 if N[j + 1] == 0.0: 220 N[j] = saved 221 saved = 0.0 222 else: 223 temp = N[j + 1] / (Uright - Uleft) 224 N[j] = saved + (Uright - knot) * temp 225 saved = (knot - Uleft) * temp 226 227 return N[0] 228 229 230def basis_functions(degree, knot_vector, spans, knots): 231 """ Computes the non-vanishing basis functions for a list of parameters. 232 233 Wrapper for :func:`.helpers.basis_function` to process multiple span 234 and knot values. Uses recurrence to compute the basis functions, also 235 known as Cox - de Boor recursion formula. 236 237 :param degree: degree, :math:`p` 238 :type degree: int 239 :param knot_vector: knot vector, :math:`U` 240 :type knot_vector: list, tuple 241 :param spans: list of knot spans 242 :type spans: list, tuple 243 :param knots: list of knots or parameters 244 :type knots: list, tuple 245 :return: basis functions 246 :rtype: list 247 """ 248 basis = [] 249 for span, knot in zip(spans, knots): 250 basis.append(basis_function(degree, knot_vector, span, knot)) 251 return basis 252 253 254def basis_function_all(degree, knot_vector, span, knot): 255 """ Computes all non-zero basis functions of all degrees from 0 up to the input degree for a single parameter. 256 257 A slightly modified version of Algorithm A2.2 from The NURBS Book by Piegl & Tiller. 258 Wrapper for :func:`.helpers.basis_function` to compute multiple basis functions. 259 Uses recurrence to compute the basis functions, also known as Cox - de Boor 260 recursion formula. 261 262 For instance; if ``degree = 2``, then this function will compute the basis function 263 values of degrees **0, 1** and **2** for the ``knot`` value at the input knot ``span`` 264 of the ``knot_vector``. 265 266 :param degree: degree, :math:`p` 267 :type degree: int 268 :param knot_vector: knot vector, :math:`U` 269 :type knot_vector: list, tuple 270 :param span: knot span, :math:`i` 271 :type span: int 272 :param knot: knot or parameter, :math:`u` 273 :type knot: float 274 :return: basis functions 275 :rtype: list 276 """ 277 N = [[None for _ in range(degree + 1)] for _ in range(degree + 1)] 278 for i in range(0, degree + 1): 279 bfuns = basis_function(i, knot_vector, span, knot) 280 for j in range(0, i + 1): 281 N[j][i] = bfuns[j] 282 return N 283 284 285def basis_function_ders(degree, knot_vector, span, knot, order): 286 """ Computes derivatives of the basis functions for a single parameter. 287 288 Implementation of Algorithm A2.3 from The NURBS Book by Piegl & Tiller. 289 290 :param degree: degree, :math:`p` 291 :type degree: int 292 :param knot_vector: knot vector, :math:`U` 293 :type knot_vector: list, tuple 294 :param span: knot span, :math:`i` 295 :type span: int 296 :param knot: knot or parameter, :math:`u` 297 :type knot: float 298 :param order: order of the derivative 299 :type order: int 300 :return: derivatives of the basis functions 301 :rtype: list 302 """ 303 # Initialize variables 304 left = [1.0 for _ in range(degree + 1)] 305 right = [1.0 for _ in range(degree + 1)] 306 ndu = [[1.0 for _ in range(degree + 1)] for _ in range(degree + 1)] # N[0][0] = 1.0 by definition 307 308 for j in range(1, degree + 1): 309 left[j] = knot - knot_vector[span + 1 - j] 310 right[j] = knot_vector[span + j] - knot 311 saved = 0.0 312 r = 0 313 for r in range(r, j): 314 # Lower triangle 315 ndu[j][r] = right[r + 1] + left[j - r] 316 temp = ndu[r][j - 1] / ndu[j][r] 317 # Upper triangle 318 ndu[r][j] = saved + (right[r + 1] * temp) 319 saved = left[j - r] * temp 320 ndu[j][j] = saved 321 322 # Load the basis functions 323 ders = [[0.0 for _ in range(degree + 1)] for _ in range((min(degree, order) + 1))] 324 for j in range(0, degree + 1): 325 ders[0][j] = ndu[j][degree] 326 327 # Start calculating derivatives 328 a = [[1.0 for _ in range(degree + 1)] for _ in range(2)] 329 # Loop over function index 330 for r in range(0, degree + 1): 331 # Alternate rows in array a 332 s1 = 0 333 s2 = 1 334 a[0][0] = 1.0 335 # Loop to compute k-th derivative 336 for k in range(1, order + 1): 337 d = 0.0 338 rk = r - k 339 pk = degree - k 340 if r >= k: 341 a[s2][0] = a[s1][0] / ndu[pk + 1][rk] 342 d = a[s2][0] * ndu[rk][pk] 343 if rk >= -1: 344 j1 = 1 345 else: 346 j1 = -rk 347 if (r - 1) <= pk: 348 j2 = k - 1 349 else: 350 j2 = degree - r 351 for j in range(j1, j2 + 1): 352 a[s2][j] = (a[s1][j] - a[s1][j - 1]) / ndu[pk + 1][rk + j] 353 d += (a[s2][j] * ndu[rk + j][pk]) 354 if r <= pk: 355 a[s2][k] = -a[s1][k - 1] / ndu[pk + 1][r] 356 d += (a[s2][k] * ndu[r][pk]) 357 ders[k][r] = d 358 359 # Switch rows 360 j = s1 361 s1 = s2 362 s2 = j 363 364 # Multiply through by the the correct factors 365 r = float(degree) 366 for k in range(1, order + 1): 367 for j in range(0, degree + 1): 368 ders[k][j] *= r 369 r *= (degree - k) 370 371 # Return the basis function derivatives list 372 return ders 373 374 375def basis_function_ders_one(degree, knot_vector, span, knot, order): 376 """ Computes the derivative of one basis functions for a single parameter. 377 378 Implementation of Algorithm A2.5 from The NURBS Book by Piegl & Tiller. 379 380 :param degree: degree, :math:`p` 381 :type degree: int 382 :param knot_vector: knot_vector, :math:`U` 383 :type knot_vector: list, tuple 384 :param span: knot span, :math:`i` 385 :type span: int 386 :param knot: knot or parameter, :math:`u` 387 :type knot: float 388 :param order: order of the derivative 389 :type order: int 390 :return: basis function derivatives 391 :rtype: list 392 """ 393 ders = [0.0 for _ in range(0, order + 1)] 394 395 # Knot is outside of span range 396 if (knot < knot_vector[span]) or (knot >= knot_vector[span + degree + 1]): 397 for k in range(0, order + 1): 398 ders[k] = 0.0 399 400 return ders 401 402 N = [[0.0 for _ in range(0, degree + 1)] for _ in range(0, degree + 1)] 403 404 # Initializing the zeroth degree basis functions 405 for j in range(0, degree + 1): 406 if knot_vector[span + j] <= knot < knot_vector[span + j + 1]: 407 N[j][0] = 1.0 408 409 # Computing all basis functions values for all degrees inside the span 410 for k in range(1, degree + 1): 411 saved = 0.0 412 # Detecting zeros saves computations 413 if N[0][k - 1] != 0.0: 414 saved = ((knot - knot_vector[span]) * N[0][k - 1]) / (knot_vector[span + k] - knot_vector[span]) 415 416 for j in range(0, degree - k + 1): 417 Uleft = knot_vector[span + j + 1] 418 Uright = knot_vector[span + j + k + 1] 419 420 # Zero detection 421 if N[j + 1][k - 1] == 0.0: 422 N[j][k] = saved 423 saved = 0.0 424 else: 425 temp = N[j + 1][k - 1] / (Uright - Uleft) 426 N[j][k] = saved + (Uright - knot) * temp 427 saved = (knot - Uleft) * temp 428 429 # The basis function value is the zeroth derivative 430 ders[0] = N[0][degree] 431 432 # Computing the basis functions derivatives 433 for k in range(1, order + 1): 434 # Buffer for computing the kth derivative 435 ND = [0.0 for _ in range(0, k + 1)] 436 437 # Basis functions values used for the derivative 438 for j in range(0, k + 1): 439 ND[j] = N[j][degree - k] 440 441 # Computing derivatives used for the kth basis function derivative 442 443 # Derivative order for the k-th basis function derivative 444 for jj in range(1, k + 1): 445 if ND[0] == 0.0: 446 saved = 0.0 447 else: 448 saved = ND[0] / (knot_vector[span + degree - k + jj] - knot_vector[span]) 449 450 # Index of the Basis function derivatives 451 for j in range(0, k - jj + 1): 452 Uleft = knot_vector[span + j + 1] 453 # Wrong in The NURBS Book: -k is missing. 454 # The right expression is the same as for saved with the added j offset 455 Uright = knot_vector[span + j + degree - k + jj + 1] 456 457 if ND[j + 1] == 0.0: 458 ND[j] = (degree - k + jj) * saved 459 saved = 0.0 460 else: 461 temp = ND[j + 1] / (Uright - Uleft) 462 463 ND[j] = (degree - k + jj) * (saved - temp) 464 saved = temp 465 466 ders[k] = ND[0] 467 468 return ders 469 470 471def basis_functions_ders(degree, knot_vector, spans, knots, order): 472 """ Computes derivatives of the basis functions for a list of parameters. 473 474 Wrapper for :func:`.helpers.basis_function_ders` to process multiple span 475 and knot values. 476 477 :param degree: degree, :math:`p` 478 :type degree: int 479 :param knot_vector: knot vector, :math:`U` 480 :type knot_vector: list, tuple 481 :param spans: list of knot spans 482 :type spans: list, tuple 483 :param knots: list of knots or parameters 484 :type knots: list, tuple 485 :param order: order of the derivative 486 :type order: int 487 :return: derivatives of the basis functions 488 :rtype: list 489 """ 490 basis_ders = [] 491 for span, knot in zip(spans, knots): 492 basis_ders.append(basis_function_ders(degree, knot_vector, span, knot, order)) 493 return basis_ders 494 495 496def knot_insertion(degree, knotvector, ctrlpts, u, **kwargs): 497 """ Computes the control points of the rational/non-rational spline after knot insertion. 498 499 Part of Algorithm A5.1 of The NURBS Book by Piegl & Tiller, 2nd Edition. 500 501 Keyword Arguments: 502 * ``num``: number of knot insertions. *Default: 1* 503 * ``s``: multiplicity of the knot. *Default: computed via :func:`.find_multiplicity`* 504 * ``span``: knot span. *Default: computed via :func:`.find_span_linear`* 505 506 :param degree: degree 507 :type degree: int 508 :param knotvector: knot vector 509 :type knotvector: list, tuple 510 :param ctrlpts: control points 511 :type ctrlpts: list 512 :param u: knot to be inserted 513 :type u: float 514 :return: updated control points 515 :rtype: list 516 """ 517 # Get keyword arguments 518 num = kwargs.get('num', 1) # number of knot insertions 519 s = kwargs.get('s', find_multiplicity(u, knotvector)) # multiplicity 520 k = kwargs.get('span', find_span_linear(degree, knotvector, len(ctrlpts), u)) # knot span 521 522 # Initialize variables 523 np = len(ctrlpts) 524 nq = np + num 525 526 # Initialize new control points array (control points may be weighted or not) 527 ctrlpts_new = [[] for _ in range(nq)] 528 529 # Initialize a local array of length p + 1 530 temp = [[] for _ in range(degree + 1)] 531 532 # Save unaltered control points 533 for i in range(0, k - degree + 1): 534 ctrlpts_new[i] = ctrlpts[i] 535 for i in range(k - s, np): 536 ctrlpts_new[i + num] = ctrlpts[i] 537 538 # Start filling the temporary local array which will be used to update control points during knot insertion 539 for i in range(0, degree - s + 1): 540 temp[i] = deepcopy(ctrlpts[k - degree + i]) 541 542 # Insert knot "num" times 543 for j in range(1, num + 1): 544 L = k - degree + j 545 for i in range(0, degree - j - s + 1): 546 alpha = knot_insertion_alpha(u, tuple(knotvector), k, i, L) 547 if isinstance(temp[i][0], float): 548 temp[i][:] = [alpha * elem2 + (1.0 - alpha) * elem1 for elem1, elem2 in zip(temp[i], temp[i + 1])] 549 else: 550 for idx in range(len(temp[i])): 551 temp[i][idx][:] = [alpha * elem2 + (1.0 - alpha) * elem1 for elem1, elem2 in 552 zip(temp[i][idx], temp[i + 1][idx])] 553 ctrlpts_new[L] = deepcopy(temp[0]) 554 ctrlpts_new[k + num - j - s] = deepcopy(temp[degree - j - s]) 555 556 # Load remaining control points 557 L = k - degree + num 558 for i in range(L + 1, k - s): 559 ctrlpts_new[i] = deepcopy(temp[i - L]) 560 561 # Return control points after knot insertion 562 return ctrlpts_new 563 564 565@lru_cache(maxsize=os.environ['GEOMDL_CACHE_SIZE'] if "GEOMDL_CACHE_SIZE" in os.environ else 128) 566def knot_insertion_alpha(u, knotvector, span, idx, leg): 567 """ Computes :math:`\\alpha` coefficient for knot insertion algorithm. 568 569 :param u: knot 570 :type u: float 571 :param knotvector: knot vector 572 :type knotvector: tuple 573 :param span: knot span 574 :type span: int 575 :param idx: index value (degree-dependent) 576 :type idx: int 577 :param leg: i-th leg of the control points polygon 578 :type leg: int 579 :return: coefficient value 580 :rtype: float 581 """ 582 return (u - knotvector[leg + idx]) / (knotvector[idx + span + 1] - knotvector[leg + idx]) 583 584 585def knot_insertion_kv(knotvector, u, span, r): 586 """ Computes the knot vector of the rational/non-rational spline after knot insertion. 587 588 Part of Algorithm A5.1 of The NURBS Book by Piegl & Tiller, 2nd Edition. 589 590 :param knotvector: knot vector 591 :type knotvector: list, tuple 592 :param u: knot 593 :type u: float 594 :param span: knot span 595 :type span: int 596 :param r: number of knot insertions 597 :type r: int 598 :return: updated knot vector 599 :rtype: list 600 """ 601 # Initialize variables 602 kv_size = len(knotvector) 603 kv_updated = [0.0 for _ in range(kv_size + r)] 604 605 # Compute new knot vector 606 for i in range(0, span + 1): 607 kv_updated[i] = knotvector[i] 608 for i in range(1, r + 1): 609 kv_updated[span + i] = u 610 for i in range(span + 1, kv_size): 611 kv_updated[i + r] = knotvector[i] 612 613 # Return the new knot vector 614 return kv_updated 615 616 617def knot_removal(degree, knotvector, ctrlpts, u, **kwargs): 618 """ Computes the control points of the rational/non-rational spline after knot removal. 619 620 Implementation based on Algorithm A5.8 and Equation 5.28 of The NURBS Book by Piegl & Tiller 621 622 Keyword Arguments: 623 * ``num``: number of knot removals 624 625 :param degree: degree 626 :type degree: int 627 :param knotvector: knot vector 628 :type knotvector: list, tuple 629 :param ctrlpts: control points 630 :type ctrlpts: list 631 :param u: knot to be removed 632 :type u: float 633 :return: updated control points 634 :rtype: list 635 """ 636 tol = kwargs.get('tol', 10e-4) # Refer to Eq 5.30 for the meaning 637 num = kwargs.get('num', 1) # number of same knot removals 638 s = kwargs.get('s', find_multiplicity(u, knotvector)) # multiplicity 639 r = kwargs.get('span', find_span_linear(degree, knotvector, len(ctrlpts), u)) # knot span 640 641 # Edge case 642 if num < 1: 643 return ctrlpts 644 645 # Initialize variables 646 first = r - degree 647 last = r - s 648 649 # Don't change input variables, prepare new ones for updating 650 ctrlpts_new = deepcopy(ctrlpts) 651 652 is_volume = True 653 if isinstance(ctrlpts_new[0][0], float): 654 is_volume = False 655 656 # Initialize temp array for storing new control points 657 if is_volume: 658 temp = [[[] for _ in range(len(ctrlpts_new[0]))] for _ in range((2 * degree) + 1)] 659 else: 660 temp = [[] for _ in range((2 * degree) + 1)] 661 662 # Loop for Eqs 5.28 & 5.29 663 for t in range(0, num): 664 temp[0] = ctrlpts[first - 1] 665 temp[last - first + 2] = ctrlpts[last + 1] 666 i = first 667 j = last 668 ii = 1 669 jj = last - first + 1 670 remflag = False 671 672 # Compute control points for one removal step 673 while j - i >= t: 674 alpha_i = knot_removal_alpha_i(u, degree, tuple(knotvector), t, i) 675 alpha_j = knot_removal_alpha_j(u, degree, tuple(knotvector), t, j) 676 if is_volume: 677 for idx in range(len(ctrlpts[0])): 678 temp[ii][idx] = [(cpt - (1.0 - alpha_i) * ti) / alpha_i for cpt, ti 679 in zip(ctrlpts[i][idx], temp[ii - 1][idx])] 680 temp[jj][idx] = [(cpt - alpha_j * tj) / (1.0 - alpha_j) for cpt, tj 681 in zip(ctrlpts[j][idx], temp[jj + 1][idx])] 682 else: 683 temp[ii] = [(cpt - (1.0 - alpha_i) * ti) / alpha_i for cpt, ti in zip(ctrlpts[i], temp[ii - 1])] 684 temp[jj] = [(cpt - alpha_j * tj) / (1.0 - alpha_j) for cpt, tj in zip(ctrlpts[j], temp[jj + 1])] 685 i += 1 686 j -= 1 687 ii += 1 688 jj -= 1 689 690 # Check if the knot is removable 691 if j - i < t: 692 if is_volume: 693 if linalg.point_distance(temp[ii - 1][0], temp[jj + 1][0]) <= tol: 694 remflag = True 695 else: 696 if linalg.point_distance(temp[ii - 1], temp[jj + 1]) <= tol: 697 remflag = True 698 else: 699 alpha_i = knot_removal_alpha_i(u, degree, tuple(knotvector), t, i) 700 if is_volume: 701 ptn = [(alpha_i * t1) + ((1.0 - alpha_i) * t2) for t1, t2 in zip(temp[ii + t + 1][0], temp[ii - 1][0])] 702 else: 703 ptn = [(alpha_i * t1) + ((1.0 - alpha_i) * t2) for t1, t2 in zip(temp[ii + t + 1], temp[ii - 1])] 704 if linalg.point_distance(ctrlpts[i], ptn) <= tol: 705 remflag = True 706 707 # Check if we can remove the knot and update new control points array 708 if remflag: 709 i = first 710 j = last 711 while j - i > t: 712 ctrlpts_new[i] = temp[i - first + 1] 713 ctrlpts_new[j] = temp[j - first + 1] 714 i += 1 715 j -= 1 716 717 # Update indices 718 first -= 1 719 last += 1 720 721 # Fix indexing 722 t += 1 723 724 # Shift control points (refer to p.183 of The NURBS Book, 2nd Edition) 725 j = int((2*r - s - degree) / 2) # first control point out 726 i = j 727 for k in range(1, t): 728 if k % 2 == 1: 729 i += 1 730 else: 731 j -= 1 732 for k in range(i+1, len(ctrlpts)): 733 ctrlpts_new[j] = ctrlpts[k] 734 j += 1 735 736 # Slice to get the new control points 737 ctrlpts_new = ctrlpts_new[0:-t] 738 739 return ctrlpts_new 740 741 742@lru_cache(maxsize=os.environ['GEOMDL_CACHE_SIZE'] if "GEOMDL_CACHE_SIZE" in os.environ else 128) 743def knot_removal_alpha_i(u, degree, knotvector, num, idx): 744 """ Computes :math:`\\alpha_{i}` coefficient for knot removal algorithm. 745 746 Please refer to Eq. 5.29 of The NURBS Book by Piegl & Tiller, 2nd Edition, p.184 for details. 747 748 :param u: knot 749 :type u: float 750 :param degree: degree 751 :type degree: int 752 :param knotvector: knot vector 753 :type knotvector: tuple 754 :param num: knot removal index 755 :type num: int 756 :param idx: iterator index 757 :type idx: int 758 :return: coefficient value 759 :rtype: float 760 """ 761 return (u - knotvector[idx]) / (knotvector[idx + degree + 1 + num] - knotvector[idx]) 762 763 764@lru_cache(maxsize=os.environ['GEOMDL_CACHE_SIZE'] if "GEOMDL_CACHE_SIZE" in os.environ else 128) 765def knot_removal_alpha_j(u, degree, knotvector, num, idx): 766 """ Computes :math:`\\alpha_{j}` coefficient for knot removal algorithm. 767 768 Please refer to Eq. 5.29 of The NURBS Book by Piegl & Tiller, 2nd Edition, p.184 for details. 769 770 :param u: knot 771 :type u: float 772 :param degree: degree 773 :type degree: int 774 :param knotvector: knot vector 775 :type knotvector: tuple 776 :param num: knot removal index 777 :type num: int 778 :param idx: iterator index 779 :type idx: int 780 :return: coefficient value 781 :rtype: float 782 """ 783 return (u - knotvector[idx - num]) / (knotvector[idx + degree + 1] - knotvector[idx - num]) 784 785 786def knot_removal_kv(knotvector, span, r): 787 """ Computes the knot vector of the rational/non-rational spline after knot removal. 788 789 Part of Algorithm A5.8 of The NURBS Book by Piegl & Tiller, 2nd Edition. 790 791 :param knotvector: knot vector 792 :type knotvector: list, tuple 793 :param span: knot span 794 :type span: int 795 :param r: number of knot removals 796 :type r: int 797 :return: updated knot vector 798 :rtype: list 799 """ 800 # Edge case 801 if r < 1: 802 return knotvector 803 804 # Create a deep copy of the input knot vector 805 kv_updated = deepcopy(knotvector) 806 807 # Shift knots 808 for k in range(span + 1, len(knotvector)): 809 kv_updated[k - r] = knotvector[k] 810 811 # Slice to get the new knot vector 812 kv_updated = kv_updated[0:-r] 813 814 # Return the new knot vector 815 return kv_updated 816 817 818def knot_refinement(degree, knotvector, ctrlpts, **kwargs): 819 """ Computes the knot vector and the control points of the rational/non-rational spline after knot refinement. 820 821 Implementation of Algorithm A5.4 of The NURBS Book by Piegl & Tiller, 2nd Edition. 822 823 The algorithm automatically find the knots to be refined, i.e. the middle knots in the knot vector, and their 824 multiplicities, i.e. number of same knots in the knot vector. This is the basis of knot refinement algorithm. 825 This operation can be overridden by providing a list of knots via ``knot_list`` argument. In addition, users can 826 provide a list of additional knots to be inserted in the knot vector via ``add_knot_list`` argument. 827 828 Moreover, a numerical ``density`` argument can be used to automate extra knot insertions. If ``density`` is bigger 829 than 1, then the algorithm finds the middle knots in each internal knot span to increase the number of knots to be 830 refined. 831 832 **Example**: Let the degree is 2 and the knot vector to be refined is ``[0, 2, 4]`` with the superfluous knots 833 from the start and end are removed. Knot vectors with the changing ``density (d)`` value will be: 834 835 * ``d = 1``, knot vector ``[0, 1, 1, 2, 2, 3, 3, 4]`` 836 * ``d = 2``, knot vector ``[0, 0.5, 0.5, 1, 1, 1.5, 1.5, 2, 2, 2.5, 2.5, 3, 3, 3.5, 3.5, 4]`` 837 838 Keyword Arguments: 839 * ``knot_list``: knot list to be refined. *Default: list of internal knots* 840 * ``add_knot_list``: additional list of knots to be refined. *Default: []* 841 * ``density``: Density of the knots. *Default: 1* 842 843 :param degree: degree 844 :type degree: int 845 :param knotvector: knot vector 846 :type knotvector: list, tuple 847 :param ctrlpts: control points 848 :return: updated control points and knot vector 849 :rtype: tuple 850 """ 851 # Get keyword arguments 852 tol = kwargs.get('tol', 10e-8) # tolerance value for zero equality checking 853 check_num = kwargs.get('check_num', True) # enables/disables input validity checking 854 knot_list = kwargs.get('knot_list', knotvector[degree:-degree]) 855 add_knot_list = kwargs.get('add_knot_list', list()) 856 density = kwargs.get('density', 1) 857 858 # Input validity checking 859 if check_num: 860 if not isinstance(density, int): 861 raise GeomdlException("Density value must be an integer", data=dict(density=density)) 862 863 if density < 1: 864 raise GeomdlException("Density value cannot be less than 1", data=dict(density=density)) 865 866 # Add additional knots to be refined 867 if add_knot_list: 868 knot_list += list(add_knot_list) 869 870 # Sort the list and convert to a set to make sure that the values are unique 871 knot_list = sorted(set(knot_list)) 872 873 # Increase knot density 874 for d in range(0, density): 875 rknots = [] 876 for i in range(len(knot_list) - 1): 877 knot_tmp = knot_list[i] + ((knot_list[i + 1] - knot_list[i]) / 2.0) 878 rknots.append(knot_list[i]) 879 rknots.append(knot_tmp) 880 rknots.append(knot_list[i + 1]) 881 knot_list = rknots 882 883 # Find how many knot insertions are necessary 884 X = [] 885 for mk in knot_list: 886 s = find_multiplicity(mk, knotvector) 887 r = degree - s 888 X += [mk for _ in range(r)] 889 890 # Check if the knot refinement is possible 891 if not X: 892 raise GeomdlException("Cannot refine knot vector on this parametric dimension") 893 894 # Initialize common variables 895 r = len(X) - 1 896 n = len(ctrlpts) - 1 897 m = n + degree + 1 898 a = find_span_linear(degree, knotvector, n, X[0]) 899 b = find_span_linear(degree, knotvector, n, X[r]) + 1 900 901 # Initialize new control points array 902 if isinstance(ctrlpts[0][0], float): 903 new_ctrlpts = [[] for _ in range(n + r + 2)] 904 else: 905 new_ctrlpts = [[[] for _ in range(len(ctrlpts[0]))] for _ in range(n + r + 2)] 906 907 # Fill unchanged control points 908 for j in range(0, a - degree + 1): 909 new_ctrlpts[j] = ctrlpts[j] 910 for j in range(b - 1, n + 1): 911 new_ctrlpts[j + r + 1] = ctrlpts[j] 912 913 # Initialize new knot vector array 914 new_kv = [0.0 for _ in range(m + r + 2)] 915 916 # Fill unchanged knots 917 for j in range(0, a + 1): 918 new_kv[j] = knotvector[j] 919 for j in range(b + degree, m + 1): 920 new_kv[j + r + 1] = knotvector[j] 921 922 # Initialize variables for knot refinement 923 i = b + degree - 1 924 k = b + degree + r 925 j = r 926 927 # Apply knot refinement 928 while j >= 0: 929 while X[j] <= knotvector[i] and i > a: 930 new_ctrlpts[k - degree - 1] = ctrlpts[i - degree - 1] 931 new_kv[k] = knotvector[i] 932 k -= 1 933 i -= 1 934 new_ctrlpts[k - degree - 1] = deepcopy(new_ctrlpts[k - degree]) 935 for l in range(1, degree + 1): 936 idx = k - degree + l 937 alpha = new_kv[k + l] - X[j] 938 if abs(alpha) < tol: 939 new_ctrlpts[idx - 1] = deepcopy(new_ctrlpts[idx]) 940 else: 941 alpha = alpha / (new_kv[k + l] - knotvector[i - degree + l]) 942 if isinstance(ctrlpts[0][0], float): 943 new_ctrlpts[idx - 1] = [alpha * p1 + (1.0 - alpha) * p2 for p1, p2 in 944 zip(new_ctrlpts[idx - 1], new_ctrlpts[idx])] 945 else: 946 for idx2 in range(len(ctrlpts[0])): 947 new_ctrlpts[idx - 1][idx2] = [alpha * p1 + (1.0 - alpha) * p2 for p1, p2 in 948 zip(new_ctrlpts[idx - 1][idx2], new_ctrlpts[idx][idx2])] 949 new_kv[k] = X[j] 950 k = k - 1 951 j -= 1 952 953 # Return control points and knot vector after refinement 954 return new_ctrlpts, new_kv 955 956 957def degree_elevation(degree, ctrlpts, **kwargs): 958 """ Computes the control points of the rational/non-rational spline after degree elevation. 959 960 Implementation of Eq. 5.36 of The NURBS Book by Piegl & Tiller, 2nd Edition, p.205 961 962 Keyword Arguments: 963 * ``num``: number of degree elevations 964 965 Please note that degree elevation algorithm can only operate on Bezier shapes, i.e. curves, surfaces, volumes. 966 967 :param degree: degree 968 :type degree: int 969 :param ctrlpts: control points 970 :type ctrlpts: list, tuple 971 :return: control points of the degree-elevated shape 972 :rtype: list 973 """ 974 # Get keyword arguments 975 num = kwargs.get('num', 1) # number of degree elevations 976 check_op = kwargs.get('check_num', True) # enable/disable input validation checks 977 978 if check_op: 979 if degree + 1 != len(ctrlpts): 980 raise GeomdlException("Degree elevation can only work with Bezier-type geometries") 981 if num <= 0: 982 raise GeomdlException("Cannot degree elevate " + str(num) + " times") 983 984 # Initialize variables 985 num_pts_elev = degree + 1 + num 986 pts_elev = [[0.0 for _ in range(len(ctrlpts[0]))] for _ in range(num_pts_elev)] 987 988 # Compute control points of degree-elevated 1-dimensional shape 989 for i in range(0, num_pts_elev): 990 start = max(0, (i - num)) 991 end = min(degree, i) 992 for j in range(start, end + 1): 993 coeff = linalg.binomial_coefficient(degree, j) * linalg.binomial_coefficient(num, (i - j)) 994 coeff /= linalg.binomial_coefficient((degree + num), i) 995 pts_elev[i] = [p1 + (coeff * p2) for p1, p2 in zip(pts_elev[i], ctrlpts[j])] 996 997 # Return computed control points after degree elevation 998 return pts_elev 999 1000 1001def degree_reduction(degree, ctrlpts, **kwargs): 1002 """ Computes the control points of the rational/non-rational spline after degree reduction. 1003 1004 Implementation of Eqs. 5.41 and 5.42 of The NURBS Book by Piegl & Tiller, 2nd Edition, p.220 1005 1006 Please note that degree reduction algorithm can only operate on Bezier shapes, i.e. curves, surfaces, volumes and 1007 this implementation does NOT compute the maximum error tolerance as described via Eqs. 5.45 and 5.46 of The NURBS 1008 Book by Piegl & Tiller, 2nd Edition, p.221 to determine whether the shape is degree reducible or not. 1009 1010 :param degree: degree 1011 :type degree: int 1012 :param ctrlpts: control points 1013 :type ctrlpts: list, tuple 1014 :return: control points of the degree-reduced shape 1015 :rtype: list 1016 """ 1017 # Get keyword arguments 1018 check_op = kwargs.get('check_num', True) # enable/disable input validation checks 1019 1020 if check_op: 1021 if degree + 1 != len(ctrlpts): 1022 raise GeomdlException("Degree reduction can only work with Bezier-type geometries") 1023 if degree < 2: 1024 raise GeomdlException("Input spline geometry must have degree > 1") 1025 1026 # Initialize variables 1027 pts_red = [[0.0 for _ in range(len(ctrlpts[0]))] for _ in range(degree)] 1028 1029 # Fix start and end control points 1030 pts_red[0] = ctrlpts[0] 1031 pts_red[-1] = ctrlpts[-1] 1032 1033 # Find if the degree is an even or an odd number 1034 p_is_odd = True if degree % 2 != 0 else False 1035 1036 # Compute control points of degree-reduced 1-dimensional shape 1037 r = int((degree - 1) / 2) 1038 # Handle a special case when degree = 2 1039 if degree == 2: 1040 r1 = r - 2 1041 else: 1042 # Determine r1 w.r.t. degree evenness 1043 r1 = r - 1 if p_is_odd else r 1044 for i in range(1, r1 + 1): 1045 alpha = float(i) / float(degree) 1046 pts_red[i] = [(c1 - (alpha * c2)) / (1 - alpha) for c1, c2 in zip(ctrlpts[i], pts_red[i - 1])] 1047 for i in range(degree - 2, r1 + 2): 1048 alpha = float(i + 1) / float(degree) 1049 pts_red[i] = [(c1 - ((1 - alpha) * c2)) / alpha for c1, c2 in zip(ctrlpts[i + 1], pts_red[i + 1])] 1050 1051 if p_is_odd: 1052 alpha = float(r) / float(degree) 1053 left = [(c1 - (alpha * c2)) / (1 - alpha) for c1, c2 in zip(ctrlpts[r], pts_red[r - 1])] 1054 alpha = float(r + 1) / float(degree) 1055 right = [(c1 - ((1 - alpha) * c2)) / alpha for c1, c2 in zip(ctrlpts[r + 1], pts_red[r + 1])] 1056 pts_red[r] = [0.5 * (pl + pr) for pl, pr in zip(left, right)] 1057 1058 # Return computed control points after degree reduction 1059 return pts_red 1060 1061 1062def curve_deriv_cpts(dim, degree, kv, cpts, rs, deriv_order=0): 1063 """ Compute control points of curve derivatives. 1064 1065 Implementation of Algorithm A3.3 from The NURBS Book by Piegl & Tiller. 1066 1067 :param dim: spatial dimension of the curve 1068 :type dim: int 1069 :param degree: degree of the curve 1070 :type degree: int 1071 :param kv: knot vector 1072 :type kv: list, tuple 1073 :param cpts: control points 1074 :type cpts: list, tuple 1075 :param rs: minimum (r1) and maximum (r2) knot spans that the curve derivative will be computed 1076 :param deriv_order: derivative order, i.e. the i-th derivative 1077 :type deriv_order: int 1078 :return: control points of the derivative curve over the input knot span range 1079 :rtype: list 1080 """ 1081 r = rs[1] - rs[0] 1082 1083 # Initialize return value (control points) 1084 PK = [[[None for _ in range(dim)] for _ in range(r + 1)] for _ in range(deriv_order + 1)] 1085 1086 # Algorithm A3.3 1087 for i in range(0, r + 1): 1088 PK[0][i][:] = [elem for elem in cpts[rs[0] + i]] 1089 1090 for k in range(1, deriv_order + 1): 1091 tmp = degree - k + 1 1092 for i in range(0, r - k + 1): 1093 PK[k][i][:] = [tmp * (elem1 - elem2) / 1094 (kv[rs[0] + i + degree + 1] - kv[rs[0] + i + k]) for elem1, elem2 1095 in zip(PK[k - 1][i + 1], PK[k - 1][i])] 1096 1097 # Return control points (as a 2-dimensional list of points) 1098 return PK 1099 1100 1101def surface_deriv_cpts(dim, degree, kv, cpts, cpsize, rs, ss, deriv_order=0): 1102 """ Compute control points of surface derivatives. 1103 1104 Implementation of Algorithm A3.7 from The NURBS Book by Piegl & Tiller. 1105 1106 :param dim: spatial dimension of the surface 1107 :type dim: int 1108 :param degree: degrees 1109 :type degree: list, tuple 1110 :param kv: knot vectors 1111 :type kv: list, tuple 1112 :param cpts: control points 1113 :type cpts: list, tuple 1114 :param cpsize: number of control points in all parametric directions 1115 :type cpsize: list, tuple 1116 :param rs: minimum (r1) and maximum (r2) knot spans for the u-direction 1117 :type rs: list, tuple 1118 :param ss: minimum (s1) and maximum (s2) knot spans for the v-direction 1119 :type ss: list, tuple 1120 :param deriv_order: derivative order, i.e. the i-th derivative 1121 :type deriv_order: int 1122 :return: control points of the derivative surface over the input knot span ranges 1123 :rtype: list 1124 """ 1125 # Initialize return value (control points) 1126 PKL = [[[[[None for _ in range(dim)] 1127 for _ in range(cpsize[1])] for _ in range(cpsize[0])] 1128 for _ in range(deriv_order + 1)] for _ in range(deriv_order + 1)] 1129 1130 du = min(degree[0], deriv_order) 1131 dv = min(degree[1], deriv_order) 1132 1133 r = rs[1] - rs[0] 1134 s = ss[1] - ss[0] 1135 1136 # Control points of the U derivatives of every U-curve 1137 for j in range(ss[0], ss[1] + 1): 1138 PKu = curve_deriv_cpts(dim=dim, 1139 degree=degree[0], 1140 kv=kv[0], 1141 cpts=[cpts[j + (cpsize[1] * i)] for i in range(cpsize[0])], 1142 rs=rs, 1143 deriv_order=du) 1144 1145 # Copy into output as the U partial derivatives 1146 for k in range(0, du + 1): 1147 for i in range(0, r - k + 1): 1148 PKL[k][0][i][j - ss[0]] = PKu[k][i] 1149 1150 # Control points of the V derivatives of every U-differentiated V-curve 1151 for k in range(0, du): 1152 for i in range(0, r - k + 1): 1153 dd = min(deriv_order - k, dv) 1154 1155 PKuv = curve_deriv_cpts(dim=dim, 1156 degree=degree[1], 1157 kv=kv[1][ss[0]:], 1158 cpts=PKL[k][0][i], 1159 rs=(0, s), 1160 deriv_order=dd) 1161 1162 # Copy into output 1163 for l in range(1, dd + 1): 1164 for j in range(0, s - l + 1): 1165 PKL[k][l][i][j] = PKuv[l][j] 1166 1167 # Return control points 1168 return PKL 1169