Skip to content
Bezier.py 33.7 KiB
Newer Older
####################################################################################################
#
Fabrice Salvaire's avatar
Fabrice Salvaire committed
# 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/>.
#
####################################################################################################

# 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

####################################################################################################

from math import log, sqrt, pow
import numpy as np

from Patro.Common.Math.Root import quadratic_root, cubic_root, fifth_root
Fabrice Salvaire's avatar
Fabrice Salvaire committed
from .Interpolation import interpolate_two_points
from .Line import Line2D
from .Primitive import Primitive3P, Primitive4P, Primitive2DMixin
from .Transformation import AffineTransformation
from .Vector import Vector2D

####################################################################################################

# Fixme:
#   max distance to the chord for linear approximation
#   fitting

####################################################################################################

class QuadraticBezier2D(Primitive2DMixin, Primitive3P):
Fabrice Salvaire's avatar
Fabrice Salvaire committed
    """Class to implements 2D Quadratic Bezier Curve."""

    LineInterpolationPrecision = 0.05

    ##############################################

    def __init__(self, p0, p1, p2):

        Primitive3P.__init__(self, p0, p1, p2)
Fabrice Salvaire's avatar
Fabrice Salvaire committed

    ##############################################

    def __repr__(self):
        return self.__class__.__name__ + '({0._p0}, {0._p1}, {0._p2})'.format(self)

    ##############################################

    @property
    def length(self):

Fabrice Salvaire's avatar
Fabrice Salvaire committed
        # Algorithm:
        #
        # http://www.gamedev.net/topic/551455-length-of-a-generalized-quadratic-bezier-curve-in-3d
        # Dave Eberly Posted October 25, 2009
        #
        # The quadratic Bezier is
        #   (x(t),y(t)) = (1-t)^2*(x0,y0) + 2*t*(1-t)*(x1,y1) + t^2*(x2,y2)
        #
        # The derivative is
        #   (x'(t),y'(t)) = -2*(1-t)*(x0,y0) + (2-4*t)*(x1,y1) + 2*t*(x2,y2)
        #
        # The length of the curve for 0 <= t <= 1 is
        #   Integral[0,1] sqrt((x'(t))^2 + (y'(t))^2) dt
        # The integrand is of the form sqrt(c*t^2 + b*t + a)
        #
        # You have three separate cases: c = 0, c > 0, or c < 0.
        # * The case c = 0 is easy.
        # * For the case c > 0, an antiderivative is
        #     (2*c*t+b)*sqrt(c*t^2+b*t+a)/(4*c) + (0.5*k)*log(2*sqrt(c*(c*t^2+b*t+a)) + 2*c*t + b)/sqrt(c)
        #   where k = 4*c/q with q = 4*a*c - b*b.
        # * For the case c < 0, an antiderivative is
        #    (2*c*t+b)*sqrt(c*t^2+b*t+a)/(4*c) - (0.5*k)*arcsin((2*c*t+b)/sqrt(-q))/sqrt(-c)

        A0 = self._p1 - self._p0
        A1 = self._p0 - self._p1 * 2 + self._p2
        if A1.magnitude_square != 0:
            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:
            return 2 * A0.magnitude

    ##############################################

    def interpolated_length(self, dt=None):

        # Length of the curve obtained via line interpolation

        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 point_at_t(self, t):
        # if 0 < t or 1 < t:
        #     raise ValueError()
Fabrice Salvaire's avatar
Fabrice Salvaire committed
        u = 1 - t
        return self._p0 * u**2 + self._p1 * 2 * t * u + self._p2 * t**2

    ##############################################

    def split_at_t(self, t):

Fabrice Salvaire's avatar
Fabrice Salvaire committed
        """Split the curve at given position"""
        if t <= 0:
            return None, self
        elif t >= 1:
            return self, None
Loading full blame...