Extend frenet_optimal_trajectory to support more scenarios (AtsushiSakai#1114)

* feat: Add third derivative and curvature rate calculations to CubicSpline classes

* feat: Extend frenet_optimal_trajectory to support more scenarios

* fix: frenet optimal trajectory type check

* fix: cubic spline planner code style check

* fix: frenet optimal trajectory review

* feat: frenet_optimal_trajectory update doc

* fix: frenet optimal trajectory review

* fix: frenet optimal trajectory

* fix: frenet optimal trajectory

Author: Aglargil

PathPlanning/CubicSpline/cubic_spline_planner.py

@@ -76,6 +76,11 @@ def calc_position(self, x):
 
         if `x` is outside the data point's `x` range, return None.
 
+        Parameters
+        ----------
+        x : float
+            x position to calculate y.
+
         Returns
         -------
         y : float
@@ -99,6 +104,11 @@ def calc_first_derivative(self, x):
 
         if x is outside the input x, return None
 
+        Parameters
+        ----------
+        x : float
+            x position to calculate first derivative.
+
         Returns
         -------
         dy : float
@@ -121,6 +131,11 @@ def calc_second_derivative(self, x):
 
         if x is outside the input x, return None
 
+        Parameters
+        ----------
+        x : float
+            x position to calculate second derivative.
+
         Returns
         -------
         ddy : float
@@ -137,6 +152,31 @@ def calc_second_derivative(self, x):
         ddy = 2.0 * self.c[i] + 6.0 * self.d[i] * dx
         return ddy
 
+    def calc_third_derivative(self, x):
+        """
+        Calc third derivative at given x.
+
+        if x is outside the input x, return None
+
+        Parameters
+        ----------
+        x : float
+            x position to calculate third derivative.
+
+        Returns
+        -------
+        dddy : float
+            third derivative for given x.
+        """
+        if x < self.x[0]:
+            return None
+        elif x > self.x[-1]:
+            return None
+
+        i = self.__search_index(x)
+        dddy = 6.0 * self.d[i]
+        return dddy
+
     def __search_index(self, x):
         """
         search data segment index
@@ -287,6 +327,33 @@ def calc_curvature(self, s):
         k = (ddy * dx - ddx * dy) / ((dx ** 2 + dy ** 2)**(3 / 2))
         return k
 
+    def calc_curvature_rate(self, s):
+        """
+        calc curvature rate
+
+        Parameters
+        ----------
+        s : float
+            distance from the start point. if `s` is outside the data point's
+            range, return None.
+
+        Returns
+        -------
+        k : float
+            curvature rate for given s.
+        """
+        dx = self.sx.calc_first_derivative(s)
+        dy = self.sy.calc_first_derivative(s)
+        ddx = self.sx.calc_second_derivative(s)
+        ddy = self.sy.calc_second_derivative(s)
+        dddx = self.sx.calc_third_derivative(s)
+        dddy = self.sy.calc_third_derivative(s)
+        a = dx * ddy - dy * ddx
+        b = dx * dddy - dy * dddx
+        c = dx * ddx + dy * ddy
+        d = dx * dx + dy * dy
+        return (b * d - 3.0 * a * c) / (d * d * d)
+
     def calc_yaw(self, s):
         """
         calc yaw

PathPlanning/FrenetOptimalTrajectory/cartesian_frenet_converter.py

@@ -0,0 +1,144 @@
+"""
+
+A converter between Cartesian and Frenet coordinate systems
+
+author: Wang Zheng (@Aglargil)
+
+Ref:
+
+- [Optimal Trajectory Generation for Dynamic Street Scenarios in a Frenet Frame]
+(https://www.researchgate.net/profile/Moritz_Werling/publication/224156269_Optimal_Trajectory_Generation_for_Dynamic_Street_Scenarios_in_a_Frenet_Frame/links/54f749df0cf210398e9277af.pdf)
+
+"""
+
+import math
+
+
+class CartesianFrenetConverter:
+    """
+    A class for converting states between Cartesian and Frenet coordinate systems
+    """
+
+    @ staticmethod
+    def cartesian_to_frenet(rs, rx, ry, rtheta, rkappa, rdkappa, x, y, v, a, theta, kappa):
+        """
+        Convert state from Cartesian coordinate to Frenet coordinate
+
+        Parameters
+        ----------
+            rs: reference line s-coordinate
+            rx, ry: reference point coordinates
+            rtheta: reference point heading
+            rkappa: reference point curvature
+            rdkappa: reference point curvature rate
+            x, y: current position
+            v: velocity
+            a: acceleration
+            theta: heading angle
+            kappa: curvature
+
+        Returns
+        -------
+            s_condition: [s(t), s'(t), s''(t)]
+            d_condition: [d(s), d'(s), d''(s)]
+        """
+        dx = x - rx
+        dy = y - ry
+
+        cos_theta_r = math.cos(rtheta)
+        sin_theta_r = math.sin(rtheta)
+
+        cross_rd_nd = cos_theta_r * dy - sin_theta_r * dx
+        d = math.copysign(math.hypot(dx, dy), cross_rd_nd)
+
+        delta_theta = theta - rtheta
+        tan_delta_theta = math.tan(delta_theta)
+        cos_delta_theta = math.cos(delta_theta)
+
+        one_minus_kappa_r_d = 1 - rkappa * d
+        d_dot = one_minus_kappa_r_d * tan_delta_theta
+
+        kappa_r_d_prime = rdkappa * d + rkappa * d_dot
+
+        d_ddot = (-kappa_r_d_prime * tan_delta_theta +
+                  one_minus_kappa_r_d / (cos_delta_theta * cos_delta_theta) *
+                  (kappa * one_minus_kappa_r_d / cos_delta_theta - rkappa))
+
+        s = rs
+        s_dot = v * cos_delta_theta / one_minus_kappa_r_d
+
+        delta_theta_prime = one_minus_kappa_r_d / cos_delta_theta * kappa - rkappa
+        s_ddot = (a * cos_delta_theta -
+                  s_dot * s_dot *
+                  (d_dot * delta_theta_prime - kappa_r_d_prime)) / one_minus_kappa_r_d
+
+        return [s, s_dot, s_ddot], [d, d_dot, d_ddot]
+
+    @ staticmethod
+    def frenet_to_cartesian(rs, rx, ry, rtheta, rkappa, rdkappa, s_condition, d_condition):
+        """
+        Convert state from Frenet coordinate to Cartesian coordinate
+
+        Parameters
+        ----------
+            rs: reference line s-coordinate
+            rx, ry: reference point coordinates
+            rtheta: reference point heading
+            rkappa: reference point curvature
+            rdkappa: reference point curvature rate
+            s_condition: [s(t), s'(t), s''(t)]
+            d_condition: [d(s), d'(s), d''(s)]
+
+        Returns
+        -------
+            x, y: position
+            theta: heading angle
+            kappa: curvature
+            v: velocity
+            a: acceleration
+        """
+        if abs(rs - s_condition[0]) >= 1.0e-6:
+            raise ValueError(
+                "The reference point s and s_condition[0] don't match")
+
+        cos_theta_r = math.cos(rtheta)
+        sin_theta_r = math.sin(rtheta)
+
+        x = rx - sin_theta_r * d_condition[0]
+        y = ry + cos_theta_r * d_condition[0]
+
+        one_minus_kappa_r_d = 1 - rkappa * d_condition[0]
+
+        tan_delta_theta = d_condition[1] / one_minus_kappa_r_d
+        delta_theta = math.atan2(d_condition[1], one_minus_kappa_r_d)
+        cos_delta_theta = math.cos(delta_theta)
+
+        theta = CartesianFrenetConverter.normalize_angle(delta_theta + rtheta)
+
+        kappa_r_d_prime = rdkappa * d_condition[0] + rkappa * d_condition[1]
+
+        kappa = (((d_condition[2] + kappa_r_d_prime * tan_delta_theta) *
+                  cos_delta_theta * cos_delta_theta) / one_minus_kappa_r_d + rkappa) * \
+            cos_delta_theta / one_minus_kappa_r_d
+
+        d_dot = d_condition[1] * s_condition[1]
+        v = math.sqrt(one_minus_kappa_r_d * one_minus_kappa_r_d *
+                      s_condition[1] * s_condition[1] + d_dot * d_dot)
+
+        delta_theta_prime = one_minus_kappa_r_d / cos_delta_theta * kappa - rkappa
+
+        a = (s_condition[2] * one_minus_kappa_r_d / cos_delta_theta +
+             s_condition[1] * s_condition[1] / cos_delta_theta *
+             (d_condition[1] * delta_theta_prime - kappa_r_d_prime))
+
+        return x, y, theta, kappa, v, a
+
+    @ staticmethod
+    def normalize_angle(angle):
+        """
+        Normalize angle to [-pi, pi]
+        """
+        a = math.fmod(angle + math.pi, 2.0 * math.pi)
+        if a < 0.0:
+            a += 2.0 * math.pi
+        return a - math.pi