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 Bézier curve see :ref:`this section <bezier-geometry-ressources-page>`.
####################################################################################################
#
# Notes: algorithm details are on bezier.rst
#
####################################################################################################
# Fixme:
# max distance to the chord for linear approximation
# fitting
# C0 = continuous
# G1 = geometric continuity
# Tangents point to the same direction
# C1 = parametric continuity
# Tangents are the same, implies G1
# C2 = curvature continuity
# Tangents and their derivatives are the same
####################################################################################################
__all__ = [
'QuadraticBezier2D',
'CubicBezier2D',
]
####################################################################################################
from Patro.Common.Math.Root import quadratic_root, cubic_root, fifth_root
from .Interpolation import interpolate_two_points
from .Primitive import Primitive3P, Primitive4P, PrimitiveNP, Primitive2DMixin
from .Transformation import AffineTransformation
from .Vector import Vector2D
####################################################################################################
_module_logger = logging.getLogger(__name__)
####################################################################################################
class BezierMixin2D(Primitive2DMixin):
"""Mixin to implements 2D Bezier Curve."""
LineInterpolationPrecision = 0.05
_logger = _module_logger.getChild('BezierMixin2D')
##############################################
def interpolated_length(self, dt=None):
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
142
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
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
if dt is None:
dt = self.LineInterpolationPrecision / (self.end_point - self.start_point).magnitude
length = 0
t = 0
while t < 1:
t0 = t
t = min(t + dt, 1)
length += (self.point_at_t(t) - self.point_at_t(t0)).magnitude
return length
##############################################
def length_at_t(self, t, cache=False):
"""Compute the length of the curve at *t*."""
if cache: # lookup cache
if not hasattr(self, '_length_cache'):
self._length_cache = {}
length = self._length_cache.get(t, None)
if length is not None:
return length
length = self.split_at_t(t).length
if cache: # save
self._length_cache[t] = length
return length
##############################################
def t_at_length(self, length, precision=1e-6):
"""Compute t for the given length. Length must lie in [0, curve length] range]. """
if length < 0:
raise ValueError('Negative length')
if length == 0:
return 0
curve_length = self.length # Fixme: cache ?
if (curve_length - length) <= precision:
return 1
if length > curve_length:
raise ValueError('Out of length')
# Search t for length using dichotomy
# convergence rate :
# 10 iterations corresponds to curve length / 1024
# 16 / 65536
# start range
inf = 0
sup = 1
while True:
middle = (sup + inf) / 2
length_at_middle = self.length_at_t(middle, cache=True) # Fixme: out of memory, use LRU ???
# exit condition
if abs(length_at_middle - length) <= precision:
return middle
elif length_at_middle < length:
inf = middle
else: # length < length_at_middle
sup = middle
##############################################
def split_at_two_t(self, t1, t2):
if t1 == t2:
return self.point_at_t(t1)
if t2 < t1:
# Fixme: raise ?
t1, t2 = t2, t1
# curve = self
# if t1 > 0:
curve = self.split_at_t(t1)[1] # right
if t2 < 1:
# Interpolate the parameter at t2 in the new curve
t = (t2 - t1) / (1 - t1)
curve = curve.split_at_t(t)[0] # left
return curve
##############################################
def _map_to_line(self, line):
transformation = AffineTransformation.Rotation(-line.v.orientation)
# Fixme: use __vector_cls__
transformation *= AffineTransformation.Translation(Vector2D(0, -line.p.y))
# Fixme: better API ?
return self.clone().transform(transformation)
##############################################
def non_parametric_curve(self, line):
"""Return the non-parametric Bezier curve D(ti, di(t)) where di(t) is the distance of the curve from
the baseline of the fat-line, ti is equally spaced in [0, 1].
"""
ts = np.arange(0, 1, 1/(self.number_of_points-1))
distances = [line.distance_to_line(p) for p in self.points]
points = [Vector2D(t, d) for t, f in zip(ts, distances)]
return self.__class__(*points)
##############################################
def distance_to_point(self, point):
p = self.closest_point(point)
if p is not None:
return (point - p).magnitude
else:
return None
####################################################################################################
class QuadraticBezier2D(BezierMixin2D, Primitive3P):
"""Class to implements 2D Quadratic Bezier Curve."""
BASIS = np.array((
(1, -2, 1),
(0, 2, -2),
(0, 0, 1),
))
INVERSE_BASIS = np.array((
(-2, 1, -2),
(-1, -3, 1),
(-1, -1, -2),
))
_logger = _module_logger.getChild('QuadraticBezier2D')
##############################################
def __init__(self, p0, p1, p2):
##############################################
def __repr__(self):
return self.__class__.__name__ + '({0._p0}, {0._p1}, {0._p2})'.format(self)
##############################################
@property
def length(self):
r"""Compute the length of the curve.
For more details see :ref:`this section <bezier-curve-length-section>`.
A0 = self._p1 - self._p0
A1 = self._p0 - self._p1 * 2 + self._p2
c = 4 * A1.dot(A1)
b = 8 * A0.dot(A1)
a = 4 * A0.dot(A0)
q = 4 * a * c - b * b
two_cb = 2 * c + b
sum_cba = c + b + a
m0 = 0.25 / c
m1 = q / (8 * c**1.5)
return (m0 * (two_cb * sqrt(sum_cba) - b * sqrt(a)) +
m1 * (log(2 * sqrt(c * sum_cba) + two_cb) - log(2 * sqrt(c * a) + b)))
else:
##############################################
def point_at_t(self, t):
# if 0 < t or 1 < t:
# raise ValueError()
u = 1 - t
return self._p0 * u**2 + self._p1 * 2 * t * u + self._p2 * t**2
##############################################
def split_at_t(self, t):
if t <= 0:
return None, self
elif t >= 1:
return self, None
else:
p01 = interpolate_two_points(self._p0, self._p1, t)
p12 = interpolate_two_points(self._p1, self._p2, t)
p = interpolate_two_points(p01, p12, t) # p = p012
# p = self.point_at_t(t)
return (QuadraticBezier2D(self._p0, p01, p), QuadraticBezier2D(p, p12, self._p2))
@property
def tangent0(self):
return (self._p1 - self._p0).normalise()
##############################################
@property
def tangent1(self):
return (self._p2 - self._p1).normalise()
##############################################
@property
def normal0(self):
return self.tangent0.normal()
##############################################
@property
def tangent1(self):
return self.tangent1.normal()
##############################################
def tangent_at(self, t):
u = 1 - t
return (self._p1 - self._p0) * u + (self._p2 - self._p1) * t
##############################################
def intersect_line(self, line):
For more details see :ref:`this section <bezier-curve-line-intersection-section>`.
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
curve = self._map_to_line(line)
p0 = curve.p0.y
p1 = curve.p1.y
p2 = curve.p2.y
return quadratic_root(
p2 - 2*p1 + p0, # t**2
2*(p1 - p0), # t
p0,
)
### a = p0 - 2*p1 + p2 # t**2
### # b = 2*(-p0 + p1) # t
### b = -p0 + p1 # was / 2 !!!
### c = p0
###
### # discriminant = b**2 - 4*a*c
### # discriminant = 4 * (p1**2 - p0*p2)
### discriminant = p1**2 - p0*p2 # was / 4 !!!
###
### if discriminant < 0:
### return None
### elif discriminant == 0:
### return -b / a # dropped 2
### else:
### # dropped 2
### y1 = (-b - sqrt(discriminant)) / a
### y2 = (-b + sqrt(discriminant)) / a
### return y1, y2
##############################################
def fat_line(self):
line = Line2D.from_two_points(self._p0, self._p3)
d1 = line.distance_to_line(self._p1)
d_min = min(0, d1 / 2)
d_max = max(0, d1 / 2)
return (line, d_min, d_max)
##############################################
def closest_point(self, point):
"""Return the closest point on the curve to the given *point*.
For more details see :ref:`this section <bezier-curve-closest-point-section>`.
A = self._p1 - self._p0
B = self._p2 - self._p1 - A
M = self._p0 - point
roots = cubic_root(
B.magnitude_square,
3*A.dot(B),
2*A.magnitude_square + M.dot(B),
M.dot(A),
)
t = [root for root in roots if 0 <= root <= 1]
if not t:
return None
elif len(t) > 1:
self._logger.warning("Found more than one root {} for {} and point {}".format(t, self, point))
return None
else:
return self.point_at_t(t)
##############################################
def to_cubic(self):
r"""Elevate the quadratic Bézier curve to a cubic Bézier cubic with the same shape.
For more details see :ref:`this section <bezier-curve-degree-elevation-section>`.
"""
p1 = (self._p0 + self._p1 * 2) / 3
p2 = (self._p2 + self._p1 * 2) / 3
####################################################################################################
_Sqrt3 = sqrt(3)
_Div18Sqrt3 = 18 / _Sqrt3
_OneThird = 1 / 3
_Sqrt3Div36 = _Sqrt3 / 36
class CubicBezier2D(BezierMixin2D, Primitive4P):
"""Class to implements 2D Cubic Bezier Curve."""
InterpolationPrecision = 0.001
BASIS = np.array((
(1, -3, 3, -1),
(0, 3, -6, 3),
(0, 0, 3, -3),
(0, 0, 0, 1),
))
INVERSE_BASIS = np.array((
(1, 1, 1, 1),
(0, 1/3, 2/3, 1),
(0, 0, 1/3, 1),
(0, 0, 0, 1),
))
_logger = _module_logger.getChild('CubicMixin2D')
#######################################
def __init__(self, p0, p1, p2, p3):
##############################################
def __repr__(self):
return self.__class__.__name__ + '({0._p0}, {0._p1}, {0._p2}, {0._p3})'.format(self)
##############################################
from .Spline import CubicUniformSpline2D
basis = np.dot(self.BASIS, CubicUniformSpline2D.INVERSE_BASIS)
points = np.dot(self.point_array, basis).transpose()
return CubicUniformSpline2D(*points)
##############################################
@property
def length(self):
return self.adaptive_length_approximation()
##############################################
def point_at_t(self, t):
# if 0 < t or 1 < t:
# raise ValueError()
return (self._p0 +
(self._p1 - self._p0) * 3 * t +
(self._p2 - self._p1*2 + self._p0) * 3 * t**2 +
(self._p3 - self._p2*3 + self._p1*3 - self._p0) * t**3)
# interpolate = point_at_t
##############################################
def _q_point(self):
"""Return the control point for mid-point quadratic approximation"""
return (self._p2*3 - self._p3 + self._p1*3 - self._p0) / 4
##############################################
def mid_point_quadratic_approximation(self):
"""Return the mid-point quadratic approximation"""
p1 = self._q_point()
return QuadraticBezier2D(self._p0, p1, self._p3)
##############################################
def split_at_t(self, t):
p01 = interpolate_two_points(self._p0, self._p1, t)
p12 = interpolate_two_points(self._p1, self._p2, t)
p23 = interpolate_two_points(self._p2, self._p3, t)
p012 = interpolate_two_points(p01, p12, t)
p123 = interpolate_two_points(p12, p23, t)
p = interpolate_two_points(p012, p123, t) # p = p0123
# p = self.point_at_t(t)
return (CubicBezier2D(self._p0, p01, p012, p), CubicBezier2D(p, p123, p23, self._p3))
##############################################
def _d01(self):
"""Return the distance between 0 and 1 quadratic aproximations"""
return (self._p3 - self._p2 * 3 + self._p1 * 3 - self._p0).magnitude / 2
##############################################
def _t_max(self):
"""Return the split point for adaptive quadratic approximation"""
return (_Div18Sqrt3 * self.InterpolationPrecision / self._d01())**_OneThird
##############################################
def q_length(self):
"""Return the length of the mid-point quadratic approximation"""
return self.mid_point_quadratic_approximation().length
##############################################
def adaptive_length_approximation(self):
"""Return the length of the adaptive quadratic approximation"""
segments = []
segment = self
t_max = segment._t_max()
while t_max < 1:
split = segment.split_at_t(t_max)
segments.append(split[0])
segment = split[1]
t_max = segment._t_max()
segments.append(segment)
return sum([segment.q_length() for segment in segments])
##############################################
@property
def tangent1(self):
return (self._p3 - self._p2).normalise()
##############################################
def tangent_at(self, t):
u = 1 - t
return (self._p1 - self._p0) * u**2 + (self._p2 - self._p1) * 2 * t * u + (self._p3 - self._p2) * t**2
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
##############################################
def intersect_line(self, line):
"""Find the intersections of the curve with a line."""
# Algorithm: same as for quadratic
# u = 1 - t
# B = p0 * u**3 +
# p1 * 3 * u**2 * t +
# p2 * 3 * u * t**2 +
# p3 * t**3
# B = p0 +
# (p1 - p0) * 3 * t +
# (p2 - p1 * 2 + p0) * 3 * t**2 +
# (p3 - p2 * 3 + p1 * 3 - p0) * t**3
# solveset(B, t)
curve = self._map_to_line(line)
p0 = curve.p0.y
p1 = curve.p1.y
p2 = curve.p2.y
p3 = curve.p3.y
return cubic_root(
p3 - 3*p2 + 3*p1 - p0,
3 * (p2 - p1 * 2 + p0),
3 * (p1 - p0),
p0,
)
##############################################
def fat_line(self):
line = Line2D.from_two_points(self._p0, self._p3)
d1 = line.distance_to_line(self._p1)
d2 = line.distance_to_line(self._p2)
if d1*d2 > 0:
factor = 3 / 4
else:
factor = 4 / 9
d_min = factor * min(0, d1, d2)
d_max = factor * max(0, d1, d2)
return (line, d_min, d_max)
##############################################
def _clipping_convex_hull(self):
line_03 = Line2D(self._p0, self._p3)
d1 = line_03.distance_to_line(self._p1)
d2 = line_03.distance_to_line(self._p2)
# Check if p1 and p2 are on the same side of the line [p0, p3]
if d1 * d2 < 0:
# p1 and p2 lie on different sides of [p0, p3].
# The hull is a quadrilateral and line [p0, p3] is not part of the hull.
# The top part includes p1, we will reverse it later if that is not the case.
hull = [
[self._p0, self._p1, self._p3], # top part
[self._p0, self._p2, self._p3] # bottom part
]
flip = d1 < 0
else:
# p1 and p2 lie on the same sides of [p0, p3]. The hull can be a triangle or a
# quadrilateral and line [p0, p3] is part of the hull. Check if the hull is a triangle
# or a quadrilateral. Also, if at least one of the distances for p1 or p2, from line
# [p0, p3] is zero then hull must at most have 3 vertices.
# Fixme: check cross product
y0, y1, y2, y3 = [p.y for p in self.points]
if abs(d1) < abs(d2):
pmax = p2;
# apex is y0 in this case, and the other apex point is y3
# vector yapex -> yapex2 or base vector which is already part of the hull
# V30xV10 * V10xV20
cross_product = ((y1 - y0) - (y3 - y0)/3) * (2*(y1 - y0) - (y2 - y0)) /3
else:
pmax = p1;
# apex is y3 and the other apex point is y0
# vector yapex -> yapex2 or base vector which is already part of the hull
# V32xV30 * V32xV31
cross_product = ((y3 - y2) - (y3 - y0)/3) * (2*(y3 - y2) - (y3 + y1)) /3
# Compare cross products of these vectors to determine if the point is in the triangle
# [p3, pmax, p0], or if it is a quadrilateral.
has_null_distance = d1 == 0 or d2 == 0 # Fixme: don't need to compute cross_product
if cross_product < 0 or has_null_distance:
# hull is a triangle
hull = [
[self._p0, pmax, self._p3], # top part is a triangle
[self._p0, self._p3], # bottom part is just an edge
]
else:
hull = [
[self._p0, self._p1, self._p2, self._p3], # top part is a quadrilateral
[self._p0, self._p3], # bottom part is just an edge
]
flip = d1 < 0 if d1 else d2 < 0
if flip:
hull.reverse()
return hull
##############################################
@staticmethod
def _clip_convex_hull(hull_top, hull_bottom, d_min, d_max) :
# Top /----
# / ---/
# / /
# d_max -------------------*---
# / / t_max
# t_min / /
# d_min -------*---------------
# / /
# / ----/ Bottom
# p0 /----
if (hull_top[0].y < d_min):
# Left of hull is below d_min,
# walk through the hull until it enters the region between d_min and d_max
return self._clip_convex_hull_part(hull_top, True, d_min);
elif (hull_bottom[0].y > d_max) :
# Left of hull is above d_max,
# walk through the hull until it enters the region between d_min and d_max
return self._clip_convex_hull_part(hull_bottom, False, d_max);
else :
# Left of hull is between d_min and d_max, no clipping possible
return hull_top[0].x; # Fixme: 0 ???
##############################################
@staticmethod
def _clip_convex_hull_part(part, top, threshold) :
"""Clip the bottom or top part of the convex hull.
*part* is a list of points, *top* is a boolean flag to indicate if it corresponds to the top
part, *threshold* is d_min if top part else d_max.
"""
# Walk on the edges
px = part[0].x;
py = part[0].y;
for i in range(1, len(part)):
qx = part[i].x;
qy = part[i].y;
if (qy >= threshold if top else qy <= threshold):
# compute a linear interpolation
# threshold = s * (t - px) + py
# t = (threshold - py) / s + px
return px + (threshold - py) * (qx - px) / (qy - py);
px = qx;
py = qy;
return None; # no intersection
##############################################
@staticmethod
def _instersect_curve(
curve1, curve2,
t_min=0, t_max=1,
u_min=0, u_max=1,
old_delta_t=1,
reverse=False, # flag to indicate that 1 <-> 2 when we store locations
recursion=0, # number of recursions
recursion_limit=32,
t_limit=0.8,
locations=[],
) :
"""Compute the intersection of two Bézier curves.
Code inspired from
* https://github.com/paperjs/paper.js/blob/master/src/path/Curve.js
* http://nbviewer.jupyter.org/gist/hkrish/0a128f21a5b9e5a7a914 The Bezier Clipping Algorithm
* https://gist.github.com/hkrish/5ef0f2da7f9882341ee5 hkrish/bezclip_manual.py
"""
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
# Note:
# see https://github.com/paperjs/paper.js/issues/565
# It was determined that more than 20 recursions are needed sometimes, depending on the
# delta_t threshold values further below when determining which curve converges the
# least. He also recommended a threshold of 0.5 instead of the initial 0.8
if recursion > recursion_limit:
return
tolerance = 1e-5
epsilon = 1e-10
# t_min_new = 0.
# t_max_new = 0.
# delta_t = 0.
# NOTE: the recursion threshold of 4 is needed to prevent this issue from occurring:
# https://github.com/paperjs/paper.js/issues/571
# when two curves share an end point
if curve1.p0.x == curve1.p3.x and u_max - u_min <= epsilon and recursion > 4:
# The fat-line of curve1 has converged to a point, the clipping is not reliable.
# Return the value we have even though we will miss the precision.
t_max_new = t_min_new = (t_max + t_min) / 2
delta_t = 0
else :
# Compute the fat-line for curve1:
# a baseline and two offsets which completely encloses the curve
fatline, d_min, d_max = curve1.fat_line()
# Calculate a non-parametric bezier curve D(ti, di(t)) where di(t) is the distance of curve2 from
# the baseline, ti is equally spaced in [0, 1]
non_parametric_curve = curve2.non_parametric_curve(fatline)
# Get the top and bottom parts of the convex-hull
top, bottom = non_parametric_curve._clip_convex_hull()
# Clip the convex-hull with d_min and d_max
t_min_clip = self.clip_convex_hull(top, bottom, d_min, d_max);
top.reverse()
bottom.reverse()
t_max_clip = clipConvexHull(top, bottom, d_min, d_max);
# No intersections if one of the t values is None
if t_min_clip is None or t_max_clip is None:
return
# Clip curve2 with the fat-line for curve1
curve2 = curve2.split_at_two_t(t_min_clip, t_max_clip)
delta_t = t_max_clip - t_min_clip
# t_min and t_max are within the range [0, 1]
# We need to project it to the original parameter range
t_min_new = t_max * t_min_clip + t_min * (1 - t_min_clip)
t_max_new = t_max * t_max_clip + t_min * (1 - t_max_clip)
delta_t_new = t_max_new - t_min_new
delta_u = u_max - u_min
# Check if we need to subdivide the curves
if old_delta_t > t_limit and delta_t > t_limit:
# Subdivide the curve which has converged the least.
args = (delta_t, not reverse, recursion+1, recursion_limit, t_limit, locations)
if delta_u < delta_t_new: # curve2 < curve1
parts = curve1.split_at_t(0.5)
t = t_min_new + delta_t_new / 2
self._intersect_curve(curve2, parts[0], u_min, u_max, t_min_new, t, *args)
self._intersect_curve(curve2, parts[1], u_min, u_max, t, t_max_new, *args)
else :
parts = curve2.split_at_t(0.5)
t = u_min + delta_u / 2
self._intersect_curve(parts[0], curve1, u_min, t, t_min_new, t_max_new, *args)
self._intersect_curve(parts[1], curve1, t, u_max, t_min_new, t_max_new, *args)
elif max(delta_u, delta_t_new) < tolerance:
# We have isolated the intersection with sufficient precision
t1 = t_min_new + delta_t_new / 2
t2 = u_min + delta_u / 2
if reverse:
t1, t2 = t2, t1
p1 = curve1.point_at_t(t1)
p2 = curve2.point_at_t(t2)
locations.append([t1, point1, t2, point2])
else:
args = (delta_t, not reverse, recursion+1, recursion_limit, t_limit)
self._intersect_curve(curve2, curve1, locations, u_min, u_max, t_min_new, t_max_new, *args)
##############################################
def is_flat_enough(self, flatness):
r"""Determines if a curve is sufficiently flat, meaning it appears as a straight line and has
curve-time that is enough linear, as specified by the given *flatness* parameter.
For more details see :ref:`this section <bezier-curve-flatness-section>`.
u = 3*self._p1 - 2*self._p0 - self._p3
v = 3*self._p2 - 2*self._p3 - self._p0
criterion = max(u.x**2, v.x**2) + max(u.y**2, v.y**2)
threshold = 16 * flatness**2
self._logger.warning("is flat {} <= {} with flatness {}".format(criterion, threshold, flatness))
return criterion <= threshold
##############################################
@property
def area(self):
"""Compute the area delimited by the curve and the segment across the start and stop point."""
# Reference: http://objectmix.com/graphics/133553-area-closed-bezier-curve.html BUT DEAD LINK
# Proof using divergence theorem ???
# Fixme: any proof !
x0, y0 = list(self._p0)
x1, y1 = list(self._p1)
x2, y2 = list(self._p2)
x3, y3 = list(self._p3)
return (3 * ((y3 - y0) * (x1 + x2) - (x3 - x0) * (y1 + y2)
+ y1 * (x0 - x2) - x1 * (y0 - y2)
+ y3 * (x2 + x0 / 3) - x3 * (y2 + y0 / 3)) / 20)
##############################################
def closest_point(self, point):
"""Return the closest point on the curve to the given *point*.
For more details see :ref:`this section <bezier-curve-closest-point-section>`.
"""
n = self._p3 - self._p2*3 + self._p1*3 - self._p0
r = (self._p2 - self._p1*2 + self._p0)*3
s = (self._p1 - self._p0)*3
v = self._p0
roots = fifth_root(
-3 * n.magnitude_square,
-5 * n.dot(r),
-2 * (2*n.dot(s) + r.magnitude_square),
3 * (point.dot(n) - n.dot(v) - r.dot(s)),
2*point.dot(r) - 2*r.dot(v) - s.magnitude_square,
point.dot(s) - s.dot(v),
)
# Fixme: to func
t = [root for root in roots if 0 <= root <= 1]
if not t:
return None
elif len(t) > 1:
# Fixme:
# Found more than one root [0, 0.516373783749732]
# for CubicBezier2D(
# Vector2D[1394.4334 1672.0004], Vector2D[1394.4334 1672.0004],
# Vector2D[1585.0004 1624.9634], Vector2D[1585.0004 1622.0004])
# and point Vector2D[1495.11502887 1649.7386517 ]
# raise NameError("Found more than one root: {}".format(t))
self._logger.warning("Found more than one root {} for {} and point {}".format(t, self, point))
# self._logger.warning("is flat {}".format(self.is_flat_enough(.1)))
if len(t) == 2 and t[0] == 0:
return self.point_at_t(t[1])
else:
return None