From 71b22e5850d529518e5f54468e47bf76ffacfb51 Mon Sep 17 00:00:00 2001 From: Adrian Mariano Date: Fri, 1 Oct 2021 21:27:52 -0400 Subject: [PATCH 1/3] helix rewrite make skin() take regions --- drawing.scad | 61 +++++++++++++++++++++++++++++++++------------------- skin.scad | 3 +++ 2 files changed, 42 insertions(+), 22 deletions(-) diff --git a/drawing.scad b/drawing.scad index 1ff77d9..b4db4fa 100644 --- a/drawing.scad +++ b/drawing.scad @@ -688,32 +688,49 @@ module arc(N, r, angle, d, cp, points, width, thickness, start, wedge=false) // Function: helix() -// Description: -// Returns a 3D helical path. // Usage: -// helix(turns, h, n, r|d, [cp], [scale]); +// helix([l|h], [turns], [angle], r|r1|r2, d|d1|d2) +// Description: +// Returns a 3D helical path on a cone, including the degerate case of flat spirals. +// You can specify start and end radii. You can give the length, the helix angle, or the number of turns: two +// of these three parameters define the helix. For a flat helix you must give length 0 and a turn count. +// Helix will be right handed if turns is positive and left handed if it is negative. +// The angle is calculateld based on the radius at the base of the helix. // Arguments: -// h = Height of spiral. -// turns = Number of turns in spiral. -// n = Number of spiral sides. -// r = Radius of spiral. -// d = Radius of spiral. -// cp = Centerpoint of spiral. Default: `[0,0]` -// scale = [X,Y] scaling factors for each axis. Default: `[1,1]` +// h|l = Height/length of helix, zero for a flat spiral +// --- +// turns = Number of turns in helix, positive for right handed +// angle = helix angle +// r = Radius of helix +// r1 = Radius of bottom of helix +// r2 = Radius of top of helix +// d = Diameter of helix +// d1 = Diameter of bottom of helix +// d2 = Diameter of top of helix // Example(3D): -// trace_path(helix(turns=2.5, h=100, n=24, r=50), N=1, showpts=true); -function helix(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 - ) [ - for (i=[0:1:cnt]) [ - rr * cos(i*360/n) * scale.x + cp.x, - rr * sin(i*360/n) * scale.y + cp.y, - i*dz - ] - ]; +// trace_path(helix(turns=2.5, h=100, r=50), N=1, showpts=true); +// Example(3D): Helix that turns the other way +// trace_path(helix(turns=-2.5, h=100, r=50), N=1, showpts=true); +// Example(3D): Flat helix (note points are still 3d) +// stroke(helix(h=0,r1=50,r2=25,l=0, turns=4)); +function helix(l,h,turns,angle, r, r1, r2, d, d1, d2)= + let( + r1=get_radius(r=r,r1=r1,d=d,d1=d1,dflt=1), + r2=get_radius(r=r,r1=r2,d=d,d1=d2,dflt=1), + length = first_defined([l,h]) + ) + assert(num_defined([length,turns,angle])==2,"Must define exactly two of l/h, turns, and angle") + assert(is_undef(angle) || length!=0, "Cannot give length 0 with an angle") + let( + // length advances dz for each turn + dz = is_def(angle) && length!=0 ? 2*PI*r1*tan(angle) : length/abs(turns), + maxtheta = is_def(turns) ? 360*turns : 360*length/dz, + N = segs(max(r1,r2)) + ) + [for(theta=lerpn(0,maxtheta, max(3,ceil(abs(maxtheta)*N/360)))) + let(R=lerp(r1,r2,theta/maxtheta)) + [R*cos(theta), R*sin(theta), abs(theta)/360 * dz]]; function _normal_segment(p1,p2) = diff --git a/skin.scad b/skin.scad index 65e0df0..7b759c4 100644 --- a/skin.scad +++ b/skin.scad @@ -395,6 +395,9 @@ module skin(profiles, slices, refine=1, method="direct", sampling, caps, closed= function skin(profiles, slices, refine=1, method="direct", sampling, caps, closed=false, z, style="min_edge") = assert(is_def(slices),"The slices argument must be specified.") assert(is_list(profiles) && len(profiles)>1, "Must provide at least two profiles") + let( + profiles = [for(p=profiles) if (is_region(p) && len(p)==1) p[0] else p] + ) let( bad = [for(i=idx(profiles)) if (!(is_path(profiles[i]) && len(profiles[i])>2)) i]) assert(len(bad)==0, str("Profiles ",bad," are not a paths or have length less than 3")) let( From a5ae4879be99c604648a58c1b54ad10283f544d5 Mon Sep 17 00:00:00 2001 From: Adrian Mariano Date: Fri, 1 Oct 2021 23:37:06 -0400 Subject: [PATCH 2/3] add full covariance gaussian random vectors add cholesky fatorization (needed for above, also useful for solving symmetric linear systems.) --- math.scad | 111 ++++++++++++++++++++++++++----------------- tests/test_math.scad | 6 +-- 2 files changed, 71 insertions(+), 46 deletions(-) diff --git a/math.scad b/math.scad index 4765983..3a9a818 100644 --- a/math.scad +++ b/math.scad @@ -498,83 +498,73 @@ function rand_int(minval, maxval, N, seed=undef) = // Function: random_points() // Usage: -// points = random_points(n, dim, scale, [seed]); +// points = random_points([N], [dim], [scale], [seed]); // See Also: random_polygon(), gaussian_random_points(), spherical_random_points() // Topics: Random, Points // Description: -// Generate `n` uniform random points of dimension `dim` with data ranging from -scale to +scale. +// Generate `N` uniform random points of dimension `dim` with data ranging from -scale to +scale. // The `scale` may be a number, in which case the random data lies in a cube, // or a vector with dimension `dim`, in which case each dimension has its own scale. // Arguments: -// n = number of points to generate. +// N = number of points to generate. Default: 1 // dim = dimension of the points. Default: 2 // scale = the scale of the point coordinates. Default: 1 // seed = an optional seed for the random generation. -function random_points(n, dim=2, scale=1, seed) = - assert( is_int(n) && n>=0, "The number of points should be a non-negative integer.") +function random_points(N, dim=2, scale=1, seed) = + assert( is_int(N) && N>=0, "The number of points should be a non-negative integer.") assert( is_int(dim) && dim>=1, "The point dimensions should be an integer greater than 1.") assert( is_finite(scale) || is_vector(scale,dim), "The scale should be a number or a vector with length equal to d.") let( rnds = is_undef(seed) - ? rands(-1,1,n*dim) - : rands(-1,1,n*dim, seed) ) + ? rands(-1,1,N*dim) + : rands(-1,1,N*dim, seed) ) is_num(scale) - ? scale*[for(i=[0:1:n-1]) [for(j=[0:dim-1]) rnds[i*dim+j] ] ] - : [for(i=[0:1:n-1]) [for(j=[0:dim-1]) scale[j]*rnds[i*dim+j] ] ]; + ? scale*[for(i=[0:1:N-1]) [for(j=[0:dim-1]) rnds[i*dim+j] ] ] + : [for(i=[0:1:N-1]) [for(j=[0:dim-1]) scale[j]*rnds[i*dim+j] ] ]; // Function: gaussian_rands() // Usage: -// arr = gaussian_rands(mean, stddev, [N], [seed]); +// arr = gaussian_rands([N],[mean], [cov], [seed]); // Description: -// Returns a random number with a gaussian/normal distribution. +// Returns a random number or vector with a Gaussian/normal distribution. // Arguments: -// mean = The average random number returned. -// stddev = The standard deviation of the numbers to be returned. -// N = Number of random numbers to return. Default: 1 +// N = the number of points to return. Default: 1 +// mean = The average of the random value (a number or vector). Default: 0 +// cov = covariance matrix of the random numbers, or variance in the 1D case. Default: 1 // seed = If given, sets the random number seed. -function gaussian_rands(mean, stddev, N=1, seed=undef) = - assert( is_finite(mean+stddev+N) && (is_undef(seed) || is_finite(seed) ), "Input must be finite numbers.") - let(nums = is_undef(seed)? rands(0,1,N*2) : rands(0,1,N*2,seed)) - [for (i = count(N,0,2)) mean + stddev*sqrt(-2*ln(nums[i]))*cos(360*nums[i+1])]; - - -// Function: gaussian_random_points() -// Usage: -// points = gaussian_random_points(n, dim, mean, stddev, [seed]); -// See Also: random_polygon(), random_points(), spherical_random_points() -// Topics: Random, Points -// Description: -// Generate `n` random points of dimension `dim` with coordinates absolute value less than `scale`. -// The gaussian distribution of all the coordinates of the points will have a mean `mean` and -// standard deviation `stddev` -// Arguments: -// n = number of points to generate. -// dim = dimension of the points. Default: 2 -// mean = the gaussian mean of the point coordinates. Default: 0 -// stddev = the gaussian standard deviation of the point coordinates. Default: 0 -// seed = an optional seed for the random generation. -function gaussian_random_points(n, dim=2, mean=0, stddev=1, seed) = - assert( is_int(n) && n>=0, "The number of points should be a non-negative integer.") - assert( is_int(dim) && dim>=1, "The point dimensions should be an integer greater than 1.") - let( rnds = gaussian_rands(mean, stddev, n*dim, seed=seed) ) - [for(i=[0:1:n-1]) [for(j=[0:dim-1]) rnds[i*dim+j] ] ]; +function gaussian_rands(N=1, mean=0, cov=1, seed=undef) = + assert(is_num(mean) || is_vector(mean)) + let( + dim = is_num(mean) ? 1 : len(mean) + ) + assert((dim==1 && is_num(cov)) || is_matrix(cov,dim,dim),"mean and covariance matrix not compatible") + assert(is_undef(seed) || is_finite(seed)) + let( + nums = is_undef(seed)? rands(0,1,dim*N*2) : rands(0,1,dim*N*2,seed), + rdata = [for (i = count(dim*N,0,2)) sqrt(-2*ln(nums[i]))*cos(360*nums[i+1])] + ) + dim==1 ? add_scalar(sqrt(cov)*rdata,mean) : + let( + L = cholesky(cov) + ) + array_group(rdata,dim)*transpose(L); // Function: spherical_random_points() // Usage: -// points = spherical_random_points(n, radius, [seed]); +// points = spherical_random_points([N], [radius], [seed]); // See Also: random_polygon(), random_points(), gaussian_random_points() // Topics: Random, Points // Description: // Generate `n` 3D uniformly distributed random points lying on a sphere centered at the origin with radius equal to `radius`. // Arguments: -// n = number of points to generate. +// n = number of points to generate. Default: 1 // radius = the sphere radius. Default: 1 // seed = an optional seed for the random generation. // See https://mathworld.wolfram.com/SpherePointPicking.html -function spherical_random_points(n, radius=1, seed) = +function spherical_random_points(N=1, radius=1, seed) = assert( is_int(n) && n>=1, "The number of points should be an integer greater than zero.") assert( is_num(radius) && radius>0, "The radius should be a non-negative number.") let( theta = is_undef(seed) @@ -1090,6 +1080,41 @@ function _back_substitute(R, b, x=[]) = _back_substitute(R, b, concat([newvalue],x)); + +// Function: cholesky() +// Usage: +// L = cholesky(A); +// Description: +// Compute the cholesky factor, L, of the symmetric positive definite matrix A. +// The matrix L is lower triangular and L * transpose(L) = A. If the A is +// not symmetric then an error is displayed. If the matrix is symmetric but +// not positive definite then undef is returned. +function cholesky(A) = + assert(is_matrix(A,square=true),"A must be a square matrix") + assert(is_matrix_symmetric(A),"Cholesky factorization requires a symmetric matrix") + echo(A=A,len=len(A)) + _cholesky(A,ident(len(A)), len(A)); + +function _cholesky(A,L,n) = let(ffee=echo(insideA=A,L,n)) + A[0][0]<0 ? undef : // Matrix not positive definite + len(A) == 1 ? submatrix_set(L,[[sqrt(A[0][0])]], n-1,n-1): + let( + i = n+1-len(A),ff=echo(i=i,lenA=len(A)) + ) + let( + sqrtAii = sqrt(A[0][0]), + Lnext = [for(j=[0:n-1]) + [for(k=[0:n-1]) + j Date: Fri, 1 Oct 2021 23:57:07 -0400 Subject: [PATCH 3/3] fix tests and bugs revealed by new tests --- math.scad | 9 +++++---- tests/test_math.scad | 11 +++++++++++ 2 files changed, 16 insertions(+), 4 deletions(-) diff --git a/math.scad b/math.scad index 3a9a818..bb5c67d 100644 --- a/math.scad +++ b/math.scad @@ -545,10 +545,12 @@ function gaussian_rands(N=1, mean=0, cov=1, seed=undef) = rdata = [for (i = count(dim*N,0,2)) sqrt(-2*ln(nums[i]))*cos(360*nums[i+1])] ) dim==1 ? add_scalar(sqrt(cov)*rdata,mean) : + assert(is_matrix_symmetric(cov),"Supplied covariance matrix is not symmetric") let( L = cholesky(cov) ) - array_group(rdata,dim)*transpose(L); + assert(is_def(L), "Supplied covariance matrix is not positive definite") + move(mean,array_group(rdata,dim)*transpose(L)); // Function: spherical_random_points() @@ -1092,14 +1094,13 @@ function _back_substitute(R, b, x=[]) = function cholesky(A) = assert(is_matrix(A,square=true),"A must be a square matrix") assert(is_matrix_symmetric(A),"Cholesky factorization requires a symmetric matrix") - echo(A=A,len=len(A)) _cholesky(A,ident(len(A)), len(A)); -function _cholesky(A,L,n) = let(ffee=echo(insideA=A,L,n)) +function _cholesky(A,L,n) = A[0][0]<0 ? undef : // Matrix not positive definite len(A) == 1 ? submatrix_set(L,[[sqrt(A[0][0])]], n-1,n-1): let( - i = n+1-len(A),ff=echo(i=i,lenA=len(A)) + i = n+1-len(A) ) let( sqrtAii = sqrt(A[0][0]), diff --git a/tests/test_math.scad b/tests/test_math.scad index 50adce3..a5a5843 100644 --- a/tests/test_math.scad +++ b/tests/test_math.scad @@ -368,6 +368,17 @@ module test_gaussian_rands() { assert_equal(len(nums3), 1000); assert_equal(nums1, nums3); assert(nums1!=nums2); + + R = [[4,2],[2,17]]; + data = gaussian_rands(100000,[0,0],R,seed=49); + assert(approx(mean(data), [0,0], eps=1e-2)); + assert(approx(transpose(data)*data/len(data), R, eps=2e-2)); + + R2 = [[4,2,-1],[2,17,4],[-1,4,11]]; + data3 = gaussian_rands(100000,[1,2,3],R2,seed=97); + assert(approx(mean(data3),[1,2,3], eps=1e-2)); + cdata = move(-mean(data3),data3); + assert(approx(transpose(cdata)*cdata/len(cdata),R2,eps=.1)); } test_gaussian_rands();