grey/BOSL2
1
0
Fork 0
mirror of https://github.com/BelfrySCAD/BOSL2.git synced 2024-12-29 16:29:40 +00:00
Code Issues Projects Releases Packages Wiki Activity Actions

Added polynomial root finding and complex math to math.scad. Changed

Browse source 
path_to_bezier to new algorithm (which uses polynomial roots).  And
updated path_smooth to use the new code.

Also removed extra (?) include<skin.scad> and
include<strings.scad> 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.
This commit is contained in:
Adrian Mariano 2020-06-18 17:50:25 -04:00
parent 5b39cefe7b
commit 7aab27266f
3 changed files with 245 additions and 68 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

95
beziers.scad
View file

@ -515,39 +515,78 @@ 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)
];

116
math.scad
View file

@ -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))<eps) z.x];
// vim: expandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap

98
rounding.scad
View file

@ -10,18 +10,10 @@
//////////////////////////////////////////////////////////////////////
include <beziers.scad>
include <strings.scad>
include <structs.scad>
include <skin.scad>
// CommonCode:
// $fn=36;
// Section: Functions
// Function: round_corners()
//
// Usage:
@ -103,18 +95,22 @@ include <skin.scad>
// 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 <skin.scad>
// 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<skin.scad>
// ten = square(50);
// cut = 5;
// linear_extrude(height=14) {
@ -184,6 +182,8 @@ include <skin.scad>
// 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<skin.scad>
// $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().
// 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.
//
// 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.
// 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(