Bezier.py

        #        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

        """

        # 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.

         *flatness* is the maximum error allowed for the straight line to deviate from the curve.

         Algorithm

         We define the flatness of the curve as the argmax of the distance from the curve to the
         line passing by the start and stop point.

         :math:`\mathrm{flatness} = argmax(d(t))` for :math:`t \in [0, 1]` where :math:`d(t) = \vert B(t) - L(t) \vert`

         The line equation is

         .. math::
             L = (1-t) \mathbf{P}_0 + t \mathbf{P}_1

         Let

         .. math::
             \begin{align}
             u &= 3\mathbf{P}_1 - 2\mathbf{P}_0 - \mathbf{P}_3 \\
             v &= 3\mathbf{P}_2 - \mathbf{P}_0 - 2\mathbf{P}_3
             \end{align}

         The distance is

         .. math::
             \begin{align}
             d(t) &= (1-t)^2 t \left(3\mathbf{P}_1 - 2\mathbf{P}_0 - \mathbf{P}_3\right)  +  (1-t) t^2 (3\mathbf{P}_2 - \mathbf{P}_0 - 2\mathbf{P}_3) \\
                  &= (1-t)^2 t u  +  (1-t) t^2 v
             \end{align}

         The square of the distance is

         .. math::
             d(t)^2 = (1 - t)^2 t^2 (((1 - t) ux + t vx)^2 + ((1 - t) uy + t vy)^2

         From

         .. math::
             \begin{align}
             argmax((1 - t)^2 t^2)   &= \frac{1}{16} \\
             argmax((1 - t) a + t b) &= argmax(a, b)
             \end{align}

         we can express a bound on the flatness

         .. math::
             \mathrm{flatness}^2 = argmax(d(t)^2) \leq \frac{1}{16} (argmax(ux^2, vx^2) + argmax(uy^2, vy^2))

         Thus an upper bound of :math:`16\,\mathrm{flatness}^2` is

         .. math::
             argmax(ux^2, vx^2) + argmax(uy^2, vy^2)

         Reference

         * Kaspar Fischer and Roger Willcocks  http://hcklbrrfnn.files.wordpress.com/2012/08/bez.pdf
         * PostScript Language Reference. Addison- Wesley, third edition, 1999

         """

         # x0, y0 = list(self._p0)
         # x1, y1 = list(self._p1)
         # x2, y2 = list(self._p2)
         # x3, y3 = list(self._p3)

         # ux = 3*x1 - 2*x0 - x3
         # uy = 3*y1 - 2*y0 - y3
         # vx = 3*x2 - 2*x3 - x0
         # vy = 3*y2 - 2*y3 - y0

         u = 3*P1 - 2*P0 - P3
         v = 3*P2 - 2*P3 - P0

         return max(u.x**2, v.x**2) + max(u.y**2, v.y**2) <= 16 * flatness**2

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

    @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):

        # n = P3 - 3*P2 + 3*P1 - P0
        # r = 3*(P2 - 2*P1 + P0
        # s = 3*(P1 - P0)
        # v = P0

        # Q(t)  = n*t**3 + r*t**2 + s*t + v
        # Q'(t) = 3*n*t**2 + 2*r*t + s

        # n, r, s, v = symbols('n r s v')
        # Q = n*t**3 + r*t**2 + s*t + v
        # Qp = simplify(Q.diff(t))
        # collect(expand((P*Qp - Q*Qp)), t)

        # -3*n**2 * t**5
        # -5*n*r * t**4
        # -2*(2*n*s + r**2) * t**3
        # 3*(P*n - n*v - r*s) * t**2
        # (2*P*r - 2*r*v - s**2) * t
        # P*s - s*v

        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:
            raise NameError("Found more than one root: {}".format(t))
        else:
            return self.point_at_t(t[0])