From 386faf9092e58a9124831714fa174e97122036ec Mon Sep 17 00:00:00 2001 From: Cameron K Brooks Date: Fri, 26 Sep 2025 18:36:20 -0400 Subject: [PATCH 01/32] fmt: format all files with topiary --- hull.scad | 516 ++++++++++++++++++++++--------------------- linalg.scad | 40 ++-- lists.scad | 20 +- mirror.scad | 25 +-- morphology.scad | 140 ++++++------ se3.scad | 118 +++++----- shapes.scad | 23 +- so3.scad | 139 ++++++------ spline.scad | 157 +++++++------ trajectory.scad | 77 ++++--- trajectory_path.scad | 172 ++++++++------- transformations.scad | 44 ++-- 12 files changed, 762 insertions(+), 709 deletions(-) diff --git a/hull.scad b/hull.scad index 5e0302e..3323ac9 100644 --- a/hull.scad +++ b/hull.scad @@ -11,314 +11,318 @@ // polygon [i1,i2,i3,...], otherwise a triangular // polyhedron [[i1,i2,i3],[i2,i3,i4],...] -function hull(points) = - !(len(points) > 0) ? [] : - len(points[0]) == 2 ? convexhull2d(points) : - len(points[0]) == 3 ? convexhull3d(points) : []; +function hull(points) = + !(len(points) > 0) ? [] + : len(points[0]) == 2 ? convexhull2d(points) + : len(points[0]) == 3 ? convexhull3d(points) : []; epsilon = 1e-9; // 2d version function convexhull2d(points) = -len(points) < 3 ? [] : let( - a=0, b=1, - - c = find_first_noncollinear([a,b], points, 2) - -) c == len(points) ? convexhull_collinear(points) : let( - - remaining = [ for (i = [2:len(points)-1]) if (i != c) i ], - - polygon = area_2d(points[a], points[b], points[c]) > 0 ? [a,b,c] : [b,a,c] - -) convex_hull_iterative_2d(points, polygon, remaining); - + len(points) < 3 ? [] + : let ( + a = 0, + b = 1, + c = find_first_noncollinear([a, b], points, 2) + ) c == len(points) ? convexhull_collinear(points) + : let ( + remaining = [for (i = [2:len(points) - 1]) if (i != c) i], + polygon = area_2d(points[a], points[b], points[c]) > 0 ? [a, b, c] : [b, a, c] + ) convex_hull_iterative_2d(points, polygon, remaining); // Adds the remaining points one by one to the convex hull -function convex_hull_iterative_2d(points, polygon, remaining, i_=0) = i_ >= len(remaining) ? polygon : - let ( - // pick a point - i = remaining[i_], - - // find the segments that are in conflict with the point (point not inside) - conflicts = find_conflicting_segments(points, polygon, points[i]) - - // no conflicts, skip point and move on - ) len(conflicts) == 0 ? convex_hull_iterative_2d(points, polygon, remaining, i_+1) : let( - - // find the first conflicting segment and the first not conflicting - // conflict will be sorted, if not wrapping around, do it the easy way - polygon = remove_conflicts_and_insert_point(polygon, conflicts, i) - ) convex_hull_iterative_2d( - points, - polygon, - remaining, - i_+1 - ); - -function find_conflicting_segments(points, polygon, point) = [ - for (i = [0:len(polygon)-1]) let(j = (i+1) % len(polygon)) - if (area_2d(points[polygon[i]], points[polygon[j]], point) < 0) - i -]; +function convex_hull_iterative_2d(points, polygon, remaining, i_ = 0) = + i_ >= len(remaining) ? polygon + : let ( + // pick a point + i = remaining[i_], + + // find the segments that are in conflict with the point (point not inside) + conflicts = find_conflicting_segments(points, polygon, points[i]) + + // no conflicts, skip point and move on + ) len(conflicts) == 0 ? convex_hull_iterative_2d(points, polygon, remaining, i_ + 1) + : let ( + + // find the first conflicting segment and the first not conflicting + // conflict will be sorted, if not wrapping around, do it the easy way + polygon = remove_conflicts_and_insert_point(polygon, conflicts, i) + ) convex_hull_iterative_2d( + points, + polygon, + remaining, + i_ + 1 + ); + +function find_conflicting_segments(points, polygon, point) = + [ + for (i = [0:len(polygon) - 1]) let (j = (i + 1) % len(polygon)) if (area_2d(points[polygon[i]], points[polygon[j]], point) < 0) i, + ]; // remove the conflicting segments from the polygon -function remove_conflicts_and_insert_point(polygon, conflicts, point) = - conflicts[0] == 0 ? let( - nonconflicting = [ for(i = [0:len(polygon)-1]) if (!contains(conflicts, i)) i ], - new_indices = concat(nonconflicting, (nonconflicting[len(nonconflicting)-1]+1) % len(polygon)), - polygon = concat([ for (i = new_indices) polygon[i] ], point) - ) polygon : let( - prior_to_first_conflict = [ for(i = [0:1:min(conflicts)]) polygon[i] ], - after_last_conflict = [ for(i = [max(conflicts)+1:1:len(polygon)-1]) polygon[i] ], - polygon = concat(prior_to_first_conflict, point, after_last_conflict) - ) polygon; - +function remove_conflicts_and_insert_point(polygon, conflicts, point) = + conflicts[0] == 0 ? let ( + nonconflicting = [for (i = [0:len(polygon) - 1]) if (!contains(conflicts, i)) i], + new_indices = concat(nonconflicting, (nonconflicting[len(nonconflicting) - 1] + 1) % len(polygon)), + polygon = concat([for (i = new_indices) polygon[i]], point) + ) polygon + : let ( + prior_to_first_conflict = [for (i = [0:1:min(conflicts)]) polygon[i]], + after_last_conflict = [for (i = [max(conflicts) + 1:1:len(polygon) - 1]) polygon[i]], + polygon = concat(prior_to_first_conflict, point, after_last_conflict) + ) polygon; // 3d version -function convexhull3d(points) = -len(points) < 3 ? [ for(i = [0:1:len(points)-1]) i ] : let ( - - // start with a single triangle - a=0, b=1, c=2, - plane = plane(points,a,b,c), - - d = find_first_noncoplanar(plane, points, 3) - -) d == len(points) ? /* all coplanar*/ let ( - - pts2d = [ for (p = points) plane_project(p, points[a], points[b], points[c]) ], - hull2d = convexhull2d(pts2d) - -) hull2d : let( - - remaining = [for (i = [3:len(points)-1]) if (i != d) i], - - // Build an initial tetrahedron - - // swap b,c if d is in front of triangle t - bc = in_front(plane, points[d]) ? [c,b] : [b,c], - b = bc[0], c = bc[1], - - triangles = [ - [a,b,c], - [d,b,a], - [c,d,a], - [b,d,c], - ], - - // calculate the plane equations - planes = [ for (t = triangles) plane(points, t[0], t[1], t[2]) ] - -) convex_hull_iterative(points, triangles, planes, remaining); +function convexhull3d(points) = + len(points) < 3 ? [for (i = [0:1:len(points) - 1]) i] + : let ( + + // start with a single triangle + a = 0, + b = 1, + c = 2, + plane = plane(points, a, b, c), + d = find_first_noncoplanar(plane, points, 3) + ) d == len(points) ? /* all coplanar*/ let ( + pts2d = [for (p = points) plane_project(p, points[a], points[b], points[c])], + hull2d = convexhull2d(pts2d) + ) hull2d + : let ( + remaining = [for (i = [3:len(points) - 1]) if (i != d) i], + + // Build an initial tetrahedron + + // swap b,c if d is in front of triangle t + bc = in_front(plane, points[d]) ? [c, b] : [b, c], + b = bc[0], + c = bc[1], + triangles = [ + [a, b, c], + [d, b, a], + [c, d, a], + [b, d, c], + ], + + // calculate the plane equations + planes = [for (t = triangles) plane(points, t[0], t[1], t[2])] + ) convex_hull_iterative(points, triangles, planes, remaining); // A plane equation (normal, offset) -function plane(points, a, b, c) = let( - normal = unit(cross(points[c]-points[a], points[b]-points[a])) -) [ - normal, - normal * points[a] -]; +function plane(points, a, b, c) = + let ( + normal = unit(cross(points[c] - points[a], points[b] - points[a])) + ) [ + normal, + normal * points[a], + ]; // Adds the remaining points one by one to the convex hull -function convex_hull_iterative(points, triangles, planes, remaining, i_=0) = i_ >= len(remaining) ? triangles : - let ( - // pick a point - i = remaining[i_], - - // find the triangles that are in conflict with the point (point not inside) - conflicts = find_conflicts(points[i], planes), - - // for all triangles that are in conflict, collect their halfedges - halfedges = [ - for(c = conflicts) - for(i = [0:2]) let(j = (i+1)%3) - [triangles[c][i], triangles[c][j]] - ], - - // find the outer perimeter of the set of conflicting triangles - horizon = remove_internal_edges(halfedges), - - // generate a new triangle for each horizon halfedge together with the picked point i - new_triangles = [ for (h = horizon) concat(h,i) ], - - // calculate the corresponding plane equations - new_planes = [ for (t = new_triangles) plane(points, t[0], t[1], t[2]) ] - - ) convex_hull_iterative( - points, - // remove the conflicting triangles and add the new ones - concat(remove_elements(triangles, conflicts), new_triangles), - concat(remove_elements(planes, conflicts), new_planes), - remaining, - i_+1 - ); - -function convexhull_collinear(points) = let( - n = points[1] - points[0], - a = points[0], - points1d = [ for(p = points) (p-a)*n ], - min_i = min_index(points1d), - max_i = max_index(points1d) -) [ min_i, max_i ]; - -function min_index(values,min_,min_i_,i_) = - i_ == undef ? min_index(values,values[0],0,1) : - i_ >= len(values) ? min_i_ : - values[i_] < min_ ? min_index(values,values[i_],i_,i_+1) - : min_index(values,min_,min_i_,i_+1); - -function max_index(values,max_,max_i_,i_) = - i_ == undef ? max_index(values,values[0],0,1) : - i_ >= len(values) ? max_i_ : - values[i_] > max_ ? max_index(values,values[i_],i_,i_+1) - : max_index(values,max_,max_i_,i_+1); - -function remove_elements(array, elements) = [ - for (i = [0:len(array)-1]) - if (!search(i, elements)) - array[i] -]; - -function remove_internal_edges(halfedges) = [ - for (h = halfedges) - if (!contains(halfedges, reverse(h))) - h -]; - -function plane_project(point, a, b, c) = let( - u = b-a, - v = c-a, - n = cross(u,v), - w = cross(n,u), - relpoint = point-a -) [relpoint * u, relpoint * w]; - -function plane_unproject(point, a, b, c) = let( - u = b-a, - v = c-a, - n = cross(u,v), - w = cross(n,u) -) a + point[0] * u + point[1] * w; - -function reverse(arr) = [ for (i = [len(arr)-1:-1:0]) arr[i] ]; - -function contains(arr, element) = search([element],arr)[0] != [] ? true : false; - -function find_conflicts(point, planes) = [ - for (i = [0:len(planes)-1]) - if (in_front(planes[i], point)) - i -]; - -function find_first_noncollinear(line, points, i) = - i >= len(points) ? len(points) : - collinear(points[line[0]], - points[line[1]], - points[i]) ? find_first_noncollinear(line, points, i+1) - : i; - -function find_first_noncoplanar(plane, points, i) = - i >= len(points) ? len(points) : - coplanar(plane, points[i]) ? find_first_noncoplanar(plane, points, i+1) - : i; +function convex_hull_iterative(points, triangles, planes, remaining, i_ = 0) = + i_ >= len(remaining) ? triangles + : let ( + // pick a point + i = remaining[i_], + + // find the triangles that are in conflict with the point (point not inside) + conflicts = find_conflicts(points[i], planes), + + // for all triangles that are in conflict, collect their halfedges + halfedges = [ + for (c = conflicts) for (i = [0:2]) let (j = (i + 1) % 3) [triangles[c][i], triangles[c][j]], + ], + + // find the outer perimeter of the set of conflicting triangles + horizon = remove_internal_edges(halfedges), + + // generate a new triangle for each horizon halfedge together with the picked point i + new_triangles = [for (h = horizon) concat(h, i)], + + // calculate the corresponding plane equations + new_planes = [for (t = new_triangles) plane(points, t[0], t[1], t[2])] + ) convex_hull_iterative( + points, + // remove the conflicting triangles and add the new ones + concat(remove_elements(triangles, conflicts), new_triangles), + concat(remove_elements(planes, conflicts), new_planes), + remaining, + i_ + 1 + ); + +function convexhull_collinear(points) = + let ( + n = points[1] - points[0], + a = points[0], + points1d = [for (p = points) (p - a) * n], + min_i = min_index(points1d), + max_i = max_index(points1d) + ) [min_i, max_i]; + +function min_index(values, min_, min_i_, i_) = + i_ == undef ? min_index(values, values[0], 0, 1) + : i_ >= len(values) ? min_i_ + : values[i_] < min_ ? min_index(values, values[i_], i_, i_ + 1) + : min_index(values, min_, min_i_, i_ + 1); + +function max_index(values, max_, max_i_, i_) = + i_ == undef ? max_index(values, values[0], 0, 1) + : i_ >= len(values) ? max_i_ + : values[i_] > max_ ? max_index(values, values[i_], i_, i_ + 1) + : max_index(values, max_, max_i_, i_ + 1); + +function remove_elements(array, elements) = + [ + for (i = [0:len(array) - 1]) if (!search(i, elements)) array[i], + ]; + +function remove_internal_edges(halfedges) = + [ + for (h = halfedges) if (!contains(halfedges, reverse(h))) h, + ]; + +function plane_project(point, a, b, c) = + let ( + u = b - a, + v = c - a, + n = cross(u, v), + w = cross(n, u), + relpoint = point - a + ) [relpoint * u, relpoint * w]; + +function plane_unproject(point, a, b, c) = + let ( + u = b - a, + v = c - a, + n = cross(u, v), + w = cross(n, u) + ) a + point[0] * u + point[1] * w; + +function reverse(arr) = [for (i = [len(arr) - 1:-1:0]) arr[i]]; + +function contains(arr, element) = search([element], arr)[0] != [] ? true : false; + +function find_conflicts(point, planes) = + [ + for (i = [0:len(planes) - 1]) if (in_front(planes[i], point)) i, + ]; + +function find_first_noncollinear(line, points, i) = + i >= len(points) ? len(points) + : collinear( + points[line[0]], + points[line[1]], + points[i] + ) ? find_first_noncollinear(line, points, i + 1) + : i; + +function find_first_noncoplanar(plane, points, i) = + i >= len(points) ? len(points) + : coplanar(plane, points[i]) ? find_first_noncoplanar(plane, points, i + 1) + : i; function distance(plane, point) = plane[0] * point - plane[1]; function in_front(plane, point) = distance(plane, point) > epsilon; -function coplanar(plane, point) = abs(distance(plane,point)) <= epsilon; +function coplanar(plane, point) = abs(distance(plane, point)) <= epsilon; -function unit(v) = v/norm(v); +function unit(v) = v / norm(v); -function area_2d(a,b,c) = ( - a[0] * (b[1] - c[1]) + - b[0] * (c[1] - a[1]) + - c[0] * (a[1] - b[1])) / 2; +function area_2d(a, b, c) = + ( + a[0] * (b[1] - c[1]) + b[0] * (c[1] - a[1]) + c[0] * (a[1] - b[1]) + ) / 2; -function collinear(a,b,c) = abs(area_2d(a,b,c)) < epsilon; +function collinear(a, b, c) = abs(area_2d(a, b, c)) < epsilon; -function spherical(cartesian) = [ +function spherical(cartesian) = + [ atan2(cartesian[1], cartesian[0]), - asin(cartesian[2]) -]; - -function cartesian(spherical) = [ - cos(spherical[1]) * cos(spherical[0]), - cos(spherical[1]) * sin(spherical[0]), - sin(spherical[1]) -]; + asin(cartesian[2]), + ]; +function cartesian(spherical) = + [ + cos(spherical[1]) * cos(spherical[0]), + cos(spherical[1]) * sin(spherical[0]), + sin(spherical[1]), + ]; /// TESTCODE - phi = 1.618033988749895; -testpoints_on_sphere = [ for(p = - [ - [1,phi,0], [-1,phi,0], [1,-phi,0], [-1,-phi,0], - [0,1,phi], [0,-1,phi], [0,1,-phi], [0,-1,-phi], - [phi,0,1], [-phi,0,1], [phi,0,-1], [-phi,0,-1] - ]) - unit(p) +testpoints_on_sphere = [ + for ( + p = [ + [1, phi, 0], + [-1, phi, 0], + [1, -phi, 0], + [-1, -phi, 0], + [0, 1, phi], + [0, -1, phi], + [0, 1, -phi], + [0, -1, -phi], + [phi, 0, 1], + [-phi, 0, 1], + [phi, 0, -1], + [-phi, 0, -1], + ] + ) unit(p), ]; -testpoints_spherical = [ for(p = testpoints_on_sphere) spherical(p) ]; -testpoints_circular = [ for(a = [0:15:360-epsilon]) [cos(a),sin(a)] ]; +testpoints_spherical = [for (p = testpoints_on_sphere) spherical(p)]; +testpoints_circular = [for (a = [0:15:360 - epsilon]) [cos(a), sin(a)]]; -testpoints_coplanar = let(u = unit([1,3,7]), v = unit([-2,1,-2])) [ for(i = [1:10]) rands(-1,1,1)[0] * u + rands(-1,1,1)[0] * v ]; +testpoints_coplanar = let (u = unit([1, 3, 7]), v = unit([-2, 1, -2])) [for (i = [1:10]) rands(-1, 1, 1)[0] * u + rands(-1, 1, 1)[0] * v]; -testpoints_collinear_2d = let(u = unit([5,3])) [ for(i = [1:20]) rands(-1,1,1)[0] * u ]; -testpoints_collinear_3d = let(u = unit([5,3,-5])) [ for(i = [1:20]) rands(-1,1,1)[0] * u ]; +testpoints_collinear_2d = let (u = unit([5, 3])) [for (i = [1:20]) rands(-1, 1, 1)[0] * u]; +testpoints_collinear_3d = let (u = unit([5, 3, -5])) [for (i = [1:20]) rands(-1, 1, 1)[0] * u]; -testpoints2d = 20 * [for (i = [1:10]) concat(rands(-1,1,2))]; -testpoints3d = 20 * [for (i = [1:50]) concat(rands(-1,1,3))]; +testpoints2d = 20 * [for (i = [1:10]) concat(rands(-1, 1, 2))]; +testpoints3d = 20 * [for (i = [1:50]) concat(rands(-1, 1, 3))]; // All points are on the sphere, no point should be red -translate([-50,0]) visualize_hull(20*testpoints_on_sphere); +translate([-50, 0]) visualize_hull(20 * testpoints_on_sphere); // 2D points -translate([50,0]) visualize_hull(testpoints2d); +translate([50, 0]) visualize_hull(testpoints2d); // All points on a circle, no point should be red -translate([0,50]) visualize_hull(20*testpoints_circular); +translate([0, 50]) visualize_hull(20 * testpoints_circular); // All points 3d but collinear -translate([0,-50]) visualize_hull(20*testpoints_coplanar); +translate([0, -50]) visualize_hull(20 * testpoints_coplanar); // Collinear -translate([50,50]) visualize_hull(20*testpoints_collinear_2d); +translate([50, 50]) visualize_hull(20 * testpoints_collinear_2d); // Collinear -translate([-50,50]) visualize_hull(20*testpoints_collinear_3d); +translate([-50, 50]) visualize_hull(20 * testpoints_collinear_3d); // 3D points visualize_hull(testpoints3d); - module visualize_hull(points) { - hull = hull(points); - - %if (len(hull) > 0 && len(hull[0]) > 0) - polyhedron(points=points, faces = hull); - else - polyhedron(points=points, faces = [hull]); - - for (i = [0:len(points)-1]) assign(p = points[i], $fn = 16) { - translate(p) { - if (hull_contains_index(hull,i)) { - color("blue") sphere(1); - } else { - color("red") sphere(1); - } - } - } - - function hull_contains_index(hull, index) = - search(index,hull,1,0) || - search(index,hull,1,1) || - search(index,hull,1,2); - + hull = hull(points); + + if (len(hull) > 0 && len(hull[0]) > 0) + %polyhedron(points=points, faces=hull); + else + %polyhedron(points=points, faces=[hull]); + + for (i = [0:len(points) - 1]) + assign (p = points[i], $fn = 16) { + translate(p) { + if (hull_contains_index(hull, i)) { + color("blue") sphere(1); + } else { + color("red") sphere(1); + } + } + } + + function hull_contains_index(hull, index) = + search(index, hull, 1, 0) || search(index, hull, 1, 1) || search(index, hull, 1, 2); } diff --git a/linalg.scad b/linalg.scad index 959a38e..573652d 100644 --- a/linalg.scad +++ b/linalg.scad @@ -4,29 +4,31 @@ //function cross(x,y) = [x[1]*y[2]-x[2]*y[1], x[2]*y[0]-x[0]*y[2], x[0]*y[1]-x[1]*y[0]]; //function norm(v) = sqrt(v*v); -function vec3(p) = len(p) < 3 ? concat(p,0) : p; -function vec4(p) = let (v3=vec3(p)) len(v3) < 4 ? concat(v3,1) : v3; -function unit(v) = v/norm(v); +function vec3(p) = len(p) < 3 ? concat(p, 0) : p; +function vec4(p) = let (v3 = vec3(p)) len(v3) < 4 ? concat(v3, 1) : v3; +function unit(v) = v / norm(v); -function identity3()=[[1,0,0],[0,1,0],[0,0,1]]; -function identity4()=[[1,0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,1]]; +function identity3() = [[1, 0, 0], [0, 1, 0], [0, 0, 1]]; +function identity4() = [[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]]; - -function take3(v) = [v[0],v[1],v[2]]; -function tail3(v) = [v[3],v[4],v[5]]; -function rotation_part(m) = [take3(m[0]),take3(m[1]),take3(m[2])]; +function take3(v) = [v[0], v[1], v[2]]; +function tail3(v) = [v[3], v[4], v[5]]; +function rotation_part(m) = [take3(m[0]), take3(m[1]), take3(m[2])]; function rot_trace(m) = m[0][0] + m[1][1] + m[2][2]; -function rot_cos_angle(m) = (rot_trace(m)-1)/2; +function rot_cos_angle(m) = (rot_trace(m) - 1) / 2; -function rotation_part(m) = [take3(m[0]),take3(m[1]),take3(m[2])]; -function translation_part(m) = [m[0][3],m[1][3],m[2][3]]; -function transpose_3(m) = [[m[0][0],m[1][0],m[2][0]],[m[0][1],m[1][1],m[2][1]],[m[0][2],m[1][2],m[2][2]]]; -function transpose_4(m) = [[m[0][0],m[1][0],m[2][0],m[3][0]], - [m[0][1],m[1][1],m[2][1],m[3][1]], - [m[0][2],m[1][2],m[2][2],m[3][2]], - [m[0][3],m[1][3],m[2][3],m[3][3]]]; +function rotation_part(m) = [take3(m[0]), take3(m[1]), take3(m[2])]; +function translation_part(m) = [m[0][3], m[1][3], m[2][3]]; +function transpose_3(m) = [[m[0][0], m[1][0], m[2][0]], [m[0][1], m[1][1], m[2][1]], [m[0][2], m[1][2], m[2][2]]]; +function transpose_4(m) = + [ + [m[0][0], m[1][0], m[2][0], m[3][0]], + [m[0][1], m[1][1], m[2][1], m[3][1]], + [m[0][2], m[1][2], m[2][2], m[3][2]], + [m[0][3], m[1][3], m[2][3], m[3][3]], + ]; function invert_rt(m) = construct_Rt(transpose_3(rotation_part(m)), -(transpose_3(rotation_part(m)) * translation_part(m))); -function construct_Rt(R,t) = [concat(R[0],t[0]),concat(R[1],t[1]),concat(R[2],t[2]),[0,0,0,1]]; +function construct_Rt(R, t) = [concat(R[0], t[0]), concat(R[1], t[1]), concat(R[2], t[2]), [0, 0, 0, 1]]; // Hadamard product of n-dimensional arrays -function hadamard(a,b) = !(len(a)>0) ? a*b : [ for(i = [0:len(a)-1]) hadamard(a[i],b[i]) ]; +function hadamard(a, b) = !(len(a) > 0) ? a * b : [for (i = [0:len(a) - 1]) hadamard(a[i], b[i])]; diff --git a/lists.scad b/lists.scad index 0d8e2b4..683e11b 100644 --- a/lists.scad +++ b/lists.scad @@ -5,22 +5,21 @@ flatten([[0,1],[2,3]]) => [0,1,2,3] */ -function flatten(list) = [ for (i = list, v = i) v ]; - +function flatten(list) = [for (i = list, v = i) v]; /*! Creates a list from a range: range([0:2:6]) => [0,2,4,6] */ -function range(r) = [ for(x=r) x ]; +function range(r) = [for (x = r) x]; /*! Reverses a list: reverse([1,2,3]) => [3,2,1] */ -function reverse(list) = [for (i = [len(list)-1:-1:0]) list[i]]; +function reverse(list) = [for (i = [len(list) - 1:-1:0]) list[i]]; /*! Extracts a subarray from index begin (inclusive) to end (exclusive) @@ -28,21 +27,20 @@ function reverse(list) = [for (i = [len(list)-1:-1:0]) list[i]]; subarray([1,2,3,4], 1, 2) => [2,3] */ -function subarray(list,begin=0,end=-1) = [ - let(end = end < 0 ? len(list) : end) - for (i = [begin : 1 : end-1]) - list[i] -]; +function subarray(list, begin = 0, end = -1) = + [ + let (end = end < 0 ? len(list) : end) for (i = [begin:1:end - 1]) list[i], + ]; /*! Returns a copy of a list with the element at index i set to x set([1,2,3,4], 2, 5) => [1,2,5,4] */ -function set(list, i, x) = [for (i_=[0:len(list)-1]) i == i_ ? x : list[i_]]; +function set(list, i, x) = [for (i_ = [0:len(list) - 1]) i == i_ ? x : list[i_]]; /*! Remove element from the list by index. remove([4,3,2,1],1) => [4,2,1] */ -function remove(list, i) = [for (i_=[0:1:len(list)-2]) list[i_ < i ? i_ : i_ + 1]]; +function remove(list, i) = [for (i_ = [0:1:len(list) - 2]) list[i_ < i ? i_ : i_ + 1]]; diff --git a/mirror.scad b/mirror.scad index d7665d6..40f4017 100644 --- a/mirror.scad +++ b/mirror.scad @@ -7,24 +7,23 @@ // mirror_y() // mirror_z() - module mirror_x() { - union() { - children(); - scale([-1,1,1]) children(); - } + union() { + children(); + scale([-1, 1, 1]) children(); + } } module mirror_y() { - union() { - children(); - scale([1,-1,1]) children(); - } + union() { + children(); + scale([1, -1, 1]) children(); + } } module mirror_z() { - union() { - children(); - scale([1,1,-1]) children(); - } + union() { + children(); + scale([1, 1, -1]) children(); + } } diff --git a/morphology.scad b/morphology.scad index 4bdd28b..53871d7 100644 --- a/morphology.scad +++ b/morphology.scad @@ -12,98 +12,100 @@ // - negative values of d places the shell on the inside // - center=true and positive d places the shell centered on the edge -module outset(d=1) { - // Bug workaround for older OpenSCAD versions - if (version_num() < 20130424) render() outset_extruded(d) children(); - else minkowski() { - circle(r=d); - children(); - } +module outset(d = 1) { + // Bug workaround for older OpenSCAD versions + if (version_num() < 20130424) + render() outset_extruded(d) children(); + else + minkowski() { + circle(r=d); + children(); + } } -module outset_extruded(d=1) { - projection(cut=true) minkowski() { - cylinder(r=d); - linear_extrude(center=true) children(); - } +module outset_extruded(d = 1) { + projection(cut=true) minkowski() { + cylinder(r=d); + linear_extrude(center=true) children(); + } } -module inset(d=1) { - render() inverse() outset(d=d) inverse() children(); +module inset(d = 1) { + render() inverse() outset(d=d) inverse() children(); } -module fillet(r=1) { - inset(d=r) render() outset(d=r) children(); +module fillet(r = 1) { + inset(d=r) render() outset(d=r) children(); } -module rounding(r=1) { - outset(d=r) inset(d=r) children(); +module rounding(r = 1) { + outset(d=r) inset(d=r) children(); } -module shell(d,center=false) { - if (center && d > 0) { - difference() { - outset(d=d/2) children(); - inset(d=d/2) children(); - } - } - if (!center && d > 0) { - difference() { - outset(d=d) children(); - children(); - } - } - if (!center && d < 0) { - difference() { - children(); - inset(d=-d) children(); - } - } - if (d == 0) children(); +module shell(d, center = false) { + if (center && d > 0) { + difference() { + outset(d=d / 2) children(); + inset(d=d / 2) children(); + } + } + if (!center && d > 0) { + difference() { + outset(d=d) children(); + children(); + } + } + if (!center && d < 0) { + difference() { + children(); + inset(d=-d) children(); + } + } + if (d == 0) children(); } - // Below are for internal use only module inverse() { - difference() { - square(1e5,center=true); - children(); - } + difference() { + square(1e5, center=true); + children(); + } } - // TEST CODE use -module arrow(l=1,w=.6,t=0.15) { - mirror_y() polygon([[0,0],[l,0],[l-w/2,w/2],[l-w/2-sqrt(2)*t,w/2],[l-t/2-sqrt(2)*t,t/2],[0,t/2]]); +module arrow(l = 1, w = .6, t = 0.15) { + mirror_y() polygon([[0, 0], [l, 0], [l - w / 2, w / 2], [l - w / 2 - sqrt(2) * t, w / 2], [l - t / 2 - sqrt(2) * t, t / 2], [0, t / 2]]); } module shape() { - polygon([[0,0],[1,0],[1.5,1],[2.5,1],[2,-1],[0,-1]]); + polygon([[0, 0], [1, 0], [1.5, 1], [2.5, 1], [2, -1], [0, -1]]); } -if(0) assign($fn=32) { - - for (p = [0:10*3-1]) assign(o=floor(p/3)) { - translate([(p%3)*2.5,-o*3]) { - //%if (p % 3 == 1) translate([0,0,1]) shape(); - if (p % 3 == 0) shape(); - if (p % 3 == 1) translate([0.6,0]) arrow(); - if (p % 3 == 2) { - if (o == 0) inset(d=0.3) shape(); - if (o == 1) outset(d=0.3) shape(); - if (o == 2) rounding(r=0.3) shape(); - if (o == 3) fillet(r=0.3) shape(); - if (o == 4) shell(d=0.3) shape(); - if (o == 5) shell(d=-0.3) shape(); - if (o == 6) shell(d=0.3,center=true) shape(); - if (o == 7) rounding(r=0.3) fillet(r=0.3) shape(); - if (o == 8) shell(d=0.3,center=true) fillet(r=0.3) rounding(r=0.3) shape(); - if (o == 9) shell(d=-0.3) fillet(r=0.3) rounding(r=0.3) shape(); - } - } - } -} +if (0) + assign ($fn = 32) { + + for (p = [0:10 * 3 - 1]) + assign (o = floor(p / 3)) { + translate([(p % 3) * 2.5, -o * 3]) { + //%if (p % 3 == 1) translate([0,0,1]) shape(); + if (p % 3 == 0) shape(); + if (p % 3 == 1) translate([0.6, 0]) arrow(); + if (p % 3 == 2) { + if (o == 0) inset(d=0.3) shape(); + if (o == 1) outset(d=0.3) shape(); + if (o == 2) rounding(r=0.3) shape(); + if (o == 3) fillet(r=0.3) shape(); + if (o == 4) shell(d=0.3) shape(); + if (o == 5) shell(d=-0.3) shape(); + if (o == 6) shell(d=0.3, center=true) shape(); + if (o == 7) rounding(r=0.3) fillet(r=0.3) shape(); + if (o == 8) shell(d=0.3, center=true) fillet(r=0.3) rounding(r=0.3) shape(); + if (o == 9) shell(d=-0.3) fillet(r=0.3) rounding(r=0.3) shape(); + } + } + } + } diff --git a/se3.scad b/se3.scad index 48b59e6..a153228 100644 --- a/se3.scad +++ b/se3.scad @@ -4,57 +4,71 @@ use function combine_se3_exp(w, ABt) = construct_Rt(rodrigues_so3_exp(w, ABt[0], ABt[1]), ABt[2]); // [A,B,t] -function se3_exp_1(t,w) = concat( - so3_exp_1(w*w), - [t + 0.5 * cross(w,t)] -); - -function se3_exp_2(t,w) = se3_exp_2_0(t,w,w*w); -function se3_exp_2_0(t,w,theta_sq) = -se3_exp_23( - so3_exp_2(theta_sq), - C = (1.0 - theta_sq/20) / 6, - t=t,w=w); - -function se3_exp_3(t,w) = se3_exp_3_0(t,w,sqrt(w*w)*180/PI,1/sqrt(w*w)); - -function se3_exp_3_0(t,w,theta_deg,inv_theta) = -se3_exp_23( - so3_exp_3_0(theta_deg = theta_deg, inv_theta = inv_theta), - C = (1 - sin(theta_deg) * inv_theta) * (inv_theta * inv_theta), - t=t,w=w); - -function se3_exp_23(AB,C,t,w) = -[AB[0], AB[1], t + AB[1] * cross(w,t) + C * cross(w,cross(w,t)) ]; - -function se3_exp(mu) = se3_exp_0(t=take3(mu),w=tail3(mu)/180*PI); - -function se3_exp_0(t,w) = -combine_se3_exp(w, -// Evaluate by Taylor expansion when near 0 - w*w < 1e-8 - ? se3_exp_1(t,w) - : w*w < 1e-6 - ? se3_exp_2(t,w) - : se3_exp_3(t,w) -); +function se3_exp_1(t, w) = + concat( + so3_exp_1(w * w), + [t + 0.5 * cross(w, t)] + ); + +function se3_exp_2(t, w) = se3_exp_2_0(t, w, w * w); +function se3_exp_2_0(t, w, theta_sq) = + se3_exp_23( + so3_exp_2(theta_sq), + C=(1.0 - theta_sq / 20) / 6, + t=t, w=w + ); + +function se3_exp_3(t, w) = se3_exp_3_0(t, w, sqrt(w * w) * 180 / PI, 1 / sqrt(w * w)); + +function se3_exp_3_0(t, w, theta_deg, inv_theta) = + se3_exp_23( + so3_exp_3_0(theta_deg=theta_deg, inv_theta=inv_theta), + C=(1 - sin(theta_deg) * inv_theta) * (inv_theta * inv_theta), + t=t, w=w + ); + +function se3_exp_23(AB, C, t, w) = + [AB[0], AB[1], t + AB[1] * cross(w, t) + C * cross(w, cross(w, t))]; + +function se3_exp(mu) = se3_exp_0(t=take3(mu), w=tail3(mu) / 180 * PI); + +function se3_exp_0(t, w) = + combine_se3_exp( + w, + // Evaluate by Taylor expansion when near 0 + w * w < 1e-8 ? se3_exp_1(t, w) + : w * w < 1e-6 ? se3_exp_2(t, w) + : se3_exp_3(t, w) + ); function se3_ln(m) = se3_ln_to_deg(se3_ln_rad(m)); -function se3_ln_to_deg(v) = concat(take3(v),tail3(v)*180/PI); - -function se3_ln_rad(m) = se3_ln_0(m, - rot = so3_ln_rad(rotation_part(m))); -function se3_ln_0(m,rot) = se3_ln_1(m,rot, - theta = sqrt(rot*rot)); -function se3_ln_1(m,rot,theta) = se3_ln_2(m,rot,theta, - shtot = theta > 0.00001 ? sin(theta/2*180/PI)/theta : 0.5, - halfrotator = so3_exp_rad(rot * -.5)); -function se3_ln_2(m,rot,theta,shtot,halfrotator) = -concat( (halfrotator * translation_part(m) - - (theta > 0.001 - ? rot * ((translation_part(m) * rot) * (1-2*shtot) / (rot*rot)) - : rot * ((translation_part(m) * rot)/24) - )) / (2 * shtot), rot); - -__se3_test = [20,-40,60,-80,100,-120]; -echo(UNITTEST_se3=norm(__se3_test-se3_ln(se3_exp(__se3_test))) < 1e-8); +function se3_ln_to_deg(v) = concat(take3(v), tail3(v) * 180 / PI); + +function se3_ln_rad(m) = + se3_ln_0( + m, + rot=so3_ln_rad(rotation_part(m)) + ); +function se3_ln_0(m, rot) = + se3_ln_1( + m, rot, + theta=sqrt(rot * rot) + ); +function se3_ln_1(m, rot, theta) = + se3_ln_2( + m, rot, theta, + shtot=theta > 0.00001 ? sin(theta / 2 * 180 / PI) / theta : 0.5, + halfrotator=so3_exp_rad(rot * -.5) + ); +function se3_ln_2(m, rot, theta, shtot, halfrotator) = + concat( + ( + halfrotator * translation_part(m) - ( + theta > 0.001 ? rot * ( (translation_part(m) * rot) * (1 - 2 * shtot) / (rot * rot)) + : rot * ( (translation_part(m) * rot) / 24) + ) + ) / (2 * shtot), rot + ); + +__se3_test = [20, -40, 60, -80, 100, -120]; +echo(UNITTEST_se3=norm(__se3_test - se3_ln(se3_exp(__se3_test))) < 1e-8); diff --git a/shapes.scad b/shapes.scad index 6d33d8a..fe6e260 100644 --- a/shapes.scad +++ b/shapes.scad @@ -1,16 +1,17 @@ -function square(size) = [[-size,-size], [-size,size], [size,size], [size,-size]] / 2; +function square(size) = [[-size, -size], [-size, size], [size, size], [size, -size]] / 2; -function circle(r) = [for (i=[0:$fn-1]) let (a=i*360/$fn) r * [cos(a), sin(a)]]; +function circle(r) = [for (i = [0:$fn - 1]) let (a = i * 360 / $fn) r * [cos(a), sin(a)]]; function regular(r, n) = circle(r, $fn=n); -function rectangle_profile(size=[1,1]) = [ - // The first point is the anchor point, put it on the point corresponding to [cos(0),sin(0)] - [ size[0]/2, 0], - [ size[0]/2, size[1]/2], - [-size[0]/2, size[1]/2], - [-size[0]/2, -size[1]/2], - [ size[0]/2, -size[1]/2], -]; +function rectangle_profile(size = [1, 1]) = + [ + // The first point is the anchor point, put it on the point corresponding to [cos(0),sin(0)] + [size[0] / 2, 0], + [size[0] / 2, size[1] / 2], + [-size[0] / 2, size[1] / 2], + [-size[0] / 2, -size[1] / 2], + [size[0] / 2, -size[1] / 2], + ]; -// FIXME: Move rectangle and rounded rectangle from extrusion \ No newline at end of file +// FIXME: Move rectangle and rounded rectangle from extrusion diff --git a/so3.scad b/so3.scad index 83308aa..66d55c3 100644 --- a/so3.scad +++ b/so3.scad @@ -2,81 +2,84 @@ use -function rodrigues_so3_exp(w, A, B) = [ -[1.0 - B*(w[1]*w[1] + w[2]*w[2]), B*(w[0]*w[1]) - A*w[2], B*(w[0]*w[2]) + A*w[1]], -[B*(w[0]*w[1]) + A*w[2], 1.0 - B*(w[0]*w[0] + w[2]*w[2]), B*(w[1]*w[2]) - A*w[0]], -[B*(w[0]*w[2]) - A*w[1], B*(w[1]*w[2]) + A*w[0], 1.0 - B*(w[0]*w[0] + w[1]*w[1])] -]; - -function so3_exp(w) = so3_exp_rad(w/180*PI); +function rodrigues_so3_exp(w, A, B) = + [ + [1.0 - B * (w[1] * w[1] + w[2] * w[2]), B * (w[0] * w[1]) - A * w[2], B * (w[0] * w[2]) + A * w[1]], + [B * (w[0] * w[1]) + A * w[2], 1.0 - B * (w[0] * w[0] + w[2] * w[2]), B * (w[1] * w[2]) - A * w[0]], + [B * (w[0] * w[2]) - A * w[1], B * (w[1] * w[2]) + A * w[0], 1.0 - B * (w[0] * w[0] + w[1] * w[1])], + ]; + +function so3_exp(w) = so3_exp_rad(w / 180 * PI); function so3_exp_rad(w) = -combine_so3_exp(w, - w*w < 1e-8 - ? so3_exp_1(w*w) - : w*w < 1e-6 - ? so3_exp_2(w*w) - : so3_exp_3(w*w)); + combine_so3_exp( + w, + w * w < 1e-8 ? so3_exp_1(w * w) + : w * w < 1e-6 ? so3_exp_2(w * w) + : so3_exp_3(w * w) + ); -function combine_so3_exp(w,AB) = rodrigues_so3_exp(w,AB[0],AB[1]); +function combine_so3_exp(w, AB) = rodrigues_so3_exp(w, AB[0], AB[1]); // Taylor series expansions close to 0 -function so3_exp_1(theta_sq) = [ - 1 - 1/6*theta_sq, - 0.5 -]; - -function so3_exp_2(theta_sq) = [ - 1.0 - theta_sq * (1.0 - theta_sq/20) / 6, - 0.5 - 0.25/6 * theta_sq -]; - -function so3_exp_3_0(theta_deg, inv_theta) = [ - sin(theta_deg) * inv_theta, - (1 - cos(theta_deg)) * (inv_theta * inv_theta) -]; - -function so3_exp_3(theta_sq) = so3_exp_3_0(sqrt(theta_sq)*180/PI, 1/sqrt(theta_sq)); - - -function rot_axis_part(m) = [m[2][1] - m[1][2], m[0][2] - m[2][0], m[1][0] - m[0][1]]*0.5; - -function so3_ln(m) = 180/PI*so3_ln_rad(m); -function so3_ln_rad(m) = so3_ln_0(m, - cos_angle = rot_cos_angle(m), - preliminary_result = rot_axis_part(m)); - -function so3_ln_0(m, cos_angle, preliminary_result) = -so3_ln_1(m, cos_angle, preliminary_result, - sin_angle_abs = sqrt(preliminary_result*preliminary_result)); - -function so3_ln_1(m, cos_angle, preliminary_result, sin_angle_abs) = - cos_angle > sqrt(1/2) - ? sin_angle_abs > 0 - ? preliminary_result * asin(sin_angle_abs)*PI/180 / sin_angle_abs - : preliminary_result - : cos_angle > -sqrt(1/2) - ? preliminary_result * acos(cos_angle)*PI/180 / sin_angle_abs - : so3_get_symmetric_part_rotation( - preliminary_result, - m, - angle = PI - asin(sin_angle_abs)*PI/180, - d0 = m[0][0] - cos_angle, - d1 = m[1][1] - cos_angle, - d2 = m[2][2] - cos_angle - ); +function so3_exp_1(theta_sq) = + [ + 1 - 1 / 6 * theta_sq, + 0.5, + ]; + +function so3_exp_2(theta_sq) = + [ + 1.0 - theta_sq * (1.0 - theta_sq / 20) / 6, + 0.5 - 0.25 / 6 * theta_sq, + ]; + +function so3_exp_3_0(theta_deg, inv_theta) = + [ + sin(theta_deg) * inv_theta, + (1 - cos(theta_deg)) * (inv_theta * inv_theta), + ]; + +function so3_exp_3(theta_sq) = so3_exp_3_0(sqrt(theta_sq) * 180 / PI, 1 / sqrt(theta_sq)); + +function rot_axis_part(m) = [m[2][1] - m[1][2], m[0][2] - m[2][0], m[1][0] - m[0][1]] * 0.5; + +function so3_ln(m) = 180 / PI * so3_ln_rad(m); +function so3_ln_rad(m) = + so3_ln_0( + m, + cos_angle=rot_cos_angle(m), + preliminary_result=rot_axis_part(m) + ); + +function so3_ln_0(m, cos_angle, preliminary_result) = + so3_ln_1( + m, cos_angle, preliminary_result, + sin_angle_abs=sqrt(preliminary_result * preliminary_result) + ); + +function so3_ln_1(m, cos_angle, preliminary_result, sin_angle_abs) = + cos_angle > sqrt(1 / 2) ? sin_angle_abs > 0 ? preliminary_result * asin(sin_angle_abs) * PI / 180 / sin_angle_abs + : preliminary_result + : cos_angle > -sqrt(1 / 2) ? preliminary_result * acos(cos_angle) * PI / 180 / sin_angle_abs + : so3_get_symmetric_part_rotation( + preliminary_result, + m, + angle=PI - asin(sin_angle_abs) * PI / 180, + d0=m[0][0] - cos_angle, + d1=m[1][1] - cos_angle, + d2=m[2][2] - cos_angle + ); function so3_get_symmetric_part_rotation(preliminary_result, m, angle, d0, d1, d2) = -so3_get_symmetric_part_rotation_0(preliminary_result,angle,so3_largest_column(m, d0, d1, d2)); + so3_get_symmetric_part_rotation_0(preliminary_result, angle, so3_largest_column(m, d0, d1, d2)); function so3_get_symmetric_part_rotation_0(preliminary_result, angle, c_max) = - angle * unit(c_max * preliminary_result < 0 ? -c_max : c_max); + angle * unit(c_max * preliminary_result < 0 ? -c_max : c_max); function so3_largest_column(m, d0, d1, d2) = - d0*d0 > d1*d1 && d0*d0 > d2*d2 - ? [d0, (m[1][0]+m[0][1])/2, (m[0][2]+m[2][0])/2] - : d1*d1 > d2*d2 - ? [(m[1][0]+m[0][1])/2, d1, (m[2][1]+m[1][2])/2] - : [(m[0][2]+m[2][0])/2, (m[2][1]+m[1][2])/2, d2]; - -__so3_test = [12,-125,110]; -echo(UNITTEST_so3=norm(__so3_test-so3_ln(so3_exp(__so3_test))) < 1e-8); + d0 * d0 > d1 * d1 && d0 * d0 > d2 * d2 ? [d0, (m[1][0] + m[0][1]) / 2, (m[0][2] + m[2][0]) / 2] + : d1 * d1 > d2 * d2 ? [(m[1][0] + m[0][1]) / 2, d1, (m[2][1] + m[1][2]) / 2] + : [(m[0][2] + m[2][0]) / 2, (m[2][1] + m[1][2]) / 2, d2]; + +__so3_test = [12, -125, 110]; +echo(UNITTEST_so3=norm(__so3_test - so3_ln(so3_exp(__so3_test))) < 1e-8); diff --git a/spline.scad b/spline.scad index dbf487b..384305c 100644 --- a/spline.scad +++ b/spline.scad @@ -11,45 +11,55 @@ use use -q1=[[1,0,0,0],[1,1,1,1],[0,1,2,3],[0,0,1,3]]; -q1inv=[[1,0,0,0],[-3,3,-2,1],[3,-3,3,-2],[-1,1,-1,1]]; -q2=[[0,0,0,0],[0,0,0,0],[0,-1,0,0],[0,0,-1,0]]; -qn1i2=-q1inv*q2; -z3=[0,0,0]; -z4=[0,0,0,0]; - -function matrix_power(m,n)= n==0? (len(m)==3?identity3():identity4()) : - n==1 ? m : (n%2==1) ? matrix_power(m*m,floor(n/2))*m : matrix_power(m*m,n/2); - -function det(m) = let(r=[for(i=[0:1:len(m)-1]) i]) det_help(m, 0, r); +q1 = [[1, 0, 0, 0], [1, 1, 1, 1], [0, 1, 2, 3], [0, 0, 1, 3]]; +q1inv = [[1, 0, 0, 0], [-3, 3, -2, 1], [3, -3, 3, -2], [-1, 1, -1, 1]]; +q2 = [[0, 0, 0, 0], [0, 0, 0, 0], [0, -1, 0, 0], [0, 0, -1, 0]]; +qn1i2 = -q1inv * q2; +z3 = [0, 0, 0]; +z4 = [0, 0, 0, 0]; + +function matrix_power(m, n) = + n == 0 ? (len(m) == 3 ? identity3() : identity4()) + : n == 1 ? m : (n % 2 == 1) ? matrix_power(m * m, floor(n / 2)) * m : matrix_power(m * m, n / 2); + +function det(m) = let (r = [for (i = [0:1:len(m) - 1]) i]) det_help(m, 0, r); // Construction indices list is inefficient, but currently there is no way to imperatively // assign to a list element -function det_help(m, i, r) = len(r) == 0 ? 1 : - m[len(m)-len(r)][r[i]]*det_help(m,0,remove(r,i)) - (i+1=n? u : u-q2*q1inv*spline_helper(i+1, n, p); +function spline_helper(i, n, p) = let (u = [p[i], p[i + 1], z3, z3]) i + 3 >= n ? u : u - q2 * q1inv * spline_helper(i + 1, n, p); // knowing s[j+1], calculate s[j]. Stop when found s[i] -function spline_si(i,n, p, sn) = i == n ? sn : q1inv*(spline_u(i,p)-q2*spline_si(i+1, n, p, sn)); +function spline_si(i, n, p, sn) = i == n ? sn : q1inv * (spline_u(i, p) - q2 * spline_si(i + 1, n, p, sn)); // Takes array of (3n+1) points or (2n + 2) points, if tangent segments are symmetric. // For non-symmetric version input is: point0, normal0, neg_normal1, point1, normal1, ... neg_normal_n, point_n @@ -58,56 +68,57 @@ function spline_si(i,n, p, sn) = i == n ? sn : q1inv*(spline_u(i,p)-q2*spline_si // I.e. if current points are p0,p1,p2 then anchor points are p0 and p2, first tangent defined by p1-p0, // second tangent defined by p3-p2. // Return array of coefficients accepted by spline(), spline_tan() and similar -function bezier3_args(p, symmetric=false) = let(step=symmetric?2:3) - [for(i=[0:step:len(p)-3]) [[1,0,0,0],[-3,3,0,0],[3,-6,3,0],[-1,3,-3,1]]* - (symmetric?[p[i],p[i]+p[i+1],p[i+2]-p[i+3],p[i+2]] : [p[i], p[i]+p[i+1], p[i+3]+p[i+2], p[i+3]])]; - +function bezier3_args(p, symmetric = false) = + let (step = symmetric ? 2 : 3) [ + for (i = [0:step:len(p) - 3]) [[1, 0, 0, 0], [-3, 3, 0, 0], [3, -6, 3, 0], [-1, 3, -3, 1]] * (symmetric ? [p[i], p[i] + p[i + 1], p[i + 2] - p[i + 3], p[i + 2]] : [p[i], p[i] + p[i + 1], p[i + 3] + p[i + 2], p[i + 3]]), + ]; + // s - spline arguments calculated by spline_args // t - defines point on curve. each segment length is 1. I.e. t= 0..1 is first segment, t=1..2 - second. -function spline(s, t)= let(i=t>=len(s)?len(s)-1: floor(t), t2=t-i) [1,t2,t2*t2,t2*t2*t2]*s[i]; - -function spline_tan(s, t)= let(i=t>=len(s)?len(s)-1: floor(t), t2=t-i) [0,1,2*t2,3*t2*t2]*s[i]; -function spline_tan_unit(s, t)= unit(spline_tan(s,t)); -function spline_d2(s,t)= let(i=t>=len(s)?len(s)-1: floor(t), t2=t-i) [0,0,2,6*t2]*s[i]; -function spline_binormal_unit(s,t)= unit(cross(spline_tan(s, t), spline_d2(s,t))); -function spline_normal_unit(s,t)= unit(cross(spline_tan(s, t), spline_binormal_unit(s,t))); - -function spline_transform(s, t)= - construct_Rt(transpose_3([spline_normal_unit(s,t), spline_binormal_unit(s,t), spline_tan_unit(s,t)]), spline(s,t)); - +function spline(s, t) = let (i = t >= len(s) ? len(s) - 1 : floor(t), t2 = t - i) [1, t2, t2 * t2, t2 * t2 * t2] * s[i]; + +function spline_tan(s, t) = let (i = t >= len(s) ? len(s) - 1 : floor(t), t2 = t - i) [0, 1, 2 * t2, 3 * t2 * t2] * s[i]; +function spline_tan_unit(s, t) = unit(spline_tan(s, t)); +function spline_d2(s, t) = let (i = t >= len(s) ? len(s) - 1 : floor(t), t2 = t - i) [0, 0, 2, 6 * t2] * s[i]; +function spline_binormal_unit(s, t) = unit(cross(spline_tan(s, t), spline_d2(s, t))); +function spline_normal_unit(s, t) = unit(cross(spline_tan(s, t), spline_binormal_unit(s, t))); + +function spline_transform(s, t) = + construct_Rt(transpose_3([spline_normal_unit(s, t), spline_binormal_unit(s, t), spline_tan_unit(s, t)]), spline(s, t)); + // Unit tests -__s = spline_args([[0,10,0], [10,0,0],[0,-5,2]], v1=[0,1,0], v2=[-1,0,0], closed=true); -for(t=[0:0.01:len(__s)]) translate(spline(__s, t)) - cube([0.2,0.2,0.2], center=true); +__s = spline_args([[0, 10, 0], [10, 0, 0], [0, -5, 2]], v1=[0, 1, 0], v2=[-1, 0, 0], closed=true); +for (t = [0:0.01:len(__s)]) + translate(spline(__s, t)) + cube([0.2, 0.2, 0.2], center=true); -__s1=spline_args([[0,0,0],[0,0,15], [26,0,26+15]], /*v1=[0,0,100],*/ v2=[40,0,0]); -for(t=[0:0.01:len(s1)]) translate(spline(__s1, t)) - cube([0.2,0.2,0.2], center=true); +__s1 = spline_args([[0, 0, 0], [0, 0, 15], [26, 0, 26 + 15]], /*v1=[0,0,100],*/ v2=[40, 0, 0]); +for (t = [0:0.01:len(s1)]) + translate(spline(__s1, t)) + cube([0.2, 0.2, 0.2], center=true); -__s2=bezier3_args([[0,0,0],[0,0,10],[0,0,15],[0,0,26*0.552284],[26,0,41],[26*0.552284,0,0]],symmetric=true); +__s2 = bezier3_args([[0, 0, 0], [0, 0, 10], [0, 0, 15], [0, 0, 26 * 0.552284], [26, 0, 41], [26 * 0.552284, 0, 0]], symmetric=true); echo(__s2); -for(t=[0:0.01:len(__s2)]) translate(spline(__s2, t)) - cube([0.2,0.2,0.2], center=true); +for (t = [0:0.01:len(__s2)]) + translate(spline(__s2, t)) + cube([0.2, 0.2, 0.2], center=true); // Rotation methods taken from list-comprehension-demos/sweep.scad to demonstrate normal and binormal // Normally spline_transform is more convenient -function __rotation_from_axis(x,y,z) = [[x[0],y[0],z[0]],[x[1],y[1],z[1]],[x[2],y[2],z[2]]]; -function __rotate_from_to(a,b,_axis=[]) = - len(_axis) == 0 - ? __rotate_from_to(a,b,unit(cross(a,b))) - : _axis*_axis >= 0.99 ? __rotation_from_axis(unit(b),_axis,cross(_axis,unit(b))) * - transpose_3(__rotation_from_axis(unit(a),_axis,cross(_axis,unit(a)))) : identity3(); - -__s3 = spline_args([[0,10,0], [6,6,0], [10,0,0],[0,-5,4]], v1=[0,1,0], v2=[-1,0,0], closed=true); -for(t=[0:0.05:len(__s3)]) translate(spline(__s3, t)) { - translate([0,0,3]) multmatrix(m=__rotate_from_to([0,0,1],spline_normal_unit(__s3,t))) +function __rotation_from_axis(x, y, z) = [[x[0], y[0], z[0]], [x[1], y[1], z[1]], [x[2], y[2], z[2]]]; +function __rotate_from_to(a, b, _axis = []) = + len(_axis) == 0 ? __rotate_from_to(a, b, unit(cross(a, b))) + : _axis * _axis >= 0.99 ? __rotation_from_axis(unit(b), _axis, cross(_axis, unit(b))) * transpose_3(__rotation_from_axis(unit(a), _axis, cross(_axis, unit(a)))) + : identity3(); + +__s3 = spline_args([[0, 10, 0], [6, 6, 0], [10, 0, 0], [0, -5, 4]], v1=[0, 1, 0], v2=[-1, 0, 0], closed=true); +for (t = [0:0.05:len(__s3)]) + translate(spline(__s3, t)) { + translate([0, 0, 3]) multmatrix(m=__rotate_from_to([0, 0, 1], spline_normal_unit(__s3, t))) cylinder(r1=0.1, r2=0, h=1, $fn=3); - translate([0,0,6]) multmatrix(m=__rotate_from_to([0,0,1],spline_binormal_unit(__s3,t))) + translate([0, 0, 6]) multmatrix(m=__rotate_from_to([0, 0, 1], spline_binormal_unit(__s3, t))) cylinder(r1=0.1, r2=0, h=1, $fn=3); -} - -translate([0,0,9]) for(t=[0:0.025:len(__s3)]) - multmatrix(spline_transform(__s3,t)) cube([1,1,0.1],center=true); + } - - +translate([0, 0, 9])for (t = [0:0.025:len(__s3)]) + multmatrix(spline_transform(__s3, t)) cube([1, 1, 0.1], center=true); diff --git a/trajectory.scad b/trajectory.scad index a7bcf81..af059f0 100644 --- a/trajectory.scad +++ b/trajectory.scad @@ -1,43 +1,48 @@ use -function val(a=undef,default=undef) = a == undef ? default : a; -function vec_is_undef(x,index_=0) = index_ >= len(x) ? true : -is_undef_or_oob(x[index_]) && vec_is_undef(x,index_+1); +function val(a = undef, default = undef) = a == undef ? default : a; +function vec_is_undef(x, index_ = 0) = + index_ >= len(x) ? true + : is_undef_or_oob(x[index_]) && vec_is_undef(x, index_ + 1); function is_undef_or_oob(x) = is_undef(x) ? true : is_list(x) ? vec_is_undef(x) : false; // Either a or b, but not both -function either(a,b,default=undef) = is_undef_or_oob(a) ? (is_undef_or_oob(b) ? default : b) : is_undef_or_oob(b) ? a : undef; - -function translationv(left=undef,right=undef,up=undef,down=undef,forward=undef,backward=undef,translation=undef) = -translationv_2( - x = either(up,is_undef(down) ? down : -down), - y = either(right,is_undef(left) ? left : -left), - z = either(forward,is_undef(backward) ? backward : -backward), - translation = translation); - -function translationv_2(x,y,z,translation) = - x == undef && y == undef && z == undef ? translation : - is_undef_or_oob(translation) ? [val(x,0),val(y,0),val(z,0)] - : undef; - -function rotationv(pitch=undef,yaw=undef,roll=undef,rotation=undef) = - rotation == undef ? [val(yaw,0),val(pitch,0),val(roll,0)] : - pitch == undef && yaw == undef && roll == undef ? rotation : - undef; +function either(a, b, default = undef) = is_undef_or_oob(a) ? (is_undef_or_oob(b) ? default : b) : is_undef_or_oob(b) ? a : undef; + +function translationv(left = undef, right = undef, up = undef, down = undef, forward = undef, backward = undef, translation = undef) = + translationv_2( + x=either(up, is_undef(down) ? down : -down), + y=either(right, is_undef(left) ? left : -left), + z=either(forward, is_undef(backward) ? backward : -backward), + translation=translation + ); + +function translationv_2(x, y, z, translation) = + x == undef && y == undef && z == undef ? translation + : is_undef_or_oob(translation) ? [val(x, 0), val(y, 0), val(z, 0)] + : undef; + +function rotationv(pitch = undef, yaw = undef, roll = undef, rotation = undef) = + rotation == undef ? [val(yaw, 0), val(pitch, 0), val(roll, 0)] + : pitch == undef && yaw == undef && roll == undef ? rotation + : undef; function trajectory( - left=undef, right=undef, - up=undef, down=undef, - forward=undef, backward=undef, - translation=undef, - - pitch=undef, - yaw=undef, - roll=undef, - rotation=undef -) = concat( - translationv(left=left,right=right,up=up,down=down,forward=forward,backward=backward,translation=translation), - rotationv(pitch=pitch,yaw=yaw,roll=roll,rotation=rotation) -); - -function rotationm(rotation=undef,pitch=undef,yaw=undef,roll=undef) = so3_exp(rotationv(rotation=rotation,pitch=pitch,yaw=yaw,roll=roll)); + left = undef, + right = undef, + up = undef, + down = undef, + forward = undef, + backward = undef, + translation = undef, + pitch = undef, + yaw = undef, + roll = undef, + rotation = undef +) = + concat( + translationv(left=left, right=right, up=up, down=down, forward=forward, backward=backward, translation=translation), + rotationv(pitch=pitch, yaw=yaw, roll=roll, rotation=rotation) + ); + +function rotationm(rotation = undef, pitch = undef, yaw = undef, roll = undef) = so3_exp(rotationv(rotation=rotation, pitch=pitch, yaw=yaw, roll=roll)); diff --git a/trajectory_path.scad b/trajectory_path.scad index f9fffe3..905e9a3 100644 --- a/trajectory_path.scad +++ b/trajectory_path.scad @@ -1,89 +1,103 @@ use use -function left_multiply(a,bs,i_=0) = i_ >= len(bs) ? [] : - concat([ - a * bs[i_] - ], left_multiply(a,bs,i_+1)); +function left_multiply(a, bs, i_ = 0) = + i_ >= len(bs) ? [] + : concat( + [ + a * bs[i_], + ], left_multiply(a, bs, i_ + 1) + ); +function right_multiply(as, b, i_ = 0) = + i_ >= len(as) ? [] + : concat( + [ + as[i_] * b, + ], right_multiply(as, b, i_ + 1) + ); -function right_multiply(as,b,i_=0) = i_ >= len(as) ? [] : - concat([ - as[i_] * b - ], right_multiply(as,b,i_+1)); +function quantize_trajectory(trajectory, step = undef, start_position = 0, steps = undef, i_ = 0, length_ = undef) = + length_ == undef ? quantize_trajectory( + trajectory=trajectory, + start_position=(step == undef ? norm(take3(trajectory)) / steps * start_position : start_position), + length_=norm(take3(trajectory)), + step=step, steps=steps, i_=i_ + ) + : (steps == undef ? start_position > length_ : i_ >= steps) ? [] + : concat( + [ + // if steps is defined, ignore start_position + se3_exp( + trajectory * ( + steps == undef ? start_position / length_ + : i_ / (steps > 1 ? steps - 1 : 1) + ) + ), + ], quantize_trajectory(trajectory=trajectory, step=step, start_position=(steps == undef ? start_position + step : start_position), steps=steps, i_=i_ + 1, length_=length_) + ); -function quantize_trajectory(trajectory,step=undef,start_position=0,steps=undef,i_=0,length_=undef) = - length_ == undef ? quantize_trajectory( - trajectory=trajectory, - start_position=(step==undef?norm(take3(trajectory))/steps*start_position:start_position), - length_=norm(take3(trajectory)), - step=step,steps=steps,i_=i_) : - (steps==undef?start_position > length_:i_>=steps) ? [] : - concat([ - // if steps is defined, ignore start_position - se3_exp(trajectory*(steps==undef ? start_position/length_ - : i_/(steps>1?steps-1:1))) - ], quantize_trajectory(trajectory=trajectory,step=step,start_position=(steps==undef?start_position+step:start_position),steps=steps,i_=i_+1,length_=length_)); +function close_trajectory_loop(trajectories) = concat(trajectories, [se3_ln(invert_rt(trajectories_end_position(trajectories)))]); -function close_trajectory_loop(trajectories) = concat(trajectories,[se3_ln(invert_rt(trajectories_end_position(trajectories)))]); +function quantize_trajectories(trajectories, step = undef, start_position = 0, steps = undef, loop = false, last_ = identity4(), i_ = 0, current_length_ = undef, j_ = 0) = + // due to quantization differences, the last step may be missed. In that case, add it: + loop == true ? quantize_trajectories( + trajectories=close_trajectory_loop(trajectories), + step=step, + start_position=start_position, + steps=steps, + loop=false, + last_=last_, + i_=i_, + current_length_=current_length_, + j_=j_ + ) + : i_ >= len(trajectories) ? (j_ < steps ? [last_] : []) + : current_length_ == undef ? + quantize_trajectories( + trajectories=trajectories, + step=(step == undef ? trajectories_length(trajectories) / steps : step), + start_position=(step == undef ? start_position * trajectories_length(trajectories) / steps : start_position), + steps=steps, + loop=loop, + last_=last_, + i_=i_, + current_length_=norm(take3(trajectories[i_])), + j_=j_ + ) + : concat( + left_multiply( + last_, quantize_trajectory( + trajectory=trajectories[i_], + start_position=start_position, + step=step + ) + ), + quantize_trajectories( + trajectories=trajectories, + step=step, + start_position=start_position > current_length_ ? start_position - current_length_ + : step - ( (current_length_ - start_position) % step), + steps=steps, + loop=loop, + last_=last_ * se3_exp(trajectories[i_]), + i_=i_ + 1, + current_length_=undef, + j_=j_ + len( -function quantize_trajectories(trajectories,step=undef,start_position=0,steps=undef,loop=false,last_=identity4(),i_=0,current_length_=undef,j_=0) = - // due to quantization differences, the last step may be missed. In that case, add it: - loop==true ? quantize_trajectories( - trajectories=close_trajectory_loop(trajectories), - step=step, - start_position = start_position, - steps=steps, - loop=false, - last_=last_, - i_=i_, - current_length_=current_length_, - j_=j_) : - i_ >= len(trajectories) ? (j_ < steps ? [last_] : []) : - current_length_ == undef ? - quantize_trajectories( - trajectories=trajectories, - step = (step == undef ? trajectories_length(trajectories) / steps : step), - start_position = (step == undef ? start_position * trajectories_length(trajectories) / steps : start_position), - steps=steps, - loop=loop, - last_=last_, - i_=i_, - current_length_=norm(take3(trajectories[i_])), - j_=j_) : - concat( - left_multiply(last_,quantize_trajectory( - trajectory=trajectories[i_], - start_position=start_position, - step=step)), - quantize_trajectories( - trajectories=trajectories, - step=step, - start_position = start_position > current_length_ - ? start_position - current_length_ - : step - ((current_length_-start_position) % step), - steps=steps, - loop=loop, - last_=last_ * se3_exp(trajectories[i_]), - i_=i_+1, - current_length_ = undef, - j_=j_+len( + quantize_trajectory( + trajectory=trajectories[i_], + start_position=start_position, + step=step + ) + ) + ) + ); - quantize_trajectory( - trajectory=trajectories[i_], - start_position=start_position, - step=step - - )) - )) -; - - -function trajectories_length(trajectories, i_=0) = i_ >= len(trajectories) ? 0 - : norm(take3(trajectories[i_])) + trajectories_length(trajectories,i_+1); - - -function trajectories_end_position(rt,i_=0,last_=identity4()) = - i_ >= len(rt) ? last_ : - trajectories_end_position(rt, i_+1, last_ * se3_exp(rt[i_])); +function trajectories_length(trajectories, i_ = 0) = + i_ >= len(trajectories) ? 0 + : norm(take3(trajectories[i_])) + trajectories_length(trajectories, i_ + 1); +function trajectories_end_position(rt, i_ = 0, last_ = identity4()) = + i_ >= len(rt) ? last_ + : trajectories_end_position(rt, i_ + 1, last_ * se3_exp(rt[i_])); diff --git a/transformations.scad b/transformations.scad index 89edff6..9659b60 100644 --- a/transformations.scad +++ b/transformations.scad @@ -8,36 +8,36 @@ use xyz = euler angles = rz * ry * rx axis = rotation_axis * rotation_angle */ -function rotation(xyz=undef, axis=undef) = - xyz != undef && axis != undef ? undef : - xyz == undef ? se3_exp([0,0,0,axis[0],axis[1],axis[2]]) : - len(xyz) == undef ? rotation(axis=[0,0,xyz]) : - (len(xyz) >= 3 ? rotation(axis=[0,0,xyz[2]]) : identity4()) * - (len(xyz) >= 2 ? rotation(axis=[0,xyz[1],0]) : identity4()) * - (len(xyz) >= 1 ? rotation(axis=[xyz[0],0,0]) : identity4()); +function rotation(xyz = undef, axis = undef) = + xyz != undef && axis != undef ? undef + : xyz == undef ? se3_exp([0, 0, 0, axis[0], axis[1], axis[2]]) + : len(xyz) == undef ? rotation(axis=[0, 0, xyz]) + : (len(xyz) >= 3 ? rotation(axis=[0, 0, xyz[2]]) : identity4()) * (len(xyz) >= 2 ? rotation(axis=[0, xyz[1], 0]) : identity4()) * (len(xyz) >= 1 ? rotation(axis=[xyz[0], 0, 0]) : identity4()); /*! Creates a scaling matrix */ -function scaling(v) = [ - [v[0],0,0,0], - [0,v[1],0,0], - [0,0,v[2],0], - [0,0,0,1], -]; +function scaling(v) = + [ + [v[0], 0, 0, 0], + [0, v[1], 0, 0], + [0, 0, v[2], 0], + [0, 0, 0, 1], + ]; /*! Creates a translation matrix */ -function translation(v) = [ - [1,0,0,v[0]], - [0,1,0,v[1]], - [0,0,1,v[2]], - [0,0,0,1], -]; +function translation(v) = + [ + [1, 0, 0, v[0]], + [0, 1, 0, v[1]], + [0, 0, 1, v[2]], + [0, 0, 0, 1], + ]; // Convert between cartesian and homogenous coordinates -function project(x) = subarray(x,end=len(x)-1) / x[len(x)-1]; +function project(x) = subarray(x, end=len(x) - 1) / x[len(x) - 1]; -function transform(m, list) = [for (p=list) project(m * vec4(p))]; -function to_3d(list) = [ for(v = list) vec3(v) ]; +function transform(m, list) = [for (p = list) project(m * vec4(p))]; +function to_3d(list) = [for (v = list) vec3(v)]; From 4b3836a1183c9316086fd7cfe7a8fca914ddf196 Mon Sep 17 00:00:00 2001 From: Cameron K Brooks Date: Fri, 26 Sep 2025 18:52:42 -0400 Subject: [PATCH 02/32] fix: hull tests, collinear still not working --- hull.scad | 45 +++++++++++++++++++++++++-------------------- 1 file changed, 25 insertions(+), 20 deletions(-) diff --git a/hull.scad b/hull.scad index 3323ac9..6143f40 100644 --- a/hull.scad +++ b/hull.scad @@ -295,34 +295,39 @@ translate([0, 50]) visualize_hull(20 * testpoints_circular); translate([0, -50]) visualize_hull(20 * testpoints_coplanar); // Collinear -translate([50, 50]) visualize_hull(20 * testpoints_collinear_2d); +//translate([50, 50]) visualize_hull(20 * testpoints_collinear_2d); // Collinear -translate([-50, 50]) visualize_hull(20 * testpoints_collinear_3d); +//translate([-50, 50]) visualize_hull(20 * testpoints_collinear_3d); // 3D points visualize_hull(testpoints3d); module visualize_hull(points) { - hull = hull(points); - - if (len(hull) > 0 && len(hull[0]) > 0) - %polyhedron(points=points, faces=hull); - else - %polyhedron(points=points, faces=[hull]); - - for (i = [0:len(points) - 1]) - assign (p = points[i], $fn = 16) { - translate(p) { - if (hull_contains_index(hull, i)) { - color("blue") sphere(1); - } else { - color("red") sphere(1); - } + // compute hull faces once + let (faces = hull(points)) { + + // faces is either [ [i,i,i], ... ] for 3D or [i,i,i,...] for 2D + // Use is_list() to branch without calling len() on a scalar. + if (len(faces) > 0 && is_list(faces[0])) + %polyhedron(points=points, faces=faces); + else + %polyhedron(points=points, faces=[faces]); + + // draw points; blue if on hull, red otherwise + for (i = [0:len(points) - 1]) { + let (p = points[i]) translate(p) { + // pass faces (not undefined) to the lookup + color(hull_contains_index(faces, i) ? "blue" : "red") + sphere(r=1, $fn=16); } } - - function hull_contains_index(hull, index) = - search(index, hull, 1, 0) || search(index, hull, 1, 1) || search(index, hull, 1, 2); + } } + +// search convenience: true if index appears anywhere in faces +function hull_contains_index(faces, index) = + search(index, faces, 1, 0) || // search depth 1 + search(index, faces, 1, 1) || // search depth 2 + search(index, faces, 1, 2); From e4c9ead8e91c50efc566211437edd74335b5edaa Mon Sep 17 00:00:00 2001 From: Cameron K Brooks Date: Fri, 26 Sep 2025 19:15:30 -0400 Subject: [PATCH 03/32] fix: hull tests, visualize_hull to properly display collinear in 2d/3d to avoid degen. polygons issue --- hull.scad | 47 ++++++++++++++++++++++++++++++++++------------- 1 file changed, 34 insertions(+), 13 deletions(-) diff --git a/hull.scad b/hull.scad index 6143f40..2cf28d7 100644 --- a/hull.scad +++ b/hull.scad @@ -295,33 +295,39 @@ translate([0, 50]) visualize_hull(20 * testpoints_circular); translate([0, -50]) visualize_hull(20 * testpoints_coplanar); // Collinear -//translate([50, 50]) visualize_hull(20 * testpoints_collinear_2d); +translate([50, 50]) visualize_hull(20 * testpoints_collinear_2d); // Collinear -//translate([-50, 50]) visualize_hull(20 * testpoints_collinear_3d); +translate([-50, 50]) visualize_hull(20 * testpoints_collinear_3d); // 3D points visualize_hull(testpoints3d); - module visualize_hull(points) { - - // compute hull faces once let (faces = hull(points)) { - // faces is either [ [i,i,i], ... ] for 3D or [i,i,i,...] for 2D - // Use is_list() to branch without calling len() on a scalar. - if (len(faces) > 0 && is_list(faces[0])) + if (len(faces) == 0) { + // nothing + } else if (is_list(faces[0])) { + // 3D hull (triangles) %polyhedron(points=points, faces=faces); - else + } else if (len(faces) >= 3) { + // 2D hull (polygon) %polyhedron(points=points, faces=[faces]); + } else if (len(faces) == 2) { + // Collinear case, draw a rod + p0 = points[faces[0]]; + p1 = points[faces[1]]; + hull() { + translate(p0) sphere(r=0.8, $fn=24); + translate(p1) sphere(r=0.8, $fn=24); + } + } - // draw points; blue if on hull, red otherwise + // Draw all points for (i = [0:len(points) - 1]) { - let (p = points[i]) translate(p) { - // pass faces (not undefined) to the lookup + translate(points[i]) color(hull_contains_index(faces, i) ? "blue" : "red") sphere(r=1, $fn=16); - } } } } @@ -331,3 +337,18 @@ function hull_contains_index(faces, index) = search(index, faces, 1, 0) || // search depth 1 search(index, faces, 1, 1) || // search depth 2 search(index, faces, 1, 2); + +// Collinear case: draw a rod between two hull points +module visualize_segment_hull(points, indices) { + p0 = points[indices[0]]; + p1 = points[indices[1]]; + hull() { + translate(p0) sphere(r=0.8, $fn=24); + translate(p1) sphere(r=0.8, $fn=24); + } +} + +// Degenerate case: only one hull point +module visualize_point_hull(points, idx) { + translate(points[idx]) color("gray") sphere(r=0.8, $fn=24); +} From 48eb4e645e535ec2e95e69cfb9eae86ee2f902b7 Mon Sep 17 00:00:00 2001 From: Cameron K Brooks Date: Fri, 26 Sep 2025 19:43:11 -0400 Subject: [PATCH 04/32] refactor: cosmetic rework of hull and split off tests to dedicated dir --- hull.scad | 379 ++++++++++++------------------------------- tests/test_hull.scad | 78 +++++++++ 2 files changed, 180 insertions(+), 277 deletions(-) create mode 100644 tests/test_hull.scad diff --git a/hull.scad b/hull.scad index 2cf28d7..32d0606 100644 --- a/hull.scad +++ b/hull.scad @@ -1,354 +1,179 @@ +// ============================================================================ +// Convex Hull Utilities (2D and 3D) +// ---------------------------------------------------------------------------- +// This implementation calculates convex hulls for 2D or 3D point sets. +// - 2D: returns polygon vertex indices [i1, i2, i3, ...] +// - 3D: returns triangular face indices [[i1,i2,i3], [i2,i3,i4], ...] +// - Collinear: returns the two extreme endpoints [i_min, i_max] +// ---------------------------------------------------------------------------- +// Notes: +// * Uses let() and list comprehensions +// * Relies on bug-fixed search() +// * Assumes vector min()/max() are available +// ============================================================================ -// NOTE: this code uses -// * experimental let() syntax -// * experimental list comprehension syntax -// * search() bugfix and feature addition -// * vector min()/max() - -// Calculates the convex hull of a set of points. -// The result is expressed in point indices. -// If the points are collinear (or 2d), the result is a convex -// polygon [i1,i2,i3,...], otherwise a triangular -// polyhedron [[i1,i2,i3],[i2,i3,i4],...] +epsilon = 1e-9; +// --- Entry point ------------------------------------------------------------ function hull(points) = - !(len(points) > 0) ? [] - : len(points[0]) == 2 ? convexhull2d(points) - : len(points[0]) == 3 ? convexhull3d(points) : []; - -epsilon = 1e-9; + (len(points) == 0) ? [] + : (len(points[0]) == 2) ? convexhull2d(points) + : (len(points[0]) == 3) ? convexhull3d(points) : []; -// 2d version +// --- 2D Convex Hull --------------------------------------------------------- function convexhull2d(points) = - len(points) < 3 ? [] + (len(points) < 3) ? [] : let ( a = 0, b = 1, c = find_first_noncollinear([a, b], points, 2) - ) c == len(points) ? convexhull_collinear(points) + ) (c == len(points)) ? convexhull_collinear(points) : let ( remaining = [for (i = [2:len(points) - 1]) if (i != c) i], - polygon = area_2d(points[a], points[b], points[c]) > 0 ? [a, b, c] : [b, a, c] + polygon = (area_2d(points[a], points[b], points[c]) > 0) ? + [a, b, c] + : [b, a, c] ) convex_hull_iterative_2d(points, polygon, remaining); -// Adds the remaining points one by one to the convex hull function convex_hull_iterative_2d(points, polygon, remaining, i_ = 0) = - i_ >= len(remaining) ? polygon + (i_ >= len(remaining)) ? polygon : let ( - // pick a point i = remaining[i_], - - // find the segments that are in conflict with the point (point not inside) conflicts = find_conflicting_segments(points, polygon, points[i]) - - // no conflicts, skip point and move on - ) len(conflicts) == 0 ? convex_hull_iterative_2d(points, polygon, remaining, i_ + 1) + ) (len(conflicts) == 0) ? + convex_hull_iterative_2d(points, polygon, remaining, i_ + 1) : let ( - - // find the first conflicting segment and the first not conflicting - // conflict will be sorted, if not wrapping around, do it the easy way polygon = remove_conflicts_and_insert_point(polygon, conflicts, i) - ) convex_hull_iterative_2d( - points, - polygon, - remaining, - i_ + 1 - ); + ) convex_hull_iterative_2d(points, polygon, remaining, i_ + 1); function find_conflicting_segments(points, polygon, point) = [ for (i = [0:len(polygon) - 1]) let (j = (i + 1) % len(polygon)) if (area_2d(points[polygon[i]], points[polygon[j]], point) < 0) i, ]; -// remove the conflicting segments from the polygon function remove_conflicts_and_insert_point(polygon, conflicts, point) = - conflicts[0] == 0 ? let ( - nonconflicting = [for (i = [0:len(polygon) - 1]) if (!contains(conflicts, i)) i], - new_indices = concat(nonconflicting, (nonconflicting[len(nonconflicting) - 1] + 1) % len(polygon)), - polygon = concat([for (i = new_indices) polygon[i]], point) - ) polygon + (conflicts[0] == 0) ? + let ( + nonconf = [for (i = [0:len(polygon) - 1]) if (!contains(conflicts, i)) i], + indices = concat(nonconf, (nonconf[len(nonconf) - 1] + 1) % len(polygon)) + ) concat([for (i = indices) polygon[i]], point) : let ( - prior_to_first_conflict = [for (i = [0:1:min(conflicts)]) polygon[i]], - after_last_conflict = [for (i = [max(conflicts) + 1:1:len(polygon) - 1]) polygon[i]], - polygon = concat(prior_to_first_conflict, point, after_last_conflict) - ) polygon; + before = [for (i = [0:1:min(conflicts)]) polygon[i]], + after = [for (i = [max(conflicts) + 1:1:len(polygon) - 1]) polygon[i]] + ) concat(before, point, after); -// 3d version +// --- 3D Convex Hull --------------------------------------------------------- function convexhull3d(points) = - len(points) < 3 ? [for (i = [0:1:len(points) - 1]) i] + (len(points) < 3) ? [for (i = [0:1:len(points) - 1]) i] : let ( - - // start with a single triangle a = 0, b = 1, c = 2, - plane = plane(points, a, b, c), - d = find_first_noncoplanar(plane, points, 3) - ) d == len(points) ? /* all coplanar*/ let ( - pts2d = [for (p = points) plane_project(p, points[a], points[b], points[c])], - hull2d = convexhull2d(pts2d) - ) hull2d + pl = plane(points, a, b, c), + d = find_first_noncoplanar(pl, points, 3) + ) (d == len(points)) ? + // all coplanar: project to 2D hull + let ( + pts2d = [for (p = points) plane_project(p, points[a], points[b], points[c])] + ) convexhull2d(pts2d) : let ( remaining = [for (i = [3:len(points) - 1]) if (i != d) i], + bc = in_front(pl, points[d]) ? [c, b] : [b, c], + b2 = bc[0], + c2 = bc[1], + tris = [[a, b2, c2], [d, b2, a], [c2, d, a], [b2, d, c2]], + planes = [for (t = tris) plane(points, t[0], t[1], t[2])] + ) convex_hull_iterative(points, tris, planes, remaining); - // Build an initial tetrahedron - - // swap b,c if d is in front of triangle t - bc = in_front(plane, points[d]) ? [c, b] : [b, c], - b = bc[0], - c = bc[1], - triangles = [ - [a, b, c], - [d, b, a], - [c, d, a], - [b, d, c], - ], - - // calculate the plane equations - planes = [for (t = triangles) plane(points, t[0], t[1], t[2])] - ) convex_hull_iterative(points, triangles, planes, remaining); - -// A plane equation (normal, offset) function plane(points, a, b, c) = - let ( - normal = unit(cross(points[c] - points[a], points[b] - points[a])) - ) [ - normal, - normal * points[a], - ]; + let (n = unit(cross(points[c] - points[a], points[b] - points[a]))) [n, n * points[a]]; -// Adds the remaining points one by one to the convex hull function convex_hull_iterative(points, triangles, planes, remaining, i_ = 0) = - i_ >= len(remaining) ? triangles + (i_ >= len(remaining)) ? triangles : let ( - // pick a point - i = remaining[i_], - - // find the triangles that are in conflict with the point (point not inside) - conflicts = find_conflicts(points[i], planes), - - // for all triangles that are in conflict, collect their halfedges + idx = remaining[i_], + conflicts = find_conflicts(points[idx], planes), halfedges = [ - for (c = conflicts) for (i = [0:2]) let (j = (i + 1) % 3) [triangles[c][i], triangles[c][j]], + for (c = conflicts) for (k = [0:2]) let (j = (k + 1) % 3) [triangles[c][k], triangles[c][j]], ], - - // find the outer perimeter of the set of conflicting triangles horizon = remove_internal_edges(halfedges), - - // generate a new triangle for each horizon halfedge together with the picked point i - new_triangles = [for (h = horizon) concat(h, i)], - - // calculate the corresponding plane equations - new_planes = [for (t = new_triangles) plane(points, t[0], t[1], t[2])] + new_tris = [for (h = horizon) concat(h, idx)], + new_pls = [for (t = new_tris) plane(points, t[0], t[1], t[2])] ) convex_hull_iterative( points, - // remove the conflicting triangles and add the new ones - concat(remove_elements(triangles, conflicts), new_triangles), - concat(remove_elements(planes, conflicts), new_planes), + concat(remove_elements(triangles, conflicts), new_tris), + concat(remove_elements(planes, conflicts), new_pls), remaining, i_ + 1 ); +// --- Collinear Special Case ------------------------------------------------- function convexhull_collinear(points) = let ( n = points[1] - points[0], a = points[0], - points1d = [for (p = points) (p - a) * n], - min_i = min_index(points1d), - max_i = max_index(points1d) + pts1d = [for (p = points) (p - a) * n], + min_i = min_index(pts1d), + max_i = max_index(pts1d) ) [min_i, max_i]; -function min_index(values, min_, min_i_, i_) = - i_ == undef ? min_index(values, values[0], 0, 1) - : i_ >= len(values) ? min_i_ - : values[i_] < min_ ? min_index(values, values[i_], i_, i_ + 1) - : min_index(values, min_, min_i_, i_ + 1); +// --- Utility Functions ------------------------------------------------------ +function min_index(values, min_ = undef, idx_min = undef, i_ = 0) = + (i_ == 0) ? min_index(values, values[0], 0, 1) + : (i_ >= len(values)) ? idx_min + : (values[i_] < min_) ? min_index(values, values[i_], i_, i_ + 1) + : min_index(values, min_, idx_min, i_ + 1); -function max_index(values, max_, max_i_, i_) = - i_ == undef ? max_index(values, values[0], 0, 1) - : i_ >= len(values) ? max_i_ - : values[i_] > max_ ? max_index(values, values[i_], i_, i_ + 1) - : max_index(values, max_, max_i_, i_ + 1); +function max_index(values, max_ = undef, idx_max = undef, i_ = 0) = + (i_ == 0) ? max_index(values, values[0], 0, 1) + : (i_ >= len(values)) ? idx_max + : (values[i_] > max_) ? max_index(values, values[i_], i_, i_ + 1) + : max_index(values, max_, idx_max, i_ + 1); -function remove_elements(array, elements) = - [ - for (i = [0:len(array) - 1]) if (!search(i, elements)) array[i], - ]; +function remove_elements(arr, to_remove) = + [for (i = [0:len(arr) - 1]) if (!search(i, to_remove)) arr[i]]; -function remove_internal_edges(halfedges) = - [ - for (h = halfedges) if (!contains(halfedges, reverse(h))) h, - ]; +function remove_internal_edges(edges) = + [for (h = edges) if (!contains(edges, reverse(h))) h]; -function plane_project(point, a, b, c) = - let ( - u = b - a, - v = c - a, - n = cross(u, v), - w = cross(n, u), - relpoint = point - a - ) [relpoint * u, relpoint * w]; +function plane_project(p, a, b, c) = + let (u = b - a, v = c - a, n = cross(u, v), w = cross(n, u), rel = p - a) [rel * u, rel * w]; -function plane_unproject(point, a, b, c) = - let ( - u = b - a, - v = c - a, - n = cross(u, v), - w = cross(n, u) - ) a + point[0] * u + point[1] * w; +function plane_unproject(p, a, b, c) = + let (u = b - a, v = c - a, n = cross(u, v), w = cross(n, u)) a + p[0] * u + p[1] * w; function reverse(arr) = [for (i = [len(arr) - 1:-1:0]) arr[i]]; -function contains(arr, element) = search([element], arr)[0] != [] ? true : false; +function contains(arr, el) = (search([el], arr)[0] != []) ? true : false; -function find_conflicts(point, planes) = - [ - for (i = [0:len(planes) - 1]) if (in_front(planes[i], point)) i, - ]; +function find_conflicts(p, planes) = + [for (i = [0:len(planes) - 1]) if (in_front(planes[i], p)) i]; -function find_first_noncollinear(line, points, i) = - i >= len(points) ? len(points) - : collinear( - points[line[0]], - points[line[1]], - points[i] - ) ? find_first_noncollinear(line, points, i + 1) +function find_first_noncollinear(line, pts, i) = + (i >= len(pts)) ? len(pts) + : collinear(pts[line[0]], pts[line[1]], pts[i]) ? + find_first_noncollinear(line, pts, i + 1) : i; -function find_first_noncoplanar(plane, points, i) = - i >= len(points) ? len(points) - : coplanar(plane, points[i]) ? find_first_noncoplanar(plane, points, i + 1) - : i; +function find_first_noncoplanar(pl, pts, i) = + (i >= len(pts)) ? len(pts) + : coplanar(pl, pts[i]) ? find_first_noncoplanar(pl, pts, i + 1) : i; -function distance(plane, point) = plane[0] * point - plane[1]; - -function in_front(plane, point) = distance(plane, point) > epsilon; - -function coplanar(plane, point) = abs(distance(plane, point)) <= epsilon; +function distance(pl, p) = pl[0] * p - pl[1]; +function in_front(pl, p) = distance(pl, p) > epsilon; +function coplanar(pl, p) = abs(distance(pl, p)) <= epsilon; function unit(v) = v / norm(v); function area_2d(a, b, c) = - ( - a[0] * (b[1] - c[1]) + b[0] * (c[1] - a[1]) + c[0] * (a[1] - b[1]) - ) / 2; + (a[0] * (b[1] - c[1]) + b[0] * (c[1] - a[1]) + c[0] * (a[1] - b[1])) / 2; function collinear(a, b, c) = abs(area_2d(a, b, c)) < epsilon; -function spherical(cartesian) = - [ - atan2(cartesian[1], cartesian[0]), - asin(cartesian[2]), - ]; - -function cartesian(spherical) = +function spherical(cart) = [atan2(cart[1], cart[0]), asin(cart[2])]; +function cartesian(sph) = [ - cos(spherical[1]) * cos(spherical[0]), - cos(spherical[1]) * sin(spherical[0]), - sin(spherical[1]), + cos(sph[1]) * cos(sph[0]), + cos(sph[1]) * sin(sph[0]), + sin(sph[1]), ]; - -/// TESTCODE - -phi = 1.618033988749895; - -testpoints_on_sphere = [ - for ( - p = [ - [1, phi, 0], - [-1, phi, 0], - [1, -phi, 0], - [-1, -phi, 0], - [0, 1, phi], - [0, -1, phi], - [0, 1, -phi], - [0, -1, -phi], - [phi, 0, 1], - [-phi, 0, 1], - [phi, 0, -1], - [-phi, 0, -1], - ] - ) unit(p), -]; - -testpoints_spherical = [for (p = testpoints_on_sphere) spherical(p)]; -testpoints_circular = [for (a = [0:15:360 - epsilon]) [cos(a), sin(a)]]; - -testpoints_coplanar = let (u = unit([1, 3, 7]), v = unit([-2, 1, -2])) [for (i = [1:10]) rands(-1, 1, 1)[0] * u + rands(-1, 1, 1)[0] * v]; - -testpoints_collinear_2d = let (u = unit([5, 3])) [for (i = [1:20]) rands(-1, 1, 1)[0] * u]; -testpoints_collinear_3d = let (u = unit([5, 3, -5])) [for (i = [1:20]) rands(-1, 1, 1)[0] * u]; - -testpoints2d = 20 * [for (i = [1:10]) concat(rands(-1, 1, 2))]; -testpoints3d = 20 * [for (i = [1:50]) concat(rands(-1, 1, 3))]; - -// All points are on the sphere, no point should be red -translate([-50, 0]) visualize_hull(20 * testpoints_on_sphere); - -// 2D points -translate([50, 0]) visualize_hull(testpoints2d); - -// All points on a circle, no point should be red -translate([0, 50]) visualize_hull(20 * testpoints_circular); - -// All points 3d but collinear -translate([0, -50]) visualize_hull(20 * testpoints_coplanar); - -// Collinear -translate([50, 50]) visualize_hull(20 * testpoints_collinear_2d); - -// Collinear -translate([-50, 50]) visualize_hull(20 * testpoints_collinear_3d); - -// 3D points -visualize_hull(testpoints3d); -module visualize_hull(points) { - let (faces = hull(points)) { - - if (len(faces) == 0) { - // nothing - } else if (is_list(faces[0])) { - // 3D hull (triangles) - %polyhedron(points=points, faces=faces); - } else if (len(faces) >= 3) { - // 2D hull (polygon) - %polyhedron(points=points, faces=[faces]); - } else if (len(faces) == 2) { - // Collinear case, draw a rod - p0 = points[faces[0]]; - p1 = points[faces[1]]; - hull() { - translate(p0) sphere(r=0.8, $fn=24); - translate(p1) sphere(r=0.8, $fn=24); - } - } - - // Draw all points - for (i = [0:len(points) - 1]) { - translate(points[i]) - color(hull_contains_index(faces, i) ? "blue" : "red") - sphere(r=1, $fn=16); - } - } -} - -// search convenience: true if index appears anywhere in faces -function hull_contains_index(faces, index) = - search(index, faces, 1, 0) || // search depth 1 - search(index, faces, 1, 1) || // search depth 2 - search(index, faces, 1, 2); - -// Collinear case: draw a rod between two hull points -module visualize_segment_hull(points, indices) { - p0 = points[indices[0]]; - p1 = points[indices[1]]; - hull() { - translate(p0) sphere(r=0.8, $fn=24); - translate(p1) sphere(r=0.8, $fn=24); - } -} - -// Degenerate case: only one hull point -module visualize_point_hull(points, idx) { - translate(points[idx]) color("gray") sphere(r=0.8, $fn=24); -} diff --git a/tests/test_hull.scad b/tests/test_hull.scad new file mode 100644 index 0000000..e24bd7a --- /dev/null +++ b/tests/test_hull.scad @@ -0,0 +1,78 @@ +include <../hull.scad> + +// --- Test Data -------------------------------------------------------------- +phi = 1.618033988749895; + +testpoints_on_sphere = [ + for ( + p = [ + [1, phi, 0], + [-1, phi, 0], + [1, -phi, 0], + [-1, -phi, 0], + [0, 1, phi], + [0, -1, phi], + [0, 1, -phi], + [0, -1, -phi], + [phi, 0, 1], + [-phi, 0, 1], + [phi, 0, -1], + [-phi, 0, -1], + ] + ) unit(p), +]; +testpoints_spherical = [for (p = testpoints_on_sphere) spherical(p)]; +testpoints_circular = [for (a = [0:15:360 - epsilon]) [cos(a), sin(a)]]; +testpoints_coplanar = let (u = unit([1, 3, 7]), v = unit([-2, 1, -2])) [for (i = [1:10]) rands(-1, 1, 1)[0] * u + rands(-1, 1, 1)[0] * v]; +testpoints_collinear_2d = let (u = unit([5, 3])) [for (i = [1:20]) rands(-1, 1, 1)[0] * u]; +testpoints_collinear_3d = let (u = unit([5, 3, -5])) [for (i = [1:20]) rands(-1, 1, 1)[0] * u]; +testpoints2d = 20 * [for (i = [1:10]) concat(rands(-1, 1, 2))]; +testpoints3d = 20 * [for (i = [1:50]) concat(rands(-1, 1, 3))]; + +// --- Visualization ---------------------------------------------------------- +echo("Test points 3d"); +visualize_hull(testpoints3d); +echo("Test points on sphere"); +translate([-50, 0]) visualize_hull(20 * testpoints_on_sphere); +echo("Test points on 2d"); +translate([50, 0]) visualize_hull(testpoints2d); +echo("Test points on circular"); +translate([0, 50]) visualize_hull(20 * testpoints_circular); +echo("Test points coplanar"); +translate([0, -50]) visualize_hull(20 * testpoints_coplanar); +echo("Test points collinear 2d"); +translate([50, 50]) visualize_hull(20 * testpoints_collinear_2d); +echo("Test points collinear 3d"); +translate([-50, 50]) visualize_hull(20 * testpoints_collinear_3d); + +// --- Visualization Helpers -------------------------------------------------- +module visualize_hull(points) { + let (faces = hull(points)) { + echo("Faces: ", faces); + if (len(faces) == 0) { + // nothing + } else if (is_list(faces[0])) { + %polyhedron(points=points, faces=faces); // 3D hull + } else if (len(faces) >= 3) { + %polyhedron(points=points, faces=[faces]); // 2D hull + } else if (len(faces) == 2) { + // collinear case: show a rod between endpoints + p0 = points[faces[0]]; + p1 = points[faces[1]]; + hull() { + translate(p0) sphere(r=0.8, $fn=24); + translate(p1) sphere(r=0.8, $fn=24); + } + } + + // Draw points: blue if on hull, red otherwise + for (i = [0:len(points) - 1]) { + translate(points[i]) + color(hull_contains_index(faces, i) ? "blue" : "red") + sphere(r=1, $fn=16); + } + } +} + +function hull_contains_index(faces, idx) = + search(idx, faces, 1, 0) || search(idx, faces, 1, 1) || search(idx, faces, 1, 2); From 8440ccf40ffe4c2efb43f15b539b05aa93fcf995 Mon Sep 17 00:00:00 2001 From: Cameron K Brooks Date: Fri, 26 Sep 2025 19:45:08 -0400 Subject: [PATCH 05/32] fix(hull): handle 3D collinear points correctly; some formatting --- hull.scad | 79 ++++++++++++++++++++++++++++++++++++++++--------------- 1 file changed, 58 insertions(+), 21 deletions(-) diff --git a/hull.scad b/hull.scad index 32d0606..2be5f33 100644 --- a/hull.scad +++ b/hull.scad @@ -1,11 +1,12 @@ // ============================================================================ // Convex Hull Utilities (2D and 3D) // ---------------------------------------------------------------------------- -// This implementation calculates convex hulls for 2D or 3D point sets. +// Computes convex hulls for 2D or 3D point sets. +// // - 2D: returns polygon vertex indices [i1, i2, i3, ...] // - 3D: returns triangular face indices [[i1,i2,i3], [i2,i3,i4], ...] // - Collinear: returns the two extreme endpoints [i_min, i_max] -// ---------------------------------------------------------------------------- +// // Notes: // * Uses let() and list comprehensions // * Relies on bug-fixed search() @@ -14,7 +15,7 @@ epsilon = 1e-9; -// --- Entry point ------------------------------------------------------------ +// --- Entry Point ------------------------------------------------------------ function hull(points) = (len(points) == 0) ? [] : (len(points[0]) == 2) ? convexhull2d(points) @@ -30,9 +31,7 @@ function convexhull2d(points) = ) (c == len(points)) ? convexhull_collinear(points) : let ( remaining = [for (i = [2:len(points) - 1]) if (i != c) i], - polygon = (area_2d(points[a], points[b], points[c]) > 0) ? - [a, b, c] - : [b, a, c] + polygon = (area_2d(points[a], points[b], points[c]) > 0) ? [a, b, c] : [b, a, c] ) convex_hull_iterative_2d(points, polygon, remaining); function convex_hull_iterative_2d(points, polygon, remaining, i_ = 0) = @@ -65,19 +64,20 @@ function remove_conflicts_and_insert_point(polygon, conflicts, point) = // --- 3D Convex Hull --------------------------------------------------------- function convexhull3d(points) = (len(points) < 3) ? [for (i = [0:1:len(points) - 1]) i] + : is_collinear3(points) ? + convexhull_collinear(points) : let ( a = 0, b = 1, - c = 2, + c = find_first_noncollinear3([a, b], points, 2), pl = plane(points, a, b, c), - d = find_first_noncoplanar(pl, points, 3) + d = find_first_noncoplanar(pl, points, 0) ) (d == len(points)) ? - // all coplanar: project to 2D hull - let ( - pts2d = [for (p = points) plane_project(p, points[a], points[b], points[c])] - ) convexhull2d(pts2d) + // all coplanar → project to 2D + let (pts2d = [for (p = points) plane_project(p, points[a], points[b], points[c])]) convexhull2d(pts2d) : let ( - remaining = [for (i = [3:len(points) - 1]) if (i != d) i], + // build initial tetrahedron + remaining = [for (i = [0:len(points) - 1]) if (i != a && i != b && i != c && i != d) i], bc = in_front(pl, points[d]) ? [c, b] : [b, c], b2 = bc[0], c2 = bc[1], @@ -85,6 +85,23 @@ function convexhull3d(points) = planes = [for (t = tris) plane(points, t[0], t[1], t[2])] ) convex_hull_iterative(points, tris, planes, remaining); +// --- Helpers for 3D Collinearity ------------------------------------------- +function is_collinear3(pts) = + (len(pts) < 3) ? true + : let (a = pts[0], dir = pts[1] - a) _all_col3(pts, a, dir, 2); + +function _all_col3(pts, a, dir, i) = + (i >= len(pts)) ? true + : (norm(cross(pts[i] - a, dir)) <= epsilon) ? + _all_col3(pts, a, dir, i + 1) + : false; + +function find_first_noncollinear3(line, pts, i) = + (i >= len(pts)) ? len(pts) + : (norm(cross(pts[line[1]] - pts[line[0]], pts[i] - pts[line[0]])) <= epsilon) ? + find_first_noncollinear3(line, pts, i + 1) + : i; + function plane(points, a, b, c) = let (n = unit(cross(points[c] - points[a], points[b] - points[a]))) [n, n * points[a]]; @@ -121,13 +138,15 @@ function convexhull_collinear(points) = function min_index(values, min_ = undef, idx_min = undef, i_ = 0) = (i_ == 0) ? min_index(values, values[0], 0, 1) : (i_ >= len(values)) ? idx_min - : (values[i_] < min_) ? min_index(values, values[i_], i_, i_ + 1) + : (values[i_] < min_) ? + min_index(values, values[i_], i_, i_ + 1) : min_index(values, min_, idx_min, i_ + 1); function max_index(values, max_ = undef, idx_max = undef, i_ = 0) = (i_ == 0) ? max_index(values, values[0], 0, 1) : (i_ >= len(values)) ? idx_max - : (values[i_] > max_) ? max_index(values, values[i_], i_, i_ + 1) + : (values[i_] > max_) ? + max_index(values, values[i_], i_, i_ + 1) : max_index(values, max_, idx_max, i_ + 1); function remove_elements(arr, to_remove) = @@ -137,10 +156,21 @@ function remove_internal_edges(edges) = [for (h = edges) if (!contains(edges, reverse(h))) h]; function plane_project(p, a, b, c) = - let (u = b - a, v = c - a, n = cross(u, v), w = cross(n, u), rel = p - a) [rel * u, rel * w]; + let ( + u = b - a, + v = c - a, + n = cross(u, v), + w = cross(n, u), + rel = p - a + ) [rel * u, rel * w]; function plane_unproject(p, a, b, c) = - let (u = b - a, v = c - a, n = cross(u, v), w = cross(n, u)) a + p[0] * u + p[1] * w; + let ( + u = b - a, + v = c - a, + n = cross(u, v), + w = cross(n, u) + ) a + p[0] * u + p[1] * w; function reverse(arr) = [for (i = [len(arr) - 1:-1:0]) arr[i]]; @@ -157,7 +187,9 @@ function find_first_noncollinear(line, pts, i) = function find_first_noncoplanar(pl, pts, i) = (i >= len(pts)) ? len(pts) - : coplanar(pl, pts[i]) ? find_first_noncoplanar(pl, pts, i + 1) : i; + : coplanar(pl, pts[i]) ? + find_first_noncoplanar(pl, pts, i + 1) + : i; function distance(pl, p) = pl[0] * p - pl[1]; function in_front(pl, p) = distance(pl, p) > epsilon; @@ -166,11 +198,16 @@ function coplanar(pl, p) = abs(distance(pl, p)) <= epsilon; function unit(v) = v / norm(v); function area_2d(a, b, c) = - (a[0] * (b[1] - c[1]) + b[0] * (c[1] - a[1]) + c[0] * (a[1] - b[1])) / 2; + ( + a[0] * (b[1] - c[1]) + b[0] * (c[1] - a[1]) + c[0] * (a[1] - b[1]) + ) / 2; + +function collinear(a, b, c) = + abs(area_2d(a, b, c)) < epsilon; -function collinear(a, b, c) = abs(area_2d(a, b, c)) < epsilon; +function spherical(cart) = + [atan2(cart[1], cart[0]), asin(cart[2])]; -function spherical(cart) = [atan2(cart[1], cart[0]), asin(cart[2])]; function cartesian(sph) = [ cos(sph[1]) * cos(sph[0]), From f9cfd5a4378b6c2a369cdd548cb8842a2b205eb6 Mon Sep 17 00:00:00 2001 From: Cameron K Brooks Date: Fri, 26 Sep 2025 19:45:39 -0400 Subject: [PATCH 06/32] remove old builtins commented --- linalg.scad | 4 ---- 1 file changed, 4 deletions(-) diff --git a/linalg.scad b/linalg.scad index 573652d..c2d5ea1 100644 --- a/linalg.scad +++ b/linalg.scad @@ -1,9 +1,5 @@ // very minimal set of linalg functions needed by so3, se3 etc. -// cross and norm are builtins -//function cross(x,y) = [x[1]*y[2]-x[2]*y[1], x[2]*y[0]-x[0]*y[2], x[0]*y[1]-x[1]*y[0]]; -//function norm(v) = sqrt(v*v); - function vec3(p) = len(p) < 3 ? concat(p, 0) : p; function vec4(p) = let (v3 = vec3(p)) len(v3) < 4 ? concat(v3, 1) : v3; function unit(v) = v / norm(v); From a3729bd50a60a49879b1fdc99dbbf0f106b1b451 Mon Sep 17 00:00:00 2001 From: Cameron K Brooks Date: Fri, 26 Sep 2025 19:50:24 -0400 Subject: [PATCH 07/32] refactor: cosmetic rework of linalg --- linalg.scad | 85 +++++++++++++++++++++++++++++++++++++++++++++-------- 1 file changed, 72 insertions(+), 13 deletions(-) diff --git a/linalg.scad b/linalg.scad index c2d5ea1..9c7f7a0 100644 --- a/linalg.scad +++ b/linalg.scad @@ -1,21 +1,65 @@ -// very minimal set of linalg functions needed by so3, se3 etc. +// ============================================================================ +// Minimal Linear Algebra Utilities (for so3, se3, etc.) +// ---------------------------------------------------------------------------- +// Provides basic vector and matrix operations needed for transformations. +// - Vector constructors and normalization +// - Identity matrices +// - Rotation/translation parts extraction +// - Transpose and inverse (rigid transform) +// - Hadamard (elementwise) product +// ============================================================================ + +epsilon = 1e-9; + +// --- Vector Constructors ---------------------------------------------------- +function vec3(p) = (len(p) < 3) ? concat(p, 0) : p; +function vec4(p) = + let (v3 = vec3(p)) (len(v3) < 4) ? concat(v3, 1) : v3; -function vec3(p) = len(p) < 3 ? concat(p, 0) : p; -function vec4(p) = let (v3 = vec3(p)) len(v3) < 4 ? concat(v3, 1) : v3; function unit(v) = v / norm(v); -function identity3() = [[1, 0, 0], [0, 1, 0], [0, 0, 1]]; -function identity4() = [[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]]; +// --- Identity Matrices ------------------------------------------------------ +function identity3() = + [ + [1, 0, 0], + [0, 1, 0], + [0, 0, 1], + ]; +function identity4() = + [ + [1, 0, 0, 0], + [0, 1, 0, 0], + [0, 0, 1, 0], + [0, 0, 0, 1], + ]; + +// --- Vector Access Helpers -------------------------------------------------- function take3(v) = [v[0], v[1], v[2]]; function tail3(v) = [v[3], v[4], v[5]]; -function rotation_part(m) = [take3(m[0]), take3(m[1]), take3(m[2])]; + +// --- Matrix Part Extraction ------------------------------------------------- +function rotation_part(m) = + [ + take3(m[0]), + take3(m[1]), + take3(m[2]), + ]; + +function translation_part(m) = [m[0][3], m[1][3], m[2][3]]; + +// --- Rotation Metrics ------------------------------------------------------- function rot_trace(m) = m[0][0] + m[1][1] + m[2][2]; function rot_cos_angle(m) = (rot_trace(m) - 1) / 2; -function rotation_part(m) = [take3(m[0]), take3(m[1]), take3(m[2])]; -function translation_part(m) = [m[0][3], m[1][3], m[2][3]]; -function transpose_3(m) = [[m[0][0], m[1][0], m[2][0]], [m[0][1], m[1][1], m[2][1]], [m[0][2], m[1][2], m[2][2]]]; +// --- Matrix Transpose ------------------------------------------------------- +function transpose_3(m) = + [ + [m[0][0], m[1][0], m[2][0]], + [m[0][1], m[1][1], m[2][1]], + [m[0][2], m[1][2], m[2][2]], + ]; + function transpose_4(m) = [ [m[0][0], m[1][0], m[2][0], m[3][0]], @@ -23,8 +67,23 @@ function transpose_4(m) = [m[0][2], m[1][2], m[2][2], m[3][2]], [m[0][3], m[1][3], m[2][3], m[3][3]], ]; -function invert_rt(m) = construct_Rt(transpose_3(rotation_part(m)), -(transpose_3(rotation_part(m)) * translation_part(m))); -function construct_Rt(R, t) = [concat(R[0], t[0]), concat(R[1], t[1]), concat(R[2], t[2]), [0, 0, 0, 1]]; -// Hadamard product of n-dimensional arrays -function hadamard(a, b) = !(len(a) > 0) ? a * b : [for (i = [0:len(a) - 1]) hadamard(a[i], b[i])]; +// --- Rigid Transform Utilities --------------------------------------------- +function invert_rt(m) = + let ( + R = transpose_3(rotation_part(m)), + t = -(R * translation_part(m)) + ) construct_Rt(R, t); + +function construct_Rt(R, t) = + [ + concat(R[0], t[0]), + concat(R[1], t[1]), + concat(R[2], t[2]), + [0, 0, 0, 1], + ]; + +// --- Elementwise Operations ------------------------------------------------- +function hadamard(a, b) = + !(len(a) > 0) ? a * b + : [for (i = [0:len(a) - 1]) hadamard(a[i], b[i])]; From 9e22e89f802dc89426a22e7f18bfc8739ad3c998 Mon Sep 17 00:00:00 2001 From: Cameron K Brooks Date: Fri, 26 Sep 2025 19:53:36 -0400 Subject: [PATCH 08/32] add tests for linalg --- linalg.scad | 4 ++-- tests/test_linalg.scad | 53 ++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 55 insertions(+), 2 deletions(-) create mode 100644 tests/test_linalg.scad diff --git a/linalg.scad b/linalg.scad index 9c7f7a0..90a0d7f 100644 --- a/linalg.scad +++ b/linalg.scad @@ -84,6 +84,6 @@ function construct_Rt(R, t) = ]; // --- Elementwise Operations ------------------------------------------------- +// Hadamard product: works recursively on arrays, multiplies scalars directly. function hadamard(a, b) = - !(len(a) > 0) ? a * b - : [for (i = [0:len(a) - 1]) hadamard(a[i], b[i])]; + is_list(a) ? [for (i = [0:len(a) - 1]) hadamard(a[i], b[i])] : a * b; diff --git a/tests/test_linalg.scad b/tests/test_linalg.scad new file mode 100644 index 0000000..62b6b10 --- /dev/null +++ b/tests/test_linalg.scad @@ -0,0 +1,53 @@ +include <../linalg.scad>; + +// --- Vector constructors ---------------------------------------------------- +echo("vec3([1,2]) =", vec3([1,2])); // expect [1,2,0] +echo("vec4([1,2,3]) =", vec4([1,2,3])); // expect [1,2,3,1] +echo("unit([3,0,0]) =", unit([3,0,0])); // expect [1,0,0] + +// --- Identity matrices ------------------------------------------------------ +echo("identity3() =", identity3()); // expect 3x3 identity +echo("identity4() =", identity4()); // expect 4x4 identity + +// --- Vector access ---------------------------------------------------------- +echo("take3([9,8,7,6]) =", take3([9,8,7,6])); // expect [9,8,7] +echo("tail3([0,1,2,3,4,5]) =", tail3([0,1,2,3,4,5])); // expect [3,4,5] + +// --- Matrix parts ----------------------------------------------------------- +M = [ + [1,0,0,5], + [0,1,0,6], + [0,0,1,7], + [0,0,0,1] +]; +echo("rotation_part(M) =", rotation_part(M)); // expect identity3() +echo("translation_part(M) =", translation_part(M)); // expect [5,6,7] + +// --- Rotation metrics ------------------------------------------------------- +R = identity3(); +echo("rot_trace(R) =", rot_trace(R)); // expect 3 +echo("rot_cos_angle(R) =", rot_cos_angle(R)); // expect 1 + +// --- Transpose -------------------------------------------------------------- +A3 = [[1,2,3],[4,5,6],[7,8,9]]; +echo("transpose_3(A3) =", transpose_3(A3)); // expect [[1,4,7],[2,5,8],[3,6,9]] + +A4 = identity4(); +echo("transpose_4(identity4) =", transpose_4(A4)); // expect identity4() + +// --- Rigid transform inverse/construct ------------------------------------- +echo("invert_rt(M) =", invert_rt(M)); +// expect transform with R=I, t=[-5,-6,-7] + +R2 = [[0,-1,0],[1,0,0],[0,0,1]]; // 90° rot about Z +t2 = [10,0,0]; +Rt = construct_Rt(R2,t2); +echo("construct_Rt(R2,t2) =", Rt); + +// --- Hadamard product ------------------------------------------------------ +echo("hadamard([1,2,3],[4,5,6]) =", hadamard([1,2,3],[4,5,6])); +// expect [4,10,18] + +echo("hadamard([[1,2],[3,4]], [[5,6],[7,8]]) =", + hadamard([[1,2],[3,4]], [[5,6],[7,8]])); +// expect [[5,12],[21,32]] From f78a5a09e072a4d62020c9308c03b259762b391c Mon Sep 17 00:00:00 2001 From: Cameron K Brooks Date: Fri, 26 Sep 2025 19:59:25 -0400 Subject: [PATCH 09/32] tests & refactor for lists --- lists.scad | 67 +++++++++++++++++++++++++++++-------------- tests/test_lists.scad | 34 ++++++++++++++++++++++ 2 files changed, 80 insertions(+), 21 deletions(-) create mode 100644 tests/test_lists.scad diff --git a/lists.scad b/lists.scad index 683e11b..cad7650 100644 --- a/lists.scad +++ b/lists.scad @@ -1,46 +1,71 @@ -// List helpers +// ============================================================================ +// List Utilities +// ---------------------------------------------------------------------------- +// Provides small helper functions for working with lists: +// - Flattening +// - Ranges +// - Reversals and subarrays +// - Element updates and removals +// ============================================================================ +// --- Flatten --------------------------------------------------------------- /*! - Flattens a list one level: + Flatten a list one level. - flatten([[0,1],[2,3]]) => [0,1,2,3] + Example: + flatten([[0,1],[2,3]]) => [0,1,2,3] */ -function flatten(list) = [for (i = list, v = i) v]; +function flatten(list) = [for (sub = list, v = sub) v]; +// --- Range ----------------------------------------------------------------- /*! - Creates a list from a range: + Create a list from a range. - range([0:2:6]) => [0,2,4,6] + Example: + range([0:2:6]) => [0,2,4,6] */ function range(r) = [for (x = r) x]; +// --- Reverse --------------------------------------------------------------- /*! - Reverses a list: + Reverse the order of elements. - reverse([1,2,3]) => [3,2,1] + Example: + reverse([1,2,3]) => [3,2,1] */ -function reverse(list) = [for (i = [len(list) - 1:-1:0]) list[i]]; +function reverse(list) = + [for (i = [len(list) - 1:-1:0]) list[i]]; +// --- Subarray -------------------------------------------------------------- /*! - Extracts a subarray from index begin (inclusive) to end (exclusive) - FIXME: Change name to use list instead of array? + Extract a subarray from index `begin` (inclusive) to `end` (exclusive). - subarray([1,2,3,4], 1, 2) => [2,3] + Notes: + - If `end < 0`, uses `len(list)`. + - FIXME: Consider renaming to `sublist` for clarity. + + Example: + subarray([1,2,3,4], 1, 3) => [2,3] */ function subarray(list, begin = 0, end = -1) = - [ - let (end = end < 0 ? len(list) : end) for (i = [begin:1:end - 1]) list[i], - ]; + let (end = (end < 0) ? len(list) : end) [for (i = [begin:end - 1]) list[i]]; +// --- Set ------------------------------------------------------------------- /*! - Returns a copy of a list with the element at index i set to x + Return a copy of a list with the element at index `i` replaced by `x`. - set([1,2,3,4], 2, 5) => [1,2,5,4] + Example: + set([1,2,3,4], 2, 5) => [1,2,5,4] */ -function set(list, i, x) = [for (i_ = [0:len(list) - 1]) i == i_ ? x : list[i_]]; +function set(list, i, x) = + [for (j = [0:len(list) - 1]) (i == j) ? x : list[j]]; +// --- Remove --------------------------------------------------------------- /*! - Remove element from the list by index. - remove([4,3,2,1],1) => [4,2,1] + Remove the element at index `i`. + + Example: + remove([4,3,2,1], 1) => [4,2,1] */ -function remove(list, i) = [for (i_ = [0:1:len(list) - 2]) list[i_ < i ? i_ : i_ + 1]]; +function remove(list, i) = + [for (j = [0:len(list) - 2]) list[ (j < i) ? j : j + 1]]; diff --git a/tests/test_lists.scad b/tests/test_lists.scad new file mode 100644 index 0000000..977ffdf --- /dev/null +++ b/tests/test_lists.scad @@ -0,0 +1,34 @@ +include <../lists.scad>; + +// --- Flatten --------------------------------------------------------------- +echo("flatten([[0,1],[2,3]]) =", flatten([[0, 1], [2, 3]])); +// expect [0,1,2,3] + +// --- Range ----------------------------------------------------------------- +echo("range([0:2:6]) =", range([0:2:6])); +// expect [0,2,4,6] + +// --- Reverse --------------------------------------------------------------- +echo("reverse([1,2,3]) =", reverse([1, 2, 3])); +// expect [3,2,1] + +// --- Subarray -------------------------------------------------------------- +echo("subarray([1,2,3,4],1,3) =", subarray([1, 2, 3, 4], 1, 3)); +// expect [2,3] + +echo("subarray([1,2,3,4],0,-1) =", subarray([1, 2, 3, 4], 0, -1)); +// expect [1,2,3,4] + +// --- Set ------------------------------------------------------------------- +echo("set([1,2,3,4],2,5) =", set([1, 2, 3, 4], 2, 5)); +// expect [1,2,5,4] + +// --- Remove --------------------------------------------------------------- +echo("remove([4,3,2,1],1) =", remove([4, 3, 2, 1], 1)); +// expect [4,2,1] + +echo("remove([10,20,30,40,50],0) =", remove([10, 20, 30, 40, 50], 0)); +// expect [20,30,40,50] + +echo("remove([10,20,30,40,50],4) =", remove([10, 20, 30, 40, 50], 4)); +// expect [10,20,30,40] From 0de46d8ed186537f15e2e36b01a7f5242b611fb1 Mon Sep 17 00:00:00 2001 From: Cameron K Brooks Date: Fri, 26 Sep 2025 20:29:41 -0400 Subject: [PATCH 10/32] enh: update and format readme --- README.md | 53 +++++++++++++++++++++-------------------------------- 1 file changed, 21 insertions(+), 32 deletions(-) diff --git a/README.md b/README.md index 8fdd4e3..66d8c1d 100644 --- a/README.md +++ b/README.md @@ -1,10 +1,8 @@ -scad-utils -========== +# scad-utils Utility libraries for OpenSCAD -Morphology ----------- +## Morphology contains basic 2D morphology operations @@ -16,9 +14,8 @@ contains basic 2D morphology operations - positive values of d places the shell on the outside - negative values of d places the shell on the inside - center=true and positive d places the shell centered on the edge - - -### Examples + +## Examples With a basic sample polygon shape, @@ -28,54 +25,46 @@ With a basic sample polygon shape, and `$fn=32;`. - -* `inset(d=0.3) shape();` +- `inset(d=0.3) shape();` ![](http://oskarlinde.github.io/scad-utils/img/morph-0.png) - -* `outset(d=0.3) shape();` +- `outset(d=0.3) shape();` ![](http://oskarlinde.github.io/scad-utils/img/morph-1.png) - -* `rounding(r=0.3) shape();` +- `rounding(r=0.3) shape();` ![](http://oskarlinde.github.io/scad-utils/img/morph-2.png) - -* `fillet(r=0.3) shape();` +- `fillet(r=0.3) shape();` ![](http://oskarlinde.github.io/scad-utils/img/morph-3.png) - -*`shell(d=0.3) shape();` +- `shell(d=0.3) shape();` ![](http://oskarlinde.github.io/scad-utils/img/morph-4.png) - -*`shell(d=-0.3) shape();` +- `shell(d=-0.3) shape();` ![](http://oskarlinde.github.io/scad-utils/img/morph-5.png) - -*`shell(d=0.3,center=true) shape();` +- `shell(d=0.3,center=true) shape();` ![](http://oskarlinde.github.io/scad-utils/img/morph-6.png) - -Mirror ------- +## Mirror contains simple mirroring functions - mirror_x() - mirror_y() - mirror_z() - -example: +- mirror_x() +- mirror_y() +- mirror_z() - module arrow(l=1,w=.6,t=0.15) { - mirror_y() polygon([[0,0],[l,0],[l-w/2,w/2],[l-w/2-sqrt(2)*t,w/2],[l-t/2-sqrt(2)*t,t/2],[0,t/2]]); - } +example: +```c +module arrow(l=1,w=.6,t=0.15) { + mirror_y() polygon([[0,0],[l,0],[l-w/2,w/2],[l-w/2-sqrt(2)*t,w/2],[l-t/2-sqrt(2)*t,t/2],[0,t/2]]); +} +``` From 025e556f8d16ccf0e1e0c0d1cb0eb60bcd051528 Mon Sep 17 00:00:00 2001 From: Cameron K Brooks Date: Fri, 26 Sep 2025 20:30:54 -0400 Subject: [PATCH 11/32] refactor and tests for mirror --- mirror.scad | 20 ++++++++++++-------- tests/test_mirror.scad | 42 ++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 54 insertions(+), 8 deletions(-) create mode 100644 tests/test_mirror.scad diff --git a/mirror.scad b/mirror.scad index 40f4017..9b739df 100644 --- a/mirror.scad +++ b/mirror.scad @@ -1,12 +1,14 @@ -// Copyright (c) 2013 Oskar Linde. All rights reserved. -// License: BSD -// -// This library contains simple mirroring functions -// -// mirror_x() -// mirror_y() -// mirror_z() +// ============================================================================ +// Mirror Utilities +// ---------------------------------------------------------------------------- +// Provides simple mirroring modules around the X, Y, and Z axes. +// Each module duplicates its children and adds a mirrored copy. +// - mirror_x(): mirror across the YZ-plane (flip X) +// - mirror_y(): mirror across the XZ-plane (flip Y) +// - mirror_z(): mirror across the XY-plane (flip Z) +// ============================================================================ +// --- Mirror across X-axis --------------------------------------------------- module mirror_x() { union() { children(); @@ -14,6 +16,7 @@ module mirror_x() { } } +// --- Mirror across Y-axis --------------------------------------------------- module mirror_y() { union() { children(); @@ -21,6 +24,7 @@ module mirror_y() { } } +// --- Mirror across Z-axis --------------------------------------------------- module mirror_z() { union() { children(); diff --git a/tests/test_mirror.scad b/tests/test_mirror.scad new file mode 100644 index 0000000..b144197 --- /dev/null +++ b/tests/test_mirror.scad @@ -0,0 +1,42 @@ +include <../mirror.scad> + +// Simple reference shape offset from the origin so the mirror is obvious +module sample_shape() { + translate([20, 10, 5]) { + cube([10, 6, 6], center=false); + translate([10, 3, 6]) sphere(r=4, $fn=24); + } +} + +// Tiny axes helper for orientation +module axes(len = 25) { + color("red") cube([len, 0.8, 0.8], center=false); // +X + color("green") cube([0.8, len, 0.8], center=false); // +Y + color("blue") cube([0.8, 0.8, len], center=false); // +Z +} + +// Arrange three demos side-by-side +translate([-60, 0, 0]) { + axes(); + mirror_x() sample_shape(); +} + +translate([0, 0, 0]) { + axes(); + mirror_y() sample_shape(); +} + +translate([60, 0, 0]) { + axes(); + mirror_z() sample_shape(); +} + +// Compare with a single plain mirror +translate([-60, -30, 0]) { + axes(); + color("orange") + union() { + sample_shape(); + mirror([1, 0, 0]) sample_shape(); + } +} From cd902b1d74da9957e64b35b88aa63b87e2dbc02d Mon Sep 17 00:00:00 2001 From: Cameron K Brooks Date: Fri, 26 Sep 2025 20:31:15 -0400 Subject: [PATCH 12/32] split morphology to tests --- morphology.scad | 37 ---------------------------------- tests/test_morphology.scad | 41 ++++++++++++++++++++++++++++++++++++++ 2 files changed, 41 insertions(+), 37 deletions(-) create mode 100644 tests/test_morphology.scad diff --git a/morphology.scad b/morphology.scad index 53871d7..3df115f 100644 --- a/morphology.scad +++ b/morphology.scad @@ -72,40 +72,3 @@ module inverse() { children(); } } - -// TEST CODE - -use - -module arrow(l = 1, w = .6, t = 0.15) { - mirror_y() polygon([[0, 0], [l, 0], [l - w / 2, w / 2], [l - w / 2 - sqrt(2) * t, w / 2], [l - t / 2 - sqrt(2) * t, t / 2], [0, t / 2]]); -} - -module shape() { - polygon([[0, 0], [1, 0], [1.5, 1], [2.5, 1], [2, -1], [0, -1]]); -} - -if (0) - assign ($fn = 32) { - - for (p = [0:10 * 3 - 1]) - assign (o = floor(p / 3)) { - translate([(p % 3) * 2.5, -o * 3]) { - //%if (p % 3 == 1) translate([0,0,1]) shape(); - if (p % 3 == 0) shape(); - if (p % 3 == 1) translate([0.6, 0]) arrow(); - if (p % 3 == 2) { - if (o == 0) inset(d=0.3) shape(); - if (o == 1) outset(d=0.3) shape(); - if (o == 2) rounding(r=0.3) shape(); - if (o == 3) fillet(r=0.3) shape(); - if (o == 4) shell(d=0.3) shape(); - if (o == 5) shell(d=-0.3) shape(); - if (o == 6) shell(d=0.3, center=true) shape(); - if (o == 7) rounding(r=0.3) fillet(r=0.3) shape(); - if (o == 8) shell(d=0.3, center=true) fillet(r=0.3) rounding(r=0.3) shape(); - if (o == 9) shell(d=-0.3) fillet(r=0.3) rounding(r=0.3) shape(); - } - } - } - } diff --git a/tests/test_morphology.scad b/tests/test_morphology.scad new file mode 100644 index 0000000..2d559f1 --- /dev/null +++ b/tests/test_morphology.scad @@ -0,0 +1,41 @@ +use <../morphology.scad> + + +// TEST CODE + +use <../mirror.scad> + +module arrow(l = 1, w = .6, t = 0.15) { + mirror_y() polygon([[0, 0], [l, 0], [l - w / 2, w / 2], [l - w / 2 - sqrt(2) * t, w / 2], [l - t / 2 - sqrt(2) * t, t / 2], [0, t / 2]]); +} + +module shape() { + polygon([[0, 0], [1, 0], [1.5, 1], [2.5, 1], [2, -1], [0, -1]]); +} + +debug = true; + +if (debug) + assign ($fn = 32) { + + for (p = [0:10 * 3 - 1]) + assign (o = floor(p / 3)) { + translate([(p % 3) * 2.5, -o * 3]) { + //%if (p % 3 == 1) translate([0,0,1]) shape(); + if (p % 3 == 0) shape(); + if (p % 3 == 1) translate([0.6, 0]) arrow(); + if (p % 3 == 2) { + if (o == 0) inset(d=0.3) shape(); + if (o == 1) outset(d=0.3) shape(); + if (o == 2) rounding(r=0.3) shape(); + if (o == 3) fillet(r=0.3) shape(); + if (o == 4) shell(d=0.3) shape(); + if (o == 5) shell(d=-0.3) shape(); + if (o == 6) shell(d=0.3, center=true) shape(); + if (o == 7) rounding(r=0.3) fillet(r=0.3) shape(); + if (o == 8) shell(d=0.3, center=true) fillet(r=0.3) rounding(r=0.3) shape(); + if (o == 9) shell(d=-0.3) fillet(r=0.3) rounding(r=0.3) shape(); + } + } + } + } From 5c18c04823025f732950c146f9c813f19e9e6d51 Mon Sep 17 00:00:00 2001 From: Cameron K Brooks Date: Fri, 26 Sep 2025 20:31:50 -0400 Subject: [PATCH 13/32] tests: morphology test update to remove assign() --- tests/test_morphology.scad | 49 ++++++++++++++++++-------------------- 1 file changed, 23 insertions(+), 26 deletions(-) diff --git a/tests/test_morphology.scad b/tests/test_morphology.scad index 2d559f1..42084ae 100644 --- a/tests/test_morphology.scad +++ b/tests/test_morphology.scad @@ -1,8 +1,4 @@ use <../morphology.scad> - - -// TEST CODE - use <../mirror.scad> module arrow(l = 1, w = .6, t = 0.15) { @@ -15,27 +11,28 @@ module shape() { debug = true; -if (debug) - assign ($fn = 32) { - - for (p = [0:10 * 3 - 1]) - assign (o = floor(p / 3)) { - translate([(p % 3) * 2.5, -o * 3]) { - //%if (p % 3 == 1) translate([0,0,1]) shape(); - if (p % 3 == 0) shape(); - if (p % 3 == 1) translate([0.6, 0]) arrow(); - if (p % 3 == 2) { - if (o == 0) inset(d=0.3) shape(); - if (o == 1) outset(d=0.3) shape(); - if (o == 2) rounding(r=0.3) shape(); - if (o == 3) fillet(r=0.3) shape(); - if (o == 4) shell(d=0.3) shape(); - if (o == 5) shell(d=-0.3) shape(); - if (o == 6) shell(d=0.3, center=true) shape(); - if (o == 7) rounding(r=0.3) fillet(r=0.3) shape(); - if (o == 8) shell(d=0.3, center=true) fillet(r=0.3) rounding(r=0.3) shape(); - if (o == 9) shell(d=-0.3) fillet(r=0.3) rounding(r=0.3) shape(); - } - } +if (debug) { + $fn = 32; + + for (p = [0:10 * 3 - 1]) { + o = floor(p / 3); + + translate([(p % 3) * 2.5, -o * 3]) { + //%if (p % 3 == 1) translate([0,0,1]) shape(); + if (p % 3 == 0) shape(); + if (p % 3 == 1) translate([0.6, 0]) arrow(); + if (p % 3 == 2) { + if (o == 0) inset(d=0.3) shape(); + if (o == 1) outset(d=0.3) shape(); + if (o == 2) rounding(r=0.3) shape(); + if (o == 3) fillet(r=0.3) shape(); + if (o == 4) shell(d=0.3) shape(); + if (o == 5) shell(d=-0.3) shape(); + if (o == 6) shell(d=0.3, center=true) shape(); + if (o == 7) rounding(r=0.3) fillet(r=0.3) shape(); + if (o == 8) shell(d=0.3, center=true) fillet(r=0.3) rounding(r=0.3) shape(); + if (o == 9) shell(d=-0.3) fillet(r=0.3) rounding(r=0.3) shape(); } + } } +} From e92a1f4612614f291d96cbf80e630cbfc5281d85 Mon Sep 17 00:00:00 2001 From: Cameron K Brooks Date: Fri, 26 Sep 2025 20:41:30 -0400 Subject: [PATCH 14/32] add arrow to mirror test --- tests/test_mirror.scad | 32 +++++++++++++++++++++++++++----- 1 file changed, 27 insertions(+), 5 deletions(-) diff --git a/tests/test_mirror.scad b/tests/test_mirror.scad index b144197..ba6e4b3 100644 --- a/tests/test_mirror.scad +++ b/tests/test_mirror.scad @@ -9,10 +9,21 @@ module sample_shape() { } // Tiny axes helper for orientation -module axes(len = 25) { - color("red") cube([len, 0.8, 0.8], center=false); // +X - color("green") cube([0.8, len, 0.8], center=false); // +Y - color("blue") cube([0.8, 0.8, len], center=false); // +Z +module axes(len = 25, thick = 0.8) { + // +X axis (red) + color("red") + translate([len / 2, 0, 0]) + cube([len, thick, thick], center=true); + + // +Y axis (green) + color("green") + translate([0, len / 2, 0]) + cube([thick, len, thick], center=true); + + // +Z axis (blue) + color("blue") + translate([0, 0, len / 2]) + cube([thick, thick, len], center=true); } // Arrange three demos side-by-side @@ -32,7 +43,7 @@ translate([60, 0, 0]) { } // Compare with a single plain mirror -translate([-60, -30, 0]) { +translate([-60, -40, 0]) { axes(); color("orange") union() { @@ -40,3 +51,14 @@ translate([-60, -30, 0]) { mirror([1, 0, 0]) sample_shape(); } } + +module arrow(l = 1, w = .6, t = 0.15) { + mirror_y() polygon([[0, 0], [l, 0], [l - w / 2, w / 2], [l - w / 2 - sqrt(2) * t, w / 2], [l - t / 2 - sqrt(2) * t, t / 2], [0, t / 2]]); +} + +// Compare with a single plain mirror +translate([60, -40, 0]) { + axes(); + color("orange") + arrow(l=20, w=10, t=2); +} From e01d66b4a9e0a4ea9b73ad770739b24e60e9f496 Mon Sep 17 00:00:00 2001 From: Cameron K Brooks Date: Fri, 26 Sep 2025 20:48:02 -0400 Subject: [PATCH 15/32] feat: add coloring function to mirror --- README.md | 2 +- mirror.scad | 28 +++++++++++++++-------- tests/test_mirror.scad | 52 ++++++++++++++++++++++++++++-------------- 3 files changed, 55 insertions(+), 27 deletions(-) diff --git a/README.md b/README.md index 66d8c1d..da9a525 100644 --- a/README.md +++ b/README.md @@ -63,7 +63,7 @@ contains simple mirroring functions example: -```c +``` module arrow(l=1,w=.6,t=0.15) { mirror_y() polygon([[0,0],[l,0],[l-w/2,w/2],[l-w/2-sqrt(2)*t,w/2],[l-t/2-sqrt(2)*t,t/2],[0,t/2]]); } diff --git a/mirror.scad b/mirror.scad index 9b739df..53945c2 100644 --- a/mirror.scad +++ b/mirror.scad @@ -3,31 +3,41 @@ // ---------------------------------------------------------------------------- // Provides simple mirroring modules around the X, Y, and Z axes. // Each module duplicates its children and adds a mirrored copy. -// - mirror_x(): mirror across the YZ-plane (flip X) -// - mirror_y(): mirror across the XZ-plane (flip Y) -// - mirror_z(): mirror across the XY-plane (flip Z) +// Optional: pass `col="colorname"` to tint the mirrored copy. +// - mirror_x([col]): mirror across the YZ-plane (flip X) +// - mirror_y([col]): mirror across the XZ-plane (flip Y) +// - mirror_z([col]): mirror across the XY-plane (flip Z) // ============================================================================ // --- Mirror across X-axis --------------------------------------------------- -module mirror_x() { +module mirror_x(col = undef) { union() { children(); - scale([-1, 1, 1]) children(); + if (is_undef(col)) + scale([-1, 1, 1]) children(); + else + color(col) scale([-1, 1, 1]) children(); } } // --- Mirror across Y-axis --------------------------------------------------- -module mirror_y() { +module mirror_y(col = undef) { union() { children(); - scale([1, -1, 1]) children(); + if (is_undef(col)) + scale([1, -1, 1]) children(); + else + color(col) scale([1, -1, 1]) children(); } } // --- Mirror across Z-axis --------------------------------------------------- -module mirror_z() { +module mirror_z(col = undef) { union() { children(); - scale([1, 1, -1]) children(); + if (is_undef(col)) + scale([1, 1, -1]) children(); + else + color(col) scale([1, 1, -1]) children(); } } diff --git a/tests/test_mirror.scad b/tests/test_mirror.scad index ba6e4b3..4993c26 100644 --- a/tests/test_mirror.scad +++ b/tests/test_mirror.scad @@ -1,6 +1,8 @@ include <../mirror.scad> -// Simple reference shape offset from the origin so the mirror is obvious +// --------------------------------------------------------------------------- +// Sample reference shape +// --------------------------------------------------------------------------- module sample_shape() { translate([20, 10, 5]) { cube([10, 6, 6], center=false); @@ -8,7 +10,9 @@ module sample_shape() { } } -// Tiny axes helper for orientation +// --------------------------------------------------------------------------- +// Axes helper for orientation +// --------------------------------------------------------------------------- module axes(len = 25, thick = 0.8) { // +X axis (red) color("red") @@ -26,39 +30,53 @@ module axes(len = 25, thick = 0.8) { cube([thick, thick, len], center=true); } -// Arrange three demos side-by-side +// --------------------------------------------------------------------------- +// Demonstrations +// --------------------------------------------------------------------------- + +// Mirror across X-axis translate([-60, 0, 0]) { axes(); - mirror_x() sample_shape(); + mirror_x("teal") sample_shape(); } +// Mirror across Y-axis translate([0, 0, 0]) { axes(); - mirror_y() sample_shape(); + mirror_y("indigo") sample_shape(); } +// Mirror across Z-axis translate([60, 0, 0]) { axes(); - mirror_z() sample_shape(); + mirror_z("salmon") sample_shape(); } -// Compare with a single plain mirror +// Compare with a plain built-in mirror translate([-60, -40, 0]) { axes(); - color("orange") - union() { - sample_shape(); - mirror([1, 0, 0]) sample_shape(); - } + union() { + sample_shape(); + mirror([1, 0, 0]) color("MediumAquamarine") sample_shape(); + } } -module arrow(l = 1, w = .6, t = 0.15) { - mirror_y() polygon([[0, 0], [l, 0], [l - w / 2, w / 2], [l - w / 2 - sqrt(2) * t, w / 2], [l - t / 2 - sqrt(2) * t, t / 2], [0, t / 2]]); +// 2D arrow example using mirror_y +module arrow(l = 1, w = 0.6, t = 0.15) { + mirror_y("orange") + polygon( + [ + [0, 0], + [l, 0], + [l - w / 2, w / 2], + [l - w / 2 - sqrt(2) * t, w / 2], + [l - t / 2 - sqrt(2) * t, t / 2], + [0, t / 2], + ] + ); } -// Compare with a single plain mirror translate([60, -40, 0]) { axes(); - color("orange") - arrow(l=20, w=10, t=2); + arrow(l=20, w=10, t=2); } From ffccd27b7cf7b6d8591cb90395ed646d149817a4 Mon Sep 17 00:00:00 2001 From: Cameron K Brooks Date: Fri, 26 Sep 2025 20:54:05 -0400 Subject: [PATCH 16/32] update test_morphology to use mirror arrow --- tests/test_morphology.scad | 17 ++++++----------- 1 file changed, 6 insertions(+), 11 deletions(-) diff --git a/tests/test_morphology.scad b/tests/test_morphology.scad index 42084ae..d00e487 100644 --- a/tests/test_morphology.scad +++ b/tests/test_morphology.scad @@ -1,9 +1,5 @@ use <../morphology.scad> -use <../mirror.scad> - -module arrow(l = 1, w = .6, t = 0.15) { - mirror_y() polygon([[0, 0], [l, 0], [l - w / 2, w / 2], [l - w / 2 - sqrt(2) * t, w / 2], [l - t / 2 - sqrt(2) * t, t / 2], [0, t / 2]]); -} +use module shape() { polygon([[0, 0], [1, 0], [1.5, 1], [2.5, 1], [2, -1], [0, -1]]); @@ -18,18 +14,17 @@ if (debug) { o = floor(p / 3); translate([(p % 3) * 2.5, -o * 3]) { - //%if (p % 3 == 1) translate([0,0,1]) shape(); if (p % 3 == 0) shape(); - if (p % 3 == 1) translate([0.6, 0]) arrow(); + if (p % 3 == 1) translate([0.6, 0]) color("grey") arrow(); if (p % 3 == 2) { if (o == 0) inset(d=0.3) shape(); if (o == 1) outset(d=0.3) shape(); if (o == 2) rounding(r=0.3) shape(); if (o == 3) fillet(r=0.3) shape(); - if (o == 4) shell(d=0.3) shape(); - if (o == 5) shell(d=-0.3) shape(); - if (o == 6) shell(d=0.3, center=true) shape(); - if (o == 7) rounding(r=0.3) fillet(r=0.3) shape(); + if (o == 4) rounding(r=0.3) fillet(r=0.3) shape(); + if (o == 5) shell(d=0.3) shape(); + if (o == 6) shell(d=-0.3) shape(); + if (o == 7) shell(d=0.3, center=true) shape(); if (o == 8) shell(d=0.3, center=true) fillet(r=0.3) rounding(r=0.3) shape(); if (o == 9) shell(d=-0.3) fillet(r=0.3) rounding(r=0.3) shape(); } From dd015d21896bf4a3c4c87787477332b750cd6fe7 Mon Sep 17 00:00:00 2001 From: Cameron K Brooks Date: Fri, 26 Sep 2025 21:02:15 -0400 Subject: [PATCH 17/32] revise morphology cosmetically --- morphology.scad | 50 ++++++++++++++++++++++++++------------ tests/test_morphology.scad | 24 ++++++++++++++---- 2 files changed, 53 insertions(+), 21 deletions(-) diff --git a/morphology.scad b/morphology.scad index 3df115f..a376249 100644 --- a/morphology.scad +++ b/morphology.scad @@ -1,17 +1,21 @@ -// Copyright (c) 2013 Oskar Linde. All rights reserved. -// License: BSD +// ============================================================================ +// Morphology Utilities (2D) +// ---------------------------------------------------------------------------- +// Provides basic 2D morphological operations: +// - outset(d=1): offset polygon outward +// - inset(d=1): offset polygon inward +// - fillet(r=1): fillet concave corners +// - rounding(r=1): round convex corners +// - shell(d, center=false): create a band (shell) around a polygon // -// This library contains basic 2D morphology operations -// -// outset(d=1) - creates a polygon at an offset d outside a 2D shape -// inset(d=1) - creates a polygon at an offset d inside a 2D shape -// fillet(r=1) - adds fillets of radius r to all concave corners of a 2D shape -// rounding(r=1) - adds rounding to all convex corners of a 2D shape -// shell(d,center=false) - makes a shell of width d along the edge of a 2D shape -// - positive values of d places the shell on the outside -// - negative values of d places the shell on the inside -// - center=true and positive d places the shell centered on the edge +// Notes: +// - Works around version-specific issues with minkowski/projection +// - Internal `_inverse()` helper used for inset/boolean operations +// ============================================================================ + +// --- Outset / Inset --------------------------------------------------------- +// Outset: grows a shape outward by distance d module outset(d = 1) { // Bug workaround for older OpenSCAD versions if (version_num() < 20130424) @@ -23,25 +27,38 @@ module outset(d = 1) { } } +// Helper for older OpenSCAD: emulate outset by extrusion + projection module outset_extruded(d = 1) { - projection(cut=true) minkowski() { + projection(cut=true) + minkowski() { cylinder(r=d); linear_extrude(center=true) children(); } } +// Inset: shrinks a shape inward by distance d module inset(d = 1) { - render() inverse() outset(d=d) inverse() children(); + render() _inverse() outset(d=d) _inverse() children(); } +// --- Corner Modifiers ------------------------------------------------------- + +// Fillet: adds arcs of radius r to concave corners module fillet(r = 1) { inset(d=r) render() outset(d=r) children(); } +// Rounding: rounds convex corners with arcs of radius r module rounding(r = 1) { outset(d=r) inset(d=r) children(); } +// --- Shells ---------------------------------------------------------------- + +// Shell: creates a band of width d around the polygon edge +// - d > 0, center=false: shell on the outside +// - d < 0, center=false: shell on the inside +// - center=true: shell straddles the edge module shell(d, center = false) { if (center && d > 0) { difference() { @@ -64,9 +81,10 @@ module shell(d, center = false) { if (d == 0) children(); } -// Below are for internal use only +// --- Internal Helpers ------------------------------------------------------- -module inverse() { +// Inverse: covers the plane and subtracts children (used in inset) +module _inverse() { difference() { square(1e5, center=true); children(); diff --git a/tests/test_morphology.scad b/tests/test_morphology.scad index d00e487..299eb53 100644 --- a/tests/test_morphology.scad +++ b/tests/test_morphology.scad @@ -1,8 +1,20 @@ use <../morphology.scad> -use +use // for arrow() +// --------------------------------------------------------------------------- +// Base test shape +// --------------------------------------------------------------------------- module shape() { - polygon([[0, 0], [1, 0], [1.5, 1], [2.5, 1], [2, -1], [0, -1]]); + polygon( + [ + [0, 0], + [1, 0], + [1.5, 1], + [2.5, 1], + [2, -1], + [0, -1], + ] + ); } debug = true; @@ -10,13 +22,15 @@ debug = true; if (debug) { $fn = 32; + // 10 groups of 3 columns: [original, arrow, transformed] for (p = [0:10 * 3 - 1]) { - o = floor(p / 3); + o = floor(p / 3); // row index (operation group) translate([(p % 3) * 2.5, -o * 3]) { - if (p % 3 == 0) shape(); - if (p % 3 == 1) translate([0.6, 0]) color("grey") arrow(); + if (p % 3 == 0) shape(); // original + if (p % 3 == 1) translate([0.6, 0]) color("grey") arrow(); // arrow if (p % 3 == 2) { + // transformed if (o == 0) inset(d=0.3) shape(); if (o == 1) outset(d=0.3) shape(); if (o == 2) rounding(r=0.3) shape(); From 09cec7c2a3ecdf105e33861ada3f8e5fa792ad4e Mon Sep 17 00:00:00 2001 From: Cameron K Brooks Date: Fri, 26 Sep 2025 21:23:37 -0400 Subject: [PATCH 18/32] revise se3 cosmetically and move test --- se3.scad | 72 ++++++++++++++++++++++++++++++++++++--------- tests/test_se3.scad | 14 +++++++++ 2 files changed, 72 insertions(+), 14 deletions(-) create mode 100644 tests/test_se3.scad diff --git a/se3.scad b/se3.scad index a153228..3d6f27b 100644 --- a/se3.scad +++ b/se3.scad @@ -1,16 +1,39 @@ +// ============================================================================ +// SE(3) Utilities (Rigid Body Transformations in 3D) +// ---------------------------------------------------------------------------- +// Provides exponential and logarithm maps for SE(3) using so3 + translation. +// - se3_exp(mu): exponential map of 6-vector [t, w] (t=translation, w=rotation) +// - se3_ln(m): logarithm map from 4x4 rigid transform matrix to 6-vector +// ---------------------------------------------------------------------------- +// Dependencies: +// - linalg.scad (matrix utilities) +// - so3.scad (SO(3) exponential/logarithm maps) +// ============================================================================ + use use -function combine_se3_exp(w, ABt) = construct_Rt(rodrigues_so3_exp(w, ABt[0], ABt[1]), ABt[2]); +// --- Core Helpers ----------------------------------------------------------- + +// Combine rotation (Rodrigues) with translated ABt components +function combine_se3_exp(w, ABt) = + construct_Rt( + rodrigues_so3_exp(w, ABt[0], ABt[1]), + ABt[2] + ); + +// --- Exponential Map (Taylor Approximations) -------------------------------- -// [A,B,t] +// Small-angle approx (1st order) function se3_exp_1(t, w) = concat( so3_exp_1(w * w), [t + 0.5 * cross(w, t)] ); +// 2nd order Taylor expansion function se3_exp_2(t, w) = se3_exp_2_0(t, w, w * w); + function se3_exp_2_0(t, w, theta_sq) = se3_exp_23( so3_exp_2(theta_sq), @@ -18,7 +41,13 @@ function se3_exp_2_0(t, w, theta_sq) = t=t, w=w ); -function se3_exp_3(t, w) = se3_exp_3_0(t, w, sqrt(w * w) * 180 / PI, 1 / sqrt(w * w)); +// 3rd order approximation +function se3_exp_3(t, w) = + se3_exp_3_0( + t, w, + theta_deg=sqrt(w * w) * 180 / PI, + inv_theta=1 / sqrt(w * w) + ); function se3_exp_3_0(t, w, theta_deg, inv_theta) = se3_exp_23( @@ -27,48 +56,63 @@ function se3_exp_3_0(t, w, theta_deg, inv_theta) = t=t, w=w ); +// Shared expansion for orders 2–3 function se3_exp_23(AB, C, t, w) = [AB[0], AB[1], t + AB[1] * cross(w, t) + C * cross(w, cross(w, t))]; -function se3_exp(mu) = se3_exp_0(t=take3(mu), w=tail3(mu) / 180 * PI); +// --- Exponential Map Wrapper ------------------------------------------------ + +// se3_exp: exponential map from 6-vector μ = [t, w] +// Converts w from degrees to radians internally +function se3_exp(mu) = + se3_exp_0( + t=take3(mu), + w=tail3(mu) / 180 * PI + ); function se3_exp_0(t, w) = combine_se3_exp( w, - // Evaluate by Taylor expansion when near 0 - w * w < 1e-8 ? se3_exp_1(t, w) - : w * w < 1e-6 ? se3_exp_2(t, w) + (w * w < 1e-8) ? se3_exp_1(t, w) + : (w * w < 1e-6) ? se3_exp_2(t, w) : se3_exp_3(t, w) ); +// --- Logarithm Map ---------------------------------------------------------- + +// Logarithm: returns [t, w] in degrees function se3_ln(m) = se3_ln_to_deg(se3_ln_rad(m)); + function se3_ln_to_deg(v) = concat(take3(v), tail3(v) * 180 / PI); +// Internal: logarithm in radians function se3_ln_rad(m) = se3_ln_0( m, rot=so3_ln_rad(rotation_part(m)) ); + function se3_ln_0(m, rot) = se3_ln_1( m, rot, theta=sqrt(rot * rot) ); + function se3_ln_1(m, rot, theta) = se3_ln_2( m, rot, theta, - shtot=theta > 0.00001 ? sin(theta / 2 * 180 / PI) / theta : 0.5, - halfrotator=so3_exp_rad(rot * -.5) + shtot=(theta > 1e-5) ? sin(theta / 2 * 180 / PI) / theta : 0.5, + halfrotator=so3_exp_rad(rot * -0.5) ); + function se3_ln_2(m, rot, theta, shtot, halfrotator) = concat( ( halfrotator * translation_part(m) - ( - theta > 0.001 ? rot * ( (translation_part(m) * rot) * (1 - 2 * shtot) / (rot * rot)) + (theta > 1e-3) ? + rot * ( (translation_part(m) * rot) * (1 - 2 * shtot) / (rot * rot)) : rot * ( (translation_part(m) * rot) / 24) ) - ) / (2 * shtot), rot + ) / (2 * shtot), + rot ); - -__se3_test = [20, -40, 60, -80, 100, -120]; -echo(UNITTEST_se3=norm(__se3_test - se3_ln(se3_exp(__se3_test))) < 1e-8); diff --git a/tests/test_se3.scad b/tests/test_se3.scad new file mode 100644 index 0000000..10760b5 --- /dev/null +++ b/tests/test_se3.scad @@ -0,0 +1,14 @@ +use <../linalg.scad> +use <../so3.scad> +use <../se3.scad> + +// --- Unit Test for SE(3) --------------------------------------------------- +// Verify se3_ln(se3_exp(mu)) ≈ mu + +mu_test = [20, -40, 60, -80, 100, -120]; +result = se3_ln(se3_exp(mu_test)); + +echo("mu_test =", mu_test); +echo("se3_ln(se3_exp(mu_test)) =", result); +echo("Error norm =", norm(mu_test - result)); +echo("PASS =", norm(mu_test - result) < 1e-8); From 998d5f78722091f7e61766093343c1965312a925 Mon Sep 17 00:00:00 2001 From: Cameron K Brooks Date: Fri, 26 Sep 2025 21:28:39 -0400 Subject: [PATCH 19/32] revise so3 cosmetically and move test --- shapes.scad | 2 +- so3.scad | 80 +++++++++++++++++++++++---------- tests/test_se3.scad | 2 - tests/test_shapes.scad | 0 tests/test_so3.scad | 12 +++++ tests/test_trajectory.scad | 0 tests/test_trajectory_path.scad | 0 tests/test_transformations.scad | 0 8 files changed, 69 insertions(+), 27 deletions(-) create mode 100644 tests/test_shapes.scad create mode 100644 tests/test_so3.scad create mode 100644 tests/test_trajectory.scad create mode 100644 tests/test_trajectory_path.scad create mode 100644 tests/test_transformations.scad diff --git a/shapes.scad b/shapes.scad index fe6e260..f5e3ed0 100644 --- a/shapes.scad +++ b/shapes.scad @@ -14,4 +14,4 @@ function rectangle_profile(size = [1, 1]) = [size[0] / 2, -size[1] / 2], ]; -// FIXME: Move rectangle and rounded rectangle from extrusion +// TODO: Move rectangle and rounded rectangle from extrusion diff --git a/so3.scad b/so3.scad index 66d55c3..f33ee85 100644 --- a/so3.scad +++ b/so3.scad @@ -1,7 +1,19 @@ -// so3 +// ============================================================================ +// SO(3) Utilities (3D Rotation Group) +// ---------------------------------------------------------------------------- +// Provides exponential and logarithm maps for SO(3). +// - so3_exp(w): exponential map from axis-angle (deg) +// - so3_ln(m): logarithm map from rotation matrix → axis-angle (deg) +// ---------------------------------------------------------------------------- +// Dependencies: +// - linalg.scad +// ============================================================================ use +// --- Rodrigues Formula (core exp) ------------------------------------------ + +// Build rotation matrix given axis w and coefficients A, B function rodrigues_so3_exp(w, A, B) = [ [1.0 - B * (w[1] * w[1] + w[2] * w[2]), B * (w[0] * w[1]) - A * w[2], B * (w[0] * w[2]) + A * w[1]], @@ -9,23 +21,26 @@ function rodrigues_so3_exp(w, A, B) = [B * (w[0] * w[2]) - A * w[1], B * (w[1] * w[2]) + A * w[0], 1.0 - B * (w[0] * w[0] + w[1] * w[1])], ]; +// --- Exponential Map -------------------------------------------------------- + +// Wrapper: exponential map with input in degrees function so3_exp(w) = so3_exp_rad(w / 180 * PI); + +// Exponential map in radians, uses Taylor expansions for small angles function so3_exp_rad(w) = combine_so3_exp( w, - w * w < 1e-8 ? so3_exp_1(w * w) - : w * w < 1e-6 ? so3_exp_2(w * w) + (w * w < 1e-8) ? so3_exp_1(w * w) + : (w * w < 1e-6) ? so3_exp_2(w * w) : so3_exp_3(w * w) ); -function combine_so3_exp(w, AB) = rodrigues_so3_exp(w, AB[0], AB[1]); +function combine_so3_exp(w, AB) = + rodrigues_so3_exp(w, AB[0], AB[1]); -// Taylor series expansions close to 0 -function so3_exp_1(theta_sq) = - [ - 1 - 1 / 6 * theta_sq, - 0.5, - ]; +// --- Taylor Expansions (near 0) -------------------------------------------- + +function so3_exp_1(theta_sq) = [1 - theta_sq / 6, 0.5]; function so3_exp_2(theta_sq) = [ @@ -39,11 +54,19 @@ function so3_exp_3_0(theta_deg, inv_theta) = (1 - cos(theta_deg)) * (inv_theta * inv_theta), ]; -function so3_exp_3(theta_sq) = so3_exp_3_0(sqrt(theta_sq) * 180 / PI, 1 / sqrt(theta_sq)); +function so3_exp_3(theta_sq) = + so3_exp_3_0(sqrt(theta_sq) * 180 / PI, 1 / sqrt(theta_sq)); + +// --- Logarithm Map ---------------------------------------------------------- -function rot_axis_part(m) = [m[2][1] - m[1][2], m[0][2] - m[2][0], m[1][0] - m[0][1]] * 0.5; +// Axis part from skew-symmetric difference +function rot_axis_part(m) = + [m[2][1] - m[1][2], m[0][2] - m[2][0], m[1][0] - m[0][1]] * 0.5; +// Logarithm: returns axis-angle in degrees function so3_ln(m) = 180 / PI * so3_ln_rad(m); + +// Logarithm in radians function so3_ln_rad(m) = so3_ln_0( m, @@ -58,28 +81,37 @@ function so3_ln_0(m, cos_angle, preliminary_result) = ); function so3_ln_1(m, cos_angle, preliminary_result, sin_angle_abs) = - cos_angle > sqrt(1 / 2) ? sin_angle_abs > 0 ? preliminary_result * asin(sin_angle_abs) * PI / 180 / sin_angle_abs - : preliminary_result - : cos_angle > -sqrt(1 / 2) ? preliminary_result * acos(cos_angle) * PI / 180 / sin_angle_abs + (cos_angle > sqrt(0.5)) ? + ( + sin_angle_abs > 0 ? + preliminary_result * asin(sin_angle_abs) * PI / 180 / sin_angle_abs + : preliminary_result + ) + : (cos_angle > -sqrt(0.5)) ? + preliminary_result * acos(cos_angle) * PI / 180 / sin_angle_abs : so3_get_symmetric_part_rotation( - preliminary_result, - m, + preliminary_result, m, angle=PI - asin(sin_angle_abs) * PI / 180, d0=m[0][0] - cos_angle, d1=m[1][1] - cos_angle, d2=m[2][2] - cos_angle ); +// --- Symmetric Part Handling (for near-π cases) ----------------------------- + function so3_get_symmetric_part_rotation(preliminary_result, m, angle, d0, d1, d2) = - so3_get_symmetric_part_rotation_0(preliminary_result, angle, so3_largest_column(m, d0, d1, d2)); + so3_get_symmetric_part_rotation_0( + preliminary_result, + angle, + so3_largest_column(m, d0, d1, d2) + ); function so3_get_symmetric_part_rotation_0(preliminary_result, angle, c_max) = - angle * unit(c_max * preliminary_result < 0 ? -c_max : c_max); + angle * unit((c_max * preliminary_result < 0) ? -c_max : c_max); function so3_largest_column(m, d0, d1, d2) = - d0 * d0 > d1 * d1 && d0 * d0 > d2 * d2 ? [d0, (m[1][0] + m[0][1]) / 2, (m[0][2] + m[2][0]) / 2] - : d1 * d1 > d2 * d2 ? [(m[1][0] + m[0][1]) / 2, d1, (m[2][1] + m[1][2]) / 2] + (d0 * d0 > d1 * d1 && d0 * d0 > d2 * d2) ? + [d0, (m[1][0] + m[0][1]) / 2, (m[0][2] + m[2][0]) / 2] + : (d1 * d1 > d2 * d2) ? + [(m[1][0] + m[0][1]) / 2, d1, (m[2][1] + m[1][2]) / 2] : [(m[0][2] + m[2][0]) / 2, (m[2][1] + m[1][2]) / 2, d2]; - -__so3_test = [12, -125, 110]; -echo(UNITTEST_so3=norm(__so3_test - so3_ln(so3_exp(__so3_test))) < 1e-8); diff --git a/tests/test_se3.scad b/tests/test_se3.scad index 10760b5..7c380bf 100644 --- a/tests/test_se3.scad +++ b/tests/test_se3.scad @@ -1,5 +1,3 @@ -use <../linalg.scad> -use <../so3.scad> use <../se3.scad> // --- Unit Test for SE(3) --------------------------------------------------- diff --git a/tests/test_shapes.scad b/tests/test_shapes.scad new file mode 100644 index 0000000..e69de29 diff --git a/tests/test_so3.scad b/tests/test_so3.scad new file mode 100644 index 0000000..5ffc804 --- /dev/null +++ b/tests/test_so3.scad @@ -0,0 +1,12 @@ +use <../so3.scad> + +// --- Unit Test for SO(3) --------------------------------------------------- +// Verify so3_ln(so3_exp(w)) ≈ w + +w_test = [12, -125, 110]; +result = so3_ln(so3_exp(w_test)); + +echo("w_test =", w_test); +echo("so3_ln(so3_exp(w_test)) =", result); +echo("Error norm =", norm(w_test - result)); +echo("PASS =", norm(w_test - result) < 1e-8); diff --git a/tests/test_trajectory.scad b/tests/test_trajectory.scad new file mode 100644 index 0000000..e69de29 diff --git a/tests/test_trajectory_path.scad b/tests/test_trajectory_path.scad new file mode 100644 index 0000000..e69de29 diff --git a/tests/test_transformations.scad b/tests/test_transformations.scad new file mode 100644 index 0000000..e69de29 From 4232ef7998a976ab22d8aeb6b8cc1396e19b995e Mon Sep 17 00:00:00 2001 From: Cameron K Brooks Date: Fri, 26 Sep 2025 21:30:36 -0400 Subject: [PATCH 20/32] revise shapes cosmetically and move test --- shapes.scad | 34 ++++++++++++++++++++++++++++++---- tests/test_shapes.scad | 14 ++++++++++++++ 2 files changed, 44 insertions(+), 4 deletions(-) diff --git a/shapes.scad b/shapes.scad index f5e3ed0..302d425 100644 --- a/shapes.scad +++ b/shapes.scad @@ -1,12 +1,38 @@ -function square(size) = [[-size, -size], [-size, size], [size, size], [size, -size]] / 2; +// ============================================================================ +// Shape Generators (2D) +// ---------------------------------------------------------------------------- +// Provides simple parametric 2D shape generators: +// - square(size): centered square polygon +// - circle(r): circle polygon using current $fn +// - regular(r, n): regular n-gon (circle with fixed $fn) +// - rectangle_profile(size=[w,h]): rectangular profile with anchor at [w/2,0] +// ---------------------------------------------------------------------------- +// Notes: +// - All outputs are polygon point lists (suitable for polygon(), etc.) +// - `circle` and `regular` depend on $fn resolution +// ============================================================================ -function circle(r) = [for (i = [0:$fn - 1]) let (a = i * 360 / $fn) r * [cos(a), sin(a)]]; +// --- Square ----------------------------------------------------------------- +// size: length of edge +function square(size) = + [[-size, -size], [-size, size], [size, size], [size, -size]] / 2; +// --- Circle ----------------------------------------------------------------- +// r: radius +// Uses global $fn for resolution +function circle(r) = + [for (i = [0:$fn - 1]) let (a = i * 360 / $fn) r * [cos(a), sin(a)]]; + +// --- Regular Polygon -------------------------------------------------------- +// r: circumradius +// n: number of sides function regular(r, n) = circle(r, $fn=n); +// --- Rectangle Profile ------------------------------------------------------ +// size = [w, h] +// Anchor point is [w/2, 0] function rectangle_profile(size = [1, 1]) = [ - // The first point is the anchor point, put it on the point corresponding to [cos(0),sin(0)] [size[0] / 2, 0], [size[0] / 2, size[1] / 2], [-size[0] / 2, size[1] / 2], @@ -14,4 +40,4 @@ function rectangle_profile(size = [1, 1]) = [size[0] / 2, -size[1] / 2], ]; -// TODO: Move rectangle and rounded rectangle from extrusion +// TODO: Move rectangle and rounded rectangle from extrusion.scad diff --git a/tests/test_shapes.scad b/tests/test_shapes.scad index e69de29..987b941 100644 --- a/tests/test_shapes.scad +++ b/tests/test_shapes.scad @@ -0,0 +1,14 @@ +include <../shapes.scad> + +$fn = 16; + +echo("square(2) =", square(2)); +echo("circle(r=5) =", circle(5)); +echo("regular(r=5, n=6) =", regular(5, 6)); +echo("rectangle_profile([4,2]) =", rectangle_profile([4, 2])); + +// Quick visualization +translate([-15, 0, 0]) polygon(square(10)); +translate([0, 0, 0]) polygon(circle(5)); +translate([15, 0, 0]) polygon(regular(5, 6)); +translate([30, 0, 0]) polygon(rectangle_profile([8, 4])); From f0b8bc9a895c0f763fdc4afed78c141547d045ff Mon Sep 17 00:00:00 2001 From: Cameron K Brooks Date: Fri, 26 Sep 2025 21:58:44 -0400 Subject: [PATCH 21/32] ensure explicit indices on iterators --- hull.scad | 10 +++++----- linalg.scad | 2 +- lists.scad | 4 ++-- 3 files changed, 8 insertions(+), 8 deletions(-) diff --git a/hull.scad b/hull.scad index 2be5f33..3278b49 100644 --- a/hull.scad +++ b/hull.scad @@ -47,13 +47,13 @@ function convex_hull_iterative_2d(points, polygon, remaining, i_ = 0) = function find_conflicting_segments(points, polygon, point) = [ - for (i = [0:len(polygon) - 1]) let (j = (i + 1) % len(polygon)) if (area_2d(points[polygon[i]], points[polygon[j]], point) < 0) i, + for (i = [0:1:len(polygon) - 1]) let (j = (i + 1) % len(polygon)) if (area_2d(points[polygon[i]], points[polygon[j]], point) < 0) i, ]; function remove_conflicts_and_insert_point(polygon, conflicts, point) = (conflicts[0] == 0) ? let ( - nonconf = [for (i = [0:len(polygon) - 1]) if (!contains(conflicts, i)) i], + nonconf = [for (i = [0:1:len(polygon) - 1]) if (!contains(conflicts, i)) i], indices = concat(nonconf, (nonconf[len(nonconf) - 1] + 1) % len(polygon)) ) concat([for (i = indices) polygon[i]], point) : let ( @@ -77,7 +77,7 @@ function convexhull3d(points) = let (pts2d = [for (p = points) plane_project(p, points[a], points[b], points[c])]) convexhull2d(pts2d) : let ( // build initial tetrahedron - remaining = [for (i = [0:len(points) - 1]) if (i != a && i != b && i != c && i != d) i], + remaining = [for (i = [0:1:len(points) - 1]) if (i != a && i != b && i != c && i != d) i], bc = in_front(pl, points[d]) ? [c, b] : [b, c], b2 = bc[0], c2 = bc[1], @@ -150,7 +150,7 @@ function max_index(values, max_ = undef, idx_max = undef, i_ = 0) = : max_index(values, max_, idx_max, i_ + 1); function remove_elements(arr, to_remove) = - [for (i = [0:len(arr) - 1]) if (!search(i, to_remove)) arr[i]]; + [for (i = [0:1:len(arr) - 1]) if (!search(i, to_remove)) arr[i]]; function remove_internal_edges(edges) = [for (h = edges) if (!contains(edges, reverse(h))) h]; @@ -177,7 +177,7 @@ function reverse(arr) = [for (i = [len(arr) - 1:-1:0]) arr[i]]; function contains(arr, el) = (search([el], arr)[0] != []) ? true : false; function find_conflicts(p, planes) = - [for (i = [0:len(planes) - 1]) if (in_front(planes[i], p)) i]; + [for (i = [0:1:len(planes) - 1]) if (in_front(planes[i], p)) i]; function find_first_noncollinear(line, pts, i) = (i >= len(pts)) ? len(pts) diff --git a/linalg.scad b/linalg.scad index 90a0d7f..955dbae 100644 --- a/linalg.scad +++ b/linalg.scad @@ -86,4 +86,4 @@ function construct_Rt(R, t) = // --- Elementwise Operations ------------------------------------------------- // Hadamard product: works recursively on arrays, multiplies scalars directly. function hadamard(a, b) = - is_list(a) ? [for (i = [0:len(a) - 1]) hadamard(a[i], b[i])] : a * b; + is_list(a) ? [for (i = [0:1:len(a) - 1]) hadamard(a[i], b[i])] : a * b; diff --git a/lists.scad b/lists.scad index cad7650..5ce01a5 100644 --- a/lists.scad +++ b/lists.scad @@ -58,7 +58,7 @@ function subarray(list, begin = 0, end = -1) = set([1,2,3,4], 2, 5) => [1,2,5,4] */ function set(list, i, x) = - [for (j = [0:len(list) - 1]) (i == j) ? x : list[j]]; + [for (j = [0:1:len(list) - 1]) (i == j) ? x : list[j]]; // --- Remove --------------------------------------------------------------- /*! @@ -68,4 +68,4 @@ function set(list, i, x) = remove([4,3,2,1], 1) => [4,2,1] */ function remove(list, i) = - [for (j = [0:len(list) - 2]) list[ (j < i) ? j : j + 1]]; + [for (j = [0:1:len(list) - 2]) list[ (j < i) ? j : j + 1]]; From 5841d93e0275f77d24e07dcb51a44e1a9a6fa131 Mon Sep 17 00:00:00 2001 From: Cameron K Brooks Date: Fri, 26 Sep 2025 21:59:05 -0400 Subject: [PATCH 22/32] revise spline cosmetrically and make tests --- spline.scad | 167 ++++++++++++++++++++++++++--------------- tests/test_spline.scad | 91 ++++++++++++++++++++++ 2 files changed, 199 insertions(+), 59 deletions(-) create mode 100644 tests/test_spline.scad diff --git a/spline.scad b/spline.scad index 384305c..20d960c 100644 --- a/spline.scad +++ b/spline.scad @@ -1,124 +1,173 @@ -// Spline module for scad-util library -// Author Sergei Kuzmin, 2014. - -// For n+1 given point and hense n intervals returns the spline coefficient matrix. -// param p defines the anchor points. -// File defines two functions: spline_args and spline. -// example usage: -// spl1 = spline_args(point, v1=[0,1,0], closed=false); -// interpolated_points = [for(t=[0:0.1:len(point)-1]) spline(spl1, t)] +// ============================================================================ +// Spline Utilities +// ---------------------------------------------------------------------------- +// Author: Sergei Kuzmin, 2014 +// License: BSD +// +// Provides cubic spline and Bezier curve utilities for OpenSCAD. +// +// Functions: +// - spline_args(p, closed=false, v1=undef, v2=undef) +// → compute spline coefficient matrices for given control points +// - spline(s, t) +// → evaluate spline at parameter t +// - spline_tan(), spline_tan_unit() +// → evaluate tangent vector +// - spline_d2(), spline_binormal_unit(), spline_normal_unit() +// → evaluate higher derivatives and Frenet frame +// - spline_transform() +// → construct SE3 transform aligned to spline +// +// - bezier3_args(p, symmetric=false) +// → construct cubic Bezier coefficient matrices +// +// Internal helpers implement determinant, matrix inversion, etc. +// ============================================================================ use use +// --- Predefined Matrices ---------------------------------------------------- q1 = [[1, 0, 0, 0], [1, 1, 1, 1], [0, 1, 2, 3], [0, 0, 1, 3]]; q1inv = [[1, 0, 0, 0], [-3, 3, -2, 1], [3, -3, 3, -2], [-1, 1, -1, 1]]; q2 = [[0, 0, 0, 0], [0, 0, 0, 0], [0, -1, 0, 0], [0, 0, -1, 0]]; qn1i2 = -q1inv * q2; + z3 = [0, 0, 0]; z4 = [0, 0, 0, 0]; +// --- Matrix Utilities ------------------------------------------------------- + +// Compute matrix power function matrix_power(m, n) = n == 0 ? (len(m) == 3 ? identity3() : identity4()) - : n == 1 ? m : (n % 2 == 1) ? matrix_power(m * m, floor(n / 2)) * m : matrix_power(m * m, n / 2); + : n == 1 ? m + : (n % 2 == 1) ? matrix_power(m * m, floor(n / 2)) * m + : matrix_power(m * m, n / 2); + +// Determinant (recursive Laplace expansion) +function det(m) = + let (r = [for (i = [0:len(m) - 1]) i]) det_help(m, 0, r); -function det(m) = let (r = [for (i = [0:1:len(m) - 1]) i]) det_help(m, 0, r); // Construction indices list is inefficient, but currently there is no way to imperatively // assign to a list element function det_help(m, i, r) = len(r) == 0 ? 1 : m[len(m) - len(r)][r[i]] * det_help(m, 0, remove(r, i)) - (i + 1 < len(r) ? det_help(m, i + 1, r) : 0); +// Matrix inversion (adjugate method) function matrix_invert(m) = let (r = [for (i = [0:len(m) - 1]) i]) [ for (i = r) [ for (j = r) ( (i + j) % 2 == 0 ? 1 : -1) * matrix_minor(m, 0, remove(r, j), remove(r, i)), ], ] / det(m); + function matrix_minor(m, k, ri, rj) = let (len_r = len(ri)) len_r == 0 ? 1 : m[ri[0]][rj[k]] * matrix_minor(m, 0, remove(ri, 0), remove(rj, k)) - (k + 1 < len_r ? matrix_minor(m, k + 1, ri, rj) : 0); +// --- Spline Construction ---------------------------------------------------- + +// Local utility function spline_u(i, p) = [p[i], p[i + 1], z3, z3]; +// Compute spline coefficient matrices function spline_args(p, closed = false, v1 = undef, v2 = undef) = len(p) < 2 ? [] : let ( - q3 = closed ? q2 : [z4, z4, v1 == undef ? [0, 0, 1, 0] : [0, 1, 0, 0], z4], - q4 = closed ? q1 : [[1, 0, 0, 0], [1, 1, 1, 1], z4, v2 == undef ? [0, 0, 1, 3] : [0, 1, 2, 3]], + q3 = closed ? q2 + : [ + z4, + z4, + v1 == undef ? [0, 0, 1, 0] : [0, 1, 0, 0], + z4, + ], + q4 = closed ? q1 + : [ + [1, 0, 0, 0], + [1, 1, 1, 1], + z4, + v2 == undef ? [0, 0, 1, 3] : [0, 1, 2, 3], + ], pcnt = closed ? len(p) + 1 : len(p), - un = [p[pcnt - 2], p[closed ? 0 : pcnt - 1], v1 == undef ? z4 : v1, v2 == undef ? z4 : v2], + un = [ + p[pcnt - 2], + p[closed ? 0 : pcnt - 1], + v1 == undef ? z4 : v1, + v2 == undef ? z4 : v2, + ], sn = matrix_invert(q4 + q3 * matrix_power(qn1i2, pcnt - 2)) * (un - q3 * q1inv * spline_helper(0, pcnt, p)) ) // result[i+1] recurrently defines result[i]. This is O(n) runtime with imperative language and // may be O(n^2) if OpenSCAD doesn't cache spline_si(i+1). [for (i = [0:pcnt - 2]) spline_si(i, pcnt - 2, p, sn)]; +// Helper recursion for spline_args // n is number of points including pseudopoint for closed contour // Weird construct cause there is no if statement for functions -function spline_helper(i, n, p) = let (u = [p[i], p[i + 1], z3, z3]) i + 3 >= n ? u : u - q2 * q1inv * spline_helper(i + 1, n, p); +function spline_helper(i, n, p) = + let (u = [p[i], p[i + 1], z3, z3]) i + 3 >= n ? u + : u - q2 * q1inv * spline_helper(i + 1, n, p); +// Recursive calculation of segment coefficients // knowing s[j+1], calculate s[j]. Stop when found s[i] -function spline_si(i, n, p, sn) = i == n ? sn : q1inv * (spline_u(i, p) - q2 * spline_si(i + 1, n, p, sn)); +function spline_si(i, n, p, sn) = + i == n ? sn + : q1inv * (spline_u(i, p) - q2 * spline_si(i + 1, n, p, sn)); + +// --- Bezier Utilities ------------------------------------------------------- // Takes array of (3n+1) points or (2n + 2) points, if tangent segments are symmetric. // For non-symmetric version input is: point0, normal0, neg_normal1, point1, normal1, ... neg_normal_n, point_n // For symmetric version: point0, normal0, point1, normal1, ... , normal_n_sub_1, point_n // In the second case second tangent is constructed from the next tangent by symmetric map. // I.e. if current points are p0,p1,p2 then anchor points are p0 and p2, first tangent defined by p1-p0, // second tangent defined by p3-p2. -// Return array of coefficients accepted by spline(), spline_tan() and similar +// Return array of coefficients accepted by spline(), spline_tan() and similar function bezier3_args(p, symmetric = false) = let (step = symmetric ? 2 : 3) [ - for (i = [0:step:len(p) - 3]) [[1, 0, 0, 0], [-3, 3, 0, 0], [3, -6, 3, 0], [-1, 3, -3, 1]] * (symmetric ? [p[i], p[i] + p[i + 1], p[i + 2] - p[i + 3], p[i + 2]] : [p[i], p[i] + p[i + 1], p[i + 3] + p[i + 2], p[i + 3]]), + for (i = [0:step:len(p) - 3]) [[1, 0, 0, 0], [-3, 3, 0, 0], [3, -6, 3, 0], [-1, 3, -3, 1]] * ( + symmetric ? [p[i], p[i] + p[i + 1], p[i + 2] - p[i + 3], p[i + 2]] + : [p[i], p[i] + p[i + 1], p[i + 3] + p[i + 2], p[i + 3]] + ), ]; +// --- Evaluation ------------------------------------------------------------- +// Evaluate spline/Bezier and derivatives // s - spline arguments calculated by spline_args // t - defines point on curve. each segment length is 1. I.e. t= 0..1 is first segment, t=1..2 - second. -function spline(s, t) = let (i = t >= len(s) ? len(s) - 1 : floor(t), t2 = t - i) [1, t2, t2 * t2, t2 * t2 * t2] * s[i]; +function spline(s, t) = + let ( + i = t >= len(s) ? len(s) - 1 : floor(t), + t2 = t - i + ) [1, t2, t2 * t2, t2 * t2 * t2] * s[i]; + +function spline_tan(s, t) = + let ( + i = t >= len(s) ? len(s) - 1 : floor(t), + t2 = t - i + ) [0, 1, 2 * t2, 3 * t2 * t2] * s[i]; -function spline_tan(s, t) = let (i = t >= len(s) ? len(s) - 1 : floor(t), t2 = t - i) [0, 1, 2 * t2, 3 * t2 * t2] * s[i]; function spline_tan_unit(s, t) = unit(spline_tan(s, t)); -function spline_d2(s, t) = let (i = t >= len(s) ? len(s) - 1 : floor(t), t2 = t - i) [0, 0, 2, 6 * t2] * s[i]; + +function spline_d2(s, t) = + let ( + i = t >= len(s) ? len(s) - 1 : floor(t), + t2 = t - i + ) [0, 0, 2, 6 * t2] * s[i]; + function spline_binormal_unit(s, t) = unit(cross(spline_tan(s, t), spline_d2(s, t))); function spline_normal_unit(s, t) = unit(cross(spline_tan(s, t), spline_binormal_unit(s, t))); function spline_transform(s, t) = - construct_Rt(transpose_3([spline_normal_unit(s, t), spline_binormal_unit(s, t), spline_tan_unit(s, t)]), spline(s, t)); - -// Unit tests -__s = spline_args([[0, 10, 0], [10, 0, 0], [0, -5, 2]], v1=[0, 1, 0], v2=[-1, 0, 0], closed=true); -for (t = [0:0.01:len(__s)]) - translate(spline(__s, t)) - cube([0.2, 0.2, 0.2], center=true); - -__s1 = spline_args([[0, 0, 0], [0, 0, 15], [26, 0, 26 + 15]], /*v1=[0,0,100],*/ v2=[40, 0, 0]); -for (t = [0:0.01:len(s1)]) - translate(spline(__s1, t)) - cube([0.2, 0.2, 0.2], center=true); - -__s2 = bezier3_args([[0, 0, 0], [0, 0, 10], [0, 0, 15], [0, 0, 26 * 0.552284], [26, 0, 41], [26 * 0.552284, 0, 0]], symmetric=true); -echo(__s2); -for (t = [0:0.01:len(__s2)]) - translate(spline(__s2, t)) - cube([0.2, 0.2, 0.2], center=true); - -// Rotation methods taken from list-comprehension-demos/sweep.scad to demonstrate normal and binormal -// Normally spline_transform is more convenient -function __rotation_from_axis(x, y, z) = [[x[0], y[0], z[0]], [x[1], y[1], z[1]], [x[2], y[2], z[2]]]; -function __rotate_from_to(a, b, _axis = []) = - len(_axis) == 0 ? __rotate_from_to(a, b, unit(cross(a, b))) - : _axis * _axis >= 0.99 ? __rotation_from_axis(unit(b), _axis, cross(_axis, unit(b))) * transpose_3(__rotation_from_axis(unit(a), _axis, cross(_axis, unit(a)))) - : identity3(); - -__s3 = spline_args([[0, 10, 0], [6, 6, 0], [10, 0, 0], [0, -5, 4]], v1=[0, 1, 0], v2=[-1, 0, 0], closed=true); -for (t = [0:0.05:len(__s3)]) - translate(spline(__s3, t)) { - translate([0, 0, 3]) multmatrix(m=__rotate_from_to([0, 0, 1], spline_normal_unit(__s3, t))) - cylinder(r1=0.1, r2=0, h=1, $fn=3); - translate([0, 0, 6]) multmatrix(m=__rotate_from_to([0, 0, 1], spline_binormal_unit(__s3, t))) - cylinder(r1=0.1, r2=0, h=1, $fn=3); - } - -translate([0, 0, 9])for (t = [0:0.025:len(__s3)]) - multmatrix(spline_transform(__s3, t)) cube([1, 1, 0.1], center=true); + construct_Rt( + transpose_3( + [ + spline_normal_unit(s, t), + spline_binormal_unit(s, t), + spline_tan_unit(s, t), + ] + ), + spline(s, t) + ); diff --git a/tests/test_spline.scad b/tests/test_spline.scad new file mode 100644 index 0000000..211d660 --- /dev/null +++ b/tests/test_spline.scad @@ -0,0 +1,91 @@ +include <../spline.scad> +include <../se3.scad> + +$fn = 24; + +// ============================================================================ +// Test Suite for spline.scad +// ---------------------------------------------------------------------------- +// Covers: +// - spline_args with open/closed curves +// - tangent / binormal / normal evaluation +// - Bezier curve generation +// - Frenet frame transform visualization +// - Regression consistency checks +// ============================================================================ + +// --- Basic Spline Example (closed curve) ----------------------------------- +p1 = [[0, 10, 0], [10, 0, 0], [0, -5, 2]]; +s1 = spline_args(p1, v1=[0, 1, 0], v2=[-1, 0, 0], closed=true); + +// Visualize spline points +for (t = [0:0.01:len(s1)]) + translate(spline(s1, t)) + color("red") sphere(r=0.1); + +// Tangent demo +for (t = [0:0.5:len(s1)]) + let (pt = spline(s1, t)) + color("blue") + translate(pt) cylinder(r=0.1, h=2, center=true, $fn=12); + +// --- Open Spline Example --------------------------------------------------- +p2 = [[0, 0, 0], [0, 0, 15], [26, 0, 41]]; +s2 = spline_args(p2, v2=[40, 0, 0]); + +for (t = [0:0.01:len(s2)]) + translate(spline(s2, t)) + color("indigo") sphere(r=0.1); + +// --- Bezier Curve Example -------------------------------------------------- +p3 = [ + [0, 0, 0], + [0, 0, 10], + [0, 0, 15], + [0, 0, 26 * 0.552], + [26, 0, 41], + [26 * 0.552, 0, 0], +]; +s3 = bezier3_args(p3, symmetric=true); + +echo("Bezier coefficients =", s3); + +for (t = [0:0.01:len(s3)]) + translate(spline(s3, t)) + color("green") sphere(r=0.1); + +// --- Frenet Frame Demo ----------------------------------------------------- +// Rotation methods (ported from list-comprehension-demos/sweep.scad) +function __rotation_from_axis(x, y, z) = + [[x[0], y[0], z[0]], [x[1], y[1], z[1]], [x[2], y[2], z[2]]]; + +function __rotate_from_to(a, b, _axis = []) = + len(_axis) == 0 ? __rotate_from_to(a, b, unit(cross(a, b))) + : _axis * _axis >= 0.99 ? __rotation_from_axis(unit(b), _axis, cross(_axis, unit(b))) * transpose_3(__rotation_from_axis(unit(a), _axis, cross(_axis, unit(a)))) + : identity3(); + +p4 = [[0, 10, 0], [6, 6, 0], [10, 0, 0], [0, -5, 4]]; +s4 = spline_args(p4, v1=[0, 1, 0], v2=[-1, 0, 0], closed=true); + +// Normal/binormal markers +for (t = [0:0.05:len(s4)]) + translate(spline(s4, t)) { + translate([0, 0, 3]) + multmatrix(m=__rotate_from_to([0, 0, 1], spline_normal_unit(s4, t))) + color("teal") + cylinder(r1=0.1, r2=0, h=1, $fn=3); + translate([0, 0, 6]) + multmatrix(m=__rotate_from_to([0, 0, 1], spline_binormal_unit(s4, t))) + color("brown") + cylinder(r1=0.1, r2=0, h=1, $fn=3); + } + +// Spline transform demo (Frenet-aligned cubes) +translate([0, 0, 9])for (t = [0:0.025:len(s4)]) + multmatrix(spline_transform(s4, t)) + color("orange") + cube([1, 1, 0.1], center=true); + +// --- Consistency Check ----------------------------------------------------- +__test = [20, -40, 60, -80, 100, -120]; +echo("SE3 consistency =", norm(__test - se3_ln(se3_exp(__test))) < 1e-8); From e40711fce1df9bfd97a91e48e3b15da261c9948c Mon Sep 17 00:00:00 2001 From: Cameron K Brooks Date: Fri, 26 Sep 2025 22:35:42 -0400 Subject: [PATCH 23/32] revise trajectory path to use is_undef and cosmetic revisio; add test file --- tests/test_trajectory_path.scad | 56 +++++++++++ trajectory_path.scad | 165 ++++++++++++++++++++++++-------- 2 files changed, 180 insertions(+), 41 deletions(-) diff --git a/tests/test_trajectory_path.scad b/tests/test_trajectory_path.scad index e69de29..6d78caa 100644 --- a/tests/test_trajectory_path.scad +++ b/tests/test_trajectory_path.scad @@ -0,0 +1,56 @@ +include <../trajectory_path.scad> +include <../linalg.scad> +include <../se3.scad> + +$fn = 24; + +// Simple frame visualizer ---------------------------------------------------- +module frame(T, s = 2) { + multmatrix(T) { + color("red") cube([s, .2, .2], center=false); + color("green") cube([.2, s, .2], center=false); + color("blue") cube([.2, .2, s], center=false); + color("gray") translate([0, 0, 0]) sphere(r=.4); + } +} + +// Example trajectories (6D twists) ------------------------------------------ +// Move +X 60, slight +Z arc; then yaw+translate; then small upward arc. +traj = [ + [60, 0, 0, 0, 0, 30], // translate 60 along X while yawing 30° + [20, 0, 0, 0, 0, -45], // short turn back + [0, 0, 30, 0, 90, 0], // arc up in pitch + [0, 0, 20, 0, 0, 0], // short straight up (relative forward to base frame) +]; + +// --- Sampling with a physical step (units of translation norm) ------------- +poses_step = quantize_trajectories(traj, step=5, start_position=0, loop=false); + +// --- Sampling with a fixed number of steps across whole path ---------------- +poses_steps = quantize_trajectories(traj, steps=30, start_position=0, loop=false); + +// --- Looping path (returns to start) --------------------------------------- +poses_loop = quantize_trajectories(traj, step=5, loop=true); + +// Render -------------------------------------------------------------------- +translate([0, -70, 0]) { + // By fixed step length + for (T = poses_step) frame(T, s=3); + echo("poses_step count =", len(poses_step)); +} + +translate([0, 0, 0]) { + // By fixed number of steps + for (T = poses_steps) frame(T, s=3); + echo("poses_steps count =", len(poses_steps)); +} + +translate([0, 70, 0]) { + // Looping variant + for (T = poses_loop) frame(T, s=3); + echo("poses_loop count =", len(poses_loop)); +} + +// Show overall end pose to verify loop closure visually +echo("End pose (open) =", trajectories_end_position(traj)); +echo("End pose (loop) =", trajectories_end_position(close_trajectory_loop(traj))); diff --git a/trajectory_path.scad b/trajectory_path.scad index 905e9a3..d28e345 100644 --- a/trajectory_path.scad +++ b/trajectory_path.scad @@ -1,47 +1,129 @@ +// ============================================================================ +// Trajectory Path Utilities +// ---------------------------------------------------------------------------- +// Purpose : Quantize SE(3) “twists” (6D vectors) into discrete transforms, +// compose segments, and optionally close loops. +// Key ideas +// - A single trajectory is a 6D vector mu = [tx,ty,tz, wx,wy,wz] (angles in deg) +// interpreted via se3_exp(mu) -> 4x4 homogeneous transform. +// - A list of trajectories is concatenated piecewise in the given order. +// - Quantization can be specified either by: +// * step: physical step length along translation part, OR +// * steps: fixed number of samples for the *entire* path +// - start_position allows skipping an initial arc-length offset. +// ---------------------------------------------------------------------------- +// Dependencies: linalg.scad, se3.scad +// ============================================================================ + use use +// --- Small helpers ---------------------------------------------------------- + +// Left-multiply a single transform `a` by each transform in array `bs`. function left_multiply(a, bs, i_ = 0) = - i_ >= len(bs) ? [] - : concat( - [ - a * bs[i_], - ], left_multiply(a, bs, i_ + 1) - ); + (i_ >= len(bs)) ? [] + : concat([a * bs[i_]], left_multiply(a, bs, i_ + 1)); +// Right-multiply each transform in array `as` by a single transform `b`. function right_multiply(as, b, i_ = 0) = - i_ >= len(as) ? [] - : concat( - [ - as[i_] * b, - ], right_multiply(as, b, i_ + 1) - ); + (i_ >= len(as)) ? [] + : concat([as[i_] * b], right_multiply(as, b, i_ + 1)); -function quantize_trajectory(trajectory, step = undef, start_position = 0, steps = undef, i_ = 0, length_ = undef) = - length_ == undef ? quantize_trajectory( +// --- Single-trajectory quantization ---------------------------------------- +// Produces an array of 4x4 transforms sampled along one 6D twist `trajectory`. +// If `steps` is provided, it overrides `step` and distributes uniformly. +// If `step` is provided, samples start at `start_position` and advance by `step`. +// `start_position` and `step` are in the units of the translation norm. +// +// Notes: +// - Length is computed from the translational part only: norm(take3(trajectory)) +// - For steps==1, returns the transform at the path end (t=1). +function quantize_trajectory( + trajectory, + step = undef, + start_position = 0, + steps = undef, + i_ = 0, + length_ = undef +) = + // Bootstrap total length if not provided + is_undef(length_) ? + quantize_trajectory( trajectory=trajectory, - start_position=(step == undef ? norm(take3(trajectory)) / steps * start_position : start_position), + start_position=is_undef(step) ? (norm(take3(trajectory)) / steps) * start_position + : start_position, length_=norm(take3(trajectory)), - step=step, steps=steps, i_=i_ + step=step, + steps=steps, + i_=i_ ) - : (steps == undef ? start_position > length_ : i_ >= steps) ? [] + // Termination: either finished by count, or past end by distance + : ( + is_undef(steps) ? (start_position > length_) + : (i_ >= steps) + ) ? [] + // Emit current sample and recurse : concat( [ - // if steps is defined, ignore start_position se3_exp( trajectory * ( - steps == undef ? start_position / length_ - : i_ / (steps > 1 ? steps - 1 : 1) + is_undef(steps) ? start_position / length_ + : i_ / (steps > 1 ? (steps - 1) : 1) ) ), - ], quantize_trajectory(trajectory=trajectory, step=step, start_position=(steps == undef ? start_position + step : start_position), steps=steps, i_=i_ + 1, length_=length_) + ], + quantize_trajectory( + trajectory=trajectory, + step=step, + start_position=is_undef(steps) ? (start_position + step) : start_position, + steps=steps, + i_=i_ + 1, + length_=length_ + ) ); -function close_trajectory_loop(trajectories) = concat(trajectories, [se3_ln(invert_rt(trajectories_end_position(trajectories)))]); +// --- Multi-trajectory helpers ---------------------------------------------- -function quantize_trajectories(trajectories, step = undef, start_position = 0, steps = undef, loop = false, last_ = identity4(), i_ = 0, current_length_ = undef, j_ = 0) = - // due to quantization differences, the last step may be missed. In that case, add it: - loop == true ? quantize_trajectories( +// Append a final segment that closes the loop back to identity when applied +// after the whole chain (computed in SE(3) log space). +function close_trajectory_loop(trajectories) = + concat( + trajectories, + [se3_ln(invert_rt(trajectories_end_position(trajectories)))] + ); + +// Total translational arc-length of a list of 6D twists. +function trajectories_length(trajectories, i_ = 0) = + (i_ >= len(trajectories)) ? 0 + : norm(take3(trajectories[i_])) + trajectories_length(trajectories, i_ + 1); + +// End pose after chaining all twists (via se3_exp and left-to-right product). +function trajectories_end_position(rt, i_ = 0, last_ = identity4()) = + (i_ >= len(rt)) ? last_ + : trajectories_end_position(rt, i_ + 1, last_ * se3_exp(rt[i_])); + +// --- Multi-trajectory quantization ----------------------------------------- +// Quantize a *list* of twists into an array of transforms along the whole path. +// Parameters mirror quantize_trajectory; additionally: +// - loop=true closes the path with an extra segment that returns to start +// Behavior: +// - If `steps` is given, we derive a uniform `step` for the *whole* path. +// - Returns an array of transforms (absolute poses), starting from identity. +function quantize_trajectories( + trajectories, + step = undef, + start_position = 0, + steps = undef, + loop = false, + last_ = identity4(), + i_ = 0, + current_length_ = undef, + j_ = 0 +) = + // Optionally close the loop by appending a final corrective twist + loop ? + quantize_trajectories( trajectories=close_trajectory_loop(trajectories), step=step, start_position=start_position, @@ -52,12 +134,16 @@ function quantize_trajectories(trajectories, step = undef, start_position = 0, s current_length_=current_length_, j_=j_ ) - : i_ >= len(trajectories) ? (j_ < steps ? [last_] : []) - : current_length_ == undef ? + // End when all segments consumed. If steps was given but rounding skipped + // the final sample, return the last pose once more. + : (i_ >= len(trajectories)) ? ( (!is_undef(steps) && (j_ < steps)) ? [last_] : []) + // Initialize global step/offset once we know total length + : is_undef(current_length_) ? quantize_trajectories( trajectories=trajectories, - step=(step == undef ? trajectories_length(trajectories) / steps : step), - start_position=(step == undef ? start_position * trajectories_length(trajectories) / steps : start_position), + step=is_undef(step) ? (trajectories_length(trajectories) / steps) : step, + start_position=is_undef(step) ? (start_position * trajectories_length(trajectories) / steps) + : start_position, steps=steps, loop=loop, last_=last_, @@ -65,26 +151,31 @@ function quantize_trajectories(trajectories, step = undef, start_position = 0, s current_length_=norm(take3(trajectories[i_])), j_=j_ ) + // Emit samples for current segment (absolute poses), then continue : concat( + // Sample the i_-th local twist, then left-multiply by `last_` to get + // absolute poses for this segment. left_multiply( - last_, quantize_trajectory( + last_, + quantize_trajectory( trajectory=trajectories[i_], start_position=start_position, step=step ) ), + // Recurse to next segment: quantize_trajectories( trajectories=trajectories, step=step, - start_position=start_position > current_length_ ? start_position - current_length_ - : step - ( (current_length_ - start_position) % step), + // Advance start_position into the next segment: + start_position=(start_position > current_length_) ? (start_position - current_length_) + : (step - ( (current_length_ - start_position) % step)), steps=steps, loop=loop, last_=last_ * se3_exp(trajectories[i_]), i_=i_ + 1, current_length_=undef, j_=j_ + len( - quantize_trajectory( trajectory=trajectories[i_], start_position=start_position, @@ -93,11 +184,3 @@ function quantize_trajectories(trajectories, step = undef, start_position = 0, s ) ) ); - -function trajectories_length(trajectories, i_ = 0) = - i_ >= len(trajectories) ? 0 - : norm(take3(trajectories[i_])) + trajectories_length(trajectories, i_ + 1); - -function trajectories_end_position(rt, i_ = 0, last_ = identity4()) = - i_ >= len(rt) ? last_ - : trajectories_end_position(rt, i_ + 1, last_ * se3_exp(rt[i_])); From 74091b3faa432ce33ce3ef5a0a8e2bb5e4796cfa Mon Sep 17 00:00:00 2001 From: Cameron K Brooks Date: Fri, 26 Sep 2025 22:44:20 -0400 Subject: [PATCH 24/32] revise trajectory to use is_undef and cosmetic revision and add test file --- tests/test_trajectory.scad | 36 ++++++++++++++++ trajectory.scad | 88 ++++++++++++++++++++++++++++++++------ 2 files changed, 111 insertions(+), 13 deletions(-) diff --git a/tests/test_trajectory.scad b/tests/test_trajectory.scad index e69de29..172638e 100644 --- a/tests/test_trajectory.scad +++ b/tests/test_trajectory.scad @@ -0,0 +1,36 @@ +include <../trajectory.scad> +include <../so3.scad> + +$fn = 24; + +// --- Regression / Echo Tests ------------------------------------------------ +echo("Forward 10 =", trajectory(forward=10)); // expect [0,0,10, 0,0,0] +echo("Up 5, Left 3 =", trajectory(up=5, left=3)); // expect [-3,5,0, 0,0,0] +echo("Yaw 45 =", rotationv(yaw=45)); // expect [45,0,0] +echo("Pitch 90 matrix =", rotationm(pitch=90)); + +// Combined twist +twist1 = trajectory(forward=20, yaw=45); +echo("Combined twist =", twist1); + +// --- Visual Demonstrations -------------------------------------------------- + +// Helper: apply trajectory to a cube +module demo_traj(traj, col = "lightblue") { + T = concat(take3(traj), [0]); // translation + R = rotationm(rotation=tail3(traj)); + multmatrix( + [ + [R[0][0], R[0][1], R[0][2], T[0]], + [R[1][0], R[1][1], R[1][2], T[1]], + [R[2][0], R[2][1], R[2][2], T[2]], + [0, 0, 0, 1], + ] + ) + color(col) cube([5, 5, 5], center=true); +} + +// Show original then two trajectory examples side by side +translate([0, 0, 0]) demo_traj(trajectory(up=0), "lightblue"); +translate([-20, 0, 0]) demo_traj(trajectory(right=10, up=5, yaw=30), "red"); +translate([20, 0, 0]) demo_traj(trajectory(translation=[5, 5, 5], rotation=[0, 90, 0]), "green"); diff --git a/trajectory.scad b/trajectory.scad index af059f0..535d64a 100644 --- a/trajectory.scad +++ b/trajectory.scad @@ -1,32 +1,88 @@ +// ============================================================================ +// Trajectory Utilities +// ---------------------------------------------------------------------------- +// Purpose : Provide a clean way to construct SE(3) “twist” vectors +// (translation + rotation) from intuitive parameters. +// Features: +// - Flexible specification of translation (left/right, up/down, fwd/back) +// - Flexible specification of rotation (pitch, yaw, roll) or direct vector +// - Helpers for dealing with undef inputs and selection logic +// - trajectory(...) returns [tx,ty,tz, yaw,pitch,roll] 6D vector +// - rotationm(...) builds a 3x3 rotation matrix from angles +// ---------------------------------------------------------------------------- +// Dependencies: so3.scad +// Author(s): adapted and cleaned up from original scad-utils +// ============================================================================ + use -function val(a = undef, default = undef) = a == undef ? default : a; +// --- Helpers ---------------------------------------------------------------- + +// Return `a` unless it is undef, then use `default`. +function val(a = undef, default = undef) = + is_undef(a) ? default : a; + +// Check if vector/list is fully undef or out-of-bounds (recursively). function vec_is_undef(x, index_ = 0) = index_ >= len(x) ? true : is_undef_or_oob(x[index_]) && vec_is_undef(x, index_ + 1); -function is_undef_or_oob(x) = is_undef(x) ? true : is_list(x) ? vec_is_undef(x) : false; -// Either a or b, but not both -function either(a, b, default = undef) = is_undef_or_oob(a) ? (is_undef_or_oob(b) ? default : b) : is_undef_or_oob(b) ? a : undef; +// Treat scalars, lists, or undef consistently. +function is_undef_or_oob(x) = + is_undef(x) ? true + : is_list(x) ? vec_is_undef(x) + : false; + +// Either return a or b (whichever is valid). If both valid → undef. +// If both invalid → default. +function either(a, b, default = undef) = + is_undef_or_oob(a) ? (is_undef_or_oob(b) ? default : b) + : is_undef_or_oob(b) ? a + : undef; + +// --- Translation builders --------------------------------------------------- -function translationv(left = undef, right = undef, up = undef, down = undef, forward = undef, backward = undef, translation = undef) = +// Accepts directional keywords (left/right, up/down, forward/backward) +// or an explicit `translation` vector. Returns [x,y,z] or undef. +function translationv( + left = undef, + right = undef, + up = undef, + down = undef, + forward = undef, + backward = undef, + translation = undef +) = translationv_2( + // X axis points "up" (OpenSCAD’s Z is height, but we’re flexible here) x=either(up, is_undef(down) ? down : -down), y=either(right, is_undef(left) ? left : -left), z=either(forward, is_undef(backward) ? backward : -backward), translation=translation ); +// Internal helper to finalize translation vector. function translationv_2(x, y, z, translation) = - x == undef && y == undef && z == undef ? translation - : is_undef_or_oob(translation) ? [val(x, 0), val(y, 0), val(z, 0)] - : undef; + (is_undef(x) && is_undef(y) && is_undef(z)) ? translation + : (is_undef_or_oob(translation) ? [val(x, 0), val(y, 0), val(z, 0)] : undef); -function rotationv(pitch = undef, yaw = undef, roll = undef, rotation = undef) = - rotation == undef ? [val(yaw, 0), val(pitch, 0), val(roll, 0)] - : pitch == undef && yaw == undef && roll == undef ? rotation +// --- Rotation builders ------------------------------------------------------ + +// Accepts yaw/pitch/roll angles (deg), or an explicit `rotation` vector. +// Returns [yaw, pitch, roll] vector or undef. +function rotationv( + pitch = undef, + yaw = undef, + roll = undef, + rotation = undef +) = + is_undef(rotation) ? [val(yaw, 0), val(pitch, 0), val(roll, 0)] + : (is_undef(pitch) && is_undef(yaw) && is_undef(roll)) ? rotation : undef; +// --- Main API --------------------------------------------------------------- + +// trajectory(...) → build a 6D twist vector = [tx,ty,tz, yaw,pitch,roll] function trajectory( left = undef, right = undef, @@ -41,8 +97,14 @@ function trajectory( rotation = undef ) = concat( - translationv(left=left, right=right, up=up, down=down, forward=forward, backward=backward, translation=translation), + translationv( + left=left, right=right, up=up, down=down, + forward=forward, backward=backward, + translation=translation + ), rotationv(pitch=pitch, yaw=yaw, roll=roll, rotation=rotation) ); -function rotationm(rotation = undef, pitch = undef, yaw = undef, roll = undef) = so3_exp(rotationv(rotation=rotation, pitch=pitch, yaw=yaw, roll=roll)); +// rotationm(...) → build 3x3 rotation matrix from angles +function rotationm(rotation = undef, pitch = undef, yaw = undef, roll = undef) = + so3_exp(rotationv(rotation=rotation, pitch=pitch, yaw=yaw, roll=roll)); From 4ee523ec2bed5881c07e4e13dc85b09c2c120e42 Mon Sep 17 00:00:00 2001 From: Cameron K Brooks Date: Fri, 26 Sep 2025 22:46:33 -0400 Subject: [PATCH 25/32] add comment revisions to transformations and test file --- tests/test_transformations.scad | 40 ++++++++++++++++++++++++ trajectory.scad | 1 - transformations.scad | 55 ++++++++++++++++++++++++++++++--- 3 files changed, 90 insertions(+), 6 deletions(-) diff --git a/tests/test_transformations.scad b/tests/test_transformations.scad index e69de29..84748b7 100644 --- a/tests/test_transformations.scad +++ b/tests/test_transformations.scad @@ -0,0 +1,40 @@ +include <../transformations.scad> + +$fn = 24; + +// ============================================================================ +// Test Suite for transformations.scad +// ---------------------------------------------------------------------------- +// Covers: +// - rotation (Euler and axis-angle) +// - scaling +// - translation +// - project / transform / to_3d +// ============================================================================ + +// --- Sample Shape ----------------------------------------------------------- +points = [[0, 0, 0], [10, 0, 0], [10, 10, 0], [0, 10, 0]]; + +// --- Rotation Tests --------------------------------------------------------- +Rz45 = rotation(axis=[0, 0, 45]); +echo("Rotation Z45 applied =", transform(Rz45, points)); + +Rx90 = rotation(xyz=[90, 0, 0]); +echo("Rotation X90 applied =", transform(Rx90, [[0, 0, 1]])); + +// --- Scaling Tests ---------------------------------------------------------- +S = scaling([2, 1, 1]); +echo("Scaling [2,1,1] =", transform(S, points)); + +// --- Translation Tests ------------------------------------------------------ +T = translation([5, 5, 0]); +echo("Translation [5,5,0] =", transform(T, points)); + +// --- Combined Transform ----------------------------------------------------- +M = T * Rz45 * S; +echo("Combined TRS =", transform(M, points)); + +// --- Project / to_3d Tests -------------------------------------------------- +p_h = [3, 4, 5, 1]; +echo("Project [3,4,5,1] =", project(p_h)); +echo("to_3d([ [1,2], [3,4,5] ]) =", to_3d([[1, 2], [3, 4, 5]])); diff --git a/trajectory.scad b/trajectory.scad index 535d64a..37616c8 100644 --- a/trajectory.scad +++ b/trajectory.scad @@ -11,7 +11,6 @@ // - rotationm(...) builds a 3x3 rotation matrix from angles // ---------------------------------------------------------------------------- // Dependencies: so3.scad -// Author(s): adapted and cleaned up from original scad-utils // ============================================================================ use diff --git a/transformations.scad b/transformations.scad index 9659b60..bd1ff20 100644 --- a/transformations.scad +++ b/transformations.scad @@ -1,21 +1,51 @@ +// ============================================================================ +// Transformations Utilities +// ---------------------------------------------------------------------------- +// Provides matrix constructors and helpers for geometric transformations: +// - rotation(xyz, axis): create rotation from Euler angles or axis-angle +// - scaling(v): scaling matrix +// - translation(v): translation matrix +// - project(x): cartesian from homogeneous coordinates +// - transform(m, list): apply matrix to list of points +// - to_3d(list): ensure vectors are 3D +// +// Notes: +// - Euler convention: R = Rz * Ry * Rx +// - Axis rotation uses se3 exponential map +// ============================================================================ + use use use +// --- Rotation --------------------------------------------------------------- + /*! Creates a rotation matrix - xyz = euler angles = rz * ry * rx - axis = rotation_axis * rotation_angle + Options: + - xyz = Euler angles (applied as Rz * Ry * Rx) + - axis = axis-angle vector (axis * angle) + + Examples: + rotation(xyz=[90,0,0]) // rotate 90° about X + rotation(axis=[0,0,45]) // rotate 45° about Z */ function rotation(xyz = undef, axis = undef) = xyz != undef && axis != undef ? undef - : xyz == undef ? se3_exp([0, 0, 0, axis[0], axis[1], axis[2]]) - : len(xyz) == undef ? rotation(axis=[0, 0, xyz]) + : // cannot define both + xyz == undef ? + se3_exp([0, 0, 0, axis[0], axis[1], axis[2]]) + : len(xyz) == undef ? + rotation(axis=[0, 0, xyz]) : (len(xyz) >= 3 ? rotation(axis=[0, 0, xyz[2]]) : identity4()) * (len(xyz) >= 2 ? rotation(axis=[0, xyz[1], 0]) : identity4()) * (len(xyz) >= 1 ? rotation(axis=[xyz[0], 0, 0]) : identity4()); +// --- Scaling ---------------------------------------------------------------- + /*! Creates a scaling matrix + + scaling([sx, sy, sz]) */ function scaling(v) = [ @@ -25,8 +55,12 @@ function scaling(v) = [0, 0, 0, 1], ]; +// --- Translation ------------------------------------------------------------ + /*! Creates a translation matrix + + translation([tx, ty, tz]) */ function translation(v) = [ @@ -36,8 +70,19 @@ function translation(v) = [0, 0, 0, 1], ]; -// Convert between cartesian and homogenous coordinates +// --- Coordinate Conversion -------------------------------------------------- + +/*! + Converts from homogeneous to cartesian coordinates +*/ function project(x) = subarray(x, end=len(x) - 1) / x[len(x) - 1]; +/*! + Applies matrix `m` to a list of points +*/ function transform(m, list) = [for (p = list) project(m * vec4(p))]; + +/*! + Ensures points are represented as 3D vectors +*/ function to_3d(list) = [for (v = list) vec3(v)]; From f07f088fd35a2ebb4738309e94f50ef9da8484fe Mon Sep 17 00:00:00 2001 From: Cameron K Brooks Date: Sat, 27 Sep 2025 01:39:32 -0400 Subject: [PATCH 26/32] some final revisions to make clean before addressing readme and docs then PR --- linalg.scad | 3 +++ lists.scad | 55 ++++++++++++-------------------------------- se3.scad | 6 +++-- shapes.scad | 2 +- spline.scad | 54 ++++++++++++++++++++++--------------------- trajectory_path.scad | 6 ++--- transformations.scad | 46 ++++++++++-------------------------- 7 files changed, 66 insertions(+), 106 deletions(-) diff --git a/linalg.scad b/linalg.scad index 955dbae..4c22992 100644 --- a/linalg.scad +++ b/linalg.scad @@ -13,6 +13,7 @@ epsilon = 1e-9; // --- Vector Constructors ---------------------------------------------------- function vec3(p) = (len(p) < 3) ? concat(p, 0) : p; + function vec4(p) = let (v3 = vec3(p)) (len(v3) < 4) ? concat(v3, 1) : v3; @@ -36,6 +37,7 @@ function identity4() = // --- Vector Access Helpers -------------------------------------------------- function take3(v) = [v[0], v[1], v[2]]; + function tail3(v) = [v[3], v[4], v[5]]; // --- Matrix Part Extraction ------------------------------------------------- @@ -50,6 +52,7 @@ function translation_part(m) = [m[0][3], m[1][3], m[2][3]]; // --- Rotation Metrics ------------------------------------------------------- function rot_trace(m) = m[0][0] + m[1][1] + m[2][2]; + function rot_cos_angle(m) = (rot_trace(m) - 1) / 2; // --- Matrix Transpose ------------------------------------------------------- diff --git a/lists.scad b/lists.scad index 5ce01a5..ee3febd 100644 --- a/lists.scad +++ b/lists.scad @@ -9,63 +9,38 @@ // ============================================================================ // --- Flatten --------------------------------------------------------------- -/*! - Flatten a list one level. - - Example: - flatten([[0,1],[2,3]]) => [0,1,2,3] -*/ +// Flatten a list one level. +// Example: flatten([[0,1],[2,3]]) => [0,1,2,3] function flatten(list) = [for (sub = list, v = sub) v]; // --- Range ----------------------------------------------------------------- -/*! - Create a list from a range. - - Example: - range([0:2:6]) => [0,2,4,6] -*/ +// Create a list from a range. +// Example: range([0:2:6]) => [0,2,4,6] function range(r) = [for (x = r) x]; // --- Reverse --------------------------------------------------------------- -/*! - Reverse the order of elements. - - Example: - reverse([1,2,3]) => [3,2,1] -*/ +// Reverse the order of elements. +// Example: reverse([1,2,3]) => [3,2,1] function reverse(list) = [for (i = [len(list) - 1:-1:0]) list[i]]; // --- Subarray -------------------------------------------------------------- -/*! - Extract a subarray from index `begin` (inclusive) to `end` (exclusive). - - Notes: - - If `end < 0`, uses `len(list)`. - - FIXME: Consider renaming to `sublist` for clarity. - - Example: - subarray([1,2,3,4], 1, 3) => [2,3] -*/ +// Extract a subarray from index `begin` (inclusive) to `end` (exclusive). +// Notes: +// - If `end < 0`, uses `len(list)`. +// - FIXME: Consider renaming to `sublist` for clarity. +// Example: subarray([1,2,3,4], 1, 3) => [2,3] function subarray(list, begin = 0, end = -1) = let (end = (end < 0) ? len(list) : end) [for (i = [begin:end - 1]) list[i]]; // --- Set ------------------------------------------------------------------- -/*! - Return a copy of a list with the element at index `i` replaced by `x`. - - Example: - set([1,2,3,4], 2, 5) => [1,2,5,4] -*/ +// Return a copy of a list with the element at index `i` replaced by `x`. +// Example: set([1,2,3,4], 2, 5) => [1,2,5,4] function set(list, i, x) = [for (j = [0:1:len(list) - 1]) (i == j) ? x : list[j]]; // --- Remove --------------------------------------------------------------- -/*! - Remove the element at index `i`. - - Example: - remove([4,3,2,1], 1) => [4,2,1] -*/ +// Remove the element at index `i`. +// Example: remove([4,3,2,1], 1) => [4,2,1] function remove(list, i) = [for (j = [0:1:len(list) - 2]) list[ (j < i) ? j : j + 1]]; diff --git a/se3.scad b/se3.scad index 3d6f27b..1493a01 100644 --- a/se3.scad +++ b/se3.scad @@ -38,7 +38,8 @@ function se3_exp_2_0(t, w, theta_sq) = se3_exp_23( so3_exp_2(theta_sq), C=(1.0 - theta_sq / 20) / 6, - t=t, w=w + t=t, + w=w ); // 3rd order approximation @@ -53,7 +54,8 @@ function se3_exp_3_0(t, w, theta_deg, inv_theta) = se3_exp_23( so3_exp_3_0(theta_deg=theta_deg, inv_theta=inv_theta), C=(1 - sin(theta_deg) * inv_theta) * (inv_theta * inv_theta), - t=t, w=w + t=t, + w=w ); // Shared expansion for orders 2–3 diff --git a/shapes.scad b/shapes.scad index 302d425..07d19e3 100644 --- a/shapes.scad +++ b/shapes.scad @@ -40,4 +40,4 @@ function rectangle_profile(size = [1, 1]) = [size[0] / 2, -size[1] / 2], ]; -// TODO: Move rectangle and rounded rectangle from extrusion.scad +// FIXME: Move rectangle and rounded rectangle from extrusion diff --git a/spline.scad b/spline.scad index 20d960c..8912f21 100644 --- a/spline.scad +++ b/spline.scad @@ -1,27 +1,26 @@ // ============================================================================ // Spline Utilities // ---------------------------------------------------------------------------- -// Author: Sergei Kuzmin, 2014 -// License: BSD -// // Provides cubic spline and Bezier curve utilities for OpenSCAD. // // Functions: -// - spline_args(p, closed=false, v1=undef, v2=undef) -// → compute spline coefficient matrices for given control points -// - spline(s, t) -// → evaluate spline at parameter t -// - spline_tan(), spline_tan_unit() -// → evaluate tangent vector -// - spline_d2(), spline_binormal_unit(), spline_normal_unit() -// → evaluate higher derivatives and Frenet frame -// - spline_transform() -// → construct SE3 transform aligned to spline -// -// - bezier3_args(p, symmetric=false) -// → construct cubic Bezier coefficient matrices +// - spline_args(p, closed=false, v1=undef, v2=undef) +// → compute spline coefficient matrices for given control points +// - spline(s, t) +// → evaluate spline at parameter t +// - spline_tan(), spline_tan_unit() +// → evaluate tangent vector +// - spline_d2(), spline_binormal_unit(), spline_normal_unit() +// → evaluate higher derivatives and Frenet frame +// - spline_transform() +// → construct SE3 transform aligned to spline +// - bezier3_args(p, symmetric=false) +// → construct cubic Bezier coefficient matrices // // Internal helpers implement determinant, matrix inversion, etc. +// +// Author: Sergei Kuzmin, 2014 +// License: BSD // ============================================================================ use @@ -49,8 +48,8 @@ function matrix_power(m, n) = function det(m) = let (r = [for (i = [0:len(m) - 1]) i]) det_help(m, 0, r); -// Construction indices list is inefficient, but currently there is no way to imperatively -// assign to a list element +// Construction indices list is inefficient, but currently there is no way to +// imperatively assign to a list element function det_help(m, i, r) = len(r) == 0 ? 1 : m[len(m) - len(r)][r[i]] * det_help(m, 0, remove(r, i)) - (i + 1 < len(r) ? det_help(m, i + 1, r) : 0); @@ -99,8 +98,8 @@ function spline_args(p, closed = false, v1 = undef, v2 = undef) = ], sn = matrix_invert(q4 + q3 * matrix_power(qn1i2, pcnt - 2)) * (un - q3 * q1inv * spline_helper(0, pcnt, p)) ) - // result[i+1] recurrently defines result[i]. This is O(n) runtime with imperative language and - // may be O(n^2) if OpenSCAD doesn't cache spline_si(i+1). + // result[i+1] recurrently defines result[i]. This is O(n) runtime with + // imperative language and may be O(n^2) if OpenSCAD doesn't cache spline_si(i+1). [for (i = [0:pcnt - 2]) spline_si(i, pcnt - 2, p, sn)]; // Helper recursion for spline_args @@ -119,12 +118,14 @@ function spline_si(i, n, p, sn) = // --- Bezier Utilities ------------------------------------------------------- // Takes array of (3n+1) points or (2n + 2) points, if tangent segments are symmetric. -// For non-symmetric version input is: point0, normal0, neg_normal1, point1, normal1, ... neg_normal_n, point_n -// For symmetric version: point0, normal0, point1, normal1, ... , normal_n_sub_1, point_n +// For non-symmetric version input is: +// point0, normal0, neg_normal1, point1, normal1, ... neg_normal_n, point_n +// For symmetric version: +// point0, normal0, point1, normal1, ... , normal_n_sub_1, point_n // In the second case second tangent is constructed from the next tangent by symmetric map. -// I.e. if current points are p0,p1,p2 then anchor points are p0 and p2, first tangent defined by p1-p0, -// second tangent defined by p3-p2. -// Return array of coefficients accepted by spline(), spline_tan() and similar +// I.e. if current points are p0,p1,p2 then anchor points are p0 and p2, first tangent +// defined by p1-p0, second tangent defined by p3-p2. +// Return array of coefficients accepted by spline(), spline_tan() and similar function bezier3_args(p, symmetric = false) = let (step = symmetric ? 2 : 3) [ for (i = [0:step:len(p) - 3]) [[1, 0, 0, 0], [-3, 3, 0, 0], [3, -6, 3, 0], [-1, 3, -3, 1]] * ( @@ -136,7 +137,8 @@ function bezier3_args(p, symmetric = false) = // --- Evaluation ------------------------------------------------------------- // Evaluate spline/Bezier and derivatives // s - spline arguments calculated by spline_args -// t - defines point on curve. each segment length is 1. I.e. t= 0..1 is first segment, t=1..2 - second. +// t - defines point on curve. Each segment length is 1. +// i.e. t= 0..1 is first segment, t=1..2 - second. function spline(s, t) = let ( i = t >= len(s) ? len(s) - 1 : floor(t), diff --git a/trajectory_path.scad b/trajectory_path.scad index d28e345..580eafc 100644 --- a/trajectory_path.scad +++ b/trajectory_path.scad @@ -31,14 +31,14 @@ function right_multiply(as, b, i_ = 0) = : concat([as[i_] * b], right_multiply(as, b, i_ + 1)); // --- Single-trajectory quantization ---------------------------------------- -// Produces an array of 4x4 transforms sampled along one 6D twist `trajectory`. +// quantize_trajectory: produces an array of 4x4 transforms sampled along one 6D twist. // If `steps` is provided, it overrides `step` and distributes uniformly. // If `step` is provided, samples start at `start_position` and advance by `step`. // `start_position` and `step` are in the units of the translation norm. // // Notes: -// - Length is computed from the translational part only: norm(take3(trajectory)) -// - For steps==1, returns the transform at the path end (t=1). +// - Length is computed from the translational part only: norm(take3(trajectory)) +// - For steps==1, returns the transform at the path end (t=1). function quantize_trajectory( trajectory, step = undef, diff --git a/transformations.scad b/transformations.scad index bd1ff20..f68d4bd 100644 --- a/transformations.scad +++ b/transformations.scad @@ -19,18 +19,13 @@ use use // --- Rotation --------------------------------------------------------------- - -/*! - Creates a rotation matrix - - Options: - - xyz = Euler angles (applied as Rz * Ry * Rx) - - axis = axis-angle vector (axis * angle) - - Examples: - rotation(xyz=[90,0,0]) // rotate 90° about X - rotation(axis=[0,0,45]) // rotate 45° about Z -*/ +// Creates a rotation matrix +// Options: +// - xyz = Euler angles (applied as Rz * Ry * Rx) +// - axis = axis-angle vector (axis * angle) +// Examples: +// rotation(xyz=[90,0,0]) // rotate 90° about X +// rotation(axis=[0,0,45]) // rotate 45° about Z function rotation(xyz = undef, axis = undef) = xyz != undef && axis != undef ? undef : // cannot define both @@ -41,12 +36,7 @@ function rotation(xyz = undef, axis = undef) = : (len(xyz) >= 3 ? rotation(axis=[0, 0, xyz[2]]) : identity4()) * (len(xyz) >= 2 ? rotation(axis=[0, xyz[1], 0]) : identity4()) * (len(xyz) >= 1 ? rotation(axis=[xyz[0], 0, 0]) : identity4()); // --- Scaling ---------------------------------------------------------------- - -/*! - Creates a scaling matrix - - scaling([sx, sy, sz]) -*/ +// Creates a scaling matrix: scaling([sx, sy, sz]) function scaling(v) = [ [v[0], 0, 0, 0], @@ -56,12 +46,7 @@ function scaling(v) = ]; // --- Translation ------------------------------------------------------------ - -/*! - Creates a translation matrix - - translation([tx, ty, tz]) -*/ +// Creates a translation matrix: translation([tx, ty, tz]) function translation(v) = [ [1, 0, 0, v[0]], @@ -71,18 +56,11 @@ function translation(v) = ]; // --- Coordinate Conversion -------------------------------------------------- - -/*! - Converts from homogeneous to cartesian coordinates -*/ +// Converts from homogeneous to cartesian coordinates function project(x) = subarray(x, end=len(x) - 1) / x[len(x) - 1]; -/*! - Applies matrix `m` to a list of points -*/ +// Applies matrix `m` to a list of points function transform(m, list) = [for (p = list) project(m * vec4(p))]; -/*! - Ensures points are represented as 3D vectors -*/ +// Ensures points are represented as 3D vectors function to_3d(list) = [for (v = list) vec3(v)]; From 18f8c93951a5db7891b7a299142f88e16dab24e7 Mon Sep 17 00:00:00 2001 From: Cameron K Brooks Date: Sat, 27 Sep 2025 01:41:46 -0400 Subject: [PATCH 27/32] cut off extraneous stuff from readme --- README.md | 29 ----------------------------- 1 file changed, 29 deletions(-) diff --git a/README.md b/README.md index da9a525..5683329 100644 --- a/README.md +++ b/README.md @@ -2,19 +2,6 @@ Utility libraries for OpenSCAD -## Morphology - -contains basic 2D morphology operations - - inset(d=1) - creates a polygon at an offset d inside a 2D shape - outset(d=1) - creates a polygon at an offset d outside a 2D shape - fillet(r=1) - adds fillets of radius r to all concave corners of a 2D shape - rounding(r=1) - adds rounding to all convex corners of a 2D shape - shell(d,center=false) - makes a shell of width d along the edge of a 2D shape - - positive values of d places the shell on the outside - - negative values of d places the shell on the inside - - center=true and positive d places the shell centered on the edge - ## Examples With a basic sample polygon shape, @@ -52,19 +39,3 @@ and `$fn=32;`. - `shell(d=0.3,center=true) shape();` ![](http://oskarlinde.github.io/scad-utils/img/morph-6.png) - -## Mirror - -contains simple mirroring functions - -- mirror_x() -- mirror_y() -- mirror_z() - -example: - -``` -module arrow(l=1,w=.6,t=0.15) { - mirror_y() polygon([[0,0],[l,0],[l-w/2,w/2],[l-w/2-sqrt(2)*t,w/2],[l-t/2-sqrt(2)*t,t/2],[0,t/2]]); -} -``` From e0237c128c47cde6a19883dd2d74ae15060856e5 Mon Sep 17 00:00:00 2001 From: Cameron K Brooks Date: Sat, 27 Sep 2025 02:08:23 -0400 Subject: [PATCH 28/32] doc: update readme to document the library --- README.md | 351 ++++++++++++++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 338 insertions(+), 13 deletions(-) diff --git a/README.md b/README.md index 5683329..72ba066 100644 --- a/README.md +++ b/README.md @@ -1,41 +1,366 @@ # scad-utils -Utility libraries for OpenSCAD +Utility libraries for OpenSCAD. +This collection provides reusable math, geometry, and transformation tools for parametric modeling. +All modules are documented, consistently structured, and tested via dedicated `.scad` test files. + +--- + +## Modules + +### `hull.scad`: Convex Hull Utilities (2D and 3D) + +Computes convex hulls for 2D or 3D point sets with special handling for collinear cases. + +**Functions:** + +- `hull(points)` - Main entry point, automatically detects 2D/3D and handles edge cases +- Returns polygon vertex indices `[i1, i2, i3, ...]` for 2D +- Returns triangular face indices `[[i1,i2,i3], [i2,i3,i4], ...]` for 3D +- Returns two extreme endpoints `[i_min, i_max]` for collinear points + +**Features:** + +- Uses optimized algorithms with incremental construction +- Handles degenerate cases (collinear, coplanar points) +- Includes utility functions for spherical/cartesian conversions +- Compatible with both clockwise and counter-clockwise orientations + +### `linalg.scad`: Linear Algebra Utilities + +Essential vector and matrix operations for 3D transformations. + +**Vector Functions:** + +- `vec3(p)` - Ensure vector is 3D (pads with 0) +- `vec4(p)` - Ensure vector is 4D homogeneous (pads with 1) +- `unit(v)` - Normalize vector to unit length +- `take3(v)` - Extract first 3 elements +- `tail3(v)` - Extract elements 3,4,5 (for 6-vectors) + +**Matrix Functions:** + +- `identity3()`, `identity4()` - Identity matrices +- `transpose_3(m)`, `transpose_4(m)` - Matrix transpose +- `rotation_part(m)`, `translation_part(m)` - Extract parts from 4×4 transform +- `rot_trace(m)`, `rot_cos_angle(m)` - Rotation analysis +- `invert_rt(m)` - Invert rigid transform (rotation + translation) +- `construct_Rt(R, t)` - Build 4×4 transform from 3×3 rotation and translation +- `hadamard(a, b)` - Element-wise multiplication (works recursively) + +### `lists.scad`: List Manipulation Utilities + +Functional programming helpers for array operations. + +**Functions:** + +- `flatten(list)` - Flatten nested arrays one level: `[[0,1],[2,3]] → [0,1,2,3]` +- `range(r)` - Convert range to list: `[0:2:6] → [0,2,4,6]` +- `reverse(list)` - Reverse element order: `[1,2,3] → [3,2,1]` +- `subarray(list, begin, end)` - Extract slice (end=-1 means full length) +- `set(list, i, x)` - Return copy with element at index i replaced by x +- `remove(list, i)` - Return copy with element at index i removed + +### `mirror.scad`: Axis Mirroring Modules + +Simple mirroring operations that duplicate and reflect geometry. + +**Modules:** + +- `mirror_x([col])` - Mirror across YZ-plane (flip X coordinate) +- `mirror_y([col])` - Mirror across XZ-plane (flip Y coordinate) +- `mirror_z([col])` - Mirror across XY-plane (flip Z coordinate) + +**Parameters:** + +- `col` - Optional color name to tint the mirrored copy (e.g., "red", "teal") + +Each module creates a union of the original children plus the mirrored copy. + +### `morphology.scad`: 2D Morphology Operations + +Advanced 2D shape modification operations for polygon processing. + +**Core Operations:** + +- `outset(d=1)` - Grow shape outward by distance d (Minkowski sum with circle) +- `inset(d=1)` - Shrink shape inward by distance d (inverse of outset) +- `fillet(r=1)` - Add rounded fillets to concave (inward) corners +- `rounding(r=1)` - Round convex (outward) corners +- `shell(d, center=false)` - Create shell/band around shape edge + +**Shell Parameters:** + +- `d > 0, center=false`: shell extends outward +- `d < 0, center=false`: shell extends inward +- `center=true`: shell straddles the original edge + +**Compatibility:** + +- Works around OpenSCAD version differences in Minkowski operations +- Uses `render()` and `projection()` fallbacks for older versions + +### `se3.scad`: SE(3) Lie Group Utilities + +Exponential and logarithm maps for rigid body transformations (translation + rotation). + +**Functions:** + +- `se3_exp(mu)` - Convert 6D twist vector `[tx,ty,tz, rx,ry,rz]` to 4×4 transform matrix +- `se3_ln(m)` - Convert 4×4 transform matrix back to 6D twist vector + +**Features:** + +- Handles small-angle approximations with Taylor series (1st, 2nd, 3rd order) +- Rotation angles specified in degrees (converted internally) +- Combines SO(3) rotations with translation using proper Lie algebra +- Numerical stability for near-identity transforms + +**Dependencies:** `linalg.scad`, `so3.scad` + +### `shapes.scad`: 2D Shape Generators + +Parametric generators for common 2D shapes returned as point arrays. + +**Functions:** + +- `square(size)` - Centered square with edge length `size` +- `circle(r)` - Circle with radius `r` (uses global `$fn` for resolution) +- `regular(r, n)` - Regular n-sided polygon with circumradius `r` +- `rectangle_profile(size=[w,h])` - Rectangle with anchor at `[w/2, 0]` + +**Output:** All functions return arrays of 2D points suitable for `polygon()`. + +### `so3.scad`: SO(3) Lie Group Utilities + +Exponential and logarithm maps for 3D rotations using Rodrigues formula. + +**Functions:** + +- `so3_exp(w)` - Convert axis-angle vector (degrees) to 3×3 rotation matrix +- `so3_ln(m)` - Convert 3×3 rotation matrix back to axis-angle vector (degrees) +- `so3_exp_rad(w)`, `so3_ln_rad(m)` - Radian versions for internal use + +**Features:** + +- Taylor expansions for small angles (1st, 2nd, 3rd order approximations) +- Handles near-π rotations with symmetric matrix decomposition +- Rodrigues formula implementation: `R = I + sin(θ)K + (1-cos(θ))K²` +- Numerical stability across full rotation range + +**Dependencies:** `linalg.scad` + +### `spline.scad`: Cubic Spline and Bezier Utilities + +Comprehensive curve generation and evaluation with Frenet frame support. + +**Spline Functions:** + +- `spline_args(points, closed=false, v1=undef, v2=undef)` - Generate spline coefficients +- `spline(s, t)` - Evaluate position at parameter t +- `spline_tan(s, t)` - Evaluate tangent vector +- `spline_tan_unit(s, t)` - Unit tangent vector +- `spline_d2(s, t)` - Second derivative (curvature) +- `spline_normal_unit(s, t)`, `spline_binormal_unit(s, t)` - Frenet frame vectors +- `spline_transform(s, t)` - SE(3) transform aligned to curve + +**Bezier Functions:** + +- `bezier3_args(points, symmetric=false)` - Generate cubic Bezier coefficients + +**Features:** + +- Supports open and closed splines with customizable end conditions +- Frenet frame calculations for swept surfaces and path following +- Matrix-based coefficient computation with automatic boundary conditions +- Parameter t scales with curve segments (t=0..1 first segment, t=1..2 second, etc.) + +**Dependencies:** `linalg.scad`, `lists.scad` + +### `trajectory.scad`: SE(3) Twist Vector Construction + +Intuitive interface for building 6D motion vectors from directional parameters. + +**Functions:** + +- `trajectory(left/right, up/down, forward/backward, translation, pitch, yaw, roll, rotation)` +- `rotationv(pitch, yaw, roll, rotation)` - Build rotation component +- `translationv(...)` - Build translation component +- `rotationm(...)` - Convert angles to 3×3 rotation matrix + +**Direction Parameters:** + +- Translation: `left/right`, `up/down`, `forward/backward` OR explicit `translation=[x,y,z]` +- Rotation: `pitch`, `yaw`, `roll` (degrees) OR explicit `rotation=[rx,ry,rz]` + +**Output:** Returns `[tx,ty,tz, yaw,pitch,roll]` 6D twist vector suitable for `se3_exp()` + +**Dependencies:** `so3.scad` + +### `trajectory_path.scad`: Multi-Segment Path Quantization + +Convert trajectory sequences into discrete transformation samples. + +**Functions:** + +- `quantize_trajectory(trajectory, step/steps, start_position)` - Sample single 6D twist +- `quantize_trajectories(trajectories, step/steps, start_position, loop)` - Sample path sequence +- `close_trajectory_loop(trajectories)` - Add segment to close loop +- `trajectories_length(trajectories)` - Compute total path length +- `trajectories_end_position(trajectories)` - Final accumulated transform + +**Parameters:** + +- `step`: Physical step length (units of translation norm) +- `steps`: Fixed number of uniform samples across entire path +- `start_position`: Arc-length offset for starting point +- `loop=true`: Automatically close path back to start + +**Output:** Arrays of 4×4 transformation matrices for each sample point + +**Dependencies:** `linalg.scad`, `se3.scad` + +### `transformations.scad`: Geometric Transformation Matrices + +High-level constructors for common 3D transformations. + +**Matrix Constructors:** + +- `rotation(xyz=[rx,ry,rz])` - Euler angles (Rz·Ry·Rx order) +- `rotation(axis=[x,y,z])` - Axis-angle rotation +- `scaling([sx,sy,sz])` - Non-uniform scaling +- `translation([tx,ty,tz])` - Translation matrix + +**Utility Functions:** + +- `project(x)` - Convert homogeneous to Cartesian coordinates +- `transform(m, points)` - Apply matrix to point list +- `to_3d(points)` - Ensure points are 3D vectors + +**Usage:** Matrices can be multiplied for composition: `T * R * S` applies scaling, then rotation, then translation. + +**Dependencies:** `se3.scad`, `linalg.scad`, `lists.scad` ## Examples -With a basic sample polygon shape, +### Morphology Operations - module shape() { - polygon([[0,0],[1,0],[1.5,1],[2.5,1],[2,-1],[0,-1]]); - } +With a basic sample polygon shape, -and `$fn=32;`. +```scad +module shape() { + polygon([[0,0],[1,0],[1.5,1],[2.5,1],[2,-1],[0,-1]]); +} +$fn = 32; +``` - `inset(d=0.3) shape();` -![](http://oskarlinde.github.io/scad-utils/img/morph-0.png) +![Inset Morphology Example](http://oskarlinde.github.io/scad-utils/img/morph-0.png) - `outset(d=0.3) shape();` -![](http://oskarlinde.github.io/scad-utils/img/morph-1.png) +![Outset Morphology Example](http://oskarlinde.github.io/scad-utils/img/morph-1.png) - `rounding(r=0.3) shape();` -![](http://oskarlinde.github.io/scad-utils/img/morph-2.png) +![Rounding Morphology Example](http://oskarlinde.github.io/scad-utils/img/morph-2.png) - `fillet(r=0.3) shape();` -![](http://oskarlinde.github.io/scad-utils/img/morph-3.png) +![Fillet Morphology Example](http://oskarlinde.github.io/scad-utils/img/morph-3.png) - `shell(d=0.3) shape();` -![](http://oskarlinde.github.io/scad-utils/img/morph-4.png) +![Shell Morphology Example Positive](http://oskarlinde.github.io/scad-utils/img/morph-4.png) - `shell(d=-0.3) shape();` -![](http://oskarlinde.github.io/scad-utils/img/morph-5.png) +![Shell Morphology Example Negative](http://oskarlinde.github.io/scad-utils/img/morph-5.png) - `shell(d=0.3,center=true) shape();` -![](http://oskarlinde.github.io/scad-utils/img/morph-6.png) +![Shell Morphology Example Centered](http://oskarlinde.github.io/scad-utils/img/morph-6.png) + +### Mirror Operations + +```scad +module arrow(l=1, w=0.6, t=0.15) { + mirror_y("orange") + polygon([[0,0], [l,0], [l-w/2,w/2], + [l-w/2-sqrt(2)*t,w/2], + [l-t/2-sqrt(2)*t,t/2], [0,t/2]]); +} + +arrow(l=20, w=10, t=2); +``` + +### Convex Hull + +```scad +use + +points_3d = [[0,0,0], [10,0,0], [5,10,0], [0,0,10]]; +faces = hull(points_3d); +polyhedron(points=points_3d, faces=faces); +``` + +### Spline Curves + +```scad +use + +points = [[0,0,0], [10,5,0], [20,0,5], [30,10,0]]; +spline_data = spline_args(points, closed=true); + +for (t = [0:0.1:len(spline_data)]) + translate(spline(spline_data, t)) + sphere(r=0.5); +``` + +### SE(3) Transformations + +```scad +use +use + +// Define a 6D twist: translate [10,0,5] + rotate 45° about Z +twist = trajectory(forward=10, up=5, yaw=45); +transform_matrix = se3_exp(twist); + +multmatrix(transform_matrix) + cube([2,2,2]); +``` + +### Multi-Segment Paths + +```scad +use + +// Define path segments +path = [ + [20, 0, 0, 0, 0, 30], // forward + yaw + [0, 15, 0, 0, 45, 0], // right + pitch + [0, 0, 10, 0, 0, -60] // up + yaw back +]; + +// Sample every 5 units +poses = quantize_trajectories(path, step=5); + +for (T = poses) + multmatrix(T) + cube([1,1,1]); +``` + +## `tests/` Directory + +Each module includes comprehensive test files in the `tests/` directory: + +**Test Coverage:** + +- **Regression Tests:** Echo-based validation of mathematical properties (e.g., `so3_ln(so3_exp(w)) ≈ w`) +- **Visual Tests:** Geometric demonstrations with color-coded results +- **Edge Cases:** Boundary conditions, degenerate inputs, numerical stability +- **Integration Tests:** Multi-module workflows (splines with SE(3), trajectory chains) + +**Running Tests:** Load test files directly in OpenSCAD to see both console output and visual results. From 0bbaf23abdbd156719fd510808e4a7373ac7be4e Mon Sep 17 00:00:00 2001 From: Cameron K Brooks Date: Sat, 27 Sep 2025 02:15:52 -0400 Subject: [PATCH 29/32] move matrix utils from spline to linalg --- linalg.scad | 33 +++++++++++++++++++++++++++++++++ spline.scad | 32 -------------------------------- tests/test_linalg.scad | 38 ++++++++++++++++++++------------------ 3 files changed, 53 insertions(+), 50 deletions(-) diff --git a/linalg.scad b/linalg.scad index 4c22992..1db0640 100644 --- a/linalg.scad +++ b/linalg.scad @@ -9,6 +9,8 @@ // - Hadamard (elementwise) product // ============================================================================ +use + epsilon = 1e-9; // --- Vector Constructors ---------------------------------------------------- @@ -50,6 +52,37 @@ function rotation_part(m) = function translation_part(m) = [m[0][3], m[1][3], m[2][3]]; +// --- Matrix Utilities ------------------------------------------------------- + +// Compute matrix power +function matrix_power(m, n) = + n == 0 ? (len(m) == 3 ? identity3() : identity4()) + : n == 1 ? m + : (n % 2 == 1) ? matrix_power(m * m, floor(n / 2)) * m + : matrix_power(m * m, n / 2); + +// Determinant (recursive Laplace expansion) +function det(m) = + let (r = [for (i = [0:len(m) - 1]) i]) det_help(m, 0, r); + +// Construction indices list is inefficient, but currently there is no way to +// imperatively assign to a list element +function det_help(m, i, r) = + len(r) == 0 ? 1 + : m[len(m) - len(r)][r[i]] * det_help(m, 0, remove(r, i)) - (i + 1 < len(r) ? det_help(m, i + 1, r) : 0); + +// Matrix inversion (adjugate method) +function matrix_invert(m) = + let (r = [for (i = [0:len(m) - 1]) i]) [ + for (i = r) [ + for (j = r) ( (i + j) % 2 == 0 ? 1 : -1) * matrix_minor(m, 0, remove(r, j), remove(r, i)), + ], + ] / det(m); + +function matrix_minor(m, k, ri, rj) = + let (len_r = len(ri)) len_r == 0 ? 1 + : m[ri[0]][rj[k]] * matrix_minor(m, 0, remove(ri, 0), remove(rj, k)) - (k + 1 < len_r ? matrix_minor(m, k + 1, ri, rj) : 0); + // --- Rotation Metrics ------------------------------------------------------- function rot_trace(m) = m[0][0] + m[1][1] + m[2][2]; diff --git a/spline.scad b/spline.scad index 8912f21..376193f 100644 --- a/spline.scad +++ b/spline.scad @@ -24,7 +24,6 @@ // ============================================================================ use -use // --- Predefined Matrices ---------------------------------------------------- q1 = [[1, 0, 0, 0], [1, 1, 1, 1], [0, 1, 2, 3], [0, 0, 1, 3]]; @@ -35,37 +34,6 @@ qn1i2 = -q1inv * q2; z3 = [0, 0, 0]; z4 = [0, 0, 0, 0]; -// --- Matrix Utilities ------------------------------------------------------- - -// Compute matrix power -function matrix_power(m, n) = - n == 0 ? (len(m) == 3 ? identity3() : identity4()) - : n == 1 ? m - : (n % 2 == 1) ? matrix_power(m * m, floor(n / 2)) * m - : matrix_power(m * m, n / 2); - -// Determinant (recursive Laplace expansion) -function det(m) = - let (r = [for (i = [0:len(m) - 1]) i]) det_help(m, 0, r); - -// Construction indices list is inefficient, but currently there is no way to -// imperatively assign to a list element -function det_help(m, i, r) = - len(r) == 0 ? 1 - : m[len(m) - len(r)][r[i]] * det_help(m, 0, remove(r, i)) - (i + 1 < len(r) ? det_help(m, i + 1, r) : 0); - -// Matrix inversion (adjugate method) -function matrix_invert(m) = - let (r = [for (i = [0:len(m) - 1]) i]) [ - for (i = r) [ - for (j = r) ( (i + j) % 2 == 0 ? 1 : -1) * matrix_minor(m, 0, remove(r, j), remove(r, i)), - ], - ] / det(m); - -function matrix_minor(m, k, ri, rj) = - let (len_r = len(ri)) len_r == 0 ? 1 - : m[ri[0]][rj[k]] * matrix_minor(m, 0, remove(ri, 0), remove(rj, k)) - (k + 1 < len_r ? matrix_minor(m, k + 1, ri, rj) : 0); - // --- Spline Construction ---------------------------------------------------- // Local utility diff --git a/tests/test_linalg.scad b/tests/test_linalg.scad index 62b6b10..9033939 100644 --- a/tests/test_linalg.scad +++ b/tests/test_linalg.scad @@ -1,24 +1,24 @@ include <../linalg.scad>; // --- Vector constructors ---------------------------------------------------- -echo("vec3([1,2]) =", vec3([1,2])); // expect [1,2,0] -echo("vec4([1,2,3]) =", vec4([1,2,3])); // expect [1,2,3,1] -echo("unit([3,0,0]) =", unit([3,0,0])); // expect [1,0,0] +echo("vec3([1,2]) =", vec3([1, 2])); // expect [1,2,0] +echo("vec4([1,2,3]) =", vec4([1, 2, 3])); // expect [1,2,3,1] +echo("unit([3,0,0]) =", unit([3, 0, 0])); // expect [1,0,0] // --- Identity matrices ------------------------------------------------------ -echo("identity3() =", identity3()); // expect 3x3 identity -echo("identity4() =", identity4()); // expect 4x4 identity +echo("identity3() =", identity3()); // expect 3x3 identity +echo("identity4() =", identity4()); // expect 4x4 identity // --- Vector access ---------------------------------------------------------- -echo("take3([9,8,7,6]) =", take3([9,8,7,6])); // expect [9,8,7] -echo("tail3([0,1,2,3,4,5]) =", tail3([0,1,2,3,4,5])); // expect [3,4,5] +echo("take3([9,8,7,6]) =", take3([9, 8, 7, 6])); // expect [9,8,7] +echo("tail3([0,1,2,3,4,5]) =", tail3([0, 1, 2, 3, 4, 5])); // expect [3,4,5] // --- Matrix parts ----------------------------------------------------------- M = [ - [1,0,0,5], - [0,1,0,6], - [0,0,1,7], - [0,0,0,1] + [1, 0, 0, 5], + [0, 1, 0, 6], + [0, 0, 1, 7], + [0, 0, 0, 1], ]; echo("rotation_part(M) =", rotation_part(M)); // expect identity3() echo("translation_part(M) =", translation_part(M)); // expect [5,6,7] @@ -29,7 +29,7 @@ echo("rot_trace(R) =", rot_trace(R)); // expect 3 echo("rot_cos_angle(R) =", rot_cos_angle(R)); // expect 1 // --- Transpose -------------------------------------------------------------- -A3 = [[1,2,3],[4,5,6],[7,8,9]]; +A3 = [[1, 2, 3], [4, 5, 6], [7, 8, 9]]; echo("transpose_3(A3) =", transpose_3(A3)); // expect [[1,4,7],[2,5,8],[3,6,9]] A4 = identity4(); @@ -39,15 +39,17 @@ echo("transpose_4(identity4) =", transpose_4(A4)); // expect identity4() echo("invert_rt(M) =", invert_rt(M)); // expect transform with R=I, t=[-5,-6,-7] -R2 = [[0,-1,0],[1,0,0],[0,0,1]]; // 90° rot about Z -t2 = [10,0,0]; -Rt = construct_Rt(R2,t2); +R2 = [[0, -1, 0], [1, 0, 0], [0, 0, 1]]; // 90° rot about Z +t2 = [10, 0, 0]; +Rt = construct_Rt(R2, t2); echo("construct_Rt(R2,t2) =", Rt); // --- Hadamard product ------------------------------------------------------ -echo("hadamard([1,2,3],[4,5,6]) =", hadamard([1,2,3],[4,5,6])); +echo("hadamard([1,2,3],[4,5,6]) =", hadamard([1, 2, 3], [4, 5, 6])); // expect [4,10,18] -echo("hadamard([[1,2],[3,4]], [[5,6],[7,8]]) =", - hadamard([[1,2],[3,4]], [[5,6],[7,8]])); +echo( + "hadamard([[1,2],[3,4]], [[5,6],[7,8]]) =", + hadamard([[1, 2], [3, 4]], [[5, 6], [7, 8]]) +); // expect [[5,12],[21,32]] From 22be1619e49284285a697033b63dcb1f9cd09286 Mon Sep 17 00:00:00 2001 From: Cameron K Brooks Date: Sat, 27 Sep 2025 02:16:43 -0400 Subject: [PATCH 30/32] add tests for moved linalg matrix utils --- tests/test_linalg.scad | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/tests/test_linalg.scad b/tests/test_linalg.scad index 9033939..d42e463 100644 --- a/tests/test_linalg.scad +++ b/tests/test_linalg.scad @@ -53,3 +53,13 @@ echo( hadamard([[1, 2], [3, 4]], [[5, 6], [7, 8]]) ); // expect [[5,12],[21,32]] + +// --- Matrix utilities ------------------------------------------------------- +B = [[2, 0], [0, 3]]; +echo("matrix_power(B,0) =", matrix_power(B, 0)); // expect identity2 (approximated as [[1,0],[0,1]]) +echo("matrix_power(B,2) =", matrix_power(B, 2)); // expect [[4,0],[0,9]] + +C = [[1, 2], [3, 4]]; +echo("det(C) =", det(C)); // expect -2 +echo("matrix_invert(C) =", matrix_invert(C)); +// expect [[-2,1],[1.5,-0.5]] From 57efdaf81c18d7f3b4dcf0ebe8680546620b9a74 Mon Sep 17 00:00:00 2001 From: Cameron K Brooks Date: Sun, 5 Oct 2025 22:02:57 -0400 Subject: [PATCH 31/32] add links to projects using scad-utils at end of readme and swap undef checks to use is_undef --- README.md | 14 ++++++++++++-- spline.scad | 8 ++++---- transformations.scad | 33 ++++++++++++++++++++++++++++----- 3 files changed, 44 insertions(+), 11 deletions(-) diff --git a/README.md b/README.md index 72ba066..9744b33 100644 --- a/README.md +++ b/README.md @@ -2,7 +2,8 @@ Utility libraries for OpenSCAD. This collection provides reusable math, geometry, and transformation tools for parametric modeling. -All modules are documented, consistently structured, and tested via dedicated `.scad` test files. +All modules are documented, consistently structured, and tested via dedicated `.scad` test files. +A list of projects making use of these utilities can be viewed in the [projects](#projects-using-scad-utils) section below. --- @@ -285,6 +286,8 @@ $fn = 32; ### Mirror Operations ```scad +use + module arrow(l=1, w=0.6, t=0.15) { mirror_y("orange") polygon([[0,0], [l,0], [l-w/2,w/2], @@ -310,7 +313,7 @@ polyhedron(points=points_3d, faces=faces); ```scad use -points = [[0,0,0], [10,5,0], [20,0,5], [30,10,0]]; +points = [[0,0,0], [10,10,0], [20,0,5], [30,10,0]]; spline_data = spline_args(points, closed=true); for (t = [0:0.1:len(spline_data)]) @@ -364,3 +367,10 @@ Each module includes comprehensive test files in the `tests/` directory: - **Integration Tests:** Multi-module workflows (splines with SE(3), trajectory chains) **Running Tests:** Load test files directly in OpenSCAD to see both console output and visual results. + +## Projects using `scad-utils` + +- [openscad/list-comprehension-demos](https://github.com/openscad/list-comprehension-demos) uses `lists.scad`, `linalg.scad`, `transformations.scad`, `shape.scad`, `trajectory.scad`, and `trajectory_path.scad` +- [adrianschlatter/threadlib](https://github.com/adrianschlatter/threadlib) uses `lists.scad` and `transformations.scad` through [MisterHW/IoP-satellite](https://github.com/MisterHW/IoP-satellite) +- [likeablob/parametric-stackable-box](https://github.com/likeablob/parametric-stackable-box) +- [likeablob/misc-printable-accessories](https://github.com/likeablob/misc-printable-accessories) diff --git a/spline.scad b/spline.scad index 376193f..71247d7 100644 --- a/spline.scad +++ b/spline.scad @@ -47,7 +47,7 @@ function spline_args(p, closed = false, v1 = undef, v2 = undef) = : [ z4, z4, - v1 == undef ? [0, 0, 1, 0] : [0, 1, 0, 0], + is_undef(v1) ? [0, 0, 1, 0] : [0, 1, 0, 0], z4, ], q4 = closed ? q1 @@ -55,14 +55,14 @@ function spline_args(p, closed = false, v1 = undef, v2 = undef) = [1, 0, 0, 0], [1, 1, 1, 1], z4, - v2 == undef ? [0, 0, 1, 3] : [0, 1, 2, 3], + is_undef(v2) ? [0, 0, 1, 3] : [0, 1, 2, 3], ], pcnt = closed ? len(p) + 1 : len(p), un = [ p[pcnt - 2], p[closed ? 0 : pcnt - 1], - v1 == undef ? z4 : v1, - v2 == undef ? z4 : v2, + is_undef(v1) ? z4 : v1, + is_undef(v2) ? z4 : v2, ], sn = matrix_invert(q4 + q3 * matrix_power(qn1i2, pcnt - 2)) * (un - q3 * q1inv * spline_helper(0, pcnt, p)) ) diff --git a/transformations.scad b/transformations.scad index f68d4bd..e6b55f1 100644 --- a/transformations.scad +++ b/transformations.scad @@ -27,13 +27,36 @@ use // rotation(xyz=[90,0,0]) // rotate 90° about X // rotation(axis=[0,0,45]) // rotate 45° about Z function rotation(xyz = undef, axis = undef) = - xyz != undef && axis != undef ? undef - : // cannot define both - xyz == undef ? + // disallow both forms together + (!is_undef(xyz) && !is_undef(axis)) ? + undef + : + // pure axis-angle exponential form + (is_undef(xyz) && !is_undef(axis)) ? se3_exp([0, 0, 0, axis[0], axis[1], axis[2]]) - : len(xyz) == undef ? + : + // shorthand for single-angle rotation about Z + (is_undef(axis) && !is_undef(xyz) && !is_list(xyz)) ? rotation(axis=[0, 0, xyz]) - : (len(xyz) >= 3 ? rotation(axis=[0, 0, xyz[2]]) : identity4()) * (len(xyz) >= 2 ? rotation(axis=[0, xyz[1], 0]) : identity4()) * (len(xyz) >= 1 ? rotation(axis=[xyz[0], 0, 0]) : identity4()); + : + // full Euler xyz case + (is_undef(axis) && is_list(xyz)) ? + ( + len(xyz) >= 3 ? + rotation(axis=[0, 0, xyz[2]]) + : identity4() + ) * ( + len(xyz) >= 2 ? + rotation(axis=[0, xyz[1], 0]) + : identity4() + ) * ( + len(xyz) >= 1 ? + rotation(axis=[xyz[0], 0, 0]) + : identity4() + ) + : + // fallback + identity4(); // --- Scaling ---------------------------------------------------------------- // Creates a scaling matrix: scaling([sx, sy, sz]) From a64094daf6134ae6ec2138ae5afd6a52e28924fe Mon Sep 17 00:00:00 2001 From: Cameron K Brooks Date: Sun, 5 Oct 2025 22:05:05 -0400 Subject: [PATCH 32/32] details of what likeablobs libs use in scad-utils --- README.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/README.md b/README.md index 9744b33..6445101 100644 --- a/README.md +++ b/README.md @@ -372,5 +372,5 @@ Each module includes comprehensive test files in the `tests/` directory: - [openscad/list-comprehension-demos](https://github.com/openscad/list-comprehension-demos) uses `lists.scad`, `linalg.scad`, `transformations.scad`, `shape.scad`, `trajectory.scad`, and `trajectory_path.scad` - [adrianschlatter/threadlib](https://github.com/adrianschlatter/threadlib) uses `lists.scad` and `transformations.scad` through [MisterHW/IoP-satellite](https://github.com/MisterHW/IoP-satellite) -- [likeablob/parametric-stackable-box](https://github.com/likeablob/parametric-stackable-box) -- [likeablob/misc-printable-accessories](https://github.com/likeablob/misc-printable-accessories) +- [likeablob/parametric-stackable-box](https://github.com/likeablob/parametric-stackable-box) uses `morphology.scad` and `mirror.scad` +- [likeablob/misc-printable-accessories](https://github.com/likeablob/misc-printable-accessories) uses `morphology.scad` and `mirror.scad`