diff --git a/pybezier/bezier_curve.py b/pybezier/bezier_curve.py index 355290a..26135b8 100644 --- a/pybezier/bezier_curve.py +++ b/pybezier/bezier_curve.py @@ -5,7 +5,7 @@ class BezierCurve(object): - def __init__(self, points : np.array, initial_time : float = 0, final_time : float = 1): + def __init__(self, points : np.ndarray, initial_time : float = 0, final_time : float = 1): if initial_time >= final_time: raise ValueError("Initial time must be smaller than final time.") self.points = points @@ -15,10 +15,10 @@ def __init__(self, points : np.array, initial_time : float = 0, final_time : flo self.final_time = final_time self.duration = final_time - initial_time - def initial_point(self) -> np.array: + def initial_point(self) -> np.ndarray: return self.points[0] - def final_point(self) -> np.array: + def final_point(self) -> np.ndarray: return self.points[-1] def _berstein(self, time : float | List[float], n : int) -> float: @@ -96,8 +96,16 @@ def __neg__(self) -> Self: def derivative(self) -> Self: points = (self.points[1:] - self.points[:-1]) * (self.degree / self.duration) return BezierCurve(points, self.initial_time, self.final_time) + + def integral(self, initial_condition : np.ndarray | None = None) -> Self: + points = self.points * self.duration / self.degree + points = np.vstack([np.zeros(self.dimension), points]) + points = np.cumsum(points, axis=0) + if initial_condition is not None: + points += initial_condition + return BezierCurve(points, self.initial_time, self.final_time) - def split(self, time : float) -> Tuple[Self, Self]: + def domain_split(self, time : float) -> Tuple[Self, Self]: if time < self.initial_time: raise ValueError("Split time must be larger than initial time.") if time > self.final_time: @@ -132,7 +140,7 @@ def l2_squared(self) -> float: a += b * self.points[i].dot(self.points[j]) return self.duration * a / (2 * self.degree + 1) - def integral_of_convex(self, f : Callable) -> float: + def integral_of_convex_function(self, f : Callable) -> float: c = self.duration / (self.degree + 1) return c * sum(f(point) for point in self.points) diff --git a/pybezier/composite_bezier_curve.py b/pybezier/composite_bezier_curve.py index 81ce391..e1cfebf 100644 --- a/pybezier/composite_bezier_curve.py +++ b/pybezier/composite_bezier_curve.py @@ -1,5 +1,5 @@ import numpy as np -from typing import List, Self +from typing import List, Callable, Self from collections.abc import Iterable from numbers import Number from pybezier.bezier_curve import BezierCurve @@ -27,14 +27,14 @@ def curve_segment(self, time : float) -> int: segment += 1 return segment - def __call__(self, time : float) -> np.array: + def __call__(self, time : float) -> np.ndarray: segment = self.curve_segment(time) return self[segment](time) - def initial_point(self) -> np.array: + def initial_point(self) -> np.ndarray: return self[0].initial_point() - def final_point(self) -> np.array: + def final_point(self) -> np.ndarray: return self[-1].final_point() def __iter__(self) -> Iterable[BezierCurve]: @@ -81,13 +81,13 @@ def __rsub__(self, composite_curve : Self | Number) -> Self: def __neg__(self) -> Self: return 0 - self - def knot_points(self) -> np.array: + def knot_points(self) -> np.ndarray: # assumes that the curve is continuous knots = [curve.initial_point() for curve in self] knots.append(self.final_point()) return np.array(knots) - def durations(self) -> np.array: + def durations(self) -> np.ndarray: return np.array([curve.duration for curve in self]) def concatenate(self, composite_curve : Self) -> Self: @@ -101,10 +101,20 @@ def concatenate(self, composite_curve : Self) -> Self: def derivative(self) -> Self: return CompositeBezierCurve([curve.derivative() for curve in self]) + + def integral(self, initial_condition : np.ndarray | None = None) -> Self: + curves = [] + for curve in self: + curves.append(curve.integral(initial_condition)) + initial_condition = curves[-1].final_point() + return CompositeBezierCurve(curves) def l2_squared(self) -> float: return sum(curve.l2_squared() for curve in self) + def integral_of_convex_function(self, f : Callable) -> float: + return sum(curve.integral_of_convex_function(f) for curve in self) + def plot_components(self, n : int = 51, legend : bool = True, **kwargs): for i, curve in enumerate(self): curve.plot_components(n, legend, **kwargs) diff --git a/tests/test_bezier_curve.py b/tests/test_bezier_curve.py index 103412e..90d4745 100644 --- a/tests/test_bezier_curve.py +++ b/tests/test_bezier_curve.py @@ -109,15 +109,26 @@ def test_neg(self): np.testing.assert_array_almost_equal(neg(time), -self.curve(time)) def test_derivative(self): - der = self.curve.derivative() - time_step = 1e-9 + derivative = self.curve.derivative() + time_step = 1e-6 for time in np.linspace(self.initial_time, self.final_time - time_step): - value = (self.curve(time + time_step) - self.curve(time)) / time_step - np.testing.assert_array_almost_equal(der(time), value) - - def test_split(self): + numerical_derivative = (self.curve(time + time_step) - self.curve(time)) / time_step + np.testing.assert_array_almost_equal(derivative(time), numerical_derivative) + + def test_integral(self): + initial_conditions = [None, np.ones(self.curve.dimension)] + for initial_condition in initial_conditions: + integral = self.curve.integral(initial_condition) + value = 0 if initial_condition is None else initial_condition + np.testing.assert_array_almost_equal(integral(self.initial_time), value) + time_step = 1e-6 + for time in np.linspace(self.initial_time, self.final_time - time_step): + numerical_integral = integral(time) + time_step * self.curve(time) + np.testing.assert_array_almost_equal(integral(time + time_step), numerical_integral) + + def test_domain_split(self): split_time = (self.initial_time + self.final_time) / 2 - curve_1, curve_2 = self.curve.split(split_time) + curve_1, curve_2 = self.curve.domain_split(split_time) for time in self.time_samples: if time < split_time: np.testing.assert_array_almost_equal(self.curve(time), curve_1(time)) @@ -131,14 +142,17 @@ def test_l2_squared(self): integral = np.trapezoid(values, times) self.assertAlmostEqual(self.curve.l2_squared(), integral) - def test_integral_of_convex(self): + def test_integral_of_convex_function(self): f = lambda point: point.dot(point) - self.assertTrue(self.curve.integral_of_convex(f) >= self.curve.l2_squared()) + value = self.curve.l2_squared() + upper_bound = self.curve.integral_of_convex_function(f) + self.assertTrue(value <= upper_bound) # upper bound for curve lenght is equal to distance of control points + diff_points = self.points[:-1] - self.points[1:] + value = sum(np.linalg.norm(diff_points, axis=1)) derivative = self.curve.derivative() - value = sum(np.linalg.norm(y - x) for x, y in zip(self.points[:-1], self.points[1:])) - - self.assertAlmostEqual(derivative.integral_of_convex(np.linalg.norm), value) + upper_bound = derivative.integral_of_convex_function(np.linalg.norm) + self.assertAlmostEqual(value, upper_bound) if __name__ == '__main__': unittest.main() \ No newline at end of file diff --git a/tests/test_composite_bezier_curve.py b/tests/test_composite_bezier_curve.py index 82f8279..cc67a9f 100644 --- a/tests/test_composite_bezier_curve.py +++ b/tests/test_composite_bezier_curve.py @@ -120,11 +120,22 @@ def test_neg(self): np.testing.assert_array_almost_equal(neg(time), -self.composite_curve(time)) def test_derivative(self): - der = self.composite_curve.derivative() - time_step = 1e-9 + derivative = self.composite_curve.derivative() + time_step = 1e-7 for time in np.linspace(self.initial_time, self.final_time - time_step): - value = (self.composite_curve(time + time_step) - self.composite_curve(time)) / time_step - np.testing.assert_array_almost_equal(der(time), value) + numerical_derivative = (self.composite_curve(time + time_step) - self.composite_curve(time)) / time_step + np.testing.assert_array_almost_equal(derivative(time), numerical_derivative) + + def test_integral(self): + initial_conditions = [None, np.ones(self.composite_curve.dimension)] + for initial_condition in initial_conditions: + integral = self.composite_curve.integral(initial_condition) + value = 0 if initial_condition is None else initial_condition + np.testing.assert_array_almost_equal(integral(self.initial_time), value) + time_step = 1e-6 + for time in np.linspace(self.initial_time, self.final_time - time_step): + numerical_integral = integral(time) + time_step * self.composite_curve(time) + np.testing.assert_array_almost_equal(integral(time + time_step), numerical_integral) def test_knot_points(self): for i, point in enumerate(self.composite_curve.knot_points()): @@ -152,5 +163,17 @@ def test_l2_squared(self): integral = np.trapezoid(values, times) self.assertAlmostEqual(self.composite_curve.l2_squared(), integral, places=4) + def test_integral_of_convex_function(self): + f = lambda point: point.dot(point) + value = self.composite_curve.l2_squared() + upper_bound = self.composite_curve.integral_of_convex_function(f) + self.assertTrue(value <= upper_bound) + # upper bound for curve lenght is equal to distance of control points + diff_points = np.vstack([curve.points[:-1] - curve.points[1:] for curve in self.composite_curve]) + value = sum(np.linalg.norm(diff_points, axis=1)) + derivative = self.composite_curve.derivative() + upper_bound = derivative.integral_of_convex_function(np.linalg.norm) + self.assertAlmostEqual(value, upper_bound) + if __name__ == '__main__': unittest.main() \ No newline at end of file