grey/BOSL2
1
0
Fork 0
mirror of https://github.com/BelfrySCAD/BOSL2.git synced 2025-01-01 09:49:45 +00:00
Code Issues Projects Releases Packages Wiki Activity Actions

Merge pull request #709 from adrianVmariano/master

Browse source 
split off linalg.scad
This commit is contained in:
Revar Desmera 2021-10-26 15:48:47 -07:00 committed by GitHub
commit 856af42d7f
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
19 changed files with 820 additions and 1000 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

417
arrays.scad
View file

@ -5,7 +5,6 @@
// include <BOSL2/std.scad>
//////////////////////////////////////////////////////////////////////
// Terminology:
// **List** = An ordered collection of zero or more items. ie: `["a", "b", "c"]`
// **Vector** = A list of numbers. ie: `[4, 5, 6]`
@ -251,7 +250,7 @@ function __find_approx(val, list, eps, i=0) =
// list = The list to get the portion of.
// 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(), columns(), last()
// See Also: slice(), column(), last()
// Example:
// l = [3,4,5,6,7,8,9];
// a = select(l, 5, 6); // Returns [8,9]
@ -292,7 +291,7 @@ function select(list, start, end) =
// list = The list to get the slice of.
// s = The index of the first item to return.
// e = The index of the last item to return.
// See Also: select(), columns(), last()
// See Also: select(), column(), last()
// Example:
// a = slice([3,4,5,6,7,8,9], 3, 5); // Returns [6,7,8]
// b = slice([3,4,5,6,7,8,9], 2, -1); // Returns [5,6,7,8,9]
@ -317,7 +316,7 @@ function slice(list,s=0,e=-1) =
// Usage:
// item = last(list);
// Topics: List Handling
// See Also: select(), slice(), columns()
// See Also: select(), slice(), column()
// Description:
// Returns the last element of a list, or undef if empty.
// Arguments:
@ -1471,63 +1470,6 @@ function permutations(l,n=2) =
: [for (i=idx(l), p=permutations([for (j=idx(l)) if (i!=j) l[j]], n=n-1)) concat([l[i]], p)];
// Function: zip()
// Usage:
// pairs = zip(a,b);
// triples = zip(a,b,c);
// quads = zip([LIST1,LIST2,LIST3,LIST4]);
// Topics: List Handling, Iteration
// See Also: zip_long()
// Description:
// Zips together two or more lists into a single list. For example, if you have two
// lists [3,4,5], and [8,7,6], and zip them together, you get [ [3,8],[4,7],[5,6] ].
// The list returned will be as long as the shortest list passed to zip().
// Arguments:
// a = The first list, or a list of lists if b and c are not given.
// b = The second list, if given.
// c = The third list, if given.
// Example:
// a = [9,8,7,6]; b = [1,2,3];
// for (p=zip(a,b)) echo(p);
// // ECHO: [9,1]
// // ECHO: [8,2]
// // ECHO: [7,3]
function zip(a,b,c) =
b!=undef? zip([a,b,if (c!=undef) c]) :
let(n = min_length(a))
[for (i=[0:1:n-1]) [for (x=a) x[i]]];
// Function: zip_long()
// Usage:
// pairs = zip_long(a,b);
// triples = zip_long(a,b,c);
// quads = zip_long([LIST1,LIST2,LIST3,LIST4]);
// Topics: List Handling, Iteration
// See Also: zip()
// Description:
// Zips together two or more lists into a single list. For example, if you have two
// lists [3,4,5], and [8,7,6], and zip them together, you get [ [3,8],[4,7],[5,6] ].
// The list returned will be as long as the longest list passed to zip_long(), with
// shorter lists padded by the value in `fill`.
// Arguments:
// a = The first list, or a list of lists if b and c are not given.
// b = The second list, if given.
// c = The third list, if given.
// fill = The value to pad shorter lists with. Default: undef
// Example:
// a = [9,8,7,6]; b = [1,2,3];
// for (p=zip_long(a,b,fill=88)) echo(p);
// // ECHO: [9,1]
// // ECHO: [8,2]
// // ECHO: [7,3]
// // ECHO: [6,88]]
function zip_long(a,b,c,fill) =
b!=undef? zip_long([a,b,if (c!=undef) c],fill=fill) :
let(n = max_length(a))
[for (i=[0:1:n-1]) [for (x=a) i<len(x)? x[i] : fill]];
// Section: Set Manipulation
@ -1616,209 +1558,7 @@ function set_intersection(a, b) =
// Section: Array Manipulation
// Function: columns()
// Usage:
// list = columns(M, idx);
// Topics: Array Handling, List Handling
// See Also: select(), slice()
// Description:
// Extracts the entries listed in idx from each entry in M. For a matrix this means
// selecting a specified set of columns. If idx is a number the return is a vector,
// otherwise it is a list of lists (the submatrix).
// This function will return `undef` at all entry positions indexed by idx not found in the input list M.
// Arguments:
// M = The given list of lists.
// idx = The index, list of indices, or range of indices to fetch.
// Example:
// M = [[1,2,3,4],[5,6,7,8],[9,10,11,12],[13,14,15,16]];
// a = columns(M,2); // Returns [3, 7, 11, 15]
// b = columns(M,[2]); // Returns [[3], [7], [11], [15]]
// c = columns(M,[2,1]); // Returns [[3, 2], [7, 6], [11, 10], [15, 14]]
// d = columns(M,[1:3]); // Returns [[2, 3, 4], [6, 7, 8], [10, 11, 12], [14, 15, 16]]
// N = [ [1,2], [3], [4,5], [6,7,8] ];
// e = columns(N,[0,1]); // Returns [ [1,2], [3,undef], [4,5], [6,7] ]
function columns(M, idx) =
assert( is_list(M), "The input is not a list." )
assert( !is_undef(idx) && _valid_idx(idx,0,1/0), "Invalid index input." )
is_finite(idx)
? [for(row=M) row[idx]]
: [for(row=M) [for(i=idx) row[i]]];
// Function: submatrix()
// Usage:
// mat = submatrix(M, idx1, idx2);
// Topics: Matrices, Array Handling
// See Also: columns(), block_matrix(), submatrix_set()
// Description:
// The input must be a list of lists (a matrix or 2d array). Returns a submatrix by selecting the rows listed in idx1 and columns listed in idx2.
// Arguments:
// M = Given list of lists
// idx1 = rows index list or range
// idx2 = column index list or range
// Example:
// M = [[ 1, 2, 3, 4, 5],
// [ 6, 7, 8, 9,10],
// [11,12,13,14,15],
// [16,17,18,19,20],
// [21,22,23,24,25]];
// submatrix(M,[1:2],[3:4]); // Returns [[9, 10], [14, 15]]
// submatrix(M,[1], [3,4])); // Returns [[9,10]]
// submatrix(M,1, [3,4])); // Returns [[9,10]]
// submatrix(M,1,3)); // Returns [[9]]
// submatrix(M, [3,4],1); // Returns [[17],[22]]);
// submatrix(M, [1,3],[2,4]); // Returns [[8,10],[18,20]]);
// A = [[true, 17, "test"],
// [[4,2], 91, false],
// [6, [3,4], undef]];
// submatrix(A,[0,2],[1,2]); // Returns [[17, "test"], [[3, 4], undef]]
function submatrix(M,idx1,idx2) =
[for(i=idx1) [for(j=idx2) M[i][j] ] ];
// Function: hstack()
// Usage:
// A = hstack(M1, M2)
// A = hstack(M1, M2, M3)
// A = hstack([M1, M2, M3, ...])
// Topics: Matrices, Array Handling
// See Also: columns(), submatrix(), block_matrix()
// Description:
// Constructs a matrix by horizontally "stacking" together compatible matrices or vectors. Vectors are treated as columsn in the stack.
// This command is the inverse of `columns`. Note: strings given in vectors are broken apart into lists of characters. Strings given
// in matrices are preserved as strings. If you need to combine vectors of strings use array_group as shown below to convert the
// vector into a column matrix. Also note that vertical stacking can be done directly with concat.
// Arguments:
// M1 = If given with other arguments, the first matrix (or vector) to stack. If given alone, a list of matrices/vectors to stack.
// M2 = Second matrix/vector to stack
// M3 = Third matrix/vector to stack.
// Example:
// M = ident(3);
// v1 = [2,3,4];
// v2 = [5,6,7];
// v3 = [8,9,10];
// a = hstack(v1,v2); // Returns [[2, 5], [3, 6], [4, 7]]
// b = hstack(v1,v2,v3); // Returns [[2, 5, 8],
// // [3, 6, 9],
// // [4, 7, 10]]
// c = hstack([M,v1,M]); // Returns [[1, 0, 0, 2, 1, 0, 0],
// // [0, 1, 0, 3, 0, 1, 0],
// // [0, 0, 1, 4, 0, 0, 1]]
// d = hstack(columns(M,0), columns(M,[1 2])); // Returns M
// strvec = ["one","two"];
// strmat = [["three","four"], ["five","six"]];
// e = hstack(strvec,strvec); // Returns [["o", "n", "e", "o", "n", "e"],
// // ["t", "w", "o", "t", "w", "o"]]
// f = hstack(array_group(strvec,1), array_group(strvec,1));
// // Returns [["one", "one"],
// // ["two", "two"]]
// g = hstack(strmat,strmat); // Returns: [["three", "four", "three", "four"],
// // [ "five", "six", "five", "six"]]
function hstack(M1, M2, M3) =
(M3!=undef)? hstack([M1,M2,M3]) :
(M2!=undef)? hstack([M1,M2]) :
assert(all([for(v=M1) is_list(v)]), "One of the inputs to hstack is not a list")
let(
minlen = min_length(M1),
maxlen = max_length(M1)
)
assert(minlen==maxlen, "Input vectors to hstack must have the same length")
[for(row=[0:1:minlen-1])
[for(matrix=M1)
each matrix[row]
]
];
// Function: block_matrix()
// Usage:
// bmat = block_matrix([[M11, M12,...],[M21, M22,...], ... ]);
// Topics: Matrices, Array Handling
// See Also: columns(), submatrix()
// Description:
// Create a block matrix by supplying a matrix of matrices, which will
// be combined into one unified matrix. Every matrix in one row
// must have the same height, and the combined width of the matrices
// in each row must be equal. Strings will stay strings.
// Example:
// A = [[1,2],
// [3,4]];
// B = ident(2);
// C = block_matrix([[A,B],[B,A],[A,B]]);
// // Returns:
// // [[1, 2, 1, 0],
// // [3, 4, 0, 1],
// // [1, 0, 1, 2],
// // [0, 1, 3, 4],
// // [1, 2, 1, 0],
// // [3, 4, 0, 1]]);
// D = block_matrix([[A,B], ident(4)]);
// // Returns:
// // [[1, 2, 1, 0],
// // [3, 4, 0, 1],
// // [1, 0, 0, 0],
// // [0, 1, 0, 0],
// // [0, 0, 1, 0],
// // [0, 0, 0, 1]]);
// E = [["one", "two"], [3,4]];
// F = block_matrix([[E,E]]);
// // Returns:
// // [["one", "two", "one", "two"],
// // [ 3, 4, 3, 4]]
function block_matrix(M) =
let(
bigM = [for(bigrow = M) each hstack(bigrow)],
len0 = len(bigM[0]),
badrows = [for(row=bigM) if (len(row)!=len0) 1]
)
assert(badrows==[], "Inconsistent or invalid input")
bigM;
// Function: diagonal_matrix()
// Usage:
// mat = diagonal_matrix(diag, [offdiag]);
// Topics: Matrices, Array Handling
// See Also: columns(), submatrix()
// Description:
// Creates a square matrix with the items in the list `diag` on
// its diagonal. The off diagonal entries are set to offdiag,
// which is zero by default.
// Arguments:
// diag = A list of items to put in the diagnal cells of the matrix.
// offdiag = Value to put in non-diagonal matrix cells.
function diagonal_matrix(diag, offdiag=0) =
assert(is_list(diag) && len(diag)>0)
[for(i=[0:1:len(diag)-1]) [for(j=[0:len(diag)-1]) i==j?diag[i] : offdiag]];
// Function: submatrix_set()
// Usage:
// mat = submatrix_set(M, A, [m], [n]);
// Topics: Matrices, Array Handling
// See Also: columns(), submatrix()
// Description:
// Sets a submatrix of M equal to the matrix A. By default the top left corner of M is set to A, but
// you can specify offset coordinates m and n. If A (as adjusted by m and n) extends beyond the bounds
// of M then the extra entries are ignored. You can pass in A=[[]], a null matrix, and M will be
// returned unchanged. Note that the input M need not be rectangular in shape.
// Arguments:
// M = Original matrix.
// A = Sub-matrix of parts to set.
// m = Row number of upper-left corner to place A at.
// n = Column number of upper-left corner to place A at.
function submatrix_set(M,A,m=0,n=0) =
assert(is_list(M))
assert(is_list(A))
assert(is_int(m))
assert(is_int(n))
let( badrows = [for(i=idx(A)) if (!is_list(A[i])) i])
assert(badrows==[], str("Input submatrix malformed rows: ",badrows))
[for(i=[0:1:len(M)-1])
assert(is_list(M[i]), str("Row ",i," of input matrix is not a list"))
[for(j=[0:1:len(M[i])-1])
i>=m && i <len(A)+m && j>=n && j<len(A[0])+n ? A[i-m][j-n] : M[i][j]]];
// Section: Changing list structure
// Function: array_group()
@ -1828,7 +1568,7 @@ function submatrix_set(M,A,m=0,n=0) =
// Takes a flat array of values, and groups items in sets of `cnt` length.
// The opposite of this is `flatten()`.
// Topics: Matrices, Array Handling
// See Also: columns(), submatrix(), hstack(), flatten(), full_flatten()
// See Also: column(), submatrix(), hstack(), flatten(), full_flatten()
// Arguments:
// v = The list of items to group.
// cnt = The number of items to put in each grouping. Default:2
@ -1847,7 +1587,7 @@ function array_group(v, cnt=2, dflt=0) =
// Usage:
// list = flatten(l);
// Topics: Matrices, Array Handling
// See Also: columns(), submatrix(), hstack(), full_flatten()
// See Also: column(), submatrix(), hstack(), full_flatten()
// Description:
// Takes a list of lists and flattens it by one level.
// Arguments:
@ -1863,7 +1603,7 @@ function flatten(l) =
// Usage:
// list = full_flatten(l);
// Topics: Matrices, Array Handling
// See Also: columns(), submatrix(), hstack(), flatten()
// See Also: column(), submatrix(), hstack(), flatten()
// Description:
// Collects in a list all elements recursively found in any level of the given list.
// The output list is ordered in depth first order.
@ -1925,110 +1665,63 @@ function array_dim(v, depth=undef) =
: let( dimlist = _array_dim_recurse(v))
(depth > len(dimlist))? 0 : dimlist[depth-1] ;
// Function: transpose()
// Function: zip()
// Usage:
// arr = transpose(arr, [reverse]);
// Topics: Matrices, Array Handling
// See Also: submatrix(), block_matrix(), hstack(), flatten()
// pairs = zip(a,b);
// triples = zip(a,b,c);
// quads = zip([LIST1,LIST2,LIST3,LIST4]);
// Topics: List Handling, Iteration
// See Also: zip_long()
// Description:
// Returns the transpose of the given input array. The input should be a list of lists that are
// all the same length. If you give a vector then transpose returns it unchanged.
// When reverse=true, the transpose is done across to the secondary diagonal. (See example below.)
// By default, reverse=false.
// Example:
// arr = [
// ["a", "b", "c"],
// ["d", "e", "f"],
// ["g", "h", "i"]
// ];
// t = transpose(arr);
// // Returns:
// // [
// // ["a", "d", "g"],
// // ["b", "e", "h"],
// // ["c", "f", "i"],
// // ]
// Example:
// arr = [
// ["a", "b", "c"],
// ["d", "e", "f"]
// ];
// t = transpose(arr);
// // Returns:
// // [
// // ["a", "d"],
// // ["b", "e"],
// // ["c", "f"],
// // ]
// Example:
// arr = [
// ["a", "b", "c"],
// ["d", "e", "f"],
// ["g", "h", "i"]
// ];
// t = transpose(arr, reverse=true);
// // Returns:
// // [
// // ["i", "f", "c"],
// // ["h", "e", "b"],
// // ["g", "d", "a"]
// // ]
// Example: Transpose on a list of numbers returns the list unchanged
// transpose([3,4,5]); // Returns: [3,4,5]
function transpose(arr, reverse=false) =
assert( is_list(arr) && len(arr)>0, "Input to transpose must be a nonempty list.")
is_list(arr[0])
? let( len0 = len(arr[0]) )
assert([for(a=arr) if(!is_list(a) || len(a)!=len0) 1 ]==[], "Input to transpose has inconsistent row lengths." )
reverse
? [for (i=[0:1:len0-1])
[ for (j=[0:1:len(arr)-1]) arr[len(arr)-1-j][len0-1-i] ] ]
: [for (i=[0:1:len0-1])
[ for (j=[0:1:len(arr)-1]) arr[j][i] ] ]
: assert( is_vector(arr), "Input to transpose must be a vector or list of lists.")
arr;
// Section: Matrices
// Function: is_matrix_symmetric()
// Usage:
// b = is_matrix_symmetric(A, [eps])
// Description:
// Returns true if the input matrix is symmetric, meaning it equals its transpose.
// Matrix should have numerical entries.
// Zips together two or more lists into a single list. For example, if you have two
// lists [3,4,5], and [8,7,6], and zip them together, you get [ [3,8],[4,7],[5,6] ].
// The list returned will be as long as the shortest list passed to zip().
// Arguments:
// A = matrix to test
// eps = epsilon for comparing equality. Default: 1e-12
function is_matrix_symmetric(A,eps=1e-12) =
approx(A,transpose(A), eps);
// a = The first list, or a list of lists if b and c are not given.
// b = The second list, if given.
// c = The third list, if given.
// Example:
// a = [9,8,7,6]; b = [1,2,3];
// for (p=zip(a,b)) echo(p);
// // ECHO: [9,1]
// // ECHO: [8,2]
// // ECHO: [7,3]
function zip(a,b,c) =
b!=undef? zip([a,b,if (c!=undef) c]) :
let(n = min_length(a))
[for (i=[0:1:n-1]) [for (x=a) x[i]]];
// Function&Module: echo_matrix()
// Function: zip_long()
// Usage:
// echo_matrix(M, [description=], [sig=], [eps=]);
// dummy = echo_matrix(M, [description=], [sig=], [eps=]),
// pairs = zip_long(a,b);
// triples = zip_long(a,b,c);
// quads = zip_long([LIST1,LIST2,LIST3,LIST4]);
// Topics: List Handling, Iteration
// See Also: zip()
// Description:
// Display a numerical matrix in a readable columnar format with `sig` significant
// digits. Values smaller than eps display as zero. If you give a description
// it is displayed at the top.
function echo_matrix(M,description,sig=4,eps=1e-9) =
let(
horiz_line = chr(8213),
matstr = matrix_strings(M,sig=sig,eps=eps),
separator = str_join(repeat(horiz_line,10)),
dummy=echo(str(separator," ",is_def(description) ? description : ""))
[for(row=matstr) echo(row)]
)
echo(separator);
// Zips together two or more lists into a single list. For example, if you have two
// lists [3,4,5], and [8,7,6], and zip them together, you get [ [3,8],[4,7],[5,6] ].
// The list returned will be as long as the longest list passed to zip_long(), with
// shorter lists padded by the value in `fill`.
// Arguments:
// a = The first list, or a list of lists if b and c are not given.
// b = The second list, if given.
// c = The third list, if given.
// fill = The value to pad shorter lists with. Default: undef
// Example:
// a = [9,8,7,6]; b = [1,2,3];
// for (p=zip_long(a,b,fill=88)) echo(p);
// // ECHO: [9,1]
// // ECHO: [8,2]
// // ECHO: [7,3]
// // ECHO: [6,88]]
function zip_long(a,b,c,fill) =
b!=undef? zip_long([a,b,if (c!=undef) c],fill=fill) :
let(n = max_length(a))
[for (i=[0:1:n-1]) [for (x=a) i<len(x)? x[i] : fill]];
module echo_matrix(M,description,sig=4,eps=1e-9)
{
dummy = echo_matrix(M,description,sig,eps);
}
// vim: expandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap

16
attachments.scad
View file

@ -1622,8 +1622,8 @@ function _find_anchor(anchor, geom) =
for (face = faces)
let(
verts = select(rpts, face),
ys = columns(verts,1),
zs = columns(verts,2)
ys = column(verts,1),
zs = column(verts,2)
)
if (max(ys) >= -eps && max(zs) >= -eps &&
min(ys) <= eps && min(zs) <= eps)
@ -1640,7 +1640,7 @@ function _find_anchor(anchor, geom) =
)
assert(len(hits)>0, "Anchor vector does not intersect with the shape. Attachment failed.")
let(
furthest = max_index(columns(hits,0)),
furthest = max_index(column(hits,0)),
dist = hits[furthest][0],
pos = hits[furthest][2],
hitnorms = [for (hit = hits) if (approx(hit[0],dist,eps=eps)) hit[1]],
@ -1665,7 +1665,7 @@ function _find_anchor(anchor, geom) =
) vnf==EMPTY_VNF? [anchor, [0,0,0], unit(anchor), 0] :
let(
rpts = apply(rot(from=anchor, to=RIGHT) * move(point3d(-cp)), vnf[0]),
maxx = max(columns(rpts,0)),
maxx = max(column(rpts,0)),
idxs = [for (i = idx(rpts)) if (approx(rpts[i].x, maxx)) i],
avep = sum(select(rpts,idxs))/len(idxs),
mpt = approx(point2d(anchor),[0,0])? [maxx,0,0] : avep,
@ -1716,7 +1716,7 @@ function _find_anchor(anchor, geom) =
)
if(!is_undef(isect) && !approx(isect,t[0])) [norm(isect), isect, n2]
],
maxidx = max_index(columns(isects,0)),
maxidx = max_index(column(isects,0)),
isect = isects[maxidx],
pos = point2d(cp) + isect[1],
vec = unit(isect[2],[0,1])
@ -1728,7 +1728,7 @@ function _find_anchor(anchor, geom) =
anchor = point2d(anchor),
m = rot(from=anchor, to=RIGHT) * move(-[cp.x, cp.y, 0]),
rpts = apply(m, flatten(rgn)),
maxx = max(columns(rpts,0)),
maxx = max(column(rpts,0)),
idxs = [for (i = idx(rpts)) if (approx(rpts[i].x, maxx)) i],
miny = min([for (i=idxs) rpts[i].y]),
maxy = max([for (i=idxs) rpts[i].y]),
@ -1757,7 +1757,7 @@ function _find_anchor(anchor, geom) =
if(!is_undef(isect) && !approx(isect,t[0]))
[norm(isect), isect, n2]
],
maxidx = max_index(columns(isects,0)),
maxidx = max_index(column(isects,0)),
isect = isects[maxidx],
pos = point3d(cp) + point3d(isect[1]) + unit([0,0,anchor.z],CENTER)*l/2,
xyvec = unit(isect[2],[0,1]),
@ -1775,7 +1775,7 @@ function _find_anchor(anchor, geom) =
rot(from=xyanch, to=RIGHT, planar=true)
) * move(-[cp.x, cp.y]),
rpts = apply(m, flatten(rgn)),
maxx = max(columns(rpts,0)),
maxx = max(column(rpts,0)),
idxs = [for (i = idx(rpts)) if (approx(rpts[i].x, maxx)) i],
ys = [for (i=idxs) rpts[i].y],
midy = (min(ys)+max(ys))/2,

20
beziers.scad
View file

@ -1043,9 +1043,9 @@ function bezier_patch_points(patch, u, v) =
assert(is_num(u) || !is_undef(u[0]))
assert(is_num(v) || !is_undef(v[0]))
let(
vbezes = [for (i = idx(patch[0])) bezier_points(columns(patch,i), is_num(u)? [u] : u)]
vbezes = [for (i = idx(patch[0])) bezier_points(column(patch,i), is_num(u)? [u] : u)]
)
[for (i = idx(vbezes[0])) bezier_points(columns(vbezes,i), is_num(v)? [v] : v)];
[for (i = idx(vbezes[0])) bezier_points(column(vbezes,i), is_num(v)? [v] : v)];
// Function: bezier_triangle_point()
@ -1358,7 +1358,7 @@ function bezier_patch_degenerate(patch, splinesteps=16, reverse=false, return_ed
all(row_degen) && all(col_degen) ? // fully degenerate case
[EMPTY_VNF, repeat([patch[0][0]],4)] :
all(row_degen) ? // degenerate to a line (top to bottom)
let(pts = bezier_points(columns(patch,0), samplepts))
let(pts = bezier_points(column(patch,0), samplepts))
[EMPTY_VNF, [pts,pts,[pts[0]],[last(pts)]]] :
all(col_degen) ? // degenerate to a line (left to right)
let(pts = bezier_points(patch[0], samplepts))
@ -1367,7 +1367,7 @@ function bezier_patch_degenerate(patch, splinesteps=16, reverse=false, return_ed
let(pts = bezier_patch_points(patch, samplepts, samplepts))
[
vnf_vertex_array(pts, reverse=!reverse),
[columns(pts,0), columns(pts,len(pts)-1), pts[0], last(pts)]
[column(pts,0), column(pts,len(pts)-1), pts[0], last(pts)]
] :
top_degen && bot_degen ?
let(
@ -1376,17 +1376,17 @@ function bezier_patch_degenerate(patch, splinesteps=16, reverse=false, return_ed
if (splinesteps%2==0) splinesteps+1,
each reverse(list([3:2:splinesteps]))
],
bpatch = [for(i=[0:1:len(patch[0])-1]) bezier_points(columns(patch,i), samplepts)],
bpatch = [for(i=[0:1:len(patch[0])-1]) bezier_points(column(patch,i), samplepts)],
pts = [
[bpatch[0][0]],
for(j=[0:splinesteps-2]) bezier_points(columns(bpatch,j+1), lerpn(0,1,rowcount[j])),
for(j=[0:splinesteps-2]) bezier_points(column(bpatch,j+1), lerpn(0,1,rowcount[j])),
[last(bpatch[0])]
],
vnf = vnf_tri_array(pts, reverse=!reverse)
) [
vnf,
[
columns(pts,0),
column(pts,0),
[for(row=pts) last(row)],
pts[0],
last(pts),
@ -1405,16 +1405,16 @@ function bezier_patch_degenerate(patch, splinesteps=16, reverse=false, return_ed
full_degen = len(patch)>=4 && all(select(row_degen,1,ceil(len(patch)/2-1))),
rowmax = full_degen ? count(splinesteps+1) :
[for(j=[0:splinesteps]) j<=splinesteps/2 ? 2*j : splinesteps],
bpatch = [for(i=[0:1:len(patch[0])-1]) bezier_points(columns(patch,i), samplepts)],
bpatch = [for(i=[0:1:len(patch[0])-1]) bezier_points(column(patch,i), samplepts)],
pts = [
[bpatch[0][0]],
for(j=[1:splinesteps]) bezier_points(columns(bpatch,j), lerpn(0,1,rowmax[j]+1))
for(j=[1:splinesteps]) bezier_points(column(bpatch,j), lerpn(0,1,rowmax[j]+1))
],
vnf = vnf_tri_array(pts, reverse=!reverse)
) [
vnf,
[
columns(pts,0),
column(pts,0),
[for(row=pts) last(row)],
pts[0],
last(pts),

2
bottlecaps.scad
View file

@ -1121,7 +1121,7 @@ module sp_neck(diam,type,wall,id,style="L",bead=false, anchor, spin, orient)
isect400 = [for(seg=pair(beadpts)) let(segisect = line_intersection([[T/2,0],[T/2,1]] , seg, LINE, SEGMENT)) if (is_def(segisect)) segisect.y];
extra_bot = type==400 && bead ? -min(columns(beadpts,1))+max(isect400) : 0;
extra_bot = type==400 && bead ? -min(column(beadpts,1))+max(isect400) : 0;
bead_shift = type==400 ? H+max(isect400) : entry[5]+W/2; // entry[5] is L
attachable(anchor,spin,orient,r=bead ? beadmax : T/2, l=H+extra_bot){

2
geometry.scad
View file

@ -1994,7 +1994,7 @@ function align_polygon(reference, poly, angles, cp, trans, return_ind=false) =
return_error=true
)
],
scores = columns(alignments,1),
scores = column(alignments,1),
minscore = min(scores),
minind = [for(i=idx(scores)) if (scores[i]<minscore+EPSILON) i],
dummy = is_def(angles) ? echo(best_angles = select(list(angles), minind)):0,

713
linalg.scad Normal file
View file

@ -0,0 +1,713 @@
//////////////////////////////////////////////////////////////////////
// LibFile: linalg.scad
// This file provides linear algebra, with support for matrix construction,
// solutions to linear systems of equations, QR and Cholesky factorizations, and
// matrix inverse.
// Includes:
// include <BOSL2/std.scad>
//////////////////////////////////////////////////////////////////////
// Section: Matrices
// The matrix, a rectangular array of numbers which represents a linear transformation,
// is the fundamental object in linear algebra. In OpenSCAD a matrix is a list of lists of numbers
// with a rectangular structure. Because OpenSCAD treats all data the same, most of the functions that
// index matrices or construct them will work on matrices (lists of lists) whose elements are not numbers but may be
// arbitrary data: strings, booleans, or even other lists. It may even be acceptable in some cases if the structure is non-rectangular.
// Of course, linear algebra computations and solutions require true matrices with rectangular structure, where all the entries are
// finite numbers.
// .
// Matrices in OpenSCAD are lists of row vectors. However, a potential source of confusion is that OpenSCAD
// treats vectors as either column vectors or row vectors as demanded by
// context. Thus both `v*M` and `M*v` are valid if `M` is square and `v` has the right length. If you want to multiply
// `M` on the left by `v` and `w` you can do this with `[v,w]*M` but if you want to multiply on the right side with `v` and `w` as
// column vectors, you now need to use {{transpose()}} because OpenSCAD doesn't adjust matrices
// contextually: `A=M*transpose([v,w])`. The solutions are now columns of A and you must extract
// them with {{column()}} or take the transpose of `A`.
// Section: Matrix testing and display
// Function: is_matrix()
// Usage:
// test = is_matrix(A, [m], [n], [square])
// Description:
// Returns true if A is a numeric matrix of height m and width n with finite entries. If m or n
// are omitted or set to undef then true is returned for any positive dimension.
// Arguments:
// A = The matrix to test.
// m = If given, requires the matrix to have this height.
// n = Is given, requires the matrix to have this width.
// square = If true, matrix must have height equal to width. Default: false
function is_matrix(A,m,n,square=false) =
is_list(A)
&& (( is_undef(m) && len(A) ) || len(A)==m)
&& (!square || len(A) == len(A[0]))
&& is_vector(A[0],n)
&& is_consistent(A);
// Function: is_matrix_symmetric()
// Usage:
// b = is_matrix_symmetric(A, [eps])
// Description:
// Returns true if the input matrix is symmetric, meaning it approximately equals its transpose.
// The matrix can have arbitrary entries.
// Arguments:
// A = matrix to test
// eps = epsilon for comparing equality. Default: 1e-12
function is_matrix_symmetric(A,eps=1e-12) =
approx(A,transpose(A), eps);
// Function&Module: echo_matrix()
// Usage:
// echo_matrix(M, [description=], [sig=], [eps=]);
// dummy = echo_matrix(M, [description=], [sig=], [eps=]),
// Description:
// Display a numerical matrix in a readable columnar format with `sig` significant
// digits. Values smaller than eps display as zero. If you give a description
// it is displayed at the top.
function echo_matrix(M,description,sig=4,eps=1e-9) =
let(
horiz_line = chr(8213),
matstr = matrix_strings(M,sig=sig,eps=eps),
separator = str_join(repeat(horiz_line,10)),
dummy=echo(str(separator," ",is_def(description) ? description : ""))
[for(row=matstr) echo(row)]
)
echo(separator);
module echo_matrix(M,description,sig=4,eps=1e-9)
{
dummy = echo_matrix(M,description,sig,eps);
}
// Section: Matrix indexing
// Function: column()
// Usage:
// list = column(M, i);
// Topics: Array Handling, List Handling
// See Also: select(), slice()
// Description:
// Extracts entry i from each list in M, or equivalently column i from the matrix M, and returns it as a vector.
// This function will return `undef` at all entry positions indexed by i not found in M.
// Arguments:
// M = The given list of lists.
// idx = The index, list of indices, or range of indices to fetch.
// Example:
// M = [[1,2,3,4],[5,6,7,8],[9,10,11,12],[13,14,15,16]];
// a = column(M,2); // Returns [3, 7, 11, 15]
// b = column(M,0); // Returns [1, 5, 9, 13]
// N = [ [1,2], [3], [4,5], [6,7,8] ];
// c = column(N,1); // Returns [1,undef,5,7]
// data = [[1,[3,4]], [3, [9,3]], [4, [3,1]]]; // Matrix with non-numeric entries
// d = column(data,0); // Returns [1,3,4]
// e = column(data,1); // Returns [[3,4],[9,3],[3,1]]
function column(M, i) =
assert( is_list(M), "The input is not a list." )
assert( is_int(i) && i>=0, "Invalid index")
[for(row=M) row[i]];
// Function: submatrix()
// Usage:
// mat = submatrix(M, idx1, idx2);
// Topics: Matrices, Array Handling
// See Also: column(), block_matrix(), submatrix_set()
// Description:
// The input must be a list of lists (a matrix or 2d array). Returns a submatrix by selecting the rows listed in idx1 and columns listed in idx2.
// Arguments:
// M = Given list of lists
// idx1 = rows index list or range
// idx2 = column index list or range
// Example:
// M = [[ 1, 2, 3, 4, 5],
// [ 6, 7, 8, 9,10],
// [11,12,13,14,15],
// [16,17,18,19,20],
// [21,22,23,24,25]];
// submatrix(M,[1:2],[3:4]); // Returns [[9, 10], [14, 15]]
// submatrix(M,[1], [3,4])); // Returns [[9,10]]
// submatrix(M,1, [3,4])); // Returns [[9,10]]
// submatrix(M,1,3)); // Returns [[9]]
// submatrix(M, [3,4],1); // Returns [[17],[22]]);
// submatrix(M, [1,3],[2,4]); // Returns [[8,10],[18,20]]);
// A = [[true, 17, "test"],
// [[4,2], 91, false],
// [6, [3,4], undef]];
// submatrix(A,[0,2],[1,2]); // Returns [[17, "test"], [[3, 4], undef]]
function submatrix(M,idx1,idx2) =
[for(i=idx1) [for(j=idx2) M[i][j] ] ];
// Section: Matrix construction and modification
// Function: ident()
// Usage:
// mat = ident(n);
// Topics: Affine, Matrices
// Description:
// Create an `n` by `n` square identity matrix.
// Arguments:
// n = The size of the identity matrix square, `n` by `n`.
// Example:
// mat = ident(3);
// // Returns:
// // [
// // [1, 0, 0],
// // [0, 1, 0],
// // [0, 0, 1]
// // ]
// Example:
// mat = ident(4);
// // Returns:
// // [
// // [1, 0, 0, 0],
// // [0, 1, 0, 0],
// // [0, 0, 1, 0],
// // [0, 0, 0, 1]
// // ]
function ident(n) = [
for (i = [0:1:n-1]) [
for (j = [0:1:n-1]) (i==j)? 1 : 0
]
];
// Function: diagonal_matrix()
// Usage:
// mat = diagonal_matrix(diag, [offdiag]);
// Topics: Matrices, Array Handling
// See Also: column(), submatrix()
// Description:
// Creates a square matrix with the items in the list `diag` on
// its diagonal. The off diagonal entries are set to offdiag,
// which is zero by default.
// Arguments:
// diag = A list of items to put in the diagnal cells of the matrix.
// offdiag = Value to put in non-diagonal matrix cells.
function diagonal_matrix(diag, offdiag=0) =
assert(is_list(diag) && len(diag)>0)
[for(i=[0:1:len(diag)-1]) [for(j=[0:len(diag)-1]) i==j?diag[i] : offdiag]];
// Function: transpose()
// Usage:
// M = transpose(M, [reverse]);
// Topics: Matrices, Array Handling
// See Also: submatrix(), block_matrix(), hstack(), flatten()
// Description:
// Returns the transpose of the given input matrix. The input can be a matrix with arbitrary entries or
// a numerical vector. If you give a vector then transpose returns it unchanged.
// When reverse=true, the transpose is done across to the secondary diagonal. (See example below.)
// By default, reverse=false.
// Example:
// M = [
// [1, 2, 3],
// [4, 5, 6],
// [7, 8, 9]
// ];
// t = transpose(M);
// // Returns:
// // [
// // [1, 4, 7],
// // [2, 5, 8],
// // [3, 6, 9]
// // ]
// Example:
// M = [
// [1, 2, 3],
// [4, 5, 6]
// ];
// t = transpose(M);
// // Returns:
// // [
// // [1, 4],
// // [2, 5],
// // [3, 6],
// // ]
// Example:
// M = [
// [1, 2, 3],
// [4, 5, 6],
// [7, 8, 9]
// ];
// t = transpose(M, reverse=true);
// // Returns:
// // [
// // [9, 6, 3],
// // [8, 5, 2],
// // [7, 4, 1]
// // ]
// Example: Transpose on a list of numbers returns the list unchanged
// transpose([3,4,5]); // Returns: [3,4,5]
// Example: Transpose on non-numeric input
// arr = [
// [ "a", "b", "c"],
// [ "d", "e", "f"],
// [[1,2],[3,4],[5,6]]
// ];
// t = transpose(arr);
// // Returns:
// // [
// // ["a", "d", [1,2]],
// // ["b", "e", [3,4]],
// // ["c", "f", [5,6]],
// // ]
function transpose(M, reverse=false) =
assert( is_list(M) && len(M)>0, "Input to transpose must be a nonempty list.")
is_list(M[0])
? let( len0 = len(M[0]) )
assert([for(a=M) if(!is_list(a) || len(a)!=len0) 1 ]==[], "Input to transpose has inconsistent row lengths." )
reverse
? [for (i=[0:1:len0-1])
[ for (j=[0:1:len(M)-1]) M[len(M)-1-j][len0-1-i] ] ]
: [for (i=[0:1:len0-1])
[ for (j=[0:1:len(M)-1]) M[j][i] ] ]
: assert( is_vector(M), "Input to transpose must be a vector or list of lists.")
M;
// Function: outer_product()
// Usage:
// x = outer_product(u,v);
// Description:
// Compute the outer product of two vectors, a matrix.
// Usage:
// M = outer_product(u,v);
function outer_product(u,v) =
assert(is_vector(u) && is_vector(v), "The inputs must be vectors.")
[for(ui=u) ui*v];
// Function: submatrix_set()
// Usage:
// mat = submatrix_set(M, A, [m], [n]);
// Topics: Matrices, Array Handling
// See Also: column(), submatrix()
// Description:
// Sets a submatrix of M equal to the matrix A. By default the top left corner of M is set to A, but
// you can specify offset coordinates m and n. If A (as adjusted by m and n) extends beyond the bounds
// of M then the extra entries are ignored. You can pass in A=[[]], a null matrix, and M will be
// returned unchanged. This function works on arbitrary lists of lists and the input M need not be rectangular in shape.
// Arguments:
// M = Original matrix.
// A = Submatrix of new values to write into M
// m = Row number of upper-left corner to place A at. Default: 0
// n = Column number of upper-left corner to place A at. Default: 0
function submatrix_set(M,A,m=0,n=0) =
assert(is_list(M))
assert(is_list(A))
assert(is_int(m))
assert(is_int(n))
let( badrows = [for(i=idx(A)) if (!is_list(A[i])) i])
assert(badrows==[], str("Input submatrix malformed rows: ",badrows))
[for(i=[0:1:len(M)-1])
assert(is_list(M[i]), str("Row ",i," of input matrix is not a list"))
[for(j=[0:1:len(M[i])-1])
i>=m && i <len(A)+m && j>=n && j<len(A[0])+n ? A[i-m][j-n] : M[i][j]]];
// Function: hstack()
// Usage:
// A = hstack(M1, M2)
// A = hstack(M1, M2, M3)
// A = hstack([M1, M2, M3, ...])
// Topics: Matrices, Array Handling
// See Also: column(), submatrix(), block_matrix()
// Description:
// Constructs a matrix by horizontally "stacking" together compatible matrices or vectors. Vectors are treated as columsn in the stack.
// This command is the inverse of `column`. Note: strings given in vectors are broken apart into lists of characters. Strings given
// in matrices are preserved as strings. If you need to combine vectors of strings use array_group as shown below to convert the
// vector into a column matrix. Also note that vertical stacking can be done directly with concat.
// Arguments:
// M1 = If given with other arguments, the first matrix (or vector) to stack. If given alone, a list of matrices/vectors to stack.
// M2 = Second matrix/vector to stack
// M3 = Third matrix/vector to stack.
// Example:
// M = ident(3);
// v1 = [2,3,4];
// v2 = [5,6,7];
// v3 = [8,9,10];
// a = hstack(v1,v2); // Returns [[2, 5], [3, 6], [4, 7]]
// b = hstack(v1,v2,v3); // Returns [[2, 5, 8],
// // [3, 6, 9],
// // [4, 7, 10]]
// c = hstack([M,v1,M]); // Returns [[1, 0, 0, 2, 1, 0, 0],
// // [0, 1, 0, 3, 0, 1, 0],
// // [0, 0, 1, 4, 0, 0, 1]]
// d = hstack(column(M,0), submatrix(M,idx(M),[1 2])); // Returns M
// strvec = ["one","two"];
// strmat = [["three","four"], ["five","six"]];
// e = hstack(strvec,strvec); // Returns [["o", "n", "e", "o", "n", "e"],
// // ["t", "w", "o", "t", "w", "o"]]
// f = hstack(array_group(strvec,1), array_group(strvec,1));
// // Returns [["one", "one"],
// // ["two", "two"]]
// g = hstack(strmat,strmat); // Returns: [["three", "four", "three", "four"],
// // [ "five", "six", "five", "six"]]
function hstack(M1, M2, M3) =
(M3!=undef)? hstack([M1,M2,M3]) :
(M2!=undef)? hstack([M1,M2]) :
assert(all([for(v=M1) is_list(v)]), "One of the inputs to hstack is not a list")
let(
minlen = min_length(M1),
maxlen = max_length(M1)
)
assert(minlen==maxlen, "Input vectors to hstack must have the same length")
[for(row=[0:1:minlen-1])
[for(matrix=M1)
each matrix[row]
]
];
// Function: block_matrix()
// Usage:
// bmat = block_matrix([[M11, M12,...],[M21, M22,...], ... ]);
// Topics: Matrices, Array Handling
// See Also: column(), submatrix()
// Description:
// Create a block matrix by supplying a matrix of matrices, which will
// be combined into one unified matrix. Every matrix in one row
// must have the same height, and the combined width of the matrices
// in each row must be equal. Strings will stay strings.
// Example:
// A = [[1,2],
// [3,4]];
// B = ident(2);
// C = block_matrix([[A,B],[B,A],[A,B]]);
// // Returns:
// // [[1, 2, 1, 0],
// // [3, 4, 0, 1],
// // [1, 0, 1, 2],
// // [0, 1, 3, 4],
// // [1, 2, 1, 0],
// // [3, 4, 0, 1]]);
// D = block_matrix([[A,B], ident(4)]);
// // Returns:
// // [[1, 2, 1, 0],
// // [3, 4, 0, 1],
// // [1, 0, 0, 0],
// // [0, 1, 0, 0],
// // [0, 0, 1, 0],
// // [0, 0, 0, 1]]);
// E = [["one", "two"], [3,4]];
// F = block_matrix([[E,E]]);
// // Returns:
// // [["one", "two", "one", "two"],
// // [ 3, 4, 3, 4]]
function block_matrix(M) =
let(
bigM = [for(bigrow = M) each hstack(bigrow)],
len0 = len(bigM[0]),
badrows = [for(row=bigM) if (len(row)!=len0) 1]
)
assert(badrows==[], "Inconsistent or invalid input")
bigM;
// Section: Solving Linear Equations and Matrix Factorizations
// Function: linear_solve()
// Usage:
// solv = linear_solve(A,b)
// Description:
// Solves the linear system Ax=b. If `A` is square and non-singular the unique solution is returned. If `A` is overdetermined
// the least squares solution is returned. If `A` is underdetermined, the minimal norm solution is returned.
// If `A` is rank deficient or singular then linear_solve returns `[]`. If `b` is a matrix that is compatible with `A`
// then the problem is solved for the matrix valued right hand side and a matrix is returned. Note that if you
// want to solve Ax=b1 and Ax=b2 that you need to form the matrix `transpose([b1,b2])` for the right hand side and then
// transpose the returned value.
function linear_solve(A,b,pivot=true) =
assert(is_matrix(A), "Input should be a matrix.")
let(
m = len(A),
n = len(A[0])
)
assert(is_vector(b,m) || is_matrix(b,m),"Invalid right hand side or incompatible with the matrix")
let (
qr = m<n? qr_factor(transpose(A),pivot) : qr_factor(A,pivot),
maxdim = max(n,m),
mindim = min(n,m),
Q = submatrix(qr[0],[0:maxdim-1], [0:mindim-1]),
R = submatrix(qr[1],[0:mindim-1], [0:mindim-1]),
P = qr[2],
zeros = [for(i=[0:mindim-1]) if (approx(R[i][i],0)) i]
)
zeros != [] ? [] :
m<n ? Q*back_substitute(R,transpose(P)*b,transpose=true) // Too messy to avoid input checks here
: P*_back_substitute(R, transpose(Q)*b); // Calling internal version skips input checks
// Function: matrix_inverse()
// Usage:
// mat = matrix_inverse(A)
// Description:
// Compute the matrix inverse of the square matrix `A`. If `A` is singular, returns `undef`.
// Note that if you just want to solve a linear system of equations you should NOT use this function.
// Instead use {{linear_solve()}}, or use {{qr_factor()}}. The computation
// will be faster and more accurate.
function matrix_inverse(A) =
assert(is_matrix(A) && len(A)==len(A[0]),"Input to matrix_inverse() must be a square matrix")
linear_solve(A,ident(len(A)));
// Function: rot_inverse()
// Usage:
// B = rot_inverse(A)
// Description:
// Inverts a 2d (3x3) or 3d (4x4) rotation matrix. The matrix can be a rotation around any center,
// so it may include a translation.
function rot_inverse(T) =
assert(is_matrix(T,square=true),"Matrix must be square")
let( n = len(T))
assert(n==3 || n==4, "Matrix must be 3x3 or 4x4")
let(
rotpart = [for(i=[0:n-2]) [for(j=[0:n-2]) T[j][i]]],
transpart = [for(row=[0:n-2]) T[row][n-1]]
)
assert(approx(determinant(T),1),"Matrix is not a rotation")
concat(hstack(rotpart, -rotpart*transpart),[[for(i=[2:n]) 0, 1]]);
// Function: null_space()
// Usage:
// x = null_space(A)
// Description:
// Returns an orthonormal basis for the null space of `A`, namely the vectors {x} such that Ax=0.
// If the null space is just the origin then returns an empty list.
function null_space(A,eps=1e-12) =
assert(is_matrix(A))
let(
Q_R = qr_factor(transpose(A),pivot=true),
R = Q_R[1],
zrows = [for(i=idx(R)) if (all_zero(R[i],eps)) i]
)
len(zrows)==0 ? [] :
select(transpose(Q_R[0]), zrows);
// Function: qr_factor()
// Usage:
// qr = qr_factor(A,[pivot]);
// Description:
// Calculates the QR factorization of the input matrix A and returns it as the list [Q,R,P]. This factorization can be
// used to solve linear systems of equations. The factorization is `A = Q*R*transpose(P)`. If pivot is false (the default)
// then P is the identity matrix and A = Q*R. If pivot is true then column pivoting results in an R matrix where the diagonal
// is non-decreasing. The use of pivoting is supposed to increase accuracy for poorly conditioned problems, and is necessary
// for rank estimation or computation of the null space, but it may be slower.
function qr_factor(A, pivot=false) =
assert(is_matrix(A), "Input must be a matrix." )
let(
m = len(A),
n = len(A[0])
)
let(
qr = _qr_factor(A, Q=ident(m),P=ident(n), pivot=pivot, col=0, m = m, n = n),
Rzero = let( R = qr[1]) [
for(i=[0:m-1]) [
let( ri = R[i] )
for(j=[0:n-1]) i>j ? 0 : ri[j]
]
]
) [qr[0], Rzero, qr[2]];
function _qr_factor(A,Q,P, pivot, col, m, n) =
col >= min(m-1,n) ? [Q,A,P] :
let(
swap = !pivot ? 1
: _swap_matrix(n,col,col+max_index([for(i=[col:n-1]) sqr([for(j=[col:m-1]) A[j][i]])])),
A = pivot ? A*swap : A,
x = [for(i=[col:1:m-1]) A[i][col]],
alpha = (x[0]<=0 ? 1 : -1) * norm(x),
u = x - concat([alpha],repeat(0,m-1)),
v = alpha==0 ? u : u / norm(u),
Qc = ident(len(x)) - 2*outer_product(v,v),
Qf = [for(i=[0:m-1]) [for(j=[0:m-1]) i<col || j<col ? (i==j ? 1 : 0) : Qc[i-col][j-col]]]
)
_qr_factor(Qf*A, Q*Qf, P*swap, pivot, col+1, m, n);
// Produces an n x n matrix that swaps column i and j (when multiplied on the right)
function _swap_matrix(n,i,j) =
assert(i<n && j<n && i>=0 && j>=0, "Swap indices out of bounds")
[for(y=[0:n-1]) [for (x=[0:n-1])
x==i ? (y==j ? 1 : 0)
: x==j ? (y==i ? 1 : 0)
: x==y ? 1 : 0]];
// Function: back_substitute()
// Usage:
// x = back_substitute(R, b, [transpose]);
// Description:
// Solves the problem Rx=b where R is an upper triangular square matrix. The lower triangular entries of R are
// ignored. If transpose==true then instead solve transpose(R)*x=b.
// You can supply a compatible matrix b and it will produce the solution for every column of b. Note that if you want to
// solve Rx=b1 and Rx=b2 you must set b to transpose([b1,b2]) and then take the transpose of the result. If the matrix
// is singular (e.g. has a zero on the diagonal) then it returns [].
function back_substitute(R, b, transpose = false) =
assert(is_matrix(R, square=true))
let(n=len(R))
assert(is_vector(b,n) || is_matrix(b,n),str("R and b are not compatible in back_substitute ",n, len(b)))
transpose
? reverse(_back_substitute(transpose(R, reverse=true), reverse(b)))
: _back_substitute(R,b);
function _back_substitute(R, b, x=[]) =
let(n=len(R))
len(x) == n ? x
: let(ind = n - len(x) - 1)
R[ind][ind] == 0 ? []
: let(
newvalue = len(x)==0
? b[ind]/R[ind][ind]
: (b[ind]-list_tail(R[ind],ind+1) * x)/R[ind][ind]
)
_back_substitute(R, b, concat([newvalue],x));
// Function: cholesky()
// Usage:
// L = cholesky(A);
// Description:
// Compute the cholesky factor, L, of the symmetric positive definite matrix A.
// The matrix L is lower triangular and `L * transpose(L) = A`. If the A is
// not symmetric then an error is displayed. If the matrix is symmetric but
// not positive definite then undef is returned.
function cholesky(A) =
assert(is_matrix(A,square=true),"A must be a square matrix")
assert(is_matrix_symmetric(A),"Cholesky factorization requires a symmetric matrix")
_cholesky(A,ident(len(A)), len(A));
function _cholesky(A,L,n) =
A[0][0]<0 ? undef : // Matrix not positive definite
len(A) == 1 ? submatrix_set(L,[[sqrt(A[0][0])]], n-1,n-1):
let(
i = n+1-len(A)
)
let(
sqrtAii = sqrt(A[0][0]),
Lnext = [for(j=[0:n-1])
[for(k=[0:n-1])
j<i-1 || k<i-1 ? (j==k ? 1 : 0)
: j==i-1 && k==i-1 ? sqrtAii
: j==i-1 ? 0
: k==i-1 ? A[j-(i-1)][0]/sqrtAii
: j==k ? 1 : 0]],
Anext = submatrix(A,[1:n-1], [1:n-1]) - outer_product(list_tail(A[0]), list_tail(A[0]))/A[0][0]
)
_cholesky(Anext,L*Lnext,n);
// Section: Matrix Properties: Determinants, Norm, Trace
// Function: det2()
// Usage:
// d = det2(M);
// Description:
// Rturns the determinant for the given 2x2 matrix.
// Arguments:
// M = The 2x2 matrix to get the determinant of.
// Example:
// M = [ [6,-2], [1,8] ];
// det = det2(M); // Returns: 50
function det2(M) =
assert(is_def(M) && M*0==[[0,0],[0,0]], "Expected square matrix (2x2)")
cross(M[0],M[1]);
// Function: det3()
// Usage:
// d = det3(M);
// Description:
// Returns the determinant for the given 3x3 matrix.
// Arguments:
// M = The 3x3 square matrix to get the determinant of.
// Example:
// M = [ [6,4,-2], [1,-2,8], [1,5,7] ];
// det = det3(M); // Returns: -334
function det3(M) =
assert(is_def(M) && M*0==[[0,0,0],[0,0,0],[0,0,0]], "Expected square matrix (3x3).")
M[0][0] * (M[1][1]*M[2][2]-M[2][1]*M[1][2]) -
M[1][0] * (M[0][1]*M[2][2]-M[2][1]*M[0][2]) +
M[2][0] * (M[0][1]*M[1][2]-M[1][1]*M[0][2]);
// Function: det4()
// Usage:
// d = det4(M);
// Description:
// Returns the determinant for the given 4x4 matrix.
// Arguments:
// M = The 4x4 square matrix to get the determinant of.
// Example:
// M = [ [6,4,-2,1], [1,-2,8,-3], [1,5,7,4], [2,3,4,7] ];
// det = det4(M); // Returns: -1773
function det4(M) =
assert(is_def(M) && M*0==[[0,0,0,0],[0,0,0,0],[0,0,0,0],[0,0,0,0]], "Expected square matrix (4x4).")
M[0][0]*M[1][1]*M[2][2]*M[3][3] + M[0][0]*M[1][2]*M[2][3]*M[3][1] + M[0][0]*M[1][3]*M[2][1]*M[3][2]
+ M[0][1]*M[1][0]*M[2][3]*M[3][2] + M[0][1]*M[1][2]*M[2][0]*M[3][3] + M[0][1]*M[1][3]*M[2][2]*M[3][0]
+ M[0][2]*M[1][0]*M[2][1]*M[3][3] + M[0][2]*M[1][1]*M[2][3]*M[3][0] + M[0][2]*M[1][3]*M[2][0]*M[3][1]
+ M[0][3]*M[1][0]*M[2][2]*M[3][1] + M[0][3]*M[1][1]*M[2][0]*M[3][2] + M[0][3]*M[1][2]*M[2][1]*M[3][0]
- M[0][0]*M[1][1]*M[2][3]*M[3][2] - M[0][0]*M[1][2]*M[2][1]*M[3][3] - M[0][0]*M[1][3]*M[2][2]*M[3][1]
- M[0][1]*M[1][0]*M[2][2]*M[3][3] - M[0][1]*M[1][2]*M[2][3]*M[3][0] - M[0][1]*M[1][3]*M[2][0]*M[3][2]
- M[0][2]*M[1][0]*M[2][3]*M[3][1] - M[0][2]*M[1][1]*M[2][0]*M[3][3] - M[0][2]*M[1][3]*M[2][1]*M[3][0]
- M[0][3]*M[1][0]*M[2][1]*M[3][2] - M[0][3]*M[1][1]*M[2][2]*M[3][0] - M[0][3]*M[1][2]*M[2][0]*M[3][1];
// Function: determinant()
// Usage:
// d = determinant(M);
// Description:
// Returns the determinant for the given square matrix.
// Arguments:
// M = The NxN square matrix to get the determinant of.
// Example:
// M = [ [6,4,-2,9], [1,-2,8,3], [1,5,7,6], [4,2,5,1] ];
// det = determinant(M); // Returns: 2267
function determinant(M) =
assert(is_list(M), "Input must be a square matrix." )
len(M)==1? M[0][0] :
len(M)==2? det2(M) :
len(M)==3? det3(M) :
len(M)==4? det4(M) :
assert(is_matrix(M, square=true), "Input must be a square matrix." )
sum(
[for (col=[0:1:len(M)-1])
((col%2==0)? 1 : -1) *
M[col][0] *
determinant(
[for (r=[1:1:len(M)-1])
[for (c=[0:1:len(M)-1])
if (c!=col) M[c][r]
]
]
)
]
);
// Function: norm_fro()
// Usage:
// norm_fro(A)
// Description:
// Computes frobenius norm of input matrix. The frobenius norm is the square root of the sum of the
// squares of all of the entries of the matrix. On vectors it is the same as the usual 2-norm.
// This is an easily computed norm that is convenient for comparing two matrices.
function norm_fro(A) =
assert(is_matrix(A) || is_vector(A))
norm(flatten(A));
// Function: matrix_trace()
// Usage:
// matrix_trace(M)
// Description:
// Computes the trace of a square matrix, the sum of the entries on the diagonal.
function matrix_trace(M) =
assert(is_matrix(M,square=true), "Input to trace must be a square matrix")
[for(i=[0:1:len(M)-1])1] * [for(i=[0:1:len(M)-1]) M[i][i]];

344
math.scad
View file

@ -5,7 +5,6 @@
// include <BOSL2/std.scad>
//////////////////////////////////////////////////////////////////////
// Section: Math Constants
// Constant: PHI
@ -669,7 +668,7 @@ function lcm(a,b=[]) =
function sum(v, dflt=0) =
v==[]? dflt :
assert(is_consistent(v), "Input to sum is non-numeric or inconsistent")
is_vector(v) || is_matrix(v) ? [for(i=v) 1]*v :
is_finite(v[0]) || is_vector(v[0]) ? [for(i=v) 1]*v :
_sum(v,v[0]*0);
function _sum(v,_total,_i=0) = _i>=len(v) ? _total : _sum(v,_total+v[_i], _i+1);
@ -799,18 +798,6 @@ function _cumprod_vec(v,_i=0,_acc=[]) =
);
// Function: outer_product()
// Usage:
// x = outer_product(u,v);
// Description:
// Compute the outer product of two vectors, a matrix.
// Usage:
// M = outer_product(u,v);
function outer_product(u,v) =
assert(is_vector(u) && is_vector(v), "The inputs must be vectors.")
[for(ui=u) ui*v];
// Function: mean()
// Usage:
// x = mean(v);
@ -877,335 +864,6 @@ function convolve(p,q) =
// Section: Matrix math
// Function: ident()
// Usage:
// mat = ident(n);
// Topics: Affine, Matrices
// Description:
// Create an `n` by `n` square identity matrix.
// Arguments:
// n = The size of the identity matrix square, `n` by `n`.
// Example:
// mat = ident(3);
// // Returns:
// // [
// // [1, 0, 0],
// // [0, 1, 0],
// // [0, 0, 1]
// // ]
// Example:
// mat = ident(4);
// // Returns:
// // [
// // [1, 0, 0, 0],
// // [0, 1, 0, 0],
// // [0, 0, 1, 0],
// // [0, 0, 0, 1]
// // ]
function ident(n) = [
for (i = [0:1:n-1]) [
for (j = [0:1:n-1]) (i==j)? 1 : 0
]
];
// Function: linear_solve()
// Usage:
// solv = linear_solve(A,b)
// Description:
// Solves the linear system Ax=b. If `A` is square and non-singular the unique solution is returned. If `A` is overdetermined
// the least squares solution is returned. If `A` is underdetermined, the minimal norm solution is returned.
// If `A` is rank deficient or singular then linear_solve returns `[]`. If `b` is a matrix that is compatible with `A`
// then the problem is solved for the matrix valued right hand side and a matrix is returned. Note that if you
// want to solve Ax=b1 and Ax=b2 that you need to form the matrix `transpose([b1,b2])` for the right hand side and then
// transpose the returned value.
function linear_solve(A,b,pivot=true) =
assert(is_matrix(A), "Input should be a matrix.")
let(
m = len(A),
n = len(A[0])
)
assert(is_vector(b,m) || is_matrix(b,m),"Invalid right hand side or incompatible with the matrix")
let (
qr = m<n? qr_factor(transpose(A),pivot) : qr_factor(A,pivot),
maxdim = max(n,m),
mindim = min(n,m),
Q = submatrix(qr[0],[0:maxdim-1], [0:mindim-1]),
R = submatrix(qr[1],[0:mindim-1], [0:mindim-1]),
P = qr[2],
zeros = [for(i=[0:mindim-1]) if (approx(R[i][i],0)) i]
)
zeros != [] ? [] :
m<n ? Q*back_substitute(R,transpose(P)*b,transpose=true) // Too messy to avoid input checks here
: P*_back_substitute(R, transpose(Q)*b); // Calling internal version skips input checks
// Function: matrix_inverse()
// Usage:
// mat = matrix_inverse(A)
// Description:
// Compute the matrix inverse of the square matrix `A`. If `A` is singular, returns `undef`.
// Note that if you just want to solve a linear system of equations you should NOT use this function.
// Instead use {{linear_solve()}}, or use {{qr_factor()}}. The computation
// will be faster and more accurate.
function matrix_inverse(A) =
assert(is_matrix(A) && len(A)==len(A[0]),"Input to matrix_inverse() must be a square matrix")
linear_solve(A,ident(len(A)));
// Function: rot_inverse()
// Usage:
// B = rot_inverse(A)
// Description:
// Inverts a 2d (3x3) or 3d (4x4) rotation matrix. The matrix can be a rotation around any center,
// so it may include a translation.
function rot_inverse(T) =
assert(is_matrix(T,square=true),"Matrix must be square")
let( n = len(T))
assert(n==3 || n==4, "Matrix must be 3x3 or 4x4")
let(
rotpart = [for(i=[0:n-2]) [for(j=[0:n-2]) T[j][i]]],
transpart = [for(row=[0:n-2]) T[row][n-1]]
)
assert(approx(determinant(T),1),"Matrix is not a rotation")
concat(hstack(rotpart, -rotpart*transpart),[[for(i=[2:n]) 0, 1]]);
// Function: null_space()
// Usage:
// x = null_space(A)
// Description:
// Returns an orthonormal basis for the null space of `A`, namely the vectors {x} such that Ax=0.
// If the null space is just the origin then returns an empty list.
function null_space(A,eps=1e-12) =
assert(is_matrix(A))
let(
Q_R = qr_factor(transpose(A),pivot=true),
R = Q_R[1],
zrow = [for(i=idx(R)) if (all_zero(R[i],eps)) i]
)
len(zrow)==0 ? [] :
transpose(columns(Q_R[0],zrow));
// Function: qr_factor()
// Usage:
// qr = qr_factor(A,[pivot]);
// Description:
// Calculates the QR factorization of the input matrix A and returns it as the list [Q,R,P]. This factorization can be
// used to solve linear systems of equations. The factorization is `A = Q*R*transpose(P)`. If pivot is false (the default)
// then P is the identity matrix and A = Q*R. If pivot is true then column pivoting results in an R matrix where the diagonal
// is non-decreasing. The use of pivoting is supposed to increase accuracy for poorly conditioned problems, and is necessary
// for rank estimation or computation of the null space, but it may be slower.
function qr_factor(A, pivot=false) =
assert(is_matrix(A), "Input must be a matrix." )
let(
m = len(A),
n = len(A[0])
)
let(
qr = _qr_factor(A, Q=ident(m),P=ident(n), pivot=pivot, column=0, m = m, n=n),
Rzero = let( R = qr[1]) [
for(i=[0:m-1]) [
let( ri = R[i] )
for(j=[0:n-1]) i>j ? 0 : ri[j]
]
]
) [qr[0], Rzero, qr[2]];
function _qr_factor(A,Q,P, pivot, column, m, n) =
column >= min(m-1,n) ? [Q,A,P] :
let(
swap = !pivot ? 1
: _swap_matrix(n,column,column+max_index([for(i=[column:n-1]) sqr([for(j=[column:m-1]) A[j][i]])])),
A = pivot ? A*swap : A,
x = [for(i=[column:1:m-1]) A[i][column]],
alpha = (x[0]<=0 ? 1 : -1) * norm(x),
u = x - concat([alpha],repeat(0,m-1)),
v = alpha==0 ? u : u / norm(u),
Qc = ident(len(x)) - 2*outer_product(v,v),
Qf = [for(i=[0:m-1]) [for(j=[0:m-1]) i<column || j<column ? (i==j ? 1 : 0) : Qc[i-column][j-column]]]
)
_qr_factor(Qf*A, Q*Qf, P*swap, pivot, column+1, m, n);
// Produces an n x n matrix that swaps column i and j (when multiplied on the right)
function _swap_matrix(n,i,j) =
assert(i<n && j<n && i>=0 && j>=0, "Swap indices out of bounds")
[for(y=[0:n-1]) [for (x=[0:n-1])
x==i ? (y==j ? 1 : 0)
: x==j ? (y==i ? 1 : 0)
: x==y ? 1 : 0]];
// Function: back_substitute()
// Usage:
// x = back_substitute(R, b, [transpose]);
// Description:
// Solves the problem Rx=b where R is an upper triangular square matrix. The lower triangular entries of R are
// ignored. If transpose==true then instead solve transpose(R)*x=b.
// You can supply a compatible matrix b and it will produce the solution for every column of b. Note that if you want to
// solve Rx=b1 and Rx=b2 you must set b to transpose([b1,b2]) and then take the transpose of the result. If the matrix
// is singular (e.g. has a zero on the diagonal) then it returns [].
function back_substitute(R, b, transpose = false) =
assert(is_matrix(R, square=true))
let(n=len(R))
assert(is_vector(b,n) || is_matrix(b,n),str("R and b are not compatible in back_substitute ",n, len(b)))
transpose
? reverse(_back_substitute(transpose(R, reverse=true), reverse(b)))
: _back_substitute(R,b);
function _back_substitute(R, b, x=[]) =
let(n=len(R))
len(x) == n ? x
: let(ind = n - len(x) - 1)
R[ind][ind] == 0 ? []
: let(
newvalue = len(x)==0
? b[ind]/R[ind][ind]
: (b[ind]-list_tail(R[ind],ind+1) * x)/R[ind][ind]
)
_back_substitute(R, b, concat([newvalue],x));
// Function: cholesky()
// Usage:
// L = cholesky(A);
// Description:
// Compute the cholesky factor, L, of the symmetric positive definite matrix A.
// The matrix L is lower triangular and `L * transpose(L) = A`. If the A is
// not symmetric then an error is displayed. If the matrix is symmetric but
// not positive definite then undef is returned.
function cholesky(A) =
assert(is_matrix(A,square=true),"A must be a square matrix")
assert(is_matrix_symmetric(A),"Cholesky factorization requires a symmetric matrix")
_cholesky(A,ident(len(A)), len(A));
function _cholesky(A,L,n) =
A[0][0]<0 ? undef : // Matrix not positive definite
len(A) == 1 ? submatrix_set(L,[[sqrt(A[0][0])]], n-1,n-1):
let(
i = n+1-len(A)
)
let(
sqrtAii = sqrt(A[0][0]),
Lnext = [for(j=[0:n-1])
[for(k=[0:n-1])
j<i-1 || k<i-1 ? (j==k ? 1 : 0)
: j==i-1 && k==i-1 ? sqrtAii
: j==i-1 ? 0
: k==i-1 ? A[j-(i-1)][0]/sqrtAii
: j==k ? 1 : 0]],
Anext = submatrix(A,[1:n-1], [1:n-1]) - outer_product(list_tail(A[0]), list_tail(A[0]))/A[0][0]
)
_cholesky(Anext,L*Lnext,n);
// Function: det2()
// Usage:
// d = det2(M);
// Description:
// Optimized function that returns the determinant for the given 2x2 square matrix.
// Arguments:
// M = The 2x2 square matrix to get the determinant of.
// Example:
// M = [ [6,-2], [1,8] ];
// det = det2(M); // Returns: 50
function det2(M) =
assert(is_matrix(M,2,2), "Matrix must be 2x2.")
M[0][0] * M[1][1] - M[0][1]*M[1][0];
// Function: det3()
// Usage:
// d = det3(M);
// Description:
// Optimized function that returns the determinant for the given 3x3 square matrix.
// Arguments:
// M = The 3x3 square matrix to get the determinant of.
// Example:
// M = [ [6,4,-2], [1,-2,8], [1,5,7] ];
// det = det3(M); // Returns: -334
function det3(M) =
assert(is_matrix(M,3,3), "Matrix must be 3x3.")
M[0][0] * (M[1][1]*M[2][2]-M[2][1]*M[1][2]) -
M[1][0] * (M[0][1]*M[2][2]-M[2][1]*M[0][2]) +
M[2][0] * (M[0][1]*M[1][2]-M[1][1]*M[0][2]);
// Function: determinant()
// Usage:
// d = determinant(M);
// Description:
// Returns the determinant for the given square matrix.
// Arguments:
// M = The NxN square matrix to get the determinant of.
// Example:
// M = [ [6,4,-2,9], [1,-2,8,3], [1,5,7,6], [4,2,5,1] ];
// det = determinant(M); // Returns: 2267
function determinant(M) =
assert(is_matrix(M, square=true), "Input should be a square matrix." )
len(M)==1? M[0][0] :
len(M)==2? det2(M) :
len(M)==3? det3(M) :
sum(
[for (col=[0:1:len(M)-1])
((col%2==0)? 1 : -1) *
M[col][0] *
determinant(
[for (r=[1:1:len(M)-1])
[for (c=[0:1:len(M)-1])
if (c!=col) M[c][r]
]
]
)
]
);
// Function: is_matrix()
// Usage:
// test = is_matrix(A, [m], [n], [square])
// Description:
// Returns true if A is a numeric matrix of height m and width n. If m or n
// are omitted or set to undef then true is returned for any positive dimension.
// Arguments:
// A = The matrix to test.
// m = Is given, requires the matrix to have the given height.
// n = Is given, requires the matrix to have the given width.
// square = If true, requires the matrix to have a width equal to its height. Default: false
function is_matrix(A,m,n,square=false) =
is_list(A)
&& (( is_undef(m) && len(A) ) || len(A)==m)
&& (!square || len(A) == len(A[0]))
&& is_vector(A[0],n)
&& is_consistent(A);
// Function: norm_fro()
// Usage:
// norm_fro(A)
// Description:
// Computes frobenius norm of input matrix. The frobenius norm is the square root of the sum of the
// squares of all of the entries of the matrix. On vectors it is the same as the usual 2-norm.
// This is an easily computed norm that is convenient for comparing two matrices.
function norm_fro(A) =
assert(is_matrix(A) || is_vector(A))
norm(flatten(A));
// Function: matrix_trace()
// Usage:
// matrix_trace(M)
// Description:
// Computes the trace of a square matrix, the sum of the entries on the diagonal.
function matrix_trace(M) =
assert(is_matrix(M,square=true), "Input to trace must be a square matrix")
[for(i=[0:1:len(M)-1])1] * [for(i=[0:1:len(M)-1]) M[i][i]];
// Section: Comparisons and Logic

4
paths.scad
View file

@ -427,7 +427,7 @@ function resample_path(path, N, spacing, closed=false) =
distlist = lerpn(0,length,N,false),
cuts = _path_cut_points(path, distlist, closed=closed)
)
[ each columns(cuts,0),
[ each column(cuts,0),
if (!closed) last(path) // Then add last point here
];
@ -1174,7 +1174,7 @@ function _assemble_path_fragments(fragments, eps=EPSILON, _finished=[]) =
len(fragments)==0? _finished :
let(
minxidx = min_index([
for (frag=fragments) min(columns(frag,0))
for (frag=fragments) min(column(frag,0))
]),
result_l = _assemble_a_path_from_fragments(
fragments=fragments,

2
regions.scad
View file

@ -267,7 +267,7 @@ function _region_region_intersections(region1, region2, closed1=true,closed2=tru
cornerpts = [for(i=[0:1])
[for(k=vector_search(points[i],eps,points[i]))
each if (len(k)>1) select(ptind[i],k)]],
risect = [for(i=[0:1]) concat(columns(intersections,i), cornerpts[i])],
risect = [for(i=[0:1]) concat(column(intersections,i), cornerpts[i])],
counts = [count(len(region1)), count(len(region2))],
pathind = [for(i=[0:1]) search(counts[i], risect[i], 0)]
)

18
rounding.scad
View file

@ -1270,7 +1270,7 @@ module convex_offset_extrude(
// The entry r[i] is [radius,z] for a given layer
r = move([0,bottom_height],p=concat(
reverse(offsets_bot), [[0,0], [0,middle]], move([0,middle], p=offsets_top)));
delta = [for(val=deltas(columns(r,0))) sign(val)];
delta = [for(val=deltas(column(r,0))) sign(val)];
below=[-thickness,0];
above=[0,thickness];
// layers is a list of pairs of the relative positions for each layer, e.g. [0,thickness]
@ -1937,8 +1937,8 @@ function rounded_prism(bottom, top, joint_bot=0, joint_top=0, joint_sides=0, k_b
verify_vert =
[for(i=[0:N-1],j=[0:4])
let(
vline = concat(select(columns(top_patch[i],j),2,4),
select(columns(bot_patch[i],j),2,4))
vline = concat(select(column(top_patch[i],j),2,4),
select(column(bot_patch[i],j),2,4))
)
if (!is_collinear(vline)) [i,j]],
//verify horiz edges
@ -1955,8 +1955,8 @@ function rounded_prism(bottom, top, joint_bot=0, joint_top=0, joint_sides=0, k_b
"Roundovers interfere with each other on bottom face: either input is self intersecting or top joint length is too large")
assert(debug || (verify_vert==[] && verify_horiz==[]), "Curvature continuity failed")
let(
vnf = vnf_merge([ each columns(top_samples,0),
each columns(bot_samples,0),
vnf = vnf_merge([ each column(top_samples,0),
each column(bot_samples,0),
for(pts=edge_points) vnf_vertex_array(pts),
debug ? vnf_from_polygons(faces)
: vnf_triangulate(vnf_from_polygons(faces))
@ -2114,7 +2114,7 @@ function _circle_mask(r) =
// ]),
// radius = [0,0,each repeat(slotradius,4),0,0], closed=false
// )
// ) apply(left(max(columns(slot,0))/2)*fwd(min(columns(slot,1))), slot);
// ) apply(left(max(column(slot,0))/2)*fwd(min(column(slot,1))), slot);
// stroke(slot(15,29,7));
// Example: A cylindrical container with rounded edges and a rounded finger slot.
// function slot(slotwidth, slotheight, slotradius) = let(
@ -2138,7 +2138,7 @@ function _circle_mask(r) =
// ]),
// radius = [0,0,each repeat(slotradius,4),0,0], closed=false
// )
// ) apply(left(max(columns(slot,0))/2)*fwd(min(columns(slot,1))), slot);
// ) apply(left(max(column(slot,0))/2)*fwd(min(column(slot,1))), slot);
// diam = 80;
// wall = 4;
// height = 40;
@ -2162,12 +2162,12 @@ module bent_cutout_mask(r, thickness, path, radius, convexity=10)
path = clockwise_polygon(path);
curvepoints = arc(d=thickness, angle = [-180,0]);
profiles = [for(pt=curvepoints) _cyl_hole(r+pt.x,apply(xscale((r+pt.x)/r), offset(path,delta=thickness/2+pt.y,check_valid=false,closed=true)))];
pathx = columns(path,0);
pathx = column(path,0);
minangle = (min(pathx)-thickness/2)*360/(2*PI*r);
maxangle = (max(pathx)+thickness/2)*360/(2*PI*r);
mindist = (r+thickness/2)/cos((maxangle-minangle)/2);
assert(maxangle-minangle<180,"Cutout angle span is too large. Must be smaller than 180.");
zmean = mean(columns(path,1));
zmean = mean(column(path,1));
innerzero = repeat([0,0,zmean], len(path));
outerpt = repeat( [1.5*mindist*cos((maxangle+minangle)/2),1.5*mindist*sin((maxangle+minangle)/2),zmean], len(path));
vnf_polyhedron(vnf_vertex_array([innerzero, each profiles, outerpt],col_wrap=true),convexity=convexity);

6
shapes2d.scad
View file

@ -363,7 +363,7 @@ function regular_ngon(n=6, r, d, or, od, ir, id, side, rounding=0, realign=false
)
each arc(N=steps, cp=p, r=rounding, start=a+180/n, angle=-360/n)
],
maxx_idx = max_index(columns(path2,0)),
maxx_idx = max_index(column(path2,0)),
path3 = polygon_shift(path2,maxx_idx)
) path3
),
@ -1005,7 +1005,7 @@ function teardrop2d(r, ang=45, cap_h, d, anchor=CENTER, spin=0) =
[-cap_w,cap_h]
], closed=true
),
maxx_idx = max_index(columns(path,0)),
maxx_idx = max_index(column(path,0)),
path2 = polygon_shift(path,maxx_idx)
) reorient(anchor,spin, two_d=true, path=path2, p=path2);
@ -1060,7 +1060,7 @@ function glued_circles(r, spread=10, tangent=30, d, anchor=CENTER, spin=0) =
[for (i=[0:1:lobesegs]) let(a=sa1+i*lobestep+180) r * [cos(a),sin(a)] + cp1],
tangent==0? [] : [for (i=[0:1:arcsegs]) let(a=ea2-i*arcstep) r2 * [cos(a),sin(a)] + cp2]
),
maxx_idx = max_index(columns(path,0)),
maxx_idx = max_index(column(path,0)),
path2 = reverse_polygon(polygon_shift(path,maxx_idx))
) reorient(anchor,spin, two_d=true, path=path2, extent=true, p=path2);

2
shapes3d.scad
View file

@ -2266,7 +2266,7 @@ module path_text(path, text, font, size, thickness, lettersize, offset=0, revers
usernorm = is_def(normal);
usetop = is_def(top);
normpts = is_undef(normal) ? (reverse?1:-1)*columns(pts,3) : _cut_interp(pts,path, normal);
normpts = is_undef(normal) ? (reverse?1:-1)*column(pts,3) : _cut_interp(pts,path, normal);
toppts = is_undef(top) ? undef : _cut_interp(pts,path,top);
for(i=idx(text))
let( tangent = pts[i][2] )

2
skin.scad
View file

@ -985,7 +985,7 @@ function path_sweep2d(shape, path, closed=false, caps, quality=1, style="min_edg
[for(pt = profile)
let(
ofs = offset(path, delta=-flip*pt.x, return_faces=true,closed=closed, quality=quality),
map = columns(_ofs_vmap(ofs,closed=closed),1)
map = column(_ofs_vmap(ofs,closed=closed),1)
)
select(path3d(ofs[0],pt.y),map)
]

1
std.scad
View file

@ -22,6 +22,7 @@ include <paths.scad>
include <edges.scad>
include <arrays.scad>
include <math.scad>
include <linalg.scad>
include <trigonometry.scad>
include <vectors.scad>
include <quaternions.scad>

100
tests/test_arrays.scad
View file

@ -451,36 +451,6 @@ module test_add_scalar() {
test_add_scalar();
module test_columns() {
v = [[1,2,3,4],[5,6,7,8],[9,10,11,12],[13,14,15,16]];
assert(columns(v,2) == [3, 7, 11, 15]);
assert(columns(v,[2]) == [[3], [7], [11], [15]]);
assert(columns(v,[2,1]) == [[3, 2], [7, 6], [11, 10], [15, 14]]);
assert(columns(v,[1:3]) == [[2, 3, 4], [6, 7, 8], [10, 11, 12], [14, 15, 16]]);
}
test_columns();
// Need decision about behavior for out of bounds ranges, empty ranges
module test_submatrix(){
M = [[1,2,3,4,5],
[6,7,8,9,10],
[11,12,13,14,15],
[16,17,18,19,20],
[21,22,23,24,25]];
assert_equal(submatrix(M,[1:2], [3:4]), [[9,10],[14,15]]);
assert_equal(submatrix(M,[1], [3,4]), [[9,10]]);
assert_equal(submatrix(M,1, [3,4]), [[9,10]]);
assert_equal(submatrix(M, [3,4],1), [[17],[22]]);
assert_equal(submatrix(M, [1,3],[2,4]), [[8,10],[18,20]]);
assert_equal(submatrix(M, 1,3), [[9]]);
A = [[true, 17, "test"],
[[4,2], 91, false],
[6, [3,4], undef]];
assert_equal(submatrix(A,[0,2],[1,2]),[[17, "test"], [[3, 4], undef]]);
}
test_submatrix();
module test_force_list() {
assert_equal(force_list([3,4,5]), [3,4,5]);
@ -535,67 +505,8 @@ module test_zip() {
v3 = [8,9,10,11];
assert(zip(v1,v3) == [[1,8],[2,9],[3,10],[4,11]]);
assert(zip([v1,v3]) == [[1,8],[2,9],[3,10],[4,11]]);
assert(zip([v1,v2],fit="short") == [[1,5],[2,6],[3,7]]);
assert(zip([v1,v2],fit="long") == [[1,5],[2,6],[3,7],[4,undef]]);
assert(zip([v1,v2],fit="long", fill=0) == [[1,5],[2,6],[3,7],[4,0]]);
assert(zip([v1,v2,v3],fit="long") == [[1,5,8],[2,6,9],[3,7,10],[4,undef,11]]);
}
//test_zip();
module test_hstack() {
M = ident(3);
v1 = [2,3,4];
v2 = [5,6,7];
v3 = [8,9,10];
a = hstack(v1,v2);
b = hstack(v1,v2,v3);
c = hstack([M,v1,M]);
d = hstack(columns(M,0), columns(M,[1, 2]));
assert_equal(a,[[2, 5], [3, 6], [4, 7]]);
assert_equal(b,[[2, 5, 8], [3, 6, 9], [4, 7, 10]]);
assert_equal(c,[[1, 0, 0, 2, 1, 0, 0], [0, 1, 0, 3, 0, 1, 0], [0, 0, 1, 4, 0, 0, 1]]);
assert_equal(d,M);
strmat = [["three","four"], ["five","six"]];
assert_equal(hstack(strmat,strmat), [["three", "four", "three", "four"], ["five", "six", "five", "six"]]);
strvec = ["one","two"];
assert_equal(hstack(strvec,strmat),[["o", "n", "e", "three", "four"], ["t", "w", "o", "five", "six"]]);
}
test_hstack();
module test_block_matrix() {
A = [[1,2],[3,4]];
B = ident(2);
assert_equal(block_matrix([[A,B],[B,A],[A,B]]), [[1,2,1,0],[3,4,0,1],[1,0,1,2],[0,1,3,4],[1,2,1,0],[3,4,0,1]]);
assert_equal(block_matrix([[A,B],ident(4)]), [[1,2,1,0],[3,4,0,1],[1,0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,1]]);
text = [["aa","bb"],["cc","dd"]];
assert_equal(block_matrix([[text,B]]), [["aa","bb",1,0],["cc","dd",0,1]]);
}
test_block_matrix();
module test_diagonal_matrix() {
assert_equal(diagonal_matrix([1,2,3]), [[1,0,0],[0,2,0],[0,0,3]]);
assert_equal(diagonal_matrix([1,"c",2]), [[1,0,0],[0,"c",0],[0,0,2]]);
assert_equal(diagonal_matrix([1,"c",2],"X"), [[1,"X","X"],["X","c","X"],["X","X",2]]);
assert_equal(diagonal_matrix([[1,1],[2,2],[3,3]], [0,0]), [[ [1,1],[0,0],[0,0]], [[0,0],[2,2],[0,0]], [[0,0],[0,0],[3,3]]]);
}
test_diagonal_matrix();
module test_submatrix_set() {
test = [[1,2,3,4,5],[6,7,8,9,10],[11,12,13,14,15], [16,17,18,19,20]];
ragged = [[1,2,3,4,5],[6,7,8,9,10],[11,12], [16,17]];
assert_equal(submatrix_set(test,[[9,8],[7,6]]), [[9,8,3,4,5],[7,6,8,9,10],[11,12,13,14,15], [16,17,18,19,20]]);
assert_equal(submatrix_set(test,[[9,7],[8,6]],1),[[1,2,3,4,5],[9,7,8,9,10],[8,6,13,14,15], [16,17,18,19,20]]);
assert_equal(submatrix_set(test,[[9,8],[7,6]],n=1), [[1,9,8,4,5],[6,7,6,9,10],[11,12,13,14,15], [16,17,18,19,20]]);
assert_equal(submatrix_set(test,[[9,8],[7,6]],1,2), [[1,2,3,4,5],[6,7,9,8,10],[11,12,7,6,15], [16,17,18,19,20]]);
assert_equal(submatrix_set(test,[[9,8],[7,6]],-1,-1), [[6,2,3,4,5],[6,7,8,9,10],[11,12,13,14,15], [16,17,18,19,20]]);
assert_equal(submatrix_set(test,[[9,8],[7,6]],n=4), [[1,2,3,4,9],[6,7,8,9,7],[11,12,13,14,15], [16,17,18,19,20]]);
assert_equal(submatrix_set(test,[[9,8],[7,6]],7,7), [[1,2,3,4,5],[6,7,8,9,10],[11,12,13,14,15], [16,17,18,19,20]]);
assert_equal(submatrix_set(ragged, [["a","b"],["c","d"]], 1, 1), [[1,2,3,4,5],[6,"a","b",9,10],[11,"c"], [16,17]]);
assert_equal(submatrix_set(test, [[]]), test);
}
test_submatrix_set();
test_zip();
module test_array_group() {
@ -645,13 +556,4 @@ module test_array_dim() {
test_array_dim();
module test_transpose() {
assert(transpose([[1,2,3],[4,5,6],[7,8,9]]) == [[1,4,7],[2,5,8],[3,6,9]]);
assert(transpose([[1,2,3],[4,5,6]]) == [[1,4],[2,5],[3,6]]);
assert(transpose([[1,2,3],[4,5,6]],reverse=true) == [[6,3], [5,2], [4,1]]);
assert(transpose([3,4,5]) == [3,4,5]);
}
test_transpose();
// vim: expandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap

155
tests/test_math.scad
View file

@ -595,13 +595,13 @@ module test_mean() {
}
test_mean();
/*
module test_median() {
assert_equal(median([2,3,7]), 4.5);
assert_equal(median([[1,2,3], [3,4,5], [8,9,10]]), [4.5,5.5,6.5]);
assert_equal(median([2,3,7]), 3);
assert_equal(median([2,4,5,8]), 4.5);
}
test_median();
*/
module test_convolve() {
@ -618,40 +618,6 @@ test_convolve();
module test_matrix_inverse() {
assert_approx(matrix_inverse(rot([20,30,40])), [[0.663413948169,0.556670399226,-0.5,0],[-0.47302145844,0.829769465589,0.296198132726,0],[0.579769465589,0.0400087565481,0.813797681349,0],[0,0,0,1]]);
}
test_matrix_inverse();
module test_det2() {
assert_equal(det2([[6,-2], [1,8]]), 50);
assert_equal(det2([[4,7], [3,2]]), -13);
assert_equal(det2([[4,3], [3,4]]), 7);
}
test_det2();
module test_det3() {
M = [ [6,4,-2], [1,-2,8], [1,5,7] ];
assert_equal(det3(M), -334);
}
test_det3();
module test_determinant() {
M = [ [6,4,-2,9], [1,-2,8,3], [1,5,7,6], [4,2,5,1] ];
assert_equal(determinant(M), 2267);
}
test_determinant();
module test_matrix_trace() {
M = [ [6,4,-2,9], [1,-2,8,3], [1,5,7,6], [4,2,5,1] ];
assert_equal(matrix_trace(M), 6-2+7+1);
}
test_matrix_trace();
// Logic
@ -975,40 +941,6 @@ test_back_substitute();
module test_norm_fro(){
assert_approx(norm_fro([[2,3,4],[4,5,6]]), 10.29563014098700);
} test_norm_fro();
module test_linear_solve(){
M = [[-2,-5,-1,3],
[3,7,6,2],
[6,5,-1,-6],
[-7,1,2,3]];
assert_approx(linear_solve(M, [-3,43,-11,13]), [1,2,3,4]);
assert_approx(linear_solve(M, [[-5,8],[18,-61],[4,7],[-1,-12]]), [[1,-2],[1,-3],[1,-4],[1,-5]]);
assert_approx(linear_solve([[2]],[4]), [2]);
assert_approx(linear_solve([[2]],[[4,8]]), [[2, 4]]);
assert_approx(linear_solve(select(M,0,2), [2,4,4]), [ 2.254871220604705e+00,
-8.378819388897780e-01,
2.330507118860985e-01,
8.511278195488737e-01]);
assert_approx(linear_solve(columns(M,[0:2]), [2,4,4,4]),
[-2.457142857142859e-01,
5.200000000000000e-01,
7.428571428571396e-02]);
assert_approx(linear_solve([[1,2,3,4]], [2]), [0.066666666666666, 0.13333333333, 0.2, 0.266666666666]);
assert_approx(linear_solve([[1],[2],[3],[4]], [4,3,2,1]), [2/3]);
rd = [[-2,-5,-1,3],
[3,7,6,2],
[3,7,6,2],
[-7,1,2,3]];
assert_equal(linear_solve(rd,[1,2,3,4]),[]);
assert_equal(linear_solve(select(rd,0,2), [2,4,4]), []);
assert_equal(linear_solve(transpose(select(rd,0,2)), [2,4,3,4]), []);
}
test_linear_solve();
module test_cumprod(){
@ -1186,85 +1118,6 @@ module test_quadratic_roots(){
test_quadratic_roots();
module test_null_space(){
assert_equal(null_space([[3,2,1],[3,6,3],[3,9,-3]]),[]);
function nullcheck(A,dim) =
let(v=null_space(A))
len(v)==dim && all_zero(A*transpose(v),eps=1e-12);
A = [[-1, 2, -5, 2],[-3,-1,3,-3],[5,0,5,0],[3,-4,11,-4]];
assert(nullcheck(A,1));
B = [
[ 4, 1, 8, 6, -2, 3],
[ 10, 5, 10, 10, 0, 5],
[ 8, 1, 8, 8, -6, 1],
[ -8, -8, 6, -1, -8, -1],
[ 2, 2, 0, 1, 2, 1],
[ 2, -3, 10, 6, -8, 1],
];
assert(nullcheck(B,3));
}
test_null_space();
module test_qr_factor() {
// Check that R is upper triangular
function is_ut(R) =
let(bad = [for(i=[1:1:len(R)-1], j=[0:min(i-1, len(R[0])-1)]) if (!approx(R[i][j],0)) 1])
bad == [];
// Test the R is upper trianglar, Q is orthogonal and qr=M
function qrok(qr,M) =
is_ut(qr[1]) && approx(qr[0]*transpose(qr[0]), ident(len(qr[0]))) && approx(qr[0]*qr[1],M) && qr[2]==ident(len(qr[2]));
// Test the R is upper trianglar, Q is orthogonal, R diagonal non-increasing and qrp=M
function qrokpiv(qr,M) =
is_ut(qr[1])
&& approx(qr[0]*transpose(qr[0]), ident(len(qr[0])))
&& approx(qr[0]*qr[1]*transpose(qr[2]),M)
&& is_decreasing([for(i=[0:1:min(len(qr[1]),len(qr[1][0]))-1]) abs(qr[1][i][i])]);
M = [[1,2,9,4,5],
[6,7,8,19,10],
[11,12,13,14,15],
[1,17,18,19,20],
[21,22,10,24,25]];
assert(qrok(qr_factor(M),M));
assert(qrok(qr_factor(select(M,0,3)),select(M,0,3)));
assert(qrok(qr_factor(transpose(select(M,0,3))),transpose(select(M,0,3))));
A = [[1,2,9,4,5],
[6,7,8,19,10],
[0,0,0,0,0],
[1,17,18,19,20],
[21,22,10,24,25]];
assert(qrok(qr_factor(A),A));
B = [[1,2,0,4,5],
[6,7,0,19,10],
[0,0,0,0,0],
[1,17,0,19,20],
[21,22,0,24,25]];
assert(qrok(qr_factor(B),B));
assert(qrok(qr_factor([[7]]), [[7]]));
assert(qrok(qr_factor([[1,2,3]]), [[1,2,3]]));
assert(qrok(qr_factor([[1],[2],[3]]), [[1],[2],[3]]));
assert(qrokpiv(qr_factor(M,pivot=true),M));
assert(qrokpiv(qr_factor(select(M,0,3),pivot=true),select(M,0,3)));
assert(qrokpiv(qr_factor(transpose(select(M,0,3)),pivot=true),transpose(select(M,0,3))));
assert(qrokpiv(qr_factor(B,pivot=true),B));
assert(qrokpiv(qr_factor([[7]],pivot=true), [[7]]));
assert(qrokpiv(qr_factor([[1,2,3]],pivot=true), [[1,2,3]]));
assert(qrokpiv(qr_factor([[1],[2],[3]],pivot=true), [[1],[2],[3]]));
}
test_qr_factor();
module test_poly_mult(){
assert_equal(poly_mult([3,2,1],[4,5,6,7]),[12,23,32,38,20,7]);

6
threading.scad
View file

@ -926,8 +926,8 @@ module generic_threaded_rod(
assert(higang1 < twist/2);
assert(higang2 < twist/2);
prof3d = path3d(profile);
pdepth = -min(columns(profile,1));
pmax = pitch * max(columns(profile,1));
pdepth = -min(column(profile,1));
pmax = pitch * max(column(profile,1));
rmax = max(_r1,_r2)+pmax;
depth = pdepth * pitch;
dummy1 = assert(_r1>depth && _r2>depth, "Screw profile deeper than rod radius");
@ -1087,7 +1087,7 @@ module generic_threaded_nut(
bevel1 = first_defined([bevel1,bevel,false]);
bevel2 = first_defined([bevel2,bevel,false]);
dummy1 = assert(is_num(pitch) && pitch>0);
depth = -pitch*min(columns(profile,1));
depth = -pitch*min(column(profile,1));
attachable(anchor,spin,orient, size=[od/cos(30),od,h]) {
difference() {
cyl(d=od/cos(30), h=h, center=true, $fn=6,chamfer1=bevel1?depth:undef,chamfer2=bevel2?depth:undef);

4
vectors.scad
View file

@ -543,8 +543,8 @@ function vector_nearest(query, k, target) =
"More results are requested than the number of points.")
tgpts
? let( tree = _bt_tree(target, count(len(target))) )
columns(_bt_nearest( query, k, target, tree),0)
: columns(_bt_nearest( query, k, target[0], target[1]),0);
column(_bt_nearest( query, k, target, tree),0)
: column(_bt_nearest( query, k, target[0], target[1]),0);
//Ball tree nearest

6
vnf.scad
View file

@ -366,7 +366,7 @@ function vnf_from_polygons(polygons) =
function _path_path_closest_vertices(path1,path2) =
let(
dists = [for (i=idx(path1)) let(j=closest_point(path1[i],path2)) [j,norm(path2[j]-path1[i])]],
i1 = min_index(columns(dists,1)),
i1 = min_index(column(dists,1)),
i2 = dists[i1][0]
) [dists[i1][1], i1, i2];
@ -396,7 +396,7 @@ function _cleave_connected_region(region) =
for (i=[1:1:len(region)-1])
_path_path_closest_vertices(region[0],region[i])
],
idxi = min_index(columns(dists,0)),
idxi = min_index(column(dists,0)),
newoline = _join_paths_at_vertices(
region[0], region[idxi+1],
dists[idxi][1], dists[idxi][2]
@ -584,7 +584,7 @@ function vnf_slice(vnf,dir,cuts) =
function _split_polygon_at_x(poly, x) =
let(
xs = columns(poly,0)
xs = column(poly,0)
) (min(xs) >= x || max(xs) <= x)? [poly] :
let(
poly2 = [