diff --git a/beziers.scad b/beziers.scad index 4217f39..adb6673 100644 --- a/beziers.scad +++ b/beziers.scad @@ -327,26 +327,40 @@ function bezier_polyline(bezier, splinesteps=16, N=3) = let( ); + // Function: path_to_bezier() // Usage: -// path_to_bezier(path,[tangent],[closed]); +// path_to_bezier(path,[tangent],[k],[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. +// 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. // 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, closed=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") + 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), lastpt = len(path) - (closed?0:1) ) - [for(i=[0:lastpt-1]) each [path[i], path[i]+tangents[i]/3, select(path,i+1)-select(tangents,i+1)/3], + [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], select(path,lastpt)]; diff --git a/geometry.scad b/geometry.scad index ec251d1..34a3a23 100644 --- a/geometry.scad +++ b/geometry.scad @@ -548,7 +548,7 @@ function triangle_area(a,b,c) = // plane3pt(p1, p2, p3); // Description: // Generates the cartesian equation of a plane from three non-collinear points on the plane. -// Returns [A,B,C,D] where Ax+By+Cz=D is the equation of a plane. +// Returns [A,B,C,D] where Ax + By + Cz = D is the equation of a plane. // Arguments: // p1 = The first point on the plane. // p2 = The second point on the plane. @@ -583,6 +583,46 @@ function plane3pt_indexed(points, i1, i2, i3) = ) plane3pt(p1,p2,p3); +// Function: plane_intersection() +// Usage: +// plane_intersection(plane1, plane2, [plane3]) +// Description: +// Compute the point which is the intersection of the three planes, or the line intersection of two planes. +// If you give three planes the intersection is returned as a point. If you give two planes the intersection +// is returned as a list of two points on the line of intersection. If any of the input planes are parallel +// then returns undef. +function plane_intersection(plane1,plane2,plane3) = + is_def(plane3) ? + let ( + matrix = [for(p=[plane1,plane2,plane3]) select(p,0,2)], + rhs = [for(p=[plane1,plane2,plane3]) p[3]] + ) + linear_solve(matrix,rhs) + : + let( + normal = cross(plane_normal(plane1), plane_normal(plane2)) + ) + approx(normal,0) ? undef : + let( + matrix = [for(p=[plane1,plane2]) select(p,0,2)], + rhs = [for(p=[plane1,plane2]) p[3]], + point = linear_solve(matrix,rhs), + dd=echo(point=point, normal=normal) + ) + [point, point+normal]; + + +// Function: plane_from_normal() +// Usage: +// plane_from_normal(normal, pt) +// Description: +// Returns a plane defined by a normal vector and a point. +// Example: +// plane_from_normal([0,0,1], [2,2,2]); // Returns the xy plane passing through the point (2,2,2) +function plane_from_normal(normal, pt) = + concat(normal, [normal*pt]); + + // Function: plane_from_pointslist() // Usage: // plane_from_pointslist(points); diff --git a/math.scad b/math.scad index 74a4a8b..f3e1b99 100644 --- a/math.scad +++ b/math.scad @@ -538,6 +538,8 @@ function mean(v) = sum(v)/len(v); // the least squares solution is returned. If A is underdetermined, the minimal norm solution is returned. // If A is rank deficient or singular then linear_solve returns `undef`. function linear_solve(A,b) = + assert(is_matrix(A)) + assert(is_vector(b)) let( dim = array_dim(A), m=dim[0], n=dim[1] @@ -569,6 +571,7 @@ function submatrix(M,ind1,ind2) = [for(i=ind1) [for(j=ind2) M[i][j] ] ]; // Calculates the QR factorization of the input matrix A and returns it as the list [Q,R]. This factorization can be // used to solve linear systems of equations. function qr_factor(A) = + assert(is_matrix(A)) let( dim = array_dim(A), m = dim[0], @@ -673,6 +676,19 @@ function determinant(M) = ); +// Function: is_matrix() +// Usage: +// is_matrix(A,[m],[n]) +// Description: +// Returns true if A is a numeric matrix of height m and width n. If m or n +// are omitted or set to undef then true is returned for any positive dimension. +function is_matrix(A,m,n) = + is_list(A) && len(A)>0 && + (is_undef(m) || len(A)==m) && + is_vector(A[0]) && + (is_undef(n) || len(A[0])==n) && + is_consistent(A); + // Section: Comparisons and Logic diff --git a/rounding.scad b/rounding.scad index 4987f21..2ccf7a2 100644 --- a/rounding.scad +++ b/rounding.scad @@ -9,10 +9,10 @@ // ``` ////////////////////////////////////////////////////////////////////// -include -include -include -include +include +include +include +include // CommonCode: @@ -409,16 +409,23 @@ function _rounding_offsets(edgespec,z_dir=1) = // Function: smooth_path() // Usage: -// smooth_path(path, [tangents], [splinesteps], [closed] +// smooth_path(path, [tangents], [k], [splinesteps], [closed] // 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(). Note that the magnitude of the tangents affects the result. See also path_to_bezier(). +// 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. // 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 // Example(2D): Original path in green, smoothed path in yellow: // color("green")stroke(square(4), width=0.1); @@ -428,21 +435,30 @@ function _rounding_offsets(edgespec,z_dir=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 can decrease or increase the amount of smoothing: -// shape = square(4); -// polygon(smooth_path(shape, tangents=0.5*deriv(shape, closed=true),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)); // 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); -function smooth_path(path, tangents, splinesteps=10, closed=false) = - let( - bez = path_to_bezier(path, tangents=tangents, closed=closed) +// 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") place_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") place_copies(path)circle(r=.1,$fn=16); +function smooth_path(path, tangents, k, splinesteps=10, closed=false) = + let ( + bez = path_to_bezier(path, tangents, k=k, closed=closed) ) bezier_polyline(bez,splinesteps=splinesteps); + // Module: offset_sweep() // // Description: diff --git a/wiring.scad b/wiring.scad index 0c97da4..8f7d1ee 100644 --- a/wiring.scad +++ b/wiring.scad @@ -9,7 +9,7 @@ ////////////////////////////////////////////////////////////////////// -include +include // Section: Functions