Bezier.py

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

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

         Reference

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

         """

         # 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.
         #
         # flatness = argmax(d(t)) for t in [0, 1] where d(t) = | B(t) - L(t) |
         #
         # L = (1-t)*P0 + t*P1
         #
         # Let
         #   u = 3*P1 - 2*P0 - P3
         #   v = 3*P2 - P0 - 2*P3
         #
         # d(t) = (1-t)**2 * t * (3*P1 - 2*P0 - P3)  +  (1-t) * t**2 * (3*P2 - P0 - 2*P3)
         #      = (1-t)**2 * t * u  +  (1-t) * t**2 * v
         #
         # d(t)**2 = (1 - t)**2 * t**2 * (((1 - t)*ux + t*vx)**2 + ((1 - t)*uy + t*vy)**2
         #
         # argmax((1 - t)**2 * t**2) = 1/16
         # argmax((1 - t)*a + t*b)   = argmax(a, b)
         #
         # flatness**2 = argmax(d(t)**2) <= 1/16 * (argmax(ux**2, vx**2) + argmax(uy**2, vy**2))
         #
         # argmax(ux**2, vx**2) + argmax(uy**2, vy**2) is thus an upper bound of 16 * flatness**2

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

        # Q(t) = (P3 - 3*P2 + 3*P1 - P0) * t**3 +
        #        3*(P2 - 2*P1 + P0) * t**2 +
        #        3*(P1 - P0) * t  +
        #        P0

        # 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

        # P0, P1, P2, P3, P, t = symbols('P0 P1 P2 P3 P t')
        # 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 on root")
        else:
            return self.point_at_t(t[0])