Newer
Older
####################################################################################################
#
# Patro - A Python library to make patterns for fashion design
# Copyright (C) 2017 Fabrice Salvaire
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
####################################################################################################
For resources on Spline curve see :ref:`this section <spline-geometry-ressources-page>`.
# The DeBoor-Cox algorithm permits to evaluate recursively a B-Spline in a similar way to the De
# Casteljaud algorithm for Bézier curves.
#
# Given `k` the degree of the B-spline, `n + 1` control points :math:`p_0, \ldots, p_n`, and an
# increasing series of scalars :math:`t_0 \le t_1 \le \ldots \le t_m` with :math:`m = n + k + 1`,
# called knots.
#
# The number of points must respect the condition :math:`n + 1 \le k`, e.g. a B-spline of degree 3
# must have 4 control points.
####################################################################################################
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
####################################################################################################
# from math import log, sqrt
import numpy as np
from .Bezier import QuadraticBezier2D, CubicBezier2D
from .Primitive import Primitive3P, Primitive4P, PrimitiveNP, Primitive2DMixin
####################################################################################################
class QuadraticUniformSpline2D(Primitive2DMixin, Primitive3P):
"""Class to implements 2D Quadratic Spline Curve."""
BASIS = np.array((
(1, -2, 1),
(1, 2, -2),
(0, 0, 1),
))
INVERSE_BASIS = np.array((
(-2, 1, -2),
(-2, -3, 1),
(-1, -1, -2),
))
#######################################
def __init__(self, p0, p1, p2):
Primitive3P.__init__(self, p0, p1, p2)
##############################################
def __repr__(self):
return self.__class__.__name__ + '({0._p0}, {0._p1}, {0._p2})'.format(self)
##############################################
def to_bezier(self):
basis = np.dot(self.BASIS, QuadraticBezier2D.INVERSE_BASIS)
points = np.dot(self.point_array, basis).transpose()
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
return QuadraticBezier2D(*points)
##############################################
def point_at_t(self, t):
# Q(t) = (
# P0 * (1-t)**3 +
# P1 * ( 3*t**3 - 6*t**2 + 4 ) +
# P2 * ( -3*t**3 + 3*t**2 + 3*t + 1 ) +
# P3 * t**3
# ) / 6
#
# = P0*(1-t)**3/6 + P1*(3*t**3 - 6*t**2 + 4)/6 + P2*(-3*t**3 + 3*t**2 + 3*t + 1)/6 + P3*t**3/6
return (self._p0/6 + self._p1*2/3 + self._p2/6 +
(-self._p0/2 + self._p2/2)*t +
(self._p0/2 - self._p1 + self._p2/2)*t**2 +
(-self._p0/6 + self._p1/2 - self._p2/2 + self._p3/6)*t**3)
####################################################################################################
class CubicUniformSpline2D(Primitive2DMixin, Primitive4P):
"""Class to implements 2D Cubic Spline Curve."""
# T = (1 t t**2 t**3)
# P = (Pi Pi+2 Pi+2 Pi+3)
# Q(t) = T M Pt
# = P Mt Tt
# Basis = Mt
BASIS = np.array((
(1, -3, 3, -1),
(4, 0, -6, 3),
(1, 3, 3, -3),
(0, 0, 0, 1),
)) / 6
INVERSE_BASIS = np.array((
( 1, 1, 1, 1),
( -1, 0, 1, 2),
(2/3, -1/3, 2/3, 11/3),
( 0, 0, 0, 6),
))
#######################################
def __init__(self, p0, p1, p2, p3):
Primitive4P.__init__(self, p0, p1, p2, p3)
##############################################
def __repr__(self):
return self.__class__.__name__ + '({0._p0}, {0._p1}, {0._p2}, {0._p3})'.format(self)
##############################################
def to_bezier(self):
basis = np.dot(self.BASIS, CubicBezier2D.INVERSE_BASIS)
points = np.dot(self.point_array, basis).transpose()
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
if self._start:
# list(self.points)[:2]
points[:2] = self._p0, self._p1
elif self._stop:
# list(self.points)[-2:]
points[-2:] = self._p2, self._p3
return CubicBezier2D(*points)
##############################################
def point_at_t(self, t):
# Q(t) = (
# P0 * (1-t)**3 +
# P1 * ( 3*t**3 - 6*t**2 + 4 ) +
# P2 * ( -3*t**3 + 3*t**2 + 3*t + 1 ) +
# P3 * t**3
# ) / 6
#
# = P0*(1-t)**3/6 + P1*(3*t**3 - 6*t**2 + 4)/6 + P2*(-3*t**3 + 3*t**2 + 3*t + 1)/6 + P3*t**3/6
return (self._p0/6 + self._p1*2/3 + self._p2/6 +
(-self._p0/2 + self._p2/2)*t +
(self._p0/2 - self._p1 + self._p2/2)*t**2 +
(-self._p0/6 + self._p1/2 - self._p2/2 + self._p3/6)*t**3)
####################################################################################################
class BSpline2D(Primitive2DMixin, PrimitiveNP):
"""Class to implement a 2D B-Spline curve.
"""
##############################################
@staticmethod
def uniform_knots(degree, number_of_points):
order = degree + 1
if number_of_points < order:
raise ValueError('Inconsistent degree and number of points')
knots = list(range(number_of_points - degree +1))
return [0]*degree + knots + [knots[-1]]*degree
##############################################
@classmethod
def check_for_unifom_knots(cls, degree, number_of_points, knots):
return np.array_equal(cls.uniform_knots(degree, number_of_points), knots)
##############################################
def __init__(self, points, degree, closed=False, knots=None):
points = self.handle_points(points)
PrimitiveNP.__init__(self, points)
self._closed = bool(closed) # Fixme: not implemented
if knots is not None:
self._knots = list(knots)
if not np.all(np.diff(self._knots) >= 0):
raise ValueError('Invalid knots {}'.format(knots))
# Fixme: check_for_unifom_knots
self._uniform = False # Fixme:
else:
self._knots = self.uniform_knots(self._degree, self.number_of_points)
self._uniform = True
# self._number_of_points = len(self._knots) - self._degree - 1
# assert (self.number_of_points >= self.order and
# (len(self._coefficients) >= self._number_of_points))
##############################################
@property
def degree(self):
return self._degree
@property
def order(self):
return self._degree +1
@property
def is_closed(self):
return self._closed
@property
def uniform(self):
return self._uniform
@property
def knots(self):
return self._knots
@property
def start_knot(self):
if self._uniform:
return 0
else:
return self._knots[0]
@property
def end_knot(self):
if self._uniform:
return self.number_of_points - self._degree
else:
return self._knots[-1]
@property
def knot_iter(self):
if self._uniform:
return range(self.end_knot)
else:
return iter(self._knots)
def number_of_spans(self):
if self._uniform:
count = self.number_of_points - self._degree
if self._closed:
count += 1
return count
else:
# multiplicity
raise NotImplementedError
##############################################
def span(self, t):
if not(self.start_knot <= t <= self.end_knot):
raise ValueError('Invalid t {}'.format(t))
if self._uniform:
return int(t) + self._degree # start padding
else:
for i, t_span in enumerate(self._knots):
if t < t_span:
return i-1
##############################################
def knot_multiplicity(self, knot):
count = 0
for t in self._knots:
if t == knot:
count += 1
elif t > knot:
break
return count
##############################################
def basis_function(self, i, k, t):
"""De Boor-Cox recursion formula"""
if k == 0:
return 1 if self._knots[i] <= t < self._knots[i+1] else 0
ki = self._knots[i]
kik = self._knots[i+k]
if kik == ki:
c1 = 0
else:
c1 = (t - ki)/(kik - ki) * self.basis_function(i, k-1, t)
ki = self._knots[i+1]
kik = self._knots[i+k+1]
if kik == ki:
c2 = 0
else:
c2 = (kik - t)/(kik - ki) * self.basis_function(i+1, k-1, t)
return c1 + c2
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
##############################################
def _deboor(self, t):
"""Compute point at t using De Boor algorithm"""
# l_minus_degree = int(t) # span index
# l = l_minus_degree + self._degree # knot index
l = self.span(t)
l_minus_degree = l - self._degree
knots = self._knots
points = [self._points[j + l_minus_degree] for j in range(self.order)]
for r in range(1, self.order):
for j in range(self._degree, r-1, -1):
k = j + l_minus_degree
d = knots[j+1+l-r] - knots[k]
if d == 0:
alpha = 0
else:
alpha = (t - knots[k]) / d
if alpha == 0:
point = points[j-1]
elif alpha == 1:
point = points[j]
else:
point = points[j-1] * (1 - alpha) + points[j] * alpha
points[j] = point
return points[self._degree]
##############################################
def _naive_point_at_t(self, t):
"""Compute point at t using a naive algorithm"""
basis = np.array([self.basis_function(i, self._degree, t)
for i in range(self.number_of_points)])
points = self.point_array
return self.__vector_cls__(np.dot(basis, points))
##############################################
def point_at_t(self, t, naive=False):
# Spline curve as a Bézier span at start and end
if self._uniform:
if t == 0:
return self.start_point
elif t == self.end_knot:
# else computation fail
return self.end_point
if naive:
return self._naive_point_at_t(t)
else:
return self._deboor(t)
##############################################
def insert_knot(self, t):
# http://pages.mtu.edu/~shene/COURSES/cs3621/NOTES/spline/B-spline/single-insertion.html
# t lie in the [t_l, t_{l+1}[ span
# this span is only affected by the control points: P_l, ..., P_{l-degree}
# (number of points = degree + 1 = order)
degree = self._degree
points = self._points
knots = self._knots
l = self.span(t)
new_points = []
for i in range(self.number_of_points +1):
# compute w
if i <= l - degree:
w = 1 # same point
elif i > l:
w = 0 # previous point
else: # blend previous and current point
ti = knots[i]
d = knots[i + degree] - ti
if d > 0:
w = (t - ti) / d
else:
w = 0
if w == 0:
point = points[i-1]
elif w == 1:
point = points[i]
else:
point = points[i-1] * (1 - w) + points[i] * w
new_points.append(point)
# l += degree
knots = self._knots[:l+1] + [t] + self._knots[l+1:]
return self.__class__(new_points, self._degree, knots=knots)
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
##############################################
def to_bezier_form(self, degree=None):
if not self._uniform:
raise ValueError('Must be uniform')
if degree is None:
degree = self._degree
new_spline = self
for t in range(1, self.end_knot):
multiplicity = self.knot_multiplicity(t)
for i in range(degree - multiplicity +1):
new_spline = new_spline.insert_knot(t)
return new_spline
##############################################
def to_bezier(self):
from . import Bezier
if self._degree == 2:
cls = Bezier.QuadraticBezier2D
elif self._degree >= 3:
cls = Bezier.CubicBezier2D
else:
not NotImplementedError
# Fixme: degree > 3
spline = self.to_bezier_form()
bezier_curves = []
order = spline.order
for i in range(self.number_of_spans):
i_inf = i * order
i_sup = i_inf + order
points = spline._points[i_inf:i_sup]
bezier = cls(*points)
bezier_curves.append(bezier)
return bezier_curves