diff --git a/Patro/GeometryEngine/Bezier.py b/Patro/GeometryEngine/Bezier.py
index 0e10b1909801204fccedb9523c10f8aa49a4e9e4..735e1314bfcc9cfe746b4b5c67e71ade2b75ef4b 100644
--- a/Patro/GeometryEngine/Bezier.py
+++ b/Patro/GeometryEngine/Bezier.py
@@ -24,6 +24,12 @@ For resources on Bézier curve see :ref:`this section <bezier-geometry-ressource
 
 """
 
+####################################################################################################
+#
+# Notes: algorithm details are on bezier.rst
+#
+####################################################################################################
+
 # Fixme:
 #   max distance to the chord for linear approximation
 #   fitting
@@ -307,18 +313,10 @@ class QuadraticBezier2D(BezierMixin2D, Primitive3P):
 
         """Find the intersections of the curve with a line.
 
-        Algorithm
-
-        * Apply a transformation to the curve that maps the line onto the X-axis.
-        * Then we only need to test the Y-values for a zero.
+        For more details see :ref:`this section <bezier-curve-line-intersection-section>`.
 
         """
 
-        # u = 1 - t
-        # B = p0 * u**2  +  p1 * 2*t*u  +  p2 * t**2
-        # collect(expand(B), t)
-        # solveset(B, t)
-
         curve = self._map_to_line(line)
 
         p0 = curve.p0.y
@@ -367,47 +365,10 @@ class QuadraticBezier2D(BezierMixin2D, Primitive3P):
 
         """Return the closest point on the curve to the given *point*.
 
-        Reference
-
-        * https://hal.archives-ouvertes.fr/inria-00518379/document
-          Improved Algebraic Algorithm On Point Projection For Bézier Curves
-          Xiao-Diao Chen, Yin Zhou, Zhenyu Shu, Hua Su, Jean-Claude Paul
+        For more details see :ref:`this section <bezier-curve-closest-point-section>`.
 
         """
 
-        # Condition:
-        #   (P - B(t)) . B'(t) = 0   where t in [0,1]
-        #
-        #   P. B'(t) - B(t). B'(t) = 0
-
-        # A = P1 - P0
-        # B = P2 - P1 - A
-        # M = P0 - P
-
-        # Q(t)  = P0*(1-t)**2 + P1*2*t*(1-t) + P2*t**2
-        # Q'(t) = -2*P0*(1 - t) + 2*P1*(1 - 2*t) + 2*P2*t
-        #       = 2*(A + B*t)
-
-        # Q = P0 * (1-t)**2  +  P1 * 2*t*(1-t)  +  P2 * t**2
-        # Qp = simplify(Q.diff(t))
-        # collect(expand((P*Qp - Q*Qp)/-2), t)
-
-        # (P0**2 - 4*P0*P1 + 2*P0*P2 + 4*P1**2 - 4*P1*P2 + P2**2) * t**3
-        # (-3*P0**2 + 9*P0*P1 - 3*P0*P2 - 6*P1**2 + 3*P1*P2) * t**2
-        # (-P*P0 + 2*P*P1 - P*P2 + 3*P0**2 - 6*P0*P1 + P0*P2 + 2*P1**2) * t
-        # P*P0 - P*P1 - P0**2 + P0*P1
-
-        # factorisation
-        # (P0 - 2*P1 + P2)**2 * t**3
-        # 3*(P1 - P0)*(P0 - 2*P1 + P2) * t**2
-        # ...
-        # (P0 - P)*(P1 - P0)
-
-        # B**2 * t**3
-        # 3*A*B * t**2
-        # (2*A**2 + M*B) * t
-        # M*A
-
         A = self._p1 - self._p0
         B = self._p2 - self._p1 - A
         M = self._p0 - point
@@ -433,15 +394,7 @@ class QuadraticBezier2D(BezierMixin2D, Primitive3P):
 
         r"""Elevate the quadratic Bézier curve to a cubic Bézier cubic with the same shape.
 
-        The new control points are
-
-        .. math::
-            \begin{align}
-            \mathbf{P'}_0 &= \mathbf{P}_0 \\
-            \mathbf{P'}_1 &= \mathbf{P}_0 + \frac{2}{3} (\mathbf{P}_1 - \mathbf{P}_0) \\
-            \mathbf{P'}_1 &= \mathbf{P}_2 + \frac{2}{3} (\mathbf{P}_1 - \mathbf{P}_2) \\
-            \mathbf{P'}_2 &= \mathbf{P}_2
-            \end{align}
+        For more details see :ref:`this section <bezier-curve-degree-elevation-section>`.
 
         """
 
@@ -876,16 +829,6 @@ class CubicBezier2D(BezierMixin2D, Primitive4P):
 
          """
 
-         # 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
 
@@ -915,26 +858,6 @@ class CubicBezier2D(BezierMixin2D, Primitive4P):
 
     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
diff --git a/Patro/GeometryEngine/Spline.py b/Patro/GeometryEngine/Spline.py
index aa9f06f1c0c2331385e5a58cc19dac46391b2bb3..c415503c3476dbc445f44378b991c8b38985890a 100644
--- a/Patro/GeometryEngine/Spline.py
+++ b/Patro/GeometryEngine/Spline.py
@@ -24,15 +24,11 @@ For resources on Spline curve see :ref:`this section <spline-geometry-ressources
 
 """
 
-# 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.
+# Notes: algorithm details are on spline.rst
 #
-# 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.
+####################################################################################################
 
 ####################################################################################################
 
diff --git a/doc/sphinx/source/resources/geometry/bezier.rst b/doc/sphinx/source/resources/geometry/bezier.rst
index e170554a26238a2fe4fa5505930c8529b5696e4c..8839380d324cab83fef3608fbb9302e3e7eb2617 100644
--- a/doc/sphinx/source/resources/geometry/bezier.rst
+++ b/doc/sphinx/source/resources/geometry/bezier.rst
@@ -8,6 +8,9 @@
 
 .. contents:: :local:
 
+
+
+
 Definitions
 -----------
 
@@ -16,22 +19,27 @@ A Bézier curve is defined by a set of control points :math:`\mathbf{P}_0` throu
 quadratic, 3 for cubic etc.). The first and last control points are always the end points of the
 curve;
 
-In the following :math:`0 \le t \le 1`.
+In the following :math:`0 \le t \le 1` and :math:`u = 1 - t`
+
+
 
 Linear Bézier Curves
----------------------
+--------------------
 
 Given distinct points :math:`\mathbf{P}_0` and :math:`\mathbf{P}_1`, a linear Bézier curve is simply
 a straight line between those two points. The curve is given by
 
 .. math::
     \begin{align}
-    \mathbf{B}(t) &= \mathbf{P}_0 + t (\mathbf{P}_1 - \mathbf{P}_0) \\
-                  &= (1-t) \mathbf{P}_0 + t \mathbf{P}_1
+    \mathbf{B}(t) &= \mathbf{P}_0 + t (\mathbf{P}_1 - \mathbf{P}_0) \\[.5em]
+                  &= (1-t) \mathbf{P}_0 + t \mathbf{P}_1 \\[.5em]
+                  &= u \mathbf{P}_0 + t \mathbf{P}_1
     \end{align}
 
 and is equivalent to linear interpolation.
 
+
+
 Quadratic Bézier Curves
 -----------------------
 
@@ -39,7 +47,10 @@ A quadratic Bézier curve is the path traced by the function :math:`\mathbf{B}(t
 :math:`\mathbf{P}_0`, :math:`\mathbf{P}_1`, and :math:`\mathbf{P}_2`,
 
 .. math::
-    \mathbf{B}(t) = (1 - t)[(1 - t) \mathbf{P}_0 + t \mathbf{P}_1] + t [(1 - t) \mathbf{P}_1 + t \mathbf{P}_2]
+    \begin{align}
+    \mathbf{B}(t) &= (1-t)[(1-t) \mathbf{P}_0 + t \mathbf{P}_1] + t [(1-t) \mathbf{P}_1 + t \mathbf{P}_2] \\[.5em]
+                  &= u[u \mathbf{P}_0 + t \mathbf{P}_1] + t [u \mathbf{P}_1 + t \mathbf{P}_2]
+    \end{align}
 
 which can be interpreted as the linear interpolant of corresponding points on the linear Bézier
 curves from :math:`\mathbf{P}_0` to :math:`\mathbf{P}_1` and from :math:`\mathbf{P}_1` to
@@ -49,21 +60,29 @@ Rearranging the preceding equation yields:
 
 .. math::
     \begin{align}
-    \mathbf{B}(t) &= (1 - t)^{2} \mathbf{P}_0 + 2(1 - t)t \mathbf{P}_1 + t^{2} \mathbf{P}_2 \\
+    \mathbf{B}(t) &= (1-t)^2 \mathbf{P}_0 + 2(1-t)t \mathbf{P}_1 + t^2 \mathbf{P}_2 \\[.5em]
+                  &= u^2 \mathbf{P}_0 + 2ut \mathbf{P}_1 + t^2 \mathbf{P}_2 \\[.5em]
                   &= (\mathbf{P}_0 - 2\mathbf{P}_1 + \mathbf{P}_2) t^2 +
-                     (-2\mathbf{P}_0 + 2\mathbf{P}_1) t +
+                     2(-\mathbf{P}_0 + \mathbf{P}_1) t +
                      \mathbf{P}_0
     \end{align}
 
 This can be written in a way that highlights the symmetry with respect to :math:`\mathbf{P}_1`:
 
 .. math::
-    \mathbf{B}(t) = \mathbf{P}_1 + (1 - t)^{2} ( \mathbf{P}_0 - \mathbf{P}_1) + t^{2} (\mathbf{P}_2 - \mathbf{P}_1)
+    \begin{align}
+    \mathbf{B}(t) &= \mathbf{P}_1 + (1-t)^2 ( \mathbf{P}_0 - \mathbf{P}_1) + t^2 (\mathbf{P}_2 - \mathbf{P}_1) \\[.5em]
+                  &= \mathbf{P}_1 + u^2 ( \mathbf{P}_0 - \mathbf{P}_1) + t^2 (\mathbf{P}_2 - \mathbf{P}_1)
+    \end{align}
 
 Which immediately gives the derivative of the Bézier curve with respect to `t`:
 
 .. math::
-    \mathbf{B}'(t) = 2(1 - t) (\mathbf{P}_1 - \mathbf{P}_0) + 2t (\mathbf{P}_2 - \mathbf{P}_1)
+    \begin{align}
+    \mathbf{B}'(t) &= 2(1-t) (\mathbf{P}_1 - \mathbf{P}_0) + 2t (\mathbf{P}_2 - \mathbf{P}_1) \\[.5em]
+                   &= 2 (\mathbf{P}_1 - \mathbf{P}_0) + 2(\mathbf{P}_0 - 2\mathbf{P}_1 + \mathbf{P}_2) t \\[.5em]
+                   &= 2u (\mathbf{P}_1 - \mathbf{P}_0) + 2t (\mathbf{P}_2 - \mathbf{P}_1)
+    \end{align}
 
 from which it can be concluded that the tangents to the curve at :math:`\mathbf{P}_0` and
 :math:`\mathbf{P}_2` intersect at :math:`\mathbf{P}_1`. As :math:`t` increases from 0 to 1, the
@@ -75,6 +94,8 @@ The second derivative of the Bézier curve with respect to :math:`t` is
 .. math::
     \mathbf{B}''(t) = 2 (\mathbf{P}_2 - 2 \mathbf{P}_1 + \mathbf{P}_0)
 
+
+
 Cubic Bézier Curves
 -------------------
 
@@ -99,7 +120,8 @@ The explicit form of the curve is:
 
 .. math::
     \begin{align}
-    \mathbf{B}(t) &= (1-t)^3 \mathbf{P}_0 + 3(1-t)^2t \mathbf{P}_1 + 3(1-t)t^2 \mathbf{P}_2 + t^3\mathbf{P}_3 \\
+    \mathbf{B}(t) &= (1-t)^3 \mathbf{P}_0 + 3(1-t)^2t \mathbf{P}_1 + 3(1-t)t^2 \mathbf{P}_2 + t^3\mathbf{P}_3 \\[.5em]
+                  &= u^3 \mathbf{P}_0 + 3u^2t \mathbf{P}_1 + 3ut^2 \mathbf{P}_2 + t^3\mathbf{P}_3 \\[.5em]
                   &= (\mathbf{P}_3 - 3\mathbf{P}_2 + 3\mathbf{P}_1 - \mathbf{P}_0) t^3 +
                      3(\mathbf{P}_2 - 2\mathbf{P}_1 + \mathbf{P}_0) t^2 +
                      3(\mathbf{P}_1 - \mathbf{P}_0) t  +
@@ -112,12 +134,20 @@ contain a cusp.
 The derivative of the cubic Bézier curve with respect to :math:`t` is
 
 .. math::
-    \mathbf{B}'(t) = 3(1-t)^2 (\mathbf{P}_1 - \mathbf{P}_0) + 6(1-t)t (\mathbf{P}_2 - \mathbf{P}_1) + 3t^2 (\mathbf{P}_3 - \mathbf{P}_2)
+    \begin{align}
+    \mathbf{B}'(t) &= 3(1-t)^2 (\mathbf{P}_1 - \mathbf{P}_0) + 6(1-t)t (\mathbf{P}_2 - \mathbf{P}_1) + 3t^2 (\mathbf{P}_3 - \mathbf{P}_2) \\[.5em]
+	           &= 3u^2 (\mathbf{P}_1 - \mathbf{P}_0) + 6ut (\mathbf{P}_2 - \mathbf{P}_1) + 3t^2 (\mathbf{P}_3 - \mathbf{P}_2)
+    \end{align}
 
 The second derivative of the Bézier curve with respect to :math:`t` is
 
 .. math::
-    \mathbf{B}''(t) = 6(1-t) (\mathbf{P}_2 - 2 \mathbf{P}_1 + \mathbf{P}_0) +  6t (\mathbf{P}_3 - 2 \mathbf{P}_2 + \mathbf{P}_1)
+    \begin{align}
+    \mathbf{B}''(t) &= 6(1-t) (\mathbf{P}_2 - 2 \mathbf{P}_1 + \mathbf{P}_0) +  6t (\mathbf{P}_3 - 2 \mathbf{P}_2 + \mathbf{P}_1) \\[.5em]
+                    &= 6u (\mathbf{P}_2 - 2 \mathbf{P}_1 + \mathbf{P}_0) +  6t (\mathbf{P}_3 - 2 \mathbf{P}_2 + \mathbf{P}_1)
+    \end{align}
+
+
 
 Recursive definition
 --------------------
@@ -134,7 +164,7 @@ The recursive definition is
 .. math::
     \begin{align}
     \mathbf{B}_{\mathbf{P}_0}(t) &= \mathbf{P}_0 \\[1em]
-    \mathbf{B}(t) &= \mathbf{B}_{\mathbf{P}_0\mathbf{P}_1\ldots\mathbf{P}_n}(t) \\
+    \mathbf{B}(t) &= \mathbf{B}_{\mathbf{P}_0\mathbf{P}_1\ldots\mathbf{P}_n}(t) \\[.5em]
                   &= (1-t) \mathbf{B}_{\mathbf{P}_0\mathbf{P}_1\ldots\mathbf{P}_{n-1}}(t) +
                          t \mathbf{B}_{\mathbf{P}_1\mathbf{P}_2\ldots\mathbf{P}_n}(t)
     \end{align}
@@ -143,18 +173,22 @@ The formula can be expressed explicitly as follows:
 
 .. math::
     \begin{align}
-    \mathbf{B}(t) &= \sum_{i=0}^n b_{i,n}(t) \mathbf{P}_i \\
-                  &= \sum_{i=0}^n {n\choose i}(1 - t)^{n - i}t^i \mathbf{P}_i \\
-                  &= (1 - t)^n \mathbf{P}_0 +
-                     {n\choose 1}(1 - t)^{n - 1}t \mathbf{P}_1 +
+    \mathbf{B}(t) &= \sum_{i=0}^n b_{i,n}(t) \mathbf{P}_i \\[.5em]
+                  &= \sum_{i=0}^n {n\choose i}(1-t)^{n - i}t^i \mathbf{P}_i \\[.5em]
+                  &= (1-t)^n \mathbf{P}_0 +
+                     {n\choose 1}(1-t)^{n - 1}t \mathbf{P}_1 +
                      \cdots +
-                     {n\choose n - 1}(1 - t)t^{n - 1} \mathbf{P}_{n - 1} +
+                     {n\choose n - 1}(1-t)t^{n - 1} \mathbf{P}_{n - 1} +
                      t^n \mathbf{P}_n
     \end{align}
 
 where :math:`b_{i,n}(t)` are the Bernstein basis polynomials of degree :math:`n` and :math:`n
 \choose i` are the binomial coefficients.
 
+
+
+.. _bezier-curve-degree-elevation-section:
+
 Degree elevation
 ----------------
 
@@ -173,13 +207,13 @@ For arbitrary :math:`n`, we have
 
 .. math::
     \begin{align}
-    \mathbf{B}(t) &= (1 - t) \sum_{i=0}^n \mathbf{b}_{i,n}(t) \mathbf{P}_i +
-                           t \sum_{i=0}^n \mathbf{b}_{i,n}(t) \mathbf{P}_i \\
+    \mathbf{B}(t) &= (1-t) \sum_{i=0}^n \mathbf{b}_{i,n}(t) \mathbf{P}_i +
+                           t \sum_{i=0}^n \mathbf{b}_{i,n}(t) \mathbf{P}_i \\[.5em]
                   &= \sum_{i=0}^n \frac{n + 1 - i}{n + 1} \mathbf{b}_{i, n + 1}(t) \mathbf{P}_i +
-                     \sum_{i=0}^n \frac{i + 1}{n + 1} \mathbf{b}_{i + 1, n + 1}(t) \mathbf{P}_i \\
+                     \sum_{i=0}^n \frac{i + 1}{n + 1} \mathbf{b}_{i + 1, n + 1}(t) \mathbf{P}_i \\[.5em]
                   &= \sum_{i=0}^{n + 1} \mathbf{b}_{i, n + 1}(t)
                                         \left(\frac{i}{n + 1}         \mathbf{P}_{i - 1} +
-                                              \frac{n + 1 - i}{n + 1} \mathbf{P}_i\right) \\
+                                              \frac{n + 1 - i}{n + 1} \mathbf{P}_i\right) \\[.5em]
                   &= \sum_{i=0}^{n + 1} \mathbf{b}_{i, n + 1}(t) \mathbf{P'}_i
     \end{align}
 
@@ -191,6 +225,18 @@ Therefore the new control points are
 It introduces two arbitrary points :math:`\mathbf{P}_{-1}` and :math:`\mathbf{P}_{n+1}` which are
 cancelled in :math:`\mathbf{P'}_i`.
 
+Example for a quadratic Bézier curve:
+
+.. math::
+   \begin{align}
+   \mathbf{P'}_0 &= \mathbf{P}_0 \\[.5em]
+   \mathbf{P'}_1 &= \mathbf{P}_0 + \frac{2}{3} (\mathbf{P}_1 - \mathbf{P}_0) \\[.5em]
+   \mathbf{P'}_1 &= \mathbf{P}_2 + \frac{2}{3} (\mathbf{P}_1 - \mathbf{P}_2) \\[.5em]
+   \mathbf{P'}_2 &= \mathbf{P}_2
+   \end{align}
+
+
+
 Matrix Forms
 ------------
 
@@ -240,6 +286,8 @@ Matrix Forms
           P1 (    2t - t^2) +
           P2           t^2
 
+
+
 Symbolic Calculation
 --------------------
 
@@ -256,9 +304,18 @@ Symbolic Calculation
     >>> B2_012 = (1-t)*((1-t)*P0 + t*P1) + t*((1-t)*P1 + t*P2)
     >>> B2_123 = (1-t)*((1-t)*P1 + t*P2) + t*((1-t)*P2 + t*P3)
     >>> B3 = (1-t)*B2_012 + t*B2_123
-    >>> collect(expand(B2), t)
+    >>> collect(expand(B3), t)
     P0 + t**3*(-P0 + 3*P1 - 3*P2 + P3) + t**2*(3*P0 - 6*P1 + 3*P2) + t*(-3*P0 + 3*P1)
 
+    # Compute derivative
+    >>> B2p = collect(simplify(B2.diff(t)), t)
+    -2*P0 + 2*P1 + t*(2*P0 - 4*P1 + 2*P2)
+
+    >>> B3p = collect(simplify(B3.diff(t)), t)
+    -3*P0 + 3*P1 + t**2*(-3*P0 + 9*P1 - 9*P2 + 3*P3) + t*(6*P0 - 12*P1 + 6*P2)
+
+
+
 .. _bezier-curve-length-section:
 
 Curve Length
@@ -302,6 +359,8 @@ For the case :math:`c < 0`, an antiderivative is
 
 where :math:`k = \frac{4c}{q}` with :math:`q = 4ac - b^2`.
 
+
+
 .. _bezier-curve-flatness-section:
 
 Determine the curve flatness
@@ -330,37 +389,150 @@ 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
+    U &= 3\mathbf{P}_1 - 2\mathbf{P}_0 - \mathbf{P}_3 \\[.5em]
+    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) \\
+    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) \\[.5em]
          &= (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
+    d(t)^2 = (1-t)^2 t^2 (((1-t) U_x + t V_x)^2 + ((1-t) U_y + t V_y)^2
 
 From
 
 .. math::
     \begin{align}
-    argmax((1 - t)^2 t^2)   &= \frac{1}{16} \\
-    argmax((1 - t) a + t b) &= argmax(a, b)
+    argmax((1-t)^2 t^2)   &= \frac{1}{16} \\[.5em]
+    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))
+    \mathrm{flatness}^2 = argmax(d(t)^2) \leq \frac{1}{16} (argmax(U_x^2, V_x^2) + argmax(U_y^2, V_y^2))
 
 Thus an upper bound of :math:`16\,\mathrm{flatness}^2` is
 
 .. math::
-    argmax(ux^2, vx^2) + argmax(uy^2, vy^2)
+    argmax(U_x^2, V_x^2) + argmax(U_y^2, V_y^2)
+
+
+
+.. _bezier-curve-line-intersection-section:
+
+Intersection of Bézier Curve with a Line
+----------------------------------------
+
+Algorithm
+
+  * Apply a transformation to the curve that maps the line onto the X-axis.
+  * Then we only need to test the Y-values for a zero.
+
+
+
+.. _bezier-curve-closest-point-section:
+
+Closest Point
+-------------
+
+Reference
+
+  * https://hal.archives-ouvertes.fr/inria-00518379/document
+    Improved Algebraic Algorithm On Point Projection For Bézier Curves
+    Xiao-Diao Chen, Yin Zhou, Zhenyu Shu, Hua Su, Jean-Claude Paul
+
+Let a point :math:`\mathbf{P}` and the closest point :math:`\mathbf{B}(t)` on the curve, we have this
+condition:
+
+.. math::
+    (\mathbf{P} - \mathbf{B}(t)) \cdot \mathbf{B}'(t) = 0 \\[.5em]
+    \mathbf{P} \cdot \mathbf{B}'(t) - \mathbf{B}(t) \cdot \mathbf{B}'(t) = 0
+
+Quadratic Bézier Curve
+~~~~~~~~~~~~~~~~~~~~~~
+
+Let
+
+.. math::
+    \begin{align}
+       \mathbf{A} &= \mathbf{P}_1 - \mathbf{P}_0 \\[.5em]
+       \mathbf{B} &= \mathbf{P}_2 - \mathbf{P}_1 -\mathbf{A} \\[.5em]
+       \mathbf{M} &= \mathbf{P}_0 - \mathbf{P}
+    \end{align}
+
+We have
+
+.. math::
+    \mathbf{B}'(t) = 2(\mathbf{A} + \mathbf{B} t)
+
+.. factorisation
+   (P0 - 2*P1 + P2)**2 * t**3
+   3*(P1 - P0)*(P0 - 2*P1 + P2) * t**2
+   ...
+   (P0 - P)*(P1 - P0)
+
+The condition can be expressed as
+
+.. math::
+    \mathbf{B}^2 t^3 + 3\mathbf{A}\mathbf{B} t^2 + (2\mathbf{A}^2 + \mathbf{M}\mathbf{B}) t + \mathbf{M}\mathbf{A} = 0
+
+.. code-block:: py3
+
+    >>> C = collect(expand((P*B2p - B2*B2p)/-2), t)
+    P*P0 - P*P1 - P0**2 + P0*P1 +
+    t**3 * (P0**2 - 4*P0*P1 + 2*P0*P2 + 4*P1**2 - 4*P1*P2 + P2**2) +
+    t**2 * (-3*P0**2 + 9*P0*P1 - 3*P0*P2 - 6*P1**2 + 3*P1*P2) +
+    t    * (-P*P0 + 2*P*P1 - P*P2 + 3*P0**2 - 6*P0*P1 + P0*P2 + 2*P1**2)
+    >>> A = P1 - P0
+        B = P2 - P1 - A
+        M = P0 - P
+    >>> C2 = B**2 * t**3 + 3*A*B * t**2 + (2*A**2 + M*B) * t + M*A
+    >>> expand(C - C2)
+    0
+
+Cubic Bézier Curve
+~~~~~~~~~~~~~~~~~~
+
+Let
+
+.. math::
+    \begin{align}
+     n &= \mathbf{P}_3 - 3\mathbf{P}_2 + 3\mathbf{P}_1 - \mathbf{P}_0 \\[.5em]
+     r &= 3(\mathbf{P}_2 - 2\mathbf{P}_1 + \mathbf{P}_0 \\[.5em]
+     s &= 3(\mathbf{P}_1 - \mathbf{P}_0) \\[.5em]
+     v &= \mathbf{P}_0
+     \end{align}
+
+We have
+
+.. math::
+    \begin{align}
+    \mathbf{B}^3(t)  &= nt^3 + rt^2 + st + v \\[.5em]
+    \mathbf{B}'(t) &= 3nt^2 + 2rt + s
+    \end{align}
+
+    .. n = P3 - 3*P2 + 3*P1 - P0
+   r = 3*(P2 - 2*P1 + P0
+   s = 3*(P1 - P0)
+   v = P0
+
+.. code-block:: py3
+
+    >>> n, r, s, v = symbols('n r s v')
+    >>> B3 = n*t**3 + r*t**2 + s*t + v
+    >>> B3p = simplify(B3.diff(t))
+    >>> C = collect(expand((P*B3p - B3*B3p)), 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
diff --git a/doc/sphinx/source/resources/geometry/spline.rst b/doc/sphinx/source/resources/geometry/spline.rst
index edb33ff908aa287ef9b0ab388e2e14f958f88983..bd8afa34003bd4489819011bc77222bdd84e2f33 100644
--- a/doc/sphinx/source/resources/geometry/spline.rst
+++ b/doc/sphinx/source/resources/geometry/spline.rst
@@ -8,6 +8,23 @@
 
 .. contents:: :local:
 
+
+..  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.
+
+References
+----------
+
+* Computer Graphics, Principle and Practice, Foley et al., Adison Wesley
+* http://web.mit.edu/hyperbook/Patrikalakis-Maekawa-Cho/node15.html
+
 B-spline Basis
 --------------
 
@@ -240,9 +257,3 @@ will become :math:`C^{k-p-1}`, where `p` is the multiplicity of the knot.
 
 The B-spline curve can be subdivided into Bézier segments by knot insertion at each internal knot
 until the multiplicity of each internal knot is equal to `k`.
-
-References
-----------
-
-* Computer Graphics, Principle and Practice, Foley et al., Adison Wesley
-* http://web.mit.edu/hyperbook/Patrikalakis-Maekawa-Cho/node15.html