diff --git a/CHANGELOG.md b/CHANGELOG.md
index 825fba2888..f0455ed43c 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -10,6 +10,10 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/).
 
 ## [Unreleased]
 
+### Added
+- Add spline joint (JointModelSplineTpl) to default joint collection ([#2784](https://github.com/stack-of-tasks/pinocchio/pull/2784))
+  - [Python example here](./examples/spline-joint.py)
+
 ## [4.1.0] - 2026-07-07
 
 ### Added
diff --git a/bindings/python/parsers/graph/expose-edges.cpp b/bindings/python/parsers/graph/expose-edges.cpp
index 91466af7ab..883c9998c4 100644
--- a/bindings/python/parsers/graph/expose-edges.cpp
+++ b/bindings/python/parsers/graph/expose-edges.cpp
@@ -53,6 +53,9 @@ namespace pinocchio
 
     const int JointMimic::nq;
     const int JointMimic::nv;
+
+    const int JointSpline::nq;
+    const int JointSpline::nv;
   } // namespace graph
 
   namespace python
@@ -160,6 +163,49 @@ namespace pinocchio
         .def_readonly("nq", &JointUniversal::nq, "Number of configuration variables.")
         .def_readonly("nv", &JointUniversal::nv, "Number of tangent variables.");
 
+      bp::class_<JointSpline>(
+        "JointSpline", "Represents a spline-based joint.",
+        bp::init<>(bp::args("self"), "Default constructor."))
+        .def(
+          bp::init<const std::vector<SE3> &, const Eigen::VectorXd &, std::size_t>(
+            bp::args("self", "ctrlFrame", "knots", "degree"),
+            "Constructor with a single control frame and degree."))
+        .def_readwrite("ctrlFrames", &JointSpline::ctrlFrames, "Control frames of the spline.")
+        .def_readwrite("knots", &JointSpline::knots, "Control frames of the spline.")
+        .def_readwrite("degree", &JointSpline::degree, "Degree of the spline.")
+        .def_readonly("nq", &JointSpline::nq, "Number of configuration variables.")
+        .def_readonly("nv", &JointSpline::nv, "Number of tangent variables.");
+
+      bp::class_<JointSplineBuilder>(
+        "JointSplineBuilder", "JointSpline builder helper.",
+        bp::init<>(bp::args("self"), "Default constructor."))
+        .def(
+          "addControlFrame", &JointSplineBuilder::addControlFrame, bp::return_self<>(),
+          (bp::arg("self"), bp::arg("frame")), "Add a B-spline control frame")
+        .def(
+          "withControlFrameVector", &JointSplineBuilder::withControlFrameVector,
+          bp::return_self<>(), (bp::arg("self"), bp::arg("frames")),
+          "Set B-spline control frame vector")
+        .def(
+          "withDegree", &JointSplineBuilder::withDegree, bp::return_self<>(),
+          (bp::arg("self"), bp::arg("degree")), "Set B-spline degree")
+        .def(
+          "withKnotVector",
+          +[](JointSplineBuilder & builder, const std::vector<double> & k) -> auto {
+            return builder.withKnotVector(k);
+          },
+          bp::return_self<>(), (bp::arg("self"), bp::arg("knots")), "Set B-spline knot vector")
+        .def(
+          "withOpenUniformKnots", &JointSplineBuilder::withOpenUniformKnots, bp::return_self<>(),
+          (bp::arg("self"), bp::arg("min"), bp::arg("max")),
+          "Set B-spline knot vector as open uniform")
+        .def(
+          "withUniformKnots", &JointSplineBuilder::withUniformKnots, bp::return_self<>(),
+          (bp::arg("self"), bp::arg("min"), bp::arg("max")), "Set B-spline knot vector as uniform")
+        .def(
+          "build", &JointSplineBuilder::build, (bp::arg("self")),
+          "Build a JointSpline from provided parameters");
+
       bp::class_<JointComposite>(
         "JointComposite", "Represents a composite joint.",
         bp::init<>(bp::args("self"), "Default constructor."))
diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt
index 0b3421200b..b123c78da2 100644
--- a/examples/CMakeLists.txt
+++ b/examples/CMakeLists.txt
@@ -155,6 +155,8 @@ if(BUILD_PYTHON_INTERFACE)
             simulation-closed-kinematic-chains
             simulation-lcaba
             model-graph-geometry
+            spline-joint
+            spline-joint-knee
         )
         if(BUILD_WITH_URDF_SUPPORT)
             list(
diff --git a/examples/spline-joint-knee.py b/examples/spline-joint-knee.py
new file mode 100644
index 0000000000..e82efc67dc
--- /dev/null
+++ b/examples/spline-joint-knee.py
@@ -0,0 +1,146 @@
+"""Knee spline joint driven by its physical flexion angle.
+
+This example uses the SE(3) spline joint (``pin.JointModelSpline`` /
+``pin.JointModelSplineBuilder``) to reproduce experimental knee kinematics: the
+joint constrains its child (tibia) frame to a cumulative B-spline on SE(3)
+defined by a list of control frames, a knot vector and a degree
+
+source: Lee, S. H., & Terzopoulos, D. (2008).
+Spline joints for multibody dynamics. In ACM SIGGRAPH 2008 papers (pp. 1-8).
+
+The 8 control frames below were fit offline (Gauss-Newton, degree=3, n_ctrl=8)
+to a planar OpenSim knee, where the tx/ty translations the rotation
+about Z (knee flexion) are driven by knee flexion itself.
+"""
+
+import time
+from pathlib import Path
+
+import coal
+import meshcat.geometry as mg
+import numpy as np
+import pinocchio as pin
+from pinocchio.visualize import MeshcatVisualizer
+
+# Bones: mesh folder, OpenSim (subject-specific) scale factors and display color.
+MESH_DIR = Path(__file__).resolve().parent.parent / "models" / "biomechanics"
+FEMUR_SCALE, TIBIA_SCALE = 1.17378, 1.12373
+BONE_COLOR = np.array([0.96, 0.94, 0.86, 0.7])  # bone ivory, semi-transparent.
+
+# Recenter so the femoral condyle (knee) sits at the world origin instead of the
+# hip: the fixed femur and the spline joint are both lifted by the knee offset.
+KNEE_FROM_HIP = np.array([-0.0046, -0.4644, 0.0])
+RECENTER = pin.SE3(np.eye(3), -KNEE_FROM_HIP)
+
+# Control frames from the offline 2D-knee regression, as (yaw [rad], tx, ty [m]).
+CONTROL_FRAMES = [
+    (-2.094400, -0.003842, -0.496205),
+    (-1.815147, 0.001278, -0.489246),
+    (-1.256640, 0.006243, -0.480186),
+    (-0.418880, -0.002220, -0.466444),
+    (0.418880, -0.008427, -0.463391),
+    (1.256640, -0.005740, -0.461613),
+    (1.815147, -0.006127, -0.464906),
+    (2.094400, -0.006285, -0.464416),
+]
+
+
+def mesh_object(name, parent_joint, stl, scale, placement):
+    """Build a mesh GeometryObject fixed to ``parent_joint`` at ``placement``."""
+    mesh_scale = scale * np.ones(3)
+    geom = coal.MeshLoader().load(str(stl), mesh_scale)
+    obj = pin.GeometryObject(name, parent_joint, placement, geom, str(stl), mesh_scale)
+    obj.meshColor = BONE_COLOR
+    return obj
+
+
+class PointTracer:
+    """Trace a fixed point in Meshcat, keeping the dotted trail."""
+
+    def __init__(self, viz, name, point=None, color=0xFF3030, size=0.004):
+        self._world_to_femur = RECENTER.inverse()
+        self._point = np.zeros(3) if point is None else np.asarray(point, float)
+        self._material = mg.PointsMaterial(size=size, color=color)
+        self._trail = []
+        self._node = viz.viewer[f"tibia_traces/{name}"]
+        self._node.set_transform(RECENTER.homogeneous)
+
+    def add(self, tibia_pose):
+        self._trail.append(self._world_to_femur.act(tibia_pose.act(self._point)))
+        pts = np.asarray(self._trail, dtype=np.float32).T
+        self._node.set_object(mg.Points(mg.PointsGeometry(pts), self._material))
+
+
+def main():
+
+    control_frames = [
+        pin.SE3(pin.rpy.rpyToMatrix(0.0, 0.0, yaw), np.array([tx, ty, 0.0]))
+        for (yaw, tx, ty) in CONTROL_FRAMES
+    ]
+
+    # Knee-angle domain of the fitted spline (in degree).
+    THETA_MIN_DEG, THETA_MAX_DEG = -120.0, 120.0
+
+    knee = (
+        pin.JointModelSplineBuilder()
+        .withDegree(
+            3
+        )  # degree of the splines, degree 3 means acceleration is continuous.
+        .withControlFrameVector(control_frames)  # list of SE3 transforms
+        .withOpenUniformKnots(
+            THETA_MIN_DEG, THETA_MAX_DEG
+        )  # Map knee flexion angle (degrees) to spline parameter q in [0, 1].
+        .build()
+    )
+    model = pin.Model()
+    joint_id = model.addJoint(0, knee, RECENTER, "knee")
+
+    # visual only code
+    try:
+        visual_model = pin.GeometryModel()
+        visual_model.addGeometryObject(
+            mesh_object("femur", 0, MESH_DIR / "femur_r.stl", FEMUR_SCALE, RECENTER)
+        )
+        visual_model.addGeometryObject(
+            mesh_object(
+                "tibia",
+                joint_id,
+                MESH_DIR / "tibia_r.stl",
+                TIBIA_SCALE,
+                pin.SE3.Identity(),
+            )
+        )
+
+        viz = MeshcatVisualizer(model, None, visual_model)
+        viz.initViewer(open=True)
+        viz.loadViewerModel()
+
+    except ImportError as e:
+        print("Error while initializing the viewer.")
+        print(e)
+        return
+
+    # Trace tibia-frame points and show the tibia frame as a triad, updated every step.
+    plateau = np.array([0.03100708, -0.03819558, -0.00347008])
+    tracers = [
+        PointTracer(viz, "origin", color=0xFF3030),
+        PointTracer(viz, "plateau", plateau, color=0x3060FF),
+    ]
+    tibia_frame = viz.viewer["tibia_frame"]
+    tibia_frame.set_object(mg.triad(0.08))
+
+    # Sweep the physical knee angle (+10 deg hyperextension to -120 deg flexion),
+    angles = np.linspace(10.0, -120.0, 120)
+    angles = np.concatenate([angles, angles[::-1]])
+    for _ in range(5):
+        for theta_deg in angles:
+            viz.display(np.array([theta_deg]))
+            tibia_pose = viz.data.oMi[joint_id]
+            for tracer in tracers:
+                tracer.add(tibia_pose)
+            tibia_frame.set_transform(tibia_pose.homogeneous)
+            time.sleep(0.015)
+
+
+if __name__ == "__main__":
+    main()
diff --git a/examples/spline-joint.py b/examples/spline-joint.py
new file mode 100644
index 0000000000..7a832e1d0c
--- /dev/null
+++ b/examples/spline-joint.py
@@ -0,0 +1,111 @@
+import time
+
+import meshcat.geometry as mg
+import numpy as np
+import pinocchio as pin
+from pinocchio.visualize import MeshcatVisualizer
+
+
+class PointTracer:
+    """Trace a fixed point in Meshcat, keeping the dotted trail."""
+
+    def __init__(self, viz):
+        self._material = mg.PointsMaterial(size=0.004, color=0xFF3030)
+        self._trail = []
+        self._node = viz.viewer["solid_pose"]
+
+    def add(self, solid_pose):
+        self._trail.append(solid_pose.translation.copy())
+        pts = np.asarray(self._trail, dtype=np.float32).T
+
+        line = mg.Line(
+            geometry=mg.PointsGeometry(pts),
+            material=mg.LineBasicMaterial(color=0xFF3030),
+        )
+        self._node.set_object(geometry=line)
+
+
+def main():
+
+    PARABOLA_CONTROL_FRAMES = [
+        (-1.000000, 1.000000),
+        (-0.888889, 0.777778),
+        (-0.666667, 0.407407),
+        (-0.333333, 0.074074),
+        (0.000000, -0.037037),
+        (0.333333, 0.074074),
+        (0.666667, 0.407407),
+        (0.888889, 0.777778),
+        (1.000000, 1.000000),
+    ]
+
+    control_frames = [
+        pin.SE3(np.eye(3), np.array([0, ty, tz]))
+        for (ty, tz) in PARABOLA_CONTROL_FRAMES
+    ]
+
+    # Create a Pinocchio model with a single free-flyer joint
+    model = pin.Model()
+    spline_joint = (
+        pin.JointModelSplineBuilder()
+        .withDegree(3)
+        .withControlFrameVector(control_frames)
+        .withOpenUniformKnots(-1, 1)
+        .build()
+    )
+    joint_id = model.addJoint(0, spline_joint, pin.SE3.Identity(), "spline-joint")
+
+    # adding some inertia to drop the solid in the parabola
+    box_size = 1.0  # edge length of each cube [m]
+    box_mass = 1.0  # mass of each cube [kg]
+    box_inertia = pin.Inertia.FromBox(box_mass, box_size, box_size, box_size)
+    model.appendBodyToJoint(joint_id, box_inertia, pin.SE3.Identity())
+
+    # visual only code
+    try:
+        visual_model = pin.GeometryModel()
+
+        viz = MeshcatVisualizer(model, None, visual_model)
+        viz.initViewer(open=True)
+        viz.loadViewerModel()
+
+    except ImportError as e:
+        print("Error while initializing the viewer.")
+        print(e)
+        return
+
+    tracers = [PointTracer(viz)]
+    solid_frame = viz.viewer["solid_frame"]
+    solid_frame.set_object(mg.triad(0.2))
+
+    dt = 0.01
+    qs, _ = sim_loop(model, dt)
+
+    for theta in qs:
+        viz.display(np.array([theta]))
+        solid_pose = viz.data.oMi[joint_id]
+        for tracer in tracers:
+            tracer.add(solid_pose)
+        solid_frame.set_transform(solid_pose.homogeneous)
+        time.sleep(dt)
+
+
+def sim_loop(model, dt=0.01, nsteps=800):
+
+    qs = [np.array([1.0])]
+    vs = [np.array([0])]
+    data = model.createData()
+    for i in range(nsteps):
+        q = qs[i]
+        v = vs[i]
+        tau = -1 * v  # a little bit of damping
+        a1 = pin.aba(model, data, q, v, tau)
+        vnext = v + dt * a1
+        qnext = pin.integrate(model, q, dt * vnext)
+        qs.append(qnext)
+        vs.append(vnext)
+    return qs, vs
+
+
+if __name__ == "__main__":
+    main()
diff --git a/include/pinocchio/autodiff/casadi.hpp b/include/pinocchio/autodiff/casadi.hpp
index e26732d6a5..182bcd819a 100644
--- a/include/pinocchio/autodiff/casadi.hpp
+++ b/include/pinocchio/autodiff/casadi.hpp
@@ -25,6 +25,7 @@
 #include "pinocchio/eigen-common.hpp"
 #include "pinocchio/math.hpp"
 #include "pinocchio/spatial.hpp"
+#include "pinocchio/src/multibody/joint/spline-utils.hxx"
 // IWYU pragma: end_keep
 
 namespace boost
@@ -349,4 +350,5 @@ namespace Eigen
 #include "pinocchio/src/autodiff/casadi/math/triangular-matrix.hxx"
 #include "pinocchio/src/autodiff/casadi/spatial/se3-tpl.hxx"
 #include "pinocchio/src/autodiff/casadi/utils/static-if.hxx"
+#include "pinocchio/src/autodiff/casadi/multibody/joint/spline-utils.hxx"
 // IWYU pragma: end_exports
diff --git a/include/pinocchio/bindings/python/context/generic.hpp b/include/pinocchio/bindings/python/context/generic.hpp
index 3426d99081..3897b1ef32 100644
--- a/include/pinocchio/bindings/python/context/generic.hpp
+++ b/include/pinocchio/bindings/python/context/generic.hpp
@@ -135,6 +135,10 @@ namespace pinocchio
       typedef JointModelUniversalTpl<Scalar> JointModelUniversal;
       typedef JointDataUniversalTpl<Scalar> JointDataUniversal;
 
+      typedef JointModelSplineTpl<Scalar> JointModelSpline;
+      typedef JointModelSplineBuilderTpl<Scalar> JointModelSplineBuilder;
+      typedef JointDataSplineTpl<Scalar> JointDataSpline;
+
       typedef JointModelTranslationTpl<Scalar> JointModelTranslation;
       typedef JointDataTranslationTpl<Scalar> JointDataTranslation;
 
diff --git a/include/pinocchio/bindings/python/multibody/joint/joints-datas.hpp b/include/pinocchio/bindings/python/multibody/joint/joints-datas.hpp
index 0aebd60835..fd505b71c7 100644
--- a/include/pinocchio/bindings/python/multibody/joint/joints-datas.hpp
+++ b/include/pinocchio/bindings/python/multibody/joint/joints-datas.hpp
@@ -83,5 +83,11 @@ namespace pinocchio
         .add_property("StU", &JointDataComposite::StU);
     }
 
+    template<>
+    inline bp::class_<JointDataSpline> &
+    expose_joint_data<JointDataSpline>(bp::class_<JointDataSpline> & cl)
+    {
+      return cl.add_property("N", &JointDataSpline::N);
+    }
   } // namespace python
 } // namespace pinocchio
diff --git a/include/pinocchio/bindings/python/multibody/joint/joints-models.hpp b/include/pinocchio/bindings/python/multibody/joint/joints-models.hpp
index 4f957968ec..f5a4649a96 100644
--- a/include/pinocchio/bindings/python/multibody/joint/joints-models.hpp
+++ b/include/pinocchio/bindings/python/multibody/joint/joints-models.hpp
@@ -311,6 +311,58 @@ namespace pinocchio
           "Second rotation axis of the JointModelUniversal.");
     }
 
+    // Specialization for JointModelSpline
+    template<>
+    bp::class_<context::JointModelSpline> &
+    expose_joint_model<context::JointModelSpline>(bp::class_<context::JointModelSpline> & cl)
+    {
+      bp::class_<context::JointModelSplineBuilder>(
+        "JointModelSplineBuilder", "JointSpline builder helper.",
+        bp::init<>(bp::args("self"), "Default constructor."))
+        .def(
+          "addControlFrame", &context::JointModelSplineBuilder::addControlFrame,
+          bp::return_self<>(), (bp::arg("self"), bp::arg("frame")), "Add a B-spline control frame")
+        .def(
+          "withControlFrameVector", &context::JointModelSplineBuilder::withControlFrameVector,
+          bp::return_self<>(), (bp::arg("self"), bp::arg("frames")),
+          "Set B-spline control frame vector")
+        .def(
+          "withDegree", &context::JointModelSplineBuilder::withDegree, bp::return_self<>(),
+          (bp::arg("self"), bp::arg("degree")), "Set B-spline degree")
+        .def(
+          "withKnotVector",
+          +[](context::JointModelSplineBuilder & builder, const std::vector<context::Scalar> & k)
+            -> auto { return builder.withKnotVector(k); },
+          bp::return_self<>(), (bp::arg("self"), bp::arg("knots")), "Set B-spline knot vector")
+        .def(
+          "withOpenUniformKnots", &context::JointModelSplineBuilder::withOpenUniformKnots,
+          bp::return_self<>(), (bp::arg("self"), bp::arg("min"), bp::arg("max")),
+          "Set B-spline knot vector as open uniform")
+        .def(
+          "withUniformKnots", &context::JointModelSplineBuilder::withUniformKnots,
+          bp::return_self<>(), (bp::arg("self"), bp::arg("min"), bp::arg("max")),
+          "Set B-spline knot vector as uniform")
+        .def(
+          "build", &context::JointModelSplineBuilder::build, (bp::arg("self")),
+          "Build a JointSpline from provided parameters");
+
+      return cl
+        .def(
+          bp::init<>(
+            bp::args("self"),
+            "Init an empty joint Spline. Default degree of spline basis function is 3."))
+        .def(
+          bp::init<const std::vector<context::SE3> &, context::VectorXs &, int>(
+            bp::args("self", "controlFrames", "knotVector", "degree"),
+            "Init a joint Spline, with a list of controlFrames, a knot vector and the degree of "
+            "the future "
+            "basis functions"))
+        .def_readwrite(
+          "degree", &context::JointModelSpline::degree, "Degree of the spline basis functions")
+        .def_readwrite("min_q", &context::JointModelSpline::min_q, "Minimum of the q entry")
+        .def_readwrite("max_q", &context::JointModelSpline::max_q, "Maximum of the q entry");
+    }
+
     // specialization for JointModelComposite
 
     struct JointModelCompositeAddJointVisitor
diff --git a/include/pinocchio/multibody/joint.hpp b/include/pinocchio/multibody/joint.hpp
index be28229c4e..142678418b 100644
--- a/include/pinocchio/multibody/joint.hpp
+++ b/include/pinocchio/multibody/joint.hpp
@@ -79,6 +79,9 @@
 
 #include "pinocchio/src/multibody/joint/joint-mimic.hxx"
 
+#include "pinocchio/src/multibody/joint/spline-utils.hxx"
+#include "pinocchio/src/multibody/joint/joint-spline.hxx"
+
 #include "pinocchio/src/multibody/joint/joint-collection.hxx"
 
 #include "pinocchio/src/multibody/joint/joint-generic.hxx"
diff --git a/include/pinocchio/src/autodiff/casadi/multibody/joint/spline-utils.hxx b/include/pinocchio/src/autodiff/casadi/multibody/joint/spline-utils.hxx
new file mode 100644
index 0000000000..fe0a559d40
--- /dev/null
+++ b/include/pinocchio/src/autodiff/casadi/multibody/joint/spline-utils.hxx
@@ -0,0 +1,50 @@
+//
+// Copyright (c) 2025 INRIA
+//
+
+#pragma once
+
+// IWYU pragma: private, include "pinocchio/autodiff/casadi.hpp"
+
+#ifdef PINOCCHIO_LSP
+  #undef PINOCCHIO_LSP
+  #include "pinocchio/autodiff/casadi.hpp"
+#endif // PINOCCHIO_LSP
+
+namespace pinocchio
+{
+  namespace internal
+  {
+    template<typename _Scalar, int Options>
+    struct SplineKinematics<::casadi::Matrix<_Scalar>, Options>
+    {
+      using Scalar = ::casadi::Matrix<_Scalar>;
+
+      static Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic>
+      allocateBasis(int degree, int knot_size)
+      {
+        return Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic>(
+          degree + 1, knot_size - degree - 1);
+      }
+
+      static void compute(
+        int degree,
+        const Eigen::Matrix<Scalar, Eigen::Dynamic, 1> & knots,
+        const std::vector<SE3Tpl<Scalar, Options>> & ctrlFrames,
+        const std::vector<MotionTpl<Scalar, Options>> & relativeMotions,
+        Scalar q,
+        const Eigen::Matrix<Scalar, 1, 1, Options> & joint_v,
+        bool computeVelocity,
+        SE3Tpl<Scalar, Options> & M,
+        MotionTpl<Scalar, Options> & v,
+        MotionTpl<Scalar, Options> & c,
+        JointMotionSubspaceTpl<1, Scalar, Options, 1> & S,
+        Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic> & basis)
+      {
+        return computeSplineKinematicsFull(
+          degree, knots, ctrlFrames, relativeMotions, q, joint_v, computeVelocity, M, v, c, S,
+          basis);
+      }
+    };
+  } // namespace internal
+} // namespace pinocchio
diff --git a/include/pinocchio/src/multibody/joint/fwd.hxx b/include/pinocchio/src/multibody/joint/fwd.hxx
index 3d53b15a9a..868dbeb500 100644
--- a/include/pinocchio/src/multibody/joint/fwd.hxx
+++ b/include/pinocchio/src/multibody/joint/fwd.hxx
@@ -159,6 +159,14 @@ namespace pinocchio
   struct JointDataTranslationTpl;
   typedef JointDataTranslationTpl<context::Scalar> JointDataTranslation;
 
+  template<typename Scalar, int Options = context::Options>
+  struct JointModelSplineTpl;
+  typedef JointModelSplineTpl<context::Scalar> JointModelSpline;
+
+  template<typename Scalar, int Options = context::Options>
+  struct JointDataSplineTpl;
+  typedef JointDataSplineTpl<context::Scalar> JointDataSpline;
+
   template<typename Scalar, int Options = context::Options>
   struct JointCollectionDefaultTpl;
   typedef JointCollectionDefaultTpl<context::Scalar> JointCollectionDefault;
@@ -205,4 +213,8 @@ namespace pinocchio
   struct JointDataTpl;
   typedef JointDataTpl<context::Scalar> JointData;
 
+  template<typename _Scalar, int _Options = context::Options>
+  struct JointModelSplineBuilderTpl;
+  typedef JointModelSplineBuilderTpl<context::Scalar> JointModelSplineBuilder;
+
 } // namespace pinocchio
diff --git a/include/pinocchio/src/multibody/joint/joint-collection.hxx b/include/pinocchio/src/multibody/joint/joint-collection.hxx
index 77e5c5400d..653dfa6e1f 100644
--- a/include/pinocchio/src/multibody/joint/joint-collection.hxx
+++ b/include/pinocchio/src/multibody/joint/joint-collection.hxx
@@ -83,6 +83,9 @@ namespace pinocchio
     // Joint Universal
     typedef JointModelUniversalTpl<Scalar, Options> JointModelUniversal;
 
+    // Joint Spline
+    typedef JointModelSplineTpl<Scalar, Options> JointModelSpline;
+
     typedef boost::variant<
       //    JointModelVoid,
       JointModelRX,
@@ -108,6 +111,7 @@ namespace pinocchio
       JointModelHz,
       JointModelHelicalUnaligned,
       JointModelUniversal,
+      JointModelSpline,
       boost::recursive_wrapper<JointModelComposite>,
       boost::recursive_wrapper<JointModelMimic>>
       JointModelVariant;
@@ -174,6 +178,9 @@ namespace pinocchio
     // Joint Universal
     typedef JointDataUniversalTpl<Scalar, Options> JointDataUniversal;
 
+    // Joint Spline
+    typedef JointDataSplineTpl<Scalar, Options> JointDataSpline;
+
     typedef boost::variant<
       //    JointDataVoid
       JointDataRX,
@@ -199,6 +206,7 @@ namespace pinocchio
       JointDataHz,
       JointDataHelicalUnaligned,
       JointDataUniversal,
+      JointDataSpline,
       boost::recursive_wrapper<JointDataComposite>,
       boost::recursive_wrapper<JointDataMimic>>
       JointDataVariant;
diff --git a/include/pinocchio/src/multibody/joint/joint-spline.hxx b/include/pinocchio/src/multibody/joint/joint-spline.hxx
new file mode 100644
index 0000000000..964db90f88
--- /dev/null
+++ b/include/pinocchio/src/multibody/joint/joint-spline.hxx
@@ -0,0 +1,504 @@
+//
+// Copyright (c) 2025 INRIA
+// Copyright (c) 2025-2026 ISIR
+//
+
+#pragma once
+
+// IWYU pragma: private, include "pinocchio/multibody/joint.hpp"
+
+#ifdef PINOCCHIO_LSP
+  #undef PINOCCHIO_LSP
+  #include "pinocchio/multibody/joint.hpp"
+#endif // PINOCCHIO_LSP
+
+namespace pinocchio
+{
+  template<typename Scalar, int Options>
+  struct JointSplineTpl;
+
+  template<typename _Scalar, int _Options>
+  struct traits<JointSplineTpl<_Scalar, _Options>>
+  {
+    enum
+    {
+      NQ = 1,
+      NV = 1,
+      NVExtended = 1
+    };
+    typedef _Scalar Scalar;
+    enum
+    {
+      Options = _Options
+    };
+    typedef JointDataSplineTpl<Scalar, Options> JointDataDerived;
+    typedef JointModelSplineTpl<Scalar, Options> JointModelDerived;
+
+    typedef JointMotionSubspaceTpl<1, Scalar, Options, 1> Constraint_t;
+    typedef SE3Tpl<Scalar, Options> Transformation_t;
+    typedef MotionTpl<Scalar, Options> Motion_t;
+    typedef MotionTpl<Scalar, Options> Bias_t;
+
+    // [ABA]
+    typedef Eigen::Matrix<Scalar, 6, NV, Options> U_t;
+    typedef Eigen::Matrix<Scalar, NV, NV, Options> D_t;
+    typedef Eigen::Matrix<Scalar, 6, NV, Options> UD_t;
+
+    typedef Eigen::Matrix<Scalar, NQ, 1, Options> ConfigVector_t;
+    typedef Eigen::Matrix<Scalar, NV, 1, Options> TangentVector_t;
+
+    typedef boost::mpl::false_ is_mimicable_t;
+
+    PINOCCHIO_JOINT_DATA_BASE_ACCESSOR_DEFAULT_RETURN_TYPE
+  };
+
+  template<typename _Scalar, int _Options>
+  struct traits<JointDataSplineTpl<_Scalar, _Options>>
+  {
+    typedef JointSplineTpl<_Scalar, _Options> JointDerived;
+    typedef _Scalar Scalar;
+  };
+
+  template<typename _Scalar, int _Options>
+  struct traits<JointModelSplineTpl<_Scalar, _Options>>
+  {
+    typedef JointSplineTpl<_Scalar, _Options> JointDerived;
+    typedef _Scalar Scalar;
+  };
+
+  template<typename _Scalar, int _Options>
+  struct JointDataSplineTpl : public JointDataBase<JointDataSplineTpl<_Scalar, _Options>>
+  {
+    EIGEN_MAKE_ALIGNED_OPERATOR_NEW
+    typedef JointSplineTpl<_Scalar, _Options> JointDerived;
+    typedef Eigen::Matrix<_Scalar, Eigen::Dynamic, Eigen::Dynamic> Matrix;
+
+    PINOCCHIO_JOINT_DATA_TYPEDEF_TEMPLATE(JointDerived);
+    PINOCCHIO_JOINT_DATA_BASE_DEFAULT_ACCESSOR
+
+    ConfigVector_t joint_q;
+    TangentVector_t joint_v;
+
+    Constraint_t S;
+    Transformation_t M;
+    Motion_t v;
+    Bias_t c;
+
+    // [ABA] specific data
+    U_t U;
+    D_t Dinv;
+    UD_t UDinv;
+    D_t StU;
+
+    // Bspline values
+    Matrix N;
+
+    JointDataSplineTpl()
+    : joint_q(ConfigVector_t::Zero())
+    , joint_v(TangentVector_t::Zero())
+    , M(Transformation_t::Identity())
+    , v(Motion_t::Zero())
+    , c(Bias_t::Zero())
+    , U(U_t::Zero())
+    , Dinv(D_t::Zero())
+    , UDinv(UD_t::Identity())
+    , StU(D_t::Zero())
+    , N(Matrix::Zero(1, 1))
+    {
+    }
+
+    JointDataSplineTpl(int degree, int knot_size)
+    : joint_q(ConfigVector_t::Zero())
+    , joint_v(TangentVector_t::Zero())
+    , M(Transformation_t::Identity())
+    , v(Motion_t::Zero())
+    , c(Bias_t::Zero())
+    , U(U_t::Zero())
+    , Dinv(D_t::Zero())
+    , UDinv(UD_t::Identity())
+    , StU(D_t::Zero())
+    , N(internal::SplineKinematics<_Scalar, _Options>::allocateBasis(degree, knot_size))
+    {
+    }
+
+    static std::string classname()
+    {
+      return std::string("JointDataSpline");
+    }
+    std::string shortname() const
+    {
+      return classname();
+    }
+
+  }; // struct JointDataSplinerTpl
+
+  PINOCCHIO_JOINT_CAST_TYPE_SPECIALIZATION(JointModelSplineTpl);
+
+  /// @brief Spline joint in \f$SE(3)\f$.
+  ///
+  /// A spline joint constrains the movement of the child frame to follow the spline defined by the
+  /// controlFrames and the degree of the spline. Implementation of the joint is based on the paper
+  /// from Lee et al. Spline Joints for Multibody Dynamics
+  /// (https://web.cs.ucla.edu/~dt/papers/siggraph08/siggraph08.pdf)
+  ///
+  template<typename _Scalar, int _Options>
+  struct JointModelSplineTpl : public JointModelBase<JointModelSplineTpl<_Scalar, _Options>>
+  {
+    EIGEN_MAKE_ALIGNED_OPERATOR_NEW
+    typedef JointSplineTpl<_Scalar, _Options> JointDerived;
+    typedef Eigen::Vector<_Scalar, Eigen::Dynamic> Vector;
+    PINOCCHIO_JOINT_TYPEDEF_TEMPLATE(JointDerived);
+
+    typedef JointModelBase<JointModelSplineTpl> Base;
+    using Base::id;
+    using Base::idx_q;
+    using Base::idx_v;
+    using Base::idx_vExtended;
+    using Base::setIndexes;
+
+    JointModelSplineTpl()
+    : degree(3)
+    , nbCtrlFrames(0)
+    , min_q(Scalar(0))
+    , max_q(Scalar(1))
+    {
+    }
+
+    JointModelSplineTpl(
+      const std::vector<Transformation_t> & controlFrames,
+      const Vector & knotVector,
+      const size_t degree)
+    : degree(degree)
+    , knots(knotVector)
+    , ctrlFrames(controlFrames)
+    {
+      if (controlFrames.size() <= degree)
+        PINOCCHIO_THROW_PRETTY(
+          std::invalid_argument,
+          "JointSpline - Number of control frames must be greater than degree of spline.");
+
+      nbCtrlFrames = controlFrames.size();
+      if (knotVector.size() != static_cast<Eigen::Index>(nbCtrlFrames + degree + 1))
+        PINOCCHIO_THROW_PRETTY(
+          std::invalid_argument,
+          "JointSpline - Size of knot vector should be nbControlFrames + degree + 1.");
+
+      size_t knot_multiplicity = 1;
+      for (Eigen::Index i = 1; i < knotVector.size(); ++i)
+      {
+        if (check_expression_if_real<Scalar>(knotVector[i] < knotVector[i - 1]))
+        {
+          PINOCCHIO_THROW_PRETTY(
+            std::invalid_argument, "JointSpline - Knot vector must be non-decreasing (knots must "
+                                   "satisfy knots[i] <= knots[i+1]).");
+        }
+
+        if (check_expression_if_real<Scalar>(knotVector[i] == knotVector[i - 1]))
+        {
+          knot_multiplicity++;
+          if (knot_multiplicity > degree + 1)
+          {
+            PINOCCHIO_THROW_PRETTY(
+              std::invalid_argument,
+              "JointSpline - Knot vector values cannot be repeated more than degree + 1 times.");
+          }
+        }
+        else
+        {
+          knot_multiplicity = 1;
+        }
+      }
+
+      min_q = knotVector[0];
+      max_q = knotVector[knots.size() - 1];
+
+      computeRelativeMotions();
+    }
+
+    JointDataDerived createData() const
+    {
+      return JointDataDerived(degree, knots.size());
+    }
+
+    const std::vector<bool> hasConfigurationLimit() const
+    {
+      return {true};
+    }
+
+    const std::vector<bool> hasConfigurationLimitInTangent() const
+    {
+      return {true};
+    }
+
+    using Base::isEqual;
+    bool isEqual(const JointModelSplineTpl & other) const
+    {
+      return Base::isEqual(other) && other.degree == degree && other.nbCtrlFrames == nbCtrlFrames
+             && other.ctrlFrames == ctrlFrames && other.knots == knots
+             && other.relativeMotions == relativeMotions;
+    }
+
+    template<typename ConfigVector>
+    void calc(JointDataDerived & data, const Eigen::MatrixBase<ConfigVector> & qs) const
+    {
+      assert(
+        check_expression_if_real<Scalar>(qs[0] >= min_q && qs[0] <= max_q)
+        && "Spline joint configuration (q) must be between min_q and max_q. ");
+
+      data.joint_q = qs.template segment<NQ>(idx_q());
+
+      internal::SplineKinematics<_Scalar, _Options>::compute(
+        degree, knots, ctrlFrames, relativeMotions, data.joint_q[0], data.joint_v, false, data.M,
+        data.v, data.c, data.S, data.N);
+    }
+
+    template<typename ConfigVector, typename TangentVector>
+    void calc(
+      JointDataDerived & data,
+      const Eigen::MatrixBase<ConfigVector> & qs,
+      const Eigen::MatrixBase<TangentVector> & vs) const
+    {
+      assert(
+        check_expression_if_real<Scalar>(qs[0] >= min_q && qs[0] <= max_q)
+        && "Spline joint configuration (q) must be between min_q and max_q. ");
+
+      data.joint_q = qs.template segment<NQ>(idx_q());
+      data.joint_v = vs.template segment<NV>(idx_v());
+
+      internal::SplineKinematics<_Scalar, _Options>::compute(
+        degree, knots, ctrlFrames, relativeMotions, data.joint_q[0], data.joint_v, true, data.M,
+        data.v, data.c, data.S, data.N);
+    }
+
+    template<typename TangentVector>
+    void
+    calc(JointDataDerived & data, const Blank, const Eigen::MatrixBase<TangentVector> & vs) const
+    {
+      data.joint_v = vs.template segment<NV>(idx_v());
+
+      internal::SplineKinematics<_Scalar, _Options>::compute(
+        degree, knots, ctrlFrames, relativeMotions, data.joint_q[0], data.joint_v, true, data.M,
+        data.v, data.c, data.S, data.N);
+    }
+
+    template<typename VectorLike, typename Matrix6Like>
+    void calc_aba(
+      JointDataDerived & data,
+      const Eigen::MatrixBase<VectorLike> & armature,
+      const Eigen::MatrixBase<Matrix6Like> & I,
+      const bool update_I) const
+    {
+      data.U.noalias() = I * data.S.matrix();
+      data.StU.noalias() = data.S.transpose() * data.U;
+      data.StU.diagonal() += armature;
+      internal::matrix_inversion(data.StU, data.Dinv);
+
+      data.UDinv.noalias() = data.U * data.Dinv;
+
+      if (update_I)
+        PINOCCHIO_EIGEN_CONST_CAST(Matrix6Like, I).noalias() -= data.UDinv * data.U.transpose();
+    }
+
+    static std::string classname()
+    {
+      return std::string("JointModelSpline");
+    }
+    std::string shortname() const
+    {
+      return classname();
+    }
+
+    template<typename NewScalar>
+    JointModelSplineTpl<NewScalar, Options> cast() const
+    {
+      typedef JointModelSplineTpl<NewScalar, Options> ReturnType;
+      ReturnType res;
+      res.degree = degree;
+      res.nbCtrlFrames = nbCtrlFrames;
+      res.ctrlFrames.reserve(ctrlFrames.size());
+      res.relativeMotions.reserve(relativeMotions.size());
+      for (const auto & cf : ctrlFrames)
+      {
+        res.ctrlFrames.push_back(cf.template cast<NewScalar>());
+      }
+      for (const auto & rm : relativeMotions)
+      {
+        res.relativeMotions.push_back(rm.template cast<NewScalar>());
+      }
+
+      res.min_q = ScalarCast<NewScalar, Scalar>::cast(min_q);
+      res.max_q = ScalarCast<NewScalar, Scalar>::cast(max_q);
+      res.knots = knots.template cast<NewScalar>();
+
+      res.setIndexes(id(), idx_q(), idx_v(), idx_vExtended());
+      return res;
+    }
+
+    // attributes
+    size_t degree;
+    size_t nbCtrlFrames;
+    Vector knots;
+    Scalar min_q;
+    Scalar max_q;
+
+    std::vector<Transformation_t> ctrlFrames;
+    std::vector<Motion_t> relativeMotions;
+
+  private:
+    void computeRelativeMotions()
+    {
+      for (size_t i = 0; i < nbCtrlFrames - 1; i++)
+        relativeMotions.push_back(log6(ctrlFrames[i].inverse() * ctrlFrames[i + 1]));
+    }
+
+  }; // struct JointModelSplineTpl
+
+  /// @brief Helper structure to specify attributes of a spline joint.
+  template<typename _Scalar, int _Options>
+  struct JointModelSplineBuilderTpl
+  {
+    using Scalar = _Scalar;
+    static constexpr int Options = _Options;
+
+    using JointModel_t = JointModelSplineTpl<Scalar, Options>;
+    using Transformation_t = typename JointModel_t::Transformation_t;
+    using Vector = Eigen::Matrix<Scalar, Eigen::Dynamic, 1>;
+
+    enum class KnotPolicy
+    {
+      OpenUniform,
+      Uniform,
+      Custom
+    };
+
+    JointModelSplineBuilderTpl()
+    : degree_(3)
+    , min_q_(Scalar(0))
+    , max_q_(Scalar(1))
+    , knot_policy_(KnotPolicy::OpenUniform)
+    {
+    }
+
+    JointModelSplineBuilderTpl & addControlFrame(const Transformation_t & frame)
+    {
+      ctrlFrames_.push_back(frame);
+      return *this;
+    }
+
+    JointModelSplineBuilderTpl &
+    withControlFrameVector(const std::vector<Transformation_t> & frames)
+    {
+      ctrlFrames_ = frames;
+      return *this;
+    }
+
+    JointModelSplineBuilderTpl & withDegree(size_t degree)
+    {
+      degree_ = degree;
+      return *this;
+    }
+
+    JointModelSplineBuilderTpl & withKnotVector(const Vector & knots)
+    {
+      knots_ = knots;
+      knot_policy_ = KnotPolicy::Custom;
+
+      return *this;
+    }
+
+    JointModelSplineBuilderTpl & withKnotVector(const std::vector<Scalar> & knots)
+    {
+      knots_.resize(knots.size());
+      for (std::size_t i = 0; i < knots.size(); ++i)
+      {
+        knots_[i] = knots[i];
+      }
+      knot_policy_ = KnotPolicy::Custom;
+
+      return *this;
+    }
+
+    JointModelSplineBuilderTpl & withOpenUniformKnots(const Scalar min_q, const Scalar max_q)
+    {
+      knot_policy_ = KnotPolicy::OpenUniform;
+      min_q_ = min_q;
+      max_q_ = max_q;
+
+      return *this;
+    }
+
+    JointModelSplineBuilderTpl & withUniformKnots(const Scalar min_q, const Scalar max_q)
+    {
+      knot_policy_ = KnotPolicy::Uniform;
+      min_q_ = min_q;
+      max_q_ = max_q;
+
+      return *this;
+    }
+
+    JointModel_t build() const
+    {
+      Vector knots;
+
+      const size_t nCtrl = ctrlFrames_.size();
+      switch (knot_policy_)
+      {
+      case KnotPolicy::OpenUniform:
+        knots = internal::generateOpenUniformKnots(min_q_, max_q_, nCtrl, degree_);
+        break;
+
+      case KnotPolicy::Uniform:
+        knots = internal::generateUniformKnots(min_q_, max_q_, nCtrl, degree_);
+        break;
+
+      case KnotPolicy::Custom:
+        knots = knots_;
+        break;
+      default:
+        break;
+      }
+      return JointModel_t(ctrlFrames_, knots, degree_);
+    }
+
+  private:
+    std::vector<Transformation_t> ctrlFrames_;
+
+    size_t degree_;
+
+    Scalar min_q_;
+    Scalar max_q_;
+
+    Vector knots_;
+    KnotPolicy knot_policy_;
+  }; // struct JointModelSplineBuilder
+
+} // namespace pinocchio
+
+#include <boost/type_traits.hpp>
+
+namespace boost
+{
+  template<typename Scalar, int Options>
+  struct has_nothrow_constructor<::pinocchio::JointModelSplineTpl<Scalar, Options>>
+  : public integral_constant<bool, true>
+  {
+  };
+
+  template<typename Scalar, int Options>
+  struct has_nothrow_copy<::pinocchio::JointModelSplineTpl<Scalar, Options>>
+  : public integral_constant<bool, true>
+  {
+  };
+
+  template<typename Scalar, int Options>
+  struct has_nothrow_constructor<::pinocchio::JointDataSplineTpl<Scalar, Options>>
+  : public integral_constant<bool, true>
+  {
+  };
+
+  template<typename Scalar, int Options>
+  struct has_nothrow_copy<::pinocchio::JointDataSplineTpl<Scalar, Options>>
+  : public integral_constant<bool, true>
+  {
+  };
+} // namespace boost
diff --git a/include/pinocchio/src/multibody/joint/spline-utils.hxx b/include/pinocchio/src/multibody/joint/spline-utils.hxx
new file mode 100644
index 0000000000..fa54be3c78
--- /dev/null
+++ b/include/pinocchio/src/multibody/joint/spline-utils.hxx
@@ -0,0 +1,677 @@
+//
+// Copyright (c) 2025 INRIA
+// Copyright (c) 2026 ISIR
+//
+
+#pragma once
+
+// IWYU pragma: private, include "pinocchio/multibody/joint.hpp
+
+#ifdef PINOCCHIO_LSP
+  #undef PINOCCHIO_LSP
+  #include "pinocchio/multibody/joint.hpp"
+#endif // PINOCCHIO_LSP
+
+namespace pinocchio
+{
+  namespace internal
+  {
+    // Compute the only positive degree 0 basis function index.
+    template<typename Scalar>
+    int findSpan(int degree, const Eigen::Matrix<Scalar, Eigen::Dynamic, 1> & knots, Scalar q)
+    {
+      assert(degree >= 0);
+      assert(knots.size() > degree + 1);
+      assert(knots[degree] <= q);
+      assert(q <= knots[knots.size() - 1 - degree]);
+
+      int low = degree;
+      for (; low < knots.size() - 1 - degree; ++low)
+      {
+        if (knots[low] <= q && q < knots[low + 1])
+        {
+          return low;
+        }
+      }
+
+      return low - 1;
+    }
+
+    template<typename Scalar>
+    Eigen::Matrix<Scalar, Eigen::Dynamic, 1>
+    generateOpenUniformKnots(const Scalar min_q, const Scalar max_q, size_t nCtrl, size_t degree)
+    {
+      using Vector = Eigen::Matrix<Scalar, Eigen::Dynamic, 1>;
+
+      const size_t n_knots = nCtrl + degree + 1;
+
+      Vector knots;
+      knots.resize(static_cast<Eigen::Index>(n_knots));
+
+      const Scalar range = max_q - min_q;
+
+      knots.head(degree + 1).setConstant(min_q);
+      const Scalar nInner = static_cast<Scalar>(nCtrl - degree - 1);
+      const Scalar denominator = static_cast<Scalar>(nInner + 1);
+
+      for (size_t i = degree + 1; i < nCtrl; i++)
+        knots[static_cast<Eigen::Index>(i)] =
+          min_q + range * static_cast<Scalar>(i - degree) / denominator;
+
+      knots.tail(degree + 1).setConstant(max_q);
+      return knots;
+    }
+
+    template<typename Scalar>
+    Eigen::Matrix<Scalar, Eigen::Dynamic, 1>
+    generateUniformKnots(const Scalar min_q, const Scalar max_q, size_t nCtrl, size_t degree)
+    {
+      using Vector = Eigen::Matrix<Scalar, Eigen::Dynamic, 1>;
+
+      const size_t n_knots = nCtrl + degree + 1;
+      Vector knots;
+      knots.resize(static_cast<Eigen::Index>(n_knots));
+
+      const Scalar step = (max_q - min_q) / static_cast<Scalar>(n_knots - 1);
+
+      for (size_t i = 0; i < n_knots; ++i)
+        knots[static_cast<Eigen::Index>(i)] = min_q + step * static_cast<Scalar>(i);
+
+      return knots;
+    }
+
+    /** De Boor algorithm modification to compute all basis involved to compute one
+     * point of the curve.
+     * \param degree Curve degree.
+     * \param knots Knot vector at least of size \p degree + 1.
+     * \param root_basis Degree 0 basis function index where knot vector contains \p x.
+     * This is the only degree 0 basis function != 0 and can be computed by \p findSpan.
+     * \param q Value to evaluate.
+     * \param basis of size (degree + 1, degree + 1).
+     * Each element i, j of this array will hold a basis function N_{i,j} where i and j
+     * are respectively the basis function index and degree.
+     */
+    template<typename Scalar>
+    void deBoorBasis(
+      int degree,
+      const Eigen::Matrix<Scalar, Eigen::Dynamic, 1> & knots,
+      int root_basis,
+      Scalar q,
+      Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic> & basis)
+    {
+      assert(degree >= 0);
+      assert(basis.rows() >= degree + 1);
+      assert(basis.cols() >= degree + 1);
+      assert(knots.size() > degree + 1);
+      assert(degree <= root_basis);
+      assert(root_basis < knots.size() - 1 - degree);
+      assert(knots[root_basis] <= q);
+      assert(q <= knots[root_basis + 1]);
+
+      /*
+       * N_{i,j} is a basis function where i and j are respectively the index and the degree.
+       * Here N_{i,2} computation scheme:
+       *
+       *  N_{0,2}       N_{1,2}       N_{2,2}
+       *   (1,3)\  (1,3)/ (2,4)\ (2,4)/
+       *         N_{1,1}       N_{2,1}
+       *        /  (2,3)\ (2,3)/     \
+       *  N_{1,0}       N_{2,0}       N_{3,0}
+       * [u1, u2[      [u2, u3[      [u3, u4[
+       *
+       * When u is in [u2, u3[ range N_{1,0} and N_{3,0} basis function value to 0.
+       * This simplify the computation scheme by allowing to only compute
+       * N_{1,1}, N_{2,1}, N_{0,2}, N_{1,2} and N_{2,2}.
+       * N_{1,1}, N_{0,2}, N_{2,1} and N_{2,2} can be computed only with one alpha value.
+       * We do it in the first pass where we compute left most and right most basis functions.
+       * Only N_{1,2} will need both alpha to be computed. We do that in the second pass.
+       */
+
+      // Compute left most and right most basis functions (first pass).
+      basis(0, 0) = Scalar(1);
+      for (int previous_degree = 0; previous_degree < degree; ++previous_degree)
+      {
+        const int current_degree = previous_degree + 1;
+        const int left_most_basis = root_basis - current_degree;
+        const int left_most_basis_start_knot = left_most_basis + 1;
+        const int left_most_basis_end_knot = left_most_basis_start_knot + current_degree;
+        const Scalar left_most_basis_alpha =
+          (knots[left_most_basis_end_knot] - q)
+          / (knots[left_most_basis_end_knot] - knots[left_most_basis_start_knot]);
+        basis(current_degree, 0) = left_most_basis_alpha * basis(previous_degree, 0);
+
+        const int right_most_basis = root_basis;
+        const int right_most_basis_start_knot = right_most_basis;
+        const int right_most_basis_end_knot = right_most_basis_start_knot + current_degree;
+        const Scalar right_most_basis_alpha =
+          (q - knots[right_most_basis_start_knot])
+          / (knots[right_most_basis_end_knot] - knots[right_most_basis_start_knot]);
+        basis(current_degree, current_degree) =
+          (right_most_basis_alpha * basis(previous_degree, previous_degree));
+      }
+
+      // Compute central basis functions (second pass).
+      for (int previous_degree = 1; previous_degree < degree; ++previous_degree)
+      {
+        const int current_degree = previous_degree + 1;
+        const int left_most_basis = root_basis - current_degree;
+        const int basis_numbers = current_degree + 1;
+        for (int i = 1; i < basis_numbers - 1; ++i)
+        {
+          const int current_basis = left_most_basis + i;
+          const int left_side_start_knot = current_basis;
+          const int left_side_end_knot = current_basis + current_degree;
+          const Scalar left_side_alpha =
+            (q - knots[left_side_start_knot])
+            / (knots[left_side_end_knot] - knots[left_side_start_knot]);
+
+          const int right_side_start_knot = left_side_start_knot + 1;
+          const int right_side_end_knot = left_side_end_knot + 1;
+          const Scalar right_side_alpha =
+            (knots[right_side_end_knot] - q)
+            / (knots[right_side_end_knot] - knots[right_side_start_knot]);
+
+          basis(current_degree, i) = left_side_alpha * basis(previous_degree, i - 1)
+                                     + right_side_alpha * basis(previous_degree, i);
+        }
+      }
+    }
+
+    /** Return num / den if den != 0 and 0 in the other case.
+     * This function help support Casadi scalar type that doesn't support
+     * comparison.
+     */
+    template<typename Scalar>
+    Scalar safeAlpha(Scalar num, Scalar den)
+    {
+      // clang-format off
+      // if(den > dummy_precision)
+      //  return (num / den)
+      // else
+      //  return 0
+      // clang-format on
+      return if_then_else(
+        GT, den, Eigen::NumTraits<Scalar>::dummy_precision(), (num / den), Scalar(0));
+    }
+
+    /** De Boor algorithm modification to compute all basis in the spline
+     * valid span (knots[degree: m - degree]).
+     * This function will compute a lot of zeros and it's implemented for
+     * Casadi scalar support. Use deBoorBasis instead.
+     * \param degree Curve degree.
+     * \param knots Knot vector of size m (at least of size \p degree + 1)
+     * \param q Value to evaluate.
+     * \param basis of size (degree + 1, m - degree - 1).
+     * Each element i, j of this array will hold a basis function N_{i,j} where i and j
+     * are respectively the basis function index and degree.
+     */
+    template<typename Scalar>
+    void deBoorFullBasis(
+      int degree,
+      const Eigen::Matrix<Scalar, Eigen::Dynamic, 1> & knots,
+      Scalar q,
+      Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic> & basis)
+    {
+      assert(degree >= 0);
+      assert(basis.rows() >= degree + 1);
+      // Number of max degree basis functions
+      assert(basis.cols() >= knots.size() - degree - 1);
+      assert(knots.size() > degree + 1);
+      basis.setZero();
+      const int first_degree0_basis = degree;
+      const int last_degree0_basis = knots.size() - degree - 2;
+      const int nb_degree0_basis = last_degree0_basis + 1 - first_degree0_basis;
+
+      // Compute degree 0 basis function values
+      for (int i = 0; i < nb_degree0_basis; ++i)
+      {
+        int current_basis = first_degree0_basis + i;
+        // clang-format off
+        // if(knots[i] <= x && x < knots[i + 1])
+        //  return 1;
+        // else
+        //  return 0;
+        // clang-format on
+        Scalar is_in_standard_range = if_then_else(
+          LE, knots[current_basis], q,
+          if_then_else(LT, q, knots[current_basis + 1], Scalar(1), Scalar(0)), Scalar(0));
+
+        // clang-format off
+        // if(x == knots.back() && x == knots[i + 1])
+        //  return 1;
+        // else
+        //  return 0;
+        // clang-format on
+        Scalar is_at_final_range = if_then_else(
+          EQ, q, knots[knots.size() - degree - 1],
+          if_then_else(EQ, q, knots[current_basis + 1], Scalar(1), Scalar(0)), Scalar(0));
+
+        basis(0, i) = is_in_standard_range + is_at_final_range;
+      }
+
+      // Compute left most and right most basis functions (first pass).
+      for (int previous_degree = 0; previous_degree < degree; ++previous_degree)
+      {
+        const int current_degree = previous_degree + 1;
+        const int basis_numbers = nb_degree0_basis + current_degree;
+        const int left_most_basis = first_degree0_basis - current_degree;
+        const int left_most_basis_start_knot = left_most_basis + 1;
+        const int left_most_basis_end_knot = left_most_basis_start_knot + current_degree;
+        const Scalar left_most_basis_alpha_num = knots[left_most_basis_end_knot] - q;
+        const Scalar left_most_basis_alpha_den =
+          knots[left_most_basis_end_knot] - knots[left_most_basis_start_knot];
+        basis(current_degree, 0) = safeAlpha(left_most_basis_alpha_num, left_most_basis_alpha_den)
+                                   * basis(previous_degree, 0);
+
+        const int right_most_basis = last_degree0_basis;
+        const int right_most_basis_start_knot = right_most_basis;
+        const int right_most_basis_end_knot = right_most_basis_start_knot + current_degree;
+        const Scalar right_most_basis_alpha_num = q - knots[right_most_basis_start_knot];
+        const Scalar right_most_basis_alpha_den =
+          knots[right_most_basis_end_knot] - knots[right_most_basis_start_knot];
+        basis(current_degree, basis_numbers - 1) =
+          safeAlpha(right_most_basis_alpha_num, right_most_basis_alpha_den)
+          * basis(previous_degree, basis_numbers - 2);
+      }
+
+      // Compute central basis functions (second pass).
+      for (int previous_degree = 0; previous_degree < degree; ++previous_degree)
+      {
+        const int current_degree = previous_degree + 1;
+        const int left_most_basis = first_degree0_basis - current_degree;
+        const int basis_numbers = nb_degree0_basis + current_degree;
+        for (int i = 1; i < basis_numbers - 1; ++i)
+        {
+          const int current_basis = left_most_basis + i;
+          const int left_side_start_knot = current_basis;
+          const int left_side_end_knot = current_basis + current_degree;
+          const Scalar left_side_alpha_num = (q - knots[left_side_start_knot]);
+          const Scalar left_side_alpha_den =
+            (knots[left_side_end_knot] - knots[left_side_start_knot]);
+
+          const int right_side_start_knot = left_side_start_knot + 1;
+          const int right_side_end_knot = left_side_end_knot + 1;
+          const Scalar right_side_alpha_num = (knots[right_side_end_knot] - q);
+          const Scalar right_side_alpha_den =
+            (knots[right_side_end_knot] - knots[right_side_start_knot]);
+
+          basis(current_degree, i) =
+            safeAlpha(left_side_alpha_num, left_side_alpha_den) * basis(previous_degree, i - 1)
+            + safeAlpha(right_side_alpha_num, right_side_alpha_den) * basis(previous_degree, i);
+        }
+      }
+    }
+
+    /** Return basis function value N_{index,degree} from basis matrix computed by \p deBoorBasis.
+     * \param root_basis Argument provided to \p deBoorBasis function.
+     * \param basis Basis matrix computed by \p deBoorBasis.
+     * \param index Index of the basis function.
+     * \param degree Degree of the basis function.
+     */
+    template<typename Scalar>
+    const Scalar & getAbsoluteBasis(
+      int root_basis,
+      const Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic> & basis,
+      int index,
+      int degree)
+    {
+      assert(0 <= root_basis);
+      assert(0 <= index);
+      assert(0 <= degree);
+      assert(degree < basis.rows());
+
+      const int offset = root_basis - degree;
+      assert(offset <= index);
+
+      return basis(degree, index - offset);
+    }
+
+    /** Return basis function value N_{index,degree} from basis matrix computed by \p
+     * deBoorBasisFull.
+     * \param spline_degree Spline degree provided to \p deBoorBasis function.
+     * \param basis Basis matrix computed by \p deBoorBasis.
+     * \param index Index of the basis function.
+     * \param degree Degree of the basis function.
+     */
+    template<typename Scalar>
+    const Scalar & getAbsoluteBasisFull(
+      int spline_degree,
+      const Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic> & basis,
+      int index,
+      int degree)
+    {
+      assert(0 <= spline_degree);
+      assert(0 <= index);
+      assert(0 <= degree);
+      assert(degree <= spline_degree);
+      assert(degree < basis.rows());
+
+      const int offset = spline_degree - degree;
+      assert(offset <= index);
+
+      return basis(degree, index - offset);
+    }
+
+    /** Compute cumulative basis first derivative for N_{index,degree}.
+     * \param root_basis Argument provided to \p deBoorBasis function.
+     * \param knots Knot vector at least of size \p degree + 1.
+     * \param basis Basis matrix computed by \p deBoorBasis.
+     * \param index Index of the basis function.
+     * \param degree Degree of the basis function.
+     */
+    template<typename Scalar>
+    Scalar cumulativeBasisDerivative(
+      int root_basis,
+      const Eigen::Matrix<Scalar, Eigen::Dynamic, 1> & knots,
+      const Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic> & basis,
+      int index,
+      int degree)
+    {
+      assert(0 <= index);
+      assert(0 <= degree);
+      assert(index + degree < knots.size());
+
+      const Scalar alpha = degree / (knots[index + degree] - knots[index]);
+      return alpha * getAbsoluteBasis(root_basis, basis, index, degree - 1);
+    }
+
+    /** Compute cumulative basis first derivative for N_{index,degree} on nodal vector computed by
+     * \p deBoorBasisFull.
+     * This method is safe to call even on basis function that span on a knot vector with same
+     * values.
+     * \param spline_degree Spline degree provided to \p deBoorBasis function.
+     * \param knots Knot vector at least of size \p degree + 1.
+     * \param basis Basis matrix computed by \p deBoorBasis.
+     * \param index Index of the basis function.
+     * \param degree Degree of the basis function.
+     */
+    template<typename Scalar>
+    Scalar cumulativeBasisDerivativeFull(
+      int spline_degree,
+      const Eigen::Matrix<Scalar, Eigen::Dynamic, 1> & knots,
+      const Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic> & basis,
+      int index,
+      int degree)
+    {
+      assert(0 <= index);
+      assert(0 <= degree);
+      assert(index + degree < knots.size());
+
+      // deBoorBasisFull only compute basis function in the spline valide span
+      // (knots[degree:m-degree]). We return zero if we try to access basis function
+      // outside of this span.
+      if (index > spline_degree - degree)
+      {
+        return safeAlpha(static_cast<Scalar>(degree), knots[index + degree] - knots[index])
+               * getAbsoluteBasisFull(spline_degree, basis, index, degree - 1);
+      }
+      else
+      {
+        return Scalar(0);
+      }
+    }
+
+    /** Compute cumulative basis second derivative for N_{index,degree}.
+     * \param root_basis Argument provided to \p deBoorBasis function.
+     * \param knots Knot vector at least of size \p degree + 1.
+     * \param basis Basis matrix computed by \p deBoorBasis.
+     * \param index Index of the basis function.
+     * \param degree Degree of the basis function.
+     */
+    template<typename Scalar>
+    Scalar cumulativeBasisDerivative2(
+      int root_basis,
+      const Eigen::Matrix<Scalar, Eigen::Dynamic, 1> & knots,
+      const Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic> & basis,
+      int index,
+      int degree)
+    {
+      assert(0 <= index);
+      assert(0 <= degree);
+      assert(index + degree + 1 < knots.size());
+
+      const int derivative1_degree = degree - 1;
+      const int derivative2_degree = degree - 2;
+      const Scalar derivative1_den = knots[index + degree] - knots[index];
+      Scalar phi_ddot_i_sum = Scalar(0);
+      // basis only contains non zero basis function.
+      // This condition prevent get an out of bound basis function on the left.
+      if (index >= root_basis - derivative2_degree)
+      {
+        const Scalar left_side_den = knots[index + derivative1_degree] - knots[index];
+        phi_ddot_i_sum =
+          getAbsoluteBasis(root_basis, basis, index, derivative2_degree) / left_side_den;
+      }
+      // basis only contains non zero basis function.
+      // This condition prevent get an out of bound basis function on the right.
+      if (index + 1 < root_basis + 1)
+      {
+        const Scalar right_side_den = knots[index + derivative1_degree + 1] - knots[index + 1];
+        phi_ddot_i_sum -=
+          getAbsoluteBasis(root_basis, basis, index + 1, derivative2_degree) / right_side_den;
+      }
+      return ((degree * derivative1_degree) / derivative1_den) * phi_ddot_i_sum;
+    }
+
+    /** Compute cumulative basis second derivative for N_{index,degree}
+     * on nodal vector computed by \p deBoorBasisFull.
+     * This method is safe to call even on basis function that span on a knot vector with same
+     * values.
+     * \param spline_degree Spline degree provided to \p deBoorBasis function.
+     * \param knots Knot vector at least of size \p degree + 1.
+     * \param basis Basis matrix computed by \p deBoorBasis.
+     * \param index Index of the basis function.
+     * \param degree Degree of the basis function.
+     */
+    template<typename Scalar>
+    Scalar cumulativeBasisDerivative2Full(
+      int spline_degree,
+      const Eigen::Matrix<Scalar, Eigen::Dynamic, 1> & knots,
+      const Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic> & basis,
+      int index,
+      int degree)
+    {
+      assert(0 <= index);
+      assert(0 <= degree);
+      assert(index + degree + 1 < knots.size());
+
+      // deBoorBasisFull only compute basis function in the spline valide span
+      // (knots[degree:m-degree]). We return zero if we try to access basis function
+      // outside of this span.
+      if (index <= spline_degree - degree)
+      {
+        return Scalar(0);
+      }
+
+      const int derivative1_degree = degree - 1;
+      const int derivative2_degree = degree - 2;
+      const Scalar derivative1_den = knots[index + degree] - knots[index];
+      Scalar phi_ddot_i_sum = Scalar(0);
+      // Like first guard, we set this part of the computation to zero if basis function
+      // needed to compute the second derivative is outside of the spline valide span.
+      if (index >= spline_degree - derivative2_degree)
+      {
+        const Scalar left_side_den = knots[index + derivative1_degree] - knots[index];
+        phi_ddot_i_sum = getAbsoluteBasisFull(spline_degree, basis, index, derivative2_degree)
+                         * safeAlpha(Scalar(1), left_side_den);
+      }
+      // All last valid basis have the same index
+      const int last_degree0_basis = knots.size() - degree - 2;
+      if (index + 1 <= last_degree0_basis)
+      {
+        const Scalar right_side_den = knots[index + derivative1_degree + 1] - knots[index + 1];
+        phi_ddot_i_sum -= getAbsoluteBasisFull(spline_degree, basis, index + 1, derivative2_degree)
+                          * safeAlpha(Scalar(1), right_side_den);
+      }
+      return ((degree * derivative1_degree) * safeAlpha(Scalar(1), derivative1_den))
+             * phi_ddot_i_sum;
+    }
+
+    /** Compute spline joint kinematics data.
+     */
+    template<typename Scalar, int Options>
+    void computeSplineKinematics(
+      int degree,
+      const Eigen::Matrix<Scalar, Eigen::Dynamic, 1> & knots,
+      const std::vector<SE3Tpl<Scalar, Options>> & ctrlFrames,
+      const std::vector<MotionTpl<Scalar, Options>> & relativeMotions,
+      Scalar q,
+      const Eigen::Matrix<Scalar, 1, 1, Options> & joint_v,
+      bool computeVelocity,
+      SE3Tpl<Scalar, Options> & M,
+      MotionTpl<Scalar, Options> & v,
+      MotionTpl<Scalar, Options> & c,
+      JointMotionSubspaceTpl<1, Scalar, Options, 1> & S,
+      Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic> & basis)
+    {
+      int root_basis_degree0 = findSpan(degree, knots, q);
+      int root_basis = root_basis_degree0 - degree;
+      deBoorBasis(degree, knots, root_basis_degree0, q, basis);
+
+      M = ctrlFrames[root_basis];
+      S.matrix().setZero();
+      if (computeVelocity)
+      {
+        v.setZero();
+        c.setZero();
+      }
+
+      for (int i = 1; i < degree + 1; i++)
+      {
+        const int current_basis = root_basis + i;
+
+        const Scalar phi_i = basis.row(degree).segment(i, degree + 1 - i).sum();
+
+        const Scalar phi_dot_i = internal::cumulativeBasisDerivative(
+          root_basis_degree0, knots, basis, current_basis, degree);
+
+        const SE3Tpl<Scalar, Options> transformation_temp(
+          exp6(relativeMotions[current_basis - 1] * phi_i));
+        M = M * transformation_temp;
+
+        if (computeVelocity)
+        {
+          const Scalar phi_ddot_i = internal::cumulativeBasisDerivative2(
+            root_basis_degree0, knots, basis, current_basis, degree);
+
+          c = relativeMotions[current_basis - 1] * phi_ddot_i
+              + transformation_temp.actInv(
+                c
+                + MotionTpl<Scalar, Options>(S.matrix()).cross(relativeMotions[current_basis - 1])
+                    * phi_dot_i);
+        }
+
+        S.matrix() =
+          transformation_temp.actInv(S) + relativeMotions[current_basis - 1].toVector() * phi_dot_i;
+      }
+      if (computeVelocity)
+      {
+        // C = Sdot * qdot = (dS/dq * qdot) * dot
+        c = c * joint_v[0] * joint_v[0];
+        v = S * joint_v;
+      }
+    }
+
+    /** Compute spline joint kinematics data on nodal vector computed by \p deBoorBasisFull.
+     * This version is less efficient than computeSplineKinematics and is dedicated to Casadi.
+     */
+    template<typename Scalar, int Options>
+    void computeSplineKinematicsFull(
+      int degree,
+      const Eigen::Matrix<Scalar, Eigen::Dynamic, 1> & knots,
+      const std::vector<SE3Tpl<Scalar, Options>> & ctrlFrames,
+      const std::vector<MotionTpl<Scalar, Options>> & relativeMotions,
+      Scalar q,
+      const Eigen::Matrix<Scalar, 1, 1, Options> & joint_v,
+      bool computeVelocity,
+      SE3Tpl<Scalar, Options> & M,
+      MotionTpl<Scalar, Options> & v,
+      MotionTpl<Scalar, Options> & c,
+      JointMotionSubspaceTpl<1, Scalar, Options, 1> & S,
+      Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic> & basis)
+    {
+      deBoorFullBasis(degree, knots, q, basis);
+
+      M = ctrlFrames[0];
+      S.matrix().setZero();
+      if (computeVelocity)
+      {
+        v.setZero();
+        c.setZero();
+      }
+
+      const int nb_basis = ctrlFrames.size();
+      for (int i = 1; i < nb_basis; i++)
+      {
+        const int current_basis = i;
+
+        const Scalar phi_i = basis.row(degree).tail(nb_basis - i).sum();
+
+        const Scalar phi_dot_i =
+          internal::cumulativeBasisDerivativeFull(degree, knots, basis, current_basis, degree);
+
+        const SE3Tpl<Scalar, Options> transformation_temp(
+          exp6(relativeMotions[current_basis - 1] * phi_i));
+        M = M * transformation_temp;
+
+        if (computeVelocity)
+        {
+          const Scalar phi_ddot_i =
+            internal::cumulativeBasisDerivative2Full(degree, knots, basis, current_basis, degree);
+
+          c = relativeMotions[current_basis - 1] * phi_ddot_i
+              + transformation_temp.actInv(
+                c
+                + MotionTpl<Scalar, Options>(S.matrix()).cross(relativeMotions[current_basis - 1])
+                    * phi_dot_i);
+        }
+
+        S.matrix() =
+          transformation_temp.actInv(S) + relativeMotions[current_basis - 1].toVector() * phi_dot_i;
+      }
+      if (computeVelocity)
+      {
+        // C = Sdot * qdot = (dS/dq * qdot) * dot
+        c = c * joint_v[0] * joint_v[0];
+        v = S * joint_v;
+      }
+    }
+
+    /** This structure allow with partial template specialization to
+     * use the right basis matrix allocation and computeSplineKinematics method for scalar type.
+     */
+    template<typename Scalar, int Options>
+    struct SplineKinematics
+    {
+      /// \return Basis vector stored in JointDataSplineTpl.
+      static Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic>
+      allocateBasis(int degree, int knot_size)
+      {
+        PINOCCHIO_UNUSED_VARIABLE(knot_size);
+        return Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic>::Zero(degree + 1, degree + 1);
+      }
+
+      /// Compute the spline kinematics.
+      static void compute(
+        int degree,
+        const Eigen::Matrix<Scalar, Eigen::Dynamic, 1> & knots,
+        const std::vector<SE3Tpl<Scalar, Options>> & ctrlFrames,
+        const std::vector<MotionTpl<Scalar, Options>> & relativeMotions,
+        Scalar q,
+        const Eigen::Matrix<Scalar, 1, 1, Options> & joint_v,
+        bool computeVelocity,
+        SE3Tpl<Scalar, Options> & M,
+        MotionTpl<Scalar, Options> & v,
+        MotionTpl<Scalar, Options> & c,
+        JointMotionSubspaceTpl<1, Scalar, Options, 1> & S,
+        Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic> & basis)
+      {
+        return computeSplineKinematics(
+          degree, knots, ctrlFrames, relativeMotions, q, joint_v, computeVelocity, M, v, c, S,
+          basis);
+      }
+    };
+
+  } // namespace internal
+
+} // namespace pinocchio
diff --git a/include/pinocchio/src/parsers/graph/fwd.hxx b/include/pinocchio/src/parsers/graph/fwd.hxx
index cd671ea45d..3595272792 100644
--- a/include/pinocchio/src/parsers/graph/fwd.hxx
+++ b/include/pinocchio/src/parsers/graph/fwd.hxx
@@ -35,6 +35,7 @@ namespace pinocchio
     struct JointUniversal;
     struct JointComposite;
     struct JointMimic;
+    struct JointSpline;
 
     // frames.hpp
     struct BodyFrame;
diff --git a/include/pinocchio/src/parsers/graph/graph-visitor.hxx b/include/pinocchio/src/parsers/graph/graph-visitor.hxx
index 705f9b7ff6..2a5b2a6bc6 100644
--- a/include/pinocchio/src/parsers/graph/graph-visitor.hxx
+++ b/include/pinocchio/src/parsers/graph/graph-visitor.hxx
@@ -121,6 +121,13 @@ namespace pinocchio
             std::invalid_argument,
             "Graph - Joint Mimic cannot have a q_ref. Please use the joint offset directly.");
         }
+
+        SE3 operator()(const JointSpline &) const
+        {
+          PINOCCHIO_THROW_PRETTY(
+            std::invalid_argument,
+            "Graph - Joint Spline cannot have a q_ref. Please use the joint offset directly.");
+        }
       };
 
       /// Return qref transformation.
diff --git a/include/pinocchio/src/parsers/graph/joints.hxx b/include/pinocchio/src/parsers/graph/joints.hxx
index fdd9236da3..cc6fb55aba 100644
--- a/include/pinocchio/src/parsers/graph/joints.hxx
+++ b/include/pinocchio/src/parsers/graph/joints.hxx
@@ -248,6 +248,141 @@ namespace pinocchio
       }
     };
 
+    struct JointSpline
+    {
+      std::vector<SE3> ctrlFrames;
+      Eigen::VectorXd knots;
+      std::size_t degree = 3;
+
+      static constexpr int nq = 1;
+      static constexpr int nv = 1;
+
+      JointSpline() = default;
+      JointSpline(std::size_t degree)
+      : degree(degree)
+      {
+      }
+
+      JointSpline(
+        const std::vector<SE3> & ctrlFrames, const Eigen::VectorXd & knots, std::size_t degree)
+      : ctrlFrames(ctrlFrames)
+      , knots(knots)
+      , degree(degree)
+      {
+      }
+
+      bool operator==(const JointSpline & other) const
+      {
+        return ctrlFrames == other.ctrlFrames && knots == other.knots && degree == other.degree;
+      }
+    };
+
+    struct JointSplineBuilder
+    {
+      using JointModelSplineBuilder = JointModelSplineBuilderTpl<double>;
+      using KnotPolicy = typename JointModelSplineBuilder::KnotPolicy;
+
+      JointSplineBuilder()
+      : degree(3)
+      , min_q(0)
+      , max_q(1)
+      , knot_policy(KnotPolicy::OpenUniform)
+      {
+      }
+
+      JointSplineBuilder & addControlFrame(const SE3 & frame)
+      {
+        ctrlFrames.push_back(frame);
+        return *this;
+      }
+
+      JointSplineBuilder & withControlFrameVector(const std::vector<SE3> & frames)
+      {
+        ctrlFrames = frames;
+        return *this;
+      }
+
+      JointSplineBuilder & withDegree(size_t p_degree)
+      {
+        degree = p_degree;
+        return *this;
+      }
+
+      JointSplineBuilder & withKnotVector(const Eigen::VectorXd & p_knots)
+      {
+        knots = p_knots;
+        knot_policy = KnotPolicy::Custom;
+
+        return *this;
+      }
+
+      JointSplineBuilder & withKnotVector(const std::vector<double> & p_knots)
+      {
+        knots.resize(static_cast<Eigen::Index>(p_knots.size()));
+        for (std::size_t i = 0; i < p_knots.size(); ++i)
+        {
+          knots[static_cast<Eigen::Index>(i)] = p_knots[i];
+        }
+        knot_policy = KnotPolicy::Custom;
+
+        return *this;
+      }
+
+      JointSplineBuilder & withOpenUniformKnots(double p_min_q, double p_max_q)
+      {
+        knot_policy = KnotPolicy::OpenUniform;
+        min_q = p_min_q;
+        max_q = p_max_q;
+
+        return *this;
+      }
+
+      JointSplineBuilder & withUniformKnots(double p_min_q, double p_max_q)
+      {
+        knot_policy = KnotPolicy::Uniform;
+        min_q = p_min_q;
+        max_q = p_max_q;
+
+        return *this;
+      }
+
+      JointSpline build() const
+      {
+        Eigen::VectorXd joint_knots;
+
+        const size_t nCtrl = ctrlFrames.size();
+        switch (knot_policy)
+        {
+        case KnotPolicy::OpenUniform:
+          joint_knots = pinocchio::internal::generateOpenUniformKnots(min_q, max_q, nCtrl, degree);
+          break;
+
+        case KnotPolicy::Uniform:
+          joint_knots = pinocchio::internal::generateUniformKnots(min_q, max_q, nCtrl, degree);
+          break;
+
+        case KnotPolicy::Custom:
+          joint_knots = knots;
+          break;
+        default:
+          break;
+        }
+
+        return JointSpline(ctrlFrames, joint_knots, degree);
+      }
+
+    private:
+      std::vector<SE3> ctrlFrames;
+
+      size_t degree;
+
+      double min_q;
+      double max_q;
+
+      Eigen::VectorXd knots;
+      KnotPolicy knot_policy;
+    };
+
     // Forward declare
     struct JointComposite;
     struct JointMimic;
@@ -265,6 +400,7 @@ namespace pinocchio
       JointPlanar,
       JointHelical,
       JointUniversal,
+      JointSpline,
       boost::recursive_wrapper<JointComposite>,
       boost::recursive_wrapper<JointMimic>>
       JointVariant;
diff --git a/include/pinocchio/src/serialization/joints-data.hxx b/include/pinocchio/src/serialization/joints-data.hxx
index e40bcdae2d..691b7e9378 100644
--- a/include/pinocchio/src/serialization/joints-data.hxx
+++ b/include/pinocchio/src/serialization/joints-data.hxx
@@ -285,5 +285,17 @@ namespace boost
       fix::serialize(ar, static_cast<pinocchio::JointDataBase<JointType> &>(joint), version);
     }
 
+    template<class Archive, typename Scalar, int Options>
+    void serialize(
+      Archive & ar,
+      pinocchio::JointDataSplineTpl<Scalar, Options> & joint,
+      const unsigned int version)
+    {
+      ar & make_nvp("N", joint.N);
+
+      typedef pinocchio::JointDataSplineTpl<Scalar, Options> JointType;
+      fix::serialize(ar, static_cast<pinocchio::JointDataBase<JointType> &>(joint), version);
+    }
+
   } // namespace serialization
 } // namespace boost
diff --git a/include/pinocchio/src/serialization/joints-model.hxx b/include/pinocchio/src/serialization/joints-model.hxx
index b82b068a5a..f587093e23 100644
--- a/include/pinocchio/src/serialization/joints-model.hxx
+++ b/include/pinocchio/src/serialization/joints-model.hxx
@@ -331,5 +331,22 @@ namespace boost
       fix::serialize(ar, *static_cast<pinocchio::JointModelBase<JointType> *>(&joint), version);
     }
 
+    template<class Archive, typename Scalar, int Options>
+    void serialize(
+      Archive & ar,
+      pinocchio::JointModelSplineTpl<Scalar, Options> & joint,
+      const unsigned int version)
+    {
+      typedef pinocchio::JointModelSplineTpl<Scalar, Options> JointType;
+      ar & make_nvp("ctrlFrames", joint.ctrlFrames);
+      ar & make_nvp("degree", joint.degree);
+      ar & make_nvp("nbCtrlFrames", joint.nbCtrlFrames);
+      ar & make_nvp("relativeMotions", joint.relativeMotions);
+      ar & make_nvp("knots", joint.knots);
+      ar & make_nvp("min_q", joint.min_q);
+      ar & make_nvp("max_q", joint.max_q);
+
+      fix::serialize(ar, *static_cast<pinocchio::JointModelBase<JointType> *>(&joint), version);
+    }
   } // namespace serialization
 } // namespace boost
diff --git a/include/pinocchio/src/spatial/act-on-set.hxx b/include/pinocchio/src/spatial/act-on-set.hxx
index 5797182aa5..d9af94d115 100644
--- a/include/pinocchio/src/spatial/act-on-set.hxx
+++ b/include/pinocchio/src/spatial/act-on-set.hxx
@@ -575,7 +575,7 @@ namespace pinocchio
         EIGEN_STATIC_ASSERT_VECTOR_ONLY(Mat);
         EIGEN_STATIC_ASSERT_VECTOR_ONLY(MatRet);
 
-        typedef MotionRef<Mat> MotionRefOnMat;
+        typedef MotionRef<const Mat> MotionRefOnMat;
         typedef MotionRef<MatRet> MotionRefOnMatRet;
 
         MotionRefOnMat min(iV.derived());
diff --git a/include/pinocchio/src/spatial/log.hxx b/include/pinocchio/src/spatial/log.hxx
index d5701ac0d1..be5bd3c1bb 100644
--- a/include/pinocchio/src/spatial/log.hxx
+++ b/include/pinocchio/src/spatial/log.hxx
@@ -29,9 +29,12 @@ namespace pinocchio
       static const long i1 = (i0 + 1) % 3;
       static const long i2 = (i0 + 2) % 3;
       Vector3 & axis = _axis.const_cast_derived();
-
+      // boost::multiprecision return an expression that can't
+      // be a math::sqrt input.
+      // To overcome this issue we force the expression to be evaluated in a Scalar.
+      Scalar sqrt_input = val + eps + eps * eps;
       const Scalar s =
-        math::sqrt(val + eps + eps * eps)
+        math::sqrt(sqrt_input)
         * if_then_else(
           GE, R.coeff(i2, i1), R.coeff(i1, i2), Scalar(1.),
           Scalar(-1.)); // Ensure value in sqrt is non negative and that s is non zero
diff --git a/models/biomechanics/femur_r.stl b/models/biomechanics/femur_r.stl
new file mode 100644
index 0000000000..d93b3d409f
Binary files /dev/null and b/models/biomechanics/femur_r.stl differ
diff --git a/models/biomechanics/tibia_r.stl b/models/biomechanics/tibia_r.stl
new file mode 100644
index 0000000000..1af887b342
Binary files /dev/null and b/models/biomechanics/tibia_r.stl differ
diff --git a/sources.cmake b/sources.cmake
index e9ff793a62..e3510c08a9 100644
--- a/sources.cmake
+++ b/sources.cmake
@@ -362,6 +362,7 @@ set(${PROJECT_NAME}_CORE_PRIVATE_HEADERS
     ${PROJECT_SOURCE_DIR}/include/pinocchio/src/multibody/joint/joint-translation.hxx
     ${PROJECT_SOURCE_DIR}/include/pinocchio/src/multibody/joint/joint-prismatic-unaligned.hxx
     ${PROJECT_SOURCE_DIR}/include/pinocchio/src/multibody/joint/joint-mimic.hxx
+    ${PROJECT_SOURCE_DIR}/include/pinocchio/src/multibody/joint/joint-spline.hxx
     ${PROJECT_SOURCE_DIR}/include/pinocchio/src/multibody/joint/fwd.hxx
     ${PROJECT_SOURCE_DIR}/include/pinocchio/src/multibody/joint/joint-free-flyer.hxx
     ${PROJECT_SOURCE_DIR}/include/pinocchio/src/multibody/joint/joint-common-operations.hxx
@@ -373,6 +374,7 @@ set(${PROJECT_NAME}_CORE_PRIVATE_HEADERS
     ${PROJECT_SOURCE_DIR}/include/pinocchio/src/multibody/joint/joint-data-base.hxx
     ${PROJECT_SOURCE_DIR}/include/pinocchio/src/multibody/joint/joint-basic-visitors.hxx
     ${PROJECT_SOURCE_DIR}/include/pinocchio/src/multibody/joint/joint-model-base.hxx
+    ${PROJECT_SOURCE_DIR}/include/pinocchio/src/multibody/joint/spline-utils.hxx
     ${PROJECT_SOURCE_DIR}/include/pinocchio/src/constraints/constraint-data-base.hxx
     ${PROJECT_SOURCE_DIR}/include/pinocchio/src/constraints/kinematics-constraint-model-base.hxx
     ${PROJECT_SOURCE_DIR}/include/pinocchio/src/constraints/baumgarte-corrector-parameters.hxx
@@ -549,6 +551,7 @@ set(${PROJECT_NAME}_CASADI_PRIVATE_HEADERS
     ${PROJECT_SOURCE_DIR}/include/pinocchio/src/autodiff/casadi/math/triangular-matrix.hxx
     ${PROJECT_SOURCE_DIR}/include/pinocchio/src/autodiff/casadi/spatial/se3-tpl.hxx
     ${PROJECT_SOURCE_DIR}/include/pinocchio/src/autodiff/casadi/utils/static-if.hxx
+    ${PROJECT_SOURCE_DIR}/include/pinocchio/src/autodiff/casadi/multibody/joint/spline-utils.hxx
 )
 
 # CPPAD
diff --git a/src/parsers/graph/model-graph-algo.cpp b/src/parsers/graph/model-graph-algo.cpp
index bb21508b7d..6a43786116 100644
--- a/src/parsers/graph/model-graph-algo.cpp
+++ b/src/parsers/graph/model-graph-algo.cpp
@@ -127,6 +127,9 @@ namespace pinocchio
         // Joint Universal
         typedef typename JointCollectionDefault::JointModelUniversal JointModelUniversal;
 
+        // JointSpline
+        typedef typename JointCollectionDefault::JointModelSpline JointModelSpline;
+
         typedef JointModel ReturnType;
 
         ReturnType operator()(const JointFixed & /*joint*/) const
@@ -243,6 +246,10 @@ namespace pinocchio
         {
           return boost::apply_visitor(*this, joint.secondary_joint);
         }
+        ReturnType operator()(const JointSpline & joint) const
+        {
+          return JointModelSpline(joint.ctrlFrames, joint.knots, joint.degree);
+        }
         ReturnType operator()(const JointComposite & joint) const
         {
           JointModelComposite jmodel;
@@ -372,6 +379,17 @@ namespace pinocchio
           addJointBetweenBodies(joint, b_f);
         }
 
+        void operator()(const JointSpline & joint, const BodyFrame & b_f)
+        {
+          if (!edge.forward)
+            PINOCCHIO_THROW_PRETTY(
+              std::invalid_argument, "Graph - JointSpline cannot be reversed.");
+          // The spline joint cannot be reversed because it's not trivial to generate
+          // a spline that will generate the inverse transform of the forward spline
+
+          addJointBetweenBodies(joint, b_f);
+        }
+
         void operator()(const JointFixed & joint, const BodyFrame & b_f)
         {
           // Need to check what's vertex the edge is coming from. If it's a body, then we add
diff --git a/src/parsers/graph/model-graph.cpp b/src/parsers/graph/model-graph.cpp
index 7e856140e0..c65bdc5ca8 100644
--- a/src/parsers/graph/model-graph.cpp
+++ b/src/parsers/graph/model-graph.cpp
@@ -96,6 +96,10 @@ namespace pinocchio
         {
           return {joint, SE3::Identity()};
         }
+        ReturnType operator()(const JointSpline & joint) const
+        {
+          return {joint, SE3::Identity()};
+        }
         ReturnType operator()(const JointComposite & joint) const
         {
           JointComposite jReturn;
diff --git a/unittest/CMakeLists.txt b/unittest/CMakeLists.txt
index 4ae98a10c2..2de25e3cab 100644
--- a/unittest/CMakeLists.txt
+++ b/unittest/CMakeLists.txt
@@ -391,6 +391,7 @@ add_pinocchio_unit_test(joint-mimic)
 add_pinocchio_unit_test(joint-helical)
 add_pinocchio_unit_test(joint-universal)
 add_pinocchio_unit_test(joint-visitors)
+add_pinocchio_unit_test(joint-spline)
 
 # Main corpus
 add_pinocchio_unit_test(model COLLISION_OPTIONAL)
diff --git a/unittest/all-joints.cpp b/unittest/all-joints.cpp
index bcd9e4386d..08c40642c9 100644
--- a/unittest/all-joints.cpp
+++ b/unittest/all-joints.cpp
@@ -3,6 +3,8 @@
 // Copyright(c) 2015 Wandercraft, 86 rue de Paris 91400 Orsay, France.
 //
 
+#include "utils/joints-init.hpp"
+
 #include "pinocchio/multibody/joint.hpp"
 #include "pinocchio/multibody/liegroup.hpp"
 
@@ -11,172 +13,6 @@
 
 using namespace pinocchio;
 
-template<typename JointModel_>
-struct init;
-
-template<typename JointModel_>
-struct init
-{
-  static JointModel_ run()
-  {
-    JointModel_ jmodel;
-    jmodel.setIndexes(0, 0, 0);
-    return jmodel;
-  }
-};
-
-template<typename Scalar, int Options>
-struct init<pinocchio::JointModelRevoluteUnalignedTpl<Scalar, Options>>
-{
-  typedef pinocchio::JointModelRevoluteUnalignedTpl<Scalar, Options> JointModel;
-
-  static JointModel run()
-  {
-    typedef typename JointModel::Vector3 Vector3;
-    JointModel jmodel(Vector3::Random().normalized());
-
-    jmodel.setIndexes(0, 0, 0);
-    return jmodel;
-  }
-};
-
-template<typename Scalar, int Options>
-struct init<pinocchio::JointModelRevoluteUnboundedUnalignedTpl<Scalar, Options>>
-{
-  typedef pinocchio::JointModelRevoluteUnboundedUnalignedTpl<Scalar, Options> JointModel;
-
-  static JointModel run()
-  {
-    typedef typename JointModel::Vector3 Vector3;
-    JointModel jmodel(Vector3::Random().normalized());
-
-    jmodel.setIndexes(0, 0, 0);
-    return jmodel;
-  }
-};
-
-template<typename Scalar, int Options>
-struct init<pinocchio::JointModelPrismaticUnalignedTpl<Scalar, Options>>
-{
-  typedef pinocchio::JointModelPrismaticUnalignedTpl<Scalar, Options> JointModel;
-
-  static JointModel run()
-  {
-    typedef typename JointModel::Vector3 Vector3;
-    JointModel jmodel(Vector3::Random().normalized());
-
-    jmodel.setIndexes(0, 0, 0);
-    return jmodel;
-  }
-};
-
-template<typename Scalar, int Options, template<typename, int> class JointCollection>
-struct init<pinocchio::JointModelTpl<Scalar, Options, JointCollection>>
-{
-  typedef pinocchio::JointModelTpl<Scalar, Options, JointCollection> JointModel;
-
-  static JointModel run()
-  {
-    typedef pinocchio::JointModelRevoluteTpl<Scalar, Options, 0> JointModelRX;
-    JointModel jmodel((JointModelRX()));
-
-    jmodel.setIndexes(0, 0, 0);
-    return jmodel;
-  }
-};
-
-template<typename Scalar, int Options, template<typename, int> class JointCollection>
-struct init<pinocchio::JointModelCompositeTpl<Scalar, Options, JointCollection>>
-{
-  typedef pinocchio::JointModelCompositeTpl<Scalar, Options, JointCollection> JointModel;
-
-  static JointModel run()
-  {
-    typedef pinocchio::JointModelRevoluteTpl<Scalar, Options, 0> JointModelRX;
-    typedef pinocchio::JointModelRevoluteTpl<Scalar, Options, 1> JointModelRY;
-    JointModel jmodel((JointModelRX()));
-    jmodel.addJoint(JointModelRY());
-
-    jmodel.setIndexes(0, 0, 0);
-    return jmodel;
-  }
-};
-
-template<typename Scalar, int Options, template<typename, int> class JointCollection>
-struct init<pinocchio::JointModelMimicTpl<Scalar, Options, JointCollection>>
-{
-  typedef pinocchio::JointModelMimicTpl<Scalar, Options, JointCollection> JointModel;
-
-  static JointModel run()
-  {
-
-    typedef pinocchio::JointModelRevoluteTpl<Scalar, Options, 0> JointModelRX;
-    JointModelRX jmodel_ref = init<JointModelRX>::run();
-
-    JointModel jmodel(jmodel_ref, 1., 0.);
-    jmodel.setIndexes(1, 0, 0, 0);
-    return jmodel;
-  }
-};
-
-template<typename Scalar, int Options>
-struct init<pinocchio::JointModelUniversalTpl<Scalar, Options>>
-{
-  typedef pinocchio::JointModelUniversalTpl<Scalar, Options> JointModel;
-
-  static JointModel run()
-  {
-    JointModel jmodel(XAxis::vector(), YAxis::vector());
-
-    jmodel.setIndexes(0, 0, 0);
-    return jmodel;
-  }
-};
-
-template<typename Scalar, int Options, int axis>
-struct init<pinocchio::JointModelHelicalTpl<Scalar, Options, axis>>
-{
-  typedef pinocchio::JointModelHelicalTpl<Scalar, Options, axis> JointModel;
-
-  static JointModel run()
-  {
-    JointModel jmodel(static_cast<Scalar>(0.5));
-
-    jmodel.setIndexes(0, 0, 0);
-    return jmodel;
-  }
-};
-
-template<typename Scalar, int Options>
-struct init<pinocchio::JointModelEllipsoidTpl<Scalar, Options>>
-{
-  typedef pinocchio::JointModelEllipsoidTpl<Scalar, Options> JointModel;
-
-  static JointModel run()
-  {
-    JointModel jmodel(
-      static_cast<Scalar>(0.01), static_cast<Scalar>(0.02), static_cast<Scalar>(0.03));
-
-    jmodel.setIndexes(0, 0, 0);
-    return jmodel;
-  }
-};
-
-template<typename Scalar, int Options>
-struct init<pinocchio::JointModelHelicalUnalignedTpl<Scalar, Options>>
-{
-  typedef pinocchio::JointModelHelicalUnalignedTpl<Scalar, Options> JointModel;
-
-  static JointModel run()
-  {
-    typedef typename JointModel::Vector3 Vector3;
-    JointModel jmodel(Vector3::Random().normalized());
-
-    jmodel.setIndexes(0, 0, 0);
-    return jmodel;
-  }
-};
-
 BOOST_AUTO_TEST_SUITE(joint_model_base_test)
 
 template<typename TestDerived>
@@ -185,7 +21,7 @@ struct TestJointModel
   template<typename JointModel>
   void operator()(const pinocchio::JointModelBase<JointModel> &) const
   {
-    JointModel jmodel = init<JointModel>::run();
+    JointModel jmodel = init<JointModel>::run().jmodel;
     return TestDerived::test(jmodel);
   }
 };
@@ -270,9 +106,6 @@ struct TestJointModelCast : TestJointModel<TestJointModelCast>
     BOOST_CHECK(jmodel == jmodel);
     BOOST_CHECK(jmodel.template cast<Scalar>().isEqual(jmodel));
     BOOST_CHECK(jmodel.template cast<Scalar>() == jmodel);
-    BOOST_CHECK_MESSAGE(
-      jmodel.template cast<long double>().template cast<double>() == jmodel,
-      std::string("Error when casting " + jmodel.shortname() + " from long double to double."));
   }
 };
 
diff --git a/unittest/casadi/joints.cpp b/unittest/casadi/joints.cpp
index 515d55f2e8..7f7ba74f59 100644
--- a/unittest/casadi/joints.cpp
+++ b/unittest/casadi/joints.cpp
@@ -11,6 +11,8 @@
 #include <boost/test/unit_test.hpp>
 #include <boost/utility/binary.hpp>
 
+#include "../utils/joints-init.hpp"
+
 BOOST_AUTO_TEST_SUITE(BOOST_TEST_MODULE)
 
 BOOST_AUTO_TEST_CASE(test_jointRX_motion_space)
@@ -142,255 +144,28 @@ BOOST_AUTO_TEST_CASE(test_jointRX_motion_space)
   }
 }
 
-template<typename JointModel_>
-struct init;
-
-template<typename JointModel_>
-struct init
-{
-  static JointModel_ run()
-  {
-    JointModel_ jmodel;
-    jmodel.setIndexes(0, 0, 0);
-    return jmodel;
-  }
-
-  static std::string name()
-  {
-    return "default " + JointModel_::classname();
-  }
-};
-
-template<typename Scalar, int Options>
-struct init<pinocchio::JointModelRevoluteUnalignedTpl<Scalar, Options>>
-{
-  typedef pinocchio::JointModelRevoluteUnalignedTpl<Scalar, Options> JointModel;
-
-  static JointModel run()
-  {
-    typedef typename JointModel::Vector3 Vector3;
-    JointModel jmodel(Vector3::Random().normalized());
-
-    jmodel.setIndexes(0, 0, 0);
-    return jmodel;
-  }
-
-  static std::string name()
-  {
-    return JointModel::classname();
-  }
-};
-
-template<typename Scalar, int Options>
-struct init<pinocchio::JointModelRevoluteUnboundedUnalignedTpl<Scalar, Options>>
-{
-  typedef pinocchio::JointModelRevoluteUnboundedUnalignedTpl<Scalar, Options> JointModel;
-
-  static JointModel run()
-  {
-    typedef typename JointModel::Vector3 Vector3;
-    JointModel jmodel(Vector3::Random().normalized());
-
-    jmodel.setIndexes(0, 0, 0);
-    return jmodel;
-  }
-
-  static std::string name()
-  {
-    return JointModel::classname();
-  }
-};
-
-template<typename Scalar, int Options>
-struct init<pinocchio::JointModelPrismaticUnalignedTpl<Scalar, Options>>
-{
-  typedef pinocchio::JointModelPrismaticUnalignedTpl<Scalar, Options> JointModel;
-
-  static JointModel run()
-  {
-    typedef typename JointModel::Vector3 Vector3;
-    JointModel jmodel(Vector3::Random().normalized());
-
-    jmodel.setIndexes(0, 0, 0);
-    return jmodel;
-  }
-
-  static std::string name()
-  {
-    return JointModel::classname();
-  }
-};
-
-template<typename Scalar, int Options, template<typename, int> class JointCollection>
-struct init<pinocchio::JointModelTpl<Scalar, Options, JointCollection>>
-{
-  typedef pinocchio::JointModelTpl<Scalar, Options, JointCollection> JointModel;
-
-  static JointModel run()
-  {
-    typedef pinocchio::JointModelRevoluteTpl<Scalar, Options, 0> JointModelRX;
-    JointModel jmodel((JointModelRX()));
-
-    jmodel.setIndexes(0, 0, 0);
-    return jmodel;
-  }
-
-  static std::string name()
-  {
-    return JointModel::classname();
-  }
-};
-
-template<typename Scalar, int Options, template<typename, int> class JointCollection>
-struct init<pinocchio::JointModelCompositeTpl<Scalar, Options, JointCollection>>
-{
-  typedef pinocchio::JointModelCompositeTpl<Scalar, Options, JointCollection> JointModel;
-
-  static JointModel run()
-  {
-    typedef pinocchio::JointModelRevoluteTpl<Scalar, Options, 0> JointModelRX;
-    typedef pinocchio::JointModelRevoluteTpl<Scalar, Options, 1> JointModelRY;
-    JointModel jmodel((JointModelRX()));
-    jmodel.addJoint(JointModelRY());
-
-    jmodel.setIndexes(0, 0, 0);
-    return jmodel;
-  }
-
-  static std::string name()
-  {
-    return JointModel::classname();
-  }
-};
-
-template<typename Scalar, int Options, template<typename, int> class JointCollection>
-struct init<pinocchio::JointModelMimicTpl<Scalar, Options, JointCollection>>
-{
-  typedef pinocchio::JointModelMimicTpl<Scalar, Options, JointCollection> JointModel;
-
-  static JointModel run()
-  {
-    typedef pinocchio::JointModelRevoluteTpl<Scalar, Options, 0> JointModelRX;
-    JointModelRX jmodel_ref = init<JointModelRX>::run();
-
-    JointModel jmodel(jmodel_ref, 1., 0.);
-    jmodel.setIndexes(1, 0, 0, 0);
-    return jmodel;
-  }
-
-  static std::string name()
-  {
-    return JointModel::classname();
-  }
-};
-
 struct TestADOnJoints
 {
   template<typename JointModel_>
   void operator()(const pinocchio::JointModelBase<JointModel_> &) const
   {
-    JointModel_ jmodel = init<JointModel_>::run();
-    test(jmodel);
-  }
+    auto jmodelWithParams = pinocchio::init<JointModel_>::run();
 
-  template<typename Scalar, int Options, int axis>
-  void operator()(const pinocchio::JointModelHelicalTpl<Scalar, Options, axis> &) const
-  {
-    typedef pinocchio::JointModelHelicalTpl<Scalar, Options, axis> JointModel;
-    JointModel jmodel(Scalar(0.4));
-    jmodel.setIndexes(0, 0, 0);
-
-    test(jmodel);
-  }
-
-  template<typename Scalar, int Options>
-  void operator()(const pinocchio::JointModelUniversalTpl<Scalar, Options> &) const
-  {
-    typedef pinocchio::JointModelUniversalTpl<Scalar, Options> JointModel;
-    typedef typename JointModel::Vector3 Vector3;
-    JointModel jmodel(Vector3::UnitX(), Vector3::UnitY());
-    jmodel.setIndexes(0, 0, 0);
-
-    test(jmodel);
-  }
-
-  template<typename Scalar, int Options>
-  void operator()(const pinocchio::JointModelHelicalUnalignedTpl<Scalar, Options> &) const
-  {
-    typedef pinocchio::JointModelHelicalUnalignedTpl<Scalar, Options> JointModel;
-    typedef typename JointModel::Vector3 Vector3;
-    JointModel jmodel(Vector3::Random().normalized());
-    jmodel.setIndexes(0, 0, 0);
-
-    test(jmodel);
-  }
-
-  template<typename Scalar, int Options>
-  void operator()(const pinocchio::JointModelRevoluteUnalignedTpl<Scalar, Options> &) const
-  {
-    typedef pinocchio::JointModelRevoluteUnalignedTpl<Scalar, Options> JointModel;
-    typedef typename JointModel::Vector3 Vector3;
-    JointModel jmodel(Vector3::Random().normalized());
-    jmodel.setIndexes(0, 0, 0);
-
-    test(jmodel);
-  }
-
-  template<typename Scalar, int Options>
-  void operator()(const pinocchio::JointModelRevoluteUnboundedUnalignedTpl<Scalar, Options> &) const
-  {
-    typedef pinocchio::JointModelRevoluteUnboundedUnalignedTpl<Scalar, Options> JointModel;
-    typedef typename JointModel::Vector3 Vector3;
-    JointModel jmodel(Vector3::Random().normalized());
-    jmodel.setIndexes(0, 0, 0);
-
-    test(jmodel);
-  }
-
-  template<typename Scalar, int Options>
-  void operator()(const pinocchio::JointModelPrismaticUnalignedTpl<Scalar, Options> &) const
-  {
-    typedef pinocchio::JointModelPrismaticUnalignedTpl<Scalar, Options> JointModel;
-    typedef typename JointModel::Vector3 Vector3;
-    JointModel jmodel(Vector3::Random().normalized());
-    jmodel.setIndexes(0, 0, 0);
-
-    test(jmodel);
-  }
-
-  template<typename Scalar, int Options, template<typename, int> class JointCollection>
-  void operator()(const pinocchio::JointModelTpl<Scalar, Options, JointCollection> &) const
-  {
-    typedef pinocchio::JointModelRevoluteTpl<Scalar, Options, 0> JointModelRX;
-    typedef pinocchio::JointModelTpl<Scalar, Options, JointCollection> JointModel;
-    JointModel jmodel((JointModelRX()));
-    jmodel.setIndexes(0, 0, 0);
-
-    test(jmodel);
-  }
-
-  template<typename Scalar, int Options, template<typename, int> class JointCollection>
-  void operator()(const pinocchio::JointModelCompositeTpl<Scalar, Options, JointCollection> &) const
-  {
-    typedef pinocchio::JointModelRevoluteTpl<Scalar, Options, 0> JointModelRX;
-    typedef pinocchio::JointModelRevoluteTpl<Scalar, Options, 1> JointModelRY;
-    typedef pinocchio::JointModelCompositeTpl<Scalar, Options, JointCollection> JointModel;
-    JointModel jmodel((JointModelRX()));
-    jmodel.addJoint(JointModelRY());
-    jmodel.setIndexes(0, 0, 0);
-
-    test(jmodel);
+    test(jmodelWithParams.jmodel, jmodelWithParams.lb, jmodelWithParams.ub);
   }
 
   // TODO: get the nq and nv quantity from LieGroups
   template<typename Scalar, int Options, template<typename, int> class JointCollection>
-  static void
-  test(const pinocchio::JointModelMimicTpl<Scalar, Options, JointCollection> & /*jmodel*/)
+  void operator()(
+    const pinocchio::JointModelMimicTpl<Scalar, Options, JointCollection> & /*jmodel*/) const
   { /* do nothing */
   }
 
   template<typename JointModel>
-  static void test(const pinocchio::JointModelBase<JointModel> & jmodel)
+  static void test(
+    const pinocchio::JointModelBase<JointModel> & jmodel,
+    const typename JointModel::ConfigVector_t & lb,
+    const typename JointModel::ConfigVector_t & ub)
   {
     std::cout << "--" << std::endl;
     std::cout << "jmodel: " << jmodel.shortname() << std::endl;
@@ -424,9 +199,6 @@ struct TestADOnJoints
 
     ConfigVector q(jmodel.nq());
 
-    ConfigVector lb(ConfigVector::Constant(jmodel.nq(), -1.));
-    ConfigVector ub(ConfigVector::Constant(jmodel.nq(), 1.));
-
     typedef pinocchio::RandomConfigurationStep<
       pinocchio::LieGroupMap, ConfigVector, ConfigVector, ConfigVector>
       RandomConfigAlgo;
diff --git a/unittest/cppad/joints.cpp b/unittest/cppad/joints.cpp
index df04a2fe84..08deda9a8a 100644
--- a/unittest/cppad/joints.cpp
+++ b/unittest/cppad/joints.cpp
@@ -186,6 +186,27 @@ struct TestADOnJoints
     test(jmodel);
   }
 
+  template<typename Scalar, int Options>
+  void operator()(const pinocchio::JointModelSplineTpl<Scalar, Options> &) const
+  {
+    typedef pinocchio::JointModelSplineTpl<Scalar, Options> JointModel;
+    typedef pinocchio::SE3Tpl<Scalar> SE3;
+    typedef pinocchio::JointModelSplineBuilderTpl<Scalar, Options> JointModelSplineBuilder;
+
+    std::vector<SE3> ctrlFrames;
+    for (int k = 0; k < 5; k++)
+      ctrlFrames.push_back(SE3::Random());
+
+    JointModel jmodel = JointModelSplineBuilder()
+                          .withControlFrameVector(ctrlFrames)
+                          .withDegree(3)
+                          .withUniformKnots(0., 1.)
+                          .build();
+    jmodel.setIndexes(0, 0, 0);
+
+    test(jmodel);
+  }
+
   // TODO: implement it
   template<typename Scalar, int Options, template<typename, int> class JointCollection>
   void operator()(
@@ -207,6 +228,101 @@ struct TestADOnJoints
     test(jmodel);
   }
 
+  template<typename Scalar, int Options>
+  static void test(const pinocchio::JointModelSplineTpl<Scalar, Options> & jmodel)
+  {
+    using CppAD::AD;
+    using CppAD::NearEqual;
+
+    typedef AD<Scalar> AD_scalar;
+
+    typedef Eigen::Matrix<AD_scalar, Eigen::Dynamic, 1> VectorXAD;
+    typedef Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic> MatrixX;
+
+    typedef pinocchio::SE3Tpl<AD_scalar> SE3AD;
+    typedef pinocchio::MotionTpl<AD_scalar> MotionAD;
+    typedef pinocchio::SE3Tpl<double> SE3;
+    typedef pinocchio::MotionTpl<double> Motion;
+    typedef pinocchio::JointMotionSubspaceTpl<Eigen::Dynamic, double> JointMotionSubspaceXd;
+
+    typedef Eigen::Matrix<AD_scalar, Eigen::Dynamic, 1> VectorXAD;
+
+    typedef pinocchio::JointModelSplineTpl<Scalar, Options> JointModel;
+    typedef typename pinocchio::CastType<AD_scalar, JointModel>::type JointModelAD;
+    typedef typename JointModelAD::JointDataDerived JointDataAD;
+
+    typedef typename JointModelAD::ConfigVector_t ConfigVectorAD;
+
+    typedef typename JointModel::JointDataDerived JointData;
+    typedef typename JointModel::ConfigVector_t ConfigVector;
+    typedef typename JointModel::TangentVector_t TangentVector;
+
+    JointData jdata(jmodel.createData());
+    pinocchio::JointDataBase<JointData> & jdata_base = jdata;
+
+    JointModelAD jmodel_ad = jmodel.template cast<AD_scalar>();
+    JointDataAD jdata_ad(jmodel_ad.createData());
+    pinocchio::JointDataBase<JointDataAD> & jdata_ad_base = jdata_ad;
+
+    ConfigVector q(jmodel.nq());
+    ConfigVector lb(ConfigVector::Constant(jmodel.nq(), 0));
+    ConfigVector ub(ConfigVector::Constant(jmodel.nq(), 1.));
+
+    typedef pinocchio::RandomConfigurationStep<
+      pinocchio::LieGroupMap, ConfigVector, ConfigVector, ConfigVector>
+      RandomConfigAlgo;
+    RandomConfigAlgo::run(jmodel.derived(), typename RandomConfigAlgo::ArgsType(q, lb, ub));
+
+    ConfigVectorAD q_ad(q.template cast<AD_scalar>());
+
+    // Zero order
+    jmodel_ad.calc(jdata_ad, q_ad);
+    jmodel.calc(jdata, q);
+
+    SE3 M1(jdata_base.M());
+    SE3AD M2(jdata_ad_base.M());
+    BOOST_CHECK(M1.isApprox(M2.template cast<Scalar>()));
+
+    // First order
+    TangentVector v(TangentVector::Random(jmodel.nv()));
+    VectorXAD X(jmodel_ad.nv());
+
+    for (Eigen::DenseIndex k = 0; k < jmodel.nv(); ++k)
+    {
+      X[k] = v[k];
+    }
+    CppAD::Independent(X);
+    jmodel_ad.calc(jdata_ad, q_ad, X);
+    jmodel.calc(jdata, q, v);
+    VectorXAD Y(6);
+    MotionAD m_ad(jdata_ad_base.v());
+    Motion m(jdata_base.v());
+    JointMotionSubspaceXd Sref(jdata_base.S().matrix());
+
+    for (Eigen::DenseIndex k = 0; k < 3; ++k)
+    {
+      Y[k + Motion::LINEAR] = m_ad.linear()[k];
+      Y[k + Motion::ANGULAR] = m_ad.angular()[k];
+    }
+
+    CppAD::ADFun<Scalar> vjoint(X, Y);
+
+    CPPAD_TESTVECTOR(Scalar) x((size_t)jmodel_ad.nv());
+    for (Eigen::DenseIndex k = 0; k < jmodel.nv(); ++k)
+    {
+      x[(size_t)k] = v[k];
+    }
+
+    CPPAD_TESTVECTOR(Scalar) jac = vjoint.Jacobian(x);
+    MatrixX S(6, jac.size() / 6);
+    S = Eigen::Map<typename PINOCCHIO_EIGEN_PLAIN_ROW_MAJOR_TYPE(MatrixX)>(
+      jac.data(), S.rows(), S.cols());
+
+    BOOST_CHECK(m.isApprox(m_ad.template cast<Scalar>()));
+
+    BOOST_CHECK(Sref.matrix().isApprox(S));
+  }
+
   template<typename JointModel>
   static void test(const pinocchio::JointModelBase<JointModel> & jmodel)
   {
diff --git a/unittest/cppadcg/basic.cpp b/unittest/cppadcg/basic.cpp
index cde394a510..78d057285f 100644
--- a/unittest/cppadcg/basic.cpp
+++ b/unittest/cppadcg/basic.cpp
@@ -101,15 +101,17 @@ BOOST_AUTO_TEST_CASE(test_eigen_support)
 
 BOOST_AUTO_TEST_CASE(test_cast)
 {
-  typedef CppAD::cg::CG<double> CGScalar;
-  typedef CppAD::AD<double> ADScalar;
-  typedef CppAD::AD<float> ADFloat;
-  typedef pinocchio::ModelTpl<CGScalar> CGModel;
+  using namespace CppAD;
+  using namespace CppAD::cg;
+
+  typedef CG<double> CGD;
+  typedef AD<CGD> ADCG;
+  typedef pinocchio::ModelTpl<ADCG> ADCGModel;
 
   pinocchio::SE3 M(pinocchio::SE3::Random());
-  typedef pinocchio::SE3Tpl<CGScalar> CGSE3;
+  typedef pinocchio::SE3Tpl<ADCG> CGSE3;
 
-  CGSE3 cg_M = M.cast<CGScalar>();
+  CGSE3 cg_M = M.cast<ADCG>();
   BOOST_CHECK(cg_M.cast<double>().isApprox(M));
 
   pinocchio::SE3::Vector3 axis(1., 1., 1.);
@@ -117,34 +119,15 @@ BOOST_AUTO_TEST_CASE(test_cast)
   BOOST_CHECK(axis.isUnitary());
 
   pinocchio::JointModelPrismaticUnaligned jmodel_prismatic(axis);
-  typedef pinocchio::JointModelPrismaticUnalignedTpl<CGScalar> CGJointModelPrismaticUnaligned;
-
-  CGScalar cg_value;
-  cg_value = -1.;
-  ADScalar ad_value;
-  ad_value = -1.;
-  ADFloat ad_float;
-  ad_float = -1.;
-  abs(ad_value);
-  abs(ad_float);
-  abs(cg_value);
+  typedef pinocchio::JointModelPrismaticUnalignedTpl<ADCG> ADCGJointModelPrismaticUnaligned;
 
-  CPPAD_TESTVECTOR(ADScalar) ad_x(3);
-  CGJointModelPrismaticUnaligned cg_jmodel_prismatic(axis.cast<CGScalar>());
+  CPPAD_TESTVECTOR(ADCG) ad_x(3);
+  ADCGJointModelPrismaticUnaligned cg_jmodel_prismatic(axis.cast<ADCG>());
 
   pinocchio::Model model;
   pinocchio::buildModels::humanoidRandom(model);
 
-  CGModel cg_model = model.cast<CGScalar>();
-
-  {
-    CppAD::AD<double> ad_value(-1.);
-    abs(ad_value); // works perfectly
-
-    CppAD::cg::CG<double> cg_value(-1.);
-    abs(cg_value); // does not compile because abs<double>(const CppAD::cg::CG<double>&) is defined
-                   // in namespace CppAD and not CppAD::cg
-  }
+  ADCGModel cg_model = model.cast<ADCG>();
 }
 
 BOOST_AUTO_TEST_CASE(test_dynamic_link)
diff --git a/unittest/finite-differences.cpp b/unittest/finite-differences.cpp
index 9e0770e2cf..7a20b4cbe3 100644
--- a/unittest/finite-differences.cpp
+++ b/unittest/finite-differences.cpp
@@ -12,6 +12,8 @@
 #include <boost/test/unit_test.hpp>
 #include <boost/utility/binary.hpp>
 
+#include "utils/joints-init.hpp"
+
 using namespace pinocchio;
 using namespace Eigen;
 
@@ -59,171 +61,6 @@ void filterValue(MatrixBase<Matrix> & mat, typename Matrix::Scalar value)
       math::fabs(mat.derived().data()[k]) <= value ? 0 : mat.derived().data()[k];
 }
 
-template<typename JointModel_>
-struct init;
-
-template<typename JointModel_>
-struct init
-{
-  static JointModel_ run()
-  {
-    JointModel_ jmodel;
-    jmodel.setIndexes(0, 0, 0);
-    return jmodel;
-  }
-};
-
-template<typename Scalar, int Options>
-struct init<pinocchio::JointModelRevoluteUnalignedTpl<Scalar, Options>>
-{
-  typedef pinocchio::JointModelRevoluteUnalignedTpl<Scalar, Options> JointModel;
-
-  static JointModel run()
-  {
-    typedef typename JointModel::Vector3 Vector3;
-    JointModel jmodel(Vector3::Random().normalized());
-
-    jmodel.setIndexes(0, 0, 0);
-    return jmodel;
-  }
-};
-
-template<typename Scalar, int Options>
-struct init<pinocchio::JointModelRevoluteUnboundedUnalignedTpl<Scalar, Options>>
-{
-  typedef pinocchio::JointModelRevoluteUnboundedUnalignedTpl<Scalar, Options> JointModel;
-
-  static JointModel run()
-  {
-    typedef typename JointModel::Vector3 Vector3;
-    JointModel jmodel(Vector3::Random().normalized());
-
-    jmodel.setIndexes(0, 0, 0);
-    return jmodel;
-  }
-};
-
-template<typename Scalar, int Options>
-struct init<pinocchio::JointModelPrismaticUnalignedTpl<Scalar, Options>>
-{
-  typedef pinocchio::JointModelPrismaticUnalignedTpl<Scalar, Options> JointModel;
-
-  static JointModel run()
-  {
-    typedef typename JointModel::Vector3 Vector3;
-    JointModel jmodel(Vector3::Random().normalized());
-
-    jmodel.setIndexes(0, 0, 0);
-    return jmodel;
-  }
-};
-
-template<typename Scalar, int Options>
-struct init<pinocchio::JointModelHelicalUnalignedTpl<Scalar, Options>>
-{
-  typedef pinocchio::JointModelHelicalUnalignedTpl<Scalar, Options> JointModel;
-
-  static JointModel run()
-  {
-    typedef typename JointModel::Vector3 Vector3;
-    JointModel jmodel(Vector3::Random().normalized());
-
-    jmodel.setIndexes(0, 0, 0);
-    return jmodel;
-  }
-};
-
-template<typename Scalar, int Options>
-struct init<pinocchio::JointModelUniversalTpl<Scalar, Options>>
-{
-  typedef pinocchio::JointModelUniversalTpl<Scalar, Options> JointModel;
-
-  static JointModel run()
-  {
-    JointModel jmodel(XAxis::vector(), YAxis::vector());
-
-    jmodel.setIndexes(0, 0, 0);
-    return jmodel;
-  }
-};
-
-template<typename Scalar, int Options, int axis>
-struct init<pinocchio::JointModelHelicalTpl<Scalar, Options, axis>>
-{
-  typedef pinocchio::JointModelHelicalTpl<Scalar, Options, axis> JointModel;
-
-  static JointModel run()
-  {
-    JointModel jmodel(static_cast<Scalar>(0.5));
-
-    jmodel.setIndexes(0, 0, 0);
-    return jmodel;
-  }
-};
-
-template<typename Scalar, int Options>
-struct init<pinocchio::JointModelEllipsoidTpl<Scalar, Options>>
-{
-  typedef pinocchio::JointModelEllipsoidTpl<Scalar, Options> JointModel;
-
-  static JointModel run()
-  {
-    JointModel jmodel(Scalar(0.01), Scalar(0.02), Scalar(0.03));
-
-    jmodel.setIndexes(0, 0, 0);
-    return jmodel;
-  }
-};
-
-template<typename Scalar, int Options, template<typename, int> class JointCollection>
-struct init<pinocchio::JointModelTpl<Scalar, Options, JointCollection>>
-{
-  typedef pinocchio::JointModelTpl<Scalar, Options, JointCollection> JointModel;
-
-  static JointModel run()
-  {
-    typedef pinocchio::JointModelRevoluteTpl<Scalar, Options, 0> JointModelRX;
-    JointModel jmodel((JointModelRX()));
-
-    jmodel.setIndexes(0, 0, 0);
-    return jmodel;
-  }
-};
-
-template<typename Scalar, int Options, template<typename, int> class JointCollection>
-struct init<pinocchio::JointModelCompositeTpl<Scalar, Options, JointCollection>>
-{
-  typedef pinocchio::JointModelCompositeTpl<Scalar, Options, JointCollection> JointModel;
-
-  static JointModel run()
-  {
-    typedef pinocchio::JointModelRevoluteTpl<Scalar, Options, 0> JointModelRX;
-    typedef pinocchio::JointModelRevoluteTpl<Scalar, Options, 1> JointModelRY;
-    JointModel jmodel((JointModelRX()));
-    jmodel.addJoint(JointModelRY());
-
-    jmodel.setIndexes(0, 0, 0);
-    return jmodel;
-  }
-};
-
-template<typename Scalar, int Options, template<typename, int> class JointCollection>
-struct init<pinocchio::JointModelMimicTpl<Scalar, Options, JointCollection>>
-{
-  typedef pinocchio::JointModelMimicTpl<Scalar, Options, JointCollection> JointModel;
-
-  static JointModel run()
-  {
-    typedef pinocchio::JointModelRevoluteTpl<Scalar, Options, 0> JointModelRX;
-    JointModelRX jmodel_ref = init<JointModelRX>::run();
-
-    JointModel jmodel(jmodel_ref, 1., 0.);
-    jmodel.setIndexes(0, 0, 0, 0);
-
-    return jmodel;
-  }
-};
-
 struct FiniteDiffJoint
 {
   void operator()(JointModelComposite & /*jmodel*/) const
@@ -241,14 +78,18 @@ struct FiniteDiffJoint
     typedef typename JointModel::TangentVector_t TV;
     typedef typename LieGroup<JointModel>::type LieGroupType;
 
-    JointModel jmodel = init<JointModel>::run();
+    auto JointModelWithParameters = init<JointModel>::run();
+    JointModel jmodel = JointModelWithParameters.jmodel;
+
     std::cout << "name: " << jmodel.classname() << std::endl;
 
     typename JointModel::JointDataDerived jdata_ = jmodel.createData();
     typedef JointDataBase<typename JointModel::JointDataDerived> DataBaseType;
     DataBaseType & jdata = static_cast<DataBaseType &>(jdata_);
 
-    CV q = LieGroupType().random();
+    CV q =
+      LieGroupType().randomConfiguration(JointModelWithParameters.lb, JointModelWithParameters.ub);
+
     jmodel.calc(jdata.derived(), q);
     SE3 M_ref(jdata.M());
 
@@ -274,8 +115,6 @@ struct FiniteDiffJoint
     }
 
     BOOST_CHECK(S.isApprox(S_ref, eps * 1e1));
-    std::cout << "S_ref:\n" << S_ref << std::endl;
-    std::cout << "S:\n" << S << std::endl;
   }
 };
 
diff --git a/unittest/joint-generic.cpp b/unittest/joint-generic.cpp
index 1c8a751bcc..5c4163a6d2 100644
--- a/unittest/joint-generic.cpp
+++ b/unittest/joint-generic.cpp
@@ -10,11 +10,16 @@
 #include <boost/test/unit_test.hpp>
 #include <boost/utility/binary.hpp>
 
+#include "utils/joints-init.hpp"
+
 using namespace pinocchio;
 
 template<typename JointModel>
 void test_joint_methods(
-  JointModelBase<JointModel> & jmodel, JointDataBase<typename JointModel::JointDataDerived> & jdata)
+  JointModelBase<JointModel> & jmodel,
+  JointDataBase<typename JointModel::JointDataDerived> & jdata,
+  const typename JointModel::ConfigVector_t & lb,
+  const typename JointModel::ConfigVector_t & ub)
 {
   typedef typename LieGroup<JointModel>::type LieGroupType;
   typedef typename JointModel::JointDataDerived JointData;
@@ -25,8 +30,8 @@ void test_joint_methods(
   Eigen::VectorXd armature =
     Eigen::VectorXd::Random(jdata.S().nv()) + Eigen::VectorXd::Ones(jdata.S().nv());
 
-  q1 = LieGroupType().random();
-  q2 = LieGroupType().random();
+  q1 = LieGroupType().randomConfiguration(lb, ub);
+  q2 = LieGroupType().randomConfiguration(lb, ub);
 
   Eigen::VectorXd v1(Eigen::VectorXd::Random(jdata.S().nv())),
     v2(Eigen::VectorXd::Random(jdata.S().nv()));
@@ -75,7 +80,7 @@ void test_joint_methods(
     jda.S().matrix().isApprox(jdata.S().matrix()),
     std::string(error_prefix + " - JointMotionSubspaceXd "));
   BOOST_CHECK_MESSAGE(
-    (jda.M()).isApprox((jdata.M())),
+    (jda.M()).isApprox((jdata.M()), 1e-6),
     std::string(error_prefix + " - Joint transforms ")); // ==  or isApprox ?
   BOOST_CHECK_MESSAGE(
     (jda.v()).isApprox((pinocchio::Motion(jdata.v()))),
@@ -128,169 +133,120 @@ void test_joint_methods(
   }
 }
 
-template<typename JointModel_>
-struct init;
-
-template<typename JointModel_>
-struct init
-{
-  static JointModel_ run()
-  {
-    JointModel_ jmodel;
-    jmodel.setIndexes(0, 0, 0);
-    return jmodel;
-  }
-};
-
 template<typename Scalar, int Options>
-struct init<pinocchio::JointModelRevoluteUnalignedTpl<Scalar, Options>>
-{
-  typedef pinocchio::JointModelRevoluteUnalignedTpl<Scalar, Options> JointModel;
-
-  static JointModel run()
-  {
-    typedef typename JointModel::Vector3 Vector3;
-    JointModel jmodel(Vector3::Random().normalized());
-
-    jmodel.setIndexes(0, 0, 0);
-    return jmodel;
-  }
-};
-
-template<typename Scalar, int Options>
-struct init<pinocchio::JointModelRevoluteUnboundedUnalignedTpl<Scalar, Options>>
+void test_joint_methods(
+  JointModelBase<JointModelSplineTpl<Scalar, Options>> & jmodel,
+  JointDataBase<typename JointModelSplineTpl<Scalar, Options>::JointDataDerived> & jdata)
 {
-  typedef pinocchio::JointModelRevoluteUnboundedUnalignedTpl<Scalar, Options> JointModel;
+  typedef typename LieGroup<JointModelSplineTpl<Scalar, Options>>::type LieGroupType;
+  typedef typename JointModelSplineTpl<Scalar, Options>::JointDataDerived JointData;
 
-  static JointModel run()
-  {
-    typedef typename JointModel::Vector3 Vector3;
-    JointModel jmodel(Vector3::Random().normalized());
+  std::cout << "Testing Joint over " << jmodel.shortname() << std::endl;
 
-    jmodel.setIndexes(0, 0, 0);
-    return jmodel;
-  }
-};
+  Eigen::VectorXd q1, q2;
+  Eigen::VectorXd armature =
+    Eigen::VectorXd::Random(jdata.S().nv()) + Eigen::VectorXd::Ones(jdata.S().nv());
 
-template<typename Scalar, int Options>
-struct init<pinocchio::JointModelPrismaticUnalignedTpl<Scalar, Options>>
-{
-  typedef pinocchio::JointModelPrismaticUnalignedTpl<Scalar, Options> JointModel;
+  q1 = LieGroupType().randomConfiguration(
+    Eigen::VectorXd::Zero(jmodel.nq()), Eigen::VectorXd::Ones(jmodel.nq()));
+  q2 = LieGroupType().randomConfiguration(
+    Eigen::VectorXd::Zero(jmodel.nq()), Eigen::VectorXd::Ones(jmodel.nq()));
 
-  static JointModel run()
-  {
-    typedef typename JointModel::Vector3 Vector3;
-    JointModel jmodel(Vector3::Random().normalized());
+  Eigen::VectorXd v1(Eigen::VectorXd::Random(jdata.S().nv())),
+    v2(Eigen::VectorXd::Random(jdata.S().nv()));
 
-    jmodel.setIndexes(0, 0, 0);
-    return jmodel;
-  }
-};
+  Inertia::Matrix6 Ia(pinocchio::Inertia::Random().matrix()),
+    Ia2(pinocchio::Inertia::Random().matrix());
+  bool update_I = false;
 
-template<typename Scalar, int Options, template<typename, int> class JointCollection>
-struct init<pinocchio::JointModelTpl<Scalar, Options, JointCollection>>
-{
-  typedef pinocchio::JointModelTpl<Scalar, Options, JointCollection> JointModel;
+  jmodel.calc(jdata.derived(), q1, v1);
+  jmodel.calc_aba(jdata.derived(), armature, Ia, update_I);
 
-  static JointModel run()
-  {
-    typedef pinocchio::JointModelRevoluteTpl<Scalar, Options, 0> JointModelRX;
-    JointModel jmodel((JointModelRX()));
+  pinocchio::JointModel jma(jmodel);
+  BOOST_CHECK(jmodel == jma);
+  BOOST_CHECK(jma == jmodel);
+  BOOST_CHECK(jma.hasSameIndexes(jmodel));
 
-    jmodel.setIndexes(0, 0, 0);
-    return jmodel;
-  }
-};
+  pinocchio::JointData jda(jdata.derived());
+  BOOST_CHECK(jda == jdata);
+  BOOST_CHECK(jdata == jda);
 
-template<typename Scalar, int Options, template<typename, int> class JointCollection>
-struct init<pinocchio::JointModelCompositeTpl<Scalar, Options, JointCollection>>
-{
-  typedef pinocchio::JointModelCompositeTpl<Scalar, Options, JointCollection> JointModel;
+  jma.calc(jda, q1, v1);
+  jma.calc_aba(jda, armature, Ia, update_I);
+  pinocchio::JointData jda_other(jdata);
 
-  static JointModel run()
-  {
-    typedef pinocchio::JointModelRevoluteTpl<Scalar, Options, 0> JointModelRX;
-    typedef pinocchio::JointModelRevoluteTpl<Scalar, Options, 1> JointModelRY;
-    JointModel jmodel((JointModelRX()));
-    jmodel.addJoint(JointModelRY());
+  jma.calc(jda_other, q2, v2);
+  jma.calc_aba(jda_other, armature, Ia2, update_I);
 
-    jmodel.setIndexes(0, 0, 0);
-    return jmodel;
-  }
-};
-template<typename Scalar, int Options, template<typename, int> class JointCollection>
-struct init<pinocchio::JointModelMimicTpl<Scalar, Options, JointCollection>>
-{
-  typedef pinocchio::JointModelMimicTpl<Scalar, Options, JointCollection> JointModel;
+  BOOST_CHECK(jda_other != jda);
+  BOOST_CHECK(jda != jda_other);
+  BOOST_CHECK(jda_other != jdata);
+  BOOST_CHECK(jdata != jda_other);
 
-  static JointModel run()
-  {
-    typedef pinocchio::JointModelRevoluteTpl<Scalar, Options, 0> JointModelRX;
-    JointModelRX jmodel_ref = init<JointModelRX>::run();
+  const std::string error_prefix("JointModel on " + jma.shortname());
+  BOOST_CHECK_MESSAGE(jmodel.nq() == jma.nq(), std::string(error_prefix + " - nq "));
+  BOOST_CHECK_MESSAGE(jmodel.nv() == jma.nv(), std::string(error_prefix + " - nv "));
 
-    JointModel jmodel(jmodel_ref, 1., 0.);
-    jmodel.setIndexes(0, 0, 0, 0);
+  BOOST_CHECK_MESSAGE(jmodel.idx_q() == jma.idx_q(), std::string(error_prefix + " - Idx_q "));
+  BOOST_CHECK_MESSAGE(jmodel.idx_v() == jma.idx_v(), std::string(error_prefix + " - Idx_v "));
+  BOOST_CHECK_MESSAGE(jmodel.id() == jma.id(), std::string(error_prefix + " - JointId "));
 
-    return jmodel;
-  }
-};
+  BOOST_CHECK_MESSAGE(
+    jda.S().matrix().isApprox(jdata.S().matrix()),
+    std::string(error_prefix + " - JointMotionSubspaceXd "));
+  BOOST_CHECK_MESSAGE(
+    (jda.M()).isApprox((jdata.M()), 1e-6),
+    std::string(error_prefix + " - Joint transforms ")); // ==  or isApprox ?
+  BOOST_CHECK_MESSAGE(
+    (jda.v()).isApprox((pinocchio::Motion(jdata.v()))),
+    std::string(error_prefix + " - Joint motions "));
+  BOOST_CHECK_MESSAGE((jda.c()) == (jdata.c()), std::string(error_prefix + " - Joint bias "));
 
-template<typename Scalar, int Options>
-struct init<pinocchio::JointModelUniversalTpl<Scalar, Options>>
-{
-  typedef pinocchio::JointModelUniversalTpl<Scalar, Options> JointModel;
+  BOOST_CHECK_MESSAGE(
+    (jda.U()).isApprox(jdata.U()),
+    std::string(error_prefix + " - Joint U inertia matrix decomposition "));
+  BOOST_CHECK_MESSAGE(
+    (jda.Dinv()).isApprox(jdata.Dinv()),
+    std::string(error_prefix + " - Joint DInv inertia matrix decomposition "));
+  BOOST_CHECK_MESSAGE(
+    (jda.UDinv()).isApprox(jdata.UDinv()),
+    std::string(error_prefix + " - Joint UDInv inertia matrix decomposition "));
 
-  static JointModel run()
-  {
-    JointModel jmodel(XAxis::vector(), YAxis::vector());
+  // Test vxS
+  typedef typename JointModel::Constraint_t Constraint_t;
+  typedef typename Constraint_t::DenseBase ConstraintDense;
 
-    jmodel.setIndexes(0, 0, 0);
-    return jmodel;
-  }
-};
+  Motion v(Motion::Random());
+  ConstraintDense vxS(v.cross(jdata.S()));
+  ConstraintDense vxS_ref = v.toActionMatrix() * jdata.S().matrix();
 
-template<typename Scalar, int Options, int axis>
-struct init<pinocchio::JointModelHelicalTpl<Scalar, Options, axis>>
-{
-  typedef pinocchio::JointModelHelicalTpl<Scalar, Options, axis> JointModel;
+  BOOST_CHECK_MESSAGE(vxS.isApprox(vxS_ref), std::string(error_prefix + "- Joint vxS operation "));
 
-  static JointModel run()
-  {
-    JointModel jmodel(static_cast<Scalar>(0.5));
+  // Test Y*S
+  const Inertia Isparse(Inertia::Random());
+  const Inertia::Matrix6 Idense(Isparse.matrix());
 
-    jmodel.setIndexes(0, 0, 0);
-    return jmodel;
-  }
-};
+  const ConstraintDense IsparseS = Isparse * jdata.S();
+  const ConstraintDense IdenseS = Idense * jdata.S();
 
-template<typename Scalar, int Options>
-struct init<pinocchio::JointModelEllipsoidTpl<Scalar, Options>>
-{
-  typedef pinocchio::JointModelEllipsoidTpl<Scalar, Options> JointModel;
+  BOOST_CHECK_MESSAGE(
+    IdenseS.isApprox(IsparseS), std::string(error_prefix + "- Joint YS operation "));
 
-  static JointModel run()
+  // Test calc
   {
-    JointModel jmodel(Scalar(0.01), Scalar(0.02), Scalar(0.03));
-
-    jmodel.setIndexes(0, 0, 0);
-    return jmodel;
-  }
-};
+    JointData jdata1(jdata.derived());
 
-template<typename Scalar, int Options>
-struct init<pinocchio::JointModelHelicalUnalignedTpl<Scalar, Options>>
-{
-  typedef pinocchio::JointModelHelicalUnalignedTpl<Scalar, Options> JointModel;
+    jmodel.calc(jdata1.derived(), q1, v1);
+    jmodel.calc(jdata1.derived(), Blank(), v2);
 
-  static JointModel run()
-  {
-    typedef typename JointModel::Vector3 Vector3;
-    JointModel jmodel(Vector3::Random().normalized());
+    JointData jdata_ref(jdata.derived());
+    jmodel.calc(jdata_ref.derived(), q1, v2);
 
-    jmodel.setIndexes(0, 0, 0);
-    return jmodel;
+    BOOST_CHECK_MESSAGE(
+      pinocchio::JointData(jdata1).v() == pinocchio::JointData(jdata_ref).v(),
+      std::string(error_prefix + "- joint.calc(jdata,*,v) "));
   }
-};
+}
 
 struct TestJoint
 {
@@ -298,11 +254,11 @@ struct TestJoint
   template<typename JointModel>
   void operator()(const JointModelBase<JointModel> &) const
   {
-    JointModel jmodel = init<JointModel>::run();
-    jmodel.setIndexes(0, 0, 0);
-    typename JointModel::JointDataDerived jdata = jmodel.createData();
+    auto jmodelParams = init<JointModel>::run();
+    JointModel jmodel = jmodelParams.jmodel;
 
-    test_joint_methods(jmodel, jdata);
+    typename JointModel::JointDataDerived jdata = jmodel.createData();
+    test_joint_methods(jmodel, jdata, jmodelParams.lb, jmodelParams.ub);
   }
 
   void operator()(const pinocchio::JointModelComposite &) const
@@ -465,13 +421,15 @@ struct TestJointOperatorEqual
   template<typename JointModel>
   void operator()(const JointModelBase<JointModel> &) const
   {
-    JointModel jmodel_init = init<JointModel>::run();
+    auto jmodelParams = init<JointModel>::run();
+    JointModel jmodel_init = jmodelParams.jmodel;
+
     typedef typename JointModel::JointDataDerived JointData;
 
     Model model;
     model.addJoint(0, jmodel_init, SE3::Random(), "toto");
-    model.lowerPositionLimit.fill(-1.);
-    model.upperPositionLimit.fill(1.);
+    model.lowerPositionLimit = jmodelParams.lb;
+    model.upperPositionLimit = jmodelParams.ub;
 
     const JointModel & jmodel = boost::get<JointModel>(model.joints[1]);
 
diff --git a/unittest/joint-motion-subspace.cpp b/unittest/joint-motion-subspace.cpp
index 5316418ee6..b71ccae565 100644
--- a/unittest/joint-motion-subspace.cpp
+++ b/unittest/joint-motion-subspace.cpp
@@ -2,6 +2,8 @@
 // Copyright (c) 2015-2019 CNRS INRIA
 //
 
+#include "utils/joints-init.hpp"
+
 #include "pinocchio/spatial.hpp"
 #include "pinocchio/multibody.hpp"
 
@@ -251,165 +253,16 @@ void test_constraint_operations(const JointModelBase<JointModel> & jmodel)
   }
 }
 
-template<typename Scalar, int Options, template<typename, int> class JointCollection>
-void test_constraint_operations(
-  const JointModelMimicTpl<Scalar, Options, JointCollection> & /*jmodel*/)
-{
-} // Disable test for JointMimic
-
-template<typename JointModel_>
-struct init;
-
-template<typename JointModel_>
-struct init
-{
-  static JointModel_ run()
-  {
-    JointModel_ jmodel;
-    jmodel.setIndexes(0, 0, 0);
-    return jmodel;
-  }
-};
-
-template<typename Scalar, int Options>
-struct init<pinocchio::JointModelRevoluteUnalignedTpl<Scalar, Options>>
-{
-  typedef pinocchio::JointModelRevoluteUnalignedTpl<Scalar, Options> JointModel;
-
-  static JointModel run()
-  {
-    typedef typename JointModel::Vector3 Vector3;
-    JointModel jmodel(Vector3::Random().normalized());
-
-    jmodel.setIndexes(0, 0, 0);
-    return jmodel;
-  }
-};
-
 template<typename Scalar, int Options>
-struct init<pinocchio::JointModelRevoluteUnboundedUnalignedTpl<Scalar, Options>>
+void test_constraint_operations(const JointModelSplineTpl<Scalar, Options> & /*jmodel*/)
 {
-  typedef pinocchio::JointModelRevoluteUnboundedUnalignedTpl<Scalar, Options> JointModel;
-
-  static JointModel run()
-  {
-    typedef typename JointModel::Vector3 Vector3;
-    JointModel jmodel(Vector3::Random().normalized());
-
-    jmodel.setIndexes(0, 0, 0);
-    return jmodel;
-  }
-};
-
-template<typename Scalar, int Options>
-struct init<pinocchio::JointModelPrismaticUnalignedTpl<Scalar, Options>>
-{
-  typedef pinocchio::JointModelPrismaticUnalignedTpl<Scalar, Options> JointModel;
-
-  static JointModel run()
-  {
-    typedef typename JointModel::Vector3 Vector3;
-    JointModel jmodel(Vector3::Random().normalized());
-
-    jmodel.setIndexes(0, 0, 0);
-    return jmodel;
-  }
-};
+} // Disable test for JointSpline, bc using generic subspace
 
 template<typename Scalar, int Options, template<typename, int> class JointCollection>
-struct init<pinocchio::JointModelTpl<Scalar, Options, JointCollection>>
-{
-  typedef pinocchio::JointModelTpl<Scalar, Options, JointCollection> JointModel;
-
-  static JointModel run()
-  {
-    typedef pinocchio::JointModelRevoluteTpl<Scalar, Options, 0> JointModelRX;
-    JointModel jmodel((JointModelRX()));
-
-    jmodel.setIndexes(0, 0, 0);
-    return jmodel;
-  }
-};
-
-template<typename Scalar, int Options>
-struct init<pinocchio::JointModelUniversalTpl<Scalar, Options>>
-{
-  typedef pinocchio::JointModelUniversalTpl<Scalar, Options> JointModel;
-
-  static JointModel run()
-  {
-    JointModel jmodel(XAxis::vector(), YAxis::vector());
-
-    jmodel.setIndexes(0, 0, 0);
-    return jmodel;
-  }
-};
-
-template<typename Scalar, int Options, template<typename, int> class JointCollection>
-struct init<pinocchio::JointModelCompositeTpl<Scalar, Options, JointCollection>>
-{
-  typedef pinocchio::JointModelCompositeTpl<Scalar, Options, JointCollection> JointModel;
-
-  static JointModel run()
-  {
-    typedef pinocchio::JointModelRevoluteTpl<Scalar, Options, 0> JointModelRX;
-    typedef pinocchio::JointModelRevoluteTpl<Scalar, Options, 1> JointModelRY;
-    typedef pinocchio::JointModelRevoluteTpl<Scalar, Options, 2> JointModelRZ;
-
-    JointModel jmodel(JointModelRX(), SE3::Random());
-    jmodel.addJoint(JointModelRY(), SE3::Random());
-    jmodel.addJoint(JointModelRZ(), SE3::Random());
-
-    jmodel.setIndexes(0, 0, 0);
-
-    return jmodel;
-  }
-};
-
-template<typename Scalar, int Options, template<typename, int> class JointCollection>
-struct init<pinocchio::JointModelMimicTpl<Scalar, Options, JointCollection>>
-{
-  typedef pinocchio::JointModelMimicTpl<Scalar, Options, JointCollection> JointModel;
-
-  static JointModel run()
-  {
-    typedef pinocchio::JointModelRevoluteTpl<Scalar, Options, 0> JointModelRX;
-    JointModelRX jmodel_ref = init<JointModelRX>::run();
-    JointModel jmodel(jmodel_ref, 1., 0.);
-    jmodel.setIndexes(1, 0, 0, 0);
-
-    return jmodel;
-  }
-};
-
-template<typename Scalar, int Options, int axis>
-struct init<pinocchio::JointModelHelicalTpl<Scalar, Options, axis>>
-{
-  typedef pinocchio::JointModelHelicalTpl<Scalar, Options, axis> JointModel;
-
-  static JointModel run()
-  {
-    JointModel jmodel(static_cast<Scalar>(0.5));
-
-    jmodel.setIndexes(0, 0, 0);
-    return jmodel;
-  }
-};
-
-template<typename Scalar, int Options>
-struct init<pinocchio::JointModelHelicalUnalignedTpl<Scalar, Options>>
+void test_constraint_operations(
+  const JointModelMimicTpl<Scalar, Options, JointCollection> & /*jmodel*/)
 {
-  typedef pinocchio::JointModelHelicalUnalignedTpl<Scalar, Options> JointModel;
-
-  static JointModel run()
-  {
-    typedef typename JointModel::Vector3 Vector3;
-    JointModel jmodel(Vector3::Random().normalized());
-
-    jmodel.setIndexes(0, 0, 0);
-    return jmodel;
-  }
-};
+} // Disable test for JointMimic
 
 struct TestJointConstraint
 {
@@ -417,7 +270,7 @@ struct TestJointConstraint
   template<typename JointModel>
   void operator()(const JointModelBase<JointModel> &) const
   {
-    JointModel jmodel = init<JointModel>::run();
+    JointModel jmodel = init<JointModel>::run().jmodel;
     test_constraint_operations(jmodel);
   }
 };
diff --git a/unittest/joint-spline.cpp b/unittest/joint-spline.cpp
new file mode 100644
index 0000000000..f82ef583d8
--- /dev/null
+++ b/unittest/joint-spline.cpp
@@ -0,0 +1,1460 @@
+//
+// Copyright (c) 2025 INRIA
+//
+
+#include "pinocchio/spatial.hpp"
+#include "pinocchio/multibody/joint/joints.hpp"
+#include "pinocchio/algorithm/rnea.hpp"
+#include "pinocchio/algorithm/aba.hpp"
+
+#include <boost/test/unit_test.hpp>
+
+using namespace pinocchio;
+
+namespace
+{
+
+  template<typename D>
+  void addJointAndBody(
+    Model & model,
+    const JointModelBase<D> & jmodel,
+    const Model::JointIndex parent_id,
+    const SE3 & joint_placement,
+    const std::string & joint_name,
+    const Inertia & Y)
+  {
+    Model::JointIndex idx;
+
+    idx = model.addJoint(parent_id, jmodel, joint_placement, joint_name);
+    model.appendBodyToJoint(idx, Y);
+  }
+
+  /// @brief Get a predefined  spline joint trajectory.
+  /// This was taken from the model opensim gait10dof18mus.osim, and correspond the knee joint.
+  void getTrajectory(std::vector<SE3> & ctrlFrames)
+  {
+    Eigen::Matrix3d rotation;
+    Eigen::Vector3d translation;
+    rotation << -0.500004, 0, -0.866023, 0, 1, 0, 0.866023, 0, -0.500004;
+    translation << -0.0032, -0.4226, 0;
+    ctrlFrames.push_back(pinocchio::SE3(rotation, translation));
+
+    rotation << -0.173649, 0, -0.984808, 0, 1, 0, 0.984808, 0, -0.173649;
+    translation << 0.00179, -0.416947, 0;
+    ctrlFrames.push_back(pinocchio::SE3(rotation, translation));
+
+    rotation << 0.173652, 0, -0.984807, 0, 1, 0, 0.984807, 0, 0.173652;
+    translation << 0.00411, -0.411057, 0;
+    ctrlFrames.push_back(pinocchio::SE3(rotation, translation));
+
+    rotation << 0.342021, 0, -0.939692, 0, 1, 0, 0.939692, 0, 0.342021;
+    translation << 0.00438827, -0.4082, 0;
+    ctrlFrames.push_back(pinocchio::SE3(rotation, translation));
+
+    rotation << 0.499998, 0, -0.866027, 0, 1, 0, 0.866027, 0, 0.499998;
+    translation << 0.0041, -0.405495, 0;
+    ctrlFrames.push_back(pinocchio::SE3(rotation, translation));
+
+    rotation << 0.766044, 0, -0.642788, 0, 1, 0, 0.642788, 0, 0.766044;
+    translation << 0.00212, -0.400825, 0;
+    ctrlFrames.push_back(pinocchio::SE3(rotation, translation));
+
+    rotation << 0.866025, 0, -0.5, 0, 1, 0, 0.5, 0, 0.866025;
+    translation << 0.000757726, -0.399, 0;
+    ctrlFrames.push_back(pinocchio::SE3(rotation, translation));
+
+    rotation << 0.939693, 0, -0.34202, 0, 1, 0, 0.34202, 0, 0.939693;
+    translation << -0.001, -0.3976, 0;
+    ctrlFrames.push_back(pinocchio::SE3(rotation, translation));
+
+    rotation << 0.984808, 0, -0.173648, 0, 1, 0, 0.173648, 0, 0.984808;
+    translation << -0.0031, -0.3966, 0;
+    ctrlFrames.push_back(pinocchio::SE3(rotation, translation));
+
+    rotation << 0.987363, 0, 0.158478, 0, 1, 0, -0.158478, 0, 0.987363;
+    translation << -0.00513719, -0.395264, 0;
+    ctrlFrames.push_back(pinocchio::SE3(rotation, translation));
+
+    rotation << 0.980591, 0, 0.196066, 0, 1, 0, -0.196066, 0, 0.980591;
+    translation << -0.005227, -0.395149, 0;
+    ctrlFrames.push_back(pinocchio::SE3(rotation, translation));
+
+    rotation << 0.94362, 0, 0.33103, 0, 1, 0, -0.33103, 0, 0.94362;
+    translation << -0.005435, -0.394792, 0;
+    ctrlFrames.push_back(pinocchio::SE3(rotation, translation));
+
+    rotation << 0.882249, 0, 0.470783, 0, 1, 0, -0.470783, 0, 0.882249;
+    translation << -0.005574, -0.394507, 0;
+    ctrlFrames.push_back(pinocchio::SE3(rotation, translation));
+
+    rotation << 0.0493163, 0, 0.998783, 0, 1, 0, -0.998783, 0, 0.0493163;
+    translation << -0.005435, -0.394812, 0;
+    ctrlFrames.push_back(pinocchio::SE3(rotation, translation));
+
+    rotation << -0.500004, 0, 0.866023, 0, 1, 0, -0.866023, 0, -0.500004;
+    translation << -0.00525, -0.396, 0;
+    ctrlFrames.push_back(pinocchio::SE3(rotation, translation));
+  }
+
+  double bsplineBasis(int index, int degree, const double x, const Eigen::VectorXd & knots)
+  {
+    assert(0 <= index);
+    assert(0 <= degree);
+    if (degree == 0)
+    {
+      if (knots[index] <= x && x < knots[index + 1])
+      {
+        return 1.;
+      }
+      else
+      {
+        if (x == knots[knots.size() - 1] && x == knots[index + 1])
+        {
+          return 1;
+        }
+        return 0.;
+      }
+    }
+
+    const double den1 = knots[index + degree] - knots[index];
+    double left = 0.;
+    if (den1 != 0.)
+    {
+      left = ((x - knots[index]) / den1) * bsplineBasis(index, degree - 1, x, knots);
+    }
+
+    const double den2 = knots[index + degree + 1] - knots[index + 1];
+    double right = 0.;
+    if (den2 != 0.)
+    {
+      right =
+        ((knots[index + degree + 1] - x) / den2) * bsplineBasis(index + 1, degree - 1, x, knots);
+    }
+
+    return left + right;
+  }
+
+  double
+  bsplineBasisDerivative(int index, int degree, const double x, const Eigen::VectorXd & knots)
+  {
+    if (degree == 0)
+    {
+      return 0.;
+    }
+
+    const double den1 = knots[index + degree] - knots[index];
+    double term1 = 0.;
+    if (den1 != 0.)
+    {
+      term1 = bsplineBasis(index, degree - 1, x, knots) / den1;
+    }
+
+    const double den2 = knots[index + degree + 1] - knots[index + 1];
+    double term2 = 0.;
+    if (den2 != 0.)
+    {
+      term2 = bsplineBasis(index + 1, degree - 1, x, knots) / den2;
+    }
+
+    return degree * (term1 - term2);
+  }
+
+  double
+  bsplineBasisDerivative2(int index, int degree, const double x, const Eigen::VectorXd & knots)
+  {
+    if (degree < 2)
+    {
+      return 0.;
+    }
+
+    const double den1 = knots[index + degree] - knots[index];
+    double term1 = 0.;
+    if (den1 != 0.)
+    {
+      term1 = bsplineBasisDerivative(index, degree - 1, x, knots) / den1;
+    }
+
+    const double den2 = knots[index + degree + 1] - knots[index + 1];
+    double term2 = 0.;
+    if (den2 != 0.)
+    {
+      term2 = bsplineBasisDerivative(index + 1, degree - 1, x, knots) / den2;
+    }
+
+    return degree * (term1 - term2);
+  }
+
+} // namespace
+
+BOOST_AUTO_TEST_SUITE(JointSpline)
+
+/// @brief Check the joint builder and the guards.
+BOOST_AUTO_TEST_CASE(jointBuilder)
+{
+  size_t degree = 3;
+
+  std::vector<SE3> ctrlFrames;
+  for (int k = 0; k < 3; k++)
+    ctrlFrames.push_back(SE3::Random());
+
+  BOOST_CHECK_THROW(
+    JointModelSplineBuilder()
+      .withControlFrameVector(ctrlFrames)
+      .withDegree(degree)
+      .withOpenUniformKnots(0., 1.)
+      .build(),
+    std::invalid_argument);
+
+  for (int k = 0; k < 3; k++)
+    ctrlFrames.push_back(SE3::Random());
+
+  // Knots vector value should not decrease (t_i <= t_{i+1})
+  Eigen::VectorXd knots_non_uniform(degree + ctrlFrames.size() + 1);
+  knots_non_uniform << 0., 0.1, 0.08, 0.15, 0.15, 0.3, 0.6, 0.6, 0.7, 1.;
+
+  BOOST_CHECK_THROW(
+    JointModelSplineBuilder()
+      .withControlFrameVector(ctrlFrames)
+      .withDegree(degree)
+      .withKnotVector(knots_non_uniform)
+      .build(),
+    std::invalid_argument);
+
+  // Knot vector value should not be repeated more than degree + 1 times
+  Eigen::VectorXd knots_too_many_repeats(degree + ctrlFrames.size() + 1);
+  knots_too_many_repeats << 0., 0., 0., 0., 0., 0.3, 0.6, 0.6, 0.7, 1.;
+
+  BOOST_CHECK_THROW(
+    JointModelSplineBuilder()
+      .withControlFrameVector(ctrlFrames)
+      .withDegree(degree)
+      .withKnotVector(knots_too_many_repeats)
+      .build(),
+    std::invalid_argument);
+}
+
+/// @brief Test on the knot vector generation
+BOOST_AUTO_TEST_CASE(makeKnots)
+{
+  size_t degree = 3;
+  size_t nbFrames = 6;
+  double min_q = 10;
+  double max_q = 40;
+
+  // Open Uniform
+  Eigen::VectorXd generated_knots =
+    internal::generateOpenUniformKnots(min_q, max_q, nbFrames, degree);
+
+  // Check size
+  BOOST_CHECK(generated_knots.size() == static_cast<Eigen::Index>(degree + nbFrames + 1));
+
+  // Check Values
+  Eigen::VectorXd open_uniform(degree + nbFrames + 1);
+  open_uniform << 10., 10., 10., 10., 20., 30., 40., 40., 40., 40.;
+  BOOST_CHECK(generated_knots.isApprox(open_uniform, 1e-5));
+
+  // Open Uniform
+  generated_knots = internal::generateUniformKnots(min_q, max_q, nbFrames, degree);
+
+  // Check size
+  BOOST_CHECK(generated_knots.size() == static_cast<Eigen::Index>(degree + nbFrames + 1));
+
+  // Check Values
+  Eigen::VectorXd uniform(degree + nbFrames + 1);
+  uniform << 10., 13.3333, 16.6667, 20.0, 23.3333, 26.6667, 30., 33.3333, 36.6667, 40.;
+  BOOST_CHECK(generated_knots.isApprox(uniform, 1e-5));
+}
+
+// Test recursive bsplineBasis function.
+// This function is only implemented in unit test, but since
+// it's used to test deBoor algorithm it worst testing it.
+// Test bsplineBasis node limit.
+// We test the degree 0 case, because other degree bound depend of it.
+BOOST_AUTO_TEST_CASE(basisFunctionsEdge)
+{
+  int degree = 0;
+
+  Eigen::VectorXd knotVector(6);
+  knotVector << 0., 0.2, 0.4, 0.6, 0.8, 1.;
+  BOOST_CHECK_EQUAL(bsplineBasis(0, degree, 0., knotVector), 1.);
+  BOOST_CHECK_EQUAL(bsplineBasis(0, degree, 0.2, knotVector), 0.);
+  BOOST_CHECK_EQUAL(bsplineBasis(1, degree, 0.2, knotVector), 1.);
+  BOOST_CHECK_EQUAL(bsplineBasis(1, degree, 0.4, knotVector), 0.);
+  BOOST_CHECK_EQUAL(bsplineBasis(2, degree, 0.4, knotVector), 1.);
+  BOOST_CHECK_EQUAL(bsplineBasis(2, degree, 0.6, knotVector), 0.);
+  BOOST_CHECK_EQUAL(bsplineBasis(3, degree, 0.6, knotVector), 1.);
+  BOOST_CHECK_EQUAL(bsplineBasis(3, degree, 0.8, knotVector), 0.);
+  BOOST_CHECK_EQUAL(bsplineBasis(4, degree, 0.8, knotVector), 1.);
+  BOOST_CHECK_EQUAL(bsplineBasis(4, degree, 1., knotVector), 1.);
+}
+
+// Test recursive bsplineBasis, bsplineDerivative and bsplineDerivative2 functions.
+// This function is only implemented in unit test, but since
+// it's used to test deBoor algorithm it worst testing it.
+// Test on an open uniform knot vector.
+BOOST_AUTO_TEST_CASE(basisFunctionsOpenUniform)
+{
+  int degree = 3;
+  int nbCtrlFrames = 6;
+
+  double min_q = 0.0;
+  double max_q = 1.;
+
+  auto knotVector = internal::generateOpenUniformKnots(
+    min_q, max_q, static_cast<size_t>(nbCtrlFrames), static_cast<size_t>(degree));
+  Eigen::Index nKnot = knotVector.size();
+  Eigen::VectorXd N = Eigen::VectorXd::Zero(nbCtrlFrames);
+  Eigen::VectorXd Nder = Eigen::VectorXd::Zero(nbCtrlFrames);
+  Eigen::VectorXd Nder2 = Eigen::VectorXd::Zero(nbCtrlFrames);
+
+  // Index Interval for unite partition [degree; nKnot - 1 - degree]
+  for (double q = knotVector[degree]; q <= knotVector[nKnot - 1 - degree]; q += 0.02)
+  {
+    for (int i = 0; i < nbCtrlFrames; i++)
+    {
+      N[i] = bsplineBasis(i, degree, q, knotVector);
+    }
+    BOOST_CHECK_CLOSE(N.sum(), 1.0, 1e-8);
+  }
+
+  // Check derivatives
+  double h = 1e-5;
+  for (double q = min_q + h; q < max_q; q += 0.02)
+  {
+    for (int i = 0; i < nbCtrlFrames; i++)
+    {
+      N[i] = bsplineBasis(i, degree, q, knotVector);
+      Nder[i] = bsplineBasisDerivative(i, degree, q, knotVector);
+      Nder2[i] = bsplineBasisDerivative2(i, degree, q, knotVector);
+
+      double n_plus = bsplineBasis(i, degree, q + h, knotVector);
+      double n_minus = bsplineBasis(i, degree, q - h, knotVector);
+
+      // First Derivative Approximation
+      double numerical_der = (n_plus - n_minus) / (2.0 * h);
+      BOOST_CHECK_SMALL(Nder[i] - numerical_der, 1e-5);
+
+      // Second Derivative Approximation
+      double n_mid = N[i];
+      double numerical_der2 = (n_plus - 2.0 * n_mid + n_minus) / (h * h);
+      BOOST_CHECK_SMALL(Nder2[i] - numerical_der2, 1e-5);
+    }
+  }
+}
+
+// Test recursive bsplineBasis, bsplineDerivative and bsplineDerivative2 functions.
+// This function is only implemented in unit test, but since
+// it's used to test deBoor algorithm it worst testing it.
+// Test on an uniform knot vector.
+BOOST_AUTO_TEST_CASE(basisFunctionsUniform)
+{
+  int degree = 3;
+  int nbCtrlFrames = 6;
+  double min_q = 0.0;
+  double max_q = 10.;
+
+  auto knotVector = internal::generateUniformKnots(
+    min_q, max_q, static_cast<size_t>(nbCtrlFrames), static_cast<size_t>(degree));
+  Eigen::Index nKnot = knotVector.size();
+
+  Eigen::VectorXd N = Eigen::VectorXd::Zero(nbCtrlFrames);
+  Eigen::VectorXd Nder = Eigen::VectorXd::Zero(nbCtrlFrames);
+  Eigen::VectorXd Nder2 = Eigen::VectorXd::Zero(nbCtrlFrames);
+  // Index Interval for unite partition [degree; len(KnotVector) - degree]
+  for (double q = knotVector[degree]; q <= knotVector[nKnot - 1 - degree]; q += 0.05)
+  {
+    for (int i = 0; i < nbCtrlFrames; i++)
+    {
+      N[i] = bsplineBasis(i, degree, q, knotVector);
+    }
+    BOOST_CHECK_CLOSE(N.sum(), 1.0, 1e-8);
+  }
+
+  // Check derivatives
+  double h = 1e-5;
+  for (double q = min_q + h; q < max_q; q += 0.02)
+  {
+    for (int i = 0; i < nbCtrlFrames; i++)
+    {
+      N[i] = bsplineBasis(i, degree, q, knotVector);
+      Nder[i] = bsplineBasisDerivative(i, degree, q, knotVector);
+      Nder2[i] = bsplineBasisDerivative2(i, degree, q, knotVector);
+
+      double n_plus = bsplineBasis(i, degree, q + h, knotVector);
+      double n_minus = bsplineBasis(i, degree, q - h, knotVector);
+
+      // First Derivative Approximation
+      double numerical_der = (n_plus - n_minus) / (2.0 * h);
+      BOOST_CHECK_SMALL(Nder[i] - numerical_der, 1e-5);
+
+      // Second Derivative Approximation
+      double n_mid = N[i];
+      double numerical_der2 = (n_plus - 2.0 * n_mid + n_minus) / (h * h);
+      BOOST_CHECK_SMALL(Nder2[i] - numerical_der2, 1e-5);
+    }
+  }
+}
+
+// Test recursive bsplineBasis, bsplineDerivative and bsplineDerivative2 functions.
+// This function is only implemented in unit test, but since
+// it's used to test deBoor algorithm it worst testing it.
+// Test on a non uniform knot vector.
+BOOST_AUTO_TEST_CASE(basisFunctionsNonUniform)
+{
+  int degree = 3;
+  int nbCtrlFrames = 5;
+  double min_q = 0.;
+  double max_q = 1.0;
+  Eigen::Index nKnot = degree + nbCtrlFrames + 1;
+  Eigen::VectorXd knotVector(nKnot);
+  knotVector << min_q, 0.1, 0.1, 0.25, 0.5, 0.5, 0.8, 0.8, max_q;
+
+  Eigen::VectorXd N = Eigen::VectorXd::Zero(nbCtrlFrames);
+  Eigen::VectorXd Nder = Eigen::VectorXd::Zero(nbCtrlFrames);
+  Eigen::VectorXd Nder2 = Eigen::VectorXd::Zero(nbCtrlFrames);
+
+  // Index Interval for unite partition [degree; len(KnotVector) -1 - degree]
+  for (double q = knotVector[degree]; q <= knotVector[nKnot - 1 - degree]; q += 0.05)
+  {
+    for (int i = 0; i < nbCtrlFrames; i++)
+    {
+      N[i] = bsplineBasis(i, degree, q, knotVector);
+    }
+    BOOST_CHECK_CLOSE(N.sum(), 1.0, 1e-8);
+  }
+
+  // Check derivatives
+  double h = 1e-5;
+  for (double q = min_q + h; q < max_q; q += 0.02)
+  {
+    for (int i = 0; i < nbCtrlFrames; i++)
+    {
+      N[i] = bsplineBasis(i, degree, q, knotVector);
+      Nder[i] = bsplineBasisDerivative(i, degree, q, knotVector);
+      Nder2[i] = bsplineBasisDerivative2(i, degree, q, knotVector);
+
+      double n_plus = bsplineBasis(i, degree, q + h, knotVector);
+      double n_minus = bsplineBasis(i, degree, q - h, knotVector);
+
+      // First Derivative Approximation
+      double numerical_der = (n_plus - n_minus) / (2.0 * h);
+      BOOST_CHECK_SMALL(Nder[i] - numerical_der, 1e-5);
+
+      // Second Derivative is not checked because the knot vector is not C2 everywhere
+    }
+  }
+}
+// Test degree 0 deBoorBasis.
+BOOST_AUTO_TEST_CASE(deBoorBasis_degree0)
+{
+  using pinocchio::internal::deBoorBasis;
+
+  const int degree = 0;
+  Eigen::VectorXd knots(3);
+  knots << 0., 2., 3.;
+  Eigen::MatrixXd basis(Eigen::MatrixXd::Zero(1, 1));
+
+  deBoorBasis(degree, knots, 0, 0., basis);
+  BOOST_CHECK_SMALL(basis(0, 0) - 1., 1e-8);
+  deBoorBasis(degree, knots, 0, 1., basis);
+  BOOST_CHECK_SMALL(basis(0, 0) - 1., 1e-8);
+  deBoorBasis(degree, knots, 0, 1.5, basis);
+  BOOST_CHECK_SMALL(basis(0, 0) - 1., 1e-8);
+  deBoorBasis(degree, knots, 1, 2., basis);
+  BOOST_CHECK_SMALL(basis(0, 0) - 1., 1e-8);
+  deBoorBasis(degree, knots, 1, 2.5, basis);
+  BOOST_CHECK_SMALL(basis(0, 0) - 1., 1e-8);
+  deBoorBasis(degree, knots, 1, 3., basis);
+  BOOST_CHECK_SMALL(basis(0, 0) - 1., 1e-8);
+}
+
+// Test degree 1 deBoorBasis.
+BOOST_AUTO_TEST_CASE(deBoorBasis_degree1)
+{
+  using pinocchio::internal::deBoorBasis;
+
+  const int degree = 1;
+  Eigen::VectorXd knots(5);
+  knots << 0., 2., 3., 5., 5.5;
+  Eigen::MatrixXd basis(Eigen::MatrixXd::Zero(2, 2));
+
+  // Evaluate between 2 and 3
+  deBoorBasis(degree, knots, 1, 2., basis);
+  BOOST_CHECK_SMALL(basis(1, 0) - 1., 1e-8);
+  BOOST_CHECK_SMALL(basis(1, 1) - 0., 1e-8);
+  deBoorBasis(degree, knots, 1, 2.5, basis);
+  BOOST_CHECK_SMALL(basis(1, 0) - 0.5, 1e-8);
+  BOOST_CHECK_SMALL(basis(1, 1) - 0.5, 1e-8);
+  deBoorBasis(degree, knots, 1, 3., basis);
+  BOOST_CHECK_SMALL(basis(1, 0) - 0., 1e-8);
+  BOOST_CHECK_SMALL(basis(1, 1) - 1., 1e-8);
+
+  // Evaluate between 3 and 5
+  deBoorBasis(degree, knots, 2, 3., basis);
+  BOOST_CHECK_SMALL(basis(1, 0) - 1., 1e-8);
+  BOOST_CHECK_SMALL(basis(1, 1) - 0., 1e-8);
+  deBoorBasis(degree, knots, 2, 4., basis);
+  BOOST_CHECK_SMALL(basis(1, 0) - 0.5, 1e-8);
+  BOOST_CHECK_SMALL(basis(1, 1) - 0.5, 1e-8);
+  deBoorBasis(degree, knots, 2, 5., basis);
+  BOOST_CHECK_SMALL(basis(1, 0) - 0., 1e-8);
+  BOOST_CHECK_SMALL(basis(1, 1) - 1., 1e-8);
+}
+
+// Test degree 3 deBoorBasis against bsplineBasis.
+BOOST_AUTO_TEST_CASE(deBoorBasis_degree3)
+{
+  using pinocchio::internal::deBoorBasis;
+  using pinocchio::internal::getAbsoluteBasis;
+
+  const int degree = 3;
+  Eigen::VectorXd knots(9);
+  knots << 0., 2., 3., 5., 5.5, 8., 8.5, 10., 11.5;
+  Eigen::MatrixXd basis(Eigen::MatrixXd::Zero(4, 4));
+
+  deBoorBasis(degree, knots, 3, 5., basis);
+  BOOST_CHECK_SMALL(
+    getAbsoluteBasis(3, basis, 0, degree) - bsplineBasis(0, degree, 5., knots), 1e-8);
+  BOOST_CHECK_SMALL(
+    getAbsoluteBasis(3, basis, 1, degree) - bsplineBasis(1, degree, 5., knots), 1e-8);
+  BOOST_CHECK_SMALL(
+    getAbsoluteBasis(3, basis, 2, degree) - bsplineBasis(2, degree, 5., knots), 1e-8);
+  BOOST_CHECK_SMALL(
+    getAbsoluteBasis(3, basis, 3, degree) - bsplineBasis(3, degree, 5., knots), 1e-8);
+
+  deBoorBasis(degree, knots, 3, 5.3, basis);
+  BOOST_CHECK_SMALL(
+    getAbsoluteBasis(3, basis, 0, degree) - bsplineBasis(0, degree, 5.3, knots), 1e-8);
+  BOOST_CHECK_SMALL(
+    getAbsoluteBasis(3, basis, 1, degree) - bsplineBasis(1, degree, 5.3, knots), 1e-8);
+  BOOST_CHECK_SMALL(
+    getAbsoluteBasis(3, basis, 2, degree) - bsplineBasis(2, degree, 5.3, knots), 1e-8);
+  BOOST_CHECK_SMALL(
+    getAbsoluteBasis(3, basis, 3, degree) - bsplineBasis(3, degree, 5.3, knots), 1e-8);
+
+  deBoorBasis(degree, knots, 3, 5.5, basis);
+  BOOST_CHECK_SMALL(
+    getAbsoluteBasis(3, basis, 0, degree) - bsplineBasis(0, degree, 5.5, knots), 1e-8);
+  BOOST_CHECK_SMALL(
+    getAbsoluteBasis(3, basis, 1, degree) - bsplineBasis(1, degree, 5.5, knots), 1e-8);
+  BOOST_CHECK_SMALL(
+    getAbsoluteBasis(3, basis, 2, degree) - bsplineBasis(2, degree, 5.5, knots), 1e-8);
+  BOOST_CHECK_SMALL(
+    getAbsoluteBasis(3, basis, 3, degree) - bsplineBasis(3, degree, 5.5, knots), 1e-8);
+
+  deBoorBasis(degree, knots, 4, 5.5, basis);
+  BOOST_CHECK_SMALL(
+    getAbsoluteBasis(4, basis, 1, degree) - bsplineBasis(1, degree, 5.5, knots), 1e-8);
+  BOOST_CHECK_SMALL(
+    getAbsoluteBasis(4, basis, 2, degree) - bsplineBasis(2, degree, 5.5, knots), 1e-8);
+  BOOST_CHECK_SMALL(
+    getAbsoluteBasis(4, basis, 3, degree) - bsplineBasis(3, degree, 5.5, knots), 1e-8);
+  BOOST_CHECK_SMALL(
+    getAbsoluteBasis(4, basis, 4, degree) - bsplineBasis(4, degree, 5.5, knots), 1e-8);
+
+  deBoorBasis(degree, knots, 4, 6., basis);
+  BOOST_CHECK_SMALL(
+    getAbsoluteBasis(4, basis, 1, degree) - bsplineBasis(1, degree, 6., knots), 1e-8);
+  BOOST_CHECK_SMALL(
+    getAbsoluteBasis(4, basis, 2, degree) - bsplineBasis(2, degree, 6., knots), 1e-8);
+  BOOST_CHECK_SMALL(
+    getAbsoluteBasis(4, basis, 3, degree) - bsplineBasis(3, degree, 6., knots), 1e-8);
+  BOOST_CHECK_SMALL(
+    getAbsoluteBasis(4, basis, 4, degree) - bsplineBasis(4, degree, 6., knots), 1e-8);
+
+  deBoorBasis(degree, knots, 4, 7.3, basis);
+  BOOST_CHECK_SMALL(
+    getAbsoluteBasis(4, basis, 1, degree) - bsplineBasis(1, degree, 7.3, knots), 1e-8);
+  BOOST_CHECK_SMALL(
+    getAbsoluteBasis(4, basis, 2, degree) - bsplineBasis(2, degree, 7.3, knots), 1e-8);
+  BOOST_CHECK_SMALL(
+    getAbsoluteBasis(4, basis, 3, degree) - bsplineBasis(3, degree, 7.3, knots), 1e-8);
+  BOOST_CHECK_SMALL(
+    getAbsoluteBasis(4, basis, 4, degree) - bsplineBasis(4, degree, 7.3, knots), 1e-8);
+
+  deBoorBasis(degree, knots, 4, 8., basis);
+  BOOST_CHECK_SMALL(
+    getAbsoluteBasis(4, basis, 1, degree) - bsplineBasis(1, degree, 8., knots), 1e-8);
+  BOOST_CHECK_SMALL(
+    getAbsoluteBasis(4, basis, 2, degree) - bsplineBasis(2, degree, 8., knots), 1e-8);
+  BOOST_CHECK_SMALL(
+    getAbsoluteBasis(4, basis, 3, degree) - bsplineBasis(3, degree, 8., knots), 1e-8);
+  BOOST_CHECK_SMALL(
+    getAbsoluteBasis(4, basis, 4, degree) - bsplineBasis(4, degree, 8., knots), 1e-8);
+}
+
+// Test deBoorBasisFull with a knot vector only defined
+// on one span.
+// This should give the same result than deBoorBasis function.
+BOOST_AUTO_TEST_CASE(deBoorBasisFull_degree2_minimal_knot)
+{
+  using pinocchio::internal::deBoorBasis;
+  using pinocchio::internal::deBoorFullBasis;
+
+  const int degree = 2;
+  Eigen::VectorXd knots(6);
+  // Defined between 3 and 5
+  knots << 0., 2., 3., 5., 5.5, 8.3;
+  Eigen::MatrixXd basis_ref(Eigen::MatrixXd::Zero(3, 3));
+  Eigen::MatrixXd basis(3, 3);
+
+  deBoorBasis(degree, knots, 2, 3., basis_ref);
+  deBoorFullBasis(degree, knots, 3., basis);
+  BOOST_CHECK(basis_ref.isApprox(basis));
+
+  deBoorBasis(degree, knots, 2, 4., basis_ref);
+  deBoorFullBasis(degree, knots, 4., basis);
+  BOOST_CHECK(basis_ref.isApprox(basis));
+
+  deBoorBasis(degree, knots, 2, 5., basis_ref);
+  deBoorFullBasis(degree, knots, 5., basis);
+  BOOST_CHECK(basis_ref.isApprox(basis));
+}
+
+// Test deBoorBasisFull with a knot vector allowing to compute
+// more basis function than deBoorBasis.
+// We compare it to deBoorBasis called on the right span.
+BOOST_AUTO_TEST_CASE(deBoorBasisFull_degree3_nominal_knot)
+{
+  using pinocchio::internal::deBoorBasis;
+  using pinocchio::internal::deBoorFullBasis;
+
+  const int degree = 2;
+  Eigen::VectorXd knots(8);
+  // Defined between 3 and 8.3
+  knots << 0., 2., 3., 5., 5.5, 8.3, 9., 11.3;
+  Eigen::MatrixXd basis_ref(Eigen::MatrixXd::Zero(3, 3));
+  Eigen::MatrixXd basis(3, 5);
+
+  deBoorFullBasis(degree, knots, 3., basis);
+  deBoorBasis(degree, knots, 2, 3., basis_ref);
+  BOOST_CHECK(basis_ref.isApprox(basis.block<3, 3>(0, 0)));
+  BOOST_CHECK((basis.block<3, 2>(0, 3).isZero()));
+
+  deBoorFullBasis(degree, knots, 4., basis);
+  deBoorBasis(degree, knots, 2, 4., basis_ref);
+  BOOST_CHECK(basis_ref.isApprox(basis.block<3, 3>(0, 0)));
+  BOOST_CHECK((basis.block<3, 2>(0, 3).isZero()));
+
+  deBoorFullBasis(degree, knots, 5., basis);
+  deBoorBasis(degree, knots, 3, 5., basis_ref);
+  BOOST_CHECK(basis_ref.isApprox(basis.block<3, 3>(0, 1)));
+  BOOST_CHECK((basis.col(0).isZero()));
+  BOOST_CHECK((basis.col(4).isZero()));
+
+  deBoorFullBasis(degree, knots, 5.3, basis);
+  deBoorBasis(degree, knots, 3, 5.3, basis_ref);
+  BOOST_CHECK(basis_ref.isApprox(basis.block<3, 3>(0, 1)));
+  BOOST_CHECK((basis.col(0).isZero()));
+  BOOST_CHECK((basis.col(4).isZero()));
+
+  deBoorFullBasis(degree, knots, 5.5, basis);
+  deBoorBasis(degree, knots, 4, 5.5, basis_ref);
+  BOOST_CHECK(basis_ref.isApprox(basis.block<3, 3>(0, 2)));
+  BOOST_CHECK((basis.block<3, 2>(0, 0).isZero()));
+
+  deBoorFullBasis(degree, knots, 8., basis);
+  deBoorBasis(degree, knots, 4, 8., basis_ref);
+  BOOST_CHECK(basis_ref.isApprox(basis.block<3, 3>(0, 2)));
+  BOOST_CHECK((basis.block<3, 2>(0, 0).isZero()));
+
+  deBoorFullBasis(degree, knots, 8.3, basis);
+  deBoorBasis(degree, knots, 4, 8.3, basis_ref);
+  BOOST_CHECK(basis_ref.isApprox(basis.block<3, 3>(0, 2)));
+  BOOST_CHECK((basis.block<3, 2>(0, 0).isZero()));
+}
+
+// Test deBoorBasisFull with a knot vector with element of multiplicity 3.
+// This should create division by 0 issue managed by the algorithm.
+BOOST_AUTO_TEST_CASE(deBoorBasisFull_degree3_multiplicity_knot)
+{
+  using pinocchio::internal::deBoorBasis;
+  using pinocchio::internal::deBoorFullBasis;
+
+  const int degree = 2;
+  Eigen::VectorXd knots(9);
+  // Defined between 0 and 2
+  knots << 0., 0., 0., 1., 1., 1., 2., 2., 2.;
+  Eigen::MatrixXd basis_ref(Eigen::MatrixXd::Zero(3, 3));
+  Eigen::MatrixXd basis(3, 6);
+
+  deBoorFullBasis(degree, knots, 0., basis);
+  deBoorBasis(degree, knots, 2, 0., basis_ref);
+  BOOST_CHECK(basis_ref.isApprox(basis.block<3, 3>(0, 0)));
+  BOOST_CHECK((basis.block<3, 3>(0, 3).isZero()));
+
+  deBoorFullBasis(degree, knots, 0.5, basis);
+  deBoorBasis(degree, knots, 2, 0.5, basis_ref);
+  BOOST_CHECK(basis_ref.isApprox(basis.block<3, 3>(0, 0)));
+  BOOST_CHECK((basis.block<3, 3>(0, 3).isZero()));
+
+  deBoorFullBasis(degree, knots, 1., basis);
+  deBoorBasis(degree, knots, 5, 1., basis_ref);
+  BOOST_CHECK(basis_ref.isApprox(basis.block<3, 3>(0, 3)));
+  BOOST_CHECK((basis.block<3, 3>(0, 0).isZero()));
+
+  deBoorFullBasis(degree, knots, 1.5, basis);
+  deBoorBasis(degree, knots, 5, 1.5, basis_ref);
+  BOOST_CHECK(basis_ref.isApprox(basis.block<3, 3>(0, 3)));
+  BOOST_CHECK((basis.block<3, 3>(0, 0).isZero()));
+
+  deBoorFullBasis(degree, knots, 2., basis);
+  deBoorBasis(degree, knots, 5, 2., basis_ref);
+  BOOST_CHECK(basis_ref.isApprox(basis.block<3, 3>(0, 3)));
+  BOOST_CHECK((basis.block<3, 3>(0, 0).isZero()));
+}
+
+// Test cumulativeBasisDerivative against bsplineBasisDerivative
+BOOST_AUTO_TEST_CASE(cumulativeBasisDerivative)
+{
+  using pinocchio::internal::cumulativeBasisDerivative;
+  using pinocchio::internal::deBoorBasis;
+
+  const int degree = 3;
+  Eigen::VectorXd knots(9);
+  knots << 0., 2., 3., 5., 5.5, 8., 8.5, 10., 11.5;
+  Eigen::MatrixXd basis(Eigen::MatrixXd::Zero(4, 4));
+
+  auto computeDerivative = [knots](int start, double q) {
+    Eigen::VectorXd res(4);
+    for (int i = 0; i < 4; ++i)
+    {
+      res(i) = bsplineBasisDerivative(i + start, degree, q, knots);
+    }
+    return res;
+  };
+
+  auto derivative_ref = computeDerivative(0, 5.);
+  deBoorBasis(degree, knots, 3, 5., basis);
+  BOOST_CHECK_SMALL(derivative_ref.sum(), 1e-8);
+  BOOST_CHECK_SMALL(
+    derivative_ref.tail(3).sum() - cumulativeBasisDerivative(3, knots, basis, 1, degree), 1e-8);
+  BOOST_CHECK_SMALL(
+    derivative_ref.tail(2).sum() - cumulativeBasisDerivative(3, knots, basis, 2, degree), 1e-8);
+  BOOST_CHECK_SMALL(
+    derivative_ref.tail(1).sum() - cumulativeBasisDerivative(3, knots, basis, 3, degree), 1e-8);
+
+  derivative_ref = computeDerivative(0, 5.3);
+  deBoorBasis(degree, knots, 3, 5.3, basis);
+  BOOST_CHECK_SMALL(derivative_ref.sum(), 1e-8);
+  BOOST_CHECK_SMALL(
+    derivative_ref.tail(3).sum() - cumulativeBasisDerivative(3, knots, basis, 1, degree), 1e-8);
+  BOOST_CHECK_SMALL(
+    derivative_ref.tail(2).sum() - cumulativeBasisDerivative(3, knots, basis, 2, degree), 1e-8);
+  BOOST_CHECK_SMALL(
+    derivative_ref.tail(1).sum() - cumulativeBasisDerivative(3, knots, basis, 3, degree), 1e-8);
+
+  derivative_ref = computeDerivative(1, 5.5);
+  deBoorBasis(degree, knots, 4, 5.5, basis);
+  BOOST_CHECK_SMALL(derivative_ref.sum(), 1e-8);
+  BOOST_CHECK_SMALL(
+    derivative_ref.tail(3).sum() - cumulativeBasisDerivative(4, knots, basis, 2, degree), 1e-8);
+  BOOST_CHECK_SMALL(
+    derivative_ref.tail(2).sum() - cumulativeBasisDerivative(4, knots, basis, 3, degree), 1e-8);
+  BOOST_CHECK_SMALL(
+    derivative_ref.tail(1).sum() - cumulativeBasisDerivative(4, knots, basis, 4, degree), 1e-8);
+
+  derivative_ref = computeDerivative(1, 6.);
+  deBoorBasis(degree, knots, 4, 6., basis);
+  BOOST_CHECK_SMALL(derivative_ref.sum(), 1e-8);
+  BOOST_CHECK_SMALL(
+    derivative_ref.tail(3).sum() - cumulativeBasisDerivative(4, knots, basis, 2, degree), 1e-8);
+  BOOST_CHECK_SMALL(
+    derivative_ref.tail(2).sum() - cumulativeBasisDerivative(4, knots, basis, 3, degree), 1e-8);
+  BOOST_CHECK_SMALL(
+    derivative_ref.tail(1).sum() - cumulativeBasisDerivative(4, knots, basis, 4, degree), 1e-8);
+
+  derivative_ref = computeDerivative(1, 7.3);
+  deBoorBasis(degree, knots, 4, 7.3, basis);
+  BOOST_CHECK_SMALL(derivative_ref.sum(), 1e-8);
+  BOOST_CHECK_SMALL(
+    derivative_ref.tail(3).sum() - cumulativeBasisDerivative(4, knots, basis, 2, degree), 1e-8);
+  BOOST_CHECK_SMALL(
+    derivative_ref.tail(2).sum() - cumulativeBasisDerivative(4, knots, basis, 3, degree), 1e-8);
+  BOOST_CHECK_SMALL(
+    derivative_ref.tail(1).sum() - cumulativeBasisDerivative(4, knots, basis, 4, degree), 1e-8);
+
+  derivative_ref = computeDerivative(1, 8.);
+  deBoorBasis(degree, knots, 4, 8., basis);
+  BOOST_CHECK_SMALL(derivative_ref.sum(), 1e-8);
+  BOOST_CHECK_SMALL(
+    derivative_ref.tail(3).sum() - cumulativeBasisDerivative(4, knots, basis, 2, degree), 1e-8);
+  BOOST_CHECK_SMALL(
+    derivative_ref.tail(2).sum() - cumulativeBasisDerivative(4, knots, basis, 3, degree), 1e-8);
+  BOOST_CHECK_SMALL(
+    derivative_ref.tail(1).sum() - cumulativeBasisDerivative(4, knots, basis, 4, degree), 1e-8);
+}
+
+// Test cumulativeBasisDerivative2 against bsplineBasisDerivative2
+BOOST_AUTO_TEST_CASE(cumulativeBasisDerivative2)
+{
+  using pinocchio::internal::cumulativeBasisDerivative2;
+  using pinocchio::internal::deBoorBasis;
+
+  const int degree = 3;
+  Eigen::VectorXd knots(9);
+  knots << 0., 2., 3., 5., 5.5, 8., 8.5, 10., 11.5;
+  Eigen::MatrixXd basis(Eigen::MatrixXd::Zero(4, 4));
+
+  auto computeDerivative2 = [knots](int start, double q) {
+    Eigen::VectorXd res(4);
+    for (int i = 0; i < 4; ++i)
+    {
+      res(i) = bsplineBasisDerivative2(i + start, degree, q, knots);
+    }
+    return res;
+  };
+
+  auto derivative2_ref = computeDerivative2(0, 5.);
+  deBoorBasis(degree, knots, 3, 5., basis);
+  BOOST_CHECK_SMALL(derivative2_ref.sum(), 1e-8);
+  BOOST_CHECK_SMALL(
+    derivative2_ref.tail(3).sum() - cumulativeBasisDerivative2(3, knots, basis, 1, degree), 1e-8);
+  BOOST_CHECK_SMALL(
+    derivative2_ref.tail(2).sum() - cumulativeBasisDerivative2(3, knots, basis, 2, degree), 1e-8);
+  BOOST_CHECK_SMALL(
+    derivative2_ref.tail(1).sum() - cumulativeBasisDerivative2(3, knots, basis, 3, degree), 1e-8);
+
+  derivative2_ref = computeDerivative2(0, 5.3);
+  deBoorBasis(degree, knots, 3, 5.3, basis);
+  BOOST_CHECK_SMALL(derivative2_ref.sum(), 1e-8);
+  BOOST_CHECK_SMALL(
+    derivative2_ref.tail(3).sum() - cumulativeBasisDerivative2(3, knots, basis, 1, degree), 1e-8);
+  BOOST_CHECK_SMALL(
+    derivative2_ref.tail(2).sum() - cumulativeBasisDerivative2(3, knots, basis, 2, degree), 1e-8);
+  BOOST_CHECK_SMALL(
+    derivative2_ref.tail(1).sum() - cumulativeBasisDerivative2(3, knots, basis, 3, degree), 1e-8);
+
+  derivative2_ref = computeDerivative2(1, 5.5);
+  deBoorBasis(degree, knots, 4, 5.5, basis);
+  BOOST_CHECK_SMALL(derivative2_ref.sum(), 1e-8);
+  BOOST_CHECK_SMALL(
+    derivative2_ref.tail(3).sum() - cumulativeBasisDerivative2(4, knots, basis, 2, degree), 1e-8);
+  BOOST_CHECK_SMALL(
+    derivative2_ref.tail(2).sum() - cumulativeBasisDerivative2(4, knots, basis, 3, degree), 1e-8);
+  BOOST_CHECK_SMALL(
+    derivative2_ref.tail(1).sum() - cumulativeBasisDerivative2(4, knots, basis, 4, degree), 1e-8);
+
+  derivative2_ref = computeDerivative2(1, 6.);
+  deBoorBasis(degree, knots, 4, 6., basis);
+  BOOST_CHECK_SMALL(derivative2_ref.sum(), 1e-8);
+  BOOST_CHECK_SMALL(
+    derivative2_ref.tail(3).sum() - cumulativeBasisDerivative2(4, knots, basis, 2, degree), 1e-8);
+  BOOST_CHECK_SMALL(
+    derivative2_ref.tail(2).sum() - cumulativeBasisDerivative2(4, knots, basis, 3, degree), 1e-8);
+  BOOST_CHECK_SMALL(
+    derivative2_ref.tail(1).sum() - cumulativeBasisDerivative2(4, knots, basis, 4, degree), 1e-8);
+
+  derivative2_ref = computeDerivative2(1, 7.3);
+  deBoorBasis(degree, knots, 4, 7.3, basis);
+  BOOST_CHECK_SMALL(derivative2_ref.sum(), 1e-8);
+  BOOST_CHECK_SMALL(
+    derivative2_ref.tail(3).sum() - cumulativeBasisDerivative2(4, knots, basis, 2, degree), 1e-8);
+  BOOST_CHECK_SMALL(
+    derivative2_ref.tail(2).sum() - cumulativeBasisDerivative2(4, knots, basis, 3, degree), 1e-8);
+  BOOST_CHECK_SMALL(
+    derivative2_ref.tail(1).sum() - cumulativeBasisDerivative2(4, knots, basis, 4, degree), 1e-8);
+
+  derivative2_ref = computeDerivative2(1, 8.);
+  deBoorBasis(degree, knots, 4, 8., basis);
+  BOOST_CHECK_SMALL(derivative2_ref.sum(), 1e-8);
+  BOOST_CHECK_SMALL(
+    derivative2_ref.tail(3).sum() - cumulativeBasisDerivative2(4, knots, basis, 2, degree), 1e-8);
+  BOOST_CHECK_SMALL(
+    derivative2_ref.tail(2).sum() - cumulativeBasisDerivative2(4, knots, basis, 3, degree), 1e-8);
+  BOOST_CHECK_SMALL(
+    derivative2_ref.tail(1).sum() - cumulativeBasisDerivative2(4, knots, basis, 4, degree), 1e-8);
+}
+
+// Test cumulativeBasisDerivativeFull against bsplineBasisDerivative.
+// This test is similar to cumulativeBasisDerivative.
+BOOST_AUTO_TEST_CASE(cumulativeBasisDerivativeFull)
+{
+  using pinocchio::internal::cumulativeBasisDerivativeFull;
+  using pinocchio::internal::deBoorFullBasis;
+
+  const int degree = 3;
+  Eigen::VectorXd knots(9);
+  knots << 0., 2., 3., 5., 5.5, 8., 8.5, 10., 11.5;
+  Eigen::MatrixXd basis(Eigen::MatrixXd::Zero(4, 5));
+
+  auto computeDerivative = [knots](double q) {
+    Eigen::VectorXd res(5);
+    for (int i = 0; i < 5; ++i)
+    {
+      res(i) = bsplineBasisDerivative(i, degree, q, knots);
+    }
+    return res;
+  };
+
+  auto derivative_ref = computeDerivative(5.);
+  deBoorFullBasis(degree, knots, 5., basis);
+  BOOST_CHECK_SMALL(
+    derivative_ref.sum() - cumulativeBasisDerivativeFull(degree, knots, basis, 0, degree), 1e-8);
+  BOOST_CHECK_SMALL(
+    derivative_ref.tail(4).sum() - cumulativeBasisDerivativeFull(degree, knots, basis, 1, degree),
+    1e-8);
+  BOOST_CHECK_SMALL(
+    derivative_ref.tail(3).sum() - cumulativeBasisDerivativeFull(degree, knots, basis, 2, degree),
+    1e-8);
+  BOOST_CHECK_SMALL(
+    derivative_ref.tail(2).sum() - cumulativeBasisDerivativeFull(degree, knots, basis, 3, degree),
+    1e-8);
+  BOOST_CHECK_SMALL(
+    derivative_ref.tail(1).sum() - cumulativeBasisDerivativeFull(degree, knots, basis, 4, degree),
+    1e-8);
+
+  derivative_ref = computeDerivative(5.3);
+  deBoorFullBasis(degree, knots, 5.3, basis);
+  BOOST_CHECK_SMALL(
+    derivative_ref.sum() - cumulativeBasisDerivativeFull(degree, knots, basis, 0, degree), 1e-8);
+  BOOST_CHECK_SMALL(
+    derivative_ref.tail(4).sum() - cumulativeBasisDerivativeFull(degree, knots, basis, 1, degree),
+    1e-8);
+  BOOST_CHECK_SMALL(
+    derivative_ref.tail(3).sum() - cumulativeBasisDerivativeFull(degree, knots, basis, 2, degree),
+    1e-8);
+  BOOST_CHECK_SMALL(
+    derivative_ref.tail(2).sum() - cumulativeBasisDerivativeFull(degree, knots, basis, 3, degree),
+    1e-8);
+  BOOST_CHECK_SMALL(
+    derivative_ref.tail(1).sum() - cumulativeBasisDerivativeFull(degree, knots, basis, 4, degree),
+    1e-8);
+
+  derivative_ref = computeDerivative(5.5);
+  deBoorFullBasis(degree, knots, 5.5, basis);
+  BOOST_CHECK_SMALL(
+    derivative_ref.sum() - cumulativeBasisDerivativeFull(degree, knots, basis, 0, degree), 1e-8);
+  BOOST_CHECK_SMALL(
+    derivative_ref.tail(4).sum() - cumulativeBasisDerivativeFull(degree, knots, basis, 1, degree),
+    1e-8);
+  BOOST_CHECK_SMALL(
+    derivative_ref.tail(3).sum() - cumulativeBasisDerivativeFull(degree, knots, basis, 2, degree),
+    1e-8);
+  BOOST_CHECK_SMALL(
+    derivative_ref.tail(2).sum() - cumulativeBasisDerivativeFull(degree, knots, basis, 3, degree),
+    1e-8);
+  BOOST_CHECK_SMALL(
+    derivative_ref.tail(1).sum() - cumulativeBasisDerivativeFull(degree, knots, basis, 4, degree),
+    1e-8);
+
+  derivative_ref = computeDerivative(6.);
+  deBoorFullBasis(degree, knots, 6., basis);
+  BOOST_CHECK_SMALL(
+    derivative_ref.sum() - cumulativeBasisDerivativeFull(degree, knots, basis, 0, degree), 1e-8);
+  BOOST_CHECK_SMALL(
+    derivative_ref.tail(4).sum() - cumulativeBasisDerivativeFull(degree, knots, basis, 1, degree),
+    1e-8);
+  BOOST_CHECK_SMALL(
+    derivative_ref.tail(3).sum() - cumulativeBasisDerivativeFull(degree, knots, basis, 2, degree),
+    1e-8);
+  BOOST_CHECK_SMALL(
+    derivative_ref.tail(2).sum() - cumulativeBasisDerivativeFull(degree, knots, basis, 3, degree),
+    1e-8);
+  BOOST_CHECK_SMALL(
+    derivative_ref.tail(1).sum() - cumulativeBasisDerivativeFull(degree, knots, basis, 4, degree),
+    1e-8);
+
+  derivative_ref = computeDerivative(7.3);
+  deBoorFullBasis(degree, knots, 7.3, basis);
+  BOOST_CHECK_SMALL(
+    derivative_ref.sum() - cumulativeBasisDerivativeFull(degree, knots, basis, 0, degree), 1e-8);
+  BOOST_CHECK_SMALL(
+    derivative_ref.tail(4).sum() - cumulativeBasisDerivativeFull(degree, knots, basis, 1, degree),
+    1e-8);
+  BOOST_CHECK_SMALL(
+    derivative_ref.tail(3).sum() - cumulativeBasisDerivativeFull(degree, knots, basis, 2, degree),
+    1e-8);
+  BOOST_CHECK_SMALL(
+    derivative_ref.tail(2).sum() - cumulativeBasisDerivativeFull(degree, knots, basis, 3, degree),
+    1e-8);
+  BOOST_CHECK_SMALL(
+    derivative_ref.tail(1).sum() - cumulativeBasisDerivativeFull(degree, knots, basis, 4, degree),
+    1e-8);
+
+  derivative_ref = computeDerivative(8.);
+  deBoorFullBasis(degree, knots, 8., basis);
+  BOOST_CHECK_SMALL(
+    derivative_ref.sum() - cumulativeBasisDerivativeFull(degree, knots, basis, 0, degree), 1e-8);
+  BOOST_CHECK_SMALL(
+    derivative_ref.tail(4).sum() - cumulativeBasisDerivativeFull(degree, knots, basis, 1, degree),
+    1e-8);
+  BOOST_CHECK_SMALL(
+    derivative_ref.tail(3).sum() - cumulativeBasisDerivativeFull(degree, knots, basis, 2, degree),
+    1e-8);
+  BOOST_CHECK_SMALL(
+    derivative_ref.tail(2).sum() - cumulativeBasisDerivativeFull(degree, knots, basis, 3, degree),
+    1e-8);
+  BOOST_CHECK_SMALL(
+    derivative_ref.tail(1).sum() - cumulativeBasisDerivativeFull(degree, knots, basis, 4, degree),
+    1e-8);
+}
+
+// Test cumulativeBasisDerivative2Full against bsplineBasisDerivative2
+// This test is similar to cumulativeBasisDerivative2.
+BOOST_AUTO_TEST_CASE(cumulativeBasisDerivativeFull2)
+{
+  using pinocchio::internal::cumulativeBasisDerivative2Full;
+  using pinocchio::internal::deBoorFullBasis;
+
+  const int degree = 3;
+  Eigen::VectorXd knots(9);
+  knots << 0., 2., 3., 5., 5.5, 8., 8.5, 10., 11.5;
+  Eigen::MatrixXd basis(Eigen::MatrixXd::Zero(4, 5));
+
+  auto computeDerivative2 = [knots](double q) {
+    Eigen::VectorXd res(5);
+    for (int i = 0; i < 5; ++i)
+    {
+      res(i) = bsplineBasisDerivative2(i, degree, q, knots);
+    }
+    return res;
+  };
+
+  auto derivative2_ref = computeDerivative2(5.);
+  deBoorFullBasis(degree, knots, 5., basis);
+  BOOST_CHECK_SMALL(
+    derivative2_ref.sum() - cumulativeBasisDerivative2Full(degree, knots, basis, 0, degree), 1e-8);
+  BOOST_CHECK_SMALL(
+    derivative2_ref.tail(4).sum() - cumulativeBasisDerivative2Full(degree, knots, basis, 1, degree),
+    1e-8);
+  BOOST_CHECK_SMALL(
+    derivative2_ref.tail(3).sum() - cumulativeBasisDerivative2Full(degree, knots, basis, 2, degree),
+    1e-8);
+  BOOST_CHECK_SMALL(
+    derivative2_ref.tail(2).sum() - cumulativeBasisDerivative2Full(degree, knots, basis, 3, degree),
+    1e-8);
+  BOOST_CHECK_SMALL(
+    derivative2_ref.tail(1).sum() - cumulativeBasisDerivative2Full(degree, knots, basis, 4, degree),
+    1e-8);
+
+  derivative2_ref = computeDerivative2(5.3);
+  deBoorFullBasis(degree, knots, 5.3, basis);
+  BOOST_CHECK_SMALL(
+    derivative2_ref.sum() - cumulativeBasisDerivative2Full(degree, knots, basis, 0, degree), 1e-8);
+  BOOST_CHECK_SMALL(
+    derivative2_ref.tail(4).sum() - cumulativeBasisDerivative2Full(degree, knots, basis, 1, degree),
+    1e-8);
+  BOOST_CHECK_SMALL(
+    derivative2_ref.tail(3).sum() - cumulativeBasisDerivative2Full(degree, knots, basis, 2, degree),
+    1e-8);
+  BOOST_CHECK_SMALL(
+    derivative2_ref.tail(2).sum() - cumulativeBasisDerivative2Full(degree, knots, basis, 3, degree),
+    1e-8);
+  BOOST_CHECK_SMALL(
+    derivative2_ref.tail(1).sum() - cumulativeBasisDerivative2Full(degree, knots, basis, 4, degree),
+    1e-8);
+
+  derivative2_ref = computeDerivative2(5.5);
+  deBoorFullBasis(degree, knots, 5.5, basis);
+  BOOST_CHECK_SMALL(
+    derivative2_ref.sum() - cumulativeBasisDerivative2Full(degree, knots, basis, 0, degree), 1e-8);
+  BOOST_CHECK_SMALL(
+    derivative2_ref.tail(4).sum() - cumulativeBasisDerivative2Full(degree, knots, basis, 1, degree),
+    1e-8);
+  BOOST_CHECK_SMALL(
+    derivative2_ref.tail(3).sum() - cumulativeBasisDerivative2Full(degree, knots, basis, 2, degree),
+    1e-8);
+  BOOST_CHECK_SMALL(
+    derivative2_ref.tail(2).sum() - cumulativeBasisDerivative2Full(degree, knots, basis, 3, degree),
+    1e-8);
+  BOOST_CHECK_SMALL(
+    derivative2_ref.tail(1).sum() - cumulativeBasisDerivative2Full(degree, knots, basis, 4, degree),
+    1e-8);
+
+  derivative2_ref = computeDerivative2(6.);
+  deBoorFullBasis(degree, knots, 6., basis);
+  BOOST_CHECK_SMALL(
+    derivative2_ref.sum() - cumulativeBasisDerivative2Full(degree, knots, basis, 0, degree), 1e-8);
+  BOOST_CHECK_SMALL(
+    derivative2_ref.tail(4).sum() - cumulativeBasisDerivative2Full(degree, knots, basis, 1, degree),
+    1e-8);
+  BOOST_CHECK_SMALL(
+    derivative2_ref.tail(3).sum() - cumulativeBasisDerivative2Full(degree, knots, basis, 2, degree),
+    1e-8);
+  BOOST_CHECK_SMALL(
+    derivative2_ref.tail(2).sum() - cumulativeBasisDerivative2Full(degree, knots, basis, 3, degree),
+    1e-8);
+  BOOST_CHECK_SMALL(
+    derivative2_ref.tail(1).sum() - cumulativeBasisDerivative2Full(degree, knots, basis, 4, degree),
+    1e-8);
+
+  derivative2_ref = computeDerivative2(7.3);
+  deBoorFullBasis(degree, knots, 7.3, basis);
+  BOOST_CHECK_SMALL(
+    derivative2_ref.sum() - cumulativeBasisDerivative2Full(degree, knots, basis, 0, degree), 1e-8);
+  BOOST_CHECK_SMALL(
+    derivative2_ref.tail(4).sum() - cumulativeBasisDerivative2Full(degree, knots, basis, 1, degree),
+    1e-8);
+  BOOST_CHECK_SMALL(
+    derivative2_ref.tail(3).sum() - cumulativeBasisDerivative2Full(degree, knots, basis, 2, degree),
+    1e-8);
+  BOOST_CHECK_SMALL(
+    derivative2_ref.tail(2).sum() - cumulativeBasisDerivative2Full(degree, knots, basis, 3, degree),
+    1e-8);
+  BOOST_CHECK_SMALL(
+    derivative2_ref.tail(1).sum() - cumulativeBasisDerivative2Full(degree, knots, basis, 4, degree),
+    1e-8);
+
+  derivative2_ref = computeDerivative2(8.);
+  deBoorFullBasis(degree, knots, 8., basis);
+  BOOST_CHECK_SMALL(
+    derivative2_ref.sum() - cumulativeBasisDerivative2Full(degree, knots, basis, 0, degree), 1e-8);
+  BOOST_CHECK_SMALL(
+    derivative2_ref.tail(4).sum() - cumulativeBasisDerivative2Full(degree, knots, basis, 1, degree),
+    1e-8);
+  BOOST_CHECK_SMALL(
+    derivative2_ref.tail(3).sum() - cumulativeBasisDerivative2Full(degree, knots, basis, 2, degree),
+    1e-8);
+  BOOST_CHECK_SMALL(
+    derivative2_ref.tail(2).sum() - cumulativeBasisDerivative2Full(degree, knots, basis, 3, degree),
+    1e-8);
+  BOOST_CHECK_SMALL(
+    derivative2_ref.tail(1).sum() - cumulativeBasisDerivative2Full(degree, knots, basis, 4, degree),
+    1e-8);
+}
+
+/// @brief Test to make sure the relative motions are correct
+BOOST_AUTO_TEST_CASE(relativeMotions)
+{
+  size_t degree = 3;
+
+  std::vector<SE3> ctrlFrames;
+  ctrlFrames.push_back(SE3::Identity());
+  std::vector<Motion> relativeMotions;
+  for (int i = 0; i < 10; i++)
+  {
+    relativeMotions.push_back(Motion::Random());
+    const SE3 & currentFrame = ctrlFrames.back();
+    ctrlFrames.push_back(currentFrame * exp6(relativeMotions.back()));
+  }
+
+  auto jmodel = JointModelSplineBuilder()
+                  .withControlFrameVector(ctrlFrames)
+                  .withDegree(degree)
+                  .withOpenUniformKnots(0., 1.)
+                  .build();
+  // Check size
+  BOOST_CHECK(jmodel.relativeMotions.size() == (ctrlFrames.size() - 1));
+
+  // check values
+  for (size_t i = 0; i < ctrlFrames.size() - 1; i++)
+    BOOST_CHECK(jmodel.relativeMotions[i].isApprox(relativeMotions[i]));
+}
+
+/// @brief Test findSpan on the simplest case (no redundant knot vector).
+BOOST_AUTO_TEST_CASE(findSpan_degree_0)
+{
+  using pinocchio::internal::findSpan;
+  const int degree = 0;
+  Eigen::VectorXd knotVector(6);
+  knotVector << 0., 0.2, 0.4, 0.6, 0.8, 1.;
+
+  BOOST_CHECK_EQUAL(findSpan(degree, knotVector, 0.), 0);
+  BOOST_CHECK_EQUAL(findSpan(degree, knotVector, 0.1), 0);
+  BOOST_CHECK_EQUAL(findSpan(degree, knotVector, 0.2), 1);
+  BOOST_CHECK_EQUAL(findSpan(degree, knotVector, 0.5), 2);
+  BOOST_CHECK_EQUAL(findSpan(degree, knotVector, 0.7), 3);
+  BOOST_CHECK_EQUAL(findSpan(degree, knotVector, 0.8), 4);
+  BOOST_CHECK_EQUAL(findSpan(degree, knotVector, 0.9), 4);
+  BOOST_CHECK_EQUAL(findSpan(degree, knotVector, 1.), 4);
+}
+
+/// @brief Test findSpan with redundant knot vector values.
+BOOST_AUTO_TEST_CASE(findSpan_degree_1)
+{
+  using pinocchio::internal::findSpan;
+  const int degree = 1;
+  Eigen::VectorXd knotVector(8);
+  knotVector << 0., 0., 0.2, 0.6, 0.6, 0.8, 1., 1.;
+
+  BOOST_CHECK_EQUAL(findSpan(degree, knotVector, 0.), 1);
+  BOOST_CHECK_EQUAL(findSpan(degree, knotVector, 0.1), 1);
+  BOOST_CHECK_EQUAL(findSpan(degree, knotVector, 0.2), 2);
+  BOOST_CHECK_EQUAL(findSpan(degree, knotVector, 0.5), 2);
+  BOOST_CHECK_EQUAL(findSpan(degree, knotVector, 0.6), 4);
+  BOOST_CHECK_EQUAL(findSpan(degree, knotVector, 0.7), 4);
+  BOOST_CHECK_EQUAL(findSpan(degree, knotVector, 0.8), 5);
+  BOOST_CHECK_EQUAL(findSpan(degree, knotVector, 0.9), 5);
+  BOOST_CHECK_EQUAL(findSpan(degree, knotVector, 1.), 5);
+}
+
+/// @brief Comparing a simple spline joint with a PZ.
+/// Make sure pose and joint subspace are the same.
+/// This test the internal::computeSplineKinematics function and
+/// we use this output to validate internal::computeSplineKinematicsFull
+/// at the same time.
+BOOST_AUTO_TEST_CASE(vsPrismaticZ)
+{
+  using namespace pinocchio;
+
+  // Spline Joint
+  std::vector<SE3> ctrlFrames;
+  ctrlFrames.push_back(SE3::Identity());
+  ctrlFrames.push_back(SE3(Eigen::Matrix3d::Identity(), Eigen::Vector3d(0., 0., 1. / 3.)));
+  ctrlFrames.push_back(SE3(Eigen::Matrix3d::Identity(), Eigen::Vector3d(0., 0., 2. / 3)));
+  ctrlFrames.push_back(SE3(Eigen::Matrix3d::Identity(), Eigen::Vector3d(0., 0., 1.)));
+
+  auto jmodel = JointModelSplineBuilder()
+                  .withControlFrameVector(ctrlFrames)
+                  .withDegree(3)
+                  .withOpenUniformKnots(0., 1.)
+                  .build();
+
+  JointDataSpline jdata = jmodel.createData();
+  jmodel.setIndexes(0, 0, 0);
+
+  // Prismatic joint
+  JointModelPZ jmodelPz;
+  JointDataPZ jdataPz = jmodelPz.createData();
+  jmodelPz.setIndexes(0, 0, 0);
+
+  // Data for computeSplineKinematicsFull
+  Eigen::MatrixXd basis_full(
+    Eigen::MatrixXd::Zero(jmodel.degree + 1, jmodel.knots.size() - jmodel.degree - 1));
+  SE3 M_full;
+  Motion v_full;
+  Motion c_full;
+  JointMotionSubspaceTpl<1, double, 0, 1> S_full;
+
+  Eigen::Matrix<double, 1, 1> q(Eigen::VectorXd::Zero(1));
+  Eigen::Matrix<double, 1, 1> q_dot(Eigen::VectorXd::Zero(1));
+
+  // -------
+  q << 0.2;
+
+  jmodel.calc(jdata, q);
+  jmodelPz.calc(jdataPz, q);
+  internal::computeSplineKinematicsFull(
+    jmodel.degree, jmodel.knots, jmodel.ctrlFrames, jmodel.relativeMotions, q[0], q_dot, false,
+    M_full, v_full, c_full, S_full, basis_full);
+
+  BOOST_CHECK(jdata.M.isApprox(jdataPz.M, 1e-12));
+  BOOST_CHECK(jdata.S.matrix().isApprox(jdataPz.S.matrix(), 1e-12));
+  BOOST_CHECK(M_full.isApprox(jdataPz.M, 1e-12));
+  BOOST_CHECK(S_full.matrix().isApprox(jdataPz.S.matrix(), 1e-12));
+
+  // -------
+  q << 0.3;
+  q_dot << 0.4;
+
+  jmodel.calc(jdata, q, q_dot);
+  jmodelPz.calc(jdataPz, q, q_dot);
+  internal::computeSplineKinematicsFull(
+    jmodel.degree, jmodel.knots, jmodel.ctrlFrames, jmodel.relativeMotions, q[0], q_dot, true,
+    M_full, v_full, c_full, S_full, basis_full);
+
+  BOOST_CHECK(jdata.M.isApprox(jdataPz.M, 1e-12));
+  BOOST_CHECK(jdata.S.matrix().isApprox(jdataPz.S.matrix(), 1e-12));
+  BOOST_CHECK(jdata.v.isApprox(jdataPz.v, 1e-12));
+  BOOST_CHECK(jdata.c.isApprox(jdataPz.c, 1e-12));
+  BOOST_CHECK(M_full.isApprox(jdataPz.M, 1e-12));
+  BOOST_CHECK(S_full.matrix().isApprox(jdataPz.S.matrix(), 1e-12));
+  BOOST_CHECK(v_full.isApprox(jdataPz.v, 1e-12));
+  BOOST_CHECK(c_full.isApprox(jdataPz.c, 1e-12));
+}
+
+/// @brief Comparing a simple spline joint with a RX.
+/// Make sure pose and joint subspace are the same.
+/// This test the internal::computeSplineKinematics function and
+/// we use this output to validate internal::computeSplineKinematicsFull
+/// at the same time.
+BOOST_AUTO_TEST_CASE(vsRevoluteX)
+{
+  using namespace pinocchio;
+
+  // Spline Joint
+  Eigen::Matrix3d rotation;
+  std::vector<SE3> ctrlFrames;
+  Eigen::AngleAxisd Rx1(1. / 3., Eigen::Vector3d::UnitX());
+  Eigen::AngleAxisd Rx2(2. / 3., Eigen::Vector3d::UnitX());
+  Eigen::AngleAxisd Rx3(1., Eigen::Vector3d::UnitX());
+
+  ctrlFrames.push_back(SE3::Identity());
+  ctrlFrames.push_back(SE3(Rx1.toRotationMatrix(), Eigen::Vector3d(0., 0., 0.)));
+  ctrlFrames.push_back(SE3(Rx2.toRotationMatrix(), Eigen::Vector3d(0., 0., 0.)));
+  ctrlFrames.push_back(SE3(Rx3.toRotationMatrix(), Eigen::Vector3d(0., 0., 0.)));
+
+  auto jmodel = JointModelSplineBuilder()
+                  .withControlFrameVector(ctrlFrames)
+                  .withDegree(3)
+                  .withOpenUniformKnots(0., 1.)
+                  .build();
+  JointDataSpline jdata = jmodel.createData();
+  jmodel.setIndexes(0, 0, 0);
+
+  // Prismatic joint
+  JointModelRX jmodelRx;
+  JointDataRX jdataRx = jmodelRx.createData();
+  jmodelRx.setIndexes(0, 0, 0);
+
+  // Data for computeSplineKinematicsFull
+  Eigen::MatrixXd basis_full(
+    Eigen::MatrixXd::Zero(jmodel.degree + 1, jmodel.knots.size() - jmodel.degree - 1));
+  SE3 M_full;
+  Motion v_full;
+  Motion c_full;
+  JointMotionSubspaceTpl<1, double, 0, 1> S_full;
+
+  Eigen::Matrix<double, 1, 1> q(Eigen::VectorXd::Zero(1));
+  Eigen::Matrix<double, 1, 1> q_dot(Eigen::VectorXd::Zero(1));
+
+  // -------
+  q << 0.2;
+
+  jmodel.calc(jdata, q);
+  jmodelRx.calc(jdataRx, q);
+  internal::computeSplineKinematicsFull(
+    jmodel.degree, jmodel.knots, jmodel.ctrlFrames, jmodel.relativeMotions, q[0], q_dot, false,
+    M_full, v_full, c_full, S_full, basis_full);
+
+  BOOST_CHECK(jdata.M.isApprox(jdataRx.M, 1e-12));
+  BOOST_CHECK(jdata.S.matrix().isApprox(jdataRx.S.matrix(), 1e-12));
+  BOOST_CHECK(M_full.isApprox(jdataRx.M, 1e-12));
+  BOOST_CHECK(S_full.matrix().isApprox(jdataRx.S.matrix(), 1e-12));
+
+  // -------
+  q << 0.3;
+  q_dot << 0.4;
+
+  jmodel.calc(jdata, q, q_dot);
+  jmodelRx.calc(jdataRx, q, q_dot);
+  internal::computeSplineKinematicsFull(
+    jmodel.degree, jmodel.knots, jmodel.ctrlFrames, jmodel.relativeMotions, q[0], q_dot, true,
+    M_full, v_full, c_full, S_full, basis_full);
+
+  BOOST_CHECK(jdata.M.isApprox(jdataRx.M, 1e-12));
+  BOOST_CHECK(jdata.S.matrix().isApprox(jdataRx.S.matrix(), 1e-12));
+  BOOST_CHECK(jdata.v.isApprox(jdataRx.v, 1e-12));
+  BOOST_CHECK(jdata.c.isApprox(jdataRx.c, 1e-12));
+  BOOST_CHECK(M_full.isApprox(jdataRx.M, 1e-12));
+  BOOST_CHECK(S_full.matrix().isApprox(jdataRx.S.matrix(), 1e-12));
+  BOOST_CHECK(v_full.isApprox(jdataRx.v, 1e-12));
+  BOOST_CHECK(c_full.isApprox(jdataRx.c, 1e-12));
+}
+
+/// @brief Test out rnea vs aba
+BOOST_AUTO_TEST_CASE(abaVSrnea)
+{
+  typedef SE3::Vector3 Vector3;
+  typedef SE3::Matrix3 Matrix3;
+
+  Model modelSpline;
+
+  Inertia inertia(1., Vector3(0., 0., 0.0), Matrix3::Identity());
+
+  std::vector<SE3> ctrlFrames;
+  getTrajectory(ctrlFrames);
+  addJointAndBody(
+    modelSpline,
+    JointModelSplineBuilder()
+      .withControlFrameVector(ctrlFrames)
+      .withDegree(3)
+      .withOpenUniformKnots(0., 1.)
+      .build(),
+    0, SE3::Identity(), "kneeSpline", inertia);
+  Data dataSplineRnea(modelSpline);
+  Data dataSplineAba(modelSpline);
+
+  Eigen::VectorXd q(modelSpline.nq);
+  pinocchio::randomConfiguration(
+    modelSpline, Eigen::VectorXd::Zero(1), Eigen::VectorXd::Ones(1), q);
+  Eigen::VectorXd vq(Eigen::VectorXd::Random(modelSpline.nv));
+  Eigen::VectorXd aq(Eigen::VectorXd::Random(modelSpline.nv));
+  Eigen::VectorXd tauRnea(1);
+  Eigen::VectorXd aAba(1);
+
+  tauRnea = pinocchio::rnea(modelSpline, dataSplineRnea, q, vq, aq);
+  aAba = pinocchio::aba(modelSpline, dataSplineAba, q, vq, tauRnea);
+
+  BOOST_CHECK(aq.isApprox(aAba));
+}
+
+/// @brief Test S and bias c computation via finite differences.
+/// This test the internal::computeSplineKinematics function and
+/// we use this output to validate internal::computeSplineKinematicsFull
+/// at the same time.
+BOOST_AUTO_TEST_CASE(vsFiniteDifference)
+{
+  using namespace pinocchio;
+
+  typedef typename JointModelSpline::ConfigVector_t CV;
+  typedef typename JointModelSpline::TangentVector_t TV;
+  typedef typename LieGroup<JointModelSpline>::type LieGroupType;
+
+  std::vector<SE3> ctrlFrames;
+  getTrajectory(ctrlFrames);
+
+  auto jmodel = JointModelSplineBuilder()
+                  .withControlFrameVector(ctrlFrames)
+                  .withDegree(3)
+                  .withOpenUniformKnots(0., 1.)
+                  .build();
+  JointDataSpline jdata = jmodel.createData();
+
+  jmodel.setIndexes(0, 0, 0);
+
+  // Data for computeSplineKinematicsFull
+  Eigen::MatrixXd basis_full(
+    Eigen::MatrixXd::Zero(jmodel.degree + 1, jmodel.knots.size() - jmodel.degree - 1));
+  SE3 M_full;
+  Motion v_full;
+  Motion c_full;
+  JointMotionSubspaceTpl<1, double, 0, 1> S_full;
+
+  const double eps = 1e-8;
+  CV q_ref(1);
+  q_ref[0] = 0.6;
+  CV q(q_ref);
+
+  const Eigen::DenseIndex nv = jdata.S.nv();
+  TV q_dot(nv);
+  TV q_dot_ref(nv);
+  q_dot.setZero();
+
+  q_dot[0] = eps;
+  q_dot_ref[0] = 0.3;
+
+  q = LieGroupType().integrate(q_ref, q_dot);
+  // Check S
+  {
+    jmodel.calc(jdata, q_ref);
+    internal::computeSplineKinematicsFull(
+      jmodel.degree, jmodel.knots, jmodel.ctrlFrames, jmodel.relativeMotions, q_ref[0], q_dot_ref,
+      false, M_full, v_full, c_full, S_full, basis_full);
+
+    SE3 M_ref(jdata.M);
+    Eigen::VectorXd S(6, JointModelSpline::NV);
+    Eigen::VectorXd S_ref(jdata.S.matrix());
+    Eigen::VectorXd S_ref_full(S_full.matrix());
+
+    jmodel.calc(jdata, q);
+    SE3 M_ = jdata.M;
+
+    S.col(0) = log6(M_ref.inverse() * M_).toVector();
+    S.col(0) /= eps;
+
+    BOOST_CHECK(S.isApprox(S_ref, 1e-6));
+    BOOST_CHECK(S.isApprox(S_ref_full, 1e-6));
+  }
+  // Check bias
+  {
+    jmodel.calc(jdata, q_ref, q_dot_ref);
+    internal::computeSplineKinematicsFull(
+      jmodel.degree, jmodel.knots, jmodel.ctrlFrames, jmodel.relativeMotions, q_ref[0], q_dot_ref,
+      true, M_full, v_full, c_full, S_full, basis_full);
+    const Motion & c_ref = jdata.c;
+    Eigen::VectorXd S_ref(jdata.S.matrix());
+
+    jmodel.calc(jdata, q);
+    Eigen::VectorXd S_(jdata.S.matrix());
+
+    Motion dSdq_fd((S_ - S_ref) / eps);
+    Motion c_fd = dSdq_fd * q_dot_ref[0] * q_dot_ref[0];
+
+    BOOST_CHECK(c_ref.isApprox(c_fd, 1e-6));
+    BOOST_CHECK(c_full.isApprox(c_fd, 1e-6));
+  }
+}
+
+BOOST_AUTO_TEST_SUITE_END()
diff --git a/unittest/model-graph.cpp b/unittest/model-graph.cpp
index 9fc94b3529..81708972f6 100644
--- a/unittest/model-graph.cpp
+++ b/unittest/model-graph.cpp
@@ -400,6 +400,53 @@ BOOST_AUTO_TEST_CASE(test_mimic_joint)
   BOOST_CHECK(d.oMf[m.getFrameId("body3", pinocchio::BODY)].isApprox(bodyPose));
 }
 
+/// @brief test spline joint, and that we cannot reverse it
+BOOST_AUTO_TEST_CASE(test_spline)
+{
+  using namespace pinocchio::graph;
+
+  ModelGraph g;
+  //////////////////////////////////////// Bodies
+  g.addBody("body1", pinocchio::Inertia::Identity());
+  g.addFrame("body2", BodyFrame(pinocchio::Inertia::Identity()));
+
+  /////////////////////////////////////// Joints
+  pinocchio::SE3 pose_body1_joint1(Eigen::Matrix3d::Identity(), Eigen::Vector3d(2., 0., 0.));
+  pinocchio::SE3 pose_body2_joint1(Eigen::Matrix3d::Identity(), Eigen::Vector3d(0., 4., 0.));
+
+  pinocchio::SE3 finalPose(pinocchio::SE3::Random());
+  std::vector<pinocchio::SE3> ctrlFrames;
+  ctrlFrames.push_back(pinocchio::SE3::Identity());
+  ctrlFrames.push_back(finalPose);
+
+  g.edgeBuilder()
+    .withName("body1_to_body2")
+    .withSourceVertex("body1")
+    .withSourcePose(pose_body1_joint1)
+    .withTargetVertex("body2")
+    .withTargetPose(pose_body2_joint1)
+    .withJointType(
+      JointSplineBuilder()
+        .withDegree(1)
+        .withControlFrameVector(ctrlFrames)
+        .withOpenUniformKnots(0., 1.)
+        .build())
+    .build();
+
+  pinocchio::SE3 pose_body1_universe = pinocchio::SE3::Identity();
+  pinocchio::Model m = buildModel(g, "body1", pose_body1_universe);
+  pinocchio::Data data(m);
+
+  Eigen::VectorXd q = Eigen::VectorXd::Ones(m.nq);
+  pinocchio::framesForwardKinematics(m, data, q);
+
+  BOOST_CHECK(data.oMf[m.getFrameId("body2", pinocchio::BODY)].isApprox(
+    pose_body1_joint1 * finalPose * pose_body2_joint1));
+
+  ///////////////// Reverse Model
+  BOOST_CHECK_THROW(buildModel(g, "body2", pinocchio::SE3::Identity()), std::invalid_argument);
+}
+
 /// @brief test out reverse joint for revolute
 BOOST_AUTO_TEST_CASE(test_reverse_revolute)
 {
diff --git a/unittest/serialization-multibody-joint.cpp b/unittest/serialization-multibody-joint.cpp
index fd02447980..17f99f9f11 100644
--- a/unittest/serialization-multibody-joint.cpp
+++ b/unittest/serialization-multibody-joint.cpp
@@ -16,181 +16,17 @@
 #include <boost/test/unit_test.hpp>
 #include <boost/utility/binary.hpp>
 
-BOOST_AUTO_TEST_SUITE(BOOST_TEST_MODULE)
-
-template<typename JointModel_>
-struct init;
-
-template<typename JointModel_>
-struct init
-{
-  static JointModel_ run()
-  {
-    JointModel_ jmodel;
-    jmodel.setIndexes(0, 0, 0);
-    return jmodel;
-  }
-};
-
-template<typename Scalar, int Options>
-struct init<pinocchio::JointModelRevoluteUnalignedTpl<Scalar, Options>>
-{
-  typedef pinocchio::JointModelRevoluteUnalignedTpl<Scalar, Options> JointModel;
-
-  static JointModel run()
-  {
-    typedef typename JointModel::Vector3 Vector3;
-    JointModel jmodel(Vector3::Random().normalized());
-
-    jmodel.setIndexes(0, 0, 0);
-    return jmodel;
-  }
-};
-
-template<typename Scalar, int Options>
-struct init<pinocchio::JointModelRevoluteUnboundedUnalignedTpl<Scalar, Options>>
-{
-  typedef pinocchio::JointModelRevoluteUnboundedUnalignedTpl<Scalar, Options> JointModel;
-
-  static JointModel run()
-  {
-    typedef typename JointModel::Vector3 Vector3;
-    JointModel jmodel(Vector3::Random().normalized());
-
-    jmodel.setIndexes(0, 0, 0);
-    return jmodel;
-  }
-};
-
-template<typename Scalar, int Options>
-struct init<pinocchio::JointModelPrismaticUnalignedTpl<Scalar, Options>>
-{
-  typedef pinocchio::JointModelPrismaticUnalignedTpl<Scalar, Options> JointModel;
-
-  static JointModel run()
-  {
-    typedef typename JointModel::Vector3 Vector3;
-    JointModel jmodel(Vector3::Random().normalized());
-
-    jmodel.setIndexes(0, 0, 0);
-    return jmodel;
-  }
-};
-
-template<typename Scalar, int Options>
-struct init<pinocchio::JointModelHelicalUnalignedTpl<Scalar, Options>>
-{
-  typedef pinocchio::JointModelHelicalUnalignedTpl<Scalar, Options> JointModel;
-
-  static JointModel run()
-  {
-    typedef typename JointModel::Vector3 Vector3;
-    JointModel jmodel(Vector3::Random().normalized());
-
-    jmodel.setIndexes(0, 0, 0);
-    return jmodel;
-  }
-};
-
-template<typename Scalar, int Options, int axis>
-struct init<pinocchio::JointModelHelicalTpl<Scalar, Options, axis>>
-{
-  typedef pinocchio::JointModelHelicalTpl<Scalar, Options, axis> JointModel;
-
-  static JointModel run()
-  {
-    JointModel jmodel(static_cast<Scalar>(0.0));
-
-    jmodel.setIndexes(0, 0, 0);
-    return jmodel;
-  }
-};
-
-template<typename Scalar, int Options>
-struct init<pinocchio::JointModelEllipsoidTpl<Scalar, Options>>
-{
-  typedef pinocchio::JointModelEllipsoidTpl<Scalar, Options> JointModel;
+#include "utils/joints-init.hpp"
 
-  static JointModel run()
-  {
-    JointModel jmodel(
-      static_cast<Scalar>(0.01), static_cast<Scalar>(0.02), static_cast<Scalar>(0.03));
-
-    jmodel.setIndexes(0, 0, 0);
-    return jmodel;
-  }
-};
-
-template<typename Scalar, int Options, template<typename, int> class JointCollection>
-struct init<pinocchio::JointModelTpl<Scalar, Options, JointCollection>>
-{
-  typedef pinocchio::JointModelTpl<Scalar, Options, JointCollection> JointModel;
-
-  static JointModel run()
-  {
-    typedef pinocchio::JointModelRevoluteTpl<Scalar, Options, 0> JointModelRX;
-    JointModel jmodel((JointModelRX()));
-
-    jmodel.setIndexes(0, 0, 0);
-    return jmodel;
-  }
-};
-
-template<typename Scalar, int Options>
-struct init<pinocchio::JointModelUniversalTpl<Scalar, Options>>
-{
-  typedef pinocchio::JointModelUniversalTpl<Scalar, Options> JointModel;
-
-  static JointModel run()
-  {
-    JointModel jmodel(pinocchio::XAxis::vector(), pinocchio::YAxis::vector());
-
-    jmodel.setIndexes(0, 0, 0);
-    return jmodel;
-  }
-};
-
-template<typename Scalar, int Options, template<typename, int> class JointCollection>
-struct init<pinocchio::JointModelCompositeTpl<Scalar, Options, JointCollection>>
-{
-  typedef pinocchio::JointModelCompositeTpl<Scalar, Options, JointCollection> JointModel;
-
-  static JointModel run()
-  {
-    typedef pinocchio::JointModelRevoluteTpl<Scalar, Options, 0> JointModelRX;
-    typedef pinocchio::JointModelRevoluteTpl<Scalar, Options, 1> JointModelRY;
-    JointModel jmodel((JointModelRX()));
-    jmodel.addJoint(JointModelRY());
-
-    jmodel.setIndexes(0, 0, 0);
-    return jmodel;
-  }
-};
-
-template<typename Scalar, int Options, template<typename, int> class JointCollection>
-struct init<pinocchio::JointModelMimicTpl<Scalar, Options, JointCollection>>
-{
-  typedef pinocchio::JointModelMimicTpl<Scalar, Options, JointCollection> JointModel;
-
-  static JointModel run()
-  {
-    typedef pinocchio::JointModelRevoluteTpl<Scalar, Options, 0> JointModelRX;
-    JointModelRX jmodel_ref = init<JointModelRX>::run();
-
-    JointModel jmodel(jmodel_ref, 1., 0.);
-    jmodel.setIndexes(1, 0, 0, 0);
-
-    return jmodel;
-  }
-};
+BOOST_AUTO_TEST_SUITE(BOOST_TEST_MODULE)
 
 struct TestJointModel
 {
   template<typename JointModel>
   void operator()(const pinocchio::JointModelBase<JointModel> &) const
   {
-    JointModel jmodel = init<JointModel>::run();
-    test(jmodel);
+    auto jmodelWithParameters = pinocchio::init<JointModel>::run();
+    test(jmodelWithParameters.jmodel);
   }
 
   template<typename JointType>
@@ -216,7 +52,8 @@ struct TestJointTransform
     typedef typename pinocchio::traits<JointDerived>::Constraint_t Constraint;
     typedef typename pinocchio::traits<JointDerived>::JointDataDerived JointData;
     typedef pinocchio::JointDataBase<JointData> JointDataBase;
-    JointModel jmodel = init<JointModel>::run();
+    auto JointModelWithParameters = pinocchio::init<JointModel>::run();
+    JointModel jmodel = JointModelWithParameters.jmodel;
 
     JointData jdata = jmodel.createData();
     JointDataBase & jdata_base = static_cast<JointDataBase &>(jdata);
@@ -224,10 +61,8 @@ struct TestJointTransform
     typedef typename pinocchio::LieGroup<JointModel>::type LieGroupType;
     LieGroupType lg;
 
-    Eigen::VectorXd lb(Eigen::VectorXd::Constant(jmodel.nq(), -1.));
-    Eigen::VectorXd ub(Eigen::VectorXd::Constant(jmodel.nq(), 1.));
-
-    Eigen::VectorXd q_random = lg.randomConfiguration(lb, ub);
+    Eigen::VectorXd q_random =
+      lg.randomConfiguration(JointModelWithParameters.lb, JointModelWithParameters.ub);
 
     jmodel.calc(jdata, q_random);
     Transform & m = jdata_base.M();
@@ -272,7 +107,8 @@ struct TestJointMotion
     typedef typename pinocchio::traits<JointDerived>::Bias_t Bias;
     typedef typename pinocchio::traits<JointDerived>::JointDataDerived JointData;
     typedef pinocchio::JointDataBase<JointData> JointDataBase;
-    JointModel jmodel = init<JointModel>::run();
+    auto JointModelWithParameters = pinocchio::init<JointModel>::run();
+    JointModel jmodel = JointModelWithParameters.jmodel;
 
     JointData jdata = jmodel.createData();
     JointDataBase & jdata_base = static_cast<JointDataBase &>(jdata);
@@ -280,10 +116,8 @@ struct TestJointMotion
     typedef typename pinocchio::LieGroup<JointModel>::type LieGroupType;
     LieGroupType lg;
 
-    Eigen::VectorXd lb(Eigen::VectorXd::Constant(jmodel.nq(), -1.));
-    Eigen::VectorXd ub(Eigen::VectorXd::Constant(jmodel.nq(), 1.));
-
-    Eigen::VectorXd q_random = lg.randomConfiguration(lb, ub);
+    Eigen::VectorXd q_random =
+      lg.randomConfiguration(JointModelWithParameters.lb, JointModelWithParameters.ub);
     Eigen::VectorXd v_random = Eigen::VectorXd::Random(jmodel.nv());
 
     jmodel.calc(jdata, q_random, v_random);
@@ -327,17 +161,16 @@ struct TestJointData
   {
     typedef typename JointModel::JointDerived JointDerived;
     typedef typename pinocchio::traits<JointDerived>::JointDataDerived JointData;
-    JointModel jmodel = init<JointModel>::run();
+    auto JointModelWithParameters = pinocchio::init<JointModel>::run();
+    JointModel jmodel = JointModelWithParameters.jmodel;
 
     JointData jdata = jmodel.createData();
 
     typedef typename pinocchio::LieGroup<JointModel>::type LieGroupType;
     LieGroupType lg;
 
-    Eigen::VectorXd lb(Eigen::VectorXd::Constant(jmodel.nq(), -1.));
-    Eigen::VectorXd ub(Eigen::VectorXd::Constant(jmodel.nq(), 1.));
-
-    Eigen::VectorXd q_random = lg.randomConfiguration(lb, ub);
+    Eigen::VectorXd q_random =
+      lg.randomConfiguration(JointModelWithParameters.lb, JointModelWithParameters.ub);
     Eigen::VectorXd v_random = Eigen::VectorXd::Random(jmodel.nv());
 
     jmodel.calc(jdata, q_random, v_random);
@@ -352,8 +185,8 @@ struct TestJointData
     typedef pinocchio::JointModelCompositeTpl<Scalar, Options, JointCollectionTpl> JointModel;
     typedef typename JointModel::JointDerived JointDerived;
     typedef typename pinocchio::traits<JointDerived>::JointDataDerived JointData;
-
-    JointModel jmodel_build = init<JointModel>::run();
+    auto JointModelWithParameters = pinocchio::init<JointModel>::run();
+    JointModel jmodel_build = JointModelWithParameters.jmodel;
 
     pinocchio::Model model;
     model.addJoint(0, jmodel_build, pinocchio::SE3::Random(), "model");
@@ -381,8 +214,9 @@ struct TestJointData
     typedef pinocchio::JointModelMimicTpl<Scalar, Options, JointCollection> JointModel;
     typedef typename JointModel::JointDerived JointDerived;
     typedef typename pinocchio::traits<JointDerived>::JointDataDerived JointData;
+    auto JointModelWithParameters = pinocchio::init<JointModel>::run();
+    JointModel jmodel = JointModelWithParameters.jmodel;
 
-    JointModel jmodel = init<JointModel>::run();
     JointData jdata = jmodel.createData();
 
     Eigen::VectorXd q_random = Eigen::VectorXd::Random(jmodel.jmodel().nq());
diff --git a/unittest/utils/joints-init.hpp b/unittest/utils/joints-init.hpp
new file mode 100644
index 0000000000..fa41a7f883
--- /dev/null
+++ b/unittest/utils/joints-init.hpp
@@ -0,0 +1,274 @@
+//
+// Copyright (c) 2025 INRIA
+//
+
+#include <vector>
+
+#include "pinocchio/multibody/joint.hpp"
+
+namespace pinocchio
+{
+  template<typename JointModel_>
+  struct JointModelWithParameters
+  {
+    typedef typename JointModel_::ConfigVector_t ConfigVector;
+    JointModel_ jmodel;
+    ConfigVector lb;
+    ConfigVector ub;
+  };
+
+  template<typename JointModel_>
+  struct init
+  {
+    static JointModelWithParameters<JointModel_> run()
+    {
+      JointModel_ jmodel;
+      jmodel.setIndexes(0, 0, 0);
+
+      // Default generic bounds
+      typename JointModel_::ConfigVector_t lb =
+        JointModel_::ConfigVector_t::Constant(jmodel.nq(), -1.0);
+      typename JointModel_::ConfigVector_t ub =
+        JointModel_::ConfigVector_t::Constant(jmodel.nq(), 1.0);
+
+      return {jmodel, lb, ub};
+    }
+  };
+
+  template<typename Scalar, int Options>
+  struct init<pinocchio::JointModelRevoluteUnalignedTpl<Scalar, Options>>
+  {
+    typedef pinocchio::JointModelRevoluteUnalignedTpl<Scalar, Options> JointModel;
+    typedef JointModelWithParameters<JointModel> ReturnType;
+
+    static ReturnType run()
+    {
+      typedef typename JointModel::Vector3 Vector3;
+      JointModel jmodel(Vector3::Random().normalized());
+
+      jmodel.setIndexes(0, 0, 0);
+      typename JointModel::ConfigVector_t lb =
+        JointModel::ConfigVector_t::Constant(jmodel.nq(), -1.0);
+      typename JointModel::ConfigVector_t ub =
+        JointModel::ConfigVector_t::Constant(jmodel.nq(), 1.0);
+
+      return {jmodel, lb, ub};
+    }
+  };
+
+  template<typename Scalar, int Options>
+  struct init<pinocchio::JointModelRevoluteUnboundedUnalignedTpl<Scalar, Options>>
+  {
+    typedef pinocchio::JointModelRevoluteUnboundedUnalignedTpl<Scalar, Options> JointModel;
+    typedef JointModelWithParameters<JointModel> ReturnType;
+
+    static ReturnType run()
+    {
+      typedef typename JointModel::Vector3 Vector3;
+      JointModel jmodel(Vector3::Random().normalized());
+
+      jmodel.setIndexes(0, 0, 0);
+      typename JointModel::ConfigVector_t lb =
+        JointModel::ConfigVector_t::Constant(jmodel.nq(), -1.0);
+      typename JointModel::ConfigVector_t ub =
+        JointModel::ConfigVector_t::Constant(jmodel.nq(), 1.0);
+
+      return {jmodel, lb, ub};
+    }
+  };
+
+  template<typename Scalar, int Options>
+  struct init<pinocchio::JointModelPrismaticUnalignedTpl<Scalar, Options>>
+  {
+    typedef pinocchio::JointModelPrismaticUnalignedTpl<Scalar, Options> JointModel;
+    typedef JointModelWithParameters<JointModel> ReturnType;
+
+    static ReturnType run()
+    {
+      typedef typename JointModel::Vector3 Vector3;
+      JointModel jmodel(Vector3::Random().normalized());
+
+      jmodel.setIndexes(0, 0, 0);
+      typename JointModel::ConfigVector_t lb =
+        JointModel::ConfigVector_t::Constant(jmodel.nq(), -1.0);
+      typename JointModel::ConfigVector_t ub =
+        JointModel::ConfigVector_t::Constant(jmodel.nq(), 1.0);
+
+      return {jmodel, lb, ub};
+    }
+  };
+
+  template<typename Scalar, int Options>
+  struct init<pinocchio::JointModelHelicalUnalignedTpl<Scalar, Options>>
+  {
+    typedef pinocchio::JointModelHelicalUnalignedTpl<Scalar, Options> JointModel;
+    typedef JointModelWithParameters<JointModel> ReturnType;
+
+    static ReturnType run()
+    {
+      typedef typename JointModel::Vector3 Vector3;
+      JointModel jmodel(Vector3::Random().normalized());
+
+      jmodel.setIndexes(0, 0, 0);
+      typename JointModel::ConfigVector_t lb =
+        JointModel::ConfigVector_t::Constant(jmodel.nq(), -1.0);
+      typename JointModel::ConfigVector_t ub =
+        JointModel::ConfigVector_t::Constant(jmodel.nq(), 1.0);
+
+      return {jmodel, lb, ub};
+    }
+  };
+
+  template<typename Scalar, int Options>
+  struct init<pinocchio::JointModelUniversalTpl<Scalar, Options>>
+  {
+    typedef pinocchio::JointModelUniversalTpl<Scalar, Options> JointModel;
+    typedef JointModelWithParameters<JointModel> ReturnType;
+
+    static ReturnType run()
+    {
+      JointModel jmodel(XAxis::vector(), YAxis::vector());
+
+      jmodel.setIndexes(0, 0, 0);
+      typename JointModel::ConfigVector_t lb =
+        JointModel::ConfigVector_t::Constant(jmodel.nq(), -1.0);
+      typename JointModel::ConfigVector_t ub =
+        JointModel::ConfigVector_t::Constant(jmodel.nq(), 1.0);
+
+      return {jmodel, lb, ub};
+    }
+  };
+
+  template<typename Scalar, int Options, int axis>
+  struct init<pinocchio::JointModelHelicalTpl<Scalar, Options, axis>>
+  {
+    typedef pinocchio::JointModelHelicalTpl<Scalar, Options, axis> JointModel;
+    typedef JointModelWithParameters<JointModel> ReturnType;
+
+    static ReturnType run()
+    {
+      JointModel jmodel(static_cast<Scalar>(0.5));
+
+      jmodel.setIndexes(0, 0, 0);
+      typename JointModel::ConfigVector_t lb =
+        JointModel::ConfigVector_t::Constant(jmodel.nq(), -1.0);
+      typename JointModel::ConfigVector_t ub =
+        JointModel::ConfigVector_t::Constant(jmodel.nq(), 1.0);
+
+      return {jmodel, lb, ub};
+    }
+  };
+
+  template<typename Scalar, int Options>
+  struct init<pinocchio::JointModelEllipsoidTpl<Scalar, Options>>
+  {
+    typedef pinocchio::JointModelEllipsoidTpl<Scalar, Options> JointModel;
+    typedef JointModelWithParameters<JointModel> ReturnType;
+
+    static ReturnType run()
+    {
+      JointModel jmodel(Scalar(0.01), Scalar(0.02), Scalar(0.03));
+
+      jmodel.setIndexes(0, 0, 0);
+      typename JointModel::ConfigVector_t lb =
+        JointModel::ConfigVector_t::Constant(jmodel.nq(), -1.0);
+      typename JointModel::ConfigVector_t ub =
+        JointModel::ConfigVector_t::Constant(jmodel.nq(), 1.0);
+
+      return {jmodel, lb, ub};
+    }
+  };
+
+  template<typename Scalar, int Options, template<typename, int> class JointCollection>
+  struct init<pinocchio::JointModelTpl<Scalar, Options, JointCollection>>
+  {
+    typedef pinocchio::JointModelTpl<Scalar, Options, JointCollection> JointModel;
+    typedef JointModelWithParameters<JointModel> ReturnType;
+
+    static ReturnType run()
+    {
+      typedef pinocchio::JointModelRevoluteTpl<Scalar, Options, 0> JointModelRX;
+      JointModel jmodel((JointModelRX()));
+
+      jmodel.setIndexes(0, 0, 0);
+      typename JointModel::ConfigVector_t lb =
+        JointModel::ConfigVector_t::Constant(jmodel.nq(), -1.0);
+      typename JointModel::ConfigVector_t ub =
+        JointModel::ConfigVector_t::Constant(jmodel.nq(), 1.0);
+
+      return {jmodel, lb, ub};
+    }
+  };
+
+  template<typename Scalar, int Options, template<typename, int> class JointCollection>
+  struct init<pinocchio::JointModelCompositeTpl<Scalar, Options, JointCollection>>
+  {
+    typedef pinocchio::JointModelCompositeTpl<Scalar, Options, JointCollection> JointModel;
+    typedef JointModelWithParameters<JointModel> ReturnType;
+
+    static ReturnType run()
+    {
+      typedef pinocchio::JointModelRevoluteTpl<Scalar, Options, 0> JointModelRX;
+      typedef pinocchio::JointModelRevoluteTpl<Scalar, Options, 1> JointModelRY;
+      JointModel jmodel((JointModelRX()));
+      jmodel.addJoint(JointModelRY());
+
+      jmodel.setIndexes(0, 0, 0);
+      typename JointModel::ConfigVector_t lb =
+        JointModel::ConfigVector_t::Constant(jmodel.nq(), -1.0);
+      typename JointModel::ConfigVector_t ub =
+        JointModel::ConfigVector_t::Constant(jmodel.nq(), 1.0);
+      return {jmodel, lb, ub};
+    }
+  };
+
+  template<typename Scalar, int Options, template<typename, int> class JointCollection>
+  struct init<pinocchio::JointModelMimicTpl<Scalar, Options, JointCollection>>
+  {
+    typedef pinocchio::JointModelMimicTpl<Scalar, Options, JointCollection> JointModel;
+    typedef JointModelWithParameters<JointModel> ReturnType;
+
+    static ReturnType run()
+    {
+      typedef pinocchio::JointModelRevoluteTpl<Scalar, Options, 0> JointModelRX;
+      auto jmodelParams_ref = init<JointModelRX>::run();
+      ;
+      JointModelRX jmodel_ref = jmodelParams_ref.jmodel;
+
+      JointModel jmodel(jmodel_ref, 1., 0.);
+      jmodel.setIndexes(1, 0, 0, 0);
+      typename JointModel::ConfigVector_t lb =
+        JointModel::ConfigVector_t::Constant(jmodel.nq(), -1.0);
+      typename JointModel::ConfigVector_t ub =
+        JointModel::ConfigVector_t::Constant(jmodel.nq(), 1.0);
+
+      return {jmodel, lb, ub};
+    }
+  };
+
+  template<typename Scalar, int Options>
+  struct init<pinocchio::JointModelSplineTpl<Scalar, Options>>
+  {
+    typedef pinocchio::JointModelSplineTpl<Scalar, Options> JointModel;
+    typedef JointModelWithParameters<JointModel> ReturnType;
+
+    static ReturnType run()
+    {
+      std::vector<pinocchio::SE3> ctrlFrames;
+      for (int k = 0; k < 5; k++)
+        ctrlFrames.push_back(SE3::Random());
+
+      JointModel jmodel = JointModelSplineBuilder()
+                            .withControlFrameVector(ctrlFrames)
+                            .withDegree(3)
+                            .withOpenUniformKnots(0., 1.)
+                            .build();
+      jmodel.setIndexes(0, 0, 0);
+
+      typename JointModel::ConfigVector_t lb = JointModel::ConfigVector_t::Constant(jmodel.nq(), 0);
+      typename JointModel::ConfigVector_t ub =
+        JointModel::ConfigVector_t::Constant(jmodel.nq(), 1.0);
+      return {jmodel, lb, ub};
+    }
+  };
+} // namespace pinocchio
diff --git a/unittest/visitor.cpp b/unittest/visitor.cpp
index 97cea16568..596f6d3031 100644
--- a/unittest/visitor.cpp
+++ b/unittest/visitor.cpp
@@ -2,6 +2,8 @@
 // Copyright (c) 2015-2019 CNRS INRIA
 //
 
+#include "utils/joints-init.hpp"
+
 #include "pinocchio/multibody.hpp"
 #include "pinocchio/multibody/visitor.hpp"
 
@@ -164,112 +166,6 @@ struct SimpleBinaryVisitor4 : public pinocchio::fusion::JointBinaryVisitorBase<S
 
 BOOST_AUTO_TEST_SUITE(BOOST_TEST_MODULE)
 
-template<typename JointModel_>
-struct init;
-
-template<typename JointModel_>
-struct init
-{
-  static JointModel_ run(const pinocchio::Model & /* model*/)
-  {
-    JointModel_ jmodel;
-    return jmodel;
-  }
-};
-
-template<typename Scalar, int Options>
-struct init<pinocchio::JointModelRevoluteUnalignedTpl<Scalar, Options>>
-{
-  typedef pinocchio::JointModelRevoluteUnalignedTpl<Scalar, Options> JointModel;
-
-  static JointModel run(const pinocchio::Model & /* model*/)
-  {
-    typedef typename JointModel::Vector3 Vector3;
-    JointModel jmodel(Vector3::Random().normalized());
-
-    return jmodel;
-  }
-};
-
-template<typename Scalar, int Options>
-struct init<pinocchio::JointModelRevoluteUnboundedUnalignedTpl<Scalar, Options>>
-{
-  typedef pinocchio::JointModelRevoluteUnboundedUnalignedTpl<Scalar, Options> JointModel;
-
-  static JointModel run(const pinocchio::Model & /* model*/)
-  {
-    typedef typename JointModel::Vector3 Vector3;
-    JointModel jmodel(Vector3::Random().normalized());
-
-    return jmodel;
-  }
-};
-
-template<typename Scalar, int Options>
-struct init<pinocchio::JointModelPrismaticUnalignedTpl<Scalar, Options>>
-{
-  typedef pinocchio::JointModelPrismaticUnalignedTpl<Scalar, Options> JointModel;
-
-  static JointModel run(const pinocchio::Model & /* model*/)
-  {
-    typedef typename JointModel::Vector3 Vector3;
-    JointModel jmodel(Vector3::Random().normalized());
-
-    return jmodel;
-  }
-};
-
-template<typename Scalar, int Options, template<typename, int> class JointCollection>
-struct init<pinocchio::JointModelTpl<Scalar, Options, JointCollection>>
-{
-  typedef pinocchio::JointModelTpl<Scalar, Options, JointCollection> JointModel;
-
-  static JointModel run(const pinocchio::Model & /* model*/)
-  {
-    typedef pinocchio::JointModelRevoluteTpl<Scalar, Options, 0> JointModelRX;
-    JointModel jmodel((JointModelRX()));
-
-    return jmodel;
-  }
-};
-
-template<typename Scalar, int Options, template<typename, int> class JointCollection>
-struct init<pinocchio::JointModelCompositeTpl<Scalar, Options, JointCollection>>
-{
-  typedef pinocchio::JointModelCompositeTpl<Scalar, Options, JointCollection> JointModel;
-
-  static JointModel run(const pinocchio::Model & /* model*/)
-  {
-    typedef pinocchio::JointModelRevoluteTpl<Scalar, Options, 0> JointModelRX;
-    typedef pinocchio::JointModelRevoluteTpl<Scalar, Options, 1> JointModelRY;
-    typedef pinocchio::JointModelRevoluteTpl<Scalar, Options, 2> JointModelRZ;
-
-    JointModel jmodel(JointModelRX(), pinocchio::SE3::Random());
-    jmodel.addJoint(JointModelRY(), pinocchio::SE3::Random());
-    jmodel.addJoint(JointModelRZ(), pinocchio::SE3::Random());
-
-    return jmodel;
-  }
-};
-
-template<typename Scalar, int Options, template<typename, int> class JointCollection>
-struct init<pinocchio::JointModelMimicTpl<Scalar, Options, JointCollection>>
-{
-  typedef pinocchio::JointModelMimicTpl<Scalar, Options, JointCollection> JointModel;
-
-  static JointModel run(const pinocchio::Model & model)
-  {
-    typedef pinocchio::JointModelRevoluteTpl<Scalar, Options, 0> JointModelRX;
-    JointModelRX jmodel_ref = init<JointModelRX>::run(model);
-    jmodel_ref.setIndexes(0, 0, 0, 0);
-
-    JointModel jmodel(jmodel_ref, 1., 0.);
-    jmodel.setIndexes(1, 0, 0, 0);
-
-    return jmodel;
-  }
-};
-
 struct AppendJointToModel
 {
   AppendJointToModel(pinocchio::Model & model)
@@ -280,7 +176,7 @@ struct AppendJointToModel
   template<typename JointModel>
   void operator()(const pinocchio::JointModelBase<JointModel> &) const
   {
-    JointModel jmodel = init<JointModel>::run(model);
+    auto jmodel = pinocchio::init<JointModel>::run().jmodel;
     model.addJoint(model.joints.size() - 1, jmodel, pinocchio::SE3::Random(), jmodel.classname());
   }