Approximating a natural function spiral with a finite series of straight lines and constant curvature circle segments

Question

I am in need of bending a wire to a spiral shape defined by

$$ r(\phi)=a + b\phi + c\phi^2+d\phi^3 $$

My bending machine is capable of bending circle segments with a constant radius $R$ and angle $\alpha=n\cdot \alpha_0$ and feeding straight wire with the length $l=n\cdot l_0$ where $n$ is an integer number and $\alpha_0,l_0$ are machine parameters.

Another constraint is that every straight line needs to be tangent to the next circle segment and that every circle segment needs to be followed by a straight segment, and vice versa.

Below is an image of a circle approximated by this method, but I am not sure how it would relate to a spiral path.

Does anybody have an idea how I could generate a series of circle segments and straight lines that approximates my curve with above constraints (i.e. $n\in \mathbb{N}$) such that the area integral between the original curve and the approximate path is minimized?

Ideally I would want a matlab or python script that outputs LRA (or YBC) coordinate data which is what a typical bending machine takes. ( LRA stands for length, rotation, angle for every segment where rotation would be required for a 3D shape but is not used in this case.)

Answer 1

I ended up writing the following script. Unfortunately, I did not get the chance to test it as my priorities shifted, and I did not figure out the springback and clearance compensations yet, but the bending commands look realistic for the spiral shape I generated for.

You can find the code and the files for my bender in my github repository https://github.com/michox/Wire-Bender/blob/main/generate_g_code.py

import sympy as smp
import numpy as np
import matplotlib.pyplot as plt
import scipy as scp
from scipy import optimize
from math import ceil, floor
# Variables
start_angle = 0
end_angle = 2 * np.pi
R_tool = 2  # radius of the bending tool
d_wire = 1  # diameter of the wire
R = 2  # R_tool + d_wire / 2
alpha_0 = (
    2 * np.pi / 200
)  # step angle in radiant. this is the value for most standard steppers with 1.8 deg
l_0 = (
    5.5 * alpha_0
)  # radius of the feeding tool multiplied with the minimum step angle of the stepper
n_min = ceil(l_0/alpha_0/R)

t, tau = smp.symbols("t, tau", nonnegative=True)
s = smp.symbols("s", nonnegative=True)
circle_start, circle_end, line_end = smp.symbols("circle_start,circle_end, line_end", nonnegative=True)

# Your functions
x_expr = R * t * smp.cos(t)
y_expr = R * t * smp.sin(t)
# Math

dtx_expr = x_expr.diff(t)
dty_expr = y_expr.diff(t)
ddtx_expr = dtx_expr.diff(t)
ddty_expr = dty_expr.diff(t)

x = smp.lambdify(t, x_expr, "numpy")
y = smp.lambdify(t, y_expr, "numpy")
dtx = smp.lambdify(t, dtx_expr, "numpy")
dty = smp.lambdify(t, dty_expr, "numpy")

curvature = (dtx_expr * ddty_expr - ddtx_expr * dty_expr) / (
    (dtx_expr**2 + dty_expr**2) ** (3 / 2)
)
curvature = curvature.simplify()

arc_length = smp.integrate(
    smp.sqrt(dtx_expr**2 + dty_expr**2), (tau, start_angle, t)
).simplify()
solution = smp.solve(smp.Eq(arc_length, s), t)
curvature_of_s_expr = curvature.subs(t, solution[0])

length: float = float(smp.N(
    arc_length.subs(t, end_angle)
 - arc_length.subs(t, start_angle)
))
n_max = ceil(length / l_0)

# calculate a general expression for the circle deviation as a function of the circle end point circle_end
circle_deviation = smp.integrate(curvature_of_s_expr - 1 / R, (s, circle_start, circle_end))
# turn the general expression into a numeric function
circle_deviation_f = smp.lambdify([circle_start, circle_end], circle_deviation, "scipy")
# calculate a general expression for the circle deviation as a function of the circle end point circle_end and the line end point line_end
line_deviation = smp.integrate(curvature_of_s_expr, (s, circle_end, line_end))
# turn the general expression into a numeric function
line_deviation_f = smp.lambdify([circle_end, line_end], line_deviation, "scipy")

# carry over any deviation of previous segments
carry_over_deviation = 0



def total_deviation(circle_start_f,circle_end_f, line_end_f):
    return (
        circle_deviation_f(circle_start_f, circle_end_f)
        + line_deviation_f(circle_end_f, line_end_f)
        + carry_over_deviation
    )


def has_root_constraint(circle_start_f, circle_end_f, line_end_f):
    # deviation for smallest possible line_end and deviation for largest possible line_end have different sign assuming there is only 1 root (decreasing curvature)
    # if the two results have a different sign, the function has a root. The constraint expects a value >0 if the constraint is fulfilled and <0 if not so we add a negative sign
    return - total_deviation(circle_start_f, circle_end_f, line_end_f) * total_deviation(circle_start_f, circle_end_f, length)
    


# start the calculation for circle_start=0
circle_start_f = 0


circle_end_p = np.linspace(0, length, n_max)
line_end_p = np.linspace(0, length, n_max)
def plot_deviation(circle_start_f):
    # calculate the total deviation for circle_start_f for all circle_end and line_end in circle_end_p and line_end_p
    total_deviation_p = np.vectorize(total_deviation)(circle_start_f, circle_end_p[:, np.newaxis], line_end_p)
    # plot the deviation as a surface over circle_end_p and line_end_p
    plt.figure()
    plt.title("deviation")
    plt.xlabel("circle_end")
    plt.ylabel("line_end")
    # surface plot
    plt.contourf(circle_end_p, line_end_p, total_deviation_p, 100)
    plt.colorbar()
    plt.show()
# Generate Segments
segments = [0] ## first segment start point is 0
circle_start_f=0
while circle_start_f < length:
    print("calculating first end point circle_end")
    
    circle_end_min_f = circle_start_f + n_min*R*alpha_0

    # find the circle end point for which the deviation is zero at the given circle_start   
    try: 
        circle_end_f = scp.optimize.root_scalar(
            lambda circle_end: circle_deviation_f(circle_start_f, circle_end)+carry_over_deviation,
            x0=circle_end_min_f,
            bracket=[circle_end_min_f, length-l_0], # boundaries for the arc end position
        ).root
    except ValueError:
        circle_end_f = circle_end_min_f
    # # round the end point to the nearest value that will be a multiple of the step angle R*alpha_0 to serve as starting point for the optimization
    circle_end_f = ceil(circle_end_f / alpha_0 / R) * alpha_0 * R

    # try different circle_end and find the shortest line_end where the total deviation will be 0
    def minimize_line_end(circle_end_f):
        # find the line end point line_end for which the total deviation of the arc and the line segment is zero
        try:
            return scp.optimize.root_scalar(
                lambda line_end_f: total_deviation(circle_start_f, circle_end_f, line_end_f),
                x0=circle_end_f + l_0, # initial guess
                bracket=[circle_end_f, length],
            ).root # returns the line end point line_end_f where the total deviation is minimized
        except ValueError:
            print("no root found")
            return length

    # find the circle end for which the total deviation is minimized
    circle_end_opt = scp.optimize.minimize(
        minimize_line_end, 
        circle_end_f, # initial guess
        constraints=({"type": "ineq", "fun": lambda circle_end : has_root_constraint(circle_start_f, circle_end,circle_end + l_0)}),
        bounds=scp.optimize.Bounds(circle_end_min_f, length, True),
    ).x[0]

    # find closest mechanically possible position circle_end for the end of the arc
    circle_end_opt = ceil(circle_end_opt / alpha_0 / R) * alpha_0 * R
    # calculate the optimal end point of the line segment line_end for the new circle_end
    line_end_opt = minimize_line_end(circle_end_opt)
    # round again
    line_end_opt = ceil(line_end_opt / l_0) * l_0

    # calculate any left over deviation resulting from mechanical limitations to consider in the next segments
    carry_over_deviation = total_deviation(circle_start_f, circle_end_opt,line_end_opt)
    # set the new start point
    circle_start_f = line_end_opt
    # add the new segments to the list of segments
    print("found new segments circle_end={}, line_end={}".format(circle_end_opt, line_end_opt))
    print("carry over deviation={}".format(carry_over_deviation))
    segments.extend([circle_end_opt, line_end_opt])
# Plotting

##### gcode generation ######   
## stepper is configured to do 200 steps per mm so that 1 mm corresponds to 1 rotation instead.
l_compensation = 1.0 # compensation factor for the length of the line segments
bend_init = 0.0 # bend angle to add to the beginning of every end to compensate for clearance in the bend tool
bend_factor = 1.0 # factor to multiply the bend angle with to compensate for spring back

def bend(n): 
    # 
    return "G91\nG0 X{:.3f}\nG0 X-{:.3f}\nG90\n".format(n*0.01,n*0.01)

def feed(n):
    return "G0 Y{:.3f}\n".format(n*0.01)

# open a new text file
file = open("bending-file.gcode", "w")
for i in range(1,len(segments)):
    if i%2 != 0:
        file.write(feed(int((segments[i]-segments[i-1])/l_0)))
        file.write(bend(int((segments[i]-segments[i-1])/R/alpha_0)))
    else:
        file.write(feed(int((segments[i]-segments[i-1])/l_0)))

Answer 2

Abstractly, we're given a regular parametrized plane curve—possibly a polar graph $r = f(\phi)$ for $a \leq \phi \leq b$, possibly a Cartesian path $(x(t), y(t))$ for $a \leq t \leq b$—and the goal is to approximate the image by a finite sequence of

1. Line segments of fixed length $\ell_{0}$ [sic],
2. Circle arcs of fixed radius $R$ and angular size $\alpha_{0}$ [sic].

Qualitatively, we have three types of curve atom (or two if an arc can bend only leftward), and we want to build a molecular chain whose shape approximates the target path.

An approximation scheme may be viewed as constructing a finite sequence of segments-and-arcs. Because successive atoms are tangent at their common endpoints (as stipulated in an edit), the position and orientation of any particular atom locks in most of the freedom: There is a two-parameter family of initial positions and a one-parameter family of initial angles; after that there are only finitely many choices in total.

For definiteness, I'll assume arcs can bend in two directions (left and right). One geometrically natural idea is to approximate the signed curvature of the target path. For a regular (non-vanishing velocity) Cartesian parametrization $(x(t), y(t))$ with $a \leq t \leq b$, the signed curvature is $$ k(t) = \frac{x'(t)y''(t) - x''(t)y'(t)}{(x'(t)^{2} + y'(t)^{2})^{3/2}}. $$ To use this, we need to normalize in terms of arc length. Define $$ s = \sigma(t) = \int_{a}^{t} \sqrt{x'(\tau)^{2} + y'(\tau)^{2}}\, d\tau, $$ let $t = \sigma^{-1}(s)$ denote the inverse functional relationship, and finally define the signed curvature as a function of arc length $$ \kappa(s) = k(t) = k(\sigma^{-1}(s)). $$ The function $k$ can generally be found analytically; the arc length parameter $s$ may involve a "non-elementary" integral, and in any case the functional inversion needed to write the path parameter $t = \sigma^{-1}(s)$ in terms of arc length generally can be found only numerically.

Those difficulties aside, our problem is fairly simple. Using the machine parameters in the question, we want to approximate the curvature function $\kappa$ of the target path with a piecewise-constant function whose values are either $0$ on intervals of length $\ell = n\ell_{0}$, or else $\pm 1/R$ on intervals of length $\alpha R = n\alpha_{0}R$. (If the machine can only bend wire to the left, discard the value $-1/R$ coming from bending to the right.) That in turn is a matter of iteratively inspecting the values of $\kappa$ along the "next prospective interval," choosing one of the three (or two) available values, and (not) bending the wire accordingly.

To carry out this scheme in practice (for "tame" target shapes in a practical sense), it suffices to

• Calculate the curvature function $k(t)$ for $a \leq t \leq b$;
• Calculate the arc length function $s = \sigma(t)$ for $a \leq t \leq b$, or numerically approximate at a resolution suitably finer that the machine's smallest parameters ($\ell_{0}$ and $\alpha_{0}R$);
• Numerically approximate the inverse function $t = \sigma^{-1}(s)$, and tabulate the resulting values of the arc length-normalized curvature $\kappa(s)$.

That leaves rather a lot of numerical and logical coding, but is mathematically well-founded and looks feasible to implement.