mirror of
https://github.com/BelfrySCAD/BOSL2.git
synced 2025-01-06 04:09:47 +00:00
Merge pull request #667 from adrianVmariano/master
new helix() and gaussian vectors with full covariance matrix
This commit is contained in:
Revar Desmera
GitHub
commit
c0ac2d7d51
4 changed files with 125 additions and 68 deletions
61
drawing.scad
61
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) =
|
||||
|
|
|
|||
112
math.scad
112
math.scad
|
|
@ -498,83 +498,75 @@ 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) :
|
||||
assert(is_matrix_symmetric(cov),"Supplied covariance matrix is not symmetric")
|
||||
let(
|
||||
L = cholesky(cov)
|
||||
)
|
||||
assert(is_def(L), "Supplied covariance matrix is not positive definite")
|
||||
move(mean,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 +1082,40 @@ 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")
|
||||
_cholesky(A,ident(len(A)), len(A));
|
||||
|
||||
function _cholesky(A,L,n) =
|
||||
A[0][0]<0 ? undef : // Matrix not positive definite
|
||||
len(A) == 1 ? submatrix_set(L,[[sqrt(A[0][0])]], n-1,n-1):
|
||||
let(
|
||||
i = n+1-len(A)
|
||||
)
|
||||
let(
|
||||
sqrtAii = sqrt(A[0][0]),
|
||||
Lnext = [for(j=[0:n-1])
|
||||
[for(k=[0:n-1])
|
||||
j<i-1 || k<i-1 ? (j==k ? 1 : 0)
|
||||
: j==i-1 && k==i-1 ? sqrtAii
|
||||
: j==i-1 ? 0
|
||||
: k==i-1 ? A[j-(i-1)][0]/sqrtAii
|
||||
: j==k ? 1 : 0]],
|
||||
Anext = submatrix(A,[1:n-1], [1:n-1]) - outer_product(list_tail(A[0]), list_tail(A[0]))/A[0][0]
|
||||
)
|
||||
_cholesky(Anext,L*Lnext,n);
|
||||
|
||||
|
||||
// Function: det2()
|
||||
// Usage:
|
||||
// d = det2(M);
|
||||
|
|
|
|||
|
|
@ -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(
|
||||
|
|
|
|||
|
|
@ -360,14 +360,25 @@ test_rand_int();
|
|||
|
||||
|
||||
module test_gaussian_rands() {
|
||||
nums1 = gaussian_rands(0,10,1000,seed=2132);
|
||||
nums2 = gaussian_rands(0,10,1000,seed=2130);
|
||||
nums3 = gaussian_rands(0,10,1000,seed=2132);
|
||||
nums1 = gaussian_rands(1000,0,10,seed=2132);
|
||||
nums2 = gaussian_rands(1000,0,10,seed=2130);
|
||||
nums3 = gaussian_rands(1000,0,10,seed=2132);
|
||||
assert_equal(len(nums1), 1000);
|
||||
assert_equal(len(nums2), 1000);
|
||||
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();
|
||||
|
||||
|
|
|
|||
Loading…
Reference in a new issue