simetri.geometry.hobby

  1import numpy as np
  2import cmath
  3
  4from ..graphics.shape import Shape
  5from ..geometry.bezier import bezier_points
  6
  7"""
  8Implementation of John Hobby's Bezier curve algorithm in Python.
  9The algorithm is very simple and efficient. Details are on page 112, 113 in Knuth's METAFONT: The Program.
 10"""
 11# Taken from https://github.com/ltrujello/Hobby_Curve_Algorithm 2/7/2025
 12
 13
 14class HobbyPoint(complex):
 15    """A class for associating numerical quantities from Hobby's algorithm with points that appear on a Hobby curve.
 16
 17    We subclass `complex` to perform complex arithmetic with points on the Hobby curve, as required in the algorithm.
 18
 19    Attributes:
 20        x (float): The x-coordinate of the point.
 21        y (float): The y-coordinate of the point.
 22        alpha (float): The reciprocal of tension for incoming segment.
 23        beta (float): The reciprocal of tension for outgoing segment.
 24        d_val (float): Distance between this point and next.
 25        theta (float): Angle of polygonal line from this point to next.
 26        phi (float): Offset angle.
 27        psi (float): Another offset angle.
 28    """
 29
 30    def __new__(cls, x: float, y: float, tension: float) -> "HobbyPoint":
 31        """Create a new instance of HobbyPoint.
 32
 33        Args:
 34            x: The x-coordinate of the point.
 35            y: The y-coordinate of the point.
 36            tension: The tension value for the curve.
 37
 38        Returns:
 39            A new HobbyPoint instance with complex value (x + y*j).
 40        """
 41        return super().__new__(cls, x, y)
 42
 43    def __init__(self, x: float, y: float, tension: float) -> None:
 44        """Initialize a HobbyPoint with coordinates and tension.
 45
 46        Args:
 47            x: The x-coordinate of the point.
 48            y: The y-coordinate of the point.
 49            tension: The tension value for the curve.
 50        """
 51        self.x = x
 52        self.y = y
 53        # In what follows, we use Knuth's notation in our variable names.
 54        self.alpha = 1 / tension
 55        self.beta = 1 / tension
 56        self.d_val = 0  # Distance between this point and next.
 57        self.theta = 0  # Angle of polygonal line from this point to next.
 58        self.phi = 0  # Offset angle.
 59        self.psi = 0  # Another offset angle.
 60
 61    def debug_info(self) -> str:
 62        """Return a string with the point's information.
 63
 64        Returns:
 65            A string containing the point's coordinates and all of its computational values.
 66        """
 67        return (
 68            f"{(self.x, self.y)} "
 69            f"alpha={self.alpha}, "
 70            f"beta={self.beta}, "
 71            f"theta={self.theta}, "
 72            f"psi={self.psi}, "
 73            f"phi={self.phi}, "
 74            f"d_val={self.d_val}"
 75        )
 76
 77    def __repr__(self) -> str:
 78        """Return a string representation of the point.
 79
 80        Returns:
 81            A string representation of the point's coordinates.
 82        """
 83        return f"{(self.x, self.y)}"
 84
 85
 86class HobbyCurve:
 87    """A class for calculating the control points required to draw a Hobby curve.
 88
 89    Attributes:
 90        points (list[HobbyPoint]): The list of points defining the curve.
 91        ctrl_pts (list[tuple]): The calculated control points.
 92        is_cyclic (bool): Whether the curve is closed.
 93        begin_curl (float): Curl value for the beginning of the curve.
 94        end_curl (float): Curl value for the end of the curve.
 95        n_points (int): Number of points in the curve.
 96        debug_mode (bool): Whether to print debug information.
 97    """
 98
 99    def __init__(
100        self,
101        points: list[tuple],
102        tension: float = 1,
103        cyclic: bool = False,
104        begin_curl: float = 1,
105        end_curl: float = 1,
106        debug: bool = False,
107    ) -> None:
108        """Initialize a HobbyCurve with the given parameters.
109
110        Args:
111            points: List of (x, y) tuples representing the curve's points.
112            tension: Tension parameter controlling the "tightness" of the curve.
113            cyclic: Whether the curve should be closed.
114            begin_curl: Curl value for the beginning of the curve.
115            end_curl: Curl value for the end of the curve.
116            debug: Whether to print debug information.
117
118        Raises:
119            ValueError: If fewer than 2 points are provided.
120        """
121        if len(points) < 2:
122            raise ValueError("Algorithm needs more than 2 points")
123        self.points = [HobbyPoint(*point, tension) for point in points]
124        self.ctrl_pts = []
125        self.is_cyclic = cyclic
126        self.begin_curl = begin_curl
127        self.end_curl = end_curl
128        self.n_points = len(points)
129        self.debug_mode = debug
130
131    def get_ctrl_points(self) -> list[tuple]:
132        """Calculate and return all of the control points of the Hobby curve.
133
134        Executes the Hobby algorithm by calculating distance values, psi values,
135        theta values, and phi values, then uses these to compute control points.
136
137        Returns:
138            A list of (x, y) tuples representing the Bezier control points.
139        """
140        self.calculate_d_vals()
141        self.calculate_psi_vals()
142        self.calculate_theta_vals()
143        self.calculate_phi_vals()
144        self.show_debug_msg()
145        self.ctrl_pts = self.calculate_ctrl_pts()
146        return self.ctrl_pts
147
148    def calculate_d_vals(self) -> None:
149        """Calculate the pairwise distances between consecutive points in the curve."""
150        # Skip last point if path is non-cyclic
151        point_inds = (
152            range(self.n_points) if self.is_cyclic else range(self.n_points - 1)
153        )
154        for i in point_inds:
155            z_i = self.points[i % self.n_points]
156            z_j = self.points[(i + 1) % self.n_points]
157            z_i.d_val = abs(z_i - z_j)
158
159    def calculate_psi_vals(self) -> None:
160        """Calculate the psi values by finding the angle of the polygonal turns.
161
162        Raises:
163            ZeroDivisionError: If consecutive points have the same coordinates.
164        """
165        # Skip first and last point if path is non-cyclic
166        point_inds = (
167            range(self.n_points) if self.is_cyclic else range(1, self.n_points - 1)
168        )
169        for i in point_inds:
170            z_h = self.points[i - 1]
171            z_i = self.points[i]
172            z_j = self.points[(i + 1) % self.n_points]
173            try:
174                polygonal_turn = (z_j - z_i) / (z_i - z_h)
175                # print(z_j - z_i, z_i - z_h)
176            except ZeroDivisionError:
177                raise ZeroDivisionError(
178                    f"Consecutive points {(z_h.x, z_h.y)} and {(z_i.x, z_i.y)} cause zero division."
179                )
180            z_i.psi = np.arctan2(polygonal_turn.imag, polygonal_turn.real)
181
182    def calculate_theta_vals(self) -> None:
183        """Calculate the theta values by solving a linear system of equations.
184
185        This is the core of Hobby's algorithm, creating and solving a system of equations
186        to find the optimal angles for smooth splines.
187        """
188        A = np.zeros(
189            self.n_points
190        )  # Inappropriate names, but they mirror Knuth's notation.
191        B = np.zeros(self.n_points)
192        C = np.zeros(self.n_points)
193        D = np.zeros(self.n_points)
194        R = np.zeros(self.n_points)
195
196        # Calculate the entries of the five vectors.
197        # Skip first and last point if path is non-cyclic.
198        point_ind = (
199            range(self.n_points) if self.is_cyclic else range(1, self.n_points - 1)
200        )
201        for i in point_ind:
202            z_h = self.points[i - 1]
203            z_i = self.points[i]
204            z_j = self.points[(i + 1) % self.n_points]
205
206            A[i] = z_h.alpha / (z_i.beta**2 * z_h.d_val)
207            B[i] = (3 - z_h.alpha) / (z_i.beta**2 * z_h.d_val)
208            C[i] = (3 - z_j.beta) / (z_i.alpha**2 * z_i.d_val)
209            D[i] = z_j.beta / (z_i.alpha**2 * z_i.d_val)
210            R[i] = -B[i] * z_i.psi - D[i] * z_j.psi
211
212        # Set up matrix M such that the soln. Mx = R are the theta values.
213        M = np.zeros((self.n_points, self.n_points))
214        for i in range(self.n_points):
215            # Fill i-th row of M
216            M[i][i - 1] = A[i]
217            M[i][i] = B[i] + C[i]
218            M[i][(i + 1) % self.n_points] = D[i]
219
220        # Special formulas for first and last rows of M with non-cyclic paths.
221        if not self.is_cyclic:
222            # First row of M
223            alpha_0 = self.points[0].alpha
224            beta_1 = self.points[1].beta
225            xi_0 = (alpha_0**2 * self.begin_curl) / beta_1**2
226            M[0][0] = alpha_0 * xi_0 + 3 - beta_1
227            M[0][1] = (3 - alpha_0) * xi_0 + beta_1
228            R[0] = -((3 - alpha_0) * xi_0 + beta_1) * self.points[1].psi
229            # Last row of M
230            alpha_n_1 = self.points[-2].alpha
231            beta_n = self.points[-1].beta
232            xi_n = (beta_n**2 * self.end_curl) / alpha_n_1**2
233            M[-1][-2] = (3 - beta_n) * xi_n + alpha_n_1
234            M[-1][-1] = beta_n * xi_n + 3 - alpha_n_1
235            R[-1] = 0
236
237        # Solve for theta values.
238        thetas = np.linalg.solve(M, R)
239        for i, point in enumerate(self.points):
240            point.theta = thetas[i]
241
242    def calculate_phi_vals(self) -> None:
243        """Calculate the phi values using the relationship theta + phi + psi = 0."""
244        for point in self.points:
245            point.phi = -(point.psi + point.theta)
246
247    def calculate_ctrl_pts(self) -> list[tuple]:
248        """Calculate the Bezier control points between consecutive points.
249
250        Returns:
251            A list of (x, y) tuples representing the control points, with two control points
252            for each curve segment between consecutive points.
253        """
254        ctrl_pts = []
255        # Skip last point if path is non-cyclic
256        point_inds = (
257            range(self.n_points) if self.is_cyclic else range(self.n_points - 1)
258        )
259        for i in point_inds:
260            z_i = self.points[i]
261            z_j = self.points[(i + 1) % self.n_points]
262            rho_coefficient = z_i.alpha * velocity(z_i.theta, z_j.phi)
263            sigma_coefficient = z_j.beta * velocity(z_j.phi, z_i.theta)
264            ctrl_pt_a = z_i + (1 / 3) * rho_coefficient * cmath.exp(
265                complex(0, z_i.theta)
266            ) * (z_j - z_i)
267            ctrl_pt_b = z_j - (1 / 3) * sigma_coefficient * cmath.exp(
268                complex(0, -z_j.phi)
269            ) * (z_j - z_i)
270            ctrl_pts.append((ctrl_pt_a.real, ctrl_pt_a.imag))
271            ctrl_pts.append((ctrl_pt_b.real, ctrl_pt_b.imag))
272        return ctrl_pts
273
274    def show_debug_msg(self) -> None:
275        """Display debug information for each point if debug mode is enabled."""
276        if self.debug_mode:
277            for point in self.points:
278                print(point.debug_info())
279
280    def __repr__(self) -> str:
281        """Return a string representation of the curve.
282
283        Returns:
284            A string representation of the curve's points in Cartesian coordinates.
285        """
286        cartesian_points = [(point.real, point.imag) for point in self.points]
287        return repr(cartesian_points)
288
289
290def hobby_ctrl_points(
291    points: list[tuple],
292    tension: float = 1,
293    cyclic: bool = False,
294    begin_curl: float = 1,
295    end_curl: float = 1,
296    debug: bool = False,
297) -> list[tuple]:
298    """Calculate cubic Bezier control points using John Hobby's algorithm.
299
300    Args:
301        points: List of (x, y) tuples representing the curve's points.
302        tension: Controls the "tightness" of the curve (lower is tighter).
303        cyclic: Whether the curve should be closed.
304        begin_curl: Curl value for the beginning of the curve.
305        end_curl: Curl value for the end of the curve.
306        debug: Whether to print debug information.
307
308    Returns:
309        A list of (x, y) tuples representing the Bezier control points.
310    """
311    curve = HobbyCurve(
312        points,
313        tension=tension,
314        cyclic=cyclic,
315        begin_curl=begin_curl,
316        end_curl=end_curl,
317        debug=debug,
318    )
319    ctrl_points = curve.get_ctrl_points()
320
321    # Calculate whitespace padding for pretty print.
322    max_pad = 0
323    for ctrl_point in ctrl_points:
324        x, y = ctrl_point[:2]
325        # Calculate number of digits in x, y before decimal, and take the
326        # max for nice padding.
327        padding = max(
328            1 if abs(x) <= 0.1 else int(np.ceil(np.log10(abs(x)))) + 1,
329            1 if abs(y) <= 0.1 else int(np.ceil(np.log10(abs(y)))) + 1,
330        )
331        if max_pad < padding:
332            max_pad = padding
333
334    if debug:
335        # Pretty print control points.
336        precision = 10
337        space = precision + max_pad + 1  # +1 for negative sign
338        i = 0
339        while i < len(ctrl_points) - 1:
340            x_1, y_1 = ctrl_points[i]
341            x_2, y_2 = ctrl_points[i + 1]
342            print(
343                f"({x_1:<{space}.{precision}f}, {y_1:<{space}.{precision}f}) "
344                f"and "
345                f"({x_2:<{space}.{precision}f}, {y_2:<{space}.{precision}f})"
346            )
347            i += 2
348
349    return ctrl_points
350
351
352def velocity(theta: float, phi: float) -> float:
353    """Calculate the "velocity" function used in Metafont's curve algorithm.
354
355    This function implements the specific velocity formula from Knuth's Metafont.
356
357    Args:
358        theta: The theta angle value.
359        phi: The phi angle value.
360
361    Returns:
362        The computed velocity value used in control point calculations.
363    """
364    numerator = 2 + np.sqrt(2) * (np.sin(theta) - (1 / 16) * np.sin(phi)) * (
365        np.sin(phi) - (1 / 16) * np.sin(theta)
366    ) * (np.cos(theta) - np.cos(phi))
367    denominator = (
368        1
369        + (1 / 2) * (np.sqrt(5) - 1) * np.cos(theta)
370        + (1 / 2) * (3 - np.sqrt(5)) * np.cos(phi)
371    )
372    return numerator / denominator
373
374def hobby_shape(points, cyclic=False, tension=1, begin_curl=1, end_curl=1, n_points=None):
375    """Create a Shape object from points using John Hobby's algorithm.
376
377    This function calculates cubic Bezier control points using Hobby's algorithm,
378    then creates a Shape object by generating points along the resulting Bezier curves.
379
380    Args:
381        points: List of (x, y) tuples representing the curve's points.
382        cyclic: Whether the curve should be closed.
383        tension: Controls the "tightness" of the curve (lower is tighter).
384        begin_curl: Curl value for the beginning of the curve.
385        end_curl: Curl value for the end of the curve.
386        debug: Whether to print debug information.
387
388    Returns:
389        A Shape object containing points along the smooth Hobby curve.
390    """
391    controls = hobby_ctrl_points(points, tension=tension, cyclic=cyclic,
392                                            begin_curl=begin_curl, end_curl=end_curl)
393    n = len(points)
394    res = []
395    if cyclic:
396        for i in range(n):
397            ind = i * 2
398            p0 = points[i]
399            p1 = controls[ind]
400            p2 = controls[ind + 1]
401            p3 = points[(i + 1)%n]
402            bez_pnts = bezier_points(p0, p1, p2, p3, 10)
403            res.extend(bez_pnts)
404    else:
405        for i in range(n - 1):
406            ind = i * 2
407            p0 = points[i]
408            p1 = controls[ind]
409            p2 = controls[ind + 1]
410            p3 = points[(i + 1)]
411            bez_pnts = bezier_points(p0, p1, p2, p3, 20)
412            res.extend(bez_pnts)
413    return Shape(res)