From 1ced82f16c80751bc1d49bdb6bc5ff425488e41d Mon Sep 17 00:00:00 2001 From: RonaldoCMP Date: Wed, 26 Aug 2020 14:07:12 +0100 Subject: [PATCH 01/23] Correction of C_times validation --- arrays.scad | 127 ++++++++++++++++++++--- math.scad | 2 +- quaternions.scad | 170 +++++++++++++++++-------------- tests/test_all.scad | 28 ++++++ tests/test_quaternions.scad | 196 ++++++++++++++++++++++++++++-------- 5 files changed, 391 insertions(+), 132 deletions(-) create mode 100644 tests/test_all.scad diff --git a/arrays.scad b/arrays.scad index 835fd65..3ae7a88 100644 --- a/arrays.scad +++ b/arrays.scad @@ -708,36 +708,108 @@ function _sort_vectors4(arr) = y ] ) concat( _sort_vectors4(lesser), equal, _sort_vectors4(greater) ); +// sort a list of vectors +function _sort_vectors(arr, _i=0) = + len(arr)<=1 || _i>=len(arr[0]) ? arr : + let( + pivot = arr[floor(len(arr)/2)][_i], + lesser = [ for (entry=arr) if (entry[_i] < pivot ) entry ], + equal = [ for (entry=arr) if (entry[_i] == pivot ) entry ], + greater = [ for (entry=arr) if (entry[_i] > pivot ) entry ] + ) + concat( + _sort_vectors(lesser, _i ), + _sort_vectors(equal, _i+1 ), + _sort_vectors(greater, _i ) ); + +// given pairs of an index and a vector, return the list of indices of the list sorted by the vectors +function _sort_vectors_indexed(arr, _i=0) = + arr==[] ? [] : + len(arr)==1 || _i>=len(arr[0][1]) ? [for(ai=arr) ai[0]] : + let( + pivot = arr[floor(len(arr)/2)][1][_i], + lesser = [ for (entry=arr) if (entry[1][_i] < pivot ) entry ], + equal = [ for (entry=arr) if (entry[1][_i] == pivot ) entry ], + greater = [ for (entry=arr) if (entry[1][_i] > pivot ) entry ] + ) + concat( + _sort_vectors_indexed(lesser, _i ), + _sort_vectors_indexed(equal, _i+1 ), + _sort_vectors_indexed(greater, _i ) ); + // when idx==undef, returns the sorted array // otherwise, returns the indices of the sorted array function _sort_general(arr, idx=undef) = - (len(arr)<=1) ? arr : + len(arr)<=1 ? arr : is_undef(idx) - ? _sort_scalar(arr) + ? _simple_sort(arr) +// : _lexical_sort(arr) : let( arrind=[for(k=[0:len(arr)-1], ark=[arr[k]]) [ k, [for (i=idx) ark[i]] ] ] ) - _indexed_sort(arrind); - -// given a list of pairs, return the first element of each pair of the list sorted by the second element of the pair -// the sorting is done using compare_vals() -function _indexed_sort(arrind) = - arrind==[] ? [] : len(arrind)==1? [arrind[0][0]] : - let( pivot = arrind[floor(len(arrind)/2)][1] ) + _indexed_sort(arrind ); + +// sort simple lists with compare_vals() +function _simple_sort(arr) = + arr==[] || len(arr)==1 ? arr : let( - lesser = [ for (entry=arrind) if (compare_vals(entry[1], pivot) <0 ) entry ], - equal = [ for (entry=arrind) if (compare_vals(entry[1], pivot)==0 ) entry[0] ], - greater = [ for (entry=arrind) if (compare_vals(entry[1], pivot) >0 ) entry ] + pivot = arr[floor(len(arr)/2)], + lesser = [ for (entry=arr) if (compare_vals(entry, pivot) <0 ) entry ], + equal = [ for (entry=arr) if (compare_vals(entry, pivot)==0 ) entry ], + greater = [ for (entry=arr) if (compare_vals(entry, pivot) >0 ) entry ] ) - concat(_indexed_sort(lesser), equal, _indexed_sort(greater)); + concat( + _simple_sort(lesser), + equal, + _simple_sort(greater) + ); + + +// given a list of pairs, return the first element of each pair of the list sorted by the second element of the pair +// it uses compare_vals() +function _lexical_sort(arr, _i=0) = + arr==[] || len(arr)==1 || _i>=len(arr[0]) ? arr : + let( + pivot = arr[floor(len(arr)/2)][_i], + lesser = [ for (entry=arr) if (compare_vals(entry[_i], pivot) <0 ) entry ], + equal = [ for (entry=arr) if (compare_vals(entry[_i], pivot)==0 ) entry ], + greater = [ for (entry=arr) if (compare_vals(entry[_i], pivot) >0 ) entry ] + ) + concat( + _lexical_sort(lesser, _i ), + _lexical_sort(equal, _i+1 ), + _lexical_sort(greater, _i ) + ); + + +// given a list of pairs, return the first element of each pair of the list sorted by the second element of the pair +// it uses compare_vals() +function _indexed_sort(arr, _i=0) = + arr==[] ? [] : + len(arr)==1? [arr[0][0]] : + _i>=len(arr[0][1]) ? [for(ai=arr) ai[0]] : + let( + pivot = arr[floor(len(arr)/2)][1][_i], + lesser = [ for (entry=arr) if (compare_vals(entry[1][_i], pivot) <0 ) entry ], + equal = [ for (entry=arr) if (compare_vals(entry[1][_i], pivot)==0 ) entry ], + greater = [ for (entry=arr) if (compare_vals(entry[1][_i], pivot) >0 ) entry ] + ) + concat( + _indexed_sort(lesser, _i ), + _indexed_sort(equal, _i+1 ), + _indexed_sort(greater, _i ) ); // returns true for valid index specifications idx in the interval [imin, imax) // note that idx can't have any value greater or EQUAL to imax +// this allows imax=INF as a bound to numerical lists function _valid_idx(idx,imin,imax) = is_undef(idx) || ( is_finite(idx) && idx>=imin && idx< imax ) || ( is_list(idx) && min(idx)>=imin && max(idx)< imax ) - || ( valid_range(idx) && idx[0]>=imin && idx[2]< imax ); + || ( ! is_list(idx) // it implicitly is a range + && (idx[1]>0 && idx[0]>=imin && idx[2]< imax) + || + (idx[0]=imin) ); // Function: sort() @@ -758,7 +830,7 @@ function _valid_idx(idx,imin,imax) = function sort(list, idx=undef) = !is_list(list) || len(list)<=1 ? list : is_def(idx) - ? assert( _valid_idx(idx,0,len(list)) , "Invalid indices.") + ? assert( _valid_idx(idx,0,len(list[0])) , "Invalid indices or out of range.") let( sarr = _sort_general(list,idx) ) [for(i=[0:len(sarr)-1]) list[sarr[i]] ] : let(size = array_dim(list)) @@ -773,6 +845,31 @@ function sort(list, idx=undef) = ) : _sort_general(list); +function sort(list, idx=undef) = + !is_list(list) || len(list)<=1 ? list : + is_vector(list) + ? assert( _valid_idx(idx,0,len(list[0])) , str("Invalid indices or out of range. ",list)) + is_def(idx) + ? sort_vector_indexed([for(i=[0:len(list)-1]) [i, [list(i)] ] ]) + : sort_scalar(list) + : is_matrix(list) + ? list==[] || list[0]==[] ? list : + assert( _valid_idx(idx,0,len(list[0])) , "Invalid indices or out of range.") + is_def(idx) + ? sort_vector_indexed([for(i=[0:len(list)-1], li=[list[i]]) [i, [for(ind=idx) li(ind)] ] ]) + : sort_vector(list) + : list==[] || list[0]==[] ? list : + let( llen = [for(li=list) !is_list(li) || is_string(li) ? 0: len(li)], + m = min(llen), + M = max(llen) + ) + M==0 ? _simple_sort(list) : + assert( m>0 && _valid_idx(idx,m-1,M) , "Invalid indices or out of range.") + is_def(idx) + ? _sort_general(list,idx) + : let( ils = _sort_general(list, [m:M]) ) + [for(i=[0:len(list)-1]) list[ils[i]] ]; + // Function: sortidx() // Description: diff --git a/math.scad b/math.scad index 888d048..e5c4625 100644 --- a/math.scad +++ b/math.scad @@ -1175,7 +1175,7 @@ function deriv3(data, h=1, closed=false) = // Description: // Multiplies two complex numbers represented by 2D vectors. function C_times(z1,z2) = - assert( is_vector(z1+z2,2), "Complex numbers should be represented by 2D vectors." ) + assert( is_matrix([z1,z2],2,2), "Complex numbers should be represented by 2D vectors" ) [ z1.x*z2.x - z1.y*z2.y, z1.x*z2.y + z1.y*z2.x ]; // Function: C_div() diff --git a/quaternions.scad b/quaternions.scad index e21801c..44b0e29 100644 --- a/quaternions.scad +++ b/quaternions.scad @@ -27,7 +27,7 @@ function _Qreal(q) = q[3]; function _Qset(v,r) = concat( v, r ); // normalizes without checking -function _Qnorm(q) = q/norm(q); +function _Qunit(q) = q/norm(q); // Function: Q_is_quat() @@ -36,7 +36,7 @@ function _Qnorm(q) = q/norm(q); // 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 Q_is_quat(q) = is_vector(q,4) ;//&& ! approx(norm(q),0) ; // Function: Quat() @@ -101,7 +101,7 @@ function QuatXYZ(a=[0,0,0]) = qy = QuatY(a[1]), qz = QuatZ(a[2]) ) - Q_Mul(qz, Q_Mul(qy, qx)); + _Qmul(qz, _Qmul(qy, qx)); // Function: Q_From_to() @@ -202,32 +202,36 @@ function Q_Sub(a, 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)") + _Qmul(a,b); + +function _Qmul(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, + a[3]*b[3] - a.x*b.x - a.y*b.y - a.z*b.z ]; - +// [ [a[3], -a.z, a.y, a.x], +// [ a.z, a[3], -a.x, a.y], +// [-a.y, a.x, a[3], a.z], +// [-a.x, -a.y, -a.z, a[3]] ]*[b.x,b.y,b.z,b[3]] + // Function: Q_Cumulative() // Usage: -// Q_Cumulative(v); +// Q_Cumulative(ql); // Description: // 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( - v, _i+1, - concat( - _acc, - [_i==0 ? v[_i] : Q_Mul(v[_i], select(_acc,-1))] - ) - ); +function Q_Cumulative(ql) = + assert( is_matrix(ql,undef,4) && len(ql)>0, "Invalid list of quaternions." ) + [for( i = 0, q = ql[0]; + i<=len(ql); + i = i+1, q = (i==len(ql))? 0: _Qmul(q,ql[i]) ) + q ]; // Function: Q_Dot() @@ -252,7 +256,7 @@ function Q_Neg(q) = // Q_Conj(q) // Description: Returns the conjugate of quaternion `q`. function Q_Conj(q) = - assert( Q_is_quat(q), "Invalid quaternion" ) + assert( Q_is_quat(q), str("Invalid quaternion",q) ) [-q.x, -q.y, -q.z, q[3]]; @@ -262,7 +266,7 @@ function Q_Conj(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) ) +// let(q = q/norm(q) ) [-q.x, -q.y, -q.z, q[3]]; @@ -283,7 +287,9 @@ function Q_Norm(q) = // Description: Normalizes quaternion `q`, so that norm([W,X,Y,Z]) == 1. function Q_Normalize(q) = assert( Q_is_quat(q) , "Invalid quaternion" ) - q/norm(q); + approx(_Qvec(q), [0,0,0]) + ? Q_Ident() + : q/norm(q); // Function: Q_Dist() @@ -318,31 +324,32 @@ function Q_Dist(q1, q2) = // Qrot(q) right(80) cube([10,10,1]); // #sphere(r=80); 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)); + 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), + q3 = dot>0.9995? q2: _Qunit(q2 - dot*q1) + ) + ! is_num(u) + ? [for (uu=u) _Qslerp(q1, q3, uu, dot)] + : _Qslerp(q1, q3, u, dot); +function _Qslerp(q1, q2, u, dot) = + dot>0.9995 + ? _Qunit(q1 + u*(q2-q1)) + : let( theta = u*acos(dot) ) + _Qunit(q1*cos(theta) + q2*sin(theta)); + -// Function: Q_Matrix3() +// Function: Q_to_matrix3() // Usage: -// Q_Matrix3(q); +// Q_to_matrix3(q); // Description: // Returns the 3x3 rotation matrix for the given normalized quaternion q. -function Q_Matrix3(q) = +function Q_to_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]], @@ -351,12 +358,12 @@ function Q_Matrix3(q) = ]; -// Function: Q_Matrix4() +// Function: Q_to_matrix4() // Usage: -// Q_Matrix4(q); +// Q_to_matrix4(q); // Description: // Returns the 4x4 rotation matrix for the given normalized quaternion q. -function Q_Matrix4(q) = +function Q_to_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], @@ -366,6 +373,35 @@ function Q_Matrix4(q) = ]; +// Function: Q_from_matrix() +// Usage: +// Q_from_matrix(M) +// Description: +// Returns a normalized quaternion corresponding to the rotation matrix M. +// M may be a 3x3 rotation matrix or a homogeneous 4x4 rotation matrix. +// The last row and last column of M are ignored for 4x4 matrices. +// It doesn't check whether M is in fact a rotation matrix. +// If M is not a rotation, the returned quaternion is unpredictable. +// +// based on https://en.wikipedia.org/wiki/Rotation_matrix +// +function Q_from_matrix(M) = + assert( is_matrix(M) && (len(M)==3 || len(M)==4) && (len(M[0])==3 || len(M[0])==4), + "The matrix should be 3x3 or 4x4") + let( tr = M[0][0]+M[1][1]+M[2][2] ) // M trace + tr>0 + ? let( r = sqrt(1+tr), s = -1/r/2 ) + _Qunit( _Qset([ M[1][2]-M[2][1], M[2][0]-M[0][2], M[0][1]-M[1][0] ]*s, r/2 ) ) + : let( + i = max_index([ M[0][0], M[1][1], M[2][2] ]), + r = sqrt(1 + 2*M[i][i] -M[0][0] -M[1][1] -M[2][2] ), + s = 1/r/2 + ) + i==0 ? _Qunit( _Qset( [ r/2, s*(M[1][0]+M[0][1]), s*(M[0][2]+M[2][0]) ], s*(M[2][1]-M[1][2])) ): + i==1 ? _Qunit( _Qset( [ s*(M[1][0]+M[0][1]), r/2, s*(M[2][1]+M[1][2]) ], s*(M[0][2]-M[2][0])) ) + : _Qunit( _Qset( [ s*(M[2][0]+M[0][2]), s*(M[1][2]+M[2][1]), r/2 ], s*(M[1][0]-M[0][1])) ) ; + + // Function: Q_Axis() // Usage: // Q_Axis(q) @@ -418,15 +454,15 @@ function Q_Angle(q1,q2) = // q = QuatXYZ([45,35,10]); // pts = Qrot(q, p=[[2,3,4], [4,5,6], [9,2,3]]); module Qrot(q) { - multmatrix(Q_Matrix4(q)) { + multmatrix(Q_to_matrix4(q)) { children(); } } function Qrot(q,p) = - is_undef(p)? Q_Matrix4(q) : + is_undef(p)? Q_to_matrix4(q) : is_vector(p)? Qrot(q,[p])[0] : - apply(Q_Matrix4(q), p); + apply(Q_to_matrix4(q), p); // Module: Qrot_copies() @@ -446,29 +482,6 @@ 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); @@ -518,11 +531,11 @@ function Q_Rotation_path(q1, n=1, q2) = 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] + ? [for( i=0, dR=Q_to_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]; + dq = Q_pow( _Qmul( q2, Q_Inverse(q1) ), 1/n ), + dR = Q_to_matrix4(dq) ) + [for( i=0, R=Q_to_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)) @@ -565,8 +578,8 @@ function Q_Nlerp(q1,q2,u) = 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 ) ]; + ? _Qunit((1-u)*q1 + u*q2 ) + : [for (ui=u) _Qunit((1-ui)*q1 + ui*q2 ) ]; // Function: Q_Squad() @@ -616,6 +629,17 @@ function Q_Squad(q1,q2,q3,q4,u) = : [for(ui=u) Q_Slerp( Q_Slerp(q1,q4,ui), Q_Slerp(q2,q3,ui), 2*ui*(1-ui) ) ]; +function Q_Scubic(q1,q2,q3,q4,u) = + assert(is_finite(u) || is_range(u) || is_vector(u) , + "Invalid interpolation coefficient(s)" ) + is_num(u) + ? let( q12 = Q_Slerp(q1,q2,u), + q23 = Q_Slerp(q2,q3,u), + q34 = Q_Slerp(q3,q4,u) ) + Q_Slerp(Q_Slerp(q12,q23,u),Q_Slerp(q23,q34,u),u) + : [for(ui=u) Q_Scubic( q1,q2,q3,q4,ui) ]; + + // Function: Q_exp() // Usage: // q2 = Q_exp(q); diff --git a/tests/test_all.scad b/tests/test_all.scad new file mode 100644 index 0000000..fe6512b --- /dev/null +++ b/tests/test_all.scad @@ -0,0 +1,28 @@ +//include +include +include +include +include +include +include +include +include +include +include +include +include +include +//include +include +include +include +include +include +include +include +include +include +include +include + + diff --git a/tests/test_quaternions.scad b/tests/test_quaternions.scad index b1e5130..af19408 100644 --- a/tests/test_quaternions.scad +++ b/tests/test_quaternions.scad @@ -1,6 +1,80 @@ include <../std.scad> include <../strings.scad> +function plane_fit(points) = + assert( is_matrix(points,undef,3) , "Improper point list." ) + len(points)< 2 ? [] : + len(points)==3 ? plane3pt(points[0],points[1],points[2]) : + let( + A = [for(pi=points) concat(pi,-1) ], + qr = qr_factor(A), + R = select(qr[1],0,3), +x=[for(ri=R) echo(R=ri)0], +//y=[for(qi=qr[0]) echo(Q=qi)0], + s = back_substitute + ( R, + sum([for(qi=qr[0])[for(j=[0:3]) qi[j] ] ]) + ) +//, s0 = echo(ls=linear_solve(A,[for(p=points) 0])) +,w=s==[]? echo("s == []")0:echo(s=s)0 + ) + s==[] + ? // points are collinear + [] + : // plane through the origin? + approx(norm([s.x, s.y, s.z]), 0) + ? let( + k = max_index([for(i=[0:2]) abs(s[i]) ]), + A = [[1,0,0], for(pi=points) pi ], + b = [1, for(i=[1:len(points)]) 0], + s = linear_solve(A,b) +, y=echo(s2=s) + ) + s==[]? []: + concat(s, 0)/norm(s) + : s/norm([s.x, s.y, s.z]) ; + +function plane_fit(points) = + assert( is_matrix(points,undef,3) , "Improper point list." ) + len(points)< 2 ? [] : + len(points)==3 ? plane3pt(points[0],points[1],points[2]) : + let( + A = [for(pi=points) concat(pi,-1) ], + B = [ for(pi=points) [0,1] ], + sB = transpose(linear_solve(A,B)), + s = sB==[] ? [] + : norm(sB[1])N/2? 10*[1,1,1]: rands(-1,1,3) ]; +pf = plane_fit(pts); +echo(pf); +//pm = concat(sum(pts)/len(pts), -1); +//echo(pmfit=pm*pf); function rec_cmp(a,b,eps=1e-9) = typeof(a)!=typeof(b)? false : @@ -8,7 +82,18 @@ 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; +function standardize(v) = + v==[] + ? [] + : sign(first_nonzero(v))*v; + +function first_nonzero(v) = + v==[] ? 0 + : is_num(v) ? v + : [for(vi=v) if(!is_list(vi) && ! approx(vi,0)) vi + else if(is_list(vi)) first_nonzero(vi), 0 ][0]; + +module assert_std(vc,ve) { assert_approx(standardize(vc),standardize(ve)); } module verify_f(actual,expected) { if (!rec_cmp(actual,expected)) { @@ -25,7 +110,7 @@ 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([0,0,0,0]),true); verify_f(Q_is_quat([1,0,2,0]),true); verify_f(Q_is_quat([1,0,2,0,0]),false); } @@ -104,7 +189,7 @@ 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(Q_to_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])); } @@ -200,7 +285,7 @@ test_Q_Mul(); module test_Q_Cumulative() { verify_f(Q_Cumulative([QuatZ(30),QuatX(57),QuatY(18)]),[[0, 0, 0.2588190451, 0.9659258263], [0.4608999698, -0.1234977747, 0.2274546059, 0.8488721457], [0.4908072659, 0.01081554785, 0.1525536221, 0.8577404293]]); } -test_Q_Cumulative(); +*test_Q_Cumulative(); module test_Q_Dot() { @@ -219,18 +304,18 @@ test_Q_Neg(); module test_Q_Conj() { + ang = rands(0,360,3); verify_f(Q_Conj([1,0,0,1]),[-1,0,0,1]); verify_f(Q_Conj([0,1,1,0]),[0,-1,-1,0]); verify_f(Q_Conj(QuatXYZ([23,45,67])),[0.0533818345, -0.4143703268, -0.4360652669, 0.7970537592]); + verify_f(Q_Mul(Q_Conj(QuatXYZ(ang)),QuatXYZ(ang)),Q_Ident()); } 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_Inverse([1,0,0,1]),[-1,0,0,1]); + verify_f(Q_Inverse([0,1,1,0]),[0,-1,-1,0]); verify_f(Q_Mul(Q_Inverse(QuatXYZ([23,45,67])),QuatXYZ([23,45,67])),Q_Ident()); } test_Q_Inverse(); @@ -247,7 +332,7 @@ test_Q_Norm(); module test_Q_Normalize() { verify_f(Q_Normalize([1,0,0,1]),[0.7071067812, 0, 0, 0.7071067812]); verify_f(Q_Normalize([0,1,1,0]),[0, 0.7071067812, 0.7071067812, 0]); - verify_f(Q_Normalize(QuatXYZ([23,45,67])),[-0.0533818345, 0.4143703268, 0.4360652669, 0.7970537592]); +// verify_f(Q_Normalize(QuatXYZ([23,45,67])),[-0.0533818345, 0.4143703268, 0.4360652669, 0.7970537592]); } test_Q_Normalize(); @@ -260,28 +345,47 @@ test_Q_Dist(); module test_Q_Slerp() { - verify_f(Q_Slerp(QuatX(45),QuatY(30),0.0),QuatX(45)); - verify_f(Q_Slerp(QuatX(45),QuatY(30),0.5),[0.1967063121, 0.1330377423, 0, 0.9713946602]); - verify_f(Q_Slerp(QuatX(45),QuatY(30),1.0),QuatY(30)); + u = rands(0,1,1)[0]; + ul = rands(0,1,5); + ul2 = [1,1,1,1,1]-ul; + a1 = rands(0,360,1)[0]; + a2 = rands(0,360,1)[0]; + a3 = rands(0,360,1)[0]; + verify_f(standardize(Q_Slerp(QuatX(a1),QuatY(a2),0.0)), standardize(QuatX(a1))); + verify_f(standardize(Q_Slerp(QuatX(45),QuatY(30),0.5)), + [0.1967063121, 0.1330377423, 0, 0.9713946602]); + verify_f(standardize(Q_Slerp(QuatX(a1),QuatY(a2),1.0)), standardize(QuatY(a2))); + verify_f(standardize(Q_Slerp(QuatXYZ([a1,a2,0]),-QuatY(a2),u)), + standardize(Q_Slerp(QuatXYZ([a1,a2,0]), QuatY(a2),u))); + verify_f(standardize(Q_Slerp(QuatX(a1),QuatX(a1),u)), standardize(QuatX(a1))); + verify_f(standardize(Q_Slerp(QuatXYZ([a1,a2,0]),QuatXYZ([a2,0,a1]),u)), + standardize(Q_Slerp(QuatXYZ([a2,0,a1]),QuatXYZ([a1,a2,0]),1-u))); + verify_f(standardize(Q_Slerp(QuatXYZ([a1,a2,0]),QuatXYZ([a2,0,a1]),ul)), + standardize(Q_Slerp(QuatXYZ([a2,0,a1]),QuatXYZ([a1,a2,0]),ul2))); } test_Q_Slerp(); -module test_Q_Matrix3() { - verify_f(Q_Matrix3(QuatZ(37)),rot(37,planar=true)); - verify_f(Q_Matrix3(QuatZ(-49)),rot(-49,planar=true)); +module test_Q_to_matrix3() { + rotZ_37 = rot(37); + rotZ_37_3 = [for(i=[0:2]) [for(j=[0:2]) rotZ_37[i][j] ] ]; + angs = [12,-49,40]; + rot4 = rot(angs); + rot3 = [for(i=[0:2]) [for(j=[0:2]) rot4[i][j] ] ]; + verify_f(Q_to_matrix3(QuatZ(37)),rotZ_37_3); + verify_f(Q_to_matrix3(QuatXYZ(angs)),rot3); } -test_Q_Matrix3(); +test_Q_to_matrix3(); -module test_Q_Matrix4() { - verify_f(Q_Matrix4(QuatZ(37)),rot(37)); - verify_f(Q_Matrix4(QuatZ(-49)),rot(-49)); - verify_f(Q_Matrix4(QuatX(37)),rot([37,0,0])); - verify_f(Q_Matrix4(QuatY(37)),rot([0,37,0])); - verify_f(Q_Matrix4(QuatXYZ([12,34,56])),rot([12,34,56])); +module test_Q_to_matrix4() { + verify_f(Q_to_matrix4(QuatZ(37)),rot(37)); + verify_f(Q_to_matrix4(QuatZ(-49)),rot(-49)); + verify_f(Q_to_matrix4(QuatX(37)),rot([37,0,0])); + verify_f(Q_to_matrix4(QuatY(37)),rot([0,37,0])); + verify_f(Q_to_matrix4(QuatXYZ([12,34,56])),rot([12,34,56])); } -test_Q_Matrix4(); +test_Q_to_matrix4(); module test_Q_Axis() { @@ -321,23 +425,23 @@ 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]))); +module test_Q_from_matrix() { + verify_f(standardize(Q_from_matrix(Q_to_matrix3(Quat([12,34,56],33)))),standardize(Quat([12,34,56],33))); + verify_f(Q_to_matrix3(Q_from_matrix(Q_to_matrix3(QuatXYZ([12,34,56])))), + Q_to_matrix3(QuatXYZ([12,34,56]))); } -test_Q_Rotation(); +test_Q_from_matrix(); 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), 5, QuatY(13.5))[0] , Q_to_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)[3] , xrot(11+(55-11)*3/4)); verify_f(Q_Rotation_path(QuatX(11), 5)[4] , xrot(55)); } @@ -347,48 +451,54 @@ 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_Rotation_path(QuatX(135), 16, QuatY(13.5))[8] , + Q_to_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)); + u = rands(0,1,5); + su = [1,1,1,1,1]-u; + verify_f(Q_Squad(QuatX(45),QuatZ(30),QuatX(32),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), + 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)); + verify_f( Q_Squad(QuatY(10),QuatZ(20),QuatX(32),QuatXYZ([210,120,30]),u[0] ), + Q_Squad(QuatXYZ([210,120,30]),QuatX(32),QuatZ(20),QuatY(10),1-u[0] ) ); + verify_f( Q_Squad(QuatY(10),QuatZ(20),QuatX(32),QuatXYZ([210,120,30]),u ), + Q_Squad(QuatXYZ([210,120,30]),QuatX(32),QuatZ(20),QuatY(10),su ) ); } test_Q_Squad(); module test_Q_exp() { + q=QuatXYZ(rands(0,360,3)); 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)))); + verify_f(Q_exp(Q_ln([1,2,3,4])), [1,2,3,4]); + verify_f(Q_exp(Q_ln(q)), q); + 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() { + q = QuatXYZ(rands(0,360,3)); 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]); + verify_f(Q_ln(Q_exp(q)), q); + verify_f(Q_ln(q)+Q_ln(Q_Conj(q)), [0,0,0,0]); } test_Q_ln(); module test_Q_pow() { - q = Quat([1,2,3],77); + q = QuatXYZ(rands(0,360,3)); verify_f(Q_pow(q,1), q); verify_f(Q_pow(q,0), Q_Ident()); verify_f(Q_pow(q,-1), Q_Inverse(q)); From 4a366967f96c29c61116125b0803a1dd405e1f91 Mon Sep 17 00:00:00 2001 From: RonaldoCMP Date: Wed, 26 Aug 2020 14:07:53 +0100 Subject: [PATCH 02/23] Revert "Correction of C_times validation" This reverts commit 1ced82f16c80751bc1d49bdb6bc5ff425488e41d. --- arrays.scad | 125 +++-------------------- math.scad | 2 +- quaternions.scad | 170 ++++++++++++++----------------- tests/test_all.scad | 28 ------ tests/test_quaternions.scad | 196 ++++++++---------------------------- 5 files changed, 131 insertions(+), 390 deletions(-) delete mode 100644 tests/test_all.scad diff --git a/arrays.scad b/arrays.scad index 3ae7a88..835fd65 100644 --- a/arrays.scad +++ b/arrays.scad @@ -708,108 +708,36 @@ function _sort_vectors4(arr) = y ] ) concat( _sort_vectors4(lesser), equal, _sort_vectors4(greater) ); -// sort a list of vectors -function _sort_vectors(arr, _i=0) = - len(arr)<=1 || _i>=len(arr[0]) ? arr : - let( - pivot = arr[floor(len(arr)/2)][_i], - lesser = [ for (entry=arr) if (entry[_i] < pivot ) entry ], - equal = [ for (entry=arr) if (entry[_i] == pivot ) entry ], - greater = [ for (entry=arr) if (entry[_i] > pivot ) entry ] - ) - concat( - _sort_vectors(lesser, _i ), - _sort_vectors(equal, _i+1 ), - _sort_vectors(greater, _i ) ); - -// given pairs of an index and a vector, return the list of indices of the list sorted by the vectors -function _sort_vectors_indexed(arr, _i=0) = - arr==[] ? [] : - len(arr)==1 || _i>=len(arr[0][1]) ? [for(ai=arr) ai[0]] : - let( - pivot = arr[floor(len(arr)/2)][1][_i], - lesser = [ for (entry=arr) if (entry[1][_i] < pivot ) entry ], - equal = [ for (entry=arr) if (entry[1][_i] == pivot ) entry ], - greater = [ for (entry=arr) if (entry[1][_i] > pivot ) entry ] - ) - concat( - _sort_vectors_indexed(lesser, _i ), - _sort_vectors_indexed(equal, _i+1 ), - _sort_vectors_indexed(greater, _i ) ); - // when idx==undef, returns the sorted array // otherwise, returns the indices of the sorted array function _sort_general(arr, idx=undef) = - len(arr)<=1 ? arr : + (len(arr)<=1) ? arr : is_undef(idx) - ? _simple_sort(arr) -// : _lexical_sort(arr) + ? _sort_scalar(arr) : let( arrind=[for(k=[0:len(arr)-1], ark=[arr[k]]) [ k, [for (i=idx) ark[i]] ] ] ) - _indexed_sort(arrind ); - -// sort simple lists with compare_vals() -function _simple_sort(arr) = - arr==[] || len(arr)==1 ? arr : - let( - pivot = arr[floor(len(arr)/2)], - lesser = [ for (entry=arr) if (compare_vals(entry, pivot) <0 ) entry ], - equal = [ for (entry=arr) if (compare_vals(entry, pivot)==0 ) entry ], - greater = [ for (entry=arr) if (compare_vals(entry, pivot) >0 ) entry ] - ) - concat( - _simple_sort(lesser), - equal, - _simple_sort(greater) - ); - - + _indexed_sort(arrind); + // given a list of pairs, return the first element of each pair of the list sorted by the second element of the pair -// it uses compare_vals() -function _lexical_sort(arr, _i=0) = - arr==[] || len(arr)==1 || _i>=len(arr[0]) ? arr : +// the sorting is done using compare_vals() +function _indexed_sort(arrind) = + arrind==[] ? [] : len(arrind)==1? [arrind[0][0]] : + let( pivot = arrind[floor(len(arrind)/2)][1] ) let( - pivot = arr[floor(len(arr)/2)][_i], - lesser = [ for (entry=arr) if (compare_vals(entry[_i], pivot) <0 ) entry ], - equal = [ for (entry=arr) if (compare_vals(entry[_i], pivot)==0 ) entry ], - greater = [ for (entry=arr) if (compare_vals(entry[_i], pivot) >0 ) entry ] + lesser = [ for (entry=arrind) if (compare_vals(entry[1], pivot) <0 ) entry ], + equal = [ for (entry=arrind) if (compare_vals(entry[1], pivot)==0 ) entry[0] ], + greater = [ for (entry=arrind) if (compare_vals(entry[1], pivot) >0 ) entry ] ) - concat( - _lexical_sort(lesser, _i ), - _lexical_sort(equal, _i+1 ), - _lexical_sort(greater, _i ) - ); - - -// given a list of pairs, return the first element of each pair of the list sorted by the second element of the pair -// it uses compare_vals() -function _indexed_sort(arr, _i=0) = - arr==[] ? [] : - len(arr)==1? [arr[0][0]] : - _i>=len(arr[0][1]) ? [for(ai=arr) ai[0]] : - let( - pivot = arr[floor(len(arr)/2)][1][_i], - lesser = [ for (entry=arr) if (compare_vals(entry[1][_i], pivot) <0 ) entry ], - equal = [ for (entry=arr) if (compare_vals(entry[1][_i], pivot)==0 ) entry ], - greater = [ for (entry=arr) if (compare_vals(entry[1][_i], pivot) >0 ) entry ] - ) - concat( - _indexed_sort(lesser, _i ), - _indexed_sort(equal, _i+1 ), - _indexed_sort(greater, _i ) ); + concat(_indexed_sort(lesser), equal, _indexed_sort(greater)); // returns true for valid index specifications idx in the interval [imin, imax) // note that idx can't have any value greater or EQUAL to imax -// this allows imax=INF as a bound to numerical lists function _valid_idx(idx,imin,imax) = is_undef(idx) || ( is_finite(idx) && idx>=imin && idx< imax ) || ( is_list(idx) && min(idx)>=imin && max(idx)< imax ) - || ( ! is_list(idx) // it implicitly is a range - && (idx[1]>0 && idx[0]>=imin && idx[2]< imax) - || - (idx[0]=imin) ); + || ( valid_range(idx) && idx[0]>=imin && idx[2]< imax ); // Function: sort() @@ -830,7 +758,7 @@ function _valid_idx(idx,imin,imax) = function sort(list, idx=undef) = !is_list(list) || len(list)<=1 ? list : is_def(idx) - ? assert( _valid_idx(idx,0,len(list[0])) , "Invalid indices or out of range.") + ? assert( _valid_idx(idx,0,len(list)) , "Invalid indices.") let( sarr = _sort_general(list,idx) ) [for(i=[0:len(sarr)-1]) list[sarr[i]] ] : let(size = array_dim(list)) @@ -845,31 +773,6 @@ function sort(list, idx=undef) = ) : _sort_general(list); -function sort(list, idx=undef) = - !is_list(list) || len(list)<=1 ? list : - is_vector(list) - ? assert( _valid_idx(idx,0,len(list[0])) , str("Invalid indices or out of range. ",list)) - is_def(idx) - ? sort_vector_indexed([for(i=[0:len(list)-1]) [i, [list(i)] ] ]) - : sort_scalar(list) - : is_matrix(list) - ? list==[] || list[0]==[] ? list : - assert( _valid_idx(idx,0,len(list[0])) , "Invalid indices or out of range.") - is_def(idx) - ? sort_vector_indexed([for(i=[0:len(list)-1], li=[list[i]]) [i, [for(ind=idx) li(ind)] ] ]) - : sort_vector(list) - : list==[] || list[0]==[] ? list : - let( llen = [for(li=list) !is_list(li) || is_string(li) ? 0: len(li)], - m = min(llen), - M = max(llen) - ) - M==0 ? _simple_sort(list) : - assert( m>0 && _valid_idx(idx,m-1,M) , "Invalid indices or out of range.") - is_def(idx) - ? _sort_general(list,idx) - : let( ils = _sort_general(list, [m:M]) ) - [for(i=[0:len(list)-1]) list[ils[i]] ]; - // Function: sortidx() // Description: diff --git a/math.scad b/math.scad index e5c4625..888d048 100644 --- a/math.scad +++ b/math.scad @@ -1175,7 +1175,7 @@ function deriv3(data, h=1, closed=false) = // Description: // Multiplies two complex numbers represented by 2D vectors. function C_times(z1,z2) = - assert( is_matrix([z1,z2],2,2), "Complex numbers should be represented by 2D vectors" ) + assert( is_vector(z1+z2,2), "Complex numbers should be represented by 2D vectors." ) [ z1.x*z2.x - z1.y*z2.y, z1.x*z2.y + z1.y*z2.x ]; // Function: C_div() diff --git a/quaternions.scad b/quaternions.scad index 44b0e29..e21801c 100644 --- a/quaternions.scad +++ b/quaternions.scad @@ -27,7 +27,7 @@ function _Qreal(q) = q[3]; function _Qset(v,r) = concat( v, r ); // normalizes without checking -function _Qunit(q) = q/norm(q); +function _Qnorm(q) = q/norm(q); // Function: Q_is_quat() @@ -36,7 +36,7 @@ function _Qunit(q) = q/norm(q); // 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 Q_is_quat(q) = is_vector(q,4) && ! approx(norm(q),0) ; // Function: Quat() @@ -101,7 +101,7 @@ function QuatXYZ(a=[0,0,0]) = qy = QuatY(a[1]), qz = QuatZ(a[2]) ) - _Qmul(qz, _Qmul(qy, qx)); + Q_Mul(qz, Q_Mul(qy, qx)); // Function: Q_From_to() @@ -202,36 +202,32 @@ function Q_Sub(a, 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)") - _Qmul(a,b); - -function _Qmul(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 + a[3]*b[3] - a.x*b.x - a.y*b.y - a.z*b.z, ]; -// [ [a[3], -a.z, a.y, a.x], -// [ a.z, a[3], -a.x, a.y], -// [-a.y, a.x, a[3], a.z], -// [-a.x, -a.y, -a.z, a[3]] ]*[b.x,b.y,b.z,b[3]] - + // Function: Q_Cumulative() // Usage: -// Q_Cumulative(ql); +// Q_Cumulative(v); // Description: // 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(ql) = - assert( is_matrix(ql,undef,4) && len(ql)>0, "Invalid list of quaternions." ) - [for( i = 0, q = ql[0]; - i<=len(ql); - i = i+1, q = (i==len(ql))? 0: _Qmul(q,ql[i]) ) - q ]; +function Q_Cumulative(v, _i=0, _acc=[]) = + _i==len(v) ? _acc : + Q_Cumulative( + v, _i+1, + concat( + _acc, + [_i==0 ? v[_i] : Q_Mul(v[_i], select(_acc,-1))] + ) + ); // Function: Q_Dot() @@ -256,7 +252,7 @@ function Q_Neg(q) = // Q_Conj(q) // Description: Returns the conjugate of quaternion `q`. function Q_Conj(q) = - assert( Q_is_quat(q), str("Invalid quaternion",q) ) + assert( Q_is_quat(q), "Invalid quaternion" ) [-q.x, -q.y, -q.z, q[3]]; @@ -266,7 +262,7 @@ function Q_Conj(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 = q/norm(q) ) + let(q = _Qnorm(q) ) [-q.x, -q.y, -q.z, q[3]]; @@ -287,9 +283,7 @@ function Q_Norm(q) = // Description: Normalizes quaternion `q`, so that norm([W,X,Y,Z]) == 1. function Q_Normalize(q) = assert( Q_is_quat(q) , "Invalid quaternion" ) - approx(_Qvec(q), [0,0,0]) - ? Q_Ident() - : q/norm(q); + q/norm(q); // Function: Q_Dist() @@ -324,32 +318,31 @@ function Q_Dist(q1, q2) = // Qrot(q) right(80) cube([10,10,1]); // #sphere(r=80); function Q_Slerp(q1, q2, u, _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), - q3 = dot>0.9995? q2: _Qunit(q2 - dot*q1) - ) - ! is_num(u) - ? [for (uu=u) _Qslerp(q1, q3, uu, dot)] - : _Qslerp(q1, q3, 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 _Qslerp(q1, q2, u, dot) = - dot>0.9995 - ? _Qunit(q1 + u*(q2-q1)) - : let( theta = u*acos(dot) ) - _Qunit(q1*cos(theta) + q2*sin(theta)); - -// Function: Q_to_matrix3() +// Function: Q_Matrix3() // Usage: -// Q_to_matrix3(q); +// Q_Matrix3(q); // Description: // Returns the 3x3 rotation matrix for the given normalized quaternion q. -function Q_to_matrix3(q) = +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]], @@ -358,12 +351,12 @@ function Q_to_matrix3(q) = ]; -// Function: Q_to_matrix4() +// Function: Q_Matrix4() // Usage: -// Q_to_matrix4(q); +// Q_Matrix4(q); // Description: // Returns the 4x4 rotation matrix for the given normalized quaternion q. -function Q_to_matrix4(q) = +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], @@ -373,35 +366,6 @@ function Q_to_matrix4(q) = ]; -// Function: Q_from_matrix() -// Usage: -// Q_from_matrix(M) -// Description: -// Returns a normalized quaternion corresponding to the rotation matrix M. -// M may be a 3x3 rotation matrix or a homogeneous 4x4 rotation matrix. -// The last row and last column of M are ignored for 4x4 matrices. -// It doesn't check whether M is in fact a rotation matrix. -// If M is not a rotation, the returned quaternion is unpredictable. -// -// based on https://en.wikipedia.org/wiki/Rotation_matrix -// -function Q_from_matrix(M) = - assert( is_matrix(M) && (len(M)==3 || len(M)==4) && (len(M[0])==3 || len(M[0])==4), - "The matrix should be 3x3 or 4x4") - let( tr = M[0][0]+M[1][1]+M[2][2] ) // M trace - tr>0 - ? let( r = sqrt(1+tr), s = -1/r/2 ) - _Qunit( _Qset([ M[1][2]-M[2][1], M[2][0]-M[0][2], M[0][1]-M[1][0] ]*s, r/2 ) ) - : let( - i = max_index([ M[0][0], M[1][1], M[2][2] ]), - r = sqrt(1 + 2*M[i][i] -M[0][0] -M[1][1] -M[2][2] ), - s = 1/r/2 - ) - i==0 ? _Qunit( _Qset( [ r/2, s*(M[1][0]+M[0][1]), s*(M[0][2]+M[2][0]) ], s*(M[2][1]-M[1][2])) ): - i==1 ? _Qunit( _Qset( [ s*(M[1][0]+M[0][1]), r/2, s*(M[2][1]+M[1][2]) ], s*(M[0][2]-M[2][0])) ) - : _Qunit( _Qset( [ s*(M[2][0]+M[0][2]), s*(M[1][2]+M[2][1]), r/2 ], s*(M[1][0]-M[0][1])) ) ; - - // Function: Q_Axis() // Usage: // Q_Axis(q) @@ -454,15 +418,15 @@ function Q_Angle(q1,q2) = // q = QuatXYZ([45,35,10]); // pts = Qrot(q, p=[[2,3,4], [4,5,6], [9,2,3]]); module Qrot(q) { - multmatrix(Q_to_matrix4(q)) { + multmatrix(Q_Matrix4(q)) { children(); } } function Qrot(q,p) = - is_undef(p)? Q_to_matrix4(q) : + is_undef(p)? Q_Matrix4(q) : is_vector(p)? Qrot(q,[p])[0] : - apply(Q_to_matrix4(q), p); + apply(Q_Matrix4(q), p); // Module: Qrot_copies() @@ -482,6 +446,29 @@ 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); @@ -531,11 +518,11 @@ function Q_Rotation_path(q1, n=1, q2) = 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_to_matrix4(q1), R=dR; i<=n; i=i+1, R=dR*R ) R] + ? [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( _Qmul( q2, Q_Inverse(q1) ), 1/n ), - dR = Q_to_matrix4(dq) ) - [for( i=0, R=Q_to_matrix4(q1); i<=n; i=i+1, R=dR*R ) R]; + 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)) @@ -578,8 +565,8 @@ function Q_Nlerp(q1,q2,u) = let( q1 = Q_Normalize(q1), q2 = Q_Normalize(q2) ) is_num(u) - ? _Qunit((1-u)*q1 + u*q2 ) - : [for (ui=u) _Qunit((1-ui)*q1 + ui*q2 ) ]; + ? _Qnorm((1-u)*q1 + u*q2 ) + : [for (ui=u) _Qnorm((1-ui)*q1 + ui*q2 ) ]; // Function: Q_Squad() @@ -629,17 +616,6 @@ function Q_Squad(q1,q2,q3,q4,u) = : [for(ui=u) Q_Slerp( Q_Slerp(q1,q4,ui), Q_Slerp(q2,q3,ui), 2*ui*(1-ui) ) ]; -function Q_Scubic(q1,q2,q3,q4,u) = - assert(is_finite(u) || is_range(u) || is_vector(u) , - "Invalid interpolation coefficient(s)" ) - is_num(u) - ? let( q12 = Q_Slerp(q1,q2,u), - q23 = Q_Slerp(q2,q3,u), - q34 = Q_Slerp(q3,q4,u) ) - Q_Slerp(Q_Slerp(q12,q23,u),Q_Slerp(q23,q34,u),u) - : [for(ui=u) Q_Scubic( q1,q2,q3,q4,ui) ]; - - // Function: Q_exp() // Usage: // q2 = Q_exp(q); diff --git a/tests/test_all.scad b/tests/test_all.scad deleted file mode 100644 index fe6512b..0000000 --- a/tests/test_all.scad +++ /dev/null @@ -1,28 +0,0 @@ -//include -include -include -include -include -include -include -include -include -include -include -include -include -include -//include -include -include -include -include -include -include -include -include -include -include -include - - diff --git a/tests/test_quaternions.scad b/tests/test_quaternions.scad index af19408..b1e5130 100644 --- a/tests/test_quaternions.scad +++ b/tests/test_quaternions.scad @@ -1,80 +1,6 @@ include <../std.scad> include <../strings.scad> -function plane_fit(points) = - assert( is_matrix(points,undef,3) , "Improper point list." ) - len(points)< 2 ? [] : - len(points)==3 ? plane3pt(points[0],points[1],points[2]) : - let( - A = [for(pi=points) concat(pi,-1) ], - qr = qr_factor(A), - R = select(qr[1],0,3), -x=[for(ri=R) echo(R=ri)0], -//y=[for(qi=qr[0]) echo(Q=qi)0], - s = back_substitute - ( R, - sum([for(qi=qr[0])[for(j=[0:3]) qi[j] ] ]) - ) -//, s0 = echo(ls=linear_solve(A,[for(p=points) 0])) -,w=s==[]? echo("s == []")0:echo(s=s)0 - ) - s==[] - ? // points are collinear - [] - : // plane through the origin? - approx(norm([s.x, s.y, s.z]), 0) - ? let( - k = max_index([for(i=[0:2]) abs(s[i]) ]), - A = [[1,0,0], for(pi=points) pi ], - b = [1, for(i=[1:len(points)]) 0], - s = linear_solve(A,b) -, y=echo(s2=s) - ) - s==[]? []: - concat(s, 0)/norm(s) - : s/norm([s.x, s.y, s.z]) ; - -function plane_fit(points) = - assert( is_matrix(points,undef,3) , "Improper point list." ) - len(points)< 2 ? [] : - len(points)==3 ? plane3pt(points[0],points[1],points[2]) : - let( - A = [for(pi=points) concat(pi,-1) ], - B = [ for(pi=points) [0,1] ], - sB = transpose(linear_solve(A,B)), - s = sB==[] ? [] - : norm(sB[1])N/2? 10*[1,1,1]: rands(-1,1,3) ]; -pf = plane_fit(pts); -echo(pf); -//pm = concat(sum(pts)/len(pts), -1); -//echo(pmfit=pm*pf); function rec_cmp(a,b,eps=1e-9) = typeof(a)!=typeof(b)? false : @@ -82,18 +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 standardize(v) = - v==[] - ? [] - : sign(first_nonzero(v))*v; - -function first_nonzero(v) = - v==[] ? 0 - : is_num(v) ? v - : [for(vi=v) if(!is_list(vi) && ! approx(vi,0)) vi - else if(is_list(vi)) first_nonzero(vi), 0 ][0]; - -module assert_std(vc,ve) { assert_approx(standardize(vc),standardize(ve)); } +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)) { @@ -110,7 +25,7 @@ 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]),true); + 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); } @@ -189,7 +104,7 @@ 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_to_matrix4(Q_From_to([1,2,3], [4,5,2])), rot(from=[1,2,3],to=[4,5,2])); + 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])); } @@ -285,7 +200,7 @@ test_Q_Mul(); module test_Q_Cumulative() { verify_f(Q_Cumulative([QuatZ(30),QuatX(57),QuatY(18)]),[[0, 0, 0.2588190451, 0.9659258263], [0.4608999698, -0.1234977747, 0.2274546059, 0.8488721457], [0.4908072659, 0.01081554785, 0.1525536221, 0.8577404293]]); } -*test_Q_Cumulative(); +test_Q_Cumulative(); module test_Q_Dot() { @@ -304,18 +219,18 @@ test_Q_Neg(); module test_Q_Conj() { - ang = rands(0,360,3); verify_f(Q_Conj([1,0,0,1]),[-1,0,0,1]); verify_f(Q_Conj([0,1,1,0]),[0,-1,-1,0]); verify_f(Q_Conj(QuatXYZ([23,45,67])),[0.0533818345, -0.4143703268, -0.4360652669, 0.7970537592]); - verify_f(Q_Mul(Q_Conj(QuatXYZ(ang)),QuatXYZ(ang)),Q_Ident()); } test_Q_Conj(); module test_Q_Inverse() { - verify_f(Q_Inverse([1,0,0,1]),[-1,0,0,1]); - verify_f(Q_Inverse([0,1,1,0]),[0,-1,-1,0]); + + 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(); @@ -332,7 +247,7 @@ test_Q_Norm(); module test_Q_Normalize() { verify_f(Q_Normalize([1,0,0,1]),[0.7071067812, 0, 0, 0.7071067812]); verify_f(Q_Normalize([0,1,1,0]),[0, 0.7071067812, 0.7071067812, 0]); -// verify_f(Q_Normalize(QuatXYZ([23,45,67])),[-0.0533818345, 0.4143703268, 0.4360652669, 0.7970537592]); + verify_f(Q_Normalize(QuatXYZ([23,45,67])),[-0.0533818345, 0.4143703268, 0.4360652669, 0.7970537592]); } test_Q_Normalize(); @@ -345,47 +260,28 @@ test_Q_Dist(); module test_Q_Slerp() { - u = rands(0,1,1)[0]; - ul = rands(0,1,5); - ul2 = [1,1,1,1,1]-ul; - a1 = rands(0,360,1)[0]; - a2 = rands(0,360,1)[0]; - a3 = rands(0,360,1)[0]; - verify_f(standardize(Q_Slerp(QuatX(a1),QuatY(a2),0.0)), standardize(QuatX(a1))); - verify_f(standardize(Q_Slerp(QuatX(45),QuatY(30),0.5)), - [0.1967063121, 0.1330377423, 0, 0.9713946602]); - verify_f(standardize(Q_Slerp(QuatX(a1),QuatY(a2),1.0)), standardize(QuatY(a2))); - verify_f(standardize(Q_Slerp(QuatXYZ([a1,a2,0]),-QuatY(a2),u)), - standardize(Q_Slerp(QuatXYZ([a1,a2,0]), QuatY(a2),u))); - verify_f(standardize(Q_Slerp(QuatX(a1),QuatX(a1),u)), standardize(QuatX(a1))); - verify_f(standardize(Q_Slerp(QuatXYZ([a1,a2,0]),QuatXYZ([a2,0,a1]),u)), - standardize(Q_Slerp(QuatXYZ([a2,0,a1]),QuatXYZ([a1,a2,0]),1-u))); - verify_f(standardize(Q_Slerp(QuatXYZ([a1,a2,0]),QuatXYZ([a2,0,a1]),ul)), - standardize(Q_Slerp(QuatXYZ([a2,0,a1]),QuatXYZ([a1,a2,0]),ul2))); + verify_f(Q_Slerp(QuatX(45),QuatY(30),0.0),QuatX(45)); + verify_f(Q_Slerp(QuatX(45),QuatY(30),0.5),[0.1967063121, 0.1330377423, 0, 0.9713946602]); + verify_f(Q_Slerp(QuatX(45),QuatY(30),1.0),QuatY(30)); } test_Q_Slerp(); -module test_Q_to_matrix3() { - rotZ_37 = rot(37); - rotZ_37_3 = [for(i=[0:2]) [for(j=[0:2]) rotZ_37[i][j] ] ]; - angs = [12,-49,40]; - rot4 = rot(angs); - rot3 = [for(i=[0:2]) [for(j=[0:2]) rot4[i][j] ] ]; - verify_f(Q_to_matrix3(QuatZ(37)),rotZ_37_3); - verify_f(Q_to_matrix3(QuatXYZ(angs)),rot3); +module test_Q_Matrix3() { + verify_f(Q_Matrix3(QuatZ(37)),rot(37,planar=true)); + verify_f(Q_Matrix3(QuatZ(-49)),rot(-49,planar=true)); } -test_Q_to_matrix3(); +test_Q_Matrix3(); -module test_Q_to_matrix4() { - verify_f(Q_to_matrix4(QuatZ(37)),rot(37)); - verify_f(Q_to_matrix4(QuatZ(-49)),rot(-49)); - verify_f(Q_to_matrix4(QuatX(37)),rot([37,0,0])); - verify_f(Q_to_matrix4(QuatY(37)),rot([0,37,0])); - verify_f(Q_to_matrix4(QuatXYZ([12,34,56])),rot([12,34,56])); +module test_Q_Matrix4() { + verify_f(Q_Matrix4(QuatZ(37)),rot(37)); + verify_f(Q_Matrix4(QuatZ(-49)),rot(-49)); + verify_f(Q_Matrix4(QuatX(37)),rot([37,0,0])); + verify_f(Q_Matrix4(QuatY(37)),rot([0,37,0])); + verify_f(Q_Matrix4(QuatXYZ([12,34,56])),rot([12,34,56])); } -test_Q_to_matrix4(); +test_Q_Matrix4(); module test_Q_Axis() { @@ -425,23 +321,23 @@ module test_Qrot() { test_Qrot(); -module test_Q_from_matrix() { - verify_f(standardize(Q_from_matrix(Q_to_matrix3(Quat([12,34,56],33)))),standardize(Quat([12,34,56],33))); - verify_f(Q_to_matrix3(Q_from_matrix(Q_to_matrix3(QuatXYZ([12,34,56])))), - Q_to_matrix3(QuatXYZ([12,34,56]))); +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_from_matrix(); +test_Q_Rotation(); module test_Q_Rotation_path() { - verify_f(Q_Rotation_path(QuatX(135), 5, QuatY(13.5))[0] , Q_to_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), 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)[3] , xrot(11+(55-11)*3/4)); verify_f(Q_Rotation_path(QuatX(11), 5)[4] , xrot(55)); } @@ -451,54 +347,48 @@ 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_to_matrix4(Q_Nlerp(QuatX(135), QuatY(13.5),0.5))); + 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() { - u = rands(0,1,5); - su = [1,1,1,1,1]-u; - verify_f(Q_Squad(QuatX(45),QuatZ(30),QuatX(32),QuatY(30),0.0),QuatX(45)); + 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), + 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(10),QuatZ(20),QuatX(32),QuatXYZ([210,120,30]),u[0] ), - Q_Squad(QuatXYZ([210,120,30]),QuatX(32),QuatZ(20),QuatY(10),1-u[0] ) ); - verify_f( Q_Squad(QuatY(10),QuatZ(20),QuatX(32),QuatXYZ([210,120,30]),u ), - Q_Squad(QuatXYZ([210,120,30]),QuatX(32),QuatZ(20),QuatY(10),su ) ); + 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() { - q=QuatXYZ(rands(0,360,3)); 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,4])), [1,2,3,4]); - verify_f(Q_exp(Q_ln(q)), q); - verify_f(Q_exp(q+Q_Inverse(q)),Q_Mul(Q_exp(q),Q_exp(Q_Inverse(q)))); + 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() { - q = QuatXYZ(rands(0,360,3)); 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(q)), q); - verify_f(Q_ln(q)+Q_ln(Q_Conj(q)), [0,0,0,0]); + 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 = QuatXYZ(rands(0,360,3)); + 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)); From e8499baba185848ff6f6f1745181a33573ee5c58 Mon Sep 17 00:00:00 2001 From: RonaldoCMP Date: Wed, 26 Aug 2020 14:09:12 +0100 Subject: [PATCH 03/23] Revert "Revert "Correction of C_times validation"" This reverts commit 4a366967f96c29c61116125b0803a1dd405e1f91. --- arrays.scad | 127 ++++++++++++++++++++--- math.scad | 2 +- quaternions.scad | 170 +++++++++++++++++-------------- tests/test_all.scad | 28 ++++++ tests/test_quaternions.scad | 196 ++++++++++++++++++++++++++++-------- 5 files changed, 391 insertions(+), 132 deletions(-) create mode 100644 tests/test_all.scad diff --git a/arrays.scad b/arrays.scad index 835fd65..3ae7a88 100644 --- a/arrays.scad +++ b/arrays.scad @@ -708,36 +708,108 @@ function _sort_vectors4(arr) = y ] ) concat( _sort_vectors4(lesser), equal, _sort_vectors4(greater) ); +// sort a list of vectors +function _sort_vectors(arr, _i=0) = + len(arr)<=1 || _i>=len(arr[0]) ? arr : + let( + pivot = arr[floor(len(arr)/2)][_i], + lesser = [ for (entry=arr) if (entry[_i] < pivot ) entry ], + equal = [ for (entry=arr) if (entry[_i] == pivot ) entry ], + greater = [ for (entry=arr) if (entry[_i] > pivot ) entry ] + ) + concat( + _sort_vectors(lesser, _i ), + _sort_vectors(equal, _i+1 ), + _sort_vectors(greater, _i ) ); + +// given pairs of an index and a vector, return the list of indices of the list sorted by the vectors +function _sort_vectors_indexed(arr, _i=0) = + arr==[] ? [] : + len(arr)==1 || _i>=len(arr[0][1]) ? [for(ai=arr) ai[0]] : + let( + pivot = arr[floor(len(arr)/2)][1][_i], + lesser = [ for (entry=arr) if (entry[1][_i] < pivot ) entry ], + equal = [ for (entry=arr) if (entry[1][_i] == pivot ) entry ], + greater = [ for (entry=arr) if (entry[1][_i] > pivot ) entry ] + ) + concat( + _sort_vectors_indexed(lesser, _i ), + _sort_vectors_indexed(equal, _i+1 ), + _sort_vectors_indexed(greater, _i ) ); + // when idx==undef, returns the sorted array // otherwise, returns the indices of the sorted array function _sort_general(arr, idx=undef) = - (len(arr)<=1) ? arr : + len(arr)<=1 ? arr : is_undef(idx) - ? _sort_scalar(arr) + ? _simple_sort(arr) +// : _lexical_sort(arr) : let( arrind=[for(k=[0:len(arr)-1], ark=[arr[k]]) [ k, [for (i=idx) ark[i]] ] ] ) - _indexed_sort(arrind); - -// given a list of pairs, return the first element of each pair of the list sorted by the second element of the pair -// the sorting is done using compare_vals() -function _indexed_sort(arrind) = - arrind==[] ? [] : len(arrind)==1? [arrind[0][0]] : - let( pivot = arrind[floor(len(arrind)/2)][1] ) + _indexed_sort(arrind ); + +// sort simple lists with compare_vals() +function _simple_sort(arr) = + arr==[] || len(arr)==1 ? arr : let( - lesser = [ for (entry=arrind) if (compare_vals(entry[1], pivot) <0 ) entry ], - equal = [ for (entry=arrind) if (compare_vals(entry[1], pivot)==0 ) entry[0] ], - greater = [ for (entry=arrind) if (compare_vals(entry[1], pivot) >0 ) entry ] + pivot = arr[floor(len(arr)/2)], + lesser = [ for (entry=arr) if (compare_vals(entry, pivot) <0 ) entry ], + equal = [ for (entry=arr) if (compare_vals(entry, pivot)==0 ) entry ], + greater = [ for (entry=arr) if (compare_vals(entry, pivot) >0 ) entry ] ) - concat(_indexed_sort(lesser), equal, _indexed_sort(greater)); + concat( + _simple_sort(lesser), + equal, + _simple_sort(greater) + ); + + +// given a list of pairs, return the first element of each pair of the list sorted by the second element of the pair +// it uses compare_vals() +function _lexical_sort(arr, _i=0) = + arr==[] || len(arr)==1 || _i>=len(arr[0]) ? arr : + let( + pivot = arr[floor(len(arr)/2)][_i], + lesser = [ for (entry=arr) if (compare_vals(entry[_i], pivot) <0 ) entry ], + equal = [ for (entry=arr) if (compare_vals(entry[_i], pivot)==0 ) entry ], + greater = [ for (entry=arr) if (compare_vals(entry[_i], pivot) >0 ) entry ] + ) + concat( + _lexical_sort(lesser, _i ), + _lexical_sort(equal, _i+1 ), + _lexical_sort(greater, _i ) + ); + + +// given a list of pairs, return the first element of each pair of the list sorted by the second element of the pair +// it uses compare_vals() +function _indexed_sort(arr, _i=0) = + arr==[] ? [] : + len(arr)==1? [arr[0][0]] : + _i>=len(arr[0][1]) ? [for(ai=arr) ai[0]] : + let( + pivot = arr[floor(len(arr)/2)][1][_i], + lesser = [ for (entry=arr) if (compare_vals(entry[1][_i], pivot) <0 ) entry ], + equal = [ for (entry=arr) if (compare_vals(entry[1][_i], pivot)==0 ) entry ], + greater = [ for (entry=arr) if (compare_vals(entry[1][_i], pivot) >0 ) entry ] + ) + concat( + _indexed_sort(lesser, _i ), + _indexed_sort(equal, _i+1 ), + _indexed_sort(greater, _i ) ); // returns true for valid index specifications idx in the interval [imin, imax) // note that idx can't have any value greater or EQUAL to imax +// this allows imax=INF as a bound to numerical lists function _valid_idx(idx,imin,imax) = is_undef(idx) || ( is_finite(idx) && idx>=imin && idx< imax ) || ( is_list(idx) && min(idx)>=imin && max(idx)< imax ) - || ( valid_range(idx) && idx[0]>=imin && idx[2]< imax ); + || ( ! is_list(idx) // it implicitly is a range + && (idx[1]>0 && idx[0]>=imin && idx[2]< imax) + || + (idx[0]=imin) ); // Function: sort() @@ -758,7 +830,7 @@ function _valid_idx(idx,imin,imax) = function sort(list, idx=undef) = !is_list(list) || len(list)<=1 ? list : is_def(idx) - ? assert( _valid_idx(idx,0,len(list)) , "Invalid indices.") + ? assert( _valid_idx(idx,0,len(list[0])) , "Invalid indices or out of range.") let( sarr = _sort_general(list,idx) ) [for(i=[0:len(sarr)-1]) list[sarr[i]] ] : let(size = array_dim(list)) @@ -773,6 +845,31 @@ function sort(list, idx=undef) = ) : _sort_general(list); +function sort(list, idx=undef) = + !is_list(list) || len(list)<=1 ? list : + is_vector(list) + ? assert( _valid_idx(idx,0,len(list[0])) , str("Invalid indices or out of range. ",list)) + is_def(idx) + ? sort_vector_indexed([for(i=[0:len(list)-1]) [i, [list(i)] ] ]) + : sort_scalar(list) + : is_matrix(list) + ? list==[] || list[0]==[] ? list : + assert( _valid_idx(idx,0,len(list[0])) , "Invalid indices or out of range.") + is_def(idx) + ? sort_vector_indexed([for(i=[0:len(list)-1], li=[list[i]]) [i, [for(ind=idx) li(ind)] ] ]) + : sort_vector(list) + : list==[] || list[0]==[] ? list : + let( llen = [for(li=list) !is_list(li) || is_string(li) ? 0: len(li)], + m = min(llen), + M = max(llen) + ) + M==0 ? _simple_sort(list) : + assert( m>0 && _valid_idx(idx,m-1,M) , "Invalid indices or out of range.") + is_def(idx) + ? _sort_general(list,idx) + : let( ils = _sort_general(list, [m:M]) ) + [for(i=[0:len(list)-1]) list[ils[i]] ]; + // Function: sortidx() // Description: diff --git a/math.scad b/math.scad index 888d048..e5c4625 100644 --- a/math.scad +++ b/math.scad @@ -1175,7 +1175,7 @@ function deriv3(data, h=1, closed=false) = // Description: // Multiplies two complex numbers represented by 2D vectors. function C_times(z1,z2) = - assert( is_vector(z1+z2,2), "Complex numbers should be represented by 2D vectors." ) + assert( is_matrix([z1,z2],2,2), "Complex numbers should be represented by 2D vectors" ) [ z1.x*z2.x - z1.y*z2.y, z1.x*z2.y + z1.y*z2.x ]; // Function: C_div() diff --git a/quaternions.scad b/quaternions.scad index e21801c..44b0e29 100644 --- a/quaternions.scad +++ b/quaternions.scad @@ -27,7 +27,7 @@ function _Qreal(q) = q[3]; function _Qset(v,r) = concat( v, r ); // normalizes without checking -function _Qnorm(q) = q/norm(q); +function _Qunit(q) = q/norm(q); // Function: Q_is_quat() @@ -36,7 +36,7 @@ function _Qnorm(q) = q/norm(q); // 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 Q_is_quat(q) = is_vector(q,4) ;//&& ! approx(norm(q),0) ; // Function: Quat() @@ -101,7 +101,7 @@ function QuatXYZ(a=[0,0,0]) = qy = QuatY(a[1]), qz = QuatZ(a[2]) ) - Q_Mul(qz, Q_Mul(qy, qx)); + _Qmul(qz, _Qmul(qy, qx)); // Function: Q_From_to() @@ -202,32 +202,36 @@ function Q_Sub(a, 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)") + _Qmul(a,b); + +function _Qmul(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, + a[3]*b[3] - a.x*b.x - a.y*b.y - a.z*b.z ]; - +// [ [a[3], -a.z, a.y, a.x], +// [ a.z, a[3], -a.x, a.y], +// [-a.y, a.x, a[3], a.z], +// [-a.x, -a.y, -a.z, a[3]] ]*[b.x,b.y,b.z,b[3]] + // Function: Q_Cumulative() // Usage: -// Q_Cumulative(v); +// Q_Cumulative(ql); // Description: // 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( - v, _i+1, - concat( - _acc, - [_i==0 ? v[_i] : Q_Mul(v[_i], select(_acc,-1))] - ) - ); +function Q_Cumulative(ql) = + assert( is_matrix(ql,undef,4) && len(ql)>0, "Invalid list of quaternions." ) + [for( i = 0, q = ql[0]; + i<=len(ql); + i = i+1, q = (i==len(ql))? 0: _Qmul(q,ql[i]) ) + q ]; // Function: Q_Dot() @@ -252,7 +256,7 @@ function Q_Neg(q) = // Q_Conj(q) // Description: Returns the conjugate of quaternion `q`. function Q_Conj(q) = - assert( Q_is_quat(q), "Invalid quaternion" ) + assert( Q_is_quat(q), str("Invalid quaternion",q) ) [-q.x, -q.y, -q.z, q[3]]; @@ -262,7 +266,7 @@ function Q_Conj(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) ) +// let(q = q/norm(q) ) [-q.x, -q.y, -q.z, q[3]]; @@ -283,7 +287,9 @@ function Q_Norm(q) = // Description: Normalizes quaternion `q`, so that norm([W,X,Y,Z]) == 1. function Q_Normalize(q) = assert( Q_is_quat(q) , "Invalid quaternion" ) - q/norm(q); + approx(_Qvec(q), [0,0,0]) + ? Q_Ident() + : q/norm(q); // Function: Q_Dist() @@ -318,31 +324,32 @@ function Q_Dist(q1, q2) = // Qrot(q) right(80) cube([10,10,1]); // #sphere(r=80); 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)); + 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), + q3 = dot>0.9995? q2: _Qunit(q2 - dot*q1) + ) + ! is_num(u) + ? [for (uu=u) _Qslerp(q1, q3, uu, dot)] + : _Qslerp(q1, q3, u, dot); +function _Qslerp(q1, q2, u, dot) = + dot>0.9995 + ? _Qunit(q1 + u*(q2-q1)) + : let( theta = u*acos(dot) ) + _Qunit(q1*cos(theta) + q2*sin(theta)); + -// Function: Q_Matrix3() +// Function: Q_to_matrix3() // Usage: -// Q_Matrix3(q); +// Q_to_matrix3(q); // Description: // Returns the 3x3 rotation matrix for the given normalized quaternion q. -function Q_Matrix3(q) = +function Q_to_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]], @@ -351,12 +358,12 @@ function Q_Matrix3(q) = ]; -// Function: Q_Matrix4() +// Function: Q_to_matrix4() // Usage: -// Q_Matrix4(q); +// Q_to_matrix4(q); // Description: // Returns the 4x4 rotation matrix for the given normalized quaternion q. -function Q_Matrix4(q) = +function Q_to_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], @@ -366,6 +373,35 @@ function Q_Matrix4(q) = ]; +// Function: Q_from_matrix() +// Usage: +// Q_from_matrix(M) +// Description: +// Returns a normalized quaternion corresponding to the rotation matrix M. +// M may be a 3x3 rotation matrix or a homogeneous 4x4 rotation matrix. +// The last row and last column of M are ignored for 4x4 matrices. +// It doesn't check whether M is in fact a rotation matrix. +// If M is not a rotation, the returned quaternion is unpredictable. +// +// based on https://en.wikipedia.org/wiki/Rotation_matrix +// +function Q_from_matrix(M) = + assert( is_matrix(M) && (len(M)==3 || len(M)==4) && (len(M[0])==3 || len(M[0])==4), + "The matrix should be 3x3 or 4x4") + let( tr = M[0][0]+M[1][1]+M[2][2] ) // M trace + tr>0 + ? let( r = sqrt(1+tr), s = -1/r/2 ) + _Qunit( _Qset([ M[1][2]-M[2][1], M[2][0]-M[0][2], M[0][1]-M[1][0] ]*s, r/2 ) ) + : let( + i = max_index([ M[0][0], M[1][1], M[2][2] ]), + r = sqrt(1 + 2*M[i][i] -M[0][0] -M[1][1] -M[2][2] ), + s = 1/r/2 + ) + i==0 ? _Qunit( _Qset( [ r/2, s*(M[1][0]+M[0][1]), s*(M[0][2]+M[2][0]) ], s*(M[2][1]-M[1][2])) ): + i==1 ? _Qunit( _Qset( [ s*(M[1][0]+M[0][1]), r/2, s*(M[2][1]+M[1][2]) ], s*(M[0][2]-M[2][0])) ) + : _Qunit( _Qset( [ s*(M[2][0]+M[0][2]), s*(M[1][2]+M[2][1]), r/2 ], s*(M[1][0]-M[0][1])) ) ; + + // Function: Q_Axis() // Usage: // Q_Axis(q) @@ -418,15 +454,15 @@ function Q_Angle(q1,q2) = // q = QuatXYZ([45,35,10]); // pts = Qrot(q, p=[[2,3,4], [4,5,6], [9,2,3]]); module Qrot(q) { - multmatrix(Q_Matrix4(q)) { + multmatrix(Q_to_matrix4(q)) { children(); } } function Qrot(q,p) = - is_undef(p)? Q_Matrix4(q) : + is_undef(p)? Q_to_matrix4(q) : is_vector(p)? Qrot(q,[p])[0] : - apply(Q_Matrix4(q), p); + apply(Q_to_matrix4(q), p); // Module: Qrot_copies() @@ -446,29 +482,6 @@ 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); @@ -518,11 +531,11 @@ function Q_Rotation_path(q1, n=1, q2) = 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] + ? [for( i=0, dR=Q_to_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]; + dq = Q_pow( _Qmul( q2, Q_Inverse(q1) ), 1/n ), + dR = Q_to_matrix4(dq) ) + [for( i=0, R=Q_to_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)) @@ -565,8 +578,8 @@ function Q_Nlerp(q1,q2,u) = 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 ) ]; + ? _Qunit((1-u)*q1 + u*q2 ) + : [for (ui=u) _Qunit((1-ui)*q1 + ui*q2 ) ]; // Function: Q_Squad() @@ -616,6 +629,17 @@ function Q_Squad(q1,q2,q3,q4,u) = : [for(ui=u) Q_Slerp( Q_Slerp(q1,q4,ui), Q_Slerp(q2,q3,ui), 2*ui*(1-ui) ) ]; +function Q_Scubic(q1,q2,q3,q4,u) = + assert(is_finite(u) || is_range(u) || is_vector(u) , + "Invalid interpolation coefficient(s)" ) + is_num(u) + ? let( q12 = Q_Slerp(q1,q2,u), + q23 = Q_Slerp(q2,q3,u), + q34 = Q_Slerp(q3,q4,u) ) + Q_Slerp(Q_Slerp(q12,q23,u),Q_Slerp(q23,q34,u),u) + : [for(ui=u) Q_Scubic( q1,q2,q3,q4,ui) ]; + + // Function: Q_exp() // Usage: // q2 = Q_exp(q); diff --git a/tests/test_all.scad b/tests/test_all.scad new file mode 100644 index 0000000..fe6512b --- /dev/null +++ b/tests/test_all.scad @@ -0,0 +1,28 @@ +//include +include +include +include +include +include +include +include +include +include +include +include +include +include +//include +include +include +include +include +include +include +include +include +include +include +include + + diff --git a/tests/test_quaternions.scad b/tests/test_quaternions.scad index b1e5130..af19408 100644 --- a/tests/test_quaternions.scad +++ b/tests/test_quaternions.scad @@ -1,6 +1,80 @@ include <../std.scad> include <../strings.scad> +function plane_fit(points) = + assert( is_matrix(points,undef,3) , "Improper point list." ) + len(points)< 2 ? [] : + len(points)==3 ? plane3pt(points[0],points[1],points[2]) : + let( + A = [for(pi=points) concat(pi,-1) ], + qr = qr_factor(A), + R = select(qr[1],0,3), +x=[for(ri=R) echo(R=ri)0], +//y=[for(qi=qr[0]) echo(Q=qi)0], + s = back_substitute + ( R, + sum([for(qi=qr[0])[for(j=[0:3]) qi[j] ] ]) + ) +//, s0 = echo(ls=linear_solve(A,[for(p=points) 0])) +,w=s==[]? echo("s == []")0:echo(s=s)0 + ) + s==[] + ? // points are collinear + [] + : // plane through the origin? + approx(norm([s.x, s.y, s.z]), 0) + ? let( + k = max_index([for(i=[0:2]) abs(s[i]) ]), + A = [[1,0,0], for(pi=points) pi ], + b = [1, for(i=[1:len(points)]) 0], + s = linear_solve(A,b) +, y=echo(s2=s) + ) + s==[]? []: + concat(s, 0)/norm(s) + : s/norm([s.x, s.y, s.z]) ; + +function plane_fit(points) = + assert( is_matrix(points,undef,3) , "Improper point list." ) + len(points)< 2 ? [] : + len(points)==3 ? plane3pt(points[0],points[1],points[2]) : + let( + A = [for(pi=points) concat(pi,-1) ], + B = [ for(pi=points) [0,1] ], + sB = transpose(linear_solve(A,B)), + s = sB==[] ? [] + : norm(sB[1])N/2? 10*[1,1,1]: rands(-1,1,3) ]; +pf = plane_fit(pts); +echo(pf); +//pm = concat(sum(pts)/len(pts), -1); +//echo(pmfit=pm*pf); function rec_cmp(a,b,eps=1e-9) = typeof(a)!=typeof(b)? false : @@ -8,7 +82,18 @@ 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; +function standardize(v) = + v==[] + ? [] + : sign(first_nonzero(v))*v; + +function first_nonzero(v) = + v==[] ? 0 + : is_num(v) ? v + : [for(vi=v) if(!is_list(vi) && ! approx(vi,0)) vi + else if(is_list(vi)) first_nonzero(vi), 0 ][0]; + +module assert_std(vc,ve) { assert_approx(standardize(vc),standardize(ve)); } module verify_f(actual,expected) { if (!rec_cmp(actual,expected)) { @@ -25,7 +110,7 @@ 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([0,0,0,0]),true); verify_f(Q_is_quat([1,0,2,0]),true); verify_f(Q_is_quat([1,0,2,0,0]),false); } @@ -104,7 +189,7 @@ 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(Q_to_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])); } @@ -200,7 +285,7 @@ test_Q_Mul(); module test_Q_Cumulative() { verify_f(Q_Cumulative([QuatZ(30),QuatX(57),QuatY(18)]),[[0, 0, 0.2588190451, 0.9659258263], [0.4608999698, -0.1234977747, 0.2274546059, 0.8488721457], [0.4908072659, 0.01081554785, 0.1525536221, 0.8577404293]]); } -test_Q_Cumulative(); +*test_Q_Cumulative(); module test_Q_Dot() { @@ -219,18 +304,18 @@ test_Q_Neg(); module test_Q_Conj() { + ang = rands(0,360,3); verify_f(Q_Conj([1,0,0,1]),[-1,0,0,1]); verify_f(Q_Conj([0,1,1,0]),[0,-1,-1,0]); verify_f(Q_Conj(QuatXYZ([23,45,67])),[0.0533818345, -0.4143703268, -0.4360652669, 0.7970537592]); + verify_f(Q_Mul(Q_Conj(QuatXYZ(ang)),QuatXYZ(ang)),Q_Ident()); } 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_Inverse([1,0,0,1]),[-1,0,0,1]); + verify_f(Q_Inverse([0,1,1,0]),[0,-1,-1,0]); verify_f(Q_Mul(Q_Inverse(QuatXYZ([23,45,67])),QuatXYZ([23,45,67])),Q_Ident()); } test_Q_Inverse(); @@ -247,7 +332,7 @@ test_Q_Norm(); module test_Q_Normalize() { verify_f(Q_Normalize([1,0,0,1]),[0.7071067812, 0, 0, 0.7071067812]); verify_f(Q_Normalize([0,1,1,0]),[0, 0.7071067812, 0.7071067812, 0]); - verify_f(Q_Normalize(QuatXYZ([23,45,67])),[-0.0533818345, 0.4143703268, 0.4360652669, 0.7970537592]); +// verify_f(Q_Normalize(QuatXYZ([23,45,67])),[-0.0533818345, 0.4143703268, 0.4360652669, 0.7970537592]); } test_Q_Normalize(); @@ -260,28 +345,47 @@ test_Q_Dist(); module test_Q_Slerp() { - verify_f(Q_Slerp(QuatX(45),QuatY(30),0.0),QuatX(45)); - verify_f(Q_Slerp(QuatX(45),QuatY(30),0.5),[0.1967063121, 0.1330377423, 0, 0.9713946602]); - verify_f(Q_Slerp(QuatX(45),QuatY(30),1.0),QuatY(30)); + u = rands(0,1,1)[0]; + ul = rands(0,1,5); + ul2 = [1,1,1,1,1]-ul; + a1 = rands(0,360,1)[0]; + a2 = rands(0,360,1)[0]; + a3 = rands(0,360,1)[0]; + verify_f(standardize(Q_Slerp(QuatX(a1),QuatY(a2),0.0)), standardize(QuatX(a1))); + verify_f(standardize(Q_Slerp(QuatX(45),QuatY(30),0.5)), + [0.1967063121, 0.1330377423, 0, 0.9713946602]); + verify_f(standardize(Q_Slerp(QuatX(a1),QuatY(a2),1.0)), standardize(QuatY(a2))); + verify_f(standardize(Q_Slerp(QuatXYZ([a1,a2,0]),-QuatY(a2),u)), + standardize(Q_Slerp(QuatXYZ([a1,a2,0]), QuatY(a2),u))); + verify_f(standardize(Q_Slerp(QuatX(a1),QuatX(a1),u)), standardize(QuatX(a1))); + verify_f(standardize(Q_Slerp(QuatXYZ([a1,a2,0]),QuatXYZ([a2,0,a1]),u)), + standardize(Q_Slerp(QuatXYZ([a2,0,a1]),QuatXYZ([a1,a2,0]),1-u))); + verify_f(standardize(Q_Slerp(QuatXYZ([a1,a2,0]),QuatXYZ([a2,0,a1]),ul)), + standardize(Q_Slerp(QuatXYZ([a2,0,a1]),QuatXYZ([a1,a2,0]),ul2))); } test_Q_Slerp(); -module test_Q_Matrix3() { - verify_f(Q_Matrix3(QuatZ(37)),rot(37,planar=true)); - verify_f(Q_Matrix3(QuatZ(-49)),rot(-49,planar=true)); +module test_Q_to_matrix3() { + rotZ_37 = rot(37); + rotZ_37_3 = [for(i=[0:2]) [for(j=[0:2]) rotZ_37[i][j] ] ]; + angs = [12,-49,40]; + rot4 = rot(angs); + rot3 = [for(i=[0:2]) [for(j=[0:2]) rot4[i][j] ] ]; + verify_f(Q_to_matrix3(QuatZ(37)),rotZ_37_3); + verify_f(Q_to_matrix3(QuatXYZ(angs)),rot3); } -test_Q_Matrix3(); +test_Q_to_matrix3(); -module test_Q_Matrix4() { - verify_f(Q_Matrix4(QuatZ(37)),rot(37)); - verify_f(Q_Matrix4(QuatZ(-49)),rot(-49)); - verify_f(Q_Matrix4(QuatX(37)),rot([37,0,0])); - verify_f(Q_Matrix4(QuatY(37)),rot([0,37,0])); - verify_f(Q_Matrix4(QuatXYZ([12,34,56])),rot([12,34,56])); +module test_Q_to_matrix4() { + verify_f(Q_to_matrix4(QuatZ(37)),rot(37)); + verify_f(Q_to_matrix4(QuatZ(-49)),rot(-49)); + verify_f(Q_to_matrix4(QuatX(37)),rot([37,0,0])); + verify_f(Q_to_matrix4(QuatY(37)),rot([0,37,0])); + verify_f(Q_to_matrix4(QuatXYZ([12,34,56])),rot([12,34,56])); } -test_Q_Matrix4(); +test_Q_to_matrix4(); module test_Q_Axis() { @@ -321,23 +425,23 @@ 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]))); +module test_Q_from_matrix() { + verify_f(standardize(Q_from_matrix(Q_to_matrix3(Quat([12,34,56],33)))),standardize(Quat([12,34,56],33))); + verify_f(Q_to_matrix3(Q_from_matrix(Q_to_matrix3(QuatXYZ([12,34,56])))), + Q_to_matrix3(QuatXYZ([12,34,56]))); } -test_Q_Rotation(); +test_Q_from_matrix(); 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), 5, QuatY(13.5))[0] , Q_to_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)[3] , xrot(11+(55-11)*3/4)); verify_f(Q_Rotation_path(QuatX(11), 5)[4] , xrot(55)); } @@ -347,48 +451,54 @@ 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_Rotation_path(QuatX(135), 16, QuatY(13.5))[8] , + Q_to_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)); + u = rands(0,1,5); + su = [1,1,1,1,1]-u; + verify_f(Q_Squad(QuatX(45),QuatZ(30),QuatX(32),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), + 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)); + verify_f( Q_Squad(QuatY(10),QuatZ(20),QuatX(32),QuatXYZ([210,120,30]),u[0] ), + Q_Squad(QuatXYZ([210,120,30]),QuatX(32),QuatZ(20),QuatY(10),1-u[0] ) ); + verify_f( Q_Squad(QuatY(10),QuatZ(20),QuatX(32),QuatXYZ([210,120,30]),u ), + Q_Squad(QuatXYZ([210,120,30]),QuatX(32),QuatZ(20),QuatY(10),su ) ); } test_Q_Squad(); module test_Q_exp() { + q=QuatXYZ(rands(0,360,3)); 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)))); + verify_f(Q_exp(Q_ln([1,2,3,4])), [1,2,3,4]); + verify_f(Q_exp(Q_ln(q)), q); + 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() { + q = QuatXYZ(rands(0,360,3)); 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]); + verify_f(Q_ln(Q_exp(q)), q); + verify_f(Q_ln(q)+Q_ln(Q_Conj(q)), [0,0,0,0]); } test_Q_ln(); module test_Q_pow() { - q = Quat([1,2,3],77); + q = QuatXYZ(rands(0,360,3)); verify_f(Q_pow(q,1), q); verify_f(Q_pow(q,0), Q_Ident()); verify_f(Q_pow(q,-1), Q_Inverse(q)); From 069432021b420aab2ce91ce1bfaf3cf73c8cef28 Mon Sep 17 00:00:00 2001 From: RonaldoCMP Date: Wed, 26 Aug 2020 14:09:32 +0100 Subject: [PATCH 04/23] Revert "Revert "Revert "Correction of C_times validation""" This reverts commit e8499baba185848ff6f6f1745181a33573ee5c58. --- arrays.scad | 125 +++-------------------- math.scad | 2 +- quaternions.scad | 170 ++++++++++++++----------------- tests/test_all.scad | 28 ------ tests/test_quaternions.scad | 196 ++++++++---------------------------- 5 files changed, 131 insertions(+), 390 deletions(-) delete mode 100644 tests/test_all.scad diff --git a/arrays.scad b/arrays.scad index 3ae7a88..835fd65 100644 --- a/arrays.scad +++ b/arrays.scad @@ -708,108 +708,36 @@ function _sort_vectors4(arr) = y ] ) concat( _sort_vectors4(lesser), equal, _sort_vectors4(greater) ); -// sort a list of vectors -function _sort_vectors(arr, _i=0) = - len(arr)<=1 || _i>=len(arr[0]) ? arr : - let( - pivot = arr[floor(len(arr)/2)][_i], - lesser = [ for (entry=arr) if (entry[_i] < pivot ) entry ], - equal = [ for (entry=arr) if (entry[_i] == pivot ) entry ], - greater = [ for (entry=arr) if (entry[_i] > pivot ) entry ] - ) - concat( - _sort_vectors(lesser, _i ), - _sort_vectors(equal, _i+1 ), - _sort_vectors(greater, _i ) ); - -// given pairs of an index and a vector, return the list of indices of the list sorted by the vectors -function _sort_vectors_indexed(arr, _i=0) = - arr==[] ? [] : - len(arr)==1 || _i>=len(arr[0][1]) ? [for(ai=arr) ai[0]] : - let( - pivot = arr[floor(len(arr)/2)][1][_i], - lesser = [ for (entry=arr) if (entry[1][_i] < pivot ) entry ], - equal = [ for (entry=arr) if (entry[1][_i] == pivot ) entry ], - greater = [ for (entry=arr) if (entry[1][_i] > pivot ) entry ] - ) - concat( - _sort_vectors_indexed(lesser, _i ), - _sort_vectors_indexed(equal, _i+1 ), - _sort_vectors_indexed(greater, _i ) ); - // when idx==undef, returns the sorted array // otherwise, returns the indices of the sorted array function _sort_general(arr, idx=undef) = - len(arr)<=1 ? arr : + (len(arr)<=1) ? arr : is_undef(idx) - ? _simple_sort(arr) -// : _lexical_sort(arr) + ? _sort_scalar(arr) : let( arrind=[for(k=[0:len(arr)-1], ark=[arr[k]]) [ k, [for (i=idx) ark[i]] ] ] ) - _indexed_sort(arrind ); - -// sort simple lists with compare_vals() -function _simple_sort(arr) = - arr==[] || len(arr)==1 ? arr : - let( - pivot = arr[floor(len(arr)/2)], - lesser = [ for (entry=arr) if (compare_vals(entry, pivot) <0 ) entry ], - equal = [ for (entry=arr) if (compare_vals(entry, pivot)==0 ) entry ], - greater = [ for (entry=arr) if (compare_vals(entry, pivot) >0 ) entry ] - ) - concat( - _simple_sort(lesser), - equal, - _simple_sort(greater) - ); - - + _indexed_sort(arrind); + // given a list of pairs, return the first element of each pair of the list sorted by the second element of the pair -// it uses compare_vals() -function _lexical_sort(arr, _i=0) = - arr==[] || len(arr)==1 || _i>=len(arr[0]) ? arr : +// the sorting is done using compare_vals() +function _indexed_sort(arrind) = + arrind==[] ? [] : len(arrind)==1? [arrind[0][0]] : + let( pivot = arrind[floor(len(arrind)/2)][1] ) let( - pivot = arr[floor(len(arr)/2)][_i], - lesser = [ for (entry=arr) if (compare_vals(entry[_i], pivot) <0 ) entry ], - equal = [ for (entry=arr) if (compare_vals(entry[_i], pivot)==0 ) entry ], - greater = [ for (entry=arr) if (compare_vals(entry[_i], pivot) >0 ) entry ] + lesser = [ for (entry=arrind) if (compare_vals(entry[1], pivot) <0 ) entry ], + equal = [ for (entry=arrind) if (compare_vals(entry[1], pivot)==0 ) entry[0] ], + greater = [ for (entry=arrind) if (compare_vals(entry[1], pivot) >0 ) entry ] ) - concat( - _lexical_sort(lesser, _i ), - _lexical_sort(equal, _i+1 ), - _lexical_sort(greater, _i ) - ); - - -// given a list of pairs, return the first element of each pair of the list sorted by the second element of the pair -// it uses compare_vals() -function _indexed_sort(arr, _i=0) = - arr==[] ? [] : - len(arr)==1? [arr[0][0]] : - _i>=len(arr[0][1]) ? [for(ai=arr) ai[0]] : - let( - pivot = arr[floor(len(arr)/2)][1][_i], - lesser = [ for (entry=arr) if (compare_vals(entry[1][_i], pivot) <0 ) entry ], - equal = [ for (entry=arr) if (compare_vals(entry[1][_i], pivot)==0 ) entry ], - greater = [ for (entry=arr) if (compare_vals(entry[1][_i], pivot) >0 ) entry ] - ) - concat( - _indexed_sort(lesser, _i ), - _indexed_sort(equal, _i+1 ), - _indexed_sort(greater, _i ) ); + concat(_indexed_sort(lesser), equal, _indexed_sort(greater)); // returns true for valid index specifications idx in the interval [imin, imax) // note that idx can't have any value greater or EQUAL to imax -// this allows imax=INF as a bound to numerical lists function _valid_idx(idx,imin,imax) = is_undef(idx) || ( is_finite(idx) && idx>=imin && idx< imax ) || ( is_list(idx) && min(idx)>=imin && max(idx)< imax ) - || ( ! is_list(idx) // it implicitly is a range - && (idx[1]>0 && idx[0]>=imin && idx[2]< imax) - || - (idx[0]=imin) ); + || ( valid_range(idx) && idx[0]>=imin && idx[2]< imax ); // Function: sort() @@ -830,7 +758,7 @@ function _valid_idx(idx,imin,imax) = function sort(list, idx=undef) = !is_list(list) || len(list)<=1 ? list : is_def(idx) - ? assert( _valid_idx(idx,0,len(list[0])) , "Invalid indices or out of range.") + ? assert( _valid_idx(idx,0,len(list)) , "Invalid indices.") let( sarr = _sort_general(list,idx) ) [for(i=[0:len(sarr)-1]) list[sarr[i]] ] : let(size = array_dim(list)) @@ -845,31 +773,6 @@ function sort(list, idx=undef) = ) : _sort_general(list); -function sort(list, idx=undef) = - !is_list(list) || len(list)<=1 ? list : - is_vector(list) - ? assert( _valid_idx(idx,0,len(list[0])) , str("Invalid indices or out of range. ",list)) - is_def(idx) - ? sort_vector_indexed([for(i=[0:len(list)-1]) [i, [list(i)] ] ]) - : sort_scalar(list) - : is_matrix(list) - ? list==[] || list[0]==[] ? list : - assert( _valid_idx(idx,0,len(list[0])) , "Invalid indices or out of range.") - is_def(idx) - ? sort_vector_indexed([for(i=[0:len(list)-1], li=[list[i]]) [i, [for(ind=idx) li(ind)] ] ]) - : sort_vector(list) - : list==[] || list[0]==[] ? list : - let( llen = [for(li=list) !is_list(li) || is_string(li) ? 0: len(li)], - m = min(llen), - M = max(llen) - ) - M==0 ? _simple_sort(list) : - assert( m>0 && _valid_idx(idx,m-1,M) , "Invalid indices or out of range.") - is_def(idx) - ? _sort_general(list,idx) - : let( ils = _sort_general(list, [m:M]) ) - [for(i=[0:len(list)-1]) list[ils[i]] ]; - // Function: sortidx() // Description: diff --git a/math.scad b/math.scad index e5c4625..888d048 100644 --- a/math.scad +++ b/math.scad @@ -1175,7 +1175,7 @@ function deriv3(data, h=1, closed=false) = // Description: // Multiplies two complex numbers represented by 2D vectors. function C_times(z1,z2) = - assert( is_matrix([z1,z2],2,2), "Complex numbers should be represented by 2D vectors" ) + assert( is_vector(z1+z2,2), "Complex numbers should be represented by 2D vectors." ) [ z1.x*z2.x - z1.y*z2.y, z1.x*z2.y + z1.y*z2.x ]; // Function: C_div() diff --git a/quaternions.scad b/quaternions.scad index 44b0e29..e21801c 100644 --- a/quaternions.scad +++ b/quaternions.scad @@ -27,7 +27,7 @@ function _Qreal(q) = q[3]; function _Qset(v,r) = concat( v, r ); // normalizes without checking -function _Qunit(q) = q/norm(q); +function _Qnorm(q) = q/norm(q); // Function: Q_is_quat() @@ -36,7 +36,7 @@ function _Qunit(q) = q/norm(q); // 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 Q_is_quat(q) = is_vector(q,4) && ! approx(norm(q),0) ; // Function: Quat() @@ -101,7 +101,7 @@ function QuatXYZ(a=[0,0,0]) = qy = QuatY(a[1]), qz = QuatZ(a[2]) ) - _Qmul(qz, _Qmul(qy, qx)); + Q_Mul(qz, Q_Mul(qy, qx)); // Function: Q_From_to() @@ -202,36 +202,32 @@ function Q_Sub(a, 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)") - _Qmul(a,b); - -function _Qmul(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 + a[3]*b[3] - a.x*b.x - a.y*b.y - a.z*b.z, ]; -// [ [a[3], -a.z, a.y, a.x], -// [ a.z, a[3], -a.x, a.y], -// [-a.y, a.x, a[3], a.z], -// [-a.x, -a.y, -a.z, a[3]] ]*[b.x,b.y,b.z,b[3]] - + // Function: Q_Cumulative() // Usage: -// Q_Cumulative(ql); +// Q_Cumulative(v); // Description: // 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(ql) = - assert( is_matrix(ql,undef,4) && len(ql)>0, "Invalid list of quaternions." ) - [for( i = 0, q = ql[0]; - i<=len(ql); - i = i+1, q = (i==len(ql))? 0: _Qmul(q,ql[i]) ) - q ]; +function Q_Cumulative(v, _i=0, _acc=[]) = + _i==len(v) ? _acc : + Q_Cumulative( + v, _i+1, + concat( + _acc, + [_i==0 ? v[_i] : Q_Mul(v[_i], select(_acc,-1))] + ) + ); // Function: Q_Dot() @@ -256,7 +252,7 @@ function Q_Neg(q) = // Q_Conj(q) // Description: Returns the conjugate of quaternion `q`. function Q_Conj(q) = - assert( Q_is_quat(q), str("Invalid quaternion",q) ) + assert( Q_is_quat(q), "Invalid quaternion" ) [-q.x, -q.y, -q.z, q[3]]; @@ -266,7 +262,7 @@ function Q_Conj(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 = q/norm(q) ) + let(q = _Qnorm(q) ) [-q.x, -q.y, -q.z, q[3]]; @@ -287,9 +283,7 @@ function Q_Norm(q) = // Description: Normalizes quaternion `q`, so that norm([W,X,Y,Z]) == 1. function Q_Normalize(q) = assert( Q_is_quat(q) , "Invalid quaternion" ) - approx(_Qvec(q), [0,0,0]) - ? Q_Ident() - : q/norm(q); + q/norm(q); // Function: Q_Dist() @@ -324,32 +318,31 @@ function Q_Dist(q1, q2) = // Qrot(q) right(80) cube([10,10,1]); // #sphere(r=80); function Q_Slerp(q1, q2, u, _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), - q3 = dot>0.9995? q2: _Qunit(q2 - dot*q1) - ) - ! is_num(u) - ? [for (uu=u) _Qslerp(q1, q3, uu, dot)] - : _Qslerp(q1, q3, 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 _Qslerp(q1, q2, u, dot) = - dot>0.9995 - ? _Qunit(q1 + u*(q2-q1)) - : let( theta = u*acos(dot) ) - _Qunit(q1*cos(theta) + q2*sin(theta)); - -// Function: Q_to_matrix3() +// Function: Q_Matrix3() // Usage: -// Q_to_matrix3(q); +// Q_Matrix3(q); // Description: // Returns the 3x3 rotation matrix for the given normalized quaternion q. -function Q_to_matrix3(q) = +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]], @@ -358,12 +351,12 @@ function Q_to_matrix3(q) = ]; -// Function: Q_to_matrix4() +// Function: Q_Matrix4() // Usage: -// Q_to_matrix4(q); +// Q_Matrix4(q); // Description: // Returns the 4x4 rotation matrix for the given normalized quaternion q. -function Q_to_matrix4(q) = +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], @@ -373,35 +366,6 @@ function Q_to_matrix4(q) = ]; -// Function: Q_from_matrix() -// Usage: -// Q_from_matrix(M) -// Description: -// Returns a normalized quaternion corresponding to the rotation matrix M. -// M may be a 3x3 rotation matrix or a homogeneous 4x4 rotation matrix. -// The last row and last column of M are ignored for 4x4 matrices. -// It doesn't check whether M is in fact a rotation matrix. -// If M is not a rotation, the returned quaternion is unpredictable. -// -// based on https://en.wikipedia.org/wiki/Rotation_matrix -// -function Q_from_matrix(M) = - assert( is_matrix(M) && (len(M)==3 || len(M)==4) && (len(M[0])==3 || len(M[0])==4), - "The matrix should be 3x3 or 4x4") - let( tr = M[0][0]+M[1][1]+M[2][2] ) // M trace - tr>0 - ? let( r = sqrt(1+tr), s = -1/r/2 ) - _Qunit( _Qset([ M[1][2]-M[2][1], M[2][0]-M[0][2], M[0][1]-M[1][0] ]*s, r/2 ) ) - : let( - i = max_index([ M[0][0], M[1][1], M[2][2] ]), - r = sqrt(1 + 2*M[i][i] -M[0][0] -M[1][1] -M[2][2] ), - s = 1/r/2 - ) - i==0 ? _Qunit( _Qset( [ r/2, s*(M[1][0]+M[0][1]), s*(M[0][2]+M[2][0]) ], s*(M[2][1]-M[1][2])) ): - i==1 ? _Qunit( _Qset( [ s*(M[1][0]+M[0][1]), r/2, s*(M[2][1]+M[1][2]) ], s*(M[0][2]-M[2][0])) ) - : _Qunit( _Qset( [ s*(M[2][0]+M[0][2]), s*(M[1][2]+M[2][1]), r/2 ], s*(M[1][0]-M[0][1])) ) ; - - // Function: Q_Axis() // Usage: // Q_Axis(q) @@ -454,15 +418,15 @@ function Q_Angle(q1,q2) = // q = QuatXYZ([45,35,10]); // pts = Qrot(q, p=[[2,3,4], [4,5,6], [9,2,3]]); module Qrot(q) { - multmatrix(Q_to_matrix4(q)) { + multmatrix(Q_Matrix4(q)) { children(); } } function Qrot(q,p) = - is_undef(p)? Q_to_matrix4(q) : + is_undef(p)? Q_Matrix4(q) : is_vector(p)? Qrot(q,[p])[0] : - apply(Q_to_matrix4(q), p); + apply(Q_Matrix4(q), p); // Module: Qrot_copies() @@ -482,6 +446,29 @@ 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); @@ -531,11 +518,11 @@ function Q_Rotation_path(q1, n=1, q2) = 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_to_matrix4(q1), R=dR; i<=n; i=i+1, R=dR*R ) R] + ? [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( _Qmul( q2, Q_Inverse(q1) ), 1/n ), - dR = Q_to_matrix4(dq) ) - [for( i=0, R=Q_to_matrix4(q1); i<=n; i=i+1, R=dR*R ) R]; + 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)) @@ -578,8 +565,8 @@ function Q_Nlerp(q1,q2,u) = let( q1 = Q_Normalize(q1), q2 = Q_Normalize(q2) ) is_num(u) - ? _Qunit((1-u)*q1 + u*q2 ) - : [for (ui=u) _Qunit((1-ui)*q1 + ui*q2 ) ]; + ? _Qnorm((1-u)*q1 + u*q2 ) + : [for (ui=u) _Qnorm((1-ui)*q1 + ui*q2 ) ]; // Function: Q_Squad() @@ -629,17 +616,6 @@ function Q_Squad(q1,q2,q3,q4,u) = : [for(ui=u) Q_Slerp( Q_Slerp(q1,q4,ui), Q_Slerp(q2,q3,ui), 2*ui*(1-ui) ) ]; -function Q_Scubic(q1,q2,q3,q4,u) = - assert(is_finite(u) || is_range(u) || is_vector(u) , - "Invalid interpolation coefficient(s)" ) - is_num(u) - ? let( q12 = Q_Slerp(q1,q2,u), - q23 = Q_Slerp(q2,q3,u), - q34 = Q_Slerp(q3,q4,u) ) - Q_Slerp(Q_Slerp(q12,q23,u),Q_Slerp(q23,q34,u),u) - : [for(ui=u) Q_Scubic( q1,q2,q3,q4,ui) ]; - - // Function: Q_exp() // Usage: // q2 = Q_exp(q); diff --git a/tests/test_all.scad b/tests/test_all.scad deleted file mode 100644 index fe6512b..0000000 --- a/tests/test_all.scad +++ /dev/null @@ -1,28 +0,0 @@ -//include -include -include -include -include -include -include -include -include -include -include -include -include -include -//include -include -include -include -include -include -include -include -include -include -include -include - - diff --git a/tests/test_quaternions.scad b/tests/test_quaternions.scad index af19408..b1e5130 100644 --- a/tests/test_quaternions.scad +++ b/tests/test_quaternions.scad @@ -1,80 +1,6 @@ include <../std.scad> include <../strings.scad> -function plane_fit(points) = - assert( is_matrix(points,undef,3) , "Improper point list." ) - len(points)< 2 ? [] : - len(points)==3 ? plane3pt(points[0],points[1],points[2]) : - let( - A = [for(pi=points) concat(pi,-1) ], - qr = qr_factor(A), - R = select(qr[1],0,3), -x=[for(ri=R) echo(R=ri)0], -//y=[for(qi=qr[0]) echo(Q=qi)0], - s = back_substitute - ( R, - sum([for(qi=qr[0])[for(j=[0:3]) qi[j] ] ]) - ) -//, s0 = echo(ls=linear_solve(A,[for(p=points) 0])) -,w=s==[]? echo("s == []")0:echo(s=s)0 - ) - s==[] - ? // points are collinear - [] - : // plane through the origin? - approx(norm([s.x, s.y, s.z]), 0) - ? let( - k = max_index([for(i=[0:2]) abs(s[i]) ]), - A = [[1,0,0], for(pi=points) pi ], - b = [1, for(i=[1:len(points)]) 0], - s = linear_solve(A,b) -, y=echo(s2=s) - ) - s==[]? []: - concat(s, 0)/norm(s) - : s/norm([s.x, s.y, s.z]) ; - -function plane_fit(points) = - assert( is_matrix(points,undef,3) , "Improper point list." ) - len(points)< 2 ? [] : - len(points)==3 ? plane3pt(points[0],points[1],points[2]) : - let( - A = [for(pi=points) concat(pi,-1) ], - B = [ for(pi=points) [0,1] ], - sB = transpose(linear_solve(A,B)), - s = sB==[] ? [] - : norm(sB[1])N/2? 10*[1,1,1]: rands(-1,1,3) ]; -pf = plane_fit(pts); -echo(pf); -//pm = concat(sum(pts)/len(pts), -1); -//echo(pmfit=pm*pf); function rec_cmp(a,b,eps=1e-9) = typeof(a)!=typeof(b)? false : @@ -82,18 +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 standardize(v) = - v==[] - ? [] - : sign(first_nonzero(v))*v; - -function first_nonzero(v) = - v==[] ? 0 - : is_num(v) ? v - : [for(vi=v) if(!is_list(vi) && ! approx(vi,0)) vi - else if(is_list(vi)) first_nonzero(vi), 0 ][0]; - -module assert_std(vc,ve) { assert_approx(standardize(vc),standardize(ve)); } +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)) { @@ -110,7 +25,7 @@ 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]),true); + 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); } @@ -189,7 +104,7 @@ 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_to_matrix4(Q_From_to([1,2,3], [4,5,2])), rot(from=[1,2,3],to=[4,5,2])); + 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])); } @@ -285,7 +200,7 @@ test_Q_Mul(); module test_Q_Cumulative() { verify_f(Q_Cumulative([QuatZ(30),QuatX(57),QuatY(18)]),[[0, 0, 0.2588190451, 0.9659258263], [0.4608999698, -0.1234977747, 0.2274546059, 0.8488721457], [0.4908072659, 0.01081554785, 0.1525536221, 0.8577404293]]); } -*test_Q_Cumulative(); +test_Q_Cumulative(); module test_Q_Dot() { @@ -304,18 +219,18 @@ test_Q_Neg(); module test_Q_Conj() { - ang = rands(0,360,3); verify_f(Q_Conj([1,0,0,1]),[-1,0,0,1]); verify_f(Q_Conj([0,1,1,0]),[0,-1,-1,0]); verify_f(Q_Conj(QuatXYZ([23,45,67])),[0.0533818345, -0.4143703268, -0.4360652669, 0.7970537592]); - verify_f(Q_Mul(Q_Conj(QuatXYZ(ang)),QuatXYZ(ang)),Q_Ident()); } test_Q_Conj(); module test_Q_Inverse() { - verify_f(Q_Inverse([1,0,0,1]),[-1,0,0,1]); - verify_f(Q_Inverse([0,1,1,0]),[0,-1,-1,0]); + + 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(); @@ -332,7 +247,7 @@ test_Q_Norm(); module test_Q_Normalize() { verify_f(Q_Normalize([1,0,0,1]),[0.7071067812, 0, 0, 0.7071067812]); verify_f(Q_Normalize([0,1,1,0]),[0, 0.7071067812, 0.7071067812, 0]); -// verify_f(Q_Normalize(QuatXYZ([23,45,67])),[-0.0533818345, 0.4143703268, 0.4360652669, 0.7970537592]); + verify_f(Q_Normalize(QuatXYZ([23,45,67])),[-0.0533818345, 0.4143703268, 0.4360652669, 0.7970537592]); } test_Q_Normalize(); @@ -345,47 +260,28 @@ test_Q_Dist(); module test_Q_Slerp() { - u = rands(0,1,1)[0]; - ul = rands(0,1,5); - ul2 = [1,1,1,1,1]-ul; - a1 = rands(0,360,1)[0]; - a2 = rands(0,360,1)[0]; - a3 = rands(0,360,1)[0]; - verify_f(standardize(Q_Slerp(QuatX(a1),QuatY(a2),0.0)), standardize(QuatX(a1))); - verify_f(standardize(Q_Slerp(QuatX(45),QuatY(30),0.5)), - [0.1967063121, 0.1330377423, 0, 0.9713946602]); - verify_f(standardize(Q_Slerp(QuatX(a1),QuatY(a2),1.0)), standardize(QuatY(a2))); - verify_f(standardize(Q_Slerp(QuatXYZ([a1,a2,0]),-QuatY(a2),u)), - standardize(Q_Slerp(QuatXYZ([a1,a2,0]), QuatY(a2),u))); - verify_f(standardize(Q_Slerp(QuatX(a1),QuatX(a1),u)), standardize(QuatX(a1))); - verify_f(standardize(Q_Slerp(QuatXYZ([a1,a2,0]),QuatXYZ([a2,0,a1]),u)), - standardize(Q_Slerp(QuatXYZ([a2,0,a1]),QuatXYZ([a1,a2,0]),1-u))); - verify_f(standardize(Q_Slerp(QuatXYZ([a1,a2,0]),QuatXYZ([a2,0,a1]),ul)), - standardize(Q_Slerp(QuatXYZ([a2,0,a1]),QuatXYZ([a1,a2,0]),ul2))); + verify_f(Q_Slerp(QuatX(45),QuatY(30),0.0),QuatX(45)); + verify_f(Q_Slerp(QuatX(45),QuatY(30),0.5),[0.1967063121, 0.1330377423, 0, 0.9713946602]); + verify_f(Q_Slerp(QuatX(45),QuatY(30),1.0),QuatY(30)); } test_Q_Slerp(); -module test_Q_to_matrix3() { - rotZ_37 = rot(37); - rotZ_37_3 = [for(i=[0:2]) [for(j=[0:2]) rotZ_37[i][j] ] ]; - angs = [12,-49,40]; - rot4 = rot(angs); - rot3 = [for(i=[0:2]) [for(j=[0:2]) rot4[i][j] ] ]; - verify_f(Q_to_matrix3(QuatZ(37)),rotZ_37_3); - verify_f(Q_to_matrix3(QuatXYZ(angs)),rot3); +module test_Q_Matrix3() { + verify_f(Q_Matrix3(QuatZ(37)),rot(37,planar=true)); + verify_f(Q_Matrix3(QuatZ(-49)),rot(-49,planar=true)); } -test_Q_to_matrix3(); +test_Q_Matrix3(); -module test_Q_to_matrix4() { - verify_f(Q_to_matrix4(QuatZ(37)),rot(37)); - verify_f(Q_to_matrix4(QuatZ(-49)),rot(-49)); - verify_f(Q_to_matrix4(QuatX(37)),rot([37,0,0])); - verify_f(Q_to_matrix4(QuatY(37)),rot([0,37,0])); - verify_f(Q_to_matrix4(QuatXYZ([12,34,56])),rot([12,34,56])); +module test_Q_Matrix4() { + verify_f(Q_Matrix4(QuatZ(37)),rot(37)); + verify_f(Q_Matrix4(QuatZ(-49)),rot(-49)); + verify_f(Q_Matrix4(QuatX(37)),rot([37,0,0])); + verify_f(Q_Matrix4(QuatY(37)),rot([0,37,0])); + verify_f(Q_Matrix4(QuatXYZ([12,34,56])),rot([12,34,56])); } -test_Q_to_matrix4(); +test_Q_Matrix4(); module test_Q_Axis() { @@ -425,23 +321,23 @@ module test_Qrot() { test_Qrot(); -module test_Q_from_matrix() { - verify_f(standardize(Q_from_matrix(Q_to_matrix3(Quat([12,34,56],33)))),standardize(Quat([12,34,56],33))); - verify_f(Q_to_matrix3(Q_from_matrix(Q_to_matrix3(QuatXYZ([12,34,56])))), - Q_to_matrix3(QuatXYZ([12,34,56]))); +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_from_matrix(); +test_Q_Rotation(); module test_Q_Rotation_path() { - verify_f(Q_Rotation_path(QuatX(135), 5, QuatY(13.5))[0] , Q_to_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), 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)[3] , xrot(11+(55-11)*3/4)); verify_f(Q_Rotation_path(QuatX(11), 5)[4] , xrot(55)); } @@ -451,54 +347,48 @@ 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_to_matrix4(Q_Nlerp(QuatX(135), QuatY(13.5),0.5))); + 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() { - u = rands(0,1,5); - su = [1,1,1,1,1]-u; - verify_f(Q_Squad(QuatX(45),QuatZ(30),QuatX(32),QuatY(30),0.0),QuatX(45)); + 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), + 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(10),QuatZ(20),QuatX(32),QuatXYZ([210,120,30]),u[0] ), - Q_Squad(QuatXYZ([210,120,30]),QuatX(32),QuatZ(20),QuatY(10),1-u[0] ) ); - verify_f( Q_Squad(QuatY(10),QuatZ(20),QuatX(32),QuatXYZ([210,120,30]),u ), - Q_Squad(QuatXYZ([210,120,30]),QuatX(32),QuatZ(20),QuatY(10),su ) ); + 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() { - q=QuatXYZ(rands(0,360,3)); 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,4])), [1,2,3,4]); - verify_f(Q_exp(Q_ln(q)), q); - verify_f(Q_exp(q+Q_Inverse(q)),Q_Mul(Q_exp(q),Q_exp(Q_Inverse(q)))); + 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() { - q = QuatXYZ(rands(0,360,3)); 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(q)), q); - verify_f(Q_ln(q)+Q_ln(Q_Conj(q)), [0,0,0,0]); + 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 = QuatXYZ(rands(0,360,3)); + 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)); From 99b18d631195c8b7a58ca7b62b487fd197ebe0a3 Mon Sep 17 00:00:00 2001 From: RonaldoCMP Date: Wed, 26 Aug 2020 14:14:26 +0100 Subject: [PATCH 05/23] Correction of C_times arg validation --- math.scad | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/math.scad b/math.scad index 888d048..e5c4625 100644 --- a/math.scad +++ b/math.scad @@ -1175,7 +1175,7 @@ function deriv3(data, h=1, closed=false) = // Description: // Multiplies two complex numbers represented by 2D vectors. function C_times(z1,z2) = - assert( is_vector(z1+z2,2), "Complex numbers should be represented by 2D vectors." ) + assert( is_matrix([z1,z2],2,2), "Complex numbers should be represented by 2D vectors" ) [ z1.x*z2.x - z1.y*z2.y, z1.x*z2.y + z1.y*z2.x ]; // Function: C_div() From dba5aab9189e8aefae685a68b09564156555a8a3 Mon Sep 17 00:00:00 2001 From: Garth Minette Date: Wed, 26 Aug 2020 20:39:45 -0700 Subject: [PATCH 06/23] Implement Issue #2. Added diameter alternates for most radius options. --- beziers.scad | 14 ++++---- distributors.scad | 6 ++-- masks.scad | 87 ++++++++++++++++++++++++++++++----------------- mutators.scad | 2 +- paths.scad | 16 ++++++--- polyhedra.scad | 3 +- shapes.scad | 14 ++++---- version.scad | 2 +- 8 files changed, 89 insertions(+), 55 deletions(-) diff --git a/beziers.scad b/beziers.scad index 13a74ab..a8484b6 100644 --- a/beziers.scad +++ b/beziers.scad @@ -356,7 +356,7 @@ function bezier_segment_length(curve, start_u=0, end_u=1, max_deflect=0.01) = // Function: fillet3pts() // Usage: -// fillet3pts(p0, p1, p2, r); +// fillet3pts(p0, p1, p2, r|d); // Description: // Takes three points, defining two line segments, and works out the // cubic (degree 3) bezier segment (and surrounding control points) @@ -368,7 +368,8 @@ function bezier_segment_length(curve, start_u=0, end_u=1, max_deflect=0.01) = // p1 = The middle point. // p2 = The ending point. // r = The radius of the fillet/rounding. -// maxerr = Max amount bezier curve should diverge from actual radius curve. Default: 0.1 +// d = The diameter of the fillet/rounding. +// maxerr = Max amount bezier curve should diverge from actual curve. Default: 0.1 // Example(2D): // p0 = [40, 0]; // p1 = [0, 0]; @@ -376,7 +377,8 @@ function bezier_segment_length(curve, start_u=0, end_u=1, max_deflect=0.01) = // trace_polyline([p0,p1,p2], showpts=true, size=0.5, color="green"); // fbez = fillet3pts(p0,p1,p2, 10); // trace_bezier(slice(fbez, 1, -2), size=1); -function fillet3pts(p0, p1, p2, r, maxerr=0.1, w=0.5, dw=0.25) = let( +function fillet3pts(p0, p1, p2, r, d, maxerr=0.1, w=0.5, dw=0.25) = let( + r = get_radius(r=r,d=d), v0 = unit(p0-p1), v1 = unit(p2-p1), midv = unit((v0+v1)/2), @@ -391,8 +393,8 @@ function fillet3pts(p0, p1, p2, r, maxerr=0.1, w=0.5, dw=0.25) = let( bp = bezier_points([tp0, cp0, cp1, tp1], 0.5), tdist = norm(cp-bp) ) (abs(tdist-cpr) <= maxerr)? [tp0, tp0, cp0, cp1, tp1, tp1] : - (tdist0 && angle<90); diff --git a/mutators.scad b/mutators.scad index b801e48..4db899b 100644 --- a/mutators.scad +++ b/mutators.scad @@ -321,7 +321,7 @@ module chain_hull() // Usage: // cylindrical_extrude(size, ir|id, or|od, [convexity]) ... // Description: -// Cylindrically extrudes all 2D children, curved around a cylidrical shape. +// Extrudes all 2D children outwards, curved around a cylindrical shape. // Arguments: // or = The outer radius to extrude to. // od = The outer diameter to extrude to. diff --git a/paths.scad b/paths.scad index 93d82fa..60bc34e 100644 --- a/paths.scad +++ b/paths.scad @@ -418,7 +418,7 @@ function path_torsion(path, closed=false) = // scale = [X,Y] scaling factors for each axis. Default: `[1,1]` // Example(3D): // trace_polyline(path3d_spiral(turns=2.5, h=100, n=24, r=50), N=1, showpts=true); -function path3d_spiral(turns=3, h=100, n=12, r=undef, d=undef, cp=[0,0], scale=[1,1]) = let( +function path3d_spiral(turns=3, h=100, n=12, r, d, cp=[0,0], scale=[1,1]) = let( rr=get_radius(r=r, d=d, dflt=100), cnt=floor(turns*n), dz=h/cnt @@ -774,15 +774,19 @@ function assemble_path_fragments(fragments, eps=EPSILON, _finished=[]) = // Module: modulated_circle() +// Usage: +// modulated_circle(r|d, sines); // Description: // Creates a 2D polygon circle, modulated by one or more superimposed sine waves. // Arguments: -// r = radius of the base circle. +// r = Radius of the base circle. Default: 40 +// d = Diameter of the base circle. // sines = array of [amplitude, frequency] pairs, where the frequency is the number of times the cycle repeats around the circle. // Example(2D): // modulated_circle(r=40, sines=[[3, 11], [1, 31]], $fn=6); -module modulated_circle(r=40, sines=[10]) +module modulated_circle(r, sines=[10], d) { + r = get_radius(r=r, d=d, dflt=40); freqs = len(sines)>0? [for (i=sines) i[1]] : [5]; points = [ for (a = [0 : (360/segs(r)/max(freqs)) : 360]) @@ -829,7 +833,8 @@ module extrude_from_to(pt1, pt2, convexity=undef, twist=undef, scale=undef, slic // Arguments: // polyline = Array of points of a polyline path, to be extruded. // h = height of the spiral to extrude along. -// r = radius of the spiral to extrude along. +// r = Radius of the spiral to extrude along. Default: 50 +// d = Diameter of the spiral to extrude along. // twist = number of degrees of rotation to spiral up along height. // anchor = Translate so anchor point is at origin (0,0,0). See [anchor](attachments.scad#anchor). Default: `CENTER` // spin = Rotate this many degrees around the Z axis after anchor. See [spin](attachments.scad#spin). Default: `0` @@ -838,7 +843,8 @@ module extrude_from_to(pt1, pt2, convexity=undef, twist=undef, scale=undef, slic // Example: // poly = [[-10,0], [-3,-5], [3,-5], [10,0], [0,-30]]; // spiral_sweep(poly, h=200, r=50, twist=1080, $fn=36); -module spiral_sweep(polyline, h, r, twist=360, center, anchor, spin=0, orient=UP) { +module spiral_sweep(polyline, h, r, twist=360, center, d, anchor, spin=0, orient=UP) { + r = get_radius(r=r, d=d, dflt=50); polyline = path3d(polyline); pline_count = len(polyline); steps = ceil(segs(r)*(twist/360)); diff --git a/polyhedra.scad b/polyhedra.scad index e845acd..3231606 100644 --- a/polyhedra.scad +++ b/polyhedra.scad @@ -730,9 +730,10 @@ function stellate_faces(scalefactor,stellate,vertices,faces_normals) = ) [newfaces, normals, allpts]; -function trapezohedron(faces, r, side, longside, h) = +function trapezohedron(faces, r, side, longside, h, d) = assert(faces%2==0, "Number of faces must be even") let( + r = get_radius(r=r, d=d, dflt=1), N = faces/2, parmcount = num_defined([r,side,longside,h]) ) diff --git a/shapes.scad b/shapes.scad index d612b0a..b010121 100644 --- a/shapes.scad +++ b/shapes.scad @@ -1498,13 +1498,14 @@ module pie_slice( // Center this part along the concave edge to be chamfered and union it in. // // Usage: -// interior_fillet(l, r, [ang], [overlap]); +// interior_fillet(l, r|d, [ang], [overlap]); // // Arguments: -// l = length of edge to fillet. -// r = radius of fillet. -// ang = angle between faces to fillet. -// overlap = overlap size for unioning with faces. +// l = Length of edge to fillet. +// r = Radius of fillet. +// d = Diameter of fillet. +// ang = Angle between faces to fillet. +// overlap = Overlap size for unioning with faces. // anchor = Translate so anchor point is at origin (0,0,0). See [anchor](attachments.scad#anchor). Default: `FRONT+LEFT` // spin = Rotate this many degrees around the Z axis after anchor. See [spin](attachments.scad#spin). Default: `0` // orient = Vector to rotate top towards, after spin. See [orient](attachments.scad#orient). Default: `UP` @@ -1526,7 +1527,8 @@ module pie_slice( // position(BOT+FRONT) // interior_fillet(l=50, r=10, spin=180, orient=RIGHT); // } -module interior_fillet(l=1.0, r=1.0, ang=90, overlap=0.01, anchor=FRONT+LEFT, spin=0, orient=UP) { +module interior_fillet(l=1.0, r, ang=90, overlap=0.01, d, anchor=FRONT+LEFT, spin=0, orient=UP) { + r = get_radius(r=r, d=d, dflt=1); dy = r/tan(ang/2); steps = ceil(segs(r)*ang/360); step = ang/steps; diff --git a/version.scad b/version.scad index c2e6cda..d10b338 100644 --- a/version.scad +++ b/version.scad @@ -8,7 +8,7 @@ ////////////////////////////////////////////////////////////////////// -BOSL_VERSION = [2,0,409]; +BOSL_VERSION = [2,0,410]; // Section: BOSL Library Version Functions From a8522854e4f4780358de8db845497a5b1bbde6e5 Mon Sep 17 00:00:00 2001 From: Adrian Mariano Date: Thu, 27 Aug 2020 17:07:12 -0400 Subject: [PATCH 07/23] small tweak to hide two functions --- polyhedra.scad | 14 +++++--------- 1 file changed, 5 insertions(+), 9 deletions(-) diff --git a/polyhedra.scad b/polyhedra.scad index 3231606..13cbe54 100644 --- a/polyhedra.scad +++ b/polyhedra.scad @@ -652,7 +652,7 @@ function regular_polyhedron_info( let( entry = ( name == "trapezohedron"? ( - trapezohedron(faces=faces, side=side, longside=longside, h=h, r=r) + _trapezohedron(faces=faces, side=side, longside=longside, h=h, r=r) ) : ( _polyhedra_[!is_undef(index)? indexlist[index] : @@ -671,7 +671,7 @@ function regular_polyhedron_info( ) / entry[edgelen] ), face_triangles = hull(entry[vertices]), - faces_normals_vertices = stellate_faces( + faces_normals_vertices = _stellate_faces( entry[edgelen], stellate, entry[vertices], entry[facevertices]==[3]? [face_triangles, [for(face=face_triangles) _facenormal(entry[vertices],face)]] : @@ -713,11 +713,7 @@ function regular_polyhedron_info( assert(false, str("Unknown info type '",info,"' requested")); - -/// hull solution fails due to roundoff -/// either cross product or just rotate to -/// -function stellate_faces(scalefactor,stellate,vertices,faces_normals) = +function _stellate_faces(scalefactor,stellate,vertices,faces_normals) = (stellate == false || stellate == 0)? concat(faces_normals,[vertices]) : let( faces = [for(face=faces_normals[0]) select(face,hull(select(vertices,face)))], @@ -730,8 +726,8 @@ function stellate_faces(scalefactor,stellate,vertices,faces_normals) = ) [newfaces, normals, allpts]; -function trapezohedron(faces, r, side, longside, h, d) = - assert(faces%2==0, "Number of faces must be even") +function _trapezohedron(faces, r, side, longside, h, d) = + assert(faces%2==0, "Must set 'faces' to an even number for trapezohedron") let( r = get_radius(r=r, d=d, dflt=1), N = faces/2, From c1782f11132431c83034fc4b4ed3859fdb7cf20c Mon Sep 17 00:00:00 2001 From: Adrian Mariano Date: Thu, 27 Aug 2020 19:25:41 -0400 Subject: [PATCH 08/23] added no_children checks, and attachable to vnf_polyhedron --- common.scad | 3 +++ strings.scad | 6 ++++-- structs.scad | 1 + version.scad | 1 + vnf.scad | 11 +++++++++-- 5 files changed, 18 insertions(+), 4 deletions(-) diff --git a/common.scad b/common.scad index db4f3bb..570ff8e 100644 --- a/common.scad +++ b/common.scad @@ -353,6 +353,7 @@ function segs(r) = // Arguments: // $children = number of children the module has. module no_children(count) { + assert($children==0, "Module no_children() does not support child modules"); assert(count==0, str("Module ",parent_module(1),"() does not support child modules")); } @@ -377,6 +378,7 @@ function _valstr(x) = // expected = The value that was expected. // info = Extra info to print out to make the error clearer. module assert_approx(got, expected, info) { + no_children($children); if (!approx(got, expected)) { echo(); echo(str("EXPECT: ", _valstr(expected))); @@ -404,6 +406,7 @@ module assert_approx(got, expected, info) { // expected = The value that was expected. // info = Extra info to print out to make the error clearer. module assert_equal(got, expected, info) { + no_children($children); if (got != expected || (is_nan(got) && is_nan(expected))) { echo(); echo(str("EXPECT: ", _valstr(expected))); diff --git a/strings.scad b/strings.scad index cc0713b..5e34fd9 100644 --- a/strings.scad +++ b/strings.scad @@ -701,7 +701,9 @@ function str_format(fmt, vals, use_nbsp=false) = // echofmt("{:-10s}{:.3f}", ["plecostamus",27.43982]); // ECHO: "plecostamus27.440" // echofmt("{:-10.9s}{:.3f}", ["plecostamus",27.43982]); // ECHO: "plecostam 27.440" function echofmt(fmt, vals, use_nbsp=false) = echo(str_format(fmt,vals,use_nbsp)); -module echofmt(fmt, vals, use_nbsp=false) echo(str_format(fmt,vals,use_nbsp)); - +module echofmt(fmt, vals, use_nbsp=false) { + no_children($children); + echo(str_format(fmt,vals,use_nbsp)); +} // vim: expandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap diff --git a/structs.scad b/structs.scad index e7effb4..1e9f321 100644 --- a/structs.scad +++ b/structs.scad @@ -101,6 +101,7 @@ function struct_echo(struct,name="") = undef; module struct_echo(struct,name="") { + no_children($children); dummy = struct_echo(struct,name); } diff --git a/version.scad b/version.scad index d10b338..b37d43f 100644 --- a/version.scad +++ b/version.scad @@ -49,6 +49,7 @@ function bosl_version_str() = version_to_str(BOSL_VERSION); // Description: // Given a version as a list, number, or string, asserts that the currently installed BOSL library is at least the given version. module bosl_required(target) { + no_children($children); assert( version_cmp(bosl_version(), target) >= 0, str( diff --git a/vnf.scad b/vnf.scad index 05fd82f..edf00e4 100644 --- a/vnf.scad +++ b/vnf.scad @@ -341,9 +341,16 @@ function vnf_vertex_array( // Arguments: // vnf = A VNF structure, or list of VNF structures. // convexity = Max number of times a line could intersect a wall of the shape. -module vnf_polyhedron(vnf, convexity=2) { +module vnf_polyhedron(vnf, + convexity=10,anchor="origin",cp, + spin=0, orient=UP, extent=false +) +{ vnf = is_vnf_list(vnf)? vnf_merge(vnf) : vnf; - polyhedron(vnf[0], vnf[1], convexity=convexity); + attachable(anchor=anchor, spin=spin, orient=orient, vnf=vnf, extent=extent, cp=is_def(cp) ? cp : vnf_centroid(vnf)){ + polyhedron(vnf[0], vnf[1], convexity=convexity); + children(); + } } From 411c1f866eea5f83b51e0569e9d852e1385f27a2 Mon Sep 17 00:00:00 2001 From: Adrian Mariano Date: Thu, 27 Aug 2020 19:48:13 -0400 Subject: [PATCH 09/23] fixed missing ; bug in test_transforms.scad --- tests/test_transforms.scad | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_transforms.scad b/tests/test_transforms.scad index a44e47c..3fb044b 100644 --- a/tests/test_transforms.scad +++ b/tests/test_transforms.scad @@ -232,7 +232,7 @@ module test_rot() { assert_equal(rot(a,p=pts2d), pts2d, info=str("rot(",a,",p=...), 2D")); assert_equal(rot(a,p=pts3d), pts3d, info=str("rot(",a,",p=...), 3D")); } - assert_equal(rot(90), [[0,-1,0,0],[1,0,0,0],[0,0,1,0],[0,0,0,1]]) + assert_equal(rot(90), [[0,-1,0,0],[1,0,0,0],[0,0,1,0],[0,0,0,1]]); for (a=angs) { assert_equal(rot(a), affine3d_zrot(a), info=str("Z angle (only) = ",a)); assert_equal(rot([a,0,0]), affine3d_xrot(a), info=str("X angle = ",a)); From e344de4aa475a67465559c7b0319b161b2824efb Mon Sep 17 00:00:00 2001 From: Garth Minette Date: Thu, 27 Aug 2020 22:29:11 -0700 Subject: [PATCH 10/23] Added cheat-sheet generation script. --- scripts/gencheat.sh | 91 +++++++++++++++++++++++++++++++++++++++++++++ version.scad | 2 +- 2 files changed, 92 insertions(+), 1 deletion(-) create mode 100755 scripts/gencheat.sh diff --git a/scripts/gencheat.sh b/scripts/gencheat.sh new file mode 100755 index 0000000..fb03406 --- /dev/null +++ b/scripts/gencheat.sh @@ -0,0 +1,91 @@ +#!/bin/bash + + +function ucase +{ + echo "$1" | tr '[:lower:]' '[:upper:]' +} + +function lcase +{ + echo "$1" | tr '[:upper:]' '[:lower:]' +} + +function columnize +{ + cols=4 + TMPFILE=$(mktemp -t $(basename $0).XXXXXX) || exit 1 + cat >>$TMPFILE + totcnt=$(wc -l $TMPFILE | awk '{print $1}') + maxrows=$((($totcnt+$cols-1)/$cols)) + maxcols=$cols + if [[ $maxcols -gt $totcnt ]] ; then + maxcols=$totcnt + fi + cnt=0 + hdrln1="| $(ucase $1) " + hdrln2='|:-----' + n=1 + while [[ $n < $maxcols ]] ; do + hdrln1+=' |  ' + hdrln2+=' |:------' + n=$(($n+1)) + done + hdrln1+=' |' + hdrln2+=' |' + n=0 + while [[ $n < $maxrows ]] ; do + lines[$n]="" + n=$(($n+1)) + done + col=0 + while IFS= read -r line; do + if [[ $col != 0 ]] ; then + lines[$cnt]+=" | " + fi + lines[$cnt]+="$line" + cnt=$(($cnt+1)) + if [[ $cnt = $maxrows ]] ; then + cnt=0 + col=$(($col+1)) + fi + done <$TMPFILE + rm -f $TMPFILE + + echo + echo $hdrln1 + echo $hdrln2 + n=0 + while [[ $n < $maxrows ]] ; do + echo "| ${lines[$n]} |" + n=$(($n+1)) + done +} + +function mkconstindex +{ + sed 's/([^)]*)//g' | sed 's/[^a-zA-Z0-9_.:$]//g' | awk -F ':' '{printf "[%s](%s#%s)\n", $3, $1, $3}' +} + +function mkotherindex +{ + sed 's/([^)]*)//g' | sed 's/[^a-zA-Z0-9_.:$]//g' | awk -F ':' '{printf "[%s()](%s#%s)\n", $3, $1, $3}' +} + +CHEAT_FILES=$(grep '^include' std.scad | sed 's/^.*<\([a-zA-Z0-9.]*\)>/\1/'|grep -v 'version.scad') + +( + echo '## Belfry OpenScad Library Cheat Sheet' + echo + echo '( [Alphabetic Index](Index.md) )' + echo + for f in $CHEAT_FILES ; do + #echo "### $f" + ( + egrep -H 'Constant: ' $f | mkconstindex + egrep -H 'Function: |Function&Module: |Module: ' $f | mkotherindex + ) | columnize $f + echo + done +) > BOSL2.wiki/CheatSheet.md + diff --git a/version.scad b/version.scad index d10b338..61bf5a8 100644 --- a/version.scad +++ b/version.scad @@ -8,7 +8,7 @@ ////////////////////////////////////////////////////////////////////// -BOSL_VERSION = [2,0,410]; +BOSL_VERSION = [2,0,412]; // Section: BOSL Library Version Functions From 7051f8ef19b6d077f49b19dd25c9fa463559c5c5 Mon Sep 17 00:00:00 2001 From: Garth Minette Date: Fri, 28 Aug 2020 00:06:28 -0700 Subject: [PATCH 11/23] Fixed cheatsheet index link. --- scripts/gencheat.sh | 2 +- version.scad | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/scripts/gencheat.sh b/scripts/gencheat.sh index fb03406..6762f0a 100755 --- a/scripts/gencheat.sh +++ b/scripts/gencheat.sh @@ -77,7 +77,7 @@ CHEAT_FILES=$(grep '^include' std.scad | sed 's/^.*<\([a-zA-Z0-9.]*\)>/\1/'|grep ( echo '## Belfry OpenScad Library Cheat Sheet' echo - echo '( [Alphabetic Index](Index.md) )' + echo '( [Alphabetic Index](Index) )' echo for f in $CHEAT_FILES ; do #echo "### $f" diff --git a/version.scad b/version.scad index 61bf5a8..55a670a 100644 --- a/version.scad +++ b/version.scad @@ -8,7 +8,7 @@ ////////////////////////////////////////////////////////////////////// -BOSL_VERSION = [2,0,412]; +BOSL_VERSION = [2,0,413]; // Section: BOSL Library Version Functions From e1b0985afcde7e4821ec19f849b425ea859e7eaf Mon Sep 17 00:00:00 2001 From: Garth Minette Date: Fri, 28 Aug 2020 19:07:10 -0700 Subject: [PATCH 12/23] Fixed broken line-plane intersections. Attachment enhanced vnf_polyhedron(). --- attachments.scad | 15 +++++++++------ geometry.scad | 26 +++++++++++++++++--------- tests/test_geometry.scad | 13 ++++++++++++- version.scad | 2 +- vnf.scad | 12 ++++++++++-- 5 files changed, 49 insertions(+), 19 deletions(-) diff --git a/attachments.scad b/attachments.scad index 48ed522..ae93412 100644 --- a/attachments.scad +++ b/attachments.scad @@ -461,13 +461,16 @@ function find_anchor(anchor, geom) = rpts = apply(rot(from=anchor, to=RIGHT) * move(point3d(-cp)), points), hits = [ for (face = faces) let( - verts = select(rpts, face) + verts = select(rpts, face), + xs = subindex(verts,0), + ys = subindex(verts,1), + zs = subindex(verts,2) ) if ( - max(subindex(verts,0)) >= -eps && - max(subindex(verts,1)) >= -eps && - max(subindex(verts,2)) >= -eps && - min(subindex(verts,1)) <= eps && - min(subindex(verts,2)) <= eps + max(xs) >= -eps && + max(ys) >= -eps && + max(zs) >= -eps && + min(ys) <= eps && + min(zs) <= eps ) let( poly = select(points, face), pt = polygon_line_intersection(poly, [cp,cp+anchor], bounded=[true,false], eps=eps) diff --git a/geometry.scad b/geometry.scad index dc86187..0ab9f73 100644 --- a/geometry.scad +++ b/geometry.scad @@ -1007,11 +1007,20 @@ function closest_point_on_plane(plane, point) = // Returns [LINE, undef] if the line is on the plane. // Returns undef if line is parallel to, but not on the given plane. function _general_plane_line_intersection(plane, line, eps=EPSILON) = - let( a = plane*[each line[0],-1], - b = plane*[each(line[1]-line[0]),-1] ) - approx(b,0,eps) - ? points_on_plane(line[0],plane,eps)? [line,undef]: undef - : [ line[0]+a/b*(line[1]-line[0]), a/b ]; + let( + l0 = line[0], // Ray start point + u = line[1] - l0, // Ray direction vector + n = plane_normal(plane), + p0 = n * plane[3], // A point on the plane + w = l0 - p0 // Vector from plane point to ray start + ) approx(n*u, 0, eps=eps) ? ( + // Line is parallel to plane. + approx(n*w, 0, eps=eps) + ? [line, undef] // Line is on the plane. + : undef // Line never intersects the plane. + ) : let( + t = (-n * w) / (n * u) // Distance ratio along ray + ) [ l0 + u*t, t ]; // Function: plane_line_angle() @@ -1098,8 +1107,8 @@ function polygon_line_intersection(poly, line, bounded=false, eps=EPSILON) = linevec = unit(line[1] - line[0]), lp1 = line[0] + (bounded[0]? 0 : -1000000) * linevec, lp2 = line[1] + (bounded[1]? 0 : 1000000) * linevec, - poly2d = clockwise_polygon(project_plane(poly, p1, p2, p3)), - line2d = project_plane([lp1,lp2], p1, p2, p3), + poly2d = clockwise_polygon(project_plane(poly, plane)), + line2d = project_plane([lp1,lp2], plane), parts = split_path_at_region_crossings(line2d, [poly2d], closed=false), inside = [for (part = parts) if (point_in_polygon(mean(part), poly2d)>0) part @@ -1107,7 +1116,7 @@ function polygon_line_intersection(poly, line, bounded=false, eps=EPSILON) = ) !inside? undef : let( - isegs = [for (seg = inside) lift_plane(seg, p1, p2, p3) ] + isegs = [for (seg = inside) lift_plane(seg, plane) ] ) isegs ) @@ -1264,7 +1273,6 @@ function find_circle_2tangents(pt1, pt2, pt3, r, d, tangents=false) = x = hyp * cos(a/2), tp1 = pt2 + x * v1, tp2 = pt2 + x * v2, -// fff=echo(tp1=tp1,cp=cp,pt2=pt2), dang1 = vector_angle(tp1-cp,pt2-cp), dang2 = vector_angle(tp2-cp,pt2-cp) ) diff --git a/tests/test_geometry.scad b/tests/test_geometry.scad index 0950e1c..e69a8ac 100644 --- a/tests/test_geometry.scad +++ b/tests/test_geometry.scad @@ -53,7 +53,7 @@ test_tri_functions(); //test__general_plane_line_intersection(); //test_plane_line_angle(); //test_plane_line_intersection(); -//test_polygon_line_intersection(); +test_polygon_line_intersection(); //test_plane_intersection(); test_coplanar(); test_points_on_plane(); @@ -542,6 +542,17 @@ module test_distance_from_plane() { *test_distance_from_plane(); +module test_polygon_line_intersection() { + poly1 = [[50,50,50], [50,-50,50], [-50,-50,50]]; + assert_approx(polygon_line_intersection(poly1, [CENTER, UP]), [0,0,50]); + assert_approx(polygon_line_intersection(poly1, [CENTER, UP+RIGHT]), [50,0,50]); + assert_approx(polygon_line_intersection(poly1, [CENTER, UP+BACK+RIGHT]), [50,50,50]); + assert_approx(polygon_line_intersection(poly1, [[0,0,50], [1,0,50]]), [[[0,0,50], [50,0,50]]]); + assert_approx(polygon_line_intersection(poly1, [[0,0,0], [1,0,0]]), undef); +} +*test_polygon_line_intersection(); + + module test_coplanar() { assert(coplanar([ [5,5,1],[0,0,1],[-1,-1,1] ]) == false); assert(coplanar([ [5,5,1],[0,0,0],[-1,-1,1] ]) == true); diff --git a/version.scad b/version.scad index 55a670a..6e4605f 100644 --- a/version.scad +++ b/version.scad @@ -8,7 +8,7 @@ ////////////////////////////////////////////////////////////////////// -BOSL_VERSION = [2,0,413]; +BOSL_VERSION = [2,0,414]; // Section: BOSL Library Version Functions diff --git a/vnf.scad b/vnf.scad index 05fd82f..b0a5607 100644 --- a/vnf.scad +++ b/vnf.scad @@ -341,9 +341,17 @@ function vnf_vertex_array( // Arguments: // vnf = A VNF structure, or list of VNF structures. // convexity = Max number of times a line could intersect a wall of the shape. -module vnf_polyhedron(vnf, convexity=2) { +// extent = If true, calculate anchors by extents, rather than intersection. Default: true. +// cp = Centerpoint of VNF to use for anchoring when `extent` is false. Default: `[0, 0, 0]` +// anchor = Translate so anchor point is at origin (0,0,0). See [anchor](attachments.scad#anchor). Default: `"origin"` +// spin = Rotate this many degrees around the Z axis after anchor. See [spin](attachments.scad#spin). Default: `0` +// orient = Vector to rotate top towards, after spin. See [orient](attachments.scad#orient). Default: `UP` +module vnf_polyhedron(vnf, convexity=2, extent=true, cp=[0,0,0], anchor="origin", spin=0, orient=UP) { vnf = is_vnf_list(vnf)? vnf_merge(vnf) : vnf; - polyhedron(vnf[0], vnf[1], convexity=convexity); + attachable(anchor,spin,orient, vnf=vnf, cp=cp, extent=extent) { + polyhedron(vnf[0], vnf[1], convexity=convexity); + children(); + } } From 2eb0ce034854e57654d45dda370440c423f10d99 Mon Sep 17 00:00:00 2001 From: Adrian Mariano Date: Fri, 28 Aug 2020 23:16:11 -0400 Subject: [PATCH 13/23] added os_mask to rounding and tweaked error messages in transpose --- arrays.scad | 26 ++++++++++++++------------ rounding.scad | 24 +++++++++++++++++++++++- 2 files changed, 37 insertions(+), 13 deletions(-) diff --git a/arrays.scad b/arrays.scad index 835fd65..a6d43d1 100644 --- a/arrays.scad +++ b/arrays.scad @@ -1286,10 +1286,12 @@ function array_dim(v, depth=undef) = // Function: transpose() -// Description: Returns the transposition of the given array. -// When reverse=true, the transposition is done in respect to the secondary diagonal, that is: -// . -// reverse(transpose(reverse(arr))) == transpose(arr, reverse=true) +// Usage: +// transpose(arr, [reverse]) +// 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 = [ @@ -1329,19 +1331,19 @@ function array_dim(v, depth=undef) = // // ["h", "e", "b"], // // ["g", "d", "a"] // // ] -// Example: +// 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, "The array is not a vector neither a matrix." ) + assert( is_list(arr) && len(arr)>0, "Input to transpose must be a nonempty list.") is_list(arr[0]) - ? let( l0 = len(arr[0]) ) - assert([for(a=arr) if(!is_list(a) || len(a)!=l0) 1 ]==[], "The array is not a vector neither a matrix." ) + ? 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:l0-1]) - [ for (j=[0:1:len(arr)-1]) arr[len(arr)-1-j][l0-1-i] ] ] - : [for (i=[0:1:l0-1]) + ? [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), "The array is not a vector neither a matrix." ) + : assert( is_vector(arr), "Input to transpose must be a vector or list of lists.") arr; diff --git a/rounding.scad b/rounding.scad index ae2aed0..c57a76f 100644 --- a/rounding.scad +++ b/rounding.scad @@ -488,6 +488,7 @@ function smooth_path(path, tangents, size, relsize, splinesteps=10, uniform=fals // - smooth: os_smooth(cut|joint). Define continuous curvature rounding, with `cut` and `joint` as for round_corners. // - teardrop: os_teardrop(r|cut). Rounding using a 1/8 circle that then changes to a 45 degree chamfer. The chamfer is at the end, and enables the object to be 3d printed without support. The radius gives the radius of the circular part. // - chamfer: os_chamfer([height], [width], [cut], [angle]). Chamfer the edge at desired angle or with desired height and width. You can specify height and width together and the angle will be ignored, or specify just one of height and width and the angle is used to determine the shape. Alternatively, specify "cut" along with angle to specify the cut back distance of the chamfer. +// - mask: os_mask(mask, [out]). Create a profile from one of the [2d masking shapes](shapes2d.scad#5-2d-masking-shapes). The `out` parameter specifies that the mask should flare outward (like crown molding or baseboard). This is set false by default. // . // The general settings that you can use with all of the helper functions are mostly used to control how offset_sweep() calls the offset() function. // - extra: Add an extra vertical step of the specified height, to be used for intersections or differences. This extra step will extend the resulting object beyond the height you specify. Default: 0 @@ -654,6 +655,15 @@ function smooth_path(path, tangents, size, relsize, splinesteps=10, uniform=fals // up(1) // offset_sweep(offset(rhex,r=-1), height=9.5, bottom=os_circle(r=2), top=os_teardrop(r=-4)); // } +// Example: Using os_mask to create ogee profiles: +// ogee = mask2d_ogee([ +// "xstep",1, "ystep",1, // Starting shoulder. +// "fillet",5, "round",5, // S-curve. +// "ystep",1, // Ending shoulder. +// ]); +// star = star(5, r=220, ir=130); +// rounded_star = round_corners(star, cut=flatten(repeat([5,0],5)), $fn=24); +// offset_sweep(rounded_star, height=100, top=os_mask(ogee), bottom=os_mask(ogee,out=true)); // This function does the actual work of repeatedly calling offset() and concatenating the resulting face and vertex lists to produce @@ -880,6 +890,18 @@ function os_profile(points, extra,check_valid, quality, offset_maxstep, offset) ]); +function os_mask(mask, out=false, extra,check_valid, quality, offset_maxstep, offset) = + let( + origin_index = [for(i=idx(mask)) if (mask[i].x<0 && mask[i].y<0) i], + xfactor = out ? -1 : 1 + ) + assert(len(origin_index)==1,"Cannot find origin in the mask") + let( + points = ([for(pt=polygon_shift(mask,origin_index[0])) [xfactor*max(pt.x,0),-max(pt.y,0)]]) + ) + os_profile(deduplicate(move(-points[1],p=select(points,1,-1))), extra,check_valid,quality,offset_maxstep,offset); + + // Module: convex_offset_extrude() // // Description: @@ -1994,4 +2016,4 @@ module bent_cutout_mask(r, thickness, path, convexity=10) } -// vim: expandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap \ No newline at end of file +// vim: expandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap From bbe4ad1467633d0da4bd3fe60af96592af709ba9 Mon Sep 17 00:00:00 2001 From: RonaldoCMP Date: Sun, 30 Aug 2020 12:12:36 +0100 Subject: [PATCH 14/23] Sort debugging and optimizing There were bugs in the previous sorting functions. They didn't check the homogeneity of the input list before calling _sort_scalars and _sort_vectors. The bug might result in wrong order and missing list elements in the output. Besides correcting the bug a recode of all sorting functions result in better performance and a enlargement of their scope. With the new functions, list of vectors of any dimension may be sorted, even with idx given, with the native comparison operators. The scope of indexed sorting is also extended. The file test_arrays has been extended to check the new funcionality. New functions: is_homogeneous - checks if a list has elements of the same type (although not distinguing booleans from numbers) up to a given depth _sort_vectors - internal function to sort homgeneous lists of vectors using native comparison operators; extends the scope of the previous _sort_vectors# functions with better performance _lexical_sort - internal function to sort non-homogeneous lists; uses compare_vals _indexed_sort - internal function to perform indexed sorting of non-homogeneous lists; uses compar_vals Changed/reviewed functions: _valid_idx - doesn't requires the input of imin and imax args sort - explores the internal functions to get better performance and an enlarged scope sortidx - explores the internal functions to get better performance and an enlarged scope _sort_general - just for sortings of non-homogeneous lists using compare_vals _array_dim_recurse - changed for bit better performance Functions eliminated: _sort_vectors1 _sort_vectors2 _sort_vectors3 _sort_vectors4 --- arrays.scad | 311 +++++++++++++++++++++-------------------- tests/test_arrays.scad | 51 +++++-- 2 files changed, 198 insertions(+), 164 deletions(-) diff --git a/arrays.scad b/arrays.scad index a6d43d1..d012c94 100644 --- a/arrays.scad +++ b/arrays.scad @@ -19,6 +19,32 @@ // Section: List Query Operations +// Function: is_homogeneous() +// Usage: +// is_homogeneous(list,depth) +// Description: +// Returns true when the list have elements of same type up to the depth `depth`. +// Booleans and numbers are not distinguinshed as of distinct types. +// Arguments: +// list - the list to check +// depth - the lowest level the check is done +// Example: +// is_homogeneous( [[1,["a"]], [2,["b"]]] ) // Returns true +// is_homogeneous( [[1,["a"]], [2,[true]]] ) // Returns false +// is_homogeneous( [[1,["a"]], [2,[true]]], 1 ) // Returns true +// is_homogeneous( [[1,["a"]], [2,[true]]], 2 ) // Returns false +// is_homogeneous( [[1,["a"]], [true,["b"]]] ) // Returns true +function is_homogeneous(l, depth) = + !is_list(l) || l==[] ? false : + let( l0=l[0] ) + [] == [for(i=[1:len(l)-1]) if( ! _same_type(l[i],l0, depth+1) ) 0 ]; + +function _same_type(a,b, depth) = + (depth==0) || (a>=b) || (a==b) || (a<=b) + || ( is_list(a) && is_list(b) && len(a)==len(b) + && []==[for(i=idx(a)) if( ! _same_type(a[i],b[i],depth-1) )0] ); + + // Function: select() // Description: // Returns a portion of a list, wrapping around past the beginning, if end=imin ) + && ( is_undef(imax) || idx< imax ) ) + || ( is_list(idx) + && ( is_undef(imin) || min(idx)>=imin ) + && ( is_undef(imax) || max(idx)< imax ) ) + || ( is_range(idx) + && ( is_undef(imin) || (idx[1]>0 && idx[0]>=imin ) || (idx[1]<0 && idx[0]<=imax ) ) + && ( is_undef(imax) || (idx[1]>0 && idx[2]<=imax ) || (idx[1]<0 && idx[2]>=imin ) ) ); + + // Function: shuffle() // Description: // Shuffles the input list into random order. @@ -611,7 +654,8 @@ function shuffle(list) = concat(shuffle(left), shuffle(right)); -// Sort a vector of scalar values +// Sort a vector of scalar values with the native comparison operator +// all elements should have the same type. function _sort_scalars(arr) = len(arr)<=1 ? arr : let( @@ -623,105 +667,67 @@ function _sort_scalars(arr) = concat( _sort_scalars(lesser), equal, _sort_scalars(greater) ); -// Sort a vector of vectors based on the first entry only of each vector -function _sort_vectors1(arr) = - len(arr)<=1 ? arr : - !(len(arr)>0) ? [] : +// lexical sort of a homogeneous list of vectors +// uses native comparison operator +function _sort_vectors(arr, _i=0) = + len(arr)<=1 || _i>=len(arr[0]) ? arr : let( - pivot = arr[floor(len(arr)/2)], - lesser = [ for (y = arr) if (y[0] < pivot[0]) y ], - equal = [ for (y = arr) if (y[0] == pivot[0]) y ], - greater = [ for (y = arr) if (y[0] > pivot[0]) y ] - ) - concat( _sort_vectors1(lesser), equal, _sort_vectors1(greater) ); + pivot = arr[floor(len(arr)/2)][_i], + lesser = [ for (entry=arr) if (entry[_i] < pivot ) entry ], + equal = [ for (entry=arr) if (entry[_i] == pivot ) entry ], + greater = [ for (entry=arr) if (entry[_i] > pivot ) entry ] + ) + concat( + _sort_vectors(lesser, _i ), + _sort_vectors(equal, _i+1 ), + _sort_vectors(greater, _i ) ); + - -// Sort a vector of vectors based on the first two entries of each vector -// Lexicographic order, remaining entries of vector ignored -function _sort_vectors2(arr) = - len(arr)<=1 ? arr : - !(len(arr)>0) ? [] : +// lexical sort of a homogeneous list of vectors by the vector components with indices in idxlist +// all idxlist indices should be in the range of the vector dimensions +// idxlist must be undef or a simple list of numbers +// uses native comparison operator +function _sort_vectors(arr, idxlist, _i=0) = + len(arr)<=1 || ( is_list(idxlist) && _i>=len(idxlist) ) || _i>=len(arr[0]) ? arr : let( - pivot = arr[floor(len(arr)/2)], - lesser = [ for (y = arr) if (y[0] < pivot[0] || (y[0]==pivot[0] && y[1] pivot[0] || (y[0]==pivot[0] && y[1]>pivot[1])) y ] - ) - concat( _sort_vectors2(lesser), equal, _sort_vectors2(greater) ); + k = is_list(idxlist) ? idxlist[_i] : _i, + pivot = arr[floor(len(arr)/2)][k], + lesser = [ for (entry=arr) if (entry[k] < pivot ) entry ], + equal = [ for (entry=arr) if (entry[k] == pivot ) entry ], + greater = [ for (entry=arr) if (entry[k] > pivot ) entry ] + ) + concat( + _sort_vectors(lesser, idxlist, _i ), + _sort_vectors(equal, idxlist, _i+1), + _sort_vectors(greater, idxlist, _i ) ); + - -// Sort a vector of vectors based on the first three entries of each vector -// Lexicographic order, remaining entries of vector ignored -function _sort_vectors3(arr) = - len(arr)<=1 ? arr : let( - pivot = arr[floor(len(arr)/2)], - lesser = [ for (y = arr) - if ( y[0] < pivot[0] - || ( y[0]==pivot[0] - && ( y[1] pivot[0] - || ( y[0]==pivot[0] - && ( y[1] > pivot[1] - || ( y[1]==pivot[1] - && y[2] > pivot[2] )))) - y ] - ) concat( _sort_vectors3(lesser), equal, _sort_vectors3(greater) ); - - -// Sort a vector of vectors based on the first four entries of each vector -// Lexicographic order, remaining entries of vector ignored -function _sort_vectors4(arr) = - len(arr)<=1 ? arr : let( - pivot = arr[floor(len(arr)/2)], - lesser = [ for (y = arr) - if ( y[0] < pivot[0] - || ( y[0]==pivot[0] - && ( y[1] pivot[0] - || ( y[0]==pivot[0] - && ( y[1]>pivot[1] - || ( y[1]==pivot[1] - && ( y[2]>pivot[2] - || ( y[2]==pivot[2] - && y[3]>pivot[3] )))))) - y ] - ) concat( _sort_vectors4(lesser), equal, _sort_vectors4(greater) ); - - -// when idx==undef, returns the sorted array -// otherwise, returns the indices of the sorted array -function _sort_general(arr, idx=undef) = +// sorting using compare_vals(); returns indexed list when `indexed==true` +function _sort_general(arr, idx=undef, indexed=false) = (len(arr)<=1) ? arr : - is_undef(idx) - ? _sort_scalar(arr) - : let( arrind=[for(k=[0:len(arr)-1], ark=[arr[k]]) [ k, [for (i=idx) ark[i]] ] ] ) - _indexed_sort(arrind); + ! indexed && is_undef(idx) + ? _lexical_sort(arr) + : let( arrind = _indexed_sort(enumerate(arr,idx)) ) + indexed + ? arrind + : [for(i=arrind) arr[i]]; +// lexical sort using compare_vals() +function _lexical_sort(arr) = + arr==[] ? [] : len(arr)==1? arr : + let( pivot = arr[floor(len(arr)/2)] ) + let( + lesser = [ for (entry=arr) if (compare_vals(entry, pivot) <0 ) entry ], + equal = [ for (entry=arr) if (compare_vals(entry, pivot)==0 ) entry ], + greater = [ for (entry=arr) if (compare_vals(entry, pivot) >0 ) entry ] + ) + concat(_lexical_sort(lesser), equal, _lexical_sort(greater)); + + // given a list of pairs, return the first element of each pair of the list sorted by the second element of the pair // the sorting is done using compare_vals() function _indexed_sort(arrind) = - arrind==[] ? [] : len(arrind)==1? [arrind[0][0]] : + arrind==[] ? [] : len(arrind)==1? [arrind[0][0]] : let( pivot = arrind[floor(len(arrind)/2)][1] ) let( lesser = [ for (entry=arrind) if (compare_vals(entry[1], pivot) <0 ) entry ], @@ -731,52 +737,46 @@ function _indexed_sort(arrind) = concat(_indexed_sort(lesser), equal, _indexed_sort(greater)); -// returns true for valid index specifications idx in the interval [imin, imax) -// note that idx can't have any value greater or EQUAL to imax -function _valid_idx(idx,imin,imax) = - is_undef(idx) - || ( is_finite(idx) && idx>=imin && idx< imax ) - || ( is_list(idx) && min(idx)>=imin && max(idx)< imax ) - || ( valid_range(idx) && idx[0]>=imin && idx[2]< imax ); - - // Function: sort() // Usage: // sort(list, [idx]) // Description: -// Sorts the given list using `compare_vals()`, sorting in lexicographic order, with types ordered according to +// Sorts the given list in lexicographic order. If the input is a homogeneous simple list or a homogeneous +// list of vectors (see function is_homogeneous), the sorting method uses the native comparison operator and is faster. +// When sorting non homogeneous list the elements are compared with `compare_vals`, with types ordered according to // `undef < boolean < number < string < list`. Comparison of lists is recursive. -// If the list is a list of vectors whose length is from 1 to 4 and the `idx` parameter is not passed, then -// `sort` uses a much more efficient method for comparisons and will run much faster. In this case, all entries -// in the data are compared using the native comparison operator, so comparisons between types will fail. +// When comparing vectors, homogeneous or not, the parameter `idx` may be used to select the components to compare. +// Note that homogeneous lists of vectors may contain mixed types provided that for any two list elements +// list[i] and list[j] satisfies type(list[i][k])==type(list[j][k]) for all k. +// Strings are allowed as any list element and are compared with the native operators although no substring +// comparison is possible. // Arguments: // list = The list to sort. // idx = If given, do the comparison based just on the specified index, range or list of indices. -// Example: -// l = [45,2,16,37,8,3,9,23,89,12,34]; -// sorted = sort(l); // Returns [2,3,8,9,12,16,23,34,37,45,89] -function sort(list, idx=undef) = +// Example: +// // Homogeneous lists +// l1 = [45,2,16,37,8,3,9,23,89,12,34]; +// sorted1 = sort(l1); // Returns [2,3,8,9,12,16,23,34,37,45,89] +// l2 = [["oat",0], ["cat",1], ["bat",3], ["bat",2], ["fat",3]]; +// sorted2 = sort(l2); // Returns: [["bat",2],["bat",3],["cat",1],["fat",3],["oat",0]] +// // Non-homegenous list +// l3 = [[4,0],[7],[3,9],20,[4],[3,1],[8]]; +// sorted3 = sort(l3); // Returns: [20,[3,1],[3,9],[4],[4,0],[7],[8]] +function sort(list, idx=undef) = !is_list(list) || len(list)<=1 ? list : - is_def(idx) - ? assert( _valid_idx(idx,0,len(list)) , "Invalid indices.") - let( sarr = _sort_general(list,idx) ) - [for(i=[0:len(sarr)-1]) list[sarr[i]] ] - : let(size = array_dim(list)) - len(size)==1 ? _sort_scalars(list) : - len(size)==2 && size[1] <=4 - ? ( - size[1]==0 ? list : - size[1]==1 ? _sort_vectors1(list) : - size[1]==2 ? _sort_vectors2(list) : - size[1]==3 ? _sort_vectors3(list) - /*size[1]==4*/ : _sort_vectors4(list) - ) - : _sort_general(list); - + is_homogeneous(list,1) + ? let(size = array_dim(list[0])) + size==0 ? _sort_scalars(list) + : len(size)!=1 ? _sort_general(list,idx) + : is_undef(idx) ? _sort_vectors(list) + : assert( _valid_idx(idx) , "Invalid indices.") + _sort_vectors(list,[for(i=idx) i]) + : _sort_general(list,idx); + // Function: sortidx() // Description: -// Given a list, calculates the sort order of the list, and returns +// Given a list, sort it as function `sort()`, and returns // a list of indexes into the original list in that sorted order. // If you iterate the returned list in order, and use the list items // to index into the original list, you will be iterating the original @@ -795,31 +795,28 @@ function sort(list, idx=undef) = // idxs1 = sortidx(lst, idx=1); // Returns: [3,0,2,1] // idxs2 = sortidx(lst, idx=0); // Returns: [1,2,0,3] // idxs3 = sortidx(lst, idx=[1,3]); // Returns: [3,0,2,1] -function sortidx(list, idx=undef) = - assert( is_list(list) || is_string(list) , "Invalid input to sort." ) - assert( _valid_idx(idx,0,len(list)) , "Invalid indices.") - list==[] ? [] : - let( - size = array_dim(list), - aug = is_undef(idx) && (len(size) == 1 || (len(size) == 2 && size[1]<=4)) - ? zip(list, list_range(len(list))) - : 0 - ) - is_undef(idx) && len(size) == 1? subindex(_sort_vectors1(aug),1) : - is_undef(idx) && len(size) == 2 && size[1] <=4 - ? ( - size[1]==0 ? list_range(len(arr)) : - size[1]==1 ? subindex(_sort_vectors1(aug),1) : - size[1]==2 ? subindex(_sort_vectors2(aug),2) : - size[1]==3 ? subindex(_sort_vectors3(aug),3) - /*size[1]==4*/ : subindex(_sort_vectors4(aug),4) - ) - : // general case - _sort_general(list,idx); - - -// sort() does not accept strings but sortidx does; isn't inconsistent ? - +function sortidx(list, idx=undef) = + !is_list(list) || len(list)<=1 ? list : + is_homogeneous(list,1) + ? let( + size = array_dim(list[0]), + aug = ! (size==0 || len(size)==1) ? 0 // for general sorting + : [for(i=[0:len(list)-1]) concat(i,list[i])], // for scalar or vector sorting + lidx = size==0? [1] : // scalar sorting + len(size)==1 + ? is_undef(idx) ? [for(i=[0:len(list[0])-1]) i+1] // vector sorting + : [for(i=idx) i+1] // vector sorting + : 0 // just to signal + ) + assert( ! ( size==0 && is_def(idx) ), + "The specification of `idx` is incompatible with scalar sorting." ) + assert( _valid_idx(idx) , "Invalid indices." ) + lidx!=0 + ? let( lsort = _sort_vectors(aug,lidx) ) + [for(li=lsort) li[0] ] + : _sort_general(list,idx,indexed=true) + : _sort_general(list,idx,indexed=true); + // Function: unique() // Usage: @@ -1241,16 +1238,22 @@ function full_flatten(l) = [for(a=l) if(is_list(a)) (each full_flatten(a)) else // Internal. Not exposed. function _array_dim_recurse(v) = !is_list(v[0]) - ? sum( [for(entry=v) is_list(entry) ? 1 : 0] ) == 0 ? [] : [undef] + ? len( [for(entry=v) if(!is_list(entry)) 0] ) == 0 ? [] : [undef] : let( - firstlen = len(v[0]), - first = sum( [for(entry = v) len(entry) == firstlen ? 0 : 1] ) == 0 ? firstlen : undef, + firstlen = is_list(v[0]) ? len(v[0]): undef, + first = len( [for(entry = v) if(! is_list(entry) || (len(entry) != firstlen)) 0 ] ) == 0 ? firstlen : undef, leveldown = flatten(v) ) is_list(leveldown[0]) ? concat([first],_array_dim_recurse(leveldown)) : [first]; +function _array_dim_recurse(v) = + let( alen = [for(vi=v) is_list(vi) ? len(vi): -1] ) + v==[] || max(alen)==-1 ? [] : + let( add = max(alen)!=min(alen) ? undef : alen[0] ) + concat( add, _array_dim_recurse(flatten(v))); + // Function: array_dim() // Usage: @@ -1281,6 +1284,8 @@ function array_dim(v, depth=undef) = ? len(v) : let( dimlist = _array_dim_recurse(v)) (depth > len(dimlist))? 0 : dimlist[depth-1] ; + + // This function may return undef! diff --git a/tests/test_arrays.scad b/tests/test_arrays.scad index 8df23a8..45cd2cd 100644 --- a/tests/test_arrays.scad +++ b/tests/test_arrays.scad @@ -3,6 +3,16 @@ include <../std.scad> // Section: List Query Operations +module test_is_homogeneous(){ + assert(is_homogeneous([[1,["a"]], [2,["b"]]])==true); + assert(is_homogeneous([[1,["a"]], [2,[true]]])==false); + assert(is_homogeneous([[1,["a"]], [2,[true]]],1)==true); + assert(is_homogeneous([[1,["a"]], [2,[true]]],2)==false); + assert(is_homogeneous([[1,["a"]], [true,["b"]]])==true); +} +test_is_homogeneous(); + + module test_select() { l = [3,4,5,6,7,8,9]; assert(select(l, 5, 6) == [8,9]); @@ -46,7 +56,6 @@ test_in_list(); module test_min_index() { assert(min_index([5,3,9,6,2,7,8,2,1])==8); assert(min_index([5,3,9,6,2,7,8,2,7],all=true)==[4,7]); -// assert(min_index([],all=true)==[]); } test_min_index(); @@ -54,7 +63,6 @@ test_min_index(); module test_max_index() { assert(max_index([5,3,9,6,2,7,8,9,1])==2); assert(max_index([5,3,9,6,2,7,8,9,7],all=true)==[2,7]); -// assert(max_index([],all=true)==[]); } test_max_index(); @@ -74,6 +82,7 @@ module test_list_decreasing() { } test_list_decreasing(); + // Section: Basic List Generation module test_repeat() { @@ -155,7 +164,6 @@ module test_list_remove() { } test_list_remove(); - module test_list_remove_values() { animals = ["bat", "cat", "rat", "dog", "bat", "rat"]; assert(list_remove_values(animals, "rat") == ["bat","cat","dog","bat","rat"]); @@ -258,15 +266,24 @@ test_shuffle(); module test_sort() { assert(sort([7,3,9,4,3,1,8]) == [1,3,3,4,7,8,9]); - assert(sort(["cat", "oat", "sat", "bat", "vat", "rat", "pat", "mat", "fat", "hat", "eat"]) == ["bat", "cat", "eat", "fat", "hat", "mat", "oat", "pat", "rat", "sat", "vat"]); + assert(sort([[4,0],[7],[3,9],20,[4],[3,1],[8]]) == [20,[3,1],[3,9],[4],[4,0],[7],[8]]); + assert(sort([[4,0],[7],[3,9],20,[4],[3,1],[8]],idx=1) == [[7],20,[4],[8],[4,0],[3,1],[3,9]]); + assert(sort([[8,6],[3,1],[9,2],[4,3],[3,4],[1,5],[8,0]]) == [[1,5],[3,1],[3,4],[4,3],[8,0],[8,6],[9,2]]); + assert(sort([[8,0],[3,1],[9,2],[4,3],[3,4],[1,5],[8,6]],idx=1) == [[8,0],[3,1],[9,2],[4,3],[3,4],[1,5],[8,6]]); + assert(sort(["cat", "oat", "sat", "bat", "vat", "rat", "pat", "mat", "fat", "hat", "eat"]) + == ["bat", "cat", "eat", "fat", "hat", "mat", "oat", "pat", "rat", "sat", "vat"]); assert(sort(enumerate([[2,3,4],[1,2,3],[2,4,3]]),idx=1)==[[1,[1,2,3]], [0,[2,3,4]], [2,[2,4,3]]]); + assert(sort([0,"1",[1,0],2,"a",[1]])== [0,2,"1","a",[1],[1,0]]); + assert(sort([["oat",0], ["cat",1], ["bat",3], ["bat",2], ["fat",3]])== [["bat",2],["bat",3],["cat",1],["fat",3],["oat",0]]); } test_sort(); module test_sortidx() { - lst1 = ["d","b","e","c"]; + lst1 = ["da","bax","eaw","cav"]; assert(sortidx(lst1) == [1,3,0,2]); + lst5 = [3,5,1,7]; + assert(sortidx(lst5) == [2,0,1,3]); lst2 = [ ["foo", 88, [0,0,1], false], ["bar", 90, [0,1,0], true], @@ -276,11 +293,22 @@ module test_sortidx() { assert(sortidx(lst2, idx=1) == [3,0,2,1]); assert(sortidx(lst2, idx=0) == [1,2,0,3]); assert(sortidx(lst2, idx=[1,3]) == [3,0,2,1]); - lst3 = [[-4, 0, 0], [0, 0, -4], [0, -4, 0], [-4, 0, 0], [0, -4, 0], [0, 0, 4], [0, 0, -4], [0, 4, 0], [4, 0, 0], [0, 0, 4], [0, 4, 0], [4, 0, 0]]; + lst3 = [[-4,0,0],[0,0,-4],[0,-4,0],[-4,0,0],[0,-4,0],[0,0,4], + [0,0,-4],[0,4,0],[4,0,0],[0,0,4],[0,4,0],[4,0,0]]; assert(sortidx(lst3)==[0,3,2,4,1,6,5,9,7,10,8,11]); + assert(sortidx([[4,0],[7],[3,9],20,[4],[3,1],[8]]) == [3,5,2,4,0,1,6]); + assert(sortidx([[4,0],[7],[3,9],20,[4],[3,1],[8]],idx=1) == [1,3,4,6,0,5,2]); + lst4=[0,"1",[1,0],2,"a",[1]]; + assert(sortidx(lst4)== [0,3,1,4,5,2]); + assert(sortidx(["cat","oat","sat","bat","vat","rat","pat","mat","fat","hat","eat"]) + == [3,0,10,8,9,7,1,6,5,2,4]); + assert(sortidx([["oat",0], ["cat",1], ["bat",3], ["bat",2], ["fat",3]])== [3,2,1,4,0]); + assert(sortidx(["Belfry", "OpenScad", "Library", "Documentation"])==[0,3,2,1]); + assert(sortidx(["x",1,[],0,"abc",true])==[5,3,1,4,0,2]); } test_sortidx(); + module test_unique() { assert(unique([]) == []); assert(unique([8]) == [8]); @@ -336,10 +364,8 @@ module test_set_intersection() { test_set_intersection(); - // Arrays - module test_add_scalar() { assert(add_scalar([1,2,3],3) == [4,5,6]); assert(add_scalar([[1,2,3],[3,4,5]],3) == [[4,5,6],[6,7,8]]); @@ -473,6 +499,12 @@ module test_array_dim() { assert(array_dim([[[1,2,3],[4,5,6]],[[7,8,9],[10,11,12]]], 0) == 2); assert(array_dim([[[1,2,3],[4,5,6]],[[7,8,9],[10,11,12]]], 2) == 3); assert(array_dim([[[1,2,3],[4,5,6]],[[7,8,9]]]) == [2,undef,3]); + assert(array_dim([1,2,3,4,5,6,7,8,9]) == [9]); + assert(array_dim([[1],[2],[3],[4],[5],[6],[7],[8],[9]]) == [9,1]); + assert(array_dim([]) == [0]); + assert(array_dim([[]]) == [1,0]); + assert(array_dim([[],[]]) == [2,0]); + assert(array_dim([[],[1]]) == [2,undef]); } test_array_dim(); @@ -486,7 +518,4 @@ module test_transpose() { test_transpose(); -cube(); - - // vim: expandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap From 42bffef35ff39e8fa0c2b0caeec77b53ad517018 Mon Sep 17 00:00:00 2001 From: RonaldoCMP Date: Sun, 30 Aug 2020 12:18:53 +0100 Subject: [PATCH 15/23] Minor doc corrections --- arrays.scad | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/arrays.scad b/arrays.scad index d012c94..e2d3ac2 100644 --- a/arrays.scad +++ b/arrays.scad @@ -26,8 +26,8 @@ // Returns true when the list have elements of same type up to the depth `depth`. // Booleans and numbers are not distinguinshed as of distinct types. // Arguments: -// list - the list to check -// depth - the lowest level the check is done +// list = the list to check +// depth = the lowest level the check is done // Example: // is_homogeneous( [[1,["a"]], [2,["b"]]] ) // Returns true // is_homogeneous( [[1,["a"]], [2,[true]]] ) // Returns false From 7c42c7f5e37114343dc34be3c4d93861f88ba581 Mon Sep 17 00:00:00 2001 From: Garth Minette Date: Mon, 31 Aug 2020 17:31:55 -0700 Subject: [PATCH 16/23] Added libfile names to index. --- scripts/genindex.sh | 36 ++++++++++++++++++++++++++++-------- version.scad | 2 +- 2 files changed, 29 insertions(+), 9 deletions(-) diff --git a/scripts/genindex.sh b/scripts/genindex.sh index 9d3e8cd..51a9075 100755 --- a/scripts/genindex.sh +++ b/scripts/genindex.sh @@ -6,16 +6,18 @@ function ucase echo "$1" | tr '[:lower:]' '[:upper:]' } + function lcase { echo "$1" | tr '[:upper:]' '[:lower:]' } -function mkindex + +function alphaindex { - TMPFILE=$(mktemp -t $(basename $0).XXXXXX) || exit 1 - sed 's/([^)]*)//g' | sed 's/[^a-zA-Z0-9_.:$]//g' | awk -F ':' '{printf "- [%s](%s#%s)\n", $3, $1, $3}' | sort -d -f -u >> $TMPFILE alpha="A B C D E F G H I J K L M N O P Q R S T U V W X Y Z" + TMPFILE=$(mktemp -t $(basename $0).XXXXXX) || exit 1 + sort -d -f >> $TMPFILE for a in $alpha; do echo -n "[$a](#$(lcase "$a")) " done @@ -33,13 +35,31 @@ function mkindex } +function constlist +{ + sed 's/([^)]*)//g' | sed 's/[^a-zA-Z0-9_.:$]//g' | awk -F ':' '{printf "- [%s](%s#%s) (in %s)\n", $3, $1, $3, $1}' +} + +function constlist2 +{ + sed 's/ *=.*$//' | sed 's/[^a-zA-Z0-9_.:$]//g' | awk -F ':' '{printf "- [%s](%s#%s) (in %s)\n", $2, $1, $2, $1}' +} + + +function funclist +{ + sed 's/([^)]*)//g' | sed 's/[^a-zA-Z0-9_.:$]//g' | awk -F ':' '{printf "- [%s()](%s#%s) (in %s)\n", $3, $1, $3, $1}' +} + + ( echo "## Belfry OpenScad Library Index" ( - grep 'Constant: ' *.scad - grep 'Function: ' *.scad - grep 'Function&Module: ' *.scad - grep 'Module: ' *.scad - ) | mkindex + ( + grep 'Constant: ' *.scad | constlist + grep '^[A-Z]* *=.*//' *.scad | constlist2 + ) | sort -u + egrep 'Function: |Function&Module: |Module: ' *.scad | sort -u | funclist + ) | sort | alphaindex ) > BOSL2.wiki/Index.md diff --git a/version.scad b/version.scad index bdbdafa..1d40e98 100644 --- a/version.scad +++ b/version.scad @@ -8,7 +8,7 @@ ////////////////////////////////////////////////////////////////////// -BOSL_VERSION = [2,0,415]; +BOSL_VERSION = [2,0,416]; // Section: BOSL Library Version Functions From 38a4c12032b43f37936345ffe0a7b3657302f6d5 Mon Sep 17 00:00:00 2001 From: Garth Minette Date: Mon, 31 Aug 2020 18:03:48 -0700 Subject: [PATCH 17/23] Cheat Sheet generation script bugfixes. --- scripts/gencheat.sh | 17 ++++++++++++----- version.scad | 2 +- 2 files changed, 13 insertions(+), 6 deletions(-) diff --git a/scripts/gencheat.sh b/scripts/gencheat.sh index 6762f0a..da801c2 100755 --- a/scripts/gencheat.sh +++ b/scripts/gencheat.sh @@ -26,7 +26,7 @@ function columnize hdrln1="| $(ucase $1) " hdrln2='|:-----' n=1 - while [[ $n < $maxcols ]] ; do + while [[ $n -lt $maxcols ]] ; do hdrln1+=' |  ' hdrln2+=' |:------' n=$(($n+1)) @@ -34,7 +34,7 @@ function columnize hdrln1+=' |' hdrln2+=' |' n=0 - while [[ $n < $maxrows ]] ; do + while [[ $n -lt $maxrows ]] ; do lines[$n]="" n=$(($n+1)) done @@ -56,7 +56,7 @@ function columnize echo $hdrln1 echo $hdrln2 n=0 - while [[ $n < $maxrows ]] ; do + while [[ $n -lt $maxrows ]] ; do echo "| ${lines[$n]} |" n=$(($n+1)) done @@ -67,6 +67,11 @@ function mkconstindex sed 's/([^)]*)//g' | sed 's/[^a-zA-Z0-9_.:$]//g' | awk -F ':' '{printf "[%s](%s#%s)\n", $3, $1, $3}' } +function mkconstindex2 +{ + sed 's/ *=.*$//' | sed 's/[^a-zA-Z0-9_.:$]//g' | awk -F ':' '{printf "[%s](%s#%s)\n", $2, $1, $2}' +} + function mkotherindex { sed 's/([^)]*)//g' | sed 's/[^a-zA-Z0-9_.:$]//g' | awk -F ':' '{printf "[%s()](%s#%s)\n", $3, $1, $3}' @@ -80,9 +85,11 @@ CHEAT_FILES=$(grep '^include' std.scad | sed 's/^.*<\([a-zA-Z0-9.]*\)>/\1/'|grep echo '( [Alphabetic Index](Index) )' echo for f in $CHEAT_FILES ; do - #echo "### $f" ( - egrep -H 'Constant: ' $f | mkconstindex + ( + grep -H 'Constant: ' $f | mkconstindex + grep -H '^[A-Z$][A-Z0-9_]* *=.*//' $f | mkconstindex2 + ) | sort -u egrep -H 'Function: |Function&Module: |Module: ' $f | mkotherindex ) | columnize $f echo diff --git a/version.scad b/version.scad index 1d40e98..ffa8d27 100644 --- a/version.scad +++ b/version.scad @@ -8,7 +8,7 @@ ////////////////////////////////////////////////////////////////////// -BOSL_VERSION = [2,0,416]; +BOSL_VERSION = [2,0,417]; // Section: BOSL Library Version Functions From cb36b5d94f539d41b31e206fba2b88f76a8527f3 Mon Sep 17 00:00:00 2001 From: Garth Minette Date: Tue, 1 Sep 2020 00:54:27 -0700 Subject: [PATCH 18/23] Bugfixes for cheat sheet generation. --- .github/workflows/docsgen.yml | 5 ++ scripts/gencheat.sh | 98 +++++++++++++++++------------------ version.scad | 2 +- 3 files changed, 55 insertions(+), 50 deletions(-) diff --git a/.github/workflows/docsgen.yml b/.github/workflows/docsgen.yml index 8ab135c..911fac2 100644 --- a/.github/workflows/docsgen.yml +++ b/.github/workflows/docsgen.yml @@ -41,6 +41,11 @@ jobs: cd $GITHUB_WORKSPACE ./scripts/genindex.sh + - name: Generate Cheat Sheet + run: | + cd $GITHUB_WORKSPACE + ./scripts/gencheat.sh + - name: Generating Docs env: GH_PAT: ${{ secrets.GH_PAT }} diff --git a/scripts/gencheat.sh b/scripts/gencheat.sh index da801c2..6fa49ec 100755 --- a/scripts/gencheat.sh +++ b/scripts/gencheat.sh @@ -16,50 +16,52 @@ function columnize cols=4 TMPFILE=$(mktemp -t $(basename $0).XXXXXX) || exit 1 cat >>$TMPFILE - totcnt=$(wc -l $TMPFILE | awk '{print $1}') - maxrows=$((($totcnt+$cols-1)/$cols)) - maxcols=$cols - if [[ $maxcols -gt $totcnt ]] ; then - maxcols=$totcnt - fi - cnt=0 - hdrln1="| $(ucase $1) " - hdrln2='|:-----' - n=1 - while [[ $n -lt $maxcols ]] ; do - hdrln1+=' |  ' - hdrln2+=' |:------' - n=$(($n+1)) - done - hdrln1+=' |' - hdrln2+=' |' - n=0 - while [[ $n -lt $maxrows ]] ; do - lines[$n]="" - n=$(($n+1)) - done - col=0 - while IFS= read -r line; do - if [[ $col != 0 ]] ; then - lines[$cnt]+=" | " + if [[ $(wc -l $TMPFILE | awk '{print $1}') -gt 1 ]] ; then + totcnt=$(wc -l $TMPFILE | awk '{print $1}') + maxrows=$((($totcnt+$cols-1)/$cols)) + maxcols=$cols + if [[ $maxcols -gt $totcnt ]] ; then + maxcols=$totcnt fi - lines[$cnt]+="$line" - cnt=$(($cnt+1)) - if [[ $cnt = $maxrows ]] ; then - cnt=0 - col=$(($col+1)) - fi - done <$TMPFILE - rm -f $TMPFILE + cnt=0 + hdrln1="| $(ucase $1) " + hdrln2='|:-----' + n=1 + while [[ $n -lt $maxcols ]] ; do + hdrln1+=' |  ' + hdrln2+=' |:------' + n=$(($n+1)) + done + hdrln1+=' |' + hdrln2+=' |' + n=0 + while [[ $n -lt $maxrows ]] ; do + lines[$n]="" + n=$(($n+1)) + done + col=0 + while IFS= read -r line; do + if [[ $col != 0 ]] ; then + lines[$cnt]+=" | " + fi + lines[$cnt]+="$line" + cnt=$(($cnt+1)) + if [[ $cnt = $maxrows ]] ; then + cnt=0 + col=$(($col+1)) + fi + done <$TMPFILE + rm -f $TMPFILE - echo - echo $hdrln1 - echo $hdrln2 - n=0 - while [[ $n -lt $maxrows ]] ; do - echo "| ${lines[$n]} |" - n=$(($n+1)) - done + echo + echo $hdrln1 + echo $hdrln2 + n=0 + while [[ $n -lt $maxrows ]] ; do + echo "| ${lines[$n]} |" + n=$(($n+1)) + done + fi } function mkconstindex @@ -84,14 +86,12 @@ CHEAT_FILES=$(grep '^include' std.scad | sed 's/^.*<\([a-zA-Z0-9.]*\)>/\1/'|grep echo echo '( [Alphabetic Index](Index) )' echo + ( + grep -H '// Constant: ' $CHEAT_FILES | mkconstindex + grep -H '^[A-Z$][A-Z0-9_]* *=.*//' $CHEAT_FILES | mkconstindex2 + ) | sort -u | columnize 'Constants' for f in $CHEAT_FILES ; do - ( - ( - grep -H 'Constant: ' $f | mkconstindex - grep -H '^[A-Z$][A-Z0-9_]* *=.*//' $f | mkconstindex2 - ) | sort -u - egrep -H 'Function: |Function&Module: |Module: ' $f | mkotherindex - ) | columnize $f + egrep -H '// Function: |// Function&Module: |// Module: ' $f | mkotherindex | sort -u | columnize $f echo done ) > BOSL2.wiki/CheatSheet.md diff --git a/version.scad b/version.scad index ffa8d27..6ef424a 100644 --- a/version.scad +++ b/version.scad @@ -8,7 +8,7 @@ ////////////////////////////////////////////////////////////////////// -BOSL_VERSION = [2,0,417]; +BOSL_VERSION = [2,0,418]; // Section: BOSL Library Version Functions From 8331d8e8032b209d5b60b8c3a40322faeaab5de9 Mon Sep 17 00:00:00 2001 From: Garth Minette Date: Tue, 1 Sep 2020 01:23:01 -0700 Subject: [PATCH 19/23] Cheat Sheet ordering tweak. --- scripts/gencheat.sh | 8 ++++---- std.scad | 28 +++++++++++++--------------- version.scad | 2 +- 3 files changed, 18 insertions(+), 20 deletions(-) diff --git a/scripts/gencheat.sh b/scripts/gencheat.sh index 6fa49ec..e07a909 100755 --- a/scripts/gencheat.sh +++ b/scripts/gencheat.sh @@ -13,7 +13,7 @@ function lcase function columnize { - cols=4 + cols=$2 TMPFILE=$(mktemp -t $(basename $0).XXXXXX) || exit 1 cat >>$TMPFILE if [[ $(wc -l $TMPFILE | awk '{print $1}') -gt 1 ]] ; then @@ -79,7 +79,7 @@ function mkotherindex sed 's/([^)]*)//g' | sed 's/[^a-zA-Z0-9_.:$]//g' | awk -F ':' '{printf "[%s()](%s#%s)\n", $3, $1, $3}' } -CHEAT_FILES=$(grep '^include' std.scad | sed 's/^.*<\([a-zA-Z0-9.]*\)>/\1/'|grep -v 'version.scad') +CHEAT_FILES=$(grep '^include' std.scad | sed 's/^.*<\([a-zA-Z0-9.]*\)>/\1/' | grep -v 'version.scad' | grep -v 'primitives.scad') ( echo '## Belfry OpenScad Library Cheat Sheet' @@ -89,9 +89,9 @@ CHEAT_FILES=$(grep '^include' std.scad | sed 's/^.*<\([a-zA-Z0-9.]*\)>/\1/'|grep ( grep -H '// Constant: ' $CHEAT_FILES | mkconstindex grep -H '^[A-Z$][A-Z0-9_]* *=.*//' $CHEAT_FILES | mkconstindex2 - ) | sort -u | columnize 'Constants' + ) | sort -u | columnize 'Constants' 6 for f in $CHEAT_FILES ; do - egrep -H '// Function: |// Function&Module: |// Module: ' $f | mkotherindex | sort -u | columnize $f + egrep -H '// Function: |// Function&Module: |// Module: ' $f | mkotherindex | columnize "[$f]($f)" 4 echo done ) > BOSL2.wiki/CheatSheet.md diff --git a/std.scad b/std.scad index 1c2e45e..685b6ab 100644 --- a/std.scad +++ b/std.scad @@ -12,21 +12,6 @@ assert(version_num()>=20190301, "BOSL2 requires OpenSCAD version 2019.03.01 or l include include -include -include -include -include -include -include - -include -include -include -include -include -include -include - include include include @@ -36,6 +21,19 @@ include include include include +include +include +include +include +include +include +include +include +include +include +include +include +include // vim: expandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap diff --git a/version.scad b/version.scad index 6ef424a..262b891 100644 --- a/version.scad +++ b/version.scad @@ -8,7 +8,7 @@ ////////////////////////////////////////////////////////////////////// -BOSL_VERSION = [2,0,418]; +BOSL_VERSION = [2,0,419]; // Section: BOSL Library Version Functions From 8f6c2e8538db046fc6c061a3ff774a507193570e Mon Sep 17 00:00:00 2001 From: Adrian Mariano Date: Tue, 1 Sep 2020 16:42:47 -0400 Subject: [PATCH 20/23] Add submatrix_set and block_matrix and tests --- arrays.scad | 36 ++++++++++++++++++++++++++++++++++++ tests/test_arrays.scad | 25 +++++++++++++++++++++++++ 2 files changed, 61 insertions(+) diff --git a/arrays.scad b/arrays.scad index e2d3ac2..4d22be2 100644 --- a/arrays.scad +++ b/arrays.scad @@ -1199,6 +1199,42 @@ function zip(vecs, v2, v3, fit=false, fill=undef) = : [for(i=[0:1:minlen-1]) [for(v=vecs) for(x=v[i]) x] ]; +// Function: block_matrix() +// Usage: +// block_matrix([[M11, M12,...],[M21, M22,...], ... ]) +// 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. +function block_matrix(M) = + let( + bigM = [for(bigrow = M) each zip(bigrow)], + len0=len(bigM[0]), + badrows = [for(row=bigM) if (len(row)!=len0) 1] + ) + assert(badrows==[], "Inconsistent or invalid input") + bigM; + + +// Function: submatrix_set() +// Usage: submatrix_set(M,A,[m],[n]) +// 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. +function submatrix_set(M,A,m=0,n=0) = + assert(is_list(M)) + assert(is_list(A)) + 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 =n && j Date: Tue, 1 Sep 2020 17:57:31 -0400 Subject: [PATCH 21/23] Added norm_fro, quadratic_roots and pivoting to qr_factor and linear_solve. Added tests. --- math.scad | 97 ++++++++++++++++++++++++++++++++---------- tests/test_arrays.scad | 1 + tests/test_math.scad | 34 ++++++++++++++- 3 files changed, 108 insertions(+), 24 deletions(-) diff --git a/math.scad b/math.scad index 509db44..a66e1cf 100644 --- a/math.scad +++ b/math.scad @@ -680,7 +680,7 @@ function convolve(p,q) = // 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) = +function linear_solve(A,b,pivot=true) = assert(is_matrix(A), "Input should be a matrix.") let( m = len(A), @@ -688,19 +688,17 @@ function linear_solve(A,b) = ) assert(is_vector(b,m) || is_matrix(b,m),"Invalid right hand side or incompatible with the matrix") let ( - qr = mj ? 0 : ri[j] ] ] - ) [qr[0],Rzero]; + ) [qr[0],Rzero,qr[2]]; -function _qr_factor(A,Q, column, m, n) = - column >= min(m-1,n) ? [Q,A] : +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]) sum_of_squares([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=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() @@ -862,6 +869,17 @@ function is_matrix(A,m,n,square=false) = && ( !square || len(A)==len(A[0])); +// 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) = + sqrt(sum([for(entry=A) sum_of_squares(entry)])); + + // Section: Comparisons and Logic // Function: is_zero() @@ -1309,6 +1327,39 @@ function C_div(z1,z2) = // Section: Polynomials +// Function: quadratic_roots() +// Usage: +// roots = quadratic_roots(a,b,c,[real]) +// Description: +// Computes roots of the quadratic equation a*x^2+b*x+c==0, where the +// coefficients are real numbers. If real is true then returns only the +// real roots. Otherwise returns a pair of complex values. + +// https://people.csail.mit.edu/bkph/articles/Quadratics.pdf + +function quadratic_roots(a,b,c,real=false) = + real ? [for(root = quadratic_roots(a,b,c,real=false)) if (root.y==0) root.x] + : + is_undef(b) && is_undef(c) && is_vector(a,3) ? quadratic_roots(a[0],a[1],a[2]) : + assert(is_num(a) && is_num(b) && is_num(c)) + assert(a!=0 || b!=0 || c!=0, "Quadratic must have a nonzero coefficient") + a==0 && b==0 ? [] : // No solutions + a==0 ? [[-c/b,0]] : + let( + descrim = b*b-4*a*c, + sqrt_des = sqrt(abs(descrim)) + ) + descrim < 0 ? // Complex case + [[-b, sqrt_des], + [-b, -sqrt_des]]/2/a : + b<0 ? // b positive + [[2*c/(-b+sqrt_des),0], + [(-b+sqrt_des)/a/2,0]] + : // b negative + [[(-b-sqrt_des)/2/a, 0], + [2*c/(-b-sqrt_des),0]]; + + // Function: polynomial() // Usage: // polynomial(p, z) diff --git a/tests/test_arrays.scad b/tests/test_arrays.scad index 90ec305..eaec5b3 100644 --- a/tests/test_arrays.scad +++ b/tests/test_arrays.scad @@ -492,6 +492,7 @@ module test_submatrix_set() { 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(); diff --git a/tests/test_math.scad b/tests/test_math.scad index c3eb601..7b53cc6 100644 --- a/tests/test_math.scad +++ b/tests/test_math.scad @@ -781,6 +781,12 @@ 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], @@ -954,6 +960,15 @@ module test_real_roots(){ test_real_roots(); + +module test_quadratic_roots(){ + assert_approx(quadratic_roots([1,4,4]),[[-2,0],[-2,0]]); + assert_approx(quadratic_roots([1,4,4],real=true),[-2,-2]); + assert_approx(quadratic_roots([1,-5,6],real=true), [2,3]); + assert_approx(quadratic_roots([1,-5,6]), [[2,0],[3,0]]); +} +test_quadratic_roots(); + module test_qr_factor() { // Check that R is upper triangular function is_ut(R) = @@ -962,7 +977,15 @@ module test_qr_factor() { // 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); + 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) + && list_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], @@ -991,6 +1014,15 @@ module test_qr_factor() { 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(); From 399c40f7a6fd9b0d48c02e2ac9e26d629beeb7c5 Mon Sep 17 00:00:00 2001 From: Adrian Mariano Date: Tue, 1 Sep 2020 18:38:31 -0400 Subject: [PATCH 22/23] Added null_space and diagonal_matrix --- arrays.scad | 10 ++++++++++ math.scad | 17 +++++++++++++++++ tests/test_arrays.scad | 8 ++++++++ tests/test_math.scad | 23 +++++++++++++++++++++++ 4 files changed, 58 insertions(+) diff --git a/arrays.scad b/arrays.scad index 4d22be2..8f44fe8 100644 --- a/arrays.scad +++ b/arrays.scad @@ -1216,6 +1216,16 @@ function block_matrix(M) = assert(badrows==[], "Inconsistent or invalid input") bigM; +// Function: diagonal_matrix() +// Usage: +// diagonal_matrix(diag, [offdiag]) +// 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. +function diagonal_matrix(diag,offdiag=0) = + [for(i=[0:1:len(diag)-1]) [for(j=[0:len(diag)-1]) i==j?diag[i] : offdiag]]; + // Function: submatrix_set() // Usage: submatrix_set(M,A,[m],[n]) diff --git a/math.scad b/math.scad index a66e1cf..6e4be5c 100644 --- a/math.scad +++ b/math.scad @@ -712,6 +712,23 @@ function matrix_inverse(A) = assert(is_matrix(A,square=true),"Input to matrix_inverse() must be a square matrix") linear_solve(A,ident(len(A))); +// Function: null_space() +// Usage: +// 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 (is_zero(R[i],eps)) i] + ) + len(zrow)==0 + ? [] + : transpose(subindex(Q_R[0],zrow)); + // Function: qr_factor() // Usage: qr = qr_factor(A,[pivot]) diff --git a/tests/test_arrays.scad b/tests/test_arrays.scad index eaec5b3..1dfe421 100644 --- a/tests/test_arrays.scad +++ b/tests/test_arrays.scad @@ -481,6 +481,14 @@ module test_block_matrix() { 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]]; diff --git a/tests/test_math.scad b/tests/test_math.scad index 7b53cc6..9a25cec 100644 --- a/tests/test_math.scad +++ b/tests/test_math.scad @@ -969,6 +969,29 @@ 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 && is_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) = From 47a1dfaa237aea221038530ea41df955b5677d0e Mon Sep 17 00:00:00 2001 From: Adrian Mariano Date: Tue, 1 Sep 2020 18:43:53 -0400 Subject: [PATCH 23/23] doc tweak --- math.scad | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/math.scad b/math.scad index 6e4be5c..599a0c4 100644 --- a/math.scad +++ b/math.scad @@ -1350,7 +1350,9 @@ function C_div(z1,z2) = // Description: // Computes roots of the quadratic equation a*x^2+b*x+c==0, where the // coefficients are real numbers. If real is true then returns only the -// real roots. Otherwise returns a pair of complex values. +// real roots. Otherwise returns a pair of complex values. This method +// may be more reliable than the general root finder at distinguishing +// real roots from complex roots. // https://people.csail.mit.edu/bkph/articles/Quadratics.pdf