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

Merge pull request #580 from RonaldoCMP/master

Browse source 
introduce convex collision and distance
This commit is contained in:
Revar Desmera 2021-06-22 14:12:43 -07:00 committed by GitHub
commit 7a3720a812
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
10 changed files with 476 additions and 162 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

112
arrays.scad
View file

@ -59,14 +59,17 @@ function _same_type(a,b, depth) =
// Returns a portion of a list, wrapping around past the beginning, if end<start.
// The first item is index 0. Negative indexes are counted back from the end.
// The last item is -1. If only the `start` index is given, returns just the value
// at that position.
// at that position when `start` is a number or the selected list of entries when `start` is
// a list of indices or a range.
// Usage:
// item = select(list,start);
// item = select(list,[s:d:e]);
// item = select(list,[i0,i1...,ik]);
// list = select(list,start,end);
// Arguments:
// list = The list to get the portion of.
// start = The index of the first item.
// end = The index of the last item.
// start = Either the index of the first item or an index range or a list of indices.
// end = The index of the last item when `start` is a number. When `start` is a list or a range, `end` should not be given.
// See Also: slice(), subindex(), last()
// Example:
// l = [3,4,5,6,7,8,9];
@ -78,7 +81,7 @@ function _same_type(a,b, depth) =
// f = select(l, 4); // Returns 7
// g = select(l, -2); // Returns 8
// h = select(l, [1:3]); // Returns [4,5,6]
// i = select(l, [1,3]); // Returns [4,6]
// i = select(l, [3,1]); // Returns [6,4]
function select(list, start, end) =
assert( is_list(list) || is_string(list), "Invalid list.")
let(l=len(list))
@ -89,7 +92,7 @@ function select(list, start, end) =
? list[ (start%l+l)%l ]
: assert( is_list(start) || is_range(start), "Invalid start parameter")
[for (i=start) list[ (i%l+l)%l ] ]
: assert(is_finite(start), "Invalid start parameter.")
: assert(is_finite(start), "When `end` is given, `start` parameter should be a number.")
assert(is_finite(end), "Invalid end parameter.")
let( s = (start%l+l)%l, e = (end%l+l)%l )
(s <= e)
@ -946,6 +949,31 @@ function shuffle(list,seed) =
)
concat(shuffle(left), shuffle(right));
// idx should be an index of the arrays l[i]
function _group_sort_by_index(l,idx) =
len(l) == 0 ? [] :
len(l) == 1 ? [l] :
let( pivot = l[floor(len(l)/2)][idx],
equal = [ for(li=l) if( li[idx]==pivot) li ] ,
lesser = [ for(li=l) if( li[idx]< pivot) li ] ,
greater = [ for(li=l) if( li[idx]> pivot) li ]
)
concat( _group_sort_by_index(lesser,idx),
[equal],
_group_sort_by_index(greater,idx) ) ;
function _group_sort(l) =
len(l) == 0 ? [] :
len(l) == 1 ? [l] :
let( pivot = l[floor(len(l)/2)],
equal = [ for(li=l) if( li==pivot) li ] ,
lesser = [ for(li=l) if( li< pivot) li ] ,
greater = [ for(li=l) if( li> pivot) li ]
)
concat( _group_sort(lesser),
[equal],
_group_sort(greater) ) ;
// Sort a vector of scalar values with the native comparison operator
// all elements should have the same type.
@ -1007,7 +1035,7 @@ function _sort_general(arr, idx=undef, indexed=false) =
// lexical sort using compare_vals()
function _lexical_sort(arr) =
arr==[] ? [] : len(arr)==1? arr :
len(arr)<=1? arr :
let( pivot = arr[floor(len(arr)/2)] )
let(
lesser = [ for (entry=arr) if (compare_vals(entry, pivot) <0 ) entry ],
@ -1034,7 +1062,7 @@ function _indexed_sort(arrind) =
// Usage:
// slist = sort(list, <idx>);
// Topics: List Handling
// See Also: shuffle(), sortidx(), unique(), unique_count()
// See Also: shuffle(), sortidx(), unique(), unique_count(), group_sort()
// Description:
// Sorts the given list in lexicographic order. If the input is a homogeneous simple list or a homogeneous
// list of vectors (see function is_homogeneous), the sorting method uses the native comparison operator and is faster.
@ -1075,7 +1103,7 @@ function sort(list, idx=undef) =
// Usage:
// idxlist = sortidx(list, <idx>);
// Topics: List Handling
// See Also: shuffle(), sort(), unique(), unique_count()
// See Also: shuffle(), sort(), group_sort(), unique(), unique_count()
// Description:
// Given a list, sort it as function `sort()`, and returns
// a list of indexes into the original list in that sorted order.
@ -1123,26 +1151,69 @@ function sortidx(list, idx=undef) =
: _sort_general(list,idx,indexed=true);
// Function: group_sort()
// Usage:
// ulist = group_sort(list);
// Topics: List Handling
// See Also: shuffle(), sort(), sortidx(), unique(), unique_count()
// Description:
// Given a list of values, returns the sorted list with all repeated items grouped in a list.
// When the list entries are themselves lists, the sorting may be done based on the `idx` entry
// of those entries, that should be numbers.
// The result is always a list of lists.
// Arguments:
// list = The list to sort.
// idx = If given, do the comparison based just on the specified index. Default: zero.
// Example:
// sorted = group_sort([5,2,8,3,1,3,8,7,5]); // Returns: [[1],[2],[3,3],[5,5],[7],[8,8]]
// sorted2 = group_sort([[5,"a"],[2,"b"], [5,"c"], [3,"d"], [2,"e"] ], idx=0); // Returns: [[[2,"b"],[2,"e"]], [[5,"a"],[5,"c"]], [[3,"d"]] ]
function group_sort(list, idx) =
assert(is_list(list), "Input should be a list." )
assert(is_undef(idx) || (is_finite(idx) && idx>=0) , "Invalid index." )
len(list)<=1 ? [list] :
is_vector(list)? _group_sort(list) :
let( idx = is_undef(idx) ? 0 : idx )
assert( [for(entry=list) if(!is_list(entry) || len(entry)<idx || !is_num(entry[idx]) ) 1]==[],
"Some entry of the list is a list shorter than `idx` or the indexed entry of it is not a number." )
_group_sort_by_index(list,idx);
// Function: unique()
// Usage:
// ulist = unique(list);
// Topics: List Handling
// See Also: shuffle(), sort(), sortidx(), unique_count()
// Description:
// Returns a sorted list with all repeated items removed.
// Given a string or a list returns the sorted string or the sorted list with all repeated items removed.
// The sorting order of non homogeneous lists is the function `sort` order.
// Arguments:
// list = The list to uniquify.
// Example:
// sorted = unique([5,2,8,3,1,3,8,7,5]); // Returns: [1,2,3,5,7,8]
// sorted = unique("axdbxxc"); // Returns: "abcdx"
// sorted = unique([true,2,"xba",[1,0],true,[0,0],3,"a",[0,0],2]); // Returns: [true,2,3,"a","xba",[0,0],[1,0]]
function unique(list) =
assert(is_list(list)||is_string(list), "Invalid input." )
is_string(list)? str_join(unique([for (x = list) x])) :
len(list)<=1? list :
let( sorted = sort(list))
[ for (i=[0:1:len(sorted)-1])
if (i==0 || (sorted[i] != sorted[i-1]))
sorted[i]
];
is_homogeneous(list,1) && ! is_list(list[0])
? _unique_sort(list)
: let( sorted = sort(list))
[ for (i=[0:1:len(sorted)-1])
if (i==0 || (sorted[i] != sorted[i-1]))
sorted[i]
];
function _unique_sort(l) =
len(l) <= 1 ? l :
let( pivot = l[floor(len(l)/2)],
equal = [ for(li=l) if( li==pivot) li ] ,
lesser = [ for(li=l) if( li<pivot ) li ] ,
greater = [ for(li=l) if( li>pivot) li ]
)
concat( _unique_sort(lesser),
equal[0],
_unique_sort(greater) ) ;
// Function: unique_count()
@ -1160,11 +1231,14 @@ function unique(list) =
function unique_count(list) =
assert(is_list(list) || is_string(list), "Invalid input." )
list == [] ? [[],[]] :
let( list=sort(list) )
let( ind = [0, for(i=[1:1:len(list)-1]) if (list[i]!=list[i-1]) i] )
[ select(list,ind), deltas( concat(ind,[len(list)]) ) ];
is_homogeneous(list,1) && ! is_list(list[0])
? let( sorted = _group_sort(list) )
[ [for(s=sorted) s[0] ], [for(s=sorted) len(s) ] ]
: let( list=sort(list) )
let( ind = [0, for(i=[1:1:len(list)-1]) if (list[i]!=list[i-1]) i] )
[ select(list,ind), deltas( concat(ind,[len(list)]) ) ];
// Section: List Iteration Helpers
// Function: idx()
@ -1692,7 +1766,7 @@ function submatrix_set(M,A,m=0,n=0) =
// Arguments:
// v = The list of items to group.
// cnt = The number of items to put in each grouping. Default:2
// dflt = The default value to fill in with is the list is not a multiple of `cnt` items long. Default: 0
// dflt = The default value to fill in with if the list is not a multiple of `cnt` items long. Default: 0
// Example:
// v = [1,2,3,4,5,6];
// a = array_group(v,2) returns [[1,2], [3,4], [5,6]]

2
coords.scad
View file

@ -237,7 +237,7 @@ function project_plane(plane,p) =
let(plane = plane_from_points(plane))
assert(is_def(plane), "Point list is not coplanar")
project_plane(plane)
: assert(is_def(p), str("Invalid plane specification",plane))
: assert(is_def(p), str("Invalid plane specification: ",plane))
is_vnf(p) ? [project_plane(plane,p[0]), p[1]]
: is_list(p) && is_list(p[0]) && is_vector(p[0][0],3) ? // bezier patch or region
[for(plist=p) project_plane(plane,plist)]

2
distributors.scad
View file

@ -34,7 +34,7 @@ module move_copies(a=[[0,0,0]])
assert(is_list(a));
for ($idx = idx(a)) {
$pos = a[$idx];
assert(is_vector($pos));
assert(is_vector($pos),"move_copies offsets should be a 2d or 3d vector.");
translate($pos) children();
}
}

402
geometry.scad
View file

@ -19,15 +19,8 @@
// edge = Array of two points forming the line segment to test against.
// eps = Tolerance in geometric comparisons. Default: `EPSILON` (1e-9)
function point_on_segment2d(point, edge, eps=EPSILON) =
assert( is_vector(point,2), "Invalid point." )
assert( is_finite(eps) && (eps>=0), "The tolerance should be a non-negative value." )
assert( _valid_line(edge,2,eps=eps), "Invalid segment." )
let( dp = point-edge[0],
de = edge[1]-edge[0],
ne = norm(de) )
( dp*de >= -eps*ne )
&& ( (dp-de)*de <= eps*ne ) // point projects on the segment
&& _dist2line(point-edge[0],unit(de))<eps; // point is on the line
point_segment_distance(point, edge)<eps;
//Internal - distance from point `d` to the line passing through the origin with unit direction n
@ -44,7 +37,7 @@ function _point_above_below_segment(point, edge) =
//Internal
function _valid_line(line,dim,eps=EPSILON) =
is_matrix(line,2,dim)
&& ! approx(norm(line[1]-line[0]), 0, eps);
&& norm(line[1]-line[0])>eps*max(norm(line[1]),norm(line[0]));
//Internal
function _valid_plane(p, eps=EPSILON) = is_vector(p,4) && ! approx(norm(p),0,eps);
@ -87,7 +80,7 @@ function collinear(a, b, c, eps=EPSILON) =
// Function: point_line_distance()
// Usage:
// point_line_distance(pt, line);
// point_line_distance(line, pt);
// Description:
// Finds the perpendicular distance of a point `pt` from the line `line`.
// Arguments:
@ -113,6 +106,8 @@ function point_line_distance(pt, line) =
// dist = point_segment_distance([3,8], [[-10,0], [10,0]]); // Returns: 8
// dist2 = point_segment_distance([14,3], [[-10,0], [10,0]]); // Returns: 5
function point_segment_distance(pt, seg) =
assert( is_matrix(concat([pt],seg),3),
"Input should be a point and a valid segment with the dimension equal to the point." )
norm(seg[0]-seg[1]) < EPSILON ? norm(pt-seg[0]) :
norm(pt-segment_closest_point(seg,pt));
@ -129,34 +124,9 @@ function point_segment_distance(pt, seg) =
// dist = segment_distance([[-14,3], [-15,9]], [[-10,0], [10,0]]); // Returns: 5
// dist2 = segment_distance([[-5,5], [5,-5]], [[-10,3], [10,-3]]); // Returns: 0
function segment_distance(seg1, seg2) =
let(
dseg1 = seg1[1]-seg1[0],
dseg2 = seg2[1]-seg2[0],
A = [ [dseg1*dseg1, -dseg1*dseg2],
[-dseg2*dseg1, dseg2*dseg2] ],
b = -[ dseg1, -dseg2 ]*(seg1[0]-seg2[0]),
uv = linear_solve(A,b)
)
!uv ?
norm(dseg1)<EPSILON
? norm(dseg2)<EPSILON
? norm(seg1[0]-seg2[0])
: point_segment_distance(seg1[0],seg2)
: point_segment_distance(seg2[0],seg1) :
uv[0]>=0 && uv[0]<=1 ?
let( p1 = seg1[0] + uv[0]*dseg1 )
uv[1]>=0 && uv[1]<=1
? norm(p1 - (seg2[0] + uv[1]*dseg2) )
: min(norm(p1-seg2[0]), norm(p1-seg2[1])) :
uv[1]>=0 && uv[1]<=1
? let( p2 = seg2[0] + uv[1]*dseg2 )
min(norm(p2-seg1[0]), norm(p2-seg1[1]))
: min(
norm(seg1[0]-seg2[0]),
norm(seg1[0]-seg2[1]),
norm(seg1[1]-seg2[0]),
norm(seg1[1]-seg2[1])
);
assert( is_matrix(concat(seg1,seg2),4),
"Inputs should be two valid segments." )
convex_distance(seg1,seg2);
// Function: line_normal()
@ -494,17 +464,9 @@ function ray_closest_point(ray,pt) =
// color("blue") translate(pt) sphere(r=1,$fn=12);
// color("red") translate(p2) sphere(r=1,$fn=12);
function segment_closest_point(seg,pt) =
assert(_valid_line(seg), "Invalid segment." )
assert(len(pt)==len(seg[0]), "Incompatible dimensions." )
approx(seg[0],seg[1])? seg[0] :
let(
seglen = norm(seg[1]-seg[0]),
segvec = (seg[1]-seg[0])/seglen,
projection = (pt-seg[0]) * segvec
)
projection<=0 ? seg[0] :
projection>=seglen ? seg[1] :
seg[0] + projection*segvec;
assert( is_matrix(concat([pt],seg),3) ,
"Invalid point or segment or incompatible dimensions." )
pt + _closest_s1([seg[0]-pt, seg[1]-pt])[0];
// Function: line_from_points()
@ -512,7 +474,7 @@ function segment_closest_point(seg,pt) =
// line_from_points(points, [fast], [eps]);
// Description:
// Given a list of 2 or more collinear points, returns a line containing them.
// If `fast` is false and the points are coincident, then `undef` is returned.
// If `fast` is false and the points are coincident or non-collinear, then `undef` is returned.
// if `fast` is true, then the collinearity test is skipped and a line passing through 2 distinct arbitrary points is returned.
// Arguments:
// points = The list of points to find the line through.
@ -522,7 +484,7 @@ function line_from_points(points, fast=false, eps=EPSILON) =
assert( is_path(points,dim=undef), "Improper point list." )
assert( is_finite(eps) && (eps>=0), "The tolerance should be a non-negative value." )
let( pb = furthest_point(points[0],points) )
approx(norm(points[pb]-points[0]),0) ? undef :
norm(points[pb]-points[0])<eps*max(norm(points[pb]),norm(points[0])) ? undef :
fast || collinear(points) ? [points[pb], points[0]] : undef;
@ -1171,7 +1133,7 @@ function _general_plane_line_intersection(plane, line, eps=EPSILON) =
// Description:
// Returns a new representation [A,B,C,D] of `plane` where norm([A,B,C]) is equal to one.
function normalize_plane(plane) =
assert( _valid_plane(plane), "Invalid plane." )
assert( _valid_plane(plane), str("Invalid plane. ",plane ) )
plane/norm(point3d(plane));
@ -1179,12 +1141,12 @@ function normalize_plane(plane) =
// Usage:
// angle = plane_line_angle(plane,line);
// Description:
// Compute the angle between a plane [A, B, C, D] and a line, specified as a pair of points [p1,p2].
// Compute the angle between a plane [A, B, C, D] and a 3d line, specified as a pair of 3d points [p1,p2].
// The resulting angle is signed, with the sign positive if the vector p2-p1 lies on
// the same side of the plane as the plane's normal vector.
function plane_line_angle(plane, line) =
assert( _valid_plane(plane), "Invalid plane." )
assert( _valid_line(line), "Invalid line." )
assert( _valid_line(line,dim=3), "Invalid 3d line." )
let(
linedir = unit(line[1]-line[0]),
normal = plane_normal(plane),
@ -1209,7 +1171,7 @@ function plane_line_angle(plane, line) =
// eps = Tolerance in geometric comparisons. Default: `EPSILON` (1e-9)
function plane_line_intersection(plane, line, bounded=false, eps=EPSILON) =
assert( is_finite(eps) && eps>=0, "The tolerance should be a positive number." )
assert(_valid_plane(plane,eps=eps) && _valid_line(line,dim=3,eps=eps), "Invalid plane and/or line.")
assert(_valid_plane(plane,eps=eps) && _valid_line(line,dim=3,eps=eps), "Invalid plane and/or 3d line.")
assert(is_bool(bounded) || is_bool_list(bounded,2), "Invalid bound condition.")
let(
bounded = is_list(bounded)? bounded : [bounded, bounded],
@ -1240,7 +1202,7 @@ function polygon_line_intersection(poly, line, bounded=false, eps=EPSILON) =
assert( is_finite(eps) && eps>=0, "The tolerance should be a positive number." )
assert(is_path(poly,dim=3), "Invalid polygon." )
assert(!is_list(bounded) || len(bounded)==2, "Invalid bound condition(s).")
assert(_valid_line(line,dim=3,eps=eps), "Invalid line." )
assert(_valid_line(line,dim=3,eps=eps), "Invalid 3D line." )
let(
bounded = is_list(bounded)? bounded : [bounded, bounded],
poly = deduplicate(poly),
@ -1694,7 +1656,7 @@ function circle_circle_tangents(c1,r1,c2,r2,d1,d2) =
// eps = epsilon used for identifying the case with one solution. Default: 1e-9
function circle_line_intersection(c,r,d,line,bounded=false,eps=EPSILON) =
let(r=get_radius(r=r,d=d,dflt=undef))
assert(_valid_line(line,2), "Input 'line' is not a valid 2d line.")
assert(_valid_line(line,2), "Invalid 2d line.")
assert(is_vector(c,2), "Circle center must be a 2-vector")
assert(is_num(r) && r>0, "Radius must be positive")
assert(is_bool(bounded) || is_bool_list(bounded,2), "Invalid bound condition")
@ -1738,14 +1700,14 @@ function noncollinear_triple(points,error=true,eps=EPSILON) =
pb = points[b],
nrm = norm(pa-pb)
)
approx(nrm, 0)
nrm <= eps*max(norm(pa),norm(pb))
? assert(!error, "Cannot find three noncollinear points in pointlist.")
[]
: let(
n = (pb-pa)/nrm,
distlist = [for(i=[0:len(points)-1]) _dist2line(points[i]-pa, n)]
)
max(distlist)<eps*nrm
max(distlist) < eps*nrm
? assert(!error, "Cannot find three noncollinear points in pointlist.")
[]
: [0,b,max_index(distlist)];
@ -1761,12 +1723,13 @@ function noncollinear_triple(points,error=true,eps=EPSILON) =
// Arguments:
// pts = List of points.
function pointlist_bounds(pts) =
assert(is_matrix(pts) && len(pts)>0 && len(pts[0])>0 , "Invalid pointlist." )
let(ptsT = transpose(pts))
[
[for(row=ptsT) min(row)],
[for(row=ptsT) max(row)]
];
assert(is_path(pts,dim=undef,fast=true) , "Invalid pointlist." )
let(
select = ident(len(pts[0])),
spread = [for(i=[0:len(pts[0])-1])
let( spreadi = pts*select[i] )
[min(spreadi), max(spreadi)] ] )
transpose(spread);
// Function: closest_point()
@ -1805,7 +1768,7 @@ function furthest_point(pt, points) =
// area = polygon_area(poly);
// Description:
// Given a 2D or 3D planar polygon, returns the area of that polygon.
// If the polygon is self-crossing, the results are undefined. For non-planar 3D polygon the result is [].
// If the polygon is self-crossing, the results are undefined. For non-planar 3D polygon the result is `undef`.
// When `signed` is true, a signed area is returned; a positive area indicates a clockwise polygon.
// Arguments:
// poly = Polygon to compute the area of.
@ -1817,53 +1780,17 @@ function polygon_area(poly, signed=false) =
? let( total = sum([for(i=[1:1:len(poly)-2]) cross(poly[i]-poly[0],poly[i+1]-poly[0]) ])/2 )
signed ? total : abs(total)
: let( plane = plane_from_polygon(poly) )
plane==[]? [] :
plane==[]? undef :
let(
n = plane_normal(plane),
total = sum([
for(i=[1:1:len(poly)-2])
let(
v1 = poly[i] - poly[0],
v2 = poly[i+1] - poly[0]
)
cross(v1,v2)
])* n/2
total =
sum([ for(i=[1:1:len(poly)-2])
cross(poly[i]-poly[0], poly[i+1]-poly[0])
]) * n/2
)
signed ? total : abs(total);
// Function: is_convex_polygon()
// Usage:
// is_convex_polygon(poly);
// Description:
// Returns true if the given 2D or 3D polygon is convex.
// The result is meaningless if the polygon is not simple (self-intersecting) or non coplanar.
// If the points are collinear an error is generated.
// Arguments:
// poly = Polygon to check.
// eps = Tolerance for the collinearity test. Default: EPSILON.
// Example:
// is_convex_polygon(circle(d=50)); // Returns: true
// is_convex_polygon(rot([50,120,30], p=path3d(circle(1,$fn=50)))); // Returns: true
// Example:
// spiral = [for (i=[0:36]) let(a=-i*10) (10+i)*[cos(a),sin(a)]];
// is_convex_polygon(spiral); // Returns: false
function is_convex_polygon(poly,eps=EPSILON) =
assert(is_path(poly), "The input should be a 2D or 3D polygon." )
let( lp = len(poly),
p0 = poly[0] )
assert( lp>=3 , "A polygon must have at least 3 points" )
let( crosses = [for(i=[0:1:lp-1]) cross(poly[(i+1)%lp]-poly[i], poly[(i+2)%lp]-poly[(i+1)%lp]) ] )
len(p0)==2
? assert( !approx(sqrt(max(max(crosses),-min(crosses))),eps), "The points are collinear" )
min(crosses) >=0 || max(crosses)<=0
: let( prod = crosses*sum(crosses),
minc = min(prod),
maxc = max(prod) )
assert( !approx(sqrt(max(maxc,-minc)),eps), "The points are collinear" )
minc>=0 || maxc<=0;
// Function: polygon_shift()
// Usage:
// polygon_shift(poly, i);
@ -2030,9 +1957,9 @@ function centroid(poly, eps=EPSILON) =
// Returns -1 if the point is outside the polygon.
// Returns 0 if the point is on the boundary.
// Returns 1 if the point lies in the interior.
// The polygon does not need to be simple: it can have self-intersections.
// The polygon does not need to be simple: it may have self-intersections.
// But the polygon cannot have holes (it must be simply connected).
// Rounding error may give mixed results for points on or near the boundary.
// Rounding errors may give mixed results for points on or near the boundary.
// Arguments:
// point = The 2D point to check position of.
// poly = The list of 2D path points forming the perimeter of the polygon.
@ -2125,7 +2052,7 @@ function ccw_polygon(poly) =
// poly = The list of the path points for the perimeter of the polygon.
function reverse_polygon(poly) =
assert(is_path(poly), "Input should be a polygon")
let(lp=len(poly)) [for (i=idx(poly)) poly[(lp-i)%lp]];
[poly[0], for(i=[len(poly)-1:-1:1]) poly[i] ];
// Function: polygon_normal()
@ -2133,7 +2060,7 @@ function reverse_polygon(poly) =
// n = polygon_normal(poly);
// Description:
// Given a 3D planar polygon, returns a unit-length normal vector for the
// clockwise orientation of the polygon. If the polygon points are collinear, returns [].
// clockwise orientation of the polygon. If the polygon points are collinear, returns `undef`.
// It doesn't check for coplanarity.
// Arguments:
// poly = The list of 3D path points for the perimeter of the polygon.
@ -2141,7 +2068,7 @@ function polygon_normal(poly) =
assert(is_path(poly,dim=3), "Invalid 3D polygon." )
len(poly)==3 ? point3d(plane3pt(poly[0],poly[1],poly[2])) :
let( triple = sort(noncollinear_triple(poly,error=false)) )
triple==[] ? [] :
triple==[] ? undef :
point3d(plane3pt(poly[triple[0]],poly[triple[1]],poly[triple[2]])) ;
@ -2295,4 +2222,255 @@ function split_polygons_at_each_z(polys, zs, _i=0) =
);
// Section: Convex Sets
// Function: is_convex_polygon()
// Usage:
// is_convex_polygon(poly);
// Description:
// Returns true if the given 2D or 3D polygon is convex.
// The result is meaningless if the polygon is not simple (self-intersecting) or non coplanar.
// If the points are collinear an error is generated.
// Arguments:
// poly = Polygon to check.
// eps = Tolerance for the collinearity test. Default: EPSILON.
// Example:
// is_convex_polygon(circle(d=50)); // Returns: true
// is_convex_polygon(rot([50,120,30], p=path3d(circle(1,$fn=50)))); // Returns: true
// Example:
// spiral = [for (i=[0:36]) let(a=-i*10) (10+i)*[cos(a),sin(a)]];
// is_convex_polygon(spiral); // Returns: false
function is_convex_polygon(poly,eps=EPSILON) =
assert(is_path(poly), "The input should be a 2D or 3D polygon." )
let( lp = len(poly),
p0 = poly[0] )
assert( lp>=3 , "A polygon must have at least 3 points" )
let( crosses = [for(i=[0:1:lp-1]) cross(poly[(i+1)%lp]-poly[i], poly[(i+2)%lp]-poly[(i+1)%lp]) ] )
len(p0)==2
? assert( !approx(sqrt(max(max(crosses),-min(crosses))),eps), "The points are collinear" )
min(crosses) >=0 || max(crosses)<=0
: let( prod = crosses*sum(crosses),
minc = min(prod),
maxc = max(prod) )
assert( !approx(sqrt(max(maxc,-minc)),eps), "The points are collinear" )
minc>=0 || maxc<=0;
// Function: convex_distance()
// Usage:
// convex_distance(pts1, pts2,<eps=>);
// See also:
// convex_collision
// Description:
// Returns the smallest distance between a point in convex hull of `points1`
// and a point in the convex hull of `points2`. All the points in the lists
// should have the same dimension, either 2D or 3D.
// A zero result means the hulls intercept whithin a tolerance `eps`.
// Arguments:
// points1 - first list of 2d or 3d points.
// points2 - second list of 2d or 3d points.
// eps - tolerance in distance evaluations. Default: EPSILON.
// Example(2D):
// pts1 = move([-3,0], p=square(3,center=true));
// pts2 = rot(a=45, p=square(2,center=true));
// pts3 = [ [2,0], [1,2],[3,2], [3,-2], [1,-2] ];
// polygon(pts1);
// polygon(pts2);
// polygon(pts3);
// echo(convex_distance(pts1,pts2)); // Returns: 0.0857864
// echo(convex_distance(pts2,pts3)); // Returns: 0
// Example(3D):
// sphr1 = sphere(2,$fn=10);
// sphr2 = move([4,0,0], p=sphr1);
// sphr3 = move([4.5,0,0], p=sphr1);
// vnf_polyhedron(sphr1);
// vnf_polyhedron(sphr2);
// echo(convex_distance(sphr1[0], sphr2[0])); // Returns: 0
// echo(convex_distance(sphr1[0], sphr3[0])); // Returns: 0.5
function convex_distance(points1, points2, eps=EPSILON) =
assert(is_matrix(points1) && is_matrix(points2,undef,len(points1[0])),
"The input list should be a consistent non empty list of points of same dimension.")
assert(len(points1[0])==2 || len(points1[0])==3 ,
"The input points should be 2d or 3d points.")
let( d = points1[0]-points2[0] )
norm(d)<eps ? 0 :
let( v = _support_diff(points1,points2,-d) )
norm(_GJK_distance(points1, points2, eps, 0, v, [v]));
// Finds the vector difference between the hulls of the two pointsets by the GJK algorithm
// Based on:
// http://www.dtecta.com/papers/jgt98convex.pdf
function _GJK_distance(points1, points2, eps=EPSILON, lbd, d, simplex=[]) =
let( nrd = norm(d) ) // distance upper bound
nrd<eps ? d :
let(
v = _support_diff(points1,points2,-d),
lbd = max(lbd, d*v/nrd), // distance lower bound
close = (nrd-lbd <= eps*nrd)
)
// v already in the simplex is a degenerence due to numerical errors
// and may produce a non-stopping loop
close || [for(nv=norm(v), s=simplex) if(norm(s-v)<=eps*nv) 1]!=[] ? d :
let( newsplx = _closest_simplex(concat(simplex,[v]),eps) )
_GJK_distance(points1, points2, eps, lbd, newsplx[0], newsplx[1]);
// Function: convex_collision()
// Usage:
// convex_collision(pts1, pts2,<eps=>);
// See also:
// convex_distance
// Description:
// Returns `true` if the convex hull of `points1` intercepts the convex hull of `points2`
// otherwise, `false`.
// All the points in the lists should have the same dimension, either 2D or 3D.
// This function is tipically faster than `convex_distance` to find a non-collision.
// Arguments:
// points1 - first list of 2d or 3d points.
// points2 - second list of 2d or 3d points.
// eps - tolerance for the intersection tests. Default: EPSILON.
// Example(2D):
// pts1 = move([-3,0], p=square(3,center=true));
// pts2 = rot(a=45, p=square(2,center=true));
// pts3 = [ [2,0], [1,2],[3,2], [3,-2], [1,-2] ];
// polygon(pts1);
// polygon(pts2);
// polygon(pts3);
// echo(convex_collision(pts1,pts2)); // Returns: false
// echo(convex_collision(pts2,pts3)); // Returns: true
// Example(3D):
// sphr1 = sphere(2,$fn=10);
// sphr2 = move([4,0,0], p=sphr1);
// sphr3 = move([4.5,0,0], p=sphr1);
// vnf_polyhedron(sphr1);
// vnf_polyhedron(sphr2);
// echo(convex_collision(sphr1[0], sphr2[0])); // Returns: true
// echo(convex_collision(sphr1[0], sphr3[0])); // Returns: false
//
function convex_collision(points1, points2, eps=EPSILON) =
assert(is_matrix(points1) && is_matrix(points2,undef,len(points1[0])),
"The input list should be a consistent non empty list of points of same dimension.")
assert(len(points1[0])==2 || len(points1[0])==3 ,
"The input points should be 2d or 3d points.")
let( d = points1[0]-points2[0] )
norm(d)<eps ? true :
let( v = _support_diff(points1,points2,-d) )
_GJK_collide(points1, points2, v, [v], eps);
// Based on the GJK collision algorithms found in:
// http://uu.diva-portal.org/smash/get/diva2/FFULLTEXT01.pdf
// or
// http://www.dtecta.com/papers/jgt98convex.pdf
function _GJK_collide(points1, points2, d, simplex, eps=EPSILON) =
norm(d) < eps ? true : // does collide
let( v = _support_diff(points1,points2,-d) )
v*d > eps ? false : // no collision
let( newsplx = _closest_simplex(concat(simplex,[v]),eps) )
_GJK_collide(points1, points2, newsplx[0], newsplx[1], eps);
// given a simplex s, returns a pair:
// - the point of the s closest to the origin
// - the smallest sub-simplex of s that contains that point
function _closest_simplex(s,eps=EPSILON) =
assert(len(s)>=2 && len(s)<=4, "Internal error.")
len(s)==2 ? _closest_s1(s,eps) :
len(s)==3 ? _closest_s2(s,eps)
: _closest_s3(s,eps);
// find the closest to a 1-simplex
// Based on: http://uu.diva-portal.org/smash/get/diva2/FFULLTEXT01.pdf
function _closest_s1(s,eps=EPSILON) =
norm(s[1]-s[0])<eps*(norm(s[0])+norm(s[1]))/2 ? [ s[0], [s[0]] ] :
let(
c = s[1]-s[0],
t = -s[0]*c/(c*c)
)
t<0 ? [ s[0], [s[0]] ] :
t>1 ? [ s[1], [s[1]] ] :
[ s[0]+t*c, s ];
// find the closest to a 2-simplex
// Based on: http://uu.diva-portal.org/smash/get/diva2/FFULLTEXT01.pdf
function _closest_s2(s,eps=EPSILON) =
let(
dim = len(s[0]),
a = dim==3 ? s[0]: [ each s[0], 0] ,
b = dim==3 ? s[1]: [ each s[1], 0] ,
c = dim==3 ? s[2]: [ each s[2], 0] ,
ab = norm(a-b),
bc = norm(b-c),
ca = norm(c-a),
nr = cross(b-a,c-a)
)
norm(nr) <= eps*max(ab,bc,ca) // degenerate case
? let( i = max_index([ab, bc, ca]) )
_closest_s1([s[i],s[(i+1)%3]],eps)
// considering that s[2] was the last inserted vertex in s,
// the only possible outcomes are :
// s, [s[0],s[2]] and [s[1],s[2]]
: let(
class = (cross(nr,a-b)*a<0 ? 1 : 0 )
+ (cross(nr,c-a)*a<0 ? 2 : 0 )
+ (cross(nr,b-c)*b<0 ? 4 : 0 )
)
assert( class!=1, "Internal error" )
class==0 ? [ nr*(nr*a)/(nr*nr), s] : // origin projects (or is) on the tri
// class==1 ? _closest_s1([s[0],s[1]]) :
class==2 ? _closest_s1([s[0],s[2]],eps) :
class==4 ? _closest_s1([s[1],s[2]],eps) :
// class==3 ? a*(a-b)> 0 ? _closest_s1([s[0],s[1]]) : _closest_s1([s[0],s[2]]) :
class==3 ? _closest_s1([s[0],s[2]],eps) :
// class==5 ? b*(b-c)<=0 ? _closest_s1([s[0],s[1]]) : _closest_s1([s[1],s[2]]) :
class==5 ? _closest_s1([s[1],s[2]],eps) :
c*(c-a)>0 ? _closest_s1([s[0],s[2]],eps) : _closest_s1([s[1],s[2]],eps);
// find the closest to a 3-simplex
// it seems that degenerate 3-simplices are correctly manage without extra code
function _closest_s3(s,eps=EPSILON) =
assert( len(s[0])==3 && len(s)==4, "Internal error." )
let( nr = cross(s[1]-s[0],s[2]-s[0]),
sz = [ norm(s[1]-s[0]), norm(s[1]-s[2]), norm(s[2]-s[0]) ] )
norm(nr)<eps*max(sz)
? let( i = max_index(sz) )
_closest_s2([ s[i], s[(i+1)%3], s[3] ], eps) // degenerate case
// considering that s[3] was the last inserted vertex in s,
// the only possible outcomes will be:
// s or some of the 3 triangles of s containing s[3]
: let(
tris = [ [s[0], s[1], s[3]],
[s[1], s[2], s[3]],
[s[2], s[0], s[3]] ],
cntr = sum(s)/4,
// indicator of the tris facing the origin
facing = [for(i=[0:2])
let( nrm = _tri_normal(tris[i]) )
if( ((nrm*(s[i]-cntr))>0)==(nrm*s[i]<0) ) i ]
)
len(facing)==0 ? [ [0,0,0], s ] : // origin is inside the simplex
len(facing)==1 ? _closest_s2(tris[facing[0]], eps) :
let( // look for the origin-facing tri closest to the origin
closest = [for(i=facing) _closest_s2(tris[i], eps) ],
dist = [for(cl=closest) norm(cl[0]) ],
nearest = min_index(dist)
)
closest[nearest];
function _tri_normal(tri) = cross(tri[1]-tri[0],tri[2]-tri[0]);
function _support_diff(p1,p2,d) =
let( p1d = p1*d, p2d = p2*d )
p1[search(max(p1d),p1d,1)[0]] - p2[search(min(p2d),p2d,1)[0]];
// vim: expandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap

14
hull.scad
View file

@ -202,17 +202,15 @@ function hull3d_faces(points) =
// Adds the remaining points one by one to the convex hull
function _hull3d_iterative(points, triangles, planes, remaining, _i=0) =
function _hull3d_iterative(points, triangles, planes, remaining, _i=0) = //let( EPSILON=1e-12 )
_i >= len(remaining) ? triangles :
let (
// pick a point
i = remaining[_i],
// evaluate the triangle plane equations at point i
// xx=[for(i=[0:len(planes)-1],p=[planes[i]]) if(len(p)!=4) echo(i=i,len(p))0 ],//echo([each points[i], -1]),
// planeq_val = [for(p=planes) p*[each points[i], -1]],
planeq_val = planes*[each points[i], -1],
// find the triangles that are in conflict with the point (point not inside)
conflicts = [ for (i = [0:1:len(planes)-1]) if (planeq_val[i]>EPSILON) i ],
conflicts = [for (i = [0:1:len(planeq_val)-1]) if (planeq_val[i]>EPSILON) i ],
// collect the halfedges of all triangles that are in conflict
halfedges = [
for(c = conflicts, i = [0:2])
@ -222,18 +220,15 @@ function _hull3d_iterative(points, triangles, planes, remaining, _i=0) =
horizon = _remove_internal_edges(halfedges),
// generate new triangles connecting point i to each horizon halfedge vertices
tri2add = [ for (h = horizon) concat(h,i) ],
// w=[for(t=tri2add) if(collinear(points[t[0]],points[t[1]],points[t[2]])) echo(t)],
// add tria2add and remove conflict triangles
new_triangles =
concat( tri2add,
[ for (i = [0:1:len(planes)-1]) if (planeq_val[i]<=EPSILON) triangles[i] ]
),
// y=[for(t=tri2add) if([]==plane3pt_indexed(points,t[0],t[1],t[2])) echo(tri2add=t,pts=[for(ti=t) points[ti]])],
// add the plane equations of new added triangles and remove the plane equations of the conflict ones
new_planes =
concat( [ for (t = tri2add) plane3pt_indexed(points, t[0], t[1], t[2]) ],
[ for (i = [0:1:len(planes)-1]) if (planeq_val[i]<=EPSILON) planes[i] ]
)
[ for (t = tri2add) plane3pt_indexed(points, t[0], t[1], t[2]) ,
for (i = [0:1:len(planes)-1]) if (planeq_val[i]<=EPSILON) planes[i] ]
) _hull3d_iterative(
points,
new_triangles,
@ -249,7 +244,6 @@ function _remove_internal_edges(halfedges) = [
h
];
function _find_first_noncoplanar(plane, points, i=0) =
(i >= len(points) || !points_on_plane([points[i]],plane))? i :
_find_first_noncoplanar(plane, points, i+1);

38
math.scad
View file

@ -640,17 +640,19 @@ function sum_of_sines(a, sines) =
// Usage:
// delts = deltas(v);
// Description:
// Returns a list with the deltas of adjacent entries in the given list.
// Returns a list with the deltas of adjacent entries in the given list, optionally wrapping back to the front.
// The list should be a consistent list of numeric components (numbers, vectors, matrix, etc).
// Given [a,b,c,d], returns [b-a,c-b,d-c].
//
// Arguments:
// v = The list to get the deltas of.
// wrap = If true, wrap back to the start from the end. ie: return the difference between the last and first items as the last delta. Default: false
// Example:
// deltas([2,5,9,17]); // returns [3,4,8].
// deltas([[1,2,3], [3,6,8], [4,8,11]]); // returns [[2,4,5], [1,2,3]]
function deltas(v) =
function deltas(v, wrap=false) =
assert( is_consistent(v) && len(v)>1 , "Inconsistent list or with length<=1.")
[for (p=pair(v)) p[1]-p[0]] ;
[for (p=pair(v,wrap)) p[1]-p[0]] ;
// Function: product()
@ -771,21 +773,31 @@ function _med3(a,b,c) =
// Usage:
// x = convolve(p,q);
// Description:
// Given two vectors, finds the convolution of them.
// The length of the returned vector is len(p)+len(q)-1 .
// Given two vectors, or one vector and a path or
// two paths of the same dimension, finds the convolution of them.
// If both parameter are vectors, returns the vector convolution.
// If one parameter is a vector and the other a path,
// convolves using products by scalars and returns a path.
// If both parameters are paths, convolve using scalar products
// and returns a vector.
// The returned vector or path has length len(p)+len(q)-1.
// Arguments:
// p = The first vector.
// q = The second vector.
// p = The first vector or path.
// q = The second vector or path.
// Example:
// a = convolve([1,1],[1,2,1]); // Returns: [1,3,3,1]
// b = convolve([1,2,3],[1,2,1])); // Returns: [1,4,8,8,3]
// c = convolve([[1,1],[2,2],[3,1]],[1,2,1])); // Returns: [[1,1],[4,4],[8,6],[8,4],[3,1]]
// d = convolve([[1,1],[2,2],[3,1]],[[1,2],[2,1]])); // Returns: [3,9,11,7]
function convolve(p,q) =
p==[] || q==[] ? [] :
assert( is_vector(p) && is_vector(q), "The inputs should be vectors.")
assert( (is_vector(p) || is_matrix(p))
&& ( is_vector(q) || (is_matrix(q) && ( !is_vector(p[0]) || (len(p[0])==len(q[0])) ) ) ) ,
"The inputs should be vectors or paths all of the same dimension.")
let( n = len(p),
m = len(q))
[for(i=[0:n+m-2], k1 = max(0,i-n+1), k2 = min(i,m-1) )
[for(j=[k1:k2]) p[i-j] ] * [for(j=[k1:k2]) q[j] ]
sum([for(j=[k1:k2]) p[i-j]*q[j] ])
];
@ -1694,7 +1706,7 @@ function polynomial(p,z,k,total) =
// x = polymult(p,q)
// x = polymult([p1,p2,p3,...])
// Description:
// Given a list of polynomials represented as real coefficient lists, with the highest degree coefficient first,
// Given a list of polynomials represented as real algebraic coefficient lists, with the highest degree coefficient first,
// computes the coefficient list of the product polynomial.
function poly_mult(p,q) =
is_undef(q) ?
@ -1714,8 +1726,8 @@ function poly_mult(p,q) =
// Description:
// Computes division of the numerator polynomial by the denominator polynomial and returns
// a list of two polynomials, [quotient, remainder]. If the division has no remainder then
// the zero polynomial [] is returned for the remainder. Similarly if the quotient is zero
// the returned quotient will be [].
// the zero polynomial [0] is returned for the remainder. Similarly if the quotient is zero
// the returned quotient will be [0].
function poly_div(n,d) =
assert( is_vector(n) && is_vector(d) , "Invalid polynomials." )
let( d = _poly_trim(d),
@ -1740,7 +1752,7 @@ function _poly_div(n,d,q) =
/// _poly_trim(p,[eps])
/// Description:
/// Removes leading zero terms of a polynomial. By default zeros must be exact,
/// or give epsilon for approximate zeros.
/// or give epsilon for approximate zeros. Returns [0] for a zero polynomial.
function _poly_trim(p,eps=0) =
let( nz = [for(i=[0:1:len(p)-1]) if ( !approx(p[i],0,eps)) i])
len(nz)==0 ? [0] : list_tail(p,nz[0]);

3
paths.scad
View file

@ -997,7 +997,7 @@ module jittered_poly(path, dist=1/512) {
// Module: extrude_from_to()
// Description:
// Extrudes a 2D shape between the points pt1 and pt2. Takes as children a set of 2D shapes to extrude.
// Extrudes a 2D shape between the 3d points pt1 and pt2. Takes as children a set of 2D shapes to extrude.
// Arguments:
// pt1 = starting point of extrusion.
// pt2 = ending point of extrusion.
@ -1010,6 +1010,7 @@ module jittered_poly(path, dist=1/512) {
// xcopies(3) circle(3, $fn=32);
// }
module extrude_from_to(pt1, pt2, convexity, twist, scale, slices) {
assert( is_path([pt1,pt2],3), "The points should be 3d points");
rtp = xyz_to_spherical(pt2-pt1);
translate(pt1) {
rotate([0, rtp[2], rtp[1]]) {

16
tests/test_arrays.scad
View file

@ -358,11 +358,21 @@ module test_sortidx() {
}
test_sortidx();
module test_group_sort() {
assert_equal(group_sort([]), [[]]);
assert_equal(group_sort([8]), [[8]]);
assert_equal(group_sort([7,3,9,4,3,1,8]), [[1], [3, 3], [4], [7], [8], [9]]);
assert_equal(group_sort([[5,"a"],[2,"b"], [5,"c"], [3,"d"], [2,"e"] ], idx=0), [[[2, "b"], [2, "e"]], [[3, "d"]], [[5, "a"], [5, "c"]]]);
assert_equal(group_sort([["a",5],["b",6], ["c",1], ["d",2], ["e",6] ], idx=1), [[["c", 1]], [["d", 2]], [["a", 5]], [["b", 6], ["e", 6]]] );
}
test_group_sort();
module test_unique() {
assert(unique([]) == []);
assert(unique([8]) == [8]);
assert(unique([7,3,9,4,3,1,8]) == [1,3,4,7,8,9]);
assert_equal(unique([]), []);
assert_equal(unique([8]), [8]);
assert_equal(unique([7,3,9,4,3,1,8]), [1,3,4,7,8,9]);
assert_equal(unique(["A","B","R","A","C","A","D","A","B","R","A"]), ["A", "B", "C", "D", "R"]);
}
test_unique();

43
tests/test_geometry.scad
View file

@ -90,7 +90,8 @@ test_cleanup_path();
test_simplify_path();
test_simplify_path_indexed();
test_is_region();
test_convex_distance();
test_convex_collision();
// to be used when there are two alternative symmetrical outcomes
// from a function like a plane output; v must be a vector
@ -301,7 +302,7 @@ module test_line_from_points() {
}
*test_line_from_points();
module test_point_on_segment2d() {
module test_point_on_segment2d() {
assert(point_on_segment2d([-15,0], [[-10,0], [10,0]]) == false);
assert(point_on_segment2d([-10,0], [[-10,0], [10,0]]) == true);
assert(point_on_segment2d([-5,0], [[-10,0], [10,0]]) == true);
@ -1075,6 +1076,44 @@ module test_is_region() {
}
*test_is_region();
module test_convex_distance() {
// 2D
c1 = circle(10,$fn=24);
c2 = move([15,0], p=c1);
assert(convex_distance(c1, c2)==0);
c3 = move([22,0],c1);
assert_approx(convex_distance(c1, c3),2);
// 3D
s1 = sphere(10,$fn=4);
s2 = move([15,0], p=s1);
assert_approx(convex_distance(s1[0], s2[0]), 0.857864376269);
s3 = move([25.3,0],s1);
assert_approx(convex_distance(s1[0], s3[0]), 11.1578643763);
s4 = move([30,25],s1);
assert_approx(convex_distance(s1[0], s4[0]), 28.8908729653);
s5 = move([10*sqrt(2),0],s1);
assert_approx(convex_distance(s1[0], s5[0]), 0);
}
*test_convex_distance();
module test_convex_collision() {
// 2D
c1 = circle(10,$fn=24);
c2 = move([15,0], p=c1);
assert(convex_collision(c1, c2));
c3 = move([22,0],c1);
assert(!convex_collision(c1, c3));
// 3D
s1 = sphere(10,$fn=4);
s2 = move([15,0], p=s1);
assert(!convex_collision(s1[0], s2[0]));
s3 = move([25.3,0],s1);
assert(!convex_collision(s1[0], s3[0]));
s4 = move([5,0],s1);
assert(convex_collision(s1[0], s4[0]));
s5 = move([10*sqrt(2),0],s1);
assert(convex_collision(s1[0], s5[0]));
}
*test_convex_distance();
// vim: expandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap

6
tests/test_math.scad
View file

@ -546,7 +546,9 @@ test_sum_of_sines();
module test_deltas() {
assert_equal(deltas([2,5,9,17]), [3,4,8]);
assert_equal(deltas([2,5,9,17],wrap=true), [3,4,8,-15]);
assert_equal(deltas([[1,2,3], [3,6,8], [4,8,11]]), [[2,4,5], [1,2,3]]);
assert_equal(deltas([[1,2,3], [3,6,8], [4,8,11]],wrap=true), [[2,4,5], [1,2,3], [-3,-6,-8]]);
}
test_deltas();
@ -582,6 +584,10 @@ module test_convolve() {
assert_equal(convolve([1,1],[]), []);
assert_equal(convolve([1,1],[1,2,1]), [1,3,3,1]);
assert_equal(convolve([1,2,3],[1,2,1]), [1,4,8,8,3]);
assert_equal(convolve([1,2,3],[[1],[2],[1]]), [[1], [4], [8], [8], [3]]);
assert_equal(convolve([[1],[2],[3]],[[1],[2],[1]]), [1,4,8,8,3]);
assert_equal(convolve([[1,0],[2,1],[3,2]],[[1,0],[2,1],[1,2]]), [1,4,9,12,7]);
assert_equal(convolve([1,2,3],[[1,0],[2,1],[1,2]]), [[1,0],[4,1],[8,4],[8,7],[3,6]]);
}
test_convolve();