From 7aab27266fbc1a0a0831151bd87a0b5ca4d19f55 Mon Sep 17 00:00:00 2001 From: Adrian Mariano Date: Thu, 18 Jun 2020 17:50:25 -0400 Subject: [PATCH] Added polynomial root finding and complex math to math.scad. Changed path_to_bezier to new algorithm (which uses polynomial roots). And updated path_smooth to use the new code. Also removed extra (?) include and include from rounding.scad and removed the common $fn=36 that was forced on all the examples. (May break something...we'll look at the examples and see.) I added $fn=36 to some examples. --- beziers.scad | 97 ++++++++++++++++++++++++++++------------- math.scad | 116 ++++++++++++++++++++++++++++++++++++++++++++++++++ rounding.scad | 100 ++++++++++++++++++++++++++----------------- 3 files changed, 245 insertions(+), 68 deletions(-) diff --git a/beziers.scad b/beziers.scad index 68c59fc..0407e49 100644 --- a/beziers.scad +++ b/beziers.scad @@ -515,40 +515,79 @@ function bezier_polyline(bezier, splinesteps=16, N=3) = // Function: path_to_bezier() // Usage: -// path_to_bezier(path,[tangent],[k],[closed]); +// path_to_bezier(path, [size|relsize], [tangents], [uniform], [closed]) // Description: -// Given an input path and optional path of tangent vectors, computes a cubic (degree 3) bezier path -// that passes through every point on the input path and matches the tangent vectors. If you do not -// supply the tangent it will be computed using path_tangents. If the path is closed specify this -// by setting closed=true. If you specify the curvature parameter k it scales the tangent vectors, -// which will increase or decrease the curvature of the interpolated bezier. Negative values of k -// create loops at the corners, so they are not allowed. Sufficiently large k values will also -// produce loops. +// Given a 2d or 3d input path and optional list of tangent vectors, computes a cubic (dgree 3) bezier +// path that passes through every poin on the input path and matches the tangent vectors. If you do +// not supply the tangent it will be computed using path_tangents. If the path is closed specify this +// by setting closed=true. The size or relsize parameter determines how far the curve can deviate from +// the input path. In the case where the curve has a single hump, the size specifies the exact distance +// between the specified path and the bezier. If you give relsize then it is relative to the segment +// length (e.g. 0.05 means 5% of the segment length). In 2d when the bezier curve makes an S-curve +// the size parameter specifies the sum of the deviations of the two peaks of the curve. In 3-space +// the bezier curve may have three extrema: two maxima and one minimum. In this case the size specifies +// the sum of the maxima minus the minimum. If you do not supply the tangents then they are +// computed using path_tangents with uniform=false by default. Tangents computed on non-uniform +// data tend to display overshoots. See smooth_path for examples. // Arguments: -// path = path of points to define the bezier -// tangents = optional list of tangent vectors at every point -// k = curvature parameter, a scalar or vector to adjust curvature at each point -// closed = set to true for a closed path. Default: false -function path_to_bezier(path, tangents, k, closed=false) = - assert(is_path(path,dim=undef),"Input path is not a valid path") - assert(is_undef(tangents) || is_path(tangents,dim=len(path[0])),"Tangents must be a path of the same dimension as the input path") +// path = 2d or 3d point list that the curve must pass through +// size = absolute size specification for the curve, a number or vector +// relsize = relative size specification for the curve, a number or vector. Default: 0.1. +// tangents = tangents constraining curve direction at each point +// uniform = set to true to compute tangents with uniform=true. Default: false +// closed = true if the curve is closed . Default: false +function path_to_bezier(path, tangents, size, relsize, uniform=false, closed=false) = + assert(is_bool(closed)) + assert(is_bool(uniform)) + assert(num_defined([size,relsize])<=1, "Can't define both size and relsize") + assert(is_path(path,[2,3]),"Input path is not a valid 2d or 3d path") + assert(is_undef(tangents) || is_path(tangents,[2,3]),"Tangents must be a 2d or 3d path") assert(is_undef(tangents) || len(path)==len(tangents), "Input tangents must be the same length as the input path") let( - k = is_undef(k) ? repeat(1, len(path)) : - is_list(k) ? k : repeat(k, len(path)), - k_bad = [for(entry=k) if (entry<0) entry] - ) - assert(len(k)==len(path), "Curvature parameter k must have the same length as the path") - assert(k_bad==[], "Curvature parameter k must be a nonnegative number or list of nonnegative numbers") - let( - tangents = is_def(tangents)? tangents : deriv(path, closed=closed), + curvesize = first_defined([size,relsize,0.1]), + relative = is_undef(size), lastpt = len(path) - (closed?0:1) - ) [ - for(i=[0:lastpt-1]) each [ - path[i], - path[i]+k[i]*tangents[i]/3, - select(path,i+1)-select(k,i+1)*select(tangents,i+1)/3 - ], + ) + assert(is_num(curvesize) || len(curvesize)==lastpt, str("Size or relsize must have length ",lastpt)) + let( + sizevect = is_num(curvesize) ? repeat(curvesize, lastpt) : curvesize, + tangents = is_def(tangents) ? [for(t=tangents) let(n=norm(t)) assert(!approx(n,0),"Zero tangent vector") t/n] : + path_tangents(path, uniform=uniform, closed=closed) + ) + assert(min(sizevect)>0, "Size and relsize must be greater than zero") + [ + for(i=[0:lastpt-1]) + let( + first = path[i], + second = select(path,i+1), + seglength = norm(second-first), + dummy = assert(seglength>0, str("Path segment has zero length from index ",i," to ",i+1)), + segdir = (second-first)/seglength, + tangent1 = tangents[i], + tangent2 = -select(tangents,i+1), // Need this to point backwards, in direction of the curve + parallel = abs(tangent1*segdir) + abs(tangent2*segdir), // Total component of tangents parallel to the segment + Lmax = seglength/parallel, // May be infinity + size = relative ? sizevect[i]*seglength : sizevect[i], + normal1 = tangent1-(tangent1*segdir)*segdir, // Components of the tangents orthogonal to the segment + normal2 = tangent2-(tangent2*segdir)*segdir, + p = [ [-3 ,6,-3 ], // polynomial in power form + [ 7,-9, 2 ], + [-5, 3, 0 ], + [ 1, 0, 0 ] ]*[normal1*normal1, normal1*normal2, normal2*normal2], + uextreme = approx(norm(p),0) ? [] + : [for(root = real_roots(p)) if (root>0 && root<1) root], + distlist = [for(d=bezier_points([normal1*0, normal1, normal2, normal2*0], uextreme)) norm(d)], + scale = len(distlist)==0 ? 0 : + len(distlist)==1 ? distlist[0] + : sum(distlist) - 2*min(distlist), + Ldesired = size/scale, // This will be infinity when the polynomial is zero + L = min(Lmax, Ldesired) + ) + each [ + first, + first + L*tangent1, + second + L*tangent2 + ], select(path,lastpt) ]; diff --git a/math.scad b/math.scad index 5ec3209..a38b2e4 100644 --- a/math.scad +++ b/math.scad @@ -1014,5 +1014,121 @@ function deriv3(data, h=1, closed=false) = ]; +// Section: Complex Numbers + +// Function: C_times() +// Usage: C_times(z1,z2) +// Description: +// Multiplies two complex numbers. +function C_times(z1,z2) = [z1.x*z2.x-z1.y*z2.y,z1.x*z2.y+z1.y*z2.x]; + +// Function: C_div() +// Usage: C_div(z1,z2) +// Description: +// Divides z1 by z2. +function C_div(z1,z2) = let(den = z2.x*z2.x + z2.y*z2.y) + [(z1.x*z2.x + z1.y*z2.y)/den, (z1.y*z2.x-z1.x*z2.y)/den]; + + +// Section: Polynomials + +// Function: polynomial() +// Usage: +// polynomial(p, z) +// Description: +// Evaluates specified real polynomial, p, at the complex or real input value, z. +// The polynomial is specified as p=[a_n, a_{n-1},...,a_1,a_0] +// where a_n is the z^n coefficient. Polynomial coefficients are real. +function polynomial(p, z, k, zk, total) = + is_undef(k) ? polynomial(p, z, len(p)-1, is_num(z)? 1 : [1,0], is_num(z) ? 0 : [0,0]) : + k==-1 ? total : + polynomial(p, z, k-1, is_num(z) ? zk*z : C_times(zk,z), total+zk*p[k]); + + +// Function: poly_roots() +// Usage: +// poly_roots(p,[tol]) +// Description: +// Returns all complex roots of the specified real polynomial p. +// The polynomial is specified as p=[a_n, a_{n-1},...,a_1,a_0] +// where a_n is the z^n coefficient. The tol parameter gives +// the stopping tolerance for the iteration. The polynomial +// must have at least one non-zero coefficient. Convergence is poor +// if the polynomial has any repeated roots other than zero. +// Arguments: +// p = polynomial coefficients with higest power coefficient first +// tol = tolerance for iteration. Default: 1e-14 + +// Uses the Aberth method https://en.wikipedia.org/wiki/Aberth_method +// +// Dario Bini. "Numerical computation of polynomial zeros by means of Aberth's Method", Numerical Algorithms, Feb 1996. +// https://www.researchgate.net/publication/225654837_Numerical_computation_of_polynomial_zeros_by_means_of_Aberth's_method + +function poly_roots(p,tol=1e-14) = + assert(p!=[], "Input polynomial must have a nonzero coefficient") + assert(is_vector(p), "Input must be a vector") + p[0] == 0 ? poly_roots(slice(p,1,-1)) : // Strip leading zero coefficients + p[len(p)-1] == 0 ? [[0,0], // Strip trailing zero coefficients + each poly_roots(select(p,0,-2))] : + len(p)==1 ? [] : // Nonzero constant case has no solutions + len(p)==2 ? [[-p[1]/p[0],0]] : // Linear case needs special handling + let( + n = len(p)-1, // polynomial degree + pderiv = [for(i=[0:n-1]) p[i]*(n-i)], + + s = [for(i=[0:1:n]) abs(p[i])*(4*(n-i)+1)], // Error bound polynomial from Bini + + // Using method from: http://www.kurims.kyoto-u.ac.jp/~kyodo/kokyuroku/contents/pdf/0915-24.pdf + beta = -p[1]/p[0]/n, + r = 1+pow(abs(polynomial(p,beta)/p[0]),1/n), + init = [for(i=[0:1:n-1]) // Initial guess for roots + let(angle = 360*i/n+270/n/PI) + [beta,0]+r*[cos(angle),sin(angle)] + ] + ) + _poly_roots(p,pderiv,s,init,tol=tol); + +// p = polynomial +// pderiv = derivative polynomial of p +// z = current guess for the roots +// tol = root tolerance +// i=iteration counter +function _poly_roots(p, pderiv, s, z, tol, i=0) = + assert(i<45, str("Polyroot exceeded iteration limit. Current solution:", z)) + let( + n = len(z), + svals = [for(zk=z) tol*polynomial(s,norm(zk))], + p_of_z = [for(zk=z) polynomial(p,zk)], + done = [for(k=[0:n-1]) norm(p_of_z[k])<=svals[k]], + newton = [for(k=[0:n-1]) C_div(p_of_z[k], polynomial(pderiv,z[k]))], + zdiff = [for(k=[0:n-1]) sum([for(j=[0:n-1]) if (j!=k) C_div([1,0], z[k]-z[j])])], + w = [for(k=[0:n-1]) done[k] ? [0,0] : C_div( newton[k], + [1,0] - C_times(newton[k], zdiff[k]))] + ) + all(done) ? z : _poly_roots(p,pderiv,s,z-w,tol,i+1); + + +// Function: real_roots() +// Usage: +// real_roots(p, [eps], [tol]) +// Description: +// Returns the real roots of the specified real polynomial p. +// The polynomial is specified as p=[a_n, a_{n-1},...,a_1,a_0] +// where a_n is the x^n coefficient. This function works by +// computing the complex roots and returning those roots where +// the imaginary part is closed to zero, specifically: z.y/(1+norm(z)) < eps. Because +// of poor convergence and higher error for repeated roots, such roots may +// be missed by the algorithm because their imaginary part is large. +// Arguments: +// p = polynomial to solve as coefficient list, highest power term first +// eps = used to determine whether imaginary parts of roots are zero +// tol = tolerance for the complex polynomial root finder + +function real_roots(p,eps=EPSILON,tol=1e-14) = + let( + roots = poly_roots(p) + ) + [for(z=roots) if (abs(z.y)/(1+norm(z)) -include include -include - - -// CommonCode: -// $fn=36; - // Section: Functions - // Function: round_corners() // // Usage: @@ -103,18 +95,22 @@ include // verbose = if true display rounding scale factors that show how close roundovers are to overlapping. Default: false // // Example(Med2D): Standard circular roundover with radius the same at every point. Compare results at the different corners. +// $fn=36; // shape = [[0,0], [10,0], [15,12], [6,6], [6, 12], [-3,7]]; // polygon(round_corners(shape, radius=1)); // color("red") down(.1) polygon(shape); // Example(Med2D): Circular roundover using the "cut" specification, the same at every corner. +// $fn=36; // shape = [[0,0], [10,0], [15,12], [6,6], [6, 12], [-3,7]]; // polygon(round_corners(shape, cut=1)); // color("red") down(.1) polygon(shape); // Example(Med2D): Continous curvature roundover using "cut", still the same at every corner. The default smoothness parameter of 0.5 was too gradual for these roundovers to fit, but 0.7 works. +// $fn=36; // shape = [[0,0], [10,0], [15,12], [6,6], [6, 12], [-3,7]]; // polygon(round_corners(shape, method="smooth", cut=1, k=0.7)); // color("red") down(.1) polygon(shape); // Example(Med2D): Continuous curvature roundover using "joint", for the last time the same at every corner. Notice how small the roundovers are. +// $fn=36; // shape = [[0,0], [10,0], [15,12], [6,6], [6, 12], [-3,7]]; // polygon(round_corners(shape, method="smooth", joint=1, k=0.7)); // color("red") down(.1) polygon(shape); @@ -130,10 +126,12 @@ include // polygon(round_corners(shape, method="smooth", cut=cuts, k=k, $fs=0.1)); // color("red") down(.1) polygon(shape); // Example(Med2D): Chamfers +// $fn=36; // shape = [[0,0], [10,0], [15,12], [6,6], [6, 12], [-3,7]]; // polygon(round_corners(shape, method="chamfer", cut=1)); // color("red") down(.1) polygon(shape); // Example(Med3D): 3D printing test pieces to display different curvature shapes. You can see the discontinuity in the curvature on the "C" piece in the rendered image. +// include // ten = square(50); // cut = 5; // linear_extrude(height=14) { @@ -184,6 +182,8 @@ include // path_sweep(regular_ngon(n=36,or=.1),round_corners(list2,closed=false, method="circle", cut = 0.75)); // Example(FlatSpin): Rounding a spiral with increased rounding along the length // // Construct a square spiral path in 3D +// include +// $fn=36; // square = [[0,0],[1,0],[1,1],[0,1]]; // spiral = flatten(repeat(concat(square,reverse(square)),5)); // Squares repeat 10 times, forward and backward // squareind = [for(i=[0:9]) each [i,i,i,i]]; // Index of the square for each point @@ -389,50 +389,71 @@ function _rounding_offsets(edgespec,z_dir=1) = // Function: smooth_path() // Usage: -// smooth_path(path, [tangents], [k], [splinesteps], [closed] +// smooth_path(path, [size|relsize], [tangents], [splinesteps], [closed], [uniform]) // Description: // Smooths the input path using a cubic spline. Every segment of the path will be replaced by a cubic curve // with `splinesteps` points. The cubic interpolation will pass through every input point on the path // and will match the tangents at every point. If you do not specify tangents they will be computed using -// deriv(). See also path_to_bezier(). -// -// Note that the magnitude of the tangents affects the result. If you increase it you will get a blunter -// corner with a larger radius of curvature. Decreasing it will produce a sharp corner. You can specify -// the curvature factor `k` to adjust the curvature. It can be a scalar or a vector the same length as -// the path and is used to scale the tangent vectors. Negative values of k create loops at the corners, -// so they are not allowed. Sufficiently large k values will also produce loops. +// path_tangents with uniform=false by default. Note that setting uniform to true with non-uniform +// sampling may be desirable in some cases but tends to produces curves that overshoot the point on the path. +// +// The size or relsize parameter determines how far the curve can bend away from +// the input path. In the case where the curve has a single hump, the size specifies the exact distance +// between the specified path and the curve. If you give relsize then it is relative to the segment +// length (e.g. 0.05 means 5% of the segment length). In 2d when the spline may make an S-curve, +// in which case the size parameter specifies the sum of the deviations of the two peaks of the curve. In 3-space +// the bezier curve may have three extrema: two maxima and one minimum. In this case the size specifies +// the sum of the maxima minus the minimum. At a given segment there is a maximum size: if your size +// value is too large it will be rounded down. See also path_to_bezier(). // Arguments: // path = path to smooth -// tangents = tangent vectors of the path -// splinesteps = number of points to insert between the path points. Default: 10 -// k = curvature parameter, a scalar or vector to adjust curvature at each point -// closed = set to true for a closed path. Default: false +// size = absolute size specification for the curve, a number or vector +// relsize = relative size specification for the curve, a number or vector. Default: 0.1 +// tangents = tangents constraining curve direction at each point +// uniform = set to true to compute tangents with uniform=true. Default: false +// closed = true if the curve is closed. Default: false. // Example(2D): Original path in green, smoothed path in yellow: // color("green")stroke(square(4), width=0.1); -// stroke(smooth_path(square(4)), width=0.1); +// stroke(smooth_path(square(4),size=0.4), width=0.1); // Example(2D): Closing the path changes the end tangents -// polygon(smooth_path(square(4), closed=true)); +// polygon(smooth_path(square(4),size=0.4,closed=true)); +// Example(2D): Turning on uniform tangent calculation also changes the end derivatives: +// color("green")stroke(square(4), width=0.1); +// stroke(smooth_path(square(4),size=0.4,uniform=true), width=0.1); +// Example(2D): Here's a wide rectangle. Using size means all edges bulge the same amount, regardless of their length. +// color("green")stroke(square([10,4]), closed=true, width=0.1); +// stroke(smooth_path(square([10,4]),size=1,closed=true),width=0.1); +// Example(2D): Here's a wide rectangle. With relsize the bulge is proportional to the side length. +// color("green")stroke(square([10,4]), closed=true, width=0.1); +// stroke(smooth_path(square([10,4]),relsize=0.1,closed=true),width=0.1); +// Example(2D): Here's a wide rectangle. Settting uniform to true biases the tangents to aline more with the line sides +// color("green")stroke(square([10,4]), closed=true, width=0.1); +// stroke(smooth_path(square([10,4]),uniform=true,relsize=0.1,closed=true),width=0.1); // Example(2D): A more interesting shape: // path = [[0,0], [4,0], [7,14], [-3,12]]; -// polygon(smooth_path(path,closed=true)); -// Example(2D): Scaling the tangent data using the curvature parameter k can decrease or increase the amount of smoothing. Note this is the same as just multiplying the deriv(square(4)) by k. -// polygon(smooth_path(square(4), k=0.5,closed=true)); -// Example(2D): Or you can specify your own tangent values to alter the shape of the curve -// polygon(smooth_path(square(4),tangents=1.25*[[-2,-1], [-2,1], [1,2], [2,-1]],closed=true)); +// polygon(smooth_path(path,size=1,closed=true)); +// Example(2D): Here's the square again with less smoothing. +// polygon(smooth_path(square(4), size=.25,closed=true)); +// Example(2D): Here's the square with a size that's too big to achieve, so you get the maximum possible curve: +// polygon(smooth_path(square(4), size=4, closed=true)); +// Example(2D): You can alter the shape of the curve by specifying your own arbitrary tangent values +// polygon(smooth_path(square(4),tangents=1.25*[[-2,-1], [-4,1], [1,2], [6,-1]],size=0.4,closed=true)); +// Example(2D): Or you can give a different size for each segment +// polygon(smooth_path(square(4),size = [.4, .05, 1, .3],closed=true)); // Example(FlatSpin): Works on 3d paths as well // path = [[0,0,0],[3,3,2],[6,0,1],[9,9,0]]; -// trace_polyline(smooth_path(path),size=.3); -// Example(2D): The curve passes through all the points, but features some unexpected wiggles. These occur because the curvature is too low to change fast enough. -// path = [[0,0], [0,3], [.5,2.8], [1,2.2], [1,0]]; -// polygon(smooth_path(path)); -// color("red") move_copies(path)circle(r=.1,$fn=16); -// Example(2D): Using the k parameter is one way to fix this problem. By allowing sharper curvature (k<1) at the two points next to the problematic point we can achieve a smoother result. The other fix is to move your points. -// path = [[0,0], [0,3], [.5,2.8], [1,2.2], [1,0]]; -// polygon(smooth_path(path,k=[1,.5,1,.5,1])); -// color("red") move_copies(path)circle(r=.1,$fn=16); -function smooth_path(path, tangents, k, splinesteps=10, closed=false) = +// stroke(smooth_path(path,relsize=.1),width=.3); +// Example(2D): This shows the type of overshoot that can occur with uniform=true. You can produce overshoots like this if you supply a tangent that is difficult to connect to the adjacent points +// pts = [[-3.3, 1.7], [-3.7, -2.2], [3.8, -4.8], [-0.9, -2.4]]; +// stroke(smooth_path(pts, uniform=true, relsize=0.1),width=.1); +// color("red")move_copies(pts)circle(r=.15,$fn=12); +// Example(2D): With the default of uniform false no overshoot occurs. Note that the shape of the curve is quite different. +// pts = [[-3.3, 1.7], [-3.7, -2.2], [3.8, -4.8], [-0.9, -2.4]]; +// stroke(smooth_path(pts, uniform=false, relsize=0.1),width=.1); +// color("red")move_copies(pts)circle(r=.15,$fn=12); +function smooth_path(path, tangents, size, relsize, splinesteps=10, uniform=false, closed=false) = let ( - bez = path_to_bezier(path, tangents, k=k, closed=closed) + bez = path_to_bezier(path, tangents=tangents, size=size, relsize=relsize, uniform=uniform, closed=closed) ) bezier_polyline(bez,splinesteps=splinesteps); @@ -1047,7 +1068,7 @@ function _remove_undefined_vals(list) = // get the expected rounding along the path, decrease `maxstep` and if the curves created by `os_round()` are too coarse, adjust $fn or $fs. // // Arguments: -// path = path that defines the stroke +// path = 2d path that defines the stroke // width = width of the stroke, a scalar or a vector of 2 values giving the offset from the path. Default: 1 // rounded = set to true to use rounded offsets, false to use sharp (delta) offsets. Default: true // chamfer = set to true to use chamfers when `rounded=false`. Default: false @@ -1141,6 +1162,7 @@ function _remove_undefined_vals(list) = // right(12) // offset_stroke(path, width=1, closed=true); function offset_stroke(path, width=1, rounded=true, start="flat", end="flat", check_valid=true, quality=1, maxstep=0.1, chamfer=false, closed=false) = + assert(is_path(path,2),"path is not a 2d path") let(closedok = !closed || (is_undef(start) && is_undef(end))) assert(closedok, "Parameters `start` and `end` not allowed with closed path") let(