openscad/BOSL2
1
0
Fork 0
mirror of https://github.com/revarbat/BOSL2.git synced 2025-09-25 10:39:10 +02:00
Code Issues Projects Releases Wiki Activity

Merge branch 'master' of https://github.com/revarbat/BOSL2

Browse Source
This commit is contained in:
Adrian Mariano
2020-07-24 22:40:20 -04:00
parent d2855718f3 65125ed57a
commit f822e4634e
40 changed files with 1329 additions and 245 deletions
Download Patch File Download Diff File Expand all files Collapse all files

6
affine.scad
View File

@@ -277,9 +277,9 @@ function affine_frame_map(x,y,z, reverse=false) =
assert(yvalid,"Input y must be a length 3 vector")
assert(zvalid,"Input z must be a length 3 vector")
let(
x = is_undef(x)? undef : unit(x),
y = is_undef(y)? undef : unit(y),
z = is_undef(z)? undef : unit(z),
x = is_undef(x)? undef : unit(x,RIGHT),
y = is_undef(y)? undef : unit(y,BACK),
z = is_undef(z)? undef : unit(z,UP),
map = is_undef(x)? [cross(y,z), y, z] :
is_undef(y)? [x, cross(z,x), z] :
is_undef(z)? [x, y, cross(x,y)] :

6
arrays.scad
View File

@@ -91,7 +91,11 @@ function slice(arr,st,end) = let(
// in_list("bar", ["foo", "bar", "baz"]); // Returns true.
// in_list("bee", ["foo", "bar", "baz"]); // Returns false.
// in_list("bar", [[2,"foo"], [4,"bar"], [3,"baz"]], idx=1); // Returns true.
function in_list(x,l,idx=undef) = search([x], l, num_returns_per_match=1, index_col_num=idx) != [[]];
function in_list(val,list,idx=undef) =
let( s = search([val], list, num_returns_per_match=1, index_col_num=idx)[0] )
s==[] ? false
: is_undef(idx) ? val==list[s]
: val==list[s][idx];
// Function: min_index()

33
attachments.scad
View File

@@ -420,12 +420,12 @@ function find_anchor(anchor, geom) =
bot = point3d(vmul(point2d(size)/2,axy),-h/2),
top = point3d(vmul(point2d(size2)/2,axy)+shift,h/2),
pos = point3d(cp) + lerp(bot,top,u) + offset,
sidevec = unit(rot(from=UP, to=top-bot, p=point3d(axy))),
vvec = anchor==CENTER? UP : unit([0,0,anchor.z]),
sidevec = unit(rot(from=UP, to=top-bot, p=point3d(axy)),UP),
vvec = anchor==CENTER? UP : unit([0,0,anchor.z],UP),
vec = anchor==CENTER? UP :
approx(axy,[0,0])? unit(anchor) :
approx(axy,[0,0])? unit(anchor,UP) :
approx(anchor.z,0)? sidevec :
unit((sidevec+vvec)/2)
unit((sidevec+vvec)/2,UP)
) [anchor, pos, vec, oang]
) : type == "cyl"? ( //r1, r2, l, shift
let(
@@ -433,24 +433,24 @@ function find_anchor(anchor, geom) =
r1 = is_num(rr1)? [rr1,rr1] : rr1,
r2 = is_num(rr2)? [rr2,rr2] : rr2,
u = (anchor.z+1)/2,
axy = unit(point2d(anchor)),
axy = unit(point2d(anchor),[0,0]),
bot = point3d(vmul(r1,axy), -l/2),
top = point3d(vmul(r2,axy)+shift, l/2),
pos = point3d(cp) + lerp(bot,top,u) + offset,
sidevec = rot(from=UP, to=top-bot, p=point3d(axy)),
vvec = anchor==CENTER? UP : unit([0,0,anchor.z]),
vvec = anchor==CENTER? UP : unit([0,0,anchor.z],UP),
vec = anchor==CENTER? UP :
approx(axy,[0,0])? unit(anchor) :
approx(axy,[0,0])? unit(anchor,UP) :
approx(anchor.z,0)? sidevec :
unit((sidevec+vvec)/2)
unit((sidevec+vvec)/2,UP)
) [anchor, pos, vec, oang]
) : type == "spheroid"? ( //r
let(
rr = geom[1],
r = is_num(rr)? [rr,rr,rr] : rr,
anchor = unit(point3d(anchor)),
anchor = unit(point3d(anchor),CENTER),
pos = point3d(cp) + vmul(r,anchor) + offset,
vec = unit(vmul(r,anchor))
vec = unit(vmul(r,anchor),UP)
) [anchor, pos, vec, oang]
) : type == "vnf_isect"? ( //vnf
let(
@@ -490,7 +490,8 @@ function find_anchor(anchor, geom) =
faceplane = plane_from_points(faceverts),
nrm = plane_normal(faceplane)
) nrm
]) / len(nfaces)
]) / len(nfaces),
UP
)
)
[anchor, pos, n, oang]
@@ -513,15 +514,15 @@ function find_anchor(anchor, geom) =
frpt = [size.x/2*anchor.x, -size.y/2],
bkpt = [size2/2*anchor.x, size.y/2],
pos = point2d(cp) + lerp(frpt, bkpt, u) + offset,
vec = unit(rot(from=BACK, to=bkpt-frpt, p=anchor))
vec = unit(rot(from=BACK, to=bkpt-frpt, p=anchor),[0,1])
) [anchor, pos, vec, 0]
) : type == "circle"? ( //r
let(
rr = geom[1],
r = is_num(rr)? [rr,rr] : rr,
pos = point2d(cp) + vmul(r,anchor) + offset,
anchor = unit(point2d(anchor)),
vec = unit(vmul([r.y,r.x],anchor))
anchor = unit(point2d(anchor),[0,0]),
vec = unit(vmul([r.y,r.x],anchor),[0,1])
) [anchor, pos, vec, 0]
) : type == "path_isect"? ( //path
let(
@@ -534,7 +535,7 @@ function find_anchor(anchor, geom) =
isect = ray_segment_intersection([[0,0],anchor], seg1),
n = is_undef(isect)? [0,1] :
!approx(isect, t[1])? line_normal(seg1) :
unit((line_normal(seg1)+line_normal(seg2))/2),
unit((line_normal(seg1)+line_normal(seg2))/2,[0,1]),
n2 = vector_angle(anchor,n)>90? -n : n
)
if(!is_undef(isect) && !approx(isect,t[0])) [norm(isect), isect, n2]
@@ -542,7 +543,7 @@ function find_anchor(anchor, geom) =
maxidx = max_index(subindex(isects,0)),
isect = isects[maxidx],
pos = point2d(cp) + isect[1],
vec = unit(isect[2])
vec = unit(isect[2],[0,1])
) [anchor, pos, vec, 0]
) : type == "path_extent"? ( //path
let(

1
coords.scad
View File

@@ -28,7 +28,6 @@ function point2d(p, fill=0) = [for (i=[0:1]) (p[i]==undef)? fill : p[i]];
// every vector has the same length.
// Arguments:
// points = A list of 2D or 3D points/vectors.
// fill = Value to fill missing values in vectors with.
function path2d(points) =
assert(is_path(points,dim=undef,fast=true),"Input to path2d is not a path")
let (result = points * concat(ident(2), repeat([0,0], len(points[0])-2)))

17
debug.scad
View File

@@ -191,16 +191,17 @@ module debug_faces(vertices, faces, size=1, disabled=false) {
if (face[0] < 0 || face[1] < 0 || face[2] < 0 || face[0] >= vlen || face[1] >= vlen || face[2] >= vlen) {
echo("BAD FACE: ", vlen=vlen, face=face);
} else {
v0 = vertices[face[0]];
v1 = vertices[face[1]];
v2 = vertices[face[2]];
c = mean(select(vertices,face));
verts = select(vertices,face);
c = mean(verts);
v0 = verts[0];
v1 = verts[1];
v2 = verts[2];
dv0 = unit(v1 - v0);
dv1 = unit(v2 - v0);
nrm0 = unit(cross(dv0, dv1));
nrm1 = [0, 0, 1];
axis = unit(cross(nrm0, nrm1));
ang = vector_angle(nrm0, nrm1);
nrm0 = cross(dv0, dv1);
nrm1 = UP;
axis = vector_axis(nrm0, nrm1);
ang = vector_angle(nrm0, nrm1);
theta = atan2(nrm0[1], nrm0[0]);
translate(c) {
rotate(a=180-ang, v=axis) {

202
geometry.scad
View File

@@ -133,8 +133,9 @@ function distance_from_line(line, pt) =
// color("green") stroke([p1,p1+10*n], endcap2="arrow2");
// color("blue") move_copies([p1,p2]) circle(d=2, $fn=12);
function line_normal(p1,p2) =
is_undef(p2)? line_normal(p1[0],p1[1]) :
unit([p1.y-p2.y,p2.x-p1.x]);
is_undef(p2)?
assert(is_path(p1,2)) line_normal(p1[0],p1[1]) :
assert(is_vector(p1,2)&&is_vector(p2,2)) unit([p1.y-p2.y,p2.x-p1.x]);
// 2D Line intersection from two segments.
@@ -252,34 +253,192 @@ function segment_intersection(s1,s2,eps=EPSILON) =
// Usage:
// line_closest_point(line,pt);
// Description:
// Returns the point on the given `line` that is closest to the given point `pt`.
// Returns the point on the given 2D or 3D `line` that is closest to the given point `pt`.
// The `line` and `pt` args should either both be 2D or both 3D.
// Arguments:
// line = A list of two points that are on the unbounded line.
// pt = The point to find the closest point on the line to.
// Example(2D):
// line = [[-30,0],[30,30]];
// pt = [-32,-10];
// p2 = line_closest_point(line,pt);
// stroke(line, endcaps="arrow2");
// color("blue") translate(pt) circle(r=1,$fn=12);
// color("red") translate(p2) circle(r=1,$fn=12);
// Example(2D):
// line = [[-30,0],[30,30]];
// pt = [-5,0];
// p2 = line_closest_point(line,pt);
// stroke(line, endcaps="arrow2");
// color("blue") translate(pt) circle(r=1,$fn=12);
// color("red") translate(p2) circle(r=1,$fn=12);
// Example(2D):
// line = [[-30,0],[30,30]];
// pt = [40,25];
// p2 = line_closest_point(line,pt);
// stroke(line, endcaps="arrow2");
// color("blue") translate(pt) circle(r=1,$fn=12);
// color("red") translate(p2) circle(r=1,$fn=12);
// Example(FlatSpin):
// line = [[-30,-15,0],[30,15,30]];
// pt = [5,5,5];
// p2 = line_closest_point(line,pt);
// stroke(line, endcaps="arrow2");
// color("blue") translate(pt) sphere(r=1,$fn=12);
// color("red") translate(p2) sphere(r=1,$fn=12);
// Example(FlatSpin):
// line = [[-30,-15,0],[30,15,30]];
// pt = [-35,-15,0];
// p2 = line_closest_point(line,pt);
// stroke(line, endcaps="arrow2");
// color("blue") translate(pt) sphere(r=1,$fn=12);
// color("red") translate(p2) sphere(r=1,$fn=12);
// Example(FlatSpin):
// line = [[-30,-15,0],[30,15,30]];
// pt = [40,15,25];
// p2 = line_closest_point(line,pt);
// stroke(line, endcaps="arrow2");
// color("blue") translate(pt) sphere(r=1,$fn=12);
// color("red") translate(p2) sphere(r=1,$fn=12);
function line_closest_point(line,pt) =
assert(is_path(line)&&len(line)==2)
assert(same_shape(pt,line[0]))
assert(!approx(line[0],line[1]))
let(
n = line_normal(line),
isect = _general_line_intersection(line,[pt,pt+n])
) isect[0];
seglen = norm(line[1]-line[0]),
segvec = (line[1]-line[0])/seglen,
projection = (pt-line[0]) * segvec
)
line[0] + projection*segvec;
// Function: ray_closest_point()
// Usage:
// ray_closest_point(seg,pt);
// Description:
// Returns the point on the given 2D or 3D ray `ray` that is closest to the given point `pt`.
// The `ray` and `pt` args should either both be 2D or both 3D.
// Arguments:
// ray = The ray, given as a list `[START,POINT]` of the start-point START, and a point POINT on the ray.
// pt = The point to find the closest point on the ray to.
// Example(2D):
// ray = [[-30,0],[30,30]];
// pt = [-32,-10];
// p2 = ray_closest_point(ray,pt);
// stroke(ray, endcap2="arrow2");
// color("blue") translate(pt) circle(r=1,$fn=12);
// color("red") translate(p2) circle(r=1,$fn=12);
// Example(2D):
// ray = [[-30,0],[30,30]];
// pt = [-5,0];
// p2 = ray_closest_point(ray,pt);
// stroke(ray, endcap2="arrow2");
// color("blue") translate(pt) circle(r=1,$fn=12);
// color("red") translate(p2) circle(r=1,$fn=12);
// Example(2D):
// ray = [[-30,0],[30,30]];
// pt = [40,25];
// p2 = ray_closest_point(ray,pt);
// stroke(ray, endcap2="arrow2");
// color("blue") translate(pt) circle(r=1,$fn=12);
// color("red") translate(p2) circle(r=1,$fn=12);
// Example(FlatSpin):
// ray = [[-30,-15,0],[30,15,30]];
// pt = [5,5,5];
// p2 = ray_closest_point(ray,pt);
// stroke(ray, endcap2="arrow2");
// color("blue") translate(pt) sphere(r=1,$fn=12);
// color("red") translate(p2) sphere(r=1,$fn=12);
// Example(FlatSpin):
// ray = [[-30,-15,0],[30,15,30]];
// pt = [-35,-15,0];
// p2 = ray_closest_point(ray,pt);
// stroke(ray, endcap2="arrow2");
// color("blue") translate(pt) sphere(r=1,$fn=12);
// color("red") translate(p2) sphere(r=1,$fn=12);
// Example(FlatSpin):
// ray = [[-30,-15,0],[30,15,30]];
// pt = [40,15,25];
// p2 = ray_closest_point(ray,pt);
// stroke(ray, endcap2="arrow2");
// color("blue") translate(pt) sphere(r=1,$fn=12);
// color("red") translate(p2) sphere(r=1,$fn=12);
function ray_closest_point(ray,pt) =
assert(is_path(ray)&&len(ray)==2)
assert(same_shape(pt,ray[0]))
assert(!approx(ray[0],ray[1]))
let(
seglen = norm(ray[1]-ray[0]),
segvec = (ray[1]-ray[0])/seglen,
projection = (pt-ray[0]) * segvec
)
projection<=0 ? ray[0] :
ray[0] + projection*segvec;
// Function: segment_closest_point()
// Usage:
// segment_closest_point(seg,pt);
// Description:
// Returns the point on the given line segment `seg` that is closest to the given point `pt`.
// Returns the point on the given 2D or 3D line segment `seg` that is closest to the given point `pt`.
// The `seg` and `pt` args should either both be 2D or both 3D.
// Arguments:
// seg = A list of two points that are the endpoints of the bounded line segment.
// pt = The point to find the closest point on the segment to.
// Example(2D):
// seg = [[-30,0],[30,30]];
// pt = [-32,-10];
// p2 = segment_closest_point(seg,pt);
// stroke(seg);
// color("blue") translate(pt) circle(r=1,$fn=12);
// color("red") translate(p2) circle(r=1,$fn=12);
// Example(2D):
// seg = [[-30,0],[30,30]];
// pt = [-5,0];
// p2 = segment_closest_point(seg,pt);
// stroke(seg);
// color("blue") translate(pt) circle(r=1,$fn=12);
// color("red") translate(p2) circle(r=1,$fn=12);
// Example(2D):
// seg = [[-30,0],[30,30]];
// pt = [40,25];
// p2 = segment_closest_point(seg,pt);
// stroke(seg);
// color("blue") translate(pt) circle(r=1,$fn=12);
// color("red") translate(p2) circle(r=1,$fn=12);
// Example(FlatSpin):
// seg = [[-30,-15,0],[30,15,30]];
// pt = [5,5,5];
// p2 = segment_closest_point(seg,pt);
// stroke(seg);
// color("blue") translate(pt) sphere(r=1,$fn=12);
// color("red") translate(p2) sphere(r=1,$fn=12);
// Example(FlatSpin):
// seg = [[-30,-15,0],[30,15,30]];
// pt = [-35,-15,0];
// p2 = segment_closest_point(seg,pt);
// stroke(seg);
// color("blue") translate(pt) sphere(r=1,$fn=12);
// color("red") translate(p2) sphere(r=1,$fn=12);
// Example(FlatSpin):
// seg = [[-30,-15,0],[30,15,30]];
// pt = [40,15,25];
// p2 = segment_closest_point(seg,pt);
// stroke(seg);
// 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(is_path(seg)&&len(seg)==2)
assert(same_shape(pt,seg[0]))
approx(seg[0],seg[1])? seg[0] :
let(
n = line_normal(seg),
isect = _general_line_intersection(seg,[pt,pt+n])
seglen = norm(seg[1]-seg[0]),
segvec = (seg[1]-seg[0])/seglen,
projection = (pt-seg[0]) * segvec
)
norm(n)==0? seg[0] :
isect[1]<=0? seg[0] :
isect[1]>=1? seg[1] :
isect[0];
projection<=0 ? seg[0] :
projection>=seglen ? seg[1] :
seg[0] + projection*segvec;
// Section: 2D Triangles
@@ -1221,15 +1380,18 @@ function find_noncollinear_points(points) =
// Usage:
// pointlist_bounds(pts);
// Description:
// Finds the bounds containing all the 2D or 3D points in `pts`.
// Returns `[[MINX, MINY, MINZ], [MAXX, MAXY, MAXZ]]`
// Finds the bounds containing all the points in `pts` which can be a list of points in any dimension.
// Returns a list of two items: a list of the minimums and a list of the maximums. For example, with
// 3d points `[[MINX, MINY, MINZ], [MAXX, MAXY, MAXZ]]`
// Arguments:
// pts = List of points.
function pointlist_bounds(pts) = [
[for (a=[0:2]) min([ for (x=pts) point3d(x)[a] ]) ],
[for (a=[0:2]) max([ for (x=pts) point3d(x)[a] ]) ]
];
function pointlist_bounds(pts) =
assert(is_matrix(pts))
let(ptsT = transpose(pts))
[
[for(row=ptsT) min(row)],
[for(row=ptsT) max(row)]
];
// Function: closest_point()
// Usage:

147
math.scad
View File

@@ -67,7 +67,7 @@ function hypot(x,y,z=0) = norm([x,y,z]);
// Usage:
// x = factorial(n,[d]);
// Description:
// Returns the factorial of the given integer value.
// Returns the factorial of the given integer value, or n!/d! if d is given.
// Arguments:
// n = The integer number to get the factorial of. (n!)
// d = If given, the returned value will be (n! / d!)
@@ -75,7 +75,10 @@ function hypot(x,y,z=0) = norm([x,y,z]);
// x = factorial(4); // Returns: 24
// y = factorial(6); // Returns: 720
// z = factorial(9); // Returns: 362880
function factorial(n,d=1) = product([for (i=[n:-1:d]) i]);
function factorial(n,d=0) =
assert(n>=0 && d>=0, "Factorial is not defined for negative numbers")
assert(d<=n, "d cannot be larger than n")
product([1,for (i=[n:-1:d+1]) i]);
// Function: lerp()
@@ -525,6 +528,17 @@ function deltas(v) = [for (p=pair(v)) p.y-p.x];
function product(v, i=0, tot=undef) = i>=len(v)? tot : product(v, i+1, ((tot==undef)? v[i] : is_vector(v[i])? vmul(tot,v[i]) : tot*v[i]));
// Function: outer_product()
// 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))
assert(len(u)==len(v))
[for(i=[0:len(u)-1]) [for(j=[0:len(u)-1]) u[i]*v[j]]];
// Function: mean()
// Description:
// Returns the arithmatic mean/average of all entries in the given array.
@@ -563,7 +577,7 @@ function median(v) =
// 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 `undef`. If b is a matrix that is compatible with A
// 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.
@@ -582,7 +596,7 @@ function linear_solve(A,b) =
R = submatrix(qr[1],[0:mindim-1], [0:mindim-1]),
zeros = [for(i=[0:mindim-1]) if (approx(R[i][i],0)) i]
)
zeros != [] ? undef :
zeros != [] ? [] :
m<n ? Q*back_substitute(R,b,transpose=true) :
back_substitute(R, transpose(Q)*b);
@@ -604,7 +618,8 @@ function matrix_inverse(A) =
// Usage: submatrix(M, ind1, ind2)
// Description:
// Returns a submatrix with the specified index ranges or index sets.
function submatrix(M,ind1,ind2) = [for(i=ind1) [for(j=ind2) M[i][j] ] ];
function submatrix(M,ind1,ind2) =
[for(i=ind1) [for(j=ind2) M[i][j] ] ];
// Function: qr_factor()
@@ -635,7 +650,7 @@ function _qr_factor(A,Q, column, m, n) =
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*transpose([v])*[v],
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, column+1, m, n);
@@ -647,11 +662,12 @@ function _qr_factor(A,Q, column, m, n) =
// Solves the problem Rx=b where R is an upper triangular square matrix. No check is made that the lower triangular entries
// are actually zero. 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.
// 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, x=[],transpose = false) =
assert(is_matrix(R, square=true))
let(n=len(R))
assert(is_vector(b,n) || is_matrix(b,n),"R and b are not compatible in back_substitute")
assert(is_vector(b,n) || is_matrix(b,n),str("R and b are not compatible in back_substitute ",n, len(b)))
!is_vector(b) ? transpose([for(i=[0:len(b[0])-1]) back_substitute(R,subindex(b,i),transpose=transpose)]) :
transpose?
reverse(back_substitute(
@@ -660,7 +676,10 @@ function back_substitute(R, b, x=[],transpose = false) =
)) :
len(x) == n ? x :
let(
ind = n - len(x) - 1,
ind = n - len(x) - 1
)
R[ind][ind] == 0 ? [] :
let(
newvalue =
len(x)==0? b[ind]/R[ind][ind] :
(b[ind]-select(R[ind],ind+1,-1) * x)/R[ind][ind]
@@ -733,7 +752,7 @@ function determinant(M) =
// m = optional height of matrix
// n = optional width of matrix
// square = set to true to require a square matrix. Default: false
function is_matrix(A,n,m,square=false) =
function is_matrix(A,m,n,square=false) =
is_vector(A[0],n) && is_vector(A*(0*A[0]),m) &&
(!square || len(A)==len(A[0]));
@@ -1037,12 +1056,82 @@ function C_div(z1,z2) = let(den = z2.x*z2.x + z2.y*z2.y)
// Evaluates specified real polynomial, p, at the complex or real input value, z.
// The polynomial is specified as p=[a_n, a_{n-1},...,a_1,a_0]
// where a_n is the z^n coefficient. Polynomial coefficients are real.
// Note: this should probably be recoded to use division by [1,-z], which is more accurate
// and avoids overflow with large coefficients, but requires poly_div to support complex coefficients.
function polynomial(p, z, k, zk, total) =
is_undef(k) ? polynomial(p, z, len(p)-1, is_num(z)? 1 : [1,0], is_num(z) ? 0 : [0,0]) :
k==-1 ? total :
polynomial(p, z, k-1, is_num(z) ? zk*z : C_times(zk,z), total+zk*p[k]);
// Function: poly_mult()
// Usage
// polymult(p,q)
// polymult([p1,p2,p3,...])
// Descriptoin:
// Given a list of polynomials represented as real coefficient lists, with the highest degree coefficient first,
// computes the coefficient list of the product polynomial.
function poly_mult(p,q) =
is_undef(q) ?
assert(is_list(p) && (is_vector(p[0]) || p[0]==[]), "Invalid arguments to poly_mult")
len(p)==2 ? poly_mult(p[0],p[1])
: poly_mult(p[0], poly_mult(select(p,1,-1)))
:
_poly_trim(
[
for(n = [len(p)+len(q)-2:-1:0])
sum( [for(i=[0:1:len(p)-1])
let(j = len(p)+len(q)- 2 - n - i)
if (j>=0 && j<len(q)) p[i]*q[j]
])
]);
// Function: poly_div()
// Usage:
// [quotient,remainder] = poly_div(n,d)
// 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 [].
function poly_div(n,d,q=[]) =
assert(len(d)>0 && d[0]!=0 , "Denominator is zero or has leading zero coefficient")
len(n)<len(d) ? [q,_poly_trim(n)] :
let(
t = n[0] / d[0],
newq = concat(q,[t]),
newn = [for(i=[1:1:len(n)-1]) i<len(d) ? n[i] - t*d[i] : n[i]]
)
poly_div(newn,d,newq);
// Internal Function: _poly_trim()
// Usage:
// _poly_trim(p,[eps])
// Description:
// Removes leading zero terms of a polynomial. By default zeros must be exact,
// or give epsilon for approximate zeros.
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 ? [] : select(p,nz[0],-1);
// Function: poly_add()
// Usage:
// sum = poly_add(p,q)
// Description:
// Computes the sum of two polynomials.
function poly_add(p,q) =
let( plen = len(p),
qlen = len(q),
long = plen>qlen ? p : q,
short = plen>qlen ? q : p
)
_poly_trim(long + concat(repeat(0,len(long)-len(short)),short));
// Function: poly_roots()
// Usage:
// poly_roots(p,[tol])
@@ -1062,14 +1151,18 @@ function polynomial(p, z, k, zk, total) =
// Dario Bini. "Numerical computation of polynomial zeros by means of Aberth's Method", Numerical Algorithms, Feb 1996.
// https://www.researchgate.net/publication/225654837_Numerical_computation_of_polynomial_zeros_by_means_of_Aberth's_method
function poly_roots(p,tol=1e-14) =
function poly_roots(p,tol=1e-14,error_bound=false) =
assert(p!=[], "Input polynomial must have a nonzero coefficient")
assert(is_vector(p), "Input must be a vector")
p[0] == 0 ? poly_roots(slice(p,1,-1)) : // Strip leading zero coefficients
p[len(p)-1] == 0 ? [[0,0], // Strip trailing zero coefficients
each poly_roots(select(p,0,-2))] :
len(p)==1 ? [] : // Nonzero constant case has no solutions
len(p)==2 ? [[-p[1]/p[0],0]] : // Linear case needs special handling
p[0] == 0 ? poly_roots(slice(p,1,-1),tol=tol,error_bound=error_bound) : // Strip leading zero coefficients
p[len(p)-1] == 0 ? // Strip trailing zero coefficients
let( solutions = poly_roots(select(p,0,-2),tol=tol, error_bound=error_bound))
(error_bound ? [ [[0,0], each solutions[0]], [0, each solutions[1]]]
: [[0,0], each solutions]) :
len(p)==1 ? (error_bound ? [[],[]] : []) : // Nonzero constant case has no solutions
len(p)==2 ? let( solution = [[-p[1]/p[0],0]]) // Linear case needs special handling
(error_bound ? [solution,[0]] : solution)
:
let(
n = len(p)-1, // polynomial degree
pderiv = [for(i=[0:n-1]) p[i]*(n-i)],
@@ -1082,9 +1175,12 @@ function poly_roots(p,tol=1e-14) =
init = [for(i=[0:1:n-1]) // Initial guess for roots
let(angle = 360*i/n+270/n/PI)
[beta,0]+r*[cos(angle),sin(angle)]
]
],
roots = _poly_roots(p,pderiv,s,init,tol=tol),
error = error_bound ? [for(xi=roots) n * (norm(polynomial(p,xi))+tol*polynomial(s,norm(xi))) /
abs(norm(polynomial(pderiv,xi))-tol*polynomial(s,norm(xi)))] : 0
)
_poly_roots(p,pderiv,s,init,tol=tol);
error_bound ? [roots, error] : roots;
// p = polynomial
// pderiv = derivative polynomial of p
@@ -1114,7 +1210,10 @@ function _poly_roots(p, pderiv, s, z, tol, i=0) =
// The polynomial is specified as p=[a_n, a_{n-1},...,a_1,a_0]
// where a_n is the x^n coefficient. This function works by
// computing the complex roots and returning those roots where
// the imaginary part is closed to zero, specifically: z.y/(1+norm(z)) < eps. Because
// the imaginary part is closed to zero. By default it uses a computed
// error bound from the polynomial solver to decide whether imaginary
// parts are zero. You can specify eps, in which case the test is
// z.y/(1+norm(z)) < eps. Because
// of poor convergence and higher error for repeated roots, such roots may
// be missed by the algorithm because their imaginary part is large.
// Arguments:
@@ -1122,11 +1221,13 @@ function _poly_roots(p, pderiv, s, z, tol, i=0) =
// eps = used to determine whether imaginary parts of roots are zero
// tol = tolerance for the complex polynomial root finder
function real_roots(p,eps=EPSILON,tol=1e-14) =
function real_roots(p,eps=undef,tol=1e-14) =
let(
roots = poly_roots(p)
roots_err = poly_roots(p,error_bound=true),
roots = roots_err[0],
err = roots_err[1]
)
[for(z=roots) if (abs(z.y)/(1+norm(z))<eps) z.x];
is_def(eps) ? [for(z=roots) if (abs(z.y)/(1+norm(z))<eps) z.x]
: [for(i=idx(roots)) if (abs(roots[i].y)<=err[i]) roots[i].x];
// vim: expandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap

483
quaternions.scad
View File

@@ -20,59 +20,107 @@
// Internal
function _Quat(a,s,w) = [a[0]*s, a[1]*s, a[2]*s, w];
function _Qvec(q) = [q.x,q.y,q.z];
function _Qreal(q) = q[3];
function _Qset(v,r) = concat( v, r );
// normalizes without checking
function _Qnorm(q) = q/norm(q);
// Function: Q_is_quat()
// Usage:
// if(Q_is_quat(q)) a=0;
// Description: Return true if q is a valid non-zero quaternion.
// Arguments:
// q = object to check.
function Q_is_quat(q) = is_vector(q,4) && ! approx(norm(q),0) ;
// Function: Quat()
// Usage:
// Quat(ax, ang);
// Description: Create a new Quaternion from axis and angle of rotation.
// Description: Create a normalized Quaternion from axis and angle of rotation.
// Arguments:
// ax = Vector of axis of rotation.
// ang = Number of degrees to rotate around the axis counter-clockwise, when facing the origin.
function Quat(ax=[0,0,1], ang=0) = _Quat(ax/norm(ax), sin(ang/2), cos(ang/2));
function Quat(ax=[0,0,1], ang=0) =
assert( is_vector(ax,3) && is_finite(ang), "Invalid input")
let( n = norm(ax) )
approx(n,0)
? _Quat([0,0,0], sin(ang/2), cos(ang/2))
: _Quat(ax/n, sin(ang/2), cos(ang/2));
// Function: QuatX()
// Usage:
// QuatX(a);
// Description: Create a new Quaternion for rotating around the X axis [1,0,0].
// Description: Create a normalized Quaternion for rotating around the X axis [1,0,0].
// Arguments:
// a = Number of degrees to rotate around the axis counter-clockwise, when facing the origin.
function QuatX(a=0) = Quat([1,0,0],a);
function QuatX(a=0) =
assert( is_finite(a), "Invalid angle" )
Quat([1,0,0],a);
// Function: QuatY()
// Usage:
// QuatY(a);
// Description: Create a new Quaternion for rotating around the Y axis [0,1,0].
// Description: Create a normalized Quaternion for rotating around the Y axis [0,1,0].
// Arguments:
// a = Number of degrees to rotate around the axis counter-clockwise, when facing the origin.
function QuatY(a=0) = Quat([0,1,0],a);
function QuatY(a=0) =
assert( is_finite(a), "Invalid angle" )
Quat([0,1,0],a);
// Function: QuatZ()
// Usage:
// QuatZ(a);
// Description: Create a new Quaternion for rotating around the Z axis [0,0,1].
// Description: Create a normalized Quaternion for rotating around the Z axis [0,0,1].
// Arguments:
// a = Number of degrees to rotate around the axis counter-clockwise, when facing the origin.
function QuatZ(a=0) = Quat([0,0,1],a);
function QuatZ(a=0) =
assert( is_finite(a), "Invalid angle" )
Quat([0,0,1],a);
// Function: QuatXYZ()
// Usage:
// QuatXYZ([X,Y,Z])
// Description:
// Creates a quaternion from standard [X,Y,Z] rotation angles in degrees.
// Creates a normalized quaternion from standard [X,Y,Z] rotation angles in degrees.
// Arguments:
// a = The triplet of rotation angles, [X,Y,Z]
function QuatXYZ(a=[0,0,0]) =
assert( is_vector(a,3), "Invalid angles")
let(
qx = QuatX(a[0]),
qy = QuatY(a[1]),
qz = QuatZ(a[2])
qx = QuatX(a[0]),
qy = QuatY(a[1]),
qz = QuatZ(a[2])
)
Q_Mul(qz, Q_Mul(qy, qx));
// Function: Q_From_to()
// Usage:
// q = Q_From_to(v1, v2);
// Description:
// Returns the normalized quaternion that rotates the non zero 3D vector v1
// to the non zero 3D vector v2.
function Q_From_to(v1, v2) =
assert( is_vector(v1,3) && is_vector(v2,3)
&& ! approx(norm(v1),0) && ! approx(norm(v2),0)
, "Invalid vector(s)")
let( ax = cross(v1,v2),
n = norm(ax) )
approx(n, 0)
? v1*v2>0 ? Q_Ident() : Quat([ v1.y, -v1.x, 0], 180)
: Quat(ax, atan2( n , v1*v2 ));
// Function: Q_Ident()
// Description: Returns the "Identity" zero-rotation Quaternion.
function Q_Ident() = [0, 0, 0, 1];
@@ -81,55 +129,85 @@ function Q_Ident() = [0, 0, 0, 1];
// Function: Q_Add_S()
// Usage:
// Q_Add_S(q, s)
// Description: Adds a scalar value `s` to the W part of a quaternion `q`.
function Q_Add_S(q, s) = q+[0,0,0,s];
// Description:
// Adds a scalar value `s` to the W part of a quaternion `q`.
// The returned quaternion is usually not normalized.
function Q_Add_S(q, s) =
assert( is_finite(s), "Invalid scalar" )
q+[0,0,0,s];
// Function: Q_Sub_S()
// Usage:
// Q_Sub_S(q, s)
// Description: Subtracts a scalar value `s` from the W part of a quaternion `q`.
function Q_Sub_S(q, s) = q-[0,0,0,s];
// Description:
// Subtracts a scalar value `s` from the W part of a quaternion `q`.
// The returned quaternion is usually not normalized.
function Q_Sub_S(q, s) =
assert( is_finite(s), "Invalid scalar" )
q-[0,0,0,s];
// Function: Q_Mul_S()
// Usage:
// Q_Mul_S(q, s)
// Description: Multiplies each part of a quaternion `q` by a scalar value `s`.
function Q_Mul_S(q, s) = q*s;
// Description:
// Multiplies each part of a quaternion `q` by a scalar value `s`.
// The returned quaternion is usually not normalized.
function Q_Mul_S(q, s) =
assert( is_finite(s), "Invalid scalar" )
q*s;
// Function: Q_Div_S()
// Usage:
// Q_Div_S(q, s)
// Description: Divides each part of a quaternion `q` by a scalar value `s`.
function Q_Div_S(q, s) = q/s;
// Description:
// Divides each part of a quaternion `q` by a scalar value `s`.
// The returned quaternion is usually not normalized.
function Q_Div_S(q, s) =
assert( is_finite(s) && ! approx(s,0) , "Invalid scalar" )
q/s;
// Function: Q_Add()
// Usage:
// Q_Add(a, b)
// Description: Adds each part of two quaternions together.
function Q_Add(a, b) = a+b;
// Description:
// Adds each part of two quaternions together.
// The returned quaternion is usually not normalized.
function Q_Add(a, b) =
assert( Q_is_quat(a) && Q_is_quat(a), "Invalid quaternion(s)")
assert( ! approx(norm(a+b),0), "Quaternions cannot be opposed" )
a+b;
// Function: Q_Sub()
// Usage:
// Q_Sub(a, b)
// Description: Subtracts each part of quaternion `b` from quaternion `a`.
function Q_Sub(a, b) = a-b;
// Description:
// Subtracts each part of quaternion `b` from quaternion `a`.
// The returned quaternion is usually not normalized.
function Q_Sub(a, b) =
assert( Q_is_quat(a) && Q_is_quat(a), "Invalid quaternion(s)")
assert( ! approx(a,b), "Quaternions cannot be equal" )
a-b;
// Function: Q_Mul()
// Usage:
// Q_Mul(a, b)
// Description: Multiplies quaternion `a` by quaternion `b`.
function Q_Mul(a, b) = [
a[3]*b.x + a.x*b[3] + a.y*b.z - a.z*b.y,
a[3]*b.y - a.x*b.z + a.y*b[3] + a.z*b.x,
a[3]*b.z + a.x*b.y - a.y*b.x + a.z*b[3],
a[3]*b[3] - a.x*b.x - a.y*b.y - a.z*b.z,
];
// Description:
// Multiplies quaternion `a` by quaternion `b`.
// The returned quaternion is normalized if both `a` and `b` are normalized
function Q_Mul(a, b) =
assert( Q_is_quat(a) && Q_is_quat(b), "Invalid quaternion(s)")
[
a[3]*b.x + a.x*b[3] + a.y*b.z - a.z*b.y,
a[3]*b.y - a.x*b.z + a.y*b[3] + a.z*b.x,
a[3]*b.z + a.x*b.y - a.y*b.x + a.z*b[3],
a[3]*b[3] - a.x*b.x - a.y*b.y - a.z*b.z,
];
// Function: Q_Cumulative()
@@ -139,6 +217,8 @@ function Q_Mul(a, b) = [
// Given a list of Quaternions, cumulatively multiplies them, returning a list
// of each cumulative Quaternion product. It starts with the first quaternion
// given in the list, and applies successive quaternion rotations in list order.
// The quaternion in the returned list are normalized if each quaternion in v
// is normalized.
function Q_Cumulative(v, _i=0, _acc=[]) =
_i==len(v) ? _acc :
Q_Cumulative(
@@ -154,42 +234,65 @@ function Q_Cumulative(v, _i=0, _acc=[]) =
// Usage:
// Q_Dot(a, b)
// Description: Calculates the dot product between quaternions `a` and `b`.
function Q_Dot(a, b) = a[0]*b[0] + a[1]*b[1] + a[2]*b[2] + a[3]*b[3];
function Q_Dot(a, b) =
assert( Q_is_quat(a) && Q_is_quat(b), "Invalid quaternion(s)" )
a*b;
// Function: Q_Neg()
// Usage:
// Q_Neg(q)
// Description: Returns the negative of quaternion `q`.
function Q_Neg(q) = -q;
function Q_Neg(q) =
assert( Q_is_quat(q), "Invalid quaternion" )
-q;
// Function: Q_Conj()
// Usage:
// Q_Conj(q)
// Description: Returns the conjugate of quaternion `q`.
function Q_Conj(q) = [-q.x, -q.y, -q.z, q[3]];
function Q_Conj(q) =
assert( Q_is_quat(q), "Invalid quaternion" )
[-q.x, -q.y, -q.z, q[3]];
// Function: Q_Inverse()
// Usage:
// qc = Q_Inverse(q)
// Description: Returns the multiplication inverse of quaternion `q` that is normalized only if `q` is normalized.
function Q_Inverse(q) =
assert( Q_is_quat(q), "Invalid quaternion" )
let(q = _Qnorm(q) )
[-q.x, -q.y, -q.z, q[3]];
// Function: Q_Norm()
// Usage:
// Q_Norm(q)
// Description: Returns the `norm()` "length" of quaternion `q`.
function Q_Norm(q) = norm(q);
// Description:
// Returns the `norm()` "length" of quaternion `q`.
// Normalized quaternions have unitary norm.
function Q_Norm(q) =
assert( Q_is_quat(q), "Invalid quaternion" )
norm(q);
// Function: Q_Normalize()
// Usage:
// Q_Normalize(q)
// Description: Normalizes quaternion `q`, so that norm([W,X,Y,Z]) == 1.
function Q_Normalize(q) = q/norm(q);
function Q_Normalize(q) =
assert( Q_is_quat(q) , "Invalid quaternion" )
q/norm(q);
// Function: Q_Dist()
// Usage:
// Q_Dist(q1, q2)
// Description: Returns the "distance" between two quaternions.
function Q_Dist(q1, q2) = norm(q2-q1);
function Q_Dist(q1, q2) =
assert( Q_is_quat(q1) && Q_is_quat(q2), "Invalid quaternion(s)" )
norm(q2-q1);
// Function: Q_Slerp()
@@ -214,24 +317,24 @@ function Q_Dist(q1, q2) = norm(q2-q1);
// for (q = Q_Slerp(a, b, [0:0.1:1]))
// Qrot(q) right(80) cube([10,10,1]);
// #sphere(r=80);
function Q_Slerp(q1, q2, u) =
assert(is_num(u) || is_num(u[0]))
!is_num(u)? [for (uu=u) Q_Slerp(q1,q2,uu)] :
let(
q1 = Q_Normalize(q1),
q2 = Q_Normalize(q2),
dot = Q_Dot(q1, q2)
) let(
q2 = dot<0? Q_Neg(q2) : q2,
dot = dot<0? -dot : dot
) (dot>0.9995)? Q_Normalize(q1 + (u * (q2-q1))) :
let(
dot = constrain(dot,-1,1),
theta_0 = acos(dot),
theta = theta_0 * u,
q3 = Q_Normalize(q2 - q1*dot),
out = q1*cos(theta) + q3*sin(theta)
) out;
function Q_Slerp(q1, q2, u, _dot) =
is_undef(_dot)
? assert(is_finite(u) || is_range(u) || is_vector(u), "Invalid interpolation coefficient(s)")
assert(Q_is_quat(q1) && Q_is_quat(q2), "Invalid quaternion(s)" )
let(
_dot = q1*q2,
q1 = q1/norm(q1),
q2 = _dot<0 ? -q2/norm(q2) : q2/norm(q2),
dot = abs(_dot)
)
! is_finite(u) ? [for (uu=u) Q_Slerp(q1, q2, uu, dot)] :
Q_Slerp(q1, q2, u, dot)
: _dot>0.9995
? _Qnorm(q1 + u*(q2-q1))
: let( theta = u*acos(_dot),
q3 = _Qnorm(q2 - _dot*q1)
)
_Qnorm(q1*cos(theta) + q3*sin(theta));
// Function: Q_Matrix3()
@@ -239,11 +342,13 @@ function Q_Slerp(q1, q2, u) =
// Q_Matrix3(q);
// Description:
// Returns the 3x3 rotation matrix for the given normalized quaternion q.
function Q_Matrix3(q) = [
[1-2*q[1]*q[1]-2*q[2]*q[2], 2*q[0]*q[1]-2*q[2]*q[3], 2*q[0]*q[2]+2*q[1]*q[3]],
[ 2*q[0]*q[1]+2*q[2]*q[3], 1-2*q[0]*q[0]-2*q[2]*q[2], 2*q[1]*q[2]-2*q[0]*q[3]],
[ 2*q[0]*q[2]-2*q[1]*q[3], 2*q[1]*q[2]+2*q[0]*q[3], 1-2*q[0]*q[0]-2*q[1]*q[1]]
];
function Q_Matrix3(q) =
let( q = Q_Normalize(q) )
[
[1-2*q[1]*q[1]-2*q[2]*q[2], 2*q[0]*q[1]-2*q[2]*q[3], 2*q[0]*q[2]+2*q[1]*q[3]],
[ 2*q[0]*q[1]+2*q[2]*q[3], 1-2*q[0]*q[0]-2*q[2]*q[2], 2*q[1]*q[2]-2*q[0]*q[3]],
[ 2*q[0]*q[2]-2*q[1]*q[3], 2*q[1]*q[2]+2*q[0]*q[3], 1-2*q[0]*q[0]-2*q[1]*q[1]]
];
// Function: Q_Matrix4()
@@ -251,12 +356,14 @@ function Q_Matrix3(q) = [
// Q_Matrix4(q);
// Description:
// Returns the 4x4 rotation matrix for the given normalized quaternion q.
function Q_Matrix4(q) = [
[1-2*q[1]*q[1]-2*q[2]*q[2], 2*q[0]*q[1]-2*q[2]*q[3], 2*q[0]*q[2]+2*q[1]*q[3], 0],
[ 2*q[0]*q[1]+2*q[2]*q[3], 1-2*q[0]*q[0]-2*q[2]*q[2], 2*q[1]*q[2]-2*q[0]*q[3], 0],
[ 2*q[0]*q[2]-2*q[1]*q[3], 2*q[1]*q[2]+2*q[0]*q[3], 1-2*q[0]*q[0]-2*q[1]*q[1], 0],
[ 0, 0, 0, 1]
];
function Q_Matrix4(q) =
let( q = Q_Normalize(q) )
[
[1-2*q[1]*q[1]-2*q[2]*q[2], 2*q[0]*q[1]-2*q[2]*q[3], 2*q[0]*q[2]+2*q[1]*q[3], 0],
[ 2*q[0]*q[1]+2*q[2]*q[3], 1-2*q[0]*q[0]-2*q[2]*q[2], 2*q[1]*q[2]-2*q[0]*q[3], 0],
[ 2*q[0]*q[2]-2*q[1]*q[3], 2*q[1]*q[2]+2*q[0]*q[3], 1-2*q[0]*q[0]-2*q[1]*q[1], 0],
[ 0, 0, 0, 1]
];
// Function: Q_Axis()
@@ -264,16 +371,28 @@ function Q_Matrix4(q) = [
// Q_Axis(q)
// Description:
// Returns the axis of rotation of a normalized quaternion `q`.
function Q_Axis(q) = let(d = sqrt(1-(q[3]*q[3]))) (d==0)? [0,0,1] : [q[0]/d, q[1]/d, q[2]/d];
// The input doesn't need to be normalized.
function Q_Axis(q) =
assert( Q_is_quat(q) , "Invalid quaternion" )
let( d = norm(_Qvec(q)) )
approx(d,0)? [0,0,1] : _Qvec(q)/d;
// Function: Q_Angle()
// Usage:
// Q_Angle(q)
// a = Q_Angle(q)
// a12 = Q_Angle(q1,q2);
// Description:
// Returns the angle of rotation (in degrees) of a normalized quaternion `q`.
function Q_Angle(q) = 2 * acos(q[3]);
// If only q1 is given, returns the angle of rotation (in degrees) of that quaternion.
// If both q1 and q2 are given, returns the angle (in degrees) between them.
// The input quaternions don't need to be normalized.
function Q_Angle(q1,q2) =
assert(Q_is_quat(q1) && (is_undef(q2) || Q_is_quat(q2)), "Invalid quaternion(s)" )
let( n1 = is_undef(q2)? norm(_Qvec(q1)): norm(q1) )
is_undef(q2)
? 2 * atan2(n1,_Qreal(q1))
: let( q1 = q1/norm(q1),
q2 = q2/norm(q2) )
4 * atan2(norm(q1 - q2), norm(q1 + q2));
// Function&Module: Qrot()
// Usage: As Module
@@ -305,9 +424,9 @@ module Qrot(q) {
}
function Qrot(q,p) =
is_undef(p)? Q_Matrix4(q) :
is_vector(p)? Qrot(q,[p])[0] :
apply(Q_Matrix4(q), p);
is_undef(p)? Q_Matrix4(q) :
is_vector(p)? Qrot(q,[p])[0] :
apply(Q_Matrix4(q), p);
// Module: Qrot_copies()
@@ -327,4 +446,214 @@ function Qrot(q,p) =
module Qrot_copies(quats) for (q=quats) Qrot(q) children();
// Function: Q_Rotation()
// Usage:
// Q_Rotation(R)
// Description:
// Returns a normalized quaternion corresponding to the rotation matrix R.
// R may be a 3x3 rotation matrix or a homogeneous 4x4 rotation matrix.
// The last row and last column of R are ignored for 4x4 matrices.
// It doesn't check whether R is in fact a rotation matrix.
// If R is not a rotation, the returned quaternion is an unpredictable quaternion .
function Q_Rotation(R) =
assert( is_matrix(R,3,3) || is_matrix(R,4,4) ,
"Matrix is neither 3x3 nor 4x4")
let( tr = R[0][0]+R[1][1]+R[2][2] ) // R trace
tr>0
? let( r = 1+tr )
_Qnorm( _Qset([ R[1][2]-R[2][1], R[2][0]-R[0][2], R[0][1]-R[1][0] ], -r ) )
: let( i = max_index([ R[0][0], R[1][1], R[2][2] ]),
r = 1 + 2*R[i][i] -R[0][0] -R[1][1] -R[2][2] )
i==0 ? _Qnorm( _Qset( [ 4*r, (R[1][0]+R[0][1]), (R[0][2]+R[2][0]) ], (R[2][1]-R[1][2])) ):
i==1 ? _Qnorm( _Qset( [ (R[1][0]+R[0][1]), 4*r, (R[2][1]+R[1][2]) ], (R[0][2]-R[2][0])) ):
_Qnorm( _Qset( [ (R[2][0]+R[0][2]), (R[1][2]+R[2][1]), 4*r ], (R[1][0]-R[0][1])) ) ;
// Function&Module: Q_Rotation_path(q1, n, [q2])
// Usage: As a function
// path = Q_Rotation_path(q1, n, q2);
// path = Q_Rotation_path(q1, n);
// Usage: As a module
// Q_Rotation_path(q1, n, q2) ...
// Description:
// If q2 is undef and it is called as a function, the path, with length n+1 (n>=1), will be the
// cumulative multiplications of the matrix rotation of q1 by itself.
// If q2 is defined and it is called as a function, returns a rotation matrix path of length n+1 (n>=1)
// that interpolates two given rotation quaternions. The first matrix of the sequence is the
// matrix rotation of q1 and the last one, the matrix rotation of q2. The intermediary matrix
// rotations are an uniform interpolation of the path extreme matrices.
// When called as a module, applies to its children() each rotation of the sequence computed
// by the function.
// The input quaternions don't need to be normalized.
// Arguments:
// q1 = The quaternion of the first rotation.
// q2 = The quaternion of the last rotation.
// n = An integer defining the path length ( path length = n+1).
// Example(3D): as a function
// a = QuatY(-135);
// b = QuatXYZ([0,-30,30]);
// for (M=Q_Rotation_path(a, 10, b))
// multmatrix(M)
// right(80) cube([10,10,1]);
// #sphere(r=80);
// Example(3D): as a module
// a = QuatY(-135);
// b = QuatXYZ([0,-30,30]);
// Q_Rotation_path(a, 10, b)
// right(80) cube([10,10,1]);
// #sphere(r=80);
// Example(3D): as a function
// a = QuatY(5);
// for (M=Q_Rotation_path(a, 10))
// multmatrix(M)
// right(80) cube([10,10,1]);
// #sphere(r=80);
// Example(3D): as a module
// a = QuatY(5);
// Q_Rotation_path(a, 10)
// right(80) cube([10,10,1]);
// #sphere(r=80);
function Q_Rotation_path(q1, n=1, q2) =
assert( Q_is_quat(q1) && (is_undef(q2) || Q_is_quat(q2) ), "Invalid quaternion(s)" )
assert( is_finite(n) && n>=1 && n==floor(n), "Invalid integer" )
assert( is_undef(q2) || ! approx(norm(q1+q2),0), "Quaternions cannot be opposed" )
is_undef(q2)
? [for( i=0, dR=Q_Matrix4(q1), R=dR; i<=n; i=i+1, R=dR*R ) R]
: let( q2 = Q_Normalize( q1*q2<0 ? -q2: q2 ),
dq = Q_pow( Q_Mul( q2, Q_Inverse(q1) ), 1/n ),
dR = Q_Matrix4(dq) )
[for( i=0, R=Q_Matrix4(q1); i<=n; i=i+1, R=dR*R ) R];
module Q_Rotation_path(q1, n=1, q2) {
for(Mi=Q_Rotation_path(q1, n, q2))
multmatrix(Mi)
children();
}
// Function: Q_Nlerp()
// Usage:
// q = Q_Nlerp(q1, q2, u);
// Description:
// Returns a quaternion that is a normalized linear interpolation between two quaternions
// when u is a number.
// If u is a list of numbers, computes the interpolations for each value in the
// list and returns the interpolated quaternions in a list.
// The input quaternions don't need to be normalized.
// Arguments:
// q1 = The first quaternion. (u=0)
// q2 = The second quaternion. (u=1)
// u = A value (or a list of values), between 0 and 1, of the proportion(s) of each quaternion in the interpolation.
// Example(3D): Giving `u` as a Scalar
// a = QuatY(-135);
// b = QuatXYZ([0,-30,30]);
// for (u=[0:0.1:1])
// Qrot(Q_Nlerp(a, b, u))
// right(80) cube([10,10,1]);
// #sphere(r=80);
// Example(3D): Giving `u` as a Range
// a = QuatZ(-135);
// b = QuatXYZ([90,0,-45]);
// for (q = Q_Nlerp(a, b, [0:0.1:1]))
// Qrot(q) right(80) cube([10,10,1]);
// #sphere(r=80);
function Q_Nlerp(q1,q2,u) =
assert(is_finite(u) || is_range(u) || is_vector(u) ,
"Invalid interpolation coefficient(s)" )
assert(Q_is_quat(q1) && Q_is_quat(q2), "Invalid quaternion(s)" )
assert( ! approx(norm(q1+q2),0), "Quaternions cannot be opposed" )
let( q1 = Q_Normalize(q1),
q2 = Q_Normalize(q2) )
is_num(u)
? _Qnorm((1-u)*q1 + u*q2 )
: [for (ui=u) _Qnorm((1-ui)*q1 + ui*q2 ) ];
// Function: Q_Squad()
// Usage:
// qn = Q_Squad(q1,q2,q3,q4,u);
// Description:
// Returns a quaternion that is a cubic spherical interpolation of the quaternions
// q1 and q4 taking the other two quaternions, q2 and q3, as parameter of a cubic
// on the sphere similar to the control points of a Bezier curve.
// If u is a number, usually between 0 and 1, returns the quaternion that results
// from the interpolation.
// If u is a list of numbers, computes the interpolations for each value in the
// list and returns the interpolated quaternions in a list.
// The input quaternions don't need to be normalized.
// Arguments:
// q1 = The start quaternion. (u=0)
// q1 = The first intermediate quaternion.
// q2 = The second intermediate quaternion.
// q4 = The end quaternion. (u=1)
// u = A value (or a list of values), of the proportion(s) of each quaternion in the cubic interpolation.
// Example(3D): Giving `u` as a Scalar
// a = QuatY(-135);
// b = QuatXYZ([-50,-50,120]);
// c = QuatXYZ([-50,-40,30]);
// d = QuatY(-45);
// color("red"){
// Qrot(b) right(80) cube([10,10,1]);
// Qrot(c) right(80) cube([10,10,1]);
// }
// for (u=[0:0.05:1])
// Qrot(Q_Squad(a, b, c, d, u))
// right(80) cube([10,10,1]);
// #sphere(r=80);
// Example(3D): Giving `u` as a Range
// a = QuatY(-135);
// b = QuatXYZ([-50,-50,120]);
// c = QuatXYZ([-50,-40,30]);
// d = QuatY(-45);
// for (q = Q_Squad(a, b, c, d, [0:0.05:1]))
// Qrot(q) right(80) cube([10,10,1]);
// #sphere(r=80);
function Q_Squad(q1,q2,q3,q4,u) =
assert(is_finite(u) || is_range(u) || is_vector(u) ,
"Invalid interpolation coefficient(s)" )
is_num(u)
? Q_Slerp( Q_Slerp(q1,q4,u), Q_Slerp(q2,q3,u), 2*u*(1-u))
: [for(ui=u) Q_Slerp( Q_Slerp(q1,q4,ui), Q_Slerp(q2,q3,ui), 2*ui*(1-ui) ) ];
// Function: Q_exp()
// Usage:
// q2 = Q_exp(q);
// Description:
// Returns the quaternion that is the exponential of the quaternion q in base e
// The returned quaternion is usually not normalized.
function Q_exp(q) =
assert( is_vector(q,4), "Input is not a valid quaternion")
let( nv = norm(_Qvec(q)) ) // q may be equal to zero here!
exp(_Qreal(q))*Quat(_Qvec(q),2*nv);
// Function: Q_ln()
// Usage:
// q2 = Q_ln(q);
// Description:
// Returns the quaternion that is the natural logarithm of the quaternion q.
// The returned quaternion is usually not normalized and may be zero.
function Q_ln(q) =
assert(Q_is_quat(q), "Input is not a valid quaternion")
let( nq = norm(q),
nv = norm(_Qvec(q)) )
approx(nv,0) ? _Qset([0,0,0] , ln(nq) ) :
_Qset(_Qvec(q)*atan2(nv,_Qreal(q))/nv, ln(nq));
// Function: Q_pow()
// Usage:
// q2 = Q_pow(q, r);
// Description:
// Returns the quaternion that is the power of the quaternion q to the real exponent r.
// The returned quaternion is normalized if `q` is normalized.
function Q_pow(q,r=1) =
assert( Q_is_quat(q) && is_finite(r),
"Invalid inputs")
let( theta = 2*atan2(norm(_Qvec(q)),_Qreal(q)) )
Quat(_Qvec(q), r*theta); // Q_exp(r*Q_ln(q));
// vim: expandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap

16
rounding.scad
View File

@@ -159,6 +159,7 @@ include <structs.scad>
// fwd(60) // Note how the different points are cut back by different amounts
// stroke(round_corners(zig,radius=1.5,closed=false),width=1);
// Example(FlatSpin): Rounding some random 3D paths
// $fn=36;
// list1= [
// [2.887360, 4.03497, 6.372090],
// [5.682210, 9.37103, 0.783548],
@@ -926,10 +927,10 @@ function os_profile(points, extra,check_valid, quality, offset_maxstep, offset)
//
// Example: Chamfered elliptical prism. If you stretch a chamfered cylinder the chamfer will be uneven.
// convex_offset_extrude(bottom = os_chamfer(height=-2), top=os_chamfer(height=1), height=7)
// xscale(4)circle(r=6);
// xscale(4)circle(r=6,$fn=64);
// Example: Elliptical prism with circular roundovers.
// convex_offset_extrude(bottom=os_circle(r=-2), top=os_circle(r=1), height=7,steps=10)
// xscale(4)circle(r=6);
// xscale(4)circle(r=6,$fn=64);
// Example: If you give a non-convex input you get a convex hull output
// right(50) linear_extrude(height=7) star(5,r=22,ir=13);
// convex_offset_extrude(bottom = os_chamfer(height=-2), top=os_chamfer(height=1), height=7)
@@ -1122,6 +1123,7 @@ function _remove_undefined_vals(list) =
// fwd(7)
// offset_stroke(arc, width=2, start=os_pointed(loc=2,dist=2),end=os_pointed(loc=.5,dist=-1));
// Example(2D): The os_round() end treatment adds roundovers to the end corners by specifying the `cut` parameter. In the first example, the cut parameter is the same at each corner. The bezier smoothness parameter `k` is given to allow a larger cut. In the second example, each corner is given a different roundover, including zero for no rounding at all. The red shows the same strokes without the roundover.
// $fn=36;
// arc = arc(points=[[1,1],[3,4],[6,3]],N=50);
// path = [[0,0],[6,2],[9,7],[8,10]];
// offset_stroke(path, width=2, rounded=false,start=os_round(angle=-20, cut=0.4,k=.9), end=os_round(angle=-35, cut=0.4,k=.9));
@@ -1133,9 +1135,9 @@ function _remove_undefined_vals(list) =
// Example(2D): Negative cut values produce a flaring end. Note how the absolute angle aligns the ends of the first example withi the axes. In the second example positive and negative cut values are combined. Note also that very different cuts are needed at the start end to produce a similar looking flare.
// arc = arc(points=[[1,1],[3,4],[6,3]],N=50);
// path = [[0,0],[6,2],[9,7],[8,10]];
// offset_stroke(path, width=2, rounded=false,start=os_round(cut=-1, abs_angle=90), end=os_round(cut=-0.5, abs_angle=0));
// offset_stroke(path, width=2, rounded=false,start=os_round(cut=-1, abs_angle=90), end=os_round(cut=-0.5, abs_angle=0),$fn=36);
// right(10)
// offset_stroke(arc, width=2, rounded=false, start=os_round(cut=[-.75,-.2], angle=-45), end=os_round(cut=[-.2,.2], angle=20));
// offset_stroke(arc, width=2, rounded=false, start=os_round(cut=[-.75,-.2], angle=-45), end=os_round(cut=[-.2,.2], angle=20),$fn=36);
// Example(2D): Setting the width to a vector allows generation of a set of parallel strokes
// path = [[0,0],[4,4],[8,4],[2,9],[10,10]];
// for(i=[0:.25:2])
@@ -1146,9 +1148,9 @@ function _remove_undefined_vals(list) =
// offset_stroke(path, rounded=true,width = [i,i+.08]);
// Example(2D): In this example a spurious triangle appears. This results from overly enthusiastic validity checking. Turning validity checking off fixes it in this case.
// path = [[0,0],[4,4],[8,4],[2,9],[10,10]];
// offset_stroke(path, check_valid=true,rounded=false,width = [1.4, 1.45]);
// offset_stroke(path, check_valid=true,rounded=false,width = [1.4, 1.5]);
// right(2)
// offset_stroke(path, check_valid=false,rounded=false,width = [1.4, 1.45]);
// offset_stroke(path, check_valid=false,rounded=false,width = [1.4, 1.5]);
// Example(2D): But in this case, disabling the validity check produces an invalid result.
// path = [[0,0],[4,4],[8,4],[2,9],[10,10]];
// offset_stroke(path, check_valid=true,rounded=false,width = [1.9, 2]);
@@ -1786,7 +1788,7 @@ function _circle_mask(r) =
// less than 180 degrees, and the path shouldn't have closely spaced points at concave points of high curvature because
// this will cause self-intersection in the mask polyhedron, resulting in CGAL failures.
// Arguments:
// r|radius = center radius of the cylindrical shell to cut a hole in
// r / radius = center radius of the cylindrical shell to cut a hole in
// thickness = thickness of cylindrical shell (may need to be slighly oversized)
// path = 2d path that defines the hole to cut
// Example: The mask as long pointed ends because this was the most efficient way to close off those ends.

2
scripts/docs_gen.py
View File

@@ -97,7 +97,7 @@ def image_compare(file1, file2):
sq = (value * (i % 256) ** 2 for i, value in enumerate(diff))
sum_squares = sum(sq)
rms = math.sqrt(sum_squares / float(img1.size[0] * img1.size[1]))
return rms<10
return rms<2
def image_resize(infile, outfile, newsize=(320,240)):

24
shapes2d.scad
View File

@@ -237,8 +237,8 @@ module stroke(
} else {
quatsums = Q_Cumulative([
for (i = idx(path2,end=-2)) let(
vec1 = i==0? UP : unit(path2[i]-path2[i-1]),
vec2 = unit(path2[i+1]-path2[i]),
vec1 = i==0? UP : unit(path2[i]-path2[i-1], UP),
vec2 = unit(path2[i+1]-path2[i], UP),
axis = vector_axis(vec1,vec2),
ang = vector_angle(vec1,vec2)
) Quat(axis,ang)
@@ -961,8 +961,8 @@ function regular_ngon(n=6, r, d, or, od, ir, id, side, rounding=0, realign=false
tipp = polar_to_xy(r-inset+rounding,a1),
pos = (p1+p2)/2
) each [
anchorpt(str("tip",i), tipp, unit(tipp), 0),
anchorpt(str("side",i), pos, unit(pos), 0),
anchorpt(str("tip",i), tipp, unit(tipp,BACK), 0),
anchorpt(str("side",i), pos, unit(pos,BACK), 0),
]
]
) reorient(anchor,spin, two_d=true, path=path, extent=false, p=path, anchors=anchors);
@@ -983,8 +983,8 @@ module regular_ngon(n=6, r, d, or, od, ir, id, side, rounding=0, realign=false,
tipp = polar_to_xy(r-inset+rounding,a1),
pos = (p1+p2)/2
) each [
anchorpt(str("tip",i), tipp, unit(tipp), 0),
anchorpt(str("side",i), pos, unit(pos), 0),
anchorpt(str("tip",i), tipp, unit(tipp,BACK), 0),
anchorpt(str("side",i), pos, unit(pos,BACK), 0),
]
];
attachable(anchor,spin, two_d=true, path=path, extent=false, anchors=anchors) {
@@ -1332,9 +1332,9 @@ function star(n, r, d, or, od, ir, id, step, realign=false, anchor=CENTER, spin=
p3 = polar_to_xy(r,a3),
pos = (p1+p3)/2
) each [
anchorpt(str("tip",i), p1, unit(p1), 0),
anchorpt(str("corner",i), p2, unit(p2), 0),
anchorpt(str("midpt",i), pos, unit(pos), 0),
anchorpt(str("tip",i), p1, unit(p1,BACK), 0),
anchorpt(str("corner",i), p2, unit(p2,BACK), 0),
anchorpt(str("midpt",i), pos, unit(pos,BACK), 0),
]
]
) reorient(anchor,spin, two_d=true, path=path, p=path, anchors=anchors);
@@ -1355,9 +1355,9 @@ module star(n, r, d, or, od, ir, id, step, realign=false, anchor=CENTER, spin=0)
p3 = polar_to_xy(r,a3),
pos = (p1+p3)/2
) each [
anchorpt(str("tip",i), p1, unit(p1), 0),
anchorpt(str("corner",i), p2, unit(p2), 0),
anchorpt(str("midpt",i), pos, unit(pos), 0),
anchorpt(str("tip",i), p1, unit(p1,BACK), 0),
anchorpt(str("corner",i), p2, unit(p2,BACK), 0),
anchorpt(str("midpt",i), pos, unit(pos,BACK), 0),
]
];
attachable(anchor,spin, two_d=true, path=path, anchors=anchors) {

4
tests/hull.scad
View File

@@ -1,5 +1,5 @@
include <BOSL2/std.scad>
include <BOSL2/hull.scad>
include <../std.scad>
include <../hull.scad>
testpoints_on_sphere = [ for(p =

4
tests/polyhedra.scad
View File

@@ -1,5 +1,5 @@
include<BOSL2/std.scad>
include<BOSL2/polyhedra.scad>
include<../std.scad>
include<../polyhedra.scad>
$fn=96;

2
tests/test_affine.scad
View File

@@ -1,4 +1,4 @@
include <BOSL2/std.scad>
include <../std.scad>
module test_ident() {

8
tests/test_arrays.scad
View File

@@ -1,5 +1,4 @@
include <BOSL2/std.scad>
include <../std.scad>
// List/Array Ops
@@ -16,6 +15,11 @@ module test_in_list() {
assert(in_list("bar", ["foo", "bar", "baz"]));
assert(!in_list("bee", ["foo", "bar", "baz"]));
assert(in_list("bar", [[2,"foo"], [4,"bar"], [3,"baz"]], idx=1));
assert(!in_list(undef, [3,4,5]));
assert(in_list(undef,[3,4,undef,5]));
assert(!in_list(3,[]));
assert(!in_list(3,[4,5,[3]]));
}
test_in_list();

2
tests/test_common.scad
View File

@@ -1,4 +1,4 @@
include <BOSL2/std.scad>
include <../std.scad>
module test_typeof() {

2
tests/test_coords.scad
View File

@@ -1,4 +1,4 @@
include <BOSL2/std.scad>
include <../std.scad>
module test_point2d() {

4
tests/test_cubetruss.scad
View File

@@ -1,5 +1,5 @@
include <BOSL2/std.scad>
include <BOSL2/cubetruss.scad>
include <../std.scad>
include <../cubetruss.scad>
module test_cubetruss_dist() {

2
tests/test_debug.scad
View File

@@ -1,4 +1,4 @@
include <BOSL2/std.scad>
include <../std.scad>
module test_standard_anchors() {

2
tests/test_edges.scad
View File

@@ -1,4 +1,4 @@
include <BOSL2/std.scad>
include <../std.scad>
module test_is_edge_array() {

2
tests/test_errors.scad
View File

@@ -1,4 +1,4 @@
include <BOSL2/std.scad>
include <../std.scad>
// Can't test echo output as yet. Include these for coverage calculations.

45
tests/test_geometry.scad
View File

@@ -1,4 +1,4 @@
include <BOSL2/std.scad>
include <../std.scad>
module test_point_on_segment2d() {
@@ -87,7 +87,7 @@ module test_line_normal() {
assert(line_normal([[0,0],[-10,0]]) == [0,-1]);
assert(line_normal([[0,0],[0,-10]]) == [1,0]);
assert(approx(line_normal([[0,0],[10,10]]), [-sqrt(2)/2,sqrt(2)/2]));
pts = [for (i=list_range(1000)) rands(-100,100,2,seed_value=4312)];
pts = [for (p=pair(rands(-100,100,1000,seed_value=4312))) p];
for (p = pair_wrap(pts)) {
p1 = p.x;
p2 = p.y;
@@ -351,18 +351,18 @@ module test_tri_functions() {
opp = p.y;
hyp = norm([opp,adj]);
ang = atan2(opp,adj);
assert(approx(hyp_opp_to_adj(hyp,opp),adj));
assert(approx(hyp_ang_to_adj(hyp,ang),adj));
assert(approx(opp_ang_to_adj(opp,ang),adj));
assert(approx(hyp_adj_to_opp(hyp,adj),opp));
assert(approx(hyp_ang_to_opp(hyp,ang),opp));
assert(approx(adj_ang_to_opp(adj,ang),opp));
assert(approx(adj_opp_to_hyp(adj,opp),hyp));
assert(approx(adj_ang_to_hyp(adj,ang),hyp));
assert(approx(opp_ang_to_hyp(opp,ang),hyp));
assert(approx(hyp_adj_to_ang(hyp,adj),ang));
assert(approx(hyp_opp_to_ang(hyp,opp),ang));
assert(approx(adj_opp_to_ang(adj,opp),ang));
assert_approx(hyp_opp_to_adj(hyp,opp), adj);
assert_approx(hyp_ang_to_adj(hyp,ang), adj);
assert_approx(opp_ang_to_adj(opp,ang), adj);
assert_approx(hyp_adj_to_opp(hyp,adj), opp);
assert_approx(hyp_ang_to_opp(hyp,ang), opp);
assert_approx(adj_ang_to_opp(adj,ang), opp);
assert_approx(adj_opp_to_hyp(adj,opp), hyp);
assert_approx(adj_ang_to_hyp(adj,ang), hyp);
assert_approx(opp_ang_to_hyp(opp,ang), hyp);
assert_approx(hyp_adj_to_ang(hyp,adj), ang);
assert_approx(hyp_opp_to_ang(hyp,opp), ang);
assert_approx(adj_opp_to_ang(adj,opp), ang);
}
}
test_tri_functions();
@@ -613,6 +613,23 @@ module test_pointlist_bounds() {
[23,57,-42]
];
assert(pointlist_bounds(pts) == [[-63,-32,-42], [84,97,42]]);
pts2d = [
[-53,12],
[-63,36],
[84,-5],
[63,42],
[23,-42]
];
assert(pointlist_bounds(pts2d) == [[-63,-42],[84,42]]);
pts5d = [
[-53,27,12,-53,12],
[-63,97,36,-63,36],
[84,-32,-5,84,-5],
[63,-24,42,63,42],
[23,57,-42,23,-42]
];
assert(pointlist_bounds(pts5d) == [[-63,-32,-42,-63,-42],[84,97,42,84,42]]);
assert(pointlist_bounds([[3,4,5,6]]), [[3,4,5,6],[3,4,5,6]]);
}
test_pointlist_bounds();

4
tests/test_linear_bearings.scad
View File

@@ -1,5 +1,5 @@
include <BOSL2/std.scad>
include <BOSL2/linear_bearings.scad>
include <../std.scad>
include <../linear_bearings.scad>
module test_get_lmXuu_bearing_diam() {

322
tests/test_math.scad
View File

@@ -1,4 +1,4 @@
include <BOSL2/std.scad>
include <../std.scad>
// Simple Calculations
@@ -76,6 +76,12 @@ module test_is_matrix() {
assert(is_matrix([[2,3,4],[5,6,7],[8,9,10]],square=false));
assert(is_matrix([[2,3],[5,6],[8,9]],m=3,n=2));
assert(is_matrix([[2,3,4],[5,6,7]],m=2,n=3));
assert(is_matrix([[2,3,4],[5,6,7]],2,3));
assert(is_matrix([[2,3,4],[5,6,7]],m=2));
assert(is_matrix([[2,3,4],[5,6,7]],2));
assert(is_matrix([[2,3,4],[5,6,7]],n=3));
assert(!is_matrix([[2,3,4],[5,6,7]],m=4));
assert(!is_matrix([[2,3,4],[5,6,7]],n=5));
assert(!is_matrix([[2,3,4],[5,6,7]],m=2,n=3,square=true));
assert(is_matrix([[2,3,4],[5,6,7],[8,9,10]],square=false));
assert(!is_matrix([[2,3],[5,6],[8,9]],m=2,n=3));
@@ -181,6 +187,9 @@ module test_sqr() {
assert_equal(sqr(2.5), 6.25);
assert_equal(sqr(3), 9);
assert_equal(sqr(16), 256);
assert_equal(sqr([2,3,4]), [4,9,16]);
assert_equal(sqr([[2,3,4],[3,5,7]]), [[4,9,16],[9,25,49]]);
assert_equal(sqr([]),[]);
}
test_sqr();
@@ -525,6 +534,9 @@ module test_any() {
assert_equal(any([1,5,true]), true);
assert_equal(any([[0,0], [0,0]]), false);
assert_equal(any([[0,0], [1,0]]), true);
assert_equal(any([[false,false],[[false,[false],[[[true]]]],false],[false,false]]), true);
assert_equal(any([[false,false],[[false,[false],[[[false]]]],false],[false,false]]), false);
assert_equal(any([]), false);
}
test_any();
@@ -536,6 +548,9 @@ module test_all() {
assert_equal(all([[0,0], [0,0]]), false);
assert_equal(all([[0,0], [1,0]]), false);
assert_equal(all([[1,1], [1,1]]), true);
assert_equal(all([[true,true],[[true,[true],[[[true]]]],true],[true,true]]), true);
assert_equal(all([[true,true],[[true,[true],[[[false]]]],true],[true,true]]), false);
assert_equal(all([]), true);
}
test_all();
@@ -554,6 +569,7 @@ test_count_true();
module test_factorial() {
assert_equal(factorial(0), 1);
assert_equal(factorial(1), 1);
assert_equal(factorial(2), 2);
assert_equal(factorial(3), 6);
@@ -562,6 +578,8 @@ module test_factorial() {
assert_equal(factorial(6), 720);
assert_equal(factorial(7), 5040);
assert_equal(factorial(8), 40320);
assert_equal(factorial(25,21), 303600);
assert_equal(factorial(25,25), 1);
}
test_factorial();
@@ -570,6 +588,11 @@ module test_gcd() {
assert_equal(gcd(15,25), 5);
assert_equal(gcd(15,27), 3);
assert_equal(gcd(270,405), 135);
assert_equal(gcd(39, 101),1);
assert_equal(gcd(15,-25), 5);
assert_equal(gcd(-15,25), 5);
assert_equal(gcd(5,0),5);
assert_equal(gcd(0,5),5);
}
test_gcd();
@@ -578,9 +601,306 @@ module test_lcm() {
assert_equal(lcm(15,25), 75);
assert_equal(lcm(15,27), 135);
assert_equal(lcm(270,405), 810);
assert_equal(lcm([3,5,15,25,35]),525);
}
test_lcm();
module test_C_times() {
assert_equal(C_times([4,5],[9,-4]), [56,29]);
assert_equal(C_times([-7,2],[24,3]), [-174, 27]);
}
test_C_times();
module test_C_div() {
assert_equal(C_div([56,29],[9,-4]), [4,5]);
assert_equal(C_div([-174,27],[-7,2]), [24,3]);
}
test_C_div();
module test_back_substitute(){
R = [[12,4,3,2],
[0,2,-4,2],
[0,0,4,5],
[0,0,0,15]];
assert_approx(back_substitute(R, [1,2,3,3]), [-0.675, 1.8, 0.5, 0.2]);
assert_approx(back_substitute(R, [6, 3, 3.5, 37], transpose=true), [0.5, 0.5, 1, 2]);
assert_approx(back_substitute(R, [[38,101],[-6,-16], [31, 71], [45, 105]]), [[1, 4],[2,3],[4,9],[3,7]]);
assert_approx(back_substitute(R, [[12,48],[8,22],[11,36],[71,164]],transpose=true), [[1, 4],[2,3],[4,9],[3,7]]);
assert_approx(back_substitute([[2]], [4]), [2]);
sing1 =[[0,4,3,2],
[0,3,-4,2],
[0,0,4,5],
[0,0,0,15]];
sing2 =[[12,4,3,2],
[0,0,-4,2],
[0,0,4,5],
[0,0,0,15]];
sing3 = [[12,4,3,2],
[0,2,-4,2],
[0,0,4,5],
[0,0,0,0]];
assert_approx(back_substitute(sing1, [1,2,3,4]), []);
assert_approx(back_substitute(sing2, [1,2,3,4]), []);
assert_approx(back_substitute(sing3, [1,2,3,4]), []);
}
test_back_substitute();
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(subindex(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_outer_product(){
assert_equal(outer_product([1,2,3],[4,5,6]), [[4,5,6],[8,10,12],[12,15,18]]);
assert_equal(outer_product([9],[7]), [[63]]);
}
test_outer_product();
module test_deriv(){
pent = [for(x=[0:70:359]) [cos(x), sin(x)]];
assert_approx(deriv(pent,closed=true),
[[-0.321393804843,0.556670399226],
[-0.883022221559,0.321393804843],
[-0.604022773555,-0.719846310393],
[0.469846310393,-0.813797681349],
[0.925416578398,0.163175911167],
[0.413175911167,0.492403876506]]);
assert_approx(deriv(pent,closed=true,h=2),
0.5*[[-0.321393804843,0.556670399226],
[-0.883022221559,0.321393804843],
[-0.604022773555,-0.719846310393],
[0.469846310393,-0.813797681349],
[0.925416578398,0.163175911167],
[0.413175911167,0.492403876506]]);
assert_approx(deriv(pent,closed=false),
[[-0.432937491789,1.55799143673],
[-0.883022221559,0.321393804843],
[-0.604022773555,-0.719846310393],
[0.469846310393,-0.813797681349],
[0.925416578398,0.163175911167],
[0.696902572292,1.45914323952]]);
spent = yscale(8,pent);
lens = path_segment_lengths(spent,closed=true);
assert_approx(deriv(spent, closed=true, h=lens),
[[-0.0381285841663,0.998065839726],
[-0.254979378104,0.0449763331253],
[-0.216850793938,-0.953089506601],
[0.123993253223,-0.982919228715],
[0.191478335034,0.0131898128456],
[0.0674850818111,0.996109041561]]);
assert_approx(deriv(spent, closed=false, h=select(lens,0,-2)),
[[-0.0871925973657,0.996191473044],
[-0.254979378104,0.0449763331253],
[-0.216850793938,-0.953089506601],
[0.123993253223,-0.982919228715],
[0.191478335034,0.0131898128456],
[0.124034734589,0.992277876714]]);
}
test_deriv();
module test_deriv2(){
oct = [for(x=[0:45:359]) [cos(x), sin(x)]];
assert_approx(deriv2(oct),
[[-0.828427124746,0.0719095841794],[-0.414213562373,-0.414213562373],[0,-0.585786437627],
[0.414213562373,-0.414213562373],[0.585786437627,0],[0.414213562373,0.414213562373],
[0,0.585786437627],[-0.636634192232,0.534938683021]]);
assert_approx(deriv2(oct,closed=false),
[[-0.828427124746,0.0719095841794],[-0.414213562373,-0.414213562373],[0,-0.585786437627],
[0.414213562373,-0.414213562373],[0.585786437627,0],[0.414213562373,0.414213562373],
[0,0.585786437627],[-0.636634192232,0.534938683021]]);
assert_approx(deriv2(oct,closed=true),
[[-0.585786437627,0],[-0.414213562373,-0.414213562373],[0,-0.585786437627],
[0.414213562373,-0.414213562373],[0.585786437627,0],[0.414213562373,0.414213562373],
[0,0.585786437627],[-0.414213562373,0.414213562373]]);
assert_approx(deriv2(oct,closed=false,h=2),
0.25*[[-0.828427124746,0.0719095841794],[-0.414213562373,-0.414213562373],[0,-0.585786437627],
[0.414213562373,-0.414213562373],[0.585786437627,0],[0.414213562373,0.414213562373],
[0,0.585786437627],[-0.636634192232,0.534938683021]]);
assert_approx(deriv2(oct,closed=true,h=2),
0.25* [[-0.585786437627,0],[-0.414213562373,-0.414213562373],[0,-0.585786437627],
[0.414213562373,-0.414213562373],[0.585786437627,0],[0.414213562373,0.414213562373],
[0,0.585786437627],[-0.414213562373,0.414213562373]]);
}
test_deriv2();
module test_deriv3(){
oct = [for(x=[0:45:359]) [cos(x), sin(x)]];
assert_approx(deriv3(oct),
[[0.414213562373,-0.686291501015],[0.414213562373,-0.343145750508],[0.414213562373,0],
[0.292893218813,0.292893218813],[0,0.414213562373],[-0.292893218813,0.292893218813],
[-0.535533905933,0.0502525316942],[-0.778174593052,-0.192388155425]]);
assert_approx(deriv3(oct,closed=false),
[[0.414213562373,-0.686291501015],[0.414213562373,-0.343145750508],[0.414213562373,0],
[0.292893218813,0.292893218813],[0,0.414213562373],[-0.292893218813,0.292893218813],
[-0.535533905933,0.0502525316942],[-0.778174593052,-0.192388155425]]);
assert_approx(deriv3(oct,closed=false,h=2),
[[0.414213562373,-0.686291501015],[0.414213562373,-0.343145750508],[0.414213562373,0],
[0.292893218813,0.292893218813],[0,0.414213562373],[-0.292893218813,0.292893218813],
[-0.535533905933,0.0502525316942],[-0.778174593052,-0.192388155425]]/8);
assert_approx(deriv3(oct,closed=true),
[[0,-0.414213562373],[0.292893218813,-0.292893218813],[0.414213562373,0],[0.292893218813,0.292893218813],
[0,0.414213562373],[-0.292893218813,0.292893218813],[-0.414213562373,0],[-0.292893218813,-0.292893218813]]);
assert_approx(deriv3(oct,closed=true,h=2),
[[0,-0.414213562373],[0.292893218813,-0.292893218813],[0.414213562373,0],[0.292893218813,0.292893218813],
[0,0.414213562373],[-0.292893218813,0.292893218813],[-0.414213562373,0],[-0.292893218813,-0.292893218813]]/8);
}
test_deriv3();
module test_polynomial(){
assert_equal(polynomial([],12),0);
assert_equal(polynomial([],[12,4]),[0,0]);
assert_equal(polynomial([1,2,3,4],3),58);
assert_equal(polynomial([1,2,3,4],[3,-1]),[47,-41]);
assert_equal(polynomial([0,0,2],4),2);
}
test_polynomial();
module test_poly_roots(){
// Fifth roots of unity
assert_approx(
poly_roots([1,0,0,0,0,-1]),
[[1,0],[0.309016994375,0.951056516295],[-0.809016994375,0.587785252292],
[-0.809016994375,-0.587785252292],[0.309016994375,-0.951056516295]]);
assert_approx(poly_roots(poly_mult([[1,-2,5],[12,-24,24],[-2, -12, -20],[1,-10,50]])),
[[1, 1], [5, 5], [1, 2], [-3, 1], [-3, -1], [1, -1], [1, -2], [5, -5]]);
assert_approx(poly_roots([.124,.231,.942, -.334]),
[[0.3242874219074053,0],[-1.093595323856930,2.666477428660098], [-1.093595323856930,-2.666477428660098]]);
}
test_poly_roots();
module test_real_roots(){
// Wilkinson polynomial is a nasty test:
assert_approx(
sort(real_roots(poly_mult([[1,-1],[1,-2],[1,-3],[1,-4],[1,-5],[1,-6],[1,-7],[1,-8],[1,-9],[1,-10]]))),
list_range(n=10,s=1));
assert_equal(real_roots([3]), []);
assert_equal(real_roots(poly_mult([[1,-2,5],[12,-24,24],[-2, -12, -20],[1,-10,50]])),[]);
assert_equal(real_roots(poly_mult([[1,-2,5],[12,-24,24],[-2, -12, -20],[1,-10,50],[1,0,0]])),[0,0]);
assert_approx(real_roots(poly_mult([[1,-2,5],[12,-24,24],[-2, -12, -20],[1,-10,50],[1,4]])),[-4]);
assert(approx(real_roots([1,-10,25]),[5,5],eps=5e-6));
assert_approx(real_roots([4,-3]), [0.75]);
assert_approx(real_roots([0,0,0,4,-3]), [0.75]);
}
test_real_roots();
// 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]]);
}
test_submatrix();
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);
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]]));
}
test_qr_factor();
module test_poly_mult(){
assert_equal(poly_mult([3,2,1],[4,5,6,7]),[12,23,32,38,20,7]);
assert_equal(poly_mult([3,2,1],[]),[]);
assert_equal(poly_mult([[1,2],[3,4],[5,6]]), [15,68,100,48]);
assert_equal(poly_mult([[1,2],[],[5,6]]), []);
assert_equal(poly_mult([[3,4,5],[0,0,0]]),[]);
}
test_poly_mult();
module test_poly_div(){
assert_equal(poly_div(poly_mult([4,3,3,2],[2,1,3]), [2,1,3]),[[4,3,3,2],[]]);
assert_equal(poly_div([1,2,3,4],[1,2,3,4,5]), [[], [1,2,3,4]]);
assert_equal(poly_div(poly_add(poly_mult([1,2,3,4],[2,0,2]), [1,1,2]), [1,2,3,4]), [[2,0,2],[1,1,2]]);
assert_equal(poly_div([1,2,3,4], [1,-3]), [[1,5,18],[58]]);
}
test_poly_div();
module test_poly_add(){
assert_equal(poly_add([2,3,4],[3,4,5,6]),[3,6,8,10]);
assert_equal(poly_add([1,2,3,4],[-1,-2,3,4]), [6,8]);
assert_equal(poly_add([1,2,3],-[1,2,3]),[]);
}
test_poly_add();
// vim: expandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap

2
tests/test_mutators.scad
View File

@@ -1,4 +1,4 @@
include <BOSL2/std.scad>
include <../std.scad>
module test_HSL() {

2
tests/test_primitives.scad
View File

@@ -1,4 +1,4 @@
include <BOSL2/std.scad>
include <../std.scad>
module test_square() {

120
tests/test_quaternions.scad
View File

@@ -1,5 +1,5 @@
include <BOSL2/std.scad>
include <BOSL2/strings.scad>
include <../std.scad>
include <../strings.scad>
function rec_cmp(a,b,eps=1e-9) =
@@ -8,6 +8,7 @@ function rec_cmp(a,b,eps=1e-9) =
is_list(a)? len(a)==len(b) && all([for (i=idx(a)) rec_cmp(a[i],b[i],eps=eps)]) :
a == b;
function Qstandard(q) = sign([for(qi=q) if( ! approx(qi,0)) qi,0 ][0])*q;
module verify_f(actual,expected) {
if (!rec_cmp(actual,expected)) {
@@ -22,6 +23,15 @@ module verify_f(actual,expected) {
}
module test_is_quat() {
verify_f(Q_is_quat([0]),false);
verify_f(Q_is_quat([0,0,0,0]),false);
verify_f(Q_is_quat([1,0,2,0]),true);
verify_f(Q_is_quat([1,0,2,0,0]),false);
}
test_is_quat();
module test_Quat() {
verify_f(Quat(UP,0),[0,0,0,1]);
verify_f(Quat(FWD,0),[0,0,0,1]);
@@ -92,6 +102,15 @@ module test_QuatXYZ() {
test_QuatXYZ();
module test_Q_From_to() {
verify_f(Q_Mul(Q_From_to([1,2,3], [4,5,2]),Q_From_to([4,5,2], [1,2,3])), Q_Ident());
verify_f(Q_Matrix4(Q_From_to([1,2,3], [4,5,2])), rot(from=[1,2,3],to=[4,5,2]));
verify_f(Qrot(Q_From_to([1,2,3], -[1,2,3]),[1,2,3]), -[1,2,3]);
verify_f(unit(Qrot(Q_From_to([1,2,3], [4,5,2]),[1,2,3])), unit([4,5,2]));
}
test_Q_From_to();
module test_Q_Ident() {
verify_f(Q_Ident(), [0,0,0,1]);
}
@@ -207,6 +226,16 @@ module test_Q_Conj() {
test_Q_Conj();
module test_Q_Inverse() {
verify_f(Q_Inverse([1,0,0,1]),[-1,0,0,1]/sqrt(2));
verify_f(Q_Inverse([0,1,1,0]),[0,-1,-1,0]/sqrt(2));
verify_f(Q_Inverse(QuatXYZ([23,45,67])),Q_Conj(QuatXYZ([23,45,67])));
verify_f(Q_Mul(Q_Inverse(QuatXYZ([23,45,67])),QuatXYZ([23,45,67])),Q_Ident());
}
test_Q_Inverse();
module test_Q_Norm() {
verify_f(Q_Norm([1,0,0,1]),1.414213562);
verify_f(Q_Norm([0,1,1,0]),1.414213562);
@@ -276,6 +305,10 @@ module test_Q_Angle() {
verify_f(Q_Angle(QuatY(-37)),37);
verify_f(Q_Angle(QuatZ(37)),37);
verify_f(Q_Angle(QuatZ(-37)),37);
verify_f(Q_Angle(QuatZ(-37),QuatZ(-37)), 0);
verify_f(Q_Angle(QuatZ( 37.123),QuatZ(-37.123)), 74.246);
verify_f(Q_Angle(QuatX( 37),QuatY(-37)), 51.86293283);
}
test_Q_Angle();
@@ -288,4 +321,87 @@ module test_Qrot() {
test_Qrot();
module test_Q_Rotation() {
verify_f(Qstandard(Q_Rotation(Q_Matrix3(Quat([12,34,56],33)))),Qstandard(Quat([12,34,56],33)));
verify_f(Q_Matrix3(Q_Rotation(Q_Matrix3(QuatXYZ([12,34,56])))),
Q_Matrix3(QuatXYZ([12,34,56])));
}
test_Q_Rotation();
module test_Q_Rotation_path() {
verify_f(Q_Rotation_path(QuatX(135), 5, QuatY(13.5))[0] , Q_Matrix4(QuatX(135)));
verify_f(Q_Rotation_path(QuatX(135), 11, QuatY(13.5))[11] , yrot(13.5));
verify_f(Q_Rotation_path(QuatX(135), 16, QuatY(13.5))[8] , Q_Rotation_path(QuatX(135), 8, QuatY(13.5))[4]);
verify_f(Q_Rotation_path(QuatX(135), 16, QuatY(13.5))[7] ,
Q_Rotation_path(QuatY(13.5),16, QuatX(135))[9]);
verify_f(Q_Rotation_path(QuatX(11), 5)[0] , xrot(11));
verify_f(Q_Rotation_path(QuatX(11), 5)[4] , xrot(55));
}
test_Q_Rotation_path();
module test_Q_Nlerp() {
verify_f(Q_Nlerp(QuatX(45),QuatY(30),0.0),QuatX(45));
verify_f(Q_Nlerp(QuatX(45),QuatY(30),0.5),[0.1967063121, 0.1330377423, 0, 0.9713946602]);
verify_f(Q_Rotation_path(QuatX(135), 16, QuatY(13.5))[8] , Q_Matrix4(Q_Nlerp(QuatX(135), QuatY(13.5),0.5)));
verify_f(Q_Nlerp(QuatX(45),QuatY(30),1.0),QuatY(30));
}
test_Q_Nlerp();
module test_Q_Squad() {
verify_f(Q_Squad(QuatX(45),QuatZ(30),QuatX(90),QuatY(30),0.0),QuatX(45));
verify_f(Q_Squad(QuatX(45),QuatZ(30),QuatX(90),QuatY(30),1.0),QuatY(30));
verify_f(Q_Squad(QuatX(0),QuatX(30),QuatX(90),QuatX(120),0.5),
Q_Slerp(QuatX(0),QuatX(120),0.5));
verify_f(Q_Squad(QuatY(0),QuatY(0),QuatX(120),QuatX(120),0.3),
Q_Slerp(QuatY(0),QuatX(120),0.3));
}
test_Q_Squad();
module test_Q_exp() {
verify_f(Q_exp(Q_Ident()), exp(1)*Q_Ident());
verify_f(Q_exp([0,0,0,33.7]), exp(33.7)*Q_Ident());
verify_f(Q_exp(Q_ln(Q_Ident())), Q_Ident());
verify_f(Q_exp(Q_ln([1,2,3,0])), [1,2,3,0]);
verify_f(Q_exp(Q_ln(QuatXYZ([31,27,34]))), QuatXYZ([31,27,34]));
let(q=QuatXYZ([12,23,34]))
verify_f(Q_exp(q+Q_Inverse(q)),Q_Mul(Q_exp(q),Q_exp(Q_Inverse(q))));
}
test_Q_exp();
module test_Q_ln() {
verify_f(Q_ln([1,2,3,0]), [24.0535117721, 48.1070235442, 72.1605353164, 1.31952866481]);
verify_f(Q_ln(Q_Ident()), [0,0,0,0]);
verify_f(Q_ln(5.5*Q_Ident()), [0,0,0,ln(5.5)]);
verify_f(Q_ln(Q_exp(QuatXYZ([13,37,43]))), QuatXYZ([13,37,43]));
verify_f(Q_ln(QuatXYZ([12,23,34]))+Q_ln(Q_Inverse(QuatXYZ([12,23,34]))), [0,0,0,0]);
}
test_Q_ln();
module test_Q_pow() {
q = Quat([1,2,3],77);
verify_f(Q_pow(q,1), q);
verify_f(Q_pow(q,0), Q_Ident());
verify_f(Q_pow(q,-1), Q_Inverse(q));
verify_f(Q_pow(q,2), Q_Mul(q,q));
verify_f(Q_pow(q,3), Q_Mul(q,Q_pow(q,2)));
verify_f(Q_Mul(Q_pow(q,0.456),Q_pow(q,0.544)), q);
verify_f(Q_Mul(Q_pow(q,0.335),Q_Mul(Q_pow(q,.552),Q_pow(q,.113))), q);
}
test_Q_pow();
// vim: expandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap

4
tests/test_queues.scad
View File

@@ -1,5 +1,5 @@
include <BOSL2/std.scad>
include <BOSL2/queues.scad>
include <../std.scad>
include <../queues.scad>
module test_queue_init() {

4
tests/test_shapes.scad
View File

@@ -1,5 +1,5 @@
include <BOSL2/std.scad>
include <BOSL2/hull.scad>
include <../std.scad>
include <../hull.scad>
module test_prismoid() {

2
tests/test_shapes2d.scad
View File

@@ -1,4 +1,4 @@
include <BOSL2/std.scad>
include <../std.scad>
module test_turtle() {

4
tests/test_skin.scad
View File

@@ -1,5 +1,5 @@
include <BOSL2/std.scad>
include <BOSL2/skin.scad>
include <../std.scad>
include <../skin.scad>
module test_skin() {

4
tests/test_stacks.scad
View File

@@ -1,5 +1,5 @@
include <BOSL2/std.scad>
include <BOSL2/stacks.scad>
include <../std.scad>
include <../stacks.scad>
module test_stack_init() {

4
tests/test_strings.scad
View File

@@ -1,5 +1,5 @@
include <BOSL2/std.scad>
include <BOSL2/strings.scad>
include <../std.scad>
include <../strings.scad>
module test_upcase() {

4
tests/test_structs.scad
View File

@@ -1,5 +1,5 @@
include <BOSL2/std.scad>
include <BOSL2/structs.scad>
include <../std.scad>
include <../structs.scad>
module test_struct_set() {

2
tests/test_transforms.scad
View File

@@ -1,4 +1,4 @@
include <BOSL2/std.scad>
include <../std.scad>
module test_translate() {

8
tests/test_vectors.scad
View File

@@ -1,4 +1,4 @@
include <BOSL2/std.scad>
include <../std.scad>
module test_is_vector() {
@@ -9,6 +9,10 @@ module test_is_vector() {
assert(is_vector(1) == false);
assert(is_vector("foo") == false);
assert(is_vector(true) == false);
assert(is_vector([0,0,0],zero=true) == true);
assert(is_vector([0,0,0],zero=false) == false);
assert(is_vector([0,1,0],zero=true) == false);
assert(is_vector([0,0,1],zero=false) == true);
}
test_is_vector();
@@ -56,7 +60,7 @@ module test_vabs() {
}
test_vabs();
include <BOSL2/strings.scad>
include <../strings.scad>
module test_vang() {
assert(vang([1,0])==0);
assert(vang([0,1])==90);

2
tests/test_version.scad
View File

@@ -1,4 +1,4 @@
include <BOSL2/std.scad>
include <../std.scad>
module test_bosl_version() {

4
tests/test_vnf.scad
View File

@@ -1,5 +1,5 @@
include <BOSL2/std.scad>
include <BOSL2/vnf.scad>
include <../std.scad>
include <../vnf.scad>
module test_is_vnf() {

66
vectors.scad
View File

@@ -19,6 +19,8 @@
// Arguments:
// v = The value to test to see if it is a vector.
// length = If given, make sure the vector is `length` items long.
// zero = If false, require that the length of the vector is not approximately zero. If true, require the length of the vector to be approx zero-length. Default: `undef` (don't check vector length.)
// eps = The minimum vector length that is considered non-zero. Default: `EPSILON` (`1e-9`)
// Example:
// is_vector(4); // Returns false
// is_vector([4,true,false]); // Returns false
@@ -28,8 +30,14 @@
// is_vector([3,4,5],3); // Returns true
// is_vector([3,4,5],4); // Returns true
// is_vector([]); // Returns false
function is_vector(v,length) =
is_list(v) && is_num(0*(v*v)) && (is_undef(length)||len(v)==length);
// is_vector([0,0,0],zero=true); // Returns true
// is_vector([0,0,0],zero=false); // Returns false
// is_vector([0,1,0],zero=true); // Returns false
// is_vector([0,0,1],zero=false); // Returns true
function is_vector(v,length,zero,eps=EPSILON) =
is_list(v) && is_num(0*(v*v))
&& (is_undef(length) || len(v)==length)
&& (is_undef(zero) || ((norm(v) >= eps) == !zero));
// Function: add_scalar()
@@ -105,18 +113,25 @@ function vceil(v) = [for (x=v) ceil(x)];
// Function: unit()
// Usage:
// unit(v, [error]);
// Description:
// Returns unit length normalized version of vector v.
// If passed a zero-length vector, returns the unchanged vector.
// Returns the unit length normalized version of vector v. If passed a zero-length vector,
// asserts an error unless `error` is given, in which case the value of `error` is returned.
// Arguments:
// v = The vector to normalize.
// error = If given, and input is a zero-length vector, this value is returned. Default: Assert error on zero-length vector.
// Examples:
// unit([10,0,0]); // Returns: [1,0,0]
// unit([0,10,0]); // Returns: [0,1,0]
// unit([0,0,10]); // Returns: [0,0,1]
// unit([0,-10,0]); // Returns: [0,-1,0]
// unit([0,0,0]); // Returns: [0,0,0]
function unit(v) = assert(is_vector(v),str(v)) norm(v)<=EPSILON? v : v/norm(v);
// unit([0,0,0],[1,2,3]); // Returns: [1,2,3]
// unit([0,0,0]); // Asserts an error.
function unit(v, error=[[["ASSERT"]]]) =
assert(is_vector(v), str("Expected a vector. Got: ",v))
norm(v)<EPSILON? (error==[[["ASSERT"]]]? assert(norm(v)>=EPSILON) : error) :
v/norm(v);
// Function: vector_angle()
@@ -181,22 +196,31 @@ function vector_angle(v1,v2,v3) =
// vector_axis([10,0,10], [0,0,0], [-10,10,0]); // Returns: [-0.57735, -0.57735, 0.57735]
// vector_axis([[10,0,10], [0,0,0], [-10,10,0]]); // Returns: [-0.57735, -0.57735, 0.57735]
function vector_axis(v1,v2=undef,v3=undef) =
(is_list(v1) && is_list(v1[0]) && is_undef(v2) && is_undef(v3))? (
assert(is_vector(v1.x))
assert(is_vector(v1.y))
len(v1)==3? assert(is_vector(v1.z)) vector_axis(v1.x, v1.y, v1.z) :
len(v1)==2? vector_axis(v1.x, v1.y) :
assert(false, "Bad arguments.")
) :
(is_vector(v1) && is_vector(v2) && is_vector(v3))? vector_axis(v1-v2, v3-v2) :
(is_vector(v1) && is_vector(v2) && is_undef(v3))? let(
is_vector(v3)
? assert(is_consistent([v3,v2,v1]), "Bad arguments.")
vector_axis(v1-v2, v3-v2)
:
assert( is_undef(v3), "Bad arguments.")
is_undef(v2)
? assert( is_list(v1), "Bad arguments.")
len(v1) == 2
? vector_axis(v1[0],v1[1])
: vector_axis(v1[0],v1[1],v1[2])
:
assert(
is_vector(v1,zero=false) &&
is_vector(v2,zero=false) &&
is_consistent([v1,v2]),
"Bad arguments."
)
let(
eps = 1e-6,
v1 = point3d(v1/norm(v1)),
v2 = point3d(v2/norm(v2)),
v3 = (norm(v1-v2) > eps && norm(v1+v2) > eps)? v2 :
(norm(vabs(v2)-UP) > eps)? UP :
RIGHT
) unit(cross(v1,v3)) : assert(false, "Bad arguments.");
w1 = point3d(v1/norm(v1)),
w2 = point3d(v2/norm(v2)),
w3 = (norm(w1-w2) > eps && norm(w1+w2) > eps) ? w2
: (norm(vabs(w2)-UP) > eps) ? UP
: RIGHT
) unit(cross(w1,w3));
// vim: expandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap

2
version.scad
View File

@@ -8,7 +8,7 @@
//////////////////////////////////////////////////////////////////////
BOSL_VERSION = [2,0,389];
BOSL_VERSION = [2,0,397];
// Section: BOSL Library Version Functions